c
      program tracsat
c
c 2/94: Main driver for tracsat code (modified from ~/uars/timessat
c       to use gettgcm.a)
c
      include 'tracsat.h'
      pointer(pgcm0,gcm0(1)),(pgcm1,gcm1(1)),(pgcmut,gcmut(1))
      pointer(pgcmht,gcmht(1)),(pcol0,col0(1)),(pcol1,col1(1))
      pointer(pcol,col(1)),(phtcol,htcol(1)),(pglb0,glb0(1)),
     +       (pglb1,glb1(1))
      dimension iasf(13)
      data iasf/13*1/, iter/0/
c 
c User input:
c
 1000 continue
      iter = iter+1
      ieof = 0
      write(6,"(' ')")
      write(6,"('BEGIN PROC ITERATION ',i3)") iter
      call getinp(ieof)
      if (ieof.gt.0) then
        write(6,"('END PROC AFTER ITER=',i3)") iter-1
        stop 'done'
      endif
c
c Prepare for graphics:
c
      call agsetf('AXIS/LEFT/CONTROL.',1.)
      call agsetf('AXIS/BOTTOM/CONTROL.',1.)
      call agsetf('AXIS/TOP/CONTROL.',1.)
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
      call cpsetup(spval)
c
c Read sat data:
c
      lensatf = lenstr(satfile)
      iuars = index(satfile,'uars')
      ishuttle = index(satfile,'shuttle')
      if (iuars.le.0.and.ishuttle.le.0) then
        write(6,"('>>> sat file name does not contain ''uars''',
     +    ' or ''shuttle''')")
        write(6,"('    --these are only two file types known')")
        stop 'satfile'
      elseif (iuars.gt.0) then
        chsat = 'UARS            '
        write(6,"(' ')")
        write(6,"('tracsat calling rduars: iseclim=',2i12)") iseclim
        call rduars(satfile,lusat,iseclim,mxpts,
     +    satut,satlat,satlon,satalt,npts,iydsat,ier)
        if (ier.ne.0) then
          write(6,"('>>> tracsat: error ',i2,' from rduars')") ier
          stop 'rduars'
        endif 
        write(6,"('tracsat after rduars: npts=',i4,' iydsat=',i5)") 
     +    npts,iydsat
      elseif (ishuttle.gt.0) then
        chsat = 'SHUTTLE         '
        call rdshuttl(satfile,lusat,iseclim,mxpts,
     +    satut,satlat,satlon,satalt,npts,iydsat,ier)
        if (ier.ne.0) then
          write(6,"('>>> tracsat: error ',i2,' from rdshuttl')") ier
          stop 'rduars'
        endif 
        write(6,"('tracsat after rdshuttl: npts=',i4,' iydsat=',i5)") 
     +    npts,iydsat
      endif
      if (npts.le.0) then
        write(6,"('>>> No data?? npts=',i1)") npts
        stop 'nodata'
      endif
c
c Bracket first sat ut with gcm if doing time interpolation:
c
      if (inttime.gt.0) then
        do it=1,ntms-1
          gcmut0 = float(mt(2,it)) + float(mt(3,it))/60.
          gcmut1 = float(mt(2,it+1)) + float(mt(3,it+1))/60.
          if (satut(1).ge.gcmut0.and.satut(1).le.gcmut1) then
            it0 = it
            it1 = it+1
            goto 100
          endif
        enddo      
        write(6,"('>>> tracsat: could not bracket first satut: ',
     +    ' satut(1)=',f10.2)") satut(1)
        stop 'satut'
 100    continue
        write(6,"('Bracketed first satut=',f9.3,' between gcmut0=',f9.3,
     +    ' and gcmut1=',f9.3,/' it0=',i3,' it1=',i3)")
     +    satut(1),gcmut0,gcmut1,it0,it1
c
c If not doing time interpolation, get only first history.
c
      else
        gcmut0 = float(mt(2,1)) + float(mt(3,1))/60.
        gcmut1 = -1.
        it0 = 1
        it1 = -1
      endif
c
      iprhdr = 1
      istart = 1
c
c Get first gcm time in gcm0:
c
      call gettgcm(histvol,nhvols,istart,mt(1,it0),pgcm0,
     +  ifget,nfhist,imx,kmx,jmx,nfget,iden,luhist,iprhdr,ivol,
     +  isdyn,istimes,ier)
      if (ier.ne.0) then
        write(6,"('>>> tracsat: error ',i3,' from gettgcm: it=',i3,
     +    ' mt=',i2,':',i2,':',i2)") ier,it,(mt(i,it0),i=1,3)
        stop 'gettgcm'
      else
        write(6,"('tracsat: mtime=',i2,':',i2,':',i2,' isdyn=',l1,
     +    ' istimes=',l1,' kmx=',i2,' nfplt=',i2)")
     +    (mt(i,it0),i=1,3),isdyn,istimes,kmx,nfget
      endif
      model = 'TIMESGCM'
      if (.not.istimes) then
        if (isdyn) then
          model = 'TIEGCM  '
        else
          model = 'TIGCM   '
        endif
      endif
c
c Allocate pressure array now that we know kmx:
c
      zp1 = zp1tigcm
      if (istimes) zp1 = zp1times
      call hpalloc(pgcmzp,kmx,ier,1)
      call hpalloc(pgcmht,kmx,ier,1)
      call hpalloc(pcol0,kmx,ier,1)
      call hpalloc(pcol1,kmx,ier,1)
      call hpalloc(pcol,kmx,ier,1)
      call hpalloc(phtcol,nhtscale,ier,1)
      call setgrid(gcmlat,glat1,dlat,jmx, gcmlon,
     +             glon1,dlon,imx,gcmzp,zp1,dzp,kmx)
c
c Global array for fields to be calculated from history fields
c (as of 5/94, these are  ui,vi,wi,ver6300,ver15m, and ver5.3m, 
c  with only the last two available:
c
      iglb = 0
      do ip=nfhist+1,nftot
        if (ifplt(ip).gt.0) iglb = 1
      enddo
      if (iglb.gt.0) then
        call hpalloc(pglb0,imx*kmx*jmx,ier,1)
        if (inttime.gt.0) call hpalloc(pglb1,imx*kmx*jmx,ier,1)
      endif
c
c ixcolz = index to gcm heights (in gcm0 and gcm1):
      nf = 0
      ixcolz = -1
      do ip=1,nfhist
        if (ifget(ip).gt.0) then
          nf = nf+1
          if (ip.eq.ixz) ixcolz = nf
        endif
      enddo
      if (ixcolz.lt.0) then
        write(6,"('>>> tracsat: could not find ixcolz <<<')")
        stop 'ixcolz'
      endif
c
c Get second time in gcm1 if doing time interp
c (make sure hist is from same model):
c
      istart = 0
      if (inttime.gt.0) then
        call gettgcm(histvol,nhvols,istart,mt(1,it1),pgcm1,
     +    ifget,nfhist,imx,kmx1,jmx,nfget1,iden,luhist,iprhdr,ivol,
     +    isdyn,istimes,ier)
        if (ier.ne.0) then
          write(6,"('>>> tracsat: error ',i3,' from gettgcm: it=',i3,
     +      ' mt=',i2,':',i2,':',i2)") 
     +      ier,it,(mt(i,it1),i=1,3)
          stop 'gettgcm'
        endif
        if (kmx1.ne.kmx.or.nfget1.ne.nfget) then
          write(6,"('>>> tracsat: different models?? kmx,kmx1=',
     +      2i3,' nfplt1,nfplt=',2i3)") kmx,kmx1,nfget1,nfget
          stop 'gettgcm'
        endif
        write(6,"('tracsat: mtime=',i2,':',i2,':',i2,' isdyn=',l1,
     +    ' istimes=',l1,' kmx=',i2,' nfget=',i2)")
     +    (mt(i,it1),i=1,3),isdyn,istimes,kmx,nfget
      endif
      iprhdr = 0
c
c Allocate space for gcmut(npts,nhtscale,nfplt):
c
      call hpalloc(pgcmut,npts*nhtscale*nfplt,ier,1)
      if (ier.ne.0) then
        write(6,"('hpalloc error = ',i6,' for pgcmut')") ier
        stop 'hpalloc'
      endif
c
c Get and interpolate tiegcm at each data point:
c
      do ipt=1,npts
c
c Get next gcm history if necessary:
c
        if (inttime.gt.0) then
 200      continue
          if (satut(ipt).gt.gcmut1.or.
     +       (ipt.gt.1.and.satut(ipt).lt.satutp)) then
            it0 = it1
            gcmut0 = gcmut1
            it1 = it1+1
            if (it1.gt.ntms) then
              write(6,"('>>> tracsat: out of gcm histories: ipt=',i4,
     +          ' satut=',f9.3,' gcmut0=',f9.3,' ntms=',i3)")
     +          ipt,satut(ipt),gcmut0,ntms
              goto 300
            endif
            gcmut1 = float(mt(2,it1)) + float(mt(3,it1))/60.
            nf = 0
            do ip=1,nftot
              if (ifplt(ip).gt.0) then
                nf = nf+1
                call xfergcm(gcm0,gcm1,imx,kmx,jmx,nfplt,nf)
              endif
            enddo
            write(6,"('tracsat getting new history ',3i4,' for satut=',
     +        f9.3,' ipt=',i4,' it0,1=',2i3,' gcmut0,1=',2f9.3)") 
     +        (mt(i,it1),i=1,3),satut(ipt),ipt,it0,it1,gcmut0,gcmut1
            call gettgcm(histvol,nhvols,istart,mt(1,it1),pgcm1,
     +        ifget,nfhist,imx,kmx1,jmx,nfget1,iden,luhist,iprhdr,ivol,
     +        isdyn,istimes,ier)
            if (ier.ne.0) then
              write(6,"('>>> tracsat: error ',i3,' from gettgcm: it=',i3,
     +          ' mt=',i2,':',i2,':',i2)") 
     +          ier,it,(mt(i,it1),i=1,3)
              stop 'gettgcm'
            endif
            if (kmx1.ne.kmx.or.nfget1.ne.nfget) then
              write(6,"('>>> tracsat: different models?? kmx,kmx1=',
     +          2i3,' nfplt1,nfplt=',2i3)") kmx,kmx1,nfget1,nfget
              stop 'gettgcm'
            endif
            if (satut(ipt).ge.gcmut0.and.satut(ipt).le.gcmut1.and.
     +        ipt.gt.1.and.satut(ipt).lt.satutp) satutp = satut(ipt)
            goto 200
          endif ! need new gcm time
        endif ! inttime > 0
c       write(6,"('ipt=',i3,' satut=',f7.3,' gcmut0,1=',2f5.1,
c    +    ' satlat,lon,alt=',f8.3,f8.2,f7.2)") ipt,satut(ipt),
c    +    gcmut0,gcmut1,satlat(ipt),satlon(ipt),satalt(ipt)
c
c mtsat = model time at sat ut for this point:
        call dhrstomt(float(mt(1,it0))*24.+satut(ipt),
     +    mtsat(1,ipt),mtsat(2,ipt),mtsat(3,ipt))
c
c Save interpolated height column in gcmht:
c
        call locint(satlat(ipt),satlon(ipt),gcm0,imx,kmx,jmx,
     +    nfplt,ixcolz,gcmlat,gcmlon,col0,ier)
        if (inttime.gt.0)
     +    call locint(satlat(ipt),satlon(ipt),gcm1,imx,kmx,jmx,
     +      nfplt,ixcolz,gcmlat,gcmlon,col1,ier)
        call utint(satut(ipt),gcmut0,gcmut1,col0,col1,gcmht,kmx)
c       write(6,"('   ixcolz=',i2,' gcmht=',/(8f8.2))") ixcolz,
c    +    (gcmht(k),k=1,kmx)
c
c Fields interpolation loop:
c     subroutine locint(slat,slon,gcmin,imx,kmx,jmx,nf,ip,gcmlat,
c    +  gcmlon,gcmout,ier)
c
        nf = 0
        do ip=1,nftot
          if (ifplt(ip).gt.0) then
            nf = nf+1
            if (ip.le.nfhist) then
              call locint(satlat(ipt),satlon(ipt),gcm0,imx,kmx,jmx,
     +          nfplt,nf,gcmlat,gcmlon,col0,ier)
              if (inttime.gt.0) 
     +          call locint(satlat(ipt),satlon(ipt),gcm1,imx,kmx,jmx,
     +            nfplt,nf,gcmlat,gcmlon,col1,ier)
            else ! calculate field from history fields
              call mkfield(gcm0,imx,kmx,jmx,nfplt,ip,glb0)
              call locint(satlat(ipt),satlon(ipt),glb0,imx,kmx,jmx,
     +          1,1,gcmlat,gcmlon,col0,ier)
              if (inttime.gt.0) then
                call mkfield(gcm1,imx,kmx,jmx,nfplt,ip,glb1)
                call locint(satlat(ipt),satlon(ipt),glb1,imx,kmx,jmx,
     +            1,1,gcmlat,gcmlon,col1,ier)
              endif
            endif
            call utint(satut(ipt),gcmut0,gcmut1,col0,col1,col,kmx)
            call htint(col,gcmht,logint(ip),htscale,nhtscale,kmx,
     +        htcol,spval,ier)
c
c plot pointer gcmut(npts,nhtscale,nfplt):
c
            do k=1,nhtscale
              gcmut((nf-1)*nhtscale*npts+(k-1)*npts+ipt) = htcol(k)
            enddo
          endif ! ifplt(ip) > 0
        enddo ! ip=1,nftot
        satutp = satut(ipt)
      enddo ! ipt=1,npts
 300  continue
c
c Make plots (contour model in top part, map sat orbit in bottom part):
c (norbfr = number of orbits per frame. If norbfr <= 0, do all points
c  on same frame)
c
      if (norbfr.gt.0) then
        irec0 = 1
        prevlon = -999.
        norb = 0
        do i=1,npts
          glon = satlon(i)
          if (glon.lt.0.) glon = glon+360.
          if (glon.lt.prevlon) then ! end of an orbit
            norb = norb+1
            if (norb.eq.norbfr) then
              call pltut(gcmut,npts,irec0,i-1)
              irec0 = i
              norb = 0
            endif
          endif
          prevlon = glon
        enddo
        if (npts-irec0.gt.1) call pltut(gcmut,npts,irec0,npts)
      else
        call pltut(gcmut,npts,1,npts)
      endif
c
c Optionally write ascii file:
c
      if (iwrdat.gt.0) call wrdat(ludat,'tracsat.dat',satut,satlat,
     +  satlon,satalt,gcmut,npts,nhtscale,nfplt)
      call hpdeallc(pgcmut,ier,1)
      call hpdeallc(pgcm0,ier,1)
      if (inttime.gt.0) call hpdeallc(pgcm1,ier,1)
      call clsgks
      if (iwrdat.gt.0) then
        close(ludat)
        if (lenstr(unixdat).gt.0)
     +    call rcpfile(ludat,'tracsat.dat',unixdat)
      endif
      if (lenstr(unixcgm).ne.0) call rcpfile(0,'gmeta',unixcgm)
      goto 1000
      stop
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine xfergcm(gcm0,gcm1,imx,kmx,jmx,nfplt,nf)
      dimension gcm0(imx,kmx,jmx,nfplt),gcm1(imx,kmx,jmx,nfplt)
      gcm0(:,:,:,nf) = gcm1(:,:,:,nf)
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine utint(sut,ut0,ut1,col0,col1,col,kmx)
      dimension col0(kmx),col1(kmx),col(kmx)
c
      if (ut0.ge.0..and.ut1.lt.0.) then
        col(:) = col0(:)
      elseif (ut1.ge.0..and.ut0.lt.0.) then
        col(:) = col1(:)
      else
        frac = (sut-ut0) / (ut1-ut0)
        col(:) = frac*col1(:) + (1.-frac)*col0(:)
      endif
      return
      end