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