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
      real*4 cbuff( sNx*nSx*nPx*sNy*nSy*nPy )

      character*(128) fname

c      integer         filei
c      integer         filej
c      integer         filek
c      integer         fileig
c      integer         filejg
c      integer         filensx
c      integer         filensy
      integer         filenopt
      _RL             fileff

cgg(
      _RL     gg
      integer igg
      integer iobcs
cgg)

c     == end of interface ==

      print *, 'pathei-lsopt in optim_readdata'

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

      read( funit ) (nWetcGlobal(k), k=1,nr)
      read( funit ) (nWetsGlobal(k), k=1,nr)
      read( funit ) (nWetwGlobal(k), k=1,nr)
#ifdef ALLOW_CTRL_WETV
      read( funit ) (nWetvGlobal(k), k=1,nr)
#endif

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