C
C-------------------------------------------------------------------------------
C
      	subroutine	bilinear_interp (zin, nx, ny, x, y, zout)
C
	Implicit	NONE
C
      	integer*4	nx, ny		! dimensions of input array
      	real*4		zin(nx,ny)	! input array
      	real*4		x, y		! input position
      	real*4		zout		! interpolated value
C
      	integer*4	ix1, iy1, ix2, iy2
      	real*4		xx, yy, z1, z2, dx, dy, dz



      	xx = x
      	if (xx .gt. nx) xx = nx
      	if (xx .lt.  1) xx = 1
      	yy = y
      	if (yy .gt. ny) yy = ny
      	if (yy .lt.  1) yy = 1

      	ix1 = xx
      	iy1 = yy
      	ix2 = ix1
      	iy2 = iy1
      	if (ix2 .lt. nx) ix2 = ix1 + 1
      	if (iy2 .lt. ny) iy2 = iy1 + 1

      	dx = xx - float (ix1)
      	dy = yy - float (iy1)

      	dz = zin(ix2,iy1) - zin(ix1,iy1)
      	z1 = zin(ix1,iy1) + (dz * dx)

      	dz = zin(ix2,iy2) - zin(ix1,iy2)
      	z2 = zin(ix1,iy2) + (dz * dx)

      	dz   = z2 - z1
      	zout = z1 + (dz * dy)


      	return
      	end



C
C*******************************************************************************
C
      	subroutine	triangular_interp (zin, nx, ny, idir, x, y, zout)
C
	Implicit	NONE
C
      	integer*4	nx, ny		! dimensions of input array
      	real*4		zin(nx,ny)	! input array
      	integer*4	idir		! direction of interpolation
      	real*4		x, y		! input position
      	real*4		zout		! interpolated value
C
      	integer*4	ix1, iy1, ix2, iy2
      	real*4		xx, yy, dx, dy, z1, z2, z3, z4, za, zb, zc

      	real*4		fuzz
      	data		fuzz	/ 1.0e-5	/




      	xx = x
      	if (xx .gt. nx) xx = nx
      	if (xx .lt.  1) xx = 1
      	yy = y
      	if (yy .gt. ny) yy = ny
      	if (yy .lt.  1) yy = 1

      	ix1 = xx
      	iy1 = yy
      	ix2 = ix1
      	iy2 = iy1
      	if (ix2 .lt. nx) ix2 = ix1 + 1
      	if (iy2 .lt. ny) iy2 = iy1 + 1

      	dx = xx - float (ix1)
      	dy = yy - float (iy1)

      	if ((dx .le. fuzz) .and. (dy .le. fuzz)) then

      	    zout = zin(ix1,iy1)

      	else

      	    if (idir .gt. 0) then
      	        z1 = zin(ix1,iy1)
      	        z2 = zin(ix2,iy1)
      	        z3 = zin(ix2,iy2)
                z4 = zin(ix1,iy2)
      	    else
      	        z1 = dx
      	        dx = 1.0 - dy
      	        dy = z1
      	        z1 = zin(ix1,iy2)
      	        z2 = zin(ix1,iy1)
      	        z3 = zin(ix2,iy1)
      	        z4 = zin(ix2,iy2)
      	    end if

      	    if (dy .le. dx) then
      	        za = z1
      	        zb = z2
      	        zc = z3
      	    else
                dx = 1.0 - dx
      	        dy = 1.0 - dy
      	        za = z3
      	        zb = z4
      	        zc = z1
      	    end if

      	    z1 = za + ((zb - za) * dx)
      	    z2 = za + ((zc - za) * dx)

      	    zout = z1 + ((z2 - z1) * (dy / dx))

      	end if


      	return
      	end