Back to home page

MITgcm
 
 

    

File indexing completed on 2026-10-08 05:08:43 UTC

view on githubraw file Latest commit 27256b2d on 2026-10-07 21:35:21 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)
                0179 #else
7236d20aeb Patr*0180       fcref = fc
ec93986742 Patr*0181 #endif
db02d8eca8 Jean*0182       WRITE(msgBuf,'(A,1PE22.14)')
                0183      &         'grdchk reference fc: fcref       =', fcref
                0184       CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0185 
                0186       DO bj = myByLo(myThid), myByHi(myThid)
                0187         DO bi = myBxLo(myThid), myBxHi(myThid)
                0188           DO k = 1, Nr
                0189             DO j = jMin, jMax
                0190               DO i = iMin, iMax
                0191                 tmpplot1(i,j,k,bi,bj) = 0. _d 0
                0192                 tmpplot2(i,j,k,bi,bj) = 0. _d 0
                0193                 tmpplot3(i,j,k,bi,bj) = 0. _d 0
                0194               ENDDO
                0195             ENDDO
                0196           ENDDO
                0197         ENDDO
                0198       ENDDO
                0199 
                0200       IF ( useCentralDiff ) THEN
b7ff4d81ac Patr*0201          grdchk_epsfac = 2. _d 0
db02d8eca8 Jean*0202       ELSE
b7ff4d81ac Patr*0203          grdchk_epsfac = 1. _d 0
db02d8eca8 Jean*0204       ENDIF
b7ff4d81ac Patr*0205 
1052c30783 Jean*0206       WRITE(standardMessageUnit,'(A)')
b41168db92 Patr*0207      &    'grad-res -------------------------------'
1052c30783 Jean*0208       WRITE(standardMessageUnit,'(2a)')
e4b263335d Patr*0209      &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
db02d8eca8 Jean*0210      &    '       fc ref            fc + eps           fc - eps'
78a0e1cce7 Patr*0211 #ifdef ALLOW_TANGENTLINEAR_RUN
1052c30783 Jean*0212       WRITE(standardMessageUnit,'(2a)')
e4b263335d Patr*0213      &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
db02d8eca8 Jean*0214      &    '      tlm grad            fd grad          1 - fd/tlm'
78a0e1cce7 Patr*0215 #else
1052c30783 Jean*0216       WRITE(standardMessageUnit,'(2a)')
e4b263335d Patr*0217      &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
db02d8eca8 Jean*0218      &    '      adj grad            fd grad          1 - fd/adj'
78a0e1cce7 Patr*0219 #endif
                0220 
db02d8eca8 Jean*0221 C--   Compute the finite difference approximations.
                0222 C--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
                0223 C--   gradient checks.
2091ce7ee7 Patr*0224 
db02d8eca8 Jean*0225       IF ( nbeg .EQ. 0 ) CALL GRDCHK_GET_POSITION( myThid )
f81d465bd0 Patr*0226 
db02d8eca8 Jean*0227       DO icomp = nbeg, nend, nstep
2091ce7ee7 Patr*0228 
db02d8eca8 Jean*0229         adxxmemo  = 0.
                0230         ichknum = (icomp - nbeg)/nstep + 1
2091ce7ee7 Patr*0231 
db02d8eca8 Jean*0232         IF ( ichknum .LE. maxgrdchecks ) THEN
                0233           WRITE(msgBuf,'(A,I4,A)')
                0234      &      '====== Starts gradient-check number', ichknum,
                0235      &                                         ' (=ichknum) ======='
                0236           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
2091ce7ee7 Patr*0237 
db02d8eca8 Jean*0238 C--       Determine the location of icomp on the grid.
                0239           IF ( myProcId .EQ. grdchkwhichproc ) THEN
                0240             CALL grdchk_loc( icomp, ichknum,
7109a141b2 Patr*0241      &              icvrec, itile, jtile, layer, obcspos,
e4b263335d Patr*0242      &              itilepos, jtilepos, icglo, itest, ierr,
db02d8eca8 Jean*0243      &              myThid )
                0244           ELSE
                0245             icvrec   = 0
                0246             itile    = 0
                0247             jtile    = 0
                0248             layer    = 0
                0249             obcspos  = 0
                0250             itilepos = 0
                0251             jtilepos = 0
                0252             icglo    = 0
                0253             itest    = 0
                0254           ENDIF
                0255 
                0256 C make sure that all procs have correct file records, so that useSingleCpuIO works
                0257           CALL GLOBAL_SUM_INT( icvrec , myThid )
                0258           CALL GLOBAL_SUM_INT( layer , myThid )
                0259           CALL GLOBAL_SUM_INT( ierr , myThid )
                0260 
                0261           WRITE(msgBuf,'(A,3I5,A,2I4,A,I3,A,I4)')
                0262      &     'grdchk pos: i,j,k=', itilepos, jtilepos, layer,
                0263      &              ' ; bi,bj=', itile, jtile, ' ; iobc=', obcspos,
                0264      &              ' ; rec=', icvrec
                0265           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0266 
                0267 C******************************************************
                0268 C--   (A): get gradient component calculated via adjoint
                0269 C******************************************************
                0270 
                0271 C--   get gradient component calculated via adjoint
                0272           CALL GRDCHK_GETADXX( icvrec,
78a0e1cce7 Patr*0273      &              itilepos, jtilepos,
1052c30783 Jean*0274      &              layer, itile, jtile,
                0275      &              grdchkwhichproc, grdchkvarindex,
                0276      &              adxxmemo, myThid )
4a63303a97 Jean*0277 C--   Add a global-sum call so that all proc will get the adjoint gradient
db02d8eca8 Jean*0278           _GLOBAL_SUM_RL( adxxmemo, myThid )
b7ff4d81ac Patr*0279 
22f0d78f5f Patr*0280 #ifdef ALLOW_TANGENTLINEAR_RUN
db02d8eca8 Jean*0281 C******************************************************
                0282 C--   (B): Get gradient component g_fc from tangent linear run:
                0283 C******************************************************
                0284 
                0285 C--   1. perturb control vector component: xx(i)=1.
                0286 
                0287           localEps = 1. _d 0
                0288           CALL GRDCHK_GETXX( icvrec, TANGENT_SIMULATION,
78a0e1cce7 Patr*0289      &              itilepos, jtilepos,
1052c30783 Jean*0290      &              layer, itile, jtile,
                0291      &              grdchkwhichproc, grdchkvarindex,
78a0e1cce7 Patr*0292      &              xxmemo_ref, xxmemo_pert, localEps,
1052c30783 Jean*0293      &              myThid )
22f0d78f5f Patr*0294 
db02d8eca8 Jean*0295 C--   2. perform tangent linear run
                0296           myTime = startTime
                0297           myIter = nIter0
b4daa24319 Shre*0298 # ifdef ALLOW_ADMTLM
db02d8eca8 Jean*0299           DO k=1,4*Nr+1
                0300             DO j=1,sNy
                0301               DO i=1,sNx
                0302                 g_objf_state_final(i,j,1,1,k) = 0.
                0303               ENDDO
                0304             ENDDO
                0305           ENDDO
b4daa24319 Shre*0306 # else
db02d8eca8 Jean*0307           g_fc = 0.
b4daa24319 Shre*0308 # endif
                0309 
                0310 # ifdef ALLOW_TAPENADE
                0311           CALL THE_MAIN_LOOP_D( myTime, myIter, myThid )
                0312 #  ifdef ALLOW_ADMTLM
                0313           ftlxxmemo = objf_state_finald(idep,jdep,1,1,1)
                0314 #  else
                0315           ftlxxmemo = fcd
                0316 #  endif
                0317 # else /* ALLOW_TAPENADE */
db02d8eca8 Jean*0318           CALL G_THE_MAIN_LOOP( myTime, myIter, myThid )
b4daa24319 Shre*0319 #  ifdef ALLOW_ADMTLM
db02d8eca8 Jean*0320           ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
b4daa24319 Shre*0321 #  else
db02d8eca8 Jean*0322           ftlxxmemo = g_fc
b4daa24319 Shre*0323 #  endif
                0324 # endif /* ALLOW_TAPENADE */
ec93986742 Patr*0325 
db02d8eca8 Jean*0326 C--   3. reset control vector
                0327           CALL GRDCHK_SETXX( icvrec, TANGENT_SIMULATION,
78a0e1cce7 Patr*0328      &              itilepos, jtilepos,
1052c30783 Jean*0329      &              layer, itile, jtile,
                0330      &              grdchkwhichproc, grdchkvarindex,
                0331      &              xxmemo_ref, myThid )
78a0e1cce7 Patr*0332 
                0333 #endif /* ALLOW_TANGENTLINEAR_RUN */
                0334 
db02d8eca8 Jean*0335 C******************************************************
                0336 C--   (C): Get gradient via finite difference perturbation
                0337 C******************************************************
22f0d78f5f Patr*0338 
db02d8eca8 Jean*0339 C--   (C-1) positive perturbation:
                0340           localEps = ABS(grdchk_eps)
22f0d78f5f Patr*0341 
db02d8eca8 Jean*0342 C--   get control vector component from file
                0343 C--   perturb it and write back to file
                0344           CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
78a0e1cce7 Patr*0345      &              itilepos, jtilepos,
1052c30783 Jean*0346      &              layer, itile, jtile,
                0347      &              grdchkwhichproc, grdchkvarindex,
78a0e1cce7 Patr*0348      &              xxmemo_ref, xxmemo_pert, localEps,
1052c30783 Jean*0349      &              myThid )
3bb97d13a4 Jean*0350 
db02d8eca8 Jean*0351 C--   forward run with perturbed control vector
                0352           myTime = startTime
                0353           myIter = nIter0
                0354           CALL THE_MAIN_LOOP( myTime, myIter, myThid )
71dd48804f Patr*0355 
b4daa24319 Shre*0356 #ifdef ALLOW_ADMTLM
db02d8eca8 Jean*0357           fcpertplus = objf_state_final(idep,jdep,1,1,1)
ec93986742 Patr*0358 #else
db02d8eca8 Jean*0359           fcpertplus = fc
ec93986742 Patr*0360 #endif
db02d8eca8 Jean*0361           WRITE(msgBuf,'(A,1PE22.14)')
                0362      &         'grdchk perturb(+)fc: fcpertplus  =', fcpertplus
                0363           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
3bb97d13a4 Jean*0364 
db02d8eca8 Jean*0365 C--   Reset control vector.
                0366           CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
78a0e1cce7 Patr*0367      &              itilepos, jtilepos,
1052c30783 Jean*0368      &              layer, itile, jtile,
                0369      &              grdchkwhichproc, grdchkvarindex,
                0370      &              xxmemo_ref, myThid )
b7ff4d81ac Patr*0371 
db02d8eca8 Jean*0372           fcpertminus = fcref
                0373           IF ( useCentralDiff ) THEN
                0374 C--   (C-2) repeat the proceedure for a negative perturbation:
                0375             localEps = - ABS(grdchk_eps)
b7ff4d81ac Patr*0376 
db02d8eca8 Jean*0377 C--   get control vector component from file
                0378 C--   perturb it and write back to file
                0379             CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
1052c30783 Jean*0380      &              itilepos, jtilepos,
                0381      &              layer, itile, jtile,
                0382      &              grdchkwhichproc, grdchkvarindex,
                0383      &              xxmemo_ref, xxmemo_pert, localEps,
                0384      &              myThid )
3bb97d13a4 Jean*0385 
db02d8eca8 Jean*0386 C--   forward run with perturbed control vector
                0387             myTime = startTime
                0388             myIter = nIter0
                0389             CALL THE_MAIN_LOOP( myTime, myIter, myThid )
71dd48804f Patr*0390 
                0391 #if (defined (ALLOW_ADMTLM))
db02d8eca8 Jean*0392             fcpertminus = objf_state_final(idep,jdep,1,1,1)
ec93986742 Patr*0393 #else
db02d8eca8 Jean*0394             fcpertminus = fc
ec93986742 Patr*0395 #endif
db02d8eca8 Jean*0396             WRITE(msgBuf,'(A,1PE22.14)')
                0397      &         'grdchk perturb(-)fc: fcpertminus =', fcpertminus
                0398             CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
3bb97d13a4 Jean*0399 
db02d8eca8 Jean*0400 C--   Reset control vector.
                0401             CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
1052c30783 Jean*0402      &              itilepos, jtilepos,
                0403      &              layer, itile, jtile,
                0404      &              grdchkwhichproc, grdchkvarindex,
                0405      &              xxmemo_ref, myThid )
db02d8eca8 Jean*0406 
                0407 C-- end of if useCentralDiff ...
                0408           ENDIF
                0409 
                0410 C******************************************************
                0411 C--   (D): calculate relative differences between gradients
                0412 C******************************************************
                0413 
                0414           IF ( grdchk_eps .EQ. 0. ) THEN
                0415             gfd = (fcpertplus-fcpertminus)
                0416           ELSE
                0417             gfd = (fcpertplus-fcpertminus) /(grdchk_epsfac*grdchk_eps)
                0418           ENDIF
                0419 
                0420           IF ( adxxmemo .EQ. 0. ) THEN
                0421             ratio_ad = ABS( adxxmemo - gfd )
                0422           ELSE
                0423             ratio_ad = 1. - gfd/adxxmemo
                0424           ENDIF
                0425 
                0426           IF ( ftlxxmemo .EQ. 0. ) THEN
                0427             ratio_ftl = ABS( ftlxxmemo - gfd )
                0428           ELSE
                0429             ratio_ftl = 1. - gfd/ftlxxmemo
                0430           ENDIF
                0431 
                0432           IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
                0433             tmpplot1(itilepos,jtilepos,layer,itile,jtile) = gfd
                0434             tmpplot2(itilepos,jtilepos,layer,itile,jtile) = ratio_ad
                0435             tmpplot3(itilepos,jtilepos,layer,itile,jtile) = ratio_ftl
                0436           ENDIF
                0437 
                0438           IF ( ierr .EQ. 0 ) THEN
                0439             fcrmem      ( ichknum ) = fcref
                0440             fcppmem     ( ichknum ) = fcpertplus
                0441             fcpmmem     ( ichknum ) = fcpertminus
                0442             xxmemref    ( ichknum ) = xxmemo_ref
                0443             xxmempert   ( ichknum ) = xxmemo_pert
                0444             gfdmem      ( ichknum ) = gfd
                0445             adxxmem     ( ichknum ) = adxxmemo
                0446             ftlxxmem    ( ichknum ) = ftlxxmemo
                0447             ratioadmem  ( ichknum ) = ratio_ad
                0448             ratioftlmem ( ichknum ) = ratio_ftl
                0449 
                0450             irecmem   ( ichknum ) = icvrec
                0451             bimem     ( ichknum ) = itile
                0452             bjmem     ( ichknum ) = jtile
                0453             ilocmem   ( ichknum ) = itilepos
                0454             jlocmem   ( ichknum ) = jtilepos
                0455             klocmem   ( ichknum ) = layer
                0456             iobcsmem  ( ichknum ) = obcspos
                0457             icompmem  ( ichknum ) = icomp
                0458             ichkmem   ( ichknum ) = ichknum
                0459             itestmem  ( ichknum ) = itest
                0460             ierrmem   ( ichknum ) = ierr
                0461             icglomem  ( ichknum ) = icglo
                0462           ENDIF
                0463 
                0464           IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
1052c30783 Jean*0465             WRITE(standardMessageUnit,'(A)')
db02d8eca8 Jean*0466      &          'grad-res -------------------------------'
1052c30783 Jean*0467             WRITE(standardMessageUnit,'(A,8I5,1x,1P3E19.11)')
db02d8eca8 Jean*0468      &           ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
                0469      &             layer,itile,jtile,obcspos,
                0470      &             fcref, fcpertplus, fcpertminus
78a0e1cce7 Patr*0471 #ifdef ALLOW_TANGENTLINEAR_RUN
1052c30783 Jean*0472             WRITE(standardMessageUnit,'(A,8I5,1x,1P3E19.11)')
db02d8eca8 Jean*0473      &           ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
                0474      &             icompmem(ichknum),itestmem(ichknum),
                0475      &             bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
                0476      &             ftlxxmemo, gfd, ratio_ftl
78a0e1cce7 Patr*0477 #else
1052c30783 Jean*0478             WRITE(standardMessageUnit,'(A,8I5,1x,1P3E19.11)')
db02d8eca8 Jean*0479      &           ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
                0480      &             icompmem(ichknum),itestmem(ichknum),
                0481      &             bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
                0482      &             adxxmemo, gfd, ratio_ad
4a63303a97 Jean*0483 #endif
db02d8eca8 Jean*0484           ENDIF
4a63303a97 Jean*0485 #ifdef ALLOW_TANGENTLINEAR_RUN
db02d8eca8 Jean*0486           WRITE(msgBuf,'(A30,1PE22.14)')
ab2fb9a209 Jean*0487      &              ' TLM  ref_cost_function      =', fcref
db02d8eca8 Jean*0488           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0489           WRITE(msgBuf,'(A30,1PE22.14)')
ab2fb9a209 Jean*0490      &              ' TLM  tangent-lin_grad       =', ftlxxmemo
db02d8eca8 Jean*0491           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0492           WRITE(msgBuf,'(A30,1PE22.14)')
ab2fb9a209 Jean*0493      &              ' TLM  finite-diff_grad       =', gfd
db02d8eca8 Jean*0494           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
4a63303a97 Jean*0495 #else
db02d8eca8 Jean*0496           WRITE(msgBuf,'(A30,1PE22.14)')
84cbbbce68 Jean*0497      &              ' ADM  ref_cost_function      =', fcref
db02d8eca8 Jean*0498           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0499           WRITE(msgBuf,'(A30,1PE22.14)')
84cbbbce68 Jean*0500      &              ' ADM  adjoint_gradient       =', adxxmemo
db02d8eca8 Jean*0501           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
                0502           WRITE(msgBuf,'(A30,1PE22.14)')
84cbbbce68 Jean*0503      &              ' ADM  finite-diff_grad       =', gfd
db02d8eca8 Jean*0504           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
78a0e1cce7 Patr*0505 #endif
                0506 
db02d8eca8 Jean*0507           WRITE(msgBuf,'(A,I4,A,I3,A)')
                0508      &      '====== End of gradient-check number', ichknum,
                0509      &                      ' (ierr=', ierr, ') ======='
                0510           CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
78a0e1cce7 Patr*0511 
db02d8eca8 Jean*0512 C-- else of if ichknum < ...
                0513         ELSE
                0514           ierr_grdchk = -1
78a0e1cce7 Patr*0515 
db02d8eca8 Jean*0516 C-- end of if ichknum < ...
                0517         ENDIF
2091ce7ee7 Patr*0518 
db02d8eca8 Jean*0519 C-- end of do icomp = ...
                0520       ENDDO
78a0e1cce7 Patr*0521 
9f85ea262e Mart*0522       CALL WRITE_REC_XYZ_RL( 'grd_findiff'  , tmpplot1, 1, 0, myThid )
                0523 # ifdef ALLOW_ADJOINT_RUN
                0524       CALL WRITE_REC_XYZ_RL( 'grd_ratio_ad' , tmpplot2, 1, 0, myThid )
                0525 # endif
                0526 # ifdef  ALLOW_TANGENTLINEAR_RUN
                0527       CALL WRITE_REC_XYZ_RL( 'grd_ratio_ftl', tmpplot3, 1, 0, myThid )
                0528 # endif
2091ce7ee7 Patr*0529 
db02d8eca8 Jean*0530 C--   Print the results of the gradient check.
                0531       CALL GRDCHK_PRINT( ichknum, ierr_grdchk, myThid )
2091ce7ee7 Patr*0532 
edd57506ae Patr*0533 #endif /* ALLOW_GRDCHK */
2091ce7ee7 Patr*0534 
db02d8eca8 Jean*0535       RETURN
                0536       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!