C********************************************************



      subroutine lint( aa, n, z, ap )
      integer n
      double precision aa(n)
      double precision z,ap
c
c     Linearly interpolate the 1D array of data points 
c     represented by aa to the 1D space point z.
c
c     REFERENCE: 
c        Numerical Recipes p.116 et seq. 3D details worked 
c        out by RM.
c
c     ASSUMPTIONS:
c        z is in units of dz, the grid spacing. 
c        This permits quick location of the
c        relevant nearest neighbour interpolating points.
c
c     INPUT ARGUMENTS:
c
c        aa -- The 1D array representing the function aa(z).
c              The relationship of aa(k) to the function
c              aa(z) is: aa(z) = aa(k*dz)
c
c        n -- The z-dimension of array aa.
c
c        z -- The point to be interpolated. The units of x
c                 y and z are assumed to be dx, dy and dz, 
c                 respectively.
c
c     OUTPUT ARGUMENTS:
c
c        ap -- The interpolated value of aa(z)
c
      integer k,kp
      double precision v,omv
c
c     Find the indices of the 3D-cube in which the point (x,y,z)
c     lies. USE THE FACT THAT (x,y,z) IS IN UNITS OF dx,dy,dz,
c     the 3D grid spacings. That is, xtrue = x * dx, 
c     ytrue = y * dy, ztrue = z * dz.
c     We apply simple chopping to get the i, j, k!
c
c      write(*,*) 'Inside lint'
      k = z
c      write(*,*)'k,z = ',k,z
c
c     Define v
c
      v = z-dble(k)
      omv = 1.d0-v
c
c     adjust k for unit offset array
c
      k = k+1
      kp = k+1

c check periodic condition 
   
      if(k.gt.n)  kp=1
c      if(k.gt.n)  k=k-n
c      if(kp.gt.n) kp=kp-n

c
c     these lines are for debug purposes
c
c      if ( k .le. 0 ) print*, 'Error: k = ',k,' <= 0; z = ',z
c      if ( kp .gt. n ) print*, 'Error: k+1 = ',kp,' > n; z = ',z
c
c     The formula for the linear interpolation...
c
      ap = aa(kp)*v+aa(k)*omv
c      write(*,*)'ap,aa(kp),kp,aa(k),k = ',ap,aa(kp),kp,aa(k),k
c
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine lintarea( aa, l, m, x, y, ap )
      implicit double precision (a-h,o-z)
      integer l,m
      dimension aa(l,m)
      double precision t,x,y,z,ap
c
c     Linearly interpolate between the 3D arrays aa1 and aa2
c     in time, to the point t.
c     Trilinearly interpolate the 3D array of data points
c     represented by aa to the 3D space point (x,y,z).
c
c     REFERENCE:
c        Numerical Recipes p.116 et seq. 3D details worked
c        out by RM.
c
c     ASSUMPTIONS:
c        The t, x, y, and z are in units of dt, dx, dy and dz,
c        respectively. This permits quick location of the
c        relevant nearest neighbour interpolating points.
c
c     INPUT ARGUMENTS:
c
c        aa -- The 2D array representing the function
c               aa(t, x, y).
c        l,m -- The x- and y-dimensions of array aa,
c                 respectively.
c
c     OUTPUT ARGUMENTS:
c
c        ap -- The interpolated value of aa(t,x,y)
c
      integer i,j,ip,jp
      real*8 s,u,oms,omu
      real*8 a1,a2,a3,a4
c
c     Find the indices of the 2D-cube in which the point (x,y)
c     lies. USE THE FACT THAT (x,y) IS IN UNITS OF dx,dy,dz,
c     the 2D grid spacings. That is, xtrue = x * dx,
c     ytrue = y * dy.
c     We apply simple chopping to get the it, i, j!
c

      i =  x - ph0/delx
      j =  y - th0/dely
c       write(*,*) 'aa',aa
c      write(*,*)'x,y',x,y,'i,j',i,j
c
c     Define q, s, u, v
c

      s = x-dfloat(i)
      u = y-dfloat(j)


      oms = 1.d0-s
      omu = 1.d0-u

c      write(*,*)'s,u,oms,omu',s,u,oms,omu

c
c     adjust i, j, k for unit offset arrays
c
      i = i+1
      j = j+1
      ip = i+1
      jp = j+1


      if(i.gt.l)  i=i-l
      if(ip.gt.l) ip=ip-l
      if(j.gt.m)  j=j-m
      if(jp.gt.m) jp=jp-m
c      write(*,*) i,ip,j,jp
c
c     these lines are for debug purposes
c
c      if ( i .le. 0 ) print*, 'Error: i = ',i,' <= 0; x = ',x
c      if ( ip .gt. l ) print*, 'Error: i+1 = ',ip,' > l; x = ',x
c      if ( j .le. 0 ) print*, 'Error: j = ',j,' <= 0; y = ',y
c      if ( jp . gt. m ) print*, 'Error: j+1 = ',jp,' > m; y = ',y
c      if ( k .le. 0 ) print*, 'Error: k = ',k,' <= 0; z = ',z
c      if ( kp .gt. n ) print*, 'Error: k+1 = ',kp,' > n; z = ',z
c
c     First linearly interpolate in time the eight cube corners
c     (eight nearest grid points) surrounding the particle
c

c      write(*,*) 'aa(ip,jp)=',aa(ip,jp)
c      write(*,*) 'aa(i,jp)=',aa(i,jp)
c      write(*,*) 'aa(ip,j)=',aa(ip,j)
c      write(*,*) 'aa(i,j)=',aa(i,j)
      a1 = aa(ip,jp)
      a2 = aa(i,jp)
      a3 = aa(ip,j)
      a4 = aa(i,j)
c      write(*,*)'a1,a2,a3,a4',a1,a2,a3,a4
c
c     Now do the trilinear interpolation in space
c
      ap = a1*s*u+a2*oms*u
     &   + a3*s*omu+a4*oms*omu
c      write(*,*)'ap',ap
c
c
c     The formula for the trilinear interpolation...
c
c      ap = aa(ip,jp,kp)*s*u*v+aa(i,jp,kp)*oms*u*v
c     &   + aa(ip,j,kp)*s*omu*v+aa(i,j,kp)*oms*omu*v
c     &   + aa(ip,jp,k)*s*u*omv+aa(i,jp,k)*oms*u*omv
c     &   + aa(ip,j,k)*s*omu*omv+aa(i,j,k)*oms*omu*omv
c
      return
      end    
c************************************************************


      subroutine trilnt( aa, l, m, n, x, y, z, ap )
      integer l,m,n
      real aa(l,m,n)
      real*8 x,y,z
      real ap
c
c     Trilinearly interpolate the 3D array of data points 
c     represented by aa to the 3D space point (x,y,z).
c
c     REFERENCE: 
c        Numerical Recipes p.116 et seq. 3D details worked 
c        out by RM.
c
c     ASSUMPTIONS:
c        The x, y, and z are in units of dx, dy and dz, 
c        respectively. This permits quick location of the
c        relevant nearest neighbour interpolating points.
c
c     INPUT ARGUMENTS:
c
c        aa -- The 3D array representing the function aa(x,y,z).
c              The relationship of aa(i,j,k) to the function
c              aa(x,y,z) is: aa(x,y,z) = aa(i*dx, j*dy, k*dz)
c
c        l,m,n -- The x-, y- and z-dimensions of array aa,
c                 respectively.
c
c        x,y,z -- The point to be interpolated. The units of x
c                 y and z are assumed to be dx, dy and dz, 
c                 respectively.
c
c     OUTPUT ARGUMENTS:
c
c        ap -- The interpolated value of aa(x,y,z)
c
      integer i,j,k,ip,jp,kp
      real*8 t,u,v,omt,omu,omv
c
c     Find the indices of the 3D-cube in which the point (x,y,z)
c     lies. USE THE FACT THAT (x,y,z) IS IN UNITS OF dx,dy,dz,
c     the 3D grid spacings. That is, xtrue = x * dx, 
c     ytrue = y * dy, ztrue = z * dz.
c     We apply simple chopping to get the i, j, k!
c
      i = x
      j = y
      k = z

c     write(*,*) 'x,y,z',x,y,z
c     stop
c
c     Define t, u, v
c
      t = x-dfloat(i)
      u = y-dfloat(j)
      v = z-dfloat(k)
      omt = 1.d0-t
      omu = 1.d0-u
      omv = 1.d0-v
c
c     adjust i, j, k for unit offset arrays
c
      i = i+1
      j = j+1
      k = k+1
      ip = i+1
      jp = j+1
      kp = k+1
      if(ip.gt.l)ip=ip-l
      if(jp.gt.m)jp=jp-m
      if(kp.gt.n)kp=kp-n
c
c     these lines are for debug purposes
c
      if ( i .le. 0 ) print*, 'Error: i = ',i,' <= 0; x = ',x
      if ( ip .gt. l ) print*, 'Error: i+1 = ',ip,' > l; x = ',x
      if ( j .le. 0 ) print*, 'Error: j = ',j,' <= 0; y = ',y
      if ( jp . gt. m ) print*, 'Error: j+1 = ',jp,' > m; y = ',y
      if ( k .le. 0 ) print*, 'Error: k = ',k,' <= 0; z = ',z
      if ( kp .gt. n ) print*, 'Error: k+1 = ',kp,' > n; z = ',z
c
c     The formula for the trilinear interpolation...
c
      ap = aa(ip,jp,kp)*t*u*v+aa(i,jp,kp)*omt*u*v
     &   + aa(ip,j,kp)*t*omu*v+aa(i,j,kp)*omt*omu*v
     &   + aa(ip,jp,k)*t*u*omv+aa(i,jp,k)*omt*u*omv
     &   + aa(ip,j,k)*t*omu*omv+aa(i,j,k)*omt*omu*omv
c
      return
      end     

      subroutine lintarea2( aa, l, m, x, y, ap )
      implicit double precision (a-h,o-z)
      integer l,m
      dimension aa(l,m)
      double precision t,x,y,z,ap
c
c     Linearly interpolate between the 3D arrays aa1 and aa2
c     in time, to the point t.
c     Trilinearly interpolate the 3D array of data points
c     represented by aa to the 3D space point (x,y,z).
c
c     REFERENCE:
c        Numerical Recipes p.116 et seq. 3D details worked
c        out by RM.
c
c     ASSUMPTIONS:
c        The t, x, y, and z are in units of dt, dx, dy and dz,
c        respectively. This permits quick location of the
c        relevant nearest neighbour interpolating points.
c
c     INPUT ARGUMENTS:
c
c        aa -- The 2D array representing the function
c               aa(t, x, y).
c        l,m -- The x- and y-dimensions of array aa,
c                 respectively.
c
c     OUTPUT ARGUMENTS:
c
c        ap -- The interpolated value of aa(t,x,y)
c
      integer i,j,ip,jp
      real*8 s,u,oms,omu
      real*8 a1,a2,a3,a4
c
c     Find the indices of the 2D-cube in which the point (x,y)
c     lies. USE THE FACT THAT (x,y) IS IN UNITS OF dx,dy,dz,
c     the 2D grid spacings. That is, xtrue = x * dx,
c     ytrue = y * dy.
c     We apply simple chopping to get the it, i, j!
c

      i =  x
      j =  y
c       write(*,*) 'aa',aa
c      write(*,*)'x,y',x,y,'i,j',i,j
c
c     Define q, s, u, v
c

      s = x-dfloat(i)
      u = y-dfloat(j)


      oms = 1.d0-s
      omu = 1.d0-u

c      write(*,*)'s,u,oms,omu',s,u,oms,omu

c
c     adjust i, j, k for unit offset arrays
c
      i = i+1
      j = (m-j)+1
      ip = i+1
      jp = j+1


      if(i.gt.l)  i=i-l
      if(ip.gt.l) ip=ip-l
      if(j.gt.m)  j=j-m
      if(jp.gt.m) jp=jp-m
c      write(*,*) i,ip,j,jp
c
c     these lines are for debug purposes
c
c      if ( i .le. 0 ) print*, 'Error: i = ',i,' <= 0; x = ',x
c      if ( ip .gt. l ) print*, 'Error: i+1 = ',ip,' > l; x = ',x
c      if ( j .le. 0 ) print*, 'Error: j = ',j,' <= 0; y = ',y
c      if ( jp . gt. m ) print*, 'Error: j+1 = ',jp,' > m; y = ',y
c      if ( k .le. 0 ) print*, 'Error: k = ',k,' <= 0; z = ',z
c      if ( kp .gt. n ) print*, 'Error: k+1 = ',kp,' > n; z = ',z
c
c     First linearly interpolate in time the eight cube corners
c     (eight nearest grid points) surrounding the particle
c

c      write(*,*) 'aa(ip,jp)=',aa(ip,jp)
c      write(*,*) 'aa(i,jp)=',aa(i,jp)
c      write(*,*) 'aa(ip,j)=',aa(ip,j)
c      write(*,*) 'aa(i,j)=',aa(i,j)
      a1 = aa(ip,jp)
      a2 = aa(i,jp)
      a3 = aa(ip,j)
      a4 = aa(i,j)
c      write(*,*)'a1,a2,a3,a4',a1,a2,a3,a4
c
c     Now do the trilinear interpolation in space
c
      ap = a1*s*u+a2*oms*u
     &   + a3*s*omu+a4*oms*omu
c      write(*,*)'ap',ap
c
c
c     The formula for the trilinear interpolation...
c
c      ap = aa(ip,jp,kp)*s*u*v+aa(i,jp,kp)*oms*u*v
c     &   + aa(ip,j,kp)*s*omu*v+aa(i,j,kp)*oms*omu*v
c     &   + aa(ip,jp,k)*s*u*omv+aa(i,jp,k)*oms*u*omv
c     &   + aa(ip,j,k)*s*omu*omv+aa(i,j,k)*oms*omu*omv
c
      return
      end    
c************************************************************

