#include <defs.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 dh2olb, etc.
! Sub solgar_bndry interpolates data to current model time,
!   defining xh2olb, etc. This routine is called once per model time-step 
!   from advance. 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,nlat,nlevp1,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
!   rwn2o, etc. Solgar_import is called from advance at each model 
!   timestep, and calls rdsolgar_import in first iteration.
!
! 7/03 btf: adapted for timegcm1 (eliminated equivalences)
!
      use params_module,only: ntndown,nlat,nlevp1
      use input_module,only: tempdir,
     |  solgar_bndry_file,solgar_import_file
      use nchist_module,only: handle_ncerr,nc_open,nc_close
      implicit none
!
#ifdef SUN
#include <netcdf3.inc>
#else
#include <netcdf.inc>
#endif
!
! Lower boundaries imported from Solomon-Garcia:
!
      real,dimension(nlat,12) :: dh2olb,dh2lb,dcolb,dco2lb,doxlb,
     |  dnozlb,dch4lb,dhoxlb
      real,dimension(nlat) :: xh2olb,xh2lb,xcolb,xco2lb,xoxlb,
     |  xnozlb,xch4lb,xhoxlb
      real :: dtndown(nlat,ntndown,12),tndown(nlat,ntndown)
!
! Fields imported by solgar_import and interpolated to current model time:
!
      real,dimension(nlat,nlevp1) :: rwn2o,rwcl,rwclo
!
! Multipliers for each species (added for Feb, 2005 co2 doubling study):
      real,parameter ::
     |  mult_h2olb  = 1.,  ! for xh2olb
     |  mult_h2lb   = 1.,  ! for xh2lb
     |  mult_colb   = 1.,  ! for xcolb
     |  mult_co2lb  = 1.,  ! for xco2lb (xco2lb currently not used in comp_co2)
     |  mult_oxlb   = 1.,  ! for xoxlb
     |  mult_nozlb  = 1.,  ! for xnozlb
     |  mult_ch4lb  = 1.,  ! for xch4lb
     |  mult_hoxlb  = 1.,  ! for xhoxlb
     |  mult_n2o    = 1.,  ! for rwn2o
     |  mult_cl     = 1.,  ! for rwcl
     |  mult_clo    = 1.   ! for rwclo
!
!    |  mult_h2olb  = 2.,  ! for xh2olb
!    |  mult_h2lb   = 2.,  ! for xh2lb
!    |  mult_colb   = 2.,  ! for xcolb
!    |  mult_co2lb  = 2.,  ! for xco2lb (xco2lb currently not used in comp_co2)
!    |  mult_oxlb   = .5,  ! for xoxlb
!    |  mult_nozlb  = 3.,  ! for xnozlb
!    |  mult_ch4lb  = 2.,  ! for xch4lb
!    |  mult_hoxlb  = 2.,  ! for xhoxlb
!    |  mult_n2o    = 1.5,  ! for rwn2o
!    |  mult_cl     = 2.,  ! for rwcl
!    |  mult_clo    = 2.   ! for rwclo
      contains
!-------------------------------------------------------------------
      subroutine rdsolgar_bndry
!
! Acquire and read netcdf file with Solomon-Garcia/zatmos lower
! boundaries, defining dlowbnd (dh2olb, dh2lb, etc).
!
      character(len=120) :: char120
      character(len=80) :: varname,dskfile
      integer,parameter :: mxdims=3
      integer :: ncid,istat,nlatrd,ntnlev,i,itype,natts,len,ndims,
     |  iddims(mxdims),idunlim,ngatts,nvars
      integer :: id_lat,id_tnlevels
      real :: glat(nlat),tnlevels(ntndown),fmin,fmax
!
! Acquire and open mss file:
!     call mkdiskflnm(solgar_bndry_file,dskfile)
!     call getms(solgar_bndry_file,dskfile,tempdir)
      dskfile = ' '
      call getfile(solgar_bndry_file,dskfile)
      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,nlatrd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting latitude dimension')
      if (nlatrd /= nlat) then
        write(6,"(/,'>>> rdsolgar_bndry: bad nlat read=',i3,
     |    ' -- should be nlat=',i3)") nlatrd,nlat
        call shutdown('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
        call shutdown('rdsolgar_bndry')
      endif
!
! Read data into dh2olb, et.al:
!
      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 = nf_get_var_double(ncid,i,glat)
!           write(6,"('glat=',/,(6e12.4))") glat
          case('tnlevels') ! tn levels dimension variable
            istat = nf_get_var_double(ncid,i,tnlevels)
!           write(6,"('tnlevels=',/,(6e12.4))") tnlevels
          case('H2O')
            istat = nf_get_var_double(ncid,i,dh2olb)
            call fminmax(dh2olb,nlat*12,fmin,fmax)
            write(6,"('H2O min,max=',2e12.4)") fmin,fmax
          case('H2')
            istat = nf_get_var_double(ncid,i,dh2lb)
            call fminmax(dh2lb,nlat*12,fmin,fmax)
            write(6,"('H2  min,max=',2e12.4)") fmin,fmax
          case('CO')
            istat = nf_get_var_double(ncid,i,dcolb)
            call fminmax(dcolb,nlat*12,fmin,fmax)
            write(6,"('CO  min,max=',2e12.4)") fmin,fmax
          case('CO2')
            istat = nf_get_var_double(ncid,i,dco2lb)
            call fminmax(dco2lb,nlat*12,fmin,fmax)
            write(6,"('CO2 min,max=',2e12.4)") fmin,fmax
          case('OX')
            istat = nf_get_var_double(ncid,i,doxlb)
            call fminmax(doxlb,nlat*12,fmin,fmax)
            write(6,"('OX  min,max=',2e12.4)") fmin,fmax
          case('NOZ')
            istat = nf_get_var_double(ncid,i,dnozlb)
            call fminmax(dnozlb,nlat*12,fmin,fmax)
            write(6,"('NOZ min,max=',2e12.4)") fmin,fmax
          case('CH4')
            istat = nf_get_var_double(ncid,i,dch4lb)
            call fminmax(dch4lb,nlat*12,fmin,fmax)
            write(6,"('CH4 min,max=',2e12.4)") fmin,fmax
          case('HOX')
            istat = nf_get_var_double(ncid,i,dhoxlb)
            call fminmax(dhoxlb,nlat*12,fmin,fmax)
            write(6,"('NOX min,max=',2e12.4)") fmin,fmax
          case('TN')
            istat = nf_get_var_double(ncid,i,dtndown)
            call fminmax(dtndown,nlat*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)
!
! Get and read netcdf file containing monthly solomon-garcia/zatmos 
!   lower boundaries (sub rdsolgar_bndry, defining dh2olb, etc.
!   and interpolate to current model time, defining xh2olb, etc.
! 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
      real fin(nlat,12,nf)
!
! 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,nlat,nf),
     |  tnmidmo(npts,nlat,ntndown),fcoeff(npts,nlat,nf),
     |  tncoeff(npts,nlat,ntndown)
      real :: tab(3),w(npts),wk(3*npts+1)
      integer :: ip,j,i,k,iyd,istat,iyr,imo,ida
      real :: sec,utsec,cursec,fmin,fmax
      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
!
! Get data file and read data for all 12 months into /dlowbnd/:
! (this is done in first call only)
!
      niter = niter+1
      if (niter.eq.1) then ! 1st call
!
! Get and read netcdf file, defining /dlowbnd/:
        call rdsolgar_bndry
!
! 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)
!
        fin(:,:,1) = dh2olb(:,:)
        fin(:,:,2) = dh2lb (:,:)
        fin(:,:,3) = dcolb (:,:)
        fin(:,:,4) = dco2lb(:,:)
        fin(:,:,5) = doxlb (:,:)
        fin(:,:,6) = dnozlb(:,:)
        fin(:,:,7) = dch4lb(:,:)
        fin(:,:,8) = dhoxlb(:,:)
        do ip=1,nf
          do j=1,nlat
            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,nlat
            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
!
! Report to stdout:
! (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_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 ',a)") 
     |    trim(solgar_bndry_file)
        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
!
! cursec = secs at current iter from beginining of year:
!
      cursec = float(inday)*86400.-86400.+isec
!
! Do cubic spline interp of dlobnd to current time, defining lobnd:
!
      do ip=1,nf
        do j=1,nlat
          call terp1(npts,secmidmo,fmidmo(1,j,ip),fcoeff(1,j,ip),
     |      cursec,1,tab,itab)
          if (ip==1) xh2olb(j) = tab(1)
          if (ip==2) xh2lb (j) = tab(1)
          if (ip==3) xcolb (j) = tab(1)
          if (ip==4) xco2lb(j) = tab(1)
          if (ip==5) xoxlb (j) = tab(1)
          if (ip==6) xnozlb(j) = tab(1)
          if (ip==7) xch4lb(j) = tab(1)
          if (ip==8) xhoxlb(j) = tab(1)
        enddo
      enddo
!     write(6,"('solgar_bndry: xoxlb=',/,(6e12.4))") xoxlb
!
! Interpolate dtndown, defining tndown:
!
      do k=1,ntndown
        do j=1,nlat
          call terp1(npts,secmidmo,tnmidmo(1,j,k),tncoeff(1,j,k),
     |      cursec,1,tab,itab)
          tndown(j,k) = tab(1)
        enddo
      enddo
!
! Modify lower boundary by multiplication factors provided 
! in module parameters above (for co2 doubling study).
!
      if (niter == 1) then
        write(6,"('Multipliers for Solomon-Garcia lower boundaries:')")
        write(6,"('  mult_h2olb  = ',f8.2)") mult_h2olb
        write(6,"('  mult_h2lb   = ',f8.2)") mult_h2lb
        write(6,"('  mult_colb   = ',f8.2)") mult_colb
        write(6,"('  mult_co2lb  = ',f8.2)") mult_co2lb
        write(6,"('  mult_oxlb   = ',f8.2)") mult_oxlb
        write(6,"('  mult_nozlb  = ',f8.2)") mult_nozlb
        write(6,"('  mult_ch4lb  = ',f8.2)") mult_ch4lb
        write(6,"('  mult_hoxlb  = ',f8.2)") mult_hoxlb
      endif
!
! Whole-array operations:
      if (mult_h2olb /= 1.) xh2olb = xh2olb * mult_h2olb
      if (mult_h2lb  /= 1.) xh2lb  = xh2lb  * mult_h2lb
      if (mult_colb  /= 1.) xcolb  = xcolb  * mult_colb
      if (mult_co2lb /= 1.) xco2lb = xco2lb * mult_co2lb
      if (mult_oxlb  /= 1.) xoxlb  = xoxlb  * mult_oxlb
      if (mult_nozlb /= 1.) xnozlb = xnozlb * mult_nozlb
      if (mult_ch4lb /= 1.) xch4lb = xch4lb * mult_ch4lb
      if (mult_hoxlb /= 1.) xhoxlb = xhoxlb * mult_hoxlb

!     call fminmax(xh2olb,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xh2olb min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(xh2lb ,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xh2lb  min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(xcolb ,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xcolb  min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(xco2lb,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xco2lb min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(xoxlb ,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xoxlb  min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(xnozlb,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xnozlb min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(xch4lb,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xch4lb min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(xhoxlb,nlat,fmin,fmax)
!     write(6,"('solgar_bndry: inyear=',i4,' inday=',i4,' isec=',
!    |  i6,' xhoxlb min,max=',2e12.4)") inyear,inday,isec,fmin,fmax

      end subroutine solgar_bndry
!-------------------------------------------------------------------
      subroutine rdsolgar_import(fimp)
!
! Read netcdf file containing solomon-garcia import fields.
!
!
      integer,parameter :: npts=13, nf=3
      real,intent(out) :: fimp(npts,nlat,nlevp1,nf)
      real :: fimp1(npts,nlat,nlevp1) ! 1 field
      character(len=120) :: char120
      character(len=80) :: varname,dskfile
      integer,parameter :: mxdims=3
      integer :: ncid,istat,nlevrd,i,itype,natts,len,ndims,
     |  iddims(mxdims),idunlim,ngatts,nvars,nlatrd
      integer :: id_lat,id_lev
      real :: glat(nlat),plev(nlevp1),fmin,fmax
!
! Get and open netcdf file:
!
      dskfile = ' '
      call getfile(solgar_import_file,dskfile)
      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)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting id of latitude dimension') 
      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,"(/,'>>> rdsolgar_import: bad nlatrd=',i3,
     |    ' -- should be nlat=',i3)") nlatrd,nlat
        call shutdown('rdsolgar_import')
      endif
!
! Get and verify levels dimension:
      istat = nf_inq_dimid(ncid,'lev',id_lev)
      istat = nf_inq_dimlen(ncid,id_lev,nlevrd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting levels dimension')
      if (nlevrd /= nlevp1) then
        write(6,"(/,'>>> rdsolgar_import: bad nlev read=',i3,
     |    ' -- should be nlevp1=',i3)") nlevrd,nlevp1
        call shutdown('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 = nf_get_var_double(ncid,i,glat)
!           write(6,"('glat=',/,(6e12.4))") glat
          case('lev')      ! levels dimension variable
            istat = nf_get_var_double(ncid,i,plev)
!           write(6,"('plev=',/,(6e12.4))") plev
          case('N2O')
            istat = nf_get_var_double(ncid,i,fimp1)
            call fminmax(fimp1,nlat*nlevp1*npts,fmin,fmax)
            write(6,"('N2O min,max=',2e12.4)") fmin,fmax
            fimp(:,:,:,1) = fimp1(:,:,:)
          case('CL')
            istat = nf_get_var_double(ncid,i,fimp1)
            call fminmax(fimp1,nlat*nlevp1*npts,fmin,fmax)
            write(6,"('CL min,max=',2e12.4)") fmin,fmax
            fimp(:,:,:,2) = fimp1(:,:,:)
          case('CLO')
            istat = nf_get_var_double(ncid,i,fimp1)
            call fminmax(fimp1,nlat*nlevp1*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)
!
! 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,nlat,nlevp1,nf),  ! coefficients for interpolation
     |  fimp(npts,nlat,nlevp1,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),fmin,fmax
      integer :: iyd,imo,ida,iyr
      real :: sec,utsec,cursec
      character(len=8) :: flab8(nf) = 
     |  (/'N2O     ','CL      ','CLO     '/)
!
      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,nlevp1
            do j=1,nlat
              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 ',a)") 
     |    trim(solgar_import_file)
        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
!
      do ip=1,nf
        do k=1,nlevp1
          do j=1,nlat
            call terp1(npts,secmidmo,fimp(1,j,k,ip),fcoeff(1,j,k,ip),
     |        cursec,1,tab,itab)
            if (ip==1) rwn2o(j,k) = tab(1)
            if (ip==2) rwcl (j,k) = tab(1)
            if (ip==3) rwclo(j,k) = tab(1)
          enddo
        enddo
      enddo
!
! Modify lower boundary by multiplication factors provided 
! in module parameters above (for co2 doubling study).
!
      if (niter == 1) then
        write(6,"('Multipliers for Solomon-Garcia import fields:')")
        write(6,"('  mult_n2o  = ',f8.2)") mult_n2o
        write(6,"('  mult_cl   = ',f8.2)") mult_cl
        write(6,"('  mult_clo  = ',f8.2)") mult_clo
      endif
      do j=1,nlat
        rwn2o(j,1) = rwn2o(j,1) * mult_n2o
        rwcl (j,1) = rwcl (j,1) * mult_cl
        rwclo(j,1) = rwclo(j,1) * mult_clo
      enddo

!     call fminmax(rwn2o,nlat*nlevp1,fmin,fmax)
!     write(6,"('solgar_import: inyear=',i4,' inday=',i4,' isec=',i6,
!    |  ' rwn2o min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(rwcl,nlat*nlevp1,fmin,fmax)
!     write(6,"('solgar_import: inyear=',i4,' inday=',i4,' isec=',i6,
!    |  ' rwcl  min,max=',2e12.4)") inyear,inday,isec,fmin,fmax
!     call fminmax(rwclo,nlat*nlevp1,fmin,fmax)
!     write(6,"('solgar_import: inyear=',i4,' inday=',i4,' isec=',i6,
!    |  ' rwclo min,max=',2e12.4)") inyear,inday,isec,fmin,fmax

      end subroutine solgar_import
      end module solgar_module