#include "unused_dummy.inc"
!> FIXME : Add documentation
module gs2_save
  use constants, only: run_name_size

  implicit none

  private

  public :: gs2_restore, gs2_save_for_restart, finish_gs2_save
  public :: read_many, save_many, gs2_save_response, gs2_restore_response
  public :: init_gs2_save, init_dt, read_t0_from_restart_file, init_ant_amp
  public :: set_restart_file, include_explicit_source_in_restart, proc_to_save_fields
  public :: init_vnm, restart_writable, EigNetcdfID
  public :: init_eigenfunc_file, finish_eigenfunc_file, add_eigenpair_to_file
  public :: gs2_save_slice

  interface gs2_restore
     module procedure gs2_restore_many
  end interface

  !> A custom type to look after the netcdf ids for the eigenvalue file
  type EigNetcdfID
     integer :: ncid            !< File handle
     integer :: nconv_count     !< Current size of `conv` dimension
  end type EigNetcdfID

  !> Read and write single or multiple restart files.
  !> Effectively ignored if not built with parallel netcdf
  logical :: read_many = .false., save_many = .false.

  !> Do we want to include the explicit source terms and
  !> related timesteps when saving/restoring from the restart file.
  logical :: include_explicit_source_in_restart = .true.

  !> Which processor should save the potentials. If <0 (default) then
  !> all processors save their potentials. If >= nproc then will be
  !> set to mod(proc_to_save_fields, nproc)
  integer :: proc_to_save_fields

  !> Used to hold the base of the restart file name.
  character (run_name_size) :: restart_file

#ifdef NETCDF
  !> The NETCDF_PARALLEL directives include code for parallel
  !> netcdf using HDF5 to write the output to a single restart file
  !> The save_many/read_many flag allows the old style multiple file output
  !> Here we set a run time flag to enable us to handle some of the different
  !> paths through the code at run time rather than with ifdef mazes
#ifdef NETCDF_PARALLEL
  logical, parameter :: has_netcdf_parallel = .true.
#else
  logical, parameter :: has_netcdf_parallel = .false.
#endif
#endif

contains
#ifdef NETCDF
  integer function get_dim_size(ncid, dim_name) result(dim_size)
    use neasyf, only: neasyf_dim
    use netcdf, only: nf90_inquire_dimension
    implicit none
    integer, intent(in) :: ncid
    character(len = *), intent(in) :: dim_name
    integer :: id, istatus
    call neasyf_dim(ncid, dim_name, dimid = id)
    istatus = nf90_inquire_dimension(ncid, id, len=dim_size)
  end function get_dim_size

  subroutine verify_compatible_dim_sizes(ncid, dim_name, expected_size)
    use mp, only: mp_abort, iproc
    integer, intent(in) :: ncid, expected_size
    character(len = *), intent(in) :: dim_name
    integer :: actual_size
    actual_size = get_dim_size(ncid, dim_name)
    if (actual_size /= expected_size) then
       write(*,'("Restart error (",A,"): got/expected = ",I0,"/",I0," : ",I0)') &
            dim_name, actual_size, expected_size, iproc
       call mp_abort('Incompatible '//dim_name//' in restart', .true.)
    end if
  end subroutine verify_compatible_dim_sizes

  !> Returns the file corresponding to restart file in current setup
  function get_file_proc(is_one_file_per_processor, fileopt) result(file_proc)
    use mp, only: iproc
    use optionals, only: get_option_with_default
    implicit none
    logical, intent(in) :: is_one_file_per_processor
    character(len=*), intent(in), optional :: fileopt
    character(run_name_size) :: file_proc
    file_proc = trim(restart_file) // trim(adjustl(get_option_with_default(fileopt, '')))
    if (is_one_file_per_processor) write(file_proc, '(A,".",I0)') trim(file_proc), iproc
  end function get_file_proc

  logical function variable_exists(ncid, var_name)
    use netcdf, only: nf90_inq_varid, NF90_NOERR
    implicit none
    integer, intent(in) :: ncid
    character(len = *), intent(in) :: var_name
    integer :: id
    variable_exists = nf90_inq_varid(ncid, var_name, id) == NF90_NOERR
  end function variable_exists
#endif

  !> Create and fill the restart files
  subroutine gs2_save_for_restart &
       (g, t0, vnm, has_phi, has_apar, has_bpar, &
       code_dt, code_dt_prev1, code_dt_prev2, code_dt_max, &
       fileopt, &
       save_glo_info_and_grids, save_velocities, header)
#ifdef NETCDF
    use fields_arrays, only: phinew, aparnew, bparnew
    use dist_fn_arrays, only: gexp_1, gexp_2, gexp_3, kx_shift, theta0_shift
    use kt_grids, only: naky, ntheta0, jtwist_out, akx, aky
    use antenna_data, only: a_ant, b_ant, ant_on
    use gs2_metadata, only: create_metadata
    use netcdf, only: nf90_collective, nf90_independent
    use neasyf, only: neasyf_open, neasyf_dim, neasyf_write, neasyf_close
    use mp, only: iproc, barrier, proc0, mp_comm, mp_info
    use theta_grid, only: theta
    use gs2_layouts, only: proc_id,ik_idx,it_idx,is_idx,ie_idx,il_idx, layout
    use layouts_type, only: g_layout_type
    use le_grids, only: energy, al, negrid, nlambda
    use species, only: nspec
    use dist_fn_arrays, only: vpa, vperp2
    use unit_tests, only: debug_message
    use optionals, only: get_option_with_default
#else
    use file_utils, only: error_unit
    use mp, only: proc0
#endif
    use theta_grid, only: ntgrid
    use gs2_layouts, only: g_lo
    use standard_header, only : standard_header_type
    implicit none
    !> The distribution function \(g\)
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in) :: g
    !> Current simulation time
    real, intent (in) :: t0
    !> FIXME: Collisionality multiplier?
    real, dimension (2), intent (in) :: vnm
    !> The current time step, the two previous time steps and the
    !> maximum allowed time step.
    real, intent(in) :: code_dt, code_dt_prev1, code_dt_prev2, code_dt_max
    !> Do we include electrostatic potential, \(\phi\)
    logical, intent (in) :: has_phi
    !> Do we include vector potential, \(A_\parallel\)
    logical, intent (in) :: has_apar
    !> Do we include parallel magnetic field, \(B_\parallel\)
    logical, intent (in) :: has_bpar
    !> Optional extra filename infix
    character (len=*), intent (in), optional :: fileopt
    !> Control if layout and diagnostic information is saved
    logical, intent (in), optional :: save_glo_info_and_grids
    !> Control if \(v_\parallel, v_\perp\) information is saved
    logical, intent (in), optional :: save_velocities
    !> Header for files with build and run information
    type(standard_header_type), optional, intent(in) :: header
#ifdef NETCDF
    !> Names of dimensions for distribution-sized variables in netCDF files
    character(len=5), dimension(*), parameter :: g_dim_names = ["theta", "sign ", "glo  "]
    character(len=5), dimension(*), parameter :: field_dim_names = ["theta", "akx  ", "aky  "]
    !> Arrays holding information about the layout for plotting
    integer, dimension (:,:,:,:,:), allocatable :: procid_kykxsel, iglo_kykxsel
    character (run_name_size) :: file_proc
    real, dimension(:), allocatable :: atmp
    real, dimension(naky) :: stmp
    integer, dimension(5) :: glo_counts
    integer, dimension(3) :: start_pos, counts, field_counts
    integer, dimension(1) :: scalar_count, exb_count
    logical :: is_one_file_per_processor, this_proc_saves_fields, has_gexp_1
    logical :: save_glo_info_and_grids_local, save_velocities_info
    logical :: should_define_single_entry_variables
    integer, parameter :: verb = 3
    integer :: ncid, iglo, n_elements
    ! Whether to use collective (parallel) or independent access to netCDF variables
    integer :: par_access
    type(standard_header_type) :: local_header

    if (present(header)) then
      local_header = header
    else
      local_header = standard_header_type()
    end if

    call debug_message(verb, 'gs2_save::gs2_save_for_restart initialising')

    ! Handle some optional inputs
    save_glo_info_and_grids_local = get_option_with_default(save_glo_info_and_grids, .false.)
    save_velocities_info = get_option_with_default(save_velocities, .false.)

    ! Decide if we want one file per processor or not
    is_one_file_per_processor = save_many .or. (.not. has_netcdf_parallel)

    ! Open the file and define dimensions and variables
    file_proc = get_file_proc(is_one_file_per_processor, fileopt)

    call debug_message(verb, 'gs2_save::gs2_save_for_restart opening file')
    call debug_message(verb, 'file proc is')
    call debug_message(verb, trim(adjustl(file_proc)))

    if (is_one_file_per_processor) then
       ncid = neasyf_open(file_proc, "w")
    else
       call barrier
       call debug_message(verb, &
            'gs2_save::gs2_save_for_restart called barrier before opening')

       ncid = neasyf_open(file_proc, "w", comm=mp_comm, info=mp_info)
    end if

    call create_metadata(ncid, "GS2 restart file", local_header)

    call debug_message(verb, 'gs2_save::gs2_save_for_restart defining dimensions')

    ! Work out the problem size
    if (is_one_file_per_processor) then
       n_elements = max(g_lo%ulim_proc-g_lo%llim_proc+1, 0)
    else
       n_elements = g_lo%ulim_world+1
    end if

    ! Note when processors have no local data we find n_elements = 0 so here we
    ! would be trying to create a zero size dimension. Netcdf interprets this as
    ! an unlimited dimension so neasyf requires unlimited = .true. in this case,
    ! hence the unlimited argument below. This will lead to empty gr/gi data being
    ! written on such processors.
    call neasyf_dim(ncid, "glo", dim_size=n_elements, unlimited = n_elements == 0)

    call neasyf_dim(ncid, "theta", theta)
    call neasyf_dim(ncid, "sign", dim_size=2)
    call neasyf_dim(ncid, "aky", aky)
    call neasyf_dim(ncid, "akx", akx)

    ! Decide if this processor should take part in defining items
    ! which should only appear once in outputs.
    should_define_single_entry_variables = proc0 .or. (.not. is_one_file_per_processor)

    !For saving distribution function
    if (save_glo_info_and_grids_local .and. should_define_single_entry_variables) then
       call neasyf_dim(ncid, "lambda", values=al, long_name="Energy/magnetic moment", units="1/(2 B_a)")
       call neasyf_dim(ncid, "species", dim_size=nspec)

       ! Energy coordinate is a function of species, but the dimension must be 1D
       call neasyf_dim(ncid, "egrid", dim_size=negrid, long_name="See 'energy' for energy grid values")
       call neasyf_write(ncid, "energy", energy, dim_names=["egrid  ", "species"], &
            long_name="Energy grid values")
    end if

    call debug_message(verb, 'gs2_save::gs2_save_for_restart defining variables')

    if (ant_on) call neasyf_dim(ncid, "nk_stir", dim_size=size(a_ant))

    call debug_message(verb, 'gs2_save::gs2_save_for_restart writing scalars')

    if (is_one_file_per_processor .or. proc0) then
      scalar_count = [1]
    else
      scalar_count = [0]
    end if

    call neasyf_write(ncid, "delt0", code_dt, count=scalar_count)
    call neasyf_write(ncid, "delt1", code_dt_prev1, count=scalar_count)
    call neasyf_write(ncid, "delt2", code_dt_prev2, count=scalar_count)
    call neasyf_write(ncid, "delt_max", code_dt_max, count=scalar_count)
    call neasyf_write(ncid, "t0", t0, count=scalar_count)
    call neasyf_write(ncid, "layout", layout, count=scalar_count)
    call neasyf_write(ncid, "ntheta0", ntheta0, count=scalar_count)
    call neasyf_write(ncid, "naky", naky, count=scalar_count)
    call neasyf_write(ncid, "nlambda", nlambda, count=scalar_count)
    call neasyf_write(ncid, "negrid", negrid, count=scalar_count)
    call neasyf_write(ncid, "nspec", nspec, count=scalar_count)
    call neasyf_write(ncid, "vnm1", vnm(1), count=scalar_count)
    call neasyf_write(ncid, "vnm2", vnm(2), count=scalar_count)

    if (ant_on) then
       atmp = real(a_ant)
       call neasyf_write(ncid, "a_ant_r", atmp, dim_names=["nk_stir"])
       atmp = aimag(a_ant)
       call neasyf_write(ncid, "a_ant_i", atmp, dim_names=["nk_stir"])

       atmp = real(b_ant)
       call neasyf_write(ncid, "b_ant_r", atmp, dim_names=["nk_stir"])
       atmp = aimag(b_ant)
       call neasyf_write(ncid, "b_ant_i", atmp, dim_names=["nk_stir"])
    end if

    call debug_message(verb, 'gs2_save::gs2_save_for_restart writing dist fn')

    if (is_one_file_per_processor) then
       start_pos = [1, 1, 1]
       par_access = nf90_independent
    else
       start_pos = [1, 1, g_lo%llim_proc+1]
       par_access = nf90_collective
    end if

    counts = shape(g)

    call neasyf_write(ncid, "gr", real(g), dim_names=g_dim_names, &
         start=start_pos, count=counts, par_access=par_access)
    call neasyf_write(ncid, "gi", aimag(g), dim_names=g_dim_names, &
         start=start_pos, count=counts, par_access=par_access)

    if (include_explicit_source_in_restart) then
       !Explicit source term gexp_1
#ifndef SHMEM
       has_gexp_1 = allocated(gexp_1)
#else
       has_gexp_1 = associated(gexp_1)
#endif

       if (has_gexp_1) then
          call neasyf_write(ncid, "gexp_1_r", real(gexp_1), dim_names=g_dim_names, &
                  start=start_pos, count=counts, par_access=par_access)
          call neasyf_write(ncid, "gexp_1_i", aimag(gexp_1), dim_names=g_dim_names, &
                  start=start_pos, count=counts, par_access=par_access)
       end if

       !Explicit source term gexp_2
       if(allocated(gexp_2)) then
          call neasyf_write(ncid, "gexp_2_r", real(gexp_2), dim_names=g_dim_names, &
                  start=start_pos, count=counts, par_access=par_access)
          call neasyf_write(ncid, "gexp_2_i", aimag(gexp_2), dim_names=g_dim_names, &
                  start=start_pos, count=counts, par_access=par_access)
       end if

       !Explicit source term gexp_3
       if(allocated(gexp_3)) then
          call neasyf_write(ncid, "gexp_3_r", real(gexp_3), dim_names=g_dim_names, &
                  start=start_pos, count=counts, par_access=par_access)
          call neasyf_write(ncid, "gexp_3_i", aimag(gexp_3), dim_names=g_dim_names, &
                  start=start_pos, count=counts, par_access=par_access)
       end if
    end if

    ! For saving distribution function
    if (save_glo_info_and_grids_local) then
       if (proc0 .or. .not. is_one_file_per_processor) then
          if(.not. allocated(iglo_kykxsel)) allocate(iglo_kykxsel(naky,ntheta0,nspec,negrid,nlambda))
          if(.not. allocated(procid_kykxsel)) allocate(procid_kykxsel(naky,ntheta0,nspec,negrid,nlambda))
       end if

       if (proc0) then
          do iglo=g_lo%llim_world,g_lo%ulim_world
             iglo_kykxsel(ik_idx(g_lo,iglo),it_idx(g_lo,iglo),is_idx(g_lo,iglo), &
                  ie_idx(g_lo,iglo),il_idx(g_lo,iglo)) = iglo
          end do

          do iglo=g_lo%llim_world,g_lo%ulim_world
             procid_kykxsel(ik_idx(g_lo,iglo),it_idx(g_lo,iglo),is_idx(g_lo,iglo), &
                  ie_idx(g_lo,iglo),il_idx(g_lo,iglo)) = proc_id(g_lo,iglo)
          end do

          glo_counts = shape(iglo_kykxsel)
       else
          glo_counts = [0, 0, 0, 0, 0]
       end if

       ! Either only write to proc0's serial file, or all processes need to call
       ! these functions in parallel
       if (proc0 .or. .not. is_one_file_per_processor) then
          call neasyf_write(ncid, "iglo_llim_proc", g_lo%llim_proc, count=scalar_count)

          ! Fill processor and layout information for plotting purposes
          call neasyf_write(ncid, "iglo_kykxsel", iglo_kykxsel, &
               dim_names=[character(len=7)::"aky", "akx", "species", "egrid", "lambda"], &
               count=glo_counts)
          call neasyf_write(ncid, "procid_kykxsel", procid_kykxsel, &
               dim_names=[character(len=7)::"aky", "akx", "species", "egrid", "lambda"], &
               count=glo_counts)

          call neasyf_write(ncid, "jtwist", jtwist_out, count=scalar_count)
       end if
    end if

    if (save_velocities_info) then
       call neasyf_write(ncid, "vpa", vpa, dim_names=["theta", "sign ", "glo  "], &
            start=start_pos, count=counts)
       call neasyf_write(ncid, "vperp2", vperp2, dim_names=["theta", "glo  "], &
            start=start_pos(1:3:2), count=counts(1:3:2))
    end if

    ! Decide if we are saving the fields
    this_proc_saves_fields = (iproc == proc_to_save_fields) .or. (proc_to_save_fields < 0)

    if (this_proc_saves_fields .or. .not. is_one_file_per_processor) then
       ! In parallel, all processes need to call `neasyf_write`, but
       ! only one should actually write any data
       if (this_proc_saves_fields) then
          field_counts = shape(phinew)
       else
          field_counts = [0, 0, 0]
       end if

       if (has_phi) then
          call neasyf_write(ncid, "phi_r", real(phinew), dim_names=field_dim_names, count=field_counts)
          call neasyf_write(ncid, "phi_i", aimag(phinew), dim_names=field_dim_names, count=field_counts)
       end if

       if (has_apar) then
          call neasyf_write(ncid, "apar_r", real(aparnew), dim_names=field_dim_names, count=field_counts)
          call neasyf_write(ncid, "apar_i", aimag(aparnew), dim_names=field_dim_names, count=field_counts)
       end if

       if (has_bpar) then
          call neasyf_write(ncid, "bpar_r", real(bparnew), dim_names=field_dim_names, count=field_counts)
          call neasyf_write(ncid, "bpar_i", aimag(bparnew), dim_names=field_dim_names, count=field_counts)
       end if
    end if

    ! Because we want to be able to restart using exb shear from a
    ! case which does not have exb shear we always add kx_shift and
    ! theta0_shift even if no exb shear present in simulation)
    if (is_one_file_per_processor .or. proc0) then
       exb_count = [naky]
    else
       exb_count = [0]
    end if

    if (allocated(kx_shift)) then
       stmp = kx_shift
    else
       stmp = 0.
    end if
    call neasyf_write(ncid, "kx_shift", stmp, dim_names=["aky"], count=exb_count)

    if (allocated(theta0_shift)) then
       stmp = theta0_shift
    else
       stmp = 0.
    end if
    call neasyf_write(ncid, "theta0_shift", stmp, dim_names=["aky"], count=exb_count)

    ! Close the file
    call neasyf_close(ncid)

    ! Deallocate local arrays
    if (allocated(iglo_kykxsel)) deallocate(iglo_kykxsel)
    if (allocated(procid_kykxsel)) deallocate(procid_kykxsel)

#else
    if (proc0) write (error_unit(),*) 'WARNING: gs2_save_for_restart is called without netcdf library'
    UNUSED_DUMMY(g); UNUSED_DUMMY(t0); UNUSED_DUMMY(vnm); UNUSED_DUMMY(has_phi)
    UNUSED_DUMMY(has_apar); UNUSED_DUMMY(has_bpar); UNUSED_DUMMY(code_dt)
    UNUSED_DUMMY(code_dt_prev1); UNUSED_DUMMY(code_dt_prev2); UNUSED_DUMMY(code_dt_max)
    UNUSED_DUMMY(fileopt); UNUSED_DUMMY(save_glo_info_and_grids)
    UNUSED_DUMMY(save_velocities); UNUSED_DUMMY(header)
#endif
  end subroutine gs2_save_for_restart

  !> Save g slices
  subroutine gs2_save_slice(g_in, it, ik, il, ie, is, time, nout, fileopt, header)
# ifdef NETCDF
    use neasyf, only: neasyf_open, neasyf_dim, neasyf_write, neasyf_close
    use gs2_metadata, only: create_metadata
    use unit_tests, only: debug_message
    use mp, only: iproc, proc0
    use optionals, only: get_option_with_default
    use convert, only: c2r
    use kt_grids, only: aky, akx, naky, ntheta0
    use le_grids, only: energy, al, negrid, nlambda
    use species, only: nspec
    use theta_grid, only: theta
    use gs2_layouts, only: proc_id, idx
# endif
    use theta_grid, only: ntgrid
    use gs2_layouts, only: g_lo
    use standard_header, only : standard_header_type
    implicit none
    !> The distribution function \(g\)
    complex, dimension (-ntgrid:, :, g_lo%llim_proc:), intent (in) :: g_in
    !> The indices of the slice to save
    integer, intent(in) :: it, ik, il, ie, is
    !> The current time
    real, intent(in) :: time
    !> The current time output index
    integer, intent(in), optional :: nout
    !> Optional extra filename infix
    character (len=*), intent (in), optional :: fileopt
    !> Header for files with build and run information
    type(standard_header_type), optional, intent(in) :: header
# ifdef NETCDF
    character (run_name_size) :: file_proc
    integer :: the_iglo, the_proc
    integer, parameter :: verb = 3
    type(standard_header_type) :: local_header
    character (len=*), parameter :: file_opt_base = '.slice'
    logical :: appending, previous_file_exists
    real, dimension(:, :, :), allocatable :: ri_gslice
    integer, save :: nout_local = 1
    logical :: error_exit
    ! File handle
    integer :: ncid

    call debug_message(verb, 'gs2_save::gs2_save_slice initialising')

    ! Sanity check the inputs
    error_exit = .false.
    if (it < 1 .or. it > ntheta0) error_exit = .true.
    if (ik < 1 .or. ik > naky) error_exit = .true.
    if (il < 1 .or. il > nlambda) error_exit = .true.
    if (ie < 1 .or. ie > negrid) error_exit = .true.
    if (is < 1 .or. is > nspec) error_exit = .true.
    if (proc0 .and. error_exit) write(6,'("Error: Invalid index based to gs2_save_slice.")')
    if (error_exit) return

    the_iglo = idx(g_lo, ik, it, il, ie, is)
    the_proc = proc_id(g_lo, the_iglo)

    ! The proc which owns this slice is responsible for writing and
    ! does the rest of the work so let us return if that's not this
    ! proc
    if (iproc /= the_proc) return

    if (present(header)) then
      local_header = header
    else
      local_header = standard_header_type()
    end if

    ! Open the file and define dimensions and variables
    file_proc = get_file_proc(.false., file_opt_base // fileopt)

    call debug_message(verb, 'gs2_save::gs2_save_slice opening file')
    call debug_message(verb, 'file proc is')
    call debug_message(verb, trim(adjustl(file_proc)))

    nout_local = get_option_with_default(nout, nout_local)

    ! See if the output file already exists, so we know how to open
    ! the file.  This would be simpler if netCDF had a "open if it
    ! exists, create if it doesn't" mode, but it doesn't so we have
    ! to do this instead
    inquire(file=trim(file_proc), exist=previous_file_exists)
    appending = (previous_file_exists .and. nout_local > 1)

    if (appending) then
       ncid = neasyf_open(trim(file_proc), 'rw')
    else
       ncid = neasyf_open(trim(file_proc), 'w')
       call create_metadata(ncid, "GS2 slice file", local_header)

       call debug_message(verb, 'gs2_save::gs2_save_slice defining dimensions')
       call neasyf_dim(ncid, "ri", dim_size = 2, long_name = "real/imaginary components")
       call neasyf_dim(ncid, "theta", values = theta, long_name="Poloidal angle", units="rad")
       call neasyf_dim(ncid, "sign",  dim_size = 2, long_name = "sign of v||")
       call neasyf_dim(ncid, "t", unlimited = .true., long_name = "Time", units = "L/vt")

       call debug_message(verb, 'gs2_save::gs2_save_slice defining constant variables')
       call neasyf_write(ncid, "it", it)
       call neasyf_write(ncid, "ik", ik)
       call neasyf_write(ncid, "il", il)
       call neasyf_write(ncid, "ie", ie)
       call neasyf_write(ncid, "is", is)

       call neasyf_write(ncid, "iglo", the_iglo)
       call neasyf_write(ncid, "iproc", the_proc)

       call neasyf_write(ncid, "kx", akx(it))
       call neasyf_write(ncid, "ky", aky(ik))
       call neasyf_write(ncid, "lambda", al(il))
       call neasyf_write(ncid, "energy", energy(ie, is))
    end if

    call neasyf_write(ncid, "t", time, start = [nout_local])

    ! Convert from complex to real
    allocate(ri_gslice(2, size(theta), 2))
    call c2r(g_in(:, :, the_iglo), ri_gslice)
    call neasyf_write(ncid, "g_slice", ri_gslice, start = [1, 1, 1, nout_local])

    ! Close the file
    call neasyf_close(ncid)

    nout_local = nout_local + 1
# else
    UNUSED_DUMMY(g_in); UNUSED_DUMMY(it); UNUSED_DUMMY(ik) ; UNUSED_DUMMY(il)
    UNUSED_DUMMY(ie); UNUSED_DUMMY(is); UNUSED_DUMMY(time); UNUSED_DUMMY(nout)
    UNUSED_DUMMY(fileopt); UNUSED_DUMMY(header)
# endif
  end subroutine gs2_save_slice

  !> FIXME : Add documentation
  subroutine gs2_restore_many (g, scale, has_phi, has_apar, has_bpar, fileopt)
!MR, 2007: restore kx_shift array if already allocated
# ifdef NETCDF
    use array_utils, only: zero_array
    use mp, only: iproc, proc0, broadcast, mp_comm, mp_info
    use fields_arrays, only: phinew, aparnew, bparnew
    use fields_arrays, only: phi, apar, bpar
    use dist_fn_arrays, only: gexp_1, gexp_2, gexp_3, kx_shift, theta0_shift
    use kt_grids, only: naky, ntheta0
    use gs2_layouts, only: layout
    use neasyf, only: neasyf_open, neasyf_read, neasyf_close
    use netcdf_utils, only: check_netcdf_file_precision, kind_nf
# endif
    use mp, only: mp_abort
    use theta_grid, only: ntgrid
    use gs2_layouts, only: g_lo
    implicit none
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (out) :: g
    real, intent (in) :: scale
    logical, intent (in) :: has_phi, has_apar, has_bpar
    character (len=*), intent(in), optional :: fileopt
# ifdef NETCDF
    logical :: this_proc_reads_fields, need_to_broadcast_fields
    integer, dimension(3) :: counts, start_pos
    real, dimension(naky) :: stmp
    character (run_name_size) :: file_proc
    character (len = 5) :: restart_layout
    real :: fac
    logical :: has_explicit_source_terms, is_one_file_per_processor, has_gexp_1
    logical, save :: explicit_warning_given = .false.
    integer(kind_nf) :: ncid
    real, allocatable, dimension(:,:,:) :: tmpr, tmpi, ftmpr, ftmpi

    is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel)
    file_proc = get_file_proc(is_one_file_per_processor, fileopt)

    if (is_one_file_per_processor) then
       ncid = neasyf_open(file_proc, "r")
    else
       ncid = neasyf_open(file_proc, "r", comm=mp_comm, info=mp_info)
    end if

    call check_netcdf_file_precision (ncid)

    ! Check dimension sizes are compatible with the current simulation inputs
    call verify_compatible_dim_sizes(ncid, "theta", 2 * ntgrid + 1)
    call verify_compatible_dim_sizes(ncid, "sign", 2)
    call verify_compatible_dim_sizes(ncid, "aky", naky)
    call verify_compatible_dim_sizes(ncid, "akx", ntheta0)
    if (is_one_file_per_processor) then
       call verify_compatible_dim_sizes(ncid, "glo", g_lo%ulim_proc - g_lo%llim_proc + 1)
    else
       call verify_compatible_dim_sizes(ncid, "glo", g_lo%ulim_world + 1)
    end if

    if (include_explicit_source_in_restart) then
       has_explicit_source_terms = .true.

#ifndef SHMEM
       has_gexp_1 = allocated(gexp_1)
#else
       has_gexp_1 = associated(gexp_1)
#endif

       if (has_gexp_1) &
          has_explicit_source_terms = has_explicit_source_terms .and. &
               (variable_exists(ncid, "gexp_1_r") .and. variable_exists(ncid, "gexp_1_i"))

       if (allocated(gexp_2)) &
          has_explicit_source_terms = has_explicit_source_terms .and. &
               (variable_exists(ncid, "gexp_2_r") .and. variable_exists(ncid, "gexp_2_i"))

       if (allocated(gexp_3))  &
          has_explicit_source_terms = has_explicit_source_terms .and. &
               (variable_exists(ncid, "gexp_3_r") .and. variable_exists(ncid, "gexp_3_i"))

       if ((.not. explicit_warning_given).and. proc0 .and. .not. has_explicit_source_terms) then
          write(*, '("Warning: At least one explicit source term absent in restart file.")')
          write(*, '("  Will not load these terms. This is probably OK unless you expect the restart file to contain these terms.")')
          write(*, '("  This warning will not be repeated.")')

          ! Update this flag to ensure we don't give this warning again
          explicit_warning_given = .true.
       end if

    else
       ! Always skip trying to get the explicit source terms if we're not
       ! including them
       has_explicit_source_terms = .false.
    end if

    call neasyf_read(ncid, "layout", restart_layout)

    ! Abort if the layouts don't match
    if(restart_layout /= layout) then
       call mp_abort("Incompatible layouts in restart file ("//restart_layout//") and simulation ("//layout//")",.true.)
    end if

    allocate (tmpr(2*ntgrid+1,2,g_lo%llim_proc:g_lo%ulim_alloc))
    allocate (tmpi(2*ntgrid+1,2,g_lo%llim_proc:g_lo%ulim_alloc))

    call zero_array(tmpr); call zero_array(tmpi)

    if (is_one_file_per_processor) then
       start_pos = [1, 1, 1]
    else
       start_pos = [1, 1, g_lo%llim_proc+1]
    end if
    counts = shape(tmpr)

    call neasyf_read(ncid, "gr", tmpr, start=start_pos, count=counts)
    call neasyf_read(ncid, "gi", tmpi, start=start_pos, count=counts)
    g = cmplx(tmpr, tmpi)

    if (include_explicit_source_in_restart .and. has_explicit_source_terms) then
       ! Explicit source term
#ifndef SHMEM
       has_gexp_1 = allocated(gexp_1)
#else
       has_gexp_1 = associated(gexp_1)
#endif

       if (has_gexp_1) then
          call neasyf_read(ncid, "gexp_1_r", tmpr, start=start_pos, count=counts)
          call neasyf_read(ncid, "gexp_1_i", tmpi, start=start_pos, count=counts)
          gexp_1 = cmplx(tmpr, tmpi)
       endif

       ! Explicit source term
       if(allocated(gexp_2)) then
          call neasyf_read(ncid, "gexp_2_r", tmpr, start=start_pos, count=counts)
          call neasyf_read(ncid, "gexp_2_i", tmpi, start=start_pos, count=counts)
          gexp_2 = cmplx(tmpr, tmpi)
       endif

       ! Explicit source term
       if(allocated(gexp_3)) then
          call neasyf_read(ncid, "gexp_3_r", tmpr, start=start_pos, count=counts)
          call neasyf_read(ncid, "gexp_3_i", tmpi, start=start_pos, count=counts)
          gexp_3 = cmplx(tmpr, tmpi)
       endif
    end if

    deallocate(tmpr, tmpi)

    if (allocated(kx_shift)) then   ! MR begin
       call neasyf_read(ncid, "kx_shift", stmp)
       kx_shift = stmp
    endif   ! MR end

    if (allocated(theta0_shift)) then
       call neasyf_read(ncid, "theta0_shift", stmp)
       theta0_shift = stmp
    endif

    ! Decide if we are reading the fields
    need_to_broadcast_fields = proc_to_save_fields >= 0
    this_proc_reads_fields = (iproc == proc_to_save_fields) .or. (.not. need_to_broadcast_fields)

    if (this_proc_reads_fields) then
       allocate (ftmpr(2*ntgrid+1,ntheta0,naky))
       allocate (ftmpi(2*ntgrid+1,ntheta0,naky))

       if (has_phi) then
          call neasyf_read(ncid, "phi_r", ftmpr)
          call neasyf_read(ncid, "phi_i", ftmpi)
          phinew = cmplx(ftmpr, ftmpi)
       end if

       if (has_apar) then
          call neasyf_read(ncid, "apar_r", ftmpr)
          call neasyf_read(ncid, "apar_i", ftmpi)
          aparnew = cmplx(ftmpr, ftmpi)
       end if

       if (has_bpar) then
          call neasyf_read(ncid, "bpar_r", ftmpr)
          call neasyf_read(ncid, "bpar_i", ftmpi)
          bparnew = cmplx(ftmpr, ftmpi)
       end if

       deallocate(ftmpr, ftmpi)
    end if

    call neasyf_close(ncid)

    if (need_to_broadcast_fields) then
       if (has_phi) call broadcast(phinew, proc_to_save_fields)
       if (has_apar) call broadcast(aparnew, proc_to_save_fields)
       if (has_bpar) call broadcast(bparnew, proc_to_save_fields)
    end if

    if (has_phi) call zero_array(phi)
    if (has_apar) call zero_array(apar)
    if (has_bpar) call zero_array(bpar)

    if (scale > 0.) then
       g = g*scale
       phinew = phinew*scale
       aparnew = aparnew*scale
       bparnew = bparnew*scale
    else
       fac = - scale/(maxval(abs(phinew)))
       g = g*fac
       phinew = phinew*fac
       aparnew = aparnew*fac
       bparnew = bparnew*fac
    end if
    
# else
    call mp_abort("Cannot restore from restart file without netcdf.", .true.)
    UNUSED_DUMMY(scale); UNUSED_DUMMY(has_phi); UNUSED_DUMMY(has_apar)
    UNUSED_DUMMY(has_bpar); UNUSED_DUMMY(fileopt) ; g = 0.0
# endif
  end subroutine gs2_restore_many

  !> This routine writes a passed square complex array to a file
  !> with passed name
  subroutine gs2_save_response(resp, fname, code_dt, condition, header)
    use file_utils, only: error_unit
    use standard_header, only: standard_header_type
#ifdef NETCDF
    use optionals, only: get_option_with_default
    use convert, only: c2r
    use gs2_metadata, only: create_metadata
    use neasyf, only: neasyf_open, neasyf_dim, neasyf_write, neasyf_close
#else
    use file_utils, only: get_unused_unit
#endif
    implicit none
    complex,dimension(:,:), intent(in) :: resp
    character(len=*), intent(in) :: fname
    real, intent(in) :: code_dt
    real, intent(in), optional :: condition
    !> Header for files with build and run information
    type(standard_header_type), optional, intent(in) :: header
    integer :: sz

#ifdef NETCDF
    integer :: ncid
    real, dimension(:,:,:), allocatable :: ri_resp
    ! Actual value for optional `header` input
    type(standard_header_type) :: local_header
#else
    integer :: unit
#endif
    !Currently only support serial writing, but could be by any proc
    !so we have to make sure only one proc calls this routine

    !Verify we have a square array
    sz=size(resp(:,1))
    if(sz/=size(resp(1,:))) then
       write(error_unit(),'("Error: gs2_save_response expects a square array input.")')
       return
    endif

#ifdef NETCDF
    if (present(header)) then
      local_header = header
    else
      local_header = standard_header_type()
    end if

    ncid = neasyf_open(fname, "w")
    call create_metadata(ncid, "GS2 response matrix", local_header)

    call neasyf_dim(ncid, "ri", dim_size=2)
    call neasyf_dim(ncid, "ax1", dim_size=sz)
    call neasyf_dim(ncid, "ax2", dim_size=sz)

    call neasyf_write(ncid, "dt", code_dt)
    call neasyf_write(ncid, "condition", get_option_with_default(condition, -1.0))

    !/Convert complex to ri and write
    allocate(ri_resp(2,sz,sz))
    call c2r(resp,ri_resp)
    call neasyf_write(ncid, "response", ri_resp, dim_names=["ri ", "ax1", "ax2"])
    deallocate(ri_resp)

    !/Now close the file
    call neasyf_close(ncid)

#else
!Fall back on binary output if no NETCDF
    !Get a free unit
    call get_unused_unit(unit)

    !Open file and write
    open(unit=unit,file=fname,form="unformatted")
    write(unit) resp
    close(unit)

    UNUSED_DUMMY(code_dt); UNUSED_DUMMY(condition); UNUSED_DUMMY(header)
#endif
  end subroutine gs2_save_response

  !> This routine reads a square complex array from a file
  !! with passed name
  subroutine gs2_restore_response(resp, fname, code_dt, condition)
    use file_utils, only: error_unit
#ifdef NETCDF
    use neasyf, only: neasyf_open, neasyf_close, neasyf_read
    use convert, only: r2c
#else
    use file_utils, only: get_unused_unit
#endif
    implicit none
    complex, dimension(:,:), intent(out) :: resp
    character(len=*), intent(in) :: fname
    real, intent(out) :: code_dt
    real, intent(out), optional :: condition
    integer :: sz
#ifdef NETCDF
    integer :: ncid
    real, dimension(:,:,:), allocatable :: ri_resp
#else
    integer :: unit
#endif
    !Currently only support serial reading, but could be by any proc
    !so we have to make sure only one proc calls this routine

    !Verify we have a square array
    sz=size(resp(:,1))
    if(sz/=size(resp(1,:))) then
       write(error_unit(),'("Error: gs2_restore_response expects a square array output.")')
       return
    endif

#ifdef NETCDF
    ncid = neasyf_open(fname, "r")

    !/Read and convert ri to complex
    allocate(ri_resp(2,sz,sz))
    call neasyf_read(ncid, "response", ri_resp)
    call r2c(resp,ri_resp)
    deallocate(ri_resp)

    !/Get the time step
    call neasyf_read(ncid, "dt", code_dt)

    !/Get the condition number
    if (present(condition)) call neasyf_read(ncid, "condition", condition)

    !/Now close the file
    call neasyf_close(ncid)
#else
!Fall back on binary output if no NETCDF
    !Get a free unit
    call get_unused_unit(unit)

    !Open file and write
    open(unit=unit,file=fname,form="unformatted")
    read(unit) resp
    close(unit)

    !Set outputs to something invalid
    code_dt = -1 ; condition = 0.0
#endif
  end subroutine gs2_restore_response

  !> Initialises a file for saving output of eigensolver to netcdf
  subroutine init_eigenfunc_file(fname, IDs)
#ifdef NETCDF
    use theta_grid, only: theta
    use neasyf, only: neasyf_open, neasyf_dim
#endif
    implicit none
    character(len=*), intent(in) :: fname
    type(EigNetcdfID), intent(inout) :: IDs

    !Set nconv counter to 0
    IDs%nconv_count=0

#ifdef NETCDF
    IDs%ncid = neasyf_open(fname, "w")

    call neasyf_dim(IDs%ncid, "ri", dim_size=2)
    call neasyf_dim(IDs%ncid, "theta", theta)
    call neasyf_dim(IDs%ncid, "conv", unlimited=.true.)
#else
    UNUSED_DUMMY(fname)
#endif
  end subroutine init_eigenfunc_file

  !> Add an eigenpairs data to file
  subroutine add_eigenpair_to_file(eval, has_phi, has_apar, has_bpar, IDs, my_conv)
#ifdef NETCDF
    use fields_arrays, only: phinew, aparnew, bparnew
    use convert, only: c2r
    use theta_grid, only: ntgrid
    use neasyf, only: neasyf_write
    use optionals, only: get_option_with_default
#endif
    complex, intent(in) :: eval !Note just use fields to get eigenvectors
    real, intent(in), optional :: my_conv
    type(EigNetcdfID), intent(inout) :: IDs
    logical, intent(in) :: has_phi, has_apar, has_bpar
#ifdef NETCDF
    real, dimension(2) :: ri_omega
    real, dimension(:,:), allocatable :: ri_field
    integer, dimension(3) :: start_field, count_field
    real :: local_conv
    character(*), dimension(*), parameter :: eigen_field_dim_names = [character(len=5)::"ri", "theta", "conv"]

    !First increment counter
    IDs%nconv_count=IDs%nconv_count+1
    start_field = [1, 1, IDs%nconv_count]
    count_field = [2, 2*ntgrid+1, 1]

    local_conv = get_option_with_default(my_conv, real(IDs%nconv_count))

    call neasyf_write(IDs%ncid, "conv", local_conv, dim_names=["conv"], &
         start=[IDs%nconv_count])

    !/Omega
    ri_omega(1)=real(eval)
    ri_omega(2)=aimag(eval)
    call neasyf_write(IDs%ncid, "omega", ri_omega, dim_names=["ri  ", "conv"], &
         start=[1, IDs%nconv_count], count=[2,1])

    !/Fields
    allocate(ri_field(2,2*ntgrid+1))
    if(has_phi)then
       call c2r(phinew(:,1,1),ri_field)
       call neasyf_write(IDs%ncid, "phi", ri_field, dim_names=eigen_field_dim_names, &
            start=start_field, count=count_field)
    endif
    if(has_apar)then
       call c2r(aparnew(:,1,1),ri_field)
       call neasyf_write(IDs%ncid, "apar", ri_field, dim_names=eigen_field_dim_names, &
            start=start_field, count=count_field)
    endif
    if(has_bpar)then
       call c2r(bparnew(:,1,1),ri_field)
       call neasyf_write(IDs%ncid, "bpar", ri_field, dim_names=eigen_field_dim_names, &
            start=start_field, count=count_field)
    endif
    deallocate(ri_field)
#else
    UNUSED_DUMMY(eval); UNUSED_DUMMY(has_phi); UNUSED_DUMMY(has_apar); UNUSED_DUMMY(has_bpar)
    UNUSED_DUMMY(IDs); UNUSED_DUMMY(my_conv)
#endif
  end subroutine add_eigenpair_to_file

  !> Close the eigenfunction file
  subroutine finish_eigenfunc_file(IDs)
#ifdef NETCDF
    use neasyf, only: neasyf_close
#endif
    implicit none
    type(EigNetcdfID), intent(inout) :: IDs
#ifdef NETCDF
    call neasyf_close(IDs%ncid)
#else
    UNUSED_DUMMY(IDs)
#endif
  end subroutine finish_eigenfunc_file

  !> This function checks to see if we can create a file with name
  !! <restart_file>//<SomeSuffix> if not then our restarts are not
  !! going to be possible and we return false. Can also be used to check
  !! that we can read from the restart file (which assumes it exists).
  function restart_writable(read_only, my_file, error_message)
    use mp, only: proc0, broadcast
    use file_utils, only: get_unused_unit
    use optionals, only: get_option_with_default
    implicit none
    !> If present and true, only check that files can be read
    logical, intent(in), optional :: read_only
    !> An optional specific filename to check
    character(len=*),intent(in),optional::my_file
    !> Error message returned from `open` if there was a problem
    character(len=:), allocatable, optional, intent(out) :: error_message

    character(len=*), parameter :: SuffixTmp = '.ThisIsATestFile'
    character(9) :: open_mode
    character(6) :: close_mode
    character(run_name_size) :: local_file, filename
    logical :: restart_writable
    integer :: unit, ierr
    character(len=200) :: io_err_msg

    ierr=-200
    local_file = trim(get_option_with_default(my_file, restart_file))

    ! On proc0 try to open tmp file for writing
    if (proc0)then
       call get_unused_unit(unit)

       ! Default open and close modes
       open_mode = 'readwrite'
       close_mode = 'delete'

       ! If we want to test write capability then do it with an unusual
       ! file name to prevent clobber
       filename = trim(local_file) // SuffixTmp

       ! Set filemode to READ if read_only=T
       if (get_option_with_default(read_only, .false.)) then
          ! If we're only checking a file can be read, then don't
          ! delete the file after we're done
          open_mode = 'read'
          close_mode = 'keep'
          ! If checking readonly then we need to make sure we try to
          ! read from an existing file
          filename = local_file
       end if

       open(unit=unit, file=trim(filename), &
            iostat=ierr, action=open_mode, &
            iomsg=io_err_msg)

       ! If open was successful then we can close the file and delete it
       if (ierr == 0) close(unit=unit, status=close_mode)
    endif

    ! Now make sure everyone knows the answer
    call broadcast(ierr)
    restart_writable = (ierr == 0)

    if (.not. present(error_message)) return

    if (restart_writable) then
      error_message = ""
    else
      error_message = trim(io_err_msg)
    end if
  end function restart_writable

  !> FIXME : Add documentation
  subroutine init_gs2_save
  end subroutine init_gs2_save

  !> Sets the base of the restart file to be used when
  !> writing/reading files.
  subroutine set_restart_file (file)
    use mp, only: proc0
    use file_utils, only: error_unit
    implicit none
    character(len=*), intent (in) :: file
    if (proc0 .and. len_trim(file) > run_name_size) then
       write( error_unit(), '("Argument to set_restart_file exceeds restart_file size so truncating")')
    end if
    restart_file = file(1:min(len_trim(file), run_name_size))
  end subroutine set_restart_file

  !> FIXME : Add documentation  
  subroutine finish_gs2_save
  end subroutine finish_gs2_save

  !> FIXME : Add documentation
  subroutine init_dt (delt0, delt1, delt2, delt_max, not_set_value)
# ifdef NETCDF
    use mp, only: proc0, broadcast
    use optionals, only: get_option_with_default
    use neasyf, only: neasyf_open, neasyf_read, neasyf_close
    use netcdf_utils, only: kind_nf
# else
    use mp, only: mp_abort
# endif
    implicit none
    real, intent (in out) :: delt0, delt1, delt2, delt_max
    real, intent(in), optional :: not_set_value
# ifdef NETCDF
    character (run_name_size) :: file_proc
    real :: not_set_value_to_use
    logical :: is_one_file_per_processor
    ! NetCDF handles
    integer(kind_nf) :: ncid
    real, dimension(0:3) :: time_values
    if (proc0) then
       is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel)
       file_proc = get_file_proc(is_one_file_per_processor)
       ncid = neasyf_open(file_proc, "r")

       not_set_value_to_use = get_option_with_default(not_set_value, -1.0)
       time_values = not_set_value_to_use
       ! Note unlike the explicit source terms all three time steps should always
       ! be available in the restart file so we don't silence the error messages here.
       ! The only situation we are likely to come across where the delt1 and delt2
       ! values aren't available is where we are trying to read an old restart file.
       ! This will then lead to the error/warning being displayed but the code should
       ! carry on as intended and the missing steps will be set to a special value to
       ! indicate they have not been set yet.
       if (variable_exists(ncid, "delt0")) call neasyf_read(ncid, "delt0", time_values(0))
       if (variable_exists(ncid, "delt1")) call neasyf_read(ncid, "delt1", time_values(1))
       if (variable_exists(ncid, "delt2")) call neasyf_read(ncid, "delt2", time_values(2))
       if (variable_exists(ncid, "delt_max")) call neasyf_read(ncid, "delt_max", time_values(3))              
       call neasyf_close(ncid)
    endif

    call broadcast(time_values)
    delt0 = time_values(0) ; delt1 = time_values(1) ; delt2 = time_values(2)
    delt_max = time_values(3)
# else
    call mp_abort("Cannot load initial dt from restart without netcdf.", .true.)
    UNUSED_DUMMY(delt0); UNUSED_DUMMY(delt1); UNUSED_DUMMY(delt2); UNUSED_DUMMY(delt_max)
    UNUSED_DUMMY(not_set_value)
# endif
  end subroutine init_dt

!> FIXME : Add documentation
  subroutine init_vnm (vnm)
# ifdef NETCDF
    use mp, only: proc0, broadcast
    use neasyf, only: neasyf_open, neasyf_read, neasyf_close
    use netcdf_utils, only: kind_nf
# else
    use mp, only: mp_abort
# endif
    implicit none
    real, dimension(2), intent (in out) :: vnm
# ifdef NETCDF
    character (run_name_size) :: file_proc
    logical :: is_one_file_per_processor
    ! NetCDF handles
    integer(kind_nf) :: ncid

    if (proc0) then
       is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel)
       file_proc = get_file_proc(is_one_file_per_processor)
       ncid = neasyf_open (file_proc, "r")
       call neasyf_read(ncid, "vnm1", vnm(1))
       call neasyf_read(ncid, "vnm2", vnm(2))
       call neasyf_close(ncid)
    endif

    call broadcast (vnm)
# else
    call mp_abort("Cannot load vnm from restart without netcdf.", .true.)
    UNUSED_DUMMY(vnm)
# endif
  end subroutine init_vnm

  !> FIXME : Add documentation
  !!
  !! @note This routine gets a_ant and b_ant for proc 0 only!
  subroutine init_ant_amp (a_ant, b_ant, nk_stir)
# ifdef NETCDF
    use constants, only: zi
    use neasyf, only: neasyf_open, neasyf_read, neasyf_close
    use mp, only: proc0, broadcast, mp_abort
    use netcdf_utils, only: kind_nf
# else
    use mp, only: mp_abort
# endif
    implicit none
    complex, dimension(:), intent (in out) :: a_ant, b_ant
    integer, intent (in) :: nk_stir
# ifdef NETCDF
    character (run_name_size) :: file_proc
    logical :: is_one_file_per_processor
    ! NetCDF handles
    integer(kind_nf) :: ncid
    real, dimension(nk_stir) :: atmp

    if (proc0) then
       is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel)
       file_proc = get_file_proc(is_one_file_per_processor)
       ncid = neasyf_open (file_proc, "r")

       call verify_compatible_dim_sizes(ncid, "nk_stir", nk_stir)

       atmp = 0.
       call neasyf_read(ncid, "a_ant_r", atmp)
       a_ant = atmp

       call neasyf_read(ncid, "a_ant_i", atmp)
       a_ant = a_ant + zi * atmp

       call neasyf_read(ncid, "b_ant_r", atmp)
       b_ant = atmp

       call neasyf_read(ncid, "b_ant_i", atmp)
       b_ant = b_ant + zi * atmp

       call neasyf_close(ncid)
    end if
    call broadcast(a_ant)
    call broadcast(b_ant)
# else
    call mp_abort("Cannot load antenna amplitudes from restart without netcdf.", .true.)
    UNUSED_DUMMY(a_ant); UNUSED_DUMMY(b_ant); UNUSED_DUMMY(nk_stir)
# endif
  end subroutine init_ant_amp

  !> FIXME : Add documentation
  subroutine read_t0_from_restart_file (tstart)
# ifdef NETCDF
    use mp, only: proc0, broadcast
    use neasyf, only: neasyf_open, neasyf_read, neasyf_close
    use netcdf_utils, only: kind_nf
# else
    use mp, only: mp_abort
# endif
    implicit none
    real, intent (in out) :: tstart
# ifdef NETCDF
    character (run_name_size) :: file_proc
    logical :: is_one_file_per_processor
    ! NetCDF handles
    integer(kind_nf) :: ncid

    if (proc0) then
       is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel)
       file_proc = get_file_proc(is_one_file_per_processor)
       ncid = neasyf_open(file_proc, "r")
       call neasyf_read(ncid, "t0", tstart)
       call neasyf_close(ncid)
    end if

    call broadcast (tstart)
# else
    call mp_abort("Cannot read t0 from restart file in build with no netcdf.", .true.)
    UNUSED_DUMMY(tstart)
# endif
  end subroutine read_t0_from_restart_file
end module gs2_save