#include "dims.h"
!
      module solgar_module
!
! Acquire and read netcdf file containing Solomon-Garcia/Zatmos 
!   lower boundary conditions, and interpolate to current model
!   time.
! Sub rdsolgar_bndry acquires and reads the data, defining /dlowbnd/
!   (dlowbnd.h). 
! Sub solgar_bndry interpolates data to current model time,
!   defining /lowbnd/ (lowbnd.h). This routine is called once
!   per model time-step from advnce. Sub rdsolgar_bndry is called
!   once per run, from solgar_bndry the first time that routine
!   is called.
!
! Also:
! Acquire and read netcdf file containing Solomon-Garcia import
!   fields dimensioned (13,zjmx,zkmxp,3), and interpolate to
!   current model time, defining fields in fimport.h.
! Sub rdsolgar_import acquires the file and reads the import data.
! Sub solgar_import interpolates to current model time, defining
!   /fimport/. Solgar_import is called from advnce at each model 
!   timestep, and calls rdsolgar_import in first iteration.
!
      use input_module,only: tempdir
      use netcdf_module,only: handle_ncerr,nc_open,nc_close,
     |  nc_get_var_real
      implicit none
#include "netcdf.inc"
#include "params.h"
!
      character(len=40) :: 
     |  solgar_bndry_mss  = '/TGCM/data/solgar_bndry.nc              '
#if (NLEV==44)
! Import file at zp -17 to +5 by 0.5:
      character(len=40) :: 
     |  solgar_import_mss = '/TGCM/data/solgar_import.nc             '
#elif (NLEV==88)
! Import file at zp -17 to +5 by 0.25:
      character(len=40) :: 
     |  solgar_import_mss = '/TGCM/data/solgar_import_dblzp5.nc      '
#elif (NLEV==96)
! Import file at zp -17 to +7 by 0.25:
      character(len=40) :: 
     |  solgar_import_mss = '/TGCM/data/solgar_import_dblz.nc        '
#endif
      contains
!-------------------------------------------------------------------
      subroutine rdsolgar_bndry
!
! Acquire and read netcdf file with Solomon-Garcia/zatmos lower
! boundaries, defining /dlowbnd/ (dlowbnd.h).
!
#include "dlowbnd.h"
      character(len=120) :: char120
      character(len=80) :: varname,dskfile
      integer,parameter :: mxdims=3
      integer :: ncid,istat,nlat,ntnlev,i,itype,natts,len,ndims,
     |  iddims(mxdims),idunlim,ngatts,nvars
      integer :: id_lat,id_tnlevels
      real :: glat(zjmx),tnlevels(ntndown),fmin,fmax
!
! Acquire and open mss file:
      call mkdiskflnm(solgar_bndry_mss,dskfile)
      call getms(solgar_bndry_mss,dskfile,tempdir,' ')
      call nc_open(ncid,dskfile,'OLD','READ')
!
! Get number of dims, vars, atts, and id of unlimited dimension:
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!
! Get and verify latitude dimension:
      istat = nf_inq_dimid(ncid,'lat',id_lat)
      istat = nf_inq_dimlen(ncid,id_lat,nlat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting latitude dimension')
      if (nlat /= zjmx) then
        write(6,"(/,'>>> rdsolgar_bndry: bad nlat=',i3,
     |    ' -- should be zjmx=',i3)") nlat,zjmx
        stop 'rdsolgar_bndry'
      endif
!
! Get and verify ntnlev dimension:
      istat = nf_inq_dimid(ncid,'tnlevels',id_tnlevels)
      istat = nf_inq_dimlen(ncid,id_tnlevels,ntnlev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting tnlevels dimension')
      if (ntnlev /= ntndown) then
        write(6,"(/,'>>> rdsolgar_bndry: bad ntnlev=',i3,
     |    ' -- should be ntndown=',i3)") ntnlev,ntndown
        stop 'rdsolgar_bndry'
      endif
!
! Read data into /dlowbnd/ (dlowbnd.h):
!     real :: dh2olb,dh2lb,dcolb,dco2lb,doxlb,dnozlb,dch4lb,
!    |  dhoxlb,dtndown
!
      do i=1,nvars
        istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
        select case(trim(varname))
          case('lat')      ! latitude dimension variable
            istat = nc_get_var_real(ncid,i,glat)
!           write(6,"('glat=',/,(6e12.4))") glat
          case('tnlevels') ! tn levels dimension variable
            istat = nc_get_var_real(ncid,i,tnlevels)
!           write(6,"('tnlevels=',/,(6e12.4))") tnlevels
          case('H2O')
            istat = nc_get_var_real(ncid,i,dh2olb)
            call fminmax(dh2olb,zjmx*12,fmin,fmax)
            write(6,"('H2O min,max=',2e12.4)") fmin,fmax
          case('H2')
            istat = nc_get_var_real(ncid,i,dh2lb)
            call fminmax(dh2lb,zjmx*12,fmin,fmax)
            write(6,"('H2  min,max=',2e12.4)") fmin,fmax
          case('CO')
            istat = nc_get_var_real(ncid,i,dcolb)
            call fminmax(dcolb,zjmx*12,fmin,fmax)
            write(6,"('CO  min,max=',2e12.4)") fmin,fmax
          case('CO2')
            istat = nc_get_var_real(ncid,i,dco2lb)
            call fminmax(dco2lb,zjmx*12,fmin,fmax)
            write(6,"('CO2 min,max=',2e12.4)") fmin,fmax
          case('OX')
            istat = nc_get_var_real(ncid,i,doxlb)
            call fminmax(doxlb,zjmx*12,fmin,fmax)
            write(6,"('OX  min,max=',2e12.4)") fmin,fmax
          case('NOZ')
            istat = nc_get_var_real(ncid,i,dnozlb)
            call fminmax(dnozlb,zjmx*12,fmin,fmax)
            write(6,"('NOZ min,max=',2e12.4)") fmin,fmax
          case('CH4')
            istat = nc_get_var_real(ncid,i,dch4lb)
            call fminmax(dch4lb,zjmx*12,fmin,fmax)
            write(6,"('CH4 min,max=',2e12.4)") fmin,fmax
          case('HOX')
            istat = nc_get_var_real(ncid,i,dhoxlb)
            call fminmax(dhoxlb,zjmx*12,fmin,fmax)
            write(6,"('NOX min,max=',2e12.4)") fmin,fmax
          case('TN')
            istat = nc_get_var_real(ncid,i,dtndown)
            call fminmax(dtndown,zjmx*ntndown*12,fmin,fmax)
            write(6,"('TN  min,max=',2e12.4)") fmin,fmax
          case default
            write(6,"('>>> rdsolgar_bndry: unknown variable ',a)") 
     |        trim(varname)
        end select
      enddo 
!
! Close the file:
      call nc_close(ncid)
      write(6,"('Completed read from solgar data file ',a)") 
     |  trim(dskfile)
      write(6,"(72('-'),/)")
      end subroutine rdsolgar_bndry
!-------------------------------------------------------------------
      subroutine solgar_bndry(inyear,inday,isec)
      use cons_module,only: jmax
!
! Get and read netcdf file containing monthly solomon-garcia/zatmos 
!   lower boundaries (sub rdsolgar_bndry, defining /dlowbnd/ in dlowbnd.h), 
!   and interpolate to current model time, defining /lowbnd/ in lowbnd.h.
! This routine is called once per model time step from advance. The file 
!   is acquired and read only in the first call.
!
! Args:
      integer,intent(in) :: inyear,inday,isec
!
      integer,parameter :: nf=8 ! number of Sol-Gar fields, excluding tndown
#include "dlowbnd.h"
      real fin(zjmx,12,nf)
      equivalence(fin,dh2olb)
#include "lowbnd.h"
      real fout(zjmx,nf)
      equivalence(fout,xh2olb)
!
! Local:
      integer,parameter :: npts=13
      integer :: iop(2) = (/3,3/) 
      integer :: itab(3) = (/1,0,0/)
      integer,save :: niter=0
      real,save :: secmidmo(npts),fmidmo(npts,zjmx,nf),
     |  tnmidmo(npts,zjmx,ntndown),fcoeff(npts,zjmx,nf),
     |  tncoeff(npts,zjmx,ntndown)
      real :: tab(3),w(npts),wk(3*npts+1)
      integer :: ip,j,i,k,iyd,istat,iyr,imo,ida
      real :: sec,utsec,cursec
      integer,save :: 
     |  ndmon(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
c                     J  F  M  A  M  J  J  A  S  O  N  D
c
c Get data file and read data for all 12 months into /dlowbnd/:
c (this is done in first call only)
c
      niter = niter+1
      if (niter.eq.1) then ! 1st call
!
! Get and read netcdf file, defining /dlowbnd/:
        call rdsolgar_bndry
c
c Define independent variable for cubspl as midpoints of each month
c   in seconds from beginning of year. Since boundaries are periodic,
c   a 13th midpoint is defined, equivalent to the first:
c
        if (mod(inyear,4).eq.0) ndmon(2) = 29
        secmidmo(1) = float(ndmon(1))*86400./2.
        do i=2,12
          secmidmo(i) = secmidmo(i-1) + float(ndmon(i-1))*86400./2. +
     +      float(ndmon(i))*86400./2.
        enddo
        secmidmo(npts) = secmidmo(12) + float(ndmon(12))*86400./2. +
     +    float(ndmon(1))*86400./2.
c
c Get coefficients of fields at month midpoints:
c (these are saved for input to terp1)
c
        do ip=1,nf
          do j=1,jmax
            do i=1,12
              fmidmo(i,j,ip) = fin(j,i,ip)
            enddo
            fmidmo(npts,j,ip) = fmidmo(1,j,ip)
            call coeff1(npts,secmidmo,fmidmo(1,j,ip),fcoeff(1,j,ip),
     +        iop,1,wk)
          enddo
        enddo
        do k=1,ntndown
          do j=1,jmax
            do i=1,12
              tnmidmo(i,j,k) = dtndown(j,k,i)
            enddo
            tnmidmo(npts,j,k) = tnmidmo(1,j,k)
            call coeff1(npts,secmidmo,tnmidmo(1,j,k),tncoeff(1,j,k),
     +        iop,1,wk)
          enddo
        enddo
c
c Report to stdout:
! (note iyd2date takes 7-digit iyd (yyyyddd), and returns 4-digit year)
c
        iyd = inyear*1000+inday
        call iyd2date(iyd,imo,ida,iyr)  ! imo,ida,iyr = current date
        iyd = (inyear-inyear/100*100)*1000+inday
        iyr = iyr-(iyr/100*100) ! back to 2-digit year
        sec = float(isec)
        if (sec.eq.0.) then
          utsec = 0.
        else
          utsec = sec/(24.*3600.)+amod(sec,sec/(24.*3600.))
        endif
        write(6,"(/72('-'))")
        write(6,"('SOLGAR_BNDRY:')")
        write(6,"('Starting at yyddd = ',i5,' (mon/day/yr = ',i2,
     |    '/',i2,'/',i2,'), ut=',f9.5)") iyd,imo,ida,iyr,utsec
        write(6,"('Lower boundaries for the following species are ',
     |    'defined from the ',/,'  Garcia-Solomon model at 30 km:')")
        write(6,"('  OX, HOX, NOX, CH4, H2O, H2, CO, CO2')")
        write(6,"('Lower boundary for TN (TNDOWN) is from ZATMOS ',
     |    ' (Forbes/msis90)',/,'  at ZP -21.5 to -17.25 by .25')")
        write(6,"('Monthly data read from mss ',a)") 
     |    trim(solgar_bndry_mss)
        write(6,"('These data will be interpolated to the current ',
     |    'time at each iteration,',/,
     |    '  using cubic spline interpolation')")
        write(6,"(72('-')/)")
      endif     ! first call only
c
c cursec = secs at current iter from beginining of year:
c
      cursec = float(inday)*86400.-86400.+isec
c
c Do cubic spline interp of dlobnd to current time, defining lobnd:
c
      do ip=1,nf
        do j=1,zjmx
          call terp1(npts,secmidmo,fmidmo(1,j,ip),fcoeff(1,j,ip),
     +      cursec,1,tab,itab)
          fout(j,ip) = tab(1)
        enddo
      enddo
c
c Interpolate dtndown, defining tndown:
c
      do k=1,ntndown
        do j=1,jmax
          call terp1(npts,secmidmo,tnmidmo(1,j,k),tncoeff(1,j,k),
     +      cursec,1,tab,itab)
          tndown(j,k) = tab(1)
        enddo
      enddo
      return
      end subroutine solgar_bndry
!-------------------------------------------------------------------
      subroutine rdsolgar_import(fimp)
!
! Read netcdf file containing solomon-garcia import fields.
!
#include "netcdf.inc"
      integer,parameter :: npts=13, nf=3
      real,intent(out) :: fimp(npts,zjmx,zkmxp,nf)
      real :: fimp1(npts,zjmx,zkmxp) ! 1 field
      character(len=120) :: char120
      character(len=80) :: varname,dskfile
      integer,parameter :: mxdims=3
      integer :: ncid,istat,nlat,nlev,i,itype,natts,len,ndims,
     |  iddims(mxdims),idunlim,ngatts,nvars
      integer :: id_lat,id_lev
      real :: glat(zjmx),plev(zkmxp),fmin,fmax
!
! Get and open netcdf file:
      call mkdiskflnm(solgar_import_mss,dskfile)
      call getms(solgar_import_mss,dskfile,tempdir,' ')
      istat = nf_open(dskfile,NF_NOWRITE,ncid)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_open for netcdf ',
     |    'file ',a)") trim(dskfile)
        call handle_ncerr(istat,char120)
        ncid = 0
      else
        write(6,"('Opened existing netcdf file ',a,' ncid=',i8)")
     |    trim(dskfile),ncid
      endif
!
! Get number of dims, vars, atts, and id of unlimited dimension:
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!
! Get and verify latitude dimension:
      istat = nf_inq_dimid(ncid,'lat',id_lat)
      istat = nf_inq_dimlen(ncid,id_lat,nlat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting latitude dimension')
      if (nlat /= zjmx) then
        write(6,"(/,'>>> rdsolgar_import: bad nlat=',i3,
     |    ' -- should be zjmx=',i3)") nlat,zjmx
        stop 'rdsolgar_import'
      endif
!
! Get and verify levels dimension:
      istat = nf_inq_dimid(ncid,'lev',id_lev)
      istat = nf_inq_dimlen(ncid,id_lev,nlev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting levels dimension')
      if (nlev /= zkmxp) then
        write(6,"(/,'>>> rdsolgar_import: bad nlev=',i3,
     |    ' -- should be zkmxp=',i3)") nlev,zkmxp
        stop 'rdsolgar'
      endif
!
! Read variables (coord vars lat,lev, data vars N2O, CL, CLO:
      do i=1,nvars
        istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
        select case(trim(varname))
          case('lat')      ! latitude dimension variable
            istat = nc_get_var_real(ncid,i,glat)
!           write(6,"('glat=',/,(6e12.4))") glat
          case('lev')      ! levels dimension variable
            istat = nc_get_var_real(ncid,i,plev)
!           write(6,"('plev=',/,(6e12.4))") plev
          case('N2O')
            istat = nc_get_var_real(ncid,i,fimp1)
            call fminmax(fimp1,zjmx*zkmxp*npts,fmin,fmax)
            write(6,"('N2O min,max=',2e12.4)") fmin,fmax
            fimp(:,:,:,1) = fimp1(:,:,:)
          case('CL')
            istat = nc_get_var_real(ncid,i,fimp1)
            call fminmax(fimp1,zjmx*zkmxp*npts,fmin,fmax)
            write(6,"('CL min,max=',2e12.4)") fmin,fmax
            fimp(:,:,:,2) = fimp1(:,:,:)
          case('CLO')
            istat = nc_get_var_real(ncid,i,fimp1)
            call fminmax(fimp1,zjmx*zkmxp*npts,fmin,fmax)
            write(6,"('CLO min,max=',2e12.4)") fmin,fmax
            fimp(:,:,:,3) = fimp1(:,:,:)
          case default
            write(6,"('rdsolgar_import: unknown variable ',a)") 
     |        trim(varname)
        end select
      enddo
      end subroutine rdsolgar_import
!-------------------------------------------------------------------
      subroutine solgar_import(inyear,inday,isec)
#include "fimport.h"
!
! Args:
      integer,intent(in) :: inyear,inday,isec
!
! Local:
      integer,parameter :: npts=13,nf=3
      integer :: niter=0
      integer :: i,j,k,ip
      integer,save :: 
     |  ndmon(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
!                     J  F  M  A  M  J  J  A  S  O  N  D
      real,save :: 
     |  secmidmo(npts),              ! midpoint of each month in seconds
     |  fcoeff(npts,zjmx,zkmxp,nf),  ! coefficients for interpolation
     |  fimp(npts,zjmx,zkmxp,nf)     ! imported data, read from mss
      integer :: iop(2) = (/3,3/)
      integer :: itab(3) = (/1,0,0/)
      real :: tab(3),w(npts),wk(3*npts+1)
      integer :: iyd,imo,ida,iyr
      real :: sec,utsec,cursec
      character(len=8) :: flab8(nf) = 
     |  (/'N2O     ','CL      ','CLO     '/)
      real fout(zjmx,zkmxp,nf)       ! equiv to /fimport/ interp fields
      equivalence(fout,rwn2o)
!
      niter = niter+1
      if (niter==1) then
        call rdsolgar_import(fimp)
!
! Define independent variable for cubspl as midpoints of each month
!   in seconds from beginning of year. Since boundaries are periodic,
!   a 13th midpoint is defined, equivalent to the first:
!
        if (mod(inyear,4).eq.0) ndmon(2) = 29
        secmidmo(1) = float(ndmon(1))*86400./2.
        do i=2,12
          secmidmo(i) = secmidmo(i-1) + float(ndmon(i-1))*86400./2. +
     +      float(ndmon(i))*86400./2.
        enddo
        secmidmo(npts) = secmidmo(12) + float(ndmon(12))*86400./2. +
     +    float(ndmon(1))*86400./2.
!
! Get coefficients of fields at month midpoints:
! (these are saved for input to terp1)
!
        do ip=1,nf
          do k=1,zkmxp
            do j=1,zjmx
              call coeff1(npts,secmidmo,fimp(1,j,k,ip),fcoeff(1,j,k,ip),
     +          iop,1,wk)
            enddo
          enddo
        enddo
!
! Report to stdout (1st iter only):
! (note iyd2date takes 7-digit iyd (yyyyddd), and returns 4-digit year)
!
        iyd = inyear*1000+inday
        call iyd2date(iyd,imo,ida,iyr)  ! imo,ida,iyr = current date
        iyd = (inyear-inyear/100*100)*1000+inday
        iyr = iyr-(iyr/100*100) ! back to 2-digit year
        sec = float(isec)
        if (sec.eq.0.) then
          utsec = 0.
        else
          utsec = sec/(24.*3600.)+amod(sec,sec/(24.*3600.))
        endif
        write(6,"(/72('-'))")
        write(6,"('SOLGAR_IMPORT:')")
        write(6,"('Starting at yyddd = ',i5,' (mon/day/yr = ',i2,
     |    '/',i2,'/',i2,'), ut=',f9.5)") iyd,imo,ida,iyr,utsec
        write(6,"('Fields (13,ZJMX,ZKMXP,NF) imported from ',
     |    'the Solomon-Garcia 2d model (NF=',i2,'):')") nf
        do ip=1,nf
          write(6,"('  ',a,$)") flab8(ip)
        enddo
        write(6,"(/'Monthly data read from mss ',a)") 
     |    trim(solgar_import_mss)
        write(6,"('These data will be interpolated to the current ',
     |    'time at each iteration,',/,
     |    '  using cubic spline interpolation')")
        write(6,"(72('-')/)")
      endif     ! niter == 1
!
! This code (to return) is executed at beginning of each iter (time step):
! cursec = secs at current iter from beginining of year:
!
      cursec = float(inday)*86400.-86400.+isec
!
! Do cubic spline interp of fimp to current time, defining fout
! (equivalenced to fields in fimport.h)
!
      do ip=1,nf
        do k=1,zkmxp
          do j=1,zjmx
            call terp1(npts,secmidmo,fimp(1,j,k,ip),fcoeff(1,j,k,ip),
     +        cursec,1,tab,itab)
            fout(j,k,ip) = tab(1)
          enddo
        enddo
      enddo
      end subroutine solgar_import
      end module solgar_module