pkg/grdchk/grdchk_main.F

C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.3.6.1 2002/05/30 19:55:17 heimbach Exp $

#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"

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

#ifdef ALLOW_GRADIENT_CHECK
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 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)

#ifdef ALLOW_TANGENTLINEAR_RUN
      _RL     g_fc
      common /g_cost_r/ g_fc
#endif

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

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.

      fcref = fc
      ierr_grdchk = 0

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

      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,
     &              itilepos, jtilepos, itest, ierr,
     &              mythid )
            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
            g_fc = 0.
            call g_the_main_loop( mytime, myiter, mythid )
            ftlxxmemo = g_fc
            _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 )
            fcpertplus = fc
            _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

            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
               icompmem  ( ichknum ) = icomp
               ichkmem   ( ichknum ) = ichknum
               itestmem  ( ichknum ) = itest
               ierrmem   ( ichknum ) = ierr
                  
               print *, 'ph-grd 3 -------------------------------'
               print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
     &              myprocid,ichknum,itilepos,jtilepos,layer,
     &              fcref, fcpertplus, fcpertminus
#ifdef ALLOW_TANGENTLINEAR_RUN
               print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
     &              myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),
     &              ftlxxmemo, gfd, ratio_ftl
#else
               print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
     &              myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),
     &              adxxmemo, gfd, ratio_ad
#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_GRADIENT_CHECK */

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