/[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.3 by heimbach, Fri Dec 6 01:47:35 2002 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          read( funit ) (nWetvGlobal(k), k=1,nr)
158    
159    cgg(    Add OBCS Mask information into the header section for optimization.
160    #ifdef ALLOW_OBCSN_CONTROL
161          read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
162    #endif
163    #ifdef ALLOW_OBCSS_CONTROL
164          read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
165    #endif
166    #ifdef ALLOW_OBCSW_CONTROL
167          read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
168    #endif
169    #ifdef ALLOW_OBCSE_CONTROL
170          read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
171    #endif
172    cgg)
173          read( funit ) (ncvarindex(i), i=1,maxcvars)
174          read( funit ) (ncvarrecs(i),  i=1,maxcvars)
175          read( funit ) (ncvarxmax(i),  i=1,maxcvars)
176          read( funit ) (ncvarymax(i),  i=1,maxcvars)
177          read( funit ) (ncvarnrmax(i), i=1,maxcvars)
178          read( funit ) (ncvargrd(i),   i=1,maxcvars)
179          read( funit )
180    
181    cph(
182    cph      if (lheaderonly) then
183             print *, 'pathei: nvartype ', nvartype
184             print *, 'pathei: nvarlength ', nvarlength
185             print *, 'pathei: expId ', expId
186             print *, 'pathei: filenopt ', filenopt
187             print *, 'pathei: fileff ', fileff
188             print *, 'pathei: fileiG ', fileiG
189             print *, 'pathei: filejG ', filejG
190             print *, 'pathei: filensx ', filensx
191             print *, 'pathei: filensy ', filensy
192            
193             print *, 'pathei: nWetcGlobal ',
194         &        (nWetcGlobal(k),  k=1,nr)
195             print *, 'pathei: nWetsGlobal ',
196         &        (nWetsGlobal(k),  k=1,nr)
197             print *, 'pathei: nWetwGlobal ',
198         &        (nWetwGlobal(k),  k=1,nr)
199             print *, 'pathei: nWetvGlobal ',
200         &        (nWetvGlobal(k),  k=1,nr)
201             print *, 'pathei: ncvarindex ',
202         &        (ncvarindex(i), i=1,maxcvars)
203             print *, 'pathei: ncvarrecs ',
204         &        (ncvarrecs(i),  i=1,maxcvars)
205             print *, 'pathei: ncvarxmax ',
206         &        (ncvarxmax(i),  i=1,maxcvars)
207             print *, 'pathei: ncvarymax ',
208         &        (ncvarymax(i),  i=1,maxcvars)
209             print *, 'pathei: ncvarnrmax ',
210         &        (ncvarnrmax(i), i=1,maxcvars)
211             print *, 'pathei: ncvargrd ',
212         &        (ncvargrd(i),   i=1,maxcvars)
213    cph      end if
214    cph)
215    c--   Check the header information for consistency.
216    
217    cph      if ( filenopt .ne. nopt ) then
218    cph         print*
219    cph         print*,' READ_HEADER: Input data belong to the wrong'
220    cph         print*,'              optimization cycle.'
221    cph         print*,'              optimization cycle = ',nopt
222    cph         print*,'              input optim  cycle = ',filenopt
223    cph         print*
224    cph         stop   ' ... stopped in READ_HEADER.'
225    cph      endif
226          
227          if ( (fileiG .ne. biG) .or. (filejG .ne. bjG) ) then
228             print*
229             print*,' READ_HEADER: Tile indices of loop and data '
230             print*,'              do not match.'
231             print*,'              loop x/y component = ',
232         &        biG,bjG
233             print*,'              data x/y component = ',
234         &        fileiG,filejG
235             print*
236             stop   ' ... stopped in READ_HEADER.'
237          endif
238          
239          if ( (filensx .ne. nsx) .or. (filensy .ne. nsy) ) then
240             print*
241             print*,' READ_HEADER: Numbers of tiles do not match.'
242             print*,'              parameter x/y no. of tiles = ',
243         &        bi,bj
244             print*,'              data      x/y no. of tiles = ',
245         &        filensx,filensy
246             print*
247             stop   ' ... stopped in READ_HEADER.'
248          endif
249    
250    ce    Add some more checks. ...
251    
252          if (.NOT. lheaderonly) then
253    c--   Read the data.
254             icvoffset = 0
255             do icvar = 1,maxcvars
256                if ( ncvarindex(icvar) .ne. -1 ) then
257                   do icvrec = 1,ncvarrecs(icvar)
258                      do bj = 1,nsy
259                         do bi = 1,nsx
260                            read( funit ) ncvarindex(icvar)
261                            read( funit ) filej
262                            read( funit ) filei
263                            do k = 1,ncvarnrmax(icvar)
264                               cbuffindex = 0
265                               if (ncvargrd(icvar) .eq. 'c') then
266                                  cbuffindex = nWetcGlobal(k)
267                               else if (ncvargrd(icvar) .eq. 's') then
268                                  cbuffindex = nWetsGlobal(k)
269                               else if (ncvargrd(icvar) .eq. 'w') then
270                                  cbuffindex = nWetwGlobal(k)
271                               else if (ncvargrd(icvar) .eq. 'v') then
272                                  cbuffindex = nWetvGlobal(k)
273    cgg(   O.B. points have the grid mask "m".
274                               else if (ncvargrd(icvar) .eq. 'm') then
275    cgg    From "icvrec", calculate what iobcs must be.
276                                 gg   = (icvrec-1)/nobcs
277                                 igg  = int(gg)
278                                 iobcs= icvrec - igg*nobcs
279    #ifdef ALLOW_OBCSN_CONTROL
280                                 if (icvar .eq. 11) then                    
281                                   cbuffindex = nWetobcsnGlo(k,iobcs)
282                                 endif
283    #endif
284    #ifdef ALLOW_OBCSS_CONTROL
285                                 if (icvar .eq. 12) then
286                                   cbuffindex = nWetobcssGlo(k,iobcs)
287                                 endif
288    #endif
289    #ifdef ALLOW_OBCSW_CONTROL
290                                 if (icvar .eq. 13) then
291                                   cbuffindex = nWetobcswGlo(k,iobcs)
292                                 endif
293    #endif
294    #ifdef ALLOW_OBCSE_CONTROL
295                                 if (icvar .eq. 14) then
296                                   cbuffindex = nWetobcseGlo(k,iobcs)
297                                 endif
298    #endif
299    cgg)
300                               endif
301                               if (cbuffindex .gt. 0) then
302                                  read( funit ) cbuffindex
303                                  read( funit ) filek
304                                  read( funit ) (cbuff(ii), ii=1,cbuffindex)
305                                  do icvcomp = 1,cbuffindex
306                                     vv(icvoffset+icvcomp) = cbuff(icvcomp)
307    cgg( Right now, the changes to the open boundary velocities are not balanced.
308    cgg( The model will crash due to physical reasons.
309    cgg( However, we can optimize with respect to just O.B. T and S if the
310    cgg( next two lines are uncommented.
311    cgg                         if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0.
312    cgg                         if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0.
313                                  enddo
314                                  icvoffset = icvoffset + cbuffindex
315                               endif
316                            enddo
317                         enddo
318                      enddo
319                   enddo
320                endif
321             enddo
322    
323          else
324    
325    c--   Assign the number of control variables.
326             nn = nvarlength
327            
328          endif
329    
330          close( funit )
331    
332    c--   Assign the cost function value in case we read the cost file.
333    
334          if      ( dfile .eq. ctrlname ) then
335            ff = 0. d 0
336          else if ( dfile .eq. costname ) then
337            ff = fileff
338          endif
339    
340          return
341          end
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.2.3
  ViewVC Help
Powered by ViewVC 1.1.22