/[MITgcm]/MITgcm_contrib/mlosch/optim_m1qn3/optim_readdata.F
ViewVC logotype

Annotation of /MITgcm_contrib/mlosch/optim_m1qn3/optim_readdata.F

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


Revision 1.2 - (hide annotations) (download)
Fri Apr 27 09:50:46 2012 UTC (11 years, 11 months ago) by mlosch
Branch: MAIN
Changes since 1.1: +3 -6 lines
get rid of the USE_POINTER flag
introduce more allocate statements for DYNAMIC
clean up a little

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

  ViewVC Help
Powered by ViewVC 1.1.22