!
      module zatmos_module
!
! This module reads a netcdf file containing t,u,z daily data, varying 
!   in latitude, from the zatmos emperical model. The file is read only
!   the first time zatmos_bndry is called (from sub lowbound, bndry.F). 
!   Then, when zatmos_bndry is called at every timestep from advance,
!   the data are interpolated to the current model time, using a cubic
!   spline (see cubspl.F). The mss file is /TGCM/data/zatmos_bndry.F.
!   See also hao:/home/tgcm/zatmos.
! If the model is not being advanced in calendar time (calendar_advance==0),
!   then uninterpolated zatmos data is used corresponding to the model
!   starting day (start_day). In this case, iday and isec are ignored.
! Note this module is used only if ncep == 0 (i.e., not an ncep run)
!
! Currently (2/05), only zatmos_tn is used in the model (see dt.F).
!   Un and ht are still read from the file by zatmos_bndry, but the
!   interpolation calls for un and ht are commented out.
!
! 2/05 btf: Module rewritten to include interpolation for timegcm1.2.
!      (Formerly, these arrays were updated only at day boundaries).
!
      use params_module,only: nlat
      implicit none
!
! Latitude varying ht, tn, and un, interpolated to current model time
! at every timestep (main output of this module):
      real,dimension(nlat) :: zatmos_ht,zatmos_tn,zatmos_un
!
! Data read from netcdf file:
      integer,parameter,private :: ndays=365
      real,dimension(ndays,nlat),private :: tn,ht,un
!
! Coefficients for cubic spline, and other needed arrays:
      real,dimension(ndays+1,nlat),private :: 
     |  tn_coef,ht_coef,un_coef,
     |  tn_ndp1,ht_ndp1,un_ndp1
      real,private :: secmidday(ndays+1)
!
      contains
!-----------------------------------------------------------------------
      subroutine zatmos_bndry(istep,iday,isec)
!
      use input_module,only: zatmos_ncfile,calendar_advance,start_day
      use nchist_module,only: handle_ncerr,nc_open,nc_close
      use params_module,only: nlevp1 ! for diags only
      use mpi_module,only: lon0,lon1 ! for diags only
      use addfld_module,only: addfld
      implicit none
!
#ifdef SUN
#include <netcdf3.inc>
#else
#include <netcdf.inc>
#endif
!
! Args:
      integer,intent(in) :: iday,isec,istep
!
! Local:
      integer,parameter :: ndays=365, mxdims=10
      integer,parameter :: nsecsperday = 24*60*60    ! 86400 secs/day
      integer,parameter :: nsecsnoon = nsecsperday/2 ! seconds at midday
      character(len=80) :: dskfile,varname
      integer,save :: ncalls=0
      integer :: istat,nlatrd,ncid,id_lat
      real :: fmin,fmax,cursec
      integer :: i,j,itype,natts,ndims,iddims(mxdims),idunlim,
     |  ngatts,nvars
      integer :: iop(2) = (/3,3/) 
      integer :: itab(3) = (/1,0,0/)
      integer(kind=8) :: isecyr
      real :: tab(3),wk(3*ndays+2)
      real,dimension(nlevp1,lon0:lon1) :: 
     |  zat_tn, zat_un, zat_ht ! diags only
!
! Exec:
      ncalls = ncalls+1
!
! Check args:
      if (iday < 0 .or. iday > 365) then
        write(6,"(/,'>>> zatmos_bndry: bad iday=',i5)") iday
        call shutdown('zatmos_bndry')
      endif
      if (isec < 0 .or. isec > nsecsperday) then
        write(6,"(/,'>>> zatmos_bndry: bad isecs=',i8)") isec
        call shutdown('zatmos_bndry')
      endif
!
! Acquire, open, and read data file (1st call only):
      if (ncalls==1) then
        dskfile = ' '
        call getfile(zatmos_ncfile,dskfile)
        write(6,"(/,72('-'))")
        write(6,"('ZATMOS_BNDRY initial read (iday=',i3,
     |    ' isec=',i6,')')") iday,isec
        call nc_open(ncid,dskfile,'OLD','READ')
        write(6,"('Opened netcdf file ',a)") trim(dskfile)
!
! Get and verify latitude dimension:
        istat = nf_inq_dimid(ncid,'lat',id_lat)
        istat = nf_inq_dimlen(ncid,id_lat,nlatrd)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'Error getting latitude dimension')
        if (nlatrd /= nlat) then
          write(6,"(/,'>>> zatmos_bndry: bad nlat read=',i3,
     |      ' -- should be nlat=',i3)") nlatrd,nlat
          write(6,"('Closing disk file ',a)") trim(dskfile)
          call nc_close(ncid)
          call shutdown('zatmos_bndry')
        endif
!
! Get number of dims, vars, atts, and id of unlimited dimension:
        istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!
! Read the variables (ht,tn,un are module data):
        write(6,"('Data for ',i3,' days, with nlat=',i3,
     |    ' for each day.')") ndays,nlat
        write(6,"('These data will be interpolated to model time',
     |    ' at every timestep.',/,'See source file zatmos.F')")
        do i=1,nvars
          istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
          select case(trim(varname))
            case('lat')      ! latitude coord variable
            case('HT')
              istat = nf_get_var_double(ncid,i,ht)
              call fminmax(ht,ndays*nlat,fmin,fmax)
              write(6,"('HT min,max=',2e12.4)") fmin,fmax
            case('TN')
              istat = nf_get_var_double(ncid,i,tn)
              call fminmax(tn,ndays*nlat,fmin,fmax)
              write(6,"('TN min,max=',2e12.4)") fmin,fmax
            case('UN')
              istat = nf_get_var_double(ncid,i,un)
              call fminmax(un,ndays*nlat,fmin,fmax)
              write(6,"('UN min,max=',2e12.4)") fmin,fmax
            case default
              write(6,"('>>> zatmos_bndry: unknown variable ',a)") 
     |          trim(varname)
          end select
        enddo 
!
! Define independent variable for cubspl as midpoints of each day
!   in seconds from beginning of the year. Since boundaries are
!   periodic, a 366th point is defined, equivalent to day 1.
! Coefficients are returned by sub coeff1 (cubspl.F), and saved
!   as module data.
!
! secmidday = seconds into the year at each midday (noon):
        secmidday(1) = float(nsecsnoon)
        do i=2,ndays+1
          secmidday(i) = secmidday(i-1)+float(nsecsperday)
        enddo
        do j=1,nlat
          do i=1,ndays
            tn_ndp1(i,j) = tn(i,j)
!           un_ndp1(i,j) = un(i,j)
!           ht_ndp1(i,j) = ht(i,j)
          enddo
          tn_ndp1(ndays+1,:) = tn_ndp1(1,:) ! periodic point
!         un_ndp1(ndays+1,:) = un_ndp1(1,:)
!         ht_ndp1(ndays+1,:) = ht_ndp1(1,:)
          call coeff1(ndays+1,secmidday,tn_ndp1(1,j),tn_coef(1,j),
     |      iop,1,wk)
!         call coeff1(ndays+1,secmidday,un_ndp1(1,j),un_coef(1,j),
!    |      iop,1,wk)
!         call coeff1(ndays+1,secmidday,ht_ndp1(1,j),ht_coef(1,j),
!    |      iop,1,wk)
        enddo
!
! Close data file:
        call nc_close(ncid)
        write(6,"('Closed netcdf file ',a)") trim(dskfile)
        write(6,"(72('-')/)")
      endif ! ncalls==1
!
! The rest of this routine is executed at every timestep 
!   (called from advance).
!
! If this is a perpetual run (calendar day not advanced), simply use
! uninterpolated zatmos data for the start_day (inputs iday and secs
! are ignored):
!
      if (calendar_advance==0) then
!       zatmos_tn(:) = tn(start_day,:)+10.
        zatmos_tn(:) = tn(start_day,:)
!       zatmos_tn(:) = tn(start_day,:)-10.
!       zatmos_un(:) = un(start_day,:)
!       zatmos_ht(:) = ht(start_day,:)
        if (ncalls==1) then
          write(6,"(/,72('-'))")
          write(6,"('ZATMOS_BNDRY first call:',/,'Using uninterpolated',
     |      ' zatmos data from start day ',i3,' (calendar_advance=',
     |      i3,')')") start_day,calendar_advance
          write(6,"('zatmos_tn=',/,(6e12.4))") zatmos_tn
!         write(6,"('zatmos_un=',/,(6e12.4))") zatmos_un
!         write(6,"('zatmos_ht=',/,(6e12.4))") zatmos_ht
          write(6,"(72('-')/)")
        endif
! Save to secondary histories (will vary in latitude only, not in time):
!       do j=1,nlat
!         zat_tn(:,:) = zatmos_tn(j)
!         zat_un(:,:) = zatmos_un(j)
!         zat_ht(:,:) = zatmos_ht(j)
!         call addfld('ZAT_TN',' ',' ',zat_tn,'lev',1,nlevp1,
!    |      'lon',lon0,lon1,j)
!         call addfld('ZAT_UN',' ',' ',zat_un,'lev',1,nlevp1,
!    |      'lon',lon0,lon1,j)
!         call addfld('ZAT_HT',' ',' ',zat_ht,'lev',1,nlevp1,
!    |      'lon',lon0,lon1,j)
!       enddo
!       return
      endif ! calendar_advance==0
!
! Interpolate to current model time using coefficients from coeff1
!   calls above (terp1 is in cubspl.F):
!
      isecyr = iday*nsecsperday-nsecsperday+isec ! 8-byte int
      cursec = float(isecyr)
      if (iday==0) cursec = float(isec)
      do j=1,nlat
        call terp1(ndays+1,secmidday,tn_ndp1(1,j),tn_coef(1,j),cursec,
     |    1,tab,itab)
!       zatmos_tn(j) = tab(1)+10.
        zatmos_tn(j) = tab(1)
!       zatmos_tn(j) = tab(1)-10.
!       call terp1(ndays+1,secmidday,un_ndp1(1,j),un_coef(1,j),cursec,
!    |    1,tab,itab)
!       zatmos_un(j) = tab(1)
!       call terp1(ndays+1,secmidday,ht_ndp1(1,j),ht_coef(1,j),cursec,
!    |    1,tab,itab)
!       zatmos_ht(j) = tab(1)
!
! Save interpolated arrays to secondary histories 
! (these will vary in latitude and time only):
!       zat_tn(:,:) = zatmos_tn(j)
!       zat_un(:,:) = zatmos_un(j)
!       zat_ht(:,:) = zatmos_ht(j)
!       call addfld('ZAT_TN',' ',' ',zat_tn,
!    |    'lev',1,nlevp1,'lon',lon0,lon1,j)
!       call addfld('ZAT_UN',' ',' ',zat_un,
!    |    'lev',1,nlevp1,'lon',lon0,lon1,j)
!       call addfld('ZAT_HT',' ',' ',zat_ht,
!    |    'lev',1,nlevp1,'lon',lon0,lon1,j)
      enddo ! j=1,nlat
!
! Print interpolated arrays to stdout at first time-step:
      if (ncalls==1) then
        write(6,"(/,72('-'))")
        write(6,"('ZATMOS_BNDRY first call: iday=',i3,' isec=',i6,
     |    ' (hours=',f6.2,') Interpolated data:')") iday,isec,
     |    float(isec)/3600.
        write(6,"('zatmos_tn(1:nlat)=',/,(6f10.4))") zatmos_tn
!       write(6,"('zatmos_un(1:nlat)=',/,(6f10.4))") zatmos_un
!       write(6,"('zatmos_ht(1:nlat)=',/,(6f10.4))") zatmos_ht
        write(6,"(72('-')/)")
      endif
      end subroutine zatmos_bndry
      end module zatmos_module