c
c------------------------------------------------------------------
c Begin file /home/sting/foster/uars/tracsat/tracsat.f
c------------------------------------------------------------------
c
      program tracsat
c
      include "tiegcm.h"
      include "tracsat.h"
      include "uars.h"
      dimension gcm0(imx,kmx,jmx,nflds),gcm1(imx,kmx,jmx,nflds),
     +  gcmcol(kmx,nflds),gcmlat(jmx),gcmlon(imx),htscale(nhtscale)
      dimension ixgcm(nflds),gcmht(kmx),loght(nflds),htcol(nhtscale),
     +  gcmut(mxpts,nhtscale,nflds)
c
c     ixt,ixu,ixv,ixo2,ixo1,ixn4s,ixno,ixop,ixn2d,ixti,ixte,ixne,
c     ixo2p,ixw,ixz,ixpot,ixn2
c
      data loght/0,0,0,1,1,1,1,1,1,0,0,1,1,0,0,1,1/
      data lusat/10/, lugcm/8/, spval/1.e36/
      equivalence(ixgcm,ixt)
c
c Get user input:
      call getinp(nvol,ntms)
      write(6,"('tracsat: nvol=',i3,' ntms=',i3)") nvol,ntms
      do i=1,nvol
        write(6,"('tracsat: i=',i3,' nvol=',i3,' histvol=',a)") i,
     +    nvol,gcmvol(i)(1:lenstr(gcmvol(i)))
      enddo
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
        ixgcm(ip) = ip
      enddo
c
c Linear height scale (y-axis):
c
      do k=1,nhtscale
        htscale(k) = ht1 + (k-1)*dht
      enddo
      write(6,"('tracsat: using linear height scale = ',/(5f9.2))")
     +  htscale
c
c Get satellite attitude:
c
      call rdapprox(satatt,lusat,satut,satlat,satlon,satalt,mxpts,npts,
     +  ier)
      write(6,"('tracsat: npts=',i5)") npts
      do i=1,npts
        if (i.eq.1.or.mod(i,200).eq.0) write(6,"('tracsat: i=',i4,
     +    ' satut=',f9.3,' lat=',f8.2,' lon=',f8.2,' alt=',f8.2)")
     +    i,satut(i),satlat(i),satlon(i),satalt(i)
      enddo
c
c Bracket first sat ut with tiegcm:
c
      do it=2,ntms
        gcmut0 = float(mt(2,it-1)) + float(mt(3,it-1))/60.
        gcmut1 = float(mt(2,it)) + float(mt(3,it))/60.
        if (satut(1).ge.gcmut0.and.satut(1).le.gcmut1) then
          it0 = it-1
          it1 = it
          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 Get first two gcm histories:
c
      istart = 1
      iprhdr = 1
      iden = 2
      call getgcm(gcmvol,nvol,istart,mt(1,it0),gcm0,imx,kmx,jmx,nflds,
     +  iden,lugcm,iprhdr,ivol,ier)
      if (ier.ne.0) then
        write(6,"('>>> tracsat: error ',i3,' from getgcm for time ',
     +    3i4)") ier,(mt(i,it0),i=1,3)
        stop 'getgcm'
      endif
      write(6,"('tracsat: obtained gcm for time ',3i4)") 
     +  (mt(i,it0),i=1,3)
      istart = 0
      call getgcm(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 getgcm for time ',
     +    3i4)") ier,(mt(i,it1),i=1,3)
        stop 'getgcm'
      endif
      write(6,"('tracsat: obtained gcm for time ',3i4)") 
     +  (mt(i,it1),i=1,3)
      iprhdr = 0
c
c Get and interpolate tiegcm at each data point:
c
c     do ipt=1,npts
      do ipt=1,5
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
 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.
          gcm0(:,:,:,:) = gcm1(:,:,:,:)
          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 getgcm(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
        call getcol(satlat(ipt),satlon(ipt),satut(ipt),gcm0,gcm1,
     +    imx,kmx,jmx,nflds,gcmut0,gcmut1,gcmcol,gcmlat,gcmlon,ier)
        gcmht(:) = gcmcol(:,ixz)
c
c Interpolate to linear height scale:
c
        do ip=1,nflds
          call htint(gcmcol(1,ip),gcmht,loght(ip),htscale,nhtscale,
     +      kmx,htcol,spval,ier)
          gcmut(ipt,:,ip) = htcol(:)
        enddo

        write(6,"('tracsat: ipt=',i4,' satut=',f9.3,' gcmut0,1=',
     +    2f9.3,' gcmcol(:,1)=',/(6e12.4))") 
     +    ipt,satut(ipt),gcmut0,gcmut1,(gcmcol(k,1),k=1,kmx)
        write(6,"('         gcmut=',/(6e12.4))")
     +    (gcmut(ipt,k,1),k=1,nhtscale)

        satutp = satut(ipt)
c
c i=1,npts:
      enddo
 300  continue
c
      stop
      end