pkg/grdchk/grdchk_main.F

C
C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.17 2006/05/12 02:17:03 heimbach Exp $
C $Name:  $

#include "AD_CONFIG.h"
#include "CPP_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"
#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 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 models' 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_grdchk = 0
cphadmtlm(
      fcref = fc
cphadmtlm      fcref = objf_state_final(45,4,1,1,1)
cphadmtlm)

      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

      print *, 'grad-res -------------------------------'
      print ('(2a)'), 
     &    ' grad-res  proc    #    i    j    k iobc',
     &    '               fc ref           fc + eps           fc - eps'
#ifdef ALLOW_TANGENTLINEAR_RUN
      print ('(2a)'), 
     &    ' grad-res  proc    #    i    j    k iobc',
     &    '             tlm grad            fd grad         1 - fd/tlm'
#else
      print ('(2a)'), 
     &    ' grad-res  proc    #    i    j    k 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

         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, itest, ierr,
     &              mythid )
cph(
               print *, 'ph-grd ----- back from loc -----',
     &             icvrec, itilepos, jtilepos, layer, obcspos
cph)
            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
            _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 )
            endif
            _BARRIER

c--
c--   2. perform tangent linear run
            mytime = starttime
            myiter = niter0
cphadmtlm(
            g_fc = 0.
cphadmtlm            do j=1,sny
cphadmtlm               do i=1,snx
cphadmtlm                  g_objf_state_final(i,j,1,1,1) = 0.
cphadmtlm                  g_objf_state_final(i,j,1,1,2) = 0.
cphadmtlm               enddo
cphadmtlm            enddo
cphadmtlm)
            call g_the_main_loop( mytime, myiter, mythid )
cphadmtlm(
            ftlxxmemo = g_fc
cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1,1)
cphadmtlm)
            _BARRIER
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 )
            endif
            _BARRIER
                     
c--   forward run with perturbed control vector
            mytime = starttime
            myiter = niter0
            call the_main_loop( mytime, myiter, mythid )
cphadmtlm(
            fcpertplus = fc
cphadmtlm            fcpertplus = objf_state_final(45,4,1,1,1)
cphadmtlm)
            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 )
               endif
               _BARRIER
                     
c--   forward run with perturbed control vector
               mytime = starttime
               myiter = niter0
               call the_main_loop( mytime, myiter, mythid )
               _BARRIER
               fcpertminus = fc
                  
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 ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then

               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
                  
               tmpplot1(itilepos,jtilepos,layer,itile,jtile)
     &              = gfd
               tmpplot2(itilepos,jtilepos,layer,itile,jtile)
     &              = ratio_ad
               tmpplot3(itilepos,jtilepos,layer,itile,jtile)
     &              = ratio_ftl

               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
                  
               print *, 'grad-res -------------------------------'
               print '(a,6I5,2x,3(1x,E18.12))', ' grad-res ',
     &              myprocid,ichknum,itilepos,jtilepos,layer,obcspos,
     &              fcref, fcpertplus, fcpertminus
#ifdef ALLOW_TANGENTLINEAR_RUN
               print '(a,6I5,2x,3(1x,E18.12))', ' grad-res ',
     &              myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),obcspos,
     &              ftlxxmemo, gfd, ratio_ftl
               WRITE(msgBuf,'(A34,2(1PE24.14,X))') 
     &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
               CALL PRINT_MESSAGE
     &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
#else
               print '(a,6I5,2x,3(1x,E18.12))', ' grad-res ',
     &              myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),obcspos,
     &              adxxmemo, gfd, ratio_ad
               WRITE(msgBuf,'(A34,2(1PE24.14,X))') 
     &              'precision_grdchk_result ADM ', fcref, adxxmemo
               CALL PRINT_MESSAGE
     &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
#endif

            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 */

      end
1	edhill	1.10	C
2	heimbach	1.18	C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.17 2006/05/12 02:17:03 heimbach Exp $
3	heimbach	1.11	C $Name: $
4	heimbach	1.2
5	edhill	1.10	#include "AD_CONFIG.h"
6	heimbach	1.2	#include "CPP_OPTIONS.h"
7
8	heimbach	1.3	CBOI
9			C
10			C !TITLE: GRADIENT CHECK
11			C !AUTHORS: mitgcm developers ( support@mitgcm.org )
12			C !AFFILIATION: Massachussetts Institute of Technology
13			C !DATE:
14			C !INTRODUCTION: gradient check package
15			c \bv
16			c Compare the gradients calculated by the adjoint model with
17			c finite difference approximations.
18			c
19			C !CALLING SEQUENCE:
20			c
21			c the_model_main
22			c \|
23			c \|-- ctrl_unpack
24			c \|-- adthe_main_loop - unperturbed cost function and
25			c \|-- ctrl_pack adjoint gradient are computed here
26			c \|
27			c \|-- grdchk_main
28			c \|
29			c \|-- grdchk_init
30			c \|-- do icomp=... - loop over control vector elements
31			c \|
32			c \|-- grdchk_loc - determine location of icomp on grid
33			c \|
34			c \|-- grdchk_getxx - get control vector component from file
35			c \| perturb it and write back to file
36			c \|-- grdchk_getadxx - get gradient component calculated
37			c \| via adjoint
38			c \|-- the_main_loop - forward run and cost evaluation
39			c \| with perturbed control vector element
40			c \|-- calculate ratio of adj. vs. finite difference gradient
41			c \|
42			c \|-- grdchk_setxx - Reset control vector element
43			c \|
44			c \|-- grdchk_print - print results
45			c \ev
46			CEOI
47
48			CBOP
49			C !ROUTINE: grdchk_main
50			C !INTERFACE:
51			subroutine grdchk_main( mythid )
52	heimbach	1.2
53	heimbach	1.3	C !DESCRIPTION: \bv
54	heimbach	1.2	c ==================================================================
55			c SUBROUTINE grdchk_main
56			c ==================================================================
57			c o Compare the gradients calculated by the adjoint model with
58			c finite difference approximations.
59			c started: Christian Eckert eckert@mit.edu 24-Feb-2000
60	heimbach	1.4	c continued&finished: heimbach@mit.edu: 13-Jun-2001
61			c changed: mlosch@ocean.mit.edu: 09-May-2002
62			c - added centered difference vs. 1-sided difference option
63			c - improved output format for readability
64			c - added control variable hFacC
65	heimbach	1.7	c heimbach@mit.edu 24-Feb-2003
66			c - added tangent linear gradient checks
67			c - fixes for multiproc. gradient checks
68			c - added more control variables
69	heimbach	1.4	c
70	heimbach	1.2	c ==================================================================
71			c SUBROUTINE grdchk_main
72			c ==================================================================
73	heimbach	1.3	C \ev
74	heimbach	1.2
75	heimbach	1.3	C !USES:
76	heimbach	1.2	implicit none
77
78			c == global variables ==
79			#include "SIZE.h"
80			#include "EEPARAMS.h"
81			#include "PARAMS.h"
82			#include "grdchk.h"
83			#include "cost.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			integer mythid
91
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			_RL mytime
97			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			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			c-- Compute the adjoint models' gradients.
152			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.9	ierr_grdchk = 0
158			cphadmtlm(
159	heimbach	1.2	fcref = fc
160	heimbach	1.15	cphadmtlm fcref = objf_state_final(45,4,1,1,1)
161	heimbach	1.9	cphadmtlm)
162
163			print *, 'ph-check fcref = ', fcref
164	heimbach	1.2
165			do bj = jtlo, jthi
166			do bi = itlo, ithi
167			do k = 1, nr
168			do j = jmin, jmax
169			do i = imin, imax
170			tmpplot1(i,j,k,bi,bj) = 0. _d 0
171			tmpplot2(i,j,k,bi,bj) = 0. _d 0
172	heimbach	1.6	tmpplot3(i,j,k,bi,bj) = 0. _d 0
173	heimbach	1.2	end do
174			end do
175			end do
176			end do
177			end do
178
179	heimbach	1.4	if ( useCentralDiff ) then
180			grdchk_epsfac = 2. _d 0
181			else
182			grdchk_epsfac = 1. _d 0
183			end if
184
185	heimbach	1.12	print *, 'grad-res -------------------------------'
186	heimbach	1.7	print ('(2a)'),
187	heimbach	1.18	& ' grad-res proc # i j k iobc',
188			& ' fc ref fc + eps fc - eps'
189	heimbach	1.7	#ifdef ALLOW_TANGENTLINEAR_RUN
190			print ('(2a)'),
191	heimbach	1.18	& ' grad-res proc # i j k iobc',
192			& ' tlm grad fd grad 1 - fd/tlm'
193	heimbach	1.7	#else
194			print ('(2a)'),
195	heimbach	1.18	& ' grad-res proc # i j k iobc',
196			& ' adj grad fd grad 1 - fd/adj'
197	heimbach	1.7	#endif
198
199	heimbach	1.2	c-- Compute the finite difference approximations.
200			c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
201			c-- gradient checks.
202	heimbach	1.14
203	heimbach	1.17	if ( nbeg .EQ. 0 )
204			& call grdchk_get_position( mythid )
205	heimbach	1.2
206	heimbach	1.7	do icomp = nbeg, nend, nstep
207	heimbach	1.2
208	heimbach	1.7	ichknum = (icomp - nbeg)/nstep + 1
209	heimbach	1.2
210	heimbach	1.7	if (ichknum .le. maxgrdchecks ) then
211	heimbach	1.2
212	heimbach	1.7	c-- Determine the location of icomp on the grid.
213			if ( myProcId .EQ. grdchkwhichproc ) then
214			call grdchk_loc( icomp, ichknum,
215	heimbach	1.8	& icvrec, itile, jtile, layer, obcspos,
216	heimbach	1.7	& itilepos, jtilepos, itest, ierr,
217			& mythid )
218	heimbach	1.18	cph(
219			print *, 'ph-grd ----- back from loc -----',
220			& icvrec, itilepos, jtilepos, layer, obcspos
221			cph)
222	heimbach	1.7	endif
223			_BARRIER
224	heimbach	1.2
225	heimbach	1.6	c******************************************************
226			c-- (A): get gradient component calculated via adjoint
227			c******************************************************
228	heimbach	1.7
229			c-- get gradient component calculated via adjoint
230			if ( myProcId .EQ. grdchkwhichproc .AND.
231			& ierr .EQ. 0 ) then
232			call grdchk_getadxx( icvrec,
233			& itile, jtile, layer,
234			& itilepos, jtilepos,
235			& adxxmemo, mythid )
236			endif
237			_BARRIER
238	heimbach	1.4
239	heimbach	1.6	#ifdef ALLOW_TANGENTLINEAR_RUN
240			c******************************************************
241			c-- (B): Get gradient component g_fc from tangent linear run:
242			c******************************************************
243			c--
244			c-- 1. perturb control vector component: xx(i)=1.
245
246	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
247			& ierr .EQ. 0 ) then
248			localEps = 1. _d 0
249			call grdchk_getxx( icvrec, TANGENT_SIMULATION,
250			& itile, jtile, layer,
251			& itilepos, jtilepos,
252			& xxmemo_ref, xxmemo_pert, localEps,
253			& mythid )
254			endif
255			_BARRIER
256	heimbach	1.6
257			c--
258			c-- 2. perform tangent linear run
259	heimbach	1.7	mytime = starttime
260			myiter = niter0
261	heimbach	1.9	cphadmtlm(
262	heimbach	1.7	g_fc = 0.
263	heimbach	1.9	cphadmtlm do j=1,sny
264			cphadmtlm do i=1,snx
265	heimbach	1.15	cphadmtlm g_objf_state_final(i,j,1,1,1) = 0.
266			cphadmtlm g_objf_state_final(i,j,1,1,2) = 0.
267	heimbach	1.9	cphadmtlm enddo
268			cphadmtlm enddo
269			cphadmtlm)
270	heimbach	1.7	call g_the_main_loop( mytime, myiter, mythid )
271	heimbach	1.9	cphadmtlm(
272	heimbach	1.7	ftlxxmemo = g_fc
273	heimbach	1.15	cphadmtlm ftlxxmemo = g_objf_state_final(45,4,1,1,1)
274	heimbach	1.9	cphadmtlm)
275	heimbach	1.7	_BARRIER
276	heimbach	1.6	c--
277			c-- 3. reset control vector
278	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
279			& ierr .EQ. 0 ) then
280			call grdchk_setxx( icvrec, TANGENT_SIMULATION,
281			& itile, jtile, layer,
282			& itilepos, jtilepos,
283			& xxmemo_ref, mythid )
284			endif
285			_BARRIER
286
287			#endif /* ALLOW_TANGENTLINEAR_RUN */
288
289	heimbach	1.6
290			c******************************************************
291			c-- (C): Get gradient via finite difference perturbation
292			c******************************************************
293
294	heimbach	1.2	c-- get control vector component from file
295	heimbach	1.7	c-- perturb it and write back to file
296	heimbach	1.6	c-- positive perturbation
297	heimbach	1.7	localEps = abs(grdchk_eps)
298			if ( myProcId .EQ. grdchkwhichproc .AND.
299			& ierr .EQ. 0 ) then
300			call grdchk_getxx( icvrec, FORWARD_SIMULATION,
301			& itile, jtile, layer,
302			& itilepos, jtilepos,
303			& xxmemo_ref, xxmemo_pert, localEps,
304			& mythid )
305			endif
306			_BARRIER
307	heimbach	1.2
308	heimbach	1.4	c-- forward run with perturbed control vector
309	heimbach	1.7	mytime = starttime
310			myiter = niter0
311			call the_main_loop( mytime, myiter, mythid )
312	heimbach	1.9	cphadmtlm(
313	heimbach	1.7	fcpertplus = fc
314	heimbach	1.15	cphadmtlm fcpertplus = objf_state_final(45,4,1,1,1)
315	heimbach	1.9	cphadmtlm)
316	heimbach	1.13	print *, 'ph-check fcpertplus = ', fcpertplus
317	heimbach	1.7	_BARRIER
318	heimbach	1.4
319			c-- Reset control vector.
320	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
321			& ierr .EQ. 0 ) then
322			call grdchk_setxx( icvrec, FORWARD_SIMULATION,
323			& itile, jtile, layer,
324			& itilepos, jtilepos,
325			& xxmemo_ref, mythid )
326			endif
327			_BARRIER
328	heimbach	1.2
329	heimbach	1.7	fcpertminus = fcref
330	heimbach	1.13	print *, 'ph-check fcpertminus = ', fcpertminus
331	heimbach	1.4
332	heimbach	1.7	if ( useCentralDiff ) then
333	heimbach	1.4
334	heimbach	1.6	c-- get control vector component from file
335	heimbach	1.7	c-- perturb it and write back to file
336	heimbach	1.4	c-- repeat the proceedure for a negative perturbation
337	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
338			& ierr .EQ. 0 ) then
339			localEps = - abs(grdchk_eps)
340			call grdchk_getxx( icvrec, FORWARD_SIMULATION,
341			& itile, jtile, layer,
342			& itilepos, jtilepos,
343			& xxmemo_ref, xxmemo_pert, localEps,
344			& mythid )
345			endif
346			_BARRIER
347	heimbach	1.4
348	heimbach	1.2	c-- forward run with perturbed control vector
349	heimbach	1.7	mytime = starttime
350			myiter = niter0
351			call the_main_loop( mytime, myiter, mythid )
352			_BARRIER
353			fcpertminus = fc
354	heimbach	1.2
355	heimbach	1.4	c-- Reset control vector.
356	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
357			& ierr .EQ. 0 ) then
358			call grdchk_setxx( icvrec, FORWARD_SIMULATION,
359			& itile, jtile, layer,
360			& itilepos, jtilepos,
361			& xxmemo_ref, mythid )
362			endif
363			_BARRIER
364
365			c-- end of if useCentralDiff ...
366			end if
367	heimbach	1.4
368	heimbach	1.6	c******************************************************
369			c-- (D): calculate relative differences between gradients
370			c******************************************************
371
372	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
373			& ierr .EQ. 0 ) then
374	heimbach	1.2
375	heimbach	1.7	if ( grdchk_eps .eq. 0. ) then
376			gfd = (fcpertplus-fcpertminus)
377			else
378			gfd = (fcpertplus-fcpertminus)
379			& /(grdchk_epsfac*grdchk_eps)
380			endif
381	heimbach	1.2
382	heimbach	1.7	if ( adxxmemo .eq. 0. ) then
383			ratio_ad = abs( adxxmemo - gfd )
384			else
385			ratio_ad = 1. - gfd/adxxmemo
386			endif
387
388			if ( ftlxxmemo .eq. 0. ) then
389			ratio_ftl = abs( ftlxxmemo - gfd )
390	heimbach	1.2	else
391	heimbach	1.7	ratio_ftl = 1. - gfd/ftlxxmemo
392	heimbach	1.2	endif
393	heimbach	1.7
394			tmpplot1(itilepos,jtilepos,layer,itile,jtile)
395			& = gfd
396			tmpplot2(itilepos,jtilepos,layer,itile,jtile)
397			& = ratio_ad
398			tmpplot3(itilepos,jtilepos,layer,itile,jtile)
399			& = ratio_ftl
400
401			fcrmem ( ichknum ) = fcref
402			fcppmem ( ichknum ) = fcpertplus
403			fcpmmem ( ichknum ) = fcpertminus
404			xxmemref ( ichknum ) = xxmemo_ref
405			xxmempert ( ichknum ) = xxmemo_pert
406			gfdmem ( ichknum ) = gfd
407			adxxmem ( ichknum ) = adxxmemo
408			ftlxxmem ( ichknum ) = ftlxxmemo
409			ratioadmem ( ichknum ) = ratio_ad
410			ratioftlmem ( ichknum ) = ratio_ftl
411
412			irecmem ( ichknum ) = icvrec
413			bimem ( ichknum ) = itile
414			bjmem ( ichknum ) = jtile
415			ilocmem ( ichknum ) = itilepos
416			jlocmem ( ichknum ) = jtilepos
417			klocmem ( ichknum ) = layer
418	heimbach	1.8	iobcsmem ( ichknum ) = obcspos
419	heimbach	1.7	icompmem ( ichknum ) = icomp
420			ichkmem ( ichknum ) = ichknum
421			itestmem ( ichknum ) = itest
422			ierrmem ( ichknum ) = ierr
423
424	heimbach	1.12	print *, 'grad-res -------------------------------'
425	heimbach	1.18	print '(a,6I5,2x,3(1x,E18.12))', ' grad-res ',
426			& myprocid,ichknum,itilepos,jtilepos,layer,obcspos,
427	heimbach	1.7	& fcref, fcpertplus, fcpertminus
428			#ifdef ALLOW_TANGENTLINEAR_RUN
429	heimbach	1.18	print '(a,6I5,2x,3(1x,E18.12))', ' grad-res ',
430	heimbach	1.7	& myprocid,ichknum,ichkmem(ichknum),
431	heimbach	1.18	& icompmem(ichknum),itestmem(ichknum),obcspos,
432	heimbach	1.7	& ftlxxmemo, gfd, ratio_ftl
433	heimbach	1.12	WRITE(msgBuf,'(A34,2(1PE24.14,X))')
434			& 'precision_grdchk_result TLM ', fcref, ftlxxmemo
435			CALL PRINT_MESSAGE
436			& (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
437	heimbach	1.7	#else
438	heimbach	1.18	print '(a,6I5,2x,3(1x,E18.12))', ' grad-res ',
439	heimbach	1.7	& myprocid,ichknum,ichkmem(ichknum),
440	heimbach	1.18	& icompmem(ichknum),itestmem(ichknum),obcspos,
441	heimbach	1.7	& adxxmemo, gfd, ratio_ad
442	heimbach	1.12	WRITE(msgBuf,'(A34,2(1PE24.14,X))')
443			& 'precision_grdchk_result ADM ', fcref, adxxmemo
444			CALL PRINT_MESSAGE
445			& (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
446	heimbach	1.7	#endif
447
448			endif
449
450			print *, 'ph-grd ierr ---------------------------'
451			print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
452			& ', ichknum = ', ichknum
453
454			_BARRIER
455
456			c-- else of if ( ichknum ...
457			else
458			ierr_grdchk = -1
459
460			c-- end of if ( ichknum ...
461	heimbach	1.2	endif
462
463	heimbach	1.7	c-- end of do icomp = ...
464			enddo
465	heimbach	1.2
466	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc ) then
467			CALL WRITE_REC_XYZ_RL(
468			& 'grd_findiff' , tmpplot1, 1, 0, myThid)
469			CALL WRITE_REC_XYZ_RL(
470			& 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
471			CALL WRITE_REC_XYZ_RL(
472			& 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
473			endif
474	heimbach	1.2
475	heimbach	1.7	c-- Everyone has to wait for the component to be reset.
476			_BARRIER
477	heimbach	1.2
478			c-- Print the results of the gradient check.
479			call grdchk_print( ichknum, ierr_grdchk, mythid )
480
481	heimbach	1.11	#endif /* ALLOW_GRDCHK */
482	heimbach	1.2
483			end