/[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.16 by heimbach, Thu Jul 28 13:54:36 2005 UTC revision 1.23 by heimbach, Sat Aug 4 18:05:23 2007 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 155  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,1)  #endif
 cphadmtlm)  
165    
166        print *, 'ph-check fcref = ', fcref        print *, 'ph-check fcref = ', fcref
167    
# Line 182  cphadmtlm) Line 185  cphadmtlm)
185           grdchk_epsfac = 1. _d 0           grdchk_epsfac = 1. _d 0
186        end if        end if
187    
188        print *, 'grad-res -------------------------------'        WRITE(standardmessageunit,'(A)')
189        print ('(2a)'),       &    'grad-res -------------------------------'
190       &     ' grad-res  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       &     ' grad-res  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       &     ' grad-res  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 ) call grdchk_get_position( mythid )        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 237  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 253  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,1) = 0.  
 cphadmtlm                  g_objf_state_final(i,j,1,1,2) = 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,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 304  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,1)  #endif
 cphadmtlm)  
331              print *, 'ph-check fcpertplus = ', fcpertplus              print *, 'ph-check fcpertplus = ', fcpertplus
332              _BARRIER              _BARRIER
333                                        
# Line 345  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 415  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 *, 'grad-res -------------------------------'  
444                 print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',                 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))', ' grad-res ',                 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))')                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')
457       &              'precision_grdchk_result TLM ', fcref, ftlxxmemo       &              'precision_grdchk_result TLM ', fcref, ftlxxmemo
458                 CALL PRINT_MESSAGE                 CALL PRINT_MESSAGE
459       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)       &              (msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
460  #else  #else
461                 print '(a,5I5,2x,3(1x,E15.9))', ' grad-res ',                 WRITE(standardmessageunit,'(a,8I5,2x,3(1x,E18.12))')
462       &              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,
465       &              adxxmemo, gfd, ratio_ad       &              adxxmemo, gfd, ratio_ad
466                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')                 WRITE(msgBuf,'(A34,2(1PE24.14,X))')
467       &              'precision_grdchk_result ADM ', fcref, adxxmemo       &              'precision_grdchk_result ADM ', fcref, adxxmemo
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22