c23456789012345678901234567890123456789012345678901234567890123456789012
      implicit none
      integer ni,nj
      parameter(ni=361,nj=181)
      real lat(nj),lon(ni),val(ni,nj)
      real interp_lat, interp_lon
      real bilinear
      real dlat,dlon,dataint
      integer  i,j

c
c enter interpolation point
c

      interp_lat=33.33
      interp_lon=40.75

c
c define latitude and longitude grid in 0.5-deg spacing
c

      dlat=0.5
      dlon=0.5

      do j=1,nj
         lat(j)=(j-1)*dlat
      enddo

      do i=1,ni
         lon(i)=(i-1)*dlon
      enddo

c
c define 2d value being interpolated
c

      do j=1,nj
         do i=1,ni
            val(i,j)=(i-1)*(j-1)
         enddo
      enddo

c
c call bilinear function to get the interpolated value
c   

      dataint = bilinear(ni, nj, dlon, dlat, lon, lat, val, 
     &                    interp_lon, interp_lat)
      print*,'interpolated value = ', dataint
      end
      
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c bilinear function
c
      real function bilinear(ni, nj, dlon, dlat, lon, lat, val, 
     &                       interp_lon, interp_lat)
      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
     
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

      bilinear=(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