c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Routines to do linear or log interpolation of gcm fields to 
c a linear height scale. Dimensions of fin, fht, and fout vary
c depending on routine called.
c
c glbhtint() is global height interpolator:
c   fin and fht (nlon,nzp,nlat) = global input and gcm heights fields
c   fout(nlon,nlat,nhts) = global interpolated field returned
c
c glbhtin() is same as glbhtint except that the intput fields fin and fht
c   are dimensioned (nlon,nlat,nzp) rather than (nlon,nzp,nlat). This is
c   more convenient for some codes, e.g. "super_procs" gcm8proc, gcm9proc,
c   etc.
c
c cuthtint() is height interpolator at some slice of the grid (e.g.,
c   along a latitude (idim1 = imx) or longitude (idim1 = jmx)
c   fin and fht (idim1,nzp) = input and gcm heights fields (at some cut)
c   fout(idim1,nhts) = interpolated field returned
c
c intloc() interpolates fields at some location (lat,lon gcm grid point)
c   
c Inputs common to all routines (dimensions vary):
c
c   fin: the gcm field at constant pressures
c   fht: gcm heights at constant pressures
c   hts: the height scale at which to do the interpolation
c   nhts: number of heights in hts
c   logint: if logint > 0, do log (exponential) interpolation, 
c           otherwise linear  
c   spval: if != 0., is special value returned in fout if 
c          ht(i) > max gcm height, or ht(i) < min gcm height. 
c          if = 0., extrapolation will be attempted if ht(i) out of range
c   iprnt: print stuff out if > 0 (e.g., if spval is used)
c
c Outputs common to all routines (dimension of fout varies):
c
c   fout: returned field, interpolated to desired height scale
c   ier: non-zero if an error has occurred
c
c Required routines:
c   fminmax(), bracket()
c
c Ben Foster (ncar/hao)
c foster@ncar.ucar.edu
c 303-497-1595
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)
c
      do j=1,nlat
        do i=1,nlon
          col(:) = fht(i,:,j)
          call fminmax(col,nzp,htmin,htmax)
          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
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      subroutine glbhtin(fin,fht,nlon,nlat,nzp,fout,hts,nhts,
     +                    logint,spval,ier,iprnt)
c
c 2/2/93: Same as glbhtint, except that fin and fht are dimensioned
c         (imx,jmx,kmx) rather than (imx,kmx,jmx). This is preferred
c         by "super-procs" (i.e., gcm8proc, gcm9proc, etc)
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,nlat,nzp) = global input field
c   fht(nlon,nlat,nzp) = global heights field
c   nlon,nlat,nzp = 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,nlat,nzp),fht(nlon,nlat,nzp),
     +  fout(nlon,nlat,nhts),hts(nhts)
      dimension col(nzp)
c
      do j=1,nlat
        do i=1,nlon
          col(:) = fht(i,j,:)
          call fminmax(col,nzp,htmin,htmax)
          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,"('glbhtin: 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,j,k1).gt.0..and.fin(i,j,k1).ne.spval.and.
     +            fin(i,j,k0).gt.0..and.fin(i,j,k0).ne.spval) then
                exparg = (alog(fin(i,j,k1) / fin(i,j,k0)) /
     +                   (col(k1)-col(k0))) * (hts(k)-col(k0))
                fout(i,j,k) = fin(i,j,k0)*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,j,k1) + (1.-f1)*fin(i,j,k0) 
            endif
 100        continue
          enddo
        enddo
      enddo
      return 
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine pglbhtint(pfin,fht,nlon,nzp,nlat,fout,hts,nhts,
     +                     logint,spval,ier,iprnt)
c
c Same as glbhtint, except that input field (1st arg) is a pointer
c rather than an array.
c
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
      pointer (pfin,fin(nlon,nzp,nlat))
      real fht(nlon,nzp,nlat),fout(nlon,nlat,nhts),hts(nhts)
      dimension col(nzp)
c
      do j=1,nlat
        do i=1,nlon
          col(:) = fht(i,:,j)
          call fminmax(col,nzp,htmin,htmax)
          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,"('pglbhtint: 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
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine cuthtint(fin,fht,idim1,nzp,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(idim1,nzp) = input field
c   fht(idim1,nzp) = heights field
c   idim1,nzp = 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(idim1,nhts) = contains interpolated field
c   ier = non-zero if an error has occurred
c
      dimension fin(idim1,nzp),fht(idim1,nzp),fout(idim1,nhts),
     +  hts(nhts)
      dimension col(nzp)
c
      do n=1,idim1
        col(:) = fht(n,:)
        call fminmax(col,nzp,htmin,htmax)
        do k=1,nhts
          if (hts(k).lt.htmin.or.hts(k).gt.htmax) then
            if (spval.ne.0.) then
              fout(n,k) = spval
              goto 100
            else
              k0 = nzp-1
              k1 = nzp
            endif
            if (iprnt.gt.0) write(6,"('cuthtint: height ',f5.1,
     +        ' out of range: htmin,max=',2f5.1,' k=',i2,' n=',i2)")
     +        hts(k),htmin,htmax,k,n
          else
            call bracket(hts(k),col,nzp,1,k0,k1,ierr)
            if (ierr.ne.0) then
              if (iprnt.gt.0) write(6,"('cuthtint: error from',
     +          ' bracket: n=',i3,' k=',i3,' hts(k)=',f8.2,
     +          ' hts col=',/(6f10.3))") n,k,hts(k),col
              fout(n,k) = spval
              ier = 1
              goto 100
            endif
          endif
c
c Do linear or log interpolation:
c
          if (logint.gt.0) then
            if (fin(n,k1).gt.0..and.fin(n,k1).ne.spval.and.
     +          fin(n,k0).gt.0..and.fin(n,k0).ne.spval) then
              exparg = (alog(fin(n,k1) / fin(n,k0)) /
     +                 (col(k1)-col(k0))) * (hts(k)-col(k0))
              fout(n,k) = fin(n,k0)*exp(exparg)
            else
              fout(n,k) = spval
            endif
          else
            if (fin(n,k0).ne.spval.and.fin(n,k1).ne.spval) then
              f1 = (hts(k)-col(k0)) / (col(k1)-col(k0))
              fout(n,k) = f1*fin(n,k1) + (1.-f1)*fin(n,k0) 
            else
              fout(n,k) = spval
            endif
          endif
 100      continue
        enddo 	! k=1,nhts
      enddo	! n=1,idim1
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine intloc(gcmin,gcmht,kmx,hts,nhts,loght,
     +  gcmout,idim1,ndim1,idim2,ier,spval,iprnt)
c
      dimension gcmin(idim1,kmx),gcmht(idim1,kmx),
     +  gcmout(idim1,idim2),hts(nhts)
c
c Interpolate idim1 values at nhts hts (all values assumed
c   to be at some tgcm grid lat/lon)
c
c On input:
c   gcmin(idim1,kmx) = Tgcm parameter to interpolate
c   gcmht(idim1,kmx) = Tgcm heights
c   hts(nhts) = heights at which to do interpolation
c   loght = 0 or 1 for linear or log interpolation
c   spval = special value is returned in gcmout if 
c           hts(ih) > max height in gcmht or < min height in gcmht
c   iprnt = 0 or 1 for print out if spval is used (height out of range)
c
c On output:
c   gcmout(idim1,idim2) is returned with interpolated data
c   ier = 0 if not error has occurred, > 0 if error has occurred
c
c
c Check idim2:
      ier = 0
      if (idim2.lt.nhts) then
        write(6,"('>>> intloc: bad idim2=',i3,' (is < nhts=',i3,')')")
     +    idim2,nhts
        ier = 1
        return
      endif
c
c First dimension loop:
      do 100 id = 1,ndim1
c
c Find max height:
        htmax = -1.e36
        htmin = 1.e36
        do 110 k=1,kmx
          htmax = amax1(htmax,gcmht(id,k))
          htmin = amin1(htmin,gcmht(id,k))
 110    continue
c
c Height loop:
        do 120 ih=1,nhts
          if (hts(ih).gt.htmax.or.hts(ih).lt.htmin) then
            if (iprnt.gt.0) then
              write(6,"('>>> intloc: id ih=',2i4,' height=',e12.4,
     +        ' outside range: htmin max=',2e12.4,' returning ',e12.4)")
     +        id,ih,hts(ih),htmin,htmax,spval
            endif
            gcmout(id,ih) = spval
            ier = 2
            goto 120
c           stop
          endif
c
c Bracket desired height:
          do 130 k=1,kmx-1
            if (hts(ih).ge.gcmht(id,k).and.hts(ih).le.gcmht(id,k+1))then
              ialt1 = k
              ialt2 = k+1
              goto 131
            endif
 130      continue
          if (iprnt.gt.0)
     +      write(6,"('>>> intloc: could not find index to height=',
     +        e12.4,' returning spval')") hts(ih)
          gcmout(id,ih) = spval
          goto 120
 131      continue
c
c Do log interpolation:
          if (loght.eq.1) then
            if (gcmin(id,ialt1).eq.spval.or.gcmin(id,ialt2).eq.spval.or.
     +          gcmin(id,ialt2).le.1.e-20.or.gcmin(id,ialt1).le.1.e-20
     +          .or.gcmht(id,ialt2)-gcmht(id,ialt1).le.1.e-20) then
              gcmout(id,ih) = spval
              goto 120
            endif
            if (gcmin(id,ialt2).gt.0.) then
              exparg = (alog(gcmin(id,ialt2) / gcmin(id,ialt1)) /
     +                 (gcmht(id,ialt2)-gcmht(id,ialt1))) *
     +                 (hts(ih)-gcmht(id,ialt1))
              gcmout(id,ih) = gcmin(id,ialt1)*exp(exparg)
            else
              gcmout(id,ih) = spval
            endif
c
c Do linear interpolation:
          else
            if (gcmht(id,ialt2).eq.spval.or.gcmht(id,ialt1).eq.spval.or.
     +          gcmht(id,ialt2)-gcmht(id,ialt1).le.1.e-20) then
              gcmout(id,ih) = spval
              goto 120
            endif
            f1 = (hts(ih)-gcmht(id,ialt1)) / 
     +           (gcmht(id,ialt2)-gcmht(id,ialt1))
            gcmout(id,ih) = f1*gcmin(id,ialt2)+(1.-f1)*gcmin(id,ialt1)
          endif
 120    continue
c
c End first dimension loop:
 100  continue
      return
      end