#include "cctk.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"

subroutine deriv_gf_2_1 ( 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(1), save :: a
  CCTK_REAL, dimension(3,2), save :: q 
  CCTK_REAL :: idel

  CCTK_INT :: il, ir, jl, jr, kl, kr

  logical, save :: first = .true.

  if ( first ) then
    a(1) = 1.0_wp/2.0_wp

    q(1,1) = -1.0_wp; q(2,1) = 1.0_wp; q(3,1) = zero
    q(1,2) = -1.0_wp/2.0_wp; q(2,2) = zero; q(3,2) = 1.0_wp/2.0_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,:,:) ) * idel
      dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(3,2) * var(3,:,:) ) * idel
      il = 3
    end if
    if ( bb(2) == 0 ) then
      ir = ni - gsize
    else
      dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + &
                           q(3,2) * var(ni-2,:,:) ) * idel
      dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + &
                         q(2,1) * var(ni-1,:,:) ) * idel
      ir = ni - 2
    end if
    dvar(il:ir,:,:) = a(1) * ( var(il+1:ir+1,:,:) - var(il-1:ir-1,:,:) ) * 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,:) ) * idel
      dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(3,2) * var(:,3,:) ) * idel
      jl = 3
    end if
    if ( bb(2) == 0 ) then
      jr = nj - gsize
    else
      dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + &
                           q(3,2) * var(:,nj-2,:) ) * idel
      dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + &
                         q(2,1) * var(:,nj-1,:) ) * idel
      jr = nj - 2
    end if
    dvar(:,jl:jr,:) = a(1) * ( var(:,jl+1:jr+1,:) - var(:,jl-1:jr-1,:) ) * 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) ) * idel
      dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(3,2) * var(:,:,3) ) * idel
      kl = 3
    end if
    if ( bb(2) == 0 ) then
      kr = nk - gsize
    else 
      dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + &
                           q(3,2) * var(:,:,nk-2) ) * idel
      dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + &
                         q(2,1) * var(:,:,nk-1) ) * idel
      kr = nk - 2
    end if
    dvar(:,:,kl:kr) = a(1) * ( var(:,:,kl+1:kr+1) - var(:,:,kl-1:kr-1) ) * idel
  end select direction
end subroutine deriv_gf_2_1