mod_slice.t

Go to the documentation of this file.

!> Writes D-1 slice, can do so in various formats, depending on slice_type
module mod_slice
  use mod_basic_types
  use mod_comm_lib, only: mpistop
  implicit none
 
  !> Slice coordinates, see @ref slices.md
  double precision :: slicecoord(1000)
 
  !> Maximum number of slices
  integer, parameter :: nslicemax=1000
 
  !> Number of slices to output
  integer :: nslices
 
  !> The slice direction for each slice
  integer :: slicedir(nslicemax)
 
  !> tag for MPI message
  integer, private :: itag
 
  !> choose data type of slice: vtu, vtuCC, dat, or csv
  character(len=std_len) :: slice_type
  
  !> igrid and ipe of the leaf blocks on the slice, with the first (1,1:2) stores the nleafs and nparents
  integer, allocatable :: sfc_sub(:,:)
 
contains
 
  subroutine write_slice
    use mod_global_parameters
    ! Writes a D-1 slice 
    ! by Oliver Porth
    ! 22.Nov 2011
    integer :: islice
 
    do islice=1,nslices
       call put_slice(slicedir(islice),slicecoord(islice))
    end do
 
    slicenext=slicenext+1
  end subroutine write_slice
 
  subroutine put_slice(dir,xslice)
    use mod_forest, only: morton_sub_start, morton_sub_stop
    use mod_global_parameters
    ! Writes a D-1 slice 
    ! For ONED simulations, output will be appended to one csv-file per slice
    ! slices are sensitive to the saveprim switch and 
    ! can contain auxiliary io variables (nwauxio)
    ! Thus csv or vtu(CC)-files with primitive variables are obtained.  
    ! when slice_type='dat', we save a D-1 dat file for potential restarts
    ! by Oliver Porth
    ! 22.Nov 2011
    integer, intent(in) :: dir
    double precision, intent(in) :: xslice
    ! .. local ..
    double precision,dimension(0:nw+nwauxio) :: normconv 
    integer :: Njgrid, jgrid
    integer, dimension(ndim-1) :: ixsubGlo, ixsubGhi
    integer, dimension(ndim-1) :: ixsubMlo, ixsubMhi
    integer :: size_subblock_io, nx^D, slice_fh, nwexpand
    integer :: type_subblock_io, type_subblockC_io, type_subblock_x_io, type_subblockC_x_io
    integer, dimension(ndim) :: sizes, subsizes, start
 
    double precision, dimension(ndim-1) :: xprobminsub, xprobmaxsub
    integer, dimension(ndim-1) :: block_nx_sub, domain_nx_sub
    logical, dimension(ndim-1) :: periodBsub
  
    ! Preamble: 
    nx^d=ixmhi^d-ixmlo^d+1;
    slice_fh=unitslice
 
    if (ndim==1) then
       nwexpand = nwauxio
    else
       if(slice_type/='dat')then
          nwexpand = nwauxio
       else
          nwexpand = 0
       endif
    end if
 
    ! Do a last consistency check:
    select case(dir)
       {case(^d)
       if(xslice<xprobmin^d.or.xslice>xprobmax^d) &
            call mpistop("slice out of bounds")
       \}
    end select
 
    ! Traverse the forest and fill nodes:
    call select_slice(dir,xslice,.false.,slice_fh,normconv)
 
    ! Create the MPI-datatype and select indices:
    {^ifthreed
    select case(dir)
    case (1)
       sizes(1) = ixghi2; sizes(2) = ixghi3;
       ixsubglo(1) = ixglo2; ixsubglo(2) = ixglo3;
       ixsubghi(1) = ixghi2; ixsubghi(2) = ixghi3;
       subsizes(1)=nx2;subsizes(2)=nx3;
       start(1)=ixmlo2-1;start(2)=ixmlo3-1;
       size_subblock_io=nx2*nx3*(nw+nwexpand)*size_double
       xprobminsub(1)=xprobmin2; xprobmaxsub(1)=xprobmax2;
       xprobminsub(2)=xprobmin3; xprobmaxsub(2)=xprobmax3;
       domain_nx_sub(1)=domain_nx2; domain_nx_sub(2)=domain_nx3;
       block_nx_sub(1)=block_nx2; block_nx_sub(2)=block_nx3;
       periodbsub(1)=periodb(2); periodbsub(2)=periodb(3);
    case (2)
       sizes(1) = ixghi1; sizes(2) = ixghi3;
       ixsubglo(1) = ixglo1; ixsubglo(2) = ixglo3;
       ixsubghi(1) = ixghi1; ixsubghi(2) = ixghi3;
       subsizes(1)=nx1;subsizes(2)=nx3;
       start(1)=ixmlo1-1;start(2)=ixmlo3-1;
       size_subblock_io=nx1*nx3*(nw+nwexpand)*size_double
       xprobminsub(1)=xprobmin1; xprobmaxsub(1)=xprobmax1;
       xprobminsub(2)=xprobmin3; xprobmaxsub(2)=xprobmax3;
       domain_nx_sub(1)=domain_nx1; domain_nx_sub(2)=domain_nx3;
       block_nx_sub(1)=block_nx1; block_nx_sub(2)=block_nx3;
       periodbsub(1)=periodb(1); periodbsub(2)=periodb(3);
    case (3)
       ixsubglo(1) = ixglo1; ixsubglo(2) = ixglo2;
       ixsubghi(1) = ixghi1; ixsubghi(2) = ixghi2;
       sizes(1) = ixghi1; sizes(2) = ixghi2;
       subsizes(1)=nx1;subsizes(2)=nx2;
       start(1)=ixmlo1-1;start(2)=ixmlo2-1;
       size_subblock_io=nx1*nx2*(nw+nwexpand)*size_double
       xprobminsub(1)=xprobmin1; xprobmaxsub(1)=xprobmax1;
       xprobminsub(2)=xprobmin2; xprobmaxsub(2)=xprobmax2;
       domain_nx_sub(1)=domain_nx1; domain_nx_sub(2)=domain_nx2;
       block_nx_sub(1)=block_nx1; block_nx_sub(2)=block_nx2;
       periodbsub(1)=periodb(1); periodbsub(2)=periodb(2);
    case default
       call mpistop("slice direction not clear in put_slice")
    end select
    }
    {^iftwod
    select case(dir)
    case (1)
       ixsubglo(1) = ixglo2; ixsubghi(1) = ixghi2;
       sizes(1) = ixghi2
       subsizes(1)=nx2
       start(1)=ixmlo2-1
       size_subblock_io=nx2*(nw+nwexpand)*size_double
       xprobminsub(1)=xprobmin2; xprobmaxsub(1)=xprobmax2;
       domain_nx_sub(1)=domain_nx2;
       block_nx_sub(1)=block_nx2;
       periodbsub(1)=periodb(2);
    case (2)
       ixsubglo(1) = ixglo1; ixsubghi(1) = ixghi1;
       sizes(1) = ixghi1
       subsizes(1)=nx1
       start(1)=ixmlo1-1
       size_subblock_io=nx1*(nw+nwexpand)*size_double
       xprobminsub(1)=xprobmin1; xprobmaxsub(1)=xprobmax1;
       domain_nx_sub(1)=domain_nx1;
       block_nx_sub(1)=block_nx1;
       periodbsub(1)=periodb(1);
    case default
       call mpistop("slice direction not clear in put_slice")
    end select
    }
    {^ifoned
    size_subblock_io=(nw+nwexpand)*size_double
    }
 
    {^nooned
    {^de&ixsubmlo(^de-1) = ixsubglo(^de-1)+nghostcells;}
    {^de&ixsubmhi(^de-1) = ixsubghi(^de-1)-nghostcells;}
    }
 
    sizes(ndim)=nw+nwexpand
    subsizes(ndim)=nw+nwexpand
    start(ndim)=0
 
    ! Types for center variables:
    call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
         mpi_order_fortran,mpi_double_precision, &
         type_subblock_io,ierrmpi)
    call mpi_type_commit(type_subblock_io,ierrmpi)
 
    sizes(ndim)=^nd
    subsizes(ndim)=^nd
    start(ndim)=0
    call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
         mpi_order_fortran,mpi_double_precision, &
         type_subblock_x_io,ierrmpi)
    call mpi_type_commit(type_subblock_x_io,ierrmpi)
 
 
    ! Types for corner variables:
    subsizes(1:ndim-1) = subsizes(1:ndim-1) + 1
    start(1:ndim-1)    = start(1:ndim-1) - 1 
    sizes(ndim)=nw+nwexpand
    subsizes(ndim)=nw+nwexpand
    start(ndim)=0
    
    call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
         mpi_order_fortran,mpi_double_precision, &
         type_subblockc_io,ierrmpi)
    call mpi_type_commit(type_subblockc_io,ierrmpi)
 
    sizes(ndim)=^nd
    subsizes(ndim)=^nd
    start(ndim)=0
    call mpi_type_create_subarray(ndim,sizes,subsizes,start, &
         mpi_order_fortran,mpi_double_precision, &
         type_subblockc_x_io,ierrmpi)
    call mpi_type_commit(type_subblockc_x_io,ierrmpi)
 
    
    ! local number of sub-grids:
    njgrid=morton_sub_stop(mype)-morton_sub_start(mype)+1
 
    ! Now output using various schemes: 
    if (ndim==1) then 
       call put_slice_zerod
    else
       select case(slice_type)
       case ('csv')
          call put_slice_csv
       case ('dat')
          call put_slice_dat
       case ('vtu', 'vtuCC')
          call put_slice_vtu
       end select
    end if
 
    do jgrid=1,njgrid
       call dealloc_subnode(jgrid)
    end do
 
    call mpi_type_free(type_subblock_io,ierrmpi)
    call mpi_type_free(type_subblock_x_io,ierrmpi)
    call mpi_type_free(type_subblockc_io,ierrmpi)
    call mpi_type_free(type_subblockc_x_io,ierrmpi)
 
    if (allocated(sfc_sub)) deallocate(sfc_sub)
 
  contains
 
    subroutine put_slice_vtu
 
      use mod_calculate_xw
      integer :: status(MPI_STATUS_SIZE), ipe
      logical             :: fileopen
      character(len=1024) :: filename
      character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
      character(len=1024) :: outfilehead
 
      if (mype==0) then
 
         inquire(slice_fh,opened=fileopen)
         if(.not.fileopen)then
            ! generate filename: 
            write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//slice_label(xslice)//'_n',slicenext,'.vtu'
            open(slice_fh,file=filename,status='unknown',form='formatted')
         end if
         ! get and write the header: 
         call getheadernames(wnamei,xandwnamei,outfilehead)
         ! generate xml header
         write(slice_fh,'(a)')'<?xml version="1.0"?>'
         write(slice_fh,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
         if(type_endian==1)then
            write(slice_fh,'(a)')' version="0.1" byte_order="LittleEndian">'
         else
            write(slice_fh,'(a)')' version="0.1" byte_order="BigEndian">'
         endif
         write(slice_fh,'(a)')'  <UnstructuredGrid>'
         write(slice_fh,'(a)')'<FieldData>'
         write(slice_fh,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
              'NumberOfTuples="1" format="ascii">'
         write(slice_fh,*) real(global_time*time_convert_factor)
         write(slice_fh,'(a)')'</DataArray>'
         write(slice_fh,'(a)')'</FieldData>'
 
         ! write to file:
         do jgrid=1, njgrid
            call write_slice_vtk(jgrid,slice_fh,wnamei)
         end do
 
         ! create a recv buffer using allocate, will be deallocated at the end of the routine:
         call alloc_subnode(njgrid+1,dir,nwauxio)
 
      end if
 
      ! Also communicate the normconv array since processor zero might not have it yet:
      if (npe>1) then
         do ipe=1,npe-1
            do jgrid=1,morton_sub_stop(ipe)-morton_sub_start(ipe)+1
               itag=morton_sub_start(ipe)+jgrid-1
               itag=itag*5
               if (ipe == mype ) then 
                  call mpi_send(ps_sub(jgrid)%x,1,type_subblock_x_io,0,itag,icomm,ierrmpi)
                  call mpi_send(ps_sub(jgrid)%w,1,type_subblock_io,0,itag+1,icomm,ierrmpi)
                  call mpi_send(ps_sub(jgrid)%xC,1,type_subblockc_x_io,0,itag+2,icomm,ierrmpi)
                  call mpi_send(ps_sub(jgrid)%wC,1,type_subblockc_io,0,itag+3,icomm,ierrmpi)
                  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag+4,icomm,ierrmpi)
               end if
               if (mype == 0) then
                  call mpi_recv(ps_sub(njgrid+1)%x,1,type_subblock_x_io,ipe,itag,icomm,status,ierrmpi)
                  call mpi_recv(ps_sub(njgrid+1)%w,1,type_subblock_io,ipe,itag+1,icomm,status,ierrmpi)
                  call mpi_recv(ps_sub(njgrid+1)%xC,1,type_subblockc_x_io,ipe,itag+2,icomm,status,ierrmpi)
                  call mpi_recv(ps_sub(njgrid+1)%wC,1,type_subblockc_io,ipe,itag+3,icomm,status,ierrmpi)
                  call mpi_recv(normconv,nw+nwauxio+1,mpi_double_precision,ipe,&
                       itag+4,icomm,status,ierrmpi)
                  call write_slice_vtk(njgrid+1,slice_fh,wnamei)
               end if
            end do
         end do
      endif
 
      if (mype==0) then
 
         write(slice_fh,'(a)')'</UnstructuredGrid>'
         write(slice_fh,'(a)')'</VTKFile>'
         close(slice_fh)
         call dealloc_subnode(njgrid+1)
 
      end if
 
    end subroutine put_slice_vtu
 
    subroutine write_slice_vtk(jgrid,slice_fh,wnamei)
 
      ! this only works for 2D and 3D, 1D reduction (line to point) not allowed
      integer, intent(in)           :: jgrid, slice_fh
      character(len=name_len), intent(in) :: wnamei(1:nw+nwauxio)
      ! This remainder part only for more than 1D, but nesting with NOONED gives problems 
      {#IFNDEF D1
      ! .. local ..
      double precision              :: x_VTK(1:3)
      double precision, parameter   :: minvalue = 1.0d-99, maxvalue = 1.0d+99
      integer                       :: ixC^L, ixCC^L, nc, np, iw
      integer                       :: nx^DM, nxC^DM, icell, ix^DM
      integer                       :: VTK_type
 
      {^dm&ixccmin^dm = ixsubmlo(^dm);}
      {^dm&ixccmax^dm = ixsubmhi(^dm);}
      {^dm&ixcmin^dm  = ixsubmlo(^dm)-1;}
      {^dm&ixcmax^dm  = ixsubmhi(^dm);}
 
      nx^dm=ixccmax^dm-ixccmin^dm+1;
      nxc^dm=nx^dm+1;
      nc={nx^dm*}      ! Number of cells per subgrid
      np={nxc^dm*}     ! Number of corner points per subgrid
 
      ! we write out every grid as one VTK PIECE
      write(slice_fh,'(a,i7,a,i7,a)') &
           '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
 
      !==============================
      ! celldata or pointdata?
      !==============================
      select case(slice_type)
 
      case('vtuCC') ! celldata
         write(slice_fh,'(a)')'<CellData>'
         do iw=1,nw+nwauxio
            if(iw<=nw) then 
               if(.not.w_write(iw)) cycle
            endif
            write(slice_fh,'(a,a,a)')&
                 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
            write(slice_fh,'(200(1pe14.6))') {^dm&(|}roundoff_minmax(ps_sub(jgrid)%w(ix^dm,iw)*normconv(iw),minvalue,maxvalue),{ix^dm=ixccmin^dm,ixccmax^dm)}
            write(slice_fh,'(a)')'</DataArray>'
         enddo
         write(slice_fh,'(a)')'</CellData>'
 
 
      case('vtu') ! pointdata
         write(slice_fh,'(a)')'<PointData>'
         do iw=1,nw+nwauxio
            if(iw<=nw) then 
               if(.not.w_write(iw)) cycle
            endif
            write(slice_fh,'(a,a,a)')&
                 '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
            write(slice_fh,'(200(1pe14.6))') {^dm&(|}roundoff_minmax(ps_sub(jgrid)%wC(ix^dm,iw)*normconv(iw),minvalue,maxvalue),{ix^dm=ixcmin^dm,ixcmax^dm)}
            write(slice_fh,'(a)')'</DataArray>'
         enddo
         write(slice_fh,'(a)')'</PointData>'
         
      
      end select
      !==============================
      ! Done: celldata or pointdata?
      !==============================
 
      !==============================
      ! output Cornerpoints
      !==============================
      write(slice_fh,'(a)')'<Points>'
      write(slice_fh,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
      ! write cell corner coordinates in a backward dimensional loop, always 3D output
      {do ix^dmb=ixcmin^dmb,ixcmax^dmb \}
            x_vtk(1:3)=zero;
            x_vtk(1:ndim)=ps_sub(jgrid)%xC(ix^dm,1:ndim)*normconv(0);
            write(slice_fh,'(3(1pe14.6))') x_vtk
      {^dm&end do \}
      write(slice_fh,'(a)')'</DataArray>'
      write(slice_fh,'(a)')'</Points>'
      !==============================
      ! Done: output Cornerpoints
      !==============================
 
      !==============================
      ! cell Metainformation
      !==============================
      write(slice_fh,'(a)')'<Cells>'
 
      ! connectivity part
      write(slice_fh,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
 
      {^dm& do ix^dmb=1,nx^dmb\}
      {^iftwod
      write(slice_fh,'(2(i7,1x))')ix1-1,ix1
      }{^ifthreed
      write(slice_fh,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
           (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
      }{^dm& end do\}
 
      write(slice_fh,'(a)')'</DataArray>'
 
      ! offsets data array
      write(slice_fh,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
      do icell=1,nc
         write(slice_fh,'(i7)') icell*(2**(^nd-1))
      end do
      write(slice_fh,'(a)')'</DataArray>'
 
      ! VTK cell type data array
      write(slice_fh,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
      ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
      {^iftwod vtk_type=3 \}
      {^ifthreed vtk_type=8 \}
      do icell=1,nc
         write(slice_fh,'(i2)') vtk_type
      enddo
      write(slice_fh,'(a)')'</DataArray>'
      
      write(slice_fh,'(a)')'</Cells>'
      !==============================
      ! Done: cell Metainformation
      !==============================
      write(slice_fh,'(a)')'</Piece>'
 
      }
    end subroutine write_slice_vtk
 
    subroutine put_slice_csv
 
      use mod_calculate_xw
      integer                       :: iw, ipe, itag
      integer                       :: status(MPI_STATUS_SIZE)
      logical                       :: fileopen
      character(len=1024)           :: filename
      character(len=name_len)       :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
      character(len=1024)           :: outfilehead
      character(len=1024)           :: line
 
      if (mype==0) then
         inquire(slice_fh,opened=fileopen)
         if(.not.fileopen)then
            ! generate filename: 
            write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//slice_label(xslice)//'_n',slicenext,'.csv'
            open(slice_fh,file=filename,status='unknown',form='formatted')
         end if
         ! get and write the header: 
         call getheadernames(wnamei,xandwnamei,outfilehead)
         line=''
         do iw=1,ndim+nw+nwauxio-1
            line = trim(line)//trim(xandwnamei(iw))//', '
         end do
         line = trim(line)//trim(xandwnamei(ndim+nw+nwauxio))
         write(slice_fh,'(a)')trim(line)
         ! create a recv buffer using allocate, will be deallocated at the end of the routine:
         call alloc_subnode(njgrid+1,dir,nwauxio)
 
         ! write to file:
         do jgrid=1, njgrid
            call put_slice_line(jgrid,slice_fh)
         end do
      end if
 
      ! Also communicate the normconv array since processor zero might not have it yet:
      if (npe>1) then
         do ipe=1,npe-1
            do jgrid=1,morton_sub_stop(ipe)-morton_sub_start(ipe)+1
               itag=morton_sub_start(ipe)+jgrid-1
               if (ipe == mype ) then 
                  call mpi_send(ps_sub(jgrid)%x,1,type_subblock_x_io,0,itag,icomm,ierrmpi)
                  call mpi_send(ps_sub(jgrid)%w,1,type_subblock_io,0,itag,icomm,ierrmpi)
                  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag,icomm,ierrmpi)
               end if
               if (mype == 0) then
                  call mpi_recv(ps_sub(njgrid+1)%x,1,type_subblock_x_io,ipe,itag,icomm,status,ierrmpi)
                  call mpi_recv(ps_sub(njgrid+1)%w,1,type_subblock_io,ipe,itag,icomm,status,ierrmpi)
                  call mpi_recv(normconv,nw+nwauxio+1,mpi_double_precision,ipe,&
                       itag,icomm,status,ierrmpi)
                  call put_slice_line(njgrid+1,slice_fh)
               end if
            end do
         end do
      endif
 
      if (mype==0) then
         close(slice_fh)
         call dealloc_subnode(njgrid+1)
      end if
 
    end subroutine put_slice_csv
 
    subroutine put_slice_line(jout,file_handle)
      integer, intent(in) :: jout, file_handle
      ! .. local ..
      double precision, parameter :: minvalue = 1.0d-99, maxvalue = 1.0d+99
      integer :: ix^D,idir,iw
      character(len=1024) ::line, data
 
      {^ifthreed
      do ix2=ixsubmlo(2),ixsubmhi(2)
         do ix1=ixsubmlo(1),ixsubmhi(1)
            \}
            {^iftwod
            do ix1=ixsubmlo(1),ixsubmhi(1)
               }
               ! Format the line:
               line = ''
               do idir=1,ndim
                  {^ifthreed
                  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%x(ix1,ix2,idir),minvalue,maxvalue)
                  }
                  {^iftwod
                  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%x(ix1,idir),minvalue,maxvalue)
                  }
                  {^ifoned
                  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%x(idir),minvalue,maxvalue)
                  }
 
                  line = trim(line)//trim(data)//', '
               end do
               do iw = 1,nw+nwauxio-1
                  {^ifthreed
                  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,ix2,iw)*normconv(iw),minvalue,maxvalue)
                  }
                  {^iftwod
                  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,iw)*normconv(iw),minvalue,maxvalue)
                  }
                  {^ifoned
                  write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(iw)*normconv(iw),minvalue,maxvalue)
                  }
                  line = trim(line)//trim(data)//', '
               end do
               {^ifthreed
               write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,ix2,nw+nwauxio)*normconv(nw+nwauxio),minvalue,maxvalue)
               }
               {^iftwod
               write(data,"(es14.6)")roundoff_minmax(ps_sub(jout)%w(ix1,nw+nwauxio)*normconv(nw+nwauxio),minvalue,maxvalue)
               }
               line = trim(line)//trim(data)
               write(file_handle,'(a)')trim(line)
               {^iftwod
            end do
            }
            {^ifthreed
         end do
      end do
      \}
 
    end subroutine put_slice_line
 
    subroutine put_slice_dat
 
      use mod_forest, only: tree_node, igrid_to_node
      use mod_input_output_helper, only: version_number, save_now
      use mod_physics, only: physics_type, phys_write_info
 
      integer, dimension(max_blocks*2) :: iorequest
      integer, dimension(MPI_STATUS_SIZE,max_blocks) :: iostatus
      integer(kind=MPI_OFFSET_KIND) :: offset, offset1
      integer :: nsubleafs
      integer :: amode, status(MPI_STATUS_SIZE), iwrite
      character(len=1024) :: filename
      character(len=name_len) :: dname
 
      integer :: nsubparents, iw, ileaf, inode
      integer(kind=MPI_OFFSET_KIND) :: block_offset, tree_offset
      integer, allocatable :: block_ig_sub(:,:)
      integer, allocatable :: block_lvl_sub(:)
      integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset_sub(:)
      logical, allocatable :: isleaf_sub(:)
      type(tree_node), pointer :: pnode
 
      nsubleafs=morton_sub_stop(npe-1)
 
      ! only on processor 0, prepare all the data for writing header
      if (mype==0) then
 
         if (nsubleafs /= sfc_sub(0,1)) error stop "nsubleafs /= sfc_sub(0,1)"
         nsubparents = sfc_sub(0,2)
         allocate(block_ig_sub(ndim-1,nsubleafs))
         allocate(block_lvl_sub(nsubleafs))
         allocate(block_offset_sub(nsubleafs))
         allocate(isleaf_sub(nsubleafs+nsubparents))
 
         block_offset = int(0,kind=mpi_offset_kind)
         ileaf = 0
         ! get block_ig_sub, block_lvl_sub, block_offset_sub
         do inode=1,nsubleafs+nsubparents
            if (sfc_sub(inode,1) == 0) then
               isleaf_sub(inode) = .false.
            else
               isleaf_sub(inode) = .true.
               pnode => igrid_to_node(sfc_sub(inode,1),sfc_sub(inode,2))%node
               ileaf = ileaf + 1
               {^ifthreed
               select case (dir)
                  case (1)
                     block_ig_sub(:,ileaf) = [ pnode%ig2, pnode%ig3 ]
                  case (2)
                     block_ig_sub(:,ileaf) = [ pnode%ig1, pnode%ig3 ]
                  case (3)
                     block_ig_sub(:,ileaf) = [ pnode%ig1, pnode%ig2 ]
               end select
               }
               {^iftwod
               select case (dir)
                  case (1)
                     block_ig_sub(:,ileaf) = [ pnode%ig2 ]
                  case (2)
                     block_ig_sub(:,ileaf) = [ pnode%ig1 ]
               end select
               }
               block_lvl_sub(ileaf) = pnode%level
               ! will not save ghost cells, start from offset_block_data
               block_offset_sub(ileaf) = block_offset
               block_offset = block_offset + int(2*(ndim-1)*size_int,kind=mpi_offset_kind)
               block_offset = block_offset + int(product(block_nx_sub)*nw*size_double,kind=mpi_offset_kind)
            end if
         end do
      end if
 
      ! generate filename
      write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//slice_label(xslice)//'_n',slicenext,'.dat'
 
      if(mype==0) then
         open(unit=slice_fh,file=filename,status='replace')
         close(unit=slice_fh)
      end if
 
      ! header writing 
      if (mype==0) then
         ! call select_slice(dir,xslice,.true.,slice_fh,normconv)
         amode=ior(mpi_mode_append,mpi_mode_wronly)
         call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null, &
              slice_fh,ierrmpi)
 
         call mpi_file_write(slice_fh, version_number, 1, mpi_integer, status, ierrmpi)
         ! two offsets for tree and block data
         call mpi_file_write(slice_fh, 0, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, 0, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, nw, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, ndir, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, ndim-1, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, levmax_sub, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, nsubleafs, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, nsubparents, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, it, 1, mpi_integer, status, ierrmpi)
         ! Note: It is nice when this double has an even number of 4 byte
         ! integers before it (for alignment)
         call mpi_file_write(slice_fh, global_time, 1, mpi_double_precision, status, ierrmpi)
 
         call mpi_file_write(slice_fh, xprobminsub, ndim-1, &
               mpi_double_precision, status, ierrmpi)
         call mpi_file_write(slice_fh, xprobmaxsub, ndim-1, &
               mpi_double_precision, status, ierrmpi)
         call mpi_file_write(slice_fh, domain_nx_sub, ndim-1, &
               mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, block_nx_sub, ndim-1, &
               mpi_integer, status, ierrmpi)
 
         ! Periodicity (assume all variables are periodic if one is)
         call mpi_file_write(slice_fh, periodbsub, ndim-1, mpi_logical, status, ierrmpi)
 
         ! Geometry
         call mpi_file_write(slice_fh, geometry_name(1:name_len), &
               name_len, mpi_character, status, ierrmpi)
 
         ! Write stagger grid mark
         call mpi_file_write(slice_fh, stagger_grid, 1, mpi_logical, status, ierrmpi)
 
         do iw = 1, nw
            ! using directly trim(adjustl((cons_wnames(iw)))) in MPI_FILE_WRITE call 
            ! does not work, there will be trailing characters
            dname = trim(adjustl((cons_wnames(iw))))
            call mpi_file_write(slice_fh, dname, name_len, mpi_character, status, ierrmpi)
         end do
 
         ! Physics related information
         call mpi_file_write(slice_fh, physics_type, name_len, mpi_character, status, ierrmpi)
 
         ! Format:
         ! integer :: n_par
         ! double precision :: values(n_par)
         ! character(n_par * name_len) :: names
         call phys_write_info(slice_fh)
 
         ! Write snapshotnext etc., which is useful for restarting.
         ! Note we add one, since snapshotnext is updated *after* this procedure
         if(pass_wall_time.or.save_now) then
            call mpi_file_write(slice_fh, snapshotnext, 1, mpi_integer, status, ierrmpi)
         else
            call mpi_file_write(slice_fh, snapshotnext+1, 1, mpi_integer, status, ierrmpi)
         end if
         call mpi_file_write(slice_fh, slicenext, 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, collapsenext, 1, mpi_integer, status, ierrmpi)
 
         call mpi_file_get_position(slice_fh, tree_offset, ierrmpi)
         call mpi_file_write(slice_fh, isleaf_sub, nsubleafs+nsubparents, mpi_logical, status, ierrmpi)
         call mpi_file_write(slice_fh, block_lvl_sub, nsubleafs, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, block_ig_sub, size(block_ig_sub), mpi_integer, status, ierrmpi)
 
         ! add the block_offset(above) to block_offset_sub
         call mpi_file_get_position(slice_fh, block_offset, ierrmpi)
         block_offset_sub(:) = block_offset_sub(:) + block_offset + &
            int(nsubleafs*2*size_int, kind=mpi_offset_kind)
 
         call mpi_file_write(slice_fh, block_offset_sub, nsubleafs, mpi_offset, status, ierrmpi)
 
         ! get the right block_offset
         call mpi_file_get_position(slice_fh, block_offset, ierrmpi)
         if (block_offset - tree_offset /= &
            (nsubleafs + nsubparents) * size_logical + &
            nsubleafs * ((ndim-1+1) * size_int + 2 * size_int)) then
            print *, "Warning: MPI_OFFSET type /= 8 bytes"
            print *, "This *could* cause problems when reading .dat files"
         end if
 
         ! come back to write offsets
         call mpi_file_seek(slice_fh, int(size_int,kind=mpi_offset_kind), mpi_seek_set, ierrmpi)
         call mpi_file_write(slice_fh, int(tree_offset), 1, mpi_integer, status, ierrmpi)
         call mpi_file_write(slice_fh, int(block_offset), 1, mpi_integer, status, ierrmpi)
 
         call mpi_file_close(slice_fh,ierrmpi)
      end if
 
      ! Broadcast nsubleafs and nsubparents to all processors
      call mpi_bcast(nsubleafs, 1, mpi_integer, 0, icomm, ierrmpi)
      call mpi_bcast(nsubparents, 1, mpi_integer, 0, icomm, ierrmpi)
 
      ! Allocate block_offset_sub on all processors (if not already allocated on processor 0)
      if (mype /= 0) then
         allocate(block_offset_sub(nsubleafs))
      end if
 
      ! Broadcast block_offset_sub from processor 0 to all processors
      call mpi_bcast(block_offset_sub, nsubleafs, mpi_offset, 0, icomm, ierrmpi)
 
      amode=ior(mpi_mode_create,mpi_mode_wronly)
      call mpi_file_open(icomm,filename,amode,mpi_info_null,slice_fh,ierrmpi)
      iorequest=mpi_request_null
      iwrite=0
 
      do jgrid=1,njgrid
         iwrite=iwrite+1
         offset=block_offset_sub(int(morton_sub_start(mype)+jgrid-1))
         call mpi_file_iwrite_at(slice_fh,offset,0,2*(ndim-1),mpi_integer,iorequest(iwrite),ierrmpi) 
         iwrite=iwrite+1
         offset1=offset+int(2*(ndim-1)*size_int,kind=mpi_offset_kind)
         call mpi_file_iwrite_at(slice_fh,offset1,ps_sub(jgrid)%w,1,type_subblock_io, &
              iorequest(iwrite),ierrmpi)
      end do
 
      if (iwrite>0) call mpi_waitall(iwrite,iorequest,iostatus,ierrmpi)
      call mpi_barrier(icomm, ierrmpi)
      call mpi_file_close(slice_fh,ierrmpi)
 
      deallocate(block_offset_sub)
      if (mype == 0) then
         deallocate(block_ig_sub)
         deallocate(block_lvl_sub)
         deallocate(isleaf_sub)
      end if
 
    end subroutine put_slice_dat
 
    subroutine put_slice_zerod
 
      use mod_calculate_xw
      integer::  iw
      integer :: amode, iwrite, status(MPI_STATUS_SIZE)
      logical, save :: opened=.false.
      character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
      character(len=1024) :: outfilehead
      character(len=1024) ::line, data
      character(len=1024) :: filename
 
      {^ifoned
      ! generate filename: 
      write(filename,"(a,i1.1,a,i4.4,a)") trim(base_filename)//'_d',dir,'_x'//slice_label(xslice)//'.csv'
 
      ! Open for header:       
      if (.not. opened .and. mype==0) then
         amode=ior(mpi_mode_create,mpi_mode_wronly)
         amode=ior(amode,mpi_mode_append)
         call mpi_file_open(mpi_comm_self,filename,amode, &
              mpi_info_null,slice_fh,ierrmpi)
         ! Create the header:
         call getheadernames(wnamei,xandwnamei,outfilehead)
         line = 'time'
         do iw=1,nw+nwauxio+ndim
            line = trim(line)//', '//trim(xandwnamei(iw))
         enddo
         ! Write header:
         call mpi_file_write(slice_fh,line,len_trim(line), &
              mpi_character,status,ierrmpi)
         call mpi_file_write(slice_fh,achar(10),1,mpi_character,status,ierrmpi)
         call mpi_file_close(slice_fh, ierrmpi)
         opened=.true.
      endif ! opened
      ! Now we let the processor, which holds the data, write - therefore barrier.
      call mpi_barrier(icomm,ierrmpi)
 
      ! There should be only one processor holding the point: 
      if (njgrid > 0) then
         amode=ior(mpi_mode_create,mpi_mode_wronly)
         amode=ior(amode,mpi_mode_append)
         call mpi_file_open(mpi_comm_self,filename,amode, &
              mpi_info_null,slice_fh,ierrmpi)
 
         ! Format the line:
         write(data,"(es14.6)")global_time
         line =  trim(data)
         write(data,"(es14.6)")ps_sub(1)%x
         line =  trim(line)//', '//trim(data)
         do iw = 1,nw+nwauxio
            write(data,"(es14.6)")ps_sub(1)%w(iw)*normconv(iw)
            line = trim(line)//', '//trim(data)
         end do
         !
         call mpi_file_write(slice_fh,trim(line),len_trim(line), &
              mpi_character,status,ierrmpi)
         call mpi_file_write(slice_fh,achar(10),1,mpi_character,status,ierrmpi)
         call mpi_file_close(slice_fh, ierrmpi)
      end if ! Njgrid > 0
      }
    end subroutine put_slice_zerod
 
  end subroutine put_slice
 
  subroutine select_slice(dir,xslice,writeonly,file_handle,normconv)
    use mod_forest, only: iglevel1_sfc, sfc_iglevel1, tree_node_ptr, tree_root, morton_sub_start, morton_sub_stop
    use mod_global_parameters
    integer, intent(in) :: dir
    double precision, intent(in) :: xslice
    integer, intent(in) :: file_handle
    logical, intent(in) :: writeonly
    double precision,dimension(0:nw+nwauxio),intent(out)       :: normconv 
    ! .. local ..
    integer :: ig^D, jgrid, slice_fh, ipe, mylevmax
    integer, dimension(nlevelshi) :: igslice
    integer, allocatable :: iglevel1_sub_sfc(:)
    integer :: ig^L, ngsub1, jglevel1, igsorted, kgrid, ngmax
 
    kgrid = 0
    jgrid = 0
    mylevmax = 0
    jglevel1 = 0
 
    ! Find the global slice index for every level:
    call get_igslice(dir,xslice,igslice)
    igmin^d=1;igmax^d=ng^d(1);
 
    ! Traverse forest to find grids indicating slice:
    select case(dir)
    {case(^d)
      igmin^d=igslice(1);
      igmax^d=igslice(1);
    \}
    case default
       call mpistop("slice direction not clear in select_slice")
    end select
 
    ! total number of level 1 blocks in the slice
    ngsub1={(igmax^d-igmin^d+1)|*}
    allocate(iglevel1_sub_sfc(ngsub1))
    if (.not. allocated(sfc_sub)) then 
      {^ifthreed
      ngmax = ngsub1*(4**refine_max_level-1)
      \}
      {^iftwod
      ngmax = ngsub1*(2**refine_max_level-1)
      \}
      {^ifoned
      ngmax = ngsub1
      \}
      allocate(sfc_sub(0:ngmax,1:3)) ! the first three are nleafs, nparents and isleaf (0 or 1)
      sfc_sub = 0
    end if
 
    ! get the sfc index of the level 1 blocks in the slice
    {do ig^db=igmin^db,igmax^db\}
       jglevel1=jglevel1+1
       iglevel1_sub_sfc(jglevel1)=iglevel1_sfc(ig^d)
    {end do \}
 
    ! Sort the array from small to big using quicksort
    call quicksort_integer(iglevel1_sub_sfc, 1, ngsub1)
 
    ! traverse level 1 blocks depending on the sub morton order
    do igsorted=1,ngsub1
       {ig^d=sfc_iglevel1(^d,iglevel1_sub_sfc(igsorted))\}
       call traverse_slice(tree_root(ig^d))
    end do
   deallocate(iglevel1_sub_sfc)
 
    if (.not.writeonly) then
       ! Synchronize the levmax_sub for output (only rank 0 needs it): 
       levmax_sub = mylevmax
       call mpi_allreduce(mpi_in_place,levmax_sub,1,mpi_integer,mpi_max,icomm,ierrmpi)
 
       ! Communicate the subgrid indices according to new Morton sub-sfc:
       morton_sub_start(:) = 1
       do ipe=0,npe-1
          call mpi_gather(jgrid,1,mpi_integer,morton_sub_stop,1,mpi_integer,ipe,icomm,ierrmpi)
       end do
 
       do ipe = 0, npe-2
          morton_sub_start(ipe+1) = morton_sub_stop(ipe)+morton_sub_start(ipe+1)
          morton_sub_stop(ipe+1)  = morton_sub_stop(ipe)+morton_sub_stop(ipe+1)
       end do
    end if
 
  contains
 
    recursive subroutine traverse_slice(tree)
      implicit none
      type(tree_node_ptr) :: tree
      integer :: ic^d
      integer, dimension(MPI_STATUS_SIZE) :: status
 
      ! if (writeonly) then
      !    call MPI_FILE_WRITE(file_handle,tree%node%leaf,1,MPI_LOGICAL, &
      !         status,ierrmpi)
      ! end if
      if (mype==0) then
         kgrid = kgrid + 1
         sfc_sub(kgrid,1) = tree%node%igrid
         sfc_sub(kgrid,2) = tree%node%ipe
         if (tree%node%leaf) then
            sfc_sub(0,1) = sfc_sub(0,1) + 1
            sfc_sub(kgrid,3) = 1
         else
            sfc_sub(0,2) = sfc_sub(0,2) + 1
            sfc_sub(kgrid,3) = 0
         end if
      end if
 
      if (tree%node%leaf) then
         if (tree%node%ipe == mype.and..not.writeonly) then
            mylevmax = max(mylevmax,tree%node%level)
            call fill_subnode(tree%node%igrid,tree%node%active,jgrid,dir,xslice,normconv)
         end if
         return
      end if
      ! We are out for leaves now, continue for branches
 
      ! Get the correct child:
      select case (dir)
         {case (^d)
         ic^d = igslice(tree%node%level+1) - 2 * tree%node%ig^d + 2
         \}
      case default
         call mpistop("slice direction not clear in traverse_slice")
      end select
 
      ! Recursively descend into the correct child branch:
      {^ifthreed
      select case(dir)
      case (1)
         do ic3=1,2
            do ic2=1,2
               call traverse_slice(tree%node%child(ic^d))
            end do
         end do
      case (2)
         do ic3=1,2
            do ic1=1,2
               call traverse_slice(tree%node%child(ic^d))
            end do
         end do
      case (3)
         do ic2=1,2
            do ic1=1,2
               call traverse_slice(tree%node%child(ic^d))
            end do
         end do
      case default
         call mpistop("slice direction not clear in traverse_slice")
      end select
      }
      {^iftwod
      select case(dir)
      case (1)
         do ic2=1,2
            call traverse_slice(tree%node%child(ic^d))
         end do
      case (2)
         do ic1=1,2
            call traverse_slice(tree%node%child(ic^d))
         end do
      case default
         call mpistop("slice direction not clear in traverse_slice")
      end select
      }
      {^ifoned
      call traverse_slice(tree%node%child(ic^d))
      }
 
    end subroutine traverse_slice
 
  end subroutine select_slice
 
  subroutine fill_subnode(igrid,active,jgrid,dir,xslice,normconv)
    use mod_global_parameters
    use mod_calculate_xw
    integer, intent(in)                                       :: igrid, dir
    integer, intent(inout)                                    :: jgrid
    logical, intent(in)                                       :: active
    double precision, intent(in)                              :: xslice
    double precision,dimension(0:nw+nwauxio),intent(out)      :: normconv 
    ! .. local ..
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim)       :: xC_TMP, xC
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)         :: xCC_TMP, xCC
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)   :: wCC_TMP
    double precision, dimension(ixG^T)                        :: x_save
    integer                :: ixslice, nwexpand, ixC^L, ixCC^L
    logical                :: mask(ixG^T)
 
    mask=.false.
    ! increase grid-count:
    jgrid=jgrid+1
    ! Allocate subdim solution array:
    if (ndim==1) then
       nwexpand = nwauxio
    else
       if(slice_type/='dat')then
          nwexpand = nwauxio
       else
          nwexpand = 0
       endif
    end if
    call alloc_subnode(jgrid,dir,nwexpand)
    call fill_subnode_info(igrid,jgrid,dir)
 
    mask(ixm^t)=.true.
 
    ! Now hunt for the index closest to the slice:
    {^ifoned
    ixslice = minloc(dabs(xslice-ps(igrid)%x(:,dir)),1,mask(:))
    }
    {^iftwod
    select case (dir)
    case (1)
       ixslice = minloc(dabs(xslice-ps(igrid)%x(:,ixmlo2,dir)),1,mask(:,ixmlo2))
    case (2)
       ixslice = minloc(dabs(xslice-ps(igrid)%x(ixmlo1,:,dir)),1,mask(ixmlo1,:))
    case default
       call mpistop("slice direction not clear in fill_subnode")
    end select
    }
    {^ifthreed
    select case (dir)
    case (1)
       ixslice = minloc(dabs(xslice-ps(igrid)%x(:,ixmlo2,ixmlo3,dir)),1,mask(:,ixmlo2,ixmlo3))
    case (2)
       ixslice = minloc(dabs(xslice-ps(igrid)%x(ixmlo1,:,ixmlo3,dir)),1,mask(ixmlo1,:,ixmlo3))
    case (3)
       ixslice = minloc(dabs(xslice-ps(igrid)%x(ixmlo1,ixmlo2,:,dir)),1,mask(ixmlo1,ixmlo2,:))
    case default
       call mpistop("slice direction not clear in fill_subnode")
    end select
    }
 
    call calc_x(igrid,xc,xcc)
    ! Set the coordinate to be exactly on the slice:
    xc(:^d&,dir)  = xslice
    xcc(:^d&,dir) = xslice
    call calc_grid(unitslice,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
         ixc^l,ixcc^l,.true.)
    ! CC stands for CellCenter
    ! C  stands for Corner
    
    ! Fill the subdimensional solution and position array:
    {^ifoned
    ps_sub(jgrid)%w(1:nw+nwexpand) = wcc_tmp(ixslice,1:nw+nwexpand)
    ps_sub(jgrid)%x(1:ndim) = xcc_tmp(ixslice,1:ndim)
    ps_sub(jgrid)%wC(1:nw+nwexpand) = wc_tmp(ixslice,1:nw+nwexpand)
    ps_sub(jgrid)%xC(1:ndim) = xc_tmp(ixslice,1:ndim)
    }
    {^iftwod
    select case (dir)
    case (1)
       ps_sub(jgrid)%w(ixccmin2:ixccmax2,1:nw+nwexpand) &
            = wcc_tmp(ixslice,ixccmin2:ixccmax2,1:nw+nwexpand)
       ps_sub(jgrid)%x(ixccmin2:ixccmax2,1:ndim) &
            = xcc_tmp(ixslice,ixccmin2:ixccmax2,1:ndim)
       ps_sub(jgrid)%wC(ixcmin2:ixcmax2,1:nw+nwexpand) &
            = wc_tmp(ixslice,ixcmin2:ixcmax2,1:nw+nwexpand)
       ps_sub(jgrid)%xC(ixcmin2:ixcmax2,1:ndim) &
            = xc_tmp(ixslice,ixcmin2:ixcmax2,1:ndim)
    case (2)
       ps_sub(jgrid)%w(ixccmin1:ixccmax1,1:nw+nwexpand) &
            = wcc_tmp(ixccmin1:ixccmax1,ixslice,1:nw+nwexpand)
       ps_sub(jgrid)%x(ixccmin1:ixccmax1,1:ndim) &
            = xcc_tmp(ixccmin1:ixccmax1,ixslice,1:ndim)
       ps_sub(jgrid)%wC(ixcmin1:ixcmax1,1:nw+nwexpand) &
            = wc_tmp(ixcmin1:ixcmax1,ixslice,1:nw+nwexpand)
       ps_sub(jgrid)%xC(ixcmin1:ixcmax1,1:ndim) &
            = xc_tmp(ixcmin1:ixcmax1,ixslice,1:ndim)
    case default
       call mpistop("slice direction not clear in fill_subnode")
    end select
    }
    {^ifthreed
    select case (dir)
    case (1)
       ps_sub(jgrid)%w(ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:nw+nwexpand) = &
            wcc_tmp(ixslice,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:nw+nwexpand)
       ps_sub(jgrid)%x(ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:ndim) = &
            xcc_tmp(ixslice,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:ndim)
       ps_sub(jgrid)%wC(ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw+nwexpand) = &
            wc_tmp(ixslice,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw+nwexpand)
       ps_sub(jgrid)%xC(ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:ndim) = &
            xc_tmp(ixslice,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:ndim)
    case (2)
       ps_sub(jgrid)%w(ixccmin1:ixccmax1,ixccmin3:ixccmax3,1:nw+nwexpand) = &
            wcc_tmp(ixccmin1:ixccmax1,ixslice,ixccmin3:ixccmax3,1:nw+nwexpand)
       ps_sub(jgrid)%x(ixccmin1:ixccmax1,ixccmin3:ixccmax3,1:ndim) = &
            xcc_tmp(ixccmin1:ixccmax1,ixslice,ixccmin3:ixccmax3,1:ndim)
       ps_sub(jgrid)%wC(ixcmin1:ixcmax1,ixcmin3:ixcmax3,1:nw+nwexpand) = &
            wc_tmp(ixcmin1:ixcmax1,ixslice,ixcmin3:ixcmax3,1:nw+nwexpand)
       ps_sub(jgrid)%xC(ixcmin1:ixcmax1,ixcmin3:ixcmax3,1:ndim) = &
            xc_tmp(ixcmin1:ixcmax1,ixslice,ixcmin3:ixcmax3,1:ndim) 
    case (3)
       ps_sub(jgrid)%w(ixccmin1:ixccmax1,ixccmin2:ixccmax2,1:nw+nwexpand) = &
            wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixslice,1:nw+nwexpand) 
       ps_sub(jgrid)%x(ixccmin1:ixccmax1,ixccmin2:ixccmax2,1:ndim) = &
            xcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixslice,1:ndim)
       ps_sub(jgrid)%wC(ixcmin1:ixcmax1,ixcmin2:ixcmax2,1:nw+nwexpand) = &
            wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixslice,1:nw+nwexpand) 
       ps_sub(jgrid)%xC(ixcmin1:ixcmax1,ixcmin2:ixcmax2,1:ndim) = &
            xc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixslice,1:ndim) 
    case default
       call mpistop("slice direction not clear in fill_subnode")
    end select
    }
 
  end subroutine fill_subnode
 
  subroutine alloc_subnode(jgrid,dir,nwexpand)
    use mod_global_parameters
    integer, intent(in) :: jgrid, dir, nwexpand
 
    ! take care, what comes out is not necessarily a right handed system!
    {^ifoned
    allocate(ps_sub(jgrid)%w(1:nw+nwexpand),ps_sub(jgrid)%x(1:ndim))
    allocate(ps_sub(jgrid)%wC(1:nw+nwexpand),ps_sub(jgrid)%xC(1:ndim))
    }
    {^iftwod
    select case (dir)
    case (1)
       allocate(ps_sub(jgrid)%w(ixglo2:ixghi2,1:nw+nwexpand),&
            ps_sub(jgrid)%x(ixglo2:ixghi2,1:ndim), &
            ps_sub(jgrid)%wC(ixglo2:ixghi2,1:nw+nwexpand),&
            ps_sub(jgrid)%xC(ixglo2:ixghi2,1:ndim))
    case (2)
       allocate(ps_sub(jgrid)%w(ixglo1:ixghi1,1:nw+nwexpand),&
            ps_sub(jgrid)%x(ixglo1:ixghi1,1:ndim), &
            ps_sub(jgrid)%wC(ixglo1:ixghi1,1:nw+nwexpand),&
            ps_sub(jgrid)%xC(ixglo1:ixghi1,1:ndim))
    case default
       call mpistop("slice direction not clear in alloc_subnode")
    end select
    }
    {^ifthreed
    select case (dir)
    case (1)
       allocate(ps_sub(jgrid)%w(ixglo2:ixghi2,ixglo3:ixghi3,1:nw+nwexpand),&
            ps_sub(jgrid)%x(ixglo2:ixghi2,ixglo3:ixghi3,1:ndim), &
            ps_sub(jgrid)%wC(ixglo2:ixghi2,ixglo3:ixghi3,1:nw+nwexpand),&
            ps_sub(jgrid)%xC(ixglo2:ixghi2,ixglo3:ixghi3,1:ndim))
    case (2)
       allocate(ps_sub(jgrid)%w(ixglo1:ixghi1,ixglo3:ixghi3,1:nw+nwexpand),&
            ps_sub(jgrid)%x(ixglo1:ixghi1,ixglo3:ixghi3,1:ndim), &
            ps_sub(jgrid)%wC(ixglo1:ixghi1,ixglo3:ixghi3,1:nw+nwexpand),&
            ps_sub(jgrid)%xC(ixglo1:ixghi1,ixglo3:ixghi3,1:ndim))
    case (3)
       allocate(ps_sub(jgrid)%w(ixglo1:ixghi1,ixglo2:ixghi2,1:nw+nwexpand),&
            ps_sub(jgrid)%x(ixglo1:ixghi1,ixglo2:ixghi2,1:ndim), &
            ps_sub(jgrid)%wC(ixglo1:ixghi1,ixglo2:ixghi2,1:nw+nwexpand),&
            ps_sub(jgrid)%xC(ixglo1:ixghi1,ixglo2:ixghi2,1:ndim))
    case default
       call mpistop("slice direction not clear in alloc_subnode")
    end select
    }
  end subroutine alloc_subnode
 
  subroutine dealloc_subnode(jgrid)
    use mod_global_parameters
    integer, intent(in) :: jgrid
 
    if (jgrid==0) then
       call mpistop("trying to delete a non-existing grid in dealloc_subnode")
    end if
 
    deallocate(ps_sub(jgrid)%w,ps_sub(jgrid)%x,ps_sub(jgrid)%wC,ps_sub(jgrid)%xC)
 
    ! reset the global node info:
    node_sub(:,jgrid)=0
    rnode_sub(:,jgrid)=zero
 
  end subroutine dealloc_subnode
 
  subroutine fill_subnode_info(igrid,jgrid,dir)
    use mod_global_parameters
    integer, intent(in) :: igrid,jgrid,dir
 
    node_sub(plevel_,jgrid)=node(plevel_,igrid)
    {^ifthreed
    select case(dir)
    case (1)
       node_sub(pig1_,jgrid)=node(pig2_,igrid)
       node_sub(pig2_,jgrid)=node(pig3_,igrid)
       rnode_sub(rpdx1_,jgrid)=rnode(rpdx2_,igrid)
       rnode_sub(rpdx2_,jgrid)=rnode(rpdx3_,igrid)
       rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin2_,igrid)
       rnode_sub(rpxmin2_,jgrid)=rnode(rpxmin3_,igrid)
    case (2)
       node_sub(pig1_,jgrid)=node(pig1_,igrid)
       node_sub(pig2_,jgrid)=node(pig3_,igrid)
       rnode_sub(rpdx1_,jgrid)=rnode(rpdx1_,igrid)
       rnode_sub(rpdx2_,jgrid)=rnode(rpdx3_,igrid)
       rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin1_,igrid)
       rnode_sub(rpxmin2_,jgrid)=rnode(rpxmin3_,igrid)
    case (3)
       node_sub(pig1_,jgrid)=node(pig1_,igrid)
       node_sub(pig2_,jgrid)=node(pig2_,igrid)
       rnode_sub(rpdx1_,jgrid)=rnode(rpdx1_,igrid)
       rnode_sub(rpdx2_,jgrid)=rnode(rpdx2_,igrid)
       rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin1_,igrid)
       rnode_sub(rpxmin2_,jgrid)=rnode(rpxmin2_,igrid)
    case default
       call mpistop("slice direction not clear in fill_subnode_info")
    end select
    }
    {^iftwod
    select case(dir)
    case (1)
       node_sub(pig1_,jgrid)=node(pig2_,igrid)
       rnode_sub(rpdx1_,jgrid)=rnode(rpdx2_,igrid)
       rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin2_,igrid)
    case (2)
       node_sub(pig1_,jgrid)=node(pig1_,igrid)
       rnode_sub(rpdx1_,jgrid)=rnode(rpdx1_,igrid)
       rnode_sub(rpxmin1_,jgrid)=rnode(rpxmin1_,igrid)
    case default
       call mpistop("slice direction not clear in fill_subnode_info")
    end select
    }
    {^ifoned
    }
 
  end subroutine fill_subnode_info
 
  subroutine get_igslice(dir,x,igslice)
    use mod_global_parameters
    integer, intent(in) :: dir
    double precision, intent(in) :: x
    integer, dimension(nlevelshi), intent(out) :: igslice
    ! .. local ..
    double precision :: distance
    integer :: level
    !double precision :: xsgrid(ndim,nlevelshi),qs(ndim),xmgrid(ndim,3),xnew,xlgrid(ndim,3)
    !integer :: nbefore
 
    if (x.ne.x) &
         call mpistop("get_igslice: your slice position is NaN!")
 
    select case (dir)
    {case (^d)
      select case (stretch_type(^d))
      case (stretch_none) ! This dimension is unstretched
        if (x<=xprobmin^d) then
          igslice=1
        else if (x>=xprobmax^d) then
          do level=1,refine_max_level
            igslice(level)=ng^d(level)
          end do
        else
          distance=x-xprobmin^d
          do level = 1, refine_max_level
            igslice(level) = ceiling(distance/dg^d(level))
          end do
        end if
 
      case (stretch_uni) ! Uniform stretching
 
        if (x<=xprobmin^d) then
          igslice=1
        else if (x>=xprobmax^d) then
          do level=1,refine_max_level
            igslice(level)=ng^d(level)
          end do
        else
          distance=x-xprobmin^d
          do level=1,refine_max_level
            igslice(level)= ceiling(dlog(distance/dxfirst(level,^d)*(qstretch(level,^d)-1.d0)+1.d0)&
              /dlog(qstretch(level,^d))/dble(block_nx^d))
          end do
        end if
 
      case (stretch_symm) ! Symmetric stretching
        ! symmetric stretch about 0.5*(xprobmin+xprobmax)
        if (x<=xprobmin^d) then
          igslice=1
        else if (x>=xprobmax^d) then
          do level=1,refine_max_level
            igslice(level)=ng^d(level)
          end do
        else if(x<xprobmin^d+xstretch^d) then
          distance=xprobmin^d+xstretch^d-x
          ! stretch to left from xprobmin+xstretch
          do level=1,refine_max_level
            igslice(level) = nstretchedblocks(level,^d)/2-int(dlog(distance*(qstretch(level,^d)-one)/&
              dxfirst(level,^d)+one)/dlog(qstretch(level,^d))/dble(block_nx^d))
          end do
        else if(x>xprobmax^d-xstretch^d) then
          distance=x-xprobmax^d+xstretch^d
          ! stretch to right from xprobmax-xstretch
          do level=1,refine_max_level
            igslice(level) = ceiling(dlog(distance*(qstretch(level,^d)-one)/&
              dxfirst(level,^d)+one)/dlog(qstretch(level,^d))/dble(block_nx^d))+ng^d(level)-nstretchedblocks(level,^d)/2
          end do
        else
          ! possible non-stretched central part
          distance=x-xprobmin^d-xstretch^d
          do level=1,refine_max_level
            igslice(level)=nstretchedblocks(level,^d)/2+ceiling(distance/dg^d(level))
          end do
        end if
        !xsgrid(^D,1)=half*(xprobmax^D-xprobmin^D)&
        !        /(half*domain_nx^D-half*nstretchedblocks_baselevel(^D)*block_nx^D &
        !        +(one-qstretch_baselevel(^D)**(half*nstretchedblocks_baselevel(^D)*block_nx^D)) &
        !        /(one-qstretch_baselevel(^D)))
        !xmgrid(^D,1)=xsgrid(^D,1)*(domain_nx^D-nstretchedblocks_baselevel(^D)*block_nx^D)
        !xlgrid(^D,1)=xsgrid(^D,1)
 
        !if (x .ge. xmgrid(^D,1)*half) then
        !  xnew=x-xmgrid(^D,1)*half
        !  do level=1,refine_max_level
        !    qs(^D)=qstretch_baselevel(^D)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(^D,level)=xsgrid(^D,level-1)/(one+qs(^D))
        !    nbefore=(domain_nx^D/block_nx^D-nstretchedblocks_baselevel(^D)/2)*2**(level-1)
        !    igslice(level)=floor((dlog(one-xnew/xsgrid(^D,level)*(one-qs(^D)))/dlog(qs(^D))-1.e-16)/block_nx^D)+1+nbefore
        !    if (x>=xprobmax^D) igslice(level)=domain_nx^D*2**(level-1)
        !  end do
        !else if(x .ge. -xmgrid(^D,1)*half) then
        !  xnew=x+xmgrid(^D,1)*half
        !  do level=1,refine_max_level
        !    nbefore=nstretchedblocks_baselevel(^D)/2*2**(level-1)
        !    if (level .gt. 1) xlgrid(^D,level)=xlgrid(^D,level-1)/2.d0
        !    igslice(level)=floor(((xnew-1.e-16)/xlgrid(^D,level))/block_nx^D)+1+nbefore
        !  end do
        !else
        !  xnew=xmgrid(^D,1)*half-x
        !  do level=1,refine_max_level
        !    qs(^D)=qstretch_baselevel(^D)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(^D,level)=xsgrid(^D,level-1)/(one+qs(^D))
        !    igslice(level)=nstretchedblocks_baselevel(^D)/2*2**(level-1) &
        !      -floor((dlog(one-xnew/xsgrid(^D,level)*(one-qs(^D)))/dlog(qs(^D)))/block_nx^D)
        !    if (x<=xprobmin^D) igslice(level)=1
        !  end do
        !end if
 
      case default
        call mpistop("stretch type not supported by get_igslice")
      end select 
    \}
    case default
      call mpistop("slice direction not clear in get_igslice")
    end select
 
  end subroutine get_igslice
 
  double precision function roundoff_minmax(val,minval,maxval)
    implicit none
    double precision,intent(in)         :: val, minval, maxval
 
    roundoff_minmax = val
 
    if (abs(roundoff_minmax) .le. minval) then
       roundoff_minmax = 0.0d0
    end if
 
    if (roundoff_minmax .gt. maxval) then
       roundoff_minmax = maxval
    else if (roundoff_minmax .lt. -maxval) then
       roundoff_minmax = -maxval
    end if
 
    ! Replace NaN with maxval (e.g. Paraview chokes on ASCII NaN): 
    if (roundoff_minmax /= roundoff_minmax) roundoff_minmax = maxval
 
  end function roundoff_minmax
 
  recursive subroutine quicksort_integer(arr, left, right)
    implicit none
    integer, intent(inout) :: arr(:)
    integer, intent(in) :: left, right
    ! .. local ..
    integer :: pivot, i, j, temp
    
    if (left < right) then
       ! Choose pivot (using middle element)
       pivot = arr((left + right) / 2)
       i = left - 1
       j = right + 1
       
       do
          ! Find elements to swap
          do
             i = i + 1
             if (arr(i) >= pivot) exit
          end do
          
          do
             j = j - 1
             if (arr(j) <= pivot) exit
          end do
          
          if (i >= j) exit
          
          ! Swap elements
          temp = arr(i)
          arr(i) = arr(j)
          arr(j) = temp
       end do
       
       ! Recursively sort subarrays
       call quicksort_integer(arr, left, j)
       call quicksort_integer(arr, j + 1, right)
    end if
    
  end subroutine quicksort_integer
 
 
  !> Signed slice coordinate for the file names, such as +0.500000D+00
  character(len=13) function slice_label(xslice)
    double precision, intent(in) :: xslice
    write(slice_label,"(D13.6)") xslice
    if(xslice >= 0.0d0) slice_label(1:1) = "+"
  end function slice_label
 
end module mod_slice