
      program rdhrdi
c
      include 'hrdi.h'
c
      character*48 datfile,ans
      logical isthere
      dimension iasf(13),un0(mxalt),vn0(mxalt)
      data iasf/13*1/
      data lu/20/, iframe/0/
      data datfile /
     + '/d/foster/uars/hrdi/uars154.dat                 '/
      data iborb,ieorb/1,1/
c
c Open data file:
c
      write(6,"('Enter data file path (<CR> = ',a,'): ',$)")
     +  datfile(1:lenstr(datfile))
      read(5,"(a)") ans
      if (ans(1:1).ne.' ') datfile = ans 
      inquire(file=datfile,exist=isthere) 
      if (.not.isthere) then
        write(6,"('Cannot find file ',a)") datfile
        stop 'nofile'
      endif
      open(lu,file=datfile,status='OLD',err=999)
      lenfile = lenstr(datfile)
      write(6,"('Opened file ',a,' (fort unit ',i2,')')") 
     +  datfile(1:lenfile),lu
c
c Get first and last orbits to plot:
c
      write(6,"(' ')")
      write(6,"('Enter first,last orbit(s) to plot (<CR> = all ',
     +  'orbits): ',$)")
      read(5,"(a)") ans
      if (ans(1:1).eq.' ') then
        iborb = 1
        ieorb = 999
      else
        read(ans,*) iborb,ieorb
      endif
      write(6,"('iborb,ieorb=',i3,',',i3)") iborb,ieorb
c
c Plot types:
c
      write(6,"(' ')")
      iplcon = 0
      write(6,"('Make contour plots? (y/n <CR>=y): ',$)")
      read(5,"(a)") ans
      if (ans(1:1).eq.' '.or.ans(1:1).eq.'y'.or.ans(1:1).eq.'Y') 
     +  iplcon = 1
      if (iplcon.gt.0) then
        write(6,"('Will make contour plots')")
      else
        write(6,"('Will NOT make contour plots')")
      endif
c
      write(6,"(' ')")
      iplvec = 0
      write(6,"('Make vector plots? (y/n <CR>=y): ',$)")
      read(5,"(a)") ans
      if (ans(1:1).eq.' '.or.ans(1:1).eq.'y'.or.ans(1:1).eq.'Y') 
     +  iplvec = 1
      if (iplvec.gt.0) then
        write(6,"('Will make vector plots')")
      else
        write(6,"('Will NOT make vector plots')")
      endif
c
c Set up gks:
c
      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
c     call cpsetup(cpspval)
c
c Read number of records and number of alts, and alt profile:
c
      read(lu,*) nrec, nalt
      if (nalt.gt.mxalt) then
        write(6,"('Too many alts=',i3,' mxalt=',i3)") nalt,mxalt
        stop 'nalt'
      endif 
      write(6,"('nrec=',i5,' nalt=',i3)") nrec,nalt
      read(lu,"(1x,10f8.1)") (alt(i),i=1,nalt)
      write(6,"('alt profile = ',/(8f8.1))") (alt(i),i=1,nalt)
c
c Read data:
c norb = current orbit number (starting at 1)
c nporb = number of points in current orbit
c orbit is defined by local time turnaround point 
c   (sltp = slt of previous point)
c
      norb = 0
      nporb = 0
      sltp = 0.
      do 100 ir = 1,nrec
        read(lu,"(1x,i8,4f12.2)",end=888) itime,rlat,rlon,rsza,rslt
        read(lu,"(1x,10f8.1)") (un0(i),i=1,nalt)
        read(lu,"(1x,10f8.1)") (vn0(i),i=1,nalt)
c
c Plot if end of current orbit or end of file:
c
        if (sltp-rslt.ge.1..or.ir.eq.nrec) then
          norb = norb+1
          if (norb.ge.iborb.and.norb.le.ieorb) then
            irb = ir-nporb
            ire = ir-1
            if (ir.eq.nrec) then
              nporb = nporb+1
              call defdat(nporb,itime,rlat,rlon,rsza,rslt,un0,vn0)
              irb = ir-nporb+1
              ire = ir
            endif
            write(6,"('Plot orbit ',i2,' nporb=',i3,' ut=',f5.2,
     +        ' to ',f5.2,' slt=',f5.2,' to ',f5.2,' ir=',i4,' to ',
     +        i4)") norb,nporb,ut(1),ut(nporb),slt(1),slt(nporb),irb,ire
            call pltorb(nporb)
          else
            write(6,"('Skipping orbit ',i2,' (iborb,ieorb=',2i3,')')")
     +        norb,iborb,ieorb
          endif
          nporb = 0
          sltp = 0.
c
c Define data for 1st point of next orbit (just read):
c
          if (ir.ne.nrec) then
            nporb = nporb + 1
            call defdat(nporb,itime,rlat,rlon,rsza,rslt,un0,vn0)
          endif
c
c Otherwise update arrays for current orbit:
c
        else
          nporb = nporb + 1
          call defdat(nporb,itime,rlat,rlon,rsza,rslt,un0,vn0)
          sltp = rslt
        endif
c       write(6,"('ir=',i4,' ut=',f8.4,' slt=',f8.4)")
c    +    ir,float(itime)/3.6e6,rslt
 100  continue
c
      call clsgks
      stop 'success!'
 999  write(6,"('Error opening file ',a)") datfile(1:lenfile)
      stop 'openerr'
 888  write(6,"('EOF encountered on file ',a,' unit ',i2)") 
     +  datfile(1:lenfile),lu 
      stop 'eof'
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine defdat(nporb,itime,rlat,rlon,rsza,rslt,un0,vn0)
      include 'hrdi.h'
      dimension un0(mxalt),vn0(mxalt)
c
      ut(nporb) = float(itime) / 3.6e+6
      glat(nporb) = rlat
      glon(nporb) = rlon
      if (glon(nporb).ge.180.) glon(nporb) = glon(nporb) - 360.
      sza(nporb) = rsza
      slt(nporb) = rslt
      do i=1,nalt
        un(nporb,i) = un0(i)
        vn(nporb,i) = vn0(i)
      enddo
      return
      end
