pkg/grdchk/grdchk_main.F

C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.8 2003/06/24 16:08:45 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"
#ifdef ALLOW_TANGENTLINEAR_RUN
#include "g_cost.h"
#endif

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 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)

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.

      ierr_grdchk = 0
cphadmtlm(
      fcref = fc
cphadmtlm      fcref = objf_state_final(45,4,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 *, '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, obcspos,
     &              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
cphadmtlm(
            g_fc = 0.
cphadmtlm            do j=1,sny
cphadmtlm               do i=1,snx
cphadmtlm                  g_objf_state_final(i,j,1,1) = 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)
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)
cphadmtlm)
            _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
               iobcsmem  ( ichknum ) = obcspos
               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.9	C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.8 2003/06/24 16:08:45 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	heimbach	1.9	#ifdef ALLOW_TANGENTLINEAR_RUN
82			#include "g_cost.h"
83			#endif
84	heimbach	1.2
85	heimbach	1.3	C !INPUT/OUTPUT PARAMETERS:
86	heimbach	1.2	c == routine arguments ==
87			integer mythid
88
89			#ifdef ALLOW_GRADIENT_CHECK
90	heimbach	1.3	C !LOCAL VARIABLES:
91	heimbach	1.2	c == local variables ==
92			integer myiter
93			_RL mytime
94			integer bi, itlo, ithi
95			integer bj, jtlo, jthi
96			integer i, imin, imax
97			integer j, jmin, jmax
98			integer k
99
100			integer icomp
101			integer ichknum
102			integer icvrec
103			integer jtile
104			integer itile
105			integer layer
106	heimbach	1.8	integer obcspos
107	heimbach	1.2	integer itilepos
108			integer jtilepos
109			integer itest
110			integer ierr
111			integer ierr_grdchk
112			_RL gfd
113			_RL fcref
114	heimbach	1.4	_RL fcpertplus, fcpertminus
115	heimbach	1.6	_RL ratio_ad
116			_RL ratio_ftl
117	heimbach	1.2	_RL xxmemo_ref
118			_RL xxmemo_pert
119			_RL adxxmemo
120	heimbach	1.6	_RL ftlxxmemo
121			_RL localEps
122			_RL grdchk_epsfac
123
124	heimbach	1.7	_RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
125			_RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
126			_RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
127
128	heimbach	1.2	c == end of interface ==
129	heimbach	1.3	CEOP
130	heimbach	1.2
131			c-- Set the loop ranges.
132			jtlo = mybylo(mythid)
133			jthi = mybyhi(mythid)
134			itlo = mybxlo(mythid)
135			ithi = mybxhi(mythid)
136			jmin = 1
137			jmax = sny
138			imin = 1
139			imax = snx
140
141			c-- initialise variables
142			call grdchk_init( mythid )
143
144			c-- Compute the adjoint models' gradients.
145			c-- Compute the unperturbed cost function.
146	heimbach	1.7	cph Gradient via adjoint has already been computed,
147			cph and so has unperturbed cost function,
148			cph assuming all xx_ fields are initialised to zero.
149	heimbach	1.2
150	heimbach	1.9	ierr_grdchk = 0
151			cphadmtlm(
152	heimbach	1.2	fcref = fc
153	heimbach	1.9	cphadmtlm fcref = objf_state_final(45,4,1,1)
154			cphadmtlm)
155
156			print *, 'ph-check fcref = ', fcref
157	heimbach	1.2
158			do bj = jtlo, jthi
159			do bi = itlo, ithi
160			do k = 1, nr
161			do j = jmin, jmax
162			do i = imin, imax
163			tmpplot1(i,j,k,bi,bj) = 0. _d 0
164			tmpplot2(i,j,k,bi,bj) = 0. _d 0
165	heimbach	1.6	tmpplot3(i,j,k,bi,bj) = 0. _d 0
166	heimbach	1.2	end do
167			end do
168			end do
169			end do
170			end do
171
172	heimbach	1.4	if ( useCentralDiff ) then
173			grdchk_epsfac = 2. _d 0
174			else
175			grdchk_epsfac = 1. _d 0
176			end if
177
178	heimbach	1.7	print *, 'ph-grd 3 -------------------------------'
179			print ('(2a)'),
180			& ' ph-grd 3 proc # i j k fc ref',
181			& ' fc + eps fc - eps'
182			#ifdef ALLOW_TANGENTLINEAR_RUN
183			print ('(2a)'),
184			& ' ph-grd 3 proc # i j k tlm grad',
185			& ' fd grad 1 - fd/tlm'
186			#else
187			print ('(2a)'),
188			& ' ph-grd 3 proc # i j k adj grad',
189			& ' fd grad 1 - fd/adj'
190			#endif
191
192	heimbach	1.2	c-- Compute the finite difference approximations.
193			c-- Cycle through all processes doing NINT(nend-nbeg+1)/nstep
194			c-- gradient checks.
195
196	heimbach	1.7	do icomp = nbeg, nend, nstep
197	heimbach	1.2
198	heimbach	1.7	ichknum = (icomp - nbeg)/nstep + 1
199	heimbach	1.2
200	heimbach	1.7	if (ichknum .le. maxgrdchecks ) then
201	heimbach	1.2
202	heimbach	1.7	c-- Determine the location of icomp on the grid.
203			if ( myProcId .EQ. grdchkwhichproc ) then
204			call grdchk_loc( icomp, ichknum,
205	heimbach	1.8	& icvrec, itile, jtile, layer, obcspos,
206	heimbach	1.7	& itilepos, jtilepos, itest, ierr,
207			& mythid )
208			endif
209			_BARRIER
210	heimbach	1.2
211	heimbach	1.6	c******************************************************
212			c-- (A): get gradient component calculated via adjoint
213			c******************************************************
214	heimbach	1.7
215			c-- get gradient component calculated via adjoint
216			if ( myProcId .EQ. grdchkwhichproc .AND.
217			& ierr .EQ. 0 ) then
218			call grdchk_getadxx( icvrec,
219			& itile, jtile, layer,
220			& itilepos, jtilepos,
221			& adxxmemo, mythid )
222			endif
223			_BARRIER
224	heimbach	1.4
225	heimbach	1.6	#ifdef ALLOW_TANGENTLINEAR_RUN
226			c******************************************************
227			c-- (B): Get gradient component g_fc from tangent linear run:
228			c******************************************************
229			c--
230			c-- 1. perturb control vector component: xx(i)=1.
231
232	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
233			& ierr .EQ. 0 ) then
234			localEps = 1. _d 0
235			call grdchk_getxx( icvrec, TANGENT_SIMULATION,
236			& itile, jtile, layer,
237			& itilepos, jtilepos,
238			& xxmemo_ref, xxmemo_pert, localEps,
239			& mythid )
240			endif
241			_BARRIER
242	heimbach	1.6
243			c--
244			c-- 2. perform tangent linear run
245	heimbach	1.7	mytime = starttime
246			myiter = niter0
247	heimbach	1.9	cphadmtlm(
248	heimbach	1.7	g_fc = 0.
249	heimbach	1.9	cphadmtlm do j=1,sny
250			cphadmtlm do i=1,snx
251			cphadmtlm g_objf_state_final(i,j,1,1) = 0.
252			cphadmtlm enddo
253			cphadmtlm enddo
254			cphadmtlm)
255	heimbach	1.7	call g_the_main_loop( mytime, myiter, mythid )
256	heimbach	1.9	cphadmtlm(
257	heimbach	1.7	ftlxxmemo = g_fc
258	heimbach	1.9	cphadmtlm ftlxxmemo = g_objf_state_final(45,4,1,1)
259			cphadmtlm)
260	heimbach	1.7	_BARRIER
261	heimbach	1.6	c--
262			c-- 3. reset control vector
263	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
264			& ierr .EQ. 0 ) then
265			call grdchk_setxx( icvrec, TANGENT_SIMULATION,
266			& itile, jtile, layer,
267			& itilepos, jtilepos,
268			& xxmemo_ref, mythid )
269			endif
270			_BARRIER
271
272			#endif /* ALLOW_TANGENTLINEAR_RUN */
273
274	heimbach	1.6
275			c******************************************************
276			c-- (C): Get gradient via finite difference perturbation
277			c******************************************************
278
279	heimbach	1.2	c-- get control vector component from file
280	heimbach	1.7	c-- perturb it and write back to file
281	heimbach	1.6	c-- positive perturbation
282	heimbach	1.7	localEps = abs(grdchk_eps)
283			if ( myProcId .EQ. grdchkwhichproc .AND.
284			& ierr .EQ. 0 ) then
285			call grdchk_getxx( icvrec, FORWARD_SIMULATION,
286			& itile, jtile, layer,
287			& itilepos, jtilepos,
288			& xxmemo_ref, xxmemo_pert, localEps,
289			& mythid )
290			endif
291			_BARRIER
292	heimbach	1.2
293	heimbach	1.4	c-- forward run with perturbed control vector
294	heimbach	1.7	mytime = starttime
295			myiter = niter0
296			call the_main_loop( mytime, myiter, mythid )
297	heimbach	1.9	cphadmtlm(
298	heimbach	1.7	fcpertplus = fc
299	heimbach	1.9	cphadmtlm fcpertplus = objf_state_final(45,4,1,1)
300			cphadmtlm)
301	heimbach	1.7	_BARRIER
302	heimbach	1.4
303			c-- Reset control vector.
304	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
305			& ierr .EQ. 0 ) then
306			call grdchk_setxx( icvrec, FORWARD_SIMULATION,
307			& itile, jtile, layer,
308			& itilepos, jtilepos,
309			& xxmemo_ref, mythid )
310			endif
311			_BARRIER
312	heimbach	1.2
313	heimbach	1.7	fcpertminus = fcref
314	heimbach	1.4
315	heimbach	1.7	if ( useCentralDiff ) then
316	heimbach	1.4
317	heimbach	1.6	c-- get control vector component from file
318	heimbach	1.7	c-- perturb it and write back to file
319	heimbach	1.4	c-- repeat the proceedure for a negative perturbation
320	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
321			& ierr .EQ. 0 ) then
322			localEps = - abs(grdchk_eps)
323			call grdchk_getxx( icvrec, FORWARD_SIMULATION,
324			& itile, jtile, layer,
325			& itilepos, jtilepos,
326			& xxmemo_ref, xxmemo_pert, localEps,
327			& mythid )
328			endif
329			_BARRIER
330	heimbach	1.4
331	heimbach	1.2	c-- forward run with perturbed control vector
332	heimbach	1.7	mytime = starttime
333			myiter = niter0
334			call the_main_loop( mytime, myiter, mythid )
335			_BARRIER
336			fcpertminus = fc
337	heimbach	1.2
338	heimbach	1.4	c-- Reset control vector.
339	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
340			& ierr .EQ. 0 ) then
341			call grdchk_setxx( icvrec, FORWARD_SIMULATION,
342			& itile, jtile, layer,
343			& itilepos, jtilepos,
344			& xxmemo_ref, mythid )
345			endif
346			_BARRIER
347
348			c-- end of if useCentralDiff ...
349			end if
350	heimbach	1.4
351	heimbach	1.6	c******************************************************
352			c-- (D): calculate relative differences between gradients
353			c******************************************************
354
355	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc .AND.
356			& ierr .EQ. 0 ) then
357	heimbach	1.2
358	heimbach	1.7	if ( grdchk_eps .eq. 0. ) then
359			gfd = (fcpertplus-fcpertminus)
360			else
361			gfd = (fcpertplus-fcpertminus)
362			& /(grdchk_epsfac*grdchk_eps)
363			endif
364	heimbach	1.2
365	heimbach	1.7	if ( adxxmemo .eq. 0. ) then
366			ratio_ad = abs( adxxmemo - gfd )
367			else
368			ratio_ad = 1. - gfd/adxxmemo
369			endif
370
371			if ( ftlxxmemo .eq. 0. ) then
372			ratio_ftl = abs( ftlxxmemo - gfd )
373	heimbach	1.2	else
374	heimbach	1.7	ratio_ftl = 1. - gfd/ftlxxmemo
375	heimbach	1.2	endif
376	heimbach	1.7
377			tmpplot1(itilepos,jtilepos,layer,itile,jtile)
378			& = gfd
379			tmpplot2(itilepos,jtilepos,layer,itile,jtile)
380			& = ratio_ad
381			tmpplot3(itilepos,jtilepos,layer,itile,jtile)
382			& = ratio_ftl
383
384			fcrmem ( ichknum ) = fcref
385			fcppmem ( ichknum ) = fcpertplus
386			fcpmmem ( ichknum ) = fcpertminus
387			xxmemref ( ichknum ) = xxmemo_ref
388			xxmempert ( ichknum ) = xxmemo_pert
389			gfdmem ( ichknum ) = gfd
390			adxxmem ( ichknum ) = adxxmemo
391			ftlxxmem ( ichknum ) = ftlxxmemo
392			ratioadmem ( ichknum ) = ratio_ad
393			ratioftlmem ( ichknum ) = ratio_ftl
394
395			irecmem ( ichknum ) = icvrec
396			bimem ( ichknum ) = itile
397			bjmem ( ichknum ) = jtile
398			ilocmem ( ichknum ) = itilepos
399			jlocmem ( ichknum ) = jtilepos
400			klocmem ( ichknum ) = layer
401	heimbach	1.8	iobcsmem ( ichknum ) = obcspos
402	heimbach	1.7	icompmem ( ichknum ) = icomp
403			ichkmem ( ichknum ) = ichknum
404			itestmem ( ichknum ) = itest
405			ierrmem ( ichknum ) = ierr
406
407			print *, 'ph-grd 3 -------------------------------'
408			print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
409			& myprocid,ichknum,itilepos,jtilepos,layer,
410			& fcref, fcpertplus, fcpertminus
411			#ifdef ALLOW_TANGENTLINEAR_RUN
412			print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
413			& myprocid,ichknum,ichkmem(ichknum),
414			& icompmem(ichknum),itestmem(ichknum),
415			& ftlxxmemo, gfd, ratio_ftl
416			#else
417			print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
418			& myprocid,ichknum,ichkmem(ichknum),
419			& icompmem(ichknum),itestmem(ichknum),
420			& adxxmemo, gfd, ratio_ad
421			#endif
422
423			endif
424
425			print *, 'ph-grd ierr ---------------------------'
426			print *, 'ph-grd ierr = ', ierr, ', icomp = ', icomp,
427			& ', ichknum = ', ichknum
428
429			_BARRIER
430
431			c-- else of if ( ichknum ...
432			else
433			ierr_grdchk = -1
434
435			c-- end of if ( ichknum ...
436	heimbach	1.2	endif
437
438	heimbach	1.7	c-- end of do icomp = ...
439			enddo
440	heimbach	1.2
441	heimbach	1.7	if ( myProcId .EQ. grdchkwhichproc ) then
442			CALL WRITE_REC_XYZ_RL(
443			& 'grd_findiff' , tmpplot1, 1, 0, myThid)
444			CALL WRITE_REC_XYZ_RL(
445			& 'grd_ratio_ad' , tmpplot2, 1, 0, myThid)
446			CALL WRITE_REC_XYZ_RL(
447			& 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
448			endif
449	heimbach	1.2
450	heimbach	1.7	c-- Everyone has to wait for the component to be reset.
451			_BARRIER
452	heimbach	1.2
453			c-- Print the results of the gradient check.
454			call grdchk_print( ichknum, ierr_grdchk, mythid )
455
456			#endif /* ALLOW_GRADIENT_CHECK */
457
458			end