mod_amr_solution_node.t

Go to the documentation of this file.

module mod_amr_solution_node
  use mod_comm_lib, only: mpistop
 
  implicit none
  private
 
  public :: getnode, putnode
  public :: alloc_node, alloc_state 
  public :: dealloc_node 
 
contains
 
 
  !> 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
    use mod_b0, only: set_b0_grid 
 
    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 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
      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(local_timestep) then
        allocate(s%dt(ixg^s))
      endif
  
      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
    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(b0fieldalloccoarse) 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)
      if(record_electric_field) deallocate(s%we)
    end if
    if(dealloc_x) then
      if(nw_extra>0) deallocate(s%wextra)
      ! deallocate coordinates
      deallocate(s%x)
      deallocate(s%dx,s%dt,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
      if(phys_trac) then
        deallocate(s%special_values)
      end if
    end if
    nullify(s%x,s%dx,s%dt,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)
    if(phys_trac) nullify(s%special_values)
  end subroutine dealloc_state
  
  subroutine dealloc_state_coarse(igrid, s)
    use mod_global_parameters
    integer, intent(in) :: igrid
    type(state) :: s
  
    deallocate(s%w)
    if(stagger_grid) then
      deallocate(s%ws)
    end if
    if(b0fieldalloccoarse) then
      deallocate(s%B0)
    end if
    ! deallocate coordinates
    deallocate(s%x)
    deallocate(s%dx,s%dt,s%ds,s%dsC)
    deallocate(s%dvolume)
    deallocate(s%surfaceC,s%surface)
    nullify(s%x,s%dx,s%dt,s%ds,s%dsC,s%dvolume,s%surfaceC,s%surface)
    nullify(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.)
      ! 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
 
end module mod_amr_solution_node