--- MITgcm/optim/optim_readdata.F 2002/02/05 20:34:35 1.1 +++ MITgcm/optim/optim_readdata.F 2002/02/05 20:34:35 1.1.2.1 @@ -0,0 +1,282 @@ + + 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" + +#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 + +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 ) (((nWetcTile(i,j,k), i=1,nsx), j=1,nsy), + & k=1,nr) + read( funit ) (((nWetsTile(i,j,k), i=1,nsx), j=1,nsy), + & k=1,nr) + read( funit ) (((nWetwTile(i,j,k), i=1,nsx), j=1,nsy), + & k=1,nr) + 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: nWetcTile ', + & (((nWetcTile(i,j,k), i=1,nsx), j=1,nsy), k=1,nr) + print *, 'pathei: nWetsTile ', + & (((nWetsTile(i,j,k), i=1,nsx), j=1,nsy), k=1,nr) + print *, 'pathei: nWetwTile ', + & (((nWetwTile(i,j,k), i=1,nsx), j=1,nsy), 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 = nwetctile(bi,bj,k) + else if (ncvargrd(icvar) .eq. 's') then + cbuffindex = nwetstile(bi,bj,k) + else if (ncvargrd(icvar) .eq. 'w') then + cbuffindex = nwetwtile(bi,bj,k) + 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) + 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