subroutine svpvf(iunit,isv,mjdt,fmjdt,  ierr,recf,vecf)
      implicit real*8 (a-h,o-z)
      dimension recf(3),vecf(3),idsv(37)
      dimension iflag(9),x(9),y(9),z(9),vx(9),vy(9),vz(9)
      data initl/1/
      data lunit /-9/
      data lastsv/-99/
      data lread1/-999/
c...............................................................
c                                                              .
c  Subroutine - svpv   Author  - Ben Remondi  Date - Jan/10/83 .
c                                                              .
c                                                              .
c  Function - svpv has the task of returning the ecf position  .
c            and velocity of a specified gps satellite (isv)   .
c            at a specified time (mjdt,fmjdt). To do this it   .
c            must read the appropriate records of the specified.
c            unit device (iunit) and interpolate.              .
c                                                              .
c            An NGS-type direct access data set is necessary   .
c            for this subroutine to work.                      .
c                                                              .
c   inputs - iunit:  the unit number of the appropriate        .
c                         orbit file                           .
c            isv:   the satellite identifier                 .
c            mjdt:    the mod. jul. day of the request time    .
c            fmjdt:   the fractional day of the request time   .
c                                                              .
c                                                              .
c  outputs:  recf(i): ecf x,y,z at mjdt+fmjdt                  .
c            vecf(i): ecf x-dot,y-dot,z-dot at mjdt+fmjdt      .
c            ierr:    appropriate error messages               .
c                     ierr=0 for normal situation              .
c                     ierr=1 for "time out of bounds early"    .
c                     ierr=2 for "time out of bounds late"     .
c                     ierr=3 for idsv not available on this    .
c                            orbit file                        .
c                     ierr=4 for position identically 0.d0     .
c                                                              .
c  recommendation: After the call to svpv, the calling program .
c                  should go as follows:                       .
c                                                              .
c                       if(ierr.ne.0) call action(ierr)        .
c                                                              .
c                                                              .
c                  where action is tailored to what specific   .
c                  actions need to be performed as a result    .
c                  of the error messages returned.             .
c                                                              .
c...............................................................
      ierr=0
      if(lunit.ne.iunit) initl=1
      if(initl.eq.0) go to 100
      initl=0
      lunit=iunit
      lastsv=-99
      lread1=-999
c
c Extract key quantities from the header block.
c
      irn=1
      read(iunit,err=9000,rec=irn) jyear,imon,iday,
     +  ihr,imin,seci,deltat,mjds,fmjds,nepoch
      irn=2
      read(iunit,err=9000,rec=irn) numsv,(idsv(k),k=1,12)
      irn=3
      read(iunit,err=9000,rec=irn) (idsv(k),k=13,25)
      irn=4
      read(iunit,err=9000,rec=irn) (idsv(k),k=26,37)
      iorder=9
c
c From this point, on, deltat is identically 1 unit of time.
c
      dnorm=1.d0
      arcl=(nepoch-1)*dnorm
      iom1=iorder-1
      iod2=iorder/2
      half=iod2*dnorm
      perday=86400.d0/deltat
c
  100 continue
c
c If the request is outside of the ephemeris span, early,
c set ierr=1.
c
c first compute request time minus start time.
c
      tfroms=((mjdt-mjds)+(fmjdt-fmjds))*perday

c     write(6,1122) tfroms
c1122 format(' "svpv" tfroms ',f20.15)
      if(tfroms.lt.(-dnorm)) ierr=1
      if(tfroms.gt.(arcl+dnorm)) ierr=2

      if(ierr.ne.0) return
c
c Itype indicates whether the request time is normal, close to
c start time, or close to end time.
c
      itype=1
      if(tfroms.lt.half) itype=2
      if(tfroms.ge.(arcl-half)) itype=3
c
c Next determine the sequential number of isv based on the header
c block ordering.
c
      ksv=0
      do 200 k=1,numsv
      if(idsv(k).eq.isv) go to 201
  200 continue
      ierr=3
      return
  201 ksv=k
      go to (350,400,500) itype
  350 continue
c
c itype=1 follows
c
c Compute the reference epoch and the epoch number of the first
c of the iorder records that will be read.
c
      irefep=tfroms+1+dnorm*0.001d0
      trefno=irefep-1
      tzero=trefno-4
      iread1=4+k+(irefep-iod2-1)*numsv
      go to 600
  400 continue
c
c itype=2 follows.
c
      iread1=4+k
      tzero=0
      go to 600
  500 continue
c
c itype=3 follows.
c
      iread1=4+k+(nepoch-iom1-1)*numsv
      tzero=nepoch-9
  600 continue
      if((ksv.eq.lastsv).and.(lread1.eq.iread1)) go to 850
c
c Now that we know which records to fetch, go fetch them.
c
      do 700 m=1,iorder
      irn=iread1+numsv*(m-1)
      read(iunit,iostat=ios,err=9000,rec=irn)  iflag(m),
     +    x(m),y(m),z(m),vx(m),vy(m),vz(m)
      if(iflag(m).eq.1) ierr=4
c
c     write(6,4444) iflag(m),x(m),y(m),z(m),irn
c4444 format('iflag,x,y,z svpv',i2,3f18.9,i6,' rec')
c
  700 continue
c
c Compute normalized requested time: tnorm=0 corresponds to
c the time of the first of the iorder values.
c
  850 tnorm=tfroms-tzero
c     write(6,9889) tnorm
c9889 format('"svpv" tnorm ',f18.15)
c
c Now perform the interpolations on x,y,z,vx,vy,vz.
c
      call lgpv9(tnorm,x,y,z,vx,vy,vz,recf(1),recf(2),recf(3),
     +   vecf(1),vecf(2),vecf(3))
c     write(6,8888) recf,vecf
c8888 format(' "svpv" recf,vecf ',3f16.8,/,1x,3f16.8)
      return
 9000 continue
      write(6,9010)irn,ios
 9010 format(' A read error in subroutine svpv at recno= ',i10,
     +  ' ios= ',i10)
      stop
      end