I forgot the actual finite difference routine the first time.

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

1 files changed, 378 insertions, 0 deletions

diff --git a/src/Derivatives_6_5.F90 b/src/Derivatives_6_5.F90
new file mode 100644
index 0000000..60b55a2
--- /dev/null
+++ b/src/Derivatives_6_5.F90
@@ -0,0 +1,378 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine deriv_gf_6_5 ( 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
+!    print*,'6-5 called'
+    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) = -137.0_wp/60.0_wp; q(2,1) = 5.0_wp; q(3,1) = -5.0_wp;
+    q(4,1) = 10.0_wp/3.0_wp; q(5,1) = -5.0_wp/4.0_wp; q(6,1) = 1.0_wp/5.0_wp;
+    q(7,1) = zero; q(8,1) = zero; q(9,1) = zero; q(10,1) = zero;
+
+    q(1,2) = -0.3209098489859668118347_wp
+    q(2,2) = -0.359441670715467075894_wp
+    q(3,2) = 0.195756852998105503889_wp
+    q(4,2) = 1.394685510250317033169_wp
+    q(5,2) = -1.448965775497476572821_wp
+    q(6,2) = 0.6519476244467816674835_wp
+    q(7,2) = -0.1115052611983591304249_wp
+    q(8,2) = -0.00156743129793461356831_wp
+    q(9,2) = zero; q(10,2) = zero
+
+    q(1,3) = 0.0980796259621024870281_wp
+    q(2,3) = -0.784752101846956387607_wp
+    q(3,3) = 0.364206221878666751300_wp
+    q(4,3) = 0.102097753636344358379_wp
+    q(5,3) = 0.377167650934576412896_wp
+    q(6,3) = -0.173241400242683302020_wp
+    q(7,3) = 0.0062120424243610783650_wp
+    q(8,3) = 0.01153111791917461507898_wp
+    q(9,3) = -0.00130091066558601341956_wp
+    q(10,3) = zero
+
+    q(1,4) = 0.0771947815148234674634_wp
+    q(2,4) = -0.4124719432622107102831_wp
+    q(3,4) = 0.654798717867622763906_wp
+    q(4,4) = -1.873474657556531216053_wp
+    q(5,4) = 2.161961221935152129759_wp
+    q(6,4) = -0.729147814207596325277_wp
+    q(7,4) = 0.1292509053479867634674_wp
+    q(8,4) = -0.01092898521375837831427_wp
+    q(9,4) = 0.00316983426828354261009_wp
+    q(10,4) = -0.000352060693772037278516_wp
+
+    q(1,5) = -0.01135836398682271847110_wp
+    q(2,5) = 0.0114974985688122284189_wp
+    q(3,5) = 0.228765565489867047381_wp
+    q(4,5) = -1.179122428698534075609_wp
+    q(5,5) = 0.774618306056840861215_wp
+    q(6,5) = 0.016612403904600041087_wp
+    q(7,5) = 0.2399304273448163041220_wp
+    q(8,5) = -0.0959180672407791626759_wp
+    q(9,5) = 0.01612460763754670343830_wp
+    q(10,5) = -0.001149949076347228906108_wp
+
+    q(1,6) = -0.01912046962023373131259_wp
+    q(2,6) = 0.1569245474161812198493_wp
+    q(3,6) = -0.567250839387340928718_wp
+    q(4,6) = 1.230411310316952671353_wp
+    q(5,6) = -2.048795717924254691221_wp
+    q(6,6) = 0.982715833107347313532_wp
+    q(7,6) = 0.2876584784911885215053_wp
+    q(8,6) = -0.0207657038198929292487_wp
+    q(9,6) = -0.00335290957957120811713_wp
+    q(10,6) = 0.001575470999623762377395_wp
+
+    q(1,7) = -0.0076072322555736921424_wp 
+    q(2,7) = 0.034503414345314189270_wp
+    q(3,7) = -0.044377233426986476366_wp
+    q(4,7) = -0.049408595803628143797_wp
+    q(5,7) = 0.304693537653175477283_wp
+    q(6,7) = -0.935863321236636595622_wp
+    q(7,7) = 0.111375967229143365899_wp 
+    q(8,7) = 0.7149322477536284813995_wp
+    q(9,7) = -0.1444768161656941519630_wp
+    q(10,7) = 0.01622803190725754603783_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,:,:) ) * 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,:,:) ) * 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,:,:) ) * 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
+!      print*,'x-'
+    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,:,:)) * 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,:,:) ) * 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,:,:) ) * idel
+      ir = ni - 7
+!      print*,'x+'
+    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,:) ) * 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,:) ) * 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,:) ) * 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
+!      print*,'y-'
+      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,:)) * 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,:) ) * 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,:) ) * idel
+!      print*,'y+'
+      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) ) * 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) ) * 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) ) * 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
+!      print*,'z-'
+      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)) * 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) ) * 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) ) * idel
+!      print*,'z+'
+      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
+!  print*,'Exiting'
+end subroutine deriv_gf_6_5