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,trueval
      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 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
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)

      trueval=(ri-1)*(rj-1)

      print*,'interpolated value = ', dataint
      print*,'true value = ',trueval

      end