m_a3_ghostcell.f90

module m_a3_ghostcell
  use m_a3_types

  implicit none
  private

  public :: a3_gc_tree
  public :: a3_gc_ids
  public :: a3_gc_box
  public :: a3_bc_dirichlet_zero
  public :: a3_bc_neumann_zero
  public :: a3_bc_continuous
  public :: a3_gc_interp
  public :: a3_gc_prolong_copy
  public :: a3_gc_interp_lim
  public :: a3_gc2_box
  public :: a3_gc2_prolong_linear
  public :: a3_bc2_neumann_zero
  public :: a3_bc2_dirichlet_zero

  ! Define interfaces so ghost cell routines can be called for a single variable
  ! or for multiple variables
  interface a3_gc_tree
     module procedure a3_gc_tree_iv, a3_gc_tree_ivs
  end interface a3_gc_tree

  interface a3_gc_ids
     module procedure a3_gc_ids_iv, a3_gc_ids_ivs, a3_gc_ids_v1
  end interface a3_gc_ids

  interface a3_gc_box
     module procedure a3_gc_box_iv, a3_gc_box_ivs, a3_gc_box_v1
  end interface

contains

  subroutine a3_gc_tree_ivs(tree, ivs, subr_rb, subr_bc, corners, leaves_only)
    type(a3_t), intent(inout)       :: tree
    integer, intent(in)              :: ivs(:)
    procedure(a3_subr_rb), optional :: subr_rb
    procedure(a3_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 "a3_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 a3_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 a3_gc_ids(tree, tree%lvls(lvl)%leaves, ivs, &
               subr_rb, subr_bc, corners)
       end do
    end if
  end subroutine a3_gc_tree_ivs

  subroutine a3_gc_tree_iv(tree, iv, subr_rb, subr_bc, corners, leaves_only)
    type(a3_t), intent(inout)       :: tree
    integer, intent(in)              :: iv
    procedure(a3_subr_rb), optional :: subr_rb
    procedure(a3_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 "a3_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 a3_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 a3_gc_ids(tree, tree%lvls(lvl)%leaves, iv, &
               subr_rb, subr_bc, corners)
       end do
    end if
  end subroutine a3_gc_tree_iv

  subroutine a3_gc_ids_ivs(tree, ids, ivs, subr_rb, subr_bc, corners)
    type(a3_t), intent(inout)       :: tree
    integer, intent(in)              :: ids(:)
    integer, intent(in)              :: ivs(:)
    procedure(a3_subr_rb), optional :: subr_rb
    procedure(a3_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    integer                          :: i

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

  subroutine a3_gc_ids_iv(tree, ids, iv, subr_rb, subr_bc, corners)
    type(a3_t), intent(inout)       :: tree
    integer, intent(in)              :: ids(:)
    integer, intent(in)              :: iv
    procedure(a3_subr_rb), optional :: subr_rb
    procedure(a3_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    integer                          :: i

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

  subroutine a3_gc_ids_v1(boxes, ids, iv, subr_rb, subr_bc, corners)
    type(box3_t), intent(inout)     :: boxes(:)
    integer, intent(in)              :: ids(:)
    integer, intent(in)              :: iv
    procedure(a3_subr_rb), optional :: subr_rb
    procedure(a3_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    integer                          :: i

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

  subroutine a3_gc_box_ivs(tree, id, ivs, subr_rb, subr_bc, corners)
    type(a3_t), intent(inout)       :: tree
    integer, intent(in)              :: id
    integer, intent(in)              :: ivs(:)
    procedure(a3_subr_rb), optional :: subr_rb
    procedure(a3_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    logical                          :: do_corners
    integer                          :: i, iv
    procedure(a3_subr_rb), pointer  :: use_rb
    procedure(a3_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 "a3_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 "a3_gc_box: no boundary condition stored"
       end if

       call a3_gc_box_sides(tree%boxes, id, iv, use_rb, use_bc)
       if (do_corners) call a3_gc_box_corner(tree%boxes, id, iv)
    end do
  end subroutine a3_gc_box_ivs

  subroutine a3_gc_box_iv(tree, id, iv, subr_rb, subr_bc, corners)
    type(a3_t), intent(inout)       :: tree
    integer, intent(in)              :: id
    integer, intent(in)              :: iv
    procedure(a3_subr_rb), optional :: subr_rb
    procedure(a3_subr_bc), optional :: subr_bc
    logical, intent(in), optional    :: corners
    logical                          :: do_corners
    procedure(a3_subr_rb), pointer  :: use_rb
    procedure(a3_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 "a3_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 "a3_gc_box: no boundary condition stored"
    end if

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

    call a3_gc_box_sides(tree%boxes, id, iv, use_rb, use_bc)
    if (do_corners) call a3_gc_box_corner(tree%boxes, id, iv)
  end subroutine a3_gc_box_iv

  subroutine a3_gc_box_v1(boxes, id, iv, subr_rb, subr_bc, corners)
    type(box3_t), intent(inout)    :: boxes(:)
    integer, intent(in)             :: id
    integer, intent(in)             :: iv
    procedure(a3_subr_rb)          :: subr_rb
    procedure(a3_subr_bc)          :: subr_bc
    logical, intent(in), optional   :: corners
    logical                         :: do_corners

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

    call a3_gc_box_sides(boxes, id, iv, subr_rb, subr_bc)
    if (do_corners) call a3_gc_box_corner(boxes, id, iv)
  end subroutine a3_gc_box_v1

  subroutine a3_gc_box_sides(boxes, id, iv, subr_rb, subr_bc)
    type(box3_t), intent(inout) :: boxes(:)
    integer, intent(in)          :: id
    integer, intent(in)          :: iv
    procedure(a3_subr_rb)       :: subr_rb
    procedure(a3_subr_bc)       :: subr_bc
    integer                      :: nb, nb_id, bc_type
    integer                      :: nb_dim, lo(3), hi(3), dnb(3)

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

       if (nb_id > af_no_box) then
          ! There is a neighbor
          nb_dim     = a3_neighb_dim(nb)
          lo(:)      = 1
          hi(:)      = boxes(id)%n_cell
          lo(nb_dim) = a3_neighb_high_01(nb) * (boxes(id)%n_cell + 1)
          hi(nb_dim) = lo(nb_dim)
          dnb        = a3_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 a3_gc_box_sides

  subroutine a3_gc_box_corner(boxes, id, iv)
    type(box3_t), intent(inout)  :: boxes(:)
    integer, intent(in)          :: id
    integer, intent(in)          :: iv
    integer                      :: n, nb_id, dnb(3), lo(3)
    integer                      :: hi(3), dim

    ! Have to do edges before corners (since extrapolation can depend on edge values)
    do n = 1, a3_num_edges
       dim = a3_edge_dim(n)

       ! Check whether there is a neighbor, and find its index
       nb_id = boxes(id)%neighbor_mat(a3_edge_dir(1, n), &
            a3_edge_dir(2, n), a3_edge_dir(3, n))

       lo = a3_edge_min_ix(:, n) * (boxes(id)%n_cell + 1)
       lo(dim) = 1

       if (nb_id > af_no_box) then
          hi      = lo
          hi(dim) = boxes(id)%n_cell
          dnb   = a3_neighb_offset(a3_nb_adj_edge(:, n))
          call copy_from_nb(boxes(id), boxes(nb_id), dnb, lo, hi, iv)
       else
          call a3_edge_gc_extrap(boxes(id), lo, dim, iv)
       end if
    end do

    do n = 1, a3_num_children
       ! Check whether there is a neighbor, and find its index
       dnb   = 2 * a3_child_dix(:, n) - 1
       nb_id = boxes(id)%neighbor_mat(dnb(1), dnb(2), dnb(3))
       lo    = a3_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 a3_corner_gc_extrap(boxes(id), lo, iv)
       end if
    end do
  end subroutine a3_gc_box_corner

  subroutine bc_to_gc(box, nb, iv, bc_type)
    type(box3_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 * a3_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 (a3_neighb_lowx)
       box%cc(0, 1:nc, 1:nc, iv) = &
            c0 * box%cc(0, 1:nc, 1:nc, iv) + &
            c1 * box%cc(1, 1:nc, 1:nc, iv) + &
            c2 * box%cc(2, 1:nc, 1:nc, iv)
    case (a3_neighb_highx)
       box%cc(nc+1, 1:nc, 1:nc, iv) = &
            c0 * box%cc(nc+1, 1:nc, 1:nc, iv) + &
            c1 * box%cc(nc, 1:nc, 1:nc, iv) + &
            c2 * box%cc(nc-1, 1:nc, 1:nc, iv)
    case (a3_neighb_lowy)
       box%cc(1:nc, 0, 1:nc, iv) = &
            c0 * box%cc(1:nc, 0, 1:nc, iv) + &
            c1 * box%cc(1:nc, 1, 1:nc, iv) + &
            c2 * box%cc(1:nc, 2, 1:nc, iv)
    case (a3_neighb_highy)
       box%cc(1:nc, nc+1, 1:nc, iv) = &
            c0 * box%cc(1:nc, nc+1, 1:nc, iv) + &
            c1 * box%cc(1:nc, nc, 1:nc, iv) + &
            c2 * box%cc(1:nc, nc-1, 1:nc, iv)
    case (a3_neighb_lowz)
       box%cc(1:nc, 1:nc, 0, iv) = &
            c0 * box%cc(1:nc, 1:nc, 0, iv) + &
            c1 * box%cc(1:nc, 1:nc, 1, iv) + &
            c2 * box%cc(1:nc, 1:nc, 2, iv)
    case (a3_neighb_highz)
       box%cc(1:nc, 1:nc, nc+1, iv) = &
            c0 * box%cc(1:nc, 1:nc, nc+1, iv) + &
            c1 * box%cc(1:nc, 1:nc, nc, iv) + &
            c2 * box%cc(1:nc, 1:nc, nc-1, iv)
    end select
  end subroutine bc_to_gc

  subroutine a3_gc_prolong_copy(boxes, id, nb, iv)
    use m_a3_prolong, only: a3_prolong_copy
    type(box3_t), intent(inout)  :: boxes(:)
    integer, intent(in)           :: id
    integer, intent(in)           :: iv
    integer, intent(in)           :: nb
    integer                       :: p_id, nb_dim, lo(3), hi(3)

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

    call a3_prolong_copy(boxes(p_id), boxes(id), iv, low=lo, high=hi)
  end subroutine a3_gc_prolong_copy

  subroutine a3_gc_interp(boxes, id, nb, iv)
    type(box3_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(3)
    real(dp), parameter         :: sixth=1/6.0_dp, third=1/3.0_dp
    integer                     :: k_c1, k_c2, k

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

    if (a3_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 (a3_neighb_dim(nb))
    case (1)
       do k = 1, nc
          k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
          k_c2 = k_c1 + 1 - 2 * iand(k, 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
             boxes(id)%cc(ix, j, k, iv) = &
                  third * boxes(p_nb_id)%cc(ix_c, j_c1, k_c1, iv) + &
                  sixth * boxes(p_nb_id)%cc(ix_c, j_c2, k_c1, iv) + &
                  sixth * boxes(p_nb_id)%cc(ix_c, j_c1, k_c2, iv) + &
                  third * boxes(id)%cc(ix_f, j, k, iv)
          end do
       end do
    case (2)
       do k = 1, nc
          k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
          k_c2 = k_c1 + 1 - 2 * iand(k, 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
             boxes(id)%cc(i, ix, k, iv) = &
                  third * boxes(p_nb_id)%cc(i_c1, ix_c, k_c1, iv) + &
                  sixth * boxes(p_nb_id)%cc(i_c2, ix_c, k_c1, iv) + &
                  sixth * boxes(p_nb_id)%cc(i_c1, ix_c, k_c2, iv) + &
                  third * boxes(id)%cc(i, ix_f, k, iv)
          end do
       end do
    case (3)
       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
          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, j, ix, iv) = &
                  third * boxes(p_nb_id)%cc(i_c1, j_c1, ix_c, iv) + &
                  sixth * boxes(p_nb_id)%cc(i_c1, j_c2, ix_c, iv) + &
                  sixth * boxes(p_nb_id)%cc(i_c2, j_c1, ix_c, iv) + &
                  third * boxes(id)%cc(i, j, ix_f, iv)
          end do
       end do
    end select

  end subroutine a3_gc_interp

  subroutine a3_gc_interp_lim(boxes, id, nb, iv)
    type(box3_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(3)
    real(dp)                    :: c1, c2
    real(dp), parameter         :: sixth=1/6.0_dp, third=1/3.0_dp
    integer                     :: k_c1, k_c2, k
    real(dp)                    :: c3

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

    if (a3_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 (a3_neighb_dim(nb))
    case (1)
       do k = 1, nc
          k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
          k_c2 = k_c1 + 1 - 2 * iand(k, 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
             c1 = boxes(p_nb_id)%cc(ix_c, j_c1, k_c1, iv)
             c2 = boxes(p_nb_id)%cc(ix_c, j_c2, k_c1, iv)
             c3 = boxes(p_nb_id)%cc(ix_c, j_c1, k_c2, iv)
             boxes(id)%cc(ix, j, k, iv) = third * c1 + sixth * c2 + &
                  sixth * c3 + third * boxes(id)%cc(ix_f, j, k, iv)
             if (boxes(id)%cc(ix, j, k, iv) > 2 * c1) &
                  boxes(id)%cc(ix, j, k, iv) = 2 * c1
          end do
       end do
    case (2)
       do k = 1, nc
          k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
          k_c2 = k_c1 + 1 - 2 * iand(k, 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
             c1 = boxes(p_nb_id)%cc(i_c1, ix_c, k_c1, iv)
             c2 = boxes(p_nb_id)%cc(i_c2, ix_c, k_c1, iv)
             c3 = boxes(p_nb_id)%cc(i_c1, ix_c, k_c2, iv)
             boxes(id)%cc(i, ix, k, iv) = third * c1 + sixth * c2 + &
                  sixth * c3 + third * boxes(id)%cc(i, ix_f, k, iv)
             if (boxes(id)%cc(i, ix, k, iv) > 2 * c1) &
                  boxes(id)%cc(i, ix, k, iv) = 2 * c1
          end do
       end do
    case (3)
       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
          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, j_c1, ix_c, iv)
             c2 = boxes(p_nb_id)%cc(i_c1, j_c2, ix_c, iv)
             c3 = boxes(p_nb_id)%cc(i_c2, j_c1, ix_c, iv)
             boxes(id)%cc(i, j, ix, iv) = third * c1 + sixth * c2 + &
                  sixth * c3 + third * boxes(id)%cc(i, j, ix_f, iv)
             if (boxes(id)%cc(i, j, ix, iv) > 2 * c1) &
                  boxes(id)%cc(i, j, ix, iv) = 2 * c1
          end do
       end do
    end select

  end subroutine a3_gc_interp_lim

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

    bc_type = af_bc_neumann
    call a3_set_box_gc(box, nb, iv, 0.0_dp)
  end subroutine a3_bc_neumann_zero

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

    bc_type = af_bc_dirichlet
    call a3_set_box_gc(box, nb, iv, 0.0_dp)
  end subroutine a3_bc_dirichlet_zero

  ! This fills ghost cells near physical boundaries using the same slope
  subroutine a3_bc_continuous(box, nb, iv, bc_type)
    type(box3_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 a3_set_box_gc(box, nb, iv, 0.0_dp)
  end subroutine a3_bc_continuous

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

    ! 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), lo(3):hi(3), iv) = &
         box_nb%cc(nlo(1):nhi(1), nlo(2):nhi(2), nlo(3):nhi(3), iv)
  end subroutine copy_from_nb

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

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

    do nb = 1, a3_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 = a3_neighb_dim(nb)
       lo(:) = 1
       hi(:) = boxes(id)%n_cell
       lo(nb_dim) = -1 + a3_neighb_high_01(nb) * (boxes(id)%n_cell + 3)
       hi(nb_dim) = lo(nb_dim)

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

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

    select case (nb)
    case (a3_neighb_lowx)
       gc_side = box_nb%cc(nc-1, 1:nc, 1:nc, iv)
    case (a3_neighb_highx)
       gc_side = box_nb%cc(2, 1:nc, 1:nc, iv)
    case (a3_neighb_lowy)
       gc_side = box_nb%cc(1:nc, nc-1, 1:nc, iv)
    case (a3_neighb_highy)
       gc_side = box_nb%cc(1:nc, 2, 1:nc, iv)
    case (a3_neighb_lowz)
       gc_side = box_nb%cc(1:nc, 1:nc, nc-1, iv)
    case (a3_neighb_highz)
       gc_side = box_nb%cc(1:nc, 1:nc, 2, iv)
    end select
  end subroutine sides2_from_nb

  subroutine a3_gc2_prolong_linear(boxes, id, nb, iv, gc_side, nc)
    type(box3_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, nc)
    integer                     :: ix, i, j
    integer                     :: i_c1, i_c2, j_c1, j_c2, p_nb_id
    integer                     :: p_id, ix_offset(3)
    integer                     :: k, k_c1, k_c2

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

    select case (a3_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 k = 1, nc
          k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
          k_c2 = k_c1 + 1 - 2 * iand(k, 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, k) = &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c1, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c2, j_c1, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c2, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c1, k_c2, iv)
          end do
       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 k = 1, nc
          k_c1 = ix_offset(3) + ishft(k+1, -1) ! (k+1)/2
          k_c2 = k_c1 + 1 - 2 * iand(k, 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, k) = &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c1, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c2, j_c1, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c2, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c1, k_c2, iv)
          end do
       end do
    case (3)
       k_c1 = ix_offset(3) + ishft(ix+1, -1) ! (k+1)/2
       k_c2 = k_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
          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, j) = &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c1, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c2, j_c1, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c2, k_c1, iv) + &
                  0.25_dp * boxes(p_nb_id)%cc(i_c1, j_c1, k_c2, iv)
          end do
       end do
    end select

  end subroutine a3_gc2_prolong_linear

  ! This fills second ghost cells near physical boundaries using Neumann zero
  subroutine a3_bc2_neumann_zero(boxes, id, nb, iv, gc_side, nc)
    type(box3_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, nc)

    select case (nb)
    case (a3_neighb_lowx)
       gc_side = boxes(id)%cc(2, 1:nc, 1:nc, iv)
    case (a3_neighb_highx)
       gc_side = boxes(id)%cc(nc-1, 1:nc, 1:nc, iv)
    case (a3_neighb_lowy)
       gc_side = boxes(id)%cc(1:nc, 2, 1:nc, iv)
    case (a3_neighb_highy)
       gc_side = boxes(id)%cc(1:nc, nc-1, 1:nc, iv)
    case (a3_neighb_lowz)
       gc_side = boxes(id)%cc(1:nc, 1:nc, 2, iv)
    case (a3_neighb_highz)
       gc_side = boxes(id)%cc(1:nc, 1:nc, nc-1, iv)
    end select
  end subroutine a3_bc2_neumann_zero

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

    select case (nb)
    case (a3_neighb_lowx)
       gc_side = -boxes(id)%cc(2, 1:nc, 1:nc, iv)
    case (a3_neighb_highx)
       gc_side = -boxes(id)%cc(nc-1, 1:nc, 1:nc, iv)
    case (a3_neighb_lowy)
       gc_side = -boxes(id)%cc(1:nc, 2, 1:nc, iv)
    case (a3_neighb_highy)
       gc_side = -boxes(id)%cc(1:nc, nc-1, 1:nc, iv)
    case (a3_neighb_lowz)
       gc_side = -boxes(id)%cc(1:nc, 1:nc, 2, iv)
    case (a3_neighb_highz)
       gc_side = -boxes(id)%cc(1:nc, 1:nc, nc-1, iv)
    end select
  end subroutine a3_bc2_dirichlet_zero

  subroutine a3_corner_gc_extrap(box, ix, iv)
    type(box3_t), intent(inout) :: box
    integer, intent(in)          :: ix(3)
    integer, intent(in)          :: iv
    integer                      :: di(3)

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

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

  subroutine a3_edge_gc_extrap(box, lo, dim, iv)
    type(box3_t), intent(inout) :: box
    integer, intent(in)          :: lo(3)
    integer, intent(in)          :: dim
    integer, intent(in)          :: iv
    integer                      :: di(3), ix(3), ia(3), ib(3), ic(3)
    integer                      :: n, o_dims(3-1)

    ! Dimensions other than/perpendicular to dim
    o_dims = [1 + mod(dim, 3), 1 + mod(dim + 1, 3)]

    ! Index offsets
    di = 1 - 2 * iand(lo, 1)    ! 0 -> di = 1, nc+1 -> di = -1
    di(dim) = 0

    ! Neighbor index in direction one
    ia = lo
    ia(o_dims(1)) = ia(o_dims(1)) + di(o_dims(1))

    ! Neighbor index in direction two
    ib = lo
    ib(o_dims(2)) = ib(o_dims(2)) + di(o_dims(2))

    ! Diagional neighbor index
    ic = lo + di
    ix = lo

    do n = 1, box%n_cell
       ia(dim) = n
       ib(dim) = n
       ic(dim) = n
       ix(dim) = n

       box%cc(ix(1), ix(2), ix(3), iv) = &
            box%cc(ia(1), ia(2), ia(3), iv) + &
            box%cc(ib(1), ib(2), ib(3), iv) - &
            box%cc(ic(1), ic(2), ic(3), iv)
    end do

  end subroutine a3_edge_gc_extrap

end module m_a3_ghostcell