pkg/grdchk/grdchk_main.F

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