/[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.4 by heimbach, Thu May 30 22:47:26 2002 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  c--   Compute the finite difference approximations.  c--   Compute the finite difference approximations.
170  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep  c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
171  c--   gradient checks.  c--   gradient checks.
# Line 173  c--         Determine the location of ic Line 188  c--         Determine the location of ic
188                                
189                    if ( ierr .eq. 0 ) then                    if ( ierr .eq. 0 ) then
190    
191    c--   positive perturbation
192                         grdchk_eps = abs(grdchk_eps)
193    
194    c--   get gradient component calculated via adjoint
195                         call grdchk_getadxx( icvrec,
196         &                    itile, jtile, layer,
197         &                    itilepos, jtilepos,
198         &                    adxxmemo, mythid )
199                         _BARRIER
200    
201  c--   get control vector component from file  c--   get control vector component from file
202  c--   perturb it and write back to file  c--   perturb it and write back to file
203                       call grdchk_getxx( icvrec,                       call grdchk_getxx( icvrec,
# Line 181  c--   perturb it and write back to file Line 206  c--   perturb it and write back to file
206       &                    xxmemo_ref, xxmemo_pert, mythid )       &                    xxmemo_ref, xxmemo_pert, mythid )
207                       _BARRIER                       _BARRIER
208                                            
209  c--   get gradient component calculated via adjoint  c--   forward run with perturbed control vector
210                       call grdchk_getadxx( icvrec,                       mytime = starttime
211                         myiter = niter0
212                         call the_main_loop( mytime, myiter, mythid )
213                         fcpertplus = fc
214                      
215    c--   Reset control vector.
216                         call grdchk_setxx( icvrec,
217       &                    itile, jtile, layer,       &                    itile, jtile, layer,
218       &                    itilepos, jtilepos,       &                    itilepos, jtilepos,
219       &                    adxxmemo, mythid )       &                    xxmemo_ref, mythid )
220                       _BARRIER                       _BARRIER
221    
222                         fcpertminus = fcref
223    
224                         if ( useCentralDiff ) then
225    
226    c--   repeat the proceedure for a negative perturbation
227                            grdchk_eps = - abs(grdchk_eps)
228    
229    c--   get control vector component from file
230    c--   perturb it and write back to file
231                            call grdchk_getxx( icvrec,
232         &                    itile, jtile, layer,
233         &                    itilepos, jtilepos,
234         &                    xxmemo_ref, xxmemo_pert, mythid )
235                            _BARRIER
236                        
237  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
238                       mytime = starttime                          mytime = starttime
239                       myiter = niter0                          myiter = niter0
240                       call the_main_loop( mytime, myiter, mythid )                          call the_main_loop( mytime, myiter, mythid )
241                       fcpert = fc                          fcpertminus = fc
242                                        
243    c--   Reset control vector.
244                            call grdchk_setxx( icvrec,
245         &                    itile, jtile, layer,
246         &                    itilepos, jtilepos,
247         &                    xxmemo_ref, mythid )
248                            _BARRIER
249    
250    c     reset grdchk_eps to a positive value
251                            grdchk_eps = abs(grdchk_eps)
252    
253                         end if
254    c--
255                       if ( grdchk_eps .eq. 0. ) then                       if ( grdchk_eps .eq. 0. ) then
256                          gfd = (fcpert-fcref)                          gfd = (fcpertplus-fcpertminus)
257                       else                       else
258                          gfd = (fcpert-fcref)/grdchk_eps                          gfd = (fcpertplus-fcpertminus)
259         &                       /(grdchk_epsfac*grdchk_eps)
260                       endif                       endif
261                                            
262                       if ( gfd .eq. 0. ) then                       if ( adxxmemo .eq. 0. ) then
263                          ratio = abs( adxxmemo - gfd )                          ratio = abs( adxxmemo - gfd )
264                       else                       else
265                          ratio = 1. - adxxmemo/gfd                          ratio = 1. - gfd/adxxmemo
266                       endif                       endif
267                                        
 c--   Reset control vector.  
                      call grdchk_setxx( icvrec,  
      &                    itile, jtile, layer,  
      &                    itilepos, jtilepos,  
      &                    xxmemo_ref, mythid )  
                      _BARRIER  
   
268  cph(  cph(
269                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =
270       &                    gfd       &                    gfd
# Line 221  cph( Line 273  cph(
273  cph)  cph)
274    
275                                            
276                       fcpmem    ( ichknum ) = fcpert                       fcrmem    ( ichknum ) = fcref
277                         fcppmem   ( ichknum ) = fcpertplus
278                         fcpmmem   ( ichknum ) = fcpertminus
279                       xxmemref  ( ichknum ) = xxmemo_ref                       xxmemref  ( ichknum ) = xxmemo_ref
280                       xxmempert ( ichknum ) = xxmemo_pert                       xxmempert ( ichknum ) = xxmemo_pert
281                       gfdmem    ( ichknum ) = gfd                       gfdmem    ( ichknum ) = gfd
# Line 241  cph) Line 295  cph)
295                                        
296  cph(  cph(
297                       print *, 'ph-grd 3 -------------------------------'                       print *, 'ph-grd 3 -------------------------------'
298                       print *, 'ph-grd 3 ',                       print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
299       &                    ichknum,itilepos,jtilepos,layer,       &                    ichknum,itilepos,jtilepos,layer,
300       &                    fcref, fcpert       &                    fcref, fcpertplus, fcpertminus
301                       print *, 'ph-grd 3 ',                       print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
302       &                    ichknum,ichkmem(ichknum),       &                    ichknum,ichkmem(ichknum),
303       &                    icompmem(ichknum),itestmem(ichknum),       &                    icompmem(ichknum),itestmem(ichknum),
304       &                    adxxmemo, gfd, ratio       &                    adxxmemo, gfd, ratio
305  cph)  cph)
306                    else                    else
307  c  c
308                         print *, 'ph-grd 3 -------------------------------'
309                         print *, 'ph-grd 3 : ierr = ', ierr,
310         &                                 ', icomp = ', icomp
311                    endif                    endif
312                 else                 else
313                    ierr_grdchk = -1                    ierr_grdchk = -1
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22