!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!   INTERPOLATE_BOUNDARY
!!!
!!!   Uses the 2-d 'simplex method' for linear interpolation on the
!!!   boundary of the input grid.  The routine coarse_to_fine_boundary
!!!   MUST be called before calling this routine, as the interpolated,
!!!   time averaged values from the parent grid are assumed to exist
!!!   on the boundary of the input CHILDGRID
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine interpolate_boundary(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) simnum
      real(SP) constcoeff,coeff(2,2)



!     Y-Z Faces

      do i = 1,nx,nx-1
      do j = 1,ny-1,refinerat
      do k = 1,nz-1,refinerat
         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 ***

         delx(2) = yy(j+refinerat) - yy(j)
         delx(3) = zz(k+refinerat) - zz(k)
         constcoeff = childgrid(i,j,k)
         coeff(1,1) = (childgrid(i,j+refinerat,k) - constcoeff)/delx(2)
         coeff(1,2) = (childgrid(i,j+refinerat,k+refinerat) - &
                       childgrid(i,j+refinerat,k))/delx(3)
         coeff(2,1) = (childgrid(i,j+refinerat,k+refinerat) - &
                       childgrid(i,j,k+refinerat))/delx(2)
         coeff(2,2) = (childgrid(i,j,k+refinerat) - constcoeff)/delx(3)

         ii = i
         do jj=j,jjmax
         do kk=k,kkmax
            ddx(2) = yy(jj) - yy(j)
            ddx(3) = zz(kk) - zz(k)
            if(ddx(3) .le. ddx(2)) then
!              SIMPLEX1 : dz <= dy 
               simnum = 1
            else
!              SIMPLEX2 : dy <= dz 
               simnum = 2
            endif
            childgrid(ii,jj,kk) = constcoeff + &
               ddx(2) * coeff(simnum,1)  + &
               ddx(3) * coeff(simnum,2) 

         enddo
         enddo

      enddo
      enddo
      enddo

!     X-Z Faces

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

!        *** put coefficient computations here ***

         delx(1) = xx(i+refinerat) - xx(i)
         delx(3) = zz(k+refinerat) - zz(k)
         constcoeff = childgrid(i,j,k)
         coeff(1,1) = (childgrid(i+refinerat,j,k) - constcoeff)/delx(1)
         coeff(1,2) = (childgrid(i+refinerat,j,k+refinerat) - &
                       childgrid(i+refinerat,j,k))/delx(3)
         coeff(2,1) = (childgrid(i+refinerat,j,k+refinerat) - &
                       childgrid(i,j,k+refinerat))/delx(1)
         coeff(2,2) = (childgrid(i,j,k+refinerat) - constcoeff)/delx(3)

         jj = j
         do ii=i,iimax
         do kk=k,kkmax
            ddx(1) = xx(ii) - xx(i)
            ddx(3) = zz(kk) - zz(k)
            if(ddx(3) .le. ddx(1)) then
!              SIMPLEX1 : dz <= dx 
               simnum = 1
            else
!              SIMPLEX2 : dx <= dz 
               simnum = 2
            endif
            childgrid(ii,jj,kk) = constcoeff + &
               ddx(1) * coeff(simnum,1)  + &
               ddx(3) * coeff(simnum,2) 

         enddo
         enddo

      enddo
      enddo
      enddo

!     X-Y Faces

      do i = 1,nx-1,refinerat
      do j = 1,ny-1,refinerat
      do k = 1,nz,nz-1
         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

!        *** put coefficient computations here ***

         delx(1) = xx(i+refinerat) - xx(i)
         delx(2) = yy(j+refinerat) - yy(j)
         constcoeff = childgrid(i,j,k)
         coeff(1,1) = (childgrid(i+refinerat,j,k) - constcoeff)/delx(1)
         coeff(1,2) = (childgrid(i+refinerat,j+refinerat,k) - &
                       childgrid(i+refinerat,j,k))/delx(2)
         coeff(2,1) = (childgrid(i+refinerat,j+refinerat,k) - &
                       childgrid(i,j+refinerat,k))/delx(1)
         coeff(2,2) = (childgrid(i,j+refinerat,k) - constcoeff)/delx(2)

         kk = k
         do ii=i,iimax
         do jj=j,jjmax
            ddx(1) = xx(ii) - xx(i)
            ddx(2) = yy(jj) - yy(j)
            if(ddx(2) .le. ddx(1)) then
!              SIMPLEX1 : dy <= dx 
               simnum = 1
            else
!              SIMPLEX2 : dx <= dy 
               simnum = 2
            endif
            childgrid(ii,jj,kk) = constcoeff + &
               ddx(1) * coeff(simnum,1)  + &
               ddx(2) * coeff(simnum,2) 
!write(*,*) 'iMMb: ',ii,jj,kk,' = ',childgrid(ii,jj,kk),'  ',simnum

         enddo
         enddo

      enddo
      enddo
      enddo


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


