#include "FIZHI_OPTIONS.h"
       subroutine update_ocean_exports (myTime, myIter, myThid)
c----------------------------------------------------------------------
c  Subroutine update_ocean_exports - 'Wrapper' routine to update
c        the fields related to the ocean surface that are needed
c        by fizhi (sst and sea ice extent).
c
c Call:  getsst  (Return the current sst field-read dataset if needed)
c        getsice (Return the current sea ice field-read data if needed)
c-----------------------------------------------------------------------
       implicit none
#include "SIZE.h"
#include "GRID.h"
#include "fizhi_ocean_coms.h"
#include "EEPARAMS.h"
#include "chronos.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_SIZE.h"
#include "W2_EXCH2_TOPOLOGY.h"
#endif /* ALLOW_EXCH2 */

       integer myIter, myThid
       _RL myTime

       INTEGER xySize
#if defined(ALLOW_EXCH2)
       PARAMETER ( xySize = W2_ioBufferSize )
#else
       PARAMETER ( xySize = Nx*Ny )
#endif
       integer i, j, bi, bj, bislot, bjslot
       integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
       integer xsize, ysize
       _RL        sstmin
       parameter ( sstmin = 273.16 )

       _RL sst1 (xySize), sst2 (xySize)
       _RL sice1(xySize), sice2(xySize)
c      _RL sst1(xsize,ysize),sst2(xsize,ysize)
c      _RL sice1(xsize,ysize),sice2(xsize,ysize)
       integer nymd1sst(nSx,nSy),nymd2sst(nSx,nSy)
       integer nymd1sice(nSx,nSy),nymd2sice(nSx,nSy)
       integer nhms1sst(nSx,nSy),nhms2sst(nSx,nSy)
       integer nhms1sice(nSx,nSy),nhms2sice(nSx,nSy)
       integer sstdates(370,nSx,nSy),sicedates(370,nSx,nSy)
       integer ssttimes(370,nSx,nSy),sicetimes(370,nSx,nSy)
       logical first(nSx,nSy)
       integer nSxnSy
       parameter(nSxnSy = nSx*nSy)
       data first/nSxnSy*.true./

       save nymd1sst,nymd2sst,nymd1sice,nymd2sice
       save nhms1sst,nhms2sst,nhms1sice,nhms2sice
       save sst1, sst2, sice1, sice2
       save sstdates, sicedates
       save ssttimes, sicetimes

#if defined(ALLOW_EXCH2)
       xsize = exch2_global_Nx
       ysize = exch2_global_Ny
#else
       xsize = Nx
       ysize = Ny
#endif
       idim1 = 1-OLx
       idim2 = sNx+OLx
       jdim1 = 1-OLy
       jdim2 = sNy+OLy
       im1 = 1
       im2 = sNx
       jm1 = 1
       jm2 = sNy

C***********************************************************************

       DO BJ = myByLo(myThid),myByHi(myThid)
       DO BI = myBxLo(myThid),myBxHi(myThid)
#if defined(ALLOW_EXCH2)
       bislot = exch2_txglobalo(W2_myTileList(bi,bj))-1
       bjslot = exch2_tyglobalo(W2_myTileList(bi,bj))-1
#else
       bislot = myXGlobalLo-1+(bi-1)*sNx
       bjslot = myYGlobalLo-1+(bj-1)*sNy
#endif

       call getsst(ksst,sstclim,idim1,idim2,jdim1,jdim2,im1,im2,
     .  jm1,jm2,nSx,nSy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
     .  sst1,sst2,first(bi,bj),nymd1sst(bi,bj),nymd2sst(bi,bj),
     .  nhms1sst(bi,bj),nhms2sst(bi,bj),sstdates(1,bi,bj),
     .  ssttimes(1,bi,bj),sst,myThid)
       call getsice(kice,siceclim,idim1,idim2,jdim1,jdim2,im1,im2,
     .  jm1,jm2,nSx,nSy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
     .  sice1,sice2,first(bi,bj),nymd1sice(bi,bj),nymd2sice(bi,bj),
     .  nhms1sice(bi,bj),nhms2sice(bi,bj),sicedates(1,bi,bj),
     .  sicetimes(1,bi,bj),sice,myThid)

c Check for Minimum Open-Water SST
c --------------------------------
       do j=jm1,jm2
       do i=im1,im2
       if(sice(i,j,bi,bj).eq.0.0 .and. sst(i,j,bi,bj).lt.sstmin)
     .                                          sst(i,j,bi,bj) = sstmin
       enddo
       enddo

       ENDDO
       ENDDO
       _EXCH_XY_RL(sst,myThid)
       _EXCH_XY_RL(sice,myThid)

       return
       end

       subroutine getsice(iunit,clim,idim1,idim2,jdim1,jdim2,im1,im2,
     .   jm1,jm2,nSumx,nSumy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
     .   sicebc1,sicebc2,first,nymdbc1,nymdbc2,nhmsbc1,nhmsbc2,
     .   nymdbc,nhmsbc,sice,mythid)
C***********************************************************************
C
C!ROUTINE: GETSICE
C!DESCRIPTION:  GETSICE returns the sea ice depth.
C!              This routine is adaptable for any frequency
C!              data upto a daily frequency.
C!              note: for diurnal data ndmax should be increased.
C
C!INPUT PARAMETERS:
C!      iunit     Unit number assigned to the sice data file
C!      idim1     Start dimension in x-direction
C!      idim2     End dimension in x-direction
C!      jdim1     Start dimension in y-direction
C!      jdim2     End dimension in y-direction
C!      im1       Begin of x-direction span for filling sice
C!      im2       End of x-direction span for filling sice
C!      jm1       Begin of y-direction span for filling sice
C!      jm2       End of y-direction span for filling sice
C!      nSumx     Number of processors in x-direction (local processor)
C!      nSumy     Number of processors in y-direction (local processor)
C!      xsize      Number of processors in x-direction (global)
C!      ysize      Number of processors in y-direction (global)
C!      bi        Processor number in x-direction (local to processor)
C!      bj        Processor number in y-direction (local to processor)
C!      bislot  Processor number in x-direction (global)
C!      bjslot  Processor number in y-direction (global)
C!      nymd    YYMMDD of the current model timestep
C!      nhms    HHMMSS of the model time
C
C!OUTPUT PARAMETERS:
C!      sice(idim1:idim2,jdim1:jdim2,nSumx,nSumy) Sea ice depth in meters
C
C!ROUTINES CALLED:
C
C!      bcdata       Reads the data for a given unit number
C!      bcheader     Reads the header info for a given unit number
C!      interp_time  Returns weights for linear interpolation
C
C--------------------------------------------------------------------------

      implicit none
#include "SIZE.h"

      integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSumx,nSumy
      integer xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,mythid
      logical clim

      _RL sicebc1(xsize,ysize)
      _RL sicebc2(xsize,ysize)
      _RL sice(idim1:idim2,jdim1:jdim2,nSumx,nSumy)
      integer nhmsbc1,nhmsbc2,nymdbc1,nymdbc2
      logical first

C Maximum number of dates in one year for the data
      integer ndmax
      parameter (ndmax = 370)
      integer nymdbc(ndmax),nhmsbc(ndmax)

      character*8  cname
      character*80 cdscrip
      character*22 sicedata
      _RL fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
      logical found, error
      integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc
      integer ndatebc,nrec
      integer nymdmod


C--------- Variable Initialization ---------------------------------

      data error /.false./

c  save header info
      save imbc,jmbc,lat0,lon0,ndatebc,undef

c  this only works for between 1950-2050
      if (nymd .lt. 500101) then
        nymdmod = 20000000 + nymd
      else if (nymd .le. 991231) then
        nymdmod = 19000000 + nymd
      else
        nymdmod = nymd
      endif

      iyear   = nymdmod/10000
      if(clim) then
       if(xsize.eq.192)sicedata='sice19232.weekly.clim'
       if(xsize.eq.612)sicedata='sice612102.weekly.clim'
      else
       if(xsize.eq.192)
     .           WRITE(sicedata,'(A,I4)')'sice19232.weekly.y',iyear
       if(xsize.eq.612)
     .           WRITE(sicedata,'(A,I4)')'sice612102.weekly.y',iyear
      endif

c  initialize so that first time through they have values for the check
c  these values make the iyear .ne. iyearbc true anyways for
c  for the first time so first isnt checked below.

      if (first) then
        nymdbc(2) = 0
        nymdbc1   = 0
        nymdbc2   = 0
        nhmsbc1   = 0
        nhmsbc2   = 0
        first = .false.
      endif

C---------- Read in Header file ----------------------------------

      iyearbc = nymdbc(2)/10000

      if( iyear.ne.iyearbc ) then

       close(iunit)
       open (iunit,file=sicedata,form='unformatted',access='direct',
     .                                         recl=xsize*ysize*4)
       nrec = 1
       call bcheader (iunit, ndmax, nrec,
     .          cname, cdscrip, imbc, jmbc, lat0, lon0,
     .          ndatebc, nymdbc, nhmsbc, undef, error)

C--------- Check data for Compatibility ------------------------------

C Check for correct data in boundary condition file
       if (.not.error .and. cname.ne.'SICE') then
        write(6,*)'Wrong data in SICE boundary condition file => ',cname
        error = .true.
       endif

C Check Horizontal Resolution
       if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
        write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
        write(6,*) ' B.C. Resolution:  ',imbc*jmbc
        write(6,*) 'Model Resolution:  ',xsize*ysize
        error = .true.
       endif

C Check Year
       iyearbc = nymdbc(2)/10000
       if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
        write(6,*)'     B.C. Year DOES NOT match REQUESTED Year!'
        write(6,*)'     B.C. Year:  ', iyearbc
        write(6,*)'Requested Year:  ', iyear
        error = .true.
       endif

       if (.not.error)   then
C if climatology, fill dates for data with current model year
        if (iyearbc.eq.0) then
         write(6,*)
         write(6,*) 'Climatological Dataset is being used.'
         write(6,*) 'Current model year to be used to fill Header Dates'
         do n = 2, ndatebc-1
          nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
         enddo
C  For the first date subtract 1 year from the current model NYMD
         n = 1
         nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
C  For the last date add 1 year to the current model NYMD
         n = ndatebc
         nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
        endif

C  Write out header info
        _BEGIN_MASTER( myThid )
        write(6,*) ' Updated boundary condition data'
        write(6,*) ' ---------------------------------'
        write(6,*) ' Variable: ',cname
        write(6,*) ' Description: ',cdscrip
        write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
     .                                       ' Undefined value = ',undef
        write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
        write(6,*) ' Data valid at these times: '
        ndby3 = ndatebc/3
        do n = 1, ndby3*3,3
         write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
 1000    format(3(2x,i3,':',i8,2x,i8))
        enddo
        write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
        _END_MASTER( myThid )
       endif

      endif

C---------- Read sice data if necessary -------------------------------

      found = .false.
      nd = 2

c  If model time is not within the times of saved sice data
c  from previous call to getsice then read new data

      timemod = dfloat(nymdmod) + dfloat(nhms)   /1000000
      timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
      timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000

      if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then

       do while (.not.found .and. nd .le. ndatebc)
        timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
        if (timebc2 .gt. timemod) then
         nymdbc1 = nymdbc(nd-1)
         nymdbc2 = nymdbc(nd)
         nhmsbc1 = nhmsbc(nd-1)
         nhmsbc2 = nhmsbc(nd)
         call bcdata (iunit,imbc,jmbc,nd,nd+1,sicebc1,sicebc2)
         found = .true.
        else
         nd = nd + 1
        endif
       enddo

c  Otherwise the data from the last time in getsice surrounds the
c  current model time.

      else
       found = .true.
      endif

      if (.not.found) then
       print *, 'STOP: Could not find SICE dates for model time.'
       call my_finalize
       call my_exit (101)
      endif

C---------- Interpolate sice data ------------------------------------

      call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
     .                                                       fac1,fac2)

      do j = jm1,jm2
      do i = im1,im2
       sice(i,j,bi,bj) = sicebc1(i+bislot,j+bjslot)*fac1
     .                 + sicebc2(i+bislot,j+bjslot)*fac2
c average to 0 or 1
c -----------------
       if (sice(i,j,bi,bj) .ge. 0.5) then
        sice(i,j,bi,bj) = 1.
       else
        sice(i,j,bi,bj) = 0.
       endif
      enddo
      enddo

C---------- Fill sice with depth of ice ------------------------------------
      do j = jm1,jm2
      do i = im1,im2
       if (sice(i,j,bi,bj) .eq. 1.) then
        sice(i,j,bi,bj) = 3.
       endif
      enddo
      enddo
C---------------------------------------------------------------------------

      return
      end
      subroutine getsst(iunit,clim,idim1,idim2,jdim1,jdim2,im1,im2,
     .   jm1,jm2,nSumx,nSumy,xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,
     .   sstbc1,sstbc2,first,nymdbc1,nymdbc2,nhmsbc1,nhmsbc2,
     .   nymdbc,nhmsbc,sst,mythid)
C***********************************************************************
C
C!ROUTINE: GETSST
C!DESCRIPTION:  GETSST gets the SST data.
C!              This routine is adaptable for any frequency
C!              data upto a daily frequency.
C!              note: for diurnal data ndmax should be increased.
C
C!INPUT PARAMETERS:
C!      iunit     Unit number assigned to the sice data file
C!      idim1     Start dimension in x-direction
C!      idim2     End dimension in x-direction
C!      jdim1     Start dimension in y-direction
C!      jdim2     End dimension in y-direction
C!      im1       Begin of x-direction span for filling sst
C!      im2       End of x-direction span for filling sst
C!      jm1       Begin of y-direction span for filling sst
C!      jm2       End of y-direction span for filling sst
C!      nSumx     Number of processors in x-direction (local processor)
C!      nSumy     Number of processors in y-direction (local processor)
C!      xsize     x-dimension of global array
C!      ysize     y-dimension of global array
C!      bi        Processor number in x-direction (local to processor)
C!      bj        Processor number in y-direction (local to processor)
C!      bislot    Slot number into global array in x-direction (global)
C!      bjslot    Slot number into global array in y-direction (global)
C!      nymd      YYMMDD of the current model timestep
C!      nhms      HHMMSS of the model time
C
C!OUTPUT PARAMETERS:
C!     sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy) Sea surface temperature (K)
C
C!ROUTINES CALLED:
C
C!     bcdata     Reads the data for a given unit number
C!     bcheader   Reads the header info for a given unit number
C!     interp_time  Returns weights for linear interpolation
C
C--------------------------------------------------------------------------

      implicit none
#include "SIZE.h"

      integer iunit,idim1,idim2,jdim1,jdim2,im1,im2,jm1,jm2,nSumx,nSumy
      integer xsize,ysize,bi,bj,bislot,bjslot,nymd,nhms,mythid
      logical clim

      _RL sstbc1(xsize,ysize)
      _RL sstbc2(xsize,ysize)
      _RL sst(idim1:idim2,jdim1:jdim2,nSumx,nSumy)
      integer nhmsbc1,nhmsbc2,nymdbc1,nymdbc2
      logical first

C Maximum number of dates in one year for the data
      integer ndmax
      parameter (ndmax = 370)
      integer nymdbc(ndmax),nhmsbc(ndmax)

      character*8  cname
      character*80 cdscrip
      character*21 sstdata
      _RL fac1, fac2, lat0, lon0, timebc1, timebc2, timemod, undef
      logical found, error
      integer i,j,n,nn,iyear,iyearbc,nd,ndby3,imbc,jmbc
      integer ndatebc,nrec
      integer nymdmod


C--------- Variable Initialization ---------------------------------

      data error /.false./

c  save header info
      save imbc,jmbc,lat0,lon0,ndatebc,undef

c  this only works for between 1950-2050
      if (nymd .lt. 500101) then
        nymdmod = 20000000 + nymd
      else if (nymd .le. 991231) then
        nymdmod = 19000000 + nymd
      else
        nymdmod = nymd
      endif

      iyear   = nymdmod/10000
      if(clim) then
       if(xsize.eq.192)sstdata='sst19232.weekly.clim'
       if(xsize.eq.612)sstdata='sst612102.weekly.clim'
      else
       if(xsize.eq.192)
     .           WRITE(sstdata,'(A,I4)')'sst19232.weekly.y',iyear
       if(xsize.eq.612)
     .           WRITE(sstdata,'(A,I4)')'sst612102.weekly.y',iyear
      endif

c  initialize so that first time through they have values for the check
c  these vaules make the iyear .ne. iyearbc true anyways for
c  for the first time so first isnt checked below.
      if (first) then
        nymdbc(2) = 0
        nymdbc1   = 0
        nymdbc2   = 0
        nhmsbc1   = 0
        nhmsbc2   = 0
        first = .false.
      endif

C---------- Read in Header file ----------------------------------

      iyearbc = nymdbc(2)/10000

      if( iyear.ne.iyearbc ) then

       close(iunit)
       open (iunit,file=sstdata,form='unformatted',access='direct',
     .                                        recl=xsize*ysize*4)
       nrec = 1
       call bcheader (iunit, ndmax, nrec,
     .          cname, cdscrip, imbc, jmbc, lat0, lon0,
     .          ndatebc, nymdbc, nhmsbc, undef, error)

C--------- Check data for Compatibility

C Check for correct data in boundary condition file
       if (.not.error .and. cname.ne.'SST') then
        write(6,*)'Wrong data in SST boundary condition file => ',cname
        error = .true.
       endif

C Check Horizontal Resolution
       if(.not.error.and.imbc*jmbc.ne.xsize*ysize)then
        write(6,*) ' B.C. Resolution DOES NOT match Model Resolution!'
        write(6,*) ' B.C. Resolution:  ',imbc*jmbc
        write(6,*) 'Model Resolution:  ',xsize*ysize
        error = .true.
       endif

C Check Year
       iyearbc = nymdbc(2)/10000
       if (.not.error .and. iyear.ne.iyearbc .and. iyearbc.ne.0) then
        write(6,*)'     B.C. Year DOES NOT match REQUESTED Year!'
        write(6,*)'     B.C. Year:  ', iyearbc
        write(6,*)'Requested Year:  ', iyear
        error = .true.
       endif

       if (.not.error)   then
C if climatology, fill dates for data with current model year
        if (iyearbc.eq.0) then
         write(6,*)
         write(6,*)'Climatological Dataset is being used.'
         write(6,*)'Current model year is used to fill Header Dates'
         do n = 2, ndatebc-1
          nymdbc(n) = nymdbc(n) +(nymdmod/10000)*10000
         enddo
C  For the first date subtract 1 year from the current model NYMD
         n = 1
         nymdbc(n) = nymdbc(n) +(nymdmod/10000-1)*10000
C  For the last date add 1 year to the current model NYMD
         n = ndatebc
         nymdbc(n) = nymdbc(n) +(nymdmod/10000+1)*10000
        endif

C  Write out header info
        _BEGIN_MASTER( myThid )
        write(6,*) ' Updated boundary condition data'
        write(6,*) ' ---------------------------------'
        write(6,*) ' Variable: ',cname
        write(6,*) ' Description: ',cdscrip
        write(6,*) ' Resolution: x= ',imbc,' y= ',jmbc,
     .                                       ' Undefined value = ',undef
        write(6,*) ' Starting latitude = ',lat0,' longitude =',lon0
        write(6,*) ' Data valid at these times: '
        ndby3 = ndatebc/3
        do n = 1, ndby3*3,3
         write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=n,n+2)
 1000    format(3(2x,i3,':',i8,2x,i8))
        enddo
        write(6,1000) (nn,nymdbc(nn),nhmsbc(nn),nn=ndby3*3+1,ndatebc)
        _END_MASTER( myThid )
       endif

       if( error ) call my_exit (101)

      endif

C---------- Read SST data if necessary -------------------------------

      found = .false.
      nd = 2

c  If model time is not within the times of saved sst data
c  from previous call to getsst then read new data

      timemod = dfloat(nymdmod) + dfloat(nhms)   /1000000
      timebc1 = dfloat(nymdbc1) + dfloat(nhmsbc1)/1000000
      timebc2 = dfloat(nymdbc2) + dfloat(nhmsbc2)/1000000

      if (timemod .lt. timebc1 .or. timemod .ge. timebc2) then

       do while (.not.found .and. nd .le. ndatebc)
        timebc2 = dfloat(nymdbc(nd)) + dfloat(nhmsbc(nd))/1000000
        if (timebc2 .gt. timemod) then
         nymdbc1 = nymdbc(nd-1)
         nymdbc2 = nymdbc(nd)
         nhmsbc1 = nhmsbc(nd-1)
         nhmsbc2 = nhmsbc(nd)
         call bcdata (iunit,imbc,jmbc,nd,nd+1,sstbc1,sstbc2)
         found = .true.
        else
         nd = nd + 1
        endif
       enddo

c  Otherwise the data from the last time in getsst surrounds the
c  current model time.

      else
       found = .true.
      endif

      if (.not.found) then
       print *, 'STOP: Could not find SST dates for model time.'
       call my_finalize
       call my_exit (101)
      endif

C---------- Interpolate SST data ------------------------------------

      call interp_time(nymdmod,nhms,nymdbc1,nhmsbc1,nymdbc2,nhmsbc2,
     .                                                        fac1,fac2)

      do j = jm1,jm2
      do i = im1,im2
       sst(i,j,bi,bj) = sstbc1(i+bislot,j+bjslot)*fac1
     .                + sstbc2(i+bislot,j+bjslot)*fac2
      enddo
      enddo


      return
      end

      subroutine bcdata (iunit,im,jm,nrec1,nrec2,field1,field2)
C************************************************************************
C
C!ROUTINE:      BCDATA
C!DESCRIPTION:  BCDATA reads the data from the file assigned to the
C!              passed unit number and returns data from the two times
C!              surrounding the current model time.  The two record
C!              postitions are not assumed to be next to each other.
C
C!INPUT PARAMETERS:
C!      im      number of x points
C!      im      number of x points
C!      nrec1   record number of the time before the model time
C!      nrec2   record number of the time after the model time
C
C!OUTPUT PARAMETERS:
C!      field1(im,jm)  data field before the model time
C!      field2(im,jm)  data field after the model time
C
C--------------------------------------------------------------------------
      implicit none

      integer iunit,im,jm,nrec1,nrec2

      _RL  field1(im,jm)
      _RL  field2(im,jm)

      integer i,j
      real*4 f1(im,jm), f2(im,jm)

C--------- Read file -----------------------------------------------
      read(iunit,rec=nrec1) f1
      read(iunit,rec=nrec2) f2

#ifdef _BYTESWAPIO
      call MDS_BYTESWAPR4( im*jm, f1)
      call MDS_BYTESWAPR4( im*jm, f2)
#endif
      do j=1,jm
      do i=1,im
       field1(i,j) = f1(i,j)
       field2(i,j) = f2(i,j)
      enddo
      enddo

      return
      end
      subroutine bcheader (iunit, ndmax, nrec,
     .           cname, cdscrip, im, jm, lat0, lon0, ndatebc,
     .           nymdbc, nhmsbc, undef, error)
C************************************************************************
C
C!ROUTINE:     BCHEADER
C!DESCRIPTION: BCHEADER reads the header from a file and returns the info.
C
C!INPUT PARAMETERS:
C!      iunit    unit number assigned to the data file
C!      ndmax    maximum number of date/times of the data
C!      nrec     record number of the header info (or assume 1 ?)
C
C!OUTPUT PARAMETERS:
C!      cname         name of the data in the file header
C!      cdscrip       description of the data in the file header
C!      im            number of x points
C!      jm            number of y points
C!      lat0          starting latitude for the data grid
C!      lon0          starting longitude for the data grid
C!      ndatebc       number of date/times of the data in the file
C!      nymdbc(ndmax) array of dates for the data including century
C!      nhmsbc(ndmax) array of times for the data
C!      undef         value for undefined values in the data
C!      error         logical TRUE if dataset problem
C
C--------------------------------------------------------------------------
      implicit none

      integer iunit, ndmax, nrec

      character*8  cname
      character*80 cdscrip
      character*112 dummy112
      integer im,jm,ndatebc,nymdbc(ndmax),nhmsbc(ndmax)
      _RL lat0,lon0,undef
      logical error

      integer i
      integer*4 im_32,jm_32
      integer*4 ndatebc_32,nhmsbc_32(ndmax),nymdbc_32(ndmax)
      real*4 lat0_32,lon0_32,undef_32

C--------- Read file -----------------------------------------------

      read(iunit,rec=nrec,err=500) cname, cdscrip,
     .     im_32, jm_32, lat0_32, lon0_32,
     .     ndatebc_32, undef_32

#ifdef _BYTESWAPIO
      call MDS_BYTESWAPI4( 1, im_32)
      call MDS_BYTESWAPI4( 1, jm_32)
      call MDS_BYTESWAPR4( 1, lat0_32)
      call MDS_BYTESWAPR4( 1, lon0_32)
      call MDS_BYTESWAPI4( 1, ndatebc_32)
      call MDS_BYTESWAPR4( 1, undef_32)
#endif

      read(iunit,rec=nrec,err=500) dummy112,
     .     (nymdbc_32(i), nhmsbc_32(i), i=1,ndatebc_32)

      im = im_32
      jm = jm_32
      lat0 = lat0_32
      lon0 = lon0_32
      undef = undef_32

      ndatebc = ndatebc_32
      do i=1,ndatebc
#ifdef _BYTESWAPIO
      call MDS_BYTESWAPI4( 1, nymdbc_32(i))
      call MDS_BYTESWAPI4( 1, nhmsbc_32(i))
#endif
      nymdbc(i) = nymdbc_32(i)
      nhmsbc(i) = nhmsbc_32(i)
      enddo

      return
  500 continue
      print *, 'Error reading boundary condition from unit ',iunit
      error = .true.
      return
      end
