New 4-3 and 6-5 operators that are optimised for accuracy. Added a new

parameter to choose between minimim bandwidth and optimized operators.


git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@31 f69c4107-0314-4c4f-9ad4-17e986b73f4a

4 files changed, 647 insertions, 1 deletions

diff --git a/param.ccl b/param.ccl
index fc70283..0866282 100644
--- a/param.ccl
+++ b/param.ccl
@@ -7,6 +7,12 @@ KEYWORD norm_type "Type of norm" STEERABLE=always
   "Full restricted" :: "Full resricted norm"
 } "Diagonal"
 
+KEYWORD operator_type "Type of operator" STEERABLE=always
+{
+  "Minimal Bandwidth" :: "Minimal bandwidth"
+  "Optimized" :: "Optimized for performance"
+} "Minimal Bandwidth"
+
 INT order "Order of accuracy" STEERABLE=always
 {
   2:8:2 :: ""
diff --git a/src/Derivatives_4_3_min_err_coeff.F90 b/src/Derivatives_4_3_min_err_coeff.F90
new file mode 100644
index 0000000..938c4bc
--- /dev/null
+++ b/src/Derivatives_4_3_min_err_coeff.F90
@@ -0,0 +1,235 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine deriv_gf_4_3_opt ( var, ni, nj, nk, dir, bb, gsize, delta, dvar )
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  CCTK_INT, intent(IN) :: ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var
+  CCTK_INT, intent(IN) :: dir
+  CCTK_INT, intent(IN) :: bb(2)
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_REAL, intent(IN) :: delta
+  CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar
+
+  CCTK_REAL, dimension(2), save :: a
+  CCTK_REAL, dimension(7,5), save :: q 
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    a(1) = 2.0_wp/3.0_wp; a(2) = -1.0_wp/12.0_wp
+
+    q(1,1) = -2.09329763466349871588733_wp
+    q(2,1) = 4.0398572053206615302160_wp
+    q(3,1) = -3.0597858079809922953240_wp
+    q(4,1) = 1.37319053865399486354933_wp
+    q(5,1) = -0.25996430133016538255400_wp
+    q(6,1) = 0_wp
+    q(7,1) = 0_wp
+    q(1,2) = -0.31641585285940445272297_wp
+    q(2,2) = -0.53930788973980422327388_wp
+    q(3,2) = 0.98517732028644343383297_wp
+    q(4,2) = -0.05264665989297578146709_wp
+    q(5,2) = -0.113807251750624235013258_wp
+    q(6,2) = 0.039879767889849911803103_wp
+    q(7,2) = -0.0028794339334846531588787_wp
+    q(1,3) = 0.13026916185021164524452_wp
+    q(2,3) = -0.87966858995059249256890_wp
+    q(3,3) = 0.38609640961100070000134_wp
+    q(4,3) = 0.31358369072435588745988_wp
+    q(5,3) = 0.085318941913678384633511_wp
+    q(6,3) = -0.039046615792734640274641_wp
+    q(7,3) = 0.0034470016440805155042908_wp
+    q(1,4) = -0.01724512193824647912172_wp
+    q(2,4) = 0.16272288227127504381134_wp
+    q(3,4) = -0.81349810248648813029217_wp
+    q(4,4) = 0.13833269266479833215645_wp
+    q(5,4) = 0.59743854328548053399616_wp
+    q(6,4) = -0.066026434346299887619324_wp
+    q(7,4) = -0.0017244594505194129307249_wp
+    q(1,5) = -0.00883569468552192965061_wp
+    q(2,5) = 0.03056074759203203857284_wp
+    q(3,5) = 0.05021168274530854232278_wp
+    q(4,5) = -0.66307364652444929534068_wp
+    q(5,5) = 0.014878787464005191116088_wp
+    q(6,5) = 0.65882706381707471953820_wp
+    q(7,5) = -0.082568940408449266558615_wp
+
+    first = .false.
+  end if
+
+  idel = 1.0_wp / delta
+
+  direction: select case (dir)
+  case (0) direction
+    if ( bb(1) == 0 ) then
+      il = 1 + gsize
+    else
+      dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + &
+                      q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) + &
+                      q(5,1) * var(5,:,:) ) * idel
+      dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + &
+                      q(3,2) * var(3,:,:) + q(4,2) * var(4,:,:) + &
+                      q(5,2) * var(5,:,:) + q(6,2) * var(6,:,:) + &
+                      q(7,2) * var(7,:,:) ) * idel
+      dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + &
+                      q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + &
+                      q(5,3) * var(5,:,:) + q(6,3) * var(6,:,:) + &
+                      q(7,3) * var(7,:,:) ) * idel
+      dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + &
+                      q(3,4) * var(3,:,:) + q(4,4) * var(4,:,:) + &
+                      q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) + &
+                      q(7,4) * var(7,:,:) ) * idel
+      dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + &
+                      q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + &
+                      q(5,5) * var(5,:,:) + q(6,5) * var(6,:,:) + &
+                      q(7,5) * var(7,:,:) ) * idel
+      il = 6
+    end if
+    if ( bb(2) == 0 ) then
+      ir = ni - gsize
+    else
+      dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + &
+                           q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + &
+                           q(5,5) * var(ni-4,:,:) + q(6,5) * var(ni-5,:,:) + &
+                           q(7,5) * var(ni-6,:,:) ) * idel
+      dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + &
+                           q(3,4) * var(ni-2,:,:) + q(4,4) * var(ni-3,:,:) + &
+                           q(5,4) * var(ni-4,:,:) + q(6,4) * var(ni-5,:,:) + &
+                           q(7,4) * var(ni-6,:,:) ) * idel
+      dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + &
+                           q(3,3) * var(ni-2,:,:) + q(4,3) * var(ni-3,:,:) + &
+                           q(5,3) * var(ni-4,:,:) + q(6,3) * var(ni-5,:,:) + &
+                           q(7,3) * var(ni-6,:,:) ) * idel
+      dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(2,2) * var(ni-1,:,:) + &
+                           q(3,2) * var(ni-2,:,:) + q(4,2) * var(ni-3,:,:) + &
+                           q(5,2) * var(ni-4,:,:) + q(6,2) * var(ni-5,:,:) + &
+                           q(7,2) * var(ni-6,:,:) ) * idel
+      dvar(ni,:,:)   = - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + &
+                           q(3,1) * var(ni-2,:,:) + q(4,1) * var(ni-3,:,:) + &
+                           q(5,1) * var(ni-4,:,:) ) * idel
+      ir = ni - 5
+    end if
+    dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - &
+                                 var(il-1:ir-1,:,:) ) + &
+                        a(2) * ( var(il+2:ir+2,:,:) - &
+                                 var(il-2:ir-2,:,:) ) ) * idel
+  case (1) direction
+    if ( bb(1) == 0 ) then
+      jl = 1 + gsize
+    else
+      dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + &
+                      q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) + &
+                      q(5,1) * var(:,5,:) ) * idel
+      dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + &
+                      q(3,2) * var(:,3,:) + q(4,2) * var(:,4,:) + &
+                      q(5,2) * var(:,5,:) + q(6,2) * var(:,6,:) + &
+                      q(7,2) * var(:,7,:) ) * idel
+      dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + &
+                      q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + &
+                      q(5,3) * var(:,5,:) + q(6,3) * var(:,6,:) + &
+                      q(7,3) * var(:,7,:) ) * idel
+      dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + &
+                      q(3,4) * var(:,3,:) + q(4,4) * var(:,4,:) + &
+                      q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) + &
+                      q(7,4) * var(:,7,:) ) * idel
+      dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + &
+                      q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + &
+                      q(5,5) * var(:,5,:) + q(6,5) * var(:,6,:) + &
+                      q(7,5) * var(:,7,:) ) * idel
+      jl = 6
+    end if
+    if ( bb(2) == 0 ) then
+      jr = nj - gsize
+    else
+      dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + &
+                           q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + &
+                           q(5,5) * var(:,nj-4,:) + q(6,5) * var(:,nj-5,:) + &
+                           q(7,5) * var(:,nj-6,:) ) * idel
+      dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + &
+                           q(3,4) * var(:,nj-2,:) + q(4,4) * var(:,nj-3,:) + &
+                           q(5,4) * var(:,nj-4,:) + q(6,4) * var(:,nj-5,:) + &
+                           q(7,4) * var(:,nj-6,:) ) * idel
+      dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + &
+                           q(3,3) * var(:,nj-2,:) + q(4,3) * var(:,nj-3,:) + &
+                           q(5,3) * var(:,nj-4,:) + q(6,3) * var(:,nj-5,:) + &
+                           q(7,3) * var(:,nj-6,:) ) * idel
+      dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(2,2) * var(:,nj-1,:) + &
+                           q(3,2) * var(:,nj-2,:) + q(4,2) * var(:,nj-3,:) + &
+                           q(5,2) * var(:,nj-4,:) + q(6,2) * var(:,nj-5,:) + &
+                           q(7,2) * var(:,nj-6,:) ) * idel
+      dvar(:,nj,:)   = - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + &
+                           q(3,1) * var(:,nj-2,:) + q(4,1) * var(:,nj-3,:) + &
+                           q(5,1) * var(:,nj-4,:) ) * idel
+      jr = nj - 5
+    end if
+    dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - &
+                                 var(:,jl-1:jr-1,:) ) + &
+                        a(2) * ( var(:,jl+2:jr+2,:) - &
+                                 var(:,jl-2:jr-2,:) ) ) * idel
+  case (2) direction
+    if ( bb(1) == 0 ) then
+      kl = 1 + gsize
+    else
+      dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + &
+                      q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) + &
+                      q(5,1) * var(:,:,5) ) * idel
+      dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + &
+                      q(3,2) * var(:,:,3) + q(4,2) * var(:,:,4) + &
+                      q(5,2) * var(:,:,5) + q(6,2) * var(:,:,6) + &
+                      q(7,2) * var(:,:,7) ) * idel
+      dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + &
+                      q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + &
+                      q(5,3) * var(:,:,5) + q(6,3) * var(:,:,6) + &
+                      q(7,3) * var(:,:,7) ) * idel
+      dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + &
+                      q(3,4) * var(:,:,3) + q(4,4) * var(:,:,4) + &
+                      q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) + &
+                      q(7,4) * var(:,:,7) ) * idel
+      dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + &
+                      q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + &
+                      q(5,5) * var(:,:,5) + q(6,5) * var(:,:,6) + &
+                      q(7,5) * var(:,:,7) ) * idel
+      kl = 6
+    end if
+    if ( bb(2) == 0 ) then
+      kr = nk - gsize
+    else 
+      dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + &
+                           q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + &
+                           q(5,5) * var(:,:,nk-4) + q(6,5) * var(:,:,nk-5) + &
+                           q(7,5) * var(:,:,nk-6) ) * idel
+      dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + &
+                           q(3,4) * var(:,:,nk-2) + q(4,4) * var(:,:,nk-3) + &
+                           q(5,4) * var(:,:,nk-4) + q(6,4) * var(:,:,nk-5) + &
+                           q(7,4) * var(:,:,nk-6) ) * idel
+      dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + &
+                           q(3,3) * var(:,:,nk-2) + q(4,3) * var(:,:,nk-3) + &
+                           q(5,3) * var(:,:,nk-4) + q(6,3) * var(:,:,nk-5) + &
+                           q(7,3) * var(:,:,nk-6) ) * idel
+      dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(2,2) * var(:,:,nk-1) + &
+                           q(3,2) * var(:,:,nk-2) + q(4,2) * var(:,:,nk-3) + &
+                           q(5,2) * var(:,:,nk-4) + q(6,2) * var(:,:,nk-5) + &
+                           q(7,2) * var(:,:,nk-6) ) * idel
+      dvar(:,:,nk)   = - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + &
+                           q(3,1) * var(:,:,nk-2) + q(4,1) * var(:,:,nk-3) + &
+                           q(5,1) * var(:,:,nk-4) ) * idel
+      kr = nk - 5
+    end if
+    dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - &
+                                 var(:,:,kl-1:kr-1) ) + &
+                        a(2) * ( var(:,:,kl+2:kr+2) - &
+                                 var(:,:,kl-2:kr-2) ) ) * idel
+  end select direction
+end subroutine deriv_gf_4_3_opt
diff --git a/src/Derivatives_6_5_min_err_coeff.F90 b/src/Derivatives_6_5_min_err_coeff.F90
new file mode 100644
index 0000000..5bc1a2f
--- /dev/null
+++ b/src/Derivatives_6_5_min_err_coeff.F90
@@ -0,0 +1,384 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine deriv_gf_6_5_opt ( var, ni, nj, nk, dir, bb, gsize, delta, dvar )
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  CCTK_INT, intent(IN) :: ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var
+  CCTK_INT, intent(IN) :: dir
+  CCTK_INT, intent(IN) :: bb(2)
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_REAL, intent(IN) :: delta
+  CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar
+
+  CCTK_REAL, dimension(3), save :: a
+  CCTK_REAL, dimension(10,7), save :: q 
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    a(1) = 3.0_wp/4.0_wp; a(2) = -3.0_wp/20.0_wp; a(3) = 1.0_wp/60.0_wp
+
+    q(1,1) = -2.465354921105726344726324628220306540433_wp
+    q(2,1) = 6.092129526634358068357947769321839242595_wp
+    q(3,1) = -7.730323816585895170894869423304598106488_wp
+    q(4,1) = 6.973765088781193561193159231072797475318_wp
+    q(5,1) = -3.980323816585895170894869423304598106488_wp
+    q(6,1) = 1.292129526634358068357947769321839242595_wp
+    q(7,1) = -0.1820215877723930113929912948869732070992_wp
+    q(8,1) = 0_wp
+    q(9,1) = 0_wp
+    q(10,1) = 0_wp
+    q(1,2) = -0.2234725650933373070210703921399196916276_wp
+    q(2,2) = -0.9329308120213103664703924027998434010527_wp
+    q(3,2) = 1.586820596323380726926380762297803979729_wp
+    q(4,2) = -0.3647002337446073060093577955438077228163_wp
+    q(5,2) = -0.2666957787028035557343386019032589667016_wp
+    q(6,2) = 0.3112949049503615973044042012025239373575_wp
+    q(7,2) = -0.1404504214992091753318473921133396942483_wp
+    q(8,2) = 0.03488568515844113812839530849342152163114_wp
+    q(9,2) = -0.004964021892185918282918047104381964571176_wp
+    q(10,2) = 0.0002126465212701664907443596108020023005085_wp
+    q(1,3) = 0.1582216736972761112130174221319685578703_wp
+    q(2,3) = -1.137049297945701138153781830545857828001_wp
+    q(3,3) = 1.212364522772435724261769516755293203384_wp
+    q(4,3) = -0.9562288727032211388608555534303520293292_wp
+    q(5,3) = 1.066548057100650347439560790660954981316_wp
+    q(6,3) = -0.3478788549812286788557041220360802499887_wp
+    q(7,3) = -0.03133923299767926953262771814410417142710_wp
+    q(8,3) = 0.04098845957876514538943748816651847063573_wp
+    q(9,3) = -0.005963188640428636310134738586972978148865_wp
+    q(10,3) = 0.0003367341191315334093187450286320436894509_wp
+    q(1,4) = 0.02915734643197329316878786487755959617414_wp
+    q(2,4) = -0.1169665090568385435258745132657230868726_wp
+    q(3,4) = -0.1112219090396860235587194333712503100976_wp
+    q(4,4) = -0.7924486264107931688501411569472360917480_wp
+    q(5,4) = 1.266650705052489006015656810787656809460_wp
+    q(6,4) = -0.2899273291614735652270415777078748545257_wp
+    q(7,4) = 0.002515684289746880247060269725567464781645_wp
+    q(8,4) = 0.01329713961045538259784924004687749664184_wp
+    q(9,4) = -0.001124464397192287407208734867324966241272_wp
+    q(10,4) = 0.00006796268131902653963123072174794242752627_wp
+    q(1,5) = -0.04582150000219473783416432201491386259525_wp
+    q(2,5) = 0.2240986548763910403084011643939472909730_wp
+    q(3,5) = -0.3246718492669968647456547557775778272008_wp
+    q(4,5) = -0.3929792922479188555634359057335494798829_wp
+    q(5,5) = 0.1166355819597526695402713402151963096924_wp
+    q(6,5) = 0.3449626905274072306980547367786542540315_wp
+    q(7,5) = 0.1430419813692381466377067274780162668130_wp
+    q(8,5) = -0.07764802500390609299919055841918594590836_wp
+    q(9,5) = 0.01332439335682098023215359404071263206533_wp
+    q(10,5) = -0.0009426355685935162741420209612996379880357_wp
+    q(1,6) = 0.003172814427534452820786948997021353279521_wp
+    q(2,6) = 0.00001061461750788952405498488693810144188215_wp
+    q(3,6) = -0.08747763621208721666864461344012255496503_wp
+    q(4,6) = 0.3975827328151903159828221701608106653957_wp
+    q(5,6) = -1.148835072887093214621630866015514532074_wp
+    q(6,6) = 0.3583006652024372613320453990260111557237_wp
+    q(7,6) = 0.5647665153537601810077289284875161241444_wp
+    q(8,6) = -0.09698196886045636272181593036009187114508_wp
+    q(9,6) = 0.008843905090944639404474662278460308660277_wp
+    q(10,6) = 0.0006174304522620539401783159789712495382110_wp
+    q(1,7) = -0.008639107520047315872218762260033429623639_wp
+    q(2,7) = 0.04722773941989840822874396849234369996233_wp
+    q(3,7) = -0.1008747534523833098882185784719780553451_wp
+    q(4,7) = 0.08043834912116582462744784937371906752091_wp
+    q(5,7) = 0.1295138677847054989778048608621896918103_wp
+    q(6,7) = -0.7909424167745257943256055480973493995702_wp
+    q(7,7) = 0.03807866849747788018644730375057882400112_wp
+    q(8,7) = 0.7367055699547762148521594995917960242549_wp
+    q(9,7) = -0.1480235854664835893774029568937961320792_wp
+    q(10,7) = 0.01651566843541618259084236365252970906860_wp
+
+    first = .false.
+  end if
+
+  idel = 1.0_wp / delta
+
+  direction: select case (dir)
+  case (0) direction
+    if ( bb(1) == 0 ) then
+      il = 1 + gsize
+    else
+      dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + &
+                      q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) + &
+                      q(5,1) * var(5,:,:) + q(6,1) * var(6,:,:) + &
+                      q(7,1) * var(7,:,:) ) * idel
+      dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + &
+                      q(3,2) * var(3,:,:) + q(4,2) * var(4,:,:) + &
+                      q(5,2) * var(5,:,:) + q(6,2) * var(6,:,:) + &
+                      q(7,2) * var(7,:,:) + q(8,2) * var(8,:,:) + &
+                      q(9,2) * var(9,:,:) + q(10,2) * var(10,:,:) ) * idel
+      dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + &
+                      q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + &
+                      q(5,3) * var(5,:,:) + q(6,3) * var(6,:,:) + &
+                      q(7,3) * var(7,:,:) + q(8,3) * var(8,:,:) + &
+                      q(9,3) * var(9,:,:) + q(10,3) * var(10,:,:) ) * idel
+      dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + &
+                      q(3,4) * var(3,:,:) + q(4,4) * var(4,:,:) + &
+                      q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) + &
+                      q(7,4) * var(7,:,:) + q(8,4) * var(8,:,:) + &
+                      q(9,4) * var(9,:,:) + q(10,4) * var(10,:,:) ) * idel
+      dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + &
+                      q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + &
+                      q(5,5) * var(5,:,:) + q(6,5) * var(6,:,:) + &
+                      q(7,5) * var(7,:,:) + q(8,5) * var(8,:,:) + &
+                      q(9,5) * var(9,:,:) + q(10,5) * var(10,:,:) ) * idel
+      dvar(6,:,:) = ( q(1,6) * var(1,:,:) + q(2,6) * var(2,:,:) + &
+                      q(3,6) * var(3,:,:) + q(4,6) * var(4,:,:) + &
+                      q(5,6) * var(5,:,:) + q(6,6) * var(6,:,:) + &
+                      q(7,6) * var(7,:,:) + q(8,6) * var(8,:,:) + &
+                      q(9,6) * var(9,:,:) + q(10,6) * var(10,:,:) ) * idel
+      dvar(7,:,:) = ( q(1,7) * var(1,:,:) + q(2,7) * var(2,:,:) + &
+                      q(3,7) * var(3,:,:) + q(4,7) * var(4,:,:) + &
+                      q(5,7) * var(5,:,:) + q(6,7) * var(6,:,:) + &
+                      q(7,7) * var(7,:,:) + q(8,7) * var(8,:,:) + &
+                      q(9,7) * var(9,:,:) + q(10,7) * var(10,:,:) ) * idel
+      il = 8
+    end if
+    if ( bb(2) == 0 ) then
+      ir = ni - gsize
+    else
+      dvar(ni-6,:,:) = - ( q(1,7) * var(ni,:,:) + q(2,7) * var(ni-1,:,:) + &
+                           q(3,7) * var(ni-2,:,:) + q(4,7) * var(ni-3,:,:) + &
+                           q(5,7) * var(ni-4,:,:) + q(6,7) * var(ni-5,:,:) + &
+                           q(7,7) * var(ni-6,:,:) + q(8,7) * var(ni-7,:,:) + &
+                           q(9,7) * var(ni-8,:,:) + &
+                           q(10,7) * var(ni-9,:,:) ) * idel
+      dvar(ni-5,:,:) = - ( q(1,6) * var(ni,:,:) + q(2,6) * var(ni-1,:,:) + &
+                           q(3,6) * var(ni-2,:,:) + q(4,6) * var(ni-3,:,:) + &
+                           q(5,6) * var(ni-4,:,:) + q(6,6) * var(ni-5,:,:) + &
+                           q(7,6) * var(ni-6,:,:) + q(8,6) * var(ni-7,:,:) + &
+                           q(9,6) * var(ni-8,:,:) + &
+                           q(10,6) * var(ni-9,:,:) ) * idel
+      dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + &
+                           q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + &
+                           q(5,5) * var(ni-4,:,:) + q(6,5) * var(ni-5,:,:) + &
+                           q(7,5) * var(ni-6,:,:) + q(8,5) * var(ni-7,:,:) + &
+                           q(9,5) * var(ni-8,:,:) + &
+                           q(10,5) * var(ni-9,:,:) ) * idel
+      dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + &
+                           q(3,4) * var(ni-2,:,:) + q(4,4) * var(ni-3,:,:) + &
+                           q(5,4) * var(ni-4,:,:) + q(6,4) * var(ni-5,:,:) + &
+                           q(7,4) * var(ni-6,:,:) + q(8,4) * var(ni-7,:,:) + &
+                           q(9,4) * var(ni-8,:,:) + &
+                           q(10,4) * var(ni-9,:,:) ) * idel
+      dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + &
+                           q(3,3) * var(ni-2,:,:) + q(4,3) * var(ni-3,:,:) + &
+                           q(5,3) * var(ni-4,:,:) + q(6,3) * var(ni-5,:,:) + &
+                           q(7,3) * var(ni-6,:,:) + q(8,3) * var(ni-7,:,:) + &
+                           q(9,3) * var(ni-8,:,:) + &
+                           q(10,3) * var(ni-9,:,:) ) * idel
+      dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(2,2) * var(ni-1,:,:) + &
+                           q(3,2) * var(ni-2,:,:) + q(4,2) * var(ni-3,:,:) + &
+                           q(5,2) * var(ni-4,:,:) + q(6,2) * var(ni-5,:,:) + &
+                           q(7,2) * var(ni-6,:,:) + q(8,2) * var(ni-7,:,:) + &
+                           q(9,2) * var(ni-8,:,:) + &
+                           q(10,2) * var(ni-9,:,:) ) * idel
+      dvar(ni,:,:)   = - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + &
+                           q(3,1) * var(ni-2,:,:) + q(4,1) * var(ni-3,:,:) + &
+                           q(5,1) * var(ni-4,:,:) + q(6,1) * var(ni-5,:,:) + &
+                           q(7,1) * var(ni-6,:,:) ) * idel
+      ir = ni - 7
+    end if
+    dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - &
+                                 var(il-1:ir-1,:,:) ) + &
+                        a(2) * ( var(il+2:ir+2,:,:) - &
+                                 var(il-2:ir-2,:,:) ) + &
+                        a(3) * ( var(il+3:ir+3,:,:) - &
+                                 var(il-3:ir-3,:,:) ) ) * idel
+  case (1) direction
+    if ( bb(1) == 0 ) then
+      jl = 1 + gsize
+    else
+      dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + &
+                      q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) + &
+                      q(5,1) * var(:,5,:) + q(6,1) * var(:,6,:) + &
+                      q(7,1) * var(:,7,:) ) * idel
+      dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + &
+                      q(3,2) * var(:,3,:) + q(4,2) * var(:,4,:) + &
+                      q(5,2) * var(:,5,:) + q(6,2) * var(:,6,:) + &
+                      q(7,2) * var(:,7,:) + q(8,2) * var(:,8,:) + &
+                      q(9,2) * var(:,9,:) + q(10,2) * var(:,10,:) ) * idel
+      dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + &
+                      q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + &
+                      q(5,3) * var(:,5,:) + q(6,3) * var(:,6,:) + &
+                      q(7,3) * var(:,7,:) + q(8,3) * var(:,8,:) + &
+                      q(9,3) * var(:,9,:) + q(10,3) * var(:,10,:) ) * idel
+      dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + &
+                      q(3,4) * var(:,3,:) + q(4,4) * var(:,4,:) + &
+                      q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) + &
+                      q(7,4) * var(:,7,:) + q(8,4) * var(:,8,:) + &
+                      q(9,4) * var(:,9,:) + q(10,4) * var(:,10,:) ) * idel
+      dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + &
+                      q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + &
+                      q(5,5) * var(:,5,:) + q(6,5) * var(:,6,:) + &
+                      q(7,5) * var(:,7,:) + q(8,5) * var(:,8,:) + &
+                      q(9,5) * var(:,9,:) + q(10,5) * var(:,10,:) ) * idel
+      dvar(:,6,:) = ( q(1,6) * var(:,1,:) + q(2,6) * var(:,2,:) + &
+                      q(3,6) * var(:,3,:) + q(4,6) * var(:,4,:) + &
+                      q(5,6) * var(:,5,:) + q(6,6) * var(:,6,:) + &
+                      q(7,6) * var(:,7,:) + q(8,6) * var(:,8,:) + &
+                      q(9,6) * var(:,9,:) + q(10,6) * var(:,10,:) ) * idel
+      dvar(:,7,:) = ( q(1,7) * var(:,1,:) + q(2,7) * var(:,2,:) + &
+                      q(3,7) * var(:,3,:) + q(4,7) * var(:,4,:) + &
+                      q(5,7) * var(:,5,:) + q(6,7) * var(:,6,:) + &
+                      q(7,7) * var(:,7,:) + q(8,7) * var(:,8,:) + &
+                      q(9,7) * var(:,9,:) + q(10,7) * var(:,10,:) ) * idel
+      jl = 8
+    end if
+    if ( bb(2) == 0 ) then
+      jr = nj - gsize
+    else
+      dvar(:,nj-6,:) = - ( q(1,7) * var(:,nj,:) + q(2,7) * var(:,nj-1,:) + &
+                           q(3,7) * var(:,nj-2,:) + q(4,7) * var(:,nj-3,:) + &
+                           q(5,7) * var(:,nj-4,:) + q(6,7) * var(:,nj-5,:) + &
+                           q(7,7) * var(:,nj-6,:) + q(8,7) * var(:,nj-7,:) + &
+                           q(9,7) * var(:,nj-8,:) + &
+                           q(10,7) * var(:,nj-9,:) ) * idel
+      dvar(:,nj-5,:) = - ( q(1,6) * var(:,nj,:) + q(2,6) * var(:,nj-1,:) + &
+                           q(3,6) * var(:,nj-2,:) + q(4,6) * var(:,nj-3,:) + &
+                           q(5,6) * var(:,nj-4,:) + q(6,6) * var(:,nj-5,:) + &
+                           q(7,6) * var(:,nj-6,:) + q(8,6) * var(:,nj-7,:) + &
+                           q(9,6) * var(:,nj-8,:) + &
+                           q(10,6) * var(:,nj-9,:) ) * idel
+      dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + &
+                           q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + &
+                           q(5,5) * var(:,nj-4,:) + q(6,5) * var(:,nj-5,:) + &
+                           q(7,5) * var(:,nj-6,:) + q(8,5) * var(:,nj-7,:) + &
+                           q(9,5) * var(:,nj-8,:) + &
+                           q(10,5) * var(:,nj-9,:) ) * idel
+      dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + &
+                           q(3,4) * var(:,nj-2,:) + q(4,4) * var(:,nj-3,:) + &
+                           q(5,4) * var(:,nj-4,:) + q(6,4) * var(:,nj-5,:) + &
+                           q(7,4) * var(:,nj-6,:) + q(8,4) * var(:,nj-7,:) + &
+                           q(9,4) * var(:,nj-8,:) + &
+                           q(10,4) * var(:,nj-9,:) ) * idel
+      dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + &
+                           q(3,3) * var(:,nj-2,:) + q(4,3) * var(:,nj-3,:) + &
+                           q(5,3) * var(:,nj-4,:) + q(6,3) * var(:,nj-5,:) + &
+                           q(7,3) * var(:,nj-6,:) + q(8,3) * var(:,nj-7,:) + &
+                           q(9,3) * var(:,nj-8,:) + &
+                           q(10,3) * var(:,nj-9,:) ) * idel
+      dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(2,2) * var(:,nj-1,:) + &
+                           q(3,2) * var(:,nj-2,:) + q(4,2) * var(:,nj-3,:) + &
+                           q(5,2) * var(:,nj-4,:) + q(6,2) * var(:,nj-5,:) + &
+                           q(7,2) * var(:,nj-6,:) + q(8,2) * var(:,nj-7,:) + &
+                           q(9,2) * var(:,nj-8,:) + &
+                           q(10,2) * var(:,nj-9,:) ) * idel
+      dvar(:,nj,:)   = - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + &
+                           q(3,1) * var(:,nj-2,:) + q(4,1) * var(:,nj-3,:) + &
+                           q(5,1) * var(:,nj-4,:) + q(6,1) * var(:,nj-5,:) + &
+                           q(7,1) * var(:,nj-6,:) ) * idel
+      jr = nj - 7
+    end if
+    dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - &
+                                 var(:,jl-1:jr-1,:) ) + &
+                        a(2) * ( var(:,jl+2:jr+2,:) - &
+                                 var(:,jl-2:jr-2,:) ) + &
+                        a(3) * ( var(:,jl+3:jr+3,:) - &
+                                 var(:,jl-3:jr-3,:) ) ) * idel
+  case (2) direction
+    if ( bb(1) == 0 ) then
+      kl = 1 + gsize
+    else
+      dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + &
+                      q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) + &
+                      q(5,1) * var(:,:,5) + q(6,1) * var(:,:,6) + &
+                      q(7,1) * var(:,:,7) ) * idel
+      dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + &
+                      q(3,2) * var(:,:,3) + q(4,2) * var(:,:,4) + &
+                      q(5,2) * var(:,:,5) + q(6,2) * var(:,:,6) + &
+                      q(7,2) * var(:,:,7) + q(8,2) * var(:,:,8) + &
+                      q(9,2) * var(:,:,9) + q(10,2) * var(:,:,10) ) * idel
+      dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + &
+                      q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + &
+                      q(5,3) * var(:,:,5) + q(6,3) * var(:,:,6) + &
+                      q(7,3) * var(:,:,7) + q(8,3) * var(:,:,8) + &
+                      q(9,3) * var(:,:,9) + q(10,3) * var(:,:,10) ) * idel
+      dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + &
+                      q(3,4) * var(:,:,3) + q(4,4) * var(:,:,4) + &
+                      q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) + &
+                      q(7,4) * var(:,:,7) + q(8,4) * var(:,:,8) + &
+                      q(9,4) * var(:,:,9) + q(10,4) * var(:,:,10) ) * idel
+      dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + &
+                      q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + &
+                      q(5,5) * var(:,:,5) + q(6,5) * var(:,:,6) + &
+                      q(7,5) * var(:,:,7) + q(8,5) * var(:,:,8) + &
+                      q(9,5) * var(:,:,9) + q(10,5) * var(:,:,10) ) * idel
+      dvar(:,:,6) = ( q(1,6) * var(:,:,1) + q(2,6) * var(:,:,2) + &
+                      q(3,6) * var(:,:,3) + q(4,6) * var(:,:,4) + &
+                      q(5,6) * var(:,:,5) + q(6,6) * var(:,:,6) + &
+                      q(7,6) * var(:,:,7) + q(8,6) * var(:,:,8) + &
+                      q(9,6) * var(:,:,9) + q(10,6) * var(:,:,10) ) * idel
+      dvar(:,:,7) = ( q(1,7) * var(:,:,1) + q(2,7) * var(:,:,2) + &
+                      q(3,7) * var(:,:,3) + q(4,7) * var(:,:,4) + &
+                      q(5,7) * var(:,:,5) + q(6,7) * var(:,:,6) + &
+                      q(7,7) * var(:,:,7) + q(8,7) * var(:,:,8) + &
+                      q(9,7) * var(:,:,9) + q(10,7) * var(:,:,10) ) * idel
+      kl = 8
+    end if
+    if ( bb(2) == 0 ) then
+      kr = nk - gsize
+    else 
+      dvar(:,:,nk-6) = - ( q(1,7) * var(:,:,nk) + q(2,7) * var(:,:,nk-1) + &
+                           q(3,7) * var(:,:,nk-2) + q(4,7) * var(:,:,nk-3) + &
+                           q(5,7) * var(:,:,nk-4) + q(6,7) * var(:,:,nk-5) + &
+                           q(7,7) * var(:,:,nk-6) + q(8,7) * var(:,:,nk-7) + &
+                           q(9,7) * var(:,:,nk-8) + &
+                           q(10,7) * var(:,:,nk-9) ) * idel
+      dvar(:,:,nk-5) = - ( q(1,6) * var(:,:,nk) + q(2,6) * var(:,:,nk-1) + &
+                           q(3,6) * var(:,:,nk-2) + q(4,6) * var(:,:,nk-3) + &
+                           q(5,6) * var(:,:,nk-4) + q(6,6) * var(:,:,nk-5) + &
+                           q(7,6) * var(:,:,nk-6) + q(8,6) * var(:,:,nk-7) + &
+                           q(9,6) * var(:,:,nk-8) + &
+                           q(10,6) * var(:,:,nk-9) ) * idel
+      dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + &
+                           q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + &
+                           q(5,5) * var(:,:,nk-4) + q(6,5) * var(:,:,nk-5) + &
+                           q(7,5) * var(:,:,nk-6) + q(8,5) * var(:,:,nk-7) + &
+                           q(9,5) * var(:,:,nk-8) + &
+                           q(10,5) * var(:,:,nk-9) ) * idel
+      dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + &
+                           q(3,4) * var(:,:,nk-2) + q(4,4) * var(:,:,nk-3) + &
+                           q(5,4) * var(:,:,nk-4) + q(6,4) * var(:,:,nk-5) + &
+                           q(7,4) * var(:,:,nk-6) + q(8,4) * var(:,:,nk-7) + &
+                           q(9,4) * var(:,:,nk-8) + &
+                           q(10,4) * var(:,:,nk-9) ) * idel
+      dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + &
+                           q(3,3) * var(:,:,nk-2) + q(4,3) * var(:,:,nk-3) + &
+                           q(5,3) * var(:,:,nk-4) + q(6,3) * var(:,:,nk-5) + &
+                           q(7,3) * var(:,:,nk-6) + q(8,3) * var(:,:,nk-7) + &
+                           q(9,3) * var(:,:,nk-8) + &
+                           q(10,3) * var(:,:,nk-9) ) * idel
+      dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(2,2) * var(:,:,nk-1) + &
+                           q(3,2) * var(:,:,nk-2) + q(4,2) * var(:,:,nk-3) + &
+                           q(5,2) * var(:,:,nk-4) + q(6,2) * var(:,:,nk-5) + &
+                           q(7,2) * var(:,:,nk-6) + q(8,2) * var(:,:,nk-7) + &
+                           q(9,2) * var(:,:,nk-8) + &
+                           q(10,2) * var(:,:,nk-9) ) * idel
+      dvar(:,:,nk)   = - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + &
+                           q(3,1) * var(:,:,nk-2) + q(4,1) * var(:,:,nk-3) + &
+                           q(5,1) * var(:,:,nk-4) + q(6,1) * var(:,:,nk-5) + &
+                           q(7,1) * var(:,:,nk-6) ) * idel
+      kr = nk - 7
+    end if
+    dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - &
+                                 var(:,:,kl-1:kr-1) ) + &
+                        a(2) * ( var(:,:,kl+2:kr+2) - &
+                                 var(:,:,kl-2:kr-2) ) + &
+                        a(3) * ( var(:,:,kl+3:kr+3) - &
+                                 var(:,:,kl-3:kr-3) ) ) * idel
+  end select direction
+end subroutine deriv_gf_6_5_opt
diff --git a/src/make.code.defn b/src/make.code.defn
index b70f869..882479e 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,28 @@
 # $Header$
 
 # Source files in this directory
-SRCS =  Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 Derivatives_8_4.F90 call_derivs.c Derivatives_4_3.F90 Derivatives_6_5.F90 call_derivs_name.c Get_Coeff.F90 set_norm_mask.F90 CheckGridSizes.F90 dissipation.c Dissipation_6_5.F90 get_coeffs.c Coefficients_2_1.F90 Coefficients_4_2.F90 Coefficients_6_3.F90 Coefficients_8_4.F90 Coefficients_4_3.F90 Coefficients_6_5.F90
+SRCS =  call_derivs.c \
+        call_derivs_name.c \
+        Derivatives_2_1.F90 \
+        Derivatives_4_2.F90 \
+        Derivatives_6_3.F90 \
+        Derivatives_8_4.F90 \
+        Derivatives_4_3.F90 \
+        Derivatives_4_3_min_err_coeff.F90 \
+        Derivatives_6_5.F90 \
+        Derivatives_6_5_min_err_coeff.F90 \
+        Get_Coeff.F90 \
+        set_norm_mask.F90 \
+        CheckGridSizes.F90 \
+        dissipation.c \
+        Dissipation_6_5.F90 \
+        get_coeffs.c \
+        Coefficients_2_1.F90 \
+        Coefficients_4_2.F90 \
+        Coefficients_6_3.F90 \
+        Coefficients_8_4.F90 \
+        Coefficients_4_3.F90 \
+        Coefficients_6_5.F90
 
 # Subdirectories containing source files
 SUBDIRS =