/[MITgcm]/MITgcm_contrib/ecco_utils/lbfgs_jpl_version/optim.2/optim_readdata.F
ViewVC logotype

Annotation of /MITgcm_contrib/ecco_utils/lbfgs_jpl_version/optim.2/optim_readdata.F

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


Revision 1.1 - (hide annotations) (download)
Wed Apr 3 23:36:08 2013 UTC (12 years, 4 months ago) by heimbach
Branch: MAIN
CVS Tags: HEAD
Add L-BFGS code adapted to ECCO Production by JPL

1 heimbach 1.1
2     subroutine optim_readdata(
3     I indic,
4     I nn,
5     I dfile,
6     I lheaderonly,
7     O ff,
8     O vv
9     & )
10    
11     c ==================================================================
12     c SUBROUTINE optim_readdata
13     c ==================================================================
14     c
15     c o Read the data written by the MITgcmUV state estimation setup and
16     c join them to one vector that is subsequently used by the minimi-
17     c zation algorithm "lsopt". Depending on the specified file name
18     c either the control vector or the gradient vector can be read.
19     c
20     c *dfile* should be the radix of the file: ecco_ctrl or ecco_cost
21     c
22     c started: Christian Eckert eckert@mit.edu 12-Apr-2000
23     c
24     c changed: Patrick Heimbach heimbach@mit.edu 19-Jun-2000
25     c - finished, revised and debugged
26     c
27     c ==================================================================
28     c SUBROUTINE optim_readdata
29     c ==================================================================
30    
31     implicit none
32    
33     c == global variables ==
34    
35     #include "EEPARAMS.h"
36     #include "SIZE.h"
37     #include "PARAMS.h"
38     #include "GRID.h"
39     cgg Include ECCO_CPPOPTIONS because the ecco_ctrl,cost files
40     cgg have headers with options for OBCS masks.
41     #include "ECCO_CPPOPTIONS.h"
42    
43     #include "ecco.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     _RS vv(nn)
55     #elif defined (USE_POINTER) || (MAX_INDEPEND == 0)
56     _RS vv
57     pointer (pvv,vv(1))
58     #else
59     integer nmax
60     parameter( nmax = MAX_INDEPEND )
61     _RS 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     integer indic
80    
81     integer cbuffindex
82     real*4 cbuff( sNx*nSx*nPx*sNy*nSy*nPy )
83    
84     character*(128) fname
85    
86     c integer filei
87     c integer filej
88     c integer filek
89     c integer fileig
90     c integer filejg
91     c integer filensx
92     c integer filensy
93     integer filenopt
94     _RL fileff
95    
96     cph(
97     _RL tmprescale
98     _RL delZnorm
99     _RL cvarlimit(maxcvars)
100     DATA delR /
101     & 10., 10., 15., 20., 20., 25., 35., 50., 75., 100.,
102     & 150., 200., 275., 350., 415., 450., 500., 500., 500., 500.,
103     & 500., 500., 500. /
104     cph)
105    
106     cgg(
107     _RL gg
108     integer igg
109     integer iobcs
110     cgg)
111    
112     c == end of interface ==
113    
114     print *, 'pathei-lsopt in optim_readdata'
115    
116     c-- The reference i/o unit.
117     funit = 20
118    
119     c-- Next optimization cycle.
120     nopt = optimcycle - indic
121    
122     if ( dfile .eq. ctrlname ) then
123     print*
124     print*,' OPTIM_READDATA: Reading control vector'
125     print*,' for optimization cycle: ',nopt
126     print*
127     else if ( dfile .eq. costname ) then
128     print*
129     print*,' OPTIM_READDATA: Reading cost function and'
130     print*,' gradient of cost function'
131     print*,' for optimization cycle: ',nopt
132     print*
133     else
134     print*
135     print*,' OPTIM_READDATA: subroutine called by a false *dfile*'
136     print*,' argument. *dfile* = ',dfile
137     print*
138     stop ' ... stopped in OPTIM_READDATA.'
139     endif
140    
141     c-- Read the data.
142    
143     bjG = 1 + (myygloballo - 1)/sny
144     biG = 1 + (myxgloballo - 1)/snx
145    
146     c-- Generate file name and open the file.
147     write(fname(1:128),'(4a,i4.4)')
148     & dfile,'_',yctrlid(1:10),'.opt', nopt
149     open( funit, file = fname,
150     & status = 'old',
151     & form = 'unformatted',
152     & access = 'sequential' )
153     print*, 'opened file ', fname
154    
155     c-- Read the header.
156     read( funit ) nvartype
157     read( funit ) nvarlength
158     read( funit ) yctrlid
159     read( funit ) filenopt
160     read( funit ) fileff
161     read( funit ) fileiG
162     read( funit ) filejG
163     read( funit ) filensx
164     read( funit ) filensy
165    
166     read( funit ) (nWetcGlobal(k), k=1,nr)
167     read( funit ) (nWetsGlobal(k), k=1,nr)
168     read( funit ) (nWetwGlobal(k), k=1,nr)
169     #ifdef ALLOW_CTRL_WETV
170     read( funit ) (nWetvGlobal(k), k=1,nr)
171     #endif
172    
173     cgg( Add OBCS Mask information into the header section for optimization.
174     #ifdef ALLOW_OBCSN_CONTROL
175     read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
176     #endif
177     #ifdef ALLOW_OBCSS_CONTROL
178     read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
179     #endif
180     #ifdef ALLOW_OBCSW_CONTROL
181     read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
182     #endif
183     #ifdef ALLOW_OBCSE_CONTROL
184     read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
185     #endif
186     cgg)
187     read( funit ) (ncvarindex(i), i=1,maxcvars)
188     read( funit ) (ncvarrecs(i), i=1,maxcvars)
189     read( funit ) (ncvarxmax(i), i=1,maxcvars)
190     read( funit ) (ncvarymax(i), i=1,maxcvars)
191     read( funit ) (ncvarnrmax(i), i=1,maxcvars)
192     read( funit ) (ncvargrd(i), i=1,maxcvars)
193     read( funit )
194    
195     cph(
196     cph if (lheaderonly) then
197     print *, 'pathei: nvartype ', nvartype
198     print *, 'pathei: nvarlength ', nvarlength
199     print *, 'pathei: yctrlid ', yctrlid
200     print *, 'pathei: filenopt ', filenopt
201     print *, 'pathei: fileff ', fileff
202     print *, 'pathei: fileiG ', fileiG
203     print *, 'pathei: filejG ', filejG
204     print *, 'pathei: filensx ', filensx
205     print *, 'pathei: filensy ', filensy
206    
207     print *, 'pathei: nWetcGlobal ',
208     & (nWetcGlobal(k), k=1,nr)
209     print *, 'pathei: nWetsGlobal ',
210     & (nWetsGlobal(k), k=1,nr)
211     print *, 'pathei: nWetwGlobal ',
212     & (nWetwGlobal(k), k=1,nr)
213     print *, 'pathei: nWetvGlobal ',
214     & (nWetvGlobal(k), k=1,nr)
215     print *, 'pathei: ncvarindex ',
216     & (ncvarindex(i), i=1,maxcvars)
217     print *, 'pathei: ncvarrecs ',
218     & (ncvarrecs(i), i=1,maxcvars)
219     print *, 'pathei: ncvarxmax ',
220     & (ncvarxmax(i), i=1,maxcvars)
221     print *, 'pathei: ncvarymax ',
222     & (ncvarymax(i), i=1,maxcvars)
223     print *, 'pathei: ncvarnrmax ',
224     & (ncvarnrmax(i), i=1,maxcvars)
225     print *, 'pathei: ncvargrd ',
226     & (ncvargrd(i), i=1,maxcvars)
227     cph end if
228     cph)
229     c-- Check the header information for consistency.
230    
231     cph if ( filenopt .ne. nopt ) then
232     cph print*
233     cph print*,' READ_HEADER: Input data belong to the wrong'
234     cph print*,' optimization cycle.'
235     cph print*,' optimization cycle = ',nopt
236     cph print*,' input optim cycle = ',filenopt
237     cph print*
238     cph stop ' ... stopped in READ_HEADER.'
239     cph endif
240    
241     if ( (fileiG .ne. biG) .or. (filejG .ne. bjG) ) then
242     print*
243     print*,' READ_HEADER: Tile indices of loop and data '
244     print*,' do not match.'
245     print*,' loop x/y component = ',
246     & biG,bjG
247     print*,' data x/y component = ',
248     & fileiG,filejG
249     print*
250     stop ' ... stopped in READ_HEADER.'
251     endif
252    
253     if ( (filensx .ne. nsx) .or. (filensy .ne. nsy) ) then
254     print*
255     print*,' READ_HEADER: Numbers of tiles do not match.'
256     print*,' parameter x/y no. of tiles = ',
257     & bi,bj
258     print*,' data x/y no. of tiles = ',
259     & filensx,filensy
260     print*
261     stop ' ... stopped in READ_HEADER.'
262     endif
263    
264     ce Add some more checks. ...
265    
266     if (.NOT. lheaderonly) then
267     c-- Read the data.
268     icvoffset = 0
269     do icvar = 1,maxcvars
270     cph(
271     cph temporary rescale to compensate for new wtheta, wsalt
272     cph and maintain scaling consistency between
273     cph ecco_ctrl, ecco_cost
274     cph if ( dfile .eq. ctrlname ) then
275     if ( ( icvar.EQ.25 .OR. icvar.EQ.26 )
276     & .AND. dfile .eq. costname ) then
277     cph tmprescale = 1./sqrt(2.)
278     tmprescale = 1.
279     else
280     tmprescale = 1.
281     endif
282     cph)
283     if ( ncvarindex(icvar) .ne. -1 ) then
284     do icvrec = 1,ncvarrecs(icvar)
285     do bj = 1,nsy
286     do bi = 1,nsx
287     read( funit ) ncvarindex(icvar)
288     read( funit ) filej
289     read( funit ) filei
290     do k = 1,ncvarnrmax(icvar)
291     cbuffindex = 0
292     if (ncvargrd(icvar) .eq. 'c') then
293     cbuffindex = nWetcGlobal(k)
294     else if (ncvargrd(icvar) .eq. 's') then
295     cbuffindex = nWetsGlobal(k)
296     else if (ncvargrd(icvar) .eq. 'w') then
297     cbuffindex = nWetwGlobal(k)
298     else if (ncvargrd(icvar) .eq. 'v') then
299     cbuffindex = nWetvGlobal(k)
300     cgg( O.B. points have the grid mask "m".
301     else if (ncvargrd(icvar) .eq. 'm') then
302     cgg From "icvrec", calculate what iobcs must be.
303     gg = (icvrec-1)/nobcs
304     igg = int(gg)
305     iobcs= icvrec - igg*nobcs
306     #ifdef ALLOW_OBCSN_CONTROL
307     if (icvar .eq. 11) then
308     cbuffindex = nWetobcsnGlo(k,iobcs)
309     endif
310     #endif
311     #ifdef ALLOW_OBCSS_CONTROL
312     if (icvar .eq. 12) then
313     cbuffindex = nWetobcssGlo(k,iobcs)
314     endif
315     #endif
316     #ifdef ALLOW_OBCSW_CONTROL
317     if (icvar .eq. 13) then
318     cbuffindex = nWetobcswGlo(k,iobcs)
319     endif
320     #endif
321     #ifdef ALLOW_OBCSE_CONTROL
322     if (icvar .eq. 14) then
323     cbuffindex = nWetobcseGlo(k,iobcs)
324     endif
325     #endif
326     cgg)
327     endif
328     if (cbuffindex .gt. 0) then
329     read( funit ) cbuffindex
330     read( funit ) filek
331     read( funit ) (cbuff(ii), ii=1,cbuffindex)
332    
333     !uniform limit for all variables
334     cvarlimit(icvar) = 1.e4
335     do icvcomp = 1,cbuffindex
336     !cph(
337     ! cvarlimit(1) = 40.
338     ! cvarlimit(2) = 80.
339     ! cvarlimit(7) = 20.
340     ! cvarlimit(8) = 20.
341     ! cvarlimit(9) = 40.
342     ! cvarlimit(10) = 40.
343     !!bnc
344     ! cvarlimit(27) = 40.
345     ! cvarlimit(28) = 40.
346     ! cvarlimit(29) = 40.
347     !!bnc
348     ! cvarlimit(32) = 20.
349     ! cvarlimit(34) = 20.
350    
351     cph temporary rescale to compensate for new scaling wtheta, wsalt
352    
353     vv(icvoffset+icvcomp) = cbuff(icvcomp)
354     & *tmprescale
355     c
356     cph goto 1234
357     c
358     if ( ( icvar.EQ.3 .OR. icvar.EQ.4 .OR.
359     & icvar.EQ.5 .OR. icvar.EQ.6 .OR.
360     & icvar.EQ.7 .OR. icvar.EQ.8 .OR.
361     & icvar.EQ.9 .OR. icvar.EQ.10 .OR.
362     !bnc
363     & icvar.EQ.27 .OR. icvar.EQ.28 .OR.icvar.EQ.29 .OR.
364     !bnc
365     & icvar.EQ.32 .OR. icvar.EQ.34 ) .AND.
366     & ( dfile .eq. costname ) ) then
367     if ( cbuff(icvcomp) .gt. cvarlimit(icvar) ) then
368     vv(icvoffset+icvcomp) = cvarlimit(icvar)
369     else if ( cbuff(icvcomp) .lt. -1.*cvarlimit(icvar) ) then
370     vv(icvoffset+icvcomp) = -1.*cvarlimit(icvar)
371     endif
372     c
373     cph(
374     cph this part to be removed
375     cph else if ( ( icvar.EQ.1 .OR. icvar.EQ.2 ) .AND.
376     cph & ( dfile .eq. ctrlname ) ) then
377     cph vv(icvoffset+icvcomp) = cbuff(icvcomp)
378     cph & *tmprescale/SQRT(delR(1)/delR(k))
379     cph)
380     c
381     else if ( ( icvar.EQ.1 .OR. icvar.EQ.2 ) .AND.
382     & ( dfile .eq. costname ) ) then
383     cno vv(icvoffset+icvcomp) = cbuff(icvcomp)
384     cph & *tmprescale*SQRT(delR(1)/delR(k))
385     c
386     if ( vv(icvoffset+icvcomp) .gt. cvarlimit(icvar) ) then
387     vv(icvoffset+icvcomp) = cvarlimit(icvar)
388     else if ( vv(icvoffset+icvcomp) .lt. -1.*cvarlimit(icvar) ) then
389     vv(icvoffset+icvcomp) = -1.*cvarlimit(icvar)
390     endif
391     endif
392     c
393     1234 continue
394     c
395     cph)
396    
397     cgg( Right now, the changes to the open boundary velocities are not balanced.
398     cgg( The model will crash due to physical reasons.
399     cgg( However, we can optimize with respect to just O.B. T and S if the
400     cgg( next two lines are uncommented.
401     cgg if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0.
402     cgg if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0.
403     enddo
404     icvoffset = icvoffset + cbuffindex
405     endif
406     enddo
407     enddo
408     enddo
409     enddo
410     endif
411     enddo
412    
413     else
414    
415     c-- Assign the number of control variables.
416     nn = nvarlength
417    
418     endif
419    
420     close( funit )
421    
422     c-- Assign the cost function value in case we read the cost file.
423    
424     if ( dfile .eq. ctrlname ) then
425     ! ff = 0. d 0
426     ff = fc
427     else if ( dfile .eq. costname ) then
428     ff = fileff
429     endif
430    
431     return
432     end

  ViewVC Help
Powered by ViewVC 1.1.22