/[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.13 - (hide annotations) (download)
Tue Jun 2 14:49:13 2015 UTC (8 years, 9 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65n, checkpoint65m, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, HEAD
Changes since 1.12: +7 -6 lines
- remove un-necessary option files includes

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

  ViewVC Help
Powered by ViewVC 1.1.22