/[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.1.2.2 - (hide annotations) (download)
Thu Apr 4 10:20:16 2002 UTC (22 years ago) by heimbach
Branch: ecco-branch
CVS Tags: ecco_ice2, ecco_ice1, ecco_c44_e25, ecco_c44_e22, ecco_c44_e23, ecco_c44_e21, ecco_c44_e24
Branch point for: c24_e25_ice
Changes since 1.1.2.1: +61 -0 lines
Added obcs control part for lsopt I/O
(by G. Gebbie).

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

  ViewVC Help
Powered by ViewVC 1.1.22