c
c------------------------------------------------------------------
c Begin file /home/sting/foster/uars/tracsat/tracsat.f
c------------------------------------------------------------------
c
      program tracsat
c
      include "timegcm.h"
      include "tracsat.h"
      include "uars.h"
      dimension gcm1(imx,kmx,jmx,nflds),gcmlat(jmx),gcmlon(imx),
     +  ixtime(nflds),gcmht(kmx),loght(nftot),ipltime(nftot),utlim(2)
      pointer(pgcm0,gcm0(imx,kmx,jmx,1)), (phtscale,htscale(1)), 
     +  (phtcol,htcol(1)), (pgcmut,gcmut(1))
      dimension iasf(13),col(kmx),col0(kmx),col1(kmx)
      dimension vtn(kmx),vte(kmx),vxo2(kmx),vxo(kmx),vxn2(kmx),
     +  vxo2p(kmx),vxne(kmx),itxver(7)
      character*8 lab(nftot)
      data iasf/13*1/
c
c     common/timeix/itxt,itxu,itxv,itxo2,itxox,itxn4s,itxnoz,itxco,
c    +              itxco2,itxh2o,itxh2,itxhox,itxop,itxch4,itxo21d,
c    +              itxno2,itxno,itxo3,itxo1,itxoh,itxho2,itxh,
c    +              itxn2d,itxti,itxte,itxne,itxo2p,itxw,itxz,itxn2
c
      data loght
c       tn   un   vn   o2    o  n4s   noz   co  co2  h2o  h2 hox
     +/  0,   0,   0,   1,   1,   1,    1,   1,   1,   1,  1,  1,
c       o+   ch4 o21d no2  no   o3 npo1   oh   ho2    h  n2d   ti
     +   1,    1,  1,   1,  1,   1,   1,   1,    1,   1,   1,   0,
c       te    ne  o2+    w   z   n2 ver63
     +   0,    1,   1,   0,  0,   1,    1/
      data lab
     + /'TN      ','UN      ','VN      ','O2      ','OX      ',  
     +  'N4S     ','NOZ     ','CO      ','CO2     ','H2O     ',  
     +  'H2      ','HOX     ','O+      ','CH4     ','O21D    ',
     +  'NO2     ','NO      ','O3      ','O1      ','OH      ',  
     +  'HO2     ','H       ','N2D     ','TI      ','TE      ',  
     +  'NE      ','O2+     ','W       ','Z       ','N2      ',
     +  'VER6300 '/
      data lusat/10/, lugcm/8/, spval/1.e36/
      equivalence(ixtime,itxt)
c
c Gcm grid:
c
      do j=1,jmx
        gcmlat(j) = rlat1 + (j-1)*dlat
      enddo
      do i=1,imx
        gcmlon(i) = rlon1 + (i-1)*dlon
      enddo
      do ip=1,nflds
        ixtime(ip) = ip
      enddo
c
c Need following fields for calculation of ver6300:
c itxt,itxte,itxo2,itxo1,itxn2,itxo2p,itxne
c
      itxver(1) = itxt
      itxver(2) = itxte
      itxver(3) = itxo2
      itxver(4) = itxo1
      itxver(5) = itxn2
      itxver(6) = itxo2p
      itxver(7) = itxne
c
c Get user input:
      call getinp(nvol,ntms,utlim,ipltime,inttime,itxver,norbfr,
     +  f107,f107a)
      write(6,"('tracsat: nvol=',i3,' ntms=',i3,' utlim=',2f9.3,
     +  ' norbfr=',i3)") nvol,ntms,utlim,norbfr
      do i=1,nvol
        write(6,"('tracsat: i=',i3,' nvol=',i3,' histvol=',a)") i,
     +    nvol,gcmvol(i)(1:lenstr(gcmvol(i)))
      enddo
c
c Prepare for graphics:
      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)
c
c Number of desired fields (for gcm1 heap allocation):
c (ixcolz will be index of heights in gcm0)
c
      nfplt = 0
c     do ip=1,nflds
      do ip=1,nftot
        if (ipltime(ip).gt.0) nfplt = nfplt+1
        if (ip.eq.itxz) ixcolz = nfplt
      enddo
      write(6,"('tracsat: nfplt=',i3,' ixcolz=',i3)") nfplt,ixcolz
c
c Linear height scale (y-axis):
c
      call hpalloc(phtscale,nhtscale,ier,1)
      if (ier.ne.0) then
        write(6,"('hpalloc error = ',i6,' for phtscale')") ier
        stop 'hpalloc'
      endif
      do k=1,nhtscale
        htscale(k) = ht1 + (k-1)*dht
      enddo
      write(6,"('tracsat: using linear height scale = ',/(5f9.2))")
     +  (htscale(k),k=1,nhtscale)
c
c Get satellite attitude:
c
      nrd = 0
 1000 continue
      nrd = nrd+1
      call rdapprox(satatt,lusat,utlim,satut,satlat,satlon,satalt,
     +  mxpts,npts,iydsat,norbfr,ier)
c
c Bracket first sat ut with tiegcm 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,"('Bracket 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
c Put first time temporarilly into gcm1 (the full array -- all fields):
c nrd > 1 means we have already read data once, and gotten first needed
c         gcm histories 
c
      if (nrd.gt.1) goto 2000
      istart = 1
      iprhdr = 1
      iden = 1	! number densities
      call gettime(gcmvol,nvol,istart,mt(1,it0),gcm1,imx,kmx,jmx,nflds,
     +  iden,lugcm,iprhdr,ivol,ier)
      if (ier.ne.0) then
        write(6,"('>>> tracsat: error ',i3,' from gettime for time ',
     +    3i4)") ier,(mt(i,it0),i=1,3)
        stop 'gettime'
      endif
      write(6,"('tracsat: obtained gcm for time ',3i4)") 
     +  (mt(i,it0),i=1,3)
c
c Allocate the partial array gcm0 (space only for desired fields),
c and transfer from gcm1:
c
      call hpalloc(pgcm0,imx*kmx*jmx*nfplt,ier,1)
      if (ier.ne.0) then
        write(6,"('hpalloc error = ',i6,' for pgcm1')") ier
        stop 'hpalloc'
      endif
      nf = 0
      do ip=1,nflds
        if (ipltime(ip).gt.0) then
          nf = nf+1
          gcm0(:,:,:,nf) = gcm1(:,:,:,ip)
        endif
      enddo
c
c Put second time into gcm1 if doing time interp:
c
      istart = 0
      if (inttime.gt.0) then
        call gettime(gcmvol,nvol,istart,mt(1,it1),gcm1,imx,kmx,jmx,
     +    nflds,iden,lugcm,iprhdr,ivol,ier)
        if (ier.ne.0) then
          write(6,"('>>> tracsat: error ',i3,' from gettime for time ',
     +      3i4)") ier,(mt(i,it1),i=1,3)
          stop 'gettime'
        endif
        write(6,"('tracsat: obtained gcm for time ',3i4)") 
     +    (mt(i,it1),i=1,3)
      endif
      iprhdr = 0
c
c Allocate space for htcol(nhtscale):
c
      call hpalloc(phtcol,nhtscale,ier,1)
      if (ier.ne.0) then
        write(6,"('hpalloc error = ',i6,' for phtcol')") ier
        stop 'hpalloc'
      endif
 2000 continue
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
      npoints = 0
      do ipt=1,npts
        npoints = npoints+1
c
c Get new gcm history if needed:
c (the second part of the conditional is meant to cover the case where
c  sat ut crosses day boundary (23 to 0), but has not yet been tested)
c Note if there are large gaps in sat ut it may be more efficient to
c  bracket the new satut, then acquire two new gcm histories if necessary,
c  rather than transferring ut0=ut1 and it1++, then getting the next history
c  in sequence -- this is the 200 loop.
c
c If we are not doing time interpolation, we'll use only the first
c   history (already acquired above), i.e., no need to acquire 
c   additional histories 
c
        if (inttime.le.0) goto 400	! no time interpolation
 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,nflds
            if (ipltime(ip).gt.0) then
              nf = nf+1
              gcm0(:,:,:,nf) = gcm1(:,:,:,ip)
            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 gettime(gcmvol,nvol,istart,mt(1,it1),gcm1,imx,kmx,jmx,
     +      nflds,iden,lugcm,iprhdr,ivol,ier)
          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
c
c Make sure satut is in between gcmut0 and gcmut1:
c
        if (satut(ipt).lt.gcmut0.or.satut(ipt).gt.gcmut1) then
          write(6,"('>>> tracsat: satut not bracketed: ipt=',i4,
     +      ' satut=',f9.3,' gcmut0,1=',2f9.3)") ipt,satut(ipt),
     +      gcmut0,gcmut1
          stop 'ut'
        endif
c
c Interpolate in lat,lon,ut to get pressure columns:
c
 400    continue
c
c First save interpolated height column in gcmht: 
c (to be used in ht interp routine htint)
c
        call locint(satlat(ipt),satlon(ipt),gcm0(1,1,1,ixcolz),
     +    imx,kmx,jmx,gcmlat,gcmlon,col0,ier)
        if (inttime.gt.0)
     +    call locint(satlat(ipt),satlon(ipt),gcm1(1,1,1,itxz),
     +      imx,kmx,jmx,gcmlat,gcmlon,col1,ier)
        call utint(satut(ipt),gcmut0,gcmut1,col0,col1,gcmht,kmx)
c
c History fields interpolation loop:
c
        nf = 0
        do ip=1,nflds
          if (ipltime(ip).gt.0) then
            nf = nf+1
c
c Do location and time interpolation:
c
            call locint(satlat(ipt),satlon(ipt),gcm0(1,1,1,nf),
     +        imx,kmx,jmx,gcmlat,gcmlon,col0,ier)
            if (inttime.gt.0) 
     +        call locint(satlat(ipt),satlon(ipt),gcm1(1,1,1,ip),
     +          imx,kmx,jmx,gcmlat,gcmlon,col1,ier)
            call utint(satut(ipt),gcmut0,gcmut1,col0,col1,col,kmx)
c
c Save certain pressure columns (interpolated to lat,lon,ut) 
c for use below (e.g., calculation of vol emission rates)
c (getinp has insured that all of these have been set in ipltime)
c
            if (ipltime(nflds+1).gt.0) then
              if (ip.eq.itxt) vtn(:) = col(:)
              if (ip.eq.itxte) vte(:) = col(:)
              if (ip.eq.itxo2) vxo2(:) = col(:)
              if (ip.eq.itxo1) vxo(:) = col(:)
              if (ip.eq.itxn2) vxn2(:) = col(:)
              if (ip.eq.itxo2p) vxo2p(:) = col(:)
              if (ip.eq.itxne) vxne(:) = col(:)
            endif
c
c Do height interpolation (linear or log):
c (gcmut is pointee, and will be final plot array)
c
            call htint(col,gcmht,loght(ip),htscale,nhtscale,kmx,
     +        htcol,spval,ier)
            do k=1,nhtscale
              gcmut((nf-1)*nhtscale*npts+(k-1)*npts+ipt) = htcol(k)
            enddo
          endif
        enddo
c
c Find additional fields not on history:
c (as of 2/20/92, this is only vol emission rate at 6300A)
c
        do ip=nflds+1,nftot
          if (ipltime(ip).gt.0) then
            iip = ip-nflds
            nf = nf+1 ! carried over from above history fields loop
            if (iip.eq.1) then ! ver6300
              slt = satut(ipt) + satlon(ipt) / 15.
              if (slt.lt.0.) slt = 24.-slt
              day = iydsat - iydsat/1000*1000
              call mkver63(gcmht,vtn,vte,vxo2,vxo,vxn2,vxo2p,vxne,
     +          col,kmx,1,slt,day,f107,f107a,satlat(ipt),satlon(ipt))
              logint = 1
c
c (to add extra fields: increase nftot (in tracsat.h), add new ipltime 
c  flag in input file, and calculate field here according to iip)
c
            endif
            call htint(col,gcmht,logint,htscale,nhtscale,kmx,htcol,
     +        spval,ier)
            do k=1,nhtscale
              gcmut((nf-1)*nhtscale*npts+(k-1)*npts+ipt) = htcol(k)
            enddo
          endif	! ipltime 
        enddo ! extra fields
c
c i=1,npts:
        satutp = satut(ipt)
      enddo
 300  continue
c
c Make plots (contour model in top part, map sat orbit in bottom part):
c Return for another orbit(s), if necessary according to input values
c   norbfr and utlim:
c
      call pltut(gcmut,htscale,npoints,nhtscale,nfplt,nftot,loght,lab,
     +  ipltime,gcmvol(ivol),inttime,cint,cmin,cmax,iydmod)
      if (satut(npts).lt.utlim(2)-.001) then
        write(6,"('tracsat: returning to rdapprox: satut(npts)=',
     +    f12.7,' utlim(2)=',f12.7)") satut(npts),utlim(2)
        call hpdeallc(pgcmut,ier,1)
        if (ier.ne.0) then
          write(6,"('>>> tracsat warning: hpdeallc error ',i3,
     +      ' for pgcmut')") ier
        endif
        goto 1000 
      endif
      write(6,"('tracsat quitting with npts=',i4,' satut(npts)=',f12.7,
     +  ' utlim=',2f12.7)") npts,satut(npts),utlim
c
c Return memory to heap:
c
      call hpdeallc(pgcm0,ier,1)
      if (ier.ne.0) then
        write(6,"('>>> tracsat warning: hpdeallc error ',i3,
     +    ' for pgcm0')") ier
      endif
      call hpdeallc(phtcol,ier,1)
      if (ier.ne.0) then
        write(6,"('>>> tracsat warning: hpdeallc error ',i3,
     +    ' for phtcol')") ier
      endif
      call hpdeallc(pgcmut,ier,1)
      if (ier.ne.0) then
        write(6,"('>>> tracsat warning: hpdeallc error ',i3,
     +    ' for pgcmut')") ier
      endif
      call hpdeallc(phtscale,ier,1)
      if (ier.ne.0) then
        write(6,"('>>> tracsat warning: hpdeallc error ',i3,
     +    ' for phtscale')") ier
      endif
c
      call clsgks
      stop
      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