/[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.21 by heimbach, Fri Apr 20 19:26:43 2007 UTC revision 1.43 by jmc, Fri Apr 4 21:39:04 2014 UTC
# Line 1  Line 1 
 C  
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "GRDCHK_OPTIONS.h"
5  #include "AD_CONFIG.h"  #include "AD_CONFIG.h"
6  #include "CPP_OPTIONS.h"  #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 12  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 125  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 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.  
161    
162          ierr      = 0
163        ierr_grdchk = 0        ierr_grdchk = 0
164  #ifdef ALLOW_ADMTLM        adxxmemo  = 0.
165          ftlxxmemo = 0.
166    #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        print *, 'grad-res -------------------------------'        WRITE(standardmessageunit,'(A)')
198        print '(2a)',       &    'grad-res -------------------------------'
199          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        print '(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        print '(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    
223  cph(          IF ( ichknum .LE. maxgrdchecks ) THEN
224  cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ',            WRITE(msgBuf,'(A,I4,A)')
225  cph-print     &        nbeg, icomp, ichknum       &      '====== Starts gradient-check number', ichknum,
226  cph)       &                                         ' (=ichknum) ======='
227           if (ichknum .le. maxgrdchecks ) then            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
228    
229  c--         Determine the location of icomp on the grid.  C--       Determine the location of icomp on the grid.
230              if ( myProcId .EQ. grdchkwhichproc ) then            IF ( myProcId .EQ. grdchkwhichproc ) THEN
231                 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              endif              layer    = 0
240              _BARRIER              obcspos  = 0
241                              itilepos = 0
242  c******************************************************              jtilepos = 0
243  c--   (A): get gradient component calculated via adjoint              icglo    = 0
244  c******************************************************              itest    = 0
245              ENDIF
246  c--   get gradient component calculated via adjoint  
247              if ( myProcId .EQ. grdchkwhichproc .AND.  C make sure that all procs have correct file records, so that useSingleCpuIO works
248       &           ierr .EQ. 0 ) then            CALL GLOBAL_SUM_INT( icvrec , myThid )
249                 call grdchk_getadxx( icvrec,            CALL GLOBAL_SUM_INT( layer , myThid )
250              CALL GLOBAL_SUM_INT( ierr , myThid )
251    
252              WRITE(msgBuf,'(A,3I5,A,2I4,A,I3,A,I4)')
253         &     'grdchk pos: i,j,k=', itilepos, jtilepos, layer,
254         &              ' ; bi,bj=', itile, jtile, ' ; iobc=', obcspos,
255         &              ' ; rec=', icvrec
256              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
257    
258    C******************************************************
259    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 )
267              endif  C--   Add a global-sum call so that all proc will get the adjoint gradient
268              _BARRIER            _GLOBAL_SUM_RL( adxxmemo, myThid )
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              endif  
284              _BARRIER  C--   2. perform tangent linear run
285              myTime = startTime
286  c--            myIter = nIter0
 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              endif  
329              _BARRIER  C--   forward run with perturbed control vector
330                                  myTime = startTime
331  c--   forward run with perturbed control vector            myIter = nIter0
332              mytime = starttime  #ifdef ALLOW_OPENAD
333              myiter = niter0            CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
             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
337              print *, 'ph-check fcpertplus = ', fcpertplus  
338              _BARRIER  #if (defined (ALLOW_ADMTLM))
339                                fcpertplus = objf_state_final(idep,jdep,1,1,1)
340  c--   Reset control vector.  #elif (defined (ALLOW_OPENAD))
341              if ( myProcId .EQ. grdchkwhichproc .AND.            fcpertplus = fc%v
342       &           ierr .EQ. 0 ) then  #else
343                 call grdchk_setxx( icvrec, FORWARD_SIMULATION,            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                 endif  
368                 _BARRIER  C--   forward run with perturbed control vector
369                                    myTime = startTime
370  c--   forward run with perturbed control vector              myIter = nIter0
371                 mytime = starttime  #ifdef ALLOW_OPENAD
372                 myiter = niter0              CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
373                 call the_main_loop( mytime, myiter, mythid )  #else
374                 _BARRIER              CALL THE_MAIN_LOOP( myTime, myIter, myThid )
375  #ifdef ALLOW_ADMTLM  #endif
376                 fcpertminus = objf_state_final(idep,jdep,1,1,1)  
377    #if (defined (ALLOW_ADMTLM))
378                fcpertminus = objf_state_final(idep,jdep,1,1,1)
379    #elif (defined (ALLOW_OPENAD))
380                fcpertminus = fc%v
381  #else  #else
382                 fcpertminus = fc              fcpertminus = fc
383  #endif  #endif
384                                  WRITE(msgBuf,'(A,1PE22.14)')
385  c--   Reset control vector.       &         'grdchk perturb(-)fc: fcpertminus =', fcpertminus
386                 if ( myProcId .EQ. grdchkwhichproc .AND.              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
387       &           ierr .EQ. 0 ) then  
388                    call grdchk_setxx( icvrec, FORWARD_SIMULATION,  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 ( myProcId .EQ. grdchkwhichproc .AND.            ELSE
404       &           ierr .EQ. 0 ) then              gfd = (fcpertplus-fcpertminus) /(grdchk_epsfac*grdchk_eps)
405              ENDIF
406                 if ( grdchk_eps .eq. 0. ) then  
407                    gfd = (fcpertplus-fcpertminus)            IF ( adxxmemo .EQ. 0. ) THEN
408                 else              ratio_ad = ABS( adxxmemo - gfd )
409                    gfd = (fcpertplus-fcpertminus)            ELSE
410       &                 /(grdchk_epsfac*grdchk_eps)              ratio_ad = 1. - gfd/adxxmemo
411                 endif            ENDIF
412                        
413                 if ( adxxmemo .eq. 0. ) then            IF ( ftlxxmemo .EQ. 0. ) THEN
414                    ratio_ad = abs( adxxmemo - gfd )              ratio_ftl = ABS( ftlxxmemo - gfd )
415                 else            ELSE
416                    ratio_ad = 1. - gfd/adxxmemo              ratio_ftl = 1. - gfd/ftlxxmemo
417                 endif            ENDIF
418    
419                 if ( ftlxxmemo .eq. 0. ) then            IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
420                    ratio_ftl = abs( ftlxxmemo - gfd )              tmpplot1(itilepos,jtilepos,layer,itile,jtile) = gfd
421                 else              tmpplot2(itilepos,jtilepos,layer,itile,jtile) = ratio_ad
422                    ratio_ftl = 1. - gfd/ftlxxmemo              tmpplot3(itilepos,jtilepos,layer,itile,jtile) = ratio_ftl
423                 endif            ENDIF
424                      
425                 tmpplot1(itilepos,jtilepos,layer,itile,jtile)            IF ( ierr .EQ. 0 ) THEN
426       &              = gfd              fcrmem      ( ichknum ) = fcref
427                 tmpplot2(itilepos,jtilepos,layer,itile,jtile)              fcppmem     ( ichknum ) = fcpertplus
428       &              = ratio_ad              fcpmmem     ( ichknum ) = fcpertminus
429                 tmpplot3(itilepos,jtilepos,layer,itile,jtile)              xxmemref    ( ichknum ) = xxmemo_ref
430       &              = ratio_ftl              xxmempert   ( ichknum ) = xxmemo_pert
431                gfdmem      ( ichknum ) = gfd
432                 fcrmem      ( ichknum ) = fcref              adxxmem     ( ichknum ) = adxxmemo
433                 fcppmem     ( ichknum ) = fcpertplus              ftlxxmem    ( ichknum ) = ftlxxmemo
434                 fcpmmem     ( ichknum ) = fcpertminus              ratioadmem  ( ichknum ) = ratio_ad
435                 xxmemref    ( ichknum ) = xxmemo_ref              ratioftlmem ( ichknum ) = ratio_ftl
436                 xxmempert   ( ichknum ) = xxmemo_pert  
437                 gfdmem      ( ichknum ) = gfd              irecmem   ( ichknum ) = icvrec
438                 adxxmem     ( ichknum ) = adxxmemo              bimem     ( ichknum ) = itile
439                 ftlxxmem    ( ichknum ) = ftlxxmemo              bjmem     ( ichknum ) = jtile
440                 ratioadmem  ( ichknum ) = ratio_ad              ilocmem   ( ichknum ) = itilepos
441                 ratioftlmem ( ichknum ) = ratio_ftl              jlocmem   ( ichknum ) = jtilepos
442                klocmem   ( ichknum ) = layer
443                 irecmem   ( ichknum ) = icvrec              iobcsmem  ( ichknum ) = obcspos
444                 bimem     ( ichknum ) = itile              icompmem  ( ichknum ) = icomp
445                 bjmem     ( ichknum ) = jtile              ichkmem   ( ichknum ) = ichknum
446                 ilocmem   ( ichknum ) = itilepos              itestmem  ( ichknum ) = itest
447                 jlocmem   ( ichknum ) = jtilepos              ierrmem   ( ichknum ) = ierr
448                 klocmem   ( ichknum ) = layer              icglomem  ( ichknum ) = icglo
449                 iobcsmem  ( ichknum ) = obcspos            ENDIF
450                 icompmem  ( ichknum ) = icomp  
451                 ichkmem   ( ichknum ) = ichknum            IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
452                 itestmem  ( ichknum ) = itest              WRITE(standardmessageunit,'(A)')
453                 ierrmem   ( ichknum ) = ierr       &          'grad-res -------------------------------'
454                 icglomem  ( ichknum ) = icglo              WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
455         &           ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
456                 print *, 'grad-res -------------------------------'       &             layer,itile,jtile,obcspos,
457                 print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',       &             fcref, fcpertplus, fcpertminus
      &              myprocid,ichknum,itilepos,jtilepos,  
      &              layer,itile,jtile,obcspos,  
      &              fcref, fcpertplus, fcpertminus  
458  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
459                 print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',              WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
460       &              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
                WRITE(msgBuf,'(A34,2(1PE24.14,X))')  
      &              'precision_grdchk_result TLM ', fcref, ftlxxmemo  
                CALL PRINT_MESSAGE  
      &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)  
464  #else  #else
465                 print '(a,8I5,2x,3(1x,E18.12))', ' grad-res ',              WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
466       &              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                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')  #endif
471       &              'precision_grdchk_result ADM ', fcref, adxxmemo            ENDIF
472                 CALL PRINT_MESSAGE  #ifdef ALLOW_TANGENTLINEAR_RUN
473       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)            WRITE(msgBuf,'(A30,1PE22.14)')
474  #endif       &              ' TLM  ref_cost_function      =', fcref
475              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
476              endif            WRITE(msgBuf,'(A30,1PE22.14)')
477         &              ' TLM  tangent-lin_grad       =', ftlxxmemo
478              print *, 'ph-grd  ierr ---------------------------'            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
479              print *, 'ph-grd  ierr = ', ierr, ', icomp = ', icomp,            WRITE(msgBuf,'(A30,1PE22.14)')
480       &           ', ichknum = ', ichknum       &              ' TLM  finite-diff_grad       =', gfd
481              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
482              _BARRIER  #else
483              WRITE(msgBuf,'(A30,1PE22.14)')
484  c-- else of if ( ichknum ...       &              ' ADM  ref_cost_function      =', fcref
485           else            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
486              ierr_grdchk = -1            WRITE(msgBuf,'(A30,1PE22.14)')
487         &              ' ADM  adjoint_gradient       =', adxxmemo
488  c-- end of if ( ichknum ...            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
489           endif            WRITE(msgBuf,'(A30,1PE22.14)')
490         &              ' ADM  finite-diff_grad       =', gfd
491  c-- end of do icomp = ...            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
492        enddo  #endif
   
       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  
493    
494  c--   Everyone has to wait for the component to be reset.            WRITE(msgBuf,'(A,I4,A,I3,A)')
495        _BARRIER       &      '====== End of gradient-check number', ichknum,
496         &                      ' (ierr=', ierr, ') ======='
497              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
498    
499    C-- else of if ichknum < ...
500            ELSE
501              ierr_grdchk = -1
502    
503    C-- end of if ichknum < ...
504            ENDIF
505    
506    C-- end of do icomp = ...
507          ENDDO
508    
509          IF (myProcId .EQ. grdchkwhichproc .AND. .NOT.useSingleCpuIO) THEN
510            CALL WRITE_REC_XYZ_RL(
511         &       'grd_findiff'   , tmpplot1, 1, 0, myThid)
512            CALL WRITE_REC_XYZ_RL(
513         &       'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
514            CALL WRITE_REC_XYZ_RL(
515         &       'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
516          ENDIF
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        end        RETURN
524          END
Legend:
Removed from v.1.21  
changed lines
  Added in v.1.43
  ViewVC Help
Powered by ViewVC 1.1.22