pkg/grdchk/grdchk_main.F

C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.5 2002/07/13 02:55:58 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     
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 jprocs
      integer proc_grdchk
      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

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

      _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.
c--   Gradient via adjoint has already been computed,
c--   and so has unperturbed cost function,
c--   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

c--   Compute the finite difference approximations.
c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
c--   gradient checks.
      do jprocs = 1,numberOfProcs
         proc_grdchk = jprocs - 1

         if ( myProcId .eq. proc_grdchk ) then

            do icomp = nbeg, nend, nstep

               ichknum = (icomp - nbeg)/nstep + 1

               if (ichknum .le. maxgrdchecks ) then

c--   Determine the location of icomp on the grid.
                  call grdchk_loc( icomp, ichknum, 
     &                 icvrec, itile, jtile, layer,
     &                 itilepos, jtilepos, itest, ierr,
     &                 mythid )
              
                  if ( ierr .eq. 0 ) then

c******************************************************
c--   (A): get gradient component calculated via adjoint
c******************************************************
                     call grdchk_getadxx( icvrec,
     &                    itile, jtile, layer,
     &                    itilepos, jtilepos,
     &                    adxxmemo, mythid )
                     _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.

                     localEps = 1.
                     call grdchk_getxx( icvrec, TANGENT_SIMULATION,
     &                    itile, jtile, layer,
     &                    itilepos, jtilepos,
     &                    xxmemo_ref, xxmemo_pert, localEps,
     &                    mythid )
                     _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
c--
c--   3. reset control vector
                     call grdchk_setxx( icvrec, TANGENT_SIMULATION,
     &                    itile, jtile, layer,
     &                    itilepos, jtilepos,
     &                    xxmemo_ref, mythid )
                     _BARRIER
#endif

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)
                     call grdchk_getxx( icvrec, FORWARD_SIMULATION,
     &                    itile, jtile, layer,
     &                    itilepos, jtilepos,
     &                    xxmemo_ref, xxmemo_pert, localEps,
     &                    mythid )
                     _BARRIER
                     
c--   forward run with perturbed control vector
                     mytime = starttime
                     myiter = niter0
                     call the_main_loop( mytime, myiter, mythid )
                     fcpertplus = fc
                  
c--   Reset control vector.
                     call grdchk_setxx( icvrec, FORWARD_SIMULATION,
     &                    itile, jtile, layer,
     &                    itilepos, jtilepos,
     &                    xxmemo_ref, mythid )
                     _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
                        localEps = - abs(grdchk_eps)
                        call grdchk_getxx( icvrec, FORWARD_SIMULATION,
     &                    itile, jtile, layer,
     &                    itilepos, jtilepos,
     &                    xxmemo_ref, xxmemo_pert, localEps,
     &                    mythid )
                        _BARRIER
                     
c--   forward run with perturbed control vector
                        mytime = starttime
                        myiter = niter0
                        call the_main_loop( mytime, myiter, mythid )
                        fcpertminus = fc
                  
c--   Reset control vector.
                        call grdchk_setxx( icvrec, FORWARD_SIMULATION,
     &                    itile, jtile, layer,
     &                    itilepos, jtilepos,
     &                    xxmemo_ref, mythid )
                        _BARRIER

                     end if
c--
c******************************************************
c--   (D): calculate relative differences between gradients
c******************************************************

                     if ( grdchk_eps .eq. 0. ) then
                        gfd = (fcpertplus-fcpertminus)
                     else
                        gfd = (fcpertplus-fcpertminus)
     &                       /(grdchk_epsfac*grdchk_eps)
                     endif

                     if ( adxxmemo .eq. 0. ) then
                        ratio_ad = abs( adxxmemo - gfd )
                     else
                        ratio_ad = 1. - gfd/adxxmemo
                     endif
                  
                     if ( ftlxxmemo .eq. 0. ) then
                        ratio_ftl = abs( ftlxxmemo - gfd )
                     else
                        ratio_ftl = 1. - gfd/ftlxxmemo
                     endif
                  
                     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
                  
cph(
                     print *, 'ph-grd 3 -------------------------------'
                     print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
     &                    ichknum,itilepos,jtilepos,layer,
     &                    fcref, fcpertplus, fcpertminus
                     print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
     &                    ichknum,ichkmem(ichknum),
     &                    icompmem(ichknum),itestmem(ichknum),
     &                    adxxmemo, gfd, ratio_ad
                     print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
     &                    ichknum,ichkmem(ichknum),
     &                    icompmem(ichknum),itestmem(ichknum),
     &                    ftlxxmemo, gfd, ratio_ftl
cph)
                  else
c
                     print *, 'ph-grd 3 -------------------------------'
                     print *, 'ph-grd 3 : ierr = ', ierr, 
     &                                 ', icomp = ', icomp
                  endif
               else
                  ierr_grdchk = -1
               endif
            
            enddo
         endif

c--   Everyone has to wait for the component to be reset.
         _BARRIER

      enddo

      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)

c--   Print the results of the gradient check.
      call grdchk_print( ichknum, ierr_grdchk, mythid )

#endif /* ALLOW_GRADIENT_CHECK */

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