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

  ViewVC Help
Powered by ViewVC 1.1.22