/[MITgcm]/MITgcm/optim/optim_readdata.F
ViewVC logotype

Diff of /MITgcm/optim/optim_readdata.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by heimbach, Tue Feb 5 20:34:35 2002 UTC revision 1.1.2.4 by heimbach, Fri Mar 7 05:21:33 2003 UTC
# Line 0  Line 1 
1    
2          subroutine optim_readdata(
3         I                      nn,
4         I                      dfile,
5         I                      lheaderonly,
6         O                      ff,
7         O                      vv
8         &                    )
9    
10    c     ==================================================================
11    c     SUBROUTINE optim_readdata
12    c     ==================================================================
13    c
14    c     o Read the data written by the MITgcmUV state estimation setup and
15    c       join them to one vector that is subsequently used by the minimi-
16    c       zation algorithm "lsopt". Depending on the specified file name
17    c       either the control vector or the gradient vector can be read.
18    c
19    c       *dfile* should be the radix of the file: ecco_ctrl or ecco_cost
20    c
21    c     started: Christian Eckert eckert@mit.edu 12-Apr-2000
22    c
23    c     changed:  Patrick Heimbach heimbach@mit.edu 19-Jun-2000
24    c               - finished, revised and debugged
25    c
26    c     ==================================================================
27    c     SUBROUTINE optim_readdata
28    c     ==================================================================
29    
30          implicit none
31    
32    c     == global variables ==
33    
34    #include "EEPARAMS.h"
35    #include "SIZE.h"
36    cgg   Include ECCO_CPPOPTIONS because the ecco_ctrl,cost files
37    cgg   have headers with options for OBCS masks.
38    #include "ECCO_CPPOPTIONS.h"
39    
40    #include "ecco.h"
41    #include "ctrl.h"
42    #include "optim.h"
43    #include "minimization.h"
44    
45    c     == routine arguments ==
46    
47          integer nn
48          _RL     ff
49    
50    #if defined (DYNAMIC)
51          _RL     vv(nn)
52    #elif defined (USE_POINTER) || (MAX_INDEPEND == 0)
53          _RL            vv
54          pointer (pvv,vv(1))
55    #else
56          integer nmax
57          parameter( nmax = MAX_INDEPEND )
58          _RL   vv(nmax)
59    #endif
60    
61          character*(9) dfile
62          logical lheaderonly
63    
64    c     == local variables ==
65    
66          integer bi,bj
67          integer biG,bjG
68          integer i,j
69          integer ii,k
70          integer icvar
71          integer icvrec
72          integer icvcomp
73          integer icvoffset
74          integer nopt
75          integer funit
76    
77          integer cbuffindex
78          _RL     cbuff( sNx*nSx*nPx*sNy*nSy*nPy )
79    
80          character*(128) fname
81    
82          integer         filei
83          integer         filej
84          integer         filek
85          integer         filenopt
86          integer         fileig
87          integer         filejg
88          integer         filensx
89          integer         filensy
90          _RL             fileff
91    
92    cgg(
93          _RL     gg
94          integer igg
95          integer iobcs
96    cgg)
97    
98    c     == end of interface ==
99    
100    c--   The reference i/o unit.
101          funit = 20
102    
103    c--   Next optimization cycle.
104          nopt = optimcycle
105    
106          if      ( dfile .eq. ctrlname ) then
107            print*
108            print*,' OPTIM_READDATA: Reading control vector'
109            print*,'            for optimization cycle: ',nopt
110            print*
111          else if ( dfile .eq. costname ) then
112            print*
113            print*,' OPTIM_READDATA: Reading cost function and'
114            print*,'            gradient of cost function'
115            print*,'            for optimization cycle: ',nopt
116            print*
117          else
118            print*
119            print*,' OPTIM_READDATA: subroutine called by a false *dfile*'
120            print*,'            argument. *dfile* = ',dfile
121            print*
122            stop   '  ...  stopped in OPTIM_READDATA.'
123          endif
124    
125    c--   Read the data.
126    
127          bjG = 1 + (myygloballo - 1)/sny
128          biG = 1 + (myxgloballo - 1)/snx
129    
130    c--   Generate file name and open the file.
131          write(fname(1:128),'(4a,i4.4)')
132         &     dfile,'_',expId(1:10),'.opt', nopt
133          open( funit, file   = fname,
134         &     status = 'old',
135         &     form   = 'unformatted',
136         &     access = 'sequential'   )
137          print*, 'opened file ', fname
138    
139    c--   Read the header.
140          read( funit ) nvartype
141          read( funit ) nvarlength
142          read( funit ) expId
143          read( funit ) filenopt
144          read( funit ) fileff
145          read( funit ) fileiG
146          read( funit ) filejG
147          read( funit ) filensx
148          read( funit ) filensy
149    
150    cph(
151          print *,'ph-opt 1 ', nvartype, nvarlength, filensx, filensy
152    cph)
153    
154          read( funit ) (nWetcGlobal(k), k=1,nr)
155          read( funit ) (nWetsGlobal(k), k=1,nr)
156          read( funit ) (nWetwGlobal(k), k=1,nr)
157    #ifdef ALLOW_CTRL_WETV
158          read( funit ) (nWetvGlobal(k), k=1,nr)
159    #endif
160    
161    cgg(    Add OBCS Mask information into the header section for optimization.
162    #ifdef ALLOW_OBCSN_CONTROL
163          read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
164    #endif
165    #ifdef ALLOW_OBCSS_CONTROL
166          read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
167    #endif
168    #ifdef ALLOW_OBCSW_CONTROL
169          read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
170    #endif
171    #ifdef ALLOW_OBCSE_CONTROL
172          read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
173    #endif
174    cgg)
175          read( funit ) (ncvarindex(i), i=1,maxcvars)
176          read( funit ) (ncvarrecs(i),  i=1,maxcvars)
177          read( funit ) (ncvarxmax(i),  i=1,maxcvars)
178          read( funit ) (ncvarymax(i),  i=1,maxcvars)
179          read( funit ) (ncvarnrmax(i), i=1,maxcvars)
180          read( funit ) (ncvargrd(i),   i=1,maxcvars)
181          read( funit )
182    
183    cph(
184    cph      if (lheaderonly) then
185             print *, 'pathei: nvartype ', nvartype
186             print *, 'pathei: nvarlength ', nvarlength
187             print *, 'pathei: expId ', expId
188             print *, 'pathei: filenopt ', filenopt
189             print *, 'pathei: fileff ', fileff
190             print *, 'pathei: fileiG ', fileiG
191             print *, 'pathei: filejG ', filejG
192             print *, 'pathei: filensx ', filensx
193             print *, 'pathei: filensy ', filensy
194            
195             print *, 'pathei: nWetcGlobal ',
196         &        (nWetcGlobal(k),  k=1,nr)
197             print *, 'pathei: nWetsGlobal ',
198         &        (nWetsGlobal(k),  k=1,nr)
199             print *, 'pathei: nWetwGlobal ',
200         &        (nWetwGlobal(k),  k=1,nr)
201             print *, 'pathei: nWetvGlobal ',
202         &        (nWetvGlobal(k),  k=1,nr)
203             print *, 'pathei: ncvarindex ',
204         &        (ncvarindex(i), i=1,maxcvars)
205             print *, 'pathei: ncvarrecs ',
206         &        (ncvarrecs(i),  i=1,maxcvars)
207             print *, 'pathei: ncvarxmax ',
208         &        (ncvarxmax(i),  i=1,maxcvars)
209             print *, 'pathei: ncvarymax ',
210         &        (ncvarymax(i),  i=1,maxcvars)
211             print *, 'pathei: ncvarnrmax ',
212         &        (ncvarnrmax(i), i=1,maxcvars)
213             print *, 'pathei: ncvargrd ',
214         &        (ncvargrd(i),   i=1,maxcvars)
215    cph      end if
216    cph)
217    c--   Check the header information for consistency.
218    
219    cph      if ( filenopt .ne. nopt ) then
220    cph         print*
221    cph         print*,' READ_HEADER: Input data belong to the wrong'
222    cph         print*,'              optimization cycle.'
223    cph         print*,'              optimization cycle = ',nopt
224    cph         print*,'              input optim  cycle = ',filenopt
225    cph         print*
226    cph         stop   ' ... stopped in READ_HEADER.'
227    cph      endif
228          
229          if ( (fileiG .ne. biG) .or. (filejG .ne. bjG) ) then
230             print*
231             print*,' READ_HEADER: Tile indices of loop and data '
232             print*,'              do not match.'
233             print*,'              loop x/y component = ',
234         &        biG,bjG
235             print*,'              data x/y component = ',
236         &        fileiG,filejG
237             print*
238             stop   ' ... stopped in READ_HEADER.'
239          endif
240          
241          if ( (filensx .ne. nsx) .or. (filensy .ne. nsy) ) then
242             print*
243             print*,' READ_HEADER: Numbers of tiles do not match.'
244             print*,'              parameter x/y no. of tiles = ',
245         &        bi,bj
246             print*,'              data      x/y no. of tiles = ',
247         &        filensx,filensy
248             print*
249             stop   ' ... stopped in READ_HEADER.'
250          endif
251    
252    ce    Add some more checks. ...
253    
254          if (.NOT. lheaderonly) then
255    c--   Read the data.
256             icvoffset = 0
257             do icvar = 1,maxcvars
258                if ( ncvarindex(icvar) .ne. -1 ) then
259                   do icvrec = 1,ncvarrecs(icvar)
260                      do bj = 1,nsy
261                         do bi = 1,nsx
262                            read( funit ) ncvarindex(icvar)
263                            read( funit ) filej
264                            read( funit ) filei
265                            do k = 1,ncvarnrmax(icvar)
266                               cbuffindex = 0
267                               if (ncvargrd(icvar) .eq. 'c') then
268                                  cbuffindex = nWetcGlobal(k)
269                               else if (ncvargrd(icvar) .eq. 's') then
270                                  cbuffindex = nWetsGlobal(k)
271                               else if (ncvargrd(icvar) .eq. 'w') then
272                                  cbuffindex = nWetwGlobal(k)
273                               else if (ncvargrd(icvar) .eq. 'v') then
274                                  cbuffindex = nWetvGlobal(k)
275    cgg(   O.B. points have the grid mask "m".
276                               else if (ncvargrd(icvar) .eq. 'm') then
277    cgg    From "icvrec", calculate what iobcs must be.
278                                 gg   = (icvrec-1)/nobcs
279                                 igg  = int(gg)
280                                 iobcs= icvrec - igg*nobcs
281    #ifdef ALLOW_OBCSN_CONTROL
282                                 if (icvar .eq. 11) then                    
283                                   cbuffindex = nWetobcsnGlo(k,iobcs)
284                                 endif
285    #endif
286    #ifdef ALLOW_OBCSS_CONTROL
287                                 if (icvar .eq. 12) then
288                                   cbuffindex = nWetobcssGlo(k,iobcs)
289                                 endif
290    #endif
291    #ifdef ALLOW_OBCSW_CONTROL
292                                 if (icvar .eq. 13) then
293                                   cbuffindex = nWetobcswGlo(k,iobcs)
294                                 endif
295    #endif
296    #ifdef ALLOW_OBCSE_CONTROL
297                                 if (icvar .eq. 14) then
298                                   cbuffindex = nWetobcseGlo(k,iobcs)
299                                 endif
300    #endif
301    cgg)
302                               endif
303                               if (cbuffindex .gt. 0) then
304                                  read( funit ) cbuffindex
305                                  read( funit ) filek
306                                  read( funit ) (cbuff(ii), ii=1,cbuffindex)
307                                  do icvcomp = 1,cbuffindex
308                                     vv(icvoffset+icvcomp) = cbuff(icvcomp)
309    cgg( Right now, the changes to the open boundary velocities are not balanced.
310    cgg( The model will crash due to physical reasons.
311    cgg( However, we can optimize with respect to just O.B. T and S if the
312    cgg( next two lines are uncommented.
313    cgg                         if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0.
314    cgg                         if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0.
315                                  enddo
316                                  icvoffset = icvoffset + cbuffindex
317                               endif
318                            enddo
319                         enddo
320                      enddo
321                   enddo
322                endif
323             enddo
324    
325          else
326    
327    c--   Assign the number of control variables.
328             nn = nvarlength
329            
330          endif
331    
332          close( funit )
333    
334    c--   Assign the cost function value in case we read the cost file.
335    
336          if      ( dfile .eq. ctrlname ) then
337            ff = 0. d 0
338          else if ( dfile .eq. costname ) then
339            ff = fileff
340          endif
341    
342          return
343          end
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.2.4
  ViewVC Help
Powered by ViewVC 1.1.22