pkg/profiles/cost_profiles.F

24462d2fa8 Patr*0001 #include "PROFILES_OPTIONS.h"
36c7a91797 Gael*0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
6a770e0a24 Patr*0005
                0006 C     o==========================================================o
98bf704dd5 Jean*0007 C     | subroutine cost_profiles                                 |
                0008 C     | o computes the cost for netcdf profiles data             |
6a770e0a24 Patr*0009 C     | started: Gael Forget 15-March-2006                       |
                0010 C     o==========================================================o
                0011
65754df434 Mart*0012       SUBROUTINE cost_profiles( myIter, myTime, myThid )
6a770e0a24 Patr*0013
                0014       IMPLICIT NONE
                0015
                0016 C     ======== Global data ============================
                0017 #include "SIZE.h"
275a56dc21 Gael*0018 #include "EEPARAMS.h"
                0019 #include "PARAMS.h"
6a770e0a24 Patr*0020 #include "GRID.h"
                0021 #include "DYNVARS.h"
24462d2fa8 Patr*0022 #ifdef ALLOW_CAL
6328b73337 Gael*0023 # include "cal.h"
24462d2fa8 Patr*0024 #endif
                0025 #ifdef ALLOW_PROFILES
6328b73337 Gael*0026 # include "PROFILES_SIZE.h"
                0027 # include "profiles.h"
                0028 # include "netcdf.inc"
6e4c90fea3 Patr*0029 #endif
36c7a91797 Gael*0030 #ifdef ALLOW_CTRL
65754df434 Mart*0031 # include "OPTIMCYCLE.h"
36c7a91797 Gael*0032 #endif
6e4c90fea3 Patr*0033
6a770e0a24 Patr*0034 c     == routine arguments ==
65754df434 Mart*0035       integer myIter
                0036       _RL     myTime
                0037       integer myThid
6e4c90fea3 Patr*0038
d28c90138c Patr*0039 #ifdef ALLOW_PROFILES
6e4c90fea3 Patr*0040
6a770e0a24 Patr*0041 C     ========= Local variables =======================
                0042       integer K,num_file,num_var,prof_num
b2a948f981 Gael*0043       integer bi,bj,iG,jG,fid
24405677c3 Gael*0044       _RL prof_traj1D(NLEVELMAX), prof_traj1D_mean(NLEVELMAX)
6a770e0a24 Patr*0045       _RL prof_data1D(NLEVELMAX), prof_weights1D(NLEVELMAX)
36c7a91797 Gael*0046 #ifndef ALLOW_CTRL
                0047       integer optimcycle
                0048 #endif
65754df434 Mart*0049       CHARACTER*(MAX_LEN_MBUF) msgbuf
                0050       CHARACTER*(MAX_LEN_FNAM) profilesfile, fnameequinc
b2a948f981 Gael*0051       integer IL, JL, err
                0052
3c8dcfdea9 Gael*0053       _RL  objf_prof_tile (nSx,nSy)
                0054       _RL  objf_prof_glo
                0055       _RL  num_prof_tile (nSx,nSy)
                0056       _RL  num_prof_glo
                0057
6b2230d510 Ou W*0058 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0059       integer iavgbin,ikzz
                0060       integer itmp
                0061       integer k2, ix9, iy9, ktmp
                0062       integer cunit
65754df434 Mart*0063       CHARACTER*(MAX_LEN_FNAM) cfile
6b2230d510 Ou W*0064
                0065       _RL prof_data1D_mean(NLEVELMAX)
                0066       _RL prof_count1D(NLEVELMAX)
                0067       _RL prof_weights1D_mean(NLEVELMAX)
                0068       _RL recip_profiles_mean_indsamples(NVARMAX)
                0069
                0070       _RL tmpr6, tmpr7, tmpr8, tmpr9
                0071       Real*4 tmp99(NAVGBINMAX)
                0072       _RL tmp11, tmp12, tmp_recip_count
                0073       LOGICAL doglbsum
                0074
                0075       _RL  objf_prof_mean_tile (nSx,nSy)
                0076       _RL  objf_prof_mean_glo
                0077       _RL  num_prof_mean_tile (nSx,nSy)
                0078       _RL  num_prof_mean_glo
                0079 #endif
                0080
b2a948f981 Gael*0081 C     !FUNCTIONS
                0082       INTEGER  ILNBLNK
                0083       EXTERNAL ILNBLNK
                0084
6a770e0a24 Patr*0085 c     == end of interface ==
98bf704dd5 Jean*0086
36c7a91797 Gael*0087 #ifndef ALLOW_CTRL
                0088       optimcycle = 0
                0089 #endif
                0090
3c8dcfdea9 Gael*0091       write(msgbuf,'(a)') ' '
                0092       call print_message( msgbuf,
65754df434 Mart*0093      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0094       write(msgbuf,'(a)') '== cost_profiles: begin =='
                0095       call print_message( msgbuf,
65754df434 Mart*0096      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0097
65754df434 Mart*0098         _BEGIN_MASTER( myThid )
6cc0cb518f Gael*0099
6b2230d510 Ou W*0100 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0101       NAVGBIN = 0
                0102 C initialize
                0103       DO iavgbin = 1, NAVGBINMAX
                0104        avgbinglbsum(iavgbin) = 0
                0105        DO ikzz = 1, NLEVELCOMBMAX
                0106         DO num_var=1,NVARMAX
                0107            prof_traj1D_all_mean(iavgbin,ikzz,num_var)
                0108      &      = 0. _d 0
                0109            prof_data1D_all_mean(iavgbin,ikzz,num_var)
                0110      &      = 0. _d 0
                0111            prof_weights1D_all_mean(iavgbin,ikzz,num_var)
                0112      &      = 0. _d 0
                0113            prof_count1D_all_mean(iavgbin,ikzz,num_var)
                0114      &      = 0. _d 0
                0115         ENDDO
                0116        ENDDO
                0117       ENDDO
                0118
                0119       DO num_var=1,NVARMAX
                0120          recip_profiles_mean_indsamples(num_var) = 0. _d 0
                0121          IF(profiles_mean_indsamples(num_var).GT. 0. _d 0) THEN
                0122           recip_profiles_mean_indsamples(num_var) = 1. _d 0 /
                0123      &     profiles_mean_indsamples(num_var)
                0124          ENDIF
                0125       ENDDO
                0126
6cc0cb518f Gael*0127       DO bj=1,nSy
                0128        DO bi=1,nSx
5058c82dba Patr*0129         do num_file=1,NFILESPROFMAX
                0130          fid=fiddata(num_file,bi,bj)
b2a948f981 Gael*0131
                0132          if ( (ProfNo(num_file,bi,bj).GT.0).AND.
f0e4bffe35 Gael*0133      &        (profilesDoNcOutput) ) then
b2a948f981 Gael*0134 c need to close the file so that the data is not lost when run finishes
                0135            err = NF_CLOSE(fidforward(num_file,bi,bj))
                0136 c then re-open it to compute cost function
                0137            iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
                0138            jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
                0139            IL  = ILNBLNK( profilesfiles(num_file) )
0e6a4460e5 Ivan*0140            write(profilesfile,'(a)')
b2a948f981 Gael*0141      &     profilesfiles(num_file)(1:IL)
                0142            IL  = ILNBLNK( profilesfile )
                0143            JL  = ILNBLNK( profilesDir )
0e6a4460e5 Ivan*0144            write(fnameequinc,'(3a,2(i3.3,a))')
b2a948f981 Gael*0145      &     profilesDir(1:JL),profilesfile(1:IL),'.',iG,'.',jG,'.equi.nc'
                0146 c
                0147            err = NF_OPEN(fnameequinc,NF_NOWRITE,
                0148      &     fidforward(num_file,bi,bj))
                0149          endif
                0150
6b2230d510 Ou W*0151 C find the vertical indices
                0152          do K=1,NLEVELMAX
                0153           prof_lev_comb(k,num_file,bi,bj) = -999
                0154           if(K.LE.ProfDepthNo(num_file,bi,bj)) then
                0155            do k2 = 1, NLEVELCOMB
                0156              if(prof_depth(num_file, k,bi,bj).EQ.
                0157      &          prof_depth_comb(k2,bi,bj).AND.
                0158      &          prof_depth_comb(k2,bi,bj).GE.0. _d 0.AND.
                0159      &          prof_lev_comb(k,num_file,bi,bj).EQ.-999) then
                0160               prof_lev_comb(k,num_file,bi,bj) = k2
                0161              endif
                0162            enddo
                0163           endif
                0164          enddo
                0165
                0166          do num_var=1,NVARMAX
                0167           if (vec_quantities(num_file,num_var,bi,bj).EQV..TRUE.) then
                0168            do prof_num=1,NOBSGLOB
                0169             if (prof_num.LE.ProfNo(num_file,bi,bj)) then
65754df434 Mart*0170
6b2230d510 Ou W*0171               do K=1,NLEVELMAX
                0172                 prof_traj1D(k)=0.
                0173 C             prof_traj1D_mean(k)=0.
                0174 C             prof_mask1D_cur(k,bi,bj)=0.
                0175                prof_data1D(k)=0.
                0176 C             prof_data1D_mean(k)=0.
                0177                prof_weights1D(k)=0.
                0178               enddo
                0179               ix9 = prof_interp_i(num_file,prof_num,1,bi,bj)
                0180               iy9 = prof_interp_j(num_file,prof_num,1,bi,bj)
65754df434 Mart*0181
6b2230d510 Ou W*0182               if(prof_ind_avgbin(num_file,prof_num,bi,bj).GT.NAVGBIN)
                0183      &         NAVGBIN = prof_ind_avgbin(num_file,prof_num,bi,bj)
                0184
                0185               if(ix9 .GE. 0. _d 0 .AND. iy9 .GE. 0. _d 0) then
                0186                itmp = prof_ind_avgbin(num_file,prof_num,bi,bj)
                0187                if(avgbinglbsum(itmp).EQ.0)
                0188      &          avgbinglbsum(itmp) = 1
                0189
                0190                call active_read_profile(num_file,
                0191      &          ProfDepthNo(num_file,bi,bj),prof_traj1D,num_var,
65754df434 Mart*0192      &          prof_num,.false.,optimcycle,bi,bj,myThid,
6b2230d510 Ou W*0193      &          profiles_dummy(num_file,num_var,bi,bj))
                0194
                0195                call profiles_readvector(num_file,num_var,
                0196      &          prof_ind_glob(num_file,prof_num,bi,bj),
                0197      &          ProfDepthNo(num_file,bi,bj),prof_data1D,bi,bj,myThid)
65754df434 Mart*0198
6b2230d510 Ou W*0199                call profiles_readvector(num_file,-num_var,
                0200      &          prof_ind_glob(num_file,prof_num,bi,bj),
                0201      &          ProfDepthNo(num_file,bi,bj),
                0202      &          prof_weights1D,bi,bj,myThid)
                0203
                0204                do K=1,ProfDepthNo(num_file,bi,bj)
                0205                 if (prof_weights1D(K).GT.0. _d 0
                0206      &           .AND. prof_mask1D_cur(K,bi,bj).NE. 0. _d 0
                0207      &             ) then
                0208                  prof_traj1D_all_mean(itmp,
                0209      &            prof_lev_comb(k,num_file,bi,bj),num_var)
                0210      &            = prof_traj1D_all_mean(itmp,
                0211      &               prof_lev_comb(k,num_file,bi,bj), num_var)
                0212      &            + prof_traj1D(k)
65754df434 Mart*0213
6b2230d510 Ou W*0214                  prof_data1D_all_mean(itmp,
                0215      &            prof_lev_comb(k,num_file,bi,bj), num_var)
                0216      &            = prof_data1D_all_mean(itmp,
                0217      &               prof_lev_comb(k,num_file,bi,bj), num_var)
                0218      &            + prof_data1D(k)
                0219
                0220                  prof_weights1D_all_mean(itmp,
                0221      &            prof_lev_comb(k,num_file,bi,bj), num_var)
                0222      &            = prof_weights1D_all_mean(itmp,
                0223      &               prof_lev_comb(k,num_file,bi,bj), num_var)
                0224      &            + 1. _d 0 /prof_weights1D(k)
                0225
                0226                  prof_count1D_all_mean(itmp,
                0227      &            prof_lev_comb(k,num_file,bi,bj), num_var)
                0228      &            = prof_count1D_all_mean(itmp,
                0229      &               prof_lev_comb(k,num_file,bi,bj), num_var)
                0230      &            + 1. _d 0
                0231                 endif
                0232                enddo !do K=1,ProfDepthNo
                0233               endif !      if(ix9 .GE. 0. _d 0 .AND. iy9 .GE. 0. _d 0) then
                0234
                0235             endif !if (prof_num.LE.ProfNo(num_file,bi,bj)) then
                0236            enddo !do prof_num=..
                0237           endif !if (vec_quantities(num_file,num_var,bi,bj).EQV..TRUE.) then
                0238          enddo !do num_var...
                0239
                0240         enddo !do num_file=1,NFILESPROFMAX
                0241
                0242        ENDDO !DO bj
                0243        ENDDO !DO bj
                0244
                0245        NAVGBINRL = NAVGBIN
                0246        _GLOBAL_MAX_RL( NAVGBINRL, myThid )
                0247        NAVGBIN = NAVGBINRL
                0248        DO iavgbin = 1, NAVGBIN
                0249           tmpr6 = avgbinglbsum(iavgbin)
                0250           _GLOBAL_SUM_RL (tmpr6, myThid)
                0251           if(tmpr6.GT.1.1) avgbinglbsum(iavgbin) = tmpr6
                0252        ENDDO
                0253
                0254 C accumulate globally
                0255        DO num_var=1,NVARMAX
                0256         doglbsum = .FALSE.
                0257         DO bj=1,nSy
                0258          DO bi=1,nSx
                0259           do num_file=1,NFILESPROFMAX
                0260             if (vec_quantities(num_file,num_var,bi,bj).EQV..TRUE.) then
                0261               doglbsum = .TRUE.
                0262             endif
                0263           enddo
                0264          ENDDO
                0265         ENDDO
                0266
                0267         if(doglbsum) then
                0268          DO iavgbin = 1, NAVGBIN
                0269            DO ikzz = 1, NLEVELCOMB
                0270             tmpr6 = prof_count1D_all_mean(iavgbin,ikzz,num_var)
                0271             _GLOBAL_SUM_RL (tmpr6, myThid)
                0272             prof_count1D_all_mean(iavgbin,ikzz,num_var) = tmpr6
                0273
                0274             tmpr9 = prof_weights1D_all_mean(iavgbin,ikzz,num_var)
                0275             _GLOBAL_SUM_RL (tmpr9, myThid)
                0276             prof_weights1D_all_mean(iavgbin,ikzz,num_var) = tmpr9
                0277
                0278             tmpr7 = prof_traj1D_all_mean(iavgbin,ikzz,num_var)
                0279             _GLOBAL_SUM_RL (tmpr7, myThid)
                0280             prof_traj1D_all_mean(iavgbin,ikzz,num_var) = tmpr7
                0281
                0282             tmpr8 = prof_data1D_all_mean(iavgbin,ikzz,num_var)
                0283             _GLOBAL_SUM_RL (tmpr8, myThid)
                0284             prof_data1D_all_mean(iavgbin,ikzz,num_var) = tmpr8
                0285            ENDDO
                0286          ENDDO
                0287         endif
                0288        ENDDO
                0289
                0290 C Now do the averaging
                0291        DO iavgbin = 1, NAVGBIN
                0292         DO ikzz = 1, NLEVELCOMB
                0293          DO num_var=1,NVARMAX
                0294             tmp_recip_count = 0. _d 0
                0295             IF(prof_count1D_all_mean(iavgbin,ikzz,num_var).GT.0)THEN
                0296              tmp_recip_count = 1. _d 0 /
                0297      &          prof_count1D_all_mean(iavgbin,ikzz,num_var)
                0298              prof_traj1D_all_mean(iavgbin,ikzz,num_var)
                0299      &        = prof_traj1D_all_mean(iavgbin,ikzz,num_var)*
65754df434 Mart*0300      &          tmp_recip_count
6b2230d510 Ou W*0301              prof_data1D_all_mean(iavgbin,ikzz,num_var)
                0302      &        = prof_data1D_all_mean(iavgbin,ikzz,num_var)*
65754df434 Mart*0303      &          tmp_recip_count
6b2230d510 Ou W*0304              prof_weights1D_all_mean(iavgbin,ikzz,num_var)
                0305      &        = prof_weights1D_all_mean(iavgbin,ikzz,num_var)*
65754df434 Mart*0306      &          tmp_recip_count
6b2230d510 Ou W*0307             ENDIF
                0308          ENDDO
                0309         ENDDO
                0310        ENDDO
                0311
                0312        DO iavgbin = 1, NAVGBIN
                0313         DO ikzz = 1, NLEVELCOMB
                0314          DO num_var=1,NVARMAX
                0315             IF(prof_count1D_all_mean(iavgbin,ikzz,num_var).GT.0)THEN
                0316 C Assuming each averaging bin has a maximum of 9 independent measurements.
                0317              tmp11 = prof_weights1D_all_mean(iavgbin,ikzz,num_var)
                0318      &             / prof_count1D_all_mean(iavgbin,ikzz,num_var)
                0319              tmp12 = prof_weights1D_all_mean(iavgbin,ikzz,num_var)
                0320      &             * recip_profiles_mean_indsamples(num_var)
                0321              prof_weights1D_all_mean(iavgbin,ikzz,num_var)
                0322      &        = max(tmp11, tmp12)
                0323
                0324 C note prof_weights1D_all_mean is still sigam^2. Need to convert to weight
                0325             if(prof_weights1D_all_mean(iavgbin,ikzz,num_var).NE.0. _d 0)
                0326      &        prof_weights1D_all_mean(iavgbin,ikzz,num_var) =
                0327      &         1. _d 0 /prof_weights1D_all_mean(iavgbin,ikzz,num_var)
                0328             ENDIF
                0329          ENDDO
                0330         ENDDO
                0331        ENDDO
                0332
                0333        IF ( myProcId .eq. 0 ) THEN
                0334
                0335         DO num_var=1,NVARMAX
                0336          iL = ILNBLNK( prof_names(1,num_var) )
                0337          write(cfile,'(2a)') prof_names(1,num_var)(1:iL),
                0338      &   '_data_mean.data'
65754df434 Mart*0339          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0340          open( cunit, file   = cfile,
                0341      &        status = 'unknown',
                0342      &        access  = 'direct',
                0343      &        recl = NAVGBINMAX*4)
65754df434 Mart*0344
6b2230d510 Ou W*0345          DO ikzz = 1, NLEVELCOMB
                0346           tmp99(1:NAVGBINMAX)=
                0347      &      prof_data1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0348           write(cunit,rec=ikzz) tmp99
                0349          ENDDO
                0350          close ( cunit )
                0351
                0352          write(cfile,'(2a)')prof_names(1,num_var)(1:iL),
                0353      &    '_model_mean.data'
65754df434 Mart*0354          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0355          open( cunit, file   = cfile,
                0356      &         status = 'unknown',
                0357 C    &         form   = 'unformatted',
                0358      &         access  = 'direct',
                0359      &         recl = NAVGBINMAX*4)
                0360 C    &         access  = 'sequential'   )
                0361
                0362          DO ikzz = 1, NLEVELCOMB
                0363           tmp99(1:NAVGBINMAX)=
                0364      &      prof_traj1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0365           write(cunit,rec=ikzz) tmp99
                0366          ENDDO
                0367          close ( cunit )
                0368
                0369          write(cfile,'(2a)')
                0370      &     prof_names(1,num_var)(1:iL),'_weight_mean.data'
65754df434 Mart*0371          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0372          open( cunit, file   = cfile,
                0373      &         status = 'unknown',
                0374      &         access  = 'direct',
                0375      &         recl = NAVGBINMAX*4)
                0376
                0377          DO ikzz = 1, NLEVELCOMB
                0378           tmp99(1:NAVGBINMAX)=
                0379      &      prof_weights1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0380           write(cunit,rec=ikzz) tmp99
                0381          ENDDO
                0382          close ( cunit )
                0383
                0384          write(cfile,'(2a)')prof_names(1,num_var)(1:iL),
                0385      &    '_count_mean.data'
65754df434 Mart*0386          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0387          open( cunit, file   = cfile,
                0388      &         status = 'unknown',
                0389      &         access  = 'direct',
                0390      &         recl = NAVGBINMAX*4)
65754df434 Mart*0391
6b2230d510 Ou W*0392          DO ikzz = 1, NLEVELCOMB
                0393           tmp99(1:NAVGBINMAX)=
                0394      &     prof_count1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0395           write(cunit,rec=ikzz) tmp99
                0396          ENDDO
                0397          close ( cunit )
                0398         ENDDO ! DO num_var=1,NVARMAX
                0399        ENDIF ! IF ( myProcId .eq. 0 ) THEN
                0400 #endif
                0401
                0402        DO bj=1,nSy
                0403         DO bi=1,nSx
                0404
                0405          do num_file=1,NFILESPROFMAX
                0406           fid=fiddata(num_file,bi,bj)
                0407
                0408           if ( (ProfNo(num_file,bi,bj).GT.0).AND.
                0409      &         (profilesDoNcOutput) ) then
                0410 c need to close the file so that the data is not lost when run finishes
                0411            err = NF_CLOSE(fidforward(num_file,bi,bj))
                0412 c then re-open it to compute cost function
                0413            iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
                0414            jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
                0415            IL  = ILNBLNK( profilesfiles(num_file) )
0e6a4460e5 Ivan*0416            write(profilesfile,'(a)')
6b2230d510 Ou W*0417      &     profilesfiles(num_file)(1:IL)
                0418            IL  = ILNBLNK( profilesfile )
                0419            JL  = ILNBLNK( profilesDir )
0e6a4460e5 Ivan*0420            write(fnameequinc,'(3a,2(i3.3,a))')
6b2230d510 Ou W*0421      &     profilesDir(1:JL),profilesfile(1:IL),'.',iG,'.',jG,'.equi.nc'
                0422 c
                0423            err = NF_OPEN(fnameequinc,NF_NOWRITE,
                0424      &     fidforward(num_file,bi,bj))
                0425           endif
                0426
                0427           do prof_num=1,NOBSGLOB
                0428            if (prof_num.LE.ProfNo(num_file,bi,bj)) then
6a770e0a24 Patr*0429
0e0f68501f Gael*0430 c would be needed to call profiles_interp to e.g. get time averages
6b2230d510 Ou W*0431 c           do k=1,NUM_INTERP_POINTS
0e0f68501f Gael*0432 c           prof_i1D(k)= prof_interp_i(num_file,prof_num,k,bi,bj)
                0433 c           prof_j1D(k)= prof_interp_j(num_file,prof_num,k,bi,bj)
                0434 c           prof_w1D(k)= prof_interp_weights(num_file,prof_num,k,bi,bj)
                0435 c          enddo
ba63501b4c Gael*0436
ea4d09597a Gael*0437            do num_var=1,NVARMAX
6a770e0a24 Patr*0438
5058c82dba Patr*0439             do K=1,NLEVELMAX
                0440              prof_traj1D(k)=0.
                0441              prof_traj1D_mean(k)=0.
                0442              prof_mask1D_cur(k,bi,bj)=0.
                0443              prof_data1D(k)=0.
                0444              prof_weights1D(k)=0.
6b2230d510 Ou W*0445 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0446              prof_data1D_mean(k)=0.
                0447              prof_weights1D_mean(k)=0.
                0448 #endif
5058c82dba Patr*0449             enddo
6a770e0a24 Patr*0450
5058c82dba Patr*0451             if (vec_quantities(num_file,num_var,bi,bj).EQV..TRUE.) then
6ebd5a1932 Patr*0452
6b2230d510 Ou W*0453 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0454              itmp = prof_ind_avgbin(num_file,prof_num,bi,bj)
                0455              if(itmp.GE. 0) then
                0456               do K=1,ProfDepthNo(num_file,bi,bj)
                0457
                0458                ktmp = prof_lev_comb(k,num_file,bi,bj)
                0459                prof_traj1D_mean(k) =
                0460      &           prof_traj1D_all_mean(itmp,ktmp,num_var)
                0461                prof_data1D_mean(k) =
                0462      &           prof_data1D_all_mean(itmp,ktmp,num_var)
                0463                prof_weights1D_mean(k) =
                0464      &           prof_weights1D_all_mean(itmp,ktmp,num_var)
                0465               enddo
                0466              endif !if(itmp.GE. 0. _d 0) then
                0467 C end of #ifndef ALLOW_PROFILES_SAMPLESPLIT_COST
                0468 #endif
                0469
5058c82dba Patr*0470              call active_read_profile(num_file,
                0471      &           ProfDepthNo(num_file,bi,bj),prof_traj1D,num_var,
65754df434 Mart*0472      &           prof_num,.false.,optimcycle,bi,bj,myThid,
5058c82dba Patr*0473      &           profiles_dummy(num_file,num_var,bi,bj))
                0474
                0475              call profiles_readvector(num_file,num_var,
                0476      &           prof_ind_glob(num_file,prof_num,bi,bj),
                0477      &           ProfDepthNo(num_file,bi,bj),prof_data1D,bi,bj,myThid)
                0478
                0479              call profiles_readvector(num_file,-num_var,
                0480      &           prof_ind_glob(num_file,prof_num,bi,bj),
                0481      &           ProfDepthNo(num_file,bi,bj),
                0482      &           prof_weights1D,bi,bj,myThid)
                0483
                0484              do K=1,ProfDepthNo(num_file,bi,bj)
6b2230d510 Ou W*0485                if (prof_weights1D(K).GT.0.
                0486 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0487      &             .AND. prof_data1D_mean(K).NE. 0. _d 0
                0488 #endif
                0489      &            ) then
4237aa7e01 Gael*0490                  objf_profiles(num_file,num_var,bi,bj)=
98bf704dd5 Jean*0491      &             objf_profiles(num_file,num_var,bi,bj)
5058c82dba Patr*0492      &             +prof_weights1D(K)*prof_mask1D_cur(K,bi,bj)
6b2230d510 Ou W*0493      &             *(prof_traj1D(K)-prof_data1D(K)-prof_traj1D_mean(K)
                0494 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0495      &               + prof_data1D_mean(K)
                0496 #endif
2df87b7069 Ou W*0497      &              )
6b2230d510 Ou W*0498      &             *(prof_traj1D(K)-prof_data1D(K)-prof_traj1D_mean(K)
                0499 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0500      &               + prof_data1D_mean(K)
                0501 #endif
2df87b7069 Ou W*0502      &              )
5058c82dba Patr*0503                  num_profiles(num_file,num_var,bi,bj)=
                0504      &               num_profiles(num_file,num_var,bi,bj)
                0505      &               +prof_mask1D_cur(K,bi,bj)
                0506                endif
                0507              enddo
                0508             endif
                0509
98bf704dd5 Jean*0510            enddo !do num_var...
5058c82dba Patr*0511           endif !if (prof_num.LE.ProfNo(num_file,bi,bj)) then
                0512          enddo !do prof_num=..
                0513
275a56dc21 Gael*0514 #ifdef ALLOW_DEBUG
                0515       IF ( debugLevel .GE. debLevD ) THEN
5058c82dba Patr*0516          if (ProfNo(num_file,bi,bj).GT.0) then
ea4d09597a Gael*0517           do num_var=1,NVARMAX
ba63501b4c Gael*0518            write(msgbuf,'(a,4I9)') 'bi,bj,prof_num,num_var ',bi,bj,
5058c82dba Patr*0519      &      ProfNo(num_file,bi,bj),num_var
ba63501b4c Gael*0520            call print_message(
65754df434 Mart*0521      &      msgbuf, standardMessageUnit, SQUEEZE_RIGHT , myThid)
0e6a4460e5 Ivan*0522            write(msgbuf,'(a,2D22.15)')
cf16ba6028 Gael*0523      &      prof_names(num_file,num_var),
5058c82dba Patr*0524      &      objf_profiles(num_file,num_var,bi,bj),
                0525      &      num_profiles(num_file,num_var,bi,bj)
                0526           enddo !do num_var...
                0527          endif
275a56dc21 Gael*0528       ENDIF
                0529 #endif /* ALLOW_DEBUG */
5058c82dba Patr*0530         enddo !do num_file=1,NFILESPROFMAX
                0531
6b2230d510 Ou W*0532 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0533       do num_var=1,NVARMAX
                0534        DO iavgbin = 1, NAVGBINMAX
                0535           do K=1,NLEVELCOMB
                0536            prof_traj1D_mean(1) =
                0537      &      prof_traj1D_all_mean(iavgbin,k,num_var)
                0538            prof_data1D_mean(1) =
                0539      &      prof_data1D_all_mean(iavgbin,k,num_var)
                0540            prof_weights1D_mean(1) =
                0541      &      prof_weights1D_all_mean(iavgbin,k,num_var)
                0542
                0543            if (prof_weights1D_mean(1).GT.0.
                0544      &         .AND. prof_data1D_mean(1).NE. 0. _d 0
                0545      &         .AND. prof_traj1D_mean(1).NE. 0. _d 0
                0546 C    &         .AND. myProcId .eq. 0
                0547      &         .AND. avgbinglbsum(iavgbin).GT.0
                0548      &        ) then
                0549              if(avgbinglbsum(iavgbin).EQ.1) then
                0550               objf_profiles_mean(num_var,bi,bj)=
                0551      &          objf_profiles_mean(num_var,bi,bj)
                0552      &          +prof_weights1D_mean(1)
                0553      &          *(prof_traj1D_mean(1)
                0554      &            - prof_data1D_mean(1)
                0555      &           )
                0556      &          *(prof_traj1D_mean(1)
                0557      &            - prof_data1D_mean(1)
                0558      &           )
                0559               num_profiles_mean(num_var,bi,bj)=
                0560      &            num_profiles_mean(num_var,bi,bj)
                0561      &            +1. _d 0
                0562              else
                0563               objf_profiles_mean(num_var,bi,bj)=
                0564      &          objf_profiles_mean(num_var,bi,bj)
                0565      &          +prof_weights1D_mean(1)
                0566      &          *(prof_traj1D_mean(1)
                0567      &            - prof_data1D_mean(1)
                0568      &           )
                0569      &          *(prof_traj1D_mean(1)
                0570      &            - prof_data1D_mean(1)
                0571      &           )/numberOfProcs
                0572               num_profiles_mean(num_var,bi,bj)=
                0573      &            num_profiles_mean(num_var,bi,bj)
                0574      &            +1. _d 0/numberOfProcs
                0575              endif ! if(avgbinglbsum(iavgbin).EQ.1) then
                0576
                0577             endif ! if (prof_weights1D_mean(1).GT.0.
                0578           enddo !do K=1,NLEVELCOMB
                0579        enddo !DO iavgbin = 1
                0580       enddo !do num_var
                0581
                0582 #ifdef ALLOW_DEBUG
                0583       IF ( debugLevel .GE. debLevD ) THEN
                0584 C        if (ProfNo(num_file,bi,bj).GT.0) then
                0585           do num_var=1,NVARMAX
                0586
                0587            write(msgbuf,'(a,4I9)') 'bi,bj,num_var ',bi,bj,
                0588      &      num_var
                0589            call print_message(
65754df434 Mart*0590      &      msgbuf, standardMessageUnit, SQUEEZE_RIGHT , myThid)
6b2230d510 Ou W*0591
0e6a4460e5 Ivan*0592            write(msgbuf,'(a,a5,2D22.15)') prof_names(1,num_var),
6b2230d510 Ou W*0593      &      '_mean',
                0594      &      objf_profiles_mean(num_var,bi,bj),
                0595      &      num_profiles_mean(num_var,bi,bj)
                0596            call print_message(
65754df434 Mart*0597      &      msgbuf, standardMessageUnit, SQUEEZE_RIGHT , myThid)
6b2230d510 Ou W*0598
                0599           enddo !do num_var...
                0600 C        endif
                0601       ENDIF
                0602 #endif /* ALLOW_DEBUG */
                0603
                0604 C        enddo !do num_file
                0605
                0606 #endif
                0607
5058c82dba Patr*0608        ENDDO
6a770e0a24 Patr*0609       ENDDO
6cc0cb518f Gael*0610
65754df434 Mart*0611       _END_MASTER( myThid )
6cc0cb518f Gael*0612
3c8dcfdea9 Gael*0613 c print cost function values
                0614       do num_file=1,NFILESPROFMAX
                0615       do num_var=1,NVARMAX
                0616 c
65754df434 Mart*0617       do bj = myByLo(myThid),myByHi(myThid)
                0618         do bi = myBxLo(myThid),myBxHi(myThid)
3c8dcfdea9 Gael*0619           objf_prof_tile(bi,bj) =
                0620      &             objf_profiles(num_file,num_var,bi,bj)
                0621           num_prof_tile(bi,bj) =
                0622      &             num_profiles(num_file,num_var,bi,bj)
                0623        enddo
                0624       enddo
                0625 c
                0626       CALL GLOBAL_SUM_TILE_RL( objf_prof_tile, objf_prof_glo, myThid )
                0627       CALL GLOBAL_SUM_TILE_RL( num_prof_tile, num_prof_glo, myThid )
                0628 c
0e6a4460e5 Ivan*0629       write(msgbuf,'(2(a,I2),a,2D12.5)')
3c8dcfdea9 Gael*0630      &  ' cost_profiles(',num_file,',',num_var,')= ',
                0631      &  objf_prof_glo,num_prof_glo
                0632
                0633       IF ( num_prof_glo .GT. 0. ) call print_message( msgbuf,
65754df434 Mart*0634      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0635 c
                0636       enddo
                0637       enddo
                0638
6b2230d510 Ou W*0639 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0640       do num_var=1,NVARMAX
                0641 c
65754df434 Mart*0642       do bj = myByLo(myThid),myByHi(myThid)
                0643         do bi = myBxLo(myThid),myBxHi(myThid)
6b2230d510 Ou W*0644           objf_prof_mean_tile(bi,bj) =
                0645      &             objf_profiles_mean(num_var,bi,bj)
                0646           num_prof_mean_tile(bi,bj) =
                0647      &             num_profiles_mean(num_var,bi,bj)
                0648        enddo
                0649       enddo
                0650 c
                0651       CALL GLOBAL_SUM_TILE_RL( objf_prof_mean_tile,
                0652      &     objf_prof_mean_glo, myThid )
                0653       CALL GLOBAL_SUM_TILE_RL( num_prof_mean_tile,
                0654      &     num_prof_mean_glo, myThid )
                0655 c
                0656       write(msgbuf,'(a,I2,a,2D12.5)')
                0657      &  ' cost_profiles_mean(',num_var,')= ',
                0658      &  objf_prof_mean_glo,num_prof_mean_glo
                0659
                0660       IF ( num_prof_mean_glo .GT. 0. ) call print_message( msgbuf,
65754df434 Mart*0661      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
6b2230d510 Ou W*0662 c
                0663       enddo
                0664
                0665 #endif
                0666 C! ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0667
3c8dcfdea9 Gael*0668       write(msgbuf,'(a)') '== cost_profiles: end   =='
                0669       call print_message( msgbuf,
65754df434 Mart*0670      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0671       write(msgbuf,'(a)') ' '
                0672       call print_message( msgbuf,
65754df434 Mart*0673      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0674
d3172737dc aver*0675       IF ( prof_make_nc ) THEN
                0676        call profiles_make_ncfile(myThid)
                0677       ENDIF
6328b73337 Gael*0678
6a770e0a24 Patr*0679 C===========================================================
6e4c90fea3 Patr*0680
                0681 #endif
6a770e0a24 Patr*0682
99621390b7 Jean*0683       RETURN
6a770e0a24 Patr*0684       END