subroutine optim_readdata( I nn, I dfile, I lheaderonly, O ff, O vv & ) c ================================================================== c SUBROUTINE optim_readdata c ================================================================== c c o Read the data written by the MITgcmUV state estimation setup and c join them to one vector that is subsequently used by the minimi- c zation algorithm "lsopt". Depending on the specified file name c either the control vector or the gradient vector can be read. c c *dfile* should be the radix of the file: ecco_ctrl or ecco_cost c c started: Christian Eckert eckert@mit.edu 12-Apr-2000 c c changed: Patrick Heimbach heimbach@mit.edu 19-Jun-2000 c - finished, revised and debugged c c ================================================================== c SUBROUTINE optim_readdata c ================================================================== implicit none c == global variables == #include "EEPARAMS.h" #include "SIZE.h" cgg Include ECCO_CPPOPTIONS because the ecco_ctrl,cost files cgg have headers with options for OBCS masks. #include "ECCO_CPPOPTIONS.h" #include "ecco.h" #include "ctrl.h" #include "optim.h" #include "minimization.h" c == routine arguments == integer nn _RL ff #if defined (DYNAMIC) _RL vv(nn) #elif defined (USE_POINTER) || (MAX_INDEPEND == 0) _RL vv pointer (pvv,vv(1)) #else integer nmax parameter( nmax = MAX_INDEPEND ) _RL vv(nmax) #endif character*(9) dfile logical lheaderonly c == local variables == integer bi,bj integer biG,bjG integer i,j integer ii,k integer icvar integer icvrec integer icvcomp integer icvoffset integer nopt integer funit integer cbuffindex _RL cbuff( sNx*nSx*nPx*sNy*nSy*nPy ) character*(128) fname integer filei integer filej integer filek integer filenopt integer fileig integer filejg integer filensx integer filensy _RL fileff cgg( _RL gg integer igg integer iobcs cgg) c == end of interface == c-- The reference i/o unit. funit = 20 c-- Next optimization cycle. nopt = optimcycle if ( dfile .eq. ctrlname ) then print* print*,' OPTIM_READDATA: Reading control vector' print*,' for optimization cycle: ',nopt print* else if ( dfile .eq. costname ) then print* print*,' OPTIM_READDATA: Reading cost function and' print*,' gradient of cost function' print*,' for optimization cycle: ',nopt print* else print* print*,' OPTIM_READDATA: subroutine called by a false *dfile*' print*,' argument. *dfile* = ',dfile print* stop ' ... stopped in OPTIM_READDATA.' endif c-- Read the data. bjG = 1 + (myygloballo - 1)/sny biG = 1 + (myxgloballo - 1)/snx c-- Generate file name and open the file. write(fname(1:128),'(4a,i4.4)') & dfile,'_',expId(1:10),'.opt', nopt open( funit, file = fname, & status = 'old', & form = 'unformatted', & access = 'sequential' ) print*, 'opened file ', fname c-- Read the header. read( funit ) nvartype read( funit ) nvarlength read( funit ) expId read( funit ) filenopt read( funit ) fileff read( funit ) fileiG read( funit ) filejG read( funit ) filensx read( funit ) filensy cph( print *,'ph-opt 1 ', nvartype, nvarlength, filensx, filensy cph) read( funit ) (nWetcGlobal(k), k=1,nr) read( funit ) (nWetsGlobal(k), k=1,nr) read( funit ) (nWetwGlobal(k), k=1,nr) read( funit ) (nWetvGlobal(k), k=1,nr) cgg( Add OBCS Mask information into the header section for optimization. #ifdef ALLOW_OBCSN_CONTROL read( funit ) ((nWetobcsnGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs) #endif #ifdef ALLOW_OBCSS_CONTROL read( funit ) ((nWetobcssGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs) #endif #ifdef ALLOW_OBCSW_CONTROL read( funit ) ((nWetobcswGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs) #endif #ifdef ALLOW_OBCSE_CONTROL read( funit ) ((nWetobcseGlo(k,iobcs), k=1,nr),iobcs= 1,nobcs) #endif cgg) read( funit ) (ncvarindex(i), i=1,maxcvars) read( funit ) (ncvarrecs(i), i=1,maxcvars) read( funit ) (ncvarxmax(i), i=1,maxcvars) read( funit ) (ncvarymax(i), i=1,maxcvars) read( funit ) (ncvarnrmax(i), i=1,maxcvars) read( funit ) (ncvargrd(i), i=1,maxcvars) read( funit ) cph( cph if (lheaderonly) then print *, 'pathei: nvartype ', nvartype print *, 'pathei: nvarlength ', nvarlength print *, 'pathei: expId ', expId print *, 'pathei: filenopt ', filenopt print *, 'pathei: fileff ', fileff print *, 'pathei: fileiG ', fileiG print *, 'pathei: filejG ', filejG print *, 'pathei: filensx ', filensx print *, 'pathei: filensy ', filensy print *, 'pathei: nWetcGlobal ', & (nWetcGlobal(k), k=1,nr) print *, 'pathei: nWetsGlobal ', & (nWetsGlobal(k), k=1,nr) print *, 'pathei: nWetwGlobal ', & (nWetwGlobal(k), k=1,nr) print *, 'pathei: nWetvGlobal ', & (nWetvGlobal(k), k=1,nr) print *, 'pathei: ncvarindex ', & (ncvarindex(i), i=1,maxcvars) print *, 'pathei: ncvarrecs ', & (ncvarrecs(i), i=1,maxcvars) print *, 'pathei: ncvarxmax ', & (ncvarxmax(i), i=1,maxcvars) print *, 'pathei: ncvarymax ', & (ncvarymax(i), i=1,maxcvars) print *, 'pathei: ncvarnrmax ', & (ncvarnrmax(i), i=1,maxcvars) print *, 'pathei: ncvargrd ', & (ncvargrd(i), i=1,maxcvars) cph end if cph) c-- Check the header information for consistency. cph if ( filenopt .ne. nopt ) then cph print* cph print*,' READ_HEADER: Input data belong to the wrong' cph print*,' optimization cycle.' cph print*,' optimization cycle = ',nopt cph print*,' input optim cycle = ',filenopt cph print* cph stop ' ... stopped in READ_HEADER.' cph endif if ( (fileiG .ne. biG) .or. (filejG .ne. bjG) ) then print* print*,' READ_HEADER: Tile indices of loop and data ' print*,' do not match.' print*,' loop x/y component = ', & biG,bjG print*,' data x/y component = ', & fileiG,filejG print* stop ' ... stopped in READ_HEADER.' endif if ( (filensx .ne. nsx) .or. (filensy .ne. nsy) ) then print* print*,' READ_HEADER: Numbers of tiles do not match.' print*,' parameter x/y no. of tiles = ', & bi,bj print*,' data x/y no. of tiles = ', & filensx,filensy print* stop ' ... stopped in READ_HEADER.' endif ce Add some more checks. ... if (.NOT. lheaderonly) then c-- Read the data. icvoffset = 0 do icvar = 1,maxcvars if ( ncvarindex(icvar) .ne. -1 ) then do icvrec = 1,ncvarrecs(icvar) do bj = 1,nsy do bi = 1,nsx read( funit ) ncvarindex(icvar) read( funit ) filej read( funit ) filei do k = 1,ncvarnrmax(icvar) cbuffindex = 0 if (ncvargrd(icvar) .eq. 'c') then cbuffindex = nWetcGlobal(k) else if (ncvargrd(icvar) .eq. 's') then cbuffindex = nWetsGlobal(k) else if (ncvargrd(icvar) .eq. 'w') then cbuffindex = nWetwGlobal(k) else if (ncvargrd(icvar) .eq. 'v') then cbuffindex = nWetvGlobal(k) cgg( O.B. points have the grid mask "m". else if (ncvargrd(icvar) .eq. 'm') then cgg From "icvrec", calculate what iobcs must be. gg = (icvrec-1)/nobcs igg = int(gg) iobcs= icvrec - igg*nobcs #ifdef ALLOW_OBCSN_CONTROL if (icvar .eq. 11) then cbuffindex = nWetobcsnGlo(k,iobcs) endif #endif #ifdef ALLOW_OBCSS_CONTROL if (icvar .eq. 12) then cbuffindex = nWetobcssGlo(k,iobcs) endif #endif #ifdef ALLOW_OBCSW_CONTROL if (icvar .eq. 13) then cbuffindex = nWetobcswGlo(k,iobcs) endif #endif #ifdef ALLOW_OBCSE_CONTROL if (icvar .eq. 14) then cbuffindex = nWetobcseGlo(k,iobcs) endif #endif cgg) endif if (cbuffindex .gt. 0) then read( funit ) cbuffindex read( funit ) filek read( funit ) (cbuff(ii), ii=1,cbuffindex) do icvcomp = 1,cbuffindex vv(icvoffset+icvcomp) = cbuff(icvcomp) cgg( Right now, the changes to the open boundary velocities are not balanced. cgg( The model will crash due to physical reasons. cgg( However, we can optimize with respect to just O.B. T and S if the cgg( next two lines are uncommented. cgg if (iobcs .eq. 3) vv(icvoffset+icvcomp)=0. cgg if (iobcs .eq. 4) vv(icvoffset+icvcomp)=0. enddo icvoffset = icvoffset + cbuffindex endif enddo enddo enddo enddo endif enddo else c-- Assign the number of control variables. nn = nvarlength endif close( funit ) c-- Assign the cost function value in case we read the cost file. if ( dfile .eq. ctrlname ) then ff = 0. d 0 else if ( dfile .eq. costname ) then ff = fileff endif return end