
      subroutine locbilin(fin,nlon,nzp,nlat,gcmlat,gcmlon,glat,glon,
     +  fout,spval,iprnt)
c
c Do bilinear interpolation of global fin to location glat,glon, 
c returning interpolated column in fout(nzp). If there are problems,
c fout is set to spval. All other parameters are input, and are
c unchanged.
c (bracket.f is in libutil.a from ~foster/lib)
c
      dimension fin(nlon,nzp,nlat), gcmlat(nlat),gcmlon(nlon),
     +  fout(nzp)
c
      call bracket(glon,gcmlon,nlon,1,lon0,lon1,ier)
      if (ier.ne.0) then
        if (iprnt.gt.0) write(6,"('locbilin: bad glon=',f12.4)") glon
        fout(:) = spval
        return
      endif
      call bracket(glat,gcmlat,nlat,1,lat0,lat1,ier)
      if (ier.ne.0) then
        if (iprnt.gt.0) write(6,"('locbilin: bad glat=',f12.4)") glat
        fout(:) = spval
        return
      endif
      dlon = gcmlon(lon1)-gcmlon(lon0)
      if (dlon.le.1.e-20) then
        if (iprnt.gt.0) write(6,"('locbilin: bad dlon=',e12.4)") dlon
        fout(:) = spval
        return
      endif
      dlat = gcmlat(lat1)-gcmlat(lat0)
      if (dlat.le.1.e-20) then
        if (iprnt.gt.0) write(6,"('locbilin: bad dlat=',e12.4)") dlat
        fout(:) = spval
        return
      endif
      flon = (glon-gcmlon(lon0)) / dlon
      flat = (glat-gcmlat(lat0)) / dlat
c
      do k=1,nzp
        fout(k) = flon*(flat*fin(lon1,k,lat1)+
     +    (1.-flat)*fin(lon1,k,lat0))+(1.-flon)*(flat*fin(lon0,k,lat1)+
     +    (1.-flat)*fin(lon0,k,lat0))
      enddo
      return
      end
