getbc Subroutine

public subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)

Uses

  • proc~~getbc~~UsesGraph proc~getbc getbc mod_physics mod_physics proc~getbc->mod_physics module~mod_boundary_conditions mod_boundary_conditions proc~getbc->module~mod_boundary_conditions module~mod_coarsen mod_coarsen proc~getbc->module~mod_coarsen module~mod_comm_lib mod_comm_lib proc~getbc->module~mod_comm_lib module~mod_global_parameters mod_global_parameters proc~getbc->module~mod_global_parameters module~mod_nvtx mod_nvtx proc~getbc->module~mod_nvtx module~mod_coarsen->mod_physics module~mod_coarsen->module~mod_global_parameters module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi iso_c_binding iso_c_binding module~mod_nvtx->iso_c_binding module~mod_variables->module~mod_basic_types

do update ghost cells of all blocks including physical boundaries

Arguments

Type IntentOptional Attributes Name
double precision, intent(in) :: time
double precision, intent(in) :: qdt
type(state), target :: psb(max_blocks)
integer, intent(in) :: nwstart
integer, intent(in) :: nwbc
logical, intent(in), optional :: req_diag

Calls

proc~~getbc~~CallsGraph proc~getbc getbc mpi_waitall mpi_waitall proc~getbc->mpi_waitall mpi_wtime mpi_wtime proc~getbc->mpi_wtime phys_boundary_adjust phys_boundary_adjust proc~getbc->phys_boundary_adjust proc~fill_boundary_after_gc fill_boundary_after_gc proc~getbc->proc~fill_boundary_after_gc proc~fill_boundary_before_gc fill_boundary_before_gc proc~getbc->proc~fill_boundary_before_gc proc~fill_coarse_boundary fill_coarse_boundary proc~getbc->proc~fill_coarse_boundary proc~getintbc getintbc proc~getbc->proc~getintbc proc~idecode~2 idecode proc~getbc->proc~idecode~2 proc~mpi_irecv_wrapper mpi_irecv_wrapper proc~getbc->proc~mpi_irecv_wrapper proc~mpi_isend_wrapper mpi_isend_wrapper proc~getbc->proc~mpi_isend_wrapper proc~nvtxendrange~2 nvtxEndRange proc~getbc->proc~nvtxendrange~2 proc~nvtxstartrange~2 nvtxStartRange proc~getbc->proc~nvtxstartrange~2 proc~skip_direction skip_direction proc~getbc->proc~skip_direction proc~bc_phys bc_phys proc~fill_boundary_after_gc->proc~bc_phys proc~fill_boundary_before_gc->proc~bc_phys proc~fill_coarse_boundary->proc~bc_phys proc~sync_positions_host sync_positions_host proc~getintbc->proc~sync_positions_host mpi_irecv mpi_irecv proc~mpi_irecv_wrapper->mpi_irecv mpi_isend mpi_isend proc~mpi_isend_wrapper->mpi_isend specialbound_usr specialbound_usr proc~bc_phys->specialbound_usr

Called by

proc~~getbc~~CalledByGraph proc~getbc getbc proc~add_split_source add_split_source proc~add_split_source->proc~getbc proc~add_split_source~2 add_split_source proc~add_split_source~2->proc~getbc proc~add_split_source~3 add_split_source proc~add_split_source~3->proc~getbc proc~add_split_source~4 add_split_source proc~add_split_source~4->proc~getbc proc~advect1 advect1 proc~advect1->proc~getbc proc~amr_coarsen_refine amr_coarsen_refine proc~amr_coarsen_refine->proc~getbc proc~global_implicit_update global_implicit_update proc~global_implicit_update->proc~getbc proc~improve_initial_condition improve_initial_condition proc~improve_initial_condition->proc~getbc proc~recalculateb recalculateB proc~improve_initial_condition->proc~recalculateb proc~initlevelone initlevelone proc~initlevelone->proc~getbc proc~main main proc~main->proc~getbc proc~main->proc~improve_initial_condition proc~main->proc~initlevelone proc~process process proc~main->proc~process proc~resettree_convert resettree_convert proc~main->proc~resettree_convert proc~settree settree proc~main->proc~settree proc~timeintegration timeintegration proc~main->proc~timeintegration proc~process->proc~getbc proc~process_advanced process_advanced proc~process_advanced->proc~getbc proc~recalculateb->proc~getbc proc~resettree_convert->proc~getbc proc~resettree_convert->proc~amr_coarsen_refine proc~saveamrfile saveamrfile proc~saveamrfile->proc~getbc proc~advect advect proc~advect->proc~advect1 proc~resettree resettree proc~resettree->proc~amr_coarsen_refine proc~settree->proc~amr_coarsen_refine proc~timeintegration->proc~process proc~timeintegration->proc~process_advanced proc~timeintegration->proc~saveamrfile proc~timeintegration->proc~resettree proc~advance advance proc~timeintegration->proc~advance program~agile agile program~agile->proc~main proc~advance->proc~advect

Source Code

  subroutine getbc(time,qdt,psb,nwstart,nwbc,req_diag)
    use mod_nvtx
    use mod_global_parameters
    use mod_physics
    use mod_coarsen, only: coarsen_grid
    use mod_boundary_conditions, only: getintbc, bc_phys
    use mod_comm_lib, only: mpistop

    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 :: nwhead, nwtail, bgstep
    integer :: iigrid, igrid, ineighbor, ipe_neighbor
    integer :: ixRmin1,ixRmin2,ixRmin3,ixRmax1,ixRmax2,ixRmax3, ixSmin1,&
         ixSmin2,ixSmin3,ixSmax1,ixSmax2,ixSmax3
    integer :: i1,i2,i3, n_i1,n_i2,n_i3, ic1,ic2,ic3, inc1,inc2,inc3, n_inc1,&
         n_inc2, n_inc3, iib1, iib2, iib3
    real(dp) :: tempval

    integer :: ix1,ix2,ix3, iw, inb, i, Nx1, Nx2, Nx3, ienc, ibuf_start
    double precision :: CoFiratio
    ! for prolongation:
    integer  :: ixCo1, ixCo2, ixCo3, ixFi1, ixFi2, ixFi3, idims
    integer  :: ixFimin1, ixFimin2, ixFimin3, ixFimax1, ixFimax2, ixFimax3
    integer  :: hxCo1, hxCo2, hxCo3, jxCo1, jxCo2, jxCo3, nigrids
    real(dp) :: slope(ndim), dxFi1, dxFi2, dxFi3, dxCo1, dxCo2, dxCo3
    real(dp) :: invdxCo1, invdxCo2, invdxCo3
    real(dp) :: xFimin1, xFimin2, xFimin3, xComin1, xComin2, xComin3
    real(dp) :: xFi1, xFi2, xFi3, xCo3, xCo2, xCo1
    real(dp) :: eta1, eta2, eta3, slopeL, slopeR, slopeC, signR, signC
    ! Polar axis.  ipole is the dimension the pole lies across (0 for an
    ! ordinary neighbour, so every expression below collapses to the plain
    ! copy), iB_pole is the physical-boundary index whose typeboundary column
    ! holds the per-variable sign, and sbase/sstep walk the source backwards
    ! along the pole dimension so the row nearest the axis lands nearest the
    ! axis.  See "The polar axis" in CLAUDE.md for the derivation.
    integer  :: ipole, iB_pole
    integer  :: sbase1, sbase2, sbase3, sstep1, sstep2, sstep3

    time_bcin=MPI_WTIME()
    call nvtxStartRange("getbc",2)

    nwhead = nwstart
    nwtail = nwstart+nwbc-1
    ! igrids(1) is undefined when this rank owns no blocks
    if (igridstail > 0) then
       bgstep = psb(igrids(1))%istep
    else
       bgstep = 1
    end if

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

    ! fill internal physical boundary
    if (internalboundary) then
       call getintbc(time,ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3)
    end if

    ! fill physical-boundary ghost cells before internal ghost-cell values exchange
    if(bcphys.and. .not.stagger_grid) then
        
       do iigrid = 1, igridstail; igrid=igrids(iigrid);
          if (.not.phyboundblock(igrid)) cycle
          call fill_boundary_before_gc(psb(igrid),igrid,time,qdt)
       end do
    end if


    ! prepare coarse values to send to coarser neighbors
     
    do iigrid = 1, igridstail; igrid=igrids(iigrid);
       if (any(neighbor_type(:,:,:,igrid)==neighbor_coarse)) then

          CoFiratio=one/dble(2**ndim)
          
          do iw = nwhead, nwtail
             do ixCo3 = ixCoMmin3,ixCoMmax3
                do ixCo2 = ixCoMmin2,ixCoMmax2
                   do ixCo1 = ixCoMmin1,ixCoMmax1
                      ixFi3=2*(ixCo3-ixCoMmin3)+ixMmin3
                      ixFi2=2*(ixCo2-ixCoMmin2)+ixMmin2
                      ixFi1=2*(ixCo1-ixCoMmin1)+ixMmin1

                      bgc(1)%w(ixCo1,ixCo2,ixCo3,iw,igrid) = 0.0d0
                      do ix3 = ixFi3,ixFi3+1
                         do ix2 = ixFi2,ixFi2+1
                            do ix1 = ixFi1,ixFi1+1
                               bgc(1)%w(ixCo1,ixCo2,ixCo3,iw,igrid) = bgc(1)%w(ixCo1,ixCo2,ixCo3,iw,igrid) &
                                    + bg(bgstep)%w(ix1,ix2,ix3,iw,igrid)
                            end do
                         end do
                      end do
                      bgc(1)%w(ixCo1,ixCo2,ixCo3,iw,igrid) = bgc(1)%w(ixCo1,ixCo2,ixCo3,iw,igrid)*coFiRatio

                   end do
                end do
             end do
          end do

          do i3 = -1, 1
             do i2 = -1, 1
                do i1 = -1, 1
                   if (skip_direction([ i1,i2,i3 ])) cycle
                   if (neighbor_type(i1,i2,i3,igrid)==neighbor_coarse) call &
                        fill_coarse_boundary(time,igrid,i1,i2,i3)
                end do
             end do
          end do
       end if
    end do

    ! MPI receive SRL
    do inb = 1, nbprocs_info%nbprocs_srl
#ifndef NOGPUDIRECT
#ifdef _CRAYFTN
      
#else
      
#endif
#endif
       call mpi_irecv_wrapper(nbprocs_info%srl_nb(inb)%rcv%buffer, &
            nbprocs_info%srl_nb(inb)%rcv%size, &
            MPI_DOUBLE_PRECISION, nbprocs_info%nbprocs_srl_list(inb), 1, icomm, recv_srl_nb(inb), ierrmpi)
       call mpi_irecv_wrapper(nbprocs_info%srl_nb(inb)%info_rcv%buffer, &
            nbprocs_info%srl_nb(inb)%info_rcv%size, &
            MPI_INTEGER, nbprocs_info%nbprocs_srl_list(inb), 2, icomm, recv_srl_nb(nbprocs_info%nbprocs_srl + inb), ierrmpi)

#ifndef NOGPUDIRECT
      
#endif
    end do

    ! MPI receive restrict (neighbor is finer)
    do inb = 1, nbprocs_info%nbprocs_f
#ifndef NOGPUDIRECT
#ifdef _CRAYFTN
      
#else
      
#endif
#endif
       call mpi_irecv_wrapper(nbprocs_info%fine_nb(inb)%rcv%buffer, &
            nbprocs_info%fine_nb(inb)%rcv%size, &
            MPI_DOUBLE_PRECISION, nbprocs_info%nbprocs_f_list(inb), 1, icomm, recv_f_nb(inb), ierrmpi)
       call mpi_irecv_wrapper(nbprocs_info%fine_nb(inb)%info_rcv%buffer, &
            nbprocs_info%fine_nb(inb)%info_rcv%size, &
            MPI_INTEGER, nbprocs_info%nbprocs_f_list(inb), 2, icomm, recv_f_nb(nbprocs_info%nbprocs_f + inb), ierrmpi)
#ifndef NOGPUDIRECT
      
#endif
    end do

    ! fill the SRL send buffers on GPU
    do inb = 1, nbprocs_info%nbprocs_srl
        
       do i = 1, nbprocs_info%srl_nb(inb)%info%nigrids
             igrid = nbprocs_info%srl_nb(inb)%info%igrid(i)
             ienc = nbprocs_info%srl_nb(inb)%info%iencode(i)
             ibuf_start = nbprocs_info%srl_nb(inb)%info%ibuf_start(i)
             call idecode( i1, i2, i3, ienc )
             iib1=idphyb(1,igrid); iib2=idphyb(2,igrid); iib3=idphyb(3,igrid)

             ! now fill the data and info buffers
             ixSmin1=ixS_srl_min1(iib1,i1); ixSmax1=ixS_srl_max1(iib1,i1)
             ixSmin2=ixS_srl_min2(iib2,i2); ixSmax2=ixS_srl_max2(iib2,i2)
             ixSmin3=ixS_srl_min3(iib3,i3); ixSmax3=ixS_srl_max3(iib3,i3)
             Nx1=ixSmax1-ixSmin1+1; Nx2=ixSmax2-ixSmin2+1; Nx3=ixSmax3-ixSmin3+1

             ! The mirror and the sign are applied here, on the send side, so
             ! that the buffer already holds what the receiver has to store and
             ! the unpack stays a plain strided copy - which is what lets the
             ! receiver get away with knowing nothing about the pole beyond the
             ! destination offset shipped in the info record below.
             ipole=neighbor_pole(i1,i2,i3,igrid)
             n_i1=-i1; n_i2=-i2; n_i3=-i3
             sbase1=ixSmin1; sstep1=1; sbase2=ixSmin2; sstep2=1
             sbase3=ixSmin3; sstep3=1; iB_pole=1
             if (ipole==1) then
                n_i1=i1; sbase1=ixSmax1; sstep1=-1; iB_pole=(i1+3)/2
             else if (ipole==2) then
                n_i2=i2; sbase2=ixSmax2; sstep2=-1; iB_pole=2+(i2+3)/2
             else if (ipole==3) then
                n_i3=i3; sbase3=ixSmax3; sstep3=-1; iB_pole=4+(i3+3)/2
             end if

              
             do iw = nwhead, nwtail
                do ix3 = ixSmin3, ixSmax3
                   do ix2 = ixSmin2, ixSmax2
                      do ix1 = ixSmin1, ixSmax1
                         nbprocs_info%srl_nb(inb)%send%buffer( &
                              ibuf_start &
                              + sstep1*(ix1-sbase1) &
                              + Nx1 * (sstep2*(ix2-sbase2)) &
                              + Nx1*Nx2 * (sstep3*(ix3-sbase3)) &
                              + Nx1*Nx2*Nx3 * (iw-nwhead) &
                              ) = merge(-one,one,ipole/=0 .and. &
                              typeboundary(iw,iB_pole)==bc_asymm) &
                              * bg(bgstep)%w( ix1, ix2, ix3, iw, igrid)
                      end do
                   end do
                end do
             end do
             ! ship the destination offset rather than our own, so that the
             ! receiver does not have to undo the mirror it cannot see
             nbprocs_info%srl_nb(inb)%info_send%buffer( 1 + 3 * (i - 1) : 3 * i ) = &
                  [neighbor(1,i1,i2,i3,igrid), &
                   1 + (n_i1+1) + 3*(n_i2+1) + 9*(n_i3+1), ibuf_start]
       end do
    end do

    ! fill the C send buffers on GPU (send_restrict)
    do inb = 1, nbprocs_info%nbprocs_c
        
       do i = 1, nbprocs_info%course_nb(inb)%info%nigrids

          igrid = nbprocs_info%course_nb(inb)%info%igrid(i)
          i1 = nbprocs_info%course_nb(inb)%info%i1(i)
          i2 = nbprocs_info%course_nb(inb)%info%i2(i)
          i3 = nbprocs_info%course_nb(inb)%info%i3(i)
          ipole=neighbor_pole(i1,i2,i3,igrid)

          ! how they will be used at the receiving end.  The plain swap is
          ! n_inc = -2*i+ic; across a pole it becomes n_inc = 2*i+(3-ic),
          ! which on this 0..3 index leaves 0 and 3 alone and swaps 1 with 2.
          inc1 = nbprocs_info%course_nb(inb)%info%inc1(i)
          if (ipole == 1) then
             if (inc1 == 1) then; inc1 = 2; else if (inc1 == 2) then; inc1 = 1; end if
          else
             if (inc1 == 0) then; inc1 = 3; else if (inc1 == 3) then; inc1 = 0; end if
          end if

          inc2 = nbprocs_info%course_nb(inb)%info%inc2(i)
          if (ipole == 2) then
             if (inc2 == 1) then; inc2 = 2; else if (inc2 == 2) then; inc2 = 1; end if
          else
             if (inc2 == 0) then; inc2 = 3; else if (inc2 == 3) then; inc2 = 0; end if
          end if

          inc3 = nbprocs_info%course_nb(inb)%info%inc3(i)
          if (ipole == 3) then
             if (inc3 == 1) then; inc3 = 2; else if (inc3 == 2) then; inc3 = 1; end if
          else
             if (inc3 == 0) then; inc3 = 3; else if (inc3 == 3) then; inc3 = 0; end if
          end if

          ibuf_start = nbprocs_info%course_nb(inb)%info%ibuf_start(i)
          iib1=idphyb(1,igrid); iib2=idphyb(2,igrid); iib3=idphyb(3,igrid)

          ! now fill the data and info buffers
          ixSmin1=ixS_r_min1(iib1,i1);ixSmin2=ixS_r_min2(iib2,i2)
          ixSmin3=ixS_r_min3(iib3,i3);ixSmax1=ixS_r_max1(iib1,i1)
          ixSmax2=ixS_r_max2(iib2,i2);ixSmax3=ixS_r_max3(iib3,i3)
          Nx1=ixSmax1-ixSmin1+1; Nx2=ixSmax2-ixSmin2+1; Nx3=ixSmax3-ixSmin3+1

          ! mirror and sign on the send side, as for srl above
          sbase1=ixSmin1; sstep1=1; sbase2=ixSmin2; sstep2=1
          sbase3=ixSmin3; sstep3=1; iB_pole=1
          if (ipole==1) then
             sbase1=ixSmax1; sstep1=-1; iB_pole=(i1+3)/2
          else if (ipole==2) then
             sbase2=ixSmax2; sstep2=-1; iB_pole=2+(i2+3)/2
          else if (ipole==3) then
             sbase3=ixSmax3; sstep3=-1; iB_pole=4+(i3+3)/2
          end if

           
          do iw = nwhead, nwtail
             do ix3 = ixSmin3, ixSmax3
                do ix2 = ixSmin2, ixSmax2
                   do ix1 = ixSmin1, ixSmax1
                      nbprocs_info%course_nb(inb)%send%buffer( &
                           ibuf_start &
                           + sstep1*(ix1-sbase1) &
                           + Nx1 * (sstep2*(ix2-sbase2)) &
                           + Nx1*Nx2 * (sstep3*(ix3-sbase3)) &
                           + Nx1*Nx2*Nx3 * (iw-nwhead) &
                           ) = merge(-one,one,ipole/=0 .and. &
                           typeboundary(iw,iB_pole)==bc_asymm) &
                           * bgc(1)%w( ix1, ix2, ix3, iw, igrid )
                   end do
                end do
             end do
          end do

          nbprocs_info%course_nb(inb)%info_send%buffer( 1 + 5 * (i - 1) : 5 * i ) = &
               [neighbor(1,i1,i2,i3,igrid), inc1, inc2, inc3, ibuf_start]
       end do
    end do

#ifdef NOGPUDIRECT
    do inb = 1, nbprocs_info%nbprocs_srl
       
       
    end do
    do inb = 1, nbprocs_info%nbprocs_c
       
       
    end do
#endif

    ! MPI send SRL
    do inb = 1, nbprocs_info%nbprocs_srl
#ifndef NOGPUDIRECT
#ifdef _CRAYFTN
      
#else
      
#endif
#endif
       call mpi_isend_wrapper(nbprocs_info%srl_nb(inb)%send%buffer, &
            nbprocs_info%srl_nb(inb)%send%size, &
            MPI_DOUBLE_PRECISION, nbprocs_info%nbprocs_srl_list(inb), 1, icomm, send_srl_nb(inb), ierrmpi)
       call mpi_isend_wrapper(nbprocs_info%srl_nb(inb)%info_send%buffer, &
            nbprocs_info%srl_nb(inb)%info_send%size, &
            MPI_INTEGER, nbprocs_info%nbprocs_srl_list(inb), 2, icomm, send_srl_nb(nbprocs_info%nbprocs_srl + inb), ierrmpi)
#ifndef NOGPUDIRECT
      
#endif
    end do

    ! MPI send C (send_restrict)
    do inb = 1, nbprocs_info%nbprocs_c
#ifndef NOGPUDIRECT
#ifdef _CRAYFTN
      
#else
      
#endif
#endif
       call mpi_isend_wrapper(nbprocs_info%course_nb(inb)%send%buffer, &
            nbprocs_info%course_nb(inb)%send%size, &
            MPI_DOUBLE_PRECISION, nbprocs_info%nbprocs_c_list(inb), 1, icomm, send_c_nb(inb), ierrmpi)
       call mpi_isend_wrapper(nbprocs_info%course_nb(inb)%info_send%buffer, &
            nbprocs_info%course_nb(inb)%info_send%size, &
            MPI_INTEGER, nbprocs_info%nbprocs_c_list(inb), 2, icomm, send_c_nb(nbprocs_info%nbprocs_c + inb), ierrmpi)
#ifndef NOGPUDIRECT
      
#endif
    end do

    ! fill ghost-cell values of sibling blocks and if neighbor is coarser (f2c)
    ! same process case
     
    do iigrid = 1, igridstail
       do i = 1, 27
          call idecode( i1, i2, i3, i)
          if (skip_direction([ i1,i2,i3 ])) cycle
          igrid=igrids(iigrid)
          iib1=idphyb(1,igrid); iib2=idphyb(2,igrid); iib3=idphyb(3,igrid)
          ipe_neighbor=neighbor(2,i1,i2,i3,igrid)

          if (ipe_neighbor==mype) then
             ineighbor=neighbor(1,i1,i2,i3,igrid)

             select case (neighbor_type(i1,i2,i3,igrid))
             case(neighbor_sibling)

                   ipole=neighbor_pole(i1,i2,i3,igrid)
                   ! across a pole the partner faces us from the same side, so
                   ! that component of the offset keeps its sign
                   n_i1=-i1; n_i2=-i2; n_i3=-i3
                   if (ipole==1) n_i1=i1
                   if (ipole==2) n_i2=i2
                   if (ipole==3) n_i3=i3

                   ixSmin1=ixS_srl_min1(iib1,i1);   ixSmin2=ixS_srl_min2(iib2,i2)
                   ixSmin3=ixS_srl_min3(iib3,i3);   ixSmax1=ixS_srl_max1(iib1,i1)
                   ixSmax2=ixS_srl_max2(iib2,i2);   ixSmax3=ixS_srl_max3(iib3,i3)
                   ixRmin1=ixR_srl_min1(iib1,n_i1); ixRmin2=ixR_srl_min2(iib2,n_i2)
                   ixRmin3=ixR_srl_min3(iib3,n_i3); ixRmax1=ixR_srl_max1(iib1,n_i1)
                   ixRmax2=ixR_srl_max2(iib2,n_i2); ixRmax3=ixR_srl_max3(iib3,n_i3)

                   sbase1=ixSmin1; sstep1=1; sbase2=ixSmin2; sstep2=1
                   sbase3=ixSmin3; sstep3=1; iB_pole=1
                   if (ipole==1) then
                      sbase1=ixSmax1; sstep1=-1; iB_pole=(i1+3)/2
                   else if (ipole==2) then
                      sbase2=ixSmax2; sstep2=-1; iB_pole=2+(i2+3)/2
                   else if (ipole==3) then
                      sbase3=ixSmax3; sstep3=-1; iB_pole=4+(i3+3)/2
                   end if

                    
                   do iw = nwhead, nwtail
                      do ix3=1,ixSmax3-ixSmin3+1
                         do ix2=1,ixSmax2-ixSmin2+1
                            do ix1=1,ixSmax1-ixSmin1+1
                               bg(bgstep)%w(ixRmin1+ix1-1,ixRmin2+ix2-1,ixRmin3+ix3-1,&
                                    iw,ineighbor) = merge(-one,one,ipole/=0 .and. &
                                    typeboundary(iw,iB_pole)==bc_asymm) &
                                    * bg(bgstep)%w(sbase1+sstep1*(ix1-1),sbase2+sstep2*(ix2-1),&
                                    sbase3+sstep3*(ix3-1),iw,igrid)
                            end do
                         end do
                      end do
                   end do

              case(neighbor_coarse)

                 ic1=1+modulo(node(pig1_,igrid)-1,2)
                 ic2=1+modulo(node(pig2_,igrid)-1,2)
                 ic3=1+modulo(node(pig3_,igrid)-1,2)

                 if(.not.(i1==0.or.i1==2*ic1-3).or..not.(i2==0.or.i2==2*ic2-3)&
                      .or..not.(i3==0.or.i3==2*ic3-3)) cycle

                   ipole=neighbor_pole(i1,i2,i3,igrid)
                   ! same rule as srl, expressed in the child-offset index:
                   ! the pole component keeps 2*i and mirrors the child index
                   n_inc1=-2*i1+ic1;n_inc2=-2*i2+ic2;n_inc3=-2*i3+ic3;
                   if (ipole==1) n_inc1=2*i1+(3-ic1)
                   if (ipole==2) n_inc2=2*i2+(3-ic2)
                   if (ipole==3) n_inc3=2*i3+(3-ic3)

                   ixSmin1=ixS_r_min1(iib1,i1);ixSmin2=ixS_r_min2(iib2,i2)
                   ixSmin3=ixS_r_min3(iib3,i3);ixSmax1=ixS_r_max1(iib1,i1)
                   ixSmax2=ixS_r_max2(iib2,i2);ixSmax3=ixS_r_max3(iib3,i3);
                   ixRmin1=ixR_r_min1(iib1,n_inc1);ixRmin2=ixR_r_min2(iib2,n_inc2)
                   ixRmin3=ixR_r_min3(iib3,n_inc3);ixRmax1=ixR_r_max1(iib1,n_inc1)
                   ixRmax2=ixR_r_max2(iib2,n_inc2);ixRmax3=ixR_r_max3(iib3,n_inc3);

                   sbase1=ixSmin1; sstep1=1; sbase2=ixSmin2; sstep2=1
                   sbase3=ixSmin3; sstep3=1; iB_pole=1
                   if (ipole==1) then
                      sbase1=ixSmax1; sstep1=-1; iB_pole=(i1+3)/2
                   else if (ipole==2) then
                      sbase2=ixSmax2; sstep2=-1; iB_pole=2+(i2+3)/2
                   else if (ipole==3) then
                      sbase3=ixSmax3; sstep3=-1; iB_pole=4+(i3+3)/2
                   end if

                    
                   do iw = nwhead, nwtail
                      do ix3=1,ixSmax3-ixSmin3+1
                         do ix2=1,ixSmax2-ixSmin2+1
                            do ix1=1,ixSmax1-ixSmin1+1
                               bg(bgstep)%w(ixRmin1+ix1-1,ixRmin2+ix2-1,&
                                    ixRmin3+ix3-1,iw,ineighbor) = merge(-one,one,&
                                    ipole/=0 .and. typeboundary(iw,iB_pole)==bc_asymm) &
                                    * bgc(1)%w(sbase1+sstep1*(ix1-1),&
                                    sbase2+sstep2*(ix2-1),sbase3+sstep3*(ix3-1),iw, igrid)
                            end do
                         end do
                      end do
                   end do

                end select

            end if

         end do
    end do


    call MPI_WAITALL(nbprocs_info%nbprocs_srl*2, recv_srl_nb, recvstatus_srl_nb, ierrmpi)
    call MPI_WAITALL(nbprocs_info%nbprocs_srl*2, send_srl_nb, sendstatus_srl_nb, ierrmpi)
#ifdef NOGPUDIRECT
    do inb = 1, nbprocs_info%nbprocs_srl
      
      
    end do
#endif

    ! unpack the MPI buffers
    do inb = 1, nbprocs_info%nbprocs_srl
        
       do i = 1, nbprocs_info%srl_nb(inb)%info%nigrids

          igrid       = nbprocs_info%srl_nb(inb)%info_rcv%buffer( 3 * (i - 1) + 1 )
          ienc        = nbprocs_info%srl_nb(inb)%info_rcv%buffer( 3 * (i - 1) + 2 )
          ibuf_start  = nbprocs_info%srl_nb(inb)%info_rcv%buffer( 3 * (i - 1) + 3 )

          ! ienc already encodes our own ghost-region offset: the sender
          ! computed it, mirrored across a pole where there is one
          call idecode( i1, i2, i3, ienc )

          iib1 = idphyb(1,igrid); iib2 = idphyb(2,igrid); iib3 = idphyb(3,igrid)

          ixRmin1=ixR_srl_min1(iib1,i1); ixRmin2=ixR_srl_min2(iib2,i2)
          ixRmin3=ixR_srl_min3(iib3,i3); ixRmax1=ixR_srl_max1(iib1,i1)
          ixRmax2=ixR_srl_max2(iib2,i2); ixRmax3=ixR_srl_max3(iib3,i3)
          Nx1=ixRmax1-ixRmin1+1; Nx2=ixRmax2-ixRmin2+1; Nx3=ixRmax3-ixRmin3+1

           
          do iw = nwhead, nwtail
             do ix3 = ixRmin3, ixRmax3
                do ix2 = ixRmin2, ixRmax2
                   do ix1 = ixRmin1, ixRmax1
                      ! going through a tempval is a workaround for Cray, which gives
                      ! a memory access fault on the GPUs otherwise
                      tempval = nbprocs_info%srl_nb(inb)%rcv%buffer( &
                           ibuf_start &
                           + (ix1-ixRmin1) &
                           + Nx1 * (ix2-ixRmin2) &
                           + Nx1*Nx2 * (ix3-ixRmin3) &
                           + Nx1*Nx2*Nx3 * (iw-nwhead) &
                           )
                      bg(bgstep)%w( ix1, ix2, ix3, iw , igrid) = tempval
                   end do
                end do
             end do
          end do

       end do
    end do

    call MPI_WAITALL(nbprocs_info%nbprocs_f*2, recv_f_nb, recvstatus_f_nb, ierrmpi)
    call MPI_WAITALL(nbprocs_info%nbprocs_c*2, send_c_nb, sendstatus_c_nb, ierrmpi)

#ifdef NOGPUDIRECT
    do inb = 1, nbprocs_info%nbprocs_f
       
       
    end do
#endif

    ! unpack the MPI buffers, fine neighbor, (f_recv), recv_restrict
    do inb = 1, nbprocs_info%nbprocs_f
        
       do i = 1,nbprocs_info%fine_nb(inb)%info%nigrids

          igrid       = nbprocs_info%fine_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 1 )
          inc1        = nbprocs_info%fine_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 2 )
          inc2        = nbprocs_info%fine_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 3 )
          inc3        = nbprocs_info%fine_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 4 )
          ibuf_start  = nbprocs_info%fine_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 5 )

          iib1 = idphyb(1,igrid); iib2 = idphyb(2,igrid); iib3 = idphyb(3,igrid)

          ixRmin1=ixR_r_min1(iib1,inc1); ixRmin2=ixR_r_min2(iib2,inc2)
          ixRmin3=ixR_r_min3(iib3,inc3); ixRmax1=ixR_r_max1(iib1,inc1)
          ixRmax2=ixR_r_max2(iib2,inc2); ixRmax3=ixR_r_max3(iib3,inc3)
          Nx1=ixRmax1-ixRmin1+1; Nx2=ixRmax2-ixRmin2+1; Nx3=ixRmax3-ixRmin3+1

           
          do iw = nwhead, nwtail
             do ix3 = ixRmin3, ixRmax3
                do ix2 = ixRmin2, ixRmax2
                   do ix1 = ixRmin1, ixRmax1
                      ! going through a tempval is a workaround for Cray, which gives
                      ! a memory access fault on the GPUs otherwise
                      tempval = nbprocs_info%fine_nb(inb)%rcv%buffer( &
                           ibuf_start &
                           + (ix1-ixRmin1) &
                           + Nx1 * (ix2-ixRmin2) &
                           + Nx1*Nx2 * (ix3-ixRmin3) &
                           + Nx1*Nx2*Nx3 * (iw-nwhead) &
                           )
                      bg(bgstep)%w( ix1, ix2, ix3, iw, igrid ) = tempval
                   end do
                end do
             end do
          end do

       end do
    end do

    ! MPI receive prolong (neighbor is coarser)
    do inb = 1, nbprocs_info%nbprocs_c
#ifndef NOGPUDIRECT
#ifdef _CRAYFTN
      
#else
      
#endif
#endif
       call mpi_irecv_wrapper(nbprocs_info%course_nb(inb)%rcv%buffer, &
            nbprocs_info%course_nb(inb)%rcv%size, &
            MPI_DOUBLE_PRECISION, nbprocs_info%nbprocs_c_list(inb), 1, icomm, recv_c_nb(inb), ierrmpi)
       call mpi_irecv_wrapper(nbprocs_info%course_nb(inb)%info_rcv%buffer, &
            nbprocs_info%course_nb(inb)%info_rcv%size, &
            MPI_INTEGER, nbprocs_info%nbprocs_c_list(inb), 2, icomm, recv_c_nb(nbprocs_info%nbprocs_c + inb), ierrmpi)
#ifndef NOGPUDIRECT
      
#endif
    end do

    ! fill the F (neighbor is finer) send buffer on GPU (send_prolong)
    do inb = 1, nbprocs_info%nbprocs_f
         
       do i = 1,nbprocs_info%fine_nb(inb)%info%nigrids

          igrid = nbprocs_info%fine_nb(inb)%info%igrid(i)
          inc1 = nbprocs_info%fine_nb(inb)%info%inc1(i)
          inc2 = nbprocs_info%fine_nb(inb)%info%inc2(i)
          inc3 = nbprocs_info%fine_nb(inb)%info%inc3(i)

          ! only inc = 2*i+ic is stored, and neighbor_pole is indexed by i;
          ! ic is 1 or 2, so inc pins i down (0 -> -1, 1 or 2 -> 0, 3 -> +1)
          if (inc1 == 0) then; i1 = -1; else if (inc1 == 3) then; i1 = 1; else; i1 = 0; end if
          if (inc2 == 0) then; i2 = -1; else if (inc2 == 3) then; i2 = 1; else; i2 = 0; end if
          if (inc3 == 0) then; i3 = -1; else if (inc3 == 3) then; i3 = 1; else; i3 = 0; end if
          ipole=neighbor_pole(i1,i2,i3,igrid)

          ibuf_start = nbprocs_info%fine_nb(inb)%info%ibuf_start(i)
          iib1=idphyb(1,igrid); iib2=idphyb(2,igrid); iib3=idphyb(3,igrid)

          ! now fill the data and info buffers
          ixSmin1=ixS_p_min1(iib1,inc1);ixSmin2=ixS_p_min2(iib2,inc2)
          ixSmin3=ixS_p_min3(iib3,inc3);ixSmax1=ixS_p_max1(iib1,inc1)
          ixSmax2=ixS_p_max2(iib2,inc2);ixSmax3=ixS_p_max3(iib3,inc3)
          Nx1=ixSmax1-ixSmin1+1; Nx2=ixSmax2-ixSmin2+1; Nx3=ixSmax3-ixSmin3+1

          ! mirror and sign on the send side, as for srl and restrict above
          sbase1=ixSmin1; sstep1=1; sbase2=ixSmin2; sstep2=1
          sbase3=ixSmin3; sstep3=1; iB_pole=1
          if (ipole==1) then
             sbase1=ixSmax1; sstep1=-1; iB_pole=(i1+3)/2
          else if (ipole==2) then
             sbase2=ixSmax2; sstep2=-1; iB_pole=2+(i2+3)/2
          else if (ipole==3) then
             sbase3=ixSmax3; sstep3=-1; iB_pole=4+(i3+3)/2
          end if

           
          do iw = nwhead, nwtail
             do ix3 = ixSmin3, ixSmax3
                do ix2 = ixSmin2, ixSmax2
                   do ix1 = ixSmin1, ixSmax1
                      nbprocs_info%fine_nb(inb)%send%buffer( &
                           ibuf_start &
                           + sstep1*(ix1-sbase1) &
                           + Nx1 * (sstep2*(ix2-sbase2)) &
                           + Nx1*Nx2 * (sstep3*(ix3-sbase3)) &
                           + Nx1*Nx2*Nx3 * (iw-nwhead) &
                           ) = merge(-one,one,ipole/=0 .and. &
                           typeboundary(iw,iB_pole)==bc_asymm) &
                           * bg(bgstep)%w( ix1, ix2, ix3, iw, igrid )
                   end do
                end do
             end do
          end do

          ! how they will be used at the receiving end.  The plain rule is
          ! n_inc = ic-i, the 0<->3 swap; across a pole n_inc = inc, no swap.
          n_inc1 = inc1; n_inc2 = inc2; n_inc3 = inc3
          if (ipole /= 1) then
             if (n_inc1 == 0) then; n_inc1 = 3; else if (n_inc1 == 3) then; n_inc1 = 0; end if
          end if
          if (ipole /= 2) then
             if (n_inc2 == 0) then; n_inc2 = 3; else if (n_inc2 == 3) then; n_inc2 = 0; end if
          end if
          if (ipole /= 3) then
             if (n_inc3 == 0) then; n_inc3 = 3; else if (n_inc3 == 3) then; n_inc3 = 0; end if
          end if

          nbprocs_info%fine_nb(inb)%info_send%buffer( 1 + 5 * (i - 1) : 5 * i ) = &
               [neighbor_child(1,inc1,inc2,inc3,igrid), n_inc1, n_inc2, n_inc3, ibuf_start]
       end do
    end do

#ifdef NOGPUDIRECT
    do inb = 1, nbprocs_info%nbprocs_f
      
      
    end do
#endif

    ! MPI send F (send_prolong)
    do inb = 1, nbprocs_info%nbprocs_f
#ifndef NOGPUDIRECT
#ifdef _CRAYFTN
      
#else
      
#endif
#endif
       call mpi_isend_wrapper(nbprocs_info%fine_nb(inb)%send%buffer, &
            nbprocs_info%fine_nb(inb)%send%size, &
            MPI_DOUBLE_PRECISION, nbprocs_info%nbprocs_f_list(inb), 1, icomm, send_f_nb(inb), ierrmpi)
       call mpi_isend_wrapper(nbprocs_info%fine_nb(inb)%info_send%buffer, &
            nbprocs_info%fine_nb(inb)%info_send%size, &
            MPI_INTEGER, nbprocs_info%nbprocs_f_list(inb), 2, icomm, send_f_nb(nbprocs_info%nbprocs_f + inb), ierrmpi)
#ifndef NOGPUDIRECT
      
#endif
    end do


    ! fill coarse ghost-cell values of finer neighbors in the same processor
     
    do iigrid=1,igridstail
       do i3=-1,1
          do i2=-1,1
             do i1=-1,1
                if (skip_direction([ i1,i2,i3 ])) cycle
                igrid = igrids(iigrid)
                if (neighbor_type(i1,i2,i3,igrid)==neighbor_fine) then
                   iib1=idphyb(1,igrid);iib2=idphyb(2,igrid);iib3=idphyb(3,igrid);
                   ipole=neighbor_pole(i1,i2,i3,igrid)
                   iB_pole=1
                   if (ipole==1) then
                      iB_pole=(i1+3)/2
                   else if (ipole==2) then
                      iB_pole=2+(i2+3)/2
                   else if (ipole==3) then
                      iB_pole=4+(i3+3)/2
                   end if
                   !  inline of call bc_fill_prolong(igrid,i1,i2,i3,iib1,iib2,iib3) :

                   do ic3 = 1+int((1-i3)/2), 2-int((1+i3)/2)
                      inc3 = 2*i3+ic3
                      do ic2 = 1+int((1-i2)/2), 2-int((1+i2)/2)
                         inc2 = 2*i2+ic2
                         do ic1 = 1+int((1-i1)/2), 2-int((1+i1)/2)
                            inc1 = 2*i1+ic1
                            ipe_neighbor = neighbor_child(2,inc1,inc2,inc3,igrid)
                            if(ipe_neighbor==mype) then
                               ixSmin1=ixS_p_min1(iib1,inc1);ixSmin2=ixS_p_min2(iib2,inc2)
                               ixSmin3=ixS_p_min3(iib3,inc3);ixSmax1=ixS_p_max1(iib1,inc1)
                               ixSmax2=ixS_p_max2(iib2,inc2);ixSmax3=ixS_p_max3(iib3,inc3)

                               ineighbor = neighbor_child(1,inc1,inc2,inc3,igrid)
                               ! plain rule n_inc = ic-i; across a pole the
                               ! child sits at the same offset, n_inc = inc
                               n_i1=-i1; n_i2=-i2; n_i3=-i3
                               n_inc1=ic1+n_i1; n_inc2=ic2+n_i2; n_inc3=ic3+n_i3
                               if (ipole==1) n_inc1=inc1
                               if (ipole==2) n_inc2=inc2
                               if (ipole==3) n_inc3=inc3

                               ixRmin1=ixR_p_min1(iib1,n_inc1)
                               ixRmin2=ixR_p_min2(iib2,n_inc2)
                               ixRmin3=ixR_p_min3(iib3,n_inc3)
                               ixRmax1=ixR_p_max1(iib1,n_inc1)
                               ixRmax2=ixR_p_max2(iib2,n_inc2)
                               ixRmax3=ixR_p_max3(iib3,n_inc3)

                               sbase1=ixSmin1; sstep1=1; sbase2=ixSmin2; sstep2=1
                               sbase3=ixSmin3; sstep3=1
                               if (ipole==1) then
                                  sbase1=ixSmax1; sstep1=-1
                               else if (ipole==2) then
                                  sbase2=ixSmax2; sstep2=-1
                               else if (ipole==3) then
                                  sbase3=ixSmax3; sstep3=-1
                               end if

                                
                               do iw = nwhead, nwtail
                                  do ix3 =0, ixRmax3-ixRmin3
                                     do ix2 = 0, ixRmax2-ixRmin2
                                        do ix1 = 0, ixRmax1-ixRmin1
                                           bgc(1)%w(ixRmin1+ix1,ixRmin2+ix2,&
                                                ixRmin3+ix3,iw,ineighbor) = merge(-one,one,&
                                                ipole/=0 .and. typeboundary(iw,iB_pole)==bc_asymm) &
                                                * bg(bgstep)%w(sbase1+sstep1*ix1,&
                                                sbase2+sstep2*ix2,sbase3+sstep3*ix3,iw, igrid)
                                        end do
                                     end do
                                  end do
                               end do

                            end if

                         end do
                      end do
                   end do

                end if

             end do
          end do
       end do

    end do

    call MPI_WAITALL(nbprocs_info%nbprocs_c*2, recv_c_nb, recvstatus_c_nb, ierrmpi)
    call MPI_WAITALL(nbprocs_info%nbprocs_f*2, send_f_nb, sendstatus_f_nb, ierrmpi)


#ifdef NOGPUDIRECT
    do inb = 1, nbprocs_info%nbprocs_c
      
      
    end do
#endif

    ! unpack the MPI buffers, coarse neighbor, (c_recv), recv_prolong
    do inb = 1, nbprocs_info%nbprocs_c
         
       do i = 1, nbprocs_info%course_nb(inb)%info%nigrids

          igrid       = nbprocs_info%course_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 1 )
          inc1        = nbprocs_info%course_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 2 )
          inc2        = nbprocs_info%course_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 3 )
          inc3        = nbprocs_info%course_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 4 )
          ibuf_start  = nbprocs_info%course_nb(inb)%info_rcv%buffer( 5 * (i - 1) + 5 )

          iib1 = idphyb(1,igrid); iib2 = idphyb(2,igrid); iib3 = idphyb(3,igrid)

          ixRmin1=ixR_p_min1(iib1,inc1); ixRmin2=ixR_p_min2(iib2,inc2)
          ixRmin3=ixR_p_min3(iib3,inc3); ixRmax1=ixR_p_max1(iib1,inc1)
          ixRmax2=ixR_p_max2(iib2,inc2); ixRmax3=ixR_p_max3(iib3,inc3)
          Nx1=ixRmax1-ixRmin1+1; Nx2=ixRmax2-ixRmin2+1; Nx3=ixRmax3-ixRmin3+1

           
          do iw = nwhead, nwtail
             do ix3 = ixRmin3, ixRmax3
                do ix2 = ixRmin2, ixRmax2
                   do ix1 = ixRmin1, ixRmax1
                      ! going through a tempval is a workaround for Cray, which gives
                      ! a memory access fault on the GPUs otherwise
                      tempval = nbprocs_info%course_nb(inb)%rcv%buffer( &
                           ibuf_start &
                           + (ix1-ixRmin1) &
                           + Nx1 * (ix2-ixRmin2) &
                           + Nx1*Nx2 * (ix3-ixRmin3) &
                           + Nx1*Nx2*Nx3 * (iw-nwhead) &
                           )
                      bgc(1)%w( ix1, ix2, ix3, iw, igrid ) = tempval
                   end do
                end do
             end do
          end do

       end do
    end do

    ! do prolongation on the ghost-cell values based on the received coarse values from coarser neighbors (f2c)
     
    do iigrid=1, igridstail
       !      inline variant of call gc_prolong(igrid)
       do i3 = -1, 1
          do i2 = -1, 1
             do i1 = -1, 1
                if (skip_direction([ i1,i2,i3 ])) cycle
                igrid = igrids(iigrid)
                if (neighbor_type(i1,i2,i3,igrid)==neighbor_coarse) then
                   iib1=idphyb(1,igrid);iib2=idphyb(2,igrid);iib3=idphyb(3,igrid)
                   !     inline variant of call bc_prolong(igrid,i1,i2,i3,iib1,iib2,iib3)

                   ixFimin1=ixR_srl_min1(iib1,i1);ixFimin2=ixR_srl_min2(iib2,i2)
                   ixFimin3=ixR_srl_min3(iib3,i3);ixFimax1=ixR_srl_max1(iib1,i1)
                   ixFimax2=ixR_srl_max2(iib2,i2);ixFimax3=ixR_srl_max3(iib3,i3);
                   dxFi1=rnode(rpdx1_,igrid);dxFi2=rnode(rpdx2_,igrid)
                   dxFi3=rnode(rpdx3_,igrid);
                   dxCo1=two*dxFi1;dxCo2=two*dxFi2;dxCo3=two*dxFi3;
                   invdxCo1=1.d0/dxCo1;invdxCo2=1.d0/dxCo2;invdxCo3=1.d0/dxCo3;

                   ! 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
                   xFimin1=rnode(rpxmin1_,igrid)-dble(nghostcells)*dxFi1
                   xFimin2=rnode(rpxmin2_,igrid)-dble(nghostcells)*dxFi2
                   xFimin3=rnode(rpxmin3_,igrid)-dble(nghostcells)*dxFi3;
                   xComin1=rnode(rpxmin1_,igrid)-dble(nghostcells)*dxCo1
                   xComin2=rnode(rpxmin2_,igrid)-dble(nghostcells)*dxCo2
                   xComin3=rnode(rpxmin3_,igrid)-dble(nghostcells)*dxCo3;

                    
                   do ixFi3 = ixFimin3,ixFimax3
                      do ixFi2 = ixFimin2,ixFimax2
                         do ixFi1 = ixFimin1,ixFimax1
                            xFi3=xFimin3+(dble(ixFi3)-half)*dxFi3
                            ixCo3=int((xFi3-xComin3)*invdxCo3)+1
                            xCo3=xComin3+(dble(ixCo3)-half)*dxCo3

                            xFi2=xFimin2+(dble(ixFi2)-half)*dxFi2
                            ixCo2=int((xFi2-xComin2)*invdxCo2)+1
                            xCo2=xComin2+(dble(ixCo2)-half)*dxCo2

                            xFi1=xFimin1+(dble(ixFi1)-half)*dxFi1
                            ixCo1=int((xFi1-xComin1)*invdxCo1)+1
                            xCo1=xComin1+(dble(ixCo1)-half)*dxCo1

                            eta1 = (xFi1-xCo1) * invdxCo1
                            eta2 = (xFi2-xCo2) * invdxCo2
                            eta3 = (xFi3-xCo3) * invdxCo3

                            do iw = nwhead, nwtail
                               do idims = 1, ndim
                                  hxCo1=ixCo1-kr(1,idims)
                                  hxCo2=ixCo2-kr(2,idims)
                                  hxCo3=ixCo3-kr(3,idims)
                                  jxCo1=ixCo1+kr(1,idims)
                                  jxCo2=ixCo2+kr(2,idims)
                                  jxCo3=ixCo3+kr(3,idims)

                                  slopeL = bgc(1)%w(ixCo1,ixCo2,ixCo3,iw,igrid) &
                                       - bgc(1)%w(hxCo1,hxCo2,hxCo3,iw,igrid)

                                  slopeR = bgc(1)%w(jxCo1,jxCo2,jxCo3,iw,igrid) &
                                       - bgc(1)%w(ixCo1,ixCo2,ixCo3,iw,igrid)

                                  slopeC = half * ( slopeR + slopeL )

                                  ! get limited slope
                                  signR = sign(one,slopeR)
                                  signC = sign(one,slopeC)

                                  slope(idims) = signC * max(zero,min(dabs(slopeC),&
                                       signC*slopeL,signC*slopeR))
                               end do

                               ! Interpolate from coarse cell using limited slopes
                               bg(bgstep)%w(ixFi1,ixFi2,ixFi3,iw,igrid) = &
                                    bgc(1)%w(ixCo1,ixCo2,ixCo3,iw, igrid) &
                                    + (slope(1)*eta1) + (slope(2)*eta2) + (slope(3)*eta3)
                            end do

                         end do
                      end do
                   end do

                end if
             end do
          end do
       end do
    end do

    ! Re-fill the physical-boundary ghost cells now that the internal
    ! ghost-cell exchange has run, mirroring upstream MPI-AMRVAC's
    ! post-exchange bc_phys pass. fill_boundary_before_gc fills only the
    ! interior transverse range of each boundary (corner info is not yet
    ! known); fill_boundary_after_gc enlarges it to the corners, which is
    ! where a physical boundary meets a pole face. No current test needs
    ! this second pass, but it keeps corner-ghost handling complete and
    ! matches upstream for setups this fork does not yet exercise.
    if(bcphys.and. .not.stagger_grid) then
        
       do iigrid = 1, igridstail; igrid=igrids(iigrid);
          if (.not.phyboundblock(igrid)) cycle
          call fill_boundary_after_gc(psb(igrid),igrid,time,qdt)
       end do
    end if

    ! modify normal component of magnetic field to fix divB=0
    if(bcphys.and.associated(phys_boundary_adjust)) then
       do iigrid=1,igridstail; igrid=igrids(iigrid);
          if(.not.phyboundblock(igrid)) cycle
          call phys_boundary_adjust(igrid,psb)
       end do
    end if

    time_bc=time_bc+(MPI_WTIME()-time_bcin)

    call nvtxEndRange

  end subroutine getbc