/[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.5 by heimbach, Sat Jul 13 02:55:58 2002 UTC revision 1.6 by heimbach, Mon Sep 16 18:11:58 2002 UTC
# Line 106  c     == local variables == Line 106  c     == local variables ==
106        _RL     gfd        _RL     gfd
107        _RL     fcref        _RL     fcref
108        _RL     fcpertplus, fcpertminus        _RL     fcpertplus, fcpertminus
109        _RL     ratio        _RL     ratio_ad
110          _RL     ratio_ftl
111        _RL     xxmemo_ref        _RL     xxmemo_ref
112        _RL     xxmemo_pert        _RL     xxmemo_pert
113        _RL     adxxmemo        _RL     adxxmemo
114          _RL     ftlxxmemo
115          _RL     localEps
116          _RL     grdchk_epsfac
117    
118    #ifdef ALLOW_TANGENTLINEAR_RUN
119          _RL     g_fc
120          common /g_cost_r/ g_fc
121    #endif
122    
 cph(  
123        _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)        _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
124        _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)
125  cph)        _RL tmpplot3(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
   
 Cml(  
       _RL     grdchk_epsfac  
 Cml)  
126    
127  c     == end of interface ==  c     == end of interface ==
128  CEOP  CEOP
# Line 138  c--   initialise variables Line 142  c--   initialise variables
142    
143  c--   Compute the adjoint models' gradients.  c--   Compute the adjoint models' gradients.
144  c--   Compute the unperturbed cost function.  c--   Compute the unperturbed cost function.
145  cph   Gradient via adjoint has already been computed,  c--   Gradient via adjoint has already been computed,
146  cph   and so has unperturbed cost function,  c--   and so has unperturbed cost function,
147  cph   assuming all xx_ fields are initialised to zero.  c--   assuming all xx_ fields are initialised to zero.
148    
149        fcref = fc        fcref = fc
150        ierr_grdchk = 0        ierr_grdchk = 0
151    
 cph(  
152        do bj = jtlo, jthi        do bj = jtlo, jthi
153           do bi = itlo, ithi           do bi = itlo, ithi
154              do k = 1, nr              do k = 1, nr
# Line 153  cph( Line 156  cph(
156                    do i = imin, imax                    do i = imin, imax
157                       tmpplot1(i,j,k,bi,bj) = 0. _d 0                       tmpplot1(i,j,k,bi,bj) = 0. _d 0
158                       tmpplot2(i,j,k,bi,bj) = 0. _d 0                       tmpplot2(i,j,k,bi,bj) = 0. _d 0
159                         tmpplot3(i,j,k,bi,bj) = 0. _d 0
160                    end do                    end do
161                 end do                 end do
162              end do              end do
163           end do           end do
164        end do        end do
 cph)  
165    
166        if ( useCentralDiff ) then        if ( useCentralDiff ) then
167           grdchk_epsfac = 2. _d 0           grdchk_epsfac = 2. _d 0
# Line 180  c--   gradient checks. Line 183  c--   gradient checks.
183    
184                 if (ichknum .le. maxgrdchecks ) then                 if (ichknum .le. maxgrdchecks ) then
185    
186  c--         Determine the location of icomp on the grid.  c--   Determine the location of icomp on the grid.
187                    call grdchk_loc( icomp, ichknum,                    call grdchk_loc( icomp, ichknum,
188       &                 icvrec, itile, jtile, layer,       &                 icvrec, itile, jtile, layer,
189       &                 itilepos, jtilepos, itest, ierr,       &                 itilepos, jtilepos, itest, ierr,
# Line 188  c--         Determine the location of ic Line 191  c--         Determine the location of ic
191                                
192                    if ( ierr .eq. 0 ) then                    if ( ierr .eq. 0 ) then
193    
194  c--   positive perturbation  c******************************************************
195                       grdchk_eps = abs(grdchk_eps)  c--   (A): get gradient component calculated via adjoint
196    c******************************************************
 c--   get gradient component calculated via adjoint  
197                       call grdchk_getadxx( icvrec,                       call grdchk_getadxx( icvrec,
198       &                    itile, jtile, layer,       &                    itile, jtile, layer,
199       &                    itilepos, jtilepos,       &                    itilepos, jtilepos,
200       &                    adxxmemo, mythid )       &                    adxxmemo, mythid )
201                       _BARRIER                       _BARRIER
202    
203    #ifdef ALLOW_TANGENTLINEAR_RUN
204    c******************************************************
205    c--   (B): Get gradient component g_fc from tangent linear run:
206    c******************************************************
207    c--
208    c--   1. perturb control vector component: xx(i)=1.
209    
210                         localEps = 1.
211                         call grdchk_getxx( icvrec, TANGENT_SIMULATION,
212         &                    itile, jtile, layer,
213         &                    itilepos, jtilepos,
214         &                    xxmemo_ref, xxmemo_pert, localEps,
215         &                    mythid )
216                         _BARRIER
217    
218    c--
219    c--   2. perform tangent linear run
220                         mytime = starttime
221                         myiter = niter0
222                         g_fc = 0.
223                         call g_the_main_loop( mytime, myiter, mythid )
224                         ftlxxmemo = g_fc
225    c--
226    c--   3. reset control vector
227                         call grdchk_setxx( icvrec, TANGENT_SIMULATION,
228         &                    itile, jtile, layer,
229         &                    itilepos, jtilepos,
230         &                    xxmemo_ref, mythid )
231                         _BARRIER
232    #endif
233    
234    c******************************************************
235    c--   (C): Get gradient via finite difference perturbation
236    c******************************************************
237    
238  c--   get control vector component from file  c--   get control vector component from file
239  c--   perturb it and write back to file  c--   perturb it and write back to file:
240                       call grdchk_getxx( icvrec,  c--   positive perturbation
241                         localEps = abs(grdchk_eps)
242                         call grdchk_getxx( icvrec, FORWARD_SIMULATION,
243       &                    itile, jtile, layer,       &                    itile, jtile, layer,
244       &                    itilepos, jtilepos,       &                    itilepos, jtilepos,
245       &                    xxmemo_ref, xxmemo_pert, mythid )       &                    xxmemo_ref, xxmemo_pert, localEps,
246         &                    mythid )
247                       _BARRIER                       _BARRIER
248                                            
249  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
# Line 213  c--   forward run with perturbed control Line 253  c--   forward run with perturbed control
253                       fcpertplus = fc                       fcpertplus = fc
254                                        
255  c--   Reset control vector.  c--   Reset control vector.
256                       call grdchk_setxx( icvrec,                       call grdchk_setxx( icvrec, FORWARD_SIMULATION,
257       &                    itile, jtile, layer,       &                    itile, jtile, layer,
258       &                    itilepos, jtilepos,       &                    itilepos, jtilepos,
259       &                    xxmemo_ref, mythid )       &                    xxmemo_ref, mythid )
# Line 223  c--   Reset control vector. Line 263  c--   Reset control vector.
263    
264                       if ( useCentralDiff ) then                       if ( useCentralDiff ) then
265    
 c--   repeat the proceedure for a negative perturbation  
                         grdchk_eps = - abs(grdchk_eps)  
   
266  c--   get control vector component from file  c--   get control vector component from file
267  c--   perturb it and write back to file  c--   perturb it and write back to file:
268                          call grdchk_getxx( icvrec,  c--   repeat the proceedure for a negative perturbation
269                            localEps = - abs(grdchk_eps)
270                            call grdchk_getxx( icvrec, FORWARD_SIMULATION,
271       &                    itile, jtile, layer,       &                    itile, jtile, layer,
272       &                    itilepos, jtilepos,       &                    itilepos, jtilepos,
273       &                    xxmemo_ref, xxmemo_pert, mythid )       &                    xxmemo_ref, xxmemo_pert, localEps,
274         &                    mythid )
275                          _BARRIER                          _BARRIER
276                                            
277  c--   forward run with perturbed control vector  c--   forward run with perturbed control vector
# Line 241  c--   forward run with perturbed control Line 281  c--   forward run with perturbed control
281                          fcpertminus = fc                          fcpertminus = fc
282                                        
283  c--   Reset control vector.  c--   Reset control vector.
284                          call grdchk_setxx( icvrec,                          call grdchk_setxx( icvrec, FORWARD_SIMULATION,
285       &                    itile, jtile, layer,       &                    itile, jtile, layer,
286       &                    itilepos, jtilepos,       &                    itilepos, jtilepos,
287       &                    xxmemo_ref, mythid )       &                    xxmemo_ref, mythid )
288                          _BARRIER                          _BARRIER
289    
 c     reset grdchk_eps to a positive value  
                         grdchk_eps = abs(grdchk_eps)  
   
290                       end if                       end if
291  c--  c--
292    c******************************************************
293    c--   (D): calculate relative differences between gradients
294    c******************************************************
295    
296                       if ( grdchk_eps .eq. 0. ) then                       if ( grdchk_eps .eq. 0. ) then
297                          gfd = (fcpertplus-fcpertminus)                          gfd = (fcpertplus-fcpertminus)
298                       else                       else
299                          gfd = (fcpertplus-fcpertminus)                          gfd = (fcpertplus-fcpertminus)
300       &                       /(grdchk_epsfac*grdchk_eps)       &                       /(grdchk_epsfac*grdchk_eps)
301                       endif                       endif
302                        
303                       if ( adxxmemo .eq. 0. ) then                       if ( adxxmemo .eq. 0. ) then
304                          ratio = abs( adxxmemo - gfd )                          ratio_ad = abs( adxxmemo - gfd )
305                       else                       else
306                          ratio = 1. - gfd/adxxmemo                          ratio_ad = 1. - gfd/adxxmemo
307                         endif
308                      
309                         if ( ftlxxmemo .eq. 0. ) then
310                            ratio_ftl = abs( ftlxxmemo - gfd )
311                         else
312                            ratio_ftl = 1. - gfd/ftlxxmemo
313                       endif                       endif
314                                        
 cph(  
315                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =                       tmpplot1(itilepos,jtilepos,layer,itile,jtile) =
316       &                    gfd       &                    gfd
317                       tmpplot2(itilepos,jtilepos,layer,itile,jtile) =                       tmpplot2(itilepos,jtilepos,layer,itile,jtile) =
318       &                    ratio       &                    ratio_ad
319  cph)                       tmpplot3(itilepos,jtilepos,layer,itile,jtile) =
320         &                    ratio_ftl
321    
322                                            
323                       fcrmem    ( ichknum ) = fcref                       fcrmem      ( ichknum ) = fcref
324                       fcppmem   ( ichknum ) = fcpertplus                       fcppmem     ( ichknum ) = fcpertplus
325                       fcpmmem   ( ichknum ) = fcpertminus                       fcpmmem     ( ichknum ) = fcpertminus
326                       xxmemref  ( ichknum ) = xxmemo_ref                       xxmemref    ( ichknum ) = xxmemo_ref
327                       xxmempert ( ichknum ) = xxmemo_pert                       xxmempert   ( ichknum ) = xxmemo_pert
328                       gfdmem    ( ichknum ) = gfd                       gfdmem      ( ichknum ) = gfd
329                       adxxmem   ( ichknum ) = adxxmemo                       adxxmem     ( ichknum ) = adxxmemo
330                       ratiomem  ( ichknum ) = ratio                       ftlxxmem    ( ichknum ) = ftlxxmemo
331                         ratioadmem  ( ichknum ) = ratio_ad
332                         ratioftlmem ( ichknum ) = ratio_ftl
333    
334                       irecmem   ( ichknum ) = icvrec                       irecmem   ( ichknum ) = icvrec
335                       bimem     ( ichknum ) = itile                       bimem     ( ichknum ) = itile
# Line 301  cph( Line 350  cph(
350                       print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',                       print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
351       &                    ichknum,ichkmem(ichknum),       &                    ichknum,ichkmem(ichknum),
352       &                    icompmem(ichknum),itestmem(ichknum),       &                    icompmem(ichknum),itestmem(ichknum),
353       &                    adxxmemo, gfd, ratio       &                    adxxmemo, gfd, ratio_ad
354                         print '(a,4I5,3(1x,E15.9))', 'ph-grd 3 ',
355         &                    ichknum,ichkmem(ichknum),
356         &                    icompmem(ichknum),itestmem(ichknum),
357         &                    ftlxxmemo, gfd, ratio_ftl
358  cph)  cph)
359                    else                    else
360  c  c
# Line 321  c--   Everyone has to wait for the compo Line 374  c--   Everyone has to wait for the compo
374    
375        enddo        enddo
376    
377  cph(        CALL WRITE_REC_XYZ_RL( 'grd_findiff'   , tmpplot1, 1, 0, myThid)
378        CALL WRITE_REC_XYZ_RL( 'grd_findiff' , tmpplot1, 1, 0, myThid)        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ad'  , tmpplot2, 1, 0, myThid)
379        CALL WRITE_REC_XYZ_RL( 'grd_ratio' , tmpplot2, 1, 0, myThid)        CALL WRITE_REC_XYZ_RL( 'grd_ratio_ftl' , tmpplot3, 1, 0, myThid)
 cph)  
380    
381  c--   Print the results of the gradient check.  c--   Print the results of the gradient check.
382        call grdchk_print( ichknum, ierr_grdchk, mythid )        call grdchk_print( ichknum, ierr_grdchk, mythid )
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22