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

  ! total size of buffer arrays
  integer :: nbuff_bc_recv_srl, nbuff_bc_send_srl, nbuff_bc_recv_r, nbuff_bc_send_r, nbuff_bc_recv_p, nbuff_bc_send_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

  ! 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 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;
    
    select case (typeghostfill)
    case ("copy")
       interpolation_order=1
    case ("linear")
       interpolation_order=2
    case default
       write (unitterm,*) "Undefined typeghostfill ",typeghostfill
       call mpistop("Undefined typeghostfill")
    end select
    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 :: my_neighbor_type, ipole, idims, iside, nwhead, nwtail
    integer :: iigrid, igrid, ineighbor, ipe_neighbor
    integer :: nrecvs, nsends, isizes
    integer :: ixG^L, ixR^L, ixS^L, ixB^L, ixI^L, k^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),bindex(ndim)
    integer :: isend_buf(npwbuf), ipwbuf, nghostcellsco,iB
    ! 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, NeedProlong(-1:1^d&)
    type(wbuffer) :: pwbuf(npwbuf)

    ! Stretching grid parameters for coarsened block of the current block

    nwhead=nwstart
    nwtail=nwstart+nwbc-1

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

    time_bcin=mpi_wtime()
    ixg^l=ixg^ll;
    
    if (internalboundary) then 
       call getintbc(time,ixg^l)
    end if
    ! fill ghost cells in physical boundaries
    if(bcphys.and. .not.stagger_grid) then
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if(.not.phyboundblock(igrid)) cycle
        saveigrid=igrid
        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=ixgmin^d+kmin^d*nghostcells
           ixbmax^d=ixgmax^d-kmax^d*nghostcells
          \}
          {^iftwod
          if(idims > 1 .and. neighbor_type(-1,0,igrid)==neighbor_boundary) ixbmin1=ixgmin1
          if(idims > 1 .and. neighbor_type( 1,0,igrid)==neighbor_boundary) ixbmax1=ixgmax1}
          {^ifthreed
          if(idims > 1 .and. neighbor_type(-1,0,0,igrid)==neighbor_boundary) ixbmin1=ixgmin1
          if(idims > 1 .and. neighbor_type( 1,0,0,igrid)==neighbor_boundary) ixbmax1=ixgmax1
          if(idims > 2 .and. neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixbmin2=ixgmin2
          if(idims > 2 .and. neighbor_type(0, 1,0,igrid)==neighbor_boundary) ixbmax2=ixgmax2}
          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^l,ixb^l)
          end do
        end do
      end do
    end if

    ! default : no singular axis
    ipole=0
    
    irecv=0
    isend=0
    isend_buf=0
    ipwbuf=1
    ! total number of times to call MPI_IRECV in each processor between sibling blocks or from finer neighbors
    nrecvs=nrecv_bc_srl+nrecv_bc_r
    ! total number of times to call MPI_ISEND in each processor between sibling blocks or to coarser neighors
    nsends=nsend_bc_srl+nsend_bc_r

    allocate(recvstatus(mpi_status_size,nrecvs),recvrequest(nrecvs))
    recvrequest=mpi_request_null

    allocate(sendstatus(mpi_status_size,nsends),sendrequest(nsends))
    sendrequest=mpi_request_null

    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
    ! receiving ghost-cell values from sibling blocks and finer neighbors
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       saveigrid=igrid
       call identifyphysbound(ps(igrid),iib^d)   
       ^d&idphyb(^d,igrid)=iib^d;
       {do i^db=-1,1\}
          if (skip_direction([ i^d ])) cycle
          my_neighbor_type=neighbor_type(i^d,igrid)
          select case (my_neighbor_type)
          case (neighbor_sibling)
             call bc_recv_srl
          case (neighbor_fine)
             call bc_recv_restrict
          end select
       {end do\}
    end do
    
    ! sending ghost-cell values to sibling blocks and coarser neighbors
    nghostcellsco=ceiling(nghostcells*0.5d0)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       saveigrid=igrid
       block=>psb(igrid)

       ! Used stored data to identify physical boundaries
       ^d&iib^d=idphyb(^d,igrid);

       if (any(neighbor_type(:^d&,igrid)==neighbor_coarse)) then
          ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
    {#IFDEF EVOLVINGBOUNDARY
          if(phyboundblock(igrid)) then
            ! coarsen finer ghost cells at physical boundaries
            ixcommin^d=ixcogmin^d+nghostcellsco;
            ixcommax^d=ixcogmax^d-nghostcellsco;
            ixmmin^d=ixgmin^d+(nghostcellsco-1);
            ixmmax^d=ixgmax^d-(nghostcellsco-1);
          else
            ixcom^l=ixcog^l^lsubnghostcells;
            ixm^l=ixg^l^lsubnghostcells;
          end if
    }
          call coarsen_grid(psb(igrid),ixg^l,ixm^l,psc(igrid),ixcog^l,ixcom^l)
       end if
    
       {do i^db=-1,1\}
          if (skip_direction([ i^d ])) cycle
          if (phi_ > 0) ipole=neighbor_pole(i^d,igrid)
          my_neighbor_type=neighbor_type(i^d,igrid)
          select case (my_neighbor_type)
          case (neighbor_sibling)
             call bc_send_srl
          case (neighbor_coarse)
             call bc_send_restrict
          end select
       {end do\}
    end do
    
    !if (irecv/=nrecvs) then
    !   call mpistop("number of recvs in phase1 in amr_ghostcells is incorrect")
    !end if
    !if (isend/=nsends) then
    !   call mpistop("number of sends in phase1 in amr_ghostcells is incorrect")
    !end if
    
    call mpi_waitall(irecv,recvrequest,recvstatus,ierrmpi)
    call mpi_waitall(isend,sendrequest,sendstatus,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 to fill ghost cells
      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
          case (neighbor_fine)
             call bc_fill_r
          end select
       {end do\}
      end do
    end if

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

    irecv=0
    isend=0
    isend_buf=0
    ipwbuf=1
    ! total number of times to call MPI_IRECV in each processor from coarser neighbors
    nrecvs=nrecv_bc_p
    ! total number of times to call MPI_ISEND in each processor to finer neighbors
    nsends=nsend_bc_p

    allocate(recvstatus(mpi_status_size,nrecvs),recvrequest(nrecvs))
    recvrequest=mpi_request_null
    allocate(sendstatus(mpi_status_size,nsends),sendrequest(nsends))
    sendrequest=mpi_request_null

    ! receiving ghost-cell values from coarser neighbors
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       saveigrid=igrid
       ^d&iib^d=idphyb(^d,igrid);
       {do i^db=-1,1\}
          if (skip_direction([ i^d ])) cycle
          my_neighbor_type=neighbor_type(i^d,igrid)
          if (my_neighbor_type==neighbor_coarse) call bc_recv_prolong
       {end do\}
    end do
    ! sending ghost-cell values to finer neighbors 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       saveigrid=igrid
       block=>psb(igrid)
       ^d&iib^d=idphyb(^d,igrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       if (any(neighbor_type(:^d&,igrid)==neighbor_fine)) then
          {do i^db=-1,1\}
             if (skip_direction([ i^d ])) cycle
             if (phi_ > 0) ipole=neighbor_pole(i^d,igrid)
             my_neighbor_type=neighbor_type(i^d,igrid)
             if (my_neighbor_type==neighbor_fine) call bc_send_prolong
          {end do\}
       end if
    end do
    
    !if (irecv/=nrecvs) then
    !   call mpistop("number of recvs in phase2 in amr_ghostcells is incorrect")
    !end if
    !if (isend/=nsends) then
    !   call mpistop("number of sends in phase2 in amr_ghostcells is incorrect")
    !end if
    
    call mpi_waitall(irecv,recvrequest,recvstatus,ierrmpi)
    call mpi_waitall(isend,sendrequest,sendstatus,ierrmpi)
    deallocate(recvstatus,recvrequest)
    deallocate(sendstatus,sendrequest)

    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 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_p
         {end do\}
      end do
    end if
    ! do prolongation on the ghost-cell values received from coarser neighbors 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&iib^d=idphyb(^d,igrid);
       if (any(neighbor_type(:^d&,igrid)==neighbor_coarse)) then
         needprolong=.false.
         {do i^db=-1,1\}
            if (skip_direction([ i^d ])) cycle
            my_neighbor_type=neighbor_type(i^d,igrid)
            if (my_neighbor_type==neighbor_coarse) then
              call bc_prolong
              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(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(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(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(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(needprolong)
          {end do\}
         end if
       end if
    end do
    
    do ipwbuf=1,npwbuf
       if (isend_buf(ipwbuf)/=0) deallocate(pwbuf(ipwbuf)%w)
    end do

    if(stagger_grid) then
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if(.not.phyboundblock(igrid)) cycle
        saveigrid=igrid
        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^l,ixb^l)
          end do
        end do
      end do
    end if

     ! modify normal component of magnetic field to fix divB=0 
    if(bcphys.and.associated(phys_boundary_adjust)) call phys_boundary_adjust()
    
    if (nwaux>0) call fix_auxiliary
    
    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

      !> Receive from sibling at same refinement level
      subroutine bc_recv_srl

        ipe_neighbor=neighbor(2,i^d,igrid)
        if (ipe_neighbor/=mype) then
           irecv=irecv+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(irecv),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=irecv+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(irecv),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

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

        if (ipole==0) then
           n_i^d=-i^d;
           if (ipe_neighbor==mype) then
              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)
           else
              isend=isend+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(isend),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
           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 (ipe_neighbor==mype) then
              ixr^l=ixr_srl_^l(iib^d,n_i^d);
              call pole_copy(psb(ineighbor)%w,ixg^l,ixr^l,psb(igrid)%w,ixg^l,ixs^l)
           else
              if (isend_buf(ipwbuf)/=0) then
                 call mpi_wait(sendrequest(isend_buf(ipwbuf)), &
                               sendstatus(:,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^l,ixs^l)
              isend=isend+1
              isend_buf(ipwbuf)=isend
              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(isend),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

      !> Send to coarser neighbor
      subroutine bc_send_restrict
        integer :: ii^D

        ic^d=1+modulo(node(pig^d_,igrid)-1,2);
        if ({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return
        if(phyboundblock(igrid).and..not.stagger_grid) 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

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

        if (ipole==0) then
           n_inc^d=-2*i^d+ic^d;
           if (ipe_neighbor==mype) then
              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)
           else
              isend=isend+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(isend),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
           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 (ipe_neighbor==mype) then
              ixr^l=ixr_r_^l(iib^d,n_inc^d);
              call pole_copy(psb(ineighbor)%w,ixg^l,ixr^l,psc(igrid)%w,ixcog^l,ixs^l)
           else
              if (isend_buf(ipwbuf)/=0) then
                 call mpi_wait(sendrequest(isend_buf(ipwbuf)), &
                               sendstatus(:,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=isend+1
              isend_buf(ipwbuf)=isend
              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(isend),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 siblings ghost cells with received data
      subroutine bc_fill_srl
        double precision :: tmp(ixgs^t)
        integer :: ixS^L,ixR^L,n_i^D,ixSsync^L,ixRsync^L
        integer :: idir, idirect

        ineighbor=neighbor(1,i^d,igrid)
        ipe_neighbor=neighbor(2,i^d,igrid)
        ipole=neighbor_pole(i^d,igrid)
        idirect={abs(i^d)|+}

        !! 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);
          if (ipe_neighbor==mype) then
            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);
              !if (idirect == 1 .and. qdt>0.d0) then
              !  ! use the same value at the face shared by two neighors
              !  call indices_for_syncing(idir,i^D,ixR^L,ixS^L,ixRsync^L,ixSsync^L)
              !   if({idir==^D .and. i^D /=0  | .or.}) write(*,*) it,'difference',maxval(abs(psb(igrid)%ws(ixRsync^S,idir)-psb(ineighbor)%ws(ixSsync^S,idir)))
              !  !psb(igrid)%ws(ixRsync^S,idir) = half*(psb(igrid)%ws(ixRsync^S,idir)+psb(ineighbor)%ws(ixSsync^S,idir))
              !end if
              psb(igrid)%ws(ixr^s,idir) = psb(ineighbor)%ws(ixs^s,idir)
            end do
          else
            !! 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)))       
              !if (idirect == 1) then
              !   ! ixR ixS maybe changed
              !   call indices_for_syncing(idir,i^D,ixR^L,ixS^L,ixRsync^L,ixSsync^L) ! Overwrites ixR, ixS
              !   if(qdt==0) psb(igrid)%ws(ixRsync^S,idir) = half*(tmp(ixSsync^S) + psb(igrid)%ws(ixRsync^S,idir))
              !   if({idir==^D .and. i^D /=0  | .or.}) write(*,*) it,'betweenpe',maxval(abs(psb(igrid)%ws(ixRsync^S,idir)-tmp(ixSsync^S)))
              !end if
              psb(igrid)%ws(ixr^s,idir) = tmp(ixs^s)
              ibuf_recv_srl=ibuf_next
            end do
          end if

        else ! There is a pole
          select case (ipole)
          {case (^d)
             n_i^d=i^d^d%n_i^dd=-i^dd;\}
          end select
          if (ipe_neighbor==mype) then
            !! Fill ghost cells
            do idir=1,ndim
              ixr^l=ixr_srl_stg_^l(idir,i^d);
              ixs^l=ixs_srl_stg_^l(idir,n_i^d);
              !! Fill ghost cells
              call pole_copy_stg(psb(igrid)%ws,ixr^l,psb(ineighbor)%ws,ixs^l,idir)
            end do
          else
             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)
             end do
          end if
        end if

      end subroutine bc_fill_srl

      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_r

        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\}
             n_i^d=-i^d;
             ineighbor=neighbor_child(1,inc^d,igrid)
             ipe_neighbor=neighbor_child(2,inc^d,igrid)
             if (ipe_neighbor==mype) then ! Same processor
               do idir=1,ndim
                  ixs^l=ixs_r_stg_^l(idir,n_i^d);
                  ixr^l=ixr_r_stg_^l(idir,inc^d);
                  psb(igrid)%ws(ixr^s,idir)=psc(ineighbor)%ws(ixs^s,idir)
               end do
             else ! Different processor
               !! 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\}
             select case(ipole)
            {case (^d)
               n_i^d=i^d^d%n_i^dd=-i^dd;\}
             end select
             ineighbor=neighbor_child(1,inc^d,igrid)
             ipe_neighbor=neighbor_child(2,inc^d,igrid)
             ixr^l=ixr_r_^l(iib^d,inc^d);
             if (ipe_neighbor==mype) then ! Same processor
               do idir=1,ndim
                 ixs^l=ixs_r_stg_^l(idir,n_i^d);
                 ixr^l=ixr_r_stg_^l(idir,inc^d);
                 !! Fill ghost cells
                 call pole_copy_stg(psb(igrid)%ws,ixr^l,psc(ineighbor)%ws,ixs^l,idir)
               end do
             else ! Different processor
               !! 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)
                 ibuf_recv_r=ibuf_next
               end do
             end if
          {end do\}
        
        end if

      end subroutine bc_fill_r

      !> 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=irecv+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(irecv),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

        {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
           inc^db=2*i^db+ic^db\}
           ixs^l=ixs_p_^l(iib^d,inc^d);

           ineighbor=neighbor_child(1,inc^d,igrid)
           ipe_neighbor=neighbor_child(2,inc^d,igrid)
        
           if (ipole==0) then
              n_i^d=-i^d;
              n_inc^d=ic^d+n_i^d;
              if (ipe_neighbor==mype) then
                ixr^l=ixr_p_^l(iib^d,n_inc^d);
                psc(ineighbor)%w(ixr^s,nwhead:nwtail) &
                   =psb(igrid)%w(ixs^s,nwhead:nwtail)
              else
                isend=isend+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(isend),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
              end if
           else
              select case (ipole)
              {case (^d)
                 n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
              end select
              if (ipe_neighbor==mype) then
                 ixr^l=ixr_p_^l(iib^d,n_inc^d);
                 call pole_copy(psc(ineighbor)%w,ixcog^l,ixr^l,psb(igrid)%w,ixg^l,ixs^l)
                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)
                  end do
                end if
              else
                if (isend_buf(ipwbuf)/=0) then
                   call mpi_wait(sendrequest(isend_buf(ipwbuf)), &
                                 sendstatus(:,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^l,ixs^l)
                isend=isend+1
                isend_buf(ipwbuf)=isend
                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(isend),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

      !> fill coarser representative with data from coarser neighbors
      subroutine bc_fill_p
        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)
        ipe_neighbor=neighbor(2,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);
          if(ipe_neighbor==mype) then !! Same processor
            n_inc^d=-2*i^d+ic^d;
            do idir=1,ndim
              ixs^l=ixs_p_stg_^l(idir,n_inc^d);
              ixr^l=ixr_p_stg_^l(idir,inc^d);
              psc(igrid)%ws(ixr^s,idir)=psb(ineighbor)%ws(ixs^s,idir)
            end do
          else !! Different processor
            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
          end if

        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
          if (ipe_neighbor==mype) then
            do idir=1,ndim
              ixs^l=ixs_p_stg_^l(idir,n_inc^d);
              ixr^l=ixr_p_stg_^l(idir,inc^d);
              call pole_copy_stg(psc(igrid)%ws,ixr^l,psb(ineighbor)%ws,ixs^l,idir)
            end do
          else
            !! 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)
              ibuf_recv_p=ibuf_next
            end do
          end if

        end if

      end subroutine bc_fill_p

      !> do prolongation for fine blocks after receipt data from coarse neighbors
      subroutine bc_prolong
        use mod_physics, only: phys_to_primitive, phys_to_conserved

        integer :: ixFi^L,ixCo^L,ii^D
        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)) then
          block=>psc(igrid)
          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;
          do idims=1,ndim
             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

        select case (typeghostfill)
        case ("linear")
           call interpolation_linear(ixfi^l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
        case ("copy")
           call interpolation_copy(ixfi^l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
        case default
           write (unitterm,*) "Undefined typeghostfill ",typeghostfill
           call mpistop("Undefined typeghostfill")
        end select

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

      end subroutine bc_prolong

      subroutine bc_prolong_stg(NeedProlong)
        use mod_amr_fct
        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;

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

        ! 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(ixFi^L,dxFi^D,xFimin^D, &
                                      dxCo^D,invdxCo^D,xComin^D)
        use mod_physics, only: phys_to_conserved
        integer, intent(in) :: 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) 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;
             {signedfactorhalf^d=(xfi^d-xco^d)*invdxco^d*two
              if(dabs(signedfactorhalf^d**2-1.0d0/4.0d0)>smalldouble) call mpistop("error in bc_prolong")
              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(typeprolonglimit)
                 case('unlimit')
                   slope(iw,idims)=slopec
                 case('minmod')
                   slope(iw,idims)=signr*max(zero,min(dabs(sloper), &
                                                     signr*slopel))
                 case('woodward')
                   slope(iw,idims)=two*signr*max(zero,min(dabs(sloper), &
                                      signr*slopel,signr*half*slopec))
                 case('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(ixFi^L,dxFi^D,xFimin^D, &
                                    dxCo^D,invdxCo^D,xComin^D)
        use mod_physics, only: phys_to_conserved
        integer, intent(in) :: 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)
      
        integer, intent(in) :: ixIR^L,ixR^L,ixIS^L,ixS^L
        double precision :: wrecv(ixir^s,1:nw), wsend(ixis^s,1:nw)

        integer :: iw, 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)
      
        integer, intent(in) :: ixR^L,ixS^L,idirs
        double precision :: wrecv(ixgs^t,1:nws), wsend(ixgs^t,1:nws)

        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, 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
        use mod_physics, only: phys_get_aux
      
        integer :: ix^L

        do iigrid=1,igridstail; igrid=igrids(iigrid);
          saveigrid=igrid
          block=>psb(igrid)
          call identifyphysbound(psb(igrid),iib^d)   
             
          {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,ps(igrid)%x,ixg^l,ix^l,"bc")
          {end do\}
        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