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

Diff of /MITgcm/optim/optim_readdata.F

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

revision 1.2 by heimbach, Fri Nov 15 04:03:25 2002 UTC revision 1.10 by mlosch, Mon Jul 25 08:58:47 2011 UTC
# Line 37  cgg   Include ECCO_CPPOPTIONS because th Line 37  cgg   Include ECCO_CPPOPTIONS because th
37  cgg   have headers with options for OBCS masks.  cgg   have headers with options for OBCS masks.
38  #include "ECCO_CPPOPTIONS.h"  #include "ECCO_CPPOPTIONS.h"
39    
 #include "ecco.h"  
40  #include "ctrl.h"  #include "ctrl.h"
41  #include "optim.h"  #include "optim.h"
42  #include "minimization.h"  #include "minimization.h"
# Line 75  c     == local variables == Line 74  c     == local variables ==
74        integer funit        integer funit
75    
76        integer cbuffindex        integer cbuffindex
77        _RL     cbuff( sNx*nSx*nPx*sNy*nSy*nPy )        real*4 cbuff( sNx*nSx*nPx*sNy*nSy*nPy )
78    
79        character*(128) fname        character*(128) fname
80    
81        integer         filei  c      integer         filei
82        integer         filej  c      integer         filej
83        integer         filek  c      integer         filek
84    c      integer         fileiG
85    c      integer         filejG
86    c      integer         filensx
87    c      integer         filensy
88        integer         filenopt        integer         filenopt
       integer         fileig  
       integer         filejg  
       integer         filensx  
       integer         filensy  
89        _RL             fileff        _RL             fileff
90    
91  cgg(  cgg(
# Line 97  cgg) Line 96  cgg)
96    
97  c     == end of interface ==  c     == end of interface ==
98    
99          print *, 'pathei-lsopt in optim_readdata'
100    
101  c--   The reference i/o unit.  c--   The reference i/o unit.
102        funit = 20        funit = 20
103    
# Line 129  c--   Read the data. Line 130  c--   Read the data.
130    
131  c--   Generate file name and open the file.  c--   Generate file name and open the file.
132        write(fname(1:128),'(4a,i4.4)')        write(fname(1:128),'(4a,i4.4)')
133       &     dfile,'_',expId(1:10),'.opt', nopt       &     dfile,'_',yctrlid(1:10),'.opt', nopt
134        open( funit, file   = fname,        open( funit, file   = fname,
135       &     status = 'old',       &     status = 'old',
136       &     form   = 'unformatted',       &     form   = 'unformatted',
# Line 139  c--   Generate file name and open the fi Line 140  c--   Generate file name and open the fi
140  c--   Read the header.  c--   Read the header.
141        read( funit ) nvartype        read( funit ) nvartype
142        read( funit ) nvarlength        read( funit ) nvarlength
143        read( funit ) expId        read( funit ) yctrlid
144        read( funit ) filenopt        read( funit ) filenopt
145        read( funit ) fileff        read( funit ) fileff
146        read( funit ) fileiG        read( funit ) fileiG
# Line 147  c--   Read the header. Line 148  c--   Read the header.
148        read( funit ) filensx        read( funit ) filensx
149        read( funit ) filensy        read( funit ) filensy
150    
151  cph(        read( funit ) (nWetcGlobal(k), k=1,nr)
152        print *,'ph-opt 1 ', nvartype, nvarlength, filensx, filensy        read( funit ) (nWetsGlobal(k), k=1,nr)
153  cph)        read( funit ) (nWetwGlobal(k), k=1,nr)
154    #ifdef ALLOW_CTRL_WETV
155        read( funit ) (((nWetcTile(i,j,k), i=1,nsx), j=1,nsy),        read( funit ) (nWetvGlobal(k), k=1,nr)
156       &     k=1,nr)  #endif
157        read( funit ) (((nWetsTile(i,j,k), i=1,nsx), j=1,nsy),  #ifdef ALLOW_SHIFWFLX_CONTROL
158       &     k=1,nr)        read(funit) (nWetiGlobal(k), k=1,nr)
159        read( funit ) (((nWetwTile(i,j,k), i=1,nsx), j=1,nsy),  c     read(funit) nWetiGlobal(1)
160       &     k=1,nr)  #endif
161    
162  cgg(    Add OBCS Mask information into the header section for optimization.  cgg(    Add OBCS Mask information into the header section for optimization.
163  #ifdef ALLOW_OBCSN_CONTROL  #ifdef ALLOW_OBCSN_CONTROL
164            read(funit) ((((nWetobcsn(i,j,k,iobcs), k=1,nr),        read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
      &          iobcs= 1,nobcs), i=1,nsx) , j=1,nsy)  
165  #endif  #endif
166  #ifdef ALLOW_OBCSS_CONTROL  #ifdef ALLOW_OBCSS_CONTROL
167            read(funit) ((((nWetobcss(i,j,k,iobcs), k=1,nr),        read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
      &          iobcs= 1,nobcs), i=1,nsx) , j=1,nsy)  
168  #endif  #endif
169  #ifdef ALLOW_OBCSW_CONTROL  #ifdef ALLOW_OBCSW_CONTROL
170            read(funit) ((((nWetobcsw(i,j,k,iobcs), k=1,nr),        read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
      &          iobcs= 1,nobcs), i=1,nsx) , j=1,nsy)  
171  #endif  #endif
172  #ifdef ALLOW_OBCSE_CONTROL  #ifdef ALLOW_OBCSE_CONTROL
173            read(funit) ((((nWetobcse(i,j,k,iobcs), k=1,nr),        read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs)
      &          iobcs= 1,nobcs), i=1,nsx) , j=1,nsy)  
174  #endif  #endif
175  cgg)  cgg)
176        read( funit ) (ncvarindex(i), i=1,maxcvars)        read( funit ) (ncvarindex(i), i=1,maxcvars)
# Line 188  cph( Line 185  cph(
185  cph      if (lheaderonly) then  cph      if (lheaderonly) then
186           print *, 'pathei: nvartype ', nvartype           print *, 'pathei: nvartype ', nvartype
187           print *, 'pathei: nvarlength ', nvarlength           print *, 'pathei: nvarlength ', nvarlength
188           print *, 'pathei: expId ', expId           print *, 'pathei: yctrlid ', yctrlid
189           print *, 'pathei: filenopt ', filenopt           print *, 'pathei: filenopt ', filenopt
190           print *, 'pathei: fileff ', fileff           print *, 'pathei: fileff ', fileff
191           print *, 'pathei: fileiG ', fileiG           print *, 'pathei: fileiG ', fileiG
# Line 196  cph      if (lheaderonly) then Line 193  cph      if (lheaderonly) then
193           print *, 'pathei: filensx ', filensx           print *, 'pathei: filensx ', filensx
194           print *, 'pathei: filensy ', filensy           print *, 'pathei: filensy ', filensy
195                    
196           print *, 'pathei: nWetcTile ',           print *, 'pathei: nWetcGlobal ',
197       &        (((nWetcTile(i,j,k), i=1,nsx), j=1,nsy), k=1,nr)       &        (nWetcGlobal(k),  k=1,nr)
198           print *, 'pathei: nWetsTile ',           print *, 'pathei: nWetsGlobal ',
199       &        (((nWetsTile(i,j,k), i=1,nsx), j=1,nsy), k=1,nr)       &        (nWetsGlobal(k),  k=1,nr)
200           print *, 'pathei: nWetwTile ',           print *, 'pathei: nWetwGlobal ',
201       &        (((nWetwTile(i,j,k), i=1,nsx), j=1,nsy), k=1,nr)       &        (nWetwGlobal(k),  k=1,nr)
202             print *, 'pathei: nWetvGlobal ',
203         &        (nWetvGlobal(k),  k=1,nr)
204    #ifdef ALLOW_SHIFWFLX_CONTROL
205             print *, 'pathei: nWetiGlobal ',
206         &        (nWetiGlobal(k), k=1,nr)
207    #endif
208    #ifdef ALLOW_OBCSN_CONTROL
209             do iobcs=1,nobcs
210              print *, 'pathei: nWetobcsnGlo (iobcs=', iobcs,')',
211         &         (nWetobcsnGlo(k,iobcs), k=1,nr)
212             enddo
213    #endif
214    #ifdef ALLOW_OBCSS_CONTROL
215             do iobcs=1,nobcs
216              print *, 'pathei: nWetobcssGlo (iobcs=', iobcs,')',
217         &         (nWetobcssGlo(k,iobcs), k=1,nr)
218             enddo
219    #endif
220    #ifdef ALLOW_OBCSW_CONTROL
221             do iobcs=1,nobcs
222              print *, 'pathei: nWetobcswGlo (iobcs=', iobcs,')',
223         &         (nWetobcswGlo(k,iobcs), k=1,nr)
224             enddo
225    #endif
226    #ifdef ALLOW_OBCSE_CONTROL
227             do iobcs=1,nobcs
228              print *, 'pathei: nWetobcseGlo (iobcs=', iobcs,')',
229         &         (nWetobcseGlo(k,iobcs), k=1,nr)
230             enddo
231    #endif
232           print *, 'pathei: ncvarindex ',           print *, 'pathei: ncvarindex ',
233       &        (ncvarindex(i), i=1,maxcvars)       &        (ncvarindex(i), i=1,maxcvars)
234           print *, 'pathei: ncvarrecs ',           print *, 'pathei: ncvarrecs ',
# Line 255  ce    Add some more checks. ... Line 282  ce    Add some more checks. ...
282    
283        if (.NOT. lheaderonly) then        if (.NOT. lheaderonly) then
284  c--   Read the data.  c--   Read the data.
285           icvoffset = 0         icvoffset = 0
286           do icvar = 1,maxcvars         do icvar = 1,maxcvars
287              if ( ncvarindex(icvar) .ne. -1 ) then          if ( ncvarindex(icvar) .ne. -1 ) then
288                 do icvrec = 1,ncvarrecs(icvar)           do icvrec = 1,ncvarrecs(icvar)
289                    do bj = 1,nsy            do bj = 1,nsy
290                       do bi = 1,nsx             do bi = 1,nsx
291                          read( funit ) ncvarindex(icvar)              read( funit ) ncvarindex(icvar)
292                          read( funit ) filej              read( funit ) filej
293                          read( funit ) filei              read( funit ) filei
294                          do k = 1,ncvarnrmax(icvar)              do k = 1,ncvarnrmax(icvar)
295                             cbuffindex = 0               cbuffindex = 0
296                             if (ncvargrd(icvar) .eq. 'c') then               if (ncvargrd(icvar) .eq. 'c') then
297                                cbuffindex = nwetctile(bi,bj,k)                cbuffindex = nWetcGlobal(k)
298                             else if (ncvargrd(icvar) .eq. 's') then               else if (ncvargrd(icvar) .eq. 's') then
299                                cbuffindex = nwetstile(bi,bj,k)                cbuffindex = nWetsGlobal(k)
300                             else if (ncvargrd(icvar) .eq. 'w') then               else if (ncvargrd(icvar) .eq. 'w') then
301                                cbuffindex = nwetwtile(bi,bj,k)                cbuffindex = nWetwGlobal(k)
302                 else if (ncvargrd(icvar) .eq. 'v') then
303                  cbuffindex = nWetvGlobal(k)
304    #ifdef ALLOW_SHIFWFLX_CONTROL
305                 else if (ncvargrd(icvar) .eq. 'i') then
306                  cbuffindex = nWetiGlobal(k)
307    #endif
308  cgg(   O.B. points have the grid mask "m".  cgg(   O.B. points have the grid mask "m".
309                             else if (ncvargrd(icvar) .eq. 'm') then               else if (ncvargrd(icvar) .eq. 'm') then
310  cgg    From "icvrec", calculate what iobcs must be.  cgg    From "icvrec", calculate what iobcs must be.
311                               gg   = (icvrec-1)/nobcs                gg   = (icvrec-1)/nobcs
312                               igg  = int(gg)                igg  = int(gg)
313                               iobcs= icvrec - igg*nobcs                iobcs= icvrec - igg*nobcs
314  #ifdef ALLOW_OBCSN_CONTROL  #ifdef ALLOW_OBCSN_CONTROL
315                               if (icvar .eq. 11) then                                    if (icvar .eq. 11) then                    
316                                 cbuffindex = nwetobcsn(bi,bj,k,iobcs)                 cbuffindex = nWetobcsnGlo(k,iobcs)
317                               endif                endif
318  #endif  #endif
319  #ifdef ALLOW_OBCSS_CONTROL  #ifdef ALLOW_OBCSS_CONTROL
320                               if (icvar .eq. 12) then                if (icvar .eq. 12) then
321                                 cbuffindex = nwetobcss(bi,bj,k,iobcs)                 cbuffindex = nWetobcssGlo(k,iobcs)
322                               endif                endif
323  #endif  #endif
324  #ifdef ALLOW_OBCSW_CONTROL  #ifdef ALLOW_OBCSW_CONTROL
325                               if (icvar .eq. 13) then                if (icvar .eq. 13) then
326                                 cbuffindex = nwetobcsw(bi,bj,k,iobcs)                 cbuffindex = nWetobcswGlo(k,iobcs)
327                               endif                endif
328  #endif  #endif
329  #ifdef ALLOW_OBCSE_CONTROL  #ifdef ALLOW_OBCSE_CONTROL
330                               if (icvar .eq. 14) then                if (icvar .eq. 14) then
331                                 cbuffindex = nwetobcse(bi,bj,k,iobcs)                 cbuffindex = nWetobcseGlo(k,iobcs)
332                               endif                endif
333  #endif  #endif
334  cgg)  cgg)
335                             endif               endif
336                             if (cbuffindex .gt. 0) then               if ( icvoffset + cbuffindex .gt. nvarlength ) then
337                                read( funit ) cbuffindex                print*
338                                read( funit ) filek                print *, ' ERROR:'
339                                read( funit ) (cbuff(ii), ii=1,cbuffindex)                print *, ' There are at least ', icvoffset+cbuffindex,
340                                do icvcomp = 1,cbuffindex       &             ' records in '//fname(1:28)//'.'
341                                   vv(icvoffset+icvcomp) = cbuff(icvcomp)                print *, ' This is more than expected from nvarlength =',
342  cgg( Right now, the changes to the open boundary velocities are not balanced.       &             nvarlength, '.'
343  cgg( The model will crash due to physical reasons.                print *, ' Something is wrong in the computation of '//
344  cgg( However, we can optimize with respect to just O.B. T and S if the       &             'the wet points or'
345  cgg( next two lines are uncommented.                print *, ' in computing the number of records in '//
346  cgg                         if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0.       &             'some variable(s).'
347  cgg                         if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0.                print *, '  ...  stopped in OPTIM_READDATA.'
348                                enddo                stop     '  ...  stopped in OPTIM_READDATA.'
349                                icvoffset = icvoffset + cbuffindex               endif
350                             endif               if (cbuffindex .gt. 0) then
351                          enddo                read( funit ) cbuffindex
352                       enddo                read( funit ) filek
353                    enddo                read( funit ) (cbuff(ii), ii=1,cbuffindex)
354                 enddo                do icvcomp = 1,cbuffindex
355              endif                 vv(icvoffset+icvcomp) = cbuff(icvcomp)
356    c     If you want to optimize with respect to just O.B. T and S
357    c     uncomment the next two lines.
358    c              if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0.
359    c              if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0.
360                  enddo
361                  icvoffset = icvoffset + cbuffindex
362                 endif
363                enddo
364               enddo
365              enddo
366           enddo           enddo
367            endif
368           enddo
369            
370        else        else
371    
372  c--   Assign the number of control variables.  c--   Assign the number of control variables.
373           nn = nvarlength         nn = nvarlength
374                    
375        endif        endif
376    
# Line 334  c--   Assign the number of control varia Line 379  c--   Assign the number of control varia
379  c--   Assign the cost function value in case we read the cost file.  c--   Assign the cost function value in case we read the cost file.
380    
381        if      ( dfile .eq. ctrlname ) then        if      ( dfile .eq. ctrlname ) then
382          ff = 0. d 0         ff = 0. d 0
383        else if ( dfile .eq. costname ) then        else if ( dfile .eq. costname ) then
384          ff = fileff         ff = fileff
385        endif        endif
386    
387        return        return
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22