/[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.10 - (hide annotations) (download)
Mon Jul 25 08:58:47 2011 UTC (12 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint63g, checkpoint63l, checkpoint63m, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c
Changes since 1.9: +74 -62 lines
change indentation so that I can fit more code in one line
add a check about nvarlength (that would have saved me a full afternoon)
adjust a comment

1 heimbach 1.2
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 "ctrl.h"
41     #include "optim.h"
42     #include "minimization.h"
43    
44     c == routine arguments ==
45    
46     integer nn
47     _RL ff
48    
49     #if defined (DYNAMIC)
50     _RL vv(nn)
51     #elif defined (USE_POINTER) || (MAX_INDEPEND == 0)
52     _RL vv
53     pointer (pvv,vv(1))
54     #else
55     integer nmax
56     parameter( nmax = MAX_INDEPEND )
57     _RL vv(nmax)
58     #endif
59    
60     character*(9) dfile
61     logical lheaderonly
62    
63     c == local variables ==
64    
65     integer bi,bj
66     integer biG,bjG
67     integer i,j
68     integer ii,k
69     integer icvar
70     integer icvrec
71     integer icvcomp
72     integer icvoffset
73     integer nopt
74     integer funit
75    
76     integer cbuffindex
77 heimbach 1.5 real*4 cbuff( sNx*nSx*nPx*sNy*nSy*nPy )
78 heimbach 1.2
79     character*(128) fname
80    
81 heimbach 1.6 c integer filei
82     c integer filej
83     c integer filek
84 dfer 1.7 c integer fileiG
85     c integer filejG
86 heimbach 1.6 c integer filensx
87     c integer filensy
88 heimbach 1.2 integer filenopt
89     _RL fileff
90    
91     cgg(
92     _RL gg
93     integer igg
94     integer iobcs
95     cgg)
96    
97     c == end of interface ==
98    
99 heimbach 1.6 print *, 'pathei-lsopt in optim_readdata'
100    
101 heimbach 1.2 c-- The reference i/o unit.
102     funit = 20
103    
104     c-- Next optimization cycle.
105     nopt = optimcycle
106    
107     if ( dfile .eq. ctrlname ) then
108     print*
109     print*,' OPTIM_READDATA: Reading control vector'
110     print*,' for optimization cycle: ',nopt
111     print*
112     else if ( dfile .eq. costname ) then
113     print*
114     print*,' OPTIM_READDATA: Reading cost function and'
115     print*,' gradient of cost function'
116     print*,' for optimization cycle: ',nopt
117     print*
118     else
119     print*
120     print*,' OPTIM_READDATA: subroutine called by a false *dfile*'
121     print*,' argument. *dfile* = ',dfile
122     print*
123     stop ' ... stopped in OPTIM_READDATA.'
124     endif
125    
126     c-- Read the data.
127    
128     bjG = 1 + (myygloballo - 1)/sny
129     biG = 1 + (myxgloballo - 1)/snx
130    
131     c-- Generate file name and open the file.
132     write(fname(1:128),'(4a,i4.4)')
133 dfer 1.7 & dfile,'_',yctrlid(1:10),'.opt', nopt
134 heimbach 1.2 open( funit, file = fname,
135     & status = 'old',
136     & form = 'unformatted',
137     & access = 'sequential' )
138     print*, 'opened file ', fname
139    
140     c-- Read the header.
141     read( funit ) nvartype
142     read( funit ) nvarlength
143 dfer 1.7 read( funit ) yctrlid
144 heimbach 1.2 read( funit ) filenopt
145     read( funit ) fileff
146     read( funit ) fileiG
147     read( funit ) filejG
148     read( funit ) filensx
149     read( funit ) filensy
150    
151 heimbach 1.3 read( funit ) (nWetcGlobal(k), k=1,nr)
152     read( funit ) (nWetsGlobal(k), k=1,nr)
153     read( funit ) (nWetwGlobal(k), k=1,nr)
154 heimbach 1.4 #ifdef ALLOW_CTRL_WETV
155 heimbach 1.3 read( funit ) (nWetvGlobal(k), k=1,nr)
156 heimbach 1.4 #endif
157 mlosch 1.8 #ifdef ALLOW_SHIFWFLX_CONTROL
158     read(funit) (nWetiGlobal(k), k=1,nr)
159     c read(funit) nWetiGlobal(1)
160     #endif
161 heimbach 1.2
162     cgg( Add OBCS Mask information into the header section for optimization.
163     #ifdef ALLOW_OBCSN_CONTROL
164 heimbach 1.3 read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
165 heimbach 1.2 #endif
166     #ifdef ALLOW_OBCSS_CONTROL
167 heimbach 1.3 read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
168 heimbach 1.2 #endif
169     #ifdef ALLOW_OBCSW_CONTROL
170 heimbach 1.3 read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
171 heimbach 1.2 #endif
172     #ifdef ALLOW_OBCSE_CONTROL
173 heimbach 1.3 read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
174 heimbach 1.2 #endif
175     cgg)
176     read( funit ) (ncvarindex(i), i=1,maxcvars)
177     read( funit ) (ncvarrecs(i), i=1,maxcvars)
178     read( funit ) (ncvarxmax(i), i=1,maxcvars)
179     read( funit ) (ncvarymax(i), i=1,maxcvars)
180     read( funit ) (ncvarnrmax(i), i=1,maxcvars)
181     read( funit ) (ncvargrd(i), i=1,maxcvars)
182     read( funit )
183    
184     cph(
185     cph if (lheaderonly) then
186     print *, 'pathei: nvartype ', nvartype
187     print *, 'pathei: nvarlength ', nvarlength
188 dfer 1.7 print *, 'pathei: yctrlid ', yctrlid
189 heimbach 1.2 print *, 'pathei: filenopt ', filenopt
190     print *, 'pathei: fileff ', fileff
191     print *, 'pathei: fileiG ', fileiG
192     print *, 'pathei: filejG ', filejG
193     print *, 'pathei: filensx ', filensx
194     print *, 'pathei: filensy ', filensy
195    
196 heimbach 1.3 print *, 'pathei: nWetcGlobal ',
197     & (nWetcGlobal(k), k=1,nr)
198     print *, 'pathei: nWetsGlobal ',
199     & (nWetsGlobal(k), k=1,nr)
200     print *, 'pathei: nWetwGlobal ',
201     & (nWetwGlobal(k), k=1,nr)
202     print *, 'pathei: nWetvGlobal ',
203     & (nWetvGlobal(k), k=1,nr)
204 mlosch 1.8 #ifdef ALLOW_SHIFWFLX_CONTROL
205     print *, 'pathei: nWetiGlobal ',
206     & (nWetiGlobal(k), k=1,nr)
207     #endif
208 mlosch 1.9 #ifdef ALLOW_OBCSN_CONTROL
209     do iobcs=1,nobcs
210     print *, 'pathei: nWetobcsnGlo (iobcs=', iobcs,')',
211     & (nWetobcsnGlo(k,iobcs), k=1,nr)
212     enddo
213     #endif
214     #ifdef ALLOW_OBCSS_CONTROL
215     do iobcs=1,nobcs
216     print *, 'pathei: nWetobcssGlo (iobcs=', iobcs,')',
217     & (nWetobcssGlo(k,iobcs), k=1,nr)
218     enddo
219     #endif
220     #ifdef ALLOW_OBCSW_CONTROL
221     do iobcs=1,nobcs
222     print *, 'pathei: nWetobcswGlo (iobcs=', iobcs,')',
223     & (nWetobcswGlo(k,iobcs), k=1,nr)
224     enddo
225     #endif
226     #ifdef ALLOW_OBCSE_CONTROL
227     do iobcs=1,nobcs
228     print *, 'pathei: nWetobcseGlo (iobcs=', iobcs,')',
229     & (nWetobcseGlo(k,iobcs), k=1,nr)
230     enddo
231     #endif
232 heimbach 1.2 print *, 'pathei: ncvarindex ',
233     & (ncvarindex(i), i=1,maxcvars)
234     print *, 'pathei: ncvarrecs ',
235     & (ncvarrecs(i), i=1,maxcvars)
236     print *, 'pathei: ncvarxmax ',
237     & (ncvarxmax(i), i=1,maxcvars)
238     print *, 'pathei: ncvarymax ',
239     & (ncvarymax(i), i=1,maxcvars)
240     print *, 'pathei: ncvarnrmax ',
241     & (ncvarnrmax(i), i=1,maxcvars)
242     print *, 'pathei: ncvargrd ',
243     & (ncvargrd(i), i=1,maxcvars)
244     cph end if
245     cph)
246     c-- Check the header information for consistency.
247    
248     cph if ( filenopt .ne. nopt ) then
249     cph print*
250     cph print*,' READ_HEADER: Input data belong to the wrong'
251     cph print*,' optimization cycle.'
252     cph print*,' optimization cycle = ',nopt
253     cph print*,' input optim cycle = ',filenopt
254     cph print*
255     cph stop ' ... stopped in READ_HEADER.'
256     cph endif
257    
258     if ( (fileiG .ne. biG) .or. (filejG .ne. bjG) ) then
259     print*
260     print*,' READ_HEADER: Tile indices of loop and data '
261     print*,' do not match.'
262     print*,' loop x/y component = ',
263     & biG,bjG
264     print*,' data x/y component = ',
265     & fileiG,filejG
266     print*
267     stop ' ... stopped in READ_HEADER.'
268     endif
269    
270     if ( (filensx .ne. nsx) .or. (filensy .ne. nsy) ) then
271     print*
272     print*,' READ_HEADER: Numbers of tiles do not match.'
273     print*,' parameter x/y no. of tiles = ',
274     & bi,bj
275     print*,' data x/y no. of tiles = ',
276     & filensx,filensy
277     print*
278     stop ' ... stopped in READ_HEADER.'
279     endif
280    
281     ce Add some more checks. ...
282    
283     if (.NOT. lheaderonly) then
284     c-- Read the data.
285 mlosch 1.10 icvoffset = 0
286     do icvar = 1,maxcvars
287     if ( ncvarindex(icvar) .ne. -1 ) then
288     do icvrec = 1,ncvarrecs(icvar)
289     do bj = 1,nsy
290     do bi = 1,nsx
291     read( funit ) ncvarindex(icvar)
292     read( funit ) filej
293     read( funit ) filei
294     do k = 1,ncvarnrmax(icvar)
295     cbuffindex = 0
296     if (ncvargrd(icvar) .eq. 'c') then
297     cbuffindex = nWetcGlobal(k)
298     else if (ncvargrd(icvar) .eq. 's') then
299     cbuffindex = nWetsGlobal(k)
300     else if (ncvargrd(icvar) .eq. 'w') then
301     cbuffindex = nWetwGlobal(k)
302     else if (ncvargrd(icvar) .eq. 'v') then
303     cbuffindex = nWetvGlobal(k)
304 mlosch 1.8 #ifdef ALLOW_SHIFWFLX_CONTROL
305 mlosch 1.10 else if (ncvargrd(icvar) .eq. 'i') then
306     cbuffindex = nWetiGlobal(k)
307 mlosch 1.8 #endif
308 heimbach 1.2 cgg( O.B. points have the grid mask "m".
309 mlosch 1.10 else if (ncvargrd(icvar) .eq. 'm') then
310 heimbach 1.2 cgg From "icvrec", calculate what iobcs must be.
311 mlosch 1.10 gg = (icvrec-1)/nobcs
312     igg = int(gg)
313     iobcs= icvrec - igg*nobcs
314 heimbach 1.2 #ifdef ALLOW_OBCSN_CONTROL
315 mlosch 1.10 if (icvar .eq. 11) then
316     cbuffindex = nWetobcsnGlo(k,iobcs)
317     endif
318 heimbach 1.2 #endif
319     #ifdef ALLOW_OBCSS_CONTROL
320 mlosch 1.10 if (icvar .eq. 12) then
321     cbuffindex = nWetobcssGlo(k,iobcs)
322     endif
323 heimbach 1.2 #endif
324     #ifdef ALLOW_OBCSW_CONTROL
325 mlosch 1.10 if (icvar .eq. 13) then
326     cbuffindex = nWetobcswGlo(k,iobcs)
327     endif
328 heimbach 1.2 #endif
329     #ifdef ALLOW_OBCSE_CONTROL
330 mlosch 1.10 if (icvar .eq. 14) then
331     cbuffindex = nWetobcseGlo(k,iobcs)
332     endif
333 heimbach 1.2 #endif
334     cgg)
335 mlosch 1.10 endif
336     if ( icvoffset + cbuffindex .gt. nvarlength ) then
337     print*
338     print *, ' ERROR:'
339     print *, ' There are at least ', icvoffset+cbuffindex,
340     & ' records in '//fname(1:28)//'.'
341     print *, ' This is more than expected from nvarlength =',
342     & nvarlength, '.'
343     print *, ' Something is wrong in the computation of '//
344     & 'the wet points or'
345     print *, ' in computing the number of records in '//
346     & 'some variable(s).'
347     print *, ' ... stopped in OPTIM_READDATA.'
348     stop ' ... stopped in OPTIM_READDATA.'
349     endif
350     if (cbuffindex .gt. 0) then
351     read( funit ) cbuffindex
352     read( funit ) filek
353     read( funit ) (cbuff(ii), ii=1,cbuffindex)
354     do icvcomp = 1,cbuffindex
355     vv(icvoffset+icvcomp) = cbuff(icvcomp)
356     c If you want to optimize with respect to just O.B. T and S
357     c uncomment the next two lines.
358     c if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0.
359     c if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0.
360     enddo
361     icvoffset = icvoffset + cbuffindex
362     endif
363     enddo
364     enddo
365     enddo
366 heimbach 1.2 enddo
367 mlosch 1.10 endif
368     enddo
369    
370 heimbach 1.2 else
371    
372     c-- Assign the number of control variables.
373 mlosch 1.10 nn = nvarlength
374 heimbach 1.2
375     endif
376    
377     close( funit )
378    
379     c-- Assign the cost function value in case we read the cost file.
380    
381     if ( dfile .eq. ctrlname ) then
382 mlosch 1.10 ff = 0. d 0
383 heimbach 1.2 else if ( dfile .eq. costname ) then
384 mlosch 1.10 ff = fileff
385 heimbach 1.2 endif
386    
387     return
388     end

  ViewVC Help
Powered by ViewVC 1.1.22