c23456789012345678901234567890123456789012345678901234567890123456789012
c
c bilinear subroutine
c
      subroutine bilinear_sub(ni, nj, dlon, dlat, lon, lat, val, 
     &                        interp_lon, interp_lat, interp_val)
      implicit none
      integer ni, nj, ileft, iright, jbot, jtop
      real dlon, dlat, interp_lon, interp_lat
      real lon(ni), lat(nj), val(ni,nj)
      real wx, wy, ri, rj
      real interp_val
     
c
c calculate floating point i and j, and interpolation weights
c

      rj=(interp_lat-lat(1))/dlat+1
      ri=(interp_lon-lon(1))/dlon+1

      ileft=int(ri)
      iright=ileft+1

      jbot=int(rj)
      jtop=jbot+1

      wx=ri-int(ri)
      wy=rj-int(rj)

c
c perform bilinear interpolation, and return the value
c

      interp_val=(1.-wy)*(1.-wx)*val(ileft,jbot) + 
     &         (1.-wy)*wx*val(iright,jbot) +
     &         wy*(1.-wx)*val(ileft,jtop) +
     &         wy*wx*val(iright,jtop)

      return
      end