diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index a55d73248..de2c2690e 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -10,7 +10,7 @@ jobs:
       matrix:
         os: [ubuntu-latest, macOS-latest]
       fail-fast: false
-    timeout-minutes: 30
+    timeout-minutes: 40
 
     steps:
       - uses: actions/checkout@v2
diff --git a/README.md b/README.md
index 731ce3061..ec27bc97a 100644
--- a/README.md
+++ b/README.md
@@ -32,13 +32,10 @@ The full documentation is online at [https://mabarnes.github.io/moment_kinetics]
     ```
    this significantly decreases the load time but prevents code changes from taking effect when `moment_kinetics.so` is used without repeating the precompilation (to use this option, add an option `-Jmoment_kinetics.so` when starting julia).
 4) To run julia with optimization, type
-    ```
-    $ julia -O3 --project run_moment_kinetics.jl
-    ```
-    Default input options are specified in `moment_kinetics_input.jl`. The defaults can be modified for a particular run by setting options in a TOML file, for example `input.toml`, which can be passed as an argument
     ```
     $ julia -O3 --project run_moment_kinetics.jl input.toml
     ```
+    Options are specified in a TOML file, e.g. `input.toml` here. The defaults are specified in `moment_kinetics_input.jl`.
     * To run in parallel, just put `mpirun -np <n>` in front of the call you would normally use, with `<n>` the number of processes to use.
     * It may be more convenient when running `moment_kinetics` more than once to work from the Julia REPL, e.g.
         ```
@@ -50,27 +47,49 @@ The full documentation is online at [https://mabarnes.github.io/moment_kinetics]
         ```
         julia> run_moment_kinetics("input.toml")
         ```
-5) To make plots and calculate frequencies/growth rates, run
+5) To restart a simulation using `input.toml` from the last time point in the existing run directory,
+    ```
+    $ julia -O3 --project run_moment_kinetics --restart input.toml
+    ```
+    or to restart from a specific output file - either from the same run or (if the settings are compatible) a different one - here `runs/example/example.dfns.h5`
+    ```
+    $ julia -O3 --project run_moment_kinetics input.toml runs/example/example.dfns.h5
+    ```
+    The output file must include distribution functions. When not using parallel I/O there will be multiple output files from different MPI ranks - any one of these can be passed.
+    * To do the same from the Julia REPL
+        ```
+        $ julia -O3 --project
+        julia> run_moment_kinetics("input.toml", restart=true)
+        ```
+        or
+        ```
+        julia> run_moment_kinetics("input.toml", restart="runs/example/example.dfns.h5")
+        ```
+    * When calling the `run_moment_kinetics()` function you can also choose a particular time index to restart from, e.g.
+        ```
+        julia> run_moment_kinetics("input.toml", restart="runs/example/example.dfns.h5", restart_time_index=42)
+        ```
+6) To make plots and calculate frequencies/growth rates, run
     ```
     $ julia --project run_post_processing.jl runs/<directory to process>
     ```
     passing the directory to process as a command line argument. Input options for post-processing can be specified in post_processing_input.jl.
 
-6) Parameter scans (see [Running parameter scans](#running-parameter-scans)) or performance tests can be performed by running
+7) Parameter scans (see [Running parameter scans](#running-parameter-scans)) or performance tests can be performed by running
     ```
     $ julia -O3 --project driver.jl
     ```
     If running a scan, it can be parallelised by passing the number of processors as an argument. Scan options are set in `scan_inputs.jl`.
 
-7) Post processing can be done for several directories at once using
+8) Post processing can be done for several directories at once using
     ```
     $ julia --project post_processing_driver.jl runs/<directory1> runs/<directory2> ...
     ```
     passing the directories to process as command line arguments. Optionally pass a number as the first argument to parallelise post processing of different directories. Input options for post-processing can be specified in `post_processing_input.jl`.
 
-8) In the course of development, it is sometimes helpful to upgrade the Julia veriosn. Upgrading the version of Julia or upgrading packages may require a fresh installation of `moment_kinetics`. To make a fresh install with the latest package versions it is necessary to remove (or rename) the `Manifest.jl` file in the main directory, and generate a new `Manifest.jl` with step 1) above. It can sometimes be necessary to remove or rename the `.julia/` folder in your root directory for this step to be successful.
+9) In the course of development, it is sometimes helpful to upgrade the Julia veriosn. Upgrading the version of Julia or upgrading packages may require a fresh installation of `moment_kinetics`. To make a fresh install with the latest package versions it is necessary to remove (or rename) the `Manifest.jl` file in the main directory, and generate a new `Manifest.jl` with step 1) above. It can sometimes be necessary to remove or rename the `.julia/` folder in your root directory for this step to be successful.
 
-9) One may have to set an environment variable to avoid error messages from the Qt library. If you execute the command
+10) One may have to set an environment variable to avoid error messages from the Qt library. If you execute the command
 
     ```
     $ julia --project run_post_processing.jl runs/your_run_dir/
diff --git a/machines/README.md b/machines/README.md
index ac977fdf3..01031eda2 100644
--- a/machines/README.md
+++ b/machines/README.md
@@ -67,7 +67,7 @@ $ ./submit-restart.sh <path to input file>.toml
 ```
 will submit a job to run and post-process a restart using input file. The
 simulation will restart from the last time point of the previous run
-(`restart_moment_kinetics.jl` supports more flexibility, but for now you would
+(`run_moment_kinetics.jl` supports more flexibility, but for now you would
 need to write your own submission script to pass the options needed for that).
 
 Default parameters for the runs (number of nodes, time limit, etc.) were set up
diff --git a/machines/archer/jobscript-restart.template b/machines/archer/jobscript-restart.template
index c4018292c..52501d169 100644
--- a/machines/archer/jobscript-restart.template
+++ b/machines/archer/jobscript-restart.template
@@ -22,6 +22,6 @@ export SRUN_CPUS_PER_TASK=$SLURM_CPUS_PER_TASK
 
 echo "running INPUTFILE $(date)"
 
-srun --distribution=block:block --hint=nomultithread --ntasks=$SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no restart_moment_kinetics.jl INPUTFILE RESTARTFROM
+srun --distribution=block:block --hint=nomultithread --ntasks=$SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no run_moment_kinetics.jl --restart INPUTFILE RESTARTFROM
 
 echo "finished INPUTFILE $(date)"
diff --git a/machines/marconi/jobscript-restart.template b/machines/marconi/jobscript-restart.template
index 53ba20059..420687c93 100644
--- a/machines/marconi/jobscript-restart.template
+++ b/machines/marconi/jobscript-restart.template
@@ -18,6 +18,6 @@ source julia.env
 
 echo "running INPUTFILE $(date)"
 
-mpirun -np $SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no restart_moment_kinetics.jl INPUTFILE RESTARTFROM
+mpirun -np $SLURM_NTASKS bin/julia -Jmoment_kinetics.so --project -O3 --check-bounds=no run_moment_kinetics.jl --restart INPUTFILE RESTARTFROM
 
 echo "finished INPUTFILE $(date)"
diff --git a/restart_moment_kinetics.jl b/restart_moment_kinetics.jl
deleted file mode 100644
index c82916954..000000000
--- a/restart_moment_kinetics.jl
+++ /dev/null
@@ -1,7 +0,0 @@
-# provide option of running from command line via 'julia run_moment_kinetics.jl'
-using Pkg
-Pkg.activate(".")
-
-using moment_kinetics
-
-restart_moment_kinetics()
diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 0908ac98e..30209ec7b 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -41,9 +41,25 @@ end
 
 """
 initialize chebyshev grid scaled to interval [-box_length/2, box_length/2]
-"""
-function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n,
-			irank, box_length, imin, imax)
+we no longer pass the box_length to this function, but instead pass precomputed
+arrays element_scale and element_shift that are needed to compute the grid.
+
+ngrid -- number of points per element (including boundary points)
+nelement_local -- number of elements in the local (distributed memory MPI) grid
+n -- total number of points in the local grid (excluding duplicate points)
+element_scale -- the scale factor in the transform from the coordinates 
+                 where the element limits are -1, 1 to the coordinate where
+                 the limits are Aj = coord.grid[imin[j]-1] and Bj = coord.grid[imax[j]]
+                 element_scale = 0.5*(Bj - Aj)
+element_shift -- the centre of the element in the extended grid coordinate
+                 element_shift = 0.5*(Aj + Bj)
+imin -- the array of minimum indices of each element on the extended grid.
+        By convention, the duplicated points are not included, so for element index j > 1
+        the lower boundary point is actually imin[j] - 1
+imax -- the array of maximum indices of each element on the extended grid.
+"""
+function scaled_chebyshev_grid(ngrid, nelement_local, n,
+			element_scale, element_shift, imin, imax)
     # initialize chebyshev grid defined on [1,-1]
     # with n grid points chosen to facilitate
     # the fast Chebyshev transform (aka the discrete cosine transform)
@@ -52,19 +68,14 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n,
     chebyshev_grid = chebyshevpoints(ngrid)
     # create array for the full grid
     grid = allocate_float(n)
-    # setup the scale factor by which the Chebyshev grid on [-1,1]
-    # is to be multiplied to account for the full domain [-L/2,L/2]
-    # and the splitting into nelement elements with ngrid grid points
-    scale_factor = 0.5*box_length/float(nelement_global)
+    
     # account for the fact that the minimum index needed for the chebyshev_grid
     # within each element changes from 1 to 2 in going from the first element
     # to the remaining elements
     k = 1
     @inbounds for j ∈ 1:nelement_local
-        #wgts[imin[j]:imax[j]] .= sqrt.(1.0 .- reverse(chebyshev_grid)[k:ngrid].^2) * scale_factor
-        # amount by which to shift the centre of this element from zero
-        iel_global = j + irank*nelement_local
-        shift = box_length*((float(iel_global)-0.5)/float(nelement_global) - 0.5)
+        scale_factor = element_scale[j]
+        shift = element_shift[j]
         # reverse the order of the original chebyshev_grid (ran from [1,-1])
         # and apply the scale factor and shift
         grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift
@@ -72,7 +83,7 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n,
         # to avoid double-counting boundary element
         k = 2
     end
-    wgts = clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, scale_factor)
+    wgts = clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_scale)
     return grid, wgts
 end
 
@@ -88,11 +99,9 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info)
     # check array bounds
     @boundscheck nelement == size(chebyshev.f,2) || throw(BoundsError(chebyshev.f))
     @boundscheck nelement == size(df,2) && coord.ngrid == size(df,1) || throw(BoundsError(df))
-    # note that one must multiply by 2*nelement/L to get derivative
-    # in scaled coordinate
-    scale_factor = 2.0*float(coord.nelement_global)/coord.L
-	# scale factor is (length of a single element/2)^{-1}
-	
+    # note that one must multiply by a coordinate transform factor 1/element_scale[j] 
+    # for each element j to get derivative on the extended grid
+    
     # variable k will be used to avoid double counting of overlapping point
     # at element boundaries (see below for further explanation)
     k = 0
@@ -112,7 +121,7 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info)
         # and multiply by scaling factor needed to go
         # from Chebyshev z coordinate to actual z
         for i ∈ 1:coord.ngrid
-            df[i,j] *= scale_factor
+            df[i,j] /= coord.element_scale[j]
         end
         k = 1
     end
@@ -323,23 +332,21 @@ end
 returns wgts array containing the integration weights associated
 with all grid points for Clenshaw-Curtis quadrature
 """
-function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, scale_factor)
+function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_scale)
     # create array containing the integration weights
     wgts = zeros(mk_float, n)
     # calculate the modified Chebshev moments of the first kind
     μ = chebyshevmoments(ngrid)
+    # calculate the raw weights for a normalised grid on [-1,1]
+    w = clenshawcurtisweights(μ)
     @inbounds begin
         # calculate the weights within a single element and
         # scale to account for modified domain (not [-1,1])
-        wgts[1:ngrid] = clenshawcurtisweights(μ)*scale_factor
+        wgts[1:ngrid] = w*element_scale[1]
         if nelement_local > 1
-            # account for double-counting of points at inner element boundaries
-            wgts[ngrid] *= 2
             for j ∈ 2:nelement_local
-                wgts[imin[j]:imax[j]] .= wgts[2:ngrid]
+                wgts[imin[j]-1:imax[j]] .+= w*element_scale[j]
             end
-            # remove double-counting of outer element boundary for last element
-            wgts[n] *= 0.5
         end
     end
     return wgts
diff --git a/src/command_line_options.jl b/src/command_line_options.jl
index ddee48f9d..430fbb8da 100644
--- a/src/command_line_options.jl
+++ b/src/command_line_options.jl
@@ -16,7 +16,7 @@ const s = ArgParseSettings()
         arg_type = String
         default = nothing
     "restartfile"
-        help = "Name of NetCDF file to restart from"
+        help = "Name of output file (HDF5 or NetCDF) to restart from"
         arg_type = String
         default = nothing
     "--debug", "-d"
@@ -24,14 +24,18 @@ const s = ArgParseSettings()
                "integer values activate more checks (and increase run time)"
         arg_type = Int
         default = 0
-    # Options for tests
-    "--long"
-        help = "Include more tests, increasing test run time."
+    "--restart"
+        help = "Restart from latest output file in run directory (ignored if " *
+               "`restartfile` is passed)"
         action = :store_true
     "--restart-time-index"
         help = "Time index in output file to restart from, defaults to final time point"
         arg_type = Int
         default = -1
+    # Options for tests
+    "--long"
+        help = "Include more tests, increasing test run time."
+        action = :store_true
     "--verbose", "-v"
         help = "Print verbose output from tests."
         action = :store_true
diff --git a/src/communication.jl b/src/communication.jl
index 8e9033e19..3259a4a06 100644
--- a/src/communication.jl
+++ b/src/communication.jl
@@ -125,8 +125,12 @@ function setup_distributed_memory_MPI(z_nelement_global,z_nelement_local,r_nelem
     end
 	# throw an error if user specified information is inconsistent
     if (nrank_per_zr_block*nblocks < nrank_global)
-        error("ERROR: You must choose global number of processes to be an integer multiple of the number of \n
-               nblocks = (r_nelement_global/r_nelement_local)*(z_nelement_global/z_nelement_local)")
+        error("ERROR: You must choose global number of processes to be an integer "
+              * "multiple of the number of\n"
+              * "nblocks($nblocks) = (r_nelement_global($r_nelement_global)/"
+              * "r_nelement_local($r_nelement_local))*"
+              * "(z_nelement_global($z_nelement_global)/"
+              * "z_nelement_local($z_nelement_local))")
     end
     
     # assign information regarding shared-memory blocks
diff --git a/src/coordinates.jl b/src/coordinates.jl
index 704f67023..6ad0164b7 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -4,6 +4,7 @@ module coordinates
 
 export define_coordinate, write_coordinate
 export equally_spaced_grid
+export set_element_boundaries
 
 using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_int
@@ -84,6 +85,12 @@ struct coordinate
     local_io_range::UnitRange{Int64}
     # global range to write into in output file
     global_io_range::UnitRange{Int64}
+    # scale for each element
+    element_scale::Array{mk_float,1}
+    # shift for each element
+    element_shift::Array{mk_float,1}
+    # option used to set up element spacing
+    element_spacing_option::String
 end
 
 """
@@ -105,10 +112,14 @@ function define_coordinate(input, parallel_io::Bool=false)
     # obtain (local) index mapping from the grid within each element
     # to the full grid
     imin, imax = elemental_to_full_grid_map(input.ngrid, input.nelement_local)
+    # initialise the data used to construct the grid
+    # boundaries for each element
+    element_boundaries = set_element_boundaries(input.nelement_global, input.L, input.element_spacing_option)
+    # shift and scale factors for each local element
+    element_scale, element_shift = set_element_scale_and_shift(input.nelement_global, input.nelement_local, input.irank, element_boundaries)
     # initialize the grid and the integration weights associated with the grid
     # also obtain the Chebyshev theta grid and spacing if chosen as discretization option
-    grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_global,
-        input.nelement_local, n_global, n_local, input.irank, input.L,
+    grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_local, n_global, n_local, input.irank, input.L, element_scale, element_shift,
         imin, imax, igrid, input.discretization, input.name)
     # calculate the widths of the cells between neighboring grid points
     cell_width = grid_spacing(grid, n_local)
@@ -126,7 +137,6 @@ function define_coordinate(input, parallel_io::Bool=false)
     # endpoints, so only two pieces of information must be shared
     send_buffer = allocate_float(1)
     receive_buffer = allocate_float(1)
-
     # Add some ranges to support parallel file io
     if !parallel_io
         # No parallel io, just write everything
@@ -149,7 +159,7 @@ function define_coordinate(input, parallel_io::Bool=false)
         cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option,
         input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch),
         scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm,
-        local_io_range, global_io_range)
+        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)
 
     if input.discretization == "chebyshev_pseudospectral" && coord.n > 1
         # create arrays needed for explicit Chebyshev pseudospectral treatment in this
@@ -167,10 +177,56 @@ function define_coordinate(input, parallel_io::Bool=false)
     return coord, spectral
 end
 
+function set_element_boundaries(nelement_global, L, element_spacing_option)
+    # set global element boundaries
+    element_boundaries = allocate_float(nelement_global+1)
+    if element_spacing_option == "sqrt" && nelement_global > 3
+        # number of boundaries of sqrt grid
+        nsqrt = floor(mk_int,(nelement_global)/2) + 1
+        if nelement_global%2 > 0 # odd
+            if nsqrt < 3
+                fac = 2.0/3.0
+            else
+                fac = 1.0/( 3.0/2.0 - 0.5*((nsqrt-2)/(nsqrt-1))^2)
+            end
+        else
+            fac = 1.0
+        end
+        
+        for j in 1:nsqrt
+            element_boundaries[j] = -(L/2.0) + fac*(L/2.0)*((j-1)/(nsqrt-1))^2
+        end
+        for j in 1:nsqrt
+            element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - fac*(L/2.0)*((j-1)/(nsqrt-1))^2
+        end
+        
+    elseif element_spacing_option == "uniform" || (element_spacing_option == "sqrt" && nelement_global < 4) # uniform spacing 
+        for j in 1:nelement_global+1
+            element_boundaries[j] = L*((j-1)/(nelement_global) - 0.5)
+        end
+    else 
+        println("ERROR: element_spacing_option: ",element_spacing_option, " not supported")
+    end
+    return element_boundaries
+end
+
+function set_element_scale_and_shift(nelement_global, nelement_local, irank, element_boundaries)
+    element_scale = allocate_float(nelement_local)
+    element_shift = allocate_float(nelement_local)
+    
+    for j in 1:nelement_local
+        iel_global = j + irank*nelement_local
+        upper_boundary = element_boundaries[iel_global+1]
+        lower_boundary = element_boundaries[iel_global]
+        element_scale[j] = 0.5*(upper_boundary-lower_boundary)
+        element_shift[j] = 0.5*(upper_boundary+lower_boundary)
+    end
+    return element_scale, element_shift
+end
 """
 setup a grid with n_global grid points on the interval [-L/2,L/2]
 """
-function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, irank, L,
+function init_grid(ngrid, nelement_local, n_global, n_local, irank, L, element_scale, element_shift,
                    imin, imax, igrid, discretization, name)
     uniform_grid = equally_spaced_grid(n_global, n_local, irank, L)
     uniform_grid_shifted = equally_spaced_grid_shifted(n_global, n_local, irank, L)
@@ -187,7 +243,7 @@ function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, ir
     elseif discretization == "chebyshev_pseudospectral"
         if name == "vperp"
             # initialize chebyshev grid defined on [-L/2,L/2]
-            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n_local, irank, L, imin, imax)
+            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax)
             grid .= grid .+ L/2.0 # shift to [0,L] appropriate to vperp variable
             wgts = 2.0 .* wgts .* grid # to include 2 vperp in jacobian of integral
                                         # see note above on normalisation
@@ -198,7 +254,7 @@ function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, ir
             # needed to obtain Chebyshev spectral coefficients
             # 'wgts' are the integration weights attached to each grid points
             # that are those associated with Clenshaw-Curtis quadrature
-            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n_local, irank, L, imin, imax)
+            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax)
         end
     elseif discretization == "finite_difference"
         if name == "vperp"
diff --git a/src/file_io.jl b/src/file_io.jl
index bb53a90f2..d77ae35dd 100644
--- a/src/file_io.jl
+++ b/src/file_io.jl
@@ -521,6 +521,10 @@ function define_io_coordinate!(parent, coord, coord_name, description, parallel_
         write_single_value!(group, "bc", coord.bc; parallel_io=parallel_io,
                             description="boundary condition for $coord_name")
 
+        # write the element spacing option for the coordinate
+        write_single_value!(group, "element_spacing_option", coord.element_spacing_option; parallel_io=parallel_io,
+                            description="element_spacing_option for $coord_name")
+
         return group
     end
 
diff --git a/src/input_structs.jl b/src/input_structs.jl
index ed3d20914..eb4707665 100644
--- a/src/input_structs.jl
+++ b/src/input_structs.jl
@@ -130,6 +130,8 @@ mutable struct grid_input_mutable
     bc::String
     # mutable struct containing advection speed options
     advection::advection_input_mutable
+    # string option determining boundary spacing
+    element_spacing_option::String
 end
 
 """
@@ -159,6 +161,8 @@ struct grid_input
     advection::advection_input
     # MPI communicator
     comm::MPI.Comm
+    # string option determining boundary spacing
+    element_spacing_option::String
 end
 
 """
diff --git a/src/load_data.jl b/src/load_data.jl
index 206873bbb..5507965b1 100644
--- a/src/load_data.jl
+++ b/src/load_data.jl
@@ -213,6 +213,7 @@ function load_coordinate_data(fid, name; printout=false)
     discretization = load_variable(coord_group, "discretization")
     fd_option = load_variable(coord_group, "fd_option")
     bc = load_variable(coord_group, "bc")
+    element_spacing_option = load_variable(coord_group, "element_spacing_option")
 
     nelement_local = nothing
     if n_local == 1 && ngrid == 1
@@ -225,11 +226,10 @@ function load_coordinate_data(fid, name; printout=false)
     else
         nelement_global = (n_global-1) ÷ (ngrid-1)
     end
-
     # Define input to create coordinate struct
     input = grid_input(name, ngrid, nelement_global, nelement_local, nrank, irank, L,
                        discretization, fd_option, bc, advection_input("", 0.0, 0.0, 0.0),
-                       MPI.COMM_NULL)
+                       MPI.COMM_NULL, element_spacing_option)
 
     coord, spectral = define_coordinate(input)
 
diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl
index 0b475f183..0611eba69 100644
--- a/src/moment_kinetics.jl
+++ b/src/moment_kinetics.jl
@@ -2,7 +2,7 @@
 """
 module moment_kinetics
 
-export run_moment_kinetics, restart_moment_kinetics
+export run_moment_kinetics
 
 using MPI
 
@@ -92,11 +92,13 @@ using .type_definitions: mk_int
 """
 main function that contains all of the content of the program
 """
-function run_moment_kinetics(to::TimerOutput, input_dict=Dict())
+function run_moment_kinetics(to::TimerOutput, input_dict=Dict(); restart=false,
+                             restart_time_index=-1)
     mk_state = nothing
     try
         # set up all the structs, etc. needed for a run
-        mk_state = setup_moment_kinetics(input_dict)
+        mk_state = setup_moment_kinetics(input_dict; restart=restart,
+                                         restart_time_index=restart_time_index)
 
         # solve the 1+1D kinetic equation to advance f in time by nstep time steps
         if run_type == performance_test
@@ -139,27 +141,38 @@ end
 """
 overload which takes a filename and loads input
 """
-function run_moment_kinetics(to::TimerOutput, input_filename::String)
-    return run_moment_kinetics(to, read_input_file(input_filename))
+function run_moment_kinetics(to::TimerOutput, input_filename::String; restart=false,
+                             restart_time_index=-1)
+    return run_moment_kinetics(to, read_input_file(input_filename); restart=restart,
+                               restart_time_index=restart_time_index)
 end
 
 """
 overload with no TimerOutput arguments
 """
-function run_moment_kinetics(input)
-    return run_moment_kinetics(TimerOutput(), input)
+function run_moment_kinetics(input; restart=false, restart_time_index=-1)
+    return run_moment_kinetics(TimerOutput(), input; restart=restart,
+                               restart_time_index=restart_time_index)
 end
 
 """
 overload which gets the input file name from command line arguments
 """
 function run_moment_kinetics()
-    inputfile = get_options()["inputfile"]
-    if inputfile == nothing
-        run_moment_kinetics(Dict())
+    options = get_options()
+    inputfile = options["inputfile"]
+    restart = options["restart"]
+    if options["restartfile"] !== nothing
+        restart = options["restartfile"]
+    end
+    restart_time_index = options["restart-time-index"]
+    if inputfile === nothing
+        this_input = Dict()
     else
-        run_moment_kinetics(inputfile)
+        this_input = inputfile
     end
+    run_moment_kinetics(this_input; restart=restart,
+                        restart_time_index=restart_time_index)
 end
 
 """
@@ -226,111 +239,9 @@ function get_backup_filename(filename)
     end
     backup_dfns_filename == "" && error("Failed to find a name for backup file.")
     backup_prefix_iblock = ("$(basename)_$(counter)", iblock)
+    original_prefix_iblock = (basename, iblock)
     return dfns_filename, backup_dfns_filename, parallel_io, moments_filename,
-           backup_moments_filename, backup_prefix_iblock
-end
-
-"""
-    restart_moment_kinetics(input_filename::String,
-                            restart_filename::Union{String,Nothing}=nothing,
-                            time_index::Int=-1)
-
-Restart moment kinetics from an existing run. Space/velocity-space resolution in the
-input must be the same as for the original run.
-
-`input_filename` is the input file to use.
-
-`restart_filename` can be used to pick a particular distribution-functions-output file to
-restart from. By default will use the most recent one (the one without the numerical
-suffix) in the run directory.
-
-`time_index` can be passed to select the time index from `restart_filename` to restart
-from. By default the latest time point is used.
-"""
-function restart_moment_kinetics(input_filename::String,
-                                 restart_filename::Union{String,Nothing}=nothing,
-                                 time_index::Int=-1)
-    restart_moment_kinetics(read_input_file(input_filename), restart_filename,
-                            time_index)
-    return nothing
-end
-function restart_moment_kinetics()
-    options = get_options()
-    inputfile = options["inputfile"]
-    if inputfile === nothing
-        error("Must pass input file as first argument to restart a run.")
-    end
-    restartfile = options["restartfile"]
-    if restartfile === nothing
-        error("Must pass output file to restart from as second argument.")
-    end
-    time_index = options["restart-time-index"]
-
-    restart_moment_kinetics(inputfile, restartfile, time_index)
-
-    return nothing
-end
-function restart_moment_kinetics(input_dict::Dict,
-                                 restart_filename::Union{String,Nothing}=nothing,
-                                 time_index::Int=-1)
-
-    if restart_filename === nothing
-        run_name = input_dict["run_name"]
-        base_directory = get(input_dict, "base_directory", "runs")
-        output_dir = joinpath(base_directory, run_name)
-        io_settings = get(input_dict, "output", Dict{String,Any}())
-        binary_format = get(io_settings, "binary_format", hdf5)
-        if binary_format === hdf5
-            ext = "h5"
-        elseif binary_format === netcdf
-            ext = "cdf"
-        else
-            error("Unrecognized binary_format '$binary_format'")
-        end
-        restart_filename = glob(joinpath(output_dir, run_name * ".dfns*." * ext))[1]
-    end
-
-    try
-        # Move the output file being restarted from to make sure it doesn't get
-        # overwritten.
-        dfns_filename, backup_dfns_filename, parallel_io, moments_filename,
-        backup_moments_filename, backup_prefix_iblock =
-            get_backup_filename(restart_filename)
-        # Ensure every process got the filenames and checked files exist before moving
-        # files
-        MPI.Barrier(comm_world)
-        if (parallel_io && global_rank[] == 0) || (!parallel_io && block_rank[] == 0)
-            mv(dfns_filename, backup_dfns_filename)
-            mv(moments_filename, backup_moments_filename)
-        end
-        # Ensure files have been moved before any process tries to read from them
-        MPI.Barrier(comm_world)
-
-        # Set up all the structs, etc. needed for a run.
-        mk_state = setup_moment_kinetics(input_dict,
-                                         restart_prefix_iblock=backup_prefix_iblock,
-                                         restart_time_index=time_index)
-
-        try
-            time_advance!(mk_state...)
-        finally
-            # clean up i/o and communications
-            # last 2 elements of mk_state are `io` and `cdf`
-            cleanup_moment_kinetics!(mk_state[end-2:end]...)
-        end
-    catch e
-        # Stop code from hanging when running on multiple processes if only one of them
-        # throws an error
-        if global_size[] > 1
-            println("Abort called on rank $(block_rank[]) due to error. Error message "
-                    * "was:\n", e)
-            MPI.Abort(comm_world, 1)
-        end
-
-        rethrow(e)
-    end
-
-    return nothing
+           backup_moments_filename, backup_prefix_iblock, original_prefix_iblock
 end
 
 """
@@ -342,14 +253,19 @@ reload data from time index given by `restart_time_index` for a restart.
 `debug_loop_type` and `debug_loop_parallel_dims` are used to force specific set ups for
 parallel loop ranges, and are only used by the tests in `debug_test/`.
 """
-function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing,
-        restart_time_index=-1,
+function setup_moment_kinetics(input_dict::Dict;
+        restart::Union{Bool,AbstractString}=false, restart_time_index::mk_int=-1,
         debug_loop_type::Union{Nothing,NTuple{N,Symbol} where N}=nothing,
         debug_loop_parallel_dims::Union{Nothing,NTuple{N,Symbol} where N}=nothing)
 
     # Set up MPI
     initialize_comms!()
 
+    if global_rank[] == 0
+        println("Starting setup   ", Dates.format(now(), dateformat"H:MM:SS"))
+        flush(stdout)
+    end
+
     input = mk_input(input_dict; save_inputs_to_txt=true, ignore_MPI=false)
     # obtain input options from moment_kinetics_input.jl
     # and check input to catch errors
@@ -383,7 +299,7 @@ function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing,
         allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz,
                                  evolve_moments, collisions, num_diss_params)
 
-    if restart_prefix_iblock === nothing
+    if restart === false
         restarting = false
         # initialize f(z,vpa) and the lowest three v-space moments (density(z), upar(z) and ppar(z)),
         # each of which may be evolved separately depending on input choices.
@@ -397,10 +313,60 @@ function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing,
     else
         restarting = true
 
+        run_name = input_dict["run_name"]
+        base_directory = get(input_dict, "base_directory", "runs")
+        output_dir = joinpath(base_directory, run_name)
+        if restart === true
+            run_name = input_dict["run_name"]
+            io_settings = get(input_dict, "output", Dict{String,Any}())
+            binary_format = get(io_settings, "binary_format", hdf5)
+            if binary_format === hdf5
+                ext = "h5"
+            elseif binary_format === netcdf
+                ext = "cdf"
+            else
+                error("Unrecognized binary_format '$binary_format'")
+            end
+            restart_filename_pattern = joinpath(output_dir, run_name * ".dfns*." * ext)
+            restart_filename_glob = glob(restart_filename_pattern)
+            if length(restart_filename_glob) == 0
+                error("No output file to restart from found matching the pattern "
+                      * "$restart_filename_pattern")
+            end
+            restart_filename = restart_filename_glob[1]
+        else
+            restart_filename = restart
+        end
+
+        # Move the output file being restarted from to make sure it doesn't get
+        # overwritten.
+        dfns_filename, backup_dfns_filename, parallel_io, moments_filename,
+        backup_moments_filename, backup_prefix_iblock, original_prefix_iblock =
+            get_backup_filename(restart_filename)
+
+        # Ensure every process got the filenames and checked files exist before moving
+        # files
+        MPI.Barrier(comm_world)
+
+        if abspath(output_dir) == abspath(dirname(dfns_filename))
+            # Only move the file if it is in our current run directory. Otherwise we are
+            # restarting from another run, and will not be overwriting the file.
+            if (parallel_io && global_rank[] == 0) || (!parallel_io && block_rank[] == 0)
+                mv(dfns_filename, backup_dfns_filename)
+                mv(moments_filename, backup_moments_filename)
+            end
+        else
+            # Reload from dfns_filename without moving the file
+            backup_prefix_iblock = original_prefix_iblock
+        end
+
+        # Ensure files have been moved before any process tries to read from them
+        MPI.Barrier(comm_world)
+
         # Reload pdf and moments from an existing output file
         code_time, previous_runs_info, restart_time_index =
             reload_evolving_fields!(pdf, moments, boundary_distributions,
-                                    restart_prefix_iblock, restart_time_index,
+                                    backup_prefix_iblock, restart_time_index,
                                     composition, r, z, vpa, vperp, vzeta, vr, vz)
         _block_synchronize()
     end
diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl
index b194b8be8..98b5d1a0b 100644
--- a/src/moment_kinetics_input.jl
+++ b/src/moment_kinetics_input.jl
@@ -230,7 +230,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # determine the boundary condition to impose in r
     # supported options are "periodic" and "Dirichlet"
     r.bc = get(scan_input, "r_bc", "periodic")
-
+    r.element_spacing_option = get(scan_input, "r_element_spacing_option", "uniform")
+    
     # overwrite some default parameters related to the z grid
     # ngrid is number of grid points per element
     z.ngrid = get(scan_input, "z_ngrid", 9)
@@ -247,7 +248,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # determine the boundary condition to impose in z
     # supported options are "constant", "periodic" and "wall"
     z.bc = get(scan_input, "z_bc", "wall")
-
+    z.element_spacing_option = get(scan_input, "z_element_spacing_option", "uniform")
+    
     # overwrite some default parameters related to the vpa grid
     # ngrid is the number of grid points per element
     vpa.ngrid = get(scan_input, "vpa_ngrid", 17)
@@ -264,7 +266,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # supported options are "chebyshev_pseudospectral" and "finite_difference"
     vpa.discretization = get(scan_input, "vpa_discretization", "chebyshev_pseudospectral")
     vpa.fd_option = get(scan_input, "vpa_finite_difference_option", "third_order_upwind")
-
+    vpa.element_spacing_option = get(scan_input, "vpa_element_spacing_option", "uniform")
+    
     num_diss_params = setup_numerical_dissipation(
         get(scan_input, "numerical_dissipation", Dict{String,Any}()), true)
 
@@ -285,7 +288,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # determine the discretization option for the vperp grid
     # supported options are "finite_difference_vperp" "chebyshev_pseudospectral_vperp"
     vperp.discretization = get(scan_input, "vperp_discretization", "chebyshev_pseudospectral_vperp")
-
+    vperp.element_spacing_option = get(scan_input, "vperp_element_spacing_option", "uniform")
+    
     # overwrite some default parameters related to the gyrophase grid
     # ngrid is the number of grid points per element
     gyrophase.ngrid = get(scan_input, "gyrophase_ngrid", 17)
@@ -311,7 +315,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
         # determine the discretization option for the vz grid
         # supported options are "chebyshev_pseudospectral" and "finite_difference"
         vz.discretization = get(scan_input, "vz_discretization", vpa.discretization)
-
+        vz.element_spacing_option = get(scan_input, "vz_element_spacing_option", "uniform")
+    
         # overwrite some default parameters related to the vr grid
         # ngrid is the number of grid points per element
         vr.ngrid = get(scan_input, "vr_ngrid", 1)
@@ -327,7 +332,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
         # determine the discretization option for the vr grid
         # supported options are "chebyshev_pseudospectral" and "finite_difference"
         vr.discretization = get(scan_input, "vr_discretization", "chebyshev_pseudospectral")
-
+        vr.element_spacing_option = get(scan_input, "vr_element_spacing_option", "uniform")
+    
         # overwrite some default parameters related to the vzeta grid
         # ngrid is the number of grid points per element
         vzeta.ngrid = get(scan_input, "vzeta_ngrid", 1)
@@ -343,6 +349,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
         # determine the discretization option for the vzeta grid
         # supported options are "chebyshev_pseudospectral" and "finite_difference"
         vzeta.discretization = get(scan_input, "vzeta_discretization", "chebyshev_pseudospectral")
+        vzeta.element_spacing_option = get(scan_input, "vzeta_element_spacing_option", "uniform")
     end
 
     is_1V = (vperp.ngrid == vperp.nelement_global == 1 && vzeta.ngrid ==
@@ -379,37 +386,37 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
 	z_advection_immutable = advection_input(z.advection.option, z.advection.constant_speed,
         z.advection.frequency, z.advection.oscillation_amplitude)
     z_immutable = grid_input("z", z.ngrid, z.nelement_global, z.nelement_local, nrank_z, irank_z, z.L, 
-        z.discretization, z.fd_option, z.bc, z_advection_immutable, comm_sub_z)
+        z.discretization, z.fd_option, z.bc, z_advection_immutable, comm_sub_z, z.element_spacing_option)
     r_advection_immutable = advection_input(r.advection.option, r.advection.constant_speed,
         r.advection.frequency, r.advection.oscillation_amplitude)
     r_immutable = grid_input("r", r.ngrid, r.nelement_global, r.nelement_local, nrank_r, irank_r, r.L,
-        r.discretization, r.fd_option, r.bc, r_advection_immutable, comm_sub_r)
+        r.discretization, r.fd_option, r.bc, r_advection_immutable, comm_sub_r, r.element_spacing_option)
 	# for dimensions below which do not currently use distributed-memory MPI
 	# assign dummy values to nrank, irank and comm of coord struct
     vpa_advection_immutable = advection_input(vpa.advection.option, vpa.advection.constant_speed,
         vpa.advection.frequency, vpa.advection.oscillation_amplitude)
     vpa_immutable = grid_input("vpa", vpa.ngrid, vpa.nelement_global, vpa.nelement_local, 1, 0, vpa.L,
-        vpa.discretization, vpa.fd_option, vpa.bc, vpa_advection_immutable, MPI.COMM_NULL)
+        vpa.discretization, vpa.fd_option, vpa.bc, vpa_advection_immutable, MPI.COMM_NULL, vpa.element_spacing_option)
     vperp_advection_immutable = advection_input(vperp.advection.option, vperp.advection.constant_speed,
         vperp.advection.frequency, vperp.advection.oscillation_amplitude)
     vperp_immutable = grid_input("vperp", vperp.ngrid, vperp.nelement_global, vperp.nelement_local, 1, 0, vperp.L,
-        vperp.discretization, vperp.fd_option, vperp.bc, vperp_advection_immutable, MPI.COMM_NULL)
+        vperp.discretization, vperp.fd_option, vperp.bc, vperp_advection_immutable, MPI.COMM_NULL, vperp.element_spacing_option)
     gyrophase_advection_immutable = advection_input(gyrophase.advection.option, gyrophase.advection.constant_speed,
         gyrophase.advection.frequency, gyrophase.advection.oscillation_amplitude)
     gyrophase_immutable = grid_input("gyrophase", gyrophase.ngrid, gyrophase.nelement_global, gyrophase.nelement_local, 1, 0, gyrophase.L,
-        gyrophase.discretization, gyrophase.fd_option, gyrophase.bc, gyrophase_advection_immutable, MPI.COMM_NULL)
+        gyrophase.discretization, gyrophase.fd_option, gyrophase.bc, gyrophase_advection_immutable, MPI.COMM_NULL, gyrophase.element_spacing_option)
     vz_advection_immutable = advection_input(vz.advection.option, vz.advection.constant_speed,
         vz.advection.frequency, vz.advection.oscillation_amplitude)
     vz_immutable = grid_input("vz", vz.ngrid, vz.nelement_global, vz.nelement_local, 1, 0, vz.L,
-        vz.discretization, vz.fd_option, vz.bc, vz_advection_immutable, MPI.COMM_NULL)
+        vz.discretization, vz.fd_option, vz.bc, vz_advection_immutable, MPI.COMM_NULL, vz.element_spacing_option)
     vr_advection_immutable = advection_input(vr.advection.option, vr.advection.constant_speed,
         vr.advection.frequency, vr.advection.oscillation_amplitude)
     vr_immutable = grid_input("vr", vr.ngrid, vr.nelement_global, vr.nelement_local, 1, 0, vr.L,
-        vr.discretization, vr.fd_option, vr.bc, vr_advection_immutable, MPI.COMM_NULL)
+        vr.discretization, vr.fd_option, vr.bc, vr_advection_immutable, MPI.COMM_NULL, vr.element_spacing_option)
     vzeta_advection_immutable = advection_input(vzeta.advection.option, vzeta.advection.constant_speed,
         vzeta.advection.frequency, vzeta.advection.oscillation_amplitude)
     vzeta_immutable = grid_input("vzeta", vzeta.ngrid, vzeta.nelement_global, vzeta.nelement_local, 1, 0, vzeta.L,
-        vzeta.discretization, vzeta.fd_option, vzeta.bc, vzeta_advection_immutable, MPI.COMM_NULL)
+        vzeta.discretization, vzeta.fd_option, vzeta.bc, vzeta_advection_immutable, MPI.COMM_NULL, vzeta.element_spacing_option)
     
     species_charged_immutable = Array{species_parameters,1}(undef,n_ion_species)
     species_neutral_immutable = Array{species_parameters,1}(undef,n_neutral_species)
@@ -569,10 +576,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_z = advection_input_mutable(advection_option_z, advection_speed_z,
         frequency_z, oscillation_amplitude_z)
+    element_spacing_option_z = "uniform"
     # create a mutable structure containing the input info related to the z grid
     z = grid_input_mutable("z", ngrid_z, nelement_global_z, nelement_local_z, L_z,
         discretization_option_z, finite_difference_option_z, boundary_option_z,
-        advection_z)
+        advection_z, element_spacing_option_z)
     #################### parameters related to the r grid ######################
     # ngrid_r is number of grid points per element
     ngrid_r = 1
@@ -609,10 +617,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for r
     advection_r = advection_input_mutable(advection_option_r, advection_speed_r,
         frequency_r, oscillation_amplitude_r)
+    element_spacing_option_r = "uniform"
     # create a mutable structure containing the input info related to the r grid
     r = grid_input_mutable("r", ngrid_r, nelement_global_r, nelement_local_r, L_r,
         discretization_option_r, finite_difference_option_r, boundary_option_r,
-        advection_r)
+        advection_r, element_spacing_option_r)
     ############################################################################
     ################### parameters related to the vpa grid #####################
     # ngrid_vpa is the number of grid points per element
@@ -647,10 +656,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vpa = advection_input_mutable(advection_option_vpa, advection_speed_vpa,
         frequency_vpa, oscillation_amplitude_vpa)
+    element_spacing_option_vpa = "uniform"
     # create a mutable structure containing the input info related to the vpa grid
     vpa = grid_input_mutable("vpa", ngrid_vpa, nelement_vpa, nelement_vpa, L_vpa,
         discretization_option_vpa, finite_difference_option_vpa, boundary_option_vpa,
-        advection_vpa)
+        advection_vpa, element_spacing_option_vpa)
     ############################################################################
     ################### parameters related to the vperp grid #####################
     # ngrid_vperp is the number of grid points per element
@@ -684,10 +694,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vperp = advection_input_mutable(advection_option_vperp, advection_speed_vperp,
         frequency_vperp, oscillation_amplitude_vperp)
+    element_spacing_option_vperp = "uniform"
     # create a mutable structure containing the input info related to the vperp grid
     vperp = grid_input_mutable("vperp", ngrid_vperp, nelement_vperp, nelement_vperp, L_vperp,
         discretization_option_vperp, finite_difference_option_vperp, boundary_option_vperp,
-        advection_vperp)
+        advection_vperp, element_spacing_option_vperp)
     ############################################################################
     ################### parameters related to the gyrophase grid #####################
     # ngrid_gyrophase is the number of grid points per element
@@ -707,10 +718,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     oscillation_amplitude_gyrophase = 1.0
     advection_gyrophase = advection_input_mutable(advection_option_gyrophase, advection_speed_gyrophase,
         frequency_gyrophase, oscillation_amplitude_gyrophase)
+    element_spacing_option_gyrophase = "uniform"
     # create a mutable structure containing the input info related to the gyrophase grid
     gyrophase = grid_input_mutable("gyrophase", ngrid_gyrophase, nelement_gyrophase, nelement_gyrophase, L_gyrophase,
         discretization_option_gyrophase, finite_difference_option_gyrophase, boundary_option_gyrophase,
-        advection_gyrophase)
+        advection_gyrophase, element_spacing_option_gyrophase)
     ############################################################################
     ################### parameters related to the vr grid #####################
     # ngrid_vr is the number of grid points per element
@@ -742,10 +754,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vr = advection_input_mutable(advection_option_vr, advection_speed_vr,
         frequency_vr, oscillation_amplitude_vr)
+    element_spacing_option_vr = "uniform"
     # create a mutable structure containing the input info related to the vr grid
     vr = grid_input_mutable("vr", ngrid_vr, nelement_vr, nelement_vr, L_vr,
         discretization_option_vr, finite_difference_option_vr, boundary_option_vr,
-        advection_vr)
+        advection_vr, element_spacing_option_vr)
     ############################################################################
     ################### parameters related to the vz grid #####################
     # ngrid_vz is the number of grid points per element
@@ -777,10 +790,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vz = advection_input_mutable(advection_option_vz, advection_speed_vz,
         frequency_vz, oscillation_amplitude_vz)
+    element_spacing_option_vz = "uniform"
     # create a mutable structure containing the input info related to the vz grid
     vz = grid_input_mutable("vz", ngrid_vz, nelement_vz, nelement_vz, L_vz,
         discretization_option_vz, finite_difference_option_vz, boundary_option_vz,
-        advection_vz)
+        advection_vz, element_spacing_option_vz)
     ############################################################################
     ################### parameters related to the vzeta grid #####################
     # ngrid_vzeta is the number of grid points per element
@@ -812,10 +826,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vzeta = advection_input_mutable(advection_option_vzeta, advection_speed_vzeta,
         frequency_vzeta, oscillation_amplitude_vzeta)
+    element_spacing_option_vzeta = "uniform"
     # create a mutable structure containing the input info related to the vzeta grid
     vzeta = grid_input_mutable("vzeta", ngrid_vzeta, nelement_vzeta, nelement_vzeta, L_vzeta,
         discretization_option_vzeta, finite_difference_option_vzeta, boundary_option_vzeta,
-        advection_vzeta)
+        advection_vzeta, element_spacing_option_vzeta)
     #############################################################################
     # define default values and create corresponding mutable structs holding
     # information about the composition of the species and their initial conditions
diff --git a/src/post_processing.jl b/src/post_processing.jl
index 3e1a7f1c8..b3d971f91 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -199,7 +199,7 @@ function construct_global_zr_coords(r_local, z_local)
         return grid_input(coord_local.name, coord_local.ngrid,
             coord_local.nelement_global, coord_local.nelement_global, 1, 0, coord_local.L,
             coord_local.discretization, coord_local.fd_option, coord_local.bc,
-            coord_local.advection, MPI.COMM_NULL)
+            coord_local.advection, MPI.COMM_NULL, coord_local.element_spacing_option)
     end
 
     r_global, r_global_spectral = define_coordinate(make_global_input(r_local))
@@ -373,12 +373,12 @@ are passed, the plots/movies are given names beginning with `compare_` and are c
 in the `comparison_plots/` subdirectory.
 
 By default plots output from all restarts in a directory. To select a single run, pass the
-`run_index` argument - the value corresponds to the `_<i>` suffix given to output files by
-`restart_moment_kinetics()`. `run_index` can be an integer (which is applied to all
-directories in `prefix...`), or a tuple of integers (which should have the same length as
-the number of directories passed to `prefix...`). Use `run_index=-1` to get the most
-recent run (which does not have a `_<i>` suffix). Note that `run_index` is only used when
-a directory (rather than the prefix of a specific output file) is passed to `prefix...`
+`run_index` argument - the value corresponds to the `_<i>` suffix given to output files
+when restarting. `run_index` can be an integer (which is applied to all directories in
+`prefix...`), or a tuple of integers (which should have the same length as the number of
+directories passed to `prefix...`). Use `run_index=-1` to get the most recent run (which
+does not have a `_<i>` suffix). Note that `run_index` is only used when a directory
+(rather than the prefix of a specific output file) is passed to `prefix...`
 """
 function analyze_and_plot_data(prefix...; run_index=nothing)
     if length(prefix) == 0
diff --git a/src/time_advance.jl b/src/time_advance.jl
index 4892dd87e..ae9ef4f28 100644
--- a/src/time_advance.jl
+++ b/src/time_advance.jl
@@ -8,8 +8,8 @@ export time_advance!
 using MPI
 using ..type_definitions: mk_float
 using ..array_allocation: allocate_float, allocate_shared_float
-using ..communication: _block_synchronize, block_rank, global_size, iblock_index,
-                       comm_world, MPISharedArray, global_rank
+using ..communication
+using ..communication: _block_synchronize
 using ..debugging
 using ..file_io: write_data_to_ascii, write_moments_data_to_binary, write_dfns_data_to_binary, debug_dump
 using ..looping
@@ -751,6 +751,7 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
     @serial_region begin
         if global_rank[] == 0
              println("beginning time advance   ", Dates.format(now(), dateformat"H:MM:SS"))
+             flush(stdout)
         end
     end
 
@@ -786,6 +787,26 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
             if isfile(t_input.stopfile)
                 # Stop cleanly if a file called 'stop' was created
                 println("Found 'stop' file $(t_input.stopfile), aborting run")
+                flush(stdout)
+                finish_now = true
+            end
+
+            # If NaNs are present, they should propagate into every field, so only need to
+            # check one. Choose phi because it is small (it has no species or velocity
+            # dimensions). If a NaN is found, stop the simulation.
+            if block_rank[] == 0
+                if any(isnan.(fields.phi))
+                    println("Found NaN, stopping simulation")
+                    found_nan = 1
+                else
+                    found_nan = 0
+                end
+                found_nan = MPI.Allreduce(found_nan, +, comm_inter_block[])
+            else
+                found_nan = 0
+            end
+            found_nan = MPI.Bcast(found_nan, 0, comm_block[])
+            if found_nan != 0
                 finish_now = true
             end
         end
@@ -844,12 +865,14 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
                 end
                 if global_rank[] == 0
                     println("    residuals:", result_string)
+                    flush(stdout)
                 end
                 if t_input.converged_residual_value > 0.0
                     if global_rank[] == 0
                         if all(r < t_input.converged_residual_value for r ∈ all_residuals)
                             println("Run converged! All tested residuals less than ",
                                     t_input.converged_residual_value)
+                            flush(stdout)
                             finish_now = true
                         end
                     end
@@ -858,6 +881,7 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
             else
                 if global_rank[] == 0
                     println()
+                    flush(stdout)
                 end
             end
 
@@ -975,6 +999,7 @@ function time_advance!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyro
                 if global_rank[] == 0
                     println("writing distribution functions at step ", i,"  ",
                                    Dates.format(now(), dateformat"H:MM:SS"))
+                    flush(stdout)
                 end
             end
             write_dfns_data_to_binary(pdf.charged.norm, pdf.neutral.norm, moments, fields,
diff --git a/submit-restart.sh b/submit-restart.sh
index 759c247c4..3186173a4 100755
--- a/submit-restart.sh
+++ b/submit-restart.sh
@@ -94,19 +94,10 @@ RUNNAME=$(util/get-run-name.jl $INPUTFILE)
 RUNDIR=runs/$RUNNAME/
 mkdir -p $RUNDIR
 
-# Get default file to restart from, which is the latest run in $RUNDIR
-if [[ -z $RESTARTFROM ]]; then
-  # "shopt -s extglob" is needed to let us use the ?() syntax within a script
-  # (it doesn't seem to be needed in an interactive shell!). See
-  # https://www.linuxjournal.com/content/pattern-matching-bash
-  shopt -s extglob
-  RESTARTFROM=$(ls $RUNDIR/$RUNNAME.dfns*.?(h5|cdf) | head -n 1)
-fi
-
 if [[ $POSTPROC -eq 0 ]]; then
-  echo "Submitting $INPUTFILE for restart from $RESTARTFROM and post-processing..."
+  echo "Submitting $INPUTFILE for restart from '$RESTARTFROM' and post-processing..."
 else
-  echo "Submitting $INPUTFILE for restart from $RESTARTFROM..."
+  echo "Submitting $INPUTFILE for restart from '$RESTARTFROM'..."
 fi
 
 # Create a submission script for the run
diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl
index c1bbdd483..ae198003c 100644
--- a/test/calculus_tests.jl
+++ b/test/calculus_tests.jl
@@ -14,7 +14,7 @@ function runtests()
         println("calculus tests")
         @testset "fundamental theorem of calculus" begin
             @testset "$discretization $ngrid $nelement" for
-                    discretization ∈ ("finite_difference", "chebyshev_pseudospectral"),
+                    (discretization, element_spacing_option, etol) ∈ (("finite_difference", "uniform", 1.0e-15), ("chebyshev_pseudospectral", "uniform", 1.0e-15), ("chebyshev_pseudospectral", "sqrt", 1.0e-2)),
                     ngrid ∈ (5,6,7,8,9,10), nelement ∈ (1, 2, 3, 4, 5)
 
                 if discretization == "finite_difference" && (ngrid - 1) * nelement % 2 == 1
@@ -26,7 +26,6 @@ function runtests()
                     continue
                 end
 
-                etol = 1.0e-15
                 # define inputs needed for the test
                 L = 6.0
                 bc = "periodic"
@@ -41,7 +40,8 @@ function runtests()
 				comm = MPI.COMM_NULL # dummy value 
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    discretization, fd_option, bc, adv_input, comm)
+                    discretization, fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
                 # create array for the function f(x) to be differentiated/integrated
@@ -86,9 +86,11 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value 
-				input = grid_input("coord", ngrid, nelement,
+				element_spacing_option = "uniform" # dummy value
+                input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -132,10 +134,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -175,10 +179,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -226,10 +232,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -440,10 +448,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -634,10 +644,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -662,7 +674,7 @@ function runtests()
         end
 
         @testset "Chebyshev pseudospectral derivatives (4 argument), Neumann" verbose=false begin
-            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"),
+            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
                     nelement ∈ (1:5), ngrid ∈ (3:33)
 
                 # define inputs needed for the test
@@ -676,10 +688,11 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
-				input = grid_input("coord", ngrid, nelement,
+				comm = MPI.COMM_NULL # dummy value
+                input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -712,7 +725,7 @@ function runtests()
         end
 
         @testset "Chebyshev pseudospectral derivatives upwinding (5 argument), Neumann" verbose=false begin
-            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"),
+            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
                     nelement ∈ (1:5), ngrid ∈ (3:33)
 
                 # define inputs needed for the test
@@ -726,10 +739,11 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
-				input = grid_input("coord", ngrid, nelement,
+				comm = MPI.COMM_NULL # dummy value
+                input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -940,9 +954,11 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index b40a55ff7..8885224cc 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -27,18 +27,21 @@ adv_input = advection_input("default", 1.0, 0.0, 0.0)
 function runtests()
     @testset "interpolation" verbose=use_verbose begin
         @testset "$discretization, $ntest, $nelement, $zlim" for
-                (discretization, rtol) ∈
-                    (("finite_difference", 1.e-5), ("chebyshev_pseudospectral", 1.e-8)),
+                (discretization, element_spacing_option, rtol) ∈
+                    (("finite_difference", "uniform", 1.e-5), ("chebyshev_pseudospectral", "uniform", 1.e-8),
+                    ("chebyshev_pseudospectral", "sqrt", 1.e-8)),
                     ntest ∈ (3, 14), nelement ∈ (2, 8), zlim ∈ (L/2.0, L/5.0)
 
             # create the 'input' struct containing input info needed to create a coordinate
 			nelement_local = nelement
 			nrank_per_block = 0 # dummy value
 			irank = 0 # dummy value
-			comm = MPI.COMM_NULL # dummy value 
+			comm = MPI.COMM_NULL # dummy value
+            #element_spacing_option = "uniform"
 			input = grid_input("coord", ngrid, nelement,
 				nelement_local, nrank_per_block, irank, L,
-				discretization, fd_option, bc, adv_input, comm)
+				discretization, fd_option, bc, adv_input, comm,
+                element_spacing_option)
             # create the coordinate struct 'z'
             z, spectral = define_coordinate(input)
 
diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl
index 9a5708921..c0b563393 100644
--- a/test/wall_bc_tests.jl
+++ b/test/wall_bc_tests.jl
@@ -82,6 +82,7 @@ test_input_finite_difference = Dict("n_ion_species" => 1,
                                     "z_nelement" => 1,
                                     "z_bc" => "wall",
                                     "z_discretization" => "finite_difference",
+                                    "z_element_spacing_option" => "uniform",
                                     "vpa_ngrid" => 400,
                                     "vpa_nelement" => 1,
                                     "vpa_L" => 8.0,
@@ -98,6 +99,32 @@ test_input_chebyshev = merge(test_input_finite_difference,
                                   "z_discretization" => "chebyshev_pseudospectral",
                                   "z_ngrid" => 9,
                                   "z_nelement" => 2,
+                                  "z_element_spacing_option" => "uniform",
+                                  "vpa_discretization" => "chebyshev_pseudospectral",
+                                  "vpa_ngrid" => 17,
+                                  "vpa_nelement" => 10,
+                                  "vz_discretization" => "chebyshev_pseudospectral",
+                                  "vz_ngrid" => 17,
+                                  "vz_nelement" => 10))
+                                  
+test_input_chebyshev_sqrt_grid_odd = merge(test_input_finite_difference,
+                             Dict("run_name" => "chebyshev_pseudospectral",
+                                  "z_discretization" => "chebyshev_pseudospectral",
+                                  "z_ngrid" => 9,
+                                  "z_nelement" => 5, # minimum nontrival nelement (odd)
+                                  "z_element_spacing_option" => "sqrt",
+                                  "vpa_discretization" => "chebyshev_pseudospectral",
+                                  "vpa_ngrid" => 17,
+                                  "vpa_nelement" => 10,
+                                  "vz_discretization" => "chebyshev_pseudospectral",
+                                  "vz_ngrid" => 17,
+                                  "vz_nelement" => 10))
+test_input_chebyshev_sqrt_grid_even = merge(test_input_finite_difference,
+                             Dict("run_name" => "chebyshev_pseudospectral",
+                                  "z_discretization" => "chebyshev_pseudospectral",
+                                  "z_ngrid" => 9,
+                                  "z_nelement" => 6, # minimum nontrival nelement (even)
+                                  "z_element_spacing_option" => "sqrt",
                                   "vpa_discretization" => "chebyshev_pseudospectral",
                                   "vpa_ngrid" => 17,
                                   "vpa_nelement" => 10,
@@ -125,7 +152,7 @@ function run_test(test_input, expected_phi, tolerance; args...)
     # update the default inputs
 
     # Convert keyword arguments to a unique name
-    name = test_input["run_name"]
+    name = test_input["run_name"] * ", with element spacing: " * test_input["z_element_spacing_option"]
     if length(args) > 0
         name = string(name, "_", (string(k, "-", v, "_") for (k, v) in args)...)
 
@@ -191,16 +218,22 @@ function run_test(test_input, expected_phi, tolerance; args...)
         adv_input = advection_input("default", 1.0, 0.0, 0.0)
 		nrank_per_block = 0 # dummy value
 		irank = 0 # dummy value
-		comm = MPI.COMM_NULL # dummy value 
+		comm = MPI.COMM_NULL # dummy value
+        element_spacing_option = "uniform"
         input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], 
 						   test_input["z_nelement"], nrank_per_block, irank, 1.0,
                            test_input["z_discretization"], "", test_input["z_bc"],
-                           adv_input, comm)
+                           adv_input, comm, test_input["z_element_spacing_option"])
         z, z_spectral = define_coordinate(input)
 
         # Cross comparison of all discretizations to same benchmark
-        phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral)
-        @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15)
+        if test_input["z_element_spacing_option"] == "uniform" 
+            # Only support this test for uniform element spacing.
+            # phi is better resolved by "sqrt" spacing grid, so disagrees with benchmark data from
+            # simulation with uniform element spacing.
+            phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral)
+            @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15)
+        end
     end
 end
 
@@ -213,12 +246,32 @@ function runtests()
             run_test(test_input_finite_difference, nothing, 2.e-3)
         end
 
-        @testset "Chebyshev" begin
+        @testset "Chebyshev uniform" begin
             run_test(test_input_chebyshev,
                      [-1.1689445031600718, -0.7479504438063098, -0.6947559936893813,
                       -0.6917252442591313, -0.7180152498764835, -0.9980114095597415],
                      2.e-3)
         end
+        
+        @testset "Chebyshev sqrt grid odd" begin
+            run_test(test_input_chebyshev_sqrt_grid_odd,
+                     [-1.2047298885671576, -0.9431378294506091, -0.8084332392927167,
+                     -0.7812620422650213, -0.7233303514000929, -0.7003878610612269,
+                     -0.69572751349158, -0.6933148921301019, -0.6992503992521327,
+                     -0.7115787972775218, -0.7596015032228407, -0.795776514029509,
+                     -0.876303297135126, -1.1471244425913258],
+                     2.e-3)
+        end
+        @testset "Chebyshev sqrt grid even" begin
+            run_test(test_input_chebyshev_sqrt_grid_even,
+                     [-1.213617049279473, -1.0054529928344382, -0.871444761913497,
+                     -0.836017699317097, -0.7552110924643832, -0.7264644073096705,
+                     -0.7149147366621806, -0.6950077192395091, -0.6923364889119271,
+                     -0.6950077192395089, -0.7149147366621814, -0.7264644073096692,
+                     -0.7552110924643836, -0.8360176993170979, -0.8714447619134948,
+                     -1.0054529928344376, -1.2136170492794727],
+                     2.e-3)
+        end
     end
 end