mod_amr_fct.t

Go to the documentation of this file.

module mod_amr_fct
  implicit none
  private
 
  type facealloc
    double precision, dimension(:^D&), pointer :: face
  end type facealloc
 
  type fake_neighbors
    integer :: igrid
    integer :: ipe
  end type fake_neighbors
 
  type(facealloc), dimension(:,:,:), allocatable, public :: pface
 
  type(fake_neighbors), dimension(:^D&,:,:), allocatable, public :: fine_neighbors
 
  integer, dimension(:,:^D&,:), allocatable, public :: old_neighbor
 
  integer :: itag, isend, irecv
  integer :: nrecv, nsend, ibuf_recv, ibuf_send, ibuf_send_next
  integer, dimension(^ND) :: isize
  integer, dimension(:), allocatable :: recvrequest, sendrequest
  integer, dimension(:,:), allocatable :: recvstatus, sendstatus
  double precision, allocatable :: recvbuffer(:), sendbuffer(:)
 
  public :: store_faces
  public :: comm_faces
  public :: end_comm_faces
  public :: deallocatebfaces
  public :: old_neighbors
  public :: prolong_2nd_stg
  public :: already_fine
 
contains
  !> This subroutine performs a 2nd order prolongation for a staggered field F,
  !> preserving the divergence of the coarse cell.
  !> This is useful for preserving DivF=0.
  !> If DivF=f(x), a different algorithm must be used.
  subroutine prolong_2nd_stg(sCo,sFi,ixCo^Lin,ixFi^Lin,dxCo^D,xComin^D,dxFi^D,xFimin^D,ghost,fine_^Lin)
    use mod_global_parameters
    use mod_physics
 
    logical, intent(in)          :: ghost
    integer, intent(in)          :: ixco^lin, ixfi^lin
    double precision, intent(in) :: dxco^d, xcomin^d, dxfi^d, xfimin^d
    type(state), intent(in)      :: sco
    type(state), intent(inout)   :: sfi
 
    logical, optional :: fine_^lin
    logical           :: fine_^l
 
    double precision :: eta^d, invdxco^d
    integer :: ixco^l,ixfi^l
    integer :: idim1,idim2,ix^de,idim3,ixfis^l,ixgs^l,ixcos^l,ixfisc^l
    integer :: ixcosv^l(1:ndim),ixfisv^l(1:ndim)
    integer :: hxcos^l,jxcos^l,ixcose^l,ixfise^l,hxfisc^l,jxfisc^l,ipxfisc^l,ixcosc^l,imxfisc^l,jpxfisc^l,jmxfisc^l,hpxfisc^l
    integer :: hxfi^l,jxfi^l,hijxfi^l,hjixfi^l,hjjxfi^l
    integer :: iihxfi^l,iijxfi^l,ijhxfi^l,ijjxfi^l,ihixfi^l,ijixfi^l,ihjxfi^l
    integer :: jihxfi^l,jijxfi^l,jjhxfi^l,jjjxfi^l,jhixfi^l,jjixfi^l,jhjxfi^l
    double precision :: bfluxco(sco%ixgs^s,nws),bfluxfi(sfi%ixgs^s,nws)
    double precision :: slopes(sco%ixgs^s,ndim),b_energy_change(ixg^t)
    {^ifthreed
    ! Directional bias, see pdf
    ! These arrays should have ranges ixO^S
    ! Since they are still not known, we give ranges ixG^T,
    ! even though it is wasteful
    double precision :: sigmau(ixg^t),sigmad(ixg^t),sigma(ixg^t,1:ndim), alpha(ixg^t,1:ndim)
    ! Auxiliary arrays for magnetic fluxes
    double precision :: f1(ixg^t),f2(ixg^t),f3(ixg^t),f4(ixg^t)
    }
 
    {^ifoned
    call mpistop("CT prolongation not implemented in 1D. But CT is not needed.")
    }
 
    {^nooned
    ! Note on the indices:
    ! ixCo  Cells where 
    !       divergence-preserving prolongation will be applied.
    ! ixFi  Fine cells which correspond to that extent.
    ! ixCoE For 'expanded', faces that need to be used to calculate
    !       the slopes for the prolongation
    ! ixCos
    ! ixFis
    ! ixCosV And
    ! ixFisV For 'variable', since the ranges depend on the direction of the faces
    ! ixCosC For 'copy',
    ! ixFisC For 'copy', to fill the information from the fine grid needed in the internal faces prolongation step.
    ! At the end, ixFisC is used to copy the assign the magnetic
    ! fields components to the state structure.
 
    fine_^l=.false.;
 
    {if(present(fine_min^din)) fine_min^d=fine_min^din;}
    {if(present(fine_max^din)) fine_max^d=fine_max^din;}
 
    ! When filling ghost cells, ixFi^L are given.
    ! When refining the block, ixCo^L are given.
    if (ghost) then
      ixfi^l=ixfi^lin;
      {ixcomin^d=int((ixfimin^d+nghostcells+1)/2);}
      {ixcomax^d=int((ixfimax^d+nghostcells+1)/2);}
    else
      ixco^l=ixco^lin;
      ixfi^l=ixm^ll;
    end if
 
    ! Expanded range for staggered variables
    ixcosmin^d=ixcomin^d-1;
    ixcosmax^d=ixcomax^d;
 
    ixfismin^d=ixfimin^d-1;
    ixfismax^d=ixfimax^d;
 
 
    associate(wcos=>sco%ws, wfis=>sfi%ws,wco=>sco%w, wfi=>sfi%w)
    ! Assemble general indices
    ixgsmin^d=sfi%ixGsmin^d;
    ixgsmax^d=sfi%ixGsmax^d;
 
    do idim1=1,ndim
      ixcosvmin^d(idim1)=ixcomin^d-kr(^d,idim1);
      ixcosvmax^d(idim1)=ixcomax^d;
      ixfisvmin^d(idim1)=ixfimin^d-kr(^d,idim1);
      ixfisvmax^d(idim1)=ixfimax^d;
    end do
 
    ! Initialize auxiliary arrays at zero
    bfluxco = zero
    bfluxfi = zero 
    slopes  = zero 
 
 
    invdxco^d=1.d0/dxco^d;
 
    ! Fill coarse magnetic flux array
    do idim1=1,ndim
      ! Fill information from parts already at the fine level
      ! First set up indices
      if (ghost) then
        ixfisemin^d=max(1-kr(^d,idim1),ixfisvmin^d(idim1)-2*(1-kr(^d,idim1)));
        ixfisemax^d=min(ixgsmax^d,ixfisvmax^d(idim1)+2*(1-kr(^d,idim1)));
        ixcose^l=int((ixfise^l+nghostcells+1)/2);
      else
        ixcose^l=ixcosv^l(idim1)^ladd(1-kr(idim1,^d));
        ixfise^l=ixfisv^l(idim1)^ladd2*(1-kr(idim1,^d));
      end if
      ! Convert fine fields to fluxes
      bfluxfi(ixfise^s,idim1)=wfis(ixfise^s,idim1)*sfi%surfaceC(ixfise^s,idim1)
      {^iftwod
      idim2=1+mod(idim1,2)
      }
      {^ifthreed
      idim2=1+mod(idim1,3)
      idim3=1+mod(idim1+1,3)
      }
      bfluxco(ixcose^s,idim1) = zero
      ! Add fine fluxes sharing the same fine face
     {do ix^de=0,1\}
        ixfisc^l=ixfise^l+ix2*kr(idim2,^d);
        {^ifthreed 
        ixfisc^l=ixfisc^l+ix3*kr(idim3,^d);
        }
        bfluxco(ixcose^s,idim1)=bfluxco(ixcose^s,idim1)+bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)
     {end do^DE&\}
    end do
 
    ! Omit indices for already refined face, if any
    {
    if (fine_min^d) then
      ixcosvmin^d(^d)=ixcosvmin^d(^d)+1
      ixfisvmin^d(^d)=ixfisvmin^d(^d)+2
    end if
    if (fine_max^d) then
      ixcosvmax^d(^d)=ixcosvmax^d(^d)-1
      ixfisvmax^d(^d)=ixfisvmax^d(^d)-2
    end if
    \}
 
    do idim1=1,ndim
      ixcose^l=ixcosv^l(idim1);
      ! Omit part already refined
      {
      if (^d/=idim1) then
        if ((.not.fine_min^d).or.(.not.ghost)) then
         ixcosemin^d=ixcosvmin^d(idim1)-1
        end if
        if ((.not.fine_max^d).or.(.not.ghost)) then
         ixcosemax^d=ixcosvmax^d(idim1)+1
        end if
      end if
      \}
      ! Fill coarse flux array from coarse field
      bfluxco(ixcose^s,idim1)=wcos(ixcose^s,idim1)*sco%surfaceC(ixcose^s,idim1)
    end do
    ! Finished filling coarse flux array
 
    ! Distribute coarse fluxes among fine fluxes
    ! There are too many loops here, perhaps optimize later
    do idim1=1,ndim
       ixcos^l=ixcosv^l(idim1);
       do idim2=1,ndim
         if(idim1==idim2) cycle
         ! Calculate slope in direction idim2
         ! Set up indices
         jxcos^l=ixcos^l+kr(idim2,^d);
         hxcos^l=ixcos^l-kr(idim2,^d);
         slopes(ixcos^s,idim2)=0.125d0*(bfluxco(jxcos^s,idim1)-bfluxco(hxcos^s,idim1))
    {^iftwod
         do ix2=0,1
           ixfiscmin^d=ixfisvmin^d(idim1)+ix2*kr(^d,idim2);
           ixfiscmax^d=ixfisvmax^d(idim1)+ix2*kr(^d,idim2);
           bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)=half*(bfluxco(ixcos^s,idim1)&
              +(2*ix2-1)*slopes(ixcos^s,idim2))
         end do
    }
    {^ifthreed
         do idim3=1,ndim
           ! Calculate slope in direction idim3
           ! Set up indices
           jxcos^l=ixcos^l+kr(idim3,^d);
           hxcos^l=ixcos^l-kr(idim3,^d);
           slopes(ixcos^s,idim3)=0.125d0*(bfluxco(jxcos^s,idim1)-bfluxco(hxcos^s,idim1))
           if(lvc(idim1,idim2,idim3)<1) cycle
           do ix2=0,1
             do ix3=0,1
              {ixfiscmin^d=ixfisvmin^d(idim1)+ix2*kr(^d,idim2)+ix3*kr(^d,idim3);}
              {ixfiscmax^d=ixfisvmax^d(idim1)+ix2*kr(^d,idim2)+ix3*kr(^d,idim3);}
               bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)=quarter*bfluxco(ixcos^s,idim1)&
                 +quarter*(2*ix2-1)*slopes(ixcos^s,idim2)&
                 +quarter*(2*ix3-1)*slopes(ixcos^s,idim3)
             end do
           end do
         end do
    }
       end do
    end do
 
    ! Calculate interior fine fluxes
    {^ifthreed
    ! Directional bias for nonlinear prolongation
    do idim1=1,ndim
      do idim2=1,ndim
        do idim3=1,ndim
          if (lvc(idim1,idim2,idim3)<1) cycle
          ! Set up indices
          hxfi^l=ixfi^l-kr(idim1,^d);
          jxfi^l=ixfi^l+kr(idim1,^d);
 
          hijxfi^l=hxfi^l+kr(idim3,^d);
          hjixfi^l=hxfi^l+kr(idim2,^d);
          hjjxfi^l=hijxfi^l+kr(idim2,^d);
 
          iihxfi^l=ixfi^l-kr(idim3,^d);
          iijxfi^l=ixfi^l+kr(idim3,^d);
          ihixfi^l=ixfi^l-kr(idim2,^d);
          ihjxfi^l=ihixfi^l+kr(idim3,^d);
          ijixfi^l=ixfi^l+kr(idim2,^d);
          ijhxfi^l=ijixfi^l-kr(idim3,^d);
          ijjxfi^l=ijixfi^l+kr(idim3,^d);
 
          jihxfi^l=jxfi^l-kr(idim3,^d);
          jijxfi^l=jxfi^l+kr(idim3,^d);
          jhixfi^l=jxfi^l-kr(idim2,^d);
          jhjxfi^l=jhixfi^l+kr(idim3,^d);
          jjixfi^l=jxfi^l+kr(idim2,^d);
          jjhxfi^l=jjixfi^l-kr(idim3,^d);
          jjjxfi^l=jjixfi^l+kr(idim3,^d);
 
          sigmau(ixco^s)=&
              abs(bfluxfi(jhixfimin^d:jhixfimax^d:2,idim2))&
            + abs(bfluxfi(jhjxfimin^d:jhjxfimax^d:2,idim2))&
            + abs(bfluxfi(jjixfimin^d:jjixfimax^d:2,idim2))&
            + abs(bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim2))&
            + abs(bfluxfi(jihxfimin^d:jihxfimax^d:2,idim3))&
            + abs(bfluxfi(jijxfimin^d:jijxfimax^d:2,idim3))&
            + abs(bfluxfi(jjhxfimin^d:jjhxfimax^d:2,idim3))&
            + abs(bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim3))
 
          sigmad(ixco^s)=&
              abs(bfluxfi(ihixfimin^d:ihixfimax^d:2,idim2))&
            + abs(bfluxfi(ihjxfimin^d:ihjxfimax^d:2,idim2))&
            + abs(bfluxfi(ijixfimin^d:ijixfimax^d:2,idim2))&
            + abs(bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim2))&
            + abs(bfluxfi(iihxfimin^d:iihxfimax^d:2,idim3))&
            + abs(bfluxfi(iijxfimin^d:iijxfimax^d:2,idim3))&
            + abs(bfluxfi(ijhxfimin^d:ijhxfimax^d:2,idim3))&
            + abs(bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim3))
 
          sigma(ixco^s,idim1)=sigmau(ixco^s)+sigmad(ixco^s)
          where(sigma(ixco^s,idim1)/=zero)
            sigma(ixco^s,idim1)=abs(sigmau(ixco^s)-sigmad(ixco^s))/sigma(ixco^s,idim1)
          elsewhere
            sigma(ixco^s,idim1)=zero
          end where
 
        end do
      end do
    end do
 
    ! Directional bias for Toth-Roe prolongation
    do idim1=1,ndim
      do idim2=1,ndim
        do idim3=1,ndim
          if (lvc(idim1,idim2,idim3)<1) cycle
    !       Nonlinear
    !        alpha(ixCo^S,idim1)=sigma(ixCo^S,idim2)-sigma(ixCo^S,idim3)
    !       Homogeneous
    !        alpha(ixCo^S,idim1)=0.d0
    !       Toth-Roe
            alpha(ixco^s,idim1)=(sco%dx(ixco^s,idim2)-sco%dx(ixco^s,idim3))/(sco%dx(ixco^s,idim2)+sco%dx(ixco^s,idim3))
        end do
      end do
    end do
    }
 
    do idim1=1,ndim
      do idim2=1,ndim
        {^iftwod
        if (idim1==idim2) cycle
        ixfiscmin^d=ixfismin^d+1;
        ixfiscmax^d=ixfismax^d-1;
        jxfisc^l=ixfisc^l+kr(idim1,^d);
        hxfisc^l=ixfisc^l-kr(idim1,^d);
        ipxfisc^l=ixfisc^l+kr(idim2,^d);
        imxfisc^l=ixfisc^l-kr(idim2,^d);
        jpxfisc^l=jxfisc^l+kr(idim2,^d);
        jmxfisc^l=jxfisc^l-kr(idim2,^d);
        hpxfisc^l=hxfisc^l+kr(idim2,^d);
 
        bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)=&
         half*(bfluxfi(jxfiscmin^d:jxfiscmax^d:2,idim1)&
              +bfluxfi(hxfiscmin^d:hxfiscmax^d:2,idim1))&
         -quarter*(bfluxfi(ipxfiscmin^d:ipxfiscmax^d:2,idim2)&
                  -bfluxfi(jpxfiscmin^d:jpxfiscmax^d:2,idim2)&
                  -bfluxfi(imxfiscmin^d:imxfiscmax^d:2,idim2)&
                  +bfluxfi(jmxfiscmin^d:jmxfiscmax^d:2,idim2))
 
        bfluxfi(ipxfiscmin^d:ipxfiscmax^d:2,idim1)=&
         half*(bfluxfi(jpxfiscmin^d:jpxfiscmax^d:2,idim1)&
              +bfluxfi(hpxfiscmin^d:hpxfiscmax^d:2,idim1))&
         -quarter*(bfluxfi(ipxfiscmin^d:ipxfiscmax^d:2,idim2)&
                  -bfluxfi(jpxfiscmin^d:jpxfiscmax^d:2,idim2)&
                  -bfluxfi(imxfiscmin^d:imxfiscmax^d:2,idim2)&
                  +bfluxfi(jmxfiscmin^d:jmxfiscmax^d:2,idim2))
        }
        {^ifthreed
        do idim3=1,ndim
          if (lvc(idim1,idim2,idim3)<1) cycle
          ! Set up indices
          hxfi^l=ixfi^l-kr(idim1,^d);
          jxfi^l=ixfi^l+kr(idim1,^d);
 
          hijxfi^l=hxfi^l+kr(idim3,^d);
          hjixfi^l=hxfi^l+kr(idim2,^d);
          hjjxfi^l=hijxfi^l+kr(idim2,^d);
 
          iihxfi^l=ixfi^l-kr(idim3,^d);
          iijxfi^l=ixfi^l+kr(idim3,^d);
          ihixfi^l=ixfi^l-kr(idim2,^d);
          ihjxfi^l=ihixfi^l+kr(idim3,^d);
          ijixfi^l=ixfi^l+kr(idim2,^d);
          ijhxfi^l=ijixfi^l-kr(idim3,^d);
          ijjxfi^l=ijixfi^l+kr(idim3,^d);
 
          jihxfi^l=jxfi^l-kr(idim3,^d);
          jijxfi^l=jxfi^l+kr(idim3,^d);
          jhixfi^l=jxfi^l-kr(idim2,^d);
          jhjxfi^l=jhixfi^l+kr(idim3,^d);
          jjixfi^l=jxfi^l+kr(idim2,^d);
          jjhxfi^l=jjixfi^l-kr(idim3,^d);
          jjjxfi^l=jjixfi^l+kr(idim3,^d);
 
          ! Prolongation formulas
          f1(ixco^s)=bfluxfi(ihixfimin^d:ihixfimax^d:2,idim2)&
                    -bfluxfi(jhixfimin^d:jhixfimax^d:2,idim2)&
                    -bfluxfi(ijixfimin^d:ijixfimax^d:2,idim2)&
                    +bfluxfi(jjixfimin^d:jjixfimax^d:2,idim2)
 
          f2(ixco^s)=bfluxfi(ihjxfimin^d:ihjxfimax^d:2,idim2)&
                    -bfluxfi(jhjxfimin^d:jhjxfimax^d:2,idim2)&
                    -bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim2)&
                    +bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim2)
 
          f3(ixco^s)=bfluxfi(iihxfimin^d:iihxfimax^d:2,idim3)&
                    -bfluxfi(jihxfimin^d:jihxfimax^d:2,idim3)&
                    -bfluxfi(iijxfimin^d:iijxfimax^d:2,idim3)&
                    +bfluxfi(jijxfimin^d:jijxfimax^d:2,idim3)
 
          f4(ixco^s)=bfluxfi(ijhxfimin^d:ijhxfimax^d:2,idim3)&
                    -bfluxfi(jjhxfimin^d:jjhxfimax^d:2,idim3)&
                    -bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim3)&
                    +bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim3)
 
          bfluxfi(ixfimin^d:ixfimax^d:2,idim1)=&
             half*(bfluxfi(hxfimin^d:hxfimax^d:2,idim1)+bfluxfi(jxfimin^d:jxfimax^d:2,idim1))&
             +6.25d-2*((3.d0+alpha(ixco^s,idim2))*f1(ixco^s)&
                      +(1.d0-alpha(ixco^s,idim2))*f2(ixco^s)&
                      +(3.d0-alpha(ixco^s,idim3))*f3(ixco^s)&
                      +(1.d0+alpha(ixco^s,idim3))*f4(ixco^s))
 
          bfluxfi(ijixfimin^d:ijixfimax^d:2,idim1)=&
             half*(bfluxfi(hjixfimin^d:hjixfimax^d:2,idim1)+bfluxfi(jjixfimin^d:jjixfimax^d:2,idim1))&
             +6.25d-2*((3.d0+alpha(ixco^s,idim2))*f1(ixco^s)&
                      +(1.d0-alpha(ixco^s,idim2))*f2(ixco^s)&
                      +(1.d0+alpha(ixco^s,idim3))*f3(ixco^s)&
                      +(3.d0-alpha(ixco^s,idim3))*f4(ixco^s))
 
          bfluxfi(iijxfimin^d:iijxfimax^d:2,idim1)=&
             half*(bfluxfi(hijxfimin^d:hijxfimax^d:2,idim1)+bfluxfi(jijxfimin^d:jijxfimax^d:2,idim1))&
             +6.25d-2*((1.d0-alpha(ixco^s,idim2))*f1(ixco^s)&
                      +(3.d0+alpha(ixco^s,idim2))*f2(ixco^s)&
                      +(3.d0-alpha(ixco^s,idim3))*f3(ixco^s)&
                      +(1.d0+alpha(ixco^s,idim3))*f4(ixco^s))
 
          bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim1)=&
             half*(bfluxfi(hjjxfimin^d:hjjxfimax^d:2,idim1)+bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim1))&
             +6.25d-2*((1.d0-alpha(ixco^s,idim2))*f1(ixco^s)&
                      +(3.d0+alpha(ixco^s,idim2))*f2(ixco^s)&
                      +(1.d0+alpha(ixco^s,idim3))*f3(ixco^s)&
                      +(3.d0-alpha(ixco^s,idim3))*f4(ixco^s))
        end do
        }
      end do
    end do
 
    ! Go back to magnetic fields
    do idim1=1,ndim
      ixfiscmax^d=ixfimax^d;
      ixfiscmin^d=ixfimin^d-kr(^d,idim1);
      where(sfi%surfaceC(ixfisc^s,idim1)/=zero)
        wfis(ixfisc^s,idim1)=bfluxfi(ixfisc^s,idim1)/sfi%surfaceC(ixfisc^s,idim1)
      elsewhere
        wfis(ixfisc^s,idim1)=zero
      end where
    end do
 
    if(phys_total_energy.and. .not.prolongprimitive) then
      b_energy_change(ixfi^s)=0.5d0*sum(wfi(ixfi^s,iw_mag(:))**2,dim=ndim+1)
    end if
    call phys_face_to_center(ixfi^l,sfi)
    if(phys_total_energy.and. .not.prolongprimitive) then
      b_energy_change(ixfi^s)=0.5d0*sum(wfi(ixfi^s,iw_mag(:))**2,dim=ndim+1)-&
        b_energy_change(ixfi^s)
      wfi(ixfi^s,iw_e)=wfi(ixfi^s,iw_e)+b_energy_change(ixfi^s)
    end if
 
    end associate
 
    } ! END NOONED
  end subroutine prolong_2nd_stg
 
  !> To achive consistency and thus conservation of divergence,
  !> when refining a block we take into account the faces of the
  !> already fine neighbours, if any. This routine stores them.
  subroutine store_faces
    use mod_forest, only: refine 
    use mod_global_parameters
    integer :: igrid, iigrid, idims, iside, ineighbor, ipe_neighbor
    integer :: nx^d, i^d, ic^d, inc^d
 
    if (npe>1) then
      nsend_fc=0
      nrecv_fc=0
    end if
 
    ! Size of the block face
    nx^d=ixmhi^d-ixmlo^d+1;
 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ! Check whether it is necessary to store any block face, i.e.
      ! if any coarser neighbour is going to be refined 
     {do iside=1,2
        i^dd=kr(^dd,^d)*(2*iside-3);
        if (neighbor_pole(i^dd,igrid)/=0) cycle
        if (neighbor_type(i^dd,igrid)==neighbor_coarse) then
          ineighbor   =neighbor(1,i^dd,igrid)
          ipe_neighbor=neighbor(2,i^dd,igrid)
          if (refine(ineighbor,ipe_neighbor)) then
            allocate(pface(iside,^d,igrid)%face(1^d%1:nx^dd))
            !! Store the faces
            if (iside==1) then !! left side
              pface(iside,^d,igrid)%face(1^d%1:nx^dd)=&
                ps(igrid)%ws(ixmlo^d-1^d%ixM^t,^d)
            else !! right side
              pface(iside,^d,igrid)%face(1^d%1:nx^dd)=&
                ps(igrid)%ws(ixmhi^d^d%ixM^t,^d)
            end if
            if (ipe_neighbor/=mype) nsend_fc(^d)=nsend_fc(^d)+1
          end if
        end if
      end do\}
 
      ! If a grid is going to be refined,
      ! remember what are its neighbours.
      if (refine(igrid,mype)) then
        ! Initialize neighbour array
        fine_neighbors(:^d&,:,igrid)%igrid=-1
        fine_neighbors(:^d&,:,igrid)%ipe=-1
        do idims=1,ndim
          do iside=1,2
            i^d=kr(^d,idims)*(2*iside-3);
            if (neighbor_pole(i^d,igrid)/=0) cycle
            if (neighbor_type(i^d,igrid)==neighbor_fine) then
             {do ic^db=1+int((1+i^db)/2),2-int((1-i^db)/2)
                inc^db=ic^db+i^db\}
                ineighbor=neighbor_child(1,inc^d,igrid)
                ipe_neighbor=neighbor_child(2,inc^d,igrid)
 
                fine_neighbors(ic^d,idims,igrid)%igrid= ineighbor
                fine_neighbors(ic^d,idims,igrid)%ipe=ipe_neighbor
 
                if (ipe_neighbor/=mype) nrecv_fc(idims)=nrecv_fc(idims)+1
             {end do\}
            end if
          end do
        end do
      end if
 
    end do
 
  end subroutine store_faces
 
  !> When refining a coarse block with fine neighbours, it is necessary
  !> prolong consistently with the already fine faces.
  !> This routine takes care of the communication of such faces.
  subroutine comm_faces
    use mod_forest, only: refine
    use mod_global_parameters
 
    integer                   :: iigrid,igrid,ineighbor,ipe_neighbor
    integer                   :: idims,iside,i^d,ic^d,inc^d,nx^d
    integer                   :: recvsize, sendsize
 
    ! Communicate the block faces to achieve consistency when refining
    ! Initialize communication structures
 
    nrecv=0
    nsend=0
    recvsize=0
    sendsize=0
 
    do idims=1,ndim
       select case (idims)
       {case (^d)
          nrecv=nrecv+nrecv_fc(^d)
          nsend=nsend+nsend_fc(^d)
          nx^d=1^d%nx^dd=ixmhi^dd-ixmlo^dd+1;
          isize(^d)={nx^dd*}
          recvsize=recvsize+nrecv_fc(^d)*isize(^d)
          sendsize=sendsize+nsend_fc(^d)*isize(^d)
       \}
       end select
    end do
 
    if (nrecv>0) then
    ! Allocate receive buffer
      allocate(recvbuffer(recvsize),recvstatus(mpi_status_size,nrecv), &
               recvrequest(nrecv))
      recvrequest=mpi_request_null
      ibuf_recv=1
      irecv=0
 
    ! Receive
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        if (refine(igrid,mype)) then
         {do ic^db=1,2\}
            ! Only one of the sides will be necessary,
            ! so we do the loop only over dimensions, instead of
            ! over dimensions and sizes as in the routines
            ! old_neighbors and already_fine.
            do idims=1,ndim
              ipe_neighbor=fine_neighbors(ic^d,idims,igrid)%ipe
              ineighbor   =fine_neighbors(ic^d,idims,igrid)%igrid
              if (ineighbor>0.and.ipe_neighbor/=mype) then
               {if (idims==^d) iside=ic^d\}
                !!! Check indices
                i^d=kr(^d,idims)*(2*iside-3);
                if (neighbor_pole(i^d,igrid)/=0) cycle
                inc^d=ic^d+i^d;
                irecv=irecv+1
                itag=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
 
                call mpi_irecv(recvbuffer(ibuf_recv),isize(idims), &
                               mpi_double_precision,ipe_neighbor,itag, &
                               icomm,recvrequest(irecv),ierrmpi)
                ibuf_recv=ibuf_recv+isize(idims)
              end if
            end do
         {end do\}
        end if
      end do
 
    end if
 
    if (nsend>0) then
    ! Allocate send buffer
      allocate(sendbuffer(sendsize),sendstatus(mpi_status_size,nsend),sendrequest(nsend))
      sendrequest=mpi_request_null
      isend=0
      ibuf_send=1
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        ! Check whether it is necessary to store any block face, i.e.
        ! if any coarser neighbour is going to be refined 
       {do iside=1,2
          i^dd=kr(^dd,^d)*(2*iside-3);
          ! When there is a pole, faces are always zero and this is not necessary
          if (neighbor_pole(i^dd,igrid)/=0) cycle
          if (neighbor_type(i^dd,igrid)==neighbor_coarse) then
            ineighbor   =neighbor(1,i^dd,igrid)
            ipe_neighbor=neighbor(2,i^dd,igrid)
            if (refine(ineighbor,ipe_neighbor)) then
              if (ipe_neighbor/=mype) then
                ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
                inc^d=-2*i^d+ic^d^d%inc^dd=ic^dd;
                itag=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
                isend=isend+1
                ibuf_send_next=ibuf_send+isize(^d)
                sendbuffer(ibuf_send:ibuf_send_next-1)=&
                reshape(pface(iside,^d,igrid)%face,(/isize(^d)/))
                call mpi_isend(sendbuffer(ibuf_send),isize(^d), &
                               mpi_double_precision,ipe_neighbor,itag, &
                               icomm,sendrequest(isend),ierrmpi)
                ibuf_send=ibuf_send_next
              end if
            end if
          end if
        end do\}
      end do
    end if
 
    ! Waitalls
    if (nrecv>0) then
       call mpi_waitall(nrecv,recvrequest,recvstatus,ierrmpi)
       deallocate(recvstatus,recvrequest)
       ibuf_recv=1
    end if
 
    if (nsend>0) then
       call mpi_waitall(nsend,sendrequest,sendstatus,ierrmpi)
       deallocate(sendbuffer,sendstatus,sendrequest)
    end if
 
  end subroutine comm_faces
 
  subroutine end_comm_faces
    use mod_global_parameters
    ! Deallocate receive buffer
    if (nrecv>0) deallocate(recvbuffer)
  end subroutine end_comm_faces
 
  subroutine deallocatebfaces
    use mod_global_parameters
    integer :: igrid, iigrid, iside
 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
     {do iside=1,2
        if (associated(pface(iside,^d,igrid)%face)) then
          deallocate(pface(iside,^d,igrid)%face)
        end if
      end do\}
    end do
 
  end subroutine deallocatebfaces
 
  subroutine old_neighbors(child_igrid,child_ipe,igrid,ipe)
    use mod_global_parameters
    integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
    integer, intent(in) :: igrid, ipe
    integer :: iside, i^d, ic^d
 
    {do ic^db=1,2\}
      old_neighbor(:,:^d&,child_igrid(ic^d))=-1
     {do iside=1,2
        if (ic^d==iside) then
          i^dd=kr(^dd,^d)*(2*iside-3);
          old_neighbor(1,i^dd,child_igrid(ic^dd))=&
             fine_neighbors(ic^dd,^d,igrid)%igrid
          old_neighbor(2,i^dd,child_igrid(ic^dd))=&
             fine_neighbors(ic^dd,^d,igrid)%ipe
        end if
      end do\}
    {end do\}
 
  end subroutine old_neighbors
 
  !> This routine fills the fine faces before prolonging.
  !> It is the face equivalent of fix_conserve 
  subroutine already_fine(sFi,ichild,fine_^L)
    use mod_forest
    use mod_global_parameters
    type(tree_node_ptr) :: tree
    type(state) :: sfi
    integer, intent(in) :: ichild
    logical :: fine_^l
 
    integer :: ineighbor,ipe_neighbor,ibufnext
    integer :: iside,iotherside,i^d,nx^d
 
    ! Size of the block face
    nx^d=ixmhi^d-ixmlo^d+1;
 
    ! Initialise everything to zero and false
    fine_min^d=.false.;
    fine_max^d=.false.;
    sfi%ws=zero
 
    {^nooned
    ! This face communication is not needed in 1D
   {do iside=1,2
      i^dd=kr(^dd,^d)*(2*iside-3);
      ! This is not necessary at the pole.
      ! We are inside a loop over the children, so the grid index is ichild
      if (neighbor_pole(i^dd,ichild)/=0) cycle
 
      ! Get old ipe and igrid of neighbour from array fake_neighbor
      ! Then plug it into the structure pfaces and get the faces 
      ineighbor   =old_neighbor(1,i^dd,ichild)
      ipe_neighbor=old_neighbor(2,i^dd,ichild)
 
      iotherside=3-iside
      if (ineighbor>0) then
        if (iside==1) then ! left side
          if (ipe_neighbor==mype) then
            sfi%ws(ixmlo^d-1^d%ixM^t,^d)=pface(iotherside,^d,ineighbor)%face(1^d%1:nx^dd)
          else
            ibufnext=ibuf_recv+isize(^d)
            sfi%ws(ixmlo^d-1^d%ixM^t,^d)=reshape(&
                 source=recvbuffer(ibuf_recv:ibufnext-1),&
                 shape=shape(sfi%ws(ixmlo^d-1^d%ixM^t,^d)))
            ibuf_recv=ibufnext
          end if
          fine_min^d=.true.
        else ! right side
          if (ipe_neighbor==mype) then
            sfi%ws(ixmhi^d^d%ixM^t,^d)=pface(iotherside,^d,ineighbor)%face(1^d%1:nx^dd)
          else
            ibufnext=ibuf_recv+isize(^d)
            sfi%ws(ixmhi^d^d%ixM^t,^d)=reshape(&
                 source=recvbuffer(ibuf_recv:ibufnext-1),&
                 shape=shape(sfi%ws(ixmhi^d^d%ixM^t,^d)))
            ibuf_recv=ibufnext
          end if
          fine_max^d=.true.
        end if
      end if
    end do\}
    }
  end subroutine already_fine
 
end module mod_amr_fct