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 ri,rj,dlat,dlon,dataint
      real wx,wy
      integer  ileft,iright,jbot,jtop
      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 search for the longitude and latitude indices
c assign indices for 4 corners
c

      do j=1,nj
         if (interp_lat < lat(j)) exit
      enddo
      jtop=j
      jbot=j-1

      do i=1,ni
         if (interp_lon < lon(i)) exit
      enddo
      iright=i
      ileft=i-1

c
c calculate interpolation weights
c

      wx=(interp_lon-lon(ileft))/dlon
      wy=(interp_lat-lat(jbot))/dlat

c
c perform bilinear interpolation
c

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

      print*,'interpolated value = ', dataint

      end