! $Header$

#include "cctk.h"
#include "cctk_Functions.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
      
subroutine SBP_SetNormMask (CCTK_ARGUMENTS)

  implicit none

  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_REAL, parameter :: zero = 0.0
  integer, parameter :: wp = kind(zero)
  CCTK_INT :: symtable, n_elements, nchar, pen_sym_handle, np
  CCTK_INT, dimension(6) :: symbnd
  CCTK_REAL, dimension(:), allocatable :: mask_x, mask_y, mask_z

! Note: The first number is twice the value from the paper, since Carpet
! Multiplies by 1/2 at the boundary when doing the sum reduction.
  CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 1.0_wp, 1.0_wp /) 
  CCTK_REAL, dimension(3), parameter :: bmask_3 = (/ 1.0_wp, 1.0_wp, 1.0_wp /) 
  CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 34.0_wp/48.0_wp, &
                                                     59.0_wp/48.0_wp, &
                                                     43.0_wp/48.0_wp, &
                                                     49.0_wp/48.0_wp /)
  CCTK_REAL, dimension(6), parameter :: bmask_6 = (/ 27298.0_wp/43200._wp, &
                                                     12013.0_wp/8640._wp, &
                                                     2711.0_wp/4320.0_wp, &
                                                     5359.0_wp/4320.0_wp, &
                                                     7877.0_wp/8640.0_wp, &
                                                     43801.0_wp/43200.0_wp /)
  CCTK_REAL, dimension(8), parameter :: bmask_8 = (/ 2996278.0_wp/5080320.0_wp,&
                                                     1107307.0_wp/725760.0_wp, &
                                                     20761.0_wp/80640.0_wp, &
                                                     1304999.0_wp/725760.0_wp, &
                                                     299527.0_wp/725760.0_wp, &
                                                     103097.0_wp/80640.0_wp, &
                                                     670091.0_wp/725760.0_wp, &
                                                     5127739.0_wp/5080320.0_wp/)
                                                     
  CCTK_REAL, dimension(8) :: bmask
  CCTK_INT :: ni, nj, nk

  ni = cctk_lsh(1); nj = cctk_lsh(2); nk = cctk_lsh(3)

  symtable = SymmetryTableHandleForGrid ( cctkGH )
  call Util_TableGetIntArray ( n_elements, symtable, 6, &
                                           symbnd, 'symmetry_handle' )
  
  pen_sym_handle = SymmetryHandleOfName ( 'multipatch' )

  if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. any ( symbnd == pen_sym_handle )) ) then
    allocate ( mask_x(ni), mask_y(nj), mask_z(nk) )
    mask_x = 1.0d0; mask_y = 1.0d0; mask_z = 1.0d0

    if ( CCTK_EQUALS(norm_type,'Diagonal') ) then
      select case (order)
      case (2)
        bmask(1:2) = bmask_2
      case (4)
        bmask(1:4) = bmask_4
      case(6)
        bmask(1:6) = bmask_6
      case(8)
        bmask(1:8) = bmask_8
      end select 
      np = order
    else
      select case (order)
      case (4,6)
        bmask(1:3) = bmask_3
        np = 3
      end select
    end if

    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(1) == pen_sym_handle) ) then
      mask_x(1:np) = bmask(1:np)
    end if
    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(2) == pen_sym_handle) ) then
      mask_x(ni:ni-np+1:-1) = bmask(1:np)
    end if
    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(3) == pen_sym_handle) ) then
      mask_y(1:np) = bmask(1:np)
    end if
    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(4) == pen_sym_handle) ) then
      mask_y(nj:nj-np+1:-1) = bmask(1:np)
    end if
    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(5) == pen_sym_handle) ) then
      mask_z(1:np) = bmask(1:np)
    end if
    if ( pen_sym_handle < 0 .or. (pen_sym_handle >= 0 .and. symbnd(6) == pen_sym_handle) ) then
      mask_z(nk:nk-np+1:-1) = bmask(1:np)
    end if

    nmask = outer_prod_3d ( mask_x, mask_y, mask_z )
  else
    nmask = 1.0_wp
  end if

contains

  function outer_prod_3d ( ax, ay, az )
    CCTK_REAL, dimension(:), intent(IN) :: ax
    CCTK_REAL, dimension(:), intent(IN) :: ay
    CCTK_REAL, dimension(:), intent(IN) :: az
    CCTK_REAL, dimension(size(ax),size(ay),size(az)) :: outer_prod_3d
    CCTK_INT :: nx, ny, nz, i, j, k
  
    nx = size(ax); ny = size(ay); nz = size(az)
    do k = 1, nz
      do j = 1, ny
        do i = 1, nx
          outer_prod_3d(i,j,k) = ax(i) * ay(j) * az(k)
        end do
      end do
    end do
  end function outer_prod_3d

end subroutine SBP_SetNormMask