pkg/grdchk/grdchk_main.F

C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.35 2012/02/09 17:55:54 jmc Exp $
C $Name:  $

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