c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timesdif/htint.f
c------------------------------------------------------------------
c
      subroutine htint(fld,pnt,htscale,nhtscale)
c
c Do height interpolation of fld (pert) and pnt (cntr) to linear
c height scale (htscale), returning interpolated values in fld and pnt 
c (2nd dimension)
c
      include 'gettime.h'
      include 'timesdif.h'
      dimension pnt(imx,kmx,jmx,ntimefld),fld(imx,kmx,jmx,nplfld)
      dimension wkglbht(imx,jmx),htscale(nhtscale)
      data iprnt/0/
c
c ixz = index to heights field in fld:
      ixz = 0
      do ip=1,ntimefld
        if (ipltime(ip).gt.0) then
          ixz = ixz+1
          if (ip.eq.itxz) goto 99
        endif
      enddo
 99   continue
c
c Field loop:
      ixp = 0
      do ip=1,ntimefld
        if (ipltime(ip).le.0.or.ip.eq.itxz) goto 100
        ixp = ixp+1
        if (ip.eq.itxz.and.iplz.le.0) goto 100
c
c Interpolate control and restore in pnt:
        do ih=1,nhtscale
          call glbhtint(pnt(1,1,1,ip),pnt(1,1,1,itxz),imx,kmx,jmx,
     +      wkglbht,htscale(ih),1,itimelog(ip),cpspval,ier,iprnt)
          pnt(:,ih,:,ip) = wkglbht(:,:)
        enddo
c
c Interpolate perturbed:
        do ih=1,nhtscale
          call glbhtint(fld(1,1,1,ixp),fld(1,1,1,ixz),imx,kmx,jmx,
     +      wkglbht,htscale(ih),1,itimelog(ip),cpspval,ier,iprnt)
c
c If glbhtint put cpspval in control or perturbed at current point,
c   (i.e., in wkglbht) then perturbed gets cpspval, and control gets 0., 
c   so that difference taken later will = cpspval, and will
c   be skipped by conpack. Otherwise, restore interpolation in fld:
c
          do j=1,jmx
            do i=1,imx
              if (pnt(i,ih,j,ip).eq.cpspval.or.
     +            wkglbht(i,j).eq.cpspval) then
                fld(i,ih,j,ixp) = cpspval
                pnt(i,ih,j,ixp) = 0.
              else
                fld(i,ih,j,ixp) = wkglbht(i,j)
              endif
            enddo
          enddo
        enddo
 100    continue
      enddo
c
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine glbhtint(fin,fht,nlon,nzp,nlat,fout,hts,nhts,
     +                    logint,spval,ier,iprnt)
c
c Do linear or log interpolation of fin to constant height surfaces
c   in hts(nhts). Interpolated field returned in fout.
c
c On input:
c   fin(nlon,nzp,nlat) = global input field
c   fht(nlon,nzp,nlat) = global heights field
c   nlon,nzp,nlat = dimensions of fin and fht
c   hts(nhts) = heights at which to interpolate
c   logint <= 0 -> do linear interpolation,
c                  otherwise do exponential (log) interpolation
c   spval: If spval is non-zero and hts(k) is not within fht range
c          at the current grid point, then fout gets spval at the
c          current grid point.
c          If spval=0 and hts(k) is out of range, then
c          fout will be extrapolated to hts(k)
c   iprnt = 1 -> some error msgs will be printed if necessary
c
c On output:
c   fout(nlon,nlat,nhts) = contains interpolated field
c   ier = non-zero if an error has occurred
c
      dimension fin(nlon,nzp,nlat),fht(nlon,nzp,nlat),
     +  fout(nlon,nlat,nhts),hts(nhts)
      dimension col(nzp)
      do j=1,nlat
        do i=1,nlon
          col(:) = fht(i,:,j)
          call fminmax(col,nzp,htmin,htmax,spval)
          do k=1,nhts
            if (hts(k).lt.htmin.or.hts(k).gt.htmax) then
              if (spval.ne.0.) then
                fout(i,j,k) = spval
                goto 100
              else
                k0 = nzp-1
                k1 = nzp
              endif
            else
              call bracket(hts(k),col,nzp,1,k0,k1,ierr)
              if (ierr.ne.0) then
                if (iprnt.gt.0) write(6,"('glbhtint: error from',
     +            ' bracket: i j=',2i3,' k=',i3,' hts(k)=',f8.2,
     +            ' hts col=',/(6f10.3))") i,j,k,hts(k),col
                fout(i,j,k) = spval
                goto 100
              endif
            endif
c
c Do linear or log interpolation:
c
            if (logint.gt.0) then
              if (fin(i,k1,j).gt.0..and.fin(i,k1,j).ne.spval.and.
     +            fin(i,k0,j).gt.0..and.fin(i,k0,j).ne.spval) then
                exparg = (alog(fin(i,k1,j) / fin(i,k0,j)) /
     +                   (col(k1)-col(k0))) * (hts(k)-col(k0))
                fout(i,j,k) = fin(i,k0,j)*exp(exparg)
              else
                fout(i,j,k) = spval
              endif
            else
              f1 = (hts(k)-col(k0)) / (col(k1)-col(k0))
              fout(i,j,k) = f1*fin(i,k1,j) + (1.-f1)*fin(i,k0,j)
            endif
 100        continue
          enddo
        enddo
      enddo
      return
      end