mod_ghostcells_update.t

Go to the documentation of this file.

!> update ghost cells of all blocks including physical boundaries 
module mod_ghostcells_update

  implicit none
  ! Special buffer for pole copy
  type wbuffer
    double precision, dimension(:^D&,:), allocatable :: w
  end type wbuffer

  ! A switch of update physical boundary or not
  logical, public :: bcphys=.true.
  integer :: ixm^l, ixcog^l, ixcom^l

  ! The number of interleaving sending buffers for ghost cells
  integer, parameter :: npwbuf=2

  ! The first index goes from -1:2, where -1 is used when a block touches the
  ! lower boundary, 1 when a block touches an upper boundary, and 0 a situation
  ! away from boundary conditions, 2 when a block touched both lower and upper
  ! boundary

  ! index ranges to send (S) to sibling blocks, receive (R) from sibling blocks
  integer, dimension(-1:2,-1:1) :: ixs_srl_^l, ixr_srl_^l

  ! index ranges of staggered variables to send (S) to sibling blocks, receive (R) from sibling blocks
  integer, dimension(^ND,-1:1) :: ixs_srl_stg_^l, ixr_srl_stg_^l

  ! index ranges to send (S) restricted (r) ghost cells to coarser blocks 
  integer, dimension(-1:1,-1:1) :: ixs_r_^l

  ! index ranges of staggered variables to send (S) restricted (r) ghost cells to coarser blocks 
  integer, dimension(^ND,-1:1) :: ixs_r_stg_^l

  ! index ranges to receive restriced ghost cells from finer blocks 
  integer, dimension(-1:1, 0:3) :: ixr_r_^l

  ! index ranges of staggered variables to receive restriced ghost cells from finer blocks 
  integer, dimension(^ND,0:3)  :: ixr_r_stg_^l

  ! send prolongated (p) ghost cells to finer blocks, receive prolongated from coarser blocks
  integer, dimension(-1:1, 0:3) :: ixs_p_^l, ixr_p_^l

  ! send prolongated (p) staggered ghost cells to finer blocks, receive prolongated from coarser blocks
  integer, dimension(^ND,0:3)  :: ixs_p_stg_^l, ixr_p_stg_^l

  ! number of MPI receive-send pairs, srl: same refinement level; r: restrict; p: prolong
  integer :: nrecv_bc_srl, nsend_bc_srl, nrecv_bc_r, nsend_bc_r, nrecv_bc_p, nsend_bc_p

  ! record index position of buffer arrays
  integer :: ibuf_send_srl, ibuf_recv_srl, ibuf_send_r, ibuf_recv_r, ibuf_send_p, ibuf_recv_p

  ! count of times of send and receive
  integer :: isend_srl, irecv_srl, isend_r, irecv_r, isend_p, irecv_p

  ! count of times of send and receive for cell center ghost cells
  integer :: isend_c, irecv_c

  ! tag of MPI send and recv
  integer, private :: itag

  ! total sizes = cell-center normal flux + stagger-grid flux of send and receive
  integer, dimension(-1:1^D&) :: sizes_srl_send_total, sizes_srl_recv_total

  ! sizes of buffer arrays for center-grid variable for siblings and restrict
  integer, dimension(:), allocatable :: recvrequest_c_sr, sendrequest_c_sr
  integer, dimension(:,:), allocatable :: recvstatus_c_sr, sendstatus_c_sr

  ! sizes of buffer arrays for center-grid variable for prolongation
  integer, dimension(:), allocatable :: recvrequest_c_p, sendrequest_c_p
  integer, dimension(:,:), allocatable :: recvstatus_c_p, sendstatus_c_p

  ! sizes of buffer arrays for stagger-grid variable
  integer, dimension(^ND,-1:1^D&) :: sizes_srl_send_stg, sizes_srl_recv_stg

  integer, dimension(:), allocatable :: recvrequest_srl, sendrequest_srl
  integer, dimension(:,:), allocatable :: recvstatus_srl, sendstatus_srl
 
  ! buffer arrays for send and receive of siblings, allocate in build_connectivity
  double precision, dimension(:), allocatable :: recvbuffer_srl, sendbuffer_srl

  integer, dimension(:), allocatable :: recvrequest_r, sendrequest_r
  integer, dimension(:,:), allocatable :: recvstatus_r, sendstatus_r

  ! buffer arrays for send and receive in restriction
  double precision, dimension(:), allocatable :: recvbuffer_r, sendbuffer_r

  integer, dimension(:), allocatable :: recvrequest_p, sendrequest_p
  integer, dimension(:,:), allocatable :: recvstatus_p, sendstatus_p

  ! buffer arrays for send and receive in prolongation
  double precision, dimension(:), allocatable :: recvbuffer_p, sendbuffer_p

  ! sizes to allocate buffer arrays for send and receive for restriction
  integer, dimension(-1:1^D&)     :: sizes_r_send_total
  integer, dimension(0:3^D&)      :: sizes_r_recv_total
  integer, dimension(^ND,-1:1^D&) :: sizes_r_send_stg
  integer, dimension(^ND,0:3^D&)  :: sizes_r_recv_stg

  ! sizes to allocate buffer arrays for send and receive for restriction
  integer, dimension(0:3^D&)      :: sizes_p_send_total, sizes_p_recv_total
  integer, dimension(^ND,0:3^D&)  :: sizes_p_send_stg, sizes_p_recv_stg

  ! There are two variants, _f indicates that all flux variables are filled,
  ! whereas _p means that part of the variables is filled 
  ! Furthermore _r_ stands for restrict, _p_ for prolongation.
  integer, dimension(-1:2^D&,-1:1^D&), target :: type_send_srl_f, type_recv_srl_f
  integer, dimension(-1:1^D&,-1:1^D&), target :: type_send_r_f
  integer, dimension(-1:1^D&, 0:3^D&), target :: type_recv_r_f, type_send_p_f, type_recv_p_f
  integer, dimension(-1:2^D&,-1:1^D&), target :: type_send_srl_p1, type_recv_srl_p1
  integer, dimension(-1:1^D&,-1:1^D&), target :: type_send_r_p1
  integer, dimension(-1:1^D&, 0:3^D&), target :: type_recv_r_p1, type_send_p_p1, type_recv_p_p1
  integer, dimension(-1:2^D&,-1:1^D&), target :: type_send_srl_p2, type_recv_srl_p2
  integer, dimension(-1:1^D&,-1:1^D&), target :: type_send_r_p2
  integer, dimension(-1:1^D&, 0:3^D&), target :: type_recv_r_p2, type_send_p_p2, type_recv_p_p2
  integer, dimension(:^D&,:^D&), pointer :: type_send_srl, type_recv_srl, type_send_r
  integer, dimension(:^D&,:^D&), pointer :: type_recv_r, type_send_p, type_recv_p

contains

  subroutine init_bc()
    use mod_global_parameters 
    use mod_physics, only: phys_req_diagonal

    integer :: nghostcellsCo, interpolation_order
    integer :: nx^D, nxCo^D, ixG^L, i^D, ic^D, inc^D, idir

    ixg^l=ixg^ll;
    ixm^l=ixg^l^lsubnghostcells;
    ixcogmin^d=1;
    !ixCoGmax^D=ixGmax^D/2+nghostcells;
    ixcogmax^d=(ixghi^d-2*nghostcells)/2+2*nghostcells;

    ixcom^l=ixcog^l^lsubnghostcells;

    nx^d=ixmmax^d-ixmmin^d+1;
    nxco^d=nx^d/2;

    if(ghost_copy) then
       interpolation_order=1
    else
       interpolation_order=2
    end if
    nghostcellsco=int((nghostcells+1)/2)
    
    if (nghostcellsco+interpolation_order-1>nghostcells) then
       call mpistop("interpolation order for prolongation in getbc too high")
    end if

    ! (iib,i) index has following meanings: iib = 0 means it is not at any physical boundary
    ! iib=-1 means it is at the minimum side of a physical boundary  
    ! iib= 1 means it is at the maximum side of a physical boundary  
    ! i=-1 means subregion prepared for the neighbor at its minimum side 
    ! i= 1 means subregion prepared for the neighbor at its maximum side 
    {
    ixs_srl_min^d(:,-1)=ixmmin^d
    ixs_srl_min^d(:, 0)=ixmmin^d
    ixs_srl_min^d(:, 1)=ixmmax^d+1-nghostcells
    ixs_srl_max^d(:,-1)=ixmmin^d-1+nghostcells
    ixs_srl_max^d(:, 0)=ixmmax^d
    ixs_srl_max^d(:, 1)=ixmmax^d
     
    ixr_srl_min^d(:,-1)=1
    ixr_srl_min^d(:, 0)=ixmmin^d
    ixr_srl_min^d(:, 1)=ixmmax^d+1
    ixr_srl_max^d(:,-1)=nghostcells
    ixr_srl_max^d(:, 0)=ixmmax^d
    ixr_srl_max^d(:, 1)=ixgmax^d
    
    ixs_r_min^d(:,-1)=ixcommin^d
    ixs_r_min^d(:, 0)=ixcommin^d
    ixs_r_min^d(:, 1)=ixcommax^d+1-nghostcells
    ixs_r_max^d(:,-1)=ixcommin^d-1+nghostcells
    ixs_r_max^d(:, 0)=ixcommax^d
    ixs_r_max^d(:, 1)=ixcommax^d
    
    ixr_r_min^d(:, 0)=1
    ixr_r_min^d(:, 1)=ixmmin^d
    ixr_r_min^d(:, 2)=ixmmin^d+nxco^d
    ixr_r_min^d(:, 3)=ixmmax^d+1
    ixr_r_max^d(:, 0)=nghostcells
    ixr_r_max^d(:, 1)=ixmmin^d-1+nxco^d
    ixr_r_max^d(:, 2)=ixmmax^d
    ixr_r_max^d(:, 3)=ixgmax^d

    ixs_p_min^d(:, 0)=ixmmin^d-(interpolation_order-1)
    ixs_p_min^d(:, 1)=ixmmin^d-(interpolation_order-1)
    ixs_p_min^d(:, 2)=ixmmin^d+nxco^d-nghostcellsco-(interpolation_order-1)
    ixs_p_min^d(:, 3)=ixmmax^d+1-nghostcellsco-(interpolation_order-1)
    ixs_p_max^d(:, 0)=ixmmin^d-1+nghostcellsco+(interpolation_order-1)
    ixs_p_max^d(:, 1)=ixmmin^d-1+nxco^d+nghostcellsco+(interpolation_order-1)
    ixs_p_max^d(:, 2)=ixmmax^d+(interpolation_order-1)
    ixs_p_max^d(:, 3)=ixmmax^d+(interpolation_order-1)

    if(.not.phys_req_diagonal) then
      ! exclude ghost-cell region when diagonal cells are unknown
      ixs_p_min^d(:, 0)=ixmmin^d
      ixs_p_max^d(:, 3)=ixmmax^d
      ixs_p_max^d(:, 1)=ixmmin^d-1+nxco^d+(interpolation_order-1)
      ixs_p_min^d(:, 2)=ixmmin^d+nxco^d-(interpolation_order-1)
    end if

    ixr_p_min^d(:, 0)=ixcommin^d-nghostcellsco-(interpolation_order-1)
    ixr_p_min^d(:, 1)=ixcommin^d-(interpolation_order-1)
    ixr_p_min^d(:, 2)=ixcommin^d-nghostcellsco-(interpolation_order-1)
    ixr_p_min^d(:, 3)=ixcommax^d+1-(interpolation_order-1)
    ixr_p_max^d(:, 0)=nghostcells+(interpolation_order-1)
    ixr_p_max^d(:, 1)=ixcommax^d+nghostcellsco+(interpolation_order-1)
    ixr_p_max^d(:, 2)=ixcommax^d+(interpolation_order-1)
    ixr_p_max^d(:, 3)=ixcommax^d+nghostcellsco+(interpolation_order-1)

    if(.not.phys_req_diagonal) then
      ixr_p_max^d(:, 0)=nghostcells
      ixr_p_min^d(:, 3)=ixcommax^d+1
      ixr_p_max^d(:, 1)=ixcommax^d+(interpolation_order-1)
      ixr_p_min^d(:, 2)=ixcommin^d-(interpolation_order-1)
    end if

    \}

    if (stagger_grid) then
      allocate(pole_buf%ws(ixgs^t,nws))
      ! Staggered (face-allocated) variables
      do idir=1,ndim
      { ixs_srl_stg_min^d(idir,-1)=ixmmin^d-kr(idir,^d)
        ixs_srl_stg_max^d(idir,-1)=ixmmin^d-1+nghostcells
        ixs_srl_stg_min^d(idir,0) =ixmmin^d-kr(idir,^d)
        ixs_srl_stg_max^d(idir,0) =ixmmax^d
        ixs_srl_stg_min^d(idir,1) =ixmmax^d-nghostcells+1-kr(idir,^d)
        ixs_srl_stg_max^d(idir,1) =ixmmax^d
        
        ixr_srl_stg_min^d(idir,-1)=1-kr(idir,^d)
        ixr_srl_stg_max^d(idir,-1)=nghostcells
        ixr_srl_stg_min^d(idir,0) =ixmmin^d-kr(idir,^d)
        ixr_srl_stg_max^d(idir,0) =ixmmax^d
        ixr_srl_stg_min^d(idir,1) =ixmmax^d+1-kr(idir,^d)
        ixr_srl_stg_max^d(idir,1) =ixgmax^d

        ixs_r_stg_min^d(idir,-1)=ixcommin^d-kr(idir,^d)
        ixs_r_stg_max^d(idir,-1)=ixcommin^d-1+nghostcells
        ixs_r_stg_min^d(idir,0) =ixcommin^d-kr(idir,^d)
        ixs_r_stg_max^d(idir,0) =ixcommax^d
        ixs_r_stg_min^d(idir,1) =ixcommax^d+1-nghostcells-kr(idir,^d)
        ixs_r_stg_max^d(idir,1) =ixcommax^d
 
        ixr_r_stg_min^d(idir,0)=1-kr(idir,^d)
        ixr_r_stg_max^d(idir,0)=nghostcells
        ixr_r_stg_min^d(idir,1)=ixmmin^d-kr(idir,^d)
        ixr_r_stg_max^d(idir,1)=ixmmin^d-1+nxco^d
        ixr_r_stg_min^d(idir,2)=ixmmin^d+nxco^d-kr(idir,^d)
        ixr_r_stg_max^d(idir,2)=ixmmax^d
        ixr_r_stg_min^d(idir,3)=ixmmax^d+1-kr(idir,^d)
        ixr_r_stg_max^d(idir,3)=ixgmax^d
        \}
       {if (idir==^d) then
          ! Parallel components
          {
          ixs_p_stg_min^d(idir,0)=ixmmin^d-1 ! -1 to make redundant 
          ixs_p_stg_max^d(idir,0)=ixmmin^d-1+nghostcellsco
          ixs_p_stg_min^d(idir,1)=ixmmin^d-1 ! -1 to make redundant 
          ixs_p_stg_max^d(idir,1)=ixmmin^d-1+nxco^d+nghostcellsco
          ixs_p_stg_min^d(idir,2)=ixmmax^d-nxco^d-nghostcellsco
          ixs_p_stg_max^d(idir,2)=ixmmax^d
          ixs_p_stg_min^d(idir,3)=ixmmax^d-nghostcellsco
          ixs_p_stg_max^d(idir,3)=ixmmax^d

          ixr_p_stg_min^d(idir,0)=ixcommin^d-1-nghostcellsco
          ixr_p_stg_max^d(idir,0)=ixcommin^d-1
          ixr_p_stg_min^d(idir,1)=ixcommin^d-1 ! -1 to make redundant 
          ixr_p_stg_max^d(idir,1)=ixcommax^d+nghostcellsco
          ixr_p_stg_min^d(idir,2)=ixcommin^d-1-nghostcellsco
          ixr_p_stg_max^d(idir,2)=ixcommax^d
          ixr_p_stg_min^d(idir,3)=ixcommax^d+1-1 ! -1 to make redundant 
          ixr_p_stg_max^d(idir,3)=ixcommax^d+nghostcellsco
          \}
        else
          {
          ! Perpendicular component
          ixs_p_stg_min^d(idir,0)=ixmmin^d
          ixs_p_stg_max^d(idir,0)=ixmmin^d-1+nghostcellsco+(interpolation_order-1)
          ixs_p_stg_min^d(idir,1)=ixmmin^d
          ixs_p_stg_max^d(idir,1)=ixmmin^d-1+nxco^d+nghostcellsco+(interpolation_order-1)
          ixs_p_stg_min^d(idir,2)=ixmmax^d+1-nxco^d-nghostcellsco-(interpolation_order-1)
          ixs_p_stg_max^d(idir,2)=ixmmax^d
          ixs_p_stg_min^d(idir,3)=ixmmax^d+1-nghostcellsco-(interpolation_order-1)
          ixs_p_stg_max^d(idir,3)=ixmmax^d
 
          ixr_p_stg_min^d(idir,0)=ixcommin^d-nghostcellsco-(interpolation_order-1)
          ixr_p_stg_max^d(idir,0)=ixcommin^d-1
          ixr_p_stg_min^d(idir,1)=ixcommin^d
          ixr_p_stg_max^d(idir,1)=ixcommax^d+nghostcellsco+(interpolation_order-1)
          ixr_p_stg_min^d(idir,2)=ixcommin^d-nghostcellsco-(interpolation_order-1)
          ixr_p_stg_max^d(idir,2)=ixcommax^d
          ixr_p_stg_min^d(idir,3)=ixcommax^d+1
          ixr_p_stg_max^d(idir,3)=ixcommax^d+nghostcellsco+(interpolation_order-1)
          \}
        end if
        }
      end do
      ! calculate sizes for buffer arrays for siblings
      {do i^db=-1,1\}
         ! Staggered (face-allocated) variables
         do idir=1,ndim
           sizes_srl_send_stg(idir,i^d)={(ixs_srl_stg_max^d(idir,i^d)-ixs_srl_stg_min^d(idir,i^d)+1)|*}
           sizes_srl_recv_stg(idir,i^d)={(ixr_srl_stg_max^d(idir,i^d)-ixr_srl_stg_min^d(idir,i^d)+1)|*}
           sizes_r_send_stg(idir,i^d)={(ixs_r_stg_max^d(idir,i^d)-ixs_r_stg_min^d(idir,i^d)+1)|*}
         end do
         sizes_srl_send_total(i^d)=sum(sizes_srl_send_stg(:,i^d))
         sizes_srl_recv_total(i^d)=sum(sizes_srl_recv_stg(:,i^d))
         sizes_r_send_total(i^d)=sum(sizes_r_send_stg(:,i^d))
      {end do\}

      {do i^db=0,3\}
         ! Staggered (face-allocated) variables
           do idir=1,ndim
             sizes_r_recv_stg(idir,i^d)={(ixr_r_stg_max^d(idir,i^d)-ixr_r_stg_min^d(idir,i^d)+1)|*}
             sizes_p_send_stg(idir,i^d)={(ixs_p_stg_max^d(idir,i^d)-ixs_p_stg_min^d(idir,i^d)+1)|*}
             sizes_p_recv_stg(idir,i^d)={(ixr_p_stg_max^d(idir,i^d)-ixr_p_stg_min^d(idir,i^d)+1)|*}
           end do
           sizes_r_recv_total(i^d)=sum(sizes_r_recv_stg(:,i^d))
           sizes_p_send_total(i^d)=sum(sizes_p_send_stg(:,i^d))
           sizes_p_recv_total(i^d)=sum(sizes_p_recv_stg(:,i^d))
      {end do\}
    else
      ! extend index range to physical boundary
      {
      ixs_srl_min^d(-1,0)=1
      ixs_srl_min^d( 1,0)=ixmmin^d
      ixs_srl_min^d( 2,0)=1
      ixs_srl_max^d(-1,0)=ixmmax^d
      ixs_srl_max^d( 1,0)=ixgmax^d
      ixs_srl_max^d( 2,0)=ixgmax^d
       
      ixr_srl_min^d(-1,0)=1
      ixr_srl_min^d( 1,0)=ixmmin^d
      ixr_srl_min^d( 2,0)=1
      ixr_srl_max^d(-1,0)=ixmmax^d
      ixr_srl_max^d( 1,0)=ixgmax^d
      ixr_srl_max^d( 2,0)=ixgmax^d
      
      ixs_r_min^d(-1,0)=1
      ixs_r_min^d( 1,0)=ixcommin^d
      ixs_r_max^d(-1,0)=ixcommax^d
      ixs_r_max^d( 1,0)=ixcogmax^d
      
      ixr_r_min^d(-1,1)=1
      ixr_r_max^d(-1,1)=ixmmin^d-1+nxco^d
      ixr_r_min^d( 1,2)=ixmmin^d+nxco^d
      ixr_r_max^d( 1,2)=ixgmax^d

      ixs_p_min^d(-1,1)=1
      ixs_p_max^d( 1,2)=ixgmax^d

      ixr_p_min^d(-1,1)=1
      ixr_p_max^d( 1,2)=ixcogmax^d
      \}
    end if

  end subroutine init_bc

  subroutine create_bc_mpi_datatype(nwstart,nwbc) 
    use mod_global_parameters 

    integer, intent(in) :: nwstart, nwbc
    integer :: i^D, ic^D, inc^D, iib^D

    {do i^db=-1,1\}
      if (i^d==0|.and.) cycle
      {do iib^db=-1,2\}
         call get_bc_comm_type(type_send_srl(iib^d,i^d),ixs_srl_^l(iib^d,i^d),ixg^ll,nwstart,nwbc)
         call get_bc_comm_type(type_recv_srl(iib^d,i^d),ixr_srl_^l(iib^d,i^d),ixg^ll,nwstart,nwbc)
         if (iib^d==2|.or.) cycle
         call get_bc_comm_type(type_send_r(iib^d,i^d),  ixs_r_^l(iib^d,i^d),ixcog^l,nwstart,nwbc)
         {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
            inc^db=2*i^db+ic^db\}
            call get_bc_comm_type(type_recv_r(iib^d,inc^d),ixr_r_^l(iib^d,inc^d), ixg^ll,nwstart,nwbc)
            call get_bc_comm_type(type_send_p(iib^d,inc^d),ixs_p_^l(iib^d,inc^d), ixg^ll,nwstart,nwbc)
            call get_bc_comm_type(type_recv_p(iib^d,inc^d),ixr_p_^l(iib^d,inc^d),ixcog^l,nwstart,nwbc)
         {end do\}
      {end do\}
    {end do\}
  
  end subroutine create_bc_mpi_datatype

  subroutine get_bc_comm_type(comm_type,ix^L,ixG^L,nwstart,nwbc)
    use mod_global_parameters 
  
    integer, intent(inout) :: comm_type
    integer, intent(in) :: ix^L, ixG^L, nwstart, nwbc
    
    integer, dimension(ndim+1) :: fullsize, subsize, start

    ^d&fullsize(^d)=ixgmax^d;
    fullsize(ndim+1)=nw
    ^d&subsize(^d)=ixmax^d-ixmin^d+1;
    subsize(ndim+1)=nwbc
    ^d&start(^d)=ixmin^d-1;
    start(ndim+1)=nwstart-1
    
    call mpi_type_create_subarray(ndim+1,fullsize,subsize,start,mpi_order_fortran, &
                                  mpi_double_precision,comm_type,ierrmpi)
    call mpi_type_commit(comm_type,ierrmpi)
    
  end subroutine get_bc_comm_type

  subroutine put_bc_comm_types()
    use mod_global_parameters 
 
    integer :: i^D, ic^D, inc^D, iib^D

    {do i^db=-1,1\}
       if (i^d==0|.and.) cycle
       {do iib^db=-1,2\}
           call mpi_type_free(type_send_srl(iib^d,i^d),ierrmpi)
           call mpi_type_free(type_recv_srl(iib^d,i^d),ierrmpi)
           if (levmin==levmax) cycle
           if (iib^d==2|.or.) cycle
           call mpi_type_free(type_send_r(iib^d,i^d),ierrmpi)
           {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
              inc^db=2*i^db+ic^db\}
              call mpi_type_free(type_recv_r(iib^d,inc^d),ierrmpi)
              call mpi_type_free(type_send_p(iib^d,inc^d),ierrmpi)
              call mpi_type_free(type_recv_p(iib^d,inc^d),ierrmpi)
           {end do\}
       {end do\}
    {end do\}
  
  end subroutine put_bc_comm_types

  !> do update ghost cells of all blocks including physical boundaries
  subroutine getbc(time,qdt,psb,nwstart,nwbc,req_diag)
    use mod_global_parameters
    use mod_physics

    double precision, intent(in)      :: time, qdt
    type(state), target               :: psb(max_blocks)
    integer, intent(in)               :: nwstart ! Fill from nwstart
    integer, intent(in)               :: nwbc    ! Number of variables to fill
    logical, intent(in), optional     :: req_diag ! If false, skip diagonal ghost cells

    double precision :: time_bcin
    integer :: ipole, nwhead, nwtail
    integer :: iigrid, igrid, ineighbor, ipe_neighbor, isizes
    integer :: ixR^L, ixS^L
    integer :: i^D, n_i^D, ic^D, inc^D, n_inc^D, iib^D, idir
    ! store physical boundary indicating index
    integer :: idphyb(ndim,max_blocks)
    integer :: isend_buf(npwbuf), ipwbuf, nghostcellsco
    ! index pointer for buffer arrays as a start for a segment
    integer :: ibuf_start, ibuf_next
    ! shapes of reshape
    integer, dimension(1) :: shapes
    logical  :: req_diagonal
    type(wbuffer) :: pwbuf(npwbuf)

    time_bcin=mpi_wtime()

    nwhead=nwstart
    nwtail=nwstart+nwbc-1

    req_diagonal = .true.
    if (present(req_diag)) req_diagonal = req_diag

    ! fill internal physical boundary
    if (internalboundary) then 
       call getintbc(time,ixg^ll)
    end if

    ! fill physical-boundary ghost cells before internal ghost-cell values exchange
    if(bcphys.and. .not.stagger_grid) then
      !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if(.not.phyboundblock(igrid)) cycle
        call fill_boundary_before_gc(igrid)
      end do
      !$OMP END PARALLEL DO
    end if

    ! default : no singular axis
    ipole=0
    irecv_c=0
    isend_c=0
    isend_buf=0
    ipwbuf=1

    if(stagger_grid) then
      ibuf_recv_srl=1
      ibuf_recv_r=1
      ibuf_recv_p=1
      ibuf_send_srl=1
      ibuf_send_r=1
      ibuf_send_p=1
      irecv_srl=0
      irecv_r=0
      irecv_p=0
      isend_srl=0
      isend_r=0
      isend_p=0
    end if

    ! prepare coarse values to send to coarser neighbors
    !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      if(any(neighbor_type(:^d&,igrid)==neighbor_coarse)) then
        call coarsen_grid(psb(igrid),ixg^ll,ixm^l,psc(igrid),ixcog^l,ixcom^l)
       {do i^db=-1,1\}
          if(skip_direction([ i^d ])) cycle
          if(neighbor_type(i^d,igrid)==neighbor_coarse) call fill_coarse_boundary(igrid,i^d)
       {end do\}
      end if
    end do
    !$OMP END PARALLEL DO

    ! MPI receive ghost-cell values from sibling blocks and finer neighbors in different processors
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      call identifyphysbound(ps(igrid),iib^d)   
      ^d&idphyb(^d,igrid)=iib^d;
      {do i^db=-1,1\}
         if (skip_direction([ i^d ])) cycle
         select case (neighbor_type(i^d,igrid))
         case (neighbor_sibling)
            call bc_recv_srl
         case (neighbor_fine)
            call bc_recv_restrict
         end select
      {end do\}
    end do

    ! MPI send ghost-cell values to sibling blocks and coarser neighbors in different processors
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ^d&iib^d=idphyb(^d,igrid);
      {do i^db=-1,1\}
         if(skip_direction([ i^d ])) cycle
         select case (neighbor_type(i^d,igrid))
         case (neighbor_sibling)
            call bc_send_srl
         case (neighbor_coarse)
            call bc_send_restrict
         end select
      {end do\}
    end do

    ! fill ghost-cell values of sibling blocks and coarser neighbors in the same processor
    !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid,iib^D)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ^d&iib^d=idphyb(^d,igrid);
      {do i^db=-1,1\}
        if(skip_direction([ i^d ])) cycle
         select case (neighbor_type(i^d,igrid))
         case(neighbor_sibling) 
           call bc_fill_srl(igrid,i^d,iib^d)
         case(neighbor_coarse)
           call bc_fill_restrict(igrid,i^d,iib^d)
         end select
      {end do\}
    end do
    !$OMP END PARALLEL DO

    call mpi_waitall(irecv_c,recvrequest_c_sr,recvstatus_c_sr,ierrmpi)
    call mpi_waitall(isend_c,sendrequest_c_sr,sendstatus_c_sr,ierrmpi)

    if(stagger_grid) then
      call mpi_waitall(nrecv_bc_srl,recvrequest_srl,recvstatus_srl,ierrmpi)
      call mpi_waitall(nsend_bc_srl,sendrequest_srl,sendstatus_srl,ierrmpi)
      call mpi_waitall(nrecv_bc_r,recvrequest_r,recvstatus_r,ierrmpi)
      call mpi_waitall(nsend_bc_r,sendrequest_r,sendstatus_r,ierrmpi)
      ! unpack the received data from sibling blocks and finer neighbors to fill ghost-cell staggered values 
      ibuf_recv_srl=1
      ibuf_recv_r=1
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        ^d&iib^d=idphyb(^d,igrid);
       {do i^db=-1,1\}
          if (skip_direction([ i^d ])) cycle
          select case (neighbor_type(i^d,igrid))
          case (neighbor_sibling)
             call bc_fill_srl_stg
          case (neighbor_fine)
             call bc_fill_restrict_stg
          end select
       {end do\}
      end do
    end if

    do ipwbuf=1,npwbuf
       if (isend_buf(ipwbuf)/=0) deallocate(pwbuf(ipwbuf)%w)
    end do

    irecv_c=0
    isend_c=0
    isend_buf=0
    ipwbuf=1

    ! MPI receive ghost-cell values from coarser neighbors in different processors
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ^d&iib^d=idphyb(^d,igrid);
      {do i^db=-1,1\}
         if (skip_direction([ i^d ])) cycle
         if (neighbor_type(i^d,igrid)==neighbor_coarse) call bc_recv_prolong
      {end do\}
    end do
    ! MPI send ghost-cell values to finer neighbors in different processors
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ^d&iib^d=idphyb(^d,igrid);
      {do i^db=-1,1\}
         if (skip_direction([ i^d ])) cycle
         if (neighbor_type(i^d,igrid)==neighbor_fine) call bc_send_prolong
      {end do\}
    end do

    ! fill coarse ghost-cell values of finer neighbors in the same processor
    !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid,iib^D)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ^d&iib^d=idphyb(^d,igrid);
      {do i^db=-1,1\}
         if (skip_direction([ i^d ])) cycle
         if (neighbor_type(i^d,igrid)==neighbor_fine) call bc_fill_prolong(igrid,i^d,iib^d)
      {end do\}
    end do
    !$OMP END PARALLEL DO

    call mpi_waitall(irecv_c,recvrequest_c_p,recvstatus_c_p,ierrmpi)
    call mpi_waitall(isend_c,sendrequest_c_p,sendstatus_c_p,ierrmpi)

    if(stagger_grid) then
      call mpi_waitall(nrecv_bc_p,recvrequest_p,recvstatus_p,ierrmpi)
      call mpi_waitall(nsend_bc_p,sendrequest_p,sendstatus_p,ierrmpi)

      ! fill coarser representative ghost cells after receipt
      ibuf_recv_p=1
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        ^d&iib^d=idphyb(^d,igrid);
        {do i^db=-1,1\}
           if (skip_direction([ i^d ])) cycle
           if(neighbor_type(i^d,igrid)==neighbor_coarse) call bc_fill_prolong_stg
        {end do\}
      end do
    end if
    ! do prolongation on the ghost-cell values based on the received coarse values from coarser neighbors 
    !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      call gc_prolong(igrid)
    end do
    !$OMP END PARALLEL DO
    
    do ipwbuf=1,npwbuf
       if (isend_buf(ipwbuf)/=0) deallocate(pwbuf(ipwbuf)%w)
    end do

    ! fill physical boundary ghost cells after internal ghost-cell values exchange
    if(bcphys.and.stagger_grid) then
      !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if(.not.phyboundblock(igrid)) cycle
        call fill_boundary_after_gc(igrid)
      end do
      !$OMP END PARALLEL DO
    end if

     ! modify normal component of magnetic field to fix divB=0 
    if(bcphys.and.associated(phys_boundary_adjust)) then
      !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if(.not.phyboundblock(igrid)) cycle
        call phys_boundary_adjust(igrid,psb)
      end do
      !$OMP END PARALLEL DO
    end if
    
    if (nwaux>0) then
      !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        call fix_auxiliary(igrid)
      end do
      !$OMP END PARALLEL DO
    end if
    
    time_bc=time_bc+(mpi_wtime()-time_bcin)
    
    contains

      logical function skip_direction(dir)
        integer, intent(in) :: dir(^nd)

        if (all(dir == 0)) then
           skip_direction = .true.
        else if (.not. req_diagonal .and. count(dir /= 0) > 1) then
           skip_direction = .true.
        else
           skip_direction = .false.
        end if
      end function skip_direction

      !> Physical boundary conditions
      subroutine fill_boundary_before_gc(igrid)

        integer, intent(in) :: igrid

        integer :: idims,iside,i^D,k^L,ixB^L

        block=>psb(igrid)
        ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
        do idims=1,ndim
          ! to avoid using as yet unknown corner info in more than 1D, we
          ! fill only interior mesh ranges of the ghost cell ranges at first,
          ! and progressively enlarge the ranges to include corners later
          {
           kmin^d=merge(0, 1, idims==^d)
           kmax^d=merge(0, 1, idims==^d)
           ixbmin^d=ixglo^d+kmin^d*nghostcells
           ixbmax^d=ixghi^d-kmax^d*nghostcells
          \}
          {^iftwod
          if(idims > 1 .and. neighbor_type(-1,0,igrid)==neighbor_boundary) ixbmin1=ixglo1
          if(idims > 1 .and. neighbor_type( 1,0,igrid)==neighbor_boundary) ixbmax1=ixghi1}
          {^ifthreed
          if(idims > 1 .and. neighbor_type(-1,0,0,igrid)==neighbor_boundary) ixbmin1=ixglo1
          if(idims > 1 .and. neighbor_type( 1,0,0,igrid)==neighbor_boundary) ixbmax1=ixghi1
          if(idims > 2 .and. neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixbmin2=ixglo2
          if(idims > 2 .and. neighbor_type(0, 1,0,igrid)==neighbor_boundary) ixbmax2=ixghi2}
          do iside=1,2
            i^d=kr(^d,idims)*(2*iside-3);
            if (aperiodb(idims)) then
              if (neighbor_type(i^d,igrid) /= neighbor_boundary .and. &
                   .not. psb(igrid)%is_physical_boundary(2*idims-2+iside)) cycle
            else
              if (neighbor_type(i^d,igrid) /= neighbor_boundary) cycle
            end if
            call bc_phys(iside,idims,time,qdt,psb(igrid),ixg^ll,ixb^l)
          end do
        end do

      end subroutine fill_boundary_before_gc

      !> Physical boundary conditions
      subroutine fill_boundary_after_gc(igrid)

        integer, intent(in) :: igrid

        integer :: idims,iside,i^D,k^L,ixB^L

        block=>psb(igrid)
        ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
        do idims=1,ndim
          ! to avoid using as yet unknown corner info in more than 1D, we
          ! fill only interior mesh ranges of the ghost cell ranges at first,
          ! and progressively enlarge the ranges to include corners later
          kmin1=0; kmax1=0;
          {^iftwod
           kmin2=merge(1, 0,  idims .lt. 2 .and. neighbor_type(0,-1,igrid)==1)
           kmax2=merge(1, 0,  idims .lt. 2 .and. neighbor_type(0, 1,igrid)==1)}
          {^ifthreed
           kmin2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0,-1,0,igrid)==1)
           kmax2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0, 1,0,igrid)==1)
           kmin3=merge(1, 0, idims .lt. 3 .and. neighbor_type(0,0,-1,igrid)==1)
           kmax3=merge(1, 0, idims .lt. 3 .and. neighbor_type(0,0, 1,igrid)==1)}
          ixbmin^d=ixglo^d+kmin^d*nghostcells;
          ixbmax^d=ixghi^d-kmax^d*nghostcells;
          do iside=1,2
            i^d=kr(^d,idims)*(2*iside-3);
            if (aperiodb(idims)) then 
              if (neighbor_type(i^d,igrid) /= neighbor_boundary .and. &
                 .not. psb(igrid)%is_physical_boundary(2*idims-2+iside)) cycle
            else 
              if (neighbor_type(i^d,igrid) /= neighbor_boundary) cycle
            end if
            call bc_phys(iside,idims,time,qdt,psb(igrid),ixg^ll,ixb^l)
          end do
        end do

      end subroutine fill_boundary_after_gc

      !> Receive from sibling at same refinement level
      subroutine bc_recv_srl

        ipe_neighbor=neighbor(2,i^d,igrid)
        if (ipe_neighbor/=mype) then
           irecv_c=irecv_c+1
           itag=(3**^nd+4**^nd)*(igrid-1)+{(i^d+1)*3**(^d-1)+}
           call mpi_irecv(psb(igrid)%w,1,type_recv_srl(iib^d,i^d), &
                          ipe_neighbor,itag,icomm,recvrequest_c_sr(irecv_c),ierrmpi)
           if(stagger_grid) then
             irecv_srl=irecv_srl+1
             call mpi_irecv(recvbuffer_srl(ibuf_recv_srl),sizes_srl_recv_total(i^d),mpi_double_precision, &
                            ipe_neighbor,itag,icomm,recvrequest_srl(irecv_srl),ierrmpi)
             ibuf_recv_srl=ibuf_recv_srl+sizes_srl_recv_total(i^d)
           end if
        end if

      end subroutine bc_recv_srl

      !> Receive from fine neighbor
      subroutine bc_recv_restrict

        {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
           inc^db=2*i^db+ic^db\}
           ipe_neighbor=neighbor_child(2,inc^d,igrid)
           if (ipe_neighbor/=mype) then
              irecv_c=irecv_c+1
              itag=(3**^nd+4**^nd)*(igrid-1)+3**^nd+{inc^d*4**(^d-1)+}
              call mpi_irecv(psb(igrid)%w,1,type_recv_r(iib^d,inc^d), &
                             ipe_neighbor,itag,icomm,recvrequest_c_sr(irecv_c),ierrmpi)
              if(stagger_grid) then
                irecv_r=irecv_r+1
                call mpi_irecv(recvbuffer_r(ibuf_recv_r),sizes_r_recv_total(inc^d), &
                               mpi_double_precision,ipe_neighbor,itag, &
                               icomm,recvrequest_r(irecv_r),ierrmpi)
                ibuf_recv_r=ibuf_recv_r+sizes_r_recv_total(inc^d)
              end if
           end if
        {end do\}

      end subroutine bc_recv_restrict

      !> Send to sibling at same refinement level
      subroutine bc_send_srl

        ipe_neighbor=neighbor(2,i^d,igrid)

        if(ipe_neighbor/=mype) then
          ineighbor=neighbor(1,i^d,igrid)
          ipole=neighbor_pole(i^d,igrid)
          if(ipole==0) then
            n_i^d=-i^d;
            isend_c=isend_c+1
            itag=(3**^nd+4**^nd)*(ineighbor-1)+{(n_i^d+1)*3**(^d-1)+}
            call mpi_isend(psb(igrid)%w,1,type_send_srl(iib^d,i^d), &
                           ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
            if(stagger_grid) then
              ibuf_start=ibuf_send_srl
              do idir=1,ndim
                ixs^l=ixs_srl_stg_^l(idir,i^d);
                ibuf_next=ibuf_start+sizes_srl_send_stg(idir,i^d)
                shapes=(/sizes_srl_send_stg(idir,i^d)/)
                sendbuffer_srl(ibuf_start:ibuf_next-1)=&
                  reshape(psb(igrid)%ws(ixs^s,idir),shapes)
                ibuf_start=ibuf_next
              end do
              isend_srl=isend_srl+1
              call mpi_isend(sendbuffer_srl(ibuf_send_srl),sizes_srl_send_total(i^d),&
                mpi_double_precision, ipe_neighbor,itag,icomm,sendrequest_srl(isend_srl),ierrmpi)
              ibuf_send_srl=ibuf_next
            end if
          else
            ixs^l=ixs_srl_^l(iib^d,i^d);
            select case (ipole)
            {case (^d)
               n_i^d=i^d^d%n_i^dd=-i^dd;\}
            end select
            if (isend_buf(ipwbuf)/=0) then
               call mpi_wait(sendrequest_c_sr(isend_buf(ipwbuf)), &
                             sendstatus_c_sr(:,isend_buf(ipwbuf)),ierrmpi)
               deallocate(pwbuf(ipwbuf)%w)
            end if
            allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
            call pole_buffer(pwbuf(ipwbuf)%w,ixs^l,ixs^l,psb(igrid)%w,ixg^ll,ixs^l)
            isend_c=isend_c+1
            isend_buf(ipwbuf)=isend_c
            itag=(3**^nd+4**^nd)*(ineighbor-1)+{(n_i^d+1)*3**(^d-1)+}
            isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
            call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
                           ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
            ipwbuf=1+modulo(ipwbuf,npwbuf)
            if(stagger_grid) then
              ibuf_start=ibuf_send_srl
              do idir=1,ndim
                ixs^l=ixs_srl_stg_^l(idir,i^d);
                ibuf_next=ibuf_start+sizes_srl_send_stg(idir,i^d)
                shapes=(/sizes_srl_send_stg(idir,i^d)/)
                sendbuffer_srl(ibuf_start:ibuf_next-1)=&
                  reshape(psb(igrid)%ws(ixs^s,idir),shapes)
                ibuf_start=ibuf_next
              end do
              isend_srl=isend_srl+1
              call mpi_isend(sendbuffer_srl(ibuf_send_srl),sizes_srl_send_total(i^d),&
                mpi_double_precision, ipe_neighbor,itag,icomm,sendrequest_srl(isend_srl),ierrmpi)
              ibuf_send_srl=ibuf_next
            end if
          end if
        end if

      end subroutine bc_send_srl

      subroutine bc_fill_srl(igrid,i^D,iib^D)
        integer, intent(in) :: igrid,i^D,iib^D
        integer :: ineighbor,ipe_neighbor,ipole,ixS^L,ixR^L,n_i^D,idir

        ipe_neighbor=neighbor(2,i^d,igrid)
        if(ipe_neighbor==mype) then
          ineighbor=neighbor(1,i^d,igrid)
          ipole=neighbor_pole(i^d,igrid)
          if(ipole==0) then
            n_i^d=-i^d;
            ixs^l=ixs_srl_^l(iib^d,i^d);
            ixr^l=ixr_srl_^l(iib^d,n_i^d);
            psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
                psb(igrid)%w(ixs^s,nwhead:nwtail)
            if(stagger_grid) then
              do idir=1,ndim
                ixs^l=ixs_srl_stg_^l(idir,i^d);
                ixr^l=ixr_srl_stg_^l(idir,n_i^d);
                psb(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
              end do
            end if
          else
            ixs^l=ixs_srl_^l(iib^d,i^d);
            select case (ipole)
            {case (^d)
              n_i^d=i^d^d%n_i^dd=-i^dd;\}
            end select
            ixr^l=ixr_srl_^l(iib^d,n_i^d);
            call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^l,psb(igrid)%w,ixg^ll,ixs^l,ipole)
            if(stagger_grid) then
              do idir=1,ndim
                ixs^l=ixs_srl_stg_^l(idir,i^d);
                ixr^l=ixr_srl_stg_^l(idir,n_i^d);
                call pole_copy_stg(psb(ineighbor)%ws,ixr^l,psb(igrid)%ws,ixs^l,idir,ipole)
              end do
            end if
          end if
        end if

      end subroutine bc_fill_srl

      subroutine fill_coarse_boundary(igrid,i^D)
        integer, intent(in) :: igrid,i^D

        integer :: idims,iside,k^L,ixB^L,ii^D

        if(phyboundblock(igrid).and..not.stagger_grid.and.bcphys) then
          ! to use block in physical boundary setup for coarse representative
          block=>psc(igrid)
          ! filling physical boundary ghost cells of a coarser representative block for
          ! sending swap region with width of nghostcells to its coarser neighbor
          do idims=1,ndim
             ! to avoid using as yet unknown corner info in more than 1D, we
             ! fill only interior mesh ranges of the ghost cell ranges at first,
             ! and progressively enlarge the ranges to include corners later
             {kmin^d=merge(0, 1, idims==^d)
             kmax^d=merge(0, 1, idims==^d)
             ixbmin^d=ixcogmin^d+kmin^d*nghostcells
             ixbmax^d=ixcogmax^d-kmax^d*nghostcells\}
             {^iftwod
             if(idims > 1 .and. neighbor_type(-1,0,igrid)==neighbor_boundary) ixbmin1=ixcogmin1
             if(idims > 1 .and. neighbor_type( 1,0,igrid)==neighbor_boundary) ixbmax1=ixcogmax1}
             {^ifthreed
             if(idims > 1 .and. neighbor_type(-1,0,0,igrid)==neighbor_boundary) ixbmin1=ixcogmin1
             if(idims > 1 .and. neighbor_type( 1,0,0,igrid)==neighbor_boundary) ixbmax1=ixcogmax1
             if(idims > 2 .and. neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixbmin2=ixcogmin2
             if(idims > 2 .and. neighbor_type(0, 1,0,igrid)==neighbor_boundary) ixbmax2=ixcogmax2}
             {if(i^d==-1) then
               ixbmin^d=ixcogmin^d+nghostcells
               ixbmax^d=ixcogmin^d+2*nghostcells-1
             else if(i^d==1) then
               ixbmin^d=ixcogmax^d-2*nghostcells+1
               ixbmax^d=ixcogmax^d-nghostcells
             end if\}
             do iside=1,2
                ii^d=kr(^d,idims)*(2*iside-3);
                if ({abs(i^d)==1.and.abs(ii^d)==1|.or.}) cycle
                if (neighbor_type(ii^d,igrid)/=neighbor_boundary) cycle
                call bc_phys(iside,idims,time,0.d0,psc(igrid),ixcog^l,ixb^l)
             end do
          end do
        end if

      end subroutine fill_coarse_boundary

      !> Send to coarser neighbor
      subroutine bc_send_restrict

        ipe_neighbor=neighbor(2,i^d,igrid)
        if(ipe_neighbor/=mype) then
          ic^d=1+modulo(node(pig^d_,igrid)-1,2);
          if({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return
          ineighbor=neighbor(1,i^d,igrid)
          ipole=neighbor_pole(i^d,igrid)
          if(ipole==0) then
            n_inc^d=-2*i^d+ic^d;
            isend_c=isend_c+1
            itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
            call mpi_isend(psc(igrid)%w,1,type_send_r(iib^d,i^d), &
                           ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
            if(stagger_grid) then 
              ibuf_start=ibuf_send_r
              do idir=1,ndim
                ixs^l=ixs_r_stg_^l(idir,i^d);
                ibuf_next=ibuf_start+sizes_r_send_stg(idir,i^d)
                shapes=(/sizes_r_send_stg(idir,i^d)/)
                sendbuffer_r(ibuf_start:ibuf_next-1)=&
                  reshape(psc(igrid)%ws(ixs^s,idir),shapes)
                ibuf_start=ibuf_next
              end do
              isend_r=isend_r+1
              call mpi_isend(sendbuffer_r(ibuf_send_r),sizes_r_send_total(i^d),&
                             mpi_double_precision,ipe_neighbor,itag, &
                             icomm,sendrequest_r(isend_r),ierrmpi)
              ibuf_send_r=ibuf_next
            end if
          else
            ixs^l=ixs_r_^l(iib^d,i^d);
            select case (ipole)
            {case (^d)
               n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
            end select
            if(isend_buf(ipwbuf)/=0) then
              call mpi_wait(sendrequest_c_sr(isend_buf(ipwbuf)), &
                            sendstatus_c_sr(:,isend_buf(ipwbuf)),ierrmpi)
              deallocate(pwbuf(ipwbuf)%w)
            end if
            allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
            call pole_buffer(pwbuf(ipwbuf)%w,ixs^l,ixs^l,psc(igrid)%w,ixcog^l,ixs^l)
            isend_c=isend_c+1
            isend_buf(ipwbuf)=isend_c
            itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
            isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
            call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
                           ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
            ipwbuf=1+modulo(ipwbuf,npwbuf)
            if(stagger_grid) then 
              ibuf_start=ibuf_send_r
              do idir=1,ndim
                ixs^l=ixs_r_stg_^l(idir,i^d);
                ibuf_next=ibuf_start+sizes_r_send_stg(idir,i^d)
                shapes=(/sizes_r_send_stg(idir,i^d)/)
                sendbuffer_r(ibuf_start:ibuf_next-1)=&
                  reshape(psc(igrid)%ws(ixs^s,idir),shapes)
                ibuf_start=ibuf_next
              end do
              isend_r=isend_r+1
              call mpi_isend(sendbuffer_r(ibuf_send_r),sizes_r_send_total(i^d),&
                             mpi_double_precision,ipe_neighbor,itag, &
                             icomm,sendrequest_r(isend_r),ierrmpi)
              ibuf_send_r=ibuf_next
            end if
          end if
        end if

      end subroutine bc_send_restrict

      !> fill coarser neighbor's ghost cells
      subroutine bc_fill_restrict(igrid,i^D,iib^D)
        integer, intent(in) :: igrid,i^D,iib^D

        integer :: ic^D,n_inc^D,ixS^L,ixR^L,ipe_neighbor,ineighbor,ipole,idir

        ipe_neighbor=neighbor(2,i^d,igrid)
        if(ipe_neighbor==mype) then
          ic^d=1+modulo(node(pig^d_,igrid)-1,2);
          if({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return
          ineighbor=neighbor(1,i^d,igrid)
          ipole=neighbor_pole(i^d,igrid)
          if(ipole==0) then
            n_inc^d=-2*i^d+ic^d;
            ixs^l=ixs_r_^l(iib^d,i^d);
            ixr^l=ixr_r_^l(iib^d,n_inc^d);
            psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
                psc(igrid)%w(ixs^s,nwhead:nwtail)
            if(stagger_grid) then
              do idir=1,ndim
                 ixs^l=ixs_r_stg_^l(idir,i^d);
                 ixr^l=ixr_r_stg_^l(idir,n_inc^d);
                 psb(ineighbor)%ws(ixr^s,idir)=psc(igrid)%ws(ixs^s,idir)
              end do
            end if
          else
            ixs^l=ixs_r_^l(iib^d,i^d);
            select case (ipole)
            {case (^d)
              n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
            end select
            ixr^l=ixr_r_^l(iib^d,n_inc^d);
            call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^l,psc(igrid)%w,ixcog^l,ixs^l,ipole)
            if(stagger_grid) then
              do idir=1,ndim
                ixs^l=ixs_r_stg_^l(idir,i^d);
                ixr^l=ixr_r_stg_^l(idir,n_inc^d);
                !! Fill ghost cells
                call pole_copy_stg(psb(ineighbor)%ws,ixr^l,psc(igrid)%ws,ixs^l,idir,ipole)
              end do
            end if
          end if
        end if

      end subroutine bc_fill_restrict

      !> fill siblings ghost cells with received data
      subroutine bc_fill_srl_stg
        double precision :: tmp(ixgs^t)
        integer :: ixS^L,ixR^L,n_i^D,ixSsync^L,ixRsync^L
        integer :: idir

        ipe_neighbor=neighbor(2,i^d,igrid)
        if(ipe_neighbor/=mype) then
          ineighbor=neighbor(1,i^d,igrid)
          ipole=neighbor_pole(i^d,igrid)

        !! Now the special treatment of the pole is done here, at the receive step
          if (ipole==0) then    
            ixr^l=ixr_srl_^l(iib^d,i^d);
            !! Unpack the buffer and fill the ghost cells
            n_i^d=-i^d;
            do idir=1,ndim
              ixs^l=ixs_srl_stg_^l(idir,n_i^d);
              ixr^l=ixr_srl_stg_^l(idir,i^d);
              ibuf_next=ibuf_recv_srl+sizes_srl_recv_stg(idir,i^d)
              tmp(ixs^s) = reshape(source=recvbuffer_srl(ibuf_recv_srl:ibuf_next-1),shape=shape(psb(igrid)%ws(ixs^s,idir)))       
              psb(igrid)%ws(ixr^s,idir) = tmp(ixs^s)
              ibuf_recv_srl=ibuf_next
            end do
          else ! There is a pole
            select case (ipole)
            {case (^d)
               n_i^d=i^d^d%n_i^dd=-i^dd;\}
            end select
            pole_buf%ws=zero
            do idir=1,ndim
             ixr^l=ixr_srl_stg_^l(idir,i^d);
             ixs^l=ixs_srl_stg_^l(idir,n_i^d);
             ibuf_next=ibuf_recv_srl+sizes_srl_recv_stg(idir,i^d)
             pole_buf%ws(ixs^s,idir)=reshape(source=recvbuffer_srl(ibuf_recv_srl:ibuf_next-1),&
               shape=shape(psb(igrid)%ws(ixs^s,idir)))
             ibuf_recv_srl=ibuf_next
             call pole_copy_stg(psb(igrid)%ws,ixr^l,pole_buf%ws,ixs^l,idir,ipole)
            end do
          end if
        end if

      end subroutine bc_fill_srl_stg

      subroutine indices_for_syncing(idir,i^D,ixR^L,ixS^L,ixRsync^L,ixSsync^L)
        integer, intent(in)       :: i^D,idir
        integer, intent(inout)    :: ixR^L,ixS^L
        integer, intent(out)      :: ixRsync^L,ixSsync^L
      
        ixrsync^l=ixr^l;
        ixssync^l=ixs^l;
        
        {
        if (i^d == -1 .and. idir == ^d) then
           ixrsyncmin^d = ixrmax^d
           ixrsyncmax^d = ixrmax^d
           ixssyncmin^d = ixsmax^d
           ixssyncmax^d = ixsmax^d
           ixrmax^d = ixrmax^d - 1
           ixsmax^d = ixsmax^d - 1
        else if (i^d == 1 .and. idir == ^d) then
           ixrsyncmin^d = ixrmin^d
           ixrsyncmax^d = ixrmin^d
           ixssyncmin^d = ixsmin^d
           ixssyncmax^d = ixsmin^d
           ixrmin^d = ixrmin^d + 1
           ixsmin^d = ixsmin^d + 1
        end if
        \}

      end subroutine indices_for_syncing

      !> fill restricted ghost cells after receipt
      subroutine bc_fill_restrict_stg

        ipole=neighbor_pole(i^d,igrid)
        if (ipole==0) then
          ! Loop over the children ic^D to and their neighbors inc^D
          {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
             inc^db=2*i^db+ic^db\}
             ipe_neighbor=neighbor_child(2,inc^d,igrid)
             if(ipe_neighbor/=mype) then
               ineighbor=neighbor_child(1,inc^d,igrid)
               n_i^d=-i^d;
               !! Unpack the buffer and fill the ghost cells
               do idir=1,ndim
                 ixr^l=ixr_r_stg_^l(idir,inc^d);
                 ibuf_next=ibuf_recv_r+sizes_r_recv_stg(idir,inc^d)
                 psb(igrid)%ws(ixr^s,idir)=reshape(source=recvbuffer_r(ibuf_recv_r:ibuf_next-1),&
                       shape=shape(psb(igrid)%ws(ixr^s,idir)))
                 ibuf_recv_r=ibuf_next
               end do
             end if
          {end do\}
        else !! There is a pole
          {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
             inc^db=2*i^db+ic^db\}
             ipe_neighbor=neighbor_child(2,inc^d,igrid)
             if(ipe_neighbor/=mype) then
               ineighbor=neighbor_child(1,inc^d,igrid)
               select case(ipole)
              {case (^d)
                 n_i^d=i^d^d%n_i^dd=-i^dd;\}
               end select
               ixr^l=ixr_r_^l(iib^d,inc^d);
               !! Unpack the buffer and fill an auxiliary array
               pole_buf%ws=zero
               do idir=1,ndim
                 ixs^l=ixs_r_stg_^l(idir,n_i^d);
                 ixr^l=ixr_r_stg_^l(idir,inc^d);
                 ibuf_next=ibuf_recv_r+sizes_r_recv_stg(idir,inc^d)
                 pole_buf%ws(ixr^s,idir)=reshape(source=recvbuffer_r(ibuf_recv_r:ibuf_next-1),&
                   shape=shape(psb(igrid)%ws(ixr^s,idir)))
                 call pole_copy_stg(psb(igrid)%ws,ixr^l,pole_buf%ws,ixr^l,idir,ipole)
                 ibuf_recv_r=ibuf_next
               end do
             end if
          {end do\}
        end if

      end subroutine bc_fill_restrict_stg

      !> Receive from coarse neighbor
      subroutine bc_recv_prolong

        ic^d=1+modulo(node(pig^d_,igrid)-1,2);
        if ({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return

        ipe_neighbor=neighbor(2,i^d,igrid)
        if (ipe_neighbor/=mype) then
           irecv_c=irecv_c+1
           inc^d=ic^d+i^d;
           itag=(3**^nd+4**^nd)*(igrid-1)+3**^nd+{inc^d*4**(^d-1)+}
           call mpi_irecv(psc(igrid)%w,1,type_recv_p(iib^d,inc^d), &
                          ipe_neighbor,itag,icomm,recvrequest_c_p(irecv_c),ierrmpi)  
           if(stagger_grid) then
             irecv_p=irecv_p+1
             call mpi_irecv(recvbuffer_p(ibuf_recv_p),sizes_p_recv_total(inc^d),& 
                            mpi_double_precision,ipe_neighbor,itag,&
                            icomm,recvrequest_p(irecv_p),ierrmpi)
             ibuf_recv_p=ibuf_recv_p+sizes_p_recv_total(inc^d)
           end if
        end if

      end subroutine bc_recv_prolong

      !> Send to finer neighbor
      subroutine bc_send_prolong
        integer :: ii^D

        ipole=neighbor_pole(i^d,igrid)

        {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
           inc^db=2*i^db+ic^db\}
           ipe_neighbor=neighbor_child(2,inc^d,igrid)
           if(ipe_neighbor/=mype) then
             ixs^l=ixs_p_^l(iib^d,inc^d);
             ineighbor=neighbor_child(1,inc^d,igrid)
             if(ipole==0) then
               n_i^d=-i^d;
               n_inc^d=ic^d+n_i^d;
               isend_c=isend_c+1
               itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
               call mpi_isend(psb(igrid)%w,1,type_send_p(iib^d,inc^d), &
                              ipe_neighbor,itag,icomm,sendrequest_c_p(isend_c),ierrmpi)
               if(stagger_grid) then 
                 ibuf_start=ibuf_send_p
                 do idir=1,ndim
                   ixs^l=ixs_p_stg_^l(idir,inc^d);
                   ibuf_next=ibuf_start+sizes_p_send_stg(idir,inc^d)
                   shapes=(/sizes_p_send_stg(idir,inc^d)/)
                   sendbuffer_p(ibuf_start:ibuf_next-1)=&
                     reshape(psb(igrid)%ws(ixs^s,idir),shapes)   
                   ibuf_start=ibuf_next
                 end do
                 isend_p=isend_p+1
                 call mpi_isend(sendbuffer_p(ibuf_send_p),sizes_p_send_total(inc^d),&
                                mpi_double_precision,ipe_neighbor,itag, &
                                icomm,sendrequest_p(isend_p),ierrmpi)
                 ibuf_send_p=ibuf_next
               end if
             else
               select case (ipole)
               {case (^d)
                 n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
               end select
               if(isend_buf(ipwbuf)/=0) then
                 call mpi_wait(sendrequest_c_p(isend_buf(ipwbuf)), &
                               sendstatus_c_p(:,isend_buf(ipwbuf)),ierrmpi)
                 deallocate(pwbuf(ipwbuf)%w)
               end if
               allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
               call pole_buffer(pwbuf(ipwbuf)%w,ixs^l,ixs^l,psb(igrid)%w,ixg^ll,ixs^l)
               isend_c=isend_c+1
               isend_buf(ipwbuf)=isend_c
               itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
               isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
               call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
                              ipe_neighbor,itag,icomm,sendrequest_c_p(isend_c),ierrmpi)
               ipwbuf=1+modulo(ipwbuf,npwbuf)
               if(stagger_grid) then 
                 ibuf_start=ibuf_send_p
                 do idir=1,ndim
                   ixs^l=ixs_p_stg_^l(idir,inc^d);
                   ibuf_next=ibuf_start+sizes_p_send_stg(idir,inc^d)
                   shapes=(/sizes_p_send_stg(idir,inc^d)/)
                   sendbuffer_p(ibuf_start:ibuf_next-1)=&
                     reshape(psb(igrid)%ws(ixs^s,idir),shapes)   
                   ibuf_start=ibuf_next
                 end do
                 isend_p=isend_p+1
                 call mpi_isend(sendbuffer_p(ibuf_send_p),sizes_p_send_total(inc^d),&
                                mpi_double_precision,ipe_neighbor,itag, &
                                icomm,sendrequest_p(isend_p),ierrmpi)
                 ibuf_send_p=ibuf_next
               end if
             end if
           end if
        {end do\}

      end subroutine bc_send_prolong

      !> Send to finer neighbor
      subroutine bc_fill_prolong(igrid,i^D,iib^D)
        integer, intent(in) :: igrid,i^D,iib^D

        integer :: ipe_neighbor,ineighbor,ixS^L,ixR^L,ic^D,inc^D,ipole,idir

        ipole=neighbor_pole(i^d,igrid)

        if(ipole==0) then
          {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
             inc^db=2*i^db+ic^db\}
             ipe_neighbor=neighbor_child(2,inc^d,igrid)
             if(ipe_neighbor==mype) then
               ixs^l=ixs_p_^l(iib^d,inc^d);
               ineighbor=neighbor_child(1,inc^d,igrid)
               ipole=neighbor_pole(i^d,igrid)
               n_i^d=-i^d;
               n_inc^d=ic^d+n_i^d;
               ixr^l=ixr_p_^l(iib^d,n_inc^d);
               psc(ineighbor)%w(ixr^s,nwhead:nwtail) &
                  =psb(igrid)%w(ixs^s,nwhead:nwtail)
               if(stagger_grid) then
                 do idir=1,ndim
                   ixs^l=ixs_p_stg_^l(idir,inc^d);
                   ixr^l=ixr_p_stg_^l(idir,n_inc^d);
                   psc(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
                 end do
               end if
             end if
          {end do\}
        else
          {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
             inc^db=2*i^db+ic^db\}
             ipe_neighbor=neighbor_child(2,inc^d,igrid)
             if(ipe_neighbor==mype) then
               ixs^l=ixs_p_^l(iib^d,inc^d);
               ineighbor=neighbor_child(1,inc^d,igrid)
               ipole=neighbor_pole(i^d,igrid)
               select case (ipole)
               {case (^d)
                  n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
               end select
               ixr^l=ixr_p_^l(iib^d,n_inc^d);
               call pole_copy(psc(ineighbor)%w,ixcog^l,ixr^l,psb(igrid)%w,ixg^ll,ixs^l,ipole)
               if(stagger_grid) then
                 do idir=1,ndim
                   ixs^l=ixs_p_stg_^l(idir,inc^d);
                   ixr^l=ixr_p_stg_^l(idir,n_inc^d);
                   call pole_copy_stg(psc(ineighbor)%ws,ixr^l,psb(igrid)%ws,ixs^l,idir,ipole)
                 end do
               end if
             end if
          {end do\}
        end if
      end subroutine bc_fill_prolong

      subroutine gc_prolong(igrid)
        integer, intent(in) :: igrid

        integer :: iib^D,i^D,idims,iside
        logical,dimension(-1:1^D&) :: NeedProlong

        ^d&iib^d=idphyb(^d,igrid);
        needprolong=.false.
        {do i^db=-1,1\}
           if (skip_direction([ i^d ])) cycle
           if (neighbor_type(i^d,igrid)==neighbor_coarse) then
             call bc_prolong(igrid,i^d,iib^d)
             needprolong(i^d)=.true.
           end if
        {end do\}
        if(stagger_grid) then
          ! Ghost cell prolongation for staggered variables
          ! must be done in a specific order.
          ! First the first neighbours, which have 2 indices=0 in 3D
          ! or one index=0 in 2D
          do idims=1,ndim
            i^d=0;
            select case(idims)
           {case(^d)
              do i^d=-1,1,2
                if (needprolong(i^dd)) call bc_prolong_stg(igrid,i^dd,iib^dd,needprolong)
              end do
            \}
            end select
          end do
          ! Then the second neighbours which have 1 index=0 in 3D
          ! (Only in 3D)
          {^ifthreed
          i1=0;
          do i2=-1,1,2
            do i3=-1,1,2
              if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
            end do
          end do
          i2=0;
          do i3=-1,1,2
            do i1=-1,1,2
              if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
            end do
          end do
          i3=0;
          do i1=-1,1,2
            do i2=-1,1,2
              if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
            end do
          end do
          }
          ! Finally, the corners, that have no index=0
         {do i^d=-1,1,2\}
            if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
         {end do\}
        end if
      end subroutine gc_prolong

      !> fill coarser representative with data from coarser neighbors
      subroutine bc_fill_prolong_stg
        ic^d=1+modulo(node(pig^d_,igrid)-1,2);
        if ({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return

        ipe_neighbor=neighbor(2,i^d,igrid)
        if(ipe_neighbor/=mype) then
          ineighbor=neighbor(1,i^d,igrid)
          ipole=neighbor_pole(i^d,igrid)

          if (ipole==0) then   !! There is no pole 
            inc^d=ic^d+i^d;
            ixr^l=ixr_p_^l(iib^d,inc^d);
            do idir=1,ndim
              ixr^l=ixr_p_stg_^l(idir,inc^d);
              ibuf_next=ibuf_recv_p+sizes_p_recv_stg(idir,inc^d)
              psc(igrid)%ws(ixr^s,idir)=reshape(source=recvbuffer_p(ibuf_recv_p:ibuf_next-1),&
                    shape=shape(psc(igrid)%ws(ixr^s,idir)))
              ibuf_recv_p=ibuf_next
            end do
          else !! There is a pole
            inc^d=ic^d+i^d;
            select case (ipole)
            {case (^d)
               n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
            end select
            !! Unpack the buffer and fill an auxiliary array
            pole_buf%ws=zero
            do idir=1,ndim
              ixr^l=ixr_p_stg_^l(idir,inc^d);
              ibuf_next=ibuf_recv_p+sizes_p_recv_stg(idir,inc^d)
              pole_buf%ws(ixr^s,idir)=reshape(source=recvbuffer_p(ibuf_recv_p:ibuf_next-1),&
                shape=shape(psc(igrid)%ws(ixr^s,idir)))
              call pole_copy_stg(psc(igrid)%ws,ixr^l,pole_buf%ws,ixr^l,idir,ipole)
              ibuf_recv_p=ibuf_next
            end do
          end if
        end if

      end subroutine bc_fill_prolong_stg

      !> do prolongation for fine blocks after receipt data from coarse neighbors
      subroutine bc_prolong(igrid,i^D,iib^D)
        use mod_physics, only: phys_to_primitive, phys_to_conserved

        integer :: i^D,iib^D,igrid
        integer :: ixFi^L,ixCo^L,ii^D, idims,iside,ixB^L
        double precision :: dxFi^D, dxCo^D, xFimin^D, xComin^D, invdxCo^D

        ixfi^l=ixr_srl_^l(iib^d,i^d);
        dxfi^d=rnode(rpdx^d_,igrid);
        dxco^d=two*dxfi^d;
        invdxco^d=1.d0/dxco^d;

        ! compute the enlarged grid lower left corner coordinates
        ! these are true coordinates for an equidistant grid, 
        ! but we can temporarily also use them for getting indices 
        ! in stretched grids
        xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
        xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;

        if(stagger_grid.and.phyboundblock(igrid).and.bcphys) then
          block=>psc(igrid)
          do idims=1,ndim
            ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
            ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
            {^ifthreed
            ! avoid using undetermined ghost cells at physical boundary edges
            if(idims == 1) then
              if(neighbor_type(-1,0,0,igrid)==neighbor_boundary .or. &
                 neighbor_type(1,0,0,igrid)==neighbor_boundary) then
                if(neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixcomin2=ixcommin2
                if(neighbor_type(0,0,-1,igrid)==neighbor_boundary) ixcomin3=ixcommin3
                if(neighbor_type(0,1,0,igrid)==neighbor_boundary) ixcomax2=ixcommax2
                if(neighbor_type(0,0,1,igrid)==neighbor_boundary) ixcomax3=ixcommax3
              end if
            else if(idims == 2) then
              if(neighbor_type(0,-1,0,igrid)==neighbor_boundary .or. &
                 neighbor_type(0,1,0,igrid)==neighbor_boundary) then
                if(neighbor_type(0,0,-1,igrid)==neighbor_boundary) ixcomin3=ixcommin3
                if(neighbor_type(0,0,1,igrid)==neighbor_boundary) ixcomax3=ixcommax3
              end if
            end if
            }
            do iside=1,2
              ii^d=kr(^d,idims)*(2*iside-3);
              if(neighbor_type(ii^d,igrid)/=neighbor_boundary) cycle
              if(( {(iside==1.and.idims==^d.and.ixcomin^d<ixcogmin^d+nghostcells)|.or.} ) &
               .or.( {(iside==2.and.idims==^d.and.ixcomax^d>ixcogmax^d-nghostcells)|.or. })) then
                {ixbmin^d=merge(ixcogmin^d,ixcomin^d,idims==^d);}
                {ixbmax^d=merge(ixcogmax^d,ixcomax^d,idims==^d);}
                call bc_phys(iside,idims,time,0.d0,psc(igrid),ixcog^l,ixb^l)
              end if
            end do
          end do
        end if

        if(prolongprimitive) then
           ! following line again assumes equidistant grid, but 
           ! just computes indices, so also ok for stretched case
           ! reason for +1-1 and +1+1: the coarse representation has 
           ! also nghostcells at each side. During
           ! prolongation, we need cells to left and right, hence -1/+1
           ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
           ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
           call phys_to_primitive(ixcog^l,ixco^l,&
             psc(igrid)%w,psc(igrid)%x)
        endif

        if(ghost_copy) then
           call interpolation_copy(igrid,ixfi^l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
        else
           call interpolation_linear(igrid,ixfi^l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
        end if

        if(prolongprimitive) call phys_to_conserved(ixcog^l,ixco^l,&
             psc(igrid)%w,psc(igrid)%x)

      end subroutine bc_prolong

      subroutine bc_prolong_stg(igrid,i^D,iib^D,NeedProlong)
        use mod_amr_fct
        integer                    :: igrid,i^D,iib^D
        logical,dimension(-1:1^D&) :: NeedProlong
        logical                    :: fine_^Lin
        integer                    :: ixFi^L,ixCo^L
        double precision           :: dxFi^D,dxCo^D,xFimin^D,xComin^D,invdxCo^D
        ! Check what is already at the desired level
        fine_^lin=.false.;
        {
        if(i^d>-1) fine_min^din=(.not.needprolong(i^dd-kr(^d,^dd)).and.neighbor_type(i^dd-kr(^d,^dd),igrid)/=1)
        if(i^d<1)  fine_max^din=(.not.needprolong(i^dd+kr(^d,^dd)).and.neighbor_type(i^dd+kr(^d,^dd),igrid)/=1)
        \}

        ixfi^l=ixr_srl_^l(iib^d,i^d);

        dxfi^d=rnode(rpdx^d_,igrid);
        dxco^d=two*dxfi^d;
        invdxco^d=1.d0/dxco^d;

        xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
        xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;

        ! moved the physical boundary filling here, to only fill the
        ! part needed

        ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
        ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;

        if(prolongprimitive) call phys_to_primitive(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)

        call prolong_2nd_stg(psc(igrid),psb(igrid),ixco^l,ixfi^l,dxco^d,xcomin^d,dxfi^d,xfimin^d,.true.,fine_^lin)

        if(prolongprimitive) call phys_to_conserved(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)

        ! The current region has already been refined, so it does not need to be prolonged again
        needprolong(i^d)=.false. 

      end subroutine bc_prolong_stg

      subroutine interpolation_linear(igrid,ixFi^L,dxFi^D,xFimin^D, &
                                      dxCo^D,invdxCo^D,xComin^D)
        use mod_physics, only: phys_to_conserved
        integer, intent(in) :: igrid, ixFi^L
        double precision, intent(in) :: dxFi^D, xFimin^D,dxCo^D, invdxCo^D, xComin^D

        integer :: ixCo^D, jxCo^D, hxCo^D, ixFi^D, ix^D, iw, idims, nwmin,nwmax
        double precision :: xCo^D, xFi^D, eta^D
        double precision :: slopeL, slopeR, slopeC, signC, signR
        double precision :: slope(1:nw,ndim)
        !!double precision :: local_invdxCo^D
        double precision :: signedfactorhalf^D
        !integer :: ixshift^D, icase

        !icase=mod(nghostcells,2)

        if(prolongprimitive) then
          nwmin=1
          nwmax=nw
        else
          nwmin=nwhead
          nwmax=nwtail
        end if

        {do ixfi^db = ixfi^lim^db
           ! cell-centered coordinates of fine grid point
           ! here we temporarily use an equidistant grid
           xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
        
           ! indices of coarse cell which contains the fine cell
           ! since we computed lower left corner earlier 
           ! in equidistant fashion: also ok for stretched case
           ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1
        
           ! cell-centered coordinates of coarse grid point
           ! here we temporarily use an equidistant grid
           xco^db=xcomin^db+(dble(ixco^db)-half)*dxco^db \}

           !if(.not.slab) then
           !   ^D&local_invdxCo^D=1.d0/psc(igrid)%dx({ixCo^DD},^D)\
           !endif

           if(slab_uniform) then
             ! actual cell-centered coordinates of fine grid point
             !!^D&xFi^D=block%x({ixFi^DD},^D)\
             ! actual cell-centered coordinates of coarse grid point
             !!^D&xCo^D=psc(igrid)%x({ixCo^DD},^D)\
             ! normalized distance between fine/coarse cell center
             ! in coarse cell: ranges from -0.5 to 0.5 in each direction
             ! (origin is coarse cell center)
             ! this is essentially +1/4 or -1/4 on cartesian mesh
             eta^d=(xfi^d-xco^d)*invdxco^d;
           else
             !select case(icase)
             ! case(0)
             !{! here we assume an even number of ghostcells!!!
             !ixshift^D=2*(mod(ixFi^D,2)-1)+1
             !if(ixshift^D>0.0d0)then
             !   ! oneven fine grid points
             !   eta^D=-0.5d0*(one-block%dvolume(ixFi^DD) &
             !     /sum(block%dvolume(ixFi^D:ixFi^D+1^D%ixFi^DD))) 
             !else
             !   ! even fine grid points
             !   eta^D=+0.5d0*(one-block%dvolume(ixFi^DD) &
             !     /sum(block%dvolume(ixFi^D-1:ixFi^D^D%ixFi^DD))) 
             !endif\}
             ! case(1)
             !{! here we assume an odd number of ghostcells!!!
             !ixshift^D=2*(mod(ixFi^D,2)-1)+1
             !if(ixshift^D>0.0d0)then
             !   ! oneven fine grid points
             !   eta^D=+0.5d0*(one-block%dvolume(ixFi^DD) &
             !     /sum(block%dvolume(ixFi^D-1:ixFi^D^D%ixFi^DD))) 
             !else
             !   ! even fine grid points
             !   eta^D=-0.5d0*(one-block%dvolume(ixFi^DD) &
             !     /sum(block%dvolume(ixFi^D:ixFi^D+1^D%ixFi^DD))) 
             !endif\}
             ! case default
             !  call mpistop("no such case")
             !end select
             ! the different cases for even/uneven number of ghost cells 
             ! are automatically handled using the relative index to ixMlo
             ! as well as the pseudo-coordinates xFi and xCo 
             ! these latter differ from actual cell centers when stretching is used
             ix^d=2*int((ixfi^d+ixmlo^d)/2)-ixmlo^d;
             {if(xfi^d>xco^d) then
                signedfactorhalf^d=0.5d0
              else
                signedfactorhalf^d=-0.5d0
              end if
              eta^d=signedfactorhalf^d*(one-psb(igrid)%dvolume(ixfi^dd) &
                   /sum(psb(igrid)%dvolume(ix^d:ix^d+1^d%ixFi^dd))) \}
             !{eta^D=(xFi^D-xCo^D)*invdxCo^D &
             !      *two*(one-block%dvolume(ixFi^DD) &
             !      /sum(block%dvolume(ix^D:ix^D+1^D%ixFi^DD))) \}
           end if
        
           do idims=1,ndim
              hxco^d=ixco^d-kr(^d,idims)\
              jxco^d=ixco^d+kr(^d,idims)\
        
              do iw=nwmin,nwmax
                 slopel=psc(igrid)%w(ixco^d,iw)-psc(igrid)%w(hxco^d,iw)
                 sloper=psc(igrid)%w(jxco^d,iw)-psc(igrid)%w(ixco^d,iw)
                 slopec=half*(sloper+slopel)
        
                 ! get limited slope
                 signr=sign(one,sloper)
                 signc=sign(one,slopec)
                 !select case(prolong_limiter)
                 !case(1)
                 !  ! unlimit
                 !  slope(iw,idims)=slopeC
                 !case(2)
                 !  ! minmod
                 !  slope(iw,idims)=signR*max(zero,min(dabs(slopeR), &
                 !                                    signR*slopeL))
                 !case(3)
                 !  ! woodward
                 !  slope(iw,idims)=two*signR*max(zero,min(dabs(slopeR), &
                 !                     signR*slopeL,signR*half*slopeC))
                 !case(4)
                 !  ! koren
                 !  slope(iw,idims)=signR*max(zero,min(two*signR*slopeL, &
                 !   (dabs(slopeR)+two*slopeL*signR)*third,two*dabs(slopeR)))
                 !case default
                   slope(iw,idims)=signc*max(zero,min(dabs(slopec), &
                                     signc*slopel,signc*sloper))
                 !end select
              end do
           end do
        
           ! Interpolate from coarse cell using limited slopes
           psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)+&
             {(slope(nwmin:nwmax,^d)*eta^d)+}
        
        {end do\}
        
        if(prolongprimitive) call phys_to_conserved(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)
      
      end subroutine interpolation_linear

      subroutine interpolation_copy(igrid, ixFi^L,dxFi^D,xFimin^D, &
                                    dxCo^D,invdxCo^D,xComin^D)
        use mod_physics, only: phys_to_conserved
        integer, intent(in) :: igrid, ixFi^L
        double precision, intent(in) :: dxFi^D, xFimin^D,dxCo^D, invdxCo^D, xComin^D

        integer :: ixCo^D, ixFi^D, nwmin,nwmax
        double precision :: xFi^D

        if(prolongprimitive) then
          nwmin=1
          nwmax=nw
        else
          nwmin=nwhead
          nwmax=nwtail
        end if

        {do ixfi^db = ixfi^lim^db
           ! cell-centered coordinates of fine grid point
           xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
        
           ! indices of coarse cell which contains the fine cell
           ! note: this also works for stretched grids
           ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1\}
        
           ! Copy from coarse cell
           psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)
        
        {end do\}
        
        if(prolongprimitive) call phys_to_conserved(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)
      
      end subroutine interpolation_copy

      subroutine pole_copy(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,ipole)
      
        integer, intent(in) :: ixIR^L,ixR^L,ixIS^L,ixS^L,ipole
        double precision :: wrecv(ixir^s,1:nw), wsend(ixis^s,1:nw)

        integer :: iw, iside, iB

        select case (ipole)
        {case (^d)
           iside=int((i^d+3)/2)
           ib=2*(^d-1)+iside
           do iw=nwhead,nwtail
             select case (typeboundary(iw,ib))
             case ("symm")
               wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
             case ("asymm")
               wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
             case default
               call mpistop("Pole boundary condition should be symm or asymm")
             end select
           end do \}
        end select
      
      end subroutine pole_copy

      subroutine pole_copy_stg(wrecv,ixR^L,wsend,ixS^L,idirs,ipole)
      
        integer, intent(in) :: ixR^L,ixS^L,idirs,ipole

        double precision :: wrecv(ixgs^t,1:nws), wsend(ixgs^t,1:nws)
        integer :: iB, iside

        select case (ipole)
        {case (^d)
           iside=int((i^d+3)/2)
           ib=2*(^d-1)+iside
           select case (typeboundary(iw_mag(idirs),ib))
           case ("symm")
             wrecv(ixr^s,idirs) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
           case ("asymm")
             wrecv(ixr^s,idirs) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
           case default
             call mpistop("Pole boundary condition should be symm or asymm")
           end select
         \}
        end select

      end subroutine pole_copy_stg

      subroutine pole_buffer(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L)
      
        integer, intent(in) :: ixIR^L,ixR^L,ixIS^L,ixS^L
        double precision :: wrecv(ixir^s,nwhead:nwtail), wsend(ixis^s,1:nw)

        integer :: iw, iside, iB

        select case (ipole)
        {case (^d)
           iside=int((i^d+3)/2)
           ib=2*(^d-1)+iside
           do iw=nwhead,nwtail
             select case (typeboundary(iw,ib))
             case ("symm")
               wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
             case ("asymm")
               wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
             case default
               call mpistop("Pole boundary condition should be symm or asymm")
             end select
           end do \}
        end select
      
      end subroutine pole_buffer

      subroutine fix_auxiliary(igrid)
        use mod_physics, only: phys_get_aux
        integer, intent(in) :: igrid

        integer :: ix^L,iib^D

        block=>psb(igrid)
        ^d&iib^d=idphyb(^d,igrid);
        {do i^db=-1,1\}
           if (skip_direction([ i^d ])) cycle
           ix^l=ixr_srl_^l(iib^d,i^d);
           call phys_get_aux(.true.,psb(igrid)%w,psb(igrid)%x,ixg^ll,ix^l,"bc")
        {end do\}

      end subroutine fix_auxiliary

  end subroutine getbc

  subroutine identifyphysbound(s,iib^D)
    use mod_global_parameters

    type(state)          :: s
    integer, intent(out) :: iib^D

    {
    if(s%is_physical_boundary(2*^d) .and. &
       s%is_physical_boundary(2*^d-1)) then
      iib^d=2
    else if(s%is_physical_boundary(2*^d-1)) then
      iib^d=-1
    else if(s%is_physical_boundary(2*^d)) then
      iib^d=1
    else
      iib^d=0
    end if
    \}

  end subroutine identifyphysbound

end module mod_ghostcells_update