c
      subroutine htint(gcmin,gcmht,loght,hts,nhts,kmx,gcmout,
     +  spval,ier)
c
      dimension gcmin(kmx),gcmht(kmx),hts(nhts),gcmout(nhts)
c
c Find min,max of gcm heights:
      htmin = 1.e36
      htmax = -1.e36
      do k=1,kmx
        if (gcmht(k).lt.htmin) htmin = gcmht(k)        
        if (gcmht(k).gt.htmax) htmax = gcmht(k)        
      enddo
c
c Height loop:
      do ih=1,nhts
        if (hts(ih).gt.htmax.or.hts(ih).lt.htmin) then
          gcmout(ih) = spval
          goto 100
        endif
c
c Bracket desired height:
c     subroutine bracket(x,xx,nx,inc,n1,n2,ier)
c
        call bracket(hts(ih),gcmht,kmx,1,iht0,iht1,ier)
        if (ier.ne.0) then
          write(6,"('>>> htscale error ',i3,' from bracket for ht=',
     +      f9.3)") hts(ih)
          stop 'hts'
        endif
c
c Log interpolation:
        if (loght.eq.1) then
          exparg = (alog(gcmin(iht1) / gcmin(iht0)) /
     +             (gcmht(iht1)-gcmht(iht0))) * (hts(ih)-gcmht(iht0))
          gcmout(ih) = gcmin(iht0)*exp(exparg)
        else
          f1 = (hts(ih)-gcmht(iht0)) /
     +         (gcmht(iht1)-gcmht(iht0))
          gcmout(ih) = f1*gcmin(iht1)+(1.-f1)*gcmin(iht0)
        endif
 100    continue
      enddo

      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine locint(slat,slon,gcmin,imx,kmx,jmx,nf,ip,gcmlat,
     +  gcmlon,gcmout,ier)
c
c Linearly interpolate gcmin to slat,slon location, 
c returning the interpolated pressure column in gcmout(kmx)
c
      dimension gcmin(imx,kmx,jmx,nf),gcmlat(jmx),gcmlon(imx),
     +  gcmout(kmx)
c
c Bracket lat and lon in gcm grid:
c
      ier = 0
      call bracket(slat,gcmlat,jmx,1,lat0,lat1,ier)
      if (ier.ne.0) then
        write(6,"('>>> locint error ',i3,' from bracket for slat=',
     +    f9.3)") slat
        ier = 1
        return
      endif
      call bracket(slon,gcmlon,imx,1,lon0,lon1,ier)
      if (ier.ne.0) then
        write(6,"('>>> locint error ',i3,' from bracket for slon=',
     +    f9.3)") slon
        ier = 2
        return
      endif
c
      do k=1,kmx
        flon = (slon-gcmlon(lon0)) / (gcmlon(lon1)-gcmlon(lon0))
        flat = (slat-gcmlat(lat0)) / (gcmlat(lat1)-gcmlat(lat0))
        f1 = gcmin(lon0,k,lat0,ip)
        f2 = gcmin(lon1,k,lat0,ip)
        f3 = gcmin(lon0,k,lat1,ip)
        f4 = gcmin(lon1,k,lat1,ip)
        gcmout(k) = flon*(flat*f4+(1.-flat)*f2) + 
     +            (1.-flon)*(flat*f3+(1.-flat)*f1)     
      enddo
      return
      end