From c92386fd0eab5c2697ef3ab8c7c5531e0c1d42b3 Mon Sep 17 00:00:00 2001
From: diener <diener@f69c4107-0314-4c4f-9ad4-17e986b73f4a>
Date: Thu, 2 Feb 2006 15:12:27 +0000
Subject: Kreiss-Oliger type dissipation operators compatible with the diagonal
 4-2 and 6-3 derivative operators.

git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@55 f69c4107-0314-4c4f-9ad4-17e986b73f4a
---
 src/Dissipation_4_2_alt.F90 | 240 +++++++++++++++++++++++++++++
 src/Dissipation_6_3_alt.F90 | 359 ++++++++++++++++++++++++++++++++++++++++++++
 src/dissipation.c           |  20 ++-
 src/make.code.defn          |   2 +
 4 files changed, 617 insertions(+), 4 deletions(-)
 create mode 100644 src/Dissipation_4_2_alt.F90
 create mode 100644 src/Dissipation_6_3_alt.F90

diff --git a/src/Dissipation_4_2_alt.F90 b/src/Dissipation_4_2_alt.F90
new file mode 100644
index 0000000..c0a68ba
--- /dev/null
+++ b/src/Dissipation_4_2_alt.F90
@@ -0,0 +1,240 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs)
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+
+  integer ::  ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
+  CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
+  CCTK_INT, dimension(6), intent(in) :: bb
+  CCTK_INT, dimension(3), intent(in) :: gsize
+  CCTK_REAL, dimension(3), intent(in) :: delta
+  CCTK_REAL, intent(in) :: epsilon 
+
+  CCTK_REAL :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  integer :: i, j, k
+  CCTK_REAL, dimension(7,4) :: a
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 64 * delta(1) )
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize(1)
+  else
+    rhs(1,:,:) = rhs(1,:,:) + &
+                 ( a(1,1) * var(1,:,:) + a(2,1) * var(2,:,:) + &
+                   a(3,1) * var(3,:,:) + a(4,1) * var(4,:,:) ) * idel
+    rhs(2,:,:) = rhs(2,:,:) + &
+                 ( a(1,2) * var(1,:,:) + a(2,2) * var(2,:,:) + &
+                   a(3,2) * var(3,:,:) + a(4,2) * var(4,:,:) + &
+                   a(5,2) * var(5,:,:) ) * idel
+    rhs(3,:,:) = rhs(3,:,:) + &
+                 ( a(1,3) * var(1,:,:) + a(2,3) * var(2,:,:) + &
+                   a(3,3) * var(3,:,:) + a(4,3) * var(4,:,:) + &
+                   a(5,3) * var(5,:,:) + a(6,3) * var(6,:,:) ) * idel
+    rhs(4,:,:) = rhs(4,:,:) + &
+                 ( a(1,4) * var(1,:,:) + a(2,4) * var(2,:,:) + &
+                   a(3,4) * var(3,:,:) + a(4,4) * var(4,:,:) + &
+                   a(5,4) * var(5,:,:) + a(6,4) * var(6,:,:) + &
+                   a(7,4) * var(7,:,:) ) * idel
+
+    il = 5
+  end if
+  if ( bb(2) == 0 ) then
+    ir = ni - gsize(1)
+  else
+    rhs(ni-3,:,:) = rhs(ni-3,:,:) + &
+                 ( a(1,4) * var(ni,:,:) + a(2,4) * var(ni-1,:,:) + &
+                   a(3,4) * var(ni-2,:,:) + a(4,4) * var(ni-3,:,:) + &
+                   a(5,4) * var(ni-4,:,:) + a(6,4) * var(ni-5,:,:) + &
+                   a(7,4) * var(ni-6,:,:) ) * idel
+    rhs(ni-2,:,:) = rhs(ni-2,:,:) + &
+                 ( a(1,3) * var(ni,:,:) + a(2,3) * var(ni-1,:,:) + &
+                   a(3,3) * var(ni-2,:,:) + a(4,3) * var(ni-3,:,:) + &
+                   a(5,3) * var(ni-4,:,:) + a(6,3) * var(ni-5,:,:) ) * idel
+    rhs(ni-1,:,:) = rhs(ni-1,:,:) + &
+                 ( a(1,2) * var(ni,:,:) + a(2,2) * var(ni-1,:,:) + &
+                   a(3,2) * var(ni-2,:,:) + a(4,2) * var(ni-3,:,:) + &
+                   a(5,2) * var(ni-4,:,:) ) * idel
+    rhs(ni,:,:)   = rhs(ni,:,:) + &
+                 ( a(1,1) * var(ni,:,:) + a(2,1) * var(ni-1,:,:) + &
+                   a(3,1) * var(ni-2,:,:) + a(4,1) * var(ni-3,:,:) ) * idel
+
+    ir = ni - 4
+  end if
+  rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+                   ( -20.0_wp * var(il:ir,:,:) + &
+                      15.0_wp * ( var(il-1:ir-1,:,:) + &
+                                  var(il+1:ir+1,:,:) ) - &
+                       6.0_wp * ( var(il-2:ir-2,:,:) + &
+                                  var(il+2:ir+2,:,:) ) + &
+                                ( var(il-3:ir-3,:,:) + &
+                                  var(il+3:ir+3,:,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 64 * delta(2) )
+
+  if ( bb(3) == 0 ) then
+    jl = 1 + gsize(2)
+  else
+    rhs(:,1,:) = rhs(:,1,:) + &
+                 ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) + &
+                   a(3,1) * var(:,3,:) + a(4,1) * var(:,4,:) ) * idel
+    rhs(:,2,:) = rhs(:,2,:) + &
+                 ( a(1,2) * var(:,1,:) + a(2,2) * var(:,2,:) + &
+                   a(3,2) * var(:,3,:) + a(4,2) * var(:,4,:) + &
+                   a(5,2) * var(:,5,:) ) * idel
+    rhs(:,3,:) = rhs(:,3,:) + &
+                 ( a(1,3) * var(:,1,:) + a(2,3) * var(:,2,:) + &
+                   a(3,3) * var(:,3,:) + a(4,3) * var(:,4,:) + &
+                   a(5,3) * var(:,5,:) + a(6,3) * var(:,6,:) ) * idel
+    rhs(:,4,:) = rhs(:,4,:) + &
+                 ( a(1,4) * var(:,1,:) + a(2,4) * var(:,2,:) + &
+                   a(3,4) * var(:,3,:) + a(4,4) * var(:,4,:) + &
+                   a(5,4) * var(:,5,:) + a(6,4) * var(:,6,:) + &
+                   a(7,4) * var(:,7,:) ) * idel
+
+    jl = 5
+  end if
+  if ( bb(4) == 0 ) then
+    jr = nj - gsize(2)
+  else
+    rhs(:,nj-3,:) = rhs(:,nj-3,:) + &
+                 ( a(1,4) * var(:,nj,:) + a(2,4) * var(:,nj-1,:) + &
+                   a(3,4) * var(:,nj-2,:) + a(4,4) * var(:,nj-3,:) + &
+                   a(5,4) * var(:,nj-4,:) + a(6,4) * var(:,nj-5,:) + &
+                   a(7,4) * var(:,nj-6,:) ) * idel
+    rhs(:,nj-2,:) = rhs(:,nj-2,:) + &
+                 ( a(1,3) * var(:,nj,:) + a(2,3) * var(:,nj-1,:) + &
+                   a(3,3) * var(:,nj-2,:) + a(4,3) * var(:,nj-3,:) + &
+                   a(5,3) * var(:,nj-4,:) + a(6,3) * var(:,nj-5,:) ) * idel
+    rhs(:,nj-1,:) = rhs(:,nj-1,:) + &
+                 ( a(1,2) * var(:,nj,:) + a(2,2) * var(:,nj-1,:) + &
+                   a(3,2) * var(:,nj-2,:) + a(4,2) * var(:,nj-3,:) + &
+                   a(5,2) * var(:,nj-4,:) ) * idel
+    rhs(:,nj,:)   = rhs(:,nj,:) + &
+                 ( a(1,1) * var(:,nj,:) + a(2,1) * var(:,nj-1,:) + &
+                   a(3,1) * var(:,nj-2,:) + a(4,1) * var(:,nj-3,:) ) * idel
+
+    jr = nj - 4
+  end if
+  rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+                   ( -20.0_wp * var(:,jl:jr,:) + &
+                      15.0_wp * ( var(:,jl-1:jr-1,:) + &
+                                  var(:,jl+1:jr+1,:) ) - &
+                       6.0_wp * ( var(:,jl-2:jr-2,:) + &
+                                  var(:,jl+2:jr+2,:) ) + &
+                                ( var(:,jl-3:jr-3,:) + &
+                                  var(:,jl+3:jr+3,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 64 * delta(3) )
+
+  if ( bb(5) == 0 ) then
+    kl = 1 + gsize(3)
+  else
+    rhs(:,:,1) = rhs(:,:,1) + &
+                 ( a(1,1) * var(:,:,1) + a(2,1) * var(:,:,2) + &
+                   a(3,1) * var(:,:,3) + a(4,1) * var(:,:,4) ) * idel
+    rhs(:,:,2) = rhs(:,:,2) + &
+                 ( a(1,2) * var(:,:,1) + a(2,2) * var(:,:,2) + &
+                   a(3,2) * var(:,:,3) + a(4,2) * var(:,:,4) + &
+                   a(5,2) * var(:,:,5) ) * idel
+    rhs(:,:,3) = rhs(:,:,3) + &
+                 ( a(1,3) * var(:,:,1) + a(2,3) * var(:,:,2) + &
+                   a(3,3) * var(:,:,3) + a(4,3) * var(:,:,4) + &
+                   a(5,3) * var(:,:,5) + a(6,3) * var(:,:,6) ) * idel
+    rhs(:,:,4) = rhs(:,:,4) + &
+                 ( a(1,4) * var(:,:,1) + a(2,4) * var(:,:,2) + &
+                   a(3,4) * var(:,:,3) + a(4,4) * var(:,:,4) + &
+                   a(5,4) * var(:,:,5) + a(6,4) * var(:,:,6) + &
+                   a(7,4) * var(:,:,7) ) * idel
+
+    kl = 5
+  end if
+  if ( bb(6) == 0 ) then
+    kr = nk - gsize(3)
+  else
+    rhs(:,:,nk-3) = rhs(:,:,nk-3) + &
+                 ( a(1,4) * var(:,:,nk) + a(2,4) * var(:,:,nk-1) + &
+                   a(3,4) * var(:,:,nk-2) + a(4,4) * var(:,:,nk-3) + &
+                   a(5,4) * var(:,:,nk-4) + a(6,4) * var(:,:,nk-5) + &
+                   a(7,4) * var(:,:,nk-6) ) * idel
+    rhs(:,:,nk-2) = rhs(:,:,nk-2) + &
+                 ( a(1,3) * var(:,:,nk) + a(2,3) * var(:,:,nk-1) + &
+                   a(3,3) * var(:,:,nk-2) + a(4,3) * var(:,:,nk-3) + &
+                   a(5,3) * var(:,:,nk-4) + a(6,3) * var(:,:,nk-5) ) * idel
+    rhs(:,:,nk-1) = rhs(:,:,nk-1) + &
+                 ( a(1,2) * var(:,:,nk) + a(2,2) * var(:,:,nk-1) + &
+                   a(3,2) * var(:,:,nk-2) + a(4,2) * var(:,:,nk-3) + &
+                   a(5,2) * var(:,:,nk-4) ) * idel
+    rhs(:,:,nk)   = rhs(:,:,nk) + &
+                 ( a(1,1) * var(:,:,nk) + a(2,1) * var(:,:,nk-1) + &
+                   a(3,1) * var(:,:,nk-2) + a(4,1) * var(:,:,nk-3) ) * idel
+
+    kr = nk - 4
+  end if
+  rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+                   ( -20.0_wp * var(:,:,kl:kr) + &
+                      15.0_wp * ( var(:,:,kl-1:kr-1) + &
+                                  var(:,:,kl+1:kr+1) ) - &
+                       6.0_wp * ( var(:,:,kl-2:kr-2) + &
+                                  var(:,:,kl+2:kr+2) ) + &
+                                ( var(:,:,kl-3:kr-3) + &
+                                  var(:,:,kl+3:kr+3) ) ) * idel
+
+contains
+
+  subroutine set_coeff ( a )
+
+    implicit none
+
+    CCTK_REAL, dimension(7,4), intent(OUT) :: a
+    CCTK_REAL :: zero = 0.0
+    integer, parameter :: wp = kind(zero)
+
+    a(1,1) = -2.8235294117647058823529411764705882352941176470588_wp
+    a(2,1) = 8.4705882352941176470588235294117647058823529411765_wp
+    a(3,1) = -8.4705882352941176470588235294117647058823529411765_wp
+    a(4,1) = 2.8235294117647058823529411764705882352941176470588_wp
+    a(5,1) = 0.0_wp
+    a(6,1) = 0.0_wp
+    a(7,1) = 0.0_wp
+    a(1,2) = 2.4406779661016949152542372881355932203389830508475_wp
+    a(2,2) = -8.1355932203389830508474576271186440677966101694915_wp
+    a(3,2) = 9.7627118644067796610169491525423728813559322033898_wp
+    a(4,2) = -4.8813559322033898305084745762711864406779661016949_wp
+    a(5,2) = 0.81355932203389830508474576271186440677966101694915_wp
+    a(6,2) = 0.0_wp
+    a(7,2) = 0.0_wp
+    a(1,3) = -3.3488372093023255813953488372093023255813953488372_wp
+    a(2,3) = 13.395348837209302325581395348837209302325581395349_wp
+    a(3,3) = -21.209302325581395348837209302325581395348837209302_wp
+    a(4,3) = 16.744186046511627906976744186046511627906976744186_wp
+    a(5,3) = -6.6976744186046511627906976744186046511627906976744_wp
+    a(6,3) = 1.1162790697674418604651162790697674418604651162791_wp
+    a(7,3) = 0.0_wp
+    a(1,4) = 0.97959183673469387755102040816326530612244897959184_wp
+    a(2,4) = -5.8775510204081632653061224489795918367346938775510_wp
+    a(3,4) = 14.693877551020408163265306122448979591836734693878_wp
+    a(4,4) = -19.591836734693877551020408163265306122448979591837_wp
+    a(5,4) = 14.693877551020408163265306122448979591836734693878_wp
+    a(6,4) = -5.8775510204081632653061224489795918367346938775510_wp
+    a(7,4) = 0.97959183673469387755102040816326530612244897959184_wp
+
+  end subroutine set_coeff
+
+end subroutine dissipation_4_2_alt
diff --git a/src/Dissipation_6_3_alt.F90 b/src/Dissipation_6_3_alt.F90
new file mode 100644
index 0000000..4a19c11
--- /dev/null
+++ b/src/Dissipation_6_3_alt.F90
@@ -0,0 +1,359 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, delta, epsilon, rhs)
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+
+  integer ::  ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
+  CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
+  CCTK_INT, dimension(6), intent(in) :: bb
+  CCTK_INT, dimension(3), intent(in) :: gsize
+  CCTK_REAL, dimension(3), intent(in) :: delta
+  CCTK_REAL, intent(in) :: epsilon 
+
+  CCTK_REAL :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  CCTK_REAL, dimension(10,6) :: a
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 256 * delta(1) )
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize(1)
+  else
+    rhs(1,:,:) = rhs(1,:,:) + &
+                 ( a(1,1) * var(1,:,:) + a(2,1) * var(2,:,:) + &
+                   a(3,1) * var(3,:,:) + a(4,1) * var(4,:,:) + &
+                   a(5,1) * var(5,:,:) ) * idel
+    rhs(2,:,:) = rhs(2,:,:) + &
+                 ( a(1,2) * var(1,:,:) + a(2,2) * var(2,:,:) + &
+                   a(3,2) * var(3,:,:) + a(4,2) * var(4,:,:) + &
+                   a(5,2) * var(5,:,:) + a(6,2) * var(6,:,:) ) * idel
+    rhs(3,:,:) = rhs(3,:,:) + &
+                 ( a(1,3) * var(1,:,:) + a(2,3) * var(2,:,:) + &
+                   a(3,3) * var(3,:,:) + a(4,3) * var(4,:,:) + &
+                   a(5,3) * var(5,:,:) + a(6,3) * var(6,:,:) + &
+                   a(7,3) * var(7,:,:) ) * idel
+    rhs(4,:,:) = rhs(4,:,:) + &
+                 ( a(1,4) * var(1,:,:) + a(2,4) * var(2,:,:) + &
+                   a(3,4) * var(3,:,:) + a(4,4) * var(4,:,:) + &
+                   a(5,4) * var(5,:,:) + a(6,4) * var(6,:,:) + &
+                   a(7,4) * var(7,:,:) + a(8,4) * var(8,:,:) ) * idel
+    rhs(5,:,:) = rhs(5,:,:) + &
+                 ( a(1,5) * var(1,:,:) + a(2,5) * var(2,:,:) + &
+                   a(3,5) * var(3,:,:) + a(4,5) * var(4,:,:) + &
+                   a(5,5) * var(5,:,:) + a(6,5) * var(6,:,:) + &
+                   a(7,5) * var(7,:,:) + a(8,5) * var(8,:,:) + &
+                   a(9,5) * var(9,:,:) ) * idel
+    rhs(6,:,:) = rhs(6,:,:) + &
+                 ( a(2,6) * var(2,:,:) + a(3,6) * var(3,:,:) + &
+                   a(4,6) * var(4,:,:) + a(5,6) * var(5,:,:) + &
+                   a(6,6) * var(6,:,:) + a(7,6) * var(7,:,:) + &
+                   a(8,6) * var(8,:,:) + a(9,6) * var(9,:,:) + &
+                   a(10,6) * var(10,:,:) ) * idel
+
+    il = 7
+  end if
+  if ( bb(2) == 0 ) then
+    ir = ni - gsize(1)
+  else
+    rhs(ni-5,:,:) = rhs(ni-5,:,:) + &
+                 ( a(2,6) * var(ni-1,:,:) + a(3,6) * var(ni-2,:,:) + &
+                   a(4,6) * var(ni-3,:,:) + a(5,6) * var(ni-4,:,:) + &
+                   a(6,6) * var(ni-5,:,:) + a(7,6) * var(ni-6,:,:) + &
+                   a(8,6) * var(ni-7,:,:) + a(9,6) * var(ni-8,:,:) + &
+                   a(10,6) * var(ni-9,:,:) ) * idel
+    rhs(ni-4,:,:) = rhs(ni-4,:,:) + &
+                 ( a(1,5) * var(ni,:,:) + a(2,5) * var(ni-1,:,:) + &
+                   a(3,5) * var(ni-2,:,:) + a(4,5) * var(ni-3,:,:) + &
+                   a(5,5) * var(ni-4,:,:) + a(6,5) * var(ni-5,:,:) + &
+                   a(7,5) * var(ni-6,:,:) + a(8,5) * var(ni-7,:,:) + &
+                   a(9,5) * var(ni-8,:,:) ) * idel
+    rhs(ni-3,:,:) = rhs(ni-3,:,:) + &
+                 ( a(1,4) * var(ni,:,:) + a(2,4) * var(ni-1,:,:) + &
+                   a(3,4) * var(ni-2,:,:) + a(4,4) * var(ni-3,:,:) + &
+                   a(5,4) * var(ni-4,:,:) + a(6,4) * var(ni-5,:,:) + &
+                   a(7,4) * var(ni-6,:,:) + a(8,4) * var(ni-7,:,:) ) * idel
+    rhs(ni-2,:,:) = rhs(ni-2,:,:) + &
+                 ( a(1,3) * var(ni,:,:) + a(2,3) * var(ni-1,:,:) + &
+                   a(3,3) * var(ni-2,:,:) + a(4,3) * var(ni-3,:,:) + &
+                   a(5,3) * var(ni-4,:,:) + a(6,3) * var(ni-5,:,:) + &
+                   a(7,3) * var(ni-6,:,:) ) * idel
+    rhs(ni-1,:,:) = rhs(ni-1,:,:) + &
+                 ( a(1,2) * var(ni,:,:) + a(2,2) * var(ni-1,:,:) + &
+                   a(3,2) * var(ni-2,:,:) + a(4,2) * var(ni-3,:,:) + &
+                   a(5,2) * var(ni-4,:,:) + a(6,2) * var(ni-5,:,:) ) * idel
+    rhs(ni,:,:)   = rhs(ni,:,:) + &
+                 ( a(1,1) * var(ni,:,:) + a(2,1) * var(ni-1,:,:) + &
+                   a(3,1) * var(ni-2,:,:) + a(4,1) * var(ni-3,:,:) + &
+                   a(5,1) * var(ni-4,:,:) ) * idel
+
+    ir = ni - 6
+  end if
+  rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+                   ( -70.0_wp * var(il:ir,:,:) + &
+                      56.0_wp * ( var(il-1:ir-1,:,:) + &
+                                  var(il+1:ir+1,:,:) ) - &
+                      28.0_wp * ( var(il-2:ir-2,:,:) + &
+                                  var(il+2:ir+2,:,:) ) + &
+                       8.0_wp * ( var(il-3:ir-3,:,:) + &
+                                  var(il+3:ir+3,:,:) ) - &
+                                ( var(il-4:ir-4,:,:) + &
+                                  var(il+4:ir+4,:,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 256 * delta(2) )
+
+  if ( bb(3) == 0 ) then
+    jl = 1 + gsize(2)
+  else
+    rhs(:,1,:) = rhs(:,1,:) + &
+                 ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) + &
+                   a(3,1) * var(:,3,:) + a(4,1) * var(:,4,:) + &
+                   a(5,1) * var(:,5,:) ) * idel
+    rhs(:,2,:) = rhs(:,2,:) + &
+                 ( a(1,2) * var(:,1,:) + a(2,2) * var(:,2,:) + &
+                   a(3,2) * var(:,3,:) + a(4,2) * var(:,4,:) + &
+                   a(5,2) * var(:,5,:) + a(6,2) * var(:,6,:) ) * idel
+    rhs(:,3,:) = rhs(:,3,:) + &
+                 ( a(1,3) * var(:,1,:) + a(2,3) * var(:,2,:) + &
+                   a(3,3) * var(:,3,:) + a(4,3) * var(:,4,:) + &
+                   a(5,3) * var(:,5,:) + a(6,3) * var(:,6,:) + &
+                   a(7,3) * var(:,7,:) ) * idel
+    rhs(:,4,:) = rhs(:,4,:) + &
+                 ( a(1,4) * var(:,1,:) + a(2,4) * var(:,2,:) + &
+                   a(3,4) * var(:,3,:) + a(4,4) * var(:,4,:) + &
+                   a(5,4) * var(:,5,:) + a(6,4) * var(:,6,:) + &
+                   a(7,4) * var(:,7,:) + a(8,4) * var(:,8,:) ) * idel
+    rhs(:,5,:) = rhs(:,5,:) + &
+                 ( a(1,5) * var(:,1,:) + a(2,5) * var(:,2,:) + &
+                   a(3,5) * var(:,3,:) + a(4,5) * var(:,4,:) + &
+                   a(5,5) * var(:,5,:) + a(6,5) * var(:,6,:) + &
+                   a(7,5) * var(:,7,:) + a(8,5) * var(:,8,:) + &
+                   a(9,5) * var(:,9,:) ) * idel
+    rhs(:,6,:) = rhs(:,6,:) + &
+                 ( a(2,6) * var(:,2,:) + a(3,6) * var(:,3,:) + &
+                   a(4,6) * var(:,4,:) + a(5,6) * var(:,5,:) + &
+                   a(6,6) * var(:,6,:) + a(7,6) * var(:,7,:) + &
+                   a(8,6) * var(:,8,:) + a(9,6) * var(:,9,:) + &
+                   a(10,6) * var(:,10,:) ) * idel
+
+    jl = 7
+  end if
+  if ( bb(4) == 0 ) then
+    jr = nj - gsize(2)
+  else
+    rhs(:,nj-5,:) = rhs(:,nj-5,:) + &
+                 ( a(2,6) * var(:,nj-1,:) + a(3,6) * var(:,nj-2,:) + &
+                   a(4,6) * var(:,nj-3,:) + a(5,6) * var(:,nj-4,:) + &
+                   a(6,6) * var(:,nj-5,:) + a(7,6) * var(:,nj-6,:) + &
+                   a(8,6) * var(:,nj-7,:) + a(9,6) * var(:,nj-8,:) + &
+                   a(10,6) * var(:,nj-9,:) ) * idel
+    rhs(:,nj-4,:) = rhs(:,nj-4,:) + &
+                 ( a(1,5) * var(:,nj,:) + a(2,5) * var(:,nj-1,:) + &
+                   a(3,5) * var(:,nj-2,:) + a(4,5) * var(:,nj-3,:) + &
+                   a(5,5) * var(:,nj-4,:) + a(6,5) * var(:,nj-5,:) + &
+                   a(7,5) * var(:,nj-6,:) + a(8,5) * var(:,nj-7,:) + &
+                   a(9,5) * var(:,nj-8,:) ) * idel
+    rhs(:,nj-3,:) = rhs(:,nj-3,:) + &
+                 ( a(1,4) * var(:,nj,:) + a(2,4) * var(:,nj-1,:) + &
+                   a(3,4) * var(:,nj-2,:) + a(4,4) * var(:,nj-3,:) + &
+                   a(5,4) * var(:,nj-4,:) + a(6,4) * var(:,nj-5,:) + &
+                   a(7,4) * var(:,nj-6,:) + a(8,4) * var(:,nj-7,:) ) * idel
+    rhs(:,nj-2,:) = rhs(:,nj-2,:) + &
+                 ( a(1,3) * var(:,nj,:) + a(2,3) * var(:,nj-1,:) + &
+                   a(3,3) * var(:,nj-2,:) + a(4,3) * var(:,nj-3,:) + &
+                   a(5,3) * var(:,nj-4,:) + a(6,3) * var(:,nj-5,:) + &
+                   a(7,3) * var(:,nj-6,:) ) * idel
+    rhs(:,nj-1,:) = rhs(:,nj-1,:) + &
+                 ( a(1,2) * var(:,nj,:) + a(2,2) * var(:,nj-1,:) + &
+                   a(3,2) * var(:,nj-2,:) + a(4,2) * var(:,nj-3,:) + &
+                   a(5,2) * var(:,nj-4,:) + a(6,2) * var(:,nj-5,:) ) * idel
+    rhs(:,nj,:)   = rhs(:,nj,:) + &
+                 ( a(1,1) * var(:,nj,:) + a(2,1) * var(:,nj-1,:) + &
+                   a(3,1) * var(:,nj-2,:) + a(4,1) * var(:,nj-3,:) + &
+                   a(5,1) * var(:,nj-4,:) ) * idel
+
+    jr = nj - 6
+  end if
+  rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+                   ( -70.0_wp * var(:,jl:jr,:) + &
+                      56.0_wp * ( var(:,jl-1:jr-1,:) + &
+                                  var(:,jl+1:jr+1,:) ) - &
+                      28.0_wp * ( var(:,jl-2:jr-2,:) + &
+                                  var(:,jl+2:jr+2,:) ) + &
+                       8.0_wp * ( var(:,jl-3:jr-3,:) + &
+                                  var(:,jl+3:jr+3,:) ) - &
+                                ( var(:,jl-4:jr-4,:) + &
+                                  var(:,jl+4:jr+4,:) ) ) * idel
+
+  call set_coeff ( a )
+
+  idel = epsilon / ( 256 * delta(3) )
+
+  if ( bb(5) == 0 ) then
+    kl = 1 + gsize(3)
+  else
+    rhs(:,:,1) = rhs(:,:,1) + &
+                 ( a(1,1) * var(:,:,1) + a(2,1) * var(:,:,2) + &
+                   a(3,1) * var(:,:,3) + a(4,1) * var(:,:,4) + &
+                   a(5,1) * var(:,:,5) ) * idel
+    rhs(:,:,2) = rhs(:,:,2) + &
+                 ( a(1,2) * var(:,:,1) + a(2,2) * var(:,:,2) + &
+                   a(3,2) * var(:,:,3) + a(4,2) * var(:,:,4) + &
+                   a(5,2) * var(:,:,5) + a(6,2) * var(:,:,6) ) * idel
+    rhs(:,:,3) = rhs(:,:,3) + &
+                 ( a(1,3) * var(:,:,1) + a(2,3) * var(:,:,2) + &
+                   a(3,3) * var(:,:,3) + a(4,3) * var(:,:,4) + &
+                   a(5,3) * var(:,:,5) + a(6,3) * var(:,:,6) + &
+                   a(7,3) * var(:,:,7) ) * idel
+    rhs(:,:,4) = rhs(:,:,4) + &
+                 ( a(1,4) * var(:,:,1) + a(2,4) * var(:,:,2) + &
+                   a(3,4) * var(:,:,3) + a(4,4) * var(:,:,4) + &
+                   a(5,4) * var(:,:,5) + a(6,4) * var(:,:,6) + &
+                   a(7,4) * var(:,:,7) + a(8,4) * var(:,:,8) ) * idel
+    rhs(:,:,5) = rhs(:,:,5) + &
+                 ( a(1,5) * var(:,:,1) + a(2,5) * var(:,:,2) + &
+                   a(3,5) * var(:,:,3) + a(4,5) * var(:,:,4) + &
+                   a(5,5) * var(:,:,5) + a(6,5) * var(:,:,6) + &
+                   a(7,5) * var(:,:,7) + a(8,5) * var(:,:,8) + &
+                   a(9,5) * var(:,:,9) ) * idel
+    rhs(:,:,6) = rhs(:,:,6) + &
+                 ( a(2,6) * var(:,:,2) + a(3,6) * var(:,:,3) + &
+                   a(4,6) * var(:,:,4) + a(5,6) * var(:,:,5) + &
+                   a(6,6) * var(:,:,6) + a(7,6) * var(:,:,7) + &
+                   a(8,6) * var(:,:,8) + a(9,6) * var(:,:,9) + &
+                   a(10,6) * var(:,:,10) ) * idel
+
+    kl = 7
+  end if
+  if ( bb(6) == 0 ) then
+    kr = nk - gsize(3)
+  else
+    rhs(:,:,nk-5) = rhs(:,:,nk-5) + &
+                 ( a(2,6) * var(:,:,nk-1) + a(3,6) * var(:,:,nk-2) + &
+                   a(4,6) * var(:,:,nk-3) + a(5,6) * var(:,:,nk-4) + &
+                   a(6,6) * var(:,:,nk-5) + a(7,6) * var(:,:,nk-6) + &
+                   a(8,6) * var(:,:,nk-7) + a(9,6) * var(:,:,nk-8) + &
+                   a(10,6) * var(:,:,nk-9) ) * idel
+    rhs(:,:,nk-4) = rhs(:,:,nk-4) + &
+                 ( a(1,5) * var(:,:,nk) + a(2,5) * var(:,:,nk-1) + &
+                   a(3,5) * var(:,:,nk-2) + a(4,5) * var(:,:,nk-3) + &
+                   a(5,5) * var(:,:,nk-4) + a(6,5) * var(:,:,nk-5) + &
+                   a(7,5) * var(:,:,nk-6) + a(8,5) * var(:,:,nk-7) + &
+                   a(9,5) * var(:,:,nk-8) ) * idel
+    rhs(:,:,nk-3) = rhs(:,:,nk-3) + &
+                 ( a(1,4) * var(:,:,nk) + a(2,4) * var(:,:,nk-1) + &
+                   a(3,4) * var(:,:,nk-2) + a(4,4) * var(:,:,nk-3) + &
+                   a(5,4) * var(:,:,nk-4) + a(6,4) * var(:,:,nk-5) + &
+                   a(7,4) * var(:,:,nk-6) + a(8,4) * var(:,:,nk-7) ) * idel
+    rhs(:,:,nk-2) = rhs(:,:,nk-2) + &
+                 ( a(1,3) * var(:,:,nk) + a(2,3) * var(:,:,nk-1) + &
+                   a(3,3) * var(:,:,nk-2) + a(4,3) * var(:,:,nk-3) + &
+                   a(5,3) * var(:,:,nk-4) + a(6,3) * var(:,:,nk-5) + &
+                   a(7,3) * var(:,:,nk-6) ) * idel
+    rhs(:,:,nk-1) = rhs(:,:,nk-1) + &
+                 ( a(1,2) * var(:,:,nk) + a(2,2) * var(:,:,nk-1) + &
+                   a(3,2) * var(:,:,nk-2) + a(4,2) * var(:,:,nk-3) + &
+                   a(5,2) * var(:,:,nk-4) + a(6,2) * var(:,:,nk-5) ) * idel
+    rhs(:,:,nk)   = rhs(:,:,nk) + &
+                 ( a(1,1) * var(:,:,nk) + a(2,1) * var(:,:,nk-1) + &
+                   a(3,1) * var(:,:,nk-2) + a(4,1) * var(:,:,nk-3) + &
+                   a(5,1) * var(:,:,nk-4) ) * idel
+    kr = nk - 6
+  end if
+  rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+                   ( -70.0_wp * var(:,:,kl:kr) + &
+                      56.0_wp * ( var(:,:,kl-1:kr-1) + &
+                                  var(:,:,kl+1:kr+1) ) - &
+                      28.0_wp * ( var(:,:,kl-2:kr-2) + &
+                                  var(:,:,kl+2:kr+2) ) + &
+                       8.0_wp * ( var(:,:,kl-3:kr-3) + &
+                                  var(:,:,kl+3:kr+3) ) - &
+                                ( var(:,:,kl-4:kr-4) + &
+                                  var(:,:,kl+4:kr+4) ) ) * idel
+contains
+
+  subroutine set_coeff ( a )
+
+    implicit none
+
+    CCTK_REAL, dimension(10,6), intent(OUT) :: a
+    CCTK_REAL :: zero = 0.0
+    integer, parameter :: wp = kind(zero)
+
+    a(1,1) = -3.1650670378782328375705179866656897941241116565316_wp
+    a(2,1) = 12.660268151512931350282071946662759176496446626126_wp
+    a(3,1) = -18.990402227269397025423107919994138764744669939190_wp
+    a(4,1) = 12.660268151512931350282071946662759176496446626126_wp
+    a(5,1) = -3.1650670378782328375705179866656897941241116565316_wp
+    a(6,1) = 0.0_wp
+    a(7,1) = 0.0_wp
+    a(8,1) = 0.0_wp
+    a(9,1) = 0.0_wp
+    a(10,1) = 0.0_wp
+    a(1,2) = 2.8768833763422958461666527928077915591442603845834_wp
+    a(2,2) = -12.226754349454757346208274369433114126363106634479_wp
+    a(3,2) = 20.138183634396070923166569549654540914009822692084_wp
+    a(4,2) = -15.822858569882627153916590360442853575293432115209_wp
+    a(5,2) = 5.7537667526845916923333055856155831182885207691667_wp
+    a(6,2) = -0.71922084408557396154166319820194788978606509614584_wp
+    a(7,2) = 0.0_wp
+    a(8,2) = 0.0_wp
+    a(9,2) = 0.0_wp
+    a(10,2) = 0.0_wp
+    a(1,3) = -9.5610475839173736628550350424197713021025451862781_wp
+    a(2,3) = 44.618222058281077093323496864625599409811877535965_wp
+    a(3,3) = -84.455920324603467355219476208041313168572482478790_wp
+    a(4,3) = 82.862412393950571744743637034304684618222058281077_wp
+    a(5,3) = -44.618222058281077093323496864625599409811877535965_wp
+    a(6,3) = 12.748063445223164883806713389893028402803393581704_wp
+    a(7,3) = -1.5935079306528956104758391737366285503504241977130_wp
+    a(8,3) = 0.0_wp
+    a(9,3) = 0.0_wp
+    a(10,3) = 0.0_wp
+    a(1,4) = 3.2244821795111028176898675125956335137152453816010_wp
+    a(2,4) = -17.734651987311065497294271319275984325433849598806_wp
+    a(3,4) = 41.918268333644336629968277663743235678298189960814_wp
+    a(4,4) = -55.622317596566523605150214592274678111587982832618_wp
+    a(5,4) = 45.142750513155439447658145176338869192013435342415_wp
+    a(6,4) = -22.571375256577719723829072588169434596006717671207_wp
+    a(7,4) = 6.4489643590222056353797350251912670274304907632021_wp
+    a(8,4) = -0.80612054487777570442246687814890837842881134540026_wp
+    a(9,4) = 0.0_wp
+    a(10,4) = 0.0_wp
+    a(1,5) = -1.0968642884346832550463374381109559476958232829758_wp
+    a(2,5) = 8.7749143074774660403706995048876475815665862638060_wp
+    a(3,5) = -30.712200076171131141297448267106766535483051923321_wp
+    a(4,5) = 61.424400152342262282594896534213533070966103846642_wp
+    a(5,5) = -76.780500190427827853243620667766916338707629808303_wp
+    a(6,5) = 61.424400152342262282594896534213533070966103846642_wp
+    a(7,5) = -30.712200076171131141297448267106766535483051923321_wp
+    a(8,5) = 8.7749143074774660403706995048876475815665862638060_wp
+    a(9,5) = -1.0968642884346832550463374381109559476958232829758_wp
+    a(10,5) = 0.0_wp
+    a(1,6) = 0.0_wp
+    a(2,6) = -0.98627885208100271683294901942878016483641925983425_wp
+    a(3,6) = 7.8902308166480217346635921554302413186913540786740_wp
+    a(4,6) = -27.615807858268076071322572544005844615419739275359_wp
+    a(5,6) = 55.231615716536152142645145088011689230839478550718_wp
+    a(6,6) = -69.039519645670190178306431360014611538549348188398_wp
+    a(7,6) = 55.231615716536152142645145088011689230839478550718_wp
+    a(8,6) = -27.615807858268076071322572544005844615419739275359_wp
+    a(9,6) = 7.8902308166480217346635921554302413186913540786740_wp
+    a(10,6) = -0.98627885208100271683294901942878016483641925983425_wp
+
+  end subroutine set_coeff
+
+end subroutine dissipation_6_3_alt
diff --git a/src/dissipation.c b/src/dissipation.c
index 29dff9f..9c8a8a9 100644
--- a/src/dissipation.c
+++ b/src/dissipation.c
@@ -192,15 +192,27 @@ apply (int const varindex, char const * const optstring, void * const arg)
       break;
     }
     case 4: {
-      CCTK_FNAME(dissipation_4_2)
-        (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+      if ( CCTK_Equals(dissipation_type,"Mattson-Svard-Nordstrom") ) {
+        CCTK_FNAME(dissipation_4_2)
+          (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+                            bbox, gsize, dx, &epsdis, rhsptr);
+      } else {
+        CCTK_FNAME(dissipation_4_2_alt)
+          (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
                             bbox, gsize, dx, &epsdis, rhsptr);
+      }
       break;
     }
     case 6: {
-      CCTK_FNAME(dissipation_6_3)
-        (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+      if ( CCTK_Equals(dissipation_type,"Mattson-Svard-Nordstrom") ) {
+        CCTK_FNAME(dissipation_6_3)
+          (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+                            bbox, gsize, dx, &epsdis, rhsptr);
+      } else {
+        CCTK_FNAME(dissipation_6_3_alt)
+          (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
                             bbox, gsize, dx, &epsdis, rhsptr);
+      }
       break;
     }
     case 8: {
diff --git a/src/make.code.defn b/src/make.code.defn
index 87f521e..50857c5 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -21,7 +21,9 @@ SRCS =  call_derivs.c \
         Dissipation_2_1.F90 \
         Dissipation_2_1_alt.F90 \
         Dissipation_4_2.F90 \
+        Dissipation_4_2_alt.F90 \
         Dissipation_6_3.F90 \
+        Dissipation_6_3_alt.F90 \
         Dissipation_8_4.F90 \
         Dissipation_8_4_alt.F90 \
         Dissipation_4_3_min_err_coeff.F90 \
-- 
cgit v1.2.3