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

Annotation of /MITgcm/optim/optim_readdata.F

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


Revision 1.1.2.4 - (hide annotations) (download)
Fri Mar 7 05:21:33 2003 UTC (21 years ago) by heimbach
Branch: ecco-branch
CVS Tags: ecco_c50_e32, ecco_c50_e33, ecco_c50_e30, ecco_c50_e31, ecco_c51_e34d, ecco_c51_e34e, ecco_c51_e34f, ecco_c51_e34g, ecco_c51_e34a, ecco_c51_e34b, ecco_c51_e34c, ecco_c50_e29, ecco_c50_e28, ecco_c50_e33a, ecco_c51_e34
Changes since 1.1.2.3: +2 -0 lines
maintaining backward compatibility of ctrl I/O without nwetv

1 heimbach 1.1.2.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 heimbach 1.1.2.2 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 heimbach 1.1.2.1
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 heimbach 1.1.2.2 cgg(
93     _RL gg
94     integer igg
95     integer iobcs
96     cgg)
97    
98 heimbach 1.1.2.1 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 heimbach 1.1.2.3 read( funit ) (nWetcGlobal(k), k=1,nr)
155     read( funit ) (nWetsGlobal(k), k=1,nr)
156     read( funit ) (nWetwGlobal(k), k=1,nr)
157 heimbach 1.1.2.4 #ifdef ALLOW_CTRL_WETV
158 heimbach 1.1.2.3 read( funit ) (nWetvGlobal(k), k=1,nr)
159 heimbach 1.1.2.4 #endif
160 heimbach 1.1.2.2
161     cgg( Add OBCS Mask information into the header section for optimization.
162     #ifdef ALLOW_OBCSN_CONTROL
163 heimbach 1.1.2.3 read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
164 heimbach 1.1.2.2 #endif
165     #ifdef ALLOW_OBCSS_CONTROL
166 heimbach 1.1.2.3 read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
167 heimbach 1.1.2.2 #endif
168     #ifdef ALLOW_OBCSW_CONTROL
169 heimbach 1.1.2.3 read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
170 heimbach 1.1.2.2 #endif
171     #ifdef ALLOW_OBCSE_CONTROL
172 heimbach 1.1.2.3 read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
173 heimbach 1.1.2.2 #endif
174     cgg)
175 heimbach 1.1.2.1 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 heimbach 1.1.2.3 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 heimbach 1.1.2.1 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 heimbach 1.1.2.3 cbuffindex = nWetcGlobal(k)
269 heimbach 1.1.2.1 else if (ncvargrd(icvar) .eq. 's') then
270 heimbach 1.1.2.3 cbuffindex = nWetsGlobal(k)
271 heimbach 1.1.2.1 else if (ncvargrd(icvar) .eq. 'w') then
272 heimbach 1.1.2.3 cbuffindex = nWetwGlobal(k)
273     else if (ncvargrd(icvar) .eq. 'v') then
274     cbuffindex = nWetvGlobal(k)
275 heimbach 1.1.2.2 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 heimbach 1.1.2.3 cbuffindex = nWetobcsnGlo(k,iobcs)
284 heimbach 1.1.2.2 endif
285     #endif
286     #ifdef ALLOW_OBCSS_CONTROL
287     if (icvar .eq. 12) then
288 heimbach 1.1.2.3 cbuffindex = nWetobcssGlo(k,iobcs)
289 heimbach 1.1.2.2 endif
290     #endif
291     #ifdef ALLOW_OBCSW_CONTROL
292     if (icvar .eq. 13) then
293 heimbach 1.1.2.3 cbuffindex = nWetobcswGlo(k,iobcs)
294 heimbach 1.1.2.2 endif
295     #endif
296     #ifdef ALLOW_OBCSE_CONTROL
297     if (icvar .eq. 14) then
298 heimbach 1.1.2.3 cbuffindex = nWetobcseGlo(k,iobcs)
299 heimbach 1.1.2.2 endif
300     #endif
301     cgg)
302 heimbach 1.1.2.1 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 heimbach 1.1.2.2 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 heimbach 1.1.2.1 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

  ViewVC Help
Powered by ViewVC 1.1.22