pkg/grdchk/grdchk_main.F

C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.33 2011/09/26 12:56:52 jmc Exp $
C $Name:  $

#include "GRDCHK_OPTIONS.h"

CBOI
C
C !TITLE: GRADIENT CHECK
C !AUTHORS: mitgcm developers ( support@mitgcm.org )
C !AFFILIATION: Massachussetts Institute of Technology
C !DATE:
C !INTRODUCTION: gradient check package
c \bv
c Compare the gradients calculated by the adjoint model with
c finite difference approximations.
c
C     !CALLING SEQUENCE:
c
c the_model_main
c |
c |-- ctrl_unpack
c |-- adthe_main_loop          - unperturbed cost function and
c |-- ctrl_pack                  adjoint gradient are computed here
c |
c |-- grdchk_main
c     |
c     |-- grdchk_init
c     |-- do icomp=...        - loop over control vector elements
c         |
c         |-- grdchk_loc      - determine location of icomp on grid
c         |
c         |-- grdchk_getxx    - get control vector component from file
c         |                     perturb it and write back to file
c         |-- grdchk_getadxx  - get gradient component calculated
c         |                     via adjoint
c         |-- the_main_loop   - forward run and cost evaluation
c         |                     with perturbed control vector element
c         |-- calculate ratio of adj. vs. finite difference gradient
c         |
c         |-- grdchk_setxx    - Reset control vector element
c         |
c         |-- grdchk_print    - print results
c \ev
CEOI

CBOP
C     !ROUTINE: grdchk_main
C     !INTERFACE:
      subroutine grdchk_main( mythid )

C     !DESCRIPTION: \bv
c     ==================================================================
c     SUBROUTINE grdchk_main
c     ==================================================================
c     o Compare the gradients calculated by the adjoint model with
c       finite difference approximations.
c     started: Christian Eckert eckert@mit.edu 24-Feb-2000
c     continued&finished: heimbach@mit.edu: 13-Jun-2001
c     changed: mlosch@ocean.mit.edu: 09-May-2002
c              - added centered difference vs. 1-sided difference option
c              - improved output format for readability
c              - added control variable hFacC
c              heimbach@mit.edu 24-Feb-2003
c              - added tangent linear gradient checks
c              - fixes for multiproc. gradient checks
c              - added more control variables
c
c     ==================================================================
c     SUBROUTINE grdchk_main
c     ==================================================================
C     \ev

C     !USES:
      implicit none

c     == global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "grdchk.h"
#include "cost.h"
#include "ctrl.h"
#ifdef ALLOW_TANGENTLINEAR_RUN
#include "g_cost.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
c     == routine arguments ==
      integer mythid

#ifdef ALLOW_GRDCHK
C     !LOCAL VARIABLES:
c     == local variables ==
      integer myiter
      _RL     mytime
      integer bi, itlo, ithi
      integer bj, jtlo, jthi
      integer i,  imin, imax
      integer j,  jmin, jmax
      integer k

      integer icomp
      integer ichknum
      integer icvrec
      integer jtile
      integer itile
      integer layer
      integer obcspos
      integer itilepos
      integer jtilepos
      integer icglo
      integer itest
      integer ierr
      integer ierr_grdchk
      _RL     gfd
      _RL     fcref
      _RL     fcpertplus, fcpertminus
      _RL     ratio_ad
      _RL     ratio_ftl
      _RL     xxmemo_ref
      _RL     xxmemo_pert
      _RL     adxxmemo
      _RL     ftlxxmemo
      _RL     localEps
      _RL     grdchk_epsfac

      _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)

      CHARACTER*(MAX_LEN_MBUF) msgBuf

c     == end of interface ==
CEOP

c--   Set the loop ranges.
      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)
      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

      print *, 'ph-check entering grdchk_main '

c--   initialise variables
      call grdchk_init( mythid )

c--   Compute the adjoint model gradients.
c--   Compute the unperturbed cost function.
cph   Gradient via adjoint has already been computed,
cph   and so has unperturbed cost function,
cph   assuming all xx_ fields are initialised to zero.

      ierr      = 0
      ierr_grdchk = 0
      adxxmemo  = 0.
      ftlxxmemo = 0.
#ifdef ALLOW_ADMTLM
      fcref = objf_state_final(idep,jdep,1,1,1)
#else
      fcref = fc
#endif

      print *, 'ph-check fcref = ', fcref

      do bj = jtlo, jthi
         do bi = itlo, ithi
            do k = 1, nr
               do j = jmin, jmax
                  do i = imin, imax
                     tmpplot1(i,j,k,bi,bj) = 0. _d 0
                     tmpplot2(i,j,k,bi,bj) = 0. _d 0
                     tmpplot3(i,j,k,bi,bj) = 0. _d 0
                  end do
               end do
            end do
         end do
      end do

      if ( useCentralDiff ) then
         grdchk_epsfac = 2. _d 0
      else
         grdchk_epsfac = 1. _d 0
      end if

      WRITE(standardmessageunit,'(A)')
     &    'grad-res -------------------------------'
      WRITE(standardmessageunit,'(2a)')
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '               fc ref           fc + eps           fc - eps'
#ifdef ALLOW_TANGENTLINEAR_RUN
      WRITE(standardmessageunit,'(2a)')
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '             tlm grad            fd grad         1 - fd/tlm'
#else
      WRITE(standardmessageunit,'(2a)')
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '             adj grad            fd grad         1 - fd/adj'
#endif

c--   Compute the finite difference approximations.
c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
c--   gradient checks.

      if ( nbeg .EQ. 0 )
     &     call grdchk_get_position( mythid )

      do icomp = nbeg, nend, nstep

         ichknum = (icomp - nbeg)/nstep + 1
         adxxmemo  = 0.

cph(
cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ',
cph-print     &        nbeg, icomp, ichknum
cph)
         if (ichknum .le. maxgrdchecks ) then

c--         Determine the location of icomp on the grid.
            if ( myProcId .EQ. grdchkwhichproc ) then
               call grdchk_loc( icomp, ichknum,
     &              icvrec, itile, jtile, layer, obcspos,
     &              itilepos, jtilepos, icglo, itest, ierr,
     &              mythid )
cph(
cph-print               print *, 'ph-grd ----- back from loc -----',
cph-print     &             icvrec, itilepos, jtilepos, layer, obcspos
cph)
            else
              itile    = 0
              jtile    = 0
              layer    = 0
              itilepos = 0
              jtilepos = 0
            endif
            _BARRIER

c******************************************************
c--   (A): get gradient component calculated via adjoint
c******************************************************

c--   get gradient component calculated via adjoint
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_getadxx( icvrec,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              adxxmemo, mythid )
            endif
C--   Add a global-sum call so that all proc will get the adjoint gradient
            _GLOBAL_SUM_RL( adxxmemo, myThid )
c           _BARRIER

#ifdef ALLOW_TANGENTLINEAR_RUN
c******************************************************
c--   (B): Get gradient component g_fc from tangent linear run:
c******************************************************
c--
c--   1. perturb control vector component: xx(i)=1.

            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               localEps = 1. _d 0
               call grdchk_getxx( icvrec, TANGENT_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, xxmemo_pert, localEps,
     &              mythid )
            else
              xxmemo_ref  = 0.
              xxmemo_pert = 0.
            endif
            _BARRIER

c--
c--   2. perform tangent linear run
            mytime = starttime
            myiter = niter0
#ifdef ALLOW_ADMTLM
            do k=1,4*Nr+1
             do j=1,sny
              do i=1,snx
               g_objf_state_final(i,j,1,1,k) = 0.
              enddo
             enddo
            enddo
#else
            g_fc = 0.
#endif

            call g_the_main_loop( mytime, myiter, mythid )
            _BARRIER
#ifdef ALLOW_ADMTLM
            ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
#else
            ftlxxmemo = g_fc
#endif

c--
c--   3. reset control vector
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_setxx( icvrec, TANGENT_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, mythid )
            endif
            _BARRIER

#endif /* ALLOW_TANGENTLINEAR_RUN */


c******************************************************
c--   (C): Get gradient via finite difference perturbation
c******************************************************

c--   get control vector component from file
c--   perturb it and write back to file
c--   positive perturbation
            localEps = abs(grdchk_eps)
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_getxx( icvrec, FORWARD_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, xxmemo_pert, localEps,
     &              mythid )
            else
              xxmemo_ref  = 0.
              xxmemo_pert = 0.
            endif
            _BARRIER

c--   forward run with perturbed control vector
            mytime = starttime
            myiter = niter0
            call the_main_loop( mytime, myiter, mythid )
#ifdef ALLOW_ADMTLM
            fcpertplus = objf_state_final(idep,jdep,1,1,1)
#else
            fcpertplus = fc
#endif
            print *, 'ph-check fcpertplus = ', fcpertplus
            _BARRIER

c--   Reset control vector.
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_setxx( icvrec, FORWARD_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, mythid )
            endif
            _BARRIER

            fcpertminus = fcref
            print *, 'ph-check fcpertminus = ', fcpertminus

            if ( useCentralDiff ) then

c--   get control vector component from file
c--   perturb it and write back to file
c--   repeat the proceedure for a negative perturbation
               if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
                  localEps = - abs(grdchk_eps)
                  call grdchk_getxx( icvrec, FORWARD_SIMULATION,
     &                 itile, jtile, layer,
     &                 itilepos, jtilepos,
     &                 xxmemo_ref, xxmemo_pert, localEps,
     &                 mythid )
               else
                 xxmemo_ref  = 0.
                 xxmemo_pert = 0.
               endif
               _BARRIER

c--   forward run with perturbed control vector
               mytime = starttime
               myiter = niter0
               call the_main_loop( mytime, myiter, mythid )
               _BARRIER
#ifdef ALLOW_ADMTLM
               fcpertminus = objf_state_final(idep,jdep,1,1,1)
#else
               fcpertminus = fc
#endif

c--   Reset control vector.
               if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
                  call grdchk_setxx( icvrec, FORWARD_SIMULATION,
     &                 itile, jtile, layer,
     &                 itilepos, jtilepos,
     &                 xxmemo_ref, mythid )
               endif
               _BARRIER

c-- end of if useCentralDiff ...
            end if

c******************************************************
c--   (D): calculate relative differences between gradients
c******************************************************

            if ( grdchk_eps .eq. 0. ) then
               gfd = (fcpertplus-fcpertminus)
            else
               gfd = (fcpertplus-fcpertminus)
     &              /(grdchk_epsfac*grdchk_eps)
            endif

            if ( adxxmemo .eq. 0. ) then
               ratio_ad = abs( adxxmemo - gfd )
            else
               ratio_ad = 1. - gfd/adxxmemo
            endif

            if ( ftlxxmemo .eq. 0. ) then
               ratio_ftl = abs( ftlxxmemo - gfd )
            else
               ratio_ftl = 1. - gfd/ftlxxmemo
            endif

            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               tmpplot1(itilepos,jtilepos,layer,itile,jtile)
     &              = gfd
               tmpplot2(itilepos,jtilepos,layer,itile,jtile)
     &              = ratio_ad
               tmpplot3(itilepos,jtilepos,layer,itile,jtile)
     &              = ratio_ftl
            endif

            if ( ierr .EQ. 0 ) then
               fcrmem      ( ichknum ) = fcref
               fcppmem     ( ichknum ) = fcpertplus
               fcpmmem     ( ichknum ) = fcpertminus
               xxmemref    ( ichknum ) = xxmemo_ref
               xxmempert   ( ichknum ) = xxmemo_pert
               gfdmem      ( ichknum ) = gfd
               adxxmem     ( ichknum ) = adxxmemo
               ftlxxmem    ( ichknum ) = ftlxxmemo
               ratioadmem  ( ichknum ) = ratio_ad
               ratioftlmem ( ichknum ) = ratio_ftl

               irecmem   ( ichknum ) = icvrec
               bimem     ( ichknum ) = itile
               bjmem     ( ichknum ) = jtile
               ilocmem   ( ichknum ) = itilepos
               jlocmem   ( ichknum ) = jtilepos
               klocmem   ( ichknum ) = layer
               iobcsmem  ( ichknum ) = obcspos
               icompmem  ( ichknum ) = icomp
               ichkmem   ( ichknum ) = ichknum
               itestmem  ( ichknum ) = itest
               ierrmem   ( ichknum ) = ierr
               icglomem  ( ichknum ) = icglo
            endif

            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then

               WRITE(standardmessageunit,'(A)')
     &             'grad-res -------------------------------'
               WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
     &              ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
     &              layer,itile,jtile,obcspos,
     &              fcref, fcpertplus, fcpertminus
#ifdef ALLOW_TANGENTLINEAR_RUN
               WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
     &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),
     &              bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
     &              ftlxxmemo, gfd, ratio_ftl
#else
               WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
     &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),
     &              bimem(ichknum),bjmem(ichknum),obcspos,
     &              adxxmemo, gfd, ratio_ad
#endif
            endif
#ifdef ALLOW_TANGENTLINEAR_RUN
            WRITE(msgBuf,'(A34,1PE24.14)')
     &              ' TLM  precision_derivative_cost =', fcref
            CALL PRINT_MESSAGE
     &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
            WRITE(msgBuf,'(A34,1PE24.14)')
     &              ' TLM  precision_derivative_grad =', ftlxxmemo
            CALL PRINT_MESSAGE
     &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
#else
            WRITE(msgBuf,'(A30,1PE22.14)')
     &              ' ADM  ref_cost_function      =', fcref
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
            WRITE(msgBuf,'(A30,1PE22.14)')
     &              ' ADM  adjoint_gradient       =', adxxmemo
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
            WRITE(msgBuf,'(A30,1PE22.14)')
     &              ' ADM  finite-diff_grad       =', gfd
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
#endif

            print *, 'ph-grd  ierr ---------------------------'
            print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
     &           ', ichknum = ', ichknum

            _BARRIER

c-- else of if ( ichknum ...
         else
            ierr_grdchk = -1

c-- end of if ( ichknum ...
         endif

c-- end of do icomp = ...
      enddo

      if ( myProcId .EQ. grdchkwhichproc ) then
         CALL WRITE_REC_XYZ_RL(
     &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
         CALL WRITE_REC_XYZ_RL(
     &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
         CALL WRITE_REC_XYZ_RL(
     &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
      endif

c--   Everyone has to wait for the component to be reset.
      _BARRIER

c--   Print the results of the gradient check.
      call grdchk_print( ichknum, ierr_grdchk, mythid )

#endif /* ALLOW_GRDCHK */

      return
      end
1	jmc	1.34	C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.33 2011/09/26 12:56:52 jmc Exp $
2	heimbach	1.11	C $Name: $
3	heimbach	1.2
4	jmc	1.31	#include "GRDCHK_OPTIONS.h"
5	heimbach	1.2
6	heimbach	1.3	CBOI
7			C
8			C !TITLE: GRADIENT CHECK
9			C !AUTHORS: mitgcm developers ( support@mitgcm.org )
10			C !AFFILIATION: Massachussetts Institute of Technology
11			C !DATE:
12			C !INTRODUCTION: gradient check package
13			c \bv
14			c Compare the gradients calculated by the adjoint model with
15			c finite difference approximations.
16			c
17			C !CALLING SEQUENCE:
18			c
19			c the_model_main
20			c \|
21			c \|-- ctrl_unpack
22			c \|-- adthe_main_loop - unperturbed cost function and
23			c \|-- ctrl_pack adjoint gradient are computed here
24			c \|
25			c \|-- grdchk_main
26			c \|
27			c \|-- grdchk_init
28			c \|-- do icomp=... - loop over control vector elements
29			c \|
30			c \|-- grdchk_loc - determine location of icomp on grid
31			c \|
32			c \|-- grdchk_getxx - get control vector component from file
33			c \| perturb it and write back to file
34	jmc	1.31	c \|-- grdchk_getadxx - get gradient component calculated
35	heimbach	1.3	c \| via adjoint
36			c \|-- the_main_loop - forward run and cost evaluation
37			c \| with perturbed control vector element
38			c \|-- calculate ratio of adj. vs. finite difference gradient
39			c \|
40			c \|-- grdchk_setxx - Reset control vector element
41			c \|
42			c \|-- grdchk_print - print results
43			c \ev
44			CEOI
45
46			CBOP
47			C !ROUTINE: grdchk_main
48			C !INTERFACE:
49			subroutine grdchk_main( mythid )
50	heimbach	1.2
51	heimbach	1.3	C !DESCRIPTION: \bv
52	heimbach	1.2	c ==================================================================
53			c SUBROUTINE grdchk_main
54			c ==================================================================
55			c o Compare the gradients calculated by the adjoint model with
56			c finite difference approximations.
57			c started: Christian Eckert eckert@mit.edu 24-Feb-2000
58	heimbach	1.4	c continued&finished: heimbach@mit.edu: 13-Jun-2001
59			c changed: mlosch@ocean.mit.edu: 09-May-2002
60			c - added centered difference vs. 1-sided difference option
61			c - improved output format for readability
62			c - added control variable hFacC
63	heimbach	1.7	c heimbach@mit.edu 24-Feb-2003
64			c - added tangent linear gradient checks
65			c - fixes for multiproc. gradient checks
66			c - added more control variables
67	jmc	1.31	c
68	heimbach	1.2	c ==================================================================
69			c SUBROUTINE grdchk_main
70			c ==================================================================
71	heimbach	1.3	C \ev
72	heimbach	1.2
73	heimbach	1.3	C !USES:
74	heimbach	1.2	implicit none
75
76			c == global variables ==
77			#include "SIZE.h"
78			#include "EEPARAMS.h"
79			#include "PARAMS.h"
80			#include "grdchk.h"
81			#include "cost.h"
82	heimbach	1.21	#include "ctrl.h"
83	heimbach	1.9	#ifdef ALLOW_TANGENTLINEAR_RUN
84			#include "g_cost.h"
85			#endif
86	heimbach	1.2
87	heimbach	1.3	C !INPUT/OUTPUT PARAMETERS:
88	heimbach	1.2	c == routine arguments ==
89	jmc	1.31	integer mythid
90	heimbach	1.2
91	heimbach	1.11	#ifdef ALLOW_GRDCHK
92	heimbach	1.3	C !LOCAL VARIABLES:
93	heimbach	1.2	c == local variables ==
94			integer myiter
95	jmc	1.31	_RL mytime
96	heimbach	1.2	integer bi, itlo, ithi
97			integer bj, jtlo, jthi
98			integer i, imin, imax
99			integer j, jmin, jmax
100			integer k
101
102			integer icomp
103			integer ichknum
104			integer icvrec
105			integer jtile
106			integer itile
107			integer layer
108	heimbach	1.8	integer obcspos
109	heimbach	1.2	integer itilepos
110			integer jtilepos
111	heimbach	1.19	integer icglo
112	heimbach	1.2	integer itest
113			integer ierr
114			integer ierr_grdchk
115			_RL gfd
116			_RL fcref
117	heimbach	1.4	_RL fcpertplus, fcpertminus
118	heimbach	1.6	_RL ratio_ad
119			_RL ratio_ftl
120	heimbach	1.2	_RL xxmemo_ref
121			_RL xxmemo_pert
122			_RL adxxmemo
123	heimbach	1.6	_RL ftlxxmemo
124			_RL localEps
125			_RL grdchk_epsfac
126
127	heimbach	1.7	_RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
128			_RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
129			_RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
130
131	heimbach	1.12	CHARACTER*(MAX_LEN_MBUF) msgBuf
132
133	heimbach	1.2	c == end of interface ==
134	heimbach	1.3	CEOP
135	heimbach	1.2
136			c-- Set the loop ranges.
137			jtlo = mybylo(mythid)
138			jthi = mybyhi(mythid)
139			itlo = mybxlo(mythid)
140			ithi = mybxhi(mythid)
141			jmin = 1
142			jmax = sny
143			imin = 1
144			imax = snx
145
146	heimbach	1.16	print *, 'ph-check entering grdchk_main '
147
148	heimbach	1.2	c-- initialise variables
149			call grdchk_init( mythid )
150
151	jmc	1.29	c-- Compute the adjoint model gradients.
152	heimbach	1.2	c-- Compute the unperturbed cost function.
153	heimbach	1.7	cph Gradient via adjoint has already been computed,
154			cph and so has unperturbed cost function,
155			cph assuming all xx_ fields are initialised to zero.
156	heimbach	1.2
157	heimbach	1.28	ierr = 0
158	heimbach	1.9	ierr_grdchk = 0
159	heimbach	1.27	adxxmemo = 0.
160			ftlxxmemo = 0.
161	heimbach	1.21	#ifdef ALLOW_ADMTLM
162			fcref = objf_state_final(idep,jdep,1,1,1)
163			#else
164	heimbach	1.2	fcref = fc
165	heimbach	1.21	#endif
166	heimbach	1.9
167			print *, 'ph-check fcref = ', fcref
168	heimbach	1.2
169			do bj = jtlo, jthi
170			do bi = itlo, ithi
171			do k = 1, nr
172			do j = jmin, jmax
173			do i = imin, imax
174			tmpplot1(i,j,k,bi,bj) = 0. _d 0
175			tmpplot2(i,j,k,bi,bj) = 0. _d 0
176	heimbach	1.6	tmpplot3(i,j,k,bi,bj) = 0. _d 0
177	heimbach	1.2	end do
178			end do
179			end do
180			end do
181			end do
182
183	heimbach	1.4	if ( useCentralDiff ) then
184			grdchk_epsfac = 2. _d 0
185			else
186			grdchk_epsfac = 1. _d 0
187			end if
188
189	heimbach	1.23	WRITE(standardmessageunit,'(A)')
190			& 'grad-res -------------------------------'
191			WRITE(standardmessageunit,'(2a)')
192	heimbach	1.19	& ' grad-res proc # i j k bi bj iobc',
193	heimbach	1.18	& ' fc ref fc + eps fc - eps'
194	heimbach	1.7	#ifdef ALLOW_TANGENTLINEAR_RUN
195	heimbach	1.23	WRITE(standardmessageunit,'(2a)')
196	heimbach	1.19	& ' grad-res proc # i j k bi bj iobc',
197	heimbach	1.18	& ' tlm grad fd grad 1 - fd/tlm'
198	heimbach	1.7	#else
199	heimbach	1.23	WRITE(standardmessageunit,'(2a)')
200	heimbach	1.19	& ' grad-res proc # i j k bi bj iobc',
201	heimbach	1.18	& ' adj grad fd grad 1 - fd/adj'
202	heimbach	1.7	#endif
203
204	heimbach	1.2	c-- Compute the finite difference approximations.
205			c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
206			c-- gradient checks.
207	heimbach	1.14
208	jmc	1.31	if ( nbeg .EQ. 0 )
209	heimbach	1.17	& call grdchk_get_position( mythid )
210	heimbach	1.2
211	heimbach	1.7	do icomp = nbeg, nend, nstep
212	heimbach	1.2
213	heimbach	1.7	ichknum = (icomp - nbeg)/nstep + 1
214	jmc	1.32	adxxmemo = 0.
215	heimbach	1.2
216	heimbach	1.19	cph(
217	jmc	1.31	cph-print print *, 'ph-grd _main: nbeg, icomp, ichknum ',
218	heimbach	1.19	cph-print & nbeg, icomp, ichknum
219			cph)
220	heimbach	1.7	if (ichknum .le. maxgrdchecks ) then
221	heimbach	1.2
222	heimbach	1.7	c-- Determine the location of icomp on the grid.
223			if ( myProcId .EQ. grdchkwhichproc ) then
224	jmc	1.31	call grdchk_loc( icomp, ichknum,
225	heimbach	1.8	& icvrec, itile, jtile, layer, obcspos,
226	heimbach	1.19	& itilepos, jtilepos, icglo, itest, ierr,
227	heimbach	1.7	& mythid )
228	heimbach	1.18	cph(
229	heimbach	1.19	cph-print print *, 'ph-grd ----- back from loc -----',
230			cph-print & icvrec, itilepos, jtilepos, layer, obcspos
231	heimbach	1.18	cph)
232	jmc	1.34	else
233			itile = 0
234			jtile = 0
235			layer = 0
236			itilepos = 0
237			jtilepos = 0
238	heimbach	1.7	endif
239			_BARRIER
240	jmc	1.31
241	heimbach	1.6	c******************************************************
242			c-- (A): get gradient component calculated via adjoint
243			c******************************************************
244	heimbach	1.7
245			c-- get gradient component calculated via adjoint
246			if ( myProcId .EQ. grdchkwhichproc .AND.
247			& ierr .EQ. 0 ) then
248			call grdchk_getadxx( icvrec,
249			& itile, jtile, layer,
250			& itilepos, jtilepos,
251			& adxxmemo, mythid )
252			endif
253	jmc	1.32	C-- Add a global-sum call so that all proc will get the adjoint gradient
254			_GLOBAL_SUM_RL( adxxmemo, myThid )
255			c _BARRIER
256	heimbach	1.4
257	heimbach	1.6	#ifdef ALLOW_TANGENTLINEAR_RUN
258			c******************************************************
259			c-- (B): Get gradient component g_fc from tangent linear run:
260			c******************************************************
261			c--
262			c-- 1. perturb control vector component: xx(i)=1.
263
264	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
265			& ierr .EQ. 0 ) then
266			localEps = 1. _d 0
267			call grdchk_getxx( icvrec, TANGENT_SIMULATION,
268			& itile, jtile, layer,
269			& itilepos, jtilepos,
270			& xxmemo_ref, xxmemo_pert, localEps,
271			& mythid )
272	jmc	1.34	else
273			xxmemo_ref = 0.
274			xxmemo_pert = 0.
275	heimbach	1.7	endif
276			_BARRIER
277	heimbach	1.6
278			c--
279			c-- 2. perform tangent linear run
280	heimbach	1.7	mytime = starttime
281			myiter = niter0
282	heimbach	1.21	#ifdef ALLOW_ADMTLM
283			do k=1,4*Nr+1
284			do j=1,sny
285			do i=1,snx
286			g_objf_state_final(i,j,1,1,k) = 0.
287			enddo
288			enddo
289			enddo
290			#else
291	heimbach	1.7	g_fc = 0.
292	heimbach	1.21	#endif
293
294	heimbach	1.7	call g_the_main_loop( mytime, myiter, mythid )
295	heimbach	1.21	_BARRIER
296			#ifdef ALLOW_ADMTLM
297			ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
298			#else
299	heimbach	1.7	ftlxxmemo = g_fc
300	heimbach	1.21	#endif
301
302	heimbach	1.6	c--
303			c-- 3. reset control vector
304	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
305			& ierr .EQ. 0 ) then
306			call grdchk_setxx( icvrec, TANGENT_SIMULATION,
307			& itile, jtile, layer,
308			& itilepos, jtilepos,
309			& xxmemo_ref, mythid )
310			endif
311			_BARRIER
312
313			#endif /* ALLOW_TANGENTLINEAR_RUN */
314
315	heimbach	1.6
316			c******************************************************
317			c-- (C): Get gradient via finite difference perturbation
318			c******************************************************
319
320	heimbach	1.2	c-- get control vector component from file
321	heimbach	1.7	c-- perturb it and write back to file
322	heimbach	1.6	c-- positive perturbation
323	heimbach	1.7	localEps = abs(grdchk_eps)
324			if ( myProcId .EQ. grdchkwhichproc .AND.
325			& ierr .EQ. 0 ) then
326			call grdchk_getxx( icvrec, FORWARD_SIMULATION,
327			& itile, jtile, layer,
328			& itilepos, jtilepos,
329			& xxmemo_ref, xxmemo_pert, localEps,
330			& mythid )
331	jmc	1.34	else
332			xxmemo_ref = 0.
333			xxmemo_pert = 0.
334	heimbach	1.7	endif
335			_BARRIER
336	jmc	1.31
337	heimbach	1.4	c-- forward run with perturbed control vector
338	heimbach	1.7	mytime = starttime
339			myiter = niter0
340			call the_main_loop( mytime, myiter, mythid )
341	heimbach	1.21	#ifdef ALLOW_ADMTLM
342			fcpertplus = objf_state_final(idep,jdep,1,1,1)
343			#else
344	heimbach	1.7	fcpertplus = fc
345	heimbach	1.21	#endif
346	heimbach	1.13	print *, 'ph-check fcpertplus = ', fcpertplus
347	heimbach	1.7	_BARRIER
348	jmc	1.31
349	heimbach	1.4	c-- Reset control vector.
350	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
351			& ierr .EQ. 0 ) then
352			call grdchk_setxx( icvrec, FORWARD_SIMULATION,
353			& itile, jtile, layer,
354			& itilepos, jtilepos,
355			& xxmemo_ref, mythid )
356			endif
357			_BARRIER
358	heimbach	1.2
359	heimbach	1.7	fcpertminus = fcref
360	heimbach	1.13	print *, 'ph-check fcpertminus = ', fcpertminus
361	heimbach	1.4
362	heimbach	1.7	if ( useCentralDiff ) then
363	heimbach	1.4
364	heimbach	1.6	c-- get control vector component from file
365	heimbach	1.7	c-- perturb it and write back to file
366	heimbach	1.4	c-- repeat the proceedure for a negative perturbation
367	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
368			& ierr .EQ. 0 ) then
369			localEps = - abs(grdchk_eps)
370			call grdchk_getxx( icvrec, FORWARD_SIMULATION,
371			& itile, jtile, layer,
372			& itilepos, jtilepos,
373			& xxmemo_ref, xxmemo_pert, localEps,
374			& mythid )
375	jmc	1.34	else
376			xxmemo_ref = 0.
377			xxmemo_pert = 0.
378	heimbach	1.7	endif
379			_BARRIER
380	jmc	1.31
381	heimbach	1.2	c-- forward run with perturbed control vector
382	heimbach	1.7	mytime = starttime
383			myiter = niter0
384			call the_main_loop( mytime, myiter, mythid )
385			_BARRIER
386	heimbach	1.21	#ifdef ALLOW_ADMTLM
387			fcpertminus = objf_state_final(idep,jdep,1,1,1)
388			#else
389	heimbach	1.7	fcpertminus = fc
390	heimbach	1.21	#endif
391	jmc	1.31
392	heimbach	1.4	c-- Reset control vector.
393	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
394			& ierr .EQ. 0 ) then
395			call grdchk_setxx( icvrec, FORWARD_SIMULATION,
396			& itile, jtile, layer,
397			& itilepos, jtilepos,
398			& xxmemo_ref, mythid )
399			endif
400			_BARRIER
401
402			c-- end of if useCentralDiff ...
403			end if
404	heimbach	1.4
405	heimbach	1.6	c******************************************************
406			c-- (D): calculate relative differences between gradients
407			c******************************************************
408
409	jmc	1.33	if ( grdchk_eps .eq. 0. ) then
410			gfd = (fcpertplus-fcpertminus)
411			else
412			gfd = (fcpertplus-fcpertminus)
413			& /(grdchk_epsfac*grdchk_eps)
414			endif
415	heimbach	1.2
416	jmc	1.33	if ( adxxmemo .eq. 0. ) then
417			ratio_ad = abs( adxxmemo - gfd )
418			else
419			ratio_ad = 1. - gfd/adxxmemo
420			endif
421	jmc	1.31
422	jmc	1.33	if ( ftlxxmemo .eq. 0. ) then
423			ratio_ftl = abs( ftlxxmemo - gfd )
424			else
425			ratio_ftl = 1. - gfd/ftlxxmemo
426			endif
427	heimbach	1.7
428	jmc	1.33	if ( myProcId .EQ. grdchkwhichproc .AND.
429			& ierr .EQ. 0 ) then
430	heimbach	1.7	tmpplot1(itilepos,jtilepos,layer,itile,jtile)
431			& = gfd
432			tmpplot2(itilepos,jtilepos,layer,itile,jtile)
433			& = ratio_ad
434			tmpplot3(itilepos,jtilepos,layer,itile,jtile)
435			& = ratio_ftl
436	jmc	1.34	endif
437	heimbach	1.7
438	jmc	1.34	if ( ierr .EQ. 0 ) then
439	heimbach	1.7	fcrmem ( ichknum ) = fcref
440			fcppmem ( ichknum ) = fcpertplus
441			fcpmmem ( ichknum ) = fcpertminus
442			xxmemref ( ichknum ) = xxmemo_ref
443			xxmempert ( ichknum ) = xxmemo_pert
444			gfdmem ( ichknum ) = gfd
445			adxxmem ( ichknum ) = adxxmemo
446			ftlxxmem ( ichknum ) = ftlxxmemo
447			ratioadmem ( ichknum ) = ratio_ad
448			ratioftlmem ( ichknum ) = ratio_ftl
449
450			irecmem ( ichknum ) = icvrec
451			bimem ( ichknum ) = itile
452			bjmem ( ichknum ) = jtile
453			ilocmem ( ichknum ) = itilepos
454			jlocmem ( ichknum ) = jtilepos
455			klocmem ( ichknum ) = layer
456	heimbach	1.8	iobcsmem ( ichknum ) = obcspos
457	heimbach	1.7	icompmem ( ichknum ) = icomp
458			ichkmem ( ichknum ) = ichknum
459			itestmem ( ichknum ) = itest
460			ierrmem ( ichknum ) = ierr
461	heimbach	1.19	icglomem ( ichknum ) = icglo
462	jmc	1.34	endif
463
464			if ( myProcId .EQ. grdchkwhichproc .AND.
465			& ierr .EQ. 0 ) then
466	heimbach	1.19
467	heimbach	1.23	WRITE(standardmessageunit,'(A)')
468			& 'grad-res -------------------------------'
469	jmc	1.32	WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
470	heimbach	1.23	& ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
471	heimbach	1.19	& layer,itile,jtile,obcspos,
472	heimbach	1.7	& fcref, fcpertplus, fcpertminus
473			#ifdef ALLOW_TANGENTLINEAR_RUN
474	jmc	1.32	WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
475	heimbach	1.23	& ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
476	heimbach	1.19	& icompmem(ichknum),itestmem(ichknum),
477			& bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
478	heimbach	1.7	& ftlxxmemo, gfd, ratio_ftl
479			#else
480	jmc	1.32	WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
481	heimbach	1.23	& ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
482	heimbach	1.19	& icompmem(ichknum),itestmem(ichknum),
483			& bimem(ichknum),bjmem(ichknum),obcspos,
484	heimbach	1.7	& adxxmemo, gfd, ratio_ad
485	jmc	1.32	#endif
486			endif
487			#ifdef ALLOW_TANGENTLINEAR_RUN
488			WRITE(msgBuf,'(A34,1PE24.14)')
489			& ' TLM precision_derivative_cost =', fcref
490			CALL PRINT_MESSAGE
491			& (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
492			WRITE(msgBuf,'(A34,1PE24.14)')
493			& ' TLM precision_derivative_grad =', ftlxxmemo
494			CALL PRINT_MESSAGE
495			& (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,myThid)
496			#else
497	jmc	1.33	WRITE(msgBuf,'(A30,1PE22.14)')
498			& ' ADM ref_cost_function =', fcref
499			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
500			& SQUEEZE_RIGHT, myThid )
501			WRITE(msgBuf,'(A30,1PE22.14)')
502			& ' ADM adjoint_gradient =', adxxmemo
503			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
504			& SQUEEZE_RIGHT, myThid )
505			WRITE(msgBuf,'(A30,1PE22.14)')
506			& ' ADM finite-diff_grad =', gfd
507			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
508			& SQUEEZE_RIGHT, myThid )
509	heimbach	1.7	#endif
510
511			print *, 'ph-grd ierr ---------------------------'
512			print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
513			& ', ichknum = ', ichknum
514
515			_BARRIER
516
517			c-- else of if ( ichknum ...
518			else
519			ierr_grdchk = -1
520
521	jmc	1.31	c-- end of if ( ichknum ...
522	heimbach	1.2	endif
523
524	jmc	1.31	c-- end of do icomp = ...
525	heimbach	1.7	enddo
526	heimbach	1.2
527	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc ) then
528	jmc	1.31	CALL WRITE_REC_XYZ_RL(
529	heimbach	1.7	& 'grd_findiff' , tmpplot1, 1, 0, myThid)
530	jmc	1.31	CALL WRITE_REC_XYZ_RL(
531	heimbach	1.7	& 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
532	jmc	1.31	CALL WRITE_REC_XYZ_RL(
533	heimbach	1.7	& 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
534			endif
535	heimbach	1.2
536	heimbach	1.7	c-- Everyone has to wait for the component to be reset.
537			_BARRIER
538	heimbach	1.2
539			c-- Print the results of the gradient check.
540			call grdchk_print( ichknum, ierr_grdchk, mythid )
541
542	heimbach	1.11	#endif /* ALLOW_GRDCHK */
543	heimbach	1.2
544	jmc	1.31	return
545	heimbach	1.2	end