c-----------------------------------------------------------------------
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996

c utility.f 
c Collection of utility routines for TCODE
c-----------------------------------------------------------------------
c  Routine to calculate the inverse 3-metric. Returns the determinant of
c  the metric also. 
c  Routines require mask function for black hole case
c
c  routine verified 
c
      subroutine inverse(g11,g12,g13,g22,g23,g33,
     *   gup11,gup12,gup13,gup22, gup23, gup33,detg,mask,nx, ny, nz)
        implicit none

        integer  nx, ny, nz

        real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx,ny,nz),
     *           gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx,ny,nz),
     *           detg(nx,ny,nz), mask(nx,ny,nz)

c     Local variables

      integer  i, j, k


		real*8   zero, value
		parameter(zero=0.0e0, value=1.0e0)

c   first calculate the determinant

        call calcdet(g11,g12,g13,g22,g23,g33,detg,nx,ny,nz)

c	 putting phony values for zero-determinant...
        
        call replace_hole_val(detg, zero, value, mask, nx,ny,nz)
c

c   construct the inverse component by component

c gup11
        do k = 1,  nz
          do j = 1,  ny
            do i = 1,  nx
            gup11(i,j,k) = g22(i,j,k)*g33(i,j,k)-g23(i,j,k)*g23(i,j,k)
            gup11(i,j,k) = gup11(i,j,k)/detg(i,j,k)
            end do
          end do
        end do

c gup12
        do k = 1,  nz
          do j = 1,  ny
            do i = 1,  nx
              gup12(i,j,k) = g13(i,j,k)*g23(i,j,k)-g12(i,j,k)*g33(i,j,k)
              gup12(i,j,k) = gup12(i,j,k)/detg(i,j,k)
            end do
          end do
        end do

c gup13
        do k = 1,  nz
          do j = 1,  ny
            do i = 1,  nx
              gup13(i,j,k) = g12(i,j,k)*g23(i,j,k)-g13(i,j,k)*g22(i,j,k)
              gup13(i,j,k) = gup13(i,j,k) / detg(i,j,k)
            end do
          end do
        end do

c gup22
        do k = 1,  nz
          do j = 1,  ny
            do i = 1,  nx
              gup22(i,j,k) = g11(i,j,k)*g33(i,j,k)-g13(i,j,k)*g13(i,j,k)
              gup22(i,j,k) = gup22(i,j,k) / detg(i,j,k)
            end do
          end do
        end do

c gup23
        do k = 1,  nz
          do j = 1,  ny
            do i = 1,  nx
              gup23(i,j,k) = g13(i,j,k)*g12(i,j,k)-g11(i,j,k)*g23(i,j,k)
              gup23(i,j,k) = gup23(i,j,k) / detg(i,j,k)
            end do
          end do
        end do

c gup33
        do k = 1,  nz
          do j = 1,  ny
            do i = 1,  nx
            gup33(i,j,k) = g11(i,j,k)*g22(i,j,k)-g12(i,j,k)*g12(i,j,k)
            gup33(i,j,k) = gup33(i,j,k) / detg(i,j,k)
            end do
          end do
        end do

        return
      end 
c-----------------------------------------------------------------------
c  Routine to calculate the determinant of the 3-metric 
c
c  routine verified 

      subroutine calcdet(g11,g12,g13,g22,g23,g33,detg,nx,ny,nz)
        implicit none
       
        integer  nx, ny, nz

        real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *           g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz),
     *           detg(nx,ny,nz)

c   Local variables

        integer i, j, k

        do k = 1, nz
          do j = 1, ny
            do i = 1, nx

             detg(i,j,k) = g11(i,j,k)*g22(i,j,k)*g33(i,j,k) -
     #               g11(i,j,k)*g23(i,j,k)*g23(i,j,k)  -
     #               g12(i,j,k)*g12(i,j,k)*g33(i,j,k)  +
     #               g12(i,j,k)*g13(i,j,k)*g23(i,j,k)  +
     #               g13(i,j,k)*g12(i,j,k)*g23(i,j,k)  -
     #               g13(i,j,k)*g13(i,j,k)*g22(i,j,k)

            end do
          end do
        end do
        
        return
      end
c-----------------------------------------------------------------------
      subroutine setval(nx,ny,nz,f,val)
      implicit none
      integer nx,ny,nz
      real*8 f(nx,ny,nz)
      real*8 val

      integer i,j,k

      do k = 1 , nz
         do j = 1 , ny
            do i = 1 , nx

               f(i,j,k) = val
            end do
         end do
      end do

      return
      end 

c-----------------------------------------------------------------------
      subroutine copy(nx, ny, nz, uin, uout)
        implicit none

        integer  nx, ny, nz

        real*8  uin(nx, ny, nz), uout(nx, ny, nz)

c    Local variables

        integer i, j, k

c    copy over uin to uout

        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              uout(i,j,k) = uin(i,j,k)
            end do
          end do
        end do
        return
      end
c-----------------------------------------------------------------------
      subroutine err(routine,reason)
      implicit none
      character*(*) routine,reason

      write(0,*) 'Error in',routine
      write(0,*) 'routine failed since',reason

      return
      end
C--------------------------------------------------------------------------
C this is not a matrix-matrix multiply but just a pointwise multiplication
      subroutine mm_mult(nx,ny,nz,u1in,u2in,uout)
      implicit none
      integer nx,ny,nz
      real*8 u1in(nx,ny,nz),u2in(nx,ny,nz)
      real*8 uout(nx,ny,nz)

      integer i,j,k

      do k = 1 , nz
         do j = 1 , ny
            do i = 1 , nx
               uout(i,j,k) = u1in(i,j,k) * u2in(i,j,k)
            end do
         end do
      end do

      return
      end
C--------------------------------------------------------------------------
      subroutine mm_add(nx,ny,nz,u1in,u2in,uout)
      implicit none
      integer nx,ny,nz
      real*8 u1in(nx,ny,nz),u2in(nx,ny,nz)
      real*8 uout(nx,ny,nz)

      integer i,j,k

      do k = 1 , nz
         do j = 1 , ny
            do i = 1 , nx
               uout(i,j,k) = u1in(i,j,k) + u2in(i,j,k)
            end do                                          
         end do
      end do

      return
      end
C--------------------------------------------------------------------------
      subroutine mm_sub(nx,ny,nz,u1in,u2in,uout)
      implicit none
      integer nx,ny,nz
      real*8 u1in(nx,ny,nz),u2in(nx,ny,nz)
      real*8 uout(nx,ny,nz)

      integer i,j,k

      do k = 1 , nz
         do j = 1 , ny
            do i = 1 , nx
               uout(i,j,k) = u1in(i,j,k) - u2in(i,j,k)
            end do                                          
         end do
      end do

      return
      end
C-----------------------------------------------------------------------------
C     scalar matrix multiply
      subroutine ms_mult(nx,ny,nz,scin,uin,uout)
      implicit none
      integer nx,ny,nz
      real*8 scin,uin(nx,ny,nz)
      real*8 uout(nx,ny,nz)

      integer i,j,k

      do k = 1 , nz
         do j = 1 , ny
            do i = 1 , nx
               uout(i,j,k) = scin * uin(i,j,k)
            end do                                          
         end do
      end do

      return
      end
C-----------------------------------------------------------------------------
      subroutine dump(nx,ny,nz,x,f,name,j,k,mf)
      implicit none
      integer nx,ny,nz
      real*8 f(nx,ny,nz),x(nx),v(513)
      character*4 name
      integer i,j,k
      integer vsrc,vsxynt
      real*8 mf
C      return

      do i = 1 , nx
         v(i) = abs(f(i,j,k)*mf)
      end do

      vsrc = vsxynt(name,0.0e0,x,v,nx)

      return
      end
c-----------------------------------------------------------------------
      real*8 function n2i(n,a)
      implicit none
      integer n
      real*8 a(n,n,n)
      integer i,j,k

      n2i = 0.0e0

      do k = 3 , n - 2
         do j = 3 , n - 2
            do i = 3 , n - 2
               n2i = n2i + a(i,j,k)**2
            end do
         end do
      end do

      n2i = sqrt( n2i / (n-4)**3 )

      return
      end
C-------------------------------------------------------------------------
      real*8 function n2ib(nx,ny,nz,a,m,nl)
      implicit none
      integer nx,ny,nz
      real*8 a(nx,ny,nz)
      integer nl
      integer m(3*nl)
		integer l,i,j,k

		n2ib = 0.0e0
      do l = 1 , nl

         i  = m(l)
         j  = m(nl+l)
         k  = m(2*nl+l)

			n2ib = n2ib + a(i,j,k)*a(i,j,k)

		end do

		n2ib = sqrt( n2ib / nl )

		return
		end
C-------------------------------------------------------------------------
      real*8 function n2b(n,a)
      implicit none
      integer n
      real*8 a(n,n,n)
      real*8 count
      integer i,j,k

      n2b = 0.0e0
      count = 0.0e0

      do k = 1 , n
         do j = 1 , n
            n2b = n2b + a(1,j,k)**2
C            n2b = n2b + a(2,j,k)**2
C            n2b = n2b + a(n,j,k)**2
C            n2b = n2b + a(n-1,j,k)**2
            count = count + 1.0e0
         end do
      end do

      do k = 1 , n
         do i = 1 , n
            n2b = n2b + a(i,1,k)**2
C            n2b = n2b + a(i,2,k)**2
C            n2b = n2b + a(i,n,k)**2
C            n2b = n2b + a(i,n-1,k)**2
            count = count + 1.0e0
         end do
      end do

      do j = 1 , n
         do i = 1 , n
            n2b = n2b + a(i,j,1)**2
C            n2b = n2b + a(i,j,2)**2
C            n2b = n2b + a(i,j,n)**2
C            n2b = n2b + a(i,j,n-1)**2
            count = count + 1.0e0
         end do
      end do


      n2b = sqrt(n2b/count)

      return
      end
c-----------------------------------------------------------------------
c  routine to raise first index of (0 2) tensor. notation->(contra cov)index
c
      subroutine rais02i1(nx, ny, nz, 
     *           gup11, gup12, gup13, gup22, gup23, gup33,
     *           t11, t12, t13, t21, t22, t23, t31, t32, t33, 
     *           temp1, temp2, temp3)
        implicit none

        integer  nx, ny, nz

        real*8  gup11(nx,ny,nz), gup12(nx,ny,nz), gup13(nx, ny,nz)
        real*8  gup22(nx,ny,nz), gup23(nx,ny,nz), gup33(nx, ny,nz)

        real*8  t11(nx,ny,nz), t12(nx,ny,nz), t13(nx, ny,nz)
        real*8  t21(nx,ny,nz), t22(nx,ny,nz), t23(nx, ny,nz)
        real*8  t31(nx,ny,nz), t32(nx,ny,nz), t33(nx, ny,nz)
        real*8  temp1(nx,ny,nz), temp2(nx,ny,nz), temp3(nx, ny,nz)

c      Local variables

        integer i, j, k
c {i}1
        call copy(nx, ny, nz, t11, temp1)
        call copy(nx, ny, nz, t21, temp2)
        call copy(nx, ny, nz, t31, temp3)

c     ^1_1
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t11(i,j,k) = gup11(i,j,k) * temp1(i,j,k) 
     *                   + gup12(i,j,k) * temp2(i,j,k)
     *                   + gup13(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c     ^2_1
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t21(i,j,k) = gup12(i,j,k) * temp1(i,j,k) 
     *                   + gup22(i,j,k) * temp2(i,j,k)
     *                   + gup23(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c     ^3_1
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t31(i,j,k) = gup13(i,j,k) * temp1(i,j,k) 
     *                   + gup23(i,j,k) * temp2(i,j,k)
     *                   + gup33(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c {i}2
        call copy(nx, ny, nz, t12, temp1)
        call copy(nx, ny, nz, t22, temp2)
        call copy(nx, ny, nz, t32, temp3)

c     ^1_2
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t12(i,j,k) = gup11(i,j,k) * temp1(i,j,k) 
     *                   + gup12(i,j,k) * temp2(i,j,k)
     *                   + gup13(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c     ^2_2
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t22(i,j,k) = gup12(i,j,k) * temp1(i,j,k) 
     *                   + gup22(i,j,k) * temp2(i,j,k)
     *                   + gup23(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c     ^3_2
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t32(i,j,k) = gup13(i,j,k) * temp1(i,j,k) 
     *                   + gup23(i,j,k) * temp2(i,j,k)
     *                   + gup33(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c {i}3
        call copy(nx, ny, nz, t13, temp1)
        call copy(nx, ny, nz, t23, temp2)
        call copy(nx, ny, nz, t33, temp3)

c     ^1_3
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t13(i,j,k) = gup11(i,j,k) * temp1(i,j,k) 
     *                   + gup12(i,j,k) * temp2(i,j,k)
     *                   + gup13(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c     ^2_3
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t23(i,j,k) = gup12(i,j,k) * temp1(i,j,k) 
     *                   + gup22(i,j,k) * temp2(i,j,k)
     *                   + gup23(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

c     ^3_3
        do k = 1, nz
          do j = 1, ny
            do i = 1, nx
              t33(i,j,k) = gup13(i,j,k) * temp1(i,j,k) 
     *                   + gup23(i,j,k) * temp2(i,j,k)
     *                   + gup33(i,j,k) * temp3(i,j,k)
              
            end do
          end do
        end do

        return
      end
c-----------------------------------------------------------------------
		subroutine replace_hole_val(array, chkval, setval, mask,
     *           nx,ny,nz)
		implicit none

		integer  nx, ny, nz

		real*8   array(nx,ny,nz), chkval, setval, mask(nx,ny,nz)

c		Local variables 
		integer	i,j,k

      do i = 1, nx
      	do j = 1, ny
      		do k = 1, nz

               array(i,j,k) = (mask(i,j,k)-chkval)*array(i,j,k)
     &                      + (1.0e0-(mask(i,j,k)-chkval))*setval

				end do
			end do
		end do

		return
		end 
************************************************************************
*  Forward smoother for inner boundary

		subroutine fwd_smooth1(nx, ny, nz, nl, m, d, mask, 
     *           g11, g12, g13, g22, g23, g33,
     *           k11, k12, k13, k21, k22, k23, k31, k32, k33,
     *           coeff)
			implicit none

			integer	nx, ny, nz, nl

			integer	m(3*nl), d(3*nl)

			real*8	mask(nx, ny, nz)

			real*8	g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *            g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz)

			real*8	k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     *            k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     *            k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)

			real*8	coeff

c		Local variables

			integer	i, j, k, l, p

			real*8	sum, c2

c		Loop over points on boundary

			do l = 1, nl

				i = m(l)
				j = m(l + nl )
				k = m(l + 2 * nl )

c			sum over defined points to get proper fraction of coeff
				sum = mask(i+1, j  , k  ) + mask(i-1, j  , k  ) +
     *            mask(i  , j+1, k  ) + mask(i  , j-1, k  ) +
     *            mask(i  , j  , k+1) + mask(i  , j  , k-1)

c			Now form other coefficient so that all add to one
				c2 = (1.0 - coeff) / sum

c			Now average

				g11(i,j,k) = coeff * g11(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g11(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g11(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g11(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g11(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g11(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g11(i,j,k-1) 

				g12(i,j,k) = coeff * g12(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g12(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g12(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g12(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g12(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g12(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g12(i,j,k-1) 

				g13(i,j,k) = coeff * g13(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g13(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g13(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g13(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g13(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g13(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g13(i,j,k-1) 

				g22(i,j,k) = coeff * g22(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g22(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g22(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g22(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g22(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g22(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g22(i,j,k-1) 

				g23(i,j,k) = coeff * g23(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g23(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g23(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g23(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g23(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g23(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g23(i,j,k-1) 

				g33(i,j,k) = coeff * g33(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g33(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g33(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g33(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g33(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g33(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g33(i,j,k-1) 


				k11(i,j,k) = coeff * k11(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k11(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k11(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k11(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k11(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k11(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k11(i,j,k-1) 

				k12(i,j,k) = coeff * k12(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k12(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k12(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k12(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k12(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k12(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k12(i,j,k-1) 

				k13(i,j,k) = coeff * k13(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k13(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k13(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k13(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k13(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k13(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k13(i,j,k-1) 

				k21(i,j,k) = coeff * k21(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k21(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k21(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k21(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k21(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k21(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k21(i,j,k-1) 

				k22(i,j,k) = coeff * k22(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k22(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k22(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k22(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k22(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k22(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k22(i,j,k-1) 

				k23(i,j,k) = coeff * k23(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k23(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k23(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k23(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k23(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k23(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k23(i,j,k-1) 

				k31(i,j,k) = coeff * k31(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k31(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k31(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k31(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k31(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k31(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k31(i,j,k-1) 

				k32(i,j,k) = coeff * k32(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k32(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k32(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k32(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k32(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k32(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k32(i,j,k-1) 

				k33(i,j,k) = coeff * k33(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k33(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k33(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k33(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k33(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k33(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k33(i,j,k-1) 
			end do

			return
		end
************************************************************************
*  Backward smoother for inner boundary

		subroutine back_smooth1(nx, ny, nz, nl, m, d, mask, 
     *           g11, g12, g13, g22, g23, g33,
     *           k11, k12, k13, k21, k22, k23, k31, k32, k33,
     *           coeff)
			implicit none

			integer	nx, ny, nz, nl

			integer	m(3*nl), d(3*nl)

			real*8	mask(nx, ny, nz)

			real*8	g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     *            g22(nx,ny,nz), g23(nx,ny,nz), g33(nx,ny,nz)

			real*8	k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     *            k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     *            k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)

			real*8	coeff

c		Local variables

			integer	i, j, k, l, p

			real*8	sum, c2

c		Loop over points on boundary

			do l = nl, 1, -1

				i = m(l)
				j = m(l + nl )
				k = m(l + 2 * nl )

c			sum over defined points to get proper fraction of coeff
				sum = mask(i+1, j  , k  ) + mask(i-1, j  , k  ) +
     *            mask(i  , j+1, k  ) + mask(i  , j-1, k  ) +
     *            mask(i  , j  , k+1) + mask(i  , j  , k-1)

c			Now form other coefficient so that all add to one
				c2 = (1.0 - coeff) / sum

c			Now average

				g11(i,j,k) = coeff * g11(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g11(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g11(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g11(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g11(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g11(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g11(i,j,k-1) 

				g12(i,j,k) = coeff * g12(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g12(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g12(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g12(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g12(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g12(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g12(i,j,k-1) 

				g13(i,j,k) = coeff * g13(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g13(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g13(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g13(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g13(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g13(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g13(i,j,k-1) 

				g22(i,j,k) = coeff * g22(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g22(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g22(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g22(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g22(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g22(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g22(i,j,k-1) 

				g23(i,j,k) = coeff * g23(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g23(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g23(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g23(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g23(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g23(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g23(i,j,k-1) 

				g33(i,j,k) = coeff * g33(i,j,k) + 
     *          c2 * mask(i+1,j,k) * g33(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * g33(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * g33(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * g33(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * g33(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * g33(i,j,k-1) 


				k11(i,j,k) = coeff * k11(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k11(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k11(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k11(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k11(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k11(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k11(i,j,k-1) 

				k12(i,j,k) = coeff * k12(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k12(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k12(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k12(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k12(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k12(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k12(i,j,k-1) 

				k13(i,j,k) = coeff * k13(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k13(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k13(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k13(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k13(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k13(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k13(i,j,k-1) 

				k21(i,j,k) = coeff * k21(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k21(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k21(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k21(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k21(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k21(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k21(i,j,k-1) 

				k22(i,j,k) = coeff * k22(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k22(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k22(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k22(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k22(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k22(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k22(i,j,k-1) 

				k23(i,j,k) = coeff * k23(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k23(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k23(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k23(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k23(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k23(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k23(i,j,k-1) 

				k31(i,j,k) = coeff * k31(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k31(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k31(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k31(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k31(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k31(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k31(i,j,k-1) 

				k32(i,j,k) = coeff * k32(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k32(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k32(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k32(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k32(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k32(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k32(i,j,k-1) 

				k33(i,j,k) = coeff * k33(i,j,k) + 
     *          c2 * mask(i+1,j,k) * k33(i+1,j,k) + 
     *          c2 * mask(i-1,j,k) * k33(i-1,j,k) + 
     *          c2 * mask(i,j+1,k) * k33(i,j+1,k) + 
     *          c2 * mask(i,j-1,k) * k33(i,j-1,k) + 
     *          c2 * mask(i,j,k+1) * k33(i,j,k+1) + 
     *          c2 * mask(i,j,k-1) * k33(i,j,k-1) 
			end do

			return
		end
************************************************************************