m_a2_ghostcell.f90

module m_a2_ghostcell
  use m_a2_types

  implicit none
  private

  public :: a2_gc_tree
  public :: a2_gc_ids
  public :: a2_gc_box
  public :: a2_bc_dirichlet_zero
  public :: a2_bc_neumann_zero
  public :: a2_bc_continuous
  public :: a2_gc_interp
  public :: a2_gc_prolong_copy
  public :: a2_gc_interp_lim
  public :: a2_gc2_box
  public :: a2_gc2_prolong_linear
  public :: a2_bc2_neumann_zero
  public :: a2_bc2_dirichlet_zero

  ! Define interfaces so ghost cell routines can be called for a single variable
  ! or for multiple variables
  interface a2_gc_tree
     module procedure a2_gc_tree_iv, a2_gc_tree_ivs
  end interface a2_gc_tree

  interface a2_gc_ids
     module procedure a2_gc_ids_iv, a2_gc_ids_ivs, a2_gc_ids_v1
  end interface a2_gc_ids

  interface a2_gc_box
     module procedure a2_gc_box_iv, a2_gc_box_ivs, a2_gc_box_v1
  end interface

contains

  subroutine a2_gc_tree_ivs(tree, ivs, subr_rb, subr_bc, corners, leaves_only)
    type(a2_t), intent(inout)       :: tree
    integer, intent(in)              :: ivs(:)
    procedure(a2_subr_rb), optional :: subr_rb
    procedure(a2_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    logical, intent(in), optional    :: leaves_only
    integer                          :: lvl
    logical                          :: all_ids

    if (.not. tree%ready) error stop "a2_gc_tree: tree not ready"
    all_ids = .true.
    if (present(leaves_only)) all_ids = .not. leaves_only

    if (all_ids) then
       do lvl = lbound(tree%lvls, 1), tree%highest_lvl
          call a2_gc_ids(tree, tree%lvls(lvl)%ids, ivs, &
               subr_rb, subr_bc, corners)
       end do
    else
       do lvl = lbound(tree%lvls, 1), tree%highest_lvl
          call a2_gc_ids(tree, tree%lvls(lvl)%leaves, ivs, &
               subr_rb, subr_bc, corners)
       end do
    end if
  end subroutine a2_gc_tree_ivs

  subroutine a2_gc_tree_iv(tree, iv, subr_rb, subr_bc, corners, leaves_only)
    type(a2_t), intent(inout)       :: tree
    integer, intent(in)              :: iv
    procedure(a2_subr_rb), optional :: subr_rb
    procedure(a2_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    logical, intent(in), optional    :: leaves_only
    integer                          :: lvl
    logical                          :: all_ids

    if (.not. tree%ready) error stop "a2_gc_tree: tree not ready"
    all_ids = .true.
    if (present(leaves_only)) all_ids = .not. leaves_only

    if (all_ids) then
       do lvl = lbound(tree%lvls, 1), tree%highest_lvl
          call a2_gc_ids(tree, tree%lvls(lvl)%ids, iv, &
               subr_rb, subr_bc, corners)
       end do
    else
       do lvl = lbound(tree%lvls, 1), tree%highest_lvl
          call a2_gc_ids(tree, tree%lvls(lvl)%leaves, iv, &
               subr_rb, subr_bc, corners)
       end do
    end if
  end subroutine a2_gc_tree_iv

  subroutine a2_gc_ids_ivs(tree, ids, ivs, subr_rb, subr_bc, corners)
    type(a2_t), intent(inout)       :: tree
    integer, intent(in)              :: ids(:)
    integer, intent(in)              :: ivs(:)
    procedure(a2_subr_rb), optional :: subr_rb
    procedure(a2_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    integer                          :: i

    !$omp parallel do
    do i = 1, size(ids)
       call a2_gc_box(tree, ids(i), ivs, subr_rb, subr_bc, corners)
    end do
    !$omp end parallel do
  end subroutine a2_gc_ids_ivs

  subroutine a2_gc_ids_iv(tree, ids, iv, subr_rb, subr_bc, corners)
    type(a2_t), intent(inout)       :: tree
    integer, intent(in)              :: ids(:)
    integer, intent(in)              :: iv
    procedure(a2_subr_rb), optional :: subr_rb
    procedure(a2_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    integer                          :: i

    !$omp parallel do
    do i = 1, size(ids)
       call a2_gc_box(tree, ids(i), iv, subr_rb, subr_bc, corners)
    end do
    !$omp end parallel do
  end subroutine a2_gc_ids_iv

  subroutine a2_gc_ids_v1(boxes, ids, iv, subr_rb, subr_bc, corners)
    type(box2_t), intent(inout)     :: boxes(:)
    integer, intent(in)              :: ids(:)
    integer, intent(in)              :: iv
    procedure(a2_subr_rb), optional :: subr_rb
    procedure(a2_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    integer                          :: i

    !$omp parallel do
    do i = 1, size(ids)
       call a2_gc_box(boxes, ids(i), iv, subr_rb, subr_bc, corners)
    end do
    !$omp end parallel do
  end subroutine a2_gc_ids_v1

  subroutine a2_gc_box_ivs(tree, id, ivs, subr_rb, subr_bc, corners)
    type(a2_t), intent(inout)       :: tree
    integer, intent(in)              :: id
    integer, intent(in)              :: ivs(:)
    procedure(a2_subr_rb), optional :: subr_rb
    procedure(a2_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    logical                          :: do_corners
    integer                          :: i, iv
    procedure(a2_subr_rb), pointer  :: use_rb
    procedure(a2_subr_bc), pointer  :: use_bc

    do_corners = .true.
    if (present(corners)) do_corners = corners

    do i = 1, size(ivs)
       iv = ivs(i)

       if (present(subr_rb)) then
          use_rb => subr_rb
       else if (tree%has_cc_method(iv)) then
          use_rb => tree%cc_methods(iv)%rb
       else
          error stop "a2_gc_box: no refinement boundary method stored"
       end if

       if (present(subr_bc)) then
          use_bc => subr_bc
       else if (tree%has_cc_method(iv)) then
          use_bc => tree%cc_methods(iv)%bc
       else
          error stop "a2_gc_box: no boundary condition stored"
       end if

       call a2_gc_box_sides(tree%boxes, id, iv, use_rb, use_bc)
       if (do_corners) call a2_gc_box_corner(tree%boxes, id, iv)
    end do
  end subroutine a2_gc_box_ivs

  subroutine a2_gc_box_iv(tree, id, iv, subr_rb, subr_bc, corners)
    type(a2_t), intent(inout)       :: tree
    integer, intent(in)              :: id
    integer, intent(in)              :: iv
    procedure(a2_subr_rb), optional :: subr_rb
    procedure(a2_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    logical                          :: do_corners
    procedure(a2_subr_rb), pointer  :: use_rb
    procedure(a2_subr_bc), pointer  :: use_bc

    if (present(subr_rb)) then
       use_rb => subr_rb
    else if (tree%has_cc_method(iv)) then
       use_rb => tree%cc_methods(iv)%rb
    else
       error stop "a2_gc_box: no refinement boundary method stored"
    end if

    if (present(subr_bc)) then
       use_bc => subr_bc
    else if (tree%has_cc_method(iv)) then
       use_bc => tree%cc_methods(iv)%bc
    else
       error stop "a2_gc_box: no boundary condition stored"
    end if

    do_corners = .true.
    if (present(corners)) do_corners = corners

    call a2_gc_box_sides(tree%boxes, id, iv, use_rb, use_bc)
    if (do_corners) call a2_gc_box_corner(tree%boxes, id, iv)
  end subroutine a2_gc_box_iv

  subroutine a2_gc_box_v1(boxes, id, iv, subr_rb, subr_bc, corners)
    type(box2_t), intent(inout)    :: boxes(:)
    integer, intent(in)             :: id
    integer, intent(in)             :: iv
    procedure(a2_subr_rb)          :: subr_rb
    procedure(a2_subr_bc)          :: subr_bc
    logical, intent(in), optional   :: corners
    logical                         :: do_corners

    do_corners = .true.
    if (present(corners)) do_corners = corners

    call a2_gc_box_sides(boxes, id, iv, subr_rb, subr_bc)
    if (do_corners) call a2_gc_box_corner(boxes, id, iv)
  end subroutine a2_gc_box_v1

  subroutine a2_gc_box_sides(boxes, id, iv, subr_rb, subr_bc)
    type(box2_t), intent(inout) :: boxes(:)
    integer, intent(in)          :: id
    integer, intent(in)          :: iv
    procedure(a2_subr_rb)       :: subr_rb
    procedure(a2_subr_bc)       :: subr_bc
    integer                      :: nb, nb_id, bc_type
    integer                      :: nb_dim, lo(2), hi(2), dnb(2)

    do nb = 1, a2_num_neighbors
       nb_id = boxes(id)%neighbors(nb)

       if (nb_id > af_no_box) then
          ! There is a neighbor
          nb_dim     = a2_neighb_dim(nb)
          lo(:)      = 1
          hi(:)      = boxes(id)%n_cell
          lo(nb_dim) = a2_neighb_high_01(nb) * (boxes(id)%n_cell + 1)
          hi(nb_dim) = lo(nb_dim)
          dnb        = a2_neighb_offset([nb])
          call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, hi, iv)
       else if (nb_id == af_no_box) then
          ! Refinement boundary
          call subr_rb(boxes, id, nb, iv)
       else
          ! Physical boundary
          call subr_bc(boxes(id), nb, iv, bc_type)
          call bc_to_gc(boxes(id), nb, iv, bc_type)
       end if
    end do

  end subroutine a2_gc_box_sides

  subroutine a2_gc_box_corner(boxes, id, iv)
    type(box2_t), intent(inout)  :: boxes(:)
    integer, intent(in)          :: id
    integer, intent(in)          :: iv
    integer                      :: n, nb_id, dnb(2), lo(2)

    do n = 1, a2_num_children
       ! Check whether there is a neighbor, and find its index
       dnb   = 2 * a2_child_dix(:, n) - 1
       nb_id = boxes(id)%neighbor_mat(dnb(1), dnb(2))
       lo    = a2_child_dix(:, n) * (boxes(id)%n_cell + 1)

       if (nb_id > af_no_box) then
          call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, lo, iv)
       else
          call a2_corner_gc_extrap(boxes(id), lo, iv)
       end if
    end do
  end subroutine a2_gc_box_corner

  subroutine bc_to_gc(box, nb, iv, bc_type)
    type(box2_t), intent(inout)  :: box
    integer, intent(in)          :: iv
    integer, intent(in)          :: nb
    integer, intent(in)          :: bc_type
    real(dp)                     :: c0, c1, c2
    integer                      :: nc

    nc = box%n_cell

    ! If we call the interior point x1, x2 and the ghost point x0, then a
    ! Dirichlet boundary value b can be imposed as:
    ! x0 = -x1 + 2*b
    ! A Neumann b.c. can be imposed as:
    ! x0 = x1 +/- dx * b
    ! A continuous boundary (same slope) as:
    ! x0 = 2 * x1 - x2
    ! Below, we set coefficients to handle these cases
    select case (bc_type)
    case (af_bc_dirichlet)
       c0 = 2
       c1 = -1
       c2 = 0
    case (af_bc_neumann)
       c0 = box%dr * a2_neighb_high_pm(nb) ! This gives a + or - sign
       c1 = 1
       c2 = 0
    case (af_bc_continuous)
       c0 = 0
       c1 = 2
       c2 = -1
    case default
       stop "fill_bc: unknown boundary condition"
    end select

    select case (nb)
    case (a2_neighb_lowx)
       box%cc(0, 1:nc, iv) = &
            c0 * box%cc(0, 1:nc, iv) + &
            c1 * box%cc(1, 1:nc, iv) + &
            c2 * box%cc(2, 1:nc, iv)
    case (a2_neighb_highx)
       box%cc(nc+1, 1:nc, iv) = &
            c0 * box%cc(nc+1, 1:nc, iv) + &
            c1 * box%cc(nc, 1:nc, iv) + &
            c2 * box%cc(nc-1, 1:nc, iv)
    case (a2_neighb_lowy)
       box%cc(1:nc, 0, iv) = &
            c0 * box%cc(1:nc, 0, iv) + &
            c1 * box%cc(1:nc, 1, iv) + &
            c2 * box%cc(1:nc, 2, iv)
    case (a2_neighb_highy)
       box%cc(1:nc, nc+1, iv) = &
            c0 * box%cc(1:nc, nc+1, iv) + &
            c1 * box%cc(1:nc, nc, iv) + &
            c2 * box%cc(1:nc, nc-1, iv)
    end select
  end subroutine bc_to_gc

  subroutine a2_gc_prolong_copy(boxes, id, nb, iv)
    use m_a2_prolong, only: a2_prolong_copy
    type(box2_t), intent(inout)  :: boxes(:)
    integer, intent(in)           :: id
    integer, intent(in)           :: iv
    integer, intent(in)           :: nb
    integer                       :: p_id, nb_dim, lo(2), hi(2)

    p_id       = boxes(id)%parent
    nb_dim     = a2_neighb_dim(nb)
    lo(:)      = 1
    hi(:)      = boxes(id)%n_cell
    lo(nb_dim) = a2_neighb_high_01(nb) * (boxes(id)%n_cell+1)
    hi(nb_dim) = a2_neighb_high_01(nb) * (boxes(id)%n_cell+1)

    call a2_prolong_copy(boxes(p_id), boxes(id), iv, low=lo, high=hi)
  end subroutine a2_gc_prolong_copy

  subroutine a2_gc_interp(boxes, id, nb, iv)
    type(box2_t), intent(inout) :: boxes(:)
    integer, intent(in)         :: id
    integer, intent(in)         :: nb
    integer, intent(in)         :: iv
    integer                     :: nc, ix, ix_c, ix_f, i, j
    integer                     :: i_c1, i_c2, j_c1, j_c2, p_nb_id
    integer                     :: p_id, ix_offset(2)
    real(dp), parameter         :: sixth=1/6.0_dp, third=1/3.0_dp

    nc        = boxes(id)%n_cell
    p_id      = boxes(id)%parent
    p_nb_id   = boxes(p_id)%neighbors(nb)
    ix_offset = a2_get_child_offset(boxes(id), nb)

    if (a2_neighb_low(nb)) then
       ix = 0
       ix_f = 1
       ix_c = nc
    else
       ix = nc+1
       ix_f = nc
       ix_c = 1
    end if

    select case (a2_neighb_dim(nb))
    case (1)
       do j = 1, nc
          j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
          j_c2 = j_c1 + 1 - 2 * iand(j, 1)          ! even: +1, odd: -1
          boxes(id)%cc(ix, j, iv) = &
               0.5_dp * boxes(p_nb_id)%cc(ix_c, j_c1, iv) + &
               sixth * boxes(p_nb_id)%cc(ix_c, j_c2, iv) + &
               third * boxes(id)%cc(ix_f, j, iv)
       end do
    case (2)
       do i = 1, nc
          i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
          i_c2 = i_c1 + 1 - 2 * iand(i, 1)          ! even: +1, odd: -1
          boxes(id)%cc(i, ix, iv) = &
               0.5_dp * boxes(p_nb_id)%cc(i_c1, ix_c, iv) + &
               sixth * boxes(p_nb_id)%cc(i_c2, ix_c, iv) + &
               third * boxes(id)%cc(i, ix_f, iv)
       end do
    end select

  end subroutine a2_gc_interp

  subroutine a2_gc_interp_lim(boxes, id, nb, iv)
    type(box2_t), intent(inout) :: boxes(:)
    integer, intent(in)         :: id
    integer, intent(in)         :: nb
    integer, intent(in)         :: iv
    integer                     :: nc, ix, ix_c, ix_f, i, j
    integer                     :: i_c1, i_c2, j_c1, j_c2, p_nb_id
    integer                     :: p_id, ix_offset(2)
    real(dp)                    :: c1, c2
    real(dp), parameter         :: sixth=1/6.0_dp, third=1/3.0_dp

    nc        = boxes(id)%n_cell
    p_id      = boxes(id)%parent
    p_nb_id   = boxes(p_id)%neighbors(nb)
    ix_offset = a2_get_child_offset(boxes(id), nb)

    if (a2_neighb_low(nb)) then
       ix = 0
       ix_f = 1
       ix_c = nc
    else
       ix = nc+1
       ix_f = nc
       ix_c = 1
    end if

    select case (a2_neighb_dim(nb))
    case (1)
       do j = 1, nc
          j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
          j_c2 = j_c1 + 1 - 2 * iand(j, 1)     ! even: +1, odd: -1
          c1 = boxes(p_nb_id)%cc(ix_c, j_c1, iv)
          c2 = boxes(p_nb_id)%cc(ix_c, j_c2, iv)
          boxes(id)%cc(ix, j, iv) = 0.5_dp * c1 + sixth * c2 + &
               third * boxes(id)%cc(ix_f, j, iv)
          if (boxes(id)%cc(ix, j, iv) > 2 * c1) boxes(id)%cc(ix, j, iv) = 2 * c1
       end do
    case (2)
       do i = 1, nc
          i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
          i_c2 = i_c1 + 1 - 2 * iand(i, 1)          ! even: +1, odd: -1
          c1 = boxes(p_nb_id)%cc(i_c1, ix_c, iv)
          c2 = boxes(p_nb_id)%cc(i_c2, ix_c, iv)
          boxes(id)%cc(i, ix, iv) = 0.5_dp * c1 + sixth * c2 + &
               third * boxes(id)%cc(i, ix_f, iv)
          if (boxes(id)%cc(i, ix, iv) > 2 * c1) boxes(id)%cc(i, ix, iv) = 2 * c1
       end do
    end select

  end subroutine a2_gc_interp_lim

  ! This fills ghost cells near physical boundaries using Neumann zero
  subroutine a2_bc_neumann_zero(box, nb, iv, bc_type)
    type(box2_t), intent(inout) :: box
    integer, intent(in)         :: nb, iv
    integer, intent(out)        :: bc_type

    bc_type = af_bc_neumann
    call a2_set_box_gc(box, nb, iv, 0.0_dp)
  end subroutine a2_bc_neumann_zero

  ! This fills ghost cells near physical boundaries using Neumann zero
  subroutine a2_bc_dirichlet_zero(box, nb, iv, bc_type)
    type(box2_t), intent(inout) :: box
    integer, intent(in)         :: nb, iv
    integer, intent(out)        :: bc_type

    bc_type = af_bc_dirichlet
    call a2_set_box_gc(box, nb, iv, 0.0_dp)
  end subroutine a2_bc_dirichlet_zero

  ! This fills ghost cells near physical boundaries using the same slope
  subroutine a2_bc_continuous(box, nb, iv, bc_type)
    type(box2_t), intent(inout) :: box
    integer, intent(in)         :: nb, iv
    integer, intent(out)        :: bc_type

    bc_type = af_bc_continuous
    ! Set values to zero (to prevent problems with NaN)
    call a2_set_box_gc(box, nb, iv, 0.0_dp)
  end subroutine a2_bc_continuous

  subroutine copy_from_nb(box, box_nb, dnb, lo, hi, iv)
    type(box2_t), intent(inout) :: box
    type(box2_t), intent(in)    :: box_nb
    integer, intent(in)          :: dnb(2)
    integer, intent(in)          :: lo(2)
    integer, intent(in)          :: hi(2)
    integer, intent(in)          :: iv
    integer                      :: nlo(2), nhi(2)

    ! Get indices on neighbor
    nlo = lo - dnb * box%n_cell
    nhi = hi - dnb * box%n_cell

    box%cc(lo(1):hi(1), lo(2):hi(2), iv) = &
         box_nb%cc(nlo(1):nhi(1), nlo(2):nhi(2), iv)
  end subroutine copy_from_nb

  subroutine a2_gc2_box(boxes, id, iv, subr_rb, subr_bc, cc, nc)
    type(box2_t), intent(inout) :: boxes(:)
    integer, intent(in)          :: id
    integer, intent(in)          :: iv
    procedure(a2_subr_egc)      :: subr_rb
    procedure(a2_subr_egc)      :: subr_bc
    integer, intent(in)          :: nc
    real(dp), intent(out)        :: cc(-1:nc+2, -1:nc+2)
    real(dp)                     :: gc(1:nc)
    integer                      :: nb, nb_id, nb_dim, lo(2), hi(2)

    cc(0:nc+1, 0:nc+1) = boxes(id)%cc(0:nc+1, 0:nc+1, iv)

    do nb = 1, a2_num_neighbors
       nb_id = boxes(id)%neighbors(nb)

       if (nb_id > af_no_box) then
          call sides2_from_nb(boxes(nb_id), nb, iv, gc, nc)
       else if (nb_id == af_no_box) then
          call subr_rb(boxes, id, nb, iv, gc, nc)
       else
          call subr_bc(boxes, id, nb, iv, gc, nc)
       end if

       ! Determine ghost cell indices
       nb_dim = a2_neighb_dim(nb)
       lo(:) = 1
       hi(:) = boxes(id)%n_cell
       lo(nb_dim) = -1 + a2_neighb_high_01(nb) * (boxes(id)%n_cell + 3)
       hi(nb_dim) = lo(nb_dim)

       cc(lo(1):hi(1), lo(2):hi(2)) = reshape(gc, 1 + hi - lo)
    end do
  end subroutine a2_gc2_box

  subroutine sides2_from_nb(box_nb, nb, iv, gc_side, nc)
    type(box2_t), intent(in) :: box_nb
    integer, intent(in)       :: nb
    integer, intent(in)       :: iv
    integer, intent(in)       :: nc
    real(dp), intent(out)     :: gc_side(nc)

    select case (nb)
    case (a2_neighb_lowx)
       gc_side = box_nb%cc(nc-1, 1:nc, iv)
    case (a2_neighb_highx)
       gc_side = box_nb%cc(2, 1:nc, iv)
    case (a2_neighb_lowy)
       gc_side = box_nb%cc(1:nc, nc-1, iv)
    case (a2_neighb_highy)
       gc_side = box_nb%cc(1:nc, 2, iv)
    end select
  end subroutine sides2_from_nb

  subroutine a2_gc2_prolong_linear(boxes, id, nb, iv, gc_side, nc)
    type(box2_t), intent(inout) :: boxes(:)
    integer, intent(in)         :: id
    integer, intent(in)         :: nb
    integer, intent(in)         :: iv
    integer, intent(in)         :: nc
    real(dp), intent(out)       :: gc_side(nc)
    integer                     :: ix, i, j
    integer                     :: i_c1, i_c2, j_c1, j_c2, p_nb_id
    integer                     :: p_id, ix_offset(2)

    p_id      = boxes(id)%parent
    p_nb_id   = boxes(p_id)%neighbors(nb)
    ix_offset = a2_get_child_offset(boxes(id), nb)
    ix        = a2_neighb_high_01(nb) * (nc+3) - 1 ! -1 or nc+2

    select case (a2_neighb_dim(nb))
    case (1)
       i_c1 = ix_offset(1) + ishft(ix+1, -1) ! (ix+1)/2
       i_c2 = i_c1 + 1 - 2 * iand(ix, 1)     ! even: +1, odd: -1
       do j = 1, nc
          j_c1 = ix_offset(2) + ishft(j+1, -1) ! (j+1)/2
          j_c2 = j_c1 + 1 - 2 * iand(j, 1)     ! even: +1, odd: -1
          gc_side(j) = &
               0.5_dp * boxes(p_nb_id)%cc(i_c1, j_c1, iv) + &
               0.25_dp * boxes(p_nb_id)%cc(i_c2, j_c1, iv) + &
               0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c2, iv)
       end do
    case (2)
       j_c1 = ix_offset(2) + ishft(ix+1, -1) ! (j+1)/2
       j_c2 = j_c1 + 1 - 2 * iand(ix, 1)     ! even: +1, odd: -1
       do i = 1, nc
          i_c1 = ix_offset(1) + ishft(i+1, -1) ! (i+1)/2
          i_c2 = i_c1 + 1 - 2 * iand(i, 1)     ! even: +1, odd: -1
          gc_side(i) = &
               0.5_dp * boxes(p_nb_id)%cc(i_c1, j_c1, iv) + &
               0.25_dp * boxes(p_nb_id)%cc(i_c2, j_c1, iv) + &
               0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c2, iv)
       end do
    end select

  end subroutine a2_gc2_prolong_linear

  ! This fills second ghost cells near physical boundaries using Neumann zero
  subroutine a2_bc2_neumann_zero(boxes, id, nb, iv, gc_side, nc)
    type(box2_t), intent(inout) :: boxes(:)
    integer, intent(in)         :: id
    integer, intent(in)         :: nb
    integer, intent(in)         :: iv
    integer, intent(in)         :: nc
    real(dp), intent(out)       :: gc_side(nc)

    select case (nb)
    case (a2_neighb_lowx)
       gc_side = boxes(id)%cc(2, 1:nc, iv)
    case (a2_neighb_highx)
       gc_side = boxes(id)%cc(nc-1, 1:nc, iv)
    case (a2_neighb_lowy)
       gc_side = boxes(id)%cc(1:nc, 2, iv)
    case (a2_neighb_highy)
       gc_side = boxes(id)%cc(1:nc, nc-1, iv)
    end select
  end subroutine a2_bc2_neumann_zero

  ! This fills second ghost cells near physical boundaries using Neumann zero
  subroutine a2_bc2_dirichlet_zero(boxes, id, nb, iv, gc_side, nc)
    type(box2_t), intent(inout) :: boxes(:)
    integer, intent(in)         :: id, nb, iv, nc
    real(dp), intent(out)       :: gc_side(nc)

    select case (nb)
    case (a2_neighb_lowx)
       gc_side = -boxes(id)%cc(2, 1:nc, iv)
    case (a2_neighb_highx)
       gc_side = -boxes(id)%cc(nc-1, 1:nc, iv)
    case (a2_neighb_lowy)
       gc_side = -boxes(id)%cc(1:nc, 2, iv)
    case (a2_neighb_highy)
       gc_side = -boxes(id)%cc(1:nc, nc-1, iv)
    end select
  end subroutine a2_bc2_dirichlet_zero

  subroutine a2_corner_gc_extrap(box, ix, iv)
    type(box2_t), intent(inout) :: box
    integer, intent(in)          :: ix(2)
    integer, intent(in)          :: iv
    integer                      :: di(2)

    di = 1 - 2 * iand(ix, 1)    ! 0 -> di = 1, nc+1 -> di = -1

    box%cc(ix(1), ix(2), iv) = box%cc(ix(1)+di(1), ix(2), iv) &
         + box%cc(ix(1), ix(2)+di(2), iv) &
         - box%cc(ix(1)+di(1), ix(2)+di(2), iv)
  end subroutine a2_corner_gc_extrap

end module m_a2_ghostcell