subroutine timeterp(d1,d2,d3,d4,n,ipos,idtime,itime,fout,
     +  iprnt,ier)
c
c Interpolate from fields d1,d2,d3,d4 (at 4 consecutive days,
c   and all at time idtime) to time itime, 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 idtime)
c   ipos = 1,2,3 for interpolation between d1 & d2, d2 & d3, or
c          d3 & d4 (usually ipos = 2)
c   idtime(3) = hr,min,isec time of fields at each of the 4 days
c   itime(3) = desired time for interpolation
c In output:
c   fout(n) is defined at desired itime
c   ier = 0 if no error
c
      dimension d1(n),d2(n),d3(n),d4(n),fout(n),idtime(3),itime(3)
      data ncalls/0/ ! for print only
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 (idtime(1).lt.0.or.idtime(1).gt.23.or.     ! hrs
     +    idtime(2).lt.0.or.idtime(2).gt.59.or.     ! min
     +    idtime(3).lt.0.or.idtime(3).gt.59) then   ! sec
        write(6,"('>>> timeterp: bad idtime=',3i4)") idtime
        ier = 2
        return
      endif
      if (itime(1).lt.0.or.itime(1).gt.23.or.     ! hrs
     +    itime(2).lt.0.or.itime(2).gt.59.or.     ! min
     +    itime(3).lt.0.or.itime(3).gt.59) then   ! sec
        write(6,"('>>> timeterp: bad itime=',3i4)") itime
        ier = 3
        return
      endif
      idelsec = (itime(1)*3600+itime(2)*60+itime(3)) -
     +          (idtime(1)*3600+idtime(2)*60+idtime(3))
      if (idelsec.lt.0) then
        write(6,"('>>> timeterp: itime must be > idtime: ',
     +    ' itime=',3i3,' idtime=',3i3)") itime,idtime
        ier = 4
        return
      endif
c
c Do interpolation:
c
      frac0 = float(idelsec) / 86400.
      if (iprnt.gt.0) write(6,"('timeterp: itime=',i2,':',i2,':',i2,
     +  ' frac0=',f10.5,' ipos=',i2)") itime,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