pkg/grdchk/grdchk_main.F

C
C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.18 2006/11/10 04:57:10 heimbach Exp $
C $Name:  $

#include "AD_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOI
C
C !TITLE: GRADIENT CHECK
C !AUTHORS: mitgcm developers ( support@mitgcm.org )
C !AFFILIATION: Massachussetts Institute of Technology
C !DATE:
C !INTRODUCTION: gradient check package
c \bv
c Compare the gradients calculated by the adjoint model with
c finite difference approximations.
c
C     !CALLING SEQUENCE:
c
c the_model_main
c |
c |-- ctrl_unpack
c |-- adthe_main_loop          - unperturbed cost function and
c |-- ctrl_pack                  adjoint gradient are computed here
c |
c |-- grdchk_main
c     |
c     |-- grdchk_init
c     |-- do icomp=...        - loop over control vector elements
c         |
c         |-- grdchk_loc      - determine location of icomp on grid
c         |
c         |-- grdchk_getxx    - get control vector component from file
c         |                     perturb it and write back to file
c         |-- grdchk_getadxx  - get gradient component calculated 
c         |                     via adjoint
c         |-- the_main_loop   - forward run and cost evaluation
c         |                     with perturbed control vector element
c         |-- calculate ratio of adj. vs. finite difference gradient
c         |
c         |-- grdchk_setxx    - Reset control vector element
c         |
c         |-- grdchk_print    - print results
c \ev
CEOI

CBOP
C     !ROUTINE: grdchk_main
C     !INTERFACE:
      subroutine grdchk_main( mythid )

C     !DESCRIPTION: \bv
c     ==================================================================
c     SUBROUTINE grdchk_main
c     ==================================================================
c     o Compare the gradients calculated by the adjoint model with
c       finite difference approximations.
c     started: Christian Eckert eckert@mit.edu 24-Feb-2000
c     continued&finished: heimbach@mit.edu: 13-Jun-2001
c     changed: mlosch@ocean.mit.edu: 09-May-2002
c              - added centered difference vs. 1-sided difference option
c              - improved output format for readability
c              - added control variable hFacC
c              heimbach@mit.edu 24-Feb-2003
c              - added tangent linear gradient checks
c              - fixes for multiproc. gradient checks
c              - added more control variables
c     
c     ==================================================================
c     SUBROUTINE grdchk_main
c     ==================================================================
C     \ev

C     !USES:
      implicit none

c     == global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "grdchk.h"
#include "cost.h"
#ifdef ALLOW_TANGENTLINEAR_RUN
#include "g_cost.h"
#endif

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

#ifdef ALLOW_GRDCHK
C     !LOCAL VARIABLES:
c     == local variables ==
      integer myiter
      _RL     mytime 
      integer bi, itlo, ithi
      integer bj, jtlo, jthi
      integer i,  imin, imax
      integer j,  jmin, jmax
      integer k

      integer icomp
      integer ichknum
      integer icvrec
      integer jtile
      integer itile
      integer layer
      integer obcspos
      integer itilepos
      integer jtilepos
      integer icglo
      integer itest
      integer ierr
      integer ierr_grdchk
      _RL     gfd
      _RL     fcref
      _RL     fcpertplus, fcpertminus
      _RL     ratio_ad
      _RL     ratio_ftl
      _RL     xxmemo_ref
      _RL     xxmemo_pert
      _RL     adxxmemo
      _RL     ftlxxmemo
      _RL     localEps
      _RL     grdchk_epsfac

      _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)

      CHARACTER*(MAX_LEN_MBUF) msgBuf

c     == end of interface ==
CEOP

c--   Set the loop ranges.
      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)
      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

      print *, 'ph-check entering grdchk_main '

c--   initialise variables
      call grdchk_init( mythid )

c--   Compute the adjoint models' gradients.
c--   Compute the unperturbed cost function.
cph   Gradient via adjoint has already been computed,
cph   and so has unperturbed cost function,
cph   assuming all xx_ fields are initialised to zero.

      ierr_grdchk = 0
cphadmtlm(
      fcref = fc
cphadmtlm      fcref = objf_state_final(45,4,1,1,1)
cphadmtlm)

      print *, 'ph-check fcref = ', fcref

      do bj = jtlo, jthi
         do bi = itlo, ithi
            do k = 1, nr
               do j = jmin, jmax
                  do i = imin, imax
                     tmpplot1(i,j,k,bi,bj) = 0. _d 0
                     tmpplot2(i,j,k,bi,bj) = 0. _d 0
                     tmpplot3(i,j,k,bi,bj) = 0. _d 0
                  end do
               end do
            end do
         end do
      end do

      if ( useCentralDiff ) then
         grdchk_epsfac = 2. _d 0
      else
         grdchk_epsfac = 1. _d 0
      end if

      print *, 'grad-res -------------------------------'
      print ('(2a)'), 
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '               fc ref           fc + eps           fc - eps'
#ifdef ALLOW_TANGENTLINEAR_RUN
      print ('(2a)'), 
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '             tlm grad            fd grad         1 - fd/tlm'
#else
      print ('(2a)'), 
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '             adj grad            fd grad         1 - fd/adj'
#endif

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

      if ( nbeg .EQ. 0 ) 
     &     call grdchk_get_position( mythid )

      do icomp = nbeg, nend, nstep

         ichknum = (icomp - nbeg)/nstep + 1

cph(
cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ', 
cph-print     &        nbeg, icomp, ichknum
cph)
         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, icglo, itest, ierr,
     &              mythid )
cph(
cph-print               print *, 'ph-grd ----- back from loc -----',
cph-print     &             icvrec, itilepos, jtilepos, layer, obcspos
cph)
            endif
            _BARRIER
              
c******************************************************
c--   (A): get gradient component calculated via adjoint
c******************************************************

c--   get gradient component calculated via adjoint
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_getadxx( icvrec,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              adxxmemo, mythid )
            endif
            _BARRIER

#ifdef ALLOW_TANGENTLINEAR_RUN
c******************************************************
c--   (B): Get gradient component g_fc from tangent linear run:
c******************************************************
c--
c--   1. perturb control vector component: xx(i)=1.

            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               localEps = 1. _d 0
               call grdchk_getxx( icvrec, TANGENT_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, xxmemo_pert, localEps,
     &              mythid )
            endif
            _BARRIER

c--
c--   2. perform tangent linear run
            mytime = starttime
            myiter = niter0
cphadmtlm(
            g_fc = 0.
cphadmtlm            do j=1,sny
cphadmtlm               do i=1,snx
cphadmtlm                  g_objf_state_final(i,j,1,1,1) = 0.
cphadmtlm                  g_objf_state_final(i,j,1,1,2) = 0.
cphadmtlm               enddo
cphadmtlm            enddo
cphadmtlm)
            call g_the_main_loop( mytime, myiter, mythid )
cphadmtlm(
            ftlxxmemo = g_fc
cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1,1)
cphadmtlm)
            _BARRIER
c--
c--   3. reset control vector
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_setxx( icvrec, TANGENT_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, mythid )
            endif
            _BARRIER

#endif /* ALLOW_TANGENTLINEAR_RUN */


c******************************************************
c--   (C): Get gradient via finite difference perturbation
c******************************************************

c--   get control vector component from file
c--   perturb it and write back to file
c--   positive perturbation
            localEps = abs(grdchk_eps)
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_getxx( icvrec, FORWARD_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, xxmemo_pert, localEps,
     &              mythid )
            endif
            _BARRIER
                     
c--   forward run with perturbed control vector
            mytime = starttime
            myiter = niter0
            call the_main_loop( mytime, myiter, mythid )
cphadmtlm(
            fcpertplus = fc
cphadmtlm            fcpertplus = objf_state_final(45,4,1,1,1)
cphadmtlm)
            print *, 'ph-check fcpertplus = ', fcpertplus
            _BARRIER
                  
c--   Reset control vector.
            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
               call grdchk_setxx( icvrec, FORWARD_SIMULATION,
     &              itile, jtile, layer,
     &              itilepos, jtilepos,
     &              xxmemo_ref, mythid )
            endif
            _BARRIER

            fcpertminus = fcref
            print *, 'ph-check fcpertminus = ', fcpertminus

            if ( useCentralDiff ) then

c--   get control vector component from file
c--   perturb it and write back to file
c--   repeat the proceedure for a negative perturbation
               if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
                  localEps = - abs(grdchk_eps)
                  call grdchk_getxx( icvrec, FORWARD_SIMULATION,
     &                 itile, jtile, layer,
     &                 itilepos, jtilepos,
     &                 xxmemo_ref, xxmemo_pert, localEps,
     &                 mythid )
               endif
               _BARRIER
                     
c--   forward run with perturbed control vector
               mytime = starttime
               myiter = niter0
               call the_main_loop( mytime, myiter, mythid )
               _BARRIER
               fcpertminus = fc
                  
c--   Reset control vector.
               if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then
                  call grdchk_setxx( icvrec, FORWARD_SIMULATION,
     &                 itile, jtile, layer,
     &                 itilepos, jtilepos,
     &                 xxmemo_ref, mythid )
               endif
               _BARRIER

c-- end of if useCentralDiff ...
            end if

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

            if ( myProcId .EQ. grdchkwhichproc .AND.
     &           ierr .EQ. 0 ) then

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

               if ( ftlxxmemo .eq. 0. ) then
                  ratio_ftl = abs( ftlxxmemo - gfd )
               else
                  ratio_ftl = 1. - gfd/ftlxxmemo
               endif
                  
               tmpplot1(itilepos,jtilepos,layer,itile,jtile)
     &              = gfd
               tmpplot2(itilepos,jtilepos,layer,itile,jtile)
     &              = ratio_ad
               tmpplot3(itilepos,jtilepos,layer,itile,jtile)
     &              = ratio_ftl

               fcrmem      ( ichknum ) = fcref
               fcppmem     ( ichknum ) = fcpertplus
               fcpmmem     ( ichknum ) = fcpertminus
               xxmemref    ( ichknum ) = xxmemo_ref
               xxmempert   ( ichknum ) = xxmemo_pert
               gfdmem      ( ichknum ) = gfd
               adxxmem     ( ichknum ) = adxxmemo
               ftlxxmem    ( ichknum ) = ftlxxmemo
               ratioadmem  ( ichknum ) = ratio_ad
               ratioftlmem ( ichknum ) = ratio_ftl

               irecmem   ( ichknum ) = icvrec
               bimem     ( ichknum ) = itile
               bjmem     ( ichknum ) = jtile
               ilocmem   ( ichknum ) = itilepos
               jlocmem   ( ichknum ) = jtilepos
               klocmem   ( ichknum ) = layer
               iobcsmem  ( ichknum ) = obcspos
               icompmem  ( ichknum ) = icomp
               ichkmem   ( ichknum ) = ichknum
               itestmem  ( ichknum ) = itest
               ierrmem   ( ichknum ) = ierr
               icglomem  ( ichknum ) = icglo

               print *, 'grad-res -------------------------------'
               print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',
     &              myprocid,ichknum,itilepos,jtilepos,
     &              layer,itile,jtile,obcspos,
     &              fcref, fcpertplus, fcpertminus
#ifdef ALLOW_TANGENTLINEAR_RUN
               print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',
     &              myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),
     &              bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
     &              ftlxxmemo, gfd, ratio_ftl
               WRITE(msgBuf,'(A34,2(1PE24.14,X))') 
     &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
               CALL PRINT_MESSAGE
     &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
#else
               print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',
     &              myprocid,ichknum,ichkmem(ichknum),
     &              icompmem(ichknum),itestmem(ichknum),
     &              bimem(ichknum),bjmem(ichknum),obcspos,
     &              adxxmemo, gfd, ratio_ad
               WRITE(msgBuf,'(A34,2(1PE24.14,X))') 
     &              'precision_grdchk_result ADM ', fcref, adxxmemo
               CALL PRINT_MESSAGE
     &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
#endif

            endif

            print *, 'ph-grd  ierr ---------------------------'
            print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,
     &           ', ichknum = ', ichknum

            _BARRIER

c-- else of if ( ichknum ...
         else
            ierr_grdchk = -1

c-- end of if ( ichknum ... 
         endif

c-- end of do icomp = ... 
      enddo

      if ( myProcId .EQ. grdchkwhichproc ) then
         CALL WRITE_REC_XYZ_RL( 
     &        'grd_findiff'   , tmpplot1, 1, 0, myThid)
         CALL WRITE_REC_XYZ_RL( 
     &        'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
         CALL WRITE_REC_XYZ_RL( 
     &        'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
      endif

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

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

#endif /* ALLOW_GRDCHK */

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