Back to home page

MITgcm
 
 

    

File indexing completed on 2024-03-02 06:10:37 UTC

view on githubraw file Latest commit 5cf43646 on 2024-03-01 18:50:49 UTC
3bb97d13a4 Jean*0001 #include "GRDCHK_OPTIONS.h"
57c22ecc45 Jean*0002 #include "AD_CONFIG.h"
e4ef0556fe Jean*0003 #ifdef ALLOW_COST
                0004 # include "COST_OPTIONS.h"
                0005 #endif
                0006 #ifdef ALLOW_CTRL
                0007 # include "CTRL_OPTIONS.h"
                0008 #endif
2091ce7ee7 Patr*0009 
d10bc38b8a Patr*0010 CBOI
                0011 C
                0012 C !TITLE: GRADIENT CHECK
                0013 C !AUTHORS: mitgcm developers ( support@mitgcm.org )
                0014 C !AFFILIATION: Massachussetts Institute of Technology
                0015 C !DATE:
                0016 C !INTRODUCTION: gradient check package
db02d8eca8 Jean*0017 C \bv
                0018 C Compare the gradients calculated by the adjoint model with
                0019 C finite difference approximations.
                0020 C
d10bc38b8a Patr*0021 C     !CALLING SEQUENCE:
db02d8eca8 Jean*0022 C
                0023 C the_model_main
                0024 C |
                0025 C |-- ctrl_unpack
                0026 C |-- adthe_main_loop          - unperturbed cost function and
                0027 C |-- ctrl_pack                  adjoint gradient are computed here
                0028 C |
                0029 C |-- grdchk_main
                0030 C     |
1052c30783 Jean*0031 C     |-- grdchk_get_mask
db02d8eca8 Jean*0032 C     |-- do icomp=...        - loop over control vector elements
                0033 C         |
                0034 C         |-- grdchk_loc      - determine location of icomp on grid
                0035 C         |
                0036 C         |-- grdchk_getadxx  - get gradient component calculated
                0037 C         |                     via adjoint
                0038 C         |-- grdchk_getxx    - get control vector component from file
                0039 C         |                     perturb it and write back to file
                0040 C         |-- the_main_loop   - forward run and cost evaluation
                0041 C         |                     with perturbed control vector element
                0042 C         |-- calculate ratio of adj. vs. finite difference gradient
                0043 C         |
                0044 C         |-- grdchk_setxx    - Reset control vector element
                0045 C     |
                0046 C     |-- grdchk_print    - print results
                0047 C \ev
d10bc38b8a Patr*0048 CEOI
2091ce7ee7 Patr*0049 
d10bc38b8a Patr*0050 CBOP
db02d8eca8 Jean*0051 C     !ROUTINE: GRDCHK_MAIN
d10bc38b8a Patr*0052 C     !INTERFACE:
db02d8eca8 Jean*0053       SUBROUTINE GRDCHK_MAIN( myThid )
2091ce7ee7 Patr*0054 
d10bc38b8a Patr*0055 C     !DESCRIPTION: \bv
db02d8eca8 Jean*0056 C     ==================================================================
                0057 C     SUBROUTINE GRDCHK_MAIN
                0058 C     ==================================================================
                0059 C     o Compare the gradients calculated by the adjoint model with
                0060 C       finite difference approximations.
                0061 C     started: Christian Eckert eckert@mit.edu 24-Feb-2000
                0062 C     continued&finished: heimbach@mit.edu: 13-Jun-2001
                0063 C     changed: mlosch@ocean.mit.edu: 09-May-2002
                0064 C              - added centered difference vs. 1-sided difference option
                0065 C              - improved output format for readability
                0066 C              - added control variable hFacC
                0067 C              heimbach@mit.edu 24-Feb-2003
                0068 C              - added tangent linear gradient checks
                0069 C              - fixes for multiproc. gradient checks
                0070 C              - added more control variables
                0071 C     ==================================================================
                0072 C     SUBROUTINE GRDCHK_MAIN
                0073 C     ==================================================================
d10bc38b8a Patr*0074 C     \ev
2091ce7ee7 Patr*0075 
d10bc38b8a Patr*0076 C     !USES:
db02d8eca8 Jean*0077       IMPLICIT NONE
2091ce7ee7 Patr*0078 
db02d8eca8 Jean*0079 C     == global variables ==
2091ce7ee7 Patr*0080 #include "SIZE.h"
                0081 #include "EEPARAMS.h"
                0082 #include "PARAMS.h"
                0083 #include "cost.h"
7236d20aeb Patr*0084 #include "g_cost.h"
5cf4364659 Mart*0085 #include "CTRL_SIZE.h"
4d72283393 Mart*0086 #include "CTRL.h"
444da61630 Mart*0087 #ifdef ALLOW_OBCS_CONTROL
                0088 C     CTRL_OBCS.h must be included before GRDCHK.h
                0089 # include "CTRL_OBCS.h"
                0090 #endif
                0091 #include "GRDCHK.h"
2091ce7ee7 Patr*0092 
b4daa24319 Shre*0093 #ifdef ALLOW_TAPENADE
                0094 # include "COST_TAP_TLM.h"
                0095 #endif
                0096 
d10bc38b8a Patr*0097 C     !INPUT/OUTPUT PARAMETERS:
db02d8eca8 Jean*0098       INTEGER myThid
2091ce7ee7 Patr*0099 
edd57506ae Patr*0100 #ifdef ALLOW_GRDCHK
d10bc38b8a Patr*0101 C     !LOCAL VARIABLES:
db02d8eca8 Jean*0102       INTEGER myIter
                0103       _RL     myTime
                0104       INTEGER bi, bj
                0105       INTEGER i, j, k
                0106       INTEGER iMin, iMax, jMin, jMax
                0107       PARAMETER( iMin = 1 , iMax = sNx , jMin = 1 , jMax = sNy )
                0108       INTEGER ioUnit
                0109       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0110 
                0111       INTEGER icomp
                0112       INTEGER ichknum
                0113       INTEGER icvrec
                0114       INTEGER jtile
                0115       INTEGER itile
                0116       INTEGER layer
                0117       INTEGER obcspos
                0118       INTEGER itilepos
                0119       INTEGER jtilepos
                0120       INTEGER icglo
                0121       INTEGER itest
                0122       INTEGER ierr
                0123       INTEGER ierr_grdchk
2091ce7ee7 Patr*0124       _RL     gfd
                0125       _RL     fcref
b7ff4d81ac Patr*0126       _RL     fcpertplus, fcpertminus
22f0d78f5f Patr*0127       _RL     ratio_ad
                0128       _RL     ratio_ftl
2091ce7ee7 Patr*0129       _RL     xxmemo_ref
                0130       _RL     xxmemo_pert
                0131       _RL     adxxmemo
22f0d78f5f Patr*0132       _RL     ftlxxmemo
                0133       _RL     localEps
                0134       _RL     grdchk_epsfac
db02d8eca8 Jean*0135       _RL tmpplot1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0136       _RL tmpplot2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0137       _RL tmpplot3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
d10bc38b8a Patr*0138 CEOP
2091ce7ee7 Patr*0139 
db02d8eca8 Jean*0140       ioUnit = standardMessageUnit
                0141       WRITE(msgBuf,'(A)')
                0142      &'// ======================================================='
                0143       CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0144       WRITE(msgBuf,'(A)') '// Gradient-check starts (grdchk_main)'
                0145       CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0146       WRITE(msgBuf,'(A)')
                0147      &'// ======================================================='
                0148       CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
2091ce7ee7 Patr*0149 
db02d8eca8 Jean*0150 #ifdef ALLOW_TANGENTLINEAR_RUN
                0151 C--   prevent writing output multiple times
                0152 C     note: already called in AD run so that only needed for TLM run
                0153       CALL TURNOFF_MODEL_IO( 0, myThid )
                0154 #endif
d7ee8fe52e Patr*0155 
db02d8eca8 Jean*0156 C--   initialise variables
1052c30783 Jean*0157       CALL GRDCHK_GET_MASK( myThid )
2091ce7ee7 Patr*0158 
db02d8eca8 Jean*0159 C--   Compute the adjoint model gradients.
                0160 C--   Compute the unperturbed cost function.
                0161 Cph   Gradient via adjoint has already been computed,
                0162 Cph   and so has unperturbed cost function,
                0163 Cph   assuming all xx_ fields are initialised to zero.
2091ce7ee7 Patr*0164 
4b1e62e6df Patr*0165       ierr      = 0
2091ce7ee7 Patr*0166       ierr_grdchk = 0
6f2c42d0ea Patr*0167       adxxmemo  = 0.
                0168       ftlxxmemo = 0.
b4daa24319 Shre*0169 
                0170 #ifdef ALLOW_TAPENADE
                0171 C--   forward run with unperturbed control vector
                0172       myTime = startTime
                0173       myIter = nIter0
                0174       CALL THE_MAIN_LOOP( myTime, myIter, myThid )
                0175 #endif
                0176 
71dd48804f Patr*0177 #if (defined  (ALLOW_ADMTLM))
ec93986742 Patr*0178       fcref = objf_state_final(idep,jdep,1,1,1)
e4ef0556fe Jean*0179 #elif (defined (ALLOW_OPENAD))
71dd48804f Patr*0180       fcref = fc%v
ec93986742 Patr*0181 #else
7236d20aeb Patr*0182       fcref = fc
ec93986742 Patr*0183 #endif
db02d8eca8 Jean*0184       WRITE(msgBuf,'(A,1PE22.14)')
                0185      &         'grdchk reference fc: fcref       =', fcref
                0186       CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0187 
                0188       DO bj = myByLo(myThid), myByHi(myThid)
                0189         DO bi = myBxLo(myThid), myBxHi(myThid)
                0190           DO k = 1, Nr
                0191             DO j = jMin, jMax
                0192               DO i = iMin, iMax
                0193                 tmpplot1(i,j,k,bi,bj) = 0. _d 0
                0194                 tmpplot2(i,j,k,bi,bj) = 0. _d 0
                0195                 tmpplot3(i,j,k,bi,bj) = 0. _d 0
                0196               ENDDO
                0197             ENDDO
                0198           ENDDO
                0199         ENDDO
                0200       ENDDO
                0201 
                0202       IF ( useCentralDiff ) THEN
b7ff4d81ac Patr*0203          grdchk_epsfac = 2. _d 0
db02d8eca8 Jean*0204       ELSE
b7ff4d81ac Patr*0205          grdchk_epsfac = 1. _d 0
db02d8eca8 Jean*0206       ENDIF
b7ff4d81ac Patr*0207 
1052c30783 Jean*0208       WRITE(standardMessageUnit,'(A)')
b41168db92 Patr*0209      &    'grad-res -------------------------------'
1052c30783 Jean*0210       WRITE(standardMessageUnit,'(2a)')
e4b263335d Patr*0211      &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
db02d8eca8 Jean*0212      &    '       fc ref            fc + eps           fc - eps'
78a0e1cce7 Patr*0213 #ifdef ALLOW_TANGENTLINEAR_RUN
1052c30783 Jean*0214       WRITE(standardMessageUnit,'(2a)')
e4b263335d Patr*0215      &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
db02d8eca8 Jean*0216      &    '      tlm grad            fd grad          1 - fd/tlm'
78a0e1cce7 Patr*0217 #else
1052c30783 Jean*0218       WRITE(standardMessageUnit,'(2a)')
e4b263335d Patr*0219      &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
db02d8eca8 Jean*0220      &    '      adj grad            fd grad          1 - fd/adj'
78a0e1cce7 Patr*0221 #endif
                0222 
db02d8eca8 Jean*0223 C--   Compute the finite difference approximations.
                0224 C--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
                0225 C--   gradient checks.
2091ce7ee7 Patr*0226 
db02d8eca8 Jean*0227       IF ( nbeg .EQ. 0 ) CALL GRDCHK_GET_POSITION( myThid )
f81d465bd0 Patr*0228 
db02d8eca8 Jean*0229       DO icomp = nbeg, nend, nstep
2091ce7ee7 Patr*0230 
db02d8eca8 Jean*0231         adxxmemo  = 0.
                0232         ichknum = (icomp - nbeg)/nstep + 1
2091ce7ee7 Patr*0233 
db02d8eca8 Jean*0234         IF ( ichknum .LE. maxgrdchecks ) THEN
                0235           WRITE(msgBuf,'(A,I4,A)')
                0236      &      '====== Starts gradient-check number', ichknum,
                0237      &                                         ' (=ichknum) ======='
                0238           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
2091ce7ee7 Patr*0239 
db02d8eca8 Jean*0240 C--       Determine the location of icomp on the grid.
                0241           IF ( myProcId .EQ. grdchkwhichproc ) THEN
                0242             CALL grdchk_loc( icomp, ichknum,
7109a141b2 Patr*0243      &              icvrec, itile, jtile, layer, obcspos,
e4b263335d Patr*0244      &              itilepos, jtilepos, icglo, itest, ierr,
db02d8eca8 Jean*0245      &              myThid )
                0246           ELSE
                0247             icvrec   = 0
                0248             itile    = 0
                0249             jtile    = 0
                0250             layer    = 0
                0251             obcspos  = 0
                0252             itilepos = 0
                0253             jtilepos = 0
                0254             icglo    = 0
                0255             itest    = 0
                0256           ENDIF
                0257 
                0258 C make sure that all procs have correct file records, so that useSingleCpuIO works
                0259           CALL GLOBAL_SUM_INT( icvrec , myThid )
                0260           CALL GLOBAL_SUM_INT( layer , myThid )
                0261           CALL GLOBAL_SUM_INT( ierr , myThid )
                0262 
                0263           WRITE(msgBuf,'(A,3I5,A,2I4,A,I3,A,I4)')
                0264      &     'grdchk pos: i,j,k=', itilepos, jtilepos, layer,
                0265      &              ' ; bi,bj=', itile, jtile, ' ; iobc=', obcspos,
                0266      &              ' ; rec=', icvrec
                0267           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0268 
                0269 C******************************************************
                0270 C--   (A): get gradient component calculated via adjoint
                0271 C******************************************************
                0272 
                0273 C--   get gradient component calculated via adjoint
                0274           CALL GRDCHK_GETADXX( icvrec,
78a0e1cce7 Patr*0275      &              itilepos, jtilepos,
1052c30783 Jean*0276      &              layer, itile, jtile,
                0277      &              grdchkwhichproc, grdchkvarindex,
                0278      &              adxxmemo, myThid )
4a63303a97 Jean*0279 C--   Add a global-sum call so that all proc will get the adjoint gradient
db02d8eca8 Jean*0280           _GLOBAL_SUM_RL( adxxmemo, myThid )
b7ff4d81ac Patr*0281 
22f0d78f5f Patr*0282 #ifdef ALLOW_TANGENTLINEAR_RUN
db02d8eca8 Jean*0283 C******************************************************
                0284 C--   (B): Get gradient component g_fc from tangent linear run:
                0285 C******************************************************
                0286 
                0287 C--   1. perturb control vector component: xx(i)=1.
                0288 
                0289           localEps = 1. _d 0
                0290           CALL GRDCHK_GETXX( icvrec, TANGENT_SIMULATION,
78a0e1cce7 Patr*0291      &              itilepos, jtilepos,
1052c30783 Jean*0292      &              layer, itile, jtile,
                0293      &              grdchkwhichproc, grdchkvarindex,
78a0e1cce7 Patr*0294      &              xxmemo_ref, xxmemo_pert, localEps,
1052c30783 Jean*0295      &              myThid )
22f0d78f5f Patr*0296 
db02d8eca8 Jean*0297 C--   2. perform tangent linear run
                0298           myTime = startTime
                0299           myIter = nIter0
b4daa24319 Shre*0300 # ifdef ALLOW_ADMTLM
db02d8eca8 Jean*0301           DO k=1,4*Nr+1
                0302             DO j=1,sNy
                0303               DO i=1,sNx
                0304                 g_objf_state_final(i,j,1,1,k) = 0.
                0305               ENDDO
                0306             ENDDO
                0307           ENDDO
b4daa24319 Shre*0308 # else
db02d8eca8 Jean*0309           g_fc = 0.
b4daa24319 Shre*0310 # endif
                0311 
                0312 # ifdef ALLOW_TAPENADE
                0313           CALL THE_MAIN_LOOP_D( myTime, myIter, myThid )
                0314 #  ifdef ALLOW_ADMTLM
                0315           ftlxxmemo = objf_state_finald(idep,jdep,1,1,1)
                0316 #  else
                0317           ftlxxmemo = fcd
                0318 #  endif
                0319 # else /* ALLOW_TAPENADE */
db02d8eca8 Jean*0320           CALL G_THE_MAIN_LOOP( myTime, myIter, myThid )
b4daa24319 Shre*0321 #  ifdef ALLOW_ADMTLM
db02d8eca8 Jean*0322           ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
b4daa24319 Shre*0323 #  else
db02d8eca8 Jean*0324           ftlxxmemo = g_fc
b4daa24319 Shre*0325 #  endif
                0326 # endif /* ALLOW_TAPENADE */
ec93986742 Patr*0327 
db02d8eca8 Jean*0328 C--   3. reset control vector
                0329           CALL GRDCHK_SETXX( icvrec, TANGENT_SIMULATION,
78a0e1cce7 Patr*0330      &              itilepos, jtilepos,
1052c30783 Jean*0331      &              layer, itile, jtile,
                0332      &              grdchkwhichproc, grdchkvarindex,
                0333      &              xxmemo_ref, myThid )
78a0e1cce7 Patr*0334 
                0335 #endif /* ALLOW_TANGENTLINEAR_RUN */
                0336 
db02d8eca8 Jean*0337 C******************************************************
                0338 C--   (C): Get gradient via finite difference perturbation
                0339 C******************************************************
22f0d78f5f Patr*0340 
db02d8eca8 Jean*0341 C--   (C-1) positive perturbation:
                0342           localEps = ABS(grdchk_eps)
22f0d78f5f Patr*0343 
db02d8eca8 Jean*0344 C--   get control vector component from file
                0345 C--   perturb it and write back to file
                0346           CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
78a0e1cce7 Patr*0347      &              itilepos, jtilepos,
1052c30783 Jean*0348      &              layer, itile, jtile,
                0349      &              grdchkwhichproc, grdchkvarindex,
78a0e1cce7 Patr*0350      &              xxmemo_ref, xxmemo_pert, localEps,
1052c30783 Jean*0351      &              myThid )
3bb97d13a4 Jean*0352 
db02d8eca8 Jean*0353 C--   forward run with perturbed control vector
                0354           myTime = startTime
                0355           myIter = nIter0
e4ef0556fe Jean*0356 #ifdef ALLOW_OPENAD
db02d8eca8 Jean*0357           CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
71dd48804f Patr*0358 #else
db02d8eca8 Jean*0359           CALL THE_MAIN_LOOP( myTime, myIter, myThid )
71dd48804f Patr*0360 #endif
                0361 
b4daa24319 Shre*0362 #ifdef ALLOW_ADMTLM
db02d8eca8 Jean*0363           fcpertplus = objf_state_final(idep,jdep,1,1,1)
e4ef0556fe Jean*0364 #elif (defined (ALLOW_OPENAD))
db02d8eca8 Jean*0365           fcpertplus = fc%v
ec93986742 Patr*0366 #else
db02d8eca8 Jean*0367           fcpertplus = fc
ec93986742 Patr*0368 #endif
db02d8eca8 Jean*0369           WRITE(msgBuf,'(A,1PE22.14)')
                0370      &         'grdchk perturb(+)fc: fcpertplus  =', fcpertplus
                0371           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
3bb97d13a4 Jean*0372 
db02d8eca8 Jean*0373 C--   Reset control vector.
                0374           CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
78a0e1cce7 Patr*0375      &              itilepos, jtilepos,
1052c30783 Jean*0376      &              layer, itile, jtile,
                0377      &              grdchkwhichproc, grdchkvarindex,
                0378      &              xxmemo_ref, myThid )
b7ff4d81ac Patr*0379 
db02d8eca8 Jean*0380           fcpertminus = fcref
                0381           IF ( useCentralDiff ) THEN
                0382 C--   (C-2) repeat the proceedure for a negative perturbation:
                0383             localEps = - ABS(grdchk_eps)
b7ff4d81ac Patr*0384 
db02d8eca8 Jean*0385 C--   get control vector component from file
                0386 C--   perturb it and write back to file
                0387             CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
1052c30783 Jean*0388      &              itilepos, jtilepos,
                0389      &              layer, itile, jtile,
                0390      &              grdchkwhichproc, grdchkvarindex,
                0391      &              xxmemo_ref, xxmemo_pert, localEps,
                0392      &              myThid )
3bb97d13a4 Jean*0393 
db02d8eca8 Jean*0394 C--   forward run with perturbed control vector
                0395             myTime = startTime
                0396             myIter = nIter0
e4ef0556fe Jean*0397 #ifdef ALLOW_OPENAD
db02d8eca8 Jean*0398             CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
71dd48804f Patr*0399 #else
db02d8eca8 Jean*0400             CALL THE_MAIN_LOOP( myTime, myIter, myThid )
71dd48804f Patr*0401 #endif
                0402 
                0403 #if (defined (ALLOW_ADMTLM))
db02d8eca8 Jean*0404             fcpertminus = objf_state_final(idep,jdep,1,1,1)
e4ef0556fe Jean*0405 #elif (defined (ALLOW_OPENAD))
db02d8eca8 Jean*0406             fcpertminus = fc%v
ec93986742 Patr*0407 #else
db02d8eca8 Jean*0408             fcpertminus = fc
ec93986742 Patr*0409 #endif
db02d8eca8 Jean*0410             WRITE(msgBuf,'(A,1PE22.14)')
                0411      &         'grdchk perturb(-)fc: fcpertminus =', fcpertminus
                0412             CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
3bb97d13a4 Jean*0413 
db02d8eca8 Jean*0414 C--   Reset control vector.
                0415             CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
1052c30783 Jean*0416      &              itilepos, jtilepos,
                0417      &              layer, itile, jtile,
                0418      &              grdchkwhichproc, grdchkvarindex,
                0419      &              xxmemo_ref, myThid )
db02d8eca8 Jean*0420 
                0421 C-- end of if useCentralDiff ...
                0422           ENDIF
                0423 
                0424 C******************************************************
                0425 C--   (D): calculate relative differences between gradients
                0426 C******************************************************
                0427 
                0428           IF ( grdchk_eps .EQ. 0. ) THEN
                0429             gfd = (fcpertplus-fcpertminus)
                0430           ELSE
                0431             gfd = (fcpertplus-fcpertminus) /(grdchk_epsfac*grdchk_eps)
                0432           ENDIF
                0433 
                0434           IF ( adxxmemo .EQ. 0. ) THEN
                0435             ratio_ad = ABS( adxxmemo - gfd )
                0436           ELSE
                0437             ratio_ad = 1. - gfd/adxxmemo
                0438           ENDIF
                0439 
                0440           IF ( ftlxxmemo .EQ. 0. ) THEN
                0441             ratio_ftl = ABS( ftlxxmemo - gfd )
                0442           ELSE
                0443             ratio_ftl = 1. - gfd/ftlxxmemo
                0444           ENDIF
                0445 
                0446           IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
                0447             tmpplot1(itilepos,jtilepos,layer,itile,jtile) = gfd
                0448             tmpplot2(itilepos,jtilepos,layer,itile,jtile) = ratio_ad
                0449             tmpplot3(itilepos,jtilepos,layer,itile,jtile) = ratio_ftl
                0450           ENDIF
                0451 
                0452           IF ( ierr .EQ. 0 ) THEN
                0453             fcrmem      ( ichknum ) = fcref
                0454             fcppmem     ( ichknum ) = fcpertplus
                0455             fcpmmem     ( ichknum ) = fcpertminus
                0456             xxmemref    ( ichknum ) = xxmemo_ref
                0457             xxmempert   ( ichknum ) = xxmemo_pert
                0458             gfdmem      ( ichknum ) = gfd
                0459             adxxmem     ( ichknum ) = adxxmemo
                0460             ftlxxmem    ( ichknum ) = ftlxxmemo
                0461             ratioadmem  ( ichknum ) = ratio_ad
                0462             ratioftlmem ( ichknum ) = ratio_ftl
                0463 
                0464             irecmem   ( ichknum ) = icvrec
                0465             bimem     ( ichknum ) = itile
                0466             bjmem     ( ichknum ) = jtile
                0467             ilocmem   ( ichknum ) = itilepos
                0468             jlocmem   ( ichknum ) = jtilepos
                0469             klocmem   ( ichknum ) = layer
                0470             iobcsmem  ( ichknum ) = obcspos
                0471             icompmem  ( ichknum ) = icomp
                0472             ichkmem   ( ichknum ) = ichknum
                0473             itestmem  ( ichknum ) = itest
                0474             ierrmem   ( ichknum ) = ierr
                0475             icglomem  ( ichknum ) = icglo
                0476           ENDIF
                0477 
                0478           IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
1052c30783 Jean*0479             WRITE(standardMessageUnit,'(A)')
db02d8eca8 Jean*0480      &          'grad-res -------------------------------'
1052c30783 Jean*0481             WRITE(standardMessageUnit,'(A,8I5,1x,1P3E19.11)')
db02d8eca8 Jean*0482      &           ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
                0483      &             layer,itile,jtile,obcspos,
                0484      &             fcref, fcpertplus, fcpertminus
78a0e1cce7 Patr*0485 #ifdef ALLOW_TANGENTLINEAR_RUN
1052c30783 Jean*0486             WRITE(standardMessageUnit,'(A,8I5,1x,1P3E19.11)')
db02d8eca8 Jean*0487      &           ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
                0488      &             icompmem(ichknum),itestmem(ichknum),
                0489      &             bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
                0490      &             ftlxxmemo, gfd, ratio_ftl
78a0e1cce7 Patr*0491 #else
1052c30783 Jean*0492             WRITE(standardMessageUnit,'(A,8I5,1x,1P3E19.11)')
db02d8eca8 Jean*0493      &           ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
                0494      &             icompmem(ichknum),itestmem(ichknum),
                0495      &             bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
                0496      &             adxxmemo, gfd, ratio_ad
4a63303a97 Jean*0497 #endif
db02d8eca8 Jean*0498           ENDIF
4a63303a97 Jean*0499 #ifdef ALLOW_TANGENTLINEAR_RUN
db02d8eca8 Jean*0500           WRITE(msgBuf,'(A30,1PE22.14)')
ab2fb9a209 Jean*0501      &              ' TLM  ref_cost_function      =', fcref
db02d8eca8 Jean*0502           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0503           WRITE(msgBuf,'(A30,1PE22.14)')
ab2fb9a209 Jean*0504      &              ' TLM  tangent-lin_grad       =', ftlxxmemo
db02d8eca8 Jean*0505           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0506           WRITE(msgBuf,'(A30,1PE22.14)')
ab2fb9a209 Jean*0507      &              ' TLM  finite-diff_grad       =', gfd
db02d8eca8 Jean*0508           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
4a63303a97 Jean*0509 #else
db02d8eca8 Jean*0510           WRITE(msgBuf,'(A30,1PE22.14)')
84cbbbce68 Jean*0511      &              ' ADM  ref_cost_function      =', fcref
db02d8eca8 Jean*0512           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0513           WRITE(msgBuf,'(A30,1PE22.14)')
84cbbbce68 Jean*0514      &              ' ADM  adjoint_gradient       =', adxxmemo
db02d8eca8 Jean*0515           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0516           WRITE(msgBuf,'(A30,1PE22.14)')
84cbbbce68 Jean*0517      &              ' ADM  finite-diff_grad       =', gfd
db02d8eca8 Jean*0518           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
78a0e1cce7 Patr*0519 #endif
                0520 
db02d8eca8 Jean*0521           WRITE(msgBuf,'(A,I4,A,I3,A)')
                0522      &      '====== End of gradient-check number', ichknum,
                0523      &                      ' (ierr=', ierr, ') ======='
                0524           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
78a0e1cce7 Patr*0525 
db02d8eca8 Jean*0526 C-- else of if ichknum < ...
                0527         ELSE
                0528           ierr_grdchk = -1
78a0e1cce7 Patr*0529 
db02d8eca8 Jean*0530 C-- end of if ichknum < ...
                0531         ENDIF
2091ce7ee7 Patr*0532 
db02d8eca8 Jean*0533 C-- end of do icomp = ...
                0534       ENDDO
78a0e1cce7 Patr*0535 
db02d8eca8 Jean*0536       IF (myProcId .EQ. grdchkwhichproc .AND. .NOT.useSingleCpuIO) THEN
                0537         CALL WRITE_REC_XYZ_RL(
                0538      &       'grd_findiff'   , tmpplot1, 1, 0, myThid)
                0539         CALL WRITE_REC_XYZ_RL(
                0540      &       'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
                0541         CALL WRITE_REC_XYZ_RL(
                0542      &       'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
                0543       ENDIF
2091ce7ee7 Patr*0544 
db02d8eca8 Jean*0545 C--   Print the results of the gradient check.
                0546       CALL GRDCHK_PRINT( ichknum, ierr_grdchk, myThid )
2091ce7ee7 Patr*0547 
edd57506ae Patr*0548 #endif /* ALLOW_GRDCHK */
2091ce7ee7 Patr*0549 
db02d8eca8 Jean*0550       RETURN
                0551       END
[ 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!