c
      program tgcmsat
c
c Read satellite data and interpolate tgcm along the satellite
c track at a given height:
c
      include "tgcmsat.h"
      pointer (pgcm0,gcm0(1)), (pgcm1,gcm1(1)), (pgcmout,gcmout(1))
      integer ifget1(nfhist),iasf(13),iframe
      real col(kmax),col0(kmax),col1(kmax),gcmht(kmax)
      character*16 xlab
      character*56 toplab
      data iasf/13*1/
c
c Get user input:
c
      do ip=1,nftot
        ixgcm(ip) = ip
      enddo
      call getinp
c
c Prepare gks:
c
      if (lineplt.gt.0) then
        call opngks
        call gsclip(0)
        call gsasf(iasf)
        call gsfais(1)
        call setag
        iframe = 0
      endif
c
c Get satellite data:
c
      call rddat
c     write(6,"('npts=',i4,' satut=',/(6f12.4))") npts,
c    +  (satut(i),i=1,npts)
c     write(6,"('npts=',i4,' satslt=',/(6f12.4))") npts,
c    +  (satslt(i),i=1,npts)
c
c Bracket first sat ut:
c
      if (inttime.gt.0) then
        do it=1,ntms-1
          gcmut0 = float(mtimes(2,it)) + float(mtimes(3,it))/60.
          gcmut1 = float(mtimes(2,it+1)) + float(mtimes(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,"('>>> tgcmsat: 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(mtimes(2,1)) + float(mtimes(3,1))/60.
        gcmut1 = -1.
        it0 = 1
        it1 = -1
      endif
c
c Put gcm fields at model time it0 into gcm0:
c
      iprhdr = 1
      istart = 1
      call gettgcm(histvols,nhvols,istart,mtimes(1,it0),pgcm0,ifget,
     +  nfhist,imx,kmx,jmx,nfget,iden,luhist,iprhdr,ivol,isdyn,
     +  istimes,ier)   
      if (ier.ne.0) then
        write(6,"('>>> Error ',i3,' from gettgcm: it0=',i3,
     +    ' mtimes=',i3,':',i2,':',i2)") ier,it0,(mtimes(i,it0),i=1,3)
        stop 'getimesf'
      endif
      nfplt = 0
      do i=1,nfhist
        if (fields(i).gt.0.and.ifget(i).gt.0) then
          nfplt = nfplt+1
          ifplt(i) = nfplt
        endif
      enddo
      write(6,"('tgcmsat: nfplt=',i2,' ifplt=',/(20i3))") nfplt,ifplt
      write(6,"('tgcmsat: nfget=',i2,' ifget=',/(20i3))") nfget,ifget
      call alloc(pgcmout,npts*nfplt)
c
c Set tgcm grid (note kmx was output from gettgcm depending on
c tiegcm or timesgcm history):
c
      call setgrid(gcmlat,glat1,dlat,jmx, gcmlon,glon1,dlon,imx, 
     +                   gcmzp,zp1,dzp,kmx)
c
c If doing time interpolation, put gcm fields at model time it1 into gcm1:
c
      istart = 0
      if (inttime.gt.0) then
        ifget1(:) = ifget(:)
        call gettgcm(histvols,nhvols,istart,mtimes(1,it1),pgcm1,ifget1,
     +    nfhist,imx,kmx,jmx,nflds,iden,luhist,iprhdr,ivol,isdyn,
     +    istimes,ier)   
        if (ier.ne.0) then
          write(6,"('>>> Error ',i3,' from gettgcm: it1=',i3,
     +      ' mtimes=',i3,':',i2,':',i2)") ier,it1,(mtimes(i,it1),i=1,3)
          stop 'getimesf'
        endif
        if (nflds.ne.nfget) then
          write(6,"('>>> nflds.ne.nfget: nflds=',i2,' nfget=',i2)")
     +      nflds,nfget
          ier = 1
        endif
        do i=1,nfhist
          if (ifget(i).ne.ifget1(i)) then
            write(6,"('>>> i=',i2,' ifget(i).ne.ifget1(i): ',
     +        ' ifget(i)=',i2,' ifget1(i)=',i2)") i,ifget(i),ifget1(i)
            ier = 1 
          endif
        enddo
        if (ier.ne.0) stop 'ifget'
      endif
      iprhdr = 0
c
c Interpolate at each sat point -- get new tgcm histories as necessary:
c
      npoints = 0
      do ipt=1,npts
        npoints = npoints+1
        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,"('>>> tgcmsat: out of gcm histories: ipt=',i4,
     +        ' satut=',f9.3,' gcmut0=',f9.3,' ntms=',i3)")
     +        ipt,satut(ipt),gcmut0,ntms
            npoints = npoints-1
            goto 300
          endif
          gcmut1 = float(mtimes(2,it1)) + float(mtimes(3,it1))/60.
          do i=1,imx*kmx*jmx*nfplt
            gcm0(i) = gcm1(i)
          enddo
          write(6,"('tracsat getting new history ',3i4,' for satut=',
     +      f9.3,' ipt=',i4,' it0,1=',2i3,' gcmut0,1=',2f9.3)") 
     +      (mtimes(i,it1),i=1,3),satut(ipt),ipt,it0,it1,gcmut0,gcmut1
          call hpdeallc(pgcm1,ier,1)
          call gettgcm(histvols,nhvols,istart,mtimes(1,it1),pgcm1,
     +      ifget1,nfhist,imx,kmx,jmx,nflds,iden,luhist,iprhdr,ivol,
     +      isdyn,istimes,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 Confirm that satut(ipt) 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       write(6,"('ipt=',i4,' satut=',f8.3,' gcmut0,1=',2f8.3)") 
c    +    ipt,satut(ipt),gcmut0,gcmut1
 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((ifget(ixz)-1)*imx*kmx*jmx+1),
     +    imx,kmx,jmx,gcmlat,gcmlon,col0,ier)
        if (inttime.gt.0)
     +    call locint(satlat(ipt),satlon(ipt),
     +      gcm1((ifget(ixz)-1)*imx*kmx*jmx+1),
     +      imx,kmx,jmx,gcmlat,gcmlon,col1,ier)
        call utint(satut(ipt),gcmut0,gcmut1,col0,col1,gcmht,kmx)
c
c Location and time interpolation:
c
        do ip=1,nftot
          if (ifplt(ip).gt.0) then
c
c Interpolate to satellite location:
c
            call locint(satlat(ipt),satlon(ipt),
     +        gcm0((ifget(ip)-1)*imx*kmx*jmx+1),imx,kmx,jmx,gcmlat,
     +          gcmlon,col0,ier)
            if (inttime.gt.0) 
     +        call locint(satlat(ipt),satlon(ipt),
     +          gcm1((ifget(ip)-1)*imx*kmx*jmx+1),imx,kmx,jmx,gcmlat,
     +          gcmlon,col1,ier)
c
c Interpolate to satellite time:
c
            call utint(satut(ipt),gcmut0,gcmut1,col0,col1,col,kmx)
c
c Interpolate to desired height:
c gcmout(npts,nfplt)
c
            call htint(col,gcmht,loght(ip),height,1,kmx,
     +        gcmout((ifplt(ip)-1)*npts+ipt),spval,ier)
          endif
        enddo	! ip=1,nftot
c
c End sat point loop:
c
        satutp = satut(ipt)
c       write(6,"('ipt=',i4,' satut=',f7.3,' satslt=',f7.3,
c    +    ' satlon=',f9.3)") ipt,satut(ipt),satslt(ipt),satlon(ipt)
      enddo	! ipt=1,npts (sat points)
      call hpdeallc(pgcm0,ier,1)
c
 300  continue
c
c Make line plots of gcmout(npts,nfplt):
c
c     write(6,"('npoints=',i4,' satut=',/(8f8.2))")
c    +  npoints,(satut(i),i=1,npoints)
c     write(6,"('npoints=',i4,' satslt=',/(8f8.2))")
c    +  npoints,(satslt(i),i=1,npoints)
c
      if (lineplt.gt.0) then
        write(toplab,"('TIMESGCM INTERPOLATED TO UARS LAT,LON: ',
     +    'HEIGHT=',f7.1,' KM')") height
        write(6,"(' ')")
        write(6,"('Making line plots: npoints=',i4)") npoints
        do ip=1,nftot
          if (ifplt(ip).gt.0) then
            iframe = iframe+1
            lflab = lenstr(flab(ip))
            call clearstr(xlab)
            xlab(1:8) = flab(ip)
            if (loght(ip).gt.0) then
              write(xlab,"('LOG10 ',a)") flab(ip)(1:lflab)
            endif
            call pltline(satut,gcmout((ifplt(ip)-1)*npts+1),npts,
     +        npoints,'SAT UT (HRS)',xlab(1:lenstr(xlab)),
     +        toplab(1:lenstr(toplab)),loght(ip),spval)
            write(6,"('Line frame ',i3,' field ',a)") iframe,
     +        xlab(1:lenstr(xlab))
          endif
        enddo
        call clsgks
      endif
c
c Write ascii file:
c
      if (iwrascii.gt.0) call wrascii(gcmout,npts,npoints,nfplt)
c
      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
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine pltline(x,y,mxh,nh,xlab,ylab,toplab,logplt,spval)
      dimension x(mxh),y(mxh)
      character*(*) toplab,xlab,ylab
c
c     write(6,"(' ')")
c     write(6,"('pltline: mxh=',i3,' nh=',i3,' logplt=',i2,' y=',
c    +  /(6e12.4))") mxh,nh,logplt,(y(i),i=1,nh)
      if (logplt.gt.0) call log10f(y,nh,1.e-20,spval)
      call anotat(xlab,ylab,0,1,-1,L)
      call ezxy(x,y,nh,toplab)
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine wrascii(fout,np,npoints,nf)
      include "tgcmsat.h"
      character*32 flnm
      character*8 lab(nf)
      data flnm/'tgcmsat.dat                     '/
      real fout(np,nf)
      integer unlink
c
      lflnm = lenstr(flnm)
      istat = unlink(flnm(1:lflnm))
      open(luascii,file=flnm(1:lflnm),status='NEW',err=900)
      write(6,"(' ')")
      write(6,"('Opened file ',a,' for ascii data')") flnm(1:lflnm)
c
c Two-line header:
c
      write(luascii,"('TIMESGCM DATA INTERPOLATED ALONG SAT ',
     +  'TRACK AT HEIGHT = ',f7.2,' KM')") height
      write(luascii,"(' SAT UT','  SAT LAT','  SAT LON',
     +  '  SAT SLT',$)")
      n = 0
      do ip=1,nftot
        if (ifplt(ip).gt.0) then
          n = n+1
          lab(n) = flab(ip)
        endif
      enddo
      do ip=1,nf-1
        write(luascii,"('  ',a8,'  ',$)") lab(ip)
      enddo
      write(luascii,"('  ',a8,'  ')") lab(nf)
c
c One line per point:
c
      do i=1,npoints
        write(luascii,"(f7.3,3f9.2,$)") satut(i),satlat(i),
     +    satlon(i),satslt(i)
        do ip=1,nfplt-1
          write(luascii,"(e12.4,$)") fout(i,ip)
        enddo
        write(luascii,"(e12.4)") fout(i,nf)
      enddo
      close(luascii)
      write(6,"('Completed writing file ',a,' (npoints=',i4,')')")
     +  flnm(1:lflnm),npoints
      return
 900  continue
      write(6,"('>>> Error opening ascii file to receive data')")
      return
      end