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

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

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

  ! MPI derived datatype to send and receive subarrays of ghost cells to
  ! neighbor blocks in a different processor.
  !
  ! 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
  !
  ! 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, iib^D

    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(:, 1)=ixmmax^d+1-nghostcells
    ixs_srl_max^d(:,-1)=ixmmin^d-1+nghostcells
    ixs_srl_max^d(:, 1)=ixmmax^d
    
    ixs_srl_min^d(-1,0)=1
    ixs_srl_min^d( 0,0)=ixmmin^d
    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( 0,0)=ixmmax^d
    ixs_srl_max^d( 1,0)=ixgmax^d
    ixs_srl_max^d( 2,0)=ixgmax^d
     
    ixr_srl_min^d(:,-1)=1
    ixr_srl_min^d(:, 1)=ixmmax^d+1
    ixr_srl_max^d(:,-1)=nghostcells
    ixr_srl_max^d(:, 1)=ixgmax^d
    
    ixr_srl_min^d(-1,0)=1
    ixr_srl_min^d( 0,0)=ixmmin^d
    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( 0,0)=ixmmax^d
    ixr_srl_max^d( 1,0)=ixgmax^d
    ixr_srl_max^d( 2,0)=ixgmax^d
    
    ixs_r_min^d(:,-1)=ixcommin^d
    ixs_r_min^d(:, 1)=ixcommax^d+1-nghostcells
    ixs_r_max^d(:,-1)=ixcommin^d-1+nghostcells
    ixs_r_max^d(:, 1)=ixcommax^d
    
    ixs_r_min^d(-1,0)=1
    ixs_r_min^d( 0,0)=ixcommin^d
    ixs_r_min^d( 1,0)=ixcommin^d
    ixs_r_max^d(-1,0)=ixcommax^d
    ixs_r_max^d( 0,0)=ixcommax^d
    ixs_r_max^d( 1,0)=ixcogmax^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

    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(:, 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

    ! extend index range to physical boundary
    ixs_p_min^d(-1,1)=1
    ixs_p_max^d( 1,2)=ixgmax^d

    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

    ! extend index range to physical boundary
    ixr_p_min^d(-1,1)=1
    ixr_p_max^d( 1,2)=ixcogmax^d
    \}

  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
    
    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

  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+1
    integer, intent(in)               :: nwbc    ! Number of variables to fill
    logical, intent(in), optional     :: req_diag ! If false, skip diagonal ghost cells

    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
    ! store physical boundary indicating index
    integer :: idphyb(ndim,max_blocks),bindex(ndim)
    integer :: isend_buf(npwbuf), ipwbuf, nghostcellsco,iB
    logical  :: req_diagonal
    type(wbuffer) :: pwbuf(npwbuf)

    double precision :: time_bcin
    ! Stretching grid parameters for coarsened block of the current block

    nwhead=nwstart+1
    nwtail=nwstart+nwbc

    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) then
      do iigrid=1,igridstail; igrid=igrids(iigrid);
         if(.not.phyboundblock(igrid)) cycle
         saveigrid=igrid
         block=>ps(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)%w,ps(igrid)%x,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

    ! 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=>ps(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)%w,ps(igrid)%x,ixg^l,ixm^l,psc(igrid)%w,psc(igrid)%x,&
                            ixcog^l,ixcom^l,igrid,igrid)
       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_coarse)
             call bc_send_restrict
          case (neighbor_sibling)
             call bc_send_srl
          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)
    deallocate(recvstatus,recvrequest)

    call mpi_waitall(isend,sendrequest,sendstatus,ierrmpi)
    do ipwbuf=1,npwbuf
       if (isend_buf(ipwbuf)/=0) deallocate(pwbuf(ipwbuf)%w)
    end do
    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=>ps(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)
    deallocate(recvstatus,recvrequest)

    ! do prolongation on the ghost-cell values received from coarser neighbors 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       saveigrid=igrid
       block=>ps(igrid)
       ^d&iib^d=idphyb(^d,igrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       if (any(neighbor_type(:^d&,igrid)==neighbor_coarse)) then
          {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_prolong
          {end do\}
       end if
    end do
    
    call mpi_waitall(isend,sendrequest,sendstatus,ierrmpi)
    do ipwbuf=1,npwbuf
       if (isend_buf(ipwbuf)/=0) deallocate(pwbuf(ipwbuf)%w)
    end do
    deallocate(sendstatus,sendrequest)

     ! modify normal component of magnetic field to fix divB=0 
    if(bcphys .and. physics_type=='mhd' .and. ndim>1) 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

      !> 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)
           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_buf(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)
           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)) then
          ! 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)%w,&
                       psc(igrid)%x,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)
           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_buf(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)
           end if
        end if

      end subroutine bc_send_restrict

      !> 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)
              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)
              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_buf(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)
              end if
           end if
        {end do\}

      end subroutine bc_send_prolong

      !> 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)
        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)
           end if
        {end do\}

      end subroutine bc_recv_restrict

      !> 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)  
        end if

      end subroutine bc_recv_prolong

      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(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 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-block%dvolume(ixfi^dd) &
                   /sum(block%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_buf(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_buf

      subroutine fix_auxiliary
        use mod_physics, only: phys_get_aux
      
        integer :: ix^L

        do iigrid=1,igridstail; igrid=igrids(iigrid);
          saveigrid=igrid
          block=>ps(igrid)
          call identifyphysbound(ps(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