/[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.12 - (hide annotations) (download)
Tue May 26 22:54:09 2015 UTC (8 years, 11 months ago) by gforget
Branch: MAIN
Changes since 1.11: +7 -4 lines
- add CTRL_OPTIONS.h that is needed to set maxcvars correctly
  when using generic controls (contributed by D. Amrhein).

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

  ViewVC Help
Powered by ViewVC 1.1.22