Back to home page

MITgcm
 
 

    

File indexing completed on 2024-03-02 06:10:02 UTC

view on githubraw file Latest commit 5cf43646 on 2024-03-01 18:50:49 UTC
65754df434 Mart*0001 C     ECCO_CPPOPTIONS used to affect maxcvars and defined ALLOW_OBCS?_CONTROL
                0002 C#include "ECCO_CPPOPTIONS.h"
                0003 C     now:
                0004 C     CTRL_OPTIONS affects maxcvars and may define ALLOW_OBCS?_CONTROL
eab6982bba Gael*0005 #include "CTRL_OPTIONS.h"
1b116cf4a9 Gael*0006 
4cee17c1be Patr*0007       subroutine optim_readdata(
                0008      I                      nn,
                0009      I                      dfile,
                0010      I                      lheaderonly,
                0011      O                      ff,
                0012      O                      vv
                0013      &                    )
                0014 
                0015 c     ==================================================================
                0016 c     SUBROUTINE optim_readdata
                0017 c     ==================================================================
                0018 c
                0019 c     o Read the data written by the MITgcmUV state estimation setup and
                0020 c       join them to one vector that is subsequently used by the minimi-
                0021 c       zation algorithm "lsopt". Depending on the specified file name
                0022 c       either the control vector or the gradient vector can be read.
                0023 c
                0024 c       *dfile* should be the radix of the file: ecco_ctrl or ecco_cost
                0025 c
                0026 c     started: Christian Eckert eckert@mit.edu 12-Apr-2000
                0027 c
                0028 c     changed:  Patrick Heimbach heimbach@mit.edu 19-Jun-2000
                0029 c               - finished, revised and debugged
                0030 c
                0031 c     ==================================================================
                0032 c     SUBROUTINE optim_readdata
                0033 c     ==================================================================
                0034 
5cf4364659 Mart*0035       IMPLICIT NONE
4cee17c1be Patr*0036 
                0037 c     == global variables ==
                0038 
                0039 #include "EEPARAMS.h"
                0040 #include "SIZE.h"
5cf4364659 Mart*0041 #include "CTRL_SIZE.h"
65754df434 Mart*0042 #ifdef ALLOW_OBCS_CONTROL
                0043 # include "CTRL_OBCS.h"
                0044 #endif
                0045 #include "CTRL.h"
4cee17c1be Patr*0046 #include "optim.h"
                0047 
                0048 c     == routine arguments ==
                0049 
                0050       integer nn
                0051       _RL     ff
                0052 
65754df434 Mart*0053 #ifdef DYNAMIC
4cee17c1be Patr*0054       _RL     vv(nn)
                0055 #else
                0056       integer nmax
                0057       parameter( nmax = MAX_INDEPEND )
                0058       _RL   vv(nmax)
                0059 #endif
                0060 
                0061       character*(9) dfile
                0062       logical lheaderonly
                0063 
                0064 c     == local variables ==
5cf4364659 Mart*0065 #ifdef READ_OLD_CTRL_PACK_FILE
                0066       INTEGER maxLocal
                0067       PARAMETER( maxLocal = old_maxcvars )
                0068       INTEGER file_varIndex(maxLocal)
                0069       INTEGER file_varRecs(maxLocal)
                0070       INTEGER file_varNxMax(maxLocal)
                0071       INTEGER file_varNyMax(maxLocal)
                0072       INTEGER file_varNrMax(maxLocal)
                0073       CHARACTER*( 1) file_varGrid(maxLocal)
                0074       CHARACTER*( 5) file_varType(maxLocal)
                0075 #endif
4cee17c1be Patr*0076 
65754df434 Mart*0077 CM      integer bi,bj
4cee17c1be Patr*0078       integer biG,bjG
                0079       integer i,j
                0080       integer ii,k
5cf4364659 Mart*0081       integer ivar
4cee17c1be Patr*0082       integer icvrec
                0083       integer icvcomp
                0084       integer icvoffset
                0085       integer nopt
                0086       integer funit
                0087 
                0088       integer cbuffindex
f50e6c1777 Patr*0089       real*4 cbuff( sNx*nSx*nPx*sNy*nSy*nPy )
4cee17c1be Patr*0090 
                0091       character*(128) fname
65754df434 Mart*0092       character*(17)  prefix
                0093       parameter ( prefix =  " OPTIM_READDATA: " )
4cee17c1be Patr*0094 
5cf4364659 Mart*0095       integer         ncvarindex_loc
                0096       integer         filei
                0097       integer         filej
                0098       integer         filek
                0099       integer         fileiG
                0100       integer         filejG
                0101       integer         filensx
                0102       integer         filensy
4cee17c1be Patr*0103       integer         filenopt
                0104       _RL             fileff
                0105 
                0106 cgg(
                0107       _RL     gg
                0108       integer igg
                0109       integer iobcs
                0110 cgg)
                0111 
                0112 c     == end of interface ==
                0113 
                0114 c--   The reference i/o unit.
                0115       funit = 20
                0116 
                0117 c--   Next optimization cycle.
                0118       nopt = optimcycle
                0119 
                0120       if      ( dfile .eq. ctrlname ) then
65754df434 Mart*0121        print*
                0122        print*,' OPTIM_READDATA: Reading control vector'
                0123        print*,'                 for optimization cycle: ',nopt
4cee17c1be Patr*0124       else if ( dfile .eq. costname ) then
65754df434 Mart*0125        print*
                0126        print*,' OPTIM_READDATA: Reading cost function and'//
                0127      &                        ' gradient of cost function'
                0128        print*,'                 for optimization cycle: ',nopt
4cee17c1be Patr*0129       else
65754df434 Mart*0130        print*
                0131        print*,' OPTIM_READDATA: subroutine called by a false *dfile*'
                0132        print*,'                 argument. *dfile* = ',dfile
                0133        print*
                0134        stop   '  ...  stopped in OPTIM_READDATA.'
4cee17c1be Patr*0135       endif
65754df434 Mart*0136       if ( lheaderonly ) then
                0137        print*,'                 ... header only'
                0138       endif
                0139       print*, ' '
4cee17c1be Patr*0140 
                0141 c--   Read the data.
                0142 
5cf4364659 Mart*0143       bjG = 1 + (myYGlobalLo - 1)/sNy
                0144       biG = 1 + (myXGlobalLo - 1)/sNx
4cee17c1be Patr*0145 
                0146 c--   Generate file name and open the file.
                0147       write(fname(1:128),'(4a,i4.4)')
91d99130a0 Davi*0148      &     dfile,'_',yctrlid(1:10),'.opt', nopt
4cee17c1be Patr*0149       open( funit, file   = fname,
                0150      &     status = 'old',
                0151      &     form   = 'unformatted',
                0152      &     access = 'sequential'   )
65754df434 Mart*0153       print*, prefix, 'opened file ', fname
4cee17c1be Patr*0154 
                0155 c--   Read the header.
                0156       read( funit ) nvartype
                0157       read( funit ) nvarlength
91d99130a0 Davi*0158       read( funit ) yctrlid
4cee17c1be Patr*0159       read( funit ) filenopt
                0160       read( funit ) fileff
65754df434 Mart*0161 C     filei/jG are dummy values and not used
4cee17c1be Patr*0162       read( funit ) fileiG
                0163       read( funit ) filejG
                0164       read( funit ) filensx
                0165       read( funit ) filensy
                0166 
5cf4364659 Mart*0167       read( funit ) (nWetcGlobal(k), k=1,Nr)
                0168       read( funit ) (nWetsGlobal(k), k=1,Nr)
                0169       read( funit ) (nWetwGlobal(k), k=1,Nr)
                0170 c#ifdef ALLOW_CTRL_WETV
                0171 c     read( funit ) (nWetvGlobal(k), k=1,Nr)
                0172 c#endif
e189f4121c Mart*0173 #ifdef ALLOW_SHIFWFLX_CONTROL
5cf4364659 Mart*0174       read(funit) (nWetiGlobal(k), k=1,Nr)
e189f4121c Mart*0175 c     read(funit) nWetiGlobal(1)
                0176 #endif
4cee17c1be Patr*0177 
                0178 cgg(    Add OBCS Mask information into the header section for optimization.
65754df434 Mart*0179 CML#if (defined (ALLOW_CTRL) && defined (ALLOW_OBCS))
4cee17c1be Patr*0180 #ifdef ALLOW_OBCSN_CONTROL
5cf4364659 Mart*0181       read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,Nr),iobcs= 1,nobcs)
4cee17c1be Patr*0182 #endif
                0183 #ifdef ALLOW_OBCSS_CONTROL
5cf4364659 Mart*0184       read( funit ) ((nWetobcssGlo(k,iobcs), k=1,Nr),iobcs= 1,nobcs)
4cee17c1be Patr*0185 #endif
                0186 #ifdef ALLOW_OBCSW_CONTROL
5cf4364659 Mart*0187       read( funit ) ((nWetobcswGlo(k,iobcs), k=1,Nr),iobcs= 1,nobcs)
4cee17c1be Patr*0188 #endif
                0189 #ifdef ALLOW_OBCSE_CONTROL
5cf4364659 Mart*0190       read( funit ) ((nWetobcseGlo(k,iobcs), k=1,Nr),iobcs= 1,nobcs)
4cee17c1be Patr*0191 #endif
65754df434 Mart*0192 CML#endif /* ALLOW_CTRL and ALLOW_OBCS */
4cee17c1be Patr*0193 cgg)
5cf4364659 Mart*0194 #ifdef READ_OLD_CTRL_PACK_FILE
                0195       DO ivar = 1,maxLocal
                0196         file_varType(ivar)  = '     '
                0197       ENDDO
                0198       read( funit ) (file_varIndex(ivar), ivar=1,maxLocal)
                0199       read( funit ) (file_varRecs(ivar),  ivar=1,maxLocal)
                0200       read( funit ) (file_varNxMax(ivar), ivar=1,maxLocal)
                0201       read( funit ) (file_varNyMax(ivar), ivar=1,maxLocal)
                0202       read( funit ) (file_varNrMax(ivar), ivar=1,maxLocal)
                0203       read( funit ) (file_varGrid(ivar),  ivar=1,maxLocal)
4cee17c1be Patr*0204       read( funit )
5cf4364659 Mart*0205       CALL CTRL_CONVERT_HEADER(
                0206      I              maxLocal, nvartype, 6,
                0207      U              file_varIndex, file_varRecs,
                0208      U              file_varNxMax, file_varNyMax, file_varNrMax,
                0209      U              file_varGrid, file_varType,
                0210      I              1 )
                0211       DO ivar = 1,nvartype
                0212         ncvarindex(ivar) = file_varIndex(ivar)
                0213         ncvarrecs(ivar)  = file_varRecs(ivar)
                0214         ncvarxmax(ivar)  = file_varNxMax(ivar)
                0215         ncvarymax(ivar)  = file_varNyMax(ivar)
                0216         ncvarnrmax(ivar) = file_varNrMax(ivar)
                0217         ncvargrd(ivar)   = file_varGrid(ivar)
                0218         ncvartype(ivar)  = file_varType(ivar)
                0219       ENDDO
                0220 #else /* READ_OLD_CTRL_PACK_FILE */
                0221       read( funit ) (ncvarindex(ivar), ivar=1,nvartype)
                0222       read( funit ) (ncvarrecs(ivar),  ivar=1,nvartype)
                0223       read( funit ) (ncvarxmax(ivar),  ivar=1,nvartype)
                0224       read( funit ) (ncvarymax(ivar),  ivar=1,nvartype)
                0225       read( funit ) (ncvarnrmax(ivar), ivar=1,nvartype)
                0226       read( funit ) (ncvargrd(ivar),   ivar=1,nvartype)
                0227       read( funit ) (ncvartype(ivar),  ivar=1,nvartype)
                0228 #endif /* READ_OLD_CTRL_PACK_FILE */
4cee17c1be Patr*0229 
65754df434 Mart*0230       print *, prefix, 'nvartype ', nvartype
                0231       print *, prefix, 'nvarlength ', nvarlength
                0232       print *, prefix, 'yctrlid ', yctrlid
                0233       print *, prefix, 'filenopt ', filenopt
                0234       print *, prefix, 'fileff ', fileff
                0235       print *, prefix, 'fileiG ', fileiG
                0236       print *, prefix, 'filejG ', filejG
                0237       print *, prefix, 'filensx ', filensx
                0238       print *, prefix, 'filensy ', filensy
                0239 
                0240       if (lheaderonly) then
5cf4364659 Mart*0241        print *, prefix, 'nWetcGlobal ', (nWetcGlobal(k), k=1,Nr)
                0242        print *, prefix, 'nWetsGlobal ', (nWetsGlobal(k), k=1,Nr)
                0243        print *, prefix, 'nWetwGlobal ', (nWetwGlobal(k), k=1,Nr)
                0244 c      print *, prefix, 'nWetvGlobal ', (nWetvGlobal(k), k=1,Nr)
e189f4121c Mart*0245 #ifdef ALLOW_SHIFWFLX_CONTROL
5cf4364659 Mart*0246        print *, prefix, 'nWetiGlobal ', (nWetiGlobal(k), k=1,Nr)
e189f4121c Mart*0247 #endif
b5856d1a66 Mart*0248 #ifdef ALLOW_OBCSN_CONTROL
65754df434 Mart*0249        do iobcs=1,nobcs
                0250         print *, prefix, 'nWetobcsnGlo (iobcs=', iobcs,')',
5cf4364659 Mart*0251      &       (nWetobcsnGlo(k,iobcs), k=1,Nr)
65754df434 Mart*0252        enddo
b5856d1a66 Mart*0253 #endif
                0254 #ifdef ALLOW_OBCSS_CONTROL
65754df434 Mart*0255        do iobcs=1,nobcs
                0256         print *, prefix, 'nWetobcssGlo (iobcs=', iobcs,')',
5cf4364659 Mart*0257      &         (nWetobcssGlo(k,iobcs), k=1,Nr)
65754df434 Mart*0258        enddo
b5856d1a66 Mart*0259 #endif
                0260 #ifdef ALLOW_OBCSW_CONTROL
65754df434 Mart*0261        do iobcs=1,nobcs
                0262         print *, prefix, 'nWetobcswGlo (iobcs=', iobcs,')',
5cf4364659 Mart*0263      &       (nWetobcswGlo(k,iobcs), k=1,Nr)
65754df434 Mart*0264        enddo
b5856d1a66 Mart*0265 #endif
                0266 #ifdef ALLOW_OBCSE_CONTROL
65754df434 Mart*0267        do iobcs=1,nobcs
                0268         print *, prefix, 'nWetobcseGlo (iobcs=', iobcs,')',
5cf4364659 Mart*0269      &       (nWetobcseGlo(k,iobcs), k=1,Nr)
65754df434 Mart*0270        enddo
b5856d1a66 Mart*0271 #endif
65754df434 Mart*0272        print *, prefix, 'ncvarindex ', (ncvarindex(i), i=1,maxcvars)
                0273        print *, prefix, 'ncvarrecs  ', (ncvarrecs(i),  i=1,maxcvars)
                0274        print *, prefix, 'ncvarxmax  ', (ncvarxmax(i),  i=1,maxcvars)
                0275        print *, prefix, 'ncvarymax  ', (ncvarymax(i),  i=1,maxcvars)
                0276        print *, prefix, 'ncvarnrmax ', (ncvarnrmax(i), i=1,maxcvars)
5cf4364659 Mart*0277        print *, prefix, 'ncvargrd   ', (ncvargrd(i), ',',i=1,maxcvars)
                0278        print *, prefix, 'ncvartype  ', (ncvartype(i),',',i=1,maxcvars)
65754df434 Mart*0279       end if
4cee17c1be Patr*0280 c--   Check the header information for consistency.
                0281 
                0282 cph      if ( filenopt .ne. nopt ) then
65754df434 Mart*0283 cph       print*
                0284 cph       print*,' READ_HEADER: Input data belong to the wrong'
                0285 cph       print*,'              optimization cycle.'
                0286 cph       print*,'              optimization cycle = ',nopt
                0287 cph       print*,'              input optim  cycle = ',filenopt
                0288 cph       print*
                0289 cph       stop   ' ... stopped in READ_HEADER.'
4cee17c1be Patr*0290 cph      endif
65754df434 Mart*0291 
4cee17c1be Patr*0292       if ( (fileiG .ne. biG) .or. (filejG .ne. bjG) ) then
65754df434 Mart*0293        print*
                0294        print*, prefix, 'Tile indices of loop and data do not match.'
                0295        print*,'                 loop x/y component = ',
                0296      &      biG,bjG
                0297        print*,'                 data x/y component = ',
                0298      &      fileiG,filejG
                0299        print*
                0300        stop   ' ... stopped in OPTIM_READDATA.'
4cee17c1be Patr*0301       endif
65754df434 Mart*0302 
5cf4364659 Mart*0303       if ( (filensx .ne. nSx) .or. (filensy .ne. nSy) ) then
65754df434 Mart*0304        print*
                0305        print*, prefix, ' Numbers of tiles do not match.'
                0306        print*,'                 parameter x/y no. of tiles = ',
5cf4364659 Mart*0307      &      nSx,nSy
65754df434 Mart*0308        print*,'                 data      x/y no. of tiles = ',
                0309      &      filensx,filensy
                0310        print*
                0311        stop   ' ... stopped in OPTIM_READDATA.'
4cee17c1be Patr*0312       endif
                0313 
                0314 ce    Add some more checks. ...
                0315 
                0316       if (.NOT. lheaderonly) then
                0317 c--   Read the data.
1a93ff8763 Mart*0318        icvoffset = 0
5cf4364659 Mart*0319 c      do ivar = 1,maxcvars
                0320        do ivar = 1,nvartype
                0321         if ( ncvarindex(ivar) .ne. -1 ) then
                0322          do icvrec = 1,ncvarrecs(ivar)
4e5349720c Patr*0323 cph          do bj = 1,nsy
                0324 cph           do bi = 1,nsx
5cf4364659 Mart*0325             read( funit ) ncvarindex_loc
1a93ff8763 Mart*0326             read( funit ) filej
                0327             read( funit ) filei
5cf4364659 Mart*0328             if ( ncvarindex_loc.NE.ncvarindex(ivar) .AND. icvrec.EQ.1 )
                0329      &       print*, prefix, ' ivar=', ivar, ' , ncvarindex(com,loc)=',
                0330      &               ncvarindex(ivar), ncvarindex_loc
                0331             do k = 1,ncvarnrmax(ivar)
1a93ff8763 Mart*0332              cbuffindex = 0
5cf4364659 Mart*0333              if (ncvargrd(ivar) .eq. 'c') then
1a93ff8763 Mart*0334               cbuffindex = nWetcGlobal(k)
5cf4364659 Mart*0335              else if (ncvargrd(ivar) .eq. 's') then
1a93ff8763 Mart*0336               cbuffindex = nWetsGlobal(k)
5cf4364659 Mart*0337              else if (ncvargrd(ivar) .eq. 'w') then
1a93ff8763 Mart*0338               cbuffindex = nWetwGlobal(k)
5cf4364659 Mart*0339 c            else if (ncvargrd(ivar) .eq. 'v') then
                0340 c             cbuffindex = nWetvGlobal(k)
e189f4121c Mart*0341 #ifdef ALLOW_SHIFWFLX_CONTROL
5cf4364659 Mart*0342              else if (ncvargrd(ivar) .eq. 'i') then
1a93ff8763 Mart*0343               cbuffindex = nWetiGlobal(k)
e189f4121c Mart*0344 #endif
65754df434 Mart*0345 #ifdef ALLOW_OBCS_CONTROL
4cee17c1be Patr*0346 cgg(   O.B. points have the grid mask "m".
5cf4364659 Mart*0347              else if (ncvargrd(ivar) .eq. 'm') then
4cee17c1be Patr*0348 cgg    From "icvrec", calculate what iobcs must be.
1a93ff8763 Mart*0349               gg   = (icvrec-1)/nobcs
                0350               igg  = int(gg)
                0351               iobcs= icvrec - igg*nobcs
65754df434 Mart*0352 # ifdef ALLOW_OBCSN_CONTROL
5cf4364659 Mart*0353               if (ncvarindex(ivar) .eq. 1)
                0354      &             cbuffindex = nWetobcsnGlo(k,iobcs)
65754df434 Mart*0355 # endif
                0356 # ifdef ALLOW_OBCSS_CONTROL
5cf4364659 Mart*0357               if (ncvarindex(ivar) .eq. 2)
                0358      &             cbuffindex = nWetobcssGlo(k,iobcs)
65754df434 Mart*0359 # endif
                0360 # ifdef ALLOW_OBCSE_CONTROL
5cf4364659 Mart*0361               if (ncvarindex(ivar) .eq. 3)
                0362      &             cbuffindex = nWetobcseGlo(k,iobcs)
                0363 # endif
                0364 # ifdef ALLOW_OBCSW_CONTROL
                0365               if (ncvarindex(ivar) .eq. 4)
                0366      &             cbuffindex = nWetobcswGlo(k,iobcs)
65754df434 Mart*0367 # endif
4cee17c1be Patr*0368 cgg)
65754df434 Mart*0369 #endif /* ALLOW_OBCS_CONTROL */
1a93ff8763 Mart*0370              endif
                0371              if ( icvoffset + cbuffindex .gt. nvarlength ) then
                0372               print*
                0373               print *, ' ERROR:'
                0374               print *, ' There are at least ', icvoffset+cbuffindex,
                0375      &             ' records in '//fname(1:28)//'.'
65754df434 Mart*0376               print *, ' This is more than expected from nvarlength =',
1a93ff8763 Mart*0377      &             nvarlength, '.'
                0378               print *, ' Something is wrong in the computation of '//
                0379      &             'the wet points or'
                0380               print *, ' in computing the number of records in '//
                0381      &             'some variable(s).'
                0382               print *, '  ...  stopped in OPTIM_READDATA.'
                0383               stop     '  ...  stopped in OPTIM_READDATA.'
                0384              endif
                0385              if (cbuffindex .gt. 0) then
                0386               read( funit ) cbuffindex
                0387               read( funit ) filek
                0388               read( funit ) (cbuff(ii), ii=1,cbuffindex)
                0389               do icvcomp = 1,cbuffindex
                0390                vv(icvoffset+icvcomp) = cbuff(icvcomp)
                0391 c     If you want to optimize with respect to just O.B. T and S
                0392 c     uncomment the next two lines.
                0393 c              if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0.
                0394 c              if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0.
                0395               enddo
                0396               icvoffset = icvoffset + cbuffindex
                0397              endif
                0398             enddo
4e5349720c Patr*0399 cph           enddo
                0400 cph          enddo
4cee17c1be Patr*0401          enddo
1a93ff8763 Mart*0402         endif
                0403        enddo
65754df434 Mart*0404 
4cee17c1be Patr*0405       else
                0406 
                0407 c--   Assign the number of control variables.
1a93ff8763 Mart*0408        nn = nvarlength
65754df434 Mart*0409 
4cee17c1be Patr*0410       endif
                0411 
                0412       close( funit )
                0413 
                0414 c--   Assign the cost function value in case we read the cost file.
                0415 
                0416       if      ( dfile .eq. ctrlname ) then
1a93ff8763 Mart*0417        ff = 0. d 0
4cee17c1be Patr*0418       else if ( dfile .eq. costname ) then
1a93ff8763 Mart*0419        ff = fileff
4cee17c1be Patr*0420       endif
65754df434 Mart*0421 c--   Always return the cost function value if lheaderonly
                0422       if ( lheaderonly) ff = fileff
                0423 
                0424       print *, prefix, 'end of optim_readdata'
                0425       print *, ' '
4cee17c1be Patr*0426 
                0427       return
                0428       end
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!