pkg/grdchk/grdchk_main.F

C
C $Header: /u/gcmpack/MITgcm/pkg/grdchk/grdchk_main.F,v 1.25 2007/09/15 17:32:41 jmc 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"
#include "ctrl.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
#ifdef ALLOW_ADMTLM
      fcref = objf_state_final(idep,jdep,1,1,1)
#else
      fcref = fc
#endif

      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

      WRITE(standardmessageunit,'(A)')
     &    'grad-res -------------------------------'
      WRITE(standardmessageunit,'(2a)')
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '               fc ref           fc + eps           fc - eps'
#ifdef ALLOW_TANGENTLINEAR_RUN
      WRITE(standardmessageunit,'(2a)')
     &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
     &    '             tlm grad            fd grad         1 - fd/tlm'
#else
      WRITE(standardmessageunit,'(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.
            ftlxxmemo = 0.

            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
#ifdef ALLOW_ADMTLM
            do k=1,4*Nr+1
             do j=1,sny
              do i=1,snx
               g_objf_state_final(i,j,1,1,k) = 0.
              enddo
             enddo
            enddo
#else
            g_fc = 0.
#endif

            call g_the_main_loop( mytime, myiter, mythid )
            _BARRIER
#ifdef ALLOW_ADMTLM
            ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
#else
            ftlxxmemo = g_fc
#endif

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 )
#ifdef ALLOW_ADMTLM
            fcpertplus = objf_state_final(idep,jdep,1,1,1)
#else
            fcpertplus = fc
#endif
            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
#ifdef ALLOW_ADMTLM
               fcpertminus = objf_state_final(idep,jdep,1,1,1)
#else
               fcpertminus = fc
#endif
                  
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

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