/[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.6 by heimbach, Mon Sep 16 18:11:58 2002 UTC revision 1.41 by jmc, Fri Aug 24 18:43:41 2012 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "GRDCHK_OPTIONS.h"
5    #include "AD_CONFIG.h"
6    
7  CBOI  CBOI
8  C  C
# Line 9  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      C              heimbach@mit.edu 24-Feb-2003
65  c     ==================================================================  C              - added tangent linear gradient checks
66  c     SUBROUTINE grdchk_main  C              - fixes for multiproc. gradient checks
67  c     ==================================================================  C              - added more control variables
68    C     ==================================================================
69    C     SUBROUTINE GRDCHK_MAIN
70    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"
80  #include "grdchk.h"  #include "grdchk.h"
81  #include "cost.h"  #include "cost.h"
82    #include "ctrl.h"
83    #ifdef ALLOW_TANGENTLINEAR_RUN
84    #include "g_cost.h"
85    #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_GRADIENT_CHECK  #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 jprocs  
103        integer proc_grdchk        INTEGER icomp
104        integer icomp        INTEGER ichknum
105        integer ichknum        INTEGER icvrec
106        integer icvrec        INTEGER jtile
107        integer jtile        INTEGER itile
108        integer itile        INTEGER layer
109        integer layer        INTEGER obcspos
110        integer itilepos        INTEGER itilepos
111        integer jtilepos        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 114  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          ioUnit = standardMessageUnit
134          WRITE(msgBuf,'(A)')
135         &'// ======================================================='
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  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
144        _RL     g_fc  C--   prevent writing output multiple times
145        common /g_cost_r/ g_fc  C     note: already called in AD run so that only needed for TLM run
146          CALL TURNOFF_MODEL_IO( 0, myThid )
147  #endif  #endif
148    
149        _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)  C--   initialise variables
150        _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        CALL GRDCHK_INIT( myThid )
       _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)  
   
 c     == end of interface ==  
 CEOP  
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  
   
 c--   initialise variables  
       call grdchk_init( mythid )  
   
 c--   Compute the adjoint models' gradients.  
 c--   Compute the unperturbed cost function.  
 c--   Gradient via adjoint has already been computed,  
 c--   and so has unperturbed cost function,  
 c--   assuming all xx_ fields are initialised to zero.  
157    
158        fcref = fc        ierr      = 0
159        ierr_grdchk = 0        ierr_grdchk = 0
160          adxxmemo  = 0.
161          ftlxxmemo = 0.
162    #if (defined  (ALLOW_ADMTLM))
163          fcref = objf_state_final(idep,jdep,1,1,1)
164    #elif (defined (ALLOW_AUTODIFF_OPENAD))
165          fcref = fc%v
166    #else
167          fcref = fc
168    #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        do bj = jtlo, jthi        IF ( useCentralDiff ) THEN
          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)')
194         &    'grad-res -------------------------------'
195          WRITE(standardmessageunit,'(2a)')
196         &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
197         &    '       fc ref            fc + eps           fc - eps'
198    #ifdef ALLOW_TANGENTLINEAR_RUN
199          WRITE(standardmessageunit,'(2a)')
200         &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
201         &    '      tlm grad            fd grad          1 - fd/tlm'
202    #else
203          WRITE(standardmessageunit,'(2a)')
204         &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
205         &    '      adj grad            fd grad          1 - fd/adj'
206    #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        do jprocs = 1,numberOfProcs  
212           proc_grdchk = jprocs - 1        IF ( nbeg .EQ. 0 ) CALL GRDCHK_GET_POSITION( myThid )
213    
214           if ( myProcId .eq. proc_grdchk ) then        DO icomp = nbeg, nend, nstep
215    
216              do icomp = nbeg, nend, nstep          adxxmemo  = 0.
217            ichknum = (icomp - nbeg)/nstep + 1
218                 ichknum = (icomp - nbeg)/nstep + 1  
219            IF ( ichknum .LE. maxgrdchecks ) THEN
220                 if (ichknum .le. maxgrdchecks ) then            WRITE(msgBuf,'(A,I4,A)')
221         &      '====== Starts gradient-check number', ichknum,
222  c--   Determine the location of icomp on the grid.       &                                         ' (=ichknum) ======='
223                    call grdchk_loc( icomp, ichknum,            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
224       &                 icvrec, itile, jtile, layer,  
225       &                 itilepos, jtilepos, itest, ierr,  C--       Determine the location of icomp on the grid.
226       &                 mythid )            IF ( myProcId .EQ. grdchkwhichproc ) THEN
227                              CALL grdchk_loc( icomp, ichknum,
228                    if ( ierr .eq. 0 ) then       &              icvrec, itile, jtile, layer, obcspos,
229         &              itilepos, jtilepos, icglo, itest, ierr,
230  c******************************************************       &              myThid )
231  c--   (A): get gradient component calculated via adjoint            ELSE
232  c******************************************************              icvrec   = 0
233                       call grdchk_getadxx( icvrec,              itile    = 0
234       &                    itile, jtile, layer,              jtile    = 0
235       &                    itilepos, jtilepos,              layer    = 0
236       &                    adxxmemo, mythid )              obcspos  = 0
237                       _BARRIER              itilepos = 0
238                jtilepos = 0
239                icglo    = 0
240                itest    = 0
241              ENDIF
242    
243    C make sure that all procs have correct file records, so that useSingleCpuIO works
244              CALL GLOBAL_SUM_INT( icvrec , myThid )
245              CALL GLOBAL_SUM_INT( layer , myThid )
246              CALL GLOBAL_SUM_INT( ierr , myThid )
247    
248              WRITE(msgBuf,'(A,3I5,A,2I4,A,I3,A,I4)')
249         &     'grdchk pos: i,j,k=', itilepos, jtilepos, layer,
250         &              ' ; bi,bj=', itile, jtile, ' ; iobc=', obcspos,
251         &              ' ; rec=', icvrec
252              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
253    
254    C******************************************************
255    C--   (A): get gradient component calculated via adjoint
256    C******************************************************
257    
258    C--   get gradient component calculated via adjoint
259              CALL GRDCHK_GETADXX( icvrec,
260         &              itile, jtile, layer,
261         &              itilepos, jtilepos,
262         &              adxxmemo, ierr, myThid )
263    C--   Add a global-sum call so that all proc will get the adjoint gradient
264              _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.            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       &                    mythid )       &              ierr, myThid )
279                       _BARRIER  
280    C--   2. perform tangent linear run
281  c--            myTime = startTime
282  c--   2. perform tangent linear run            myIter = nIter0
283                       mytime = starttime  #ifdef ALLOW_ADMTLM
284                       myiter = niter0            DO k=1,4*Nr+1
285                       g_fc = 0.              DO j=1,sNy
286                       call g_the_main_loop( mytime, myiter, mythid )                DO i=1,sNx
287                       ftlxxmemo = g_fc                  g_objf_state_final(i,j,1,1,k) = 0.
288  c--                ENDDO
289  c--   3. reset control vector              ENDDO
290                       call grdchk_setxx( icvrec, TANGENT_SIMULATION,            ENDDO
291       &                    itile, jtile, layer,  #else
292       &                    itilepos, jtilepos,            g_fc = 0.
293       &                    xxmemo_ref, mythid )  #endif
294                       _BARRIER  
295  #endif            CALL G_THE_MAIN_LOOP( myTime, myIter, myThid )
296    #ifdef ALLOW_ADMTLM
297  c******************************************************            ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
298  c--   (C): Get gradient via finite difference perturbation  #else
299  c******************************************************            ftlxxmemo = g_fc
300    #endif
301  c--   get control vector component from file  
302  c--   perturb it and write back to file:  C--   3. reset control vector
303  c--   positive perturbation            CALL GRDCHK_SETXX( icvrec, TANGENT_SIMULATION,
304                       localEps = abs(grdchk_eps)       &              itile, jtile, layer,
305                       call grdchk_getxx( icvrec, FORWARD_SIMULATION,       &              itilepos, jtilepos,
306       &                    itile, jtile, layer,       &              xxmemo_ref, ierr, myThid )
307       &                    itilepos, jtilepos,  
308       &                    xxmemo_ref, xxmemo_pert, localEps,  #endif /* ALLOW_TANGENTLINEAR_RUN */
309       &                    mythid )  
310                       _BARRIER  C******************************************************
311                        C--   (C): Get gradient via finite difference perturbation
312  c--   forward run with perturbed control vector  C******************************************************
313                       mytime = starttime  
314                       myiter = niter0  C--   (C-1) positive perturbation:
315                       call the_main_loop( mytime, myiter, mythid )            localEps = ABS(grdchk_eps)
316                       fcpertplus = fc  
317                      C--   get control vector component from file
318  c--   Reset control vector.  C--   perturb it and write back to file
319                       call grdchk_setxx( icvrec, FORWARD_SIMULATION,            CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
320       &                    itile, jtile, layer,       &              itile, jtile, layer,
321       &                    itilepos, jtilepos,       &              itilepos, jtilepos,
322       &                    xxmemo_ref, mythid )       &              xxmemo_ref, xxmemo_pert, localEps,
323                       _BARRIER       &              ierr, myThid )
324    
325                       fcpertminus = fcref  C--   forward run with perturbed control vector
326              myTime = startTime
327                       if ( useCentralDiff ) then            myIter = nIter0
328    #ifdef ALLOW_AUTODIFF_OPENAD
329  c--   get control vector component from file            CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
330  c--   perturb it and write back to file:  #else
331  c--   repeat the proceedure for a negative perturbation            CALL THE_MAIN_LOOP( myTime, myIter, myThid )
332                          localEps = - abs(grdchk_eps)  #endif
333                          call grdchk_getxx( icvrec, FORWARD_SIMULATION,  
334       &                    itile, jtile, layer,  #if (defined (ALLOW_ADMTLM))
335       &                    itilepos, jtilepos,            fcpertplus = objf_state_final(idep,jdep,1,1,1)
336       &                    xxmemo_ref, xxmemo_pert, localEps,  #elif (defined (ALLOW_AUTODIFF_OPENAD))
337       &                    mythid )            fcpertplus = fc%v
338                          _BARRIER  #else
339                                  fcpertplus = fc
340  c--   forward run with perturbed control vector  #endif
341                          mytime = starttime            WRITE(msgBuf,'(A,1PE22.14)')
342                          myiter = niter0       &         'grdchk perturb(+)fc: fcpertplus  =', fcpertplus
343                          call the_main_loop( mytime, myiter, mythid )            CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
344                          fcpertminus = fc  
345                      C--   Reset control vector.
346  c--   Reset control vector.            CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
347                          call grdchk_setxx( icvrec, FORWARD_SIMULATION,       &              itile, jtile, layer,
348       &                    itile, jtile, layer,       &              itilepos, jtilepos,
349       &                    itilepos, jtilepos,       &              xxmemo_ref, ierr, myThid )
350       &                    xxmemo_ref, mythid )  
351                          _BARRIER            fcpertminus = fcref
352              IF ( useCentralDiff ) THEN
353                       end if  C--   (C-2) repeat the proceedure for a negative perturbation:
354  c--              localEps = - ABS(grdchk_eps)
355  c******************************************************  
356  c--   (D): calculate relative differences between gradients  C--   get control vector component from file
357  c******************************************************  C--   perturb it and write back to file
358                CALL GRDCHK_GETXX( icvrec, FORWARD_SIMULATION,
359                       if ( grdchk_eps .eq. 0. ) then       &                 itile, jtile, layer,
360                          gfd = (fcpertplus-fcpertminus)       &                 itilepos, jtilepos,
361                       else       &                 xxmemo_ref, xxmemo_pert, localEps,
362                          gfd = (fcpertplus-fcpertminus)       &                 ierr, myThid )
363       &                       /(grdchk_epsfac*grdchk_eps)  
364                       endif  C--   forward run with perturbed control vector
365                myTime = startTime
366                       if ( adxxmemo .eq. 0. ) then              myIter = nIter0
367                          ratio_ad = abs( adxxmemo - gfd )  #ifdef ALLOW_AUTODIFF_OPENAD
368                       else              CALL OpenAD_THE_MAIN_LOOP( myTime, myIter, myThid )
369                          ratio_ad = 1. - gfd/adxxmemo  #else
370                       endif              CALL THE_MAIN_LOOP( myTime, myIter, myThid )
371                      #endif
372                       if ( ftlxxmemo .eq. 0. ) then  
373                          ratio_ftl = abs( ftlxxmemo - gfd )  #if (defined (ALLOW_ADMTLM))
374                       else              fcpertminus = objf_state_final(idep,jdep,1,1,1)
375                          ratio_ftl = 1. - gfd/ftlxxmemo  #elif (defined (ALLOW_AUTODIFF_OPENAD))
376                       endif              fcpertminus = fc%v
377                      #else
378                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =              fcpertminus = fc
379       &                    gfd  #endif
380                       tmpplot2(itilepos,jtilepos,layer,itile,jtile) =              WRITE(msgBuf,'(A,1PE22.14)')
381       &                    ratio_ad       &         'grdchk perturb(-)fc: fcpertminus =', fcpertminus
382                       tmpplot3(itilepos,jtilepos,layer,itile,jtile) =              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
383       &                    ratio_ftl  
384    C--   Reset control vector.
385                                    CALL GRDCHK_SETXX( icvrec, FORWARD_SIMULATION,
386                       fcrmem      ( ichknum ) = fcref       &                 itile, jtile, layer,
387                       fcppmem     ( ichknum ) = fcpertplus       &                 itilepos, jtilepos,
388                       fcpmmem     ( ichknum ) = fcpertminus       &                 xxmemo_ref, ierr, myThid )
389                       xxmemref    ( ichknum ) = xxmemo_ref  
390                       xxmempert   ( ichknum ) = xxmemo_pert  C-- end of if useCentralDiff ...
391                       gfdmem      ( ichknum ) = gfd            ENDIF
392                       adxxmem     ( ichknum ) = adxxmemo  
393                       ftlxxmem    ( ichknum ) = ftlxxmemo  C******************************************************
394                       ratioadmem  ( ichknum ) = ratio_ad  C--   (D): calculate relative differences between gradients
395                       ratioftlmem ( ichknum ) = ratio_ftl  C******************************************************
396    
397                       irecmem   ( ichknum ) = icvrec            IF ( grdchk_eps .EQ. 0. ) THEN
398                       bimem     ( ichknum ) = itile              gfd = (fcpertplus-fcpertminus)
399                       bjmem     ( ichknum ) = jtile            ELSE
400                       ilocmem   ( ichknum ) = itilepos              gfd = (fcpertplus-fcpertminus) /(grdchk_epsfac*grdchk_eps)
401                       jlocmem   ( ichknum ) = jtilepos            ENDIF
402                       klocmem   ( ichknum ) = layer  
403                       icompmem  ( ichknum ) = icomp            IF ( adxxmemo .EQ. 0. ) THEN
404                       ichkmem   ( ichknum ) = ichknum              ratio_ad = ABS( adxxmemo - gfd )
405                       itestmem  ( ichknum ) = itest            ELSE
406                       ierrmem   ( ichknum ) = ierr              ratio_ad = 1. - gfd/adxxmemo
407                                ENDIF
408  cph(  
409                       print *, 'ph-grd 3 -------------------------------'            IF ( ftlxxmemo .EQ. 0. ) THEN
410                       print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',              ratio_ftl = ABS( ftlxxmemo - gfd )
411       &                    ichknum,itilepos,jtilepos,layer,            ELSE
412       &                    fcref, fcpertplus, fcpertminus              ratio_ftl = 1. - gfd/ftlxxmemo
413                       print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',            ENDIF
414       &                    ichknum,ichkmem(ichknum),  
415       &                    icompmem(ichknum),itestmem(ichknum),            IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
416       &                    adxxmemo, gfd, ratio_ad              tmpplot1(itilepos,jtilepos,layer,itile,jtile) = gfd
417                       print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',              tmpplot2(itilepos,jtilepos,layer,itile,jtile) = ratio_ad
418       &                    ichknum,ichkmem(ichknum),              tmpplot3(itilepos,jtilepos,layer,itile,jtile) = ratio_ftl
419       &                    icompmem(ichknum),itestmem(ichknum),            ENDIF
420       &                    ftlxxmemo, gfd, ratio_ftl  
421  cph)            IF ( ierr .EQ. 0 ) THEN
422                    else              fcrmem      ( ichknum ) = fcref
423  c              fcppmem     ( ichknum ) = fcpertplus
424                       print *, 'ph-grd 3 -------------------------------'              fcpmmem     ( ichknum ) = fcpertminus
425                       print *, 'ph-grd 3 : ierr = ', ierr,              xxmemref    ( ichknum ) = xxmemo_ref
426       &                                 ', icomp = ', icomp              xxmempert   ( ichknum ) = xxmemo_pert
427                    endif              gfdmem      ( ichknum ) = gfd
428                 else              adxxmem     ( ichknum ) = adxxmemo
429                    ierr_grdchk = -1              ftlxxmem    ( ichknum ) = ftlxxmemo
430                 endif              ratioadmem  ( ichknum ) = ratio_ad
431                            ratioftlmem ( ichknum ) = ratio_ftl
432              enddo  
433           endif              irecmem   ( ichknum ) = icvrec
434                bimem     ( ichknum ) = itile
435  c--   Everyone has to wait for the component to be reset.              bjmem     ( ichknum ) = jtile
436           _BARRIER              ilocmem   ( ichknum ) = itilepos
437                jlocmem   ( ichknum ) = jtilepos
438        enddo              klocmem   ( ichknum ) = layer
439                iobcsmem  ( ichknum ) = obcspos
440        CALL WRITE_REC_XYZ_RL( 'grd_findiff'   , tmpplot1, 1, 0, myThid)              icompmem  ( ichknum ) = icomp
441        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)              ichkmem   ( ichknum ) = ichknum
442        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)              itestmem  ( ichknum ) = itest
443                ierrmem   ( ichknum ) = ierr
444                icglomem  ( ichknum ) = icglo
445              ENDIF
446    
447              IF ( myProcId .EQ. grdchkwhichproc .AND. ierr .EQ. 0 ) THEN
448                WRITE(standardmessageunit,'(A)')
449         &          'grad-res -------------------------------'
450                WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
451         &           ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
452         &             layer,itile,jtile,obcspos,
453         &             fcref, fcpertplus, fcpertminus
454    #ifdef ALLOW_TANGENTLINEAR_RUN
455                WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
456         &           ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
457         &             icompmem(ichknum),itestmem(ichknum),
458         &             bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
459         &             ftlxxmemo, gfd, ratio_ftl
460    #else
461                WRITE(standardmessageunit,'(A,8I5,1x,1P3E19.11)')
462         &           ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
463         &             icompmem(ichknum),itestmem(ichknum),
464         &             bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
465         &             adxxmemo, gfd, ratio_ad
466    #endif
467              ENDIF
468    #ifdef ALLOW_TANGENTLINEAR_RUN
469              WRITE(msgBuf,'(A30,1PE22.14)')
470         &              ' TLM  ref_cost_function      =', fcref
471              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
472              WRITE(msgBuf,'(A30,1PE22.14)')
473         &              ' TLM  tangent-lin_grad       =', ftlxxmemo
474              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
475              WRITE(msgBuf,'(A30,1PE22.14)')
476         &              ' TLM  finite-diff_grad       =', gfd
477              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
478    #else
479              WRITE(msgBuf,'(A30,1PE22.14)')
480         &              ' ADM  ref_cost_function      =', fcref
481              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
482              WRITE(msgBuf,'(A30,1PE22.14)')
483         &              ' ADM  adjoint_gradient       =', adxxmemo
484              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
485              WRITE(msgBuf,'(A30,1PE22.14)')
486         &              ' ADM  finite-diff_grad       =', gfd
487              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
488    #endif
489    
490              WRITE(msgBuf,'(A,I4,A,I3,A)')
491         &      '====== End of gradient-check number', ichknum,
492         &                      ' (ierr=', ierr, ') ======='
493              CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
494    
495    C-- else of if ichknum < ...
496            ELSE
497              ierr_grdchk = -1
498    
499    C-- end of if ichknum < ...
500            ENDIF
501    
502    C-- end of do icomp = ...
503          ENDDO
504    
505          IF (myProcId .EQ. grdchkwhichproc .AND. .NOT.useSingleCpuIO) THEN
506            CALL WRITE_REC_XYZ_RL(
507         &       'grd_findiff'   , tmpplot1, 1, 0, myThid)
508            CALL WRITE_REC_XYZ_RL(
509         &       'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
510            CALL WRITE_REC_XYZ_RL(
511         &       'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
512          ENDIF
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_GRADIENT_CHECK */  #endif /* ALLOW_GRDCHK */
518    
519        end        RETURN
520          END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.41
  ViewVC Help
Powered by ViewVC 1.1.22