/[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.1 by heimbach, Fri Jul 13 13:08:17 2001 UTC revision 1.2 by heimbach, Fri Jul 13 14:50:46 2001 UTC
# Line 0  Line 1 
1    C $Header$
2    
3    #include "CPP_OPTIONS.h"
4    
5    
6          subroutine grdchk_main(
7         I                        mythid
8         &                      )
9    
10    c     ==================================================================
11    c     SUBROUTINE grdchk_main
12    c     ==================================================================
13    c
14    c     o Compare the gradients calculated by the adjoint model with
15    c       finite difference approximations.
16    c
17    c     started: Christian Eckert eckert@mit.edu 24-Feb-2000
18    c     continued: heimbach@mit.edu: 13-Jun-2001
19    c
20    c     ==================================================================
21    c     SUBROUTINE grdchk_main
22    c     ==================================================================
23    
24          implicit none
25    
26    c     == global variables ==
27    
28    #include "SIZE.h"
29    #include "EEPARAMS.h"
30    #include "PARAMS.h"
31    
32    #include "grdchk.h"
33    #include "cost.h"
34    
35    c     == routine arguments ==
36    
37          integer mythid
38    
39    #ifdef ALLOW_GRADIENT_CHECK
40    c     == local variables ==
41          integer myiter
42          _RL     mytime
43          integer bi, itlo, ithi
44          integer bj, jtlo, jthi
45          integer i,  imin, imax
46          integer j,  jmin, jmax
47          integer k
48    
49          integer jprocs
50          integer proc_grdchk
51          integer icomp
52          integer ichknum
53          integer icvrec
54          integer jtile
55          integer itile
56          integer layer
57          integer itilepos
58          integer jtilepos
59          integer itest
60          integer ierr
61          integer ierr_grdchk
62          _RL     gfd
63          _RL     fcref
64          _RL     fcpert
65          _RL     ratio
66          _RL     xxmemo_ref
67          _RL     xxmemo_pert
68          _RL     adxxmemo
69    
70    cph(
71          _RL tmpplot1(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
72          _RL tmpplot2(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
73    cph)
74    
75    c     == end of interface ==
76    
77    c--   Set the loop ranges.
78          jtlo = mybylo(mythid)
79          jthi = mybyhi(mythid)
80          itlo = mybxlo(mythid)
81          ithi = mybxhi(mythid)
82          jmin = 1
83          jmax = sny
84          imin = 1
85          imax = snx
86    
87    c--   initialise variables
88          call grdchk_init( mythid )
89    
90    c--   Compute the adjoint models' gradients.
91    c--   Compute the unperturbed cost function.
92    cph   Gradient via adjoint has already been computed,
93    cph   and so has unperturbed cost function,
94    cph   assuming all xx_ fields are initialised to zero.
95    
96          fcref = fc
97          ierr_grdchk = 0
98    
99    cph(
100          do bj = jtlo, jthi
101             do bi = itlo, ithi
102                do k = 1, nr
103                   do j = jmin, jmax
104                      do i = imin, imax
105                         tmpplot1(i,j,k,bi,bj) = 0. _d 0
106                         tmpplot2(i,j,k,bi,bj) = 0. _d 0
107                      end do
108                   end do
109                end do
110             end do
111          end do
112    cph)
113    
114    c--   Compute the finite difference approximations.
115    c--   Cycle through all processes doing NINT(nend-nbeg+1)/nstep
116    c--   gradient checks.
117          do jprocs = 1,numberOfProcs
118             proc_grdchk = jprocs - 1
119    
120             if ( myProcId .eq. proc_grdchk ) then
121    
122                do icomp = nbeg, nend, nstep
123    
124                   ichknum = (icomp - nbeg)/nstep + 1
125    
126                   if (ichknum .le. maxgrdchecks ) then
127    
128    c--         Determine the location of icomp on the grid.
129                      call grdchk_loc( icomp, ichknum,
130         &                 icvrec, itile, jtile, layer,
131         &                 itilepos, jtilepos, itest, ierr,
132         &                 mythid )
133                  
134                      if ( ierr .eq. 0 ) then
135    
136    c--   get control vector component from file
137    c--   perturb it and write back to file
138                         call grdchk_getxx( icvrec,
139         &                    itile, jtile, layer,
140         &                    itilepos, jtilepos,
141         &                    xxmemo_ref, xxmemo_pert, mythid )
142                         _BARRIER
143                        
144    c--   get gradient component calculated via adjoint
145                         call grdchk_getadxx( icvrec,
146         &                    itile, jtile, layer,
147         &                    itilepos, jtilepos,
148         &                    adxxmemo, mythid )
149                         _BARRIER
150    
151    c--   forward run with perturbed control vector
152                         mytime = starttime
153                         myiter = niter0
154                         call the_main_loop( mytime, myiter, mythid )
155                         fcpert = fc
156                      
157                         if ( grdchk_eps .eq. 0. ) then
158                            gfd = (fcpert-fcref)
159                         else
160                            gfd = (fcpert-fcref)/grdchk_eps
161                         endif
162                        
163                         if ( gfd .eq. 0. ) then
164                            ratio = abs( adxxmemo - gfd )
165                         else
166                            ratio = 1. - adxxmemo/gfd
167                         endif
168                      
169    c--   Reset control vector.
170                         call grdchk_setxx( icvrec,
171         &                    itile, jtile, layer,
172         &                    itilepos, jtilepos,
173         &                    xxmemo_ref, mythid )
174                         _BARRIER
175    
176    cph(
177                         tmpplot1(itilepos,jtilepos,layer,itile,jtile) =
178         &                    gfd
179                         tmpplot2(itilepos,jtilepos,layer,itile,jtile) =
180         &                    ratio
181    cph)
182    
183                        
184                         fcpmem    ( ichknum ) = fcpert
185                         xxmemref  ( ichknum ) = xxmemo_ref
186                         xxmempert ( ichknum ) = xxmemo_pert
187                         gfdmem    ( ichknum ) = gfd
188                         adxxmem   ( ichknum ) = adxxmemo
189                         ratiomem  ( ichknum ) = ratio
190    
191                         irecmem   ( ichknum ) = icvrec
192                         bimem     ( ichknum ) = itile
193                         bjmem     ( ichknum ) = jtile
194                         ilocmem   ( ichknum ) = itilepos
195                         jlocmem   ( ichknum ) = jtilepos
196                         klocmem   ( ichknum ) = layer
197                         icompmem  ( ichknum ) = icomp
198                         ichkmem   ( ichknum ) = ichknum
199                         itestmem  ( ichknum ) = itest
200                         ierrmem   ( ichknum ) = ierr
201                      
202    cph(
203                         print *, 'ph-grd 3 -------------------------------'
204                         print *, 'ph-grd 3 ',
205         &                    ichknum,itilepos,jtilepos,layer,
206         &                    fcref, fcpert
207                         print *, 'ph-grd 3 ',
208         &                    ichknum,ichkmem(ichknum),
209         &                    icompmem(ichknum),itestmem(ichknum),
210         &                    adxxmemo, gfd, ratio
211    cph)
212                      else
213    c
214                      endif
215                   else
216                      ierr_grdchk = -1
217                   endif
218                
219                enddo
220             endif
221    
222    c--   Everyone has to wait for the component to be reset.
223             _BARRIER
224    
225          enddo
226    
227    cph(
228          CALL WRITE_REC_XYZ_RL( 'grd_findiff' , tmpplot1, 1, 0, myThid)
229          CALL WRITE_REC_XYZ_RL( 'grd_ratio' , tmpplot2, 1, 0, myThid)
230    cph)
231    
232    c--   Print the results of the gradient check.
233          call grdchk_print( ichknum, ierr_grdchk, mythid )
234    
235    #endif /* ALLOW_GRADIENT_CHECK */
236    
237          end
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22