c
      program tst
c
      iyear = 1990
      isec = 0
      lastcall = 0
      niter = 365
      do iday=1,niter
        call rpt(iyear,iday,isec)
        if (iday.eq.niter) lastcall = 1
        call mksolgar(iyear,iday,isec,lastcall)
      enddo
      stop
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine mksolgar(iyear,iday,isec,lastcall)
c
      integer zjmx
      parameter(zjmx=36,ntndown=18,jmax=zjmx)
      parameter(nf=8)	! number of fields excluding dtndown
      parameter(npts=13)
c
c /dlobnd/ is in *CDK DLOWBND (*CA DLOWBND in mods):
c
      common/dlobnd/ dh2olb(zjmx,12),dh2lb(zjmx,12),dcolb(zjmx,12),
     +  dco2lb(zjmx,12),doxlb(zjmx,12),dnozlb(zjmx,12),dch4lb(zjmx,12),
     +  dhoxlb(zjmx,12),dtndown(zjmx,ntndown,12)
      real fin(zjmx,12,nf)
      equivalence(fin,dh2olb)
c
c /lobnd/ is in *CDK LOWBND (*CA LOWBND in mods):
c
      common/lobnd/ xh2olb(zjmx),xh2lb(zjmx),xcolb(zjmx),xco2lb(zjmx),
     +  xoxlb(zjmx),xnozlb(zjmx),xch4lb(zjmx),xhoxlb(zjmx),
     +  tndown(zjmx,ntndown)
      real fout(zjmx,nf)
      equivalence(fout,xh2olb)
c
c For test plots only (do not include in model version):
c
      parameter (mxiter=400)
      real pltout(mxiter,zjmx,nf),plt(12,zjmx),xiter(mxiter),
     +  xmon(12),ylat(zjmx),vp(4)
      integer iasf(13)
      character*72 toplab
      character*8 flab(nf)
      data iasf/13*1/
      data flab
     +  /'H2O     ','H2      ','CO      ','CO2     ','OX      ',
     +   'NOZ     ','CH4     ','HOX     '/
      data xmon/1,2,3,4,5,6,7,8,9,10,11,12/
      data vp /.13,.87,.26,.91/
      data niter/0/
      save pltout,ylat,niter,xmid
c
c Below dimensions must be in model version:
c
      integer ndmon(12),iop(2),itab(3)
      real secmidmo(npts),tab(3),w(npts),wk(3*npts+1),
     +  fmidmo(npts,zjmx,nf),tnmidmo(npts,zjmx,ntndown),
     +  fcoeff(npts,zjmx,nf),tncoeff(npts,zjmx,ntndown)
      character*16 diskfile
      character*40 mssfile
      character*80 errmsg
      data mssfile/'/FOSTER/lowbound/sol_gar_zat.dat        '/
      data diskfile/'sol_gar.dat     '/, lu/95/
      data iop/3,3/, itab/1,0,0/
c                 J  F  M  A  M  J  J  A  S  O  N  D
      data ndmon/31,28,31,30,31,30,31,31,30,31,30,31/
      data ifile/0/
      save ifile,secmidmo,fcoeff,tncoeff,fmidmo,tnmidmo
c
c Set up for ncarg (test version only):
c
      if (ifile.eq.0) then
        call opngks
        call gsclip(0)
        call gsasf(iasf)
        call gsfais(1)
        xmid = 0.5*(vp(1)+vp(2))
        do j=1,jmax
          ylat(j) = -87.5+(j-1)*5.
        enddo
        call cpsetr('ILX',xmid)
        call cpsetr('ILY',-.16)
        call cpseti('ILP',0)
        call cpsetr('ILS',.016)
        call cpseti('SET',0)
        call cpsetr('YC1',ylat(1))
        call cpsetr('YCN',ylat(jmax))
      endif
c
c Get data file and read data for all 12 months into /dlowbnd/:
c (this is done in first call only)
c   
      if (ifile.eq.0) then
        lmss = lenstr(mssfile)
        ldisk = lenstr(diskfile)
        call msread(ier,diskfile(1:ldisk),mssfile(1:lmss),' ',' ')
        if (ier.ne.3.and.ier.ne.0) then
          write(6,"(' ')")
          write(6,"('>>> mksolgar: error from msread: mssfile=',a)")
     +      mssfile(1:lmss)
          call mserror(errmsg)
          write(6,"('    mss error: ',a)") errmsg
          write(6,"(' ')")
          stop 'mksolgar'
        endif
        write(6,"('Acquired file ',a,' to disk file ',a)")
     +    mssfile(1:lmss),diskfile(1:ldisk)
        open(lu,file=diskfile(1:ldisk),status='OLD',
     +    form='UNFORMATTED',err=900,iostat=istat)
        read(lu) dh2olb,dh2lb,dcolb,dco2lb,doxlb,dnozlb,dch4lb,
     +    dhoxlb,dtndown
        ifile = 1
        write(6,"('Read /dlowbnd/ from disk:')")
        do ip=1,nf
          call fminmax(fin(1,1,ip),zjmx*12,rmin,rmax,1.e36)
          write(6,"('ip=',i2,' min,max=',2e12.4)") ip,rmin,rmax
        enddo
        call fminmax(dtndown,zjmx*ntndown,12,rmin,rmax,1.e36)
        write(6,"('dtndown min,max=',2e12.4)") rmin,rmax
c
c Define data points as midpoints of each month in seconds from
c   beginning of year (independent variable for cubspl). Since
c   boundaries are periodic, a 13th midpoint is defined, equivalent
c   to the first:
c
        if (mod(iyear,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       write(6,"('Midpoint secs for each month +1 = ',/(5f12.2))")
c    +    secmidmo
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 Make plots of input data (test version only):
c     subroutine pltline(plt,x,mx,nx,xlab,ylab,toplab,logplt,spval)
c
        call set(vp(1),vp(2),vp(3),vp(4),1.,12.,-90.,90.,1)
        call cpsetr('XC1',1.)
        call cpsetr('XCM',12.)
        do ip=1,nf
          do j=1,jmax
            plt(:,j) = fin(j,:,ip)
          enddo
          call clearstr(toplab)
          lflab = lenstr(flab(ip))
          write(toplab,"('Monthly Sol-Gar for ',a)") flab(ip)(1:lflab)
          call cpcnrc(plt,12,12,jmax,0.,0.,0.,1,-1,-1)
          call labrect(xmon,12,ylat,jmax,'MONTH','LATITUDE',0.)
          call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+.03,0.)
          call frame
          write(6,"('Contoured monthly Sol-Gar data for field ',a)")
     +      flab(ip)(1:lflab)
        enddo        
      endif	! first call only
c   
c Find current secs from beginining of year:
c
      cursec = float(iday)*86400.-86400.+isec
c
c Do cubic spline interp of dlobnd to current time, defining lobnd:
c
c
      do ip=1,nf
        do j=1,jmax
          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   dtndown(zjmx,ntndown,12),tndown(zjmx,ntndown)
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
c
c Plot interpolated values (test version only):
c     real pltout(mxiter,zjmx,nf),plt(12,zjmx),xiter(mxiter),
c    +  xmon(12),ylat(zjmx),vp(4)
c
      niter = niter+1
      if (niter.gt.mxiter) write(6,"('>>> WARNING: niter > mxiter:',
     +  ' niter=',i3,' mxiter=',i3)") niter,mxiter
      xiter(niter) = float(niter)
      do ip=1,nf
        pltout(niter,:,ip) = fout(:,ip)
      enddo
      if (lastcall.gt.0) then
        if (niter.ge.3) then
          call set(vp(1),vp(2),vp(3),vp(4),1.,float(niter),-90.,90.,1)
          call cpsetr('XC1',1.)
          call cpsetr('XCM',float(niter))
          do ip=1,nf
            call clearstr(toplab)
            lflab = lenstr(flab(ip))
            write(toplab,"('Interpolated Sol-Gar for ',a,' niter=',i4)") 
     +        flab(ip)(1:lflab),niter
            call cpcnrc(pltout(1,1,ip),mxiter,niter,jmax,0.,0.,0.,1,
     +        -1,-1)
            call labrect(xiter,niter,ylat,jmax,'ITER','LATITUDE',0.)
            call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+.03,0.)
            call frame
            write(6,"('Contoured interpolated Sol-Gar data for field ',
     +        a,' niter=',i4)") flab(ip),niter
          enddo
        endif
        call clsgks
      endif
      return
c
c Error opening disk file:
c
 900  continue
      write(6,"('>>> mksolgar: Error opening disk file ',a)")
     +  diskfile(1:ldisk)
      stop 'mksolgar'
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      function lenstr(str)
      character*(*) str
c
c Return index to last non-blank char in str
c
      length = len(str)
      do i=length,1,-1
        if (str(i:i).ne.' ') then
          lenstr = i
          return
        endif
      enddo
      lenstr = 0
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine fminmax(f,n,rmin,rmax,spval)
      dimension f(n)
c
      rmin = 1.e36
      rmax = -1.e36
      do i=1,n
        if (f(i).ne.spval.and.f(i).gt.rmax) rmax = f(i)         
        if (f(i).ne.spval.and.f(i).lt.rmin) rmin = f(i)         
      enddo
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine rpt(iyear,iday,isec)
c
c Temp version only (not in model):
c
      iyd = (iyear-iyear/100*100)*1000+iday
      call iyd2date(iyd,imo,ida,iyr)	! imo,ida,iyr = current date
      sec = float(isec)
      if (sec.eq.0.) then 
        utsec = 0.
      else
        utsec = sec/(24.*3600.)+amod(sec,sec/(24.*3600.))
      endif
      write(6,"('yyddd=',i5,' month/day/yr = ',i2,'/',i2,'/',
     +  i2,' isec=',i5,' (ut=',f8.4,')')") iyd,imo,ida,iyr,isec,utsec
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine clearstr(str)
c 
c Set given string to all blanks 
c 
      character*(*) str
      length = len(str)
      do i=1,length
        str(i:i) = ' '
      enddo
      return
      end