/[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.9 by heimbach, Thu Oct 2 21:30:22 2003 UTC revision 1.25 by jmc, Sat Sep 15 17:32:41 2007 UTC
# Line 1  Line 1 
1    C
2  C $Header$  C $Header$
3    C $Name$
4    
5    #include "AD_CONFIG.h"
6  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
7    
8  CBOI  CBOI
# Line 78  c     == global variables == Line 81  c     == global variables ==
81  #include "PARAMS.h"  #include "PARAMS.h"
82  #include "grdchk.h"  #include "grdchk.h"
83  #include "cost.h"  #include "cost.h"
84    #include "ctrl.h"
85  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
86  #include "g_cost.h"  #include "g_cost.h"
87  #endif  #endif
# Line 86  C     !INPUT/OUTPUT PARAMETERS: Line 90  C     !INPUT/OUTPUT PARAMETERS:
90  c     == routine arguments ==  c     == routine arguments ==
91        integer mythid        integer mythid
92    
93  #ifdef ALLOW_GRADIENT_CHECK  #ifdef ALLOW_GRDCHK
94  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
95  c     == local variables ==  c     == local variables ==
96        integer myiter        integer myiter
# Line 106  c     == local variables == Line 110  c     == local variables ==
110        integer obcspos        integer obcspos
111        integer itilepos        integer itilepos
112        integer jtilepos        integer jtilepos
113          integer icglo
114        integer itest        integer itest
115        integer ierr        integer ierr
116        integer ierr_grdchk        integer ierr_grdchk
# Line 125  c     == local variables == Line 130  c     == local variables ==
130        _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
131        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
132    
133          CHARACTER*(MAX_LEN_MBUF) msgBuf
134    
135  c     == end of interface ==  c     == end of interface ==
136  CEOP  CEOP
137    
# Line 138  c--   Set the loop ranges. Line 145  c--   Set the loop ranges.
145        imin = 1        imin = 1
146        imax = snx        imax = snx
147    
148          print *, 'ph-check entering grdchk_main '
149    
150  c--   initialise variables  c--   initialise variables
151        call grdchk_init( mythid )        call grdchk_init( mythid )
152    
# Line 148  cph   and so has unperturbed cost functi Line 157  cph   and so has unperturbed cost functi
157  cph   assuming all xx_ fields are initialised to zero.  cph   assuming all xx_ fields are initialised to zero.
158    
159        ierr_grdchk = 0        ierr_grdchk = 0
160  cphadmtlm(  #ifdef ALLOW_ADMTLM
161          fcref = objf_state_final(idep,jdep,1,1,1)
162    #else
163        fcref = fc        fcref = fc
164  cphadmtlm      fcref = objf_state_final(45,4,1,1)  #endif
 cphadmtlm)  
165    
166        print *, 'ph-check fcref = ', fcref        print *, 'ph-check fcref = ', fcref
167    
# Line 175  cphadmtlm) Line 185  cphadmtlm)
185           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
186        end if        end if
187    
188        print *, 'ph-grd 3 -------------------------------'        WRITE(standardmessageunit,'(A)')
189        print ('(2a)'),       &    'grad-res -------------------------------'
190       &     ' ph-grd 3  proc    #    i    j    k            fc ref',        WRITE(standardmessageunit,'(2a)')
191       &     '        fc + eps        fc - eps'       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
192         &    '               fc ref           fc + eps           fc - eps'
193  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
194        print ('(2a)'),        WRITE(standardmessageunit,'(2a)')
195       &     ' ph-grd 3  proc    #    i    j    k          tlm grad',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
196       &     '         fd grad      1 - fd/tlm'       &    '             tlm grad            fd grad         1 - fd/tlm'
197  #else  #else
198        print ('(2a)'),        WRITE(standardmessageunit,'(2a)')
199       &     ' ph-grd 3  proc    #    i    j    k          adj grad',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
200       &     '         fd grad      1 - fd/adj'       &    '             adj grad            fd grad         1 - fd/adj'
201  #endif  #endif
202    
203  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
204  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
205  c--   gradient checks.  c--   gradient checks.
206    
207          if ( nbeg .EQ. 0 )
208         &     call grdchk_get_position( mythid )
209    
210        do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
211    
212           ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
213    
214    cph(
215    cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ',
216    cph-print     &        nbeg, icomp, ichknum
217    cph)
218           if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
219    
220  c--         Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
221              if ( myProcId .EQ. grdchkwhichproc ) then              if ( myProcId .EQ. grdchkwhichproc ) then
222                 call grdchk_loc( icomp, ichknum,                 call grdchk_loc( icomp, ichknum,
223       &              icvrec, itile, jtile, layer, obcspos,       &              icvrec, itile, jtile, layer, obcspos,
224       &              itilepos, jtilepos, itest, ierr,       &              itilepos, jtilepos, icglo, itest, ierr,
225       &              mythid )       &              mythid )
226    cph(
227    cph-print               print *, 'ph-grd ----- back from loc -----',
228    cph-print     &             icvrec, itilepos, jtilepos, layer, obcspos
229    cph)
230              endif              endif
231              _BARRIER              _BARRIER
232                                
# Line 228  c--   (B): Get gradient component g_fc f Line 250  c--   (B): Get gradient component g_fc f
250  c******************************************************  c******************************************************
251  c--  c--
252  c--   1. perturb control vector component: xx(i)=1.  c--   1. perturb control vector component: xx(i)=1.
253                ftlxxmemo = 0.
254    
255              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
256       &           ierr .EQ. 0 ) then       &           ierr .EQ. 0 ) then
# Line 244  c-- Line 267  c--
267  c--   2. perform tangent linear run  c--   2. perform tangent linear run
268              mytime = starttime              mytime = starttime
269              myiter = niter0              myiter = niter0
270  cphadmtlm(  #ifdef ALLOW_ADMTLM
271                do k=1,4*Nr+1
272                 do j=1,sny
273                  do i=1,snx
274                   g_objf_state_final(i,j,1,1,k) = 0.
275                  enddo
276                 enddo
277                enddo
278    #else
279              g_fc = 0.              g_fc = 0.
280  cphadmtlm            do j=1,sny  #endif
281  cphadmtlm               do i=1,snx  
 cphadmtlm                  g_objf_state_final(i,j,1,1) = 0.  
 cphadmtlm               enddo  
 cphadmtlm            enddo  
 cphadmtlm)  
282              call g_the_main_loop( mytime, myiter, mythid )              call g_the_main_loop( mytime, myiter, mythid )
 cphadmtlm(  
             ftlxxmemo = g_fc  
 cphadmtlm            ftlxxmemo = g_objf_state_final(45,4,1,1)  
 cphadmtlm)  
283              _BARRIER              _BARRIER
284    #ifdef ALLOW_ADMTLM
285                ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
286    #else
287                ftlxxmemo = g_fc
288    #endif
289    
290  c--  c--
291  c--   3. reset control vector  c--   3. reset control vector
292              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
# Line 294  c--   forward run with perturbed control Line 323  c--   forward run with perturbed control
323              mytime = starttime              mytime = starttime
324              myiter = niter0              myiter = niter0
325              call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
326  cphadmtlm(  #ifdef ALLOW_ADMTLM
327                fcpertplus = objf_state_final(idep,jdep,1,1,1)
328    #else
329              fcpertplus = fc              fcpertplus = fc
330  cphadmtlm            fcpertplus = objf_state_final(45,4,1,1)  #endif
331  cphadmtlm)              print *, 'ph-check fcpertplus = ', fcpertplus
332              _BARRIER              _BARRIER
333                                        
334  c--   Reset control vector.  c--   Reset control vector.
# Line 311  c--   Reset control vector. Line 342  c--   Reset control vector.
342              _BARRIER              _BARRIER
343    
344              fcpertminus = fcref              fcpertminus = fcref
345                print *, 'ph-check fcpertminus = ', fcpertminus
346    
347              if ( useCentralDiff ) then              if ( useCentralDiff ) then
348    
# Line 333  c--   forward run with perturbed control Line 365  c--   forward run with perturbed control
365                 myiter = niter0                 myiter = niter0
366                 call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
367                 _BARRIER                 _BARRIER
368    #ifdef ALLOW_ADMTLM
369                   fcpertminus = objf_state_final(idep,jdep,1,1,1)
370    #else
371                 fcpertminus = fc                 fcpertminus = fc
372    #endif
373                                        
374  c--   Reset control vector.  c--   Reset control vector.
375                 if ( myProcId .EQ. grdchkwhichproc .AND.                 if ( myProcId .EQ. grdchkwhichproc .AND.
# Line 403  c*************************************** Line 439  c***************************************
439                 ichkmem   ( ichknum ) = ichknum                 ichkmem   ( ichknum ) = ichknum
440                 itestmem  ( ichknum ) = itest                 itestmem  ( ichknum ) = itest
441                 ierrmem   ( ichknum ) = ierr                 ierrmem   ( ichknum ) = ierr
442                                     icglomem  ( ichknum ) = icglo
443                 print *, 'ph-grd 3 -------------------------------'  
444                 print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',                 WRITE(standardmessageunit,'(A)')
445       &              myprocid,ichknum,itilepos,jtilepos,layer,       &             'grad-res -------------------------------'
446                   WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
447         &              ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
448         &              layer,itile,jtile,obcspos,
449       &              fcref, fcpertplus, fcpertminus       &              fcref, fcpertplus, fcpertminus
450  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
451                 print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
452       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
453       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
454         &              bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
455       &              ftlxxmemo, gfd, ratio_ftl       &              ftlxxmemo, gfd, ratio_ftl
456                   WRITE(msgBuf,'(A34,2(1PE24.14,X))')
457         &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
458                   CALL PRINT_MESSAGE
459         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
460    c        
461                   WRITE(msgBuf,'(A34,1PE24.14)')
462         &              ' TLM  precision_derivative_cost =', fcref
463                   CALL PRINT_MESSAGE
464         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
465                   WRITE(msgBuf,'(A34,1PE24.14)')
466         &              ' TLM  precision_derivative_grad =', ftlxxmemo
467                   CALL PRINT_MESSAGE
468         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
469  #else  #else
470                 print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
471       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
472       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
473         &              bimem(ichknum),bjmem(ichknum),obcspos,
474       &              adxxmemo, gfd, ratio_ad       &              adxxmemo, gfd, ratio_ad
475                   WRITE(msgBuf,'(A34,2(1PE24.14,X))')
476         &              'precision_grdchk_result ADM ', fcref, adxxmemo
477                   CALL PRINT_MESSAGE
478         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
479    c
480                   WRITE(msgBuf,'(A34,1PE24.14)')
481         &              ' ADM  precision_derivative_cost =', fcref
482                   CALL PRINT_MESSAGE
483         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
484                   WRITE(msgBuf,'(A34,1PE24.14)')
485         &              ' ADM  precision_derivative_grad =', adxxmemo
486                   CALL PRINT_MESSAGE
487         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
488  #endif  #endif
489    
490              endif              endif
# Line 453  c--   Everyone has to wait for the compo Line 520  c--   Everyone has to wait for the compo
520  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
521        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
522    
523  #endif /* ALLOW_GRADIENT_CHECK */  #endif /* ALLOW_GRDCHK */
524    
525        end        end
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22