/[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.15 by heimbach, Tue Oct 26 20:10:25 2004 UTC revision 1.28 by heimbach, Fri May 29 06:14:47 2009 UTC
# Line 81  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 109  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 143  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 152  cph   Gradient via adjoint has already b Line 156  cph   Gradient via adjoint has already b
156  cph   and so has unperturbed cost function,  cph   and so has unperturbed cost function,
157  cph   assuming all xx_ fields are initialised to zero.  cph   assuming all xx_ fields are initialised to zero.
158    
159          ierr      = 0
160        ierr_grdchk = 0        ierr_grdchk = 0
161  cphadmtlm(        adxxmemo  = 0.
162          ftlxxmemo = 0.
163    #ifdef ALLOW_ADMTLM
164          fcref = objf_state_final(idep,jdep,1,1,1)
165    #else
166        fcref = fc        fcref = fc
167  cphadmtlm      fcref = objf_state_final(45,4,1,1,1)  #endif
 cphadmtlm)  
168    
169        print *, 'ph-check fcref = ', fcref        print *, 'ph-check fcref = ', fcref
170    
# Line 180  cphadmtlm) Line 188  cphadmtlm)
188           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
189        end if        end if
190    
191        print *, 'grad-res -------------------------------'        WRITE(standardmessageunit,'(A)')
192        print ('(2a)'),       &    'grad-res -------------------------------'
193       &     ' grad-res  proc    #    i    j    k            fc ref',        WRITE(standardmessageunit,'(2a)')
194       &     '        fc + eps        fc - eps'       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
195         &    '               fc ref           fc + eps           fc - eps'
196  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
197        print ('(2a)'),        WRITE(standardmessageunit,'(2a)')
198       &     ' grad-res  proc    #    i    j    k          tlm grad',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
199       &     '         fd grad      1 - fd/tlm'       &    '             tlm grad            fd grad         1 - fd/tlm'
200  #else  #else
201        print ('(2a)'),        WRITE(standardmessageunit,'(2a)')
202       &     ' grad-res  proc    #    i    j    k          adj grad',       &    ' grad-res  proc    #    i    j    k   bi   bj iobc',
203       &     '         fd grad      1 - fd/adj'       &    '             adj grad            fd grad         1 - fd/adj'
204  #endif  #endif
205    
206  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
207  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
208  c--   gradient checks.  c--   gradient checks.
209    
210        if ( nbeg .EQ. 0 ) call grdchk_get_position( mythid )        if ( nbeg .EQ. 0 )
211         &     call grdchk_get_position( mythid )
212    
213        do icomp = nbeg, nend, nstep        do icomp = nbeg, nend, nstep
214    
215           ichknum = (icomp - nbeg)/nstep + 1           ichknum = (icomp - nbeg)/nstep + 1
216    
217    cph(
218    cph-print         print *, 'ph-grd _main: nbeg, icomp, ichknum ',
219    cph-print     &        nbeg, icomp, ichknum
220    cph)
221           if (ichknum .le. maxgrdchecks ) then           if (ichknum .le. maxgrdchecks ) then
222    
223  c--         Determine the location of icomp on the grid.  c--         Determine the location of icomp on the grid.
224              if ( myProcId .EQ. grdchkwhichproc ) then              if ( myProcId .EQ. grdchkwhichproc ) then
225                 call grdchk_loc( icomp, ichknum,                 call grdchk_loc( icomp, ichknum,
226       &              icvrec, itile, jtile, layer, obcspos,       &              icvrec, itile, jtile, layer, obcspos,
227       &              itilepos, jtilepos, itest, ierr,       &              itilepos, jtilepos, icglo, itest, ierr,
228       &              mythid )       &              mythid )
229    cph(
230    cph-print               print *, 'ph-grd ----- back from loc -----',
231    cph-print     &             icvrec, itilepos, jtilepos, layer, obcspos
232    cph)
233              endif              endif
234              _BARRIER              _BARRIER
235                                
# Line 251  c-- Line 269  c--
269  c--   2. perform tangent linear run  c--   2. perform tangent linear run
270              mytime = starttime              mytime = starttime
271              myiter = niter0              myiter = niter0
272  cphadmtlm(  #ifdef ALLOW_ADMTLM
273                do k=1,4*Nr+1
274                 do j=1,sny
275                  do i=1,snx
276                   g_objf_state_final(i,j,1,1,k) = 0.
277                  enddo
278                 enddo
279                enddo
280    #else
281              g_fc = 0.              g_fc = 0.
282  cphadmtlm            do j=1,sny  #endif
283  cphadmtlm               do i=1,snx  
 cphadmtlm                  g_objf_state_final(i,j,1,1,1) = 0.  
 cphadmtlm                  g_objf_state_final(i,j,1,1,2) = 0.  
 cphadmtlm               enddo  
 cphadmtlm            enddo  
 cphadmtlm)  
284              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,1)  
 cphadmtlm)  
285              _BARRIER              _BARRIER
286    #ifdef ALLOW_ADMTLM
287                ftlxxmemo = g_objf_state_final(idep,jdep,1,1,1)
288    #else
289                ftlxxmemo = g_fc
290    #endif
291    
292  c--  c--
293  c--   3. reset control vector  c--   3. reset control vector
294              if ( myProcId .EQ. grdchkwhichproc .AND.              if ( myProcId .EQ. grdchkwhichproc .AND.
# Line 302  c--   forward run with perturbed control Line 325  c--   forward run with perturbed control
325              mytime = starttime              mytime = starttime
326              myiter = niter0              myiter = niter0
327              call the_main_loop( mytime, myiter, mythid )              call the_main_loop( mytime, myiter, mythid )
328  cphadmtlm(  #ifdef ALLOW_ADMTLM
329                fcpertplus = objf_state_final(idep,jdep,1,1,1)
330    #else
331              fcpertplus = fc              fcpertplus = fc
332  cphadmtlm            fcpertplus = objf_state_final(45,4,1,1,1)  #endif
 cphadmtlm)  
333              print *, 'ph-check fcpertplus = ', fcpertplus              print *, 'ph-check fcpertplus = ', fcpertplus
334              _BARRIER              _BARRIER
335                                        
# Line 343  c--   forward run with perturbed control Line 367  c--   forward run with perturbed control
367                 myiter = niter0                 myiter = niter0
368                 call the_main_loop( mytime, myiter, mythid )                 call the_main_loop( mytime, myiter, mythid )
369                 _BARRIER                 _BARRIER
370    #ifdef ALLOW_ADMTLM
371                   fcpertminus = objf_state_final(idep,jdep,1,1,1)
372    #else
373                 fcpertminus = fc                 fcpertminus = fc
374    #endif
375                                        
376  c--   Reset control vector.  c--   Reset control vector.
377                 if ( myProcId .EQ. grdchkwhichproc .AND.                 if ( myProcId .EQ. grdchkwhichproc .AND.
# Line 413  c*************************************** Line 441  c***************************************
441                 ichkmem   ( ichknum ) = ichknum                 ichkmem   ( ichknum ) = ichknum
442                 itestmem  ( ichknum ) = itest                 itestmem  ( ichknum ) = itest
443                 ierrmem   ( ichknum ) = ierr                 ierrmem   ( ichknum ) = ierr
444                                     icglomem  ( ichknum ) = icglo
445                 print *, 'grad-res -------------------------------'  
446                 print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',                 WRITE(standardmessageunit,'(A)')
447       &              myprocid,ichknum,itilepos,jtilepos,layer,       &             'grad-res -------------------------------'
448                   WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
449         &              ' grad-res ',myprocid,ichknum,itilepos,jtilepos,
450         &              layer,itile,jtile,obcspos,
451       &              fcref, fcpertplus, fcpertminus       &              fcref, fcpertplus, fcpertminus
452  #ifdef ALLOW_TANGENTLINEAR_RUN  #ifdef ALLOW_TANGENTLINEAR_RUN
453                 print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
454       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
455       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
456         &              bimem(ichknum),bjmem(ichknum),iobcsmem(ichknum),
457       &              ftlxxmemo, gfd, ratio_ftl       &              ftlxxmemo, gfd, ratio_ftl
458                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')
459       &              'precision_grdchk_result TLM ', fcref, ftlxxmemo       &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
460                 CALL PRINT_MESSAGE                 CALL PRINT_MESSAGE
461       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
462    c        
463                   WRITE(msgBuf,'(A34,1PE24.14)')
464         &              ' TLM  precision_derivative_cost =', fcref
465                   CALL PRINT_MESSAGE
466         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
467                   WRITE(msgBuf,'(A34,1PE24.14)')
468         &              ' TLM  precision_derivative_grad =', ftlxxmemo
469                   CALL PRINT_MESSAGE
470         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
471  #else  #else
472                 print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
473       &              myprocid,ichknum,ichkmem(ichknum),       &              ' grad-res ',myprocid,ichknum,ichkmem(ichknum),
474       &              icompmem(ichknum),itestmem(ichknum),       &              icompmem(ichknum),itestmem(ichknum),
475         &              bimem(ichknum),bjmem(ichknum),obcspos,
476       &              adxxmemo, gfd, ratio_ad       &              adxxmemo, gfd, ratio_ad
477                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')  c              WRITE(msgBuf,'(A34,2(1PE24.14,X))')
478       &              'precision_grdchk_result ADM ', fcref, adxxmemo  c    &              'precision_grdchk_result ADM ', fcref, adxxmemo
479    c              CALL PRINT_MESSAGE
480    c    &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
481                   WRITE(msgBuf,'(A34,1PE24.14)')
482         &              ' ADM  precision_derivative_cost =', fcref
483                   CALL PRINT_MESSAGE
484         &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
485                   WRITE(msgBuf,'(A34,1PE24.14)')
486         &              ' ADM  precision_derivative_grad =', adxxmemo
487                 CALL PRINT_MESSAGE                 CALL PRINT_MESSAGE
488       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
489  #endif  #endif
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22