c
      subroutine timeterp(d1,d2,d3,d4,n,ipos,isec_dat,isec_req,fout,
     +  iprnt,ier)
      implicit none
c
c Interpolate from fields d1,d2,d3,d4 (at 4 consecutive days,
c   and all at time isec_dat) to time isec_req, returning fout(n)
c On input:
c   d1(n),d2(n),d3(n),d4(n) = the field at 4 consecutive days
c                             (all are at time isec_dat)
c   ipos = 1,2,3 for interpolation between d1 & d2, d2 & d3, or
c          d3 & d4 (usually ipos = 2)
c   isec_dat = time of fields at each of the 4 days (secs)
c   isec_req = desired time for interpolation (secs)
c In output:
c   fout(n) is defined at desired isec_req
c   ier = 0 if no error
c
      integer,intent(in)  :: n,ipos,isec_dat,isec_req,iprnt
      integer,intent(out) :: ier
      real,intent(in)  :: d1(n),d2(n),d3(n),d4(n)
      real,intent(out) :: fout(n)
!
! Local:
      integer :: i,itim_dat(3),itim_req(3)
      integer :: ncalls=0 ! for print only
      real :: frac0,frac1
c
c Check position and times:
c
      ncalls = ncalls+1
      ier = 0
      if (ipos.lt.1.or.ipos.gt.3) then
        write(6,"('>>> timeterp: bad ipos =',i3,
     +    ' (must be 1, 2, or 3)')") ipos
        ier = 1
        return
      endif
      if (isec_dat.lt.0.or.isec_dat.gt.86399) then
        write(6,"('>>> timeterp: bad isec_dat=',3i4)") isec_dat
        ier = 2
        return
      endif
      if (isec_req.lt.0.or.isec_req.gt.86399) then
        write(6,"('>>> timeterp: bad isec_req=',3i4)") isec_req
        ier = 3
        return
      endif
      call isec2hms(isec_dat,itim_dat(1),itim_dat(2),itim_dat(3))
      call isec2hms(isec_req,itim_req(1),itim_req(2),itim_req(3))
c
c Do interpolation:
c (note isec_dat may be 12*3600 or 0*3600, i.e., data can be at
c  0 or 12 ut)
c
      frac0 = float(mod(isec_req-isec_dat+86400,86400))/86400.
      if (iprnt.gt.0) write(6,"('timeterp: ih:im:is (req)=',i2,':',
     +  i2,':',i2,' ih:im:is (dat)=',i2,':',i2,':',i2,' frac0=',
     +  f10.5,' ipos=',i2)") itim_req,itim_dat,frac0,ipos
      if (ipos.eq.1) then ! between d1 and d2
        do i=1,n
          fout(i) = -(((frac0-1.)*(frac0-2.)*(frac0-3.)*d1(i))/6.) +
     +               ((frac0*(frac0-2.)*(frac0-3.)*d2(i))/2.) -
     +               ((frac0*(frac0-1.)*(frac0-3.)*d3(i))/2.) +
     +               ((frac0*(frac0-1.)*(frac0-2.)*d4(i))/6.)
        enddo
      elseif (ipos.eq.2) then ! between d2 and d3
        frac1 = 1.-frac0
        do i=1,n
          fout(i) = -((frac0*(1.-frac0)*(2.-frac0)*d1(i))/6.) +
     +               (((1.+frac0)*(1.-frac0)*(2.-frac0)*d2(i))/2.) +
     +               (((1.+frac1)*(1.-frac1)*(2.-frac1)*d3(i))/2.) -
     +               ((frac1*(1.-frac1)*(2.-frac1)*d4(i))/6.)
        enddo
      elseif (ipos.eq.3) then ! between d3 and d4
        do i=1,n
          fout(i) = -(((frac0-1.)*(frac0-2.)*(frac0-3.)*d4(i))/6.) +
     +               ((frac0*(frac0-2.)*(frac0-3.)*d3(i))/2.) -
     +               ((frac0*(frac0-1.)*(frac0-3.)*d2(i))/2.) +
     +               ((frac0*(frac0-1.)*(frac0-2.)*d1(i))/6.)
        enddo
      endif
      return
      end