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) )
de57a2ec4b Mart*0140            write(profilesfile,'(1a)')
b2a948f981 Gael*0141      &     profilesfiles(num_file)(1:IL)
                0142            IL  = ILNBLNK( profilesfile )
                0143            JL  = ILNBLNK( profilesDir )
de57a2ec4b Mart*0144            write(fnameequinc,'(3a,i3.3,a,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
65754df434 Mart*0333 C     _BEGIN_MASTER( myThid )
6b2230d510 Ou W*0334        IF ( myProcId .eq. 0 ) THEN
                0335
                0336         DO num_var=1,NVARMAX
                0337          iL = ILNBLNK( prof_names(1,num_var) )
                0338          write(cfile,'(2a)') prof_names(1,num_var)(1:iL),
                0339      &   '_data_mean.data'
65754df434 Mart*0340          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0341          open( cunit, file   = cfile,
                0342      &        status = 'unknown',
                0343      &        access  = 'direct',
                0344      &        recl = NAVGBINMAX*4)
65754df434 Mart*0345
6b2230d510 Ou W*0346          DO ikzz = 1, NLEVELCOMB
                0347           tmp99(1:NAVGBINMAX)=
                0348      &      prof_data1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0349           write(cunit,rec=ikzz) tmp99
                0350          ENDDO
                0351          close ( cunit )
                0352
                0353          write(cfile,'(2a)')prof_names(1,num_var)(1:iL),
                0354      &    '_model_mean.data'
65754df434 Mart*0355          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0356          open( cunit, file   = cfile,
                0357      &         status = 'unknown',
                0358 C    &         form   = 'unformatted',
                0359      &         access  = 'direct',
                0360      &         recl = NAVGBINMAX*4)
                0361 C    &         access  = 'sequential'   )
                0362
                0363          DO ikzz = 1, NLEVELCOMB
                0364           tmp99(1:NAVGBINMAX)=
                0365      &      prof_traj1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0366           write(cunit,rec=ikzz) tmp99
                0367          ENDDO
                0368          close ( cunit )
                0369
                0370          write(cfile,'(2a)')
                0371      &     prof_names(1,num_var)(1:iL),'_weight_mean.data'
65754df434 Mart*0372          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0373          open( cunit, file   = cfile,
                0374      &         status = 'unknown',
                0375      &         access  = 'direct',
                0376      &         recl = NAVGBINMAX*4)
                0377
                0378          DO ikzz = 1, NLEVELCOMB
                0379           tmp99(1:NAVGBINMAX)=
                0380      &      prof_weights1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0381           write(cunit,rec=ikzz) tmp99
                0382          ENDDO
                0383          close ( cunit )
                0384
                0385          write(cfile,'(2a)')prof_names(1,num_var)(1:iL),
                0386      &    '_count_mean.data'
65754df434 Mart*0387          call mdsfindunit( cunit, myThid )
6b2230d510 Ou W*0388          open( cunit, file   = cfile,
                0389      &         status = 'unknown',
                0390      &         access  = 'direct',
                0391      &         recl = NAVGBINMAX*4)
65754df434 Mart*0392
6b2230d510 Ou W*0393          DO ikzz = 1, NLEVELCOMB
                0394           tmp99(1:NAVGBINMAX)=
                0395      &     prof_count1D_all_mean(1:NAVGBINMAX,ikzz,num_var)
                0396           write(cunit,rec=ikzz) tmp99
                0397          ENDDO
                0398          close ( cunit )
                0399         ENDDO ! DO num_var=1,NVARMAX
                0400        ENDIF ! IF ( myProcId .eq. 0 ) THEN
65754df434 Mart*0401 C      _END_MASTER( myThid )
6b2230d510 Ou W*0402        _BARRIER
                0403 #endif
                0404
                0405        DO bj=1,nSy
                0406         DO bi=1,nSx
                0407
                0408          do num_file=1,NFILESPROFMAX
                0409           fid=fiddata(num_file,bi,bj)
                0410
                0411           if ( (ProfNo(num_file,bi,bj).GT.0).AND.
                0412      &         (profilesDoNcOutput) ) then
                0413 c need to close the file so that the data is not lost when run finishes
                0414            err = NF_CLOSE(fidforward(num_file,bi,bj))
                0415 c then re-open it to compute cost function
                0416            iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles
                0417            jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles
                0418            IL  = ILNBLNK( profilesfiles(num_file) )
de57a2ec4b Mart*0419            write(profilesfile,'(1a)')
6b2230d510 Ou W*0420      &     profilesfiles(num_file)(1:IL)
                0421            IL  = ILNBLNK( profilesfile )
                0422            JL  = ILNBLNK( profilesDir )
de57a2ec4b Mart*0423            write(fnameequinc,'(3a,i3.3,a,i3.3,a)')
6b2230d510 Ou W*0424      &     profilesDir(1:JL),profilesfile(1:IL),'.',iG,'.',jG,'.equi.nc'
                0425 c
                0426            err = NF_OPEN(fnameequinc,NF_NOWRITE,
                0427      &     fidforward(num_file,bi,bj))
                0428           endif
                0429
                0430           do prof_num=1,NOBSGLOB
                0431            if (prof_num.LE.ProfNo(num_file,bi,bj)) then
6a770e0a24 Patr*0432
0e0f68501f Gael*0433 c would be needed to call profiles_interp to e.g. get time averages
6b2230d510 Ou W*0434 c           do k=1,NUM_INTERP_POINTS
0e0f68501f Gael*0435 c           prof_i1D(k)= prof_interp_i(num_file,prof_num,k,bi,bj)
                0436 c           prof_j1D(k)= prof_interp_j(num_file,prof_num,k,bi,bj)
                0437 c           prof_w1D(k)= prof_interp_weights(num_file,prof_num,k,bi,bj)
                0438 c          enddo
ba63501b4c Gael*0439
ea4d09597a Gael*0440            do num_var=1,NVARMAX
6a770e0a24 Patr*0441
5058c82dba Patr*0442             do K=1,NLEVELMAX
                0443              prof_traj1D(k)=0.
                0444              prof_traj1D_mean(k)=0.
                0445              prof_mask1D_cur(k,bi,bj)=0.
                0446              prof_data1D(k)=0.
                0447              prof_weights1D(k)=0.
6b2230d510 Ou W*0448 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0449              prof_data1D_mean(k)=0.
                0450              prof_weights1D_mean(k)=0.
                0451 #endif
5058c82dba Patr*0452             enddo
6a770e0a24 Patr*0453
5058c82dba Patr*0454             if (vec_quantities(num_file,num_var,bi,bj).EQV..TRUE.) then
6ebd5a1932 Patr*0455
6b2230d510 Ou W*0456 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0457              itmp = prof_ind_avgbin(num_file,prof_num,bi,bj)
                0458              if(itmp.GE. 0) then
                0459               do K=1,ProfDepthNo(num_file,bi,bj)
                0460
                0461                ktmp = prof_lev_comb(k,num_file,bi,bj)
                0462                prof_traj1D_mean(k) =
                0463      &           prof_traj1D_all_mean(itmp,ktmp,num_var)
                0464                prof_data1D_mean(k) =
                0465      &           prof_data1D_all_mean(itmp,ktmp,num_var)
                0466                prof_weights1D_mean(k) =
                0467      &           prof_weights1D_all_mean(itmp,ktmp,num_var)
                0468               enddo
                0469              endif !if(itmp.GE. 0. _d 0) then
                0470 C end of #ifndef ALLOW_PROFILES_SAMPLESPLIT_COST
                0471 #endif
                0472
5058c82dba Patr*0473              call active_read_profile(num_file,
                0474      &           ProfDepthNo(num_file,bi,bj),prof_traj1D,num_var,
65754df434 Mart*0475      &           prof_num,.false.,optimcycle,bi,bj,myThid,
5058c82dba Patr*0476      &           profiles_dummy(num_file,num_var,bi,bj))
                0477
                0478              call profiles_readvector(num_file,num_var,
                0479      &           prof_ind_glob(num_file,prof_num,bi,bj),
                0480      &           ProfDepthNo(num_file,bi,bj),prof_data1D,bi,bj,myThid)
                0481
                0482              call profiles_readvector(num_file,-num_var,
                0483      &           prof_ind_glob(num_file,prof_num,bi,bj),
                0484      &           ProfDepthNo(num_file,bi,bj),
                0485      &           prof_weights1D,bi,bj,myThid)
                0486
                0487              do K=1,ProfDepthNo(num_file,bi,bj)
6b2230d510 Ou W*0488                if (prof_weights1D(K).GT.0.
                0489 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0490      &             .AND. prof_data1D_mean(K).NE. 0. _d 0
                0491 #endif
                0492      &            ) then
4237aa7e01 Gael*0493                  objf_profiles(num_file,num_var,bi,bj)=
98bf704dd5 Jean*0494      &             objf_profiles(num_file,num_var,bi,bj)
5058c82dba Patr*0495      &             +prof_weights1D(K)*prof_mask1D_cur(K,bi,bj)
6b2230d510 Ou W*0496      &             *(prof_traj1D(K)-prof_data1D(K)-prof_traj1D_mean(K)
                0497 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0498      &               + prof_data1D_mean(K)
                0499 #endif
2df87b7069 Ou W*0500      &              )
6b2230d510 Ou W*0501      &             *(prof_traj1D(K)-prof_data1D(K)-prof_traj1D_mean(K)
                0502 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0503      &               + prof_data1D_mean(K)
                0504 #endif
2df87b7069 Ou W*0505      &              )
5058c82dba Patr*0506                  num_profiles(num_file,num_var,bi,bj)=
                0507      &               num_profiles(num_file,num_var,bi,bj)
                0508      &               +prof_mask1D_cur(K,bi,bj)
                0509                endif
                0510              enddo
                0511             endif
                0512
98bf704dd5 Jean*0513            enddo !do num_var...
5058c82dba Patr*0514           endif !if (prof_num.LE.ProfNo(num_file,bi,bj)) then
                0515          enddo !do prof_num=..
                0516
275a56dc21 Gael*0517 #ifdef ALLOW_DEBUG
                0518       IF ( debugLevel .GE. debLevD ) THEN
5058c82dba Patr*0519          if (ProfNo(num_file,bi,bj).GT.0) then
ea4d09597a Gael*0520           do num_var=1,NVARMAX
ba63501b4c Gael*0521            write(msgbuf,'(a,4I9)') 'bi,bj,prof_num,num_var ',bi,bj,
5058c82dba Patr*0522      &      ProfNo(num_file,bi,bj),num_var
ba63501b4c Gael*0523            call print_message(
65754df434 Mart*0524      &      msgbuf, standardMessageUnit, SQUEEZE_RIGHT , myThid)
2767dff983 Jean*0525            write(msgbuf,'(a,D22.15,D22.15)')
cf16ba6028 Gael*0526      &      prof_names(num_file,num_var),
5058c82dba Patr*0527      &      objf_profiles(num_file,num_var,bi,bj),
                0528      &      num_profiles(num_file,num_var,bi,bj)
                0529           enddo !do num_var...
                0530          endif
275a56dc21 Gael*0531       ENDIF
                0532 #endif /* ALLOW_DEBUG */
5058c82dba Patr*0533         enddo !do num_file=1,NFILESPROFMAX
                0534
6b2230d510 Ou W*0535 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0536       do num_var=1,NVARMAX
                0537        DO iavgbin = 1, NAVGBINMAX
                0538           do K=1,NLEVELCOMB
                0539            prof_traj1D_mean(1) =
                0540      &      prof_traj1D_all_mean(iavgbin,k,num_var)
                0541            prof_data1D_mean(1) =
                0542      &      prof_data1D_all_mean(iavgbin,k,num_var)
                0543            prof_weights1D_mean(1) =
                0544      &      prof_weights1D_all_mean(iavgbin,k,num_var)
                0545
                0546            if (prof_weights1D_mean(1).GT.0.
                0547      &         .AND. prof_data1D_mean(1).NE. 0. _d 0
                0548      &         .AND. prof_traj1D_mean(1).NE. 0. _d 0
                0549 C    &         .AND. myProcId .eq. 0
                0550      &         .AND. avgbinglbsum(iavgbin).GT.0
                0551      &        ) then
                0552              if(avgbinglbsum(iavgbin).EQ.1) then
                0553               objf_profiles_mean(num_var,bi,bj)=
                0554      &          objf_profiles_mean(num_var,bi,bj)
                0555      &          +prof_weights1D_mean(1)
                0556      &          *(prof_traj1D_mean(1)
                0557      &            - prof_data1D_mean(1)
                0558      &           )
                0559      &          *(prof_traj1D_mean(1)
                0560      &            - prof_data1D_mean(1)
                0561      &           )
                0562               num_profiles_mean(num_var,bi,bj)=
                0563      &            num_profiles_mean(num_var,bi,bj)
                0564      &            +1. _d 0
                0565              else
                0566               objf_profiles_mean(num_var,bi,bj)=
                0567      &          objf_profiles_mean(num_var,bi,bj)
                0568      &          +prof_weights1D_mean(1)
                0569      &          *(prof_traj1D_mean(1)
                0570      &            - prof_data1D_mean(1)
                0571      &           )
                0572      &          *(prof_traj1D_mean(1)
                0573      &            - prof_data1D_mean(1)
                0574      &           )/numberOfProcs
                0575               num_profiles_mean(num_var,bi,bj)=
                0576      &            num_profiles_mean(num_var,bi,bj)
                0577      &            +1. _d 0/numberOfProcs
                0578              endif ! if(avgbinglbsum(iavgbin).EQ.1) then
                0579
                0580             endif ! if (prof_weights1D_mean(1).GT.0.
                0581           enddo !do K=1,NLEVELCOMB
                0582        enddo !DO iavgbin = 1
                0583       enddo !do num_var
                0584
                0585 #ifdef ALLOW_DEBUG
                0586       IF ( debugLevel .GE. debLevD ) THEN
                0587 C        if (ProfNo(num_file,bi,bj).GT.0) then
                0588           do num_var=1,NVARMAX
                0589
                0590            write(msgbuf,'(a,4I9)') 'bi,bj,num_var ',bi,bj,
                0591      &      num_var
                0592            call print_message(
65754df434 Mart*0593      &      msgbuf, standardMessageUnit, SQUEEZE_RIGHT , myThid)
6b2230d510 Ou W*0594
                0595            write(msgbuf,'(a,a5,D22.15,D22.15)') prof_names(1,num_var),
                0596      &      '_mean',
                0597      &      objf_profiles_mean(num_var,bi,bj),
                0598      &      num_profiles_mean(num_var,bi,bj)
                0599            call print_message(
65754df434 Mart*0600      &      msgbuf, standardMessageUnit, SQUEEZE_RIGHT , myThid)
6b2230d510 Ou W*0601
                0602           enddo !do num_var...
                0603 C        endif
                0604       ENDIF
                0605 #endif /* ALLOW_DEBUG */
                0606
                0607 C        enddo !do num_file
                0608
                0609 #endif
                0610
5058c82dba Patr*0611        ENDDO
6a770e0a24 Patr*0612       ENDDO
6cc0cb518f Gael*0613
65754df434 Mart*0614       _END_MASTER( myThid )
6cc0cb518f Gael*0615
3c8dcfdea9 Gael*0616 c print cost function values
                0617       do num_file=1,NFILESPROFMAX
                0618       do num_var=1,NVARMAX
                0619 c
65754df434 Mart*0620       do bj = myByLo(myThid),myByHi(myThid)
                0621         do bi = myBxLo(myThid),myBxHi(myThid)
3c8dcfdea9 Gael*0622           objf_prof_tile(bi,bj) =
                0623      &             objf_profiles(num_file,num_var,bi,bj)
                0624           num_prof_tile(bi,bj) =
                0625      &             num_profiles(num_file,num_var,bi,bj)
                0626        enddo
                0627       enddo
                0628 c
                0629       CALL GLOBAL_SUM_TILE_RL( objf_prof_tile, objf_prof_glo, myThid )
                0630       CALL GLOBAL_SUM_TILE_RL( num_prof_tile, num_prof_glo, myThid )
                0631 c
                0632       write(msgbuf,'(a,I2,a,I2,a,2D12.5)')
                0633      &  ' cost_profiles(',num_file,',',num_var,')= ',
                0634      &  objf_prof_glo,num_prof_glo
                0635
                0636       IF ( num_prof_glo .GT. 0. ) call print_message( msgbuf,
65754df434 Mart*0637      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0638 c
                0639       enddo
                0640       enddo
                0641
6b2230d510 Ou W*0642 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0643       do num_var=1,NVARMAX
                0644 c
65754df434 Mart*0645       do bj = myByLo(myThid),myByHi(myThid)
                0646         do bi = myBxLo(myThid),myBxHi(myThid)
6b2230d510 Ou W*0647           objf_prof_mean_tile(bi,bj) =
                0648      &             objf_profiles_mean(num_var,bi,bj)
                0649           num_prof_mean_tile(bi,bj) =
                0650      &             num_profiles_mean(num_var,bi,bj)
                0651        enddo
                0652       enddo
                0653 c
                0654       CALL GLOBAL_SUM_TILE_RL( objf_prof_mean_tile,
                0655      &     objf_prof_mean_glo, myThid )
                0656       CALL GLOBAL_SUM_TILE_RL( num_prof_mean_tile,
                0657      &     num_prof_mean_glo, myThid )
                0658 c
                0659       write(msgbuf,'(a,I2,a,2D12.5)')
                0660      &  ' cost_profiles_mean(',num_var,')= ',
                0661      &  objf_prof_mean_glo,num_prof_mean_glo
                0662
                0663       IF ( num_prof_mean_glo .GT. 0. ) call print_message( msgbuf,
65754df434 Mart*0664      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
6b2230d510 Ou W*0665 c
                0666       enddo
                0667
                0668 #endif
                0669 C! ifdef ALLOW_PROFILES_SAMPLESPLIT_COST
                0670
3c8dcfdea9 Gael*0671       write(msgbuf,'(a)') '== cost_profiles: end   =='
                0672       call print_message( msgbuf,
65754df434 Mart*0673      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0674       write(msgbuf,'(a)') ' '
                0675       call print_message( msgbuf,
65754df434 Mart*0676      &  standardMessageUnit,SQUEEZE_RIGHT , myThid)
3c8dcfdea9 Gael*0677
65754df434 Mart*0678 c     call profiles_make_ncfile(myThid)
6328b73337 Gael*0679
6a770e0a24 Patr*0680 C===========================================================
6e4c90fea3 Patr*0681
                0682 #endif
6a770e0a24 Patr*0683
99621390b7 Jean*0684       RETURN
6a770e0a24 Patr*0685       END