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

Contents of /MITgcm/optim/optim_readdata.F

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


Revision 1.13 - (show 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
2 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
10 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 real*4 cbuff( sNx*nSx*nPx*sNy*nSy*nPy )
82
83 character*(128) fname
84
85 c integer filei
86 c integer filej
87 c integer filek
88 c integer fileiG
89 c integer filejG
90 c integer filensx
91 c integer filensy
92 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 print *, 'pathei-lsopt in optim_readdata'
104
105 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 & dfile,'_',yctrlid(1:10),'.opt', nopt
138 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 read( funit ) yctrlid
148 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 read( funit ) (nWetcGlobal(k), k=1,nr)
156 read( funit ) (nWetsGlobal(k), k=1,nr)
157 read( funit ) (nWetwGlobal(k), k=1,nr)
158 #ifdef ALLOW_CTRL_WETV
159 read( funit ) (nWetvGlobal(k), k=1,nr)
160 #endif
161 #ifdef ALLOW_SHIFWFLX_CONTROL
162 read(funit) (nWetiGlobal(k), k=1,nr)
163 c read(funit) nWetiGlobal(1)
164 #endif
165
166 cgg( Add OBCS Mask information into the header section for optimization.
167 #ifdef ALLOW_OBCSN_CONTROL
168 read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
169 #endif
170 #ifdef ALLOW_OBCSS_CONTROL
171 read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
172 #endif
173 #ifdef ALLOW_OBCSW_CONTROL
174 read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
175 #endif
176 #ifdef ALLOW_OBCSE_CONTROL
177 read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
178 #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 print *, 'pathei: yctrlid ', yctrlid
193 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 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 #ifdef ALLOW_SHIFWFLX_CONTROL
209 print *, 'pathei: nWetiGlobal ',
210 & (nWetiGlobal(k), k=1,nr)
211 #endif
212 #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 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 icvoffset = 0
290 do icvar = 1,maxcvars
291 if ( ncvarindex(icvar) .ne. -1 ) then
292 do icvrec = 1,ncvarrecs(icvar)
293 cph do bj = 1,nsy
294 cph do bi = 1,nsx
295 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 #ifdef ALLOW_SHIFWFLX_CONTROL
309 else if (ncvargrd(icvar) .eq. 'i') then
310 cbuffindex = nWetiGlobal(k)
311 #endif
312 cgg( O.B. points have the grid mask "m".
313 else if (ncvargrd(icvar) .eq. 'm') then
314 cgg From "icvrec", calculate what iobcs must be.
315 gg = (icvrec-1)/nobcs
316 igg = int(gg)
317 iobcs= icvrec - igg*nobcs
318 #ifdef ALLOW_OBCSN_CONTROL
319 if (icvar .eq. 11) then
320 cbuffindex = nWetobcsnGlo(k,iobcs)
321 endif
322 #endif
323 #ifdef ALLOW_OBCSS_CONTROL
324 if (icvar .eq. 12) then
325 cbuffindex = nWetobcssGlo(k,iobcs)
326 endif
327 #endif
328 #ifdef ALLOW_OBCSW_CONTROL
329 if (icvar .eq. 13) then
330 cbuffindex = nWetobcswGlo(k,iobcs)
331 endif
332 #endif
333 #ifdef ALLOW_OBCSE_CONTROL
334 if (icvar .eq. 14) then
335 cbuffindex = nWetobcseGlo(k,iobcs)
336 endif
337 #endif
338 cgg)
339 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 cph enddo
369 cph enddo
370 enddo
371 endif
372 enddo
373
374 else
375
376 c-- Assign the number of control variables.
377 nn = nvarlength
378
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 ff = 0. d 0
387 else if ( dfile .eq. costname ) then
388 ff = fileff
389 endif
390
391 return
392 end

  ViewVC Help
Powered by ViewVC 1.1.22