!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   INTERPOLATE
!!!
!!!   Uses the 3-d 'simplex method' for linear interpolation.  The routine
!!!   coarse_to_fine MUST be called before calling this routine, as the
!!!   interpolated values from the parent grid are assumed to exist
!!!   within the input CHILDGRID.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine interpolate(childgrid,order,dim,refinerat,xx,yy,zz,nx,ny,nz)

      use types
      use bh02

      real(SP), dimension(:,:,:)  :: childgrid
      integer(I4B)                :: order,dim,refinerat
      real(SP), dimension(:)      :: xx,yy,zz
      integer(I4B)                :: nx,ny,nz
!     LOCAL VARIABLES
      integer(I4B)                :: i,j,k,ii,jj,kk
      integer(I4B)                :: iimax,jjmax,kkmax
      real(SP), dimension(3)      :: ddx,delx
      integer(I4B) svert(6,6,3),nlow(3),simnum
      real(SP) constcoeff(6),coeff(6,3)


      do i = 1,nx-1,refinerat
      do j = 1,ny-1,refinerat
      do k = 1,nz-1,refinerat
         if(i + refinerat .gt. nx -1) then
            iimax = i + refinerat
         else
            iimax = i + refinerat - 1
         endif
         if(j + refinerat .gt. ny -1) then
            jjmax = j + refinerat
         else
            jjmax = j + refinerat - 1
         endif
         if(k + refinerat .gt. nz -1) then
            kkmax = k + refinerat
         else
            kkmax = k + refinerat - 1
         endif

!        *** put coefficient computations here ***

         nlow(1) = i
         nlow(2) = j
         nlow(3) = k
         Forall(ii=1:6,jj=1:6,kk=1:3) svert(ii,jj,kk) = nlow(kk)
!        SIMPLEX 1
         svert(1,1,1) = svert(1,1,1) + refinerat
         svert(1,1,2) = svert(1,1,2) + refinerat
         svert(1,1,3) = svert(1,1,3) + refinerat
         svert(1,2,2) = svert(1,2,2) + refinerat
         svert(1,2,3) = svert(1,2,3) + refinerat
         svert(1,3,2) = svert(1,3,2) + refinerat
         svert(1,3,3) = svert(1,3,3) + refinerat
         svert(1,4,3) = svert(1,4,3) + refinerat
         svert(1,5,3) = svert(1,5,3) + refinerat
!        SIMPLEX 2
         svert(2,1,1) = svert(2,1,1) + refinerat
         svert(2,1,2) = svert(2,1,2) + refinerat
         svert(2,1,3) = svert(2,1,3) + refinerat
         svert(2,2,2) = svert(2,2,2) + refinerat
         svert(2,2,3) = svert(2,2,3) + refinerat
         svert(2,3,2) = svert(2,3,2) + refinerat
         svert(2,5,2) = svert(2,5,2) + refinerat
         svert(2,5,3) = svert(2,5,3) + refinerat
         svert(2,6,2) = svert(2,6,2) + refinerat
!        SIMPLEX3
         svert(3,1,1) = svert(3,1,1) + refinerat
         svert(3,1,3) = svert(3,1,3) + refinerat
         svert(3,2,3) = svert(3,2,3) + refinerat
         svert(3,3,1) = svert(3,3,1) + refinerat
         svert(3,3,2) = svert(3,3,2) + refinerat
         svert(3,3,3) = svert(3,3,3) + refinerat
         svert(3,4,1) = svert(3,4,1) + refinerat
         svert(3,4,3) = svert(3,4,3) + refinerat
         svert(3,5,3) = svert(3,5,3) + refinerat
!        SIMPLEX4
         svert(4,1,1) = svert(4,1,1) + refinerat
         svert(4,3,1) = svert(4,3,1) + refinerat
         svert(4,3,2) = svert(4,3,2) + refinerat
         svert(4,3,3) = svert(4,3,3) + refinerat
         svert(4,4,1) = svert(4,4,1) + refinerat
         svert(4,4,3) = svert(4,4,3) + refinerat
         svert(4,5,1) = svert(4,5,1) + refinerat
         svert(4,5,3) = svert(4,5,3) + refinerat
         svert(4,6,1) = svert(4,6,1) + refinerat
!        SIMPLEX5
         svert(5,1,1) = svert(5,1,1) + refinerat
         svert(5,1,2) = svert(5,1,2) + refinerat
         svert(5,2,2) = svert(5,2,2) + refinerat
         svert(5,3,2) = svert(5,3,2) + refinerat
         svert(5,5,1) = svert(5,5,1) + refinerat
         svert(5,5,2) = svert(5,5,2) + refinerat
         svert(5,5,3) = svert(5,5,3) + refinerat
         svert(5,6,1) = svert(5,6,1) + refinerat
         svert(5,6,2) = svert(5,6,2) + refinerat
!        SIMPLEX6
         svert(6,1,1) = svert(6,1,1) + refinerat
         svert(6,3,1) = svert(6,3,1) + refinerat
         svert(6,3,2) = svert(6,3,2) + refinerat
         svert(6,4,1) = svert(6,4,1) + refinerat
         svert(6,5,1) = svert(6,5,1) + refinerat
         svert(6,5,2) = svert(6,5,2) + refinerat
         svert(6,5,3) = svert(6,5,3) + refinerat
         svert(6,6,1) = svert(6,6,1) + refinerat
         svert(6,6,2) = svert(6,6,2) + refinerat

         delx(1) = xx(i+refinerat) - xx(i)
         delx(2) = yy(j+refinerat) - yy(j)
         delx(3) = zz(k+refinerat) - zz(k)
         do simnum = 1,6
            constcoeff(simnum) = childgrid(i,j,k)
            coeff(simnum,1) = &
          childgrid(svert(simnum,1,1),svert(simnum,1,2),svert(simnum,1,3)) - &
          childgrid(svert(simnum,2,1),svert(simnum,2,2),svert(simnum,2,3))
            coeff(simnum,2) = &
          childgrid(svert(simnum,3,1),svert(simnum,3,2),svert(simnum,3,3)) - &
          childgrid(svert(simnum,4,1),svert(simnum,4,2),svert(simnum,4,3))
            coeff(simnum,3) = &
          childgrid(svert(simnum,5,1),svert(simnum,5,2),svert(simnum,5,3)) - &
          childgrid(svert(simnum,6,1),svert(simnum,6,2),svert(simnum,6,3))
         enddo

         do ii=i,iimax
         do jj=j,jjmax
         do kk=k,kkmax
            ddx(1) = xx(ii) - xx(i)
            ddx(2) = yy(jj) - yy(j)
            ddx(3) = zz(kk) - zz(k)
            if(ddx(1) .le. ddx(2)) then
               if(ddx(2) .le. ddx(3)) then
!                 SIMPLEX1 : dx <= dy <= dz
                  simnum = 1
               else
                  if(ddx(3) .le. ddx(1)) then
!                    SIMPLEX5 : dz <= dx <= dy
                     simnum = 5
                  else
!                    SIMPLEX2 : dx <= dz <= dy
                     simnum = 2
                  endif
               endif
            else
               if(ddx(3) .le. ddx(2)) then
!                 SIMPLEX6 : dz <= dy <= dx
                  simnum = 6
               else
                  if(ddx(3) .le. ddx(1)) then
!                    SIMPLEX4: dy <= dz <= dx
                     simnum = 4
                  else
!                    SIMPLEX3: dy <= dx <= dz
                     simnum = 3
                  endif
               endif
            endif
            childgrid(ii,jj,kk) = constcoeff(simnum) + &
               ddx(1) * coeff(simnum,1) / delx(1) + &
               ddx(2) * coeff(simnum,2) / delx(2) + &
               ddx(3) * coeff(simnum,3) / delx(3)

         enddo
         enddo
         enddo

      enddo
      enddo
      enddo

!     write(6,*) '------- INTERPOLATION DONE ------------'
      return
      end subroutine interpolate


