subroutine locint(slat,slon,gcmin,imx,kmx,jmx,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),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)") ier,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)") ier,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)
        f2 = gcmin(lon1,k,lat0)
        f3 = gcmin(lon0,k,lat1)
        f4 = gcmin(lon1,k,lat1)
        gcmout(k) = flon*(flat*f4+(1.-flat)*f2) + 
     +            (1.-flon)*(flat*f3+(1.-flat)*f1)     
      enddo
      return
      end