/[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.3 by heimbach, Fri Sep 28 16:53:46 2001 UTC revision 1.3.6.2 by heimbach, Thu Feb 13 23:36:18 2003 UTC
# Line 54  c     ================================== Line 54  c     ==================================
54  c     o Compare the gradients calculated by the adjoint model with  c     o Compare the gradients calculated by the adjoint model with
55  c       finite difference approximations.  c       finite difference approximations.
56  c     started: Christian Eckert eckert@mit.edu 24-Feb-2000  c     started: Christian Eckert eckert@mit.edu 24-Feb-2000
57  c     continued: heimbach@mit.edu: 13-Jun-2001  c     continued&finished: heimbach@mit.edu: 13-Jun-2001
58    c     changed: mlosch@ocean.mit.edu: 09-May-2002
59    c              - added centered difference vs. 1-sided difference option
60    c              - improved output format for readability
61    c              - added control variable hFacC
62    c    
63  c     ==================================================================  c     ==================================================================
64  c     SUBROUTINE grdchk_main  c     SUBROUTINE grdchk_main
65  c     ==================================================================  c     ==================================================================
# Line 100  c     == local variables == Line 105  c     == local variables ==
105        integer ierr_grdchk        integer ierr_grdchk
106        _RL     gfd        _RL     gfd
107        _RL     fcref        _RL     fcref
108        _RL     fcpert        _RL     fcpertplus, fcpertminus
109        _RL     ratio        _RL     ratio
110        _RL     xxmemo_ref        _RL     xxmemo_ref
111        _RL     xxmemo_pert        _RL     xxmemo_pert
# Line 111  cph( Line 116  cph(
116        _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)
117  cph)  cph)
118    
119    Cml(
120          _RL     grdchk_epsfac
121    Cml)
122    
123  c     == end of interface ==  c     == end of interface ==
124  CEOP  CEOP
125    
# Line 151  cph( Line 160  cph(
160        end do        end do
161  cph)  cph)
162    
163          if ( useCentralDiff ) then
164             grdchk_epsfac = 2. _d 0
165          else
166             grdchk_epsfac = 1. _d 0
167          end if
168    
169          print *, 'ph-grd 3 -------------------------------'
170          print ('(2a)'),
171         &     ' ph-grd 3  proc    #    i    j    k            fc ref',
172         &     '        fc + eps        fc - eps'
173          print ('(2a)'),
174         &     ' ph-grd 3  proc    #    i    j    k          adj grad',
175         &     '         fd grad      1 - fd/adj'
176    
177  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
178  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
179  c--   gradient checks.  c--   gradient checks.
# Line 173  c--         Determine the location of ic Line 196  c--         Determine the location of ic
196                                
197                    if ( ierr .eq. 0 ) then                    if ( ierr .eq. 0 ) then
198    
199    c--   positive perturbation
200                         grdchk_eps = abs(grdchk_eps)
201    
202    c--   get gradient component calculated via adjoint
203                         call grdchk_getadxx( icvrec,
204         &                    itile, jtile, layer,
205         &                    itilepos, jtilepos,
206         &                    adxxmemo, mythid )
207                         _BARRIER
208    
209  c--   get control vector component from file  c--   get control vector component from file
210  c--   perturb it and write back to file  c--   perturb it and write back to file
211                       call grdchk_getxx( icvrec,                       call grdchk_getxx( icvrec,
# Line 181  c--   perturb it and write back to file Line 214  c--   perturb it and write back to file
214       &                    xxmemo_ref, xxmemo_pert, mythid )       &                    xxmemo_ref, xxmemo_pert, mythid )
215                       _BARRIER                       _BARRIER
216                                            
217  c--   get gradient component calculated via adjoint  c--   forward run with perturbed control vector
218                       call grdchk_getadxx( icvrec,                       mytime = starttime
219                         myiter = niter0
220                         call the_main_loop( mytime, myiter, mythid )
221                         fcpertplus = fc
222                      
223    c--   Reset control vector.
224                         call grdchk_setxx( icvrec,
225       &                    itile, jtile, layer,       &                    itile, jtile, layer,
226       &                    itilepos, jtilepos,       &                    itilepos, jtilepos,
227       &                    adxxmemo, mythid )       &                    xxmemo_ref, mythid )
228                       _BARRIER                       _BARRIER
229    
230                         fcpertminus = fcref
231    
232                         if ( useCentralDiff ) then
233    
234    c--   repeat the proceedure for a negative perturbation
235                            grdchk_eps = - abs(grdchk_eps)
236    
237    c--   get control vector component from file
238    c--   perturb it and write back to file
239                            call grdchk_getxx( icvrec,
240         &                    itile, jtile, layer,
241         &                    itilepos, jtilepos,
242         &                    xxmemo_ref, xxmemo_pert, mythid )
243                            _BARRIER
244                        
245  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
246                       mytime = starttime                          mytime = starttime
247                       myiter = niter0                          myiter = niter0
248                       call the_main_loop( mytime, myiter, mythid )                          call the_main_loop( mytime, myiter, mythid )
249                       fcpert = fc                          fcpertminus = fc
250                                        
251    c--   Reset control vector.
252                            call grdchk_setxx( icvrec,
253         &                    itile, jtile, layer,
254         &                    itilepos, jtilepos,
255         &                    xxmemo_ref, mythid )
256                            _BARRIER
257    
258    c     reset grdchk_eps to a positive value
259                            grdchk_eps = abs(grdchk_eps)
260    
261                         end if
262    c--
263                       if ( grdchk_eps .eq. 0. ) then                       if ( grdchk_eps .eq. 0. ) then
264                          gfd = (fcpert-fcref)                          gfd = (fcpertplus-fcpertminus)
265                       else                       else
266                          gfd = (fcpert-fcref)/grdchk_eps                          gfd = (fcpertplus-fcpertminus)
267         &                       /(grdchk_epsfac*grdchk_eps)
268                       endif                       endif
269                                            
270                       if ( gfd .eq. 0. ) then                       if ( adxxmemo .eq. 0. ) then
271                          ratio = abs( adxxmemo - gfd )                          ratio = abs( adxxmemo - gfd )
272                       else                       else
273                          ratio = 1. - adxxmemo/gfd                          ratio = 1. - gfd/adxxmemo
274                       endif                       endif
275                                        
 c--   Reset control vector.  
                      call grdchk_setxx( icvrec,  
      &                    itile, jtile, layer,  
      &                    itilepos, jtilepos,  
      &                    xxmemo_ref, mythid )  
                      _BARRIER  
   
276  cph(  cph(
277                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =
278       &                    gfd       &                    gfd
# Line 221  cph( Line 281  cph(
281  cph)  cph)
282    
283                                            
284                       fcpmem    ( ichknum ) = fcpert                       fcrmem    ( ichknum ) = fcref
285                         fcppmem   ( ichknum ) = fcpertplus
286                         fcpmmem   ( ichknum ) = fcpertminus
287                       xxmemref  ( ichknum ) = xxmemo_ref                       xxmemref  ( ichknum ) = xxmemo_ref
288                       xxmempert ( ichknum ) = xxmemo_pert                       xxmempert ( ichknum ) = xxmemo_pert
289                       gfdmem    ( ichknum ) = gfd                       gfdmem    ( ichknum ) = gfd
# Line 241  cph) Line 303  cph)
303                                        
304  cph(  cph(
305                       print *, 'ph-grd 3 -------------------------------'                       print *, 'ph-grd 3 -------------------------------'
306                       print *, 'ph-grd 3 ',                       print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
307       &                    ichknum,itilepos,jtilepos,layer,       &                    myprocid,ichknum,itilepos,jtilepos,layer,
308       &                    fcref, fcpert       &                    fcref, fcpertplus, fcpertminus
309                       print *, 'ph-grd 3 ',                       print '(a,5I5,2x,3(1x,E15.9))', ' ph-grd 3 ',
310       &                    ichknum,ichkmem(ichknum),       &                    myprocid,ichknum,ichkmem(ichknum),
311       &                    icompmem(ichknum),itestmem(ichknum),       &                    icompmem(ichknum),itestmem(ichknum),
312       &                    adxxmemo, gfd, ratio       &                    adxxmemo, gfd, ratio
313  cph)  cph)
314                    else                    else
315  c  c
316                         print *, 'ph-grd 3 -------------------------------'
317                         print *, 'ph-grd 3 : ierr = ', ierr,
318         &                                 ', icomp = ', icomp
319                    endif                    endif
320                 else                 else
321                    ierr_grdchk = -1                    ierr_grdchk = -1
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.3.6.2
  ViewVC Help
Powered by ViewVC 1.1.22