



























































































































 


































































































































































 
































































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                

















































































subroutine deriv2_gf_2_1 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar2 )

  use All_Coeffs_mod

  implicit none

  interface
  integer*8 function CCTK_PointerTo (var)
  implicit none
  type(*),dimension(..),TARGET :: var
  end function CCTK_PointerTo
  end interface
  interface
  integer function CCTK_Equals (arg1, arg2)
  implicit none
  integer*8 arg1
  character(*) arg2
  end function CCTK_Equals
  integer function CCTK_MyProc (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_MyProc
  integer function CCTK_nProcs (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_nProcs
  integer function CCTK_IsThornActive (name)
  implicit none
  character(*) name
  end function CCTK_IsThornActive
  integer*8 function CCTK_NullPointer ()
  implicit none
  end function CCTK_NullPointer
  end interface
  interface
  subroutine Diff2_coeff (cctkGH, dir, nsize, imin, imax, q, table_handle)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir
  INTEGER*4 nsize
  INTEGER*4 imin(*)
  INTEGER*4 imax(*)
  REAL*8 q(*)
  INTEGER*4 table_handle
  end subroutine Diff2_coeff
  end interface
  interface
  subroutine Diff2_gv (cctkGH, dir1, dir2, var, dvar, table_handle)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir1
  INTEGER*4 dir2
  REAL*8 var(*)
  REAL*8 dvar(*)
  INTEGER*4 table_handle
  end subroutine Diff2_gv
  end interface
  interface
  subroutine Diff_coeff (cctkGH, dir, nsize, imin, imax, q, table_handle)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir
  INTEGER*4 nsize
  INTEGER*4 imin(*)
  INTEGER*4 imax(*)
  REAL*8 q(*)
  INTEGER*4 table_handle
  end subroutine Diff_coeff
  end interface
  interface
  subroutine Diff_gf (cctkGH, dir, var_name, dvar_name)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir
  character(*) var_name
  character(*) dvar_name
  end subroutine Diff_gf
  end interface
  interface
  subroutine Diff_gv (cctkGH, dir, var, dvar, table_handle)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir
  REAL*8 var(*)
  REAL*8 dvar(*)
  INTEGER*4 table_handle
  end subroutine Diff_gv
  end interface
  interface
  subroutine Diff_up_coeff (cctkGH, dir, nsize, imin, imax, q, up, table_handle)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir
  INTEGER*4 nsize
  INTEGER*4 imin(*)
  INTEGER*4 imax(*)
  REAL*8 q(*)
  INTEGER*4 up
  INTEGER*4 table_handle
  end subroutine Diff_up_coeff
  end interface
  interface
  subroutine Diff_up_gv (cctkGH, dir, var, dvar, up, table_handle)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir
  REAL*8 var(*)
  REAL*8 dvar(*)
  REAL*8 up(*)
  INTEGER*4 table_handle
  end subroutine Diff_up_gv
  end interface
  interface
  subroutine GetBoundWidth (cctkGH, bsize, table_handle)
  implicit none
  integer*8 cctkGH
  INTEGER*4 bsize(*)
  INTEGER*4 table_handle
  end subroutine GetBoundWidth
  end interface
  interface
  INTEGER*4 function GetBoundarySpecification (size, nboundaryzones, is_internal, is_staggered, shiftout)
  implicit none
  INTEGER*4 size
  INTEGER*4 nboundaryzones(*)
  INTEGER*4 is_internal(*)
  INTEGER*4 is_staggered(*)
  INTEGER*4 shiftout(*)
  end function GetBoundarySpecification
  end interface
  interface
  INTEGER*4 function GetDomainSpecification (size, physical_min, physical_max, interior_min, interior_max, exterior_min, exterior_m&
  &ax, spacing)
  implicit none
  INTEGER*4 size
  REAL*8 physical_min(*)
  REAL*8 physical_max(*)
  REAL*8 interior_min(*)
  REAL*8 interior_max(*)
  REAL*8 exterior_min(*)
  REAL*8 exterior_max(*)
  REAL*8 spacing(*)
  end function GetDomainSpecification
  end interface
  interface
  subroutine GetLshIndexRanges (cctkGH, imin, imax)
  implicit none
  integer*8 cctkGH
  INTEGER*4 imin(*)
  INTEGER*4 imax(*)
  end subroutine GetLshIndexRanges
  end interface
  interface
  REAL*8 function GetScalProdCoeff ()
  implicit none
  end function GetScalProdCoeff
  end interface
  interface
  subroutine GetScalProdDiag (cctkGH, dir, nsize, sigmad)
  implicit none
  integer*8 cctkGH
  INTEGER*4 dir
  INTEGER*4 nsize
  REAL*8 sigmad(*)
  end subroutine GetScalProdDiag
  end interface
  interface
  INTEGER*4 function MoLQueryEvolvedRHS (EvolvedIndex)
  implicit none
  INTEGER*4 EvolvedIndex
  end function MoLQueryEvolvedRHS
  end interface
  interface
  INTEGER*4 function MultiPatch_GetBbox (cctkGH, size, bbox)
  implicit none
  integer*8 cctkGH
  INTEGER*4 size
  INTEGER*4 bbox(*)
  end function MultiPatch_GetBbox
  end interface
  interface
  INTEGER*4 function MultiPatch_GetBoundarySpecification (map, size, nboundaryzones, is_internal, is_staggered, shiftout)
  implicit none
  INTEGER*4 map
  INTEGER*4 size
  INTEGER*4 nboundaryzones(*)
  INTEGER*4 is_internal(*)
  INTEGER*4 is_staggered(*)
  INTEGER*4 shiftout(*)
  end function MultiPatch_GetBoundarySpecification
  end interface
  interface
  INTEGER*4 function MultiPatch_GetMap (cctkGH)
  implicit none
  integer*8 cctkGH
  end function MultiPatch_GetMap
  end interface
  interface
  INTEGER*4 function MultiPatch_GetMaps (cctkGH)
  implicit none
  integer*8 cctkGH
  end function MultiPatch_GetMaps
  end interface
  interface
  INTEGER*4 function SymmetryHandleOfName (sym_name)
  implicit none
  character(*) sym_name
  end function SymmetryHandleOfName
  end interface
  interface
  INTEGER*4 function SymmetryTableHandleForGrid (cctkGH)
  implicit none
  integer*8 cctkGH
  end function SymmetryTableHandleForGrid
  end interface
  REAL*8 diss_fraction (3)
  integer, parameter :: cctki_use_diss_fraction = kind(diss_fraction)
  REAL*8 epsdis
  integer, parameter :: cctki_use_epsdis = kind(epsdis)
  REAL*8 h_scaling (3)
  integer, parameter :: cctki_use_h_scaling = kind(h_scaling)
  REAL*8 poison_value
  integer, parameter :: cctki_use_poison_value = kind(poison_value)
  integer*8 dissipation_type
  integer, parameter :: cctki_use_dissipation_type = kind(dissipation_type)
  integer*8 norm_type
  integer, parameter :: cctki_use_norm_type = kind(norm_type)
  integer*8 operator_type
  integer, parameter :: cctki_use_operator_type = kind(operator_type)
  integer*8 vars
  integer, parameter :: cctki_use_vars = kind(vars)
  INTEGER*4 check_grid_sizes
  integer, parameter :: cctki_use_check_grid_sizes = kind(check_grid_sizes)
  INTEGER*4 onesided_interpatch_boundaries
  integer, parameter :: cctki_use_onesided_interpatch_boundaries = kind(onesided_interpatch_boundaries)
  INTEGER*4 onesided_outer_boundaries
  integer, parameter :: cctki_use_onesided_outer_boundaries = kind(onesided_outer_boundaries)
  INTEGER*4 order
  integer, parameter :: cctki_use_order = kind(order)
  INTEGER*4 poison_derivatives
  integer, parameter :: cctki_use_poison_derivatives = kind(poison_derivatives)
  INTEGER*4 poison_dissipation
  integer, parameter :: cctki_use_poison_dissipation = kind(poison_dissipation)
  INTEGER*4 sbp_1st_deriv
  integer, parameter :: cctki_use_sbp_1st_deriv = kind(sbp_1st_deriv)
  INTEGER*4 sbp_2nd_deriv
  integer, parameter :: cctki_use_sbp_2nd_deriv = kind(sbp_2nd_deriv)
  INTEGER*4 sbp_upwind_deriv
  integer, parameter :: cctki_use_sbp_upwind_deriv = kind(sbp_upwind_deriv)
  INTEGER*4 scale_with_h
  integer, parameter :: cctki_use_scale_with_h = kind(scale_with_h)
  INTEGER*4 use_dissipation
  integer, parameter :: cctki_use_use_dissipation = kind(use_dissipation)
  INTEGER*4 use_shiftout
  integer, parameter :: cctki_use_use_shiftout = kind(use_shiftout)
  INTEGER*4 use_variable_deltas
  integer, parameter :: cctki_use_use_variable_deltas = kind(use_variable_deltas)
  INTEGER*4 zero_derivs_y
  integer, parameter :: cctki_use_zero_derivs_y = kind(zero_derivs_y)
  INTEGER*4 zero_derivs_z
  integer, parameter :: cctki_use_zero_derivs_z = kind(zero_derivs_z)
  COMMON /SummationByPartsrest/diss_fraction, epsdis, h_scaling, poison_value, dissipation_type, norm_type, operator_type, vars, ch&
  &eck_grid_sizes, onesided_interpatch_boundaries, onesided_outer_boundaries, order, poison_derivatives, poison_dissipation, sbp_1s&
  &t_deriv, sbp_2nd_deriv, sbp_upwind_deriv, scale_with_h, use_dissipation, use_shiftout, use_variable_deltas, zero_derivs_y, zero_&
  &derivs_z

  INTEGER*4, intent(IN) :: ni, nj, nk
  REAL*8, dimension(ni,nj,nk), intent(IN) :: var
  INTEGER*4, intent(IN) :: dir
  INTEGER*4, intent(IN) :: bb(2)
  INTEGER*4, intent(IN) :: gsize
  REAL*8, intent(IN) :: delta
  REAL*8, dimension(ni,nj,nk), intent(OUT) :: dvar2

  REAL*8, dimension(2), save :: a
  REAL*8, dimension(3,2), save :: q 
  REAL*8 :: idel2

  INTEGER*4 :: il, ir, jl, jr, kl, kr

  logical, save :: first = .true.

  if ( first ) then
    call coeffs_2_2_1 ( a, q )
    first = .false.
  end if

  idel2 = 1.0_wp / delta**2

  direction: select case (dir)
  case (0) direction
    if ( bb(1) == 0 ) then
      il = 1 + gsize
    else
      dvar2(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + &
                       q(3,1) * var(3,:,:) ) * idel2
      dvar2(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + &
                       q(3,2) * var(3,:,:) ) * idel2
      il = 3
    end if
    if ( bb(2) == 0 ) then
      ir = ni - gsize
    else
      dvar2(ni-1,:,:) = ( q(1,2) * var(ni,:,:) + &
                          q(2,2) * var(ni-1,:,:) + &
                          q(3,2) * var(ni-2,:,:) ) * idel2
      dvar2(ni,:,:) = ( q(1,1) * var(ni,:,:) + &
                        q(2,1) * var(ni-1,:,:) + &
                        q(3,1) * var(ni-2,:,:) ) * idel2
      ir = ni - 2
    end if
    dvar2(il:ir,:,:) = ( a(2) * ( var(il+1:ir+1,:,:) + &
                                  var(il-1:ir-1,:,:) ) + &
                         a(1) * var(il:ir,:,:) ) * idel2
  case (1) direction
    if ( zero_derivs_y /= 0 ) then
      dvar2 = zero
    else
      if ( bb(1) == 0 ) then
        jl = 1 + gsize
      else
        dvar2(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + &
                         q(3,1) * var(:,3,:) ) * idel2
        dvar2(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + &
                         q(3,2) * var(:,3,:) ) * idel2
        jl = 3
      end if
      if ( bb(2) == 0 ) then
        jr = nj - gsize
      else
        dvar2(:,nj-1,:) = ( q(1,2) * var(:,nj,:) + &
                            q(2,2) * var(:,nj-1,:) + &
                            q(3,2) * var(:,nj-2,:) ) * idel2
        dvar2(:,nj,:) = ( q(1,1) * var(:,nj,:) + &
                          q(2,1) * var(:,nj-1,:) + &
                          q(3,1) * var(:,nj-2,:) ) * idel2
        jr = nj - 2
      end if
      dvar2(:,jl:jr,:) = ( a(2) * ( var(:,jl+1:jr+1,:) + &
                                    var(:,jl-1:jr-1,:) ) + &
                           a(1) * var(:,jl:jr,:) ) * idel2
    end if
  case (2) direction
    if ( zero_derivs_z /= 0 ) then
      dvar2 = zero
    else
      if ( bb(1) == 0 ) then
        kl = 1 + gsize
      else
        dvar2(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + &
                         q(3,1) * var(:,:,3) ) * idel2
        dvar2(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + &
                         q(3,2) * var(:,:,3) ) * idel2
        kl = 3
      end if
      if ( bb(2) == 0 ) then
        kr = nk - gsize
      else 
        dvar2(:,:,nk-1) = ( q(1,2) * var(:,:,nk) + &
                            q(2,2) * var(:,:,nk-1) + &
                            q(3,2) * var(:,:,nk-2) ) * idel2
        dvar2(:,:,nk) = ( q(1,1) * var(:,:,nk) + &
                          q(2,1) * var(:,:,nk-1) + &
                          q(3,1) * var(:,:,nk-2) ) * idel2
        kr = nk - 2
      end if
      dvar2(:,:,kl:kr) = ( a(2) * ( var(:,:,kl+1:kr+1) +  &
                                    var(:,:,kl-1:kr-1) ) + &
                           a(1) * var(:,:,kl:kr) ) * idel2
    end if
  end select direction
end subroutine deriv2_gf_2_1
