#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
    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

  if (gsize < 3) call CCTK_WARN (0, "not enough ghostzones")

  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
    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
    end if
    if (il > ir+1) call CCTK_WARN (0, "domain too small")
    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 ( zero_derivs_y /= 0 ) then
      dvar = zero
    else
      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
        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
        jr = nj - 7
      end if
      if (jl > jr+1) call CCTK_WARN (0, "domain too small")
      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
    end if
  case (2) direction
    if ( zero_derivs_z /= 0 ) then
      dvar = zero
    else
      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
        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
        kr = nk - 7
      end if
      if (kl > kr+1) call CCTK_WARN (0, "domain too small")
      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 if
  end select direction
end subroutine deriv_gf_6_5