/[MITgcm]/MITgcm/pkg/grdchk/grdchk_main.F
ViewVC logotype

Diff of /MITgcm/pkg/grdchk/grdchk_main.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

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