amr_solution_node.t

Go to the documentation of this file.

!> Get first available igrid on processor ipe
integer function getnode(ipe)
  use mod_forest, only: igrid_inuse
  use mod_global_parameters
 
  integer, intent(in) :: ipe
  integer :: igrid, igrid_available
 
  igrid_available=0
 
  do igrid=1,max_blocks
     if (igrid_inuse(igrid,ipe)) cycle
 
     igrid_available=igrid
     exit
  end do
 
  if (igrid_available == 0) then
     getnode = -1
     print *, "Current maximum number of grid blocks:", max_blocks
     call mpistop("Insufficient grid blocks; increase max_blocks in meshlist")
  else
     getnode=igrid_available
     igrid_inuse(igrid,ipe)=.true.
  end if
 
  if (ipe==mype) then
     ! initialize nodal block
     node(1:nodehi,getnode) = 0
     rnode(1:rnodehi,getnode) = zero
  end if
 
end function getnode
 
! put igrid on processor ipe to be not in use
subroutine putnode(igrid,ipe)
  use mod_forest
  implicit none
 
  integer, intent(in) :: igrid, ipe
 
  igrid_inuse(igrid,ipe)=.false.
 
end subroutine putnode
 
!> allocate arrays on igrid node
subroutine alloc_node(igrid)
  use mod_forest
  use mod_global_parameters
  use mod_geometry
  use mod_usr_methods, only: usr_set_surface
  use mod_physics, only: phys_set_equi_vars
 
  integer, intent(in) :: igrid
 
  integer :: level, ig^D, ign^D, ixCoG^L, ix, i^D
  integer :: imin, imax, index, igCo^D, ixshift, offset, ifirst
  integer :: icase, ixGext^L
  double precision :: dx^D, summeddx, sizeuniformpart^D
  double precision :: xext(ixGlo^D-1:ixGhi^D+1,1:ndim)
  double precision :: delx_ext(ixGlo1-1:ixGhi1+1)
  double precision :: exp_factor_ext(ixGlo1-1:ixGhi1+1),del_exp_factor_ext(ixGlo1-1:ixGhi1+1),exp_factor_primitive_ext(ixGlo1-1:ixGhi1+1)
  double precision :: xc(ixGlo1:ixGhi1),delxc(ixGlo1:ixGhi1)
  double precision :: exp_factor_coarse(ixGlo1:ixGhi1),del_exp_factor_coarse(ixGlo1:ixGhi1),exp_factor_primitive_coarse(ixGlo1:ixGhi1)
 
  ixcogmin^d=1;
  ixcogmax^d=(ixghi^d-2*nghostcells)/2+2*nghostcells;
 
  icase=mod(nghostcells,2)
  if(stagger_grid) icase=1
  select case(icase)
    case(0)
      ixgext^l=ixg^ll;
    case(1)
      ! for ghost cell related prolongations, we need
      ! an extra layer with known volumes and dx-intervals
      ! in case the number of ghost cells is odd
      ixgext^l=ixg^ll^ladd1;
    case default
      call mpistop("no such case")
  end select
 
  ! set level information
  level=igrid_to_node(igrid,mype)%node%level
 
  if(.not. allocated(ps(igrid)%w)) then
 
    ! allocate arrays for solution and space
    call alloc_state(igrid, ps(igrid), ixg^ll, ixgext^l, .true.)
    ! allocate arrays for one level coarser solution
    call alloc_state_coarse(igrid, psc(igrid), ixcog^l, ixcog^l)
    if(.not.convert) then
      ! allocate arrays for old solution
      call alloc_state(igrid, pso(igrid), ixg^ll, ixgext^l, .false.)
      ! allocate arrays for temp solution 1
      call alloc_state(igrid, ps1(igrid), ixg^ll, ixgext^l, .false.)
 
      ! allocate temporary solution space
      select case (t_integrator)
      case(ssprk3,ssprk4,imex_midpoint,imex_trapezoidal,imex_222)
        call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
      case(rk3_bt,rk4,ssprk5,imex_cb3a)
        call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
        call alloc_state(igrid, ps3(igrid), ixg^ll, ixgext^l, .false.)
      case(imex_ars3,imex_232)
        call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
        call alloc_state(igrid, ps3(igrid), ixg^ll, ixgext^l, .false.)
        call alloc_state(igrid, ps4(igrid), ixg^ll, ixgext^l, .false.)
      end select
    end if
 
  end if
 
  ! avoid dividing by zero rho in skipped corner ghostcells when phys_req_diagonal=F
  ps(igrid)%w(:^d&,1)=1.d0
  ps(igrid)%level=level
  psc(igrid)%level=level-1
  ! avoid dividing by zero rho in skipped corner ghostcells when phys_req_diagonal=F
  psc(igrid)%w(:^d&,1)=1.d0
  if(phys_trac) ps(igrid)%special_values=0.d0
  if(.not.convert) then
    pso(igrid)%level=level
    ps1(igrid)%level=level
    select case (t_integrator)
    case(ssprk3,ssprk4,imex_midpoint,imex_trapezoidal,imex_222)
      ps2(igrid)%level=level
    case(rk3_bt,rk4,ssprk5,imex_cb3a)
      ps2(igrid)%level=level
      ps3(igrid)%level=level
    case(imex_ars3,imex_232)
      ps2(igrid)%level=level
      ps3(igrid)%level=level
      ps4(igrid)%level=level
    end select
  end if
 
  ! block pointer to current block
  block=>ps(igrid)
 
  ig^d=igrid_to_node(igrid,mype)%node%ig^d;
 
  node(plevel_,igrid)=level
  ^d&node(pig^d_,igrid)=ig^d\
 
  ! set dx information
  ^d&rnode(rpdx^d_,igrid)=dx(^d,level)\
  dxlevel(:)=dx(:,level)
 
  ! uniform cartesian case as well as all unstretched coordinates
  ! determine the minimal and maximal corners
  ^d&rnode(rpxmin^d_,igrid)=xprobmin^d+dble(ig^d-1)*dg^d(level)\
  ^d&rnode(rpxmax^d_,igrid)=xprobmin^d+dble(ig^d)*dg^d(level)\
 {if(rnode(rpxmax^d_,igrid)>xprobmax^d) rnode(rpxmax^d_,igrid)=xprobmax^d\}
 
  ^d&dx^d=rnode(rpdx^d_,igrid)\
 {do ix=ixglo^d,ixmhi^d-nghostcells
    ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)+(dble(ix-nghostcells)-half)*dx^d
  end do\}
 ! update overlap cells of neighboring blocks in the same way to get the same values
 {do ix=ixmhi^d-nghostcells+1,ixghi^d
    ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmax^d_,igrid)+(dble(ix-ixmhi^d)-half)*dx^d
  end do\}
 
  ^d&dx^d=2.0d0*rnode(rpdx^d_,igrid)\
 {do ix=ixcogmin^d,ixcogmax^d
    psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmin^d_,igrid)+(dble(ix-nghostcells)-half)*dx^d
  end do\}
 
  ^d&ps(igrid)%dx(ixgext^s,^d)=rnode(rpdx^d_,igrid);
  ^d&psc(igrid)%dx(ixcog^s,^d)=2.0d0*rnode(rpdx^d_,igrid);
  ^d&dx^d=rnode(rpdx^d_,igrid)\
 {do ix=ixgextmin^d,ixgextmax^d
    xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)+(dble(ix-nghostcells)-half)*dx^d
  end do\}
 
  if(any(stretched_dim)) then
   {if(stretch_type(^d) == stretch_uni)then
      imin=(ig^d-1)*block_nx^d
      imax=ig^d*block_nx^d
      rnode(rpxmin^d_,igrid)=xprobmin^d+dxfirst_1mq(level,^d) &
                   *(1.0d0-qstretch(level,^d)**imin)
      rnode(rpxmax^d_,igrid)=xprobmin^d+dxfirst_1mq(level,^d) &
                   *(1.0d0-qstretch(level,^d)**imax)
      ! fix possible out of bound due to precision
      if(rnode(rpxmax^d_,igrid)>xprobmax^d) rnode(rpxmax^d_,igrid)=xprobmax^d
      ixshift=(ig^d-1)*block_nx^d-nghostcells
      do ix=ixgextmin^d,ixgextmax^d
        index=ixshift+ix
        ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(index-1)
      enddo
      igco^d=(ig^d-1)/2
      ixshift=igco^d*block_nx^d+(1-modulo(ig^d,2))*block_nx^d/2-nghostcells
      do ix=ixcogmin^d,ixcogmax^d
        index=ixshift+ix
        psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(index-1)
        psc(igrid)%x(ix^d%ixCoG^s,^d)=xprobmin^d+dxfirst_1mq(level-1,^d)&
                                *(1.0d0-qstretch(level-1,^d)**(index-1))&
                    + 0.5d0*dxfirst(level-1,^d)*qstretch(level-1,^d)**(index-1)
      end do
      ! now that dx and grid boundaries are known: fill cell centers
      ifirst=nghostcells+1
      ! first fill the mesh
      summeddx=0.0d0
      do ix=ixmlo^d,ixmhi^d
        ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
        summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
      enddo
      ! then ghost cells to left
      summeddx=0.0d0
      do ix=nghostcells,1,-1
        ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
        summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
      enddo
      ! then ghost cells to right
      summeddx=0.0d0
      do ix=ixghi^d-nghostcells+1,ixghi^d
        ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
        summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
      enddo
      select case(icase)
        case(0)
          ! if even number of ghost cells: xext is just copy of local x
          xext(ixgext^s,^d)=ps(igrid)%x(ixgext^s,^d)
        case(1)
          ! if uneven number of ghost cells: extra layer left/right
          summeddx=0.0d0
          do ix=ixmlo^d,ixmhi^d
            xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
           summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
          enddo
          ! then ghost cells to left
          summeddx=0.0d0
          do ix=nghostcells,ixgextmin^d,-1
            xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
            summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
          enddo
          ! then ghost cells to right
          summeddx=0.0d0
          do ix=ixghi^d-nghostcells+1,ixgextmax^d
             xext(ix^d%ixGext^s,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
             summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
          enddo
        case default
          call mpistop("no such case")
      end select
     endif\}
    {if(stretch_type(^d) == stretch_symm)then
       ! here we distinguish three kinds of grid blocks
       ! depending on their ig-index, set per level
       !      the first n_stretchedblocks/2  will stretch to the left
       !      the middle ntotal-n_stretchedblocks will be uniform
       !      the last  n_stretchedblocks/2  will stretch to the right
       if(ig^d<=nstretchedblocks(level,^d)/2)then
         ! stretch to the left
         offset=block_nx^d*nstretchedblocks(level,^d)/2
         imin=(ig^d-1)*block_nx^d
         imax=ig^d*block_nx^d
         rnode(rpxmin^d_,igrid)=xprobmin^d+xstretch^d-dxfirst_1mq(level,^d) &
                                    *(1.0d0-qstretch(level,^d)**(offset-imin))
         rnode(rpxmax^d_,igrid)=xprobmin^d+xstretch^d-dxfirst_1mq(level,^d) &
                                    *(1.0d0-qstretch(level,^d)**(offset-imax))
         ! fix possible out of bound due to precision
         if(rnode(rpxmin^d_,igrid)<xprobmin^d) rnode(rpxmin^d_,igrid)=xprobmin^d
         ixshift=(ig^d-1)*block_nx^d-nghostcells
         do ix=ixgextmin^d,ixgextmax^d
           index=ixshift+ix
           ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(offset-index)
         enddo
         ixshift=(nstretchedblocks(level,^d)/2-ig^d)*(block_nx^d/2)+block_nx^d/2+nghostcells
         do ix=ixcogmin^d,ixcogmax^d
           index=ixshift-ix
           psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**index
         enddo
         ! last block: to modify ghost cells!!!
         if(ig^d==nstretchedblocks(level,^d)/2)then
           if(ng^d(level)==nstretchedblocks(level,^d))then
             ! if middle blocks do not exist then use symmetry
             do ix=ixghi^d-nghostcells+1,ixgextmax^d
                ps(igrid)%dx(ix^d%ixGext^s,^d)= &
                ps(igrid)%dx(2*(ixghi^d-nghostcells)+1-ix^d%ixGext^s,^d)
             enddo
             do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
                psc(igrid)%dx(ix^d%ixCoG^s,^d)= &
                psc(igrid)%dx(2*(ixcogmax^d-nghostcells)+1-ix^d%ixCoG^s,^d)
             enddo
           else
             ! if middle blocks exist then use same as middle blocks:
             do ix=ixghi^d-nghostcells+1,ixgextmax^d
                ps(igrid)%dx(ix^d%ixGext^s,^d)=dxmid(level,^d)
             enddo
             do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
                psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxmid(level-1,^d)
             enddo
           endif
         endif
         ! first block: make ghost cells symmetric (to allow periodicity)
         if(ig^d==1)then
           do ix=ixgextmin^d,nghostcells
             ps(igrid)%dx(ix^d%ixGext^s,^d)=ps(igrid)%dx(2*nghostcells+1-ix^d%ixGext^s,^d)
           enddo
           do ix=1,nghostcells
             psc(igrid)%dx(ix^d%ixCoG^s,^d)=psc(igrid)%dx(2*nghostcells+1-ix^d%ixCoG^s,^d)
           enddo
         endif
       else
         if(ig^d<=ng^d(level)-nstretchedblocks(level,^d)/2) then
           ! keep uniform
           ps(igrid)%dx(ixgext^s,^d)=dxmid(level,^d)
           psc(igrid)%dx(ixcog^s,^d)=dxmid(level-1,^d)
           rnode(rpxmin^d_,igrid)=xprobmin^d+xstretch^d+(ig^d-nstretchedblocks(level,^d)/2-1)*block_nx^d*dxmid(level,^d)
           rnode(rpxmax^d_,igrid)=xprobmin^d+xstretch^d+(ig^d-nstretchedblocks(level,^d)/2)  *block_nx^d*dxmid(level,^d)
           ! first and last block: to modify the ghost cells!!!
           if(ig^d==nstretchedblocks(level,^d)/2+1)then
             do ix=ixgextmin^d,nghostcells
               ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(nghostcells-ix)
             enddo
             do ix=1,nghostcells
               psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(nghostcells-ix)
             enddo
           endif
           if(ig^d==ng^d(level)-nstretchedblocks(level,^d))then
             do ix=ixghi^d-nghostcells+1,ixgextmax^d
               ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(ix-block_nx^d-nghostcells-1)
             enddo
             do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
               psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(ix-ixcogmax^d+nghostcells-1)
             enddo
           endif
         else
           ! stretch to the right
           offset=block_nx^d*(ng^d(level)-nstretchedblocks(level,^d)/2)
           sizeuniformpart^d=dxmid(1,^d)*(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d)
           imin=(ig^d-1)*block_nx^d-offset
           imax=ig^d*block_nx^d-offset
           rnode(rpxmin^d_,igrid)=xprobmin^d+xstretch^d+sizeuniformpart^d+dxfirst_1mq(level,^d) &
                                   *(1.0d0-qstretch(level,^d)**imin)
           rnode(rpxmax^d_,igrid)=xprobmin^d+xstretch^d+sizeuniformpart^d+dxfirst_1mq(level,^d) &
                                   *(1.0d0-qstretch(level,^d)**imax)
           ! fix possible out of bound due to precision
           if(rnode(rpxmax^d_,igrid)>xprobmax^d) rnode(rpxmax^d_,igrid)=xprobmax^d
           ixshift=(ig^d-1)*block_nx^d-nghostcells-offset
           do ix=ixgextmin^d,ixgextmax^d
             index=ixshift+ix
             ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(index-1)
           enddo
           ixshift=(ig^d+nstretchedblocks(level,^d)/2-ng^d(level)-1)*(block_nx^d/2)-nghostcells
           do ix=ixcogmin^d,ixcogmax^d
             index=ixshift+ix
             psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(index-1)
           enddo
           ! first block: modify ghost cells!!!
           if(ig^d==ng^d(level)-nstretchedblocks(level,^d)+1)then
             if(ng^d(level)==nstretchedblocks(level,^d))then
               ! if middle blocks do not exist then use symmetry
               do ix=ixgextmin^d,nghostcells
                 ps(igrid)%dx(ix^d%ixGext^s,^d)=ps(igrid)%dx(2*nghostcells+1-ix^d%ixGext^s,^d)
               enddo
               do ix=1,nghostcells
                 psc(igrid)%dx(ix^d%ixCoG^s,^d)=psc(igrid)%dx(2*nghostcells+1-ix^d%ixCoG^s,^d)
               enddo
             else
               ! if middle blocks exist then use same as middle blocks:
               do ix=ixgextmin^d,nghostcells
                 ps(igrid)%dx(ix^d%ixGext^s,^d)=dxmid(level,^d)
               enddo
               do ix=1,nghostcells
                 psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxmid(level-1,^d)
               enddo
             endif
           endif
           ! last block: make ghost cells symmetric (to allow periodicity)
           if(ig^d==ng^d(level))then
             do ix=ixghi^d-nghostcells+1,ixgextmax^d
               ps(igrid)%dx(ix^d%ixGext^s,^d)=ps(igrid)%dx(2*(ixghi^d-nghostcells)+1-ix^d%ixGext^s,^d)
             enddo
             do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
               psc(igrid)%dx(ix^d%ixCoG^s,^d)=psc(igrid)%dx(2*(ixcogmax^d-nghostcells)+1-ix^d%ixCoG^s,^d)
             enddo
           endif
         endif
       endif
       ! now that dx and grid boundaries are known: fill cell centers
       ifirst=nghostcells+1
       ! first fill the mesh
       summeddx=0.0d0
       do ix=ixmlo^d,ixmhi^d
         ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
         summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
       enddo
       ! then ghost cells to left
       summeddx=0.0d0
       do ix=nghostcells,1,-1
         ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
         summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
       enddo
       ! then ghost cells to right
       summeddx=0.0d0
       do ix=ixghi^d-nghostcells+1,ixghi^d
         ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
         summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
       enddo
       ! and next for the coarse representation
       ! first fill the mesh
       summeddx=0.0d0
       do ix=nghostcells+1,ixcogmax^d-nghostcells
         psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*psc(igrid)%dx(ix^d%ixCoG^s,^d)
         summeddx=summeddx+psc(igrid)%dx(ix^d%ifirst,^d)
       enddo
       ! then ghost cells to left
       summeddx=0.0d0
       do ix=nghostcells,1,-1
         psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*psc(igrid)%dx(ix^d%ixCoG^s,^d)
         summeddx=summeddx+psc(igrid)%dx(ix^d%ifirst,^d)
       enddo
       ! then ghost cells to right
       summeddx=0.0d0
       do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
         psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*psc(igrid)%dx(ix^d%ixCoG^s,^d)
         summeddx=summeddx+psc(igrid)%dx(ix^d%ifirst,^d)
       enddo
       select case(icase)
         case(0)
           ! if even number of ghost cells: xext is just copy of local x
           xext(ixgext^s,^d)=ps(igrid)%x(ixgext^s,^d)
         case(1)
           ! if uneven number of ghost cells: extra layer left/right
           summeddx=0.0d0
           do ix=ixmlo^d,ixmhi^d
             xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
            summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
           enddo
           ! then ghost cells to left
           summeddx=0.0d0
           do ix=nghostcells,ixgextmin^d,-1
             xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
             summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
           enddo
          ! then ghost cells to right
          summeddx=0.0d0
          do ix=ixghi^d-nghostcells+1,ixgextmax^d
             xext(ix^d%ixGext^s,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
             summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
          enddo
         case default
           call mpistop("no such case")
       end select
     endif\}
  endif
 
  ! calculate area of cell surfaces for standard block
  call get_surface_area(ps(igrid),ixg^ll)
  ! calculate area of cell surfaces for coarser representative block
  call get_surface_area(psc(igrid),ixcog^l)
  ! calculate volume and distance of cells
  ps(igrid)%dsC=1.d0
  select case (coordinate)
    case (cartesian)
      ps(igrid)%dvolume(ixgext^s)= {^d&rnode(rpdx^d_,igrid)|*}
      ps(igrid)%ds(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
      ps(igrid)%dsC(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
      psc(igrid)%dvolume(ixcog^s)= {^d&2.d0*rnode(rpdx^d_,igrid)|*}
      psc(igrid)%ds(ixcog^s,1:ndim)=psc(igrid)%dx(ixcog^s,1:ndim)
    case (cartesian_stretched)
      ps(igrid)%dvolume(ixgext^s)= {^d&ps(igrid)%dx(ixgext^s,^d)|*}
      ps(igrid)%ds(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
      ps(igrid)%dsC(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
      psc(igrid)%dvolume(ixcog^s)= {^d&psc(igrid)%dx(ixcog^s,^d)|*}
      psc(igrid)%ds(ixcog^s,1:ndim)=psc(igrid)%dx(ixcog^s,1:ndim)
    case (cartesian_expansion)
      {^ifoned
      delx_ext(ixgext^s) = ps(igrid)%dx(ixgext^s,1)
      if(associated(usr_set_surface)) call usr_set_surface(ixgext^l,xext(ixgext^s,1),delx_ext(ixgext^s),exp_factor_ext(ixgext^s),del_exp_factor_ext(ixgext^s),exp_factor_primitive_ext(ixgext^s))
      ps(igrid)%dvolume(ixgext^s)= exp_factor_primitive_ext(ixgext^s)
      ps(igrid)%ds(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
      ps(igrid)%dsC(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
      xc(ixcog^s) = psc(igrid)%x(ixcog^s,1)
      delxc(ixcog^s) = psc(igrid)%dx(ixcog^s,1)
      if(associated(usr_set_surface)) call usr_set_surface(ixcog^l,xc(ixcog^s),delxc(ixcog^s),exp_factor_coarse(ixcog^s),del_exp_factor_coarse(ixcog^s),exp_factor_primitive_coarse(ixcog^s))
      psc(igrid)%dvolume(ixcog^s)= exp_factor_primitive_coarse(ixcog^s)
      psc(igrid)%ds(ixcog^s,1)=psc(igrid)%dx(ixcog^s,1)
      }
    case (spherical)
      ps(igrid)%dvolume(ixgext^s)=(xext(ixgext^s,1)**2 &
                                +ps(igrid)%dx(ixgext^s,1)**2/12.0d0)*&
              ps(igrid)%dx(ixgext^s,1){^nooned &
             *two*dabs(dsin(xext(ixgext^s,2))) &
             *dsin(half*ps(igrid)%dx(ixgext^s,2))}{^ifthreed*ps(igrid)%dx(ixgext^s,3)}
      psc(igrid)%dvolume(ixcog^s)=(psc(igrid)%x(ixcog^s,1)**2 &
                                       +psc(igrid)%dx(ixcog^s,1)**2/12.0d0)*&
              psc(igrid)%dx(ixcog^s,1){^nooned &
             *two*dabs(dsin(psc(igrid)%x(ixcog^s,2))) &
             *dsin(half*psc(igrid)%dx(ixcog^s,2))}{^ifthreed*psc(igrid)%dx(ixcog^s,3)}
      ps(igrid)%ds(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
      {^nooned   ps(igrid)%ds(ixgext^s,2)=xext(ixgext^s,1)*ps(igrid)%dx(ixgext^s,2)}
      {^ifthreed ps(igrid)%ds(ixgext^s,3)=xext(ixgext^s,1)*dsin(xext(ixgext^s,2))*&
                                          ps(igrid)%dx(ixgext^s,3)}
      ps(igrid)%dsC(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
      {^nooned   ps(igrid)%dsC(ixgext^s,2)=(xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1))*&
                                          ps(igrid)%dx(ixgext^s,2)
      if(ndir>ndim) then
        ps(igrid)%dsC(ixgext^s,3)=(xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1))*&
                                       dsin(xext(ixgext^s,2)+half*ps(igrid)%dx(ixgext^s,2))
      end if
      }
      {^ifthreed ps(igrid)%dsC(ixgext^s,3)=(xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1))*&
                                       dsin(xext(ixgext^s,2)+half*ps(igrid)%dx(ixgext^s,2))*&
                                       ps(igrid)%dx(ixgext^s,3)}
    case (cylindrical)
      ps(igrid)%dvolume(ixgext^s)=dabs(xext(ixgext^s,1)) &
           *ps(igrid)%dx(ixgext^s,1){^de&*ps(igrid)%dx(ixgext^s,^de) }
      psc(igrid)%dvolume(ixcog^s)=dabs(psc(igrid)%x(ixcog^s,1)) &
           *psc(igrid)%dx(ixcog^s,1){^de&*psc(igrid)%dx(ixcog^s,^de) }
      ps(igrid)%ds(ixgext^s,r_)=ps(igrid)%dx(ixgext^s,r_)
      ps(igrid)%dsC(ixgext^s,r_)=ps(igrid)%dx(ixgext^s,r_)
      if(z_>0.and.z_<=ndim) then
        ps(igrid)%ds(ixgext^s,z_)=ps(igrid)%dx(ixgext^s,z_)
        ps(igrid)%dsC(ixgext^s,z_)=ps(igrid)%dx(ixgext^s,z_)
        if(phi_>z_.and.ndir>ndim) then
          ps(igrid)%dsC(ixgext^s,phi_)=xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1)
        end if
      end if
      if(phi_>0.and.phi_<=ndim) then
        ps(igrid)%ds(ixgext^s,phi_)=xext(ixgext^s,1)*ps(igrid)%dx(ixgext^s,phi_)
        ps(igrid)%dsC(ixgext^s,phi_)=(xext(ixgext^s,1)+&
                   half*ps(igrid)%dx(ixgext^s,1))*ps(igrid)%dx(ixgext^s,phi_)
        if(z_>phi_.and.ndir>ndim) ps(igrid)%dsC(ixgext^s,z_)=1.d0
      end if
    case default
      call mpistop("Sorry, coordinate unknown")
  end select
 
  ! initialize background non-evolving solution
  if (b0field) call set_b0_grid(igrid)
  if (number_equi_vars>0) call phys_set_equi_vars(igrid)
 
  ! find the blocks on the boundaries
  ps(igrid)%is_physical_boundary=.false.
  {
  do i^d=-1,1
    if(i^d==0) cycle
    ign^d=ig^d+i^d
    ! blocks at periodic boundary have neighbors in the physical domain
    ! thus threated at internal blocks with no physical boundary
    if (periodb(^d)) ign^d=1+modulo(ign^d-1,ng^d(level))
    if (ign^d > ng^d(level)) then
       if(phi_ > 0 .and. poleb(2,^d)) then
         ! if at a pole, the boundary is not physical boundary
         ps(igrid)%is_physical_boundary(2*^d)=.false.
       else
         ps(igrid)%is_physical_boundary(2*^d)=.true.
       end if
    else if (ign^d < 1) then
       if(phi_ > 0 .and. poleb(1,^d)) then
         ! if at a pole, the boundary is not physical boundary
         ps(igrid)%is_physical_boundary(2*^d-1)=.false.
       else
         ps(igrid)%is_physical_boundary(2*^d-1)=.true.
       end if
    end if
  end do
  \}
  if(any(ps(igrid)%is_physical_boundary)) then
    phyboundblock(igrid)=.true.
  else
    phyboundblock(igrid)=.false.
  end if
 
end subroutine alloc_node
 
!> allocate memory to physical state of igrid node
subroutine alloc_state(igrid, s, ixG^L, ixGext^L, alloc_once_for_ps)
  use mod_global_parameters
  type(state) :: s
  integer, intent(in) :: igrid, ixG^L, ixGext^L
  logical, intent(in) :: alloc_once_for_ps
  integer             :: ixGs^L
 
  allocate(s%w(ixg^s,1:nw))
  s%igrid=igrid
  s%w=0.d0
  s%ixG^l=ixg^l;
  {^d& ixgsmin^d = ixgmin^d-1; ixgsmax^d = ixgmax^d|;}
  if(stagger_grid) then
    allocate(s%ws(ixgs^s,1:nws))
    s%ws=0.d0
    if(record_electric_field) then
      allocate(s%we(ixgs^s,1:nws))
      s%we=0.d0
    end if
    s%ixGs^l=ixgs^l;
  end if
  if(alloc_once_for_ps) then
    ! allocate extra variables for ps state
    if(nw_extra>0) allocate(s%wextra(ixg^s,1:nw_extra))
    ! allocate coordinates
    allocate(s%x(ixg^s,1:ndim))
    allocate(s%dx(ixgext^s,1:ndim), &
               s%ds(ixgext^s,1:ndim),s%dsC(ixgext^s,1:3))
    allocate(s%dvolume(ixgext^s))
    allocate(s%surfaceC(ixgs^s,1:ndim), &
             s%surface(ixg^s,1:ndim))
    ! allocate physical boundary flag
    allocate(s%is_physical_boundary(2*ndim))
    if(b0field) then
      allocate(s%B0(ixg^s,1:ndir,0:ndim))
      allocate(s%J0(ixg^s,7-2*ndir:3))
    end if
    if(number_equi_vars > 0) then
      allocate(s%equi_vars(ixg^s,1:number_equi_vars,0:ndim))
    endif
 
    ! allocate space for special values for each block state
    if(phys_trac) then
      ! special_values(1) Tcoff local
      ! special_values(2) Tmax local
      ! special_values(3:2+ndim) Bdir local
      allocate(s%special_values(ndim+2))
    end if
  else
    ! share common info from ps states to save memory
    if(nw_extra>0) s%wextra=>ps(igrid)%wextra
    s%x=>ps(igrid)%x
    s%dx=>ps(igrid)%dx
    s%ds=>ps(igrid)%ds
    s%dsC=>ps(igrid)%dsC
    s%dvolume=>ps(igrid)%dvolume
    s%surfaceC=>ps(igrid)%surfaceC
    s%surface=>ps(igrid)%surface
    s%is_physical_boundary=>ps(igrid)%is_physical_boundary
    if(b0field) then
      s%B0=>ps(igrid)%B0
      s%J0=>ps(igrid)%J0
    end if
    if(number_equi_vars > 0) then
      s%equi_vars=>ps(igrid)%equi_vars
    endif
    if(phys_trac) s%special_values=>ps(igrid)%special_values
  end if
end subroutine alloc_state
 
!> allocate memory to one-level coarser physical state of igrid node
subroutine alloc_state_coarse(igrid, s, ixG^L, ixGext^L)
  use mod_global_parameters 
  use mod_mhd_phys, only: mhd_semirelativistic
  type(state) :: s
  integer, intent(in) :: igrid, ixG^L, ixGext^L
  integer             :: ixGs^L
 
  allocate(s%w(ixg^s,1:nw))
  s%igrid=igrid
  s%w=0.d0
  s%ixG^l=ixg^l;
  {^d& ixgsmin^d = ixgmin^d-1; ixgsmax^d = ixgmax^d|;}
  if(stagger_grid) then
    allocate(s%ws(ixgs^s,1:nws))
    s%ws=0.d0
    s%ixGs^l=ixgs^l;
  end if
  if(b0field.and.mhd_semirelativistic) then
    allocate(s%B0(ixg^s,1:ndir,0:ndim))
  end if
  ! allocate coordinates
  allocate(s%x(ixg^s,1:ndim))
  allocate(s%dx(ixgext^s,1:ndim), &
             s%ds(ixgext^s,1:ndim),s%dsC(ixgext^s,1:3))
  allocate(s%dvolume(ixgext^s))
  allocate(s%surfaceC(ixgs^s,1:ndim), &
           s%surface(ixg^s,1:ndim))
  ! allocate physical boundary flag
  allocate(s%is_physical_boundary(2*ndim))
end subroutine alloc_state_coarse
 
subroutine dealloc_state(igrid, s,dealloc_x)
  use mod_global_parameters
  integer, intent(in) :: igrid
  type(state) :: s
  logical, intent(in) :: dealloc_x
 
  deallocate(s%w)
  if(stagger_grid) then
    deallocate(s%ws)
  end if
  if(dealloc_x) then
    if(nw_extra>0) deallocate(s%wextra)
    ! deallocate coordinates
    deallocate(s%x)
    deallocate(s%dx,s%ds,s%dsC)
    deallocate(s%dvolume)
    deallocate(s%surfaceC,s%surface)
    deallocate(s%is_physical_boundary)
    if(b0field) then
      deallocate(s%B0)
      deallocate(s%J0)
    end if
    if(number_equi_vars > 0) then
      deallocate(s%equi_vars)
    end if
  else
    nullify(s%x,s%dx,s%ds,s%dsC,s%dvolume,s%surfaceC,s%surface)
    nullify(s%is_physical_boundary)
    if(b0field) nullify(s%B0,s%J0)
    if(number_equi_vars > 0) then
      nullify(s%equi_vars)
    end if
    if(nw_extra>0) nullify(s%wextra)
  end if
end subroutine dealloc_state
 
subroutine dealloc_state_coarse(igrid, s)
  use mod_global_parameters
  use mod_mhd_phys, only: mhd_semirelativistic
  integer, intent(in) :: igrid
  type(state) :: s
 
  deallocate(s%w)
  if(stagger_grid) then
    deallocate(s%ws)
  end if
  if(b0field.and.mhd_semirelativistic) then
    deallocate(s%B0)
  end if
  ! deallocate coordinates
  deallocate(s%x)
  deallocate(s%dx,s%ds,s%dsC)
  deallocate(s%dvolume)
  deallocate(s%surfaceC,s%surface)
  deallocate(s%is_physical_boundary)
end subroutine dealloc_state_coarse
 
subroutine dealloc_node(igrid)
  use mod_global_parameters
 
  integer, intent(in) :: igrid
 
  if (igrid==0) then
     call mpistop("trying to delete a non-existing grid in dealloc_node")
  end if
 
  call dealloc_state(igrid, ps(igrid),.true.)
  call dealloc_state_coarse(igrid, psc(igrid))
  if(.not.convert) then
    call dealloc_state(igrid, ps1(igrid),.false.)
    call dealloc_state(igrid, pso(igrid),.false.)
    ! deallocate temporary solution space
    select case (t_integrator)
    case(ssprk3,ssprk4,imex_midpoint,imex_trapezoidal,imex_222)
      call dealloc_state(igrid, ps2(igrid),.false.)
    case(rk3_bt,rk4,ssprk5,imex_cb3a)
      call dealloc_state(igrid, ps2(igrid),.false.)
      call dealloc_state(igrid, ps3(igrid),.false.)
    case(imex_ars3,imex_232)
      call dealloc_state(igrid, ps2(igrid),.false.)
      call dealloc_state(igrid, ps3(igrid),.false.)
      call dealloc_state(igrid, ps4(igrid),.false.)
    end select
  end if
 
end subroutine dealloc_node