mod_input_output.t

Go to the documentation of this file.

!> Module for reading input and writing output
module mod_input_output
  use mod_comm_lib, only: mpistop
  use mod_basic_types, only: name_len
 
  implicit none
  public
 
 
  !> coefficient for rk2_alfa
  double precision, private :: rk2_alfa
 
  !> debug field-dump variable names (set by debug_alloc, used by save_wdebug)
  character(len=name_len), allocatable :: wdebug_names(:)
 
  !> List of compatible versions
  integer, parameter :: compatible_versions(3) = [3, 4, 5]
 
  !> number of w found in dat files
  integer :: nw_found
 
  !> tag for MPI message
  integer, private :: itag
 
  !> whether staggered field is in dat
  logical, private :: stagger_mark_dat=.false.
 
  ! Formats used in output
  character(len=*), parameter :: fmt_r  = 'es16.8' ! Default precision
  character(len=*), parameter :: fmt_r2 = 'es10.2' ! Two digits
  character(len=*), parameter :: fmt_i  = 'i8'     ! Integer format
 
  ! public methods
  public :: snapshot_write_header  
 
contains
 
  !> Read the command line arguments passed to amrvac
  subroutine read_arguments()
    use mod_global_parameters
 
    integer                          :: len, stat, n, i, ipars
    integer, parameter               :: max_files = 20 ! Maximum number of par files
    integer                          :: n_par_files
    logical                          :: unknown_arg, help, morepars
    character(len=max_files*std_len) :: all_par_files
    character(len=std_len)           :: tmp_files(max_files), arg
 
    if (mype == 0) then
       print *, '-----------------------------------------------------------------------------'
       print *, '-----------------------------------------------------------------------------'
       print *, '|         __  __ ____ ___        _    __  __ ______     ___    ____         |'
       print *, '|        |  \/  |  _ \_ _|      / \  |  \/  |  _ \ \   / / \  / ___|        |'
       print *, '|        | |\/| | |_) | |_____ / _ \ | |\/| | |_) \ \ / / _ \| |            |'
       print *, '|        | |  | |  __/| |_____/ ___ \| |  | |  _ < \ V / ___ \ |___         |'
       print *, '|        |_|  |_|_|  |___|   /_/   \_\_|  |_|_| \_\ \_/_/   \_\____|        |'
       print *, '-----------------------------------------------------------------------------'
       print *, '-----------------------------------------------------------------------------'
    end if
 
    ! =============== Fortran 2003 command line reading ================
 
    ! Default command line arguments
    all_par_files="amrvac.par"
    restart_from_file=undefined
    snapshotnext=-1
    slicenext=-1
    collapsenext=-1
    help=.false.
    convert=.false.
    resume_previous_run=.false.
 
    ! Argument 0 is program name, so we start from one
    i = 1
    unknown_arg=.false.
    DO
       CALL get_command_argument(i, arg)
       IF (len_trim(arg) == 0) EXIT
       select case(arg)
         case("-i")
           i = i+1
           CALL get_command_argument(i, arg)
           !if(mype==0)print *,'found argument=',arg
           all_par_files=trim(arg)
           morepars=.true.
           ipars=1
           do while (ipars<max_files.and.morepars)
              CALL get_command_argument(i+ipars,arg)
              !if(mype==0)print *,'found argument=',arg
              if (index(trim(arg),"-")==1.or.len_trim(arg)==0) then
                 morepars=.false.
              else
                 ipars=ipars+1
                 all_par_files=trim(all_par_files)//" "//trim(arg)
                 !if(mype==0)print *,'now all_par_files=',TRIM(all_par_files)
              endif
           end do
           !if(mype==0)print *,'-i identified ',ipars, ' par-file arguments'
           i=i+ipars-1
           !if(mype==0)print *,'-i arguments passed to all_par_files=',TRIM(all_par_files)
         case("-if")
           i = i+1
           CALL get_command_argument(i, arg)
           restart_from_file=trim(arg)
           !if(mype==0)print *,'-if has argument=',TRIM(arg),' passed to restart_from_file=',TRIM(restart_from_file)
         case("-slicenext")
           i = i+1
           CALL get_command_argument(i, arg)
           read(arg,*,iostat=stat) slicenext
           !if(mype==0)print *,'-slicenext has argument=',arg,' passed to slicenext=',slicenext
         case("-collapsenext")
           i = i+1
           CALL get_command_argument(i, arg)
           read(arg,*,iostat=stat) collapsenext
           !if(mype==0)print *,'-collapsenext has argument=',arg,' passed to collapsenext=',collapsenext
         case("-snapshotnext")
           i = i+1
           CALL get_command_argument(i, arg)
           read(arg,*,iostat=stat) snapshotnext
           !if(mype==0)print *,'-snapshotnext has argument=',arg,' passed to snapshotnext=',snapshotnext
         case("-resume")
           resume_previous_run=.true.
           !if(mype==0)print *,'resume specified: resume_previous_run=T'
         case("-convert")
           convert=.true.
           if(mype==0)print *,'convert specified: convert=T'
         case("--help","-help")
           help=.true.
           EXIT
         case default
           unknown_arg=.true.
           help=.true.
           EXIT
       end select
       i = i+1
    END DO
 
    if (unknown_arg) then
       print*,"======================================="
       print*,"Error: Command line argument ' ",trim(arg)," ' not recognized"
       print*,"======================================="
       help=.true.
    end if
 
    ! Show the usage if the help flag was given, or no par file was specified
    if (help) then
       if (mype == 0) then
          print *, 'Usage example:'
          print *, 'mpirun -np 4 ./amrvac -i file.par [file2.par ...]'
          print *, '         (later .par files override earlier ones)'
          print *, ''
          print *, 'Optional arguments:'
          print *, '-convert             Convert snapshot files'
          print *, '-if file0001.dat     Use this snapshot to restart from'
          print *, '                     (you can modify e.g. output names)'
          print *, '-resume              Automatically resume previous run'
          print *, '                     (useful for long runs on HPC systems)'
          print *, '-snapshotnext N      Manual index for next snapshot'
          print *, '-slicenext N         Manual index for next slice output'
          print *, '-collapsenext N      Manual index for next collapsed output'
          print *, ''
       end if
       call mpi_finalize(ierrmpi)
       stop
    end if
 
    ! Split the input files, in case multiple were given
    call get_fields_string(all_par_files, " ,'"""//char(9), max_files, &
         tmp_files, n_par_files)
 
    allocate(par_files(n_par_files))
    par_files = tmp_files(1:n_par_files)
 
  end subroutine read_arguments
 
  !> Read in the user-supplied parameter-file
  subroutine read_par_files()
    use mod_global_parameters
    use mod_physics, only: phys_energy, phys_total_energy, physics_type, phys_wider_stencil
    use mod_small_values
    use mod_limiter
    use mod_slice
    use mod_geometry
    use mod_source
    use mod_input_output_helper, only: get_names_from_string
    use mod_timing, only: write_timing_log, timing_log_interval
 
    double precision, dimension(nsavehi) :: tsave_log, tsave_dat, tsave_slice, &
         tsave_collapsed, tsave_custom
    double precision :: dtsave_log, dtsave_dat, dtsave_slice, &
         dtsave_collapsed, dtsave_custom
    double precision :: tsavestart_log, tsavestart_dat, tsavestart_slice, &
         tsavestart_collapsed, tsavestart_custom
    double precision :: sizeuniformpart^D
    double precision :: im_delta,im_nu,rka54,rka51,rkb54,rka55
    double precision :: dx_vec(^ND)
    integer :: ditsave_log, ditsave_dat, ditsave_slice, &
         ditsave_collapsed, ditsave_custom
    integer :: windex, ipower
    integer :: i, j, k, ifile, io_state
    integer :: iB, isave, iw, level, idim, islice
    integer :: nx_vec(^ND), block_nx_vec(^ND)
    integer :: my_unit, iostate
    integer :: ilev
    logical :: fileopen, file_exists
 
    character              :: c_ndim
    character(len=80)      :: fmt_string
    character(len=std_len) :: err_msg, restart_from_file_arg
    character(len=std_len) :: basename_full, basename_prev, dummy_file
    character(len=std_len), dimension(:), allocatable :: &
         typeboundary_min^D, typeboundary_max^D
    character(len=std_len), allocatable :: limiter(:)
    character(len=std_len), allocatable :: gradient_limiter(:)
    character(len=name_len) :: stretch_dim(ndim)
    !> How to apply dimensional splitting to the source terms, see
    !> @ref discretization.md
    character(len=std_len) :: typesourcesplit
    !> Which flux scheme of spatial discretization to use (per grid level)
    character(len=std_len), allocatable :: flux_scheme(:)
    !> Which type of the maximal bound speed of Riemann fan to use
    character(len=std_len) :: typeboundspeed
    !> Which time stepper to use
    character(len=std_len) :: time_stepper
    !> Which time integrator to use
    character(len=std_len) :: time_integrator
    !> type of curl operator
    character(len=std_len) :: typecurl
    !> Limiter used for prolongation to refined grids and ghost cells
    character(len=std_len) :: typeprolonglimit
    !> How to compute the CFL-limited time step.
    !> Options are 'maxsum': max(sum(c/dx)); 'summax': sum(max(c/dx)) and
    !> 'minimum: max(c/dx), where the summations loop over the grid dimensions and
    !> c is the velocity. The default 'maxsum' is the conventiontal way of
    !> computing CFL-limited time steps.
    character(len=std_len) :: typecourant
 
 
    namelist /filelist/ base_filename,restart_from_file, &
         typefilelog,firstprocess,reset_grid,snapshotnext, &
         allow_ndir_change, &
         convert,convert_type,saveprim,usr_filename,&
         nwauxio,nocartesian, w_write,writelevel,&
         writespshift,length_convert_factor, w_convert_factor, &
         time_convert_factor,level_io,level_io_min, level_io_max, &
         autoconvert,slice_type,slicenext,collapsenext,collapse_type, &
         type_endian
 
    namelist /savelist/ tsave,itsave,dtsave,ditsave,nslices,slicedir, &
         slicecoord,collapse,collapselevel, time_between_print,&
         tsave_log, tsave_dat, tsave_slice, tsave_collapsed, tsave_custom, &
         dtsave_log, dtsave_dat, dtsave_slice, dtsave_collapsed, dtsave_custom, &
         ditsave_log, ditsave_dat, ditsave_slice, ditsave_collapsed, ditsave_custom,&
         tsavestart_log, tsavestart_dat, tsavestart_slice, tsavestart_collapsed,&
         tsavestart_custom, tsavestart, &
         write_timing_log, timing_log_interval
 
    namelist /stoplist/ it_init,time_init,it_max,time_max,dtmin,reset_it,reset_time,&
         wall_time_max,final_dt_reduction
 
    namelist /methodlist/ time_stepper, time_integrator, &
         source_split_usr, typesourcesplit, local_timestep, &
         dimsplit, typedimsplit, flux_scheme, limit_on_wold, &
         limiter, gradient_limiter, cada3_radius, &
         loglimit, typeboundspeed, h_correction, &
         typetvd, typeentropy, entropycoef, typeaverage, &
         typegrad, typediv, typecurl, &
         nxdiffusehllc, flathllc, tvdlfeps, flux_adaptive_diffusion, &
         flux_adaptive_diffusion_min, flux_adaptive_diffusion_scale, &
         flatcd, flatsh, ppm_avisc, ppm_rjv, &
         rk2_alfa, imex222_lambda, ssprk_order, rk3_switch, imex_switch, &
         small_temperature, small_pressure, small_density, &
         small_values_method, small_values_daverage, fix_small_values, check_small_values, &
         trace_small_values, small_values_fix_iw, minfip, maxfip, &
         ghostcell_comm_batched, ghostcell_comm_batch_size
 
    namelist /boundlist/ nghostcells,ghost_copy,&
         internalboundary, typeboundary_^l, save_physical_boundary
 
    namelist /meshlist/ refine_max_level,nbufferx^d,refine_threshold,&
         derefine_ratio, refine_criterion, &
         stretch_dim, stretch_uncentered, &
         qstretch_baselevel, nstretchedblocks_baselevel, &
         amr_wavefilter,max_blocks,block_nx^d,domain_nx^d,iprob,xprob^l, &
         w_refine_weight, prolongprimitive,coarsenprimitive, &
         typeprolonglimit, &
         logflag,tfixgrid,itfixgrid,ditregrid, &
         lb_diagnose,lb_automatic,lb_interval,lb_alpha,lb_max_block_ratio
    namelist /paramlist/  courantpar, dtpar, dtdiffpar, &
         typecourant, slowsteps
 
    namelist /emissionlist/ filename_euv,wavelength,&
          filename_sxr,emin_sxr,emax_sxr,&
          los_theta,los_phi,image_rotate,x_origin,big_image,&
          spectrum_wl,location_slit,filename_spectrum,&
          radiation_transfer,ray_method,dat_resolution_mode,emission_model,&
          instrument_postprocess,&
          radio_frequency,radio_beam_fwhm,radio_beam_pixel_size,&
          radsyn_pixel_batch,radsyn_segment_batch_factor,radsyn_segment_memory_mb,&
          radsyn_segment_comm_factor,&
          radsyn_verbose,&
          spectrum_window_min,spectrum_window_max,&
          instrument_resolution_factor,activate_unit_arcsec,&
          filename_whitelight,whitelight_instrument,r_occultor,r_opt_thick,&
          dat_resolution,direction_slit,output_tau,output_absorption_fraction
 
    ! default maximum number of grid blocks in a processor
    max_blocks=4000
 
    ! allocate cell size of all levels
    allocate(dx(ndim,nlevelshi))
    {allocate(dg^d(nlevelshi))\}
    {allocate(ng^d(nlevelshi))\}
 
    ! default block size excluding ghost cells
    {block_nx^d = 16\}
 
    ! default resolution of level-1 mesh (full domain)
    {domain_nx^d = 32\}
 
    !! default number of ghost-cell layers at each boundary of a block
    ! this is now done when the variable is defined in mod_global_parameters
    ! the physics modules might set this variable in their init subroutine called earlier
    nghostcells = 2
 
    ! Allocate boundary conditions arrays in new and old style
    {
    allocate(typeboundary_min^d(nwfluxbc))
    allocate(typeboundary_max^d(nwfluxbc))
    typeboundary_min^d = undefined
    typeboundary_max^d = undefined
    }
 
    allocate(typeboundary(nwflux+nwaux,2*ndim))
    typeboundary=0
 
    ! not save physical boundary in dat files by default
    save_physical_boundary = .false.
 
    internalboundary   = .false.
 
    ! defaults for specific options
    typegrad   = 'central'
    typediv    = 'central'
    typecurl   = 'central'
 
    ! defaults for smallest physical values allowed
    small_temperature = 0.d0
    small_pressure    = 0.d0
    small_density     = 0.d0
 
    allocate(small_values_fix_iw(nw))
    small_values_fix_iw(:) = .true.
 
    ! defaults for convert behavior
 
    ! store the -if value from argument in command line
    restart_from_file_arg    = restart_from_file
    nwauxio                  = 0
    nocartesian              = .false.
    saveprim                 = .false.
    autoconvert              = .false.
    convert_type             = 'vtuBCCmpi'
    slice_type               = 'vtuCC'
    collapse_type            = 'vti'
    allocate(w_write(nw))
    w_write(1:nw)             = .true.
    allocate(writelevel(nlevelshi))
    writelevel(1:nlevelshi)  = .true.
    writespshift(1:ndim,1:2) = zero
    level_io                 = -1
    level_io_min             = 1
    level_io_max             = nlevelshi
    ! endianness: littleendian (default) is 1, bigendian otherwise
    type_endian              = 1
 
    ! normalization of primitive variables: only for output
    ! note that length_convert_factor is for length
    ! this scaling is optional, and must be set consistently if used
    allocate(w_convert_factor(nw))
    w_convert_factor(:)   = 1.0d0
    time_convert_factor   = 1.0d0
    length_convert_factor = 1.0d0
 
    ! AMR related defaults
    refine_max_level      = 1
    {nbufferx^d                 = 0\}
    allocate(refine_threshold(nlevelshi))
    refine_threshold(1:nlevelshi)            = 0.1d0
    allocate(derefine_ratio(nlevelshi))
    derefine_ratio(1:nlevelshi)       = 1.0d0/8.0d0
    typeprolonglimit            = 'default'
    refine_criterion               = 3
    allocate(w_refine_weight(nw+1))
    w_refine_weight                    = 0.d0
    allocate(logflag(nw+1))
    logflag                     = .false.
    allocate(amr_wavefilter(nlevelshi))
    amr_wavefilter(1:nlevelshi) = 1.0d-2
    tfixgrid                    = bigdouble
    itfixgrid                   = biginteger
    ditregrid                   = 1
 
    ! Grid stretching defaults
    stretch_uncentered = .true.
    stretch_dim(1:ndim) = undefined
    qstretch_baselevel(1:ndim) = bigdouble
    nstretchedblocks_baselevel(1:ndim)=0
 
    ! IO defaults
    it_init       = 0
    it_max        = biginteger
    time_init     = 0.d0
    time_max      = bigdouble
    wall_time_max = bigdouble
    if(local_timestep) then
      final_dt_reduction=.false.
    else
      final_dt_reduction=.true.
    endif
    final_dt_exit=.false.
    dtmin         = 1.0d-10
    nslices       = 0
    collapse      = .false.
    collapselevel = 1
    time_between_print = 30.0d0 ! Print status every 30 seconds
 
    do ifile=1,nfile
      do isave=1,nsavehi
        tsave(isave,ifile)  = bigdouble  ! global_time  of saves into the output files
        itsave(isave,ifile) = biginteger ! it of saves into the output files
      end do
      dtsave(ifile)  = bigdouble  ! time between saves
      ditsave(ifile) = biginteger ! timesteps between saves
      isavet(ifile)  = 1          ! index for saves by global_time
      isaveit(ifile) = 1          ! index for saves by it
      tsavestart(ifile) = 0.0d0
    end do
 
    tsave_log       = bigdouble
    tsave_dat       = bigdouble
    tsave_slice     = bigdouble
    tsave_collapsed = bigdouble
    tsave_custom    = bigdouble
 
    dtsave_log       = bigdouble
    dtsave_dat       = bigdouble
    dtsave_slice     = bigdouble
    dtsave_collapsed = bigdouble
    dtsave_custom    = bigdouble
 
    ditsave_log       = biginteger
    ditsave_dat       = biginteger
    ditsave_slice     = biginteger
    ditsave_collapsed = biginteger
    ditsave_custom    = biginteger
 
    tsavestart_log       = bigdouble
    tsavestart_dat       = bigdouble
    tsavestart_slice     = bigdouble
    tsavestart_collapsed = bigdouble
    tsavestart_custom    = bigdouble
 
    typefilelog = 'default'
 
    ! defaults for input
    reset_time = .false.
    reset_it = .false.
    firstprocess  = .false.
    reset_grid     = .false.
    allow_ndir_change = .false.
    base_filename   = 'data'
    usr_filename    = ''
 
 
    ! Defaults for discretization methods
    typeaverage     = 'default'
    tvdlfeps        = one
    flux_adaptive_diffusion = .false.
    flux_adaptive_diffusion_min = 0.d0
    flux_adaptive_diffusion_scale = 1.d0
    nxdiffusehllc   = 0
    flathllc        = .false.
    slowsteps       = -1
    courantpar      = 0.8d0
    typecourant     = 'maxsum'
    dimsplit        = .false.
    typedimsplit    = 'default'
    if(physics_type=='mhd') then
      cada3_radius  = 0.1d0
    else
      cada3_radius  = 0.1d0
    end if
    typetvd         = 'roe'
    typeboundspeed  = 'Einfeldt'
    source_split_usr= .false.
    time_stepper    = 'twostep'
    time_integrator = 'default'
    limit_on_wold   = .false.
    ! default PC or explicit midpoint, hence alfa=0.5
    rk2_alfa        = half 
    ! default IMEX-RK22Ln hence lambda = 1 - 1/sqrt(2)
    imex222_lambda  = 1.0d0 - 1.0d0 / dsqrt(2.0d0)
    ! default SSPRK(3,3) or Gottlieb-Shu 1998 for threestep
    ! default SSPRK(4,3) or Spireti-Ruuth for fourstep
    ! default SSPRK(5,4) using Gottlieb coeffs
    ssprk_order    = 3
    ! default RK3 butcher table: Heun 3rd order
    rk3_switch      = 3
    ! default IMEX threestep is IMEX_ARK(232)
    imex_switch     = 1
 
    ! Defaults for synthesing emission
    los_theta = 0.d0
    los_phi = 0.d0
    image_rotate = 0.d0
    x_origin = 0.d0
    big_image = .false.
    location_slit = 0.d0
    radiation_transfer = 'thin'
    ray_method = 'auto'
    dat_resolution_mode = 'nominal'
    emission_model = 'auto'
    instrument_postprocess=.false.
    radio_frequency = 17.d9
    radio_beam_fwhm = 0.d0
    radio_beam_pixel_size = 0.d0
    radsyn_pixel_batch = 128
    radsyn_segment_batch_factor = 256
    radsyn_segment_memory_mb = 8.d0
    radsyn_segment_comm_factor = 16
    radsyn_verbose = .false.
    direction_slit = -1
    instrument_resolution_factor=1.d0
    activate_unit_arcsec=.true.
    whitelight_instrument='LASCO/C2'
    r_occultor=-1.d0
    r_opt_thick=1.d0
    dat_resolution=.false.
    output_tau=.false.
    output_absorption_fraction=.false.
 
    allocate(flux_scheme(nlevelshi),typepred1(nlevelshi),flux_method(nlevelshi))
    allocate(limiter(nlevelshi),gradient_limiter(nlevelshi))
    do level=1,nlevelshi
       flux_scheme(level) = 'tvdlf'
       typepred1(level)   = 0
       limiter(level)     = 'minmod'
       gradient_limiter(level) = 'minmod'
    end do
 
    ppm_avisc       = 0.0d0
    ppm_rjv         = 0.0d0
    flatcd          = .false.
    flatsh          = .false.
    typesourcesplit = 'sfs'
    allocate(loglimit(nw))
    loglimit(1:nw)  = .false.
 
    allocate(typeentropy(nw))
 
    do iw=1,nw
       typeentropy(iw)='nul'      ! Entropy fix type
    end do
 
    dtdiffpar     = 0.5d0
    dtpar         = -1.d0
 
    ! problem setup defaults
    iprob    = 1
 
    ! end defaults
 
    ! Initialize Kronecker delta, and Levi-Civita tensor
    do i=1,3
       do j=1,3
          if(i==j)then
             kr(i,j)=1
          else
             kr(i,j)=0
          endif
          do k=1,3
             if(i==j.or.j==k.or.k==i)then
                lvc(i,j,k)=0
             else if(i+1==j.or.i-2==j)then
                lvc(i,j,k)=1
             else
                lvc(i,j,k)=-1
             endif
          enddo
       enddo
    enddo
 
    ! These are used to construct file and log names from multiple par files
    basename_full = ''
    basename_prev = ''
 
    do i = 1, size(par_files)
       if (mype == 0) print *, "Reading " // trim(par_files(i))
 
       ! Check whether the file exists
       inquire(file=trim(par_files(i)), exist=file_exists)
 
       if (.not. file_exists) then
          write(err_msg, *) "The parameter file " // trim(par_files(i)) // &
               " does not exist"
          call mpistop(trim(err_msg))
       end if
 
       open(unitpar, file=trim(par_files(i)), status='old')
 
       ! Try to read in the namelists. They can be absent or in a different
       ! order, since we rewind before each read.
       rewind(unitpar)
       read(unitpar, filelist, end=101)
 
101    rewind(unitpar)
       read(unitpar, savelist, end=102)
 
102    rewind(unitpar)
       read(unitpar, stoplist, end=103)
 
103    rewind(unitpar)
       read(unitpar, methodlist, end=104)
 
104    rewind(unitpar)
       read(unitpar, boundlist, end=105)
 
105    rewind(unitpar)
       read(unitpar, meshlist, end=106)
 
106    rewind(unitpar)
       read(unitpar, paramlist, end=107)
 
107    rewind(unitpar)
       read(unitpar, emissionlist, end=108)
 
108    close(unitpar)
 
       ! Append the log and file names given in the par files
       if (base_filename /= basename_prev) &
            basename_full = trim(basename_full) // trim(base_filename)
       basename_prev = base_filename
    end do
 
    base_filename = basename_full
 
    ! Check whether output directory is writable
    if(mype==0) then
      dummy_file = trim(base_filename)//"DUMMY"
      open(newunit=my_unit, file=trim(dummy_file), iostat=iostate)
      if (iostate /= 0) then
         call mpistop("Can't write to output directory (" // &
              trim(base_filename) // ")")
      else
         close(my_unit, status='delete')
      end if
    end if
 
    if(source_split_usr) any_source_split=.true.
 
    ! restart filename from command line overwrites the one in par file
    if(restart_from_file_arg /= undefined) &
      restart_from_file=restart_from_file_arg
 
    ! Root process will search snapshot
    if (mype == 0) then
      if(restart_from_file == undefined) then
        ! search file from highest index
        file_exists=.false.
        do index_latest_data = 9999, 0, -1
          if(snapshot_exists(index_latest_data)) then
            file_exists=.true.
            exit
          end if
        end do
        if(.not.file_exists) index_latest_data=-1
      else
        ! get index of the given data restarted from
        index_latest_data=get_snapshot_index(trim(restart_from_file))
      end if
    end if
    call mpi_bcast(index_latest_data, 1, mpi_integer, 0, icomm, ierrmpi)
 
    if (resume_previous_run) then
      if (index_latest_data == -1) then
        if(mype==0) write(*,*) "No snapshots found to resume from, start a new run..."
      else
        ! Set file name to restart from
        write(restart_from_file, "(a,i4.4,a)") trim(base_filename),index_latest_data, ".dat"
      end if
    end if
 
    if (restart_from_file == undefined) then
      snapshotnext = 0
      slicenext    = 0
      collapsenext = 0
      if (firstprocess) &
           call mpistop("Please restart from a snapshot when firstprocess=T")
      if (convert) &
           call mpistop('Change convert to .false. for a new run!')
    end if
 
    if (small_pressure < 0.d0) call mpistop("small_pressure should be positive.")
    if (small_density < 0.d0) call mpistop("small_density should be positive.")
    if (ghostcell_comm_batched .and. ghostcell_comm_batch_size < 1) &
         call mpistop("ghostcell_comm_batch_size should be positive.")
    ! Give priority to non-zero small temperature
    if (small_temperature>0.d0) small_pressure=small_density*small_temperature
 
    if(convert) autoconvert=.false.
 
    where (tsave_log < bigdouble) tsave(:, 1) = tsave_log
    where (tsave_dat < bigdouble) tsave(:, 2) = tsave_dat
    where (tsave_slice < bigdouble) tsave(:, 3) = tsave_slice
    where (tsave_collapsed < bigdouble) tsave(:, 4) = tsave_collapsed
    where (tsave_custom < bigdouble) tsave(:, 5) = tsave_custom
 
    if (dtsave_log < bigdouble) dtsave(1) = dtsave_log
    if (dtsave_dat < bigdouble) dtsave(2) = dtsave_dat
    if (dtsave_slice < bigdouble) dtsave(3) = dtsave_slice
    if (dtsave_collapsed < bigdouble) dtsave(4) = dtsave_collapsed
    if (dtsave_custom < bigdouble) dtsave(5) = dtsave_custom
 
    if (tsavestart_log < bigdouble) tsavestart(1) = tsavestart_log
    if (tsavestart_dat < bigdouble) tsavestart(2) = tsavestart_dat
    if (tsavestart_slice < bigdouble) tsavestart(3) = tsavestart_slice
    if (tsavestart_collapsed < bigdouble) tsavestart(4) = tsavestart_collapsed
    if (tsavestart_custom < bigdouble) tsavestart(5) = tsavestart_custom
 
    if (ditsave_log < bigdouble) ditsave(1) = ditsave_log
    if (ditsave_dat < bigdouble) ditsave(2) = ditsave_dat
    if (ditsave_slice < bigdouble) ditsave(3) = ditsave_slice
    if (ditsave_collapsed < bigdouble) ditsave(4) = ditsave_collapsed
    if (ditsave_custom < bigdouble) ditsave(5) = ditsave_custom
    ! convert hours to seconds for ending wall time
    if (wall_time_max < bigdouble) wall_time_max=wall_time_max*3600.d0
 
    if (mype == 0) then
       write(unitterm, *) ''
       write(unitterm, *) 'Output type | tsavestart |  dtsave   | ditsave | itsave(1) | tsave(1)'
       write(fmt_string, *) '(A12," | ",E9.3E2,"  | ",E9.3E2," | ",I6,"  | "'//&
            ',I6, "    | ",E9.3E2)'
    end if
 
    do ifile = 1, nfile
       if (mype == 0) write(unitterm, fmt_string) trim(output_names(ifile)), &
            tsavestart(ifile), dtsave(ifile), ditsave(ifile), itsave(1, ifile), tsave(1, ifile)
    end do
 
    if (mype == 0) write(unitterm, *) ''
 
    do islice=1,nslices
       if(slicedir(islice) > ndim) &
            write(uniterr,*)'Warning in read_par_files: ', &
            'Slice ', islice,' direction',slicedir(islice),'larger than ndim=',ndim
       if(slicedir(islice) < 1) &
            write(uniterr,*)'Warning in read_par_files: ', &
            'Slice ', islice,' direction',slicedir(islice),'too small, should be [',1,ndim,']'
    end do
 
    if(it_max==biginteger .and. time_max==bigdouble.and.mype==0) write(uniterr,*) &
         'Warning in read_par_files: it_max or time_max not given!'
 
    select case (typecourant)
    case ('maxsum')
      type_courant=type_maxsum
    case ('summax')
      type_courant=type_summax
      if (local_timestep) then
       call mpistop("Type courant summax incompatible with local_timestep")
      endif
    case ('minimum')
      type_courant=type_minimum
      if (local_timestep) then
       call mpistop("Type courant minimum incompatible with local_timestep")
      endif
    case default
       write(unitterm,*)'Unknown typecourant=',typecourant
       call mpistop("Error from read_par_files: no such typecourant!")
    end select
 
 
    do level=1,nlevelshi
       select case (flux_scheme(level))
       case ('hll')
          flux_method(level)=fs_hll
       case ('hllc')
          flux_method(level)=fs_hllc
       case ('hlld')
          flux_method(level)=fs_hlld
       case ('hllcd')
          flux_method(level)=fs_hllcd
       case ('tvdlf')
          flux_method(level)=fs_tvdlf
       case ('tvdmu')
          flux_method(level)=fs_tvdmu
       case ('tvd')
          flux_method(level)=fs_tvd
       case ('cd')
          flux_method(level)=fs_cd
       case ('cd4')
          flux_method(level)=fs_cd4
       case ('fd')
          flux_method(level)=fs_fd
       case ('source')
          flux_method(level)=fs_source
       case ('nul','null')
          flux_method(level)=fs_nul
       case default
          call mpistop("unkown or bad flux scheme")
       end select
       if(flux_scheme(level)=='tvd'.and.time_stepper/='onestep') &
            call mpistop(" tvd is onestep method, reset time_stepper='onestep'")
       if(flux_scheme(level)=='tvd')then
          if(mype==0.and.(.not.dimsplit)) write(unitterm,*) &
               'Warning: setting dimsplit=T for tvd, as used for level=',level
          dimsplit=.true.
       endif
       if(flux_scheme(level)=='hlld'.and.physics_type/='mhd' .and. physics_type/='twofl') &
          call mpistop("Cannot use hlld flux if not using MHD or 2FL only charges physics!")
 
       if(flux_scheme(level)=='hlld'.and. .not.phys_total_energy) &
          call mpistop("Cannot use hlld flux if not solving total energy equation!")
 
       if(flux_scheme(level)=='hllc'.and.physics_type=='mf') &
          call mpistop("Cannot use hllc flux if using magnetofriction physics!")
 
       if(flux_scheme(level)=='tvd'.and.physics_type=='mf') &
          call mpistop("Cannot use tvd flux if using magnetofriction physics!")
 
       if(flux_scheme(level)=='tvdmu'.and.physics_type=='mf') &
          call mpistop("Cannot use tvdmu flux if using magnetofriction physics!")
 
       if (typepred1(level)==0) then
          select case (flux_scheme(level))
          case ('cd')
             typepred1(level)=fs_cd
          case ('cd4')
             typepred1(level)=fs_cd4
          case ('fd')
             typepred1(level)=fs_fd
          case ('tvdlf','tvdmu')
             typepred1(level)=fs_hancock
          case ('hll')
             typepred1(level)=fs_hll
          case ('hllc')
             typepred1(level)=fs_hllc
          case ('hllcd')
             typepred1(level)=fs_hllcd
          case ('hlld')
             typepred1(level)=fs_hlld
          case ('nul','source','tvd')
             typepred1(level)=fs_nul
          case default
             call mpistop("No default predictor for this full step")
          end select
       end if
    end do
 
    ! finite difference scheme fd need global maximal speed
    if(any(flux_scheme=='fd')) need_global_cmax=.true.
 
    ! initialize type_curl
     select case (typecurl)
     case ("central")
        type_curl=central
     case ("Gaussbased")
        type_curl=gaussbased
     case ("Stokesbased")
        type_curl=stokesbased
     case default
        write(unitterm,*) "typecurl=",typecurl
        call mpistop("unkown type of curl operator in read_par_files")
     end select
 
    ! initialize types of time stepper and time integrator
    select case (time_stepper)
    case ("onestep")
       t_stepper=onestep
       nstep=1
       if (time_integrator=='default') then
          time_integrator="Forward_Euler"
       end if
       select case (time_integrator)
       case ("Forward_Euler")
          t_integrator=forward_euler
       case ("IMEX_Euler")
          t_integrator=imex_euler
       case ("IMEX_SP")
          t_integrator=imex_sp
       case default
          write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
          call mpistop("unkown onestep time_integrator in read_par_files")
       end select
       use_imex_scheme=(t_integrator==imex_euler.or.t_integrator==imex_sp)
    case ("twostep")
       t_stepper=twostep
       nstep=2
       if (time_integrator=='default') then
          time_integrator="Predictor_Corrector"
       endif
       select case (time_integrator)
       case ("Predictor_Corrector")
          t_integrator=predictor_corrector
       case ("RK2_alfa")
          t_integrator=rk2_alf
       case ("ssprk2")
          t_integrator=ssprk2
       case ("IMEX_Midpoint")
          t_integrator=imex_midpoint
       case ("IMEX_Trapezoidal")
          t_integrator=imex_trapezoidal
       case ("IMEX_222")
          t_integrator=imex_222
       case default
          write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
          call mpistop("unkown twostep time_integrator in read_par_files")
       end select
       use_imex_scheme=(t_integrator==imex_midpoint.or.t_integrator==imex_trapezoidal&
            .or.t_integrator==imex_222)
       if (t_integrator==rk2_alf) then
          if(rk2_alfa<smalldouble.or.rk2_alfa>one)call mpistop("set rk2_alfa within [0,1]")
          rk_a21=rk2_alfa 
          rk_b2=1.0d0/(2.0d0*rk2_alfa)
          rk_b1=1.0d0-rk_b2
       endif
    case ("threestep")
       t_stepper=threestep
       nstep=3
       if (time_integrator=='default') then
          time_integrator='ssprk3'
       endif
       select case (time_integrator)
       case ("ssprk3")
          t_integrator=ssprk3
       case ("RK3_BT")
          t_integrator=rk3_bt
       case ("IMEX_ARS3")
          t_integrator=imex_ars3
       case ("IMEX_232")
          t_integrator=imex_232
       case ("IMEX_CB3a")
          t_integrator=imex_cb3a
       case default
          write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
          call mpistop("unkown threestep time_integrator in read_par_files")
       end select
       if(t_integrator==rk3_bt) then
           select case(rk3_switch)
             case(1) 
              ! we code up Ralston 3rd order here
              rk3_a21=1.0d0/2.0d0
              rk3_a31=0.0d0
              rk3_a32=3.0d0/4.0d0
              rk3_b1=2.0d0/9.0d0
              rk3_b2=1.0d0/3.0d0
             case(2) 
              ! we code up RK-Wray 3rd order here
              rk3_a21=8.0d0/15.0d0
              rk3_a31=1.0d0/4.0d0
              rk3_a32=5.0d0/12.0d0
              rk3_b1=1.0d0/4.0d0
              rk3_b2=0.0d0
             case(3) 
              ! we code up Heun 3rd order here
              rk3_a21=1.0d0/3.0d0
              rk3_a31=0.0d0
              rk3_a32=2.0d0/3.0d0
              rk3_b1=1.0d0/4.0d0
              rk3_b2=0.0d0
             case(4) 
              ! we code up Nystrom 3rd order here
              rk3_a21=2.0d0/3.0d0
              rk3_a31=0.0d0
              rk3_a32=2.0d0/3.0d0
              rk3_b1=1.0d0/4.0d0
              rk3_b2=3.0d0/8.0d0
             case default
                call mpistop("Unknown rk3_switch")
            end select
           ! the rest is fixed from above
           rk3_b3=1.0d0-rk3_b1-rk3_b2
           rk3_c2=rk3_a21
           rk3_c3=rk3_a31+rk3_a32
       endif
       if(t_integrator==ssprk3) then
         select case(ssprk_order)
             case(3) ! this is SSPRK(3,3) Gottlieb-Shu
                rk_beta11=1.0d0
                rk_beta22=1.0d0/4.0d0
                rk_beta33=2.0d0/3.0d0
                rk_alfa21=3.0d0/4.0d0
                rk_alfa31=1.0d0/3.0d0
                rk_c2=1.0d0
                rk_c3=1.0d0/2.0d0
             case(2) ! this is SSP(3,2)
                rk_beta11=1.0d0/2.0d0
                rk_beta22=1.0d0/2.0d0
                rk_beta33=1.0d0/3.0d0
                rk_alfa21=0.0d0
                rk_alfa31=1.0d0/3.0d0
                rk_c2=1.0d0/2.0d0
                rk_c3=1.0d0
             case default
                call mpistop("Unknown ssprk3_order")
         end select
         rk_alfa22=1.0d0-rk_alfa21
         rk_alfa33=1.0d0-rk_alfa31
       endif
       if(t_integrator==imex_ars3) then
          ars_gamma=(3.0d0+dsqrt(3.0d0))/6.0d0
       endif
       if(t_integrator==imex_232) then
           select case(imex_switch)
             case(1) ! this is IMEX_ARK(232)
              im_delta=1.0d0-1.0d0/dsqrt(2.0d0)
              im_nu=(3.0d0+2.0d0*dsqrt(2.0d0))/6.0d0
              imex_a21=2.0d0*im_delta
              imex_a31=1.0d0-im_nu
              imex_a32=im_nu
              imex_b1=1.0d0/(2.0d0*dsqrt(2.0d0))
              imex_b2=1.0d0/(2.0d0*dsqrt(2.0d0))
              imex_ha21=im_delta
              imex_ha22=im_delta
             case(2) ! this is IMEX_SSP(232)
              ! doi 10.1002/2017MS001065 Rokhzadi et al
              imex_a21=0.711664700366941d0
              imex_a31=0.077338168947683d0
              imex_a32=0.917273367886007d0
              imex_b1=0.398930808264688d0
              imex_b2=0.345755244189623d0
              imex_ha21=0.353842865099275d0
              imex_ha22=0.353842865099275d0
             case default
                call mpistop("Unknown imex_siwtch")
           end select
           imex_c2=imex_a21
           imex_c3=imex_a31+imex_a32
           imex_b3=1.0d0-imex_b1-imex_b2
       endif
       if(t_integrator==imex_cb3a) then
          imex_c2   = 0.8925502329346865
          imex_a22  = imex_c2
          imex_ha21 = imex_c2
          imex_c3   = imex_c2 / (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0)
          imex_ha32 = imex_c3
          imex_b2   = (3.0d0*imex_c2 - 1.0d0) / (6.0d0*imex_c2**2)
          imex_b3   = (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0) / (6.0d0*imex_c2**2)
          imex_a33  = (1.0d0/6.0d0 - imex_b2*imex_c2**2 - imex_b3*imex_c2*imex_c3) / (imex_b3*(imex_c3-imex_c2))
          imex_a32  = imex_c3 - imex_a33
          ! if (mype == 0) then
          !    write(*,*) "================================="
          !    write(*,*) "Asserting the order conditions..."
          !    ! First order convergence: OK
          !    ! Second order convergence
          !    write(*,*) -1.0d0/2.0d0 + imex_b2*imex_c2 + imex_b3*imex_c3
          !    ! Third order convergence
          !    write(*,*) -1.0d0/3.0d0 + imex_b2*imex_c2**2 + imex_b3*imex_c3**2
          !    write(*,*) -1.0d0/6.0d0 + imex_b3*imex_ha32*imex_c2
          !    write(*,*) -1.0d0/6.0d0 + imex_b2*imex_a22*imex_c2 + imex_b3*imex_a32*imex_c2 + imex_b3*imex_a33*imex_c3
          !    write(*,*) "================================="
          ! end if
       end if
       use_imex_scheme=(t_integrator==imex_ars3.or.t_integrator==imex_232.or.t_integrator==imex_cb3a)
    case ("fourstep")
       t_stepper=fourstep
       nstep=4
       if (time_integrator=='default') then
          time_integrator="ssprk4"
       endif
       select case (time_integrator)
       case ("ssprk4")
          t_integrator=ssprk4
       case ("rk4")
          t_integrator=rk4
       case default
          write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
          call mpistop("unkown fourstep time_integrator in read_par_files")
       end select
       if(t_integrator==ssprk4) then
         select case(ssprk_order)
             case(3) ! this is SSPRK(4,3) Spireti-Ruuth
                rk_beta11=1.0d0/2.0d0
                rk_beta22=1.0d0/2.0d0
                rk_beta33=1.0d0/6.0d0
                rk_beta44=1.0d0/2.0d0
                rk_alfa21=0.0d0
                rk_alfa31=2.0d0/3.0d0
                rk_alfa41=0.0d0
                rk_c2=1.0d0/2.0d0
                rk_c3=1.0d0
                rk_c4=1.0d0/2.0d0
             case(2) ! this is SSP(4,2)
                rk_beta11=1.0d0/3.0d0
                rk_beta22=1.0d0/3.0d0
                rk_beta33=1.0d0/3.0d0
                rk_beta44=1.0d0/4.0d0
                rk_alfa21=0.0d0
                rk_alfa31=0.0d0
                rk_alfa41=1.0d0/4.0d0
                rk_c2=1.0d0/3.0d0
                rk_c3=2.0d0/3.0d0
                rk_c4=1.0d0
             case default
                call mpistop("Unknown ssprk_order")
         end select
         rk_alfa22=1.0d0-rk_alfa21
         rk_alfa33=1.0d0-rk_alfa31
         rk_alfa44=1.0d0-rk_alfa41
       endif
    case ("fivestep")
       t_stepper=fivestep
       nstep=5
       if (time_integrator=='default') then
          time_integrator="ssprk5"
       end if
       select case (time_integrator)
       case ("ssprk5")
          t_integrator=ssprk5
       case default
          write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
          call mpistop("unkown fivestep time_integrator in read_par_files")
       end select
       if(t_integrator==ssprk5) then
         select case(ssprk_order)
           ! we use ssprk_order to intercompare the different coefficient choices 
           case(3) ! From Gottlieb 2005
            rk_beta11=0.391752226571890d0
            rk_beta22=0.368410593050371d0
            rk_beta33=0.251891774271694d0
            rk_beta44=0.544974750228521d0
            rk_beta54=0.063692468666290d0
            rk_beta55=0.226007483236906d0
            rk_alfa21=0.444370493651235d0
            rk_alfa31=0.620101851488403d0
            rk_alfa41=0.178079954393132d0
            rk_alfa53=0.517231671970585d0
            rk_alfa54=0.096059710526147d0
 
            rk_alfa22=0.555629506348765d0
            rk_alfa33=0.379898148511597d0
            rk_alfa44=0.821920045606868d0
            rk_alfa55=0.386708617503269d0
            rk_alfa22=1.0d0-rk_alfa21
            rk_alfa33=1.0d0-rk_alfa31
            rk_alfa44=1.0d0-rk_alfa41
            rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
            rk_c2=rk_beta11
            rk_c3=rk_alfa22*rk_c2+rk_beta22
            rk_c4=rk_alfa33*rk_c3+rk_beta33
            rk_c5=rk_alfa44*rk_c4+rk_beta44
           case(2) ! From Spireti-Ruuth
            rk_beta11=0.39175222700392d0
            rk_beta22=0.36841059262959d0
            rk_beta33=0.25189177424738d0
            rk_beta44=0.54497475021237d0
            !rk_beta54=0.06369246925946d0
            !rk_beta55=0.22600748319395d0
            rk_alfa21=0.44437049406734d0
            rk_alfa31=0.62010185138540d0
            rk_alfa41=0.17807995410773d0
            rk_alfa53=0.51723167208978d0
            !rk_alfa54=0.09605971145044d0
 
            rk_alfa22=1.0d0-rk_alfa21
            rk_alfa33=1.0d0-rk_alfa31
            rk_alfa44=1.0d0-rk_alfa41
 
            rka51=0.00683325884039d0
            rka54=0.12759831133288d0
            rkb54=0.08460416338212d0
            rk_beta54=rkb54-rk_beta44*rka51/rk_alfa41
            rk_alfa54=rka54-rk_alfa44*rka51/rk_alfa41
 
            rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
            rk_c2=rk_beta11
            rk_c3=rk_alfa22*rk_c2+rk_beta22
            rk_c4=rk_alfa33*rk_c3+rk_beta33
            rk_c5=rk_alfa44*rk_c4+rk_beta44
            rk_beta55=1.0d0-rk_beta54-rk_alfa53*rk_c3-rk_alfa54*rk_c4-rk_alfa55*rk_c5
           case default
                call mpistop("Unknown ssprk_order")
         end select
         ! the following combinations must be unity
         !print *,rk_beta55+rk_beta54+rk_alfa53*rk_c3+rk_alfa54*rk_c4+rk_alfa55*rk_c5
         !print *,rk_alfa22+rk_alfa21
         !print *,rk_alfa33+rk_alfa31
         !print *,rk_alfa44+rk_alfa41
         !print *,rk_alfa55+rk_alfa53+rk_alfa54
       endif
       use_imex_scheme=.false.
    case default
       call mpistop("Unknown time_stepper in read_par_files")
    end select
 
    do i = 1, ndim
       select case (stretch_dim(i))
       case (undefined, 'none')
          stretch_type(i) = stretch_none
          stretched_dim(i) = .false.
       case ('uni','uniform')
          stretch_type(i) = stretch_uni
          stretched_dim(i) = .true.
       case ('symm', 'symmetric')
          stretch_type(i) = stretch_symm
          stretched_dim(i) = .true.
       case default
          stretch_type(i) = stretch_none
          stretched_dim(i) = .false.
          if (mype == 0) print *, 'Got stretch_type = ', stretch_type(i)
          call mpistop('Unknown stretch type')
       end select
    end do
 
    ! Harmonize the parameters for dimensional splitting and source splitting
    if(typedimsplit   =='default'.and.     dimsplit)   typedimsplit='xyyx'
    if(typedimsplit   =='default'.and..not.dimsplit)   typedimsplit='unsplit'
    dimsplit   = typedimsplit   /='unsplit'
 
    ! initialize types of split-source addition
    select case (typesourcesplit)
    case ('sfs')
      sourcesplit=sourcesplit_sfs
    case ('sf')
      sourcesplit=sourcesplit_sf
    case ('ssf')
      sourcesplit=sourcesplit_ssf
    case ('ssfss')
      sourcesplit=sourcesplit_ssfss
    case default
       write(unitterm,*)'No such typesourcesplit=',typesourcesplit
       call mpistop("Error: Unknown typesourcesplit!")
    end select
 
    if(coordinate==-1) then
      coordinate=cartesian
      if(mype==0) then
        write(*,*) 'Warning: coordinate system is not specified!'
        write(*,*) 'call set_coordinate_system in usr_init in mod_usr.t'
        write(*,*) 'Now use Cartesian coordinate'
      end if
    end if
 
    if(coordinate==cartesian) then
      slab=.true.
      slab_uniform=.true.
      if(any(stretched_dim)) then
        coordinate=cartesian_stretched
        slab_uniform=.false.
      end if
    else
      slab=.false.
      slab_uniform=.false.
    end if
 
    if(coordinate==spherical) then
      if(dimsplit) then
        if(mype==0)print *,'Warning: spherical symmetry needs dimsplit=F, resetting'
        dimsplit=.false.
      end if
    end if
 
    if (ndim==1) dimsplit=.false.
 
    ! type limiter of prolongation
    select case(typeprolonglimit)
    case('unlimit')
      ! unlimited
      prolong_limiter=1
    case('minmod')
      prolong_limiter=2
    case('woodward')
      prolong_limiter=3
    case('koren')
      prolong_limiter=4
    case default
      prolong_limiter=0
    end select
 
    ! Type limiter is of integer type for performance
    allocate(type_limiter(nlevelshi))
    allocate(type_gradient_limiter(nlevelshi))
 
    do level=1,nlevelshi
       type_limiter(level) = limiter_type(limiter(level))
       type_gradient_limiter(level) = limiter_type(gradient_limiter(level))
    end do
 
    if (any(limiter(1:nlevelshi)== 'ppm')&
         .and.(flatsh.and.physics_type=='rho')) then
       call mpistop(" PPM with flatsh=.true. can not be used with physics_type='rho'!")
    end if
 
    ! Copy boundary conditions to typeboundary, which is used internally
    {
    do iw=1,nwfluxbc
      select case(typeboundary_min^d(iw))
      case("special")
        typeboundary(iw,2*^d-1)=bc_special
      case("cont")
        typeboundary(iw,2*^d-1)=bc_cont
      case("symm")
        typeboundary(iw,2*^d-1)=bc_symm
      case("asymm")
        typeboundary(iw,2*^d-1)=bc_asymm
      case("periodic")
        typeboundary(iw,2*^d-1)=bc_periodic
      case("aperiodic")
        typeboundary(iw,2*^d-1)=bc_aperiodic
      case("noinflow")
        typeboundary(iw,2*^d-1)=bc_noinflow
      case("pole")
        typeboundary(iw,2*^d-1)=12
      case("bc_data")
        typeboundary(iw,2*^d-1)=bc_data
      case("bc_icarus")
        typeboundary(iw,2*^d-1)=bc_icarus
      case("character")
        typeboundary(iw,2*^d-1)=bc_character
      case default
         write (unitterm,*) "Undefined boundarytype found in read_par_files", &
           typeboundary_min^d(iw),"for variable iw=",iw," and side iB=",2*^d-1
      end select
    end do
    do iw=1,nwfluxbc
      select case(typeboundary_max^d(iw))
      case("special")
        typeboundary(iw,2*^d)=bc_special
      case("cont")
        typeboundary(iw,2*^d)=bc_cont
      case("symm")
        typeboundary(iw,2*^d)=bc_symm
      case("asymm")
        typeboundary(iw,2*^d)=bc_asymm
      case("periodic")
        typeboundary(iw,2*^d)=bc_periodic
      case("aperiodic")
        typeboundary(iw,2*^d)=bc_aperiodic
      case("noinflow")
        typeboundary(iw,2*^d)=bc_noinflow
      case("pole")
        typeboundary(iw,2*^d)=12
      case("bc_data")
        typeboundary(iw,2*^d)=bc_data
      case("bc_icarus")
        typeboundary(iw,2*^d-1)=bc_icarus
      case("bc_character")
        typeboundary(iw,2*^d)=bc_character
      case default
         write (unitterm,*) "Undefined boundarytype found in read_par_files", &
           typeboundary_max^d(iw),"for variable iw=",iw," and side iB=",2*^d
      end select
    end do
    }
 
    ! psi, tracers take the same boundary type as the first variable
    if (nwfluxbc<nwflux) then
      do iw=nwfluxbc+1,nwflux
        typeboundary(iw,:) = typeboundary(1, :)
      end do
    end if
    ! auxiliary variables take the same boundary type as the first variable
    if (nwaux>0) then
      do iw=nwflux+1, nwflux+nwaux
        typeboundary(iw,:) = typeboundary(1, :)
      end do
    end if
 
    if (any(typeboundary == 0)) then
      call mpistop("Not all boundary conditions have been defined")
    end if
 
    do idim=1,ndim
       periodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_periodic))
       aperiodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_aperiodic))
       if (periodb(idim).or.aperiodb(idim)) then
          do iw=1,nwflux
             if (typeboundary(iw,2*idim-1) .ne. typeboundary(iw,2*idim)) &
                  call mpistop("Wrong counterpart in periodic boundary")
 
             if (typeboundary(iw,2*idim-1) /= bc_periodic .and. &
                  typeboundary(iw,2*idim-1) /= bc_aperiodic) then
               call mpistop("Each dimension should either have all "//&
                    "or no variables periodic, some can be aperiodic")
             end if
          end do
       end if
    end do
    {^nooned
    do idim=1,ndim
      if(any(typeboundary(:,2*idim-1)==12)) then
        if(any(typeboundary(:,2*idim-1)/=12)) typeboundary(:,2*idim-1)=12
        select case(physics_type)
        case ('rho','ard','nonlinear','ffhd')
           ! all symmetric at pole
           typeboundary(:,2*idim-1)=bc_symm
           if(mype==0) print *,'symmetric minimal pole'
        case ('hd','rhd','srhd','mhd','rmhd')
           typeboundary(:,2*idim-1)=bc_symm
           ! here we assume the ordering of variables is fixed to rho-mom-[e]-B
           if(phys_energy) then
            windex=2
           else
            windex=1
           end if
           select case(coordinate)
           case(cylindrical)
            ! across r=0 both the radial and azimuthal components flip sign
            typeboundary(r_+1,2*idim-1)=bc_asymm
            typeboundary(phi_+1,2*idim-1)=bc_asymm
            if(physics_type=='mhd'.or.physics_type=='rmhd') then
              typeboundary(ndir+windex+r_,2*idim-1)=bc_asymm
              typeboundary(ndir+windex+phi_,2*idim-1)=bc_asymm
            end if
           case(spherical)
            typeboundary(3:ndir+1,2*idim-1)=bc_asymm
            if(physics_type=='mhd'.or.physics_type=='rmhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim-1)=bc_asymm
           case default
            call mpistop('Pole is in cylindrical, polar, spherical coordinates!')
           end select
        case ('twofl')
           call mpistop('Pole treatment for twofl not implemented yet')
        case ('mf')
           if (ndim /= 3 .or. coordinate /= spherical .or. idim /= 2) &
             call mpistop("MF pole requires a 3D spherical theta boundary")
           {^ifthreed
           ! Spherical angles are still in turns at this point in the parser.
           if (.not. periodb(3) .or. abs(xprobmax3-xprobmin3-one) >= smalldouble) &
             call mpistop("MF pole requires full periodic longitude")
           if (abs(xprobmin2) >= smalldouble) &
             call mpistop("MF pole requires a true theta pole")
           ! Radial components and scalars (including GLM psi) are even.
           typeboundary(:,2*idim-1)=bc_symm
           typeboundary(iw_mom(2:3),2*idim-1)=bc_asymm
           typeboundary(iw_mag(2:3),2*idim-1)=bc_asymm
           }
        case default
           call mpistop('unknown physics type for setting minimal pole boundary treatment')
        end select
      end if
      if(any(typeboundary(:,2*idim)==12)) then
        if(any(typeboundary(:,2*idim)/=12)) typeboundary(:,2*idim)=12
        select case(physics_type)
        case ('rho','ard','nonlinear','ffhd')
           ! all symmetric at pole
           typeboundary(:,2*idim)=bc_symm
           if(mype==0) print *,'symmetric maximal pole'
        case ('hd','rhd','srhd','mhd','rmhd')
           typeboundary(:,2*idim)=bc_symm
           ! here we assume the ordering of variables is fixed to rho-mom-[e]-B
           if(phys_energy) then
            windex=2
           else
            windex=1
           end if
           select case(coordinate)
           case(cylindrical)
            ! across r=0 both the radial and azimuthal components flip sign
            typeboundary(r_+1,2*idim)=bc_asymm
            typeboundary(phi_+1,2*idim)=bc_asymm
            if(physics_type=='mhd'.or.physics_type=='rmhd') then
              typeboundary(ndir+windex+r_,2*idim)=bc_asymm
              typeboundary(ndir+windex+phi_,2*idim)=bc_asymm
            end if
           case(spherical)
            typeboundary(3:ndir+1,2*idim)=bc_asymm
            if(physics_type=='mhd'.or.physics_type=='rmhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim)=bc_asymm
           case default
            call mpistop('Pole is in cylindrical, polar, spherical coordinates!')
           end select
        case ('twofl')
           call mpistop('Pole treatment for twofl not implemented yet')
        case ('mf')
           if (ndim /= 3 .or. coordinate /= spherical .or. idim /= 2) &
             call mpistop("MF pole requires a 3D spherical theta boundary")
           {^ifthreed
           ! Spherical angles are still in turns at this point in the parser.
           if (.not. periodb(3) .or. abs(xprobmax3-xprobmin3-one) >= smalldouble) &
             call mpistop("MF pole requires full periodic longitude")
           if (abs(xprobmax2-half) >= smalldouble) &
             call mpistop("MF pole requires a true theta pole")
           ! Radial components and scalars (including GLM psi) are even.
           typeboundary(:,2*idim)=bc_symm
           typeboundary(iw_mom(2:3),2*idim)=bc_asymm
           typeboundary(iw_mag(2:3),2*idim)=bc_asymm
           }
        case default
           call mpistop('unknown physics type for setting maximal pole boundary treatment')
        end select
      end if
    end do
    }
 
    if(.not.phys_energy) then
      flatcd=.false.
      flatsh=.false.
    end if
 
    if(any(limiter(1:nlevelshi)=='mp5')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='weno5')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='weno5nm')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='wenoz5')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='wenoz5nm')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='wenozp5')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='wenozp5nm')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='teno5ad')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='weno5cu6')) then
      nghostcells=max(nghostcells,3)
    end if
 
    if(any(limiter(1:nlevelshi)=='ppm')) then
      if(flatsh .or. flatcd) then
        nghostcells=max(nghostcells,4)
      else
        nghostcells=max(nghostcells,3)
      end if
    end if
 
    if(any(limiter(1:nlevelshi)=='weno7')) then
      nghostcells=max(nghostcells,4)
    end if
 
    if(any(limiter(1:nlevelshi)=='mpweno7')) then
      nghostcells=max(nghostcells,4)
    end if
 
    ! If a wider stencil is used, extend the number of ghost cells
    nghostcells = nghostcells + phys_wider_stencil
 
    ! prolongation in AMR for constrained transport MHD needs even number ghosts
    if(stagger_grid .and. refine_max_level>1 .and. mod(nghostcells,2)/=0) then
      nghostcells=nghostcells+1
    end if
 
      select case (coordinate)
         {^nooned
      case (spherical)
         xprob^lim^de=xprob^lim^de*two*dpi;
         \}
      case (cylindrical)
         {
         if (^d==phi_) then
            xprob^lim^d=xprob^lim^d*two*dpi;
         end if
         \}
      end select
 
    ! full block size including ghostcells
    {ixghi^d = block_nx^d + 2*nghostcells\}
    {ixgshi^d = ixghi^d\}
 
    nx_vec = [{domain_nx^d|, }]
    block_nx_vec = [{block_nx^d|, }]
 
    if (any(nx_vec < 4) .or. any(mod(nx_vec, 2) == 1)) &
         call mpistop('Grid size (domain_nx^D) has to be even and >= 4')
 
    if (any(block_nx_vec < 4) .or. any(mod(block_nx_vec, 2) == 1)) &
         call mpistop('Block size (block_nx^D) has to be even and >= 4')
 
    { if(mod(domain_nx^d,block_nx^d)/=0) &
       call mpistop('Grid (domain_nx^D) and block (block_nx^D) must be consistent') \}
 
    if(refine_max_level>nlevelshi.or.refine_max_level<1)then
       write(unitterm,*)'Error: refine_max_level',refine_max_level,'>nlevelshi ',nlevelshi
       call mpistop("Reset nlevelshi and recompile!")
    endif
 
    if (any(stretched_dim)) then
       allocate(qstretch(0:nlevelshi,1:ndim),dxfirst(0:nlevelshi,1:ndim),&
             dxfirst_1mq(0:nlevelshi,1:ndim),dxmid(0:nlevelshi,1:ndim))
       allocate(nstretchedblocks(1:nlevelshi,1:ndim))
       qstretch(0:nlevelshi,1:ndim)=0.0d0
       dxfirst(0:nlevelshi,1:ndim)=0.0d0
       nstretchedblocks(1:nlevelshi,1:ndim)=0
       {if (stretch_type(^d) == stretch_uni) then
           ! first some sanity checks
           if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble) then
             if(mype==0) then
               write(*,*) 'stretched grid needs finite qstretch_baselevel>1'
               write(*,*) 'will try default value for qstretch_baselevel in dimension', ^d
             endif
             if(xprobmin^d>smalldouble)then
               qstretch_baselevel(^d)=(xprobmax^d/xprobmin^d)**(1.d0/dble(domain_nx^d))
             else
               call mpistop("can not set qstretch_baselevel automatically")
             endif
           endif
           if(mod(block_nx^d,2)==1) &
             call mpistop("stretched grid needs even block size block_nxD")
           if(mod(domain_nx^d/block_nx^d,2)/=0) &
             call mpistop("number level 1 blocks in D must be even")
           qstretch(1,^d)=qstretch_baselevel(^d)
           dxfirst(1,^d)=(xprobmax^d-xprobmin^d) &
                *(1.0d0-qstretch(1,^d))/(1.0d0-qstretch(1,^d)**domain_nx^d)
           qstretch(0,^d)=qstretch(1,^d)**2
           dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
           if(refine_max_level>1)then
              do ilev=2,refine_max_level
                 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
                 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
                       /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
              enddo
           endif
        endif \}
        if(mype==0) then
           {if(stretch_type(^d) == stretch_uni) then
              write(*,*) 'Stretched dimension ', ^d
              write(*,*) 'Using stretched grid with qs=',qstretch(0:refine_max_level,^d)
              write(*,*) '        and first cell sizes=',dxfirst(0:refine_max_level,^d)
           endif\}
        end if
       {if(stretch_type(^d) == stretch_symm) then
           if(mype==0) then
               write(*,*) 'will apply symmetric stretch in dimension', ^d
           endif
           if(mod(block_nx^d,2)==1) &
             call mpistop("stretched grid needs even block size block_nxD")
           ! checks on the input variable nstretchedblocks_baselevel
           if(nstretchedblocks_baselevel(^d)==0) &
             call mpistop("need finite even number of stretched blocks at baselevel")
           if(mod(nstretchedblocks_baselevel(^d),2)==1) &
             call mpistop("need even number of stretched blocks at baselevel")
           if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble) &
             call mpistop('stretched grid needs finite qstretch_baselevel>1')
           ! compute stretched part to ensure uniform center
           ipower=(nstretchedblocks_baselevel(^d)/2)*block_nx^d
           if(nstretchedblocks_baselevel(^d)==domain_nx^d/block_nx^d)then
              xstretch^d=0.5d0*(xprobmax^d-xprobmin^d)
           else
              xstretch^d=(xprobmax^d-xprobmin^d) &
                /(2.0d0+dble(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d) &
               *(1.0d0-qstretch_baselevel(^d))/(1.0d0-qstretch_baselevel(^d)**ipower))
           endif
           if(xstretch^d>(xprobmax^d-xprobmin^d)*0.5d0) &
             call mpistop(" stretched grid part should not exceed full domain")
           dxfirst(1,^d)=xstretch^d*(1.0d0-qstretch_baselevel(^d)) &
                                   /(1.0d0-qstretch_baselevel(^d)**ipower)
           nstretchedblocks(1,^d)=nstretchedblocks_baselevel(^d)
           qstretch(1,^d)=qstretch_baselevel(^d)
           qstretch(0,^d)=qstretch(1,^d)**2
           dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
           dxmid(1,^d)=dxfirst(1,^d)
           dxmid(0,^d)=dxfirst(1,^d)*2.0d0
           if(refine_max_level>1)then
              do ilev=2,refine_max_level
                 nstretchedblocks(ilev,^d)=2*nstretchedblocks(ilev-1,^d)
                 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
                 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
                         /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
                 dxmid(ilev,^d)=dxmid(ilev-1,^d)/2.0d0
              enddo
           endif
           ! sanity check on total domain size:
           sizeuniformpart^d=dxfirst(1,^d) &
             *(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d)
           if(mype==0) then
              print *,'uniform part of size=',sizeuniformpart^d
              print *,'setting of domain is then=',2*xstretch^d+sizeuniformpart^d
              print *,'versus=',xprobmax^d-xprobmin^d
           endif
           if(dabs(xprobmax^d-xprobmin^d-2*xstretch^d-sizeuniformpart^d)>smalldouble) then
              call mpistop('mismatch in domain size!')
           endif
        endif \}
        dxfirst_1mq(0:refine_max_level,1:ndim)=dxfirst(0:refine_max_level,1:ndim) &
                              /(1.0d0-qstretch(0:refine_max_level,1:ndim))
    end if
 
    dx_vec = [{xprobmax^d-xprobmin^d|, }] / nx_vec
 
    if (mype==0) then
       write(c_ndim, '(I1)') ^nd
       write(unitterm, '(A30,' // c_ndim // '(I0," "))') &
            ' Domain size (cells): ', nx_vec
       write(unitterm, '(A30,' // c_ndim // '(E9.3," "))') &
            ' Level one dx: ', dx_vec
    end if
 
    if (any(dx_vec < smalldouble)) &
         call mpistop("Incorrect domain size (too small grid spacing)")
 
    dx(:, 1) = dx_vec
 
    if(sum(w_refine_weight(:))==0) w_refine_weight(1) = 1.d0
    if(dabs(sum(w_refine_weight(:))-1.d0)>smalldouble) then
      write(unitterm,*) "Sum of all elements in w_refine_weight be 1.d0"
      call mpistop("Reset w_refine_weight so the sum is 1.d0")
    end if
 
    select case (typeboundspeed)
    case('Einfeldt')
      boundspeed=1
    case('cmaxmean')
      boundspeed=2
    case('cmaxleftright')
      boundspeed=3
    case('pvrs')
      boundspeed=4
    case default
      call mpistop("set typeboundspeed='Einfeldt' or 'cmaxmean' or 'cmaxleftright' or 'pvrs'")
    end select
 
    if (mype==0) write(unitterm, '(A30)', advance='no') 'Refine estimation: '
 
    select case (refine_criterion)
    case (0)
       if (mype==0) write(unitterm, '(A)') "user defined"
    case (1)
       if (mype==0) write(unitterm, '(A)') "relative error"
    case (2)
       if (mype==0) write(unitterm, '(A)') "Lohner's original scheme"
    case (3)
       if (mype==0) write(unitterm, '(A)') "Lohner's scheme"
    case default
       call mpistop("Unknown error estimator, change refine_criterion")
    end select
 
    if (tfixgrid<bigdouble/2.0d0) then
       if(mype==0)print*,'Warning, at time=',tfixgrid,'the grid will be fixed'
    end if
    if (itfixgrid<biginteger/2) then
       if(mype==0)print*,'Warning, at iteration=',itfixgrid,'the grid will be fixed'
    end if
    if (ditregrid>1) then
       if(mype==0)print*,'Note, Grid is reconstructed once every',ditregrid,'iterations'
    end if
 
 
    do islice=1,nslices
       select case(slicedir(islice))
          {case(^d)
          if(slicecoord(islice)<xprobmin^d.or.slicecoord(islice)>xprobmax^d) &
               write(uniterr,*)'Warning in read_par_files: ', &
               'Slice ', islice, ' coordinate',slicecoord(islice),'out of bounds for dimension ',slicedir(islice)
          \}
       end select
    end do
 
    if (mype==0) then
       write(unitterm, '(A30,A,A)') 'restart_from_file: ', ' ', trim(restart_from_file)
       write(unitterm, '(A30,L1)') 'converting: ', convert
       write(unitterm, '(A)') ''
    endif
 
    deallocate(flux_scheme)
 
  end subroutine read_par_files
 
  !> Routine to find entries in a string
  subroutine get_fields_string(line, delims, n_max, fields, n_found, fully_read)
    !> The line from which we want to read
    character(len=*), intent(in)    :: line
    !> A string with delimiters. For example delims = " ,'"""//char(9)
    character(len=*), intent(in)    :: delims
    !> Maximum number of entries to read in
    integer, intent(in)             :: n_max
    !> Number of entries found
    integer, intent(inout)          :: n_found
    !> Fields in the strings
    character(len=*), intent(inout) :: fields(n_max)
    logical, intent(out), optional  :: fully_read
 
    integer :: ixs_start(n_max)
    integer :: ixs_end(n_max)
    integer :: ix, ix_prev
 
    ix_prev = 0
    n_found = 0
 
    do while (n_found < n_max)
       ! Find the starting point of the next entry (a non-delimiter value)
       ix = verify(line(ix_prev+1:), delims)
       if (ix == 0) exit
 
       n_found            = n_found + 1
       ixs_start(n_found) = ix_prev + ix ! This is the absolute position in 'line'
 
       ! Get the end point of the current entry (next delimiter index minus one)
       ix = scan(line(ixs_start(n_found)+1:), delims) - 1
 
       if (ix == -1) then              ! If there is no last delimiter,
          ixs_end(n_found) = len(line) ! the end of the line is the endpoint
       else
          ixs_end(n_found) = ixs_start(n_found) + ix
       end if
 
       fields(n_found) = line(ixs_start(n_found):ixs_end(n_found))
       ix_prev = ixs_end(n_found) ! We continue to search from here
    end do
 
    if (present(fully_read)) then
      ix = verify(line(ix_prev+1:), delims)
      fully_read = (ix == 0)    ! Are there only delimiters?
    end if
 
  end subroutine get_fields_string
 
  subroutine saveamrfile(ifile)
 
    use mod_usr_methods, only: usr_print_log, usr_write_analysis
    use mod_global_parameters
    use mod_particles, only: write_particles_snapshot
    use mod_slice, only: write_slice
    use mod_collapse, only: write_collapsed
    use mod_convert_files, only: generate_plotfile
    integer:: ifile
 
    select case (ifile)
    case (fileout_)
       ! Write .dat snapshot
       call write_snapshot()
 
       ! Generate formatted output (e.g., VTK)
       if(autoconvert) call generate_plotfile
 
       if(use_particles) call write_particles_snapshot()
 
       snapshotnext = snapshotnext + 1
    case (fileslice_)
       call write_slice
    case (filecollapse_)
       call write_collapsed
    case (filelog_)
       select case (typefilelog)
       case ('default')
          call printlog_default
       case ('regression_test')
          call printlog_regression_test()
       case ('special')
          if (.not. associated(usr_print_log)) then
             call mpistop("usr_print_log not defined")
          else
             call usr_print_log()
          end if
       case default
          call mpistop("Error in SaveFile: Unknown typefilelog")
       end select
    case (fileanalysis_)
       if (associated(usr_write_analysis)) then
          call usr_write_analysis()
       end if
    case default
       write(*,*) 'No save method is defined for ifile=',ifile
       call mpistop("")
    end select
 
    ! opedit: Flush stdout and stderr from time to time.
    flush(unit=unitterm)
 
  end subroutine saveamrfile
 
 
  ! Check if a snapshot exists
  logical function snapshot_exists(ix)
    use mod_global_parameters
    integer, intent(in)    :: ix !< Index of snapshot
    character(len=std_len) :: filename
 
    write(filename, "(a,i4.4,a)") trim(base_filename), ix, ".dat"
    inquire(file=trim(filename), exist=snapshot_exists)
  end function snapshot_exists
 
  integer function get_snapshot_index(filename)
    character(len=*), intent(in) :: filename
    integer                      :: i
 
    ! Try to parse index in restart_from_file string (e.g. basename0000.dat)
    i = len_trim(filename) - 7
    read(filename(i:i+3), '(I4)') get_snapshot_index
  end function get_snapshot_index
 
 
 
  !> Write header for a snapshot
  !>
  !> If you edit the header, don't forget to update: snapshot_write_header(),
  !> snapshot_read_header(), doc/fileformat.md, tools/python/dat_reader.py
  subroutine snapshot_write_header(fh, offset_tree, offset_block)
    use mod_forest
    use mod_physics
    use mod_global_parameters
    use mod_input_output_helper, only: snapshot_write_header1
    integer, intent(in)                       :: fh           !< File handle
    integer(kind=MPI_OFFSET_KIND), intent(in) :: offset_tree  !< Offset of tree info
    integer(kind=MPI_OFFSET_KIND), intent(in) :: offset_block !< Offset of block data
    call snapshot_write_header1(fh, offset_tree, offset_block, cons_wnames, nw)
  end subroutine snapshot_write_header
 
  !> Read header for a snapshot
  !>
  !> If you edit the header, don't forget to update: snapshot_write_header(),
  !> snapshot_read_header(), doc/fileformat.md, tools/python/dat_reader.py
  subroutine snapshot_read_header(fh, offset_tree, offset_block)
    use mod_forest
    use mod_global_parameters
    use mod_physics, only: physics_type
    integer, intent(in)                   :: fh           !< File handle
    integer(MPI_OFFSET_KIND), intent(out) :: offset_tree  !< Offset of tree info
    integer(MPI_OFFSET_KIND), intent(out) :: offset_block !< Offset of block data
 
    double precision                      :: rbuf(ndim)
    double precision, allocatable         :: params(:)
    integer                               :: i, version
    integer                               :: ibuf(ndim), iw
    integer                               :: er, n_par, tmp_int
    integer, dimension(MPI_STATUS_SIZE)   :: st
    logical                               :: periodic(ndim)
    character(len=name_len), allocatable  :: var_names(:), param_names(:)
    character(len=name_len)               :: phys_name, geom_name
 
    ! Version number
    call mpi_file_read(fh, version, 1, mpi_integer, st, er)
    if (all(compatible_versions /= version)) then
      call mpistop("Incompatible file version (maybe old format?)")
    end if
 
    ! offset_tree
    call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
    offset_tree = ibuf(1)
 
    ! offset_block
    call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
    offset_block = ibuf(1)
 
    ! nw
    call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
    nw_found=ibuf(1)
    if (nw /= ibuf(1)) then
      write(*,*) "nw=",nw," and nw found in restart file=",ibuf(1)
      write(*,*) "Please be aware of changes in w at restart."
      !call mpistop("currently, changing nw at restart is not allowed")
    end if
 
    ! ndir
    call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
    if (ibuf(1) /= ndir) then
      if (allow_ndir_change) then
        if (mype==0) then
          write(*,*) "WARNING: ndir in restart file = ",ibuf(1)," but current ndir = ",ndir
          write(*,*) "allow_ndir_change=T: loading anyway (block I/O is ndim-based)."
          write(*,*) "Ensure usr_transform_w maps the source vars into the right slots."
        end if
      else
        write(*,*) "ndir in restart file = ",ibuf(1)
        write(*,*) "ndir = ",ndir
        call mpistop("reset ndir to ndir in restart file (or set allow_ndir_change=T)")
      end if
    end if
 
    ! ndim
    call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
    if (ibuf(1) /= ndim) then
      write(*,*) "ndim in restart file = ",ibuf(1)
      write(*,*) "ndim = ",ndim
      call mpistop("reset ndim to ndim in restart file")
    end if
 
    ! levmax
    call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
    if (ibuf(1) > refine_max_level) then
      write(*,*) "number of levels in restart file = ",ibuf(1)
      write(*,*) "refine_max_level = ",refine_max_level
      call mpistop("refine_max_level < num. levels in restart file")
    end if
 
    ! nleafs
    call mpi_file_read(fh, nleafs, 1, mpi_integer, st, er)
 
    ! nparents
    call mpi_file_read(fh, nparents, 1, mpi_integer, st, er)
 
    ! it
    call mpi_file_read(fh, it, 1, mpi_integer, st, er)
 
    ! global time
    call mpi_file_read(fh, global_time, 1, mpi_double_precision, st, er)
 
    ! xprobmin^D
    call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
    if (maxval(abs(rbuf(1:ndim) - [ xprobmin^d ])) > 0) then
      write(*,*) "Error: xprobmin differs from restart data: ", rbuf(1:ndim)
      call mpistop("change xprobmin^D in par file")
    end if
 
    ! xprobmax^D
    call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
    if (maxval(abs(rbuf(1:ndim) - [ xprobmax^d ])) > 0) then
      write(*,*) "Error: xprobmax differs from restart data: ", rbuf(1:ndim)
      call mpistop("change xprobmax^D in par file")
    end if
 
    ! domain_nx^D
    call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
    if (any(ibuf(1:ndim) /= [ domain_nx^d ])) then
      write(*,*) "Error: mesh size differs from restart data: ", ibuf(1:ndim)
      call mpistop("change domain_nx^D in par file")
    end if
 
    ! block_nx^D
    call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
    if (any(ibuf(1:ndim) /= [ block_nx^d ])) then
      write(*,*) "Error: block size differs from restart data:", ibuf(1:ndim)
      call mpistop("change block_nx^D in par file")
    end if
 
    ! From version 5, read more info about the grid
    if (version > 4) then
      call mpi_file_read(fh, periodic, ndim, mpi_logical, st, er)
      if ({periodic(^d) .and. .not.periodb(^d) .or. .not.periodic(^d) .and. periodb(^d)| .or. }) &
           call mpistop("change in periodicity in par file")
 
      call mpi_file_read(fh, geom_name, name_len, mpi_character, st, er)
 
      if (geom_name /= geometry_name(1:name_len)) then
        if (allow_ndir_change) then
          if (mype==0) write(*,*) "WARNING: coordinates in data = ",trim(geom_name), &
               " vs current ",trim(geometry_name),"; allow_ndir_change=T (e.g. Cartesian_2D->2.5D), loading anyway."
        else
          write(*,*) "type of coordinates in data is: ", geom_name
          call mpistop("select the correct coordinates in mod_usr.t file")
        end if
      end if
 
      call mpi_file_read(fh, stagger_mark_dat, 1, mpi_logical, st, er)
      if (stagger_grid .and. .not. stagger_mark_dat .or. .not.stagger_grid.and.stagger_mark_dat) then
        write(*,*) "Warning: stagger grid flag differs from restart data:", stagger_mark_dat
        !call mpistop("change parameter to use stagger grid")
      end if
    end if
 
    ! From version 4 onwards, the later parts of the header must be present
    if (version > 3) then
      ! w_names (not used here)
      allocate(var_names(nw_found))
      do iw = 1, nw_found
        call mpi_file_read(fh, var_names(iw), name_len, mpi_character, st, er)
      end do
 
      ! Physics related information
      call mpi_file_read(fh, phys_name, name_len, mpi_character, st, er)
 
      if (phys_name /= physics_type) then
!        call mpistop("Cannot restart with a different physics type")
      end if
 
      call mpi_file_read(fh, n_par, 1, mpi_integer, st, er)
      allocate(params(n_par))
      allocate(param_names(n_par))
      call mpi_file_read(fh, params, n_par, mpi_double_precision, st, er)
      call mpi_file_read(fh, param_names, name_len * n_par, mpi_character, st, er)
 
      ! Read snapshotnext etc. for restarting
      call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
 
      ! Only set snapshotnext if the user hasn't specified it
      if (snapshotnext == -1) snapshotnext = tmp_int
 
      call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
      if (slicenext == -1) slicenext = tmp_int
 
      call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
      if (collapsenext == -1) collapsenext = tmp_int
    else
      ! Guess snapshotnext from file name if not set
      if (snapshotnext == -1) &
           snapshotnext = get_snapshot_index(trim(restart_from_file)) + 1
      ! Set slicenext and collapsenext if not set
      if (slicenext == -1) slicenext = 0
      if (collapsenext == -1) collapsenext = 0
    end if
 
    ! Still used in convert
    snapshotini = snapshotnext-1
 
  end subroutine snapshot_read_header
 
  subroutine write_snapshot
    use mod_forest
    use mod_global_parameters
    use mod_physics
    use mod_input_output_helper, only: count_ix,block_shape_io,create_output_file
    use mod_functions_forest, only: write_forest
 
    double precision, allocatable :: w_buffer(:)
    integer                       :: file_handle, igrid, Morton_no, iwrite
    integer                       :: ipe, ix_buffer(2*ndim+1), n_values
    integer                       :: ixO^L, n_ghost(2*ndim)
    integer                       :: ixOs^L,n_values_stagger
    integer                       :: iorecvstatus(MPI_STATUS_SIZE)
    integer                       :: ioastatus(MPI_STATUS_SIZE)
    integer                       :: igrecvstatus(MPI_STATUS_SIZE)
    integer                       :: istatus(MPI_STATUS_SIZE)
    integer(kind=MPI_OFFSET_KIND) :: offset_tree_info
    integer(kind=MPI_OFFSET_KIND) :: offset_block_data
    integer(kind=MPI_OFFSET_KIND) :: offset_offsets
    integer, allocatable                       :: block_ig(:, :)
    integer, allocatable                       :: block_lvl(:)
    integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset(:)
    type(tree_node), pointer      :: pnode
 
    call mpi_barrier(icomm, ierrmpi)
 
    ! Allocate send/receive buffer
    n_values = count_ix(ixg^ll) * nw
    if(stagger_grid) then
      n_values = n_values + count_ix(ixgs^ll) * nws
    end if
    allocate(w_buffer(n_values))
 
    ! Allocate arrays with information about grid blocks
    allocate(block_ig(ndim, nleafs))
    allocate(block_lvl(nleafs))
    allocate(block_offset(nleafs+1))
 
    ! master processor
    if (mype==0) then
      call create_output_file(file_handle, snapshotnext, ".dat")
 
      ! Don't know offsets yet, we will write header again later
      offset_tree_info = -1
      offset_block_data = -1
      call snapshot_write_header(file_handle, offset_tree_info, &
           offset_block_data)
 
      call mpi_file_get_position(file_handle, offset_tree_info, ierrmpi)
 
      call write_forest(file_handle)
 
      ! Collect information about the spatial index (ig^D) and refinement level
      ! of leaves
      do morton_no = morton_start(0), morton_stop(npe-1)
         igrid = sfc(1, morton_no)
         ipe = sfc(2, morton_no)
         pnode => igrid_to_node(igrid, ipe)%node
 
         block_ig(:, morton_no) = [ pnode%ig^d ]
         block_lvl(morton_no) = pnode%level
         block_offset(morton_no) = 0 ! Will be determined later
      end do
 
      call mpi_file_write(file_handle, block_lvl, size(block_lvl), &
           mpi_integer, istatus, ierrmpi)
 
      call mpi_file_write(file_handle, block_ig, size(block_ig), &
           mpi_integer, istatus, ierrmpi)
 
      ! Block offsets are currently unknown, but will be overwritten later
      call mpi_file_get_position(file_handle, offset_offsets, ierrmpi)
      call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
           mpi_offset, istatus, ierrmpi)
 
      call mpi_file_get_position(file_handle, offset_block_data, ierrmpi)
 
      ! Check whether data was written as expected
      if (offset_block_data - offset_tree_info /= &
           (nleafs + nparents) * size_logical + &
           nleafs * ((1+ndim) * size_int + 2 * size_int)) then
        if (mype == 0) then
          print *, "Warning: MPI_OFFSET type /= 8 bytes"
          print *, "This *could* cause problems when reading .dat files"
        end if
      end if
 
      block_offset(1) = offset_block_data
      iwrite = 0
    end if
 
    do morton_no=morton_start(mype), morton_stop(mype)
      igrid  = sfc_to_igrid(morton_no)
      itag   = morton_no
 
      call block_shape_io(igrid, n_ghost, ixo^l, n_values)
      if(stagger_grid) then
        w_buffer(1:n_values) = pack(ps(igrid)%w(ixo^s, 1:nw), .true.)
        {ixosmin^d = ixomin^d -1\}
        {ixosmax^d = ixomax^d \}
        n_values_stagger= count_ix(ixos^l)*nws
        w_buffer(n_values+1:n_values+n_values_stagger) = pack(ps(igrid)%ws(ixos^s, 1:nws), .true.)
        n_values=n_values+n_values_stagger
      else
        w_buffer(1:n_values) = pack(ps(igrid)%w(ixo^s, 1:nw), .true.)
      end if
      ix_buffer(1) = n_values
      ix_buffer(2:) = n_ghost
 
      if (mype /= 0) then
        call mpi_send(ix_buffer, 2*ndim+1, &
             mpi_integer, 0, itag, icomm, ierrmpi)
        call mpi_send(w_buffer, n_values, &
             mpi_double_precision, 0, itag, icomm, ierrmpi)
      else
        iwrite = iwrite+1
        call mpi_file_write(file_handle, ix_buffer(2:), &
             2*ndim, mpi_integer, istatus, ierrmpi)
        call mpi_file_write(file_handle, w_buffer, &
             n_values, mpi_double_precision, istatus, ierrmpi)
 
        ! Set offset of next block
        block_offset(iwrite+1) = block_offset(iwrite) + &
             int(n_values, mpi_offset_kind) * size_double + &
             2 * ndim * size_int
      end if
    end do
 
    ! Write data communicated from other processors
    if (mype == 0) then
      do ipe = 1, npe-1
        do morton_no=morton_start(ipe), morton_stop(ipe)
          iwrite=iwrite+1
          itag=morton_no
 
          call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, ipe, itag, icomm,&
               igrecvstatus, ierrmpi)
          n_values = ix_buffer(1)
 
          call mpi_recv(w_buffer, n_values, mpi_double_precision,&
               ipe, itag, icomm, iorecvstatus, ierrmpi)
 
          call mpi_file_write(file_handle, ix_buffer(2:), &
             2*ndim, mpi_integer, istatus, ierrmpi)
          call mpi_file_write(file_handle, w_buffer, &
             n_values, mpi_double_precision, istatus, ierrmpi)
 
          ! Set offset of next block
          block_offset(iwrite+1) = block_offset(iwrite) + &
               int(n_values, mpi_offset_kind) * size_double + &
               2 * ndim * size_int
        end do
      end do
 
      ! Write block offsets (now we know them)
      call mpi_file_seek(file_handle, offset_offsets, mpi_seek_set, ierrmpi)
      call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
           mpi_offset, istatus, ierrmpi)
 
      ! Write header again, now with correct offsets
      call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set, ierrmpi)
      call snapshot_write_header(file_handle, offset_tree_info, &
           offset_block_data)
 
      call mpi_file_close(file_handle, ierrmpi)
    end if
 
    call mpi_barrier(icomm, ierrmpi)
  end subroutine write_snapshot
 
  !> Enable the debug field dump: register n named slots. Capture fields with
  !> ps(igrid)%wdebug(...,islot) anywhere in a block loop (lazy per-block
  !> allocation at the capture site), then flush with save_wdebug.
  subroutine debug_alloc(n, names)
    use mod_global_parameters
    integer, intent(in)          :: n
    character(len=*), intent(in) :: names(n)
    integer :: i
    n_wdebug = n
    if (allocated(wdebug_names)) deallocate(wdebug_names)
    allocate(wdebug_names(n))
    do i = 1, n
      wdebug_names(i) = trim(adjustl(names(i)))
    end do
    wdebug_on = .true.
  end subroutine debug_alloc
 
  !> Flush ps(:)%wdebug to a standalone cell-centred .dat (no staggered),
  !> readable by the standard AMRVAC .dat readers. Collective; call at a
  !> barrier-safe point (end of a step), NOT inside a block loop.
  subroutine save_wdebug(suffix)
    use mod_forest
    use mod_global_parameters
    use mod_physics
    use mod_input_output_helper, only: count_ix, block_shape_io, create_output_file, snapshot_write_header1
    use mod_functions_forest, only: write_forest
    character(len=*), intent(in)  :: suffix
    logical                       :: stagger_save
    double precision, allocatable :: w_buffer(:)
    integer                       :: file_handle, igrid, Morton_no, iwrite
    integer                       :: ipe, ix_buffer(2*ndim+1), n_values
    integer                       :: ixO^L, n_ghost(2*ndim)
    integer                       :: iorecvstatus(MPI_STATUS_SIZE)
    integer                       :: igrecvstatus(MPI_STATUS_SIZE)
    integer                       :: istatus(MPI_STATUS_SIZE)
    integer(kind=MPI_OFFSET_KIND) :: offset_tree_info, offset_block_data, offset_offsets
    integer, allocatable          :: block_ig(:, :), block_lvl(:)
    integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset(:)
    type(tree_node), pointer      :: pnode
 
    if (n_wdebug <= 0) return
    call mpi_barrier(icomm, ierrmpi)
 
    n_values = count_ix(ixg^ll) * n_wdebug
    allocate(w_buffer(n_values))
    allocate(block_ig(ndim, nleafs), block_lvl(nleafs), block_offset(nleafs+1))
 
    if (mype == 0) then
      call create_output_file(file_handle, it, ".dat", trim(suffix))
      offset_tree_info = -1; offset_block_data = -1
      stagger_save = stagger_grid; stagger_grid = .false.
      call snapshot_write_header1(file_handle, offset_tree_info, offset_block_data, wdebug_names, n_wdebug)
      stagger_grid = stagger_save
      call mpi_file_get_position(file_handle, offset_tree_info, ierrmpi)
      call write_forest(file_handle)
      do morton_no = morton_start(0), morton_stop(npe-1)
        igrid = sfc(1, morton_no); ipe = sfc(2, morton_no)
        pnode => igrid_to_node(igrid, ipe)%node
        block_ig(:, morton_no) = [ pnode%ig^d ]
        block_lvl(morton_no) = pnode%level
        block_offset(morton_no) = 0
      end do
      call mpi_file_write(file_handle, block_lvl, size(block_lvl), mpi_integer, istatus, ierrmpi)
      call mpi_file_write(file_handle, block_ig, size(block_ig), mpi_integer, istatus, ierrmpi)
      call mpi_file_get_position(file_handle, offset_offsets, ierrmpi)
      call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, mpi_offset, istatus, ierrmpi)
      call mpi_file_get_position(file_handle, offset_block_data, ierrmpi)
      block_offset(1) = offset_block_data
      iwrite = 0
    end if
 
    do morton_no = morton_start(mype), morton_stop(mype)
      igrid = sfc_to_igrid(morton_no); itag = morton_no
      call block_shape_io(igrid, n_ghost, ixo^l, n_values)
      n_values = count_ix(ixo^l) * n_wdebug
      if (allocated(ps(igrid)%wdebug)) then
        w_buffer(1:n_values) = pack(ps(igrid)%wdebug(ixo^s, 1:n_wdebug), .true.)
      else
        w_buffer(1:n_values) = 0.0d0
      end if
      ix_buffer(1) = n_values
      ix_buffer(2:) = n_ghost
      if (mype /= 0) then
        call mpi_send(ix_buffer, 2*ndim+1, mpi_integer, 0, itag, icomm, ierrmpi)
        call mpi_send(w_buffer, n_values, mpi_double_precision, 0, itag, icomm, ierrmpi)
      else
        iwrite = iwrite+1
        call mpi_file_write(file_handle, ix_buffer(2:), 2*ndim, mpi_integer, istatus, ierrmpi)
        call mpi_file_write(file_handle, w_buffer, n_values, mpi_double_precision, istatus, ierrmpi)
        block_offset(iwrite+1) = block_offset(iwrite) + &
             int(n_values, mpi_offset_kind) * size_double + 2 * ndim * size_int
      end if
    end do
 
    if (mype == 0) then
      do ipe = 1, npe-1
        do morton_no = morton_start(ipe), morton_stop(ipe)
          iwrite = iwrite+1; itag = morton_no
          call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, ipe, itag, icomm, igrecvstatus, ierrmpi)
          n_values = ix_buffer(1)
          call mpi_recv(w_buffer, n_values, mpi_double_precision, ipe, itag, icomm, iorecvstatus, ierrmpi)
          call mpi_file_write(file_handle, ix_buffer(2:), 2*ndim, mpi_integer, istatus, ierrmpi)
          call mpi_file_write(file_handle, w_buffer, n_values, mpi_double_precision, istatus, ierrmpi)
          block_offset(iwrite+1) = block_offset(iwrite) + &
               int(n_values, mpi_offset_kind) * size_double + 2 * ndim * size_int
        end do
      end do
      call mpi_file_seek(file_handle, offset_offsets, mpi_seek_set, ierrmpi)
      call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, mpi_offset, istatus, ierrmpi)
      call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set, ierrmpi)
      stagger_save = stagger_grid; stagger_grid = .false.
      call snapshot_write_header1(file_handle, offset_tree_info, offset_block_data, wdebug_names, n_wdebug)
      stagger_grid = stagger_save
      call mpi_file_close(file_handle, ierrmpi)
    end if
    deallocate(w_buffer, block_ig, block_lvl, block_offset)
    call mpi_barrier(icomm, ierrmpi)
    if (mype==0) write(*,*) 'save_wdebug: wrote ', n_wdebug, ' field(s) suffix=', trim(suffix), ' it=', it
  end subroutine save_wdebug
 
 
  !> Routine to read in snapshots (.dat files). When it cannot recognize the
  !> file version, it will automatically try the 'old' reader.
  subroutine read_snapshot
    use mod_usr_methods, only: usr_transform_w
    use mod_input_output_helper, only: count_ix
    use mod_forest
    use mod_global_parameters
    use mod_amr_solution_node, only: alloc_node
    use mod_functions_forest, only: read_forest
 
    double precision :: ws(ixGs^T,1:ndim)
    double precision, allocatable :: w_buffer(:)
    double precision, dimension(:^D&,:), allocatable :: w
    integer                       :: ix_buffer(2*ndim+1), n_values, n_values_stagger
    integer                       :: ixO^L, ixOs^L
    integer                       :: file_handle, amode, igrid, Morton_no, iread
    integer                       :: istatus(MPI_STATUS_SIZE)
    integer                       :: iorecvstatus(MPI_STATUS_SIZE)
    integer                       :: ipe,inrecv,nrecv, file_version
    integer(MPI_OFFSET_KIND)      :: offset_tree_info
    integer(MPI_OFFSET_KIND)      :: offset_block_data
    logical                       :: fexist
 
    if (mype==0) then
      inquire(file=trim(restart_from_file), exist=fexist)
      if (.not.fexist) call mpistop(trim(restart_from_file)//" not found!")
 
      call mpi_file_open(mpi_comm_self,restart_from_file,mpi_mode_rdonly, &
           mpi_info_null,file_handle,ierrmpi)
      call mpi_file_read(file_handle, file_version, 1, mpi_integer, &
           istatus, ierrmpi)
    end if
 
    call mpi_bcast(file_version,1,mpi_integer,0,icomm,ierrmpi)
 
    if (all(compatible_versions /= file_version)) then
      if (mype == 0) print *, "Unknown version, trying old snapshot reader..."
      call mpi_file_close(file_handle,ierrmpi)
      call read_snapshot_old()
 
      ! Guess snapshotnext from file name if not set
      if (snapshotnext == -1) &
           snapshotnext = get_snapshot_index(trim(restart_from_file)) + 1
      ! Set slicenext and collapsenext if not set
      if (slicenext == -1) slicenext = 0
      if (collapsenext == -1) collapsenext = 0
 
      ! Still used in convert
      snapshotini = snapshotnext-1
 
      return ! Leave this routine
    else if (mype == 0) then
      call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set, ierrmpi)
      call snapshot_read_header(file_handle, offset_tree_info, &
           offset_block_data)
    end if
 
    ! Share information about restart file
    call mpi_bcast(nw_found,1,mpi_integer,0,icomm,ierrmpi)
    call mpi_bcast(nleafs,1,mpi_integer,0,icomm,ierrmpi)
    call mpi_bcast(nparents,1,mpi_integer,0,icomm,ierrmpi)
    call mpi_bcast(it,1,mpi_integer,0,icomm,ierrmpi)
    call mpi_bcast(global_time,1,mpi_double_precision,0,icomm,ierrmpi)
 
    call mpi_bcast(snapshotnext,1,mpi_integer,0,icomm,ierrmpi)
    call mpi_bcast(slicenext,1,mpi_integer,0,icomm,ierrmpi)
    call mpi_bcast(collapsenext,1,mpi_integer,0,icomm,ierrmpi)
    call mpi_bcast(stagger_mark_dat,1,mpi_logical,0,icomm,ierrmpi)
 
    ! Allocate send/receive buffer
    n_values = count_ix(ixg^ll) * nw_found
    if(stagger_mark_dat) then
      n_values = n_values + count_ix(ixgs^ll) * nws
    end if
    allocate(w_buffer(n_values))
    allocate(w(ixg^t,1:nw_found))
 
    nleafs_active = nleafs
 
    if (mype == 0) then
      call mpi_file_seek(file_handle, offset_tree_info, &
           mpi_seek_set, ierrmpi)
    end if
 
    call read_forest(file_handle)
 
    do morton_no=morton_start(mype),morton_stop(mype)
      igrid=sfc_to_igrid(morton_no)
      call alloc_node(igrid)
    end do
 
    if (mype==0) then
      call mpi_file_seek(file_handle, offset_block_data, mpi_seek_set, ierrmpi)
 
      iread = 0
      do ipe = 0, npe-1
        do morton_no=morton_start(ipe),morton_stop(ipe)
          iread=iread+1
          itag=morton_no
 
          call mpi_file_read(file_handle,ix_buffer(1:2*ndim), 2*ndim, &
               mpi_integer, istatus,ierrmpi)
 
          ! Construct ixO^L array from number of ghost cells
          {ixomin^d = ixmlo^d - ix_buffer(^d)\}
          {ixomax^d = ixmhi^d + ix_buffer(ndim+^d)\}
          n_values = count_ix(ixo^l) * nw_found
          if(stagger_mark_dat) then
            {ixosmin^d = ixomin^d - 1\}
            {ixosmax^d = ixomax^d\}
            n_values_stagger = n_values
            n_values = n_values + count_ix(ixos^l) * nws
          end if
 
          call mpi_file_read(file_handle, w_buffer, n_values, &
               mpi_double_precision, istatus, ierrmpi)
 
          if (mype == ipe) then ! Root task
            igrid=sfc_to_igrid(morton_no)
            block=>ps(igrid)
            if(stagger_mark_dat) then
              w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values_stagger), &
                   shape(w(ixo^s, 1:nw_found)))
              if(stagger_grid) &
              ps(igrid)%ws(ixos^s,1:nws)=reshape(w_buffer(n_values_stagger+1:n_values), &
                   shape(ws(ixos^s, 1:nws)))
            else
              w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values), &
                   shape(w(ixo^s, 1:nw_found)))
            end if
            if (nw_found<nw) then
              if (associated(usr_transform_w)) then
                call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
              else
                ps(igrid)%w(ixo^s,1:nw_found)=w(ixo^s,1:nw_found)
              end if
            else if (nw_found>nw) then
              if (associated(usr_transform_w)) then
                call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
              else
                ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
              end if
            else
              ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
            end if
          else
            call mpi_send([ ixo^l, n_values ], 2*ndim+1, &
                 mpi_integer, ipe, itag, icomm, ierrmpi)
            call mpi_send(w_buffer, n_values, &
                 mpi_double_precision, ipe, itag, icomm, ierrmpi)
          end if
        end do
      end do
 
      call mpi_file_close(file_handle,ierrmpi)
 
    else ! mype > 0
 
      do morton_no=morton_start(mype),morton_stop(mype)
        igrid=sfc_to_igrid(morton_no)
        block=>ps(igrid)
        itag=morton_no
 
        call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, 0, itag, icomm,&
             iorecvstatus, ierrmpi)
        {ixomin^d = ix_buffer(^d)\}
        {ixomax^d = ix_buffer(ndim+^d)\}
        n_values = ix_buffer(2*ndim+1)
 
        call mpi_recv(w_buffer, n_values, mpi_double_precision,&
             0, itag, icomm, iorecvstatus, ierrmpi)
 
        if(stagger_mark_dat) then
          n_values_stagger = count_ix(ixo^l) * nw_found
          {ixosmin^d = ixomin^d - 1\}
          {ixosmax^d = ixomax^d\}
          w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values_stagger), &
               shape(w(ixo^s, 1:nw_found)))
          if(stagger_grid) &
          ps(igrid)%ws(ixos^s,1:nws)=reshape(w_buffer(n_values_stagger+1:n_values), &
               shape(ws(ixos^s, 1:nws)))
        else
          w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values), &
               shape(w(ixo^s, 1:nw_found)))
        end if
        if (nw_found<nw) then
          if (associated(usr_transform_w)) then
            call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
          else
            ps(igrid)%w(ixo^s,1:nw_found)=w(ixo^s,1:nw_found)
          end if
        else if (nw_found>nw) then
          if (associated(usr_transform_w)) then
            call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
          else
            ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
          end if
        else
          ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
        end if
      end do
    end if
 
    call mpi_barrier(icomm,ierrmpi)
 
  end subroutine read_snapshot
 
  subroutine read_snapshot_old()
    use mod_forest
    use mod_global_parameters
    use mod_amr_solution_node, only: alloc_node
    use mod_functions_forest, only: read_forest
 
    double precision              :: wio(ixG^T,1:nw)
    double precision              :: eqpar_dummy(100)
    integer                       :: fh, igrid, Morton_no, iread
    integer                       :: levmaxini, ndimini, ndirini
    integer                       :: nwini, neqparini, nxini^D
    integer(kind=MPI_OFFSET_KIND) :: offset
    integer                       :: istatus(MPI_STATUS_SIZE)
    integer, allocatable          :: iorecvstatus(:,:)
    integer                       :: ipe,inrecv,nrecv
    integer                       :: sendini(7+^ND)
    logical                       :: fexist
    character(len=80)             :: filename
 
    if (mype==0) then
      call mpi_file_open(mpi_comm_self,trim(restart_from_file), &
           mpi_mode_rdonly,mpi_info_null,fh,ierrmpi)
 
      offset=-int(7*size_int+size_double,kind=mpi_offset_kind)
      call mpi_file_seek(fh,offset,mpi_seek_end,ierrmpi)
 
      call mpi_file_read(fh,nleafs,1,mpi_integer,istatus,ierrmpi)
      nleafs_active = nleafs
      call mpi_file_read(fh,levmaxini,1,mpi_integer,istatus,ierrmpi)
      call mpi_file_read(fh,ndimini,1,mpi_integer,istatus,ierrmpi)
      call mpi_file_read(fh,ndirini,1,mpi_integer,istatus,ierrmpi)
      call mpi_file_read(fh,nwini,1,mpi_integer,istatus,ierrmpi)
      call mpi_file_read(fh,neqparini,1,mpi_integer,istatus,ierrmpi)
      call mpi_file_read(fh,it,1,mpi_integer,istatus,ierrmpi)
      call mpi_file_read(fh,global_time,1,mpi_double_precision,istatus,ierrmpi)
 
      ! check if settings are suitable for restart
      if (levmaxini>refine_max_level) then
        write(*,*) "number of levels in restart file = ",levmaxini
        write(*,*) "refine_max_level = ",refine_max_level
        call mpistop("refine_max_level < number of levels in restart file")
      end if
      if (ndimini/=ndim) then
        write(*,*) "ndim in restart file = ",ndimini
        write(*,*) "ndim = ",ndim
        call mpistop("reset ndim to ndim in restart file")
      end if
      if (ndirini/=ndir) then
        write(*,*) "ndir in restart file = ",ndirini
        write(*,*) "ndir = ",ndir
        call mpistop("reset ndir to ndir in restart file")
      end if
      if (nw/=nwini) then
        write(*,*) "nw=",nw," and nw in restart file=",nwini
        call mpistop("currently, changing nw at restart is not allowed")
      end if
 
      offset=offset-int(ndimini*size_int+neqparini*size_double,kind=mpi_offset_kind)
      call mpi_file_seek(fh,offset,mpi_seek_end,ierrmpi)
 
      {call mpi_file_read(fh,nxini^d,1,mpi_integer,istatus,ierrmpi)\}
      if (ixghi^d/=nxini^d+2*nghostcells|.or.) then
        write(*,*) "Error: reset resolution to ",nxini^d+2*nghostcells
        call mpistop("change with setamrvac")
      end if
 
      call mpi_file_read(fh,eqpar_dummy,neqparini, &
           mpi_double_precision,istatus,ierrmpi)
    end if
 
    ! broadcast the global parameters first
    if (npe>1) then
      if (mype==0) then
        sendini=(/nleafs,levmaxini,ndimini,ndirini,nwini,neqparini,it ,^d&nxini^d /)
      end if
      call mpi_bcast(sendini,7+^nd,mpi_integer,0,icomm,ierrmpi)
      nleafs=sendini(1);levmaxini=sendini(2);ndimini=sendini(3);
      ndirini=sendini(4);nwini=sendini(5);
      neqparini=sendini(6);it=sendini(7);
      nxini^d=sendini(7+^d);
      nleafs_active = nleafs
      call mpi_bcast(global_time,1,mpi_double_precision,0,icomm,ierrmpi)
    end if
 
    if (mype == 0) then
      offset = int(size_block_io,kind=mpi_offset_kind) * &
           int(nleafs,kind=mpi_offset_kind)
      call mpi_file_seek(fh,offset,mpi_seek_set,ierrmpi)
    end if
 
    call read_forest(fh)
 
    do morton_no=morton_start(mype),morton_stop(mype)
      igrid=sfc_to_igrid(morton_no)
      call alloc_node(igrid)
    end do
 
    if (mype==0)then
      iread=0
 
      do morton_no=morton_start(0),morton_stop(0)
        igrid=sfc_to_igrid(morton_no)
        iread=iread+1
        offset=int(size_block_io,kind=mpi_offset_kind) &
             *int(morton_no-1,kind=mpi_offset_kind)
        call mpi_file_read_at(fh,offset,ps(igrid)%w,1,type_block_io, &
             istatus,ierrmpi)
      end do
      if (npe>1) then
        do ipe=1,npe-1
          do morton_no=morton_start(ipe),morton_stop(ipe)
            iread=iread+1
            itag=morton_no
            offset=int(size_block_io,kind=mpi_offset_kind)&
                 *int(morton_no-1,kind=mpi_offset_kind)
            call mpi_file_read_at(fh,offset,wio,1,type_block_io,&
                 istatus,ierrmpi)
            call mpi_send(wio,1,type_block_io,ipe,itag,icomm,ierrmpi)
          end do
        end do
      end if
      call mpi_file_close(fh,ierrmpi)
    else
      nrecv=(morton_stop(mype)-morton_start(mype)+1)
      allocate(iorecvstatus(mpi_status_size,nrecv))
      inrecv=0
      do morton_no=morton_start(mype),morton_stop(mype)
        igrid=sfc_to_igrid(morton_no)
        itag=morton_no
        inrecv=inrecv+1
        call mpi_recv(ps(igrid)%w,1,type_block_io,0,itag,icomm,&
             iorecvstatus(:,inrecv),ierrmpi)
      end do
      deallocate(iorecvstatus)
    end if
 
    call mpi_barrier(icomm,ierrmpi)
 
  end subroutine read_snapshot_old
 
  !> Write volume-averaged values and other information to the log file
  subroutine printlog_default
 
    use mod_timing
    use mod_forest, only: nleafs, nleafs_active, nleafs_level
    use mod_global_parameters
 
    double precision     :: dtTimeLast, now, cellupdatesPerSecond
    double precision     :: activeBlocksPerCore, wctPerCodeTime, timeToFinish
    double precision     :: wmean(1:nw), total_volume
    double precision     :: volume_coverage(refine_max_level)
    integer              :: i, iw, level
    integer              :: nx^D, nc, ncells, dit
    integer              :: amode, istatus(MPI_STATUS_SIZE)
    integer, parameter   :: my_unit = 20
    logical, save        :: opened  = .false.
    logical              :: fileopen
    character(len=40)    :: fmt_string
    character(len=80)    :: filename
    character(len=2048)  :: line
 
    ! Compute the volume-average of w**1 = w
    call get_volume_average(1, wmean, total_volume)
 
    ! Compute the volume coverage
    call get_volume_coverage(volume_coverage)
 
    if (mype == 0) then
 
       ! To compute cell updates per second, we do the following:
       nx^d=ixmhi^d-ixmlo^d+1;
       nc={nx^d*}
       ncells = nc * nleafs_active
 
       ! assumes the number of active leafs haven't changed since last compute.
       now        = mpi_wtime()
       dit        = it - ittimelast
       dttimelast = now - timelast
       ittimelast = it
       timelast   = now
       cellupdatespersecond = dble(ncells) * dble(nstep) * &
            dble(dit) / (dttimelast * dble(npe))
 
       ! blocks per core:
       activeblockspercore = dble(nleafs_active) / dble(npe)
 
       ! Wall clock time per code time unit in seconds:
       wctpercodetime = dttimelast / max(dit * dt, epsilon(1.0d0))
 
       ! Wall clock time to finish in hours:
       timetofinish = (time_max - global_time) * wctpercodetime / 3600.0d0
 
       ! On first entry, open the file and generate the header
       if (.not. opened) then
 
          filename = trim(base_filename) // ".log"
 
          ! Delete the log when not doing a restart run
          if (restart_from_file == undefined) then
             open(unit=my_unit,file=trim(filename),status='replace')
             close(my_unit, status='delete')
          end if
 
          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, log_fh, ierrmpi)
 
          opened = .true.
 
          ! Start of file headern
          line = "it global_time dt"
          do level=1,nw
             i = len_trim(line) + 2
             write(line(i:),"(a,a)") trim(cons_wnames(level)), " "
          end do
 
          ! Volume coverage per level
          do level = 1, refine_max_level
             i = len_trim(line) + 2
             write(line(i:), "(a,i0)") "c", level
          end do
 
          ! Cell counts per level
          do level=1,refine_max_level
             i = len_trim(line) + 2
             write(line(i:), "(a,i0)") "n", level
          end do
 
          ! Rest of file header
          line = trim(line) // " | Xload Xmemory 'Cell_Updates /second/core'"
          line = trim(line) // " 'Active_Blocks/Core' 'Wct Per Code Time [s]'"
          line = trim(line) // " 'TimeToFinish [hrs]'"
 
          ! Only write header if not restarting
          if (restart_from_file == undefined .or. reset_time) then
            call mpi_file_write(log_fh, trim(line) // new_line('a'), &
                 len_trim(line)+1, mpi_character, istatus, ierrmpi)
          end if
       end if
 
       ! Construct the line to be added to the log
 
       fmt_string = '(' // fmt_i // ',2' // fmt_r // ')'
       write(line, fmt_string) it, global_time, dt
       i = len_trim(line) + 2
 
       write(fmt_string, '(a,i0,a)') '(', nw, fmt_r // ')'
       write(line(i:), fmt_string) wmean(1:nw)
       i = len_trim(line) + 2
 
       write(fmt_string, '(a,i0,a)') '(', refine_max_level, fmt_r // ')'
       write(line(i:), fmt_string) volume_coverage(1:refine_max_level)
       i = len_trim(line) + 2
 
       write(fmt_string, '(a,i0,a)') '(', refine_max_level, fmt_i // ')'
       write(line(i:), fmt_string) nleafs_level(1:refine_max_level)
       i = len_trim(line) + 2
 
       fmt_string = '(a,6' // fmt_r2 // ')'
       write(line(i:), fmt_string) '| ', xload, xmemory, cellupdatespersecond, &
            activeblockspercore, wctpercodetime, timetofinish
 
       call mpi_file_write(log_fh, trim(line) // new_line('a') , &
            len_trim(line)+1, mpi_character, istatus, ierrmpi)
    end if
 
  end subroutine printlog_default
 
  !> Print a log that can be used to check whether the code still produces the
  !> same output (regression test)
  subroutine printlog_regression_test()
    use mod_global_parameters
 
    double precision   :: modes(nw, 2), volume
    integer, parameter :: n_modes = 2
    integer            :: power
    integer              :: amode, istatus(MPI_STATUS_SIZE)
    logical, save      :: file_open = .false.
    character(len=40)  :: fmt_string
    character(len=2048)  :: line
    character(len=80)    :: filename
 
    do power = 1, n_modes
       call get_volume_average(power, modes(:, power), volume)
    end do
 
    if (mype == 0) then
       if (.not. file_open) then
          filename = trim(base_filename) // ".log"
          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, log_fh, ierrmpi)
          file_open = .true.
 
          line= "# time mean(w) mean(w**2)"
          call mpi_file_write(log_fh, trim(line) // new_line('a'), &
                 len_trim(line)+1, mpi_character, istatus, ierrmpi)
       end if
 
       write(fmt_string, "(a,i0,a)") "(", nw * n_modes + 1, fmt_r // ")"
       write(line, fmt_string) global_time, modes
       call mpi_file_write(log_fh, trim(line) // new_line('a') , &
            len_trim(line)+1, mpi_character, istatus, ierrmpi)
    end if
  end subroutine printlog_regression_test
 
  !> Compute mean(w**power) over the leaves of the grid. The first mode
  !> (power=1) corresponds to the mean, the second to the mean squared values
  !> and so on.
  subroutine get_volume_average(power, mode, volume)
    use mod_global_parameters
 
    integer, intent(in)           :: power     !< Which mode to compute
    double precision, intent(out) :: mode(nw)  !< The computed mode
    double precision, intent(out) :: volume    !< The total grid volume
 
    double precision              :: wsum(nw+1)
    double precision              :: dsum_recv(1:nw+1)
    integer                       :: iigrid, igrid, iw
 
    wsum(:) = 0
 
    ! Loop over all the grids
    do iigrid = 1, igridstail
       igrid = igrids(iigrid)
 
       ! Store total volume in last element
       wsum(nw+1) = wsum(nw+1) + sum(ps(igrid)%dvolume(ixm^t))
 
       ! Compute the modes of the cell-centered variables, weighted by volume
       do iw = 1, nw
          wsum(iw) = wsum(iw) + &
               sum(ps(igrid)%dvolume(ixm^t)*ps(igrid)%w(ixm^t,iw)**power)
       end do
    end do
 
    ! Make the information available on all tasks
    call mpi_allreduce(wsum, dsum_recv, nw+1, mpi_double_precision, &
         mpi_sum, icomm, ierrmpi)
 
    ! Set the volume and the average
    volume = dsum_recv(nw+1)
    mode   = dsum_recv(1:nw) / volume
 
  end subroutine get_volume_average
 
  !> Compute how much of the domain is covered by each grid level. This routine
  !> does not take a non-Cartesian geometry into account.
  subroutine get_volume_coverage(vol_cov)
    use mod_global_parameters
 
    double precision, intent(out) :: vol_cov(1:refine_max_level)
    double precision              :: dsum_recv(1:refine_max_level)
    integer                       :: iigrid, igrid, iw, level
 
    ! First determine the total 'flat' volume in each level
    vol_cov(1:refine_max_level)=zero
 
    do iigrid = 1, igridstail
       igrid          = igrids(iigrid);
       level          = node(plevel_,igrid)
       vol_cov(level) = vol_cov(level)+ &
            {(rnode(rpxmax^d_,igrid)-rnode(rpxmin^d_,igrid))|*}
    end do
 
    ! Make the information available on all tasks
    call mpi_allreduce(vol_cov, dsum_recv, refine_max_level, mpi_double_precision, &
         mpi_sum, icomm, ierrmpi)
 
    ! Normalize
    vol_cov = dsum_recv / sum(dsum_recv)
  end subroutine get_volume_coverage
 
  !> Compute the volume average of func(w) over the leaves of the grid.
  subroutine get_volume_average_func(func, f_avg, volume)
    use mod_global_parameters
 
    interface
       pure function func(w_vec, w_size) result(val)
         integer, intent(in)          :: w_size
         double precision, intent(in) :: w_vec(w_size)
         double precision             :: val
       end function func
    end interface
    double precision, intent(out) :: f_avg  !< The volume average of func
    double precision, intent(out) :: volume    !< The total grid volume
    double precision              :: wsum(2)
    double precision              :: dsum_recv(2)
    integer                       :: iigrid, igrid, i^D
 
    wsum(:) = 0
 
    ! Loop over all the grids
    do iigrid = 1, igridstail
       igrid = igrids(iigrid)
 
       ! Store total volume in last element
       wsum(2) = wsum(2) + sum(ps(igrid)%dvolume(ixm^t))
 
       ! Compute the modes of the cell-centered variables, weighted by volume
       {do i^d = ixmlo^d, ixmhi^d\}
       wsum(1) = wsum(1) + ps(igrid)%dvolume(i^d) * &
            func(ps(igrid)%w(i^d, :), nw)
       {end do\}
    end do
 
    ! Make the information available on all tasks
    call mpi_allreduce(wsum, dsum_recv, 2, mpi_double_precision, &
         mpi_sum, icomm, ierrmpi)
 
    ! Set the volume and the average
    volume = dsum_recv(2)
    f_avg  = dsum_recv(1) / volume
 
  end subroutine get_volume_average_func
 
  !> Compute global maxima of iw variables over the leaves of the grid.
  subroutine get_global_maxima(wmax,psa)
    use mod_global_parameters
 
    double precision, intent(out) :: wmax(nw)  !< The global maxima
    type(state), target :: psa(max_blocks)
 
    double precision              :: wmax_mype(nw),wmax_recv(nw)
    integer                       :: iigrid, igrid, iw
 
    wmax_mype(1:nw) = -bigdouble
 
    ! Loop over all the grids
    do iigrid = 1, igridstail
       igrid = igrids(iigrid)
       do iw = 1, nw
          wmax_mype(iw)=max(wmax_mype(iw),maxval(psa(igrid)%w(ixm^t,iw)))
       end do
    end do
 
    ! Make the information available on all tasks
    call mpi_allreduce(wmax_mype, wmax_recv, nw, mpi_double_precision, &
         mpi_max, icomm, ierrmpi)
 
    wmax(1:nw)=wmax_recv(1:nw)
 
  end subroutine get_global_maxima
 
  !> Compute global minima of iw variables over the leaves of the grid.
  subroutine get_global_minima(wmin,psa)
    use mod_global_parameters
 
    double precision, intent(out) :: wmin(nw)  !< The global maxima
    type(state), target :: psa(max_blocks)
 
    double precision              :: wmin_mype(nw),wmin_recv(nw)
    integer                       :: iigrid, igrid, iw
 
    wmin_mype(1:nw) = bigdouble
 
    ! Loop over all the grids
    do iigrid = 1, igridstail
       igrid = igrids(iigrid)
       do iw = 1, nw
          wmin_mype(iw)=min(wmin_mype(iw),minval(psa(igrid)%w(ixm^t,iw)))
       end do
    end do
 
    ! Make the information available on all tasks
    call mpi_allreduce(wmin_mype, wmin_recv, nw, mpi_double_precision, &
         mpi_min, icomm, ierrmpi)
 
    wmin(1:nw)=wmin_recv(1:nw)
 
  end subroutine get_global_minima
 
end module mod_input_output