Back to home page

MITgcm
 
 

    

File indexing completed on 2025-12-15 06:14:27 UTC

view on githubraw file Latest commit ad59256d on 2025-12-15 00:05:36 UTC
ad59256d7d aver*0001 #include "OBSFIT_OPTIONS.h"
                0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
                0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
                0008 
                0009 CBOP
                0010 C     !ROUTINE: OBSFIT_COST
                0011 
                0012 C     !INTERFACE:
                0013       SUBROUTINE OBSFIT_COST( myTime, myIter, myThid )
                0014 
                0015 C     !DESCRIPTION:
                0016 C     ==================================================================
                0017 C     | Computes the cost function for ObsFit data
                0018 C     ==================================================================
                0019 
                0020 C     !USES:
                0021       IMPLICIT NONE
                0022 C     == Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "DYNVARS.h"
                0028 #ifdef ALLOW_CAL
                0029 # include "cal.h"
                0030 #endif
                0031 #ifdef ALLOW_OBSFIT
                0032 # include "OBSFIT_SIZE.h"
                0033 # include "OBSFIT.h"
                0034 # include "netcdf.inc"
                0035 #endif
                0036 #ifdef ALLOW_CTRL
                0037 # include "OPTIMCYCLE.h"
                0038 #endif
                0039 #ifdef ALLOW_AUTODIFF
                0040 # include "tamc.h"
                0041 #endif
                0042 
                0043 C     !INPUT PARAMETERS:
                0044 C     myTime  :: Current time in simulation
                0045 C     myIter  :: Current iteration number in simulation
                0046 C     myThid  :: my thread ID number
                0047       _RL     myTime
                0048       INTEGER myIter
                0049       INTEGER myThid
                0050 CEOP
                0051 
                0052 #ifdef ALLOW_OBSFIT
                0053 C     !FUNCTIONS
                0054       INTEGER  ILNBLNK
                0055       EXTERNAL ILNBLNK
                0056 
                0057 C     !LOCAL VARIABLES:
                0058       INTEGER num_file,sample_num
                0059       INTEGER bi,bj
                0060       _RL sample_modval
                0061       _RL obs_modval, obs_modmask
                0062       _RL obs_data, obs_uncert, obs_weight
                0063 #ifndef ALLOW_CTRL
                0064       INTEGER optimcycle
                0065 #endif
                0066       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0067       CHARACTER*(MAX_LEN_FNAM) obsfitfile
                0068       CHARACTER*(MAX_LEN_FNAM) fnamemisfit
                0069       INTEGER IL, JL, err
                0070       INTEGER irec, ii, varID
                0071       INTEGER obs_num, num_valid_samples
                0072       _RL sample_mask_sum
                0073       _RL objf_obsfit_glo
                0074       _RL num_obsfit_glo
                0075       _RL samples_buff(NSAMPLES_MAX_GLO)
                0076       _RL samples_mask_buff(NSAMPLES_MAX_GLO)
                0077       _RL samples_modval_glob(NSAMPLES_MAX_GLO)
                0078       _RL samples_mask_glob(NSAMPLES_MAX_GLO)
                0079       _RL tmpgs
                0080       INTEGER nobsmean
                0081       _RL offset, mod_mean, obs_mean, misfit
                0082       _RL spval
                0083       PARAMETER ( spval = -9999. _d 0 )
                0084 
                0085 #ifndef ALLOW_CTRL
                0086       optimcycle = 0
                0087 #endif
                0088 
                0089       WRITE( msgBuf,'(A)' ) ' '
                0090       CALL PRINT_MESSAGE( msgBuf,
                0091      &     standardMessageUnit,SQUEEZE_RIGHT, myThid )
                0092       WRITE( msgBuf,'(A)' ) '== obsfit_cost: begin =='
                0093       CALL PRINT_MESSAGE( msgBuf,
                0094      &     standardMessageUnit,SQUEEZE_RIGHT, myThid )
                0095 
                0096 C Initialise local storage
                0097 #ifdef ALLOW_AUTODIFF_TAMC
                0098 CADJ INIT tapelev_obsfit = COMMON, NFILESMAX_OBS
                0099 #endif
                0100 
                0101       _BEGIN_MASTER( myThid )
                0102 
                0103       DO num_file = 1, NFILESMAX_OBS
                0104 
                0105 C File maintenance
                0106         DO bj = 1, nSy
                0107           DO bi= 1, nSx
                0108 
                0109             IF ( ( sampleNo(num_file,bi,bj).GT.0 ).AND.
                0110      &           (obsfitDoNcOutput) ) THEN
                0111 C Need to sync the file so that the data is not lost when run finishes
                0112               err = NF_SYNC( fidfwd_obs(num_file,bi,bj) )
                0113               CALL OBSFIT_NF_ERROR( 'COST: NF_SYNC fidfwd_obs',
                0114      &             err,bi,bj,myThid )
                0115             ENDIF
                0116 
                0117           ENDDO
                0118         ENDDO
                0119 
                0120 C Loop over samples
                0121         DO ii = 1, NSAMP_PER_TILE_MAX
                0122           samples_buff(ii)      = zeroRL
                0123           samples_mask_buff(ii) = zeroRL
                0124         ENDDO
                0125 
                0126         DO bj = 1, nSy
                0127           DO bi = 1, nSx
                0128 
                0129 C Open tiled files and read to buffer
                0130             DO sample_num = 1, NSAMP_PER_TILE_MAX
                0131               IF ( sample_num.LE.sampleNo(num_file,bi,bj) ) THEN
                0132 
                0133                 sample_modval = zeroRL
                0134 
                0135                 CALL ACTIVE_READ_OBS_TILE(
                0136      I               num_file,
                0137      O               sample_modval,
                0138      I               sample_num,.false.,
                0139      I               optimcycle,bi,bj,myThid,
                0140      I               obsfit_dummy(num_file,bi,bj) )
                0141 
                0142 C Save model equi (of samples) and masks in buffer
                0143 C Combine all threads here
                0144                 irec = sample_ind_glob(num_file,sample_num,bi,bj)
                0145                 samples_buff(irec) = samples_buff(irec)+sample_modval
                0146                 samples_mask_buff(irec) = samples_mask_buff(irec)
                0147      &           +sample_modmask(bi,bj)
                0148 
                0149               ENDIF !IF (sample_num.LE.sampleNo(num_file,bi,bj))
                0150             ENDDO !DO sample_num
                0151 
                0152           ENDDO !DO bj
                0153         ENDDO !DO bi
                0154 
                0155 C Combine all processes
                0156         DO ii = 1, NSAMP_PER_TILE_MAX
                0157           tmpgs = samples_buff(ii)
                0158           _GLOBAL_SUM_RL(tmpgs, myThid)
                0159           samples_modval_glob(ii) = tmpgs
                0160           tmpgs = samples_mask_buff(ii)
                0161           _GLOBAL_SUM_RL(tmpgs, myThid)
                0162           samples_mask_glob(ii) = tmpgs
                0163         ENDDO
                0164 
                0165 #ifdef ALLOW_AUTODIFF_TAMC
                0166 C     Not really necessary, but cheap (for small NFILESMAX_OBS)
                0167 C     and avoids some hidden recomputations.
                0168 CADJ STORE samples_modval_glob, samples_mask_glob
                0169 CADJ &     = tapelev_obsfit, key = num_file, kind = isbyte
                0170 #endif
                0171 
                0172         IF ( myProcId .EQ. 0 ) THEN
                0173 
                0174 C Loop over obs
                0175           DO obs_num = 1, NOBSMAX_OBS
                0176             IF ( obs_num.LE.ObsNo(num_file) ) THEN
                0177 
                0178               obs_modval        = zeroRL
                0179               sample_mask_sum   = zeroRL
                0180               num_valid_samples = 0
                0181 
                0182 C Calculate model equi of each obs by averaging NP samples
                0183               DO sample_num = 1, NSAMP_PER_OBS_MAX
                0184                 IF ( sample_num.LE.obs_np(num_file,obs_num) ) THEN
                0185                   irec = obs_sample1_ind(num_file,obs_num)+sample_num-1
                0186                   obs_modval = obs_modval+samples_modval_glob(irec)
                0187      &             *samples_mask_glob(irec)
                0188                   sample_mask_sum = sample_mask_sum
                0189      &             +samples_mask_glob(irec)
                0190                   IF ( samples_mask_glob(irec).GT.zeroRL ) THEN
                0191                     num_valid_samples = num_valid_samples+1
                0192                   ENDIF
                0193 
                0194                 ENDIF
                0195               ENDDO
                0196 
                0197 C Time averaging
                0198               IF ( obsfitOperation(num_file).EQ.1 ) THEN
                0199                 obs_modval = obs_modval/obs_delT(num_file, obs_num)
                0200               ENDIF
                0201 
                0202 C Spatial averaging
                0203               IF (sample_mask_sum.GT.zeroRL) THEN
                0204                 obs_modval = obs_modval/sample_mask_sum
                0205                 obs_modmask = oneRL
                0206               ELSE
                0207                 obs_modval = spval
                0208                 obs_modmask = zeroRL
                0209               ENDIF
                0210 
                0211 C Write to global netcdf file
                0212               CALL ACTIVE_WRITE_OBS_GLOB(
                0213      I             num_file, obs_modval, obs_modmask, obs_num,
                0214      I             optimcycle, myThid, obsfit_globaldummy(num_file) )
                0215 
                0216             ENDIF !IF (obs_num.LE.ObsNo(num_file))
                0217           ENDDO !DO obs_num
                0218 #ifdef ALLOW_AUTODIFF_TAMC
                0219 CADJ STORE obsfit_globaldummy(num_file)
                0220 CADJ &     = tapelev_obsfit, key = num_file, kind = isbyte
                0221 #endif
                0222 
                0223 C Need to sync the file so that the data is not lost when run finishes
                0224           IF ( ObsNo(num_file).GT.0 ) THEN
                0225             err = NF_SYNC( fidglobal(num_file) )
                0226             CALL OBSFIT_NF_ERROR(
                0227      &           'COST: NF_SYNC fidglobal',err,0,0,myThid )
                0228           ENDIF
                0229           IF ( obs_is_ssh(num_file).GT.0 ) THEN
                0230 C Read data to calculate the mean offset between model and obs
                0231             offset   = zeroRL
                0232             mod_mean = zeroRL
                0233             obs_mean = zeroRL
                0234             nobsmean = 0
                0235 
                0236 C Loop over obs
                0237             DO obs_num = 1, NOBSMAX_OBS
                0238               IF ( obs_num.LE.ObsNo(num_file) ) THEN
                0239 
                0240                 obs_modval = zeroRL
                0241                 obs_data   = zeroRL
                0242                 obs_uncert = zeroRL
                0243 
                0244 C Read model equivalent from global file
                0245                 CALL ACTIVE_READ_OBS_GLOB(
                0246      I               num_file,
                0247      O               obs_modval, obs_modmask,
                0248      I               obs_num, .FALSE., optimcycle, myThid,
                0249      I               obsfit_globaldummy(num_file) )
                0250 
                0251 C Read observation and uncertainty
                0252                 CALL OBSFIT_READ_OBS(
                0253      I               num_file,1,obs_ind_glob(num_file,obs_num),
                0254      O               obs_data,
                0255      I               myThid )
                0256 
                0257                 IF ( obs_data.GT.spval .AND. obs_modval.GT.spval ) THEN
                0258                   obs_mean = obs_mean + obs_data
                0259                   mod_mean = mod_mean + obs_modval
                0260                   nobsmean = nobsmean + 1
                0261                 ENDIF
                0262 
                0263               ENDIF !IF (obs_num.LE.ObsNo(num_file))
                0264             ENDDO !DO obs_num
                0265 
                0266             obs_mean = obs_mean/nobsmean
                0267             mod_mean = mod_mean/nobsmean
                0268             offset   = mod_mean-obs_mean
                0269 
                0270           ELSE
                0271             offset = zeroRL
                0272           ENDIF
                0273 
                0274 C Read data to calculate the cost
                0275 C and write misfits to global file
                0276 
                0277 C Global file for misfits
                0278           IF (ObsNo(num_file).GT.0) THEN
                0279             IL = ILNBLNK( obsfitfiles(num_file) )
                0280             WRITE( obsfitfile,'(1A)' ) obsfitfiles(num_file)(1:IL)
                0281             IL = ILNBLNK( obsfitfile )
                0282             JL = ILNBLNK( obsfitDir )
                0283             WRITE( fnamemisfit,'(3A)' )
                0284      &       obsfitDir(1:JL),obsfitfile(1:IL),'.misfit.nc'
                0285             err = NF_OPEN( fnamemisfit,NF_WRITE,fidmisfit(num_file) )
                0286             err = NF_INQ_VARID( fidmisfit(num_file),'misfit',varID )
                0287             CALL OBSFIT_NF_ERROR(
                0288      &           'COST: NF_INQ_VARID misfit',err,0,0,myThid )
                0289           ENDIF
                0290 C Loop over obs
                0291           DO obs_num = 1, NOBSMAX_OBS
                0292             IF ( obs_num.LE.ObsNo(num_file) ) THEN
                0293 
                0294               obs_modval = zeroRL
                0295               obs_data   = zeroRL
                0296               obs_uncert = zeroRL
                0297 
                0298 C Read observation and uncertainty
                0299               CALL OBSFIT_READ_OBS(
                0300      I             num_file, 1, obs_ind_glob(num_file,obs_num),
                0301      O             obs_data,
                0302      I             myThid )
                0303 
                0304               CALL OBSFIT_READ_OBS(
                0305      I             num_file, -1, obs_ind_glob(num_file,obs_num),
                0306      O             obs_uncert,
                0307      I             myThid )
                0308 
                0309               IF ( obs_data .EQ. spval ) obs_uncert = zeroRL
                0310 
                0311 C Read model equivalent from global file
                0312               CALL ACTIVE_READ_OBS_GLOB(
                0313      I             num_file,
                0314      O             obs_modval,obs_modmask,
                0315      I             obs_num,.FALSE.,optimcycle,myThid,
                0316      I             obsfit_globaldummy(num_file) )
                0317 
                0318               IF ( obs_uncert.GT.zeroRL ) THEN
                0319                 obs_weight = 1. _d 0 / (obs_uncert*obs_uncert)
                0320                 misfit = obs_modval-obs_data-offset
                0321                 objf_obsfit(num_file) = objf_obsfit(num_file)
                0322      &           + obs_modmask* obs_weight * misfit * misfit
                0323                 num_obsfit(num_file)  = num_obsfit(num_file)
                0324      &           + obs_modmask
                0325 C Write misfit to global netcdf file
                0326                 err = NF_PUT_VARA_DOUBLE( fidmisfit(num_file),varID,
                0327      &                obs_ind_glob(num_file,obs_num),1,misfit )
                0328                 CALL OBSFIT_NF_ERROR(
                0329      &               'COST: NF_PUT_VARA_DOUBLE misfit',
                0330      &               err,0,0,myThid )
                0331               ENDIF
                0332 
                0333             ENDIF !IF (obs_num.LE.ObsNo(num_file))
                0334           ENDDO !DO obs_num
                0335 
                0336         ENDIF !IF myprocid = 0
                0337       ENDDO !DO num_file
                0338 
                0339       _END_MASTER( myThid )
                0340 
                0341 C Print cost function values
                0342       DO num_file = 1, NFILESMAX_OBS
                0343 
                0344         objf_obsfit_glo = objf_obsfit(num_file)
                0345         num_obsfit_glo = num_obsfit(num_file)
                0346 
                0347         WRITE( msgBuf,'(A,I2,A,E20.13,1X,D12.5)' )
                0348      &       'obsfit_cost(',num_file,') = ',
                0349      &       objf_obsfit_glo, num_obsfit_glo
                0350 
                0351         IF ( num_obsfit_glo.GT.zeroRL )
                0352      &    CALL PRINT_MESSAGE( msgBuf,
                0353      &         standardMessageUnit,SQUEEZE_RIGHT, myThid )
                0354 
                0355       ENDDO
                0356 
                0357       WRITE( msgBuf,'(A)') '== obsfit_cost: end   =='
                0358       CALL PRINT_MESSAGE( msgBuf,
                0359      &     standardMessageUnit,SQUEEZE_RIGHT, myThid )
                0360       WRITE( msgBuf,'(A)' ) ' '
                0361       CALL PRINT_MESSAGE( msgBuf,
                0362      &     standardMessageUnit,SQUEEZE_RIGHT, myThid )
                0363 
                0364 #endif /* ALLOW_OBSFIT */
                0365 
                0366       RETURN
                0367       END
                0368 
                0369 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!