geometry.f90 Source File


Contents

Source Code


Source Code

!> Deals with most types of equilibria to generate GS2's expected geometrical parameters
!> Input:
!> rhoc :: minor radius of surface of interest (see irho below)
!> use_large_aspect :: default false.  Choose true to get "s-alpha" equilibria (see local_eq)
!> Equilibrium options : Choose 1 of:
!>          local_eq : Analytic parameterization of local equilibrium
!>          ppl_eq = .true. :: use PPL style NetCDF equilibrium (Menard)
!>          transp_eq = .true. :: use TRXPL style NetCDF equilibrium (Menard)
!>          gen_eq = .true. :: use GA style NetCDF equilibrium (TOQ)
!>          chs_eq = .true. ::  Use CHEASE generated equilibrium
!>          efit_eq = .true. :: use EFIT  equilibrium
!>          dfit_eq = .true. :: use dipole equilibrium
!> -------------------> [note: if any of the above except local_eq are true set eqfile = input file name]
!> irho :: choose flux surface label
!>         case (1) :: sqrt(toroidal flux)/sqrt(toroidal flux of LCFS)
!>         case (2) :: diameter/(diameter of LCFS).  recommended
!>         case (3) :: poloidal flux/(poloidal flux of LCFS)
!>         case (4) :: rho_mid? Only used for the vacuum ring dipole, dfit_eq = T
!>
!> equal_arc :: .true. makes grad_par coefficient a constant.
!>
!> New feature:  Given an actual equilibrium, allow s_hat and d pressure/d rho
!> to vary self-consistently (a la Greene and Chance, Bishop).  We use the
!> Bishop formalism [C. M. Bishop, et al., NF, Vol. 24, 1579 (1984).]
!>
!> bishop :: uses Bishop relations to generate coefficients.
!>        case(0) :: do not use Bishop relations -- primarily for testing
!>        case(1) :: use Bishop relations with actual equilibrium shat, p'
!>        case(2) :: use Bishop relations with new shat, and p' from p_prime_input
!>        case(3) :: use Bishop relations with new shat, L_p from below
!> s_hat_input :: new magnetic shear, rho/q*dq/drho
!> invLp_input :: new -1/pressure * d pressure/d rho
!>        case(4) :: use Bishop relations with new shat, beta' from below
!> s_hat_input :: new magnetic shear, rho/q*dq/drho
!> beta_prime_input :: new d beta/d rho
!>        case(5) :: use Bishop relations with new shat, alpha from below
!> s_hat_input :: new magnetic shear, rho/q*dq/drho
!> alpha_input :: new alpha = -q**2 * rmaj * d pressure/d rho
!>        case(6) :: use Bishop relations with equilibrium shat, beta' from below
!> beta_prime_input :: new d beta/d rho
!>        case(7) :: use Bishop relations with equilibrium shat, beta' from below
!> dp_mult :: d pressure/d rho = equilibrium gradient * dp_mult
!>         (In each case, the definition of rho is determined by irho above)
!>
!> nperiod :: number of cells of width 2 pi in theta.
!>
!> geoType :: determines the analytic geometry specification
!>         case (0) :: the tilted Miller geometry: generalization of Miller et al. PoP 1998
!>           (see Chapter 3.2.1 of J. Ball MIT Masters thesis)
!>         case (1) :: the global MHD solution geometry: calculated by inverting the large
!>           aspect ratio, constant-current MHD solution
!>           (see J. Ball "Optimized up-down asymmetry to drive fast intrinsic rotation in tokamak reactors")
!>         case (2) :: the generalized ellipticity geometry: the polar formula for an ellipse generalized
!>           to arbitrary shaping effect mode number (see Chapter 5.1.2 of J. Ball Oxford PhD thesis)
!>         case (3) :: the Fourier series geometry: Fourier decomposition of the minor radius as a function
!>           of poloidal angle (see Chapter 5.1.1 of J. Ball Oxford PhD thesis)
!> rmaj :: major radius of magnetic axis (in units of aminor) if not local_eq
!>      or major radius of center of flux surface of interest (units of aminor) otherwise
!> r_geo :: major radius of last closed flux surface or
!>       if local_eq: sets the reference magnetic field to be the value of the toroidal
!>       field at R=r_geo on the flux surface of interest
!> aSurf :: the real-space flux surface label of the surface with deltam and deltan shaping as well as
!>   a major radius of Rmaj (only used when geoType={1})
!> shift :: d rmaj/drho.  (Typically < 0.) (only used when geoType={0,2,3})
!> -> Alternative definitions:
!>   Raxis :: major radial location of the magnetic axis. (only used when geoType={1})
!> shiftVert :: d zmaj/drho. (< 0 for a vertical shift) (only used when geoType={0,2,3})
!> -> Alternative definitions:
!>   Zaxis :: axial location of the magnetic axis. (only used when geoType={1})
!> mMode :: the first poloidal mode number (only used when geoType={2,3})
!> nMode :: the second poloidal mode number (only used when geoType={2,3})
!> deltam :: magnitude of the first poloidal shaping effect (only used when geoType={2,3})
!> -> Alternative definitions:
!>   delta2 :: identical to akappa (only used when geoType={1})
!>   akappa :: flux surface elongation (only used when geoType={0})
!> deltan :: magnitude of the second poloidal shaping effect (only used when geoType={2,3})
!> -> Alternative definitions:
!>   delta3 :: magnitude of the triangulaity-like shaping effect (only used when geoType={1})
!>   tri :: flux surface triangulaity (only used when geoType={0})
!> deltampri :: d(deltam)/d(rhoc) (only used when geoType={2,3})
!> -> Alternative definitions:
!>   akappapri :: d(akappa)/d(rhoc) (only used when geoType={0})
!> deltanpri :: d(deltan)/d(rhoc) (only used when geoType={2,3})
!> -> Alternative definitions:
!>   tripri :: d(tri)/d(rhoc) (only used when geoType={0})
!> thetam :: tilt angle of the first poloidal shaping effect (only used when geoType={2,3})
!> -> Alternative definitions:
!>   theta2 :: tilt angle of the flux surface elongation (only used when geoType={1})
!>   thetak :: tilt angle of the flux surface elongation (only used when geoType={0})
!> thetan :: tilt angle of the first poloidal shaping effect (only used when geoType={2,3})
!> -> Alternative definitions:
!>   theta3 :: tilt angle of the triangularity-like shaping effect (only used when geoType={1})
!>   thetad :: tilt angle of the flux surface triangularity (only used when geoType={0})
!> qinp :: safety factor q at surface of interest (local_eq = T only)
!> shat :: d q/drho at surface of interest
!>
!> delrho :: numerical parameter.  Should be "small enough".  Typically 0.001 ok.
!>
!> force_sym :: true -> assume up-down symmetric equilibrium.
!> writelots :: .true. for more output
!>
!>  Use:
!> Include this module, set input variables appropriately and call eikcoefs.
!>
!> I do not recommend calling other routines, but a few are available for
!> diagnostic purposes.
module geometry
  use geo_utils, only: abstract_geo_type, EQFILE_LENGTH, flux_surface_type
  implicit none

  private
  public :: eikcoefs, pbarfun, pbarofrho, rpofrho, diameter, run_eikcoefs_and_resample
  public :: eqfile, eqnormfile, irho, surf, verb, big
  public :: ppl_eq, transp_eq, chs_eq, efit_eq, gs2d_eq, idfit_eq, dfit_eq, local_eq, gen_eq
  public :: force_sym, writelots, equal_arc, bishop, use_large_aspect
  public :: p_prime_input, invLp_input, beta_prime_input, alpha_input, dp_mult
  public :: geom, eikcoefs_output_type, finish_geometry

  class(abstract_geo_type), allocatable :: geom
  ! Start of inputs
  type(flux_surface_type) :: surf
  real :: p_prime_input, invLp_input, beta_prime_input, alpha_input, dp_mult
  integer :: irho, bishop
  logical :: gen_eq, efit_eq, ppl_eq, local_eq, chs_eq, dfit_eq, idfit_eq, gs2d_eq, transp_eq
  logical :: writelots, equal_arc, force_sym, use_large_aspect
  character(len = EQFILE_LENGTH) :: eqfile, eqnormfile
  ! End of inputs
  real :: rpmin, rpmax
  integer :: big = 8

  !> The single-period theta domain "half-grid" size. That is, the number of points in
  !> \((0, \pi]\). `theta(-nth:nth)` is the single central period \([-\pi, \pi]\).  For
  !> EFIT and local equilbria, this is set from
  !> [[theta_grid_parameters_config_type:ntheta]]; for all other numerical equilbria, it
  !> is determined by the data file. Always even.
  integer :: nth
  !> The full theta domain "half-grid" size. That is, the total size of the
  !> [[geometry:theta]] grid is `(2 * ntgrid) + 1`, with array bounds
  !> `theta(-ntgrid:ntgrid)`. Equal to `(2*nperiod - 1)*nth`
  integer :: ntgrid
  !> Verbosity of print statements | Should we use runtime_tests:verbosity?
  integer :: verb = 2
  !> If debug then print out lots of debug info.
  !> debug is true if (verb >= 3)
  logical :: debug = .false.

  !> A type to hold a collection of geometrical quantities as calculated by
  !> [[eikcoefs]]
  type eikcoefs_output_type
     integer :: ntheta
     real :: drhodpsi
     real :: kxfac
     real :: qsf
     real :: surfarea
     real :: dvdrhon
     real :: shat
     real :: dbetadrho
     real :: rhoc
     !> theta :: the theta grid of the results.  The final grid has
     !>          pi - 2*pi*nperiod <= theta <= 2*pi*nperiod - pi
     real, dimension(:), allocatable :: theta
     !> bmag :: |B|
     real, dimension(:), allocatable :: bmag
     !> gradpar :: coefficient of b_hat.grad operator
     real, dimension(:), allocatable :: gradpar
     !> grho :: grad(rho)
     real, dimension(:), allocatable :: grho
     !> cdrift :: part of coriolis drift operator independent of theta_0
     real, dimension(:), allocatable :: cdrift
     !> cvdrift  :: coefficient of v_par**2/Omega b_hat x (b_hat.grad b_hat).grad
     real, dimension(:), allocatable :: cvdrift
     !> gbdrift  :: coefficient of mu/Omega b_hat x (grad B).grad operator
     real, dimension(:), allocatable :: gbdrift
     !> cdrift0 :: part of coriolis drift operator proportional to theta_0
     real, dimension(:), allocatable :: cdrift0
     !> cvdrift0 :: part of total curvature drift operator proportional to theta_0
     real, dimension(:), allocatable :: cvdrift0
     !> gbdrift0 :: part of total grad B drift operator proportional to theta_0
     real, dimension(:), allocatable :: gbdrift0
     !> gds2 :: part of grad_perp**2 operator independent of theta_0
     real, dimension(:), allocatable :: gds2
     !> gds21 :: part of grad_perp**2 operator proportional to theta_0
     real, dimension(:), allocatable :: gds21
     !> gds22 :: part of grad_perp**2 operator proportional to theta_0**2
     real, dimension(:), allocatable :: gds22
     !> gds23 :: part of v_E . grad theta operator independent of theta_0
     real, dimension(:), allocatable :: gds23
     !> gds24 :: part of v_E . grad theta operator proportional to theta_0
     real, dimension(:), allocatable :: gds24
     !> gds24_noq :: gds24/dqdrp
     real, dimension(:), allocatable :: gds24_noq
     !> rplot :: R(theta) - major radius on theta grid
     real, dimension(:), allocatable :: rplot
     !> zplot :: Z(theta) - height on theta grid
     real, dimension(:), allocatable :: zplot
     !> aplot :: a(theta) - alpha-phi on theta grid
     real, dimension(:), allocatable :: aplot
     !> rprime :: derivative wrt rho of R(theta) - major radius on theta grid
     real, dimension(:), allocatable :: rprime
     !> zprime :: derivative  wrt rho Z(theta) - height on theta grid
     real, dimension(:), allocatable :: zprime
     !> aprime :: derivative  wrt rho a(theta) - alpha-phi on theta grid
     real, dimension(:), allocatable :: aprime
     !> bpol :: Poloidal magnetic field strength
     real, dimension(:), allocatable :: bpol
  end type eikcoefs_output_type

contains

  !> Run eikcoefs with `ntheta = ntheta_input` and then resample onto
  !> a theta grid with resolution `ntheta_target`.
  subroutine run_eikcoefs_and_resample(ntheta_input, ntheta_target, nperiod, results, job_id)
    implicit none
    !> This is the number of theta grid points per two pi period used in the
    !> geometry calculation (actually one less than this).
    integer, intent(in out) :: ntheta_input
    !> This is the desired number of theta grid points per two pi period
    !> in the resulting geometrical coefficients used elsewhere.
    integer, intent(in) :: ntheta_target
    integer, intent(in) :: nperiod
    type(eikcoefs_output_type), intent(out) :: results
    integer, intent(in), optional :: job_id
    type(eikcoefs_output_type) :: full_results
    real, dimension(:), allocatable :: input_index_space, target_index_space
    integer :: ig, ntgrid_target, ntgrid_input
    real :: spacing

    ! First run eikcoefs using requested resolution
    call eikcoefs(ntheta_input, nperiod, full_results, job_id)

    ntgrid_input = (size(full_results%theta) - 1) / 2
    ntgrid_target =(ntheta_target*(2*nperiod-1)) / 2

    ! Set scalar outputs
    results%ntheta = ntheta_target
    results%drhodpsi = full_results%drhodpsi
    results%kxfac = full_results%kxfac
    results%qsf = full_results%qsf
    results%surfarea = full_results%surfarea
    results%dvdrhon = full_results%dvdrhon
    results%shat = full_results%shat
    results%dbetadrho = full_results%dbetadrho
    results%rhoc = full_results%rhoc

    ! Make an index space array for the original theta
    ! grid. Spans [0,1] in size(theta)=2*ntgrid_input points
    allocate(input_index_space(-ntgrid_input:ntgrid_input))
    spacing = 1.0 / (2 * ntgrid_input)
    do ig = -ntgrid_input, ntgrid_input
       input_index_space(ig) = (ig + ntgrid_input) * spacing
    end do

    ! Make new theta grid
    allocate(target_index_space(-ntgrid_target:ntgrid_target))

    ! First make target index space grid
    spacing = 1.0 / (2 * ntgrid_target)
    do ig = -ntgrid_target, ntgrid_target
       target_index_space(ig) = (ig + ntgrid_target) * spacing
    end do

    ! Now interpolate theta grid onto target index space grid
    call resample( input_index_space, full_results%theta, target_index_space, results%theta)

    ! Now interpolate other output quantities onto target theta grid
    call resample(full_results%theta, full_results%bmag, results%theta, results%bmag)
    call resample(full_results%theta, full_results%gradpar, results%theta, results%gradpar)
    call resample(full_results%theta, full_results%grho, results%theta, results%grho)
    call resample(full_results%theta, full_results%cdrift, results%theta, results%cdrift)
    call resample(full_results%theta, full_results%cvdrift, results%theta, results%cvdrift)
    call resample(full_results%theta, full_results%gbdrift, results%theta, results%gbdrift)
    call resample(full_results%theta, full_results%cdrift0, results%theta, results%cdrift0)
    call resample(full_results%theta, full_results%cvdrift0, results%theta, results%cvdrift0)
    call resample(full_results%theta, full_results%gbdrift0, results%theta, results%gbdrift0)
    call resample(full_results%theta, full_results%gds2, results%theta, results%gds2)
    call resample(full_results%theta, full_results%gds21, results%theta, results%gds21)
    call resample(full_results%theta, full_results%gds22, results%theta, results%gds22)
    call resample(full_results%theta, full_results%gds23, results%theta, results%gds23)
    call resample(full_results%theta, full_results%gds24, results%theta, results%gds24)
    call resample(full_results%theta, full_results%gds24_noq, results%theta, results%gds24_noq)
    call resample(full_results%theta, full_results%rplot, results%theta, results%rplot)
    call resample(full_results%theta, full_results%zplot, results%theta, results%zplot)
    call resample(full_results%theta, full_results%aplot, results%theta, results%aplot)
    call resample(full_results%theta, full_results%rprime, results%theta, results%rprime)
    call resample(full_results%theta, full_results%zprime, results%theta, results%zprime)
    call resample(full_results%theta, full_results%aprime, results%theta, results%aprime)
    call resample(full_results%theta, full_results%bpol, results%theta, results%bpol)

  contains
    subroutine resample(x_orig, y_orig, x_new, y_new)
      use splines, only: new_spline, delete_spline, splint, spline
      real, dimension(:), intent(in) :: x_orig, y_orig
      real, dimension(-ntgrid_target:), intent(in) :: x_new
      real, dimension(:), allocatable, intent(out) ::  y_new
      type(spline) :: the_spline
      integer :: i, lower, upper

      ! Spline the old data
      the_spline = new_spline(x_orig, y_orig)
      lower = lbound(x_new, 1) ; upper = ubound(x_new, 1)

      allocate(y_new(lower:upper))

      do i = lower, upper
         y_new(i) = splint(x_new(i), the_spline)
      end do

      ! Get rid of the spline
      call delete_spline(the_spline)
    end subroutine resample
  end subroutine run_eikcoefs_and_resample

  !> FIXME : Add documentation
  subroutine eikcoefs (ntheta_input, nperiod, outputs, job_id_in)
    use constants, only: twopi
    use file_utils, only: open_output_file, close_output_file, finish_file_utils
    use file_utils, only: futils_initialized => initialized, init_file_utils
    use optionals, only: get_option_with_default
    implicit none
    !> Sets the number of theta grid points. If odd will be updated to be ntheta_input - 1
    integer, intent(in out) :: ntheta_input
    integer, intent(in) :: nperiod
    !> A type containing all of the outputs of eikcoefs
    type(eikcoefs_output_type), intent(out) :: outputs
    !> job_id is the job id in a Trinity or list mode run. Only for debug
    integer, intent(in), optional :: job_id_in
    logical :: start_file_utils, dummy_flag
    integer :: job_id, i, iBmax, fort11_unit
    real, allocatable, dimension(:) :: seik, rgrid, invR, trip, dsdthet, rmajor, dummy
    real, allocatable, dimension(:,:) :: thgrad, rpgrad, crpgrad, grads
    real, allocatable, dimension(:,:) :: gradrptot, gradstot, gradztot, gradthtot, gradbtot
    real :: shat, dbetadrho, dqdrp, qval, dum, bi
    real :: rp, pbar, drhodrp, drhodrhod, drhodpsi, dpsidrp, dpsidrho
    ! Used to translate definition of the poloidal angle such that the location
    ! of the maximium magnetic field is at pi (see section 3.2.3 of Ball MIT
    ! Masters thesis or section 2.4 of "GS2 analytic geometry specification")
    real :: thetaShift
    job_id = get_option_with_default(job_id_in, -1)
    start_file_utils = .not. futils_initialized

    if (start_file_utils) then
      ! Initialize file utils with a run name we provide
      ! and by setting trin_run prevent it from searching
      ! the command line arguments.
      write(*,*) 'WARNING: Automatically initializing file_utils&
      & in eikcoefs. This is a hack to &
      & support older programs like eiktest. You should&
      & initialize file utils before calling eikcoefs'
      call init_file_utils(dummy_flag, name='geometry_module', trin_run=.true., input=.false.)
    end if
    !Create module level output file and get unit
    call open_output_file(fort11_unit,".fort.11")

    !Compute the initial constants
    debug = (verb > 2)
    if (debug) writelots = .true.

    if (debug) write(6,*) "eikcoefs: local_eq=",local_eq, 'jid=',job_id

    if (local_eq) then
       if (irho /= 2) write(*,*) 'Forcing irho = 2 for local_eq'
       irho = 2
       if (bishop == 0) then
          write(*,*) 'Bishop = 0 not compatible with local_eq. Forcing bishop = 1'
          bishop = 1 !Not clear that this is a good choice for Miller
       end if
    end if

    if (debug) write(6,*) "eikcoefs: call check ",'jid=',job_id
    call check(gen_eq, efit_eq, ppl_eq, local_eq, dfit_eq, idfit_eq, chs_eq, gs2d_eq, transp_eq, bishop)

    call init_uniform_theta_grid(ntheta_input, outputs%theta, nperiod, ntgrid, nth)
    surf%nt = ntheta_input
    if (debug) write(6,*) 'eikcoefs: alloc_local_arrays'
    if (.not. allocated(seik)) call alloc_local_arrays(ntgrid)
    thetaShift = 0.0
    allocate(outputs%bmag(-ntgrid:ntgrid), outputs%bpol(-ntgrid:ntgrid))
    if (.not. allocated(geom)) call allocate_geom(geom)
    call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol)

    ! Start theta shift calc. This is used to modify the definition of
    ! theta such that the maxiumum magnetic field is at pi. This is
    ! necessary to correctly calculate bounce points. Can only shift with local_eq
    if (local_eq) then
       ! Locate the index of the maximum magnetic field in the range [-pi,pi)
       iBmax = maxloc(outputs%bmag(-nth:nth-1), dim = 1) - (nth + 1)

       ! If max is not at -pi and theta then shift theta.
       if (iBmax /= -nth) then
          thetaShift = outputs%theta(iBmax) - outputs%theta(-nth)
          call geom%finalise !Clean up geometry to allow reinitialisation
          write(*,*) "Bmax is not located at +/-pi, translating theta by",thetaShift," so that it is..."
          call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol)
       end if
    end if

    !     compute |grad rho|:
    allocate(outputs%grho(-ntgrid:ntgrid))
    outputs%grho(-nth:nth) = hypot(rpgrad(-nth:nth,1), rpgrad(-nth:nth,2)) * abs(drhodrp)
    call periodic_copy(outputs%grho, 0., nperiod)
    if (force_sym) call sym(outputs%grho, 0, ntgrid)

    if (writelots) write(fort11_unit,*) 'dpsidrp= ',dpsidrp
    drhodpsi = drhodrp / dpsidrp ; dpsidrho = 1.0 / drhodpsi

    !     compute derivatives of eikonal, grad B and grad Z
    if (bishop /= 0) then
       call bishop_grads(surf%r, surf%rmaj, outputs%bmag, outputs%bpol)
    else
       call bishop_large_aspect_ratio_approximation(surf%r, surf%dr, ntheta_input)
    end if

    if (writelots) then
       write(fort11_unit,*) 'i,theta,grads1,2=  '
       do i = -ntgrid, ntgrid
          write(fort11_unit,*) i, outputs%theta(i), grads(i, 1), grads(i, 2)
       end do
    end if
    !     compute total grad S including toroidal part
    call gradstottgrid(invR, grads, gradstot)

    gradrptot(-nth:nth, 1:2) = rpgrad(-nth:nth, :)
    call periodic_copy(gradrptot(:, 1), 0., nperiod)
    call periodic_copy(gradrptot(:, 2), 0., nperiod)
    gradrptot(:, 3) = 0.

    ! gradthtot is grad theta since ballooning, theta goes from
    ! -ntheta*nperiod -> ntheta*nperiod, but grad theta is periodic in theta
    gradthtot(-nth:nth, 1:2) = thgrad(-nth:nth, :)
    call periodic_copy(gradthtot(:, 1), 0., nperiod)
    call periodic_copy(gradthtot(:, 2), 0., nperiod)
    gradthtot(:, 3) = 0.

    !     compute the magnitude squared of grad S
    allocate(outputs%gds2(-ntgrid:ntgrid), outputs%gds21(-ntgrid:ntgrid), &
         outputs%gds22(-ntgrid:ntgrid), outputs%gds23(-ntgrid:ntgrid), &
         outputs%gds24(-ntgrid:ntgrid), outputs%gds24_noq(-ntgrid:ntgrid))
    call calculate_metric_terms(gradstot, gradrptot, gradthtot, dpsidrho, dqdrp, &
         qval, surf%r, outputs%bmag, force_sym, outputs%gds2, outputs%gds21, outputs%gds22, &
         outputs%gds23, outputs%gds24, outputs%gds24_noq)

    !     compute curvature, grad b, and coriolis drift terms
    allocate(outputs%gbdrift(-ntgrid:ntgrid), outputs%gbdrift0(-ntgrid:ntgrid))
    allocate(outputs%cvdrift(-ntgrid:ntgrid), outputs%cvdrift0(-ntgrid:ntgrid))
    allocate(outputs%cdrift(-ntgrid:ntgrid), outputs%cdrift0(-ntgrid:ntgrid))
    call drift(rgrid, rp, gradstot, gradrptot, gradztot, gradbtot, dqdrp, dpsidrho, dpsidrp,&
         outputs%bmag, invR, outputs%theta, dbetadrho, bi, force_sym, nperiod, &
         outputs%gbdrift, outputs%gbdrift0, outputs%cvdrift, outputs%cvdrift0, &
         outputs%cdrift, outputs%cdrift0)

    !Recompute the magnetic field along theta if use_large_aspect - should this come earlier?
    if (use_large_aspect) then
       outputs%bmag(-nth:nth)= 1 / (1 + rgrid(-nth:nth) * cos(outputs%theta(-nth:nth)) / surf%rmaj)
       call periodic_copy(outputs%bmag, 0., nperiod)
    end if
    if (force_sym) call sym(outputs%bmag, 0, ntgrid)

    !     compute the grad-parallel coefficient
    allocate(outputs%gradpar(-ntgrid:ntgrid))
    if (dfit_eq) then
       outputs%gradpar = trip / outputs%bmag
    else
       outputs%gradpar(-nth:nth) = -bi / (rmajor(-nth:nth)**2 * outputs%bmag(-nth:nth) * dsdthet(-nth:nth))
       call periodic_copy(outputs%gradpar, 0., nperiod)
    end if
    if (force_sym) call sym(outputs%gradpar, 0, ntgrid)

    ! Note : We should call this before we potentially redefine the theta grid
    ! due to equal_arc below as some specific equilibrium types rely on theta
    ! having a consistent definition with that used to calculate rgrid initially.
    ! In particular this is important for Efit and GS2D equilibria (in eeq). Other
    ! equilibrium types use interpolation so shouldn't depend on the specific theta grid
    ! (although it's possible there could still be some interaction, given that one of the
    ! outputs of the below 'aplot' is given by `seik` which is evaluated on the original
    ! grid rather than the "new" one defined below).
    ! @note Check if restriction on cartesian geo types still holds.
    allocate(outputs%rplot(-ntgrid:ntgrid), outputs%rprime(-ntgrid:ntgrid))
    allocate(outputs%zplot(-ntgrid:ntgrid), outputs%zprime(-ntgrid:ntgrid))
    allocate(outputs%aplot(-ntgrid:ntgrid), outputs%aprime(-ntgrid:ntgrid))
    call plotdata (rgrid, seik, grads, dpsidrho, surf%dr, outputs%theta, &
         outputs%rplot, outputs%zplot, outputs%aplot, &
         outputs%rprime, outputs%zprime, outputs%aprime)

    ! In the lines below we redefine theta to be theta', and gradpar
    ! to be gradpar'. This means that the values of theta in the gs2
    ! output file are not values of the poloidal angle but values of theta_prime
    if (equal_arc) then
       allocate(dummy(-ntgrid:ntgrid))
       call arclength (ntheta_input, nperiod, outputs%theta, outputs%gradpar, dummy)
       outputs%theta = dummy
    end if

    call calculate_surface_area(outputs%theta, outputs%gradpar, outputs%bmag, outputs%grho, &
         dpsidrho, surf%r, nth, outputs%surfarea, outputs%dvdrhon)
    if (writelots) write (fort11_unit,*) 'surfarea=', outputs%surfarea

    !Close output file
    call close_output_file(fort11_unit)

    outputs%kxfac = abs(qval * dpsidrho / surf%r)

    if (qval == 0. .and. (dfit_eq .or. idfit_eq)) then
       shat = max(1.e-7, shat)
       !RN> I think kxfac should include the sign. However, from GS2 convention,
       !RN> eventually kxfac = abs(kxfac). See get_grids in theta_grid. Even if the sign is
       !RN> effective, probably it has no effect due to symmetry. I just leave the sign here.
       outputs%kxfac = - dpsidrho
       outputs%gds22 = (outputs%grho * dpsidrho)**2
    end if

    if (start_file_utils) call finish_file_utils

    ! Pack data into output type
    outputs%ntheta = ntheta_input ; outputs%drhodpsi = drhodpsi ; outputs%rhoc = surf%r
    outputs%qsf = qval ; outputs%shat = shat ; outputs%dbetadrho = dbetadrho

  contains
    !> Initialises the specific geometry backend
    subroutine initialise_geometry_and_get_bmag(geom, rhoc, thetaShift, theta, bmag, bpol)
      use geo_utils, only: geo_input_type
      class(abstract_geo_type), intent(in out) :: geom
      real, intent(in) :: rhoc, thetaShift
      real, dimension(-ntgrid:), intent(in) :: theta
      real, dimension(-ntgrid:), intent(out) :: bmag, bpol
      type(geo_input_type) :: geom_inputs
      real :: avgrmid, rbar, B_T0
      integer :: fort78_unit, fort99_unit
      ! Initialise the specific geometry types (mostly reading data
      ! from file and calculating derived quantities).
      geom_inputs%surf = surf ; geom_inputs%thShift = thetaShift ; geom_inputs%big = big
      geom_inputs%eqfile = eqfile ; geom_inputs%eqnormfile = eqnormfile ; geom_inputs%theta = theta
      call geom%initialise(geom_inputs, rpmin, rpmax, surf%rmaj, B_T0, avgrmid)
      if (debug) write(6,*) "eikcoefs: rpmin,rpmax=", rpmin, rpmax
      if (.not. local_eq) surf%r_geo = geom%rcenter(rpmax) !Only used for writelots/checks
      if (surf%r_geo > surf%rmaj + 1.e-5) then
         write(6,*)
         write(6,*) 'Warning: the center of the LCFS is outboard of the '
         write(6,*) 'center of the reference surface.  This is probably wrong.'
         write(6,*)
      end if

      if (writelots) then
         write(fort11_unit,*) 'All lengths normalized to the avg midplane minor radius'
         if (irho /= 2) then
            write(fort11_unit,*) 'except rho, which always goes from 0 to 1.'
            write(fort11_unit,*) 'The rho normalization is determined by irho.'
         end if
         write(fort11_unit,*) 'avgrmid = ',avgrmid,' meters'
         if (.not. local_eq) write(fort11_unit,*) 'R_mag = ',surf%rmaj
         if (.not. local_eq) write(fort11_unit,*) 'B_T0 = ',B_T0
      end if

      if (debug) write(6,*) "eikcoefs: find rp ", 'jid=',job_id
      rp = rpofrho(rhoc) ; drhodrp = drho_drp(rp, surf%dr)
      if (debug) write(6,*) "eikcoefs: rp=",rp, 'jid=',job_id

      if (debug) write(6,*) "eikcoefs: find rgrid"
      call rtgrid(rgrid, rp, theta)
      call periodic_copy(rgrid, 0.0, nperiod)

      pbar = pbarfun(rgrid(0), theta(0))
      bi = geom%btori(pbar) ! Get the toroidal field function (i.e. fp = R B_T)

      if (writelots) then
         write(fort11_unit,*) 'rp of rho = ',rp
         write(fort11_unit,*) 'Rcenter(rho) = ',geom%rcenter(rp)
         write(fort11_unit,*) 'R_geo = ',surf%r_geo,' !or ',bi

         call open_output_file(fort78_unit,".fort.78")
         call open_output_file(fort99_unit,".fort.99")
         write(fort78_unit,*) '# theta_pol,  R,   Z'
         do i = -nth, nth
            write(fort78_unit,*) theta(i), geom%Rpos(rgrid(i), theta(i)), &
                 geom%Zpos(rgrid(i), theta(i))
            write(fort99_unit,*) 'i=', i, ' thet= ', theta(i), ' rgrid= ', rgrid(i)
         end do
         call close_output_file(fort78_unit) ; call close_output_file(fort99_unit)

         rbar = 0.
         do i = -nth, nth - 1
            rbar = rbar + rgrid(i) * (theta(i + 1) - theta(i))
         end do
         rbar = rbar / twopi
         write (*,*) 'r_bar = ', rbar
         write (*,*)

         write(fort11_unit,*) 'drhodrp=', drhodrp
         if (irho == 1) then
            drhodrhod = drho_drhod(rp, drhodrp, surf%dr)
            write(fort11_unit,*) 'drho_drhod=', drhodrhod
         end if
      end if

      if (debug) write(fort11_unit,*) 'eikcoefs: various logs', gen_eq, ppl_eq, chs_eq, efit_eq, dfit_eq, idfit_eq, gs2d_eq, transp_eq

      if (debug) write(6,*) "eikcoefs: call rmajortgrid"
      call rmajortgrid(rgrid, theta, rmajor) ; call periodic_copy(rmajor, 0., nperiod)
      if (use_large_aspect) then
         invR = 1.0 / surf%rmaj
      else
         call invRtgrid(rgrid, theta, invR) ; call periodic_copy(invR, 0., nperiod)
      end if

      !     compute gradient of rp in R, Z, phi (cartesian - crpgrad)
      ! Note that rpgrad is given in (R,Z,phi) coordinates if bishop == 0
      ! and (rho,l,phi) coordinates otherwise
      ! Also, rpgrad is grad psi if not local_eq and grad rho if local_eq -- MAB
      ! When bishop > 0
      ! rpgrad(:,1) is |grad rp| ; rpgrad(:,2) is 0
      ! crpgrad(:,1) is drp/dR ; crpgrad(:,2) is drp/dZ
      !
      ! bishop = 0 is high aspect ratio coefficients for debugging -- EGH
      call geom%gradient(rgrid, theta, crpgrad, 'P', dum, nth, ntgrid, .false.)
      call periodic_copy(crpgrad(:,1), 0., nperiod)
      call periodic_copy(crpgrad(:,2), 0., nperiod)

      !     compute coordinate gradients
      call geom%gradient(rgrid, theta, rpgrad, 'P', dum, nth, ntgrid, bishop /= 0)
      call geom%gradient(rgrid, theta, thgrad, 'T', dum, nth, ntgrid, bishop /= 0)
      call tripprod2dtgrid(rpgrad, thgrad, invR, trip)
      call periodic_copy(trip, 0., nperiod)

      !     compute eikonal S
      if (local_eq) qval = geom%qfun(pbar) !leq:qfun ignores input pbar
      if (debug) write(6,*) "eikcoefs: call eikonal"
      call eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp)
      if (debug) write(6,*) "eikcoefs: done eikonal"
      if (writelots) write(fort11_unit,*) 'q= ', qval
      if (local_eq) trip = trip * dpsidrp

      if (bishop > 0) then !Following bits _could_ be calculated for all bishop
         call B_pol(invR, crpgrad, dpsidrp, bpol, nth) !<Should this use rpgrad for local?
         call periodic_copy(bpol, 0., nperiod)
         call B_mod(invR, bpol, bmag, bi, nth)
         call periodic_copy(bmag, 0., nperiod)
      end if !Note bmag not set for bishop == 0
    end subroutine initialise_geometry_and_get_bmag

    subroutine bishop_grads(rhoc, rmaj, bmag, bpol)
      real, intent(in) :: rhoc, rmaj
      real, dimension(-ntgrid:), intent(in) :: bmag, bpol
      logical, parameter :: save_bishop = .true.
      integer :: bishop_unit
      real :: a_b, b_b, c_b, dp, dp_new, di, di_new, s_hat, s_hat_new, pressure, tmp
      real, allocatable, dimension(:) :: a_bish, b_bish, c_bish, da_bish, db_bish, dc_bish, &
           zoftheta, dzdl, dsdl, ltheta, th_bish, invRpol
      allocate(a_bish(-ntgrid:ntgrid), b_bish(-ntgrid:ntgrid), c_bish(-ntgrid:ntgrid), &
           da_bish(-ntgrid:ntgrid), db_bish(-ntgrid:ntgrid), dc_bish(-ntgrid:ntgrid), &
           zoftheta(-ntgrid:ntgrid), dzdl(-ntgrid:ntgrid), dsdl(-ntgrid:ntgrid), &
           ltheta(-ntgrid:ntgrid), th_bish(-ntgrid:ntgrid), invRpol(-ntgrid:ntgrid))

      call th_bishop(crpgrad, th_bish, nth)
      call periodic_copy(th_bish, 0., nperiod)
      call loftheta(rgrid, outputs%theta, ltheta)
      call periodic_copy(ltheta, ltheta(nth) - ltheta(-nth), nperiod)
      call inv_Rpol(outputs%theta, th_bish, ltheta, invRpol, nth)
      call periodic_copy(invRpol, 0., nperiod)
      if (debug) write(6,*) "eikcoefs:  gradients"
      if (debug) write(6,*) 'drhodpsi', drhodpsi
      if (debug) write(6,*) 'bi', bi
      if (verb > 5) write(6,*) 'rmajor', rmajor(:)
      if (verb > 5) write(6,*) 'bpol', bpol(:)
      if (verb > 5) write(6,*) rpgrad(-nth:nth,1)
      if (verb > 5) write(6,*) rpgrad(-nth:nth,2)
      if (verb > 5) write(6,*) thgrad(-nth:nth,1)
      if (verb > 5) write(6,*) thgrad(-nth:nth,2)
      if (verb > 5) write(6,*) 1.0 / trip(-nth:nth)
      if (verb > 5) write(6,*) -1.0 / (bpol(-nth:nth) * invRpol(-nth:nth))
      if (debug) write(6,*) "eikcoefs: end gradients"

      da_bish = (1 / rmajor**2) + (bi / (bpol * rmajor**2))**2
      db_bish = bi / (bpol * rmajor)**2
      dc_bish = bi * (sin(th_bish) + rmajor * invRpol) / (-bpol * rmajor**4)

      da_bish = da_bish / trip
      db_bish = db_bish / trip
      dc_bish = dc_bish / trip

      call integrate(da_bish, outputs%theta, a_bish, ntgrid)
      call integrate(db_bish, outputs%theta, b_bish, ntgrid)
      call integrate(dc_bish, outputs%theta, c_bish, ntgrid)

      a_b = a_bish(nth) - a_bish(-nth)
      b_b = b_bish(nth) - b_bish(-nth)
      c_b = c_bish(nth) - c_bish(-nth)

      a_bish = - a_bish * bpol * rmajor
      b_bish = - b_bish * bpol * rmajor
      c_bish = - c_bish * bpol * rmajor

      ! find shat, pressure' and report them
      dp = geom%dpfun(pbar) ; di = geom%dbtori(pbar)

      if (abs(qval) > epsilon(0.)) then
         tmp = rhoc / (qval * twopi * drhodpsi)
         s_hat = tmp * (a_b * di + b_b * dp + c_b * 2) !< Calculate the shear from equib.
         if (local_eq) s_hat = surf%shat !Can't determine this by calculation in local_eq
         if (debug) write(6,*) "a_b (f'), b_b (p'), c_b terms =", a_b*di, b_b*dp, c_b*2
         if (debug) write(6,*) 'qval, rho, drhodpsi, s_hat =', qval, rhoc, drhodpsi, s_hat
      else
         tmp = 1. ; s_hat = 0.
      end if

      write(fort11_unit,*)
      write(fort11_unit,*) 'Quantity, equilibrium value, value used'

      ! Decide how we want to set s_hat_new in the default case
      s_hat_new = s_hat
      select case (bishop)
      case (2, 3, 4, 5, 8, 9) !1, 6, 7 keep the initial value
         s_hat_new = surf%shat
      end select

      ! Decide how we want to set dp_new
      dp_new = dp
      select case (bishop)
      case (2)
         dp_new = p_prime_input
      case (3)
         pressure = geom%pfun(pbar)
         dp_new = -invLp_input * pressure * drhodpsi
      case (4, 6, 9)
         dp_new = 0.5 * beta_prime_input * drhodpsi
      case (5)
         dp_new = -alpha_input * drhodpsi / (qval**2 * rmaj)
      case (7)
         dp_new = dp * dp_mult
      case (8)
         dp_new = 0.5 * beta_prime_input * drhodpsi * dp_mult
      end select

      !Report on updated dp
      select case(bishop)
      case (2)
         write(fort11_unit,*) 'p_prime = ',dp,', ',dp_new
         if (writelots) write(*,*) 'p_prime = ',dp,', ',dp_new
      case (3)
         write(fort11_unit,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', &
              -dp_new / (drhodpsi *pressure)
         if (writelots) write(*,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', &
              -dp_new / (drhodpsi * pressure)
      case (5)
         write(fort11_unit,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', &
              -qval**2 * rmaj * dp_new / drhodpsi
         if (writelots) write(*,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', &
              -qval**2 * rmaj * dp_new / drhodpsi
      case default
         write(fort11_unit,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi
         if (writelots) write(*,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi
      end select

      write(fort11_unit,*) 's_hat ',s_hat,', ',s_hat_new
      if (writelots) write(*,*) 's_hat ',s_hat,', ',s_hat_new
      write(fort11_unit,*)

      if ((.not. local_eq) .and. (bishop == 1 .or. bishop > 9)) then
         di_new = di
      else
         di_new = (s_hat_new / tmp - 2 * c_b - dp_new * b_b) / a_b
      end if

      ! Calculate grad(S)
      ! here the coordinates are rho and l
      call gradl(ltheta, seik, dSdl, twopi * qval, nth)
      call periodic_copy(dSdl, 0., nperiod)
      grads(:,1) = a_bish * di_new + b_bish * dp_new + c_bish * 2
      grads(:,2) = dSdl

      ! Calculate grad(Z)
      ! obtain Z(theta) and hence dZ/dl
      call Ztgrid(rgrid, outputs%theta, Zoftheta)
      call periodic_copy(Zoftheta, 0., nperiod)
      call gradl(ltheta, Zoftheta, dZdl, 0., nth)
      call periodic_copy(dZdl, 0., nperiod)

      ! gradztot(:,1) is grad Z . rhohat
      gradztot(:,1) = crpgrad(:,2) / outputs%grho
      ! gradztot(:,2) is grad Z . lhat
      gradztot(:,2) = dZdl
      ! gradztot(:,3) is grad Z . phihat
      gradztot(:,3) = 0.0

      shat = s_hat_new ! Output
      dqdrp = s_hat_new * qval * drhodrp / rhoc ! Output
      dbetadrho = 2 * dp_new * dpsidrho ! Output

      call bishop_gradB(bmag, bpol, invRpol, th_bish, ltheta, rmajor, dp_new, bi, gradbtot)

      if (save_bishop) then
         call open_output_file(bishop_unit, '.bishop')
         write(bishop_unit, fmt="('#di, di_new = ',F12.7,' ',F12.7)") di, di_new
         write(bishop_unit, fmt="('#dp, dp_new = ',F12.7,' ',F12.7)") dp, dp_new
         write(bishop_unit, fmt="('#bi =', F12.7)") bi
         write(bishop_unit, fmt="('#drdpsi =', F12.7)") drhodpsi
         write(bishop_unit, fmt="('#a_b (fp)=', F20.10, ' ,b_b (pp)=', F20.10,' ,c_b  =', F20.10,' ,s_hat =', F20.10)") a_b*di, b_b*dp, c_b, s_hat
         write(bishop_unit, fmt ="('#  theta, a_bish,  b_bish,  c_bish, rmajor, Z, Bpol, invRpol, th_bish')")
         do i = -ntgrid, ntgrid
            write(bishop_unit, '( 10(2X, F25.17) )') outputs%theta(i), a_bish(i), b_bish(i), &
                 c_bish(i), rmajor(i), Zoftheta(i),  bpol(i), invRpol(i), th_bish(i)
         end do
         call close_output_file(bishop_unit)
      end if
    end subroutine bishop_grads

    !> @Note local_eq _must_ be false if we end in this routine. We make use of this.
    subroutine bishop_large_aspect_ratio_approximation(rhoc, dr, ntheta)
      real, intent(in) :: rhoc, dr
      integer, intent(in) :: ntheta
      real :: delrp, rp1, rp2, qval1, qval2, bi1, bi2, dqdrho, pbar_1, pbar_2
      real, dimension(:), allocatable :: seik1, seik2, dsdrp
      integer :: itot, k
      !     rp derivative
      delrp = dr / drhodrp
      rp1 = rp - delrp ; rp2 = rp + delrp
      pbar_1 = pbarfun(rp1, outputs%theta(0)) ; pbar_2 = pbarfun(rp2, outputs%theta(0))
      ! Calculate modified btori -- note we use the fact that rgrid == rp in
      ! seikon to avoid having to form rgrid1/rgrid2 here.
      bi1 = geom%btori(pbar_1) ; bi2 = geom%btori(pbar_2)

      allocate(seik1(-ntgrid:ntgrid), seik2(-ntgrid:ntgrid), dsdrp(-ntgrid:ntgrid))
      call seikon(rp1, outputs%theta, invR, qval1, bi1, nperiod, seik1)
      call seikon(rp2, outputs%theta, invR, qval2, bi2, nperiod, seik2)

      dqdrho = (qval2 - qval1) / (2 * dr)
      shat = dqdrho * rhoc / qval ! Output
      dqdrp = dqdrho * drhodrp ! Output
      dbetadrho = 2 * geom%dpfun(pbar) / drhodrp ! Output

      if (writelots) then
         write (fort11_unit,*) 'dqdrho=',dqdrho,' qval1,2= ',qval1, qval2
         write (fort11_unit,*) 's_hat= ',shat
         write (fort11_unit,*) 'i,rgrid,rgrid1,rgrid2:'
         do i = -nth, nth
            write (fort11_unit,*) i, rgrid(i), rp1, rp2
         end do
      end if

      do i = -nth, nth
         dsdrp(i) = (seik2(i) - seik1(i)) / (2 * delrp)
      end do
      call periodic_copy(dsdrp, -twopi * dqdrp, nperiod)

      !     compute gradient of S
      ! this is really gradient of alpha
      ! note that coordinates are R and Z (not rho and l, as in bishop case)
      do k = -nperiod+1, nperiod-1
         do i = -nth, nth
            itot = i + k * ntheta
            grads(itot, :) = rpgrad(i, :) * dsdrp(itot) + thgrad(i, :) * dsdthet(i)  ! Output
         end do
      end do
      ! gradztot(:,1) is grad Z . Rhat ! Output
      gradztot(:,1) = 0.
      ! gradztot(:,2) is grad Z . Zhat ! Output
      gradztot(:,2) = 1.
      ! gradztot(:,3) is grad Z . phihat ! Output
      gradztot(:,3) = 0.

      call geom%gradient(rgrid, outputs%theta, gradbtot(:, 1:2), 'B', rp, nth, ntgrid, bishop /= 0)
      gradbtot(:, 3) = 0.0
    end subroutine bishop_large_aspect_ratio_approximation

    !> FIXME : Add documentation
    subroutine alloc_local_arrays(n)
      implicit none
      integer, intent(in) :: n

      allocate(seik   (-n:n), &
           dsdthet    (-n:n), &
           invR       (-n:n), &
           trip       (-n:n), &
           rgrid      (-n:n), &
           rmajor     (-n:n))

      allocate(thgrad (-n:n,2), &
           rpgrad     (-n:n,2), &
           crpgrad    (-n:n,2), &
           grads      (-n:n,2))

      allocate(gradrptot (-n:n,3), &
           gradstot   (-n:n,3), &
           gradztot   (-n:n,3), &
           gradthtot  (-n:n,3), & ! MAB
           gradbtot   (-n:n,3))
    end subroutine alloc_local_arrays
  end subroutine eikcoefs

  !> Allocate the geom instance to the correct type
  pure subroutine allocate_geom(geom)
    use geq, only: geq_type
    use leq, only: leq_type
    use ceq, only: ceq_type
    use peq, only: peq_type, teq_type
    use eeq, only: eeq_type, gs2d_type
    use deq, only: deq_type
    use ideq, only: ideq_type
    use geo_utils, only: abstract_geo_type
    implicit none
    class(abstract_geo_type), allocatable, intent(in out) :: geom
    if (gen_eq) then
       allocate(geq_type :: geom)
    else if (chs_eq) then
       allocate(ceq_type :: geom)
    else if (gs2d_eq) then
       allocate(gs2d_type :: geom)
    else if (efit_eq) then
       allocate(eeq_type :: geom)
    else if (transp_eq) then
       allocate(teq_type :: geom)
    else if (ppl_eq) then
       allocate(peq_type :: geom)
    else if (dfit_eq) then
       allocate(deq_type :: geom)
    else if (idfit_eq) then
       allocate(ideq_type :: geom)
    else if (local_eq) then
       allocate(leq_type :: geom)
    end if
  end subroutine allocate_geom

  pure subroutine calculate_metric_terms(gradstot, gradrptot, gradthtot, dpsidrho, dqdrp, &
       qval, rhoc, bmag, force_sym, gds2, gds21, gds22, gds23, gds24, gds24_noq)
    real, dimension(-ntgrid:, :), intent(in) :: gradstot, gradrptot, gradthtot
    real, intent(in) :: dpsidrho, dqdrp, qval, rhoc
    real, dimension(-ntgrid:), intent(in) :: bmag
    logical, intent(in) :: force_sym
    real, dimension(-ntgrid:), intent(out) :: gds2, gds21, gds22, gds23, gds24, gds24_noq
    real, dimension(:), allocatable :: gdsdum1, gdsdum2, gdsdum3, gdsdum4, gdsdum5
    allocate(gdsdum1(-ntgrid:ntgrid), gdsdum2(-ntgrid:ntgrid), gdsdum3(-ntgrid:ntgrid))
    call dottgridf(gradstot, gradstot, gdsdum1)
    call dottgridf(gradstot, gradrptot, gdsdum2)
    call dottgridf(gradrptot, gradrptot, gdsdum3)
    gds2  = gdsdum1            *dpsidrho**2
    gds21 = gdsdum2*dqdrp      *dpsidrho**2
    gds22 = gdsdum3*dqdrp*dqdrp*dpsidrho**2

    ! v_E . grad theta = i*ky*phi*(rho_{r}*v_{t,r}/a**2)*(gds23 + gsd24*theta0)
    ! calculation of gds23 and gds24 assumes rp=rho, which is not necessarily true
    ! for numerical equilibria
    allocate(gdsdum4(-ntgrid:ntgrid), gdsdum5(-ntgrid:ntgrid))
    call dottgridf(gradrptot, gradthtot, gdsdum4) !grad rho dot grad theta
    call dottgridf(gradstot, gradthtot, gdsdum5) !grad alpha dot grad theta
    gds23 = (gdsdum4 * gds2 - gdsdum5 * gdsdum2 * dpsidrho**2) / bmag**2
    gds24_noq = (gdsdum4 * gdsdum2 - gdsdum5 * gdsdum3) * (dpsidrho / bmag)**2
    gds24 = gds24_noq * dqdrp
    gds24_noq = gds24_noq * qval / rhoc

    if (force_sym) then
       call sym(gds2, 0, ntgrid)
       call sym(gds21, 1, ntgrid)
       call sym(gds22, 0, ntgrid)
       ! not sure about the following yet -- MAB
       ! call sym(gds23, 0, ntgrid); call sym(gds24, 1, ntgrid)
    end if
  end subroutine calculate_metric_terms

  subroutine calculate_surface_area(theta, gradpar, bmag, grho, dpsidrho, rhoc, nth, &
       surfarea, dvdrhon)
    use constants, only: twopi
    use file_utils, only: open_output_file, close_output_file
    real, dimension(-ntgrid:), intent(in) :: theta, gradpar, bmag, grho
    real, intent(in) :: dpsidrho, rhoc
    integer, intent(in) :: nth
    real, intent(out) :: surfarea, dvdrhon
    real, dimension(:), allocatable :: jacob, ds, dummy
    integer :: i, fort25_unit
    ! compute grad rho*Jacobian of the transformation to the new coordinates:
    ! Axisymmetry assumed to get this (nu indep. of phi):
    allocate(jacob(-ntgrid:ntgrid), ds(-ntgrid:ntgrid), dummy(-ntgrid:ntgrid))
    do i = -nth, nth
       jacob(i) = abs(dpsidrho / (gradpar(i)*bmag(i)))
       ds(i) = grho(i) * jacob(i)
    end do

    ! compute surface area from A=Int(R sqrt(r**2+drdtheta**2) dtheta dphi)
    !  does not work with equal_arc grid
    ! compute the surface area from A=Int(J |grad rho| dtheta dalpha)
    call integrate(ds, theta, dummy, nth)
    surfarea = twopi * (dummy(nth) - dummy(-nth))

    ! compute dV/drhon = 2 pi Int(J dtheta)
    call integrate(jacob, theta, dummy, nth)
    dvdrhon = twopi * (dummy(nth) - dummy(-nth))

    call open_output_file(fort25_unit,".fort.25")
    write(fort25_unit,*) rhoc, f_trap(bmag(-nth:nth), theta(-nth:nth), jacob(-nth:nth))
    call close_output_file(fort25_unit)
  end subroutine calculate_surface_area

  !> Computes something like x cross y dot grad zeta.
  pure subroutine tripprod2dtgrid(x, y, invR, val)
    implicit none
    real, dimension(-ntgrid:), intent (in) :: invR
    real, dimension(-ntgrid:), intent (out) :: val
    real, dimension(-ntgrid:, :), intent (in) :: x, y
    integer :: i
    ! factor of 1/R from grad zeta
    do i = -nth, nth
       val(i) = (x(i, 1) * y(i, 2) - x(i, 2) * y(i, 1)) * invR(i)
    end do
  end subroutine tripprod2dtgrid

  !> Calculate rgrid for the case of eeq and deq: in this case, rgrid
  !> is the distance to the magnetic axis as a function of theta.
  subroutine rtgrid(rgrid, rp, theta)
    implicit none
    real, dimension(-ntgrid:), intent (out) :: rgrid
    real, intent (in) :: rp
    real, dimension(-ntgrid:), intent (in) :: theta
    integer :: i
    do i = -nth, nth
       rgrid(i) = geom%rfun(rp, theta(i))
    end do
  end subroutine rtgrid

  !> Calculates rmajor = R(theta), where R is the distance
  !> to the central axis (rmajor not to be confused with rmaj, r_geo etc)
  !> rgrid is the submodule radial grid (r_circ for EFIT & deq, rp otherwise)
  subroutine rmajortgrid(rgrid, theta, rmajor)
    implicit none
    real, dimension(-ntgrid:), intent(in) :: rgrid, theta
    real, dimension(-ntgrid:), intent(out) :: rmajor
    integer :: i
    do i = -nth, nth
       rmajor(i) = geom%Rpos(rgrid(i), theta(i))
    end do
  end subroutine rmajortgrid

  !> Calculates invR = 1/R(theta), where R is the distance
  !> to the central axis (rmajor not to be confused with rmaj, r_geo etc)
  !> rgrid is the submodule radial grid (r_circ for EFIT & deq, rp otherwise)
  subroutine invRtgrid(rgrid, theta, invR)
    implicit none
    real, dimension(-ntgrid:), intent(in) :: rgrid, theta
    real, dimension(-ntgrid:), intent(out) :: invR
    integer :: i
    do i = -nth, nth
       invR(i) = geom%invR(rgrid(i), theta(i))
    end do
  end subroutine invRtgrid

  !> Calculates the vertical height, Z, on the theta grid
  subroutine Ztgrid(rgrid, theta, Zoftheta)
    implicit none
    real, dimension(-ntgrid:), intent(in) :: rgrid, theta
    real, dimension(-ntgrid:), intent(out) :: Zoftheta
    integer :: i
    do i = -nth, nth
       Zoftheta(i) = geom%Zpos(rgrid(i), theta(i))
    end do
  end subroutine Ztgrid

  !> FIXME : Add documentation
  pure subroutine bvectortgrid(invR, nth, gradrptot, dpsidrp, bi, bvector)
    integer, intent (in) :: nth
    real, dimension(-ntgrid:  ), intent (in) :: invR
    real, dimension(-ntgrid:,:), intent (in) :: gradrptot
    real, dimension(-ntgrid:,:), intent (out) :: bvector
    real, intent (in) :: dpsidrp, bi
    integer :: i
    do i = -nth, nth
       bvector(i,1) = -dpsidrp * gradrptot(i,2) * invR(i)
       bvector(i,2) = dpsidrp * gradrptot(i,1) * invR(i)
       bvector(i,3) = bi * invR(i)
    end do
  end subroutine bvectortgrid

  !> FIXME : Add documentation
  pure subroutine gradstottgrid (invR, grads, gradstot)
    real, dimension(-ntgrid:), intent (in) :: invR
    real, dimension(-ntgrid:, :), intent (in) :: grads
    real, dimension(-ntgrid:, :), intent (out) :: gradstot
    gradstot(:, 1:2) = grads
    if (use_large_aspect) then
       gradstot(:, 3) = 0.
    else
       gradstot(:, 3) = invR
    end if
  end subroutine gradstottgrid

  !> Calculate the cross product c = a x b on the `[-nth, nth]` domain
  pure subroutine crosstgrid(a, b, c)
    real, dimension(-ntgrid:, :), intent (in) ::  a, b
    real, dimension(-ntgrid:, :), intent (out) ::  c
    integer :: i
    do i = -nth, nth
       c(i, 1) = a(i, 2) * b(i, 3) - b(i, 2) * a(i, 3)
       c(i, 2) = a(i, 3) * b(i, 1) - b(i, 3) * a(i, 1)
       c(i, 3) = a(i, 1) * b(i, 2) - b(i, 1) * a(i, 2)
    end do
  end subroutine crosstgrid

  !> Calculate the dot product c = a.b where a, b and c are on the `[-ntgrid, ntgrid]` domain
  pure subroutine dottgridf(a, b, c)
    real, dimension(-ntgrid:,:), intent(in) :: a, b
    real, dimension(-ntgrid:), intent(out) :: c
    c = a(:, 1) * b(:, 1) + a(:, 2) * b(:, 2) + a(:, 3) * b(:, 3)
  end subroutine dottgridf

  !> Returns the diameter of the surface labelled by rp.
  real function diameter(rp)
    implicit none
    real, intent(in) :: rp
    diameter = geom%diameter(rp) !Used to allow rp<rpmin here by just returning zero
  end function diameter

  !> FIXME : Add documentation
  !> Only used with ~undocumented `irho = 4` option. Ideally would update to take pbar
  real function rhofun (rp)
    implicit none
    real, intent (in) :: rp
    real :: pbar
    pbar = psifun(rp)
    rhofun = geom%rhofun(pbar)
  end function rhofun

  !> Calculates the normalised radial coordinate
  pure real function psifun (rp)
    implicit none
    real, intent (in) :: rp
    psifun = min(1.0, max(0.0, (rp - rpmin) / (rpmax - rpmin)))
  end function psifun

  real function pbarfun(r, theta)
    implicit none
    real, intent(in) :: r, theta
    pbarfun = psifun(geom%psi(r, theta))
  end function pbarfun

  !> FIXME : Add documentation
  real function pbarofrho(rho)
    implicit none
    real, intent (in) :: rho
    pbarofrho = psifun(rpofrho(rho))
  end function pbarofrho

  !> FIXME : Add documentation
  subroutine drift(rgrid, rp, gradstot, gradrptot, gradztot, gradbtot, dqdrp, dpsidrho, &
       dpsidrp, bmag, invR, theta, dbetadrho, bi, force_sym, nperiod, &
       gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0)
    implicit none
    real, dimension (-ntgrid:), intent (in) :: rgrid, bmag, invR, theta
    real, dimension (-ntgrid:, :), intent (in) :: gradstot, gradrptot, gradztot, gradbtot
    real, intent (in) :: rp, dqdrp, dpsidrho, dpsidrp, dbetadrho, bi
    logical, intent(in) :: force_sym
    integer, intent(in) :: nperiod
    real, dimension(-ntgrid:ntgrid), intent(out) :: gbdrift, gbdrift0, cvdrift, cvdrift0, cdrift, cdrift0
    real, dimension (-ntgrid:ntgrid, 3) :: pgradtot, dummy, curve, bvector
    real, dimension (-ntgrid:ntgrid) :: gbdrift1, gbdrift2, cvdrift1, cvdrift2, cdrift1, cdrift2
    real :: dp
    if (use_large_aspect) then
       gbdrift = (-sin(theta) * gradstot(:,1) - cos(theta) * gradstot(:,2)) * invR
       gbdrift = 2 * dpsidrho * gbdrift

       gbdrift0 = (-sin(theta) * gradrptot(:,1) - cos(theta) * gradrptot(:,2)) * invR
       gbdrift0 = 2 * dpsidrho * gbdrift0 * dqdrp

       cvdrift = gbdrift ; cvdrift0 = gbdrift0
    else
       !     compute magnetic field
       call bvectortgrid(invR, nth, gradrptot, dpsidrp, bi, bvector)
       dp = dbetadrho / (2 * dpsidrho)

       call crosstgrid(bvector, gradbtot, dummy) !B cross Grad B
       call periodic_copy(dummy(:, 1), 0.0, nperiod)
       call periodic_copy(dummy(:, 2), 0.0, nperiod)
       call periodic_copy(dummy(:, 3), 0.0, nperiod)
       call dottgridf(dummy, gradstot, gbdrift1)
       call dottgridf(dummy, gradrptot, gbdrift2)

       !     create curvature
       if (bishop == 0) then
          call geom%gradient(rgrid, theta, pgradtot(:, 1:2), 'R', rp, nth, ntgrid, bishop /= 0)
          pgradtot(-nth:nth, 3) = 0.0
       else
          pgradtot = gradrptot * dpsidrp * dp !Grad rho * dp/drho
       end if

       curve = gradbtot
       curve(-nth:nth, 1) = curve(-nth:nth, 1) + pgradtot(-nth:nth, 1) / bmag(-nth:nth)
       curve(-nth:nth, 2) = curve(-nth:nth, 2) + pgradtot(-nth:nth, 2) / bmag(-nth:nth)

       call crosstgrid(bvector, curve, dummy) !B cross [Grad B + Grad P / B]
       call periodic_copy(dummy(:, 1), 0.0, nperiod)
       call periodic_copy(dummy(:, 2), 0.0, nperiod)
       call periodic_copy(dummy(:, 3), 0.0, nperiod)
       call dottgridf(dummy, gradstot, cvdrift1)
       call dottgridf(dummy, gradrptot, cvdrift2)

       ! factor of 2 appears below because will later multiply
       ! by wunits, which has a factor of 1/2 built-in
       gbdrift = 2 * dpsidrho * gbdrift1 / bmag**3
       gbdrift0 = 2 * dpsidrho * gbdrift2 * dqdrp / bmag**3
       cvdrift = 2 * dpsidrho * cvdrift1 / bmag**3
       cvdrift0 = 2 * dpsidrho * cvdrift2 * dqdrp / bmag**3
    end if

    ! coriolis drift term is -2*vpa*om_phi/Omega * (bhat x (zhat x bhat)) . grad (g+q<phi>F_0/T)
    ! here we calculate 2*(zhat . grad (alpha+q*theta0))/(B/Bref)

    ! this is the Zhat . grad(alpha) part of Zhat . grad(S)
    ! note that a factor of k_y = (n0/a) (drho/dpsiN) has been factored out
    call dottgridf(gradztot, gradstot, cdrift1)
    ! this is the Zhat . grad(q) part of Zhat . grad(S)
    ! note that a factor of k_y * theta0 has been factored out
    call dottgridf(gradztot, gradrptot, cdrift2)

    ! factor of 4 (instead of 2) below because will later mutliply
    ! by wunits, which has a built-in factor of 1/2
    ! cdrift0 is part of coriolis drift dependent on theta0
    cdrift  = -4 * dpsidrho * cdrift1 / bmag
    cdrift0 = -4 * dpsidrho * cdrift2 * dqdrp / bmag

    if (force_sym) then
       call sym(gbdrift, 0, ntgrid)
       call sym(cvdrift, 0, ntgrid)
       call sym(gbdrift0, 1, ntgrid)
       call sym(cvdrift0, 1, ntgrid)
       ! should maybe add in cdrift and cdrift0 here -- MAB
    end if
  end subroutine drift

  !> Forces symmetry in theta in the passed array. If isign == 0 then
  !> `a(+/- theta)` is replaced by the average of `a(+/-theta)`, otherwise
  !> `a(+/- theta)` is replaced by `+/- (a(theta)-a(-theta))`.
  pure subroutine sym(a, isign, ntgrid)
    implicit none
    integer, intent(in) :: isign, ntgrid
    real, dimension(-ntgrid:), intent (in out) :: a
    integer :: i
    if (isign == 0) then
       do i = 1, ntgrid
          a(i) = 0.5 * (a(i) + a(-i))
          a(-i) = a(i)
       end do
    else
       do i = 1, ntgrid
          a(i) = 0.5 * (a(i) - a(-i))
          a(i) = -a(-i)
       end do
       a(0) = 0.
    end if
  end subroutine sym

  !> FIXME: Add documentation. Only used for bishop == 0 testing
  subroutine seikon(rp, theta, invR, qval, bi, nperiod, seik)
    real, intent(in) :: rp
    real, dimension(-ntgrid:), intent(in) :: theta, invR
    real, intent(in out) :: qval
    real, intent(in) :: bi
    integer, intent(in) :: nperiod
    real, dimension(-ntgrid:), intent(out) :: seik
    real :: dpsidrp
    real, dimension(-ntgrid:ntgrid, 2) :: thgrad, rpgrad
    real, dimension(-ntgrid:ntgrid) :: trip, dsdthet, rgrid
    rgrid = rp !<This is not right for cart types, but ok for those compatible with bishop=0.
    call geom%gradient(rgrid, theta, rpgrad, 'P', rp, nth, ntgrid, bishop /= 0)
    call geom%gradient(rgrid, theta, thgrad, 'T', rp, nth, ntgrid, bishop /= 0)
    call tripprod2dtgrid(rpgrad, thgrad, invR, trip) ; call periodic_copy(trip, 0., nperiod)
    call eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp)
  end subroutine seikon

  !> FIXME : Add documentation
  pure subroutine eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp)
    use constants, only: twopi
    implicit none
    real, dimension(-ntgrid:), intent(in) :: theta, invR, trip
    real, dimension(-ntgrid:), intent(out) :: seik, dsdthet
    real, intent(in out) :: qval
    real, intent(in) :: bi
    real, intent(out) :: dpsidrp

    dsdthet = -bi * invR**2 / trip
    call integrate(dsdthet, theta, seik, ntgrid)

    if (local_eq) then
       dpsidrp = -(seik(nth) - seik(-nth)) / (twopi * qval)
       seik = seik / dpsidrp
       dsdthet = dsdthet / dpsidrp
    else
       dpsidrp = 1.
       qval = -(seik(nth) - seik(-nth)) / (twopi)
    end if
  end subroutine eikonal

  !> FIXME : Add documentation
  real function rpofrho(rho)
    implicit none
    real, intent(in) :: rho
    real, parameter :: xerrbi = 1.e-4, xerrsec = 1.e-8
    interface
      real function f(x_)
        real, intent(in) :: x_
      end function f
    end interface
    procedure(f), pointer :: func
    real :: a, b, fval
    integer :: ier

    ier = 123456 ! Default to indicating error
    a = rpmin ; b = rpmax

    if (irho == 1) then
       func => phi ; fval = rho * rho * func(b)
    else if (irho == 2) then
       func => diameter ; fval = rho * func(b)
    else if (irho == 3) then
       func => psifun ; fval = rho * func(b)
    else if (irho == 4) then
       func => rhofun ; fval = rho * func(b)
    end if

    call root(func, fval, a, b, xerrbi, xerrsec, ier, rpofrho)

    if (debug) then
       write (*,*) "Values in root: "
       write (*,*) 'rho:  ', rho
       write (*,*) 'fval: ', fval
       write (*,*) 'rpofrho: ', rpofrho
       write (*,*) 'rpmin: ', rpmin
       write (*,*) 'rpmax: ', rpmax
    end if
    if (ier > 0) write(*, *) 'error in rpofrho, rho=', rho, rpofrho, ier
  end function rpofrho

  !> Integrates arg(grid) from i = 0 (generally theta = 0) to either
  !> ends of the grid (usually +/- pi) storing the cumulative integral
  !> in the output `ans`.
  pure subroutine integrate(arg, grid, ans, n)
    real, dimension(-ntgrid:), intent (in) :: arg, grid
    real, dimension(-ntgrid:), intent (out) :: ans
    integer, intent(in) :: n
    integer :: i
    ans = 0.
    do i = 1, n
       ans(i) = ans(i-1) + 0.5 * (grid(i) - grid(i-1)) * (arg(i) + arg(i-1))
    end do
    do i = -1, -n, -1
       ans(i) = ans(i+1) + 0.5 * (grid(i) - grid(i+1)) * (arg(i) + arg(i+1))
    end do
  end subroutine integrate

  !> Find soln where f(soln) ~= fval. Works best if a and b bracket the root.
  !> We can then first bisect the bracket to narrow in on the solution.
  !> Following this we use the secant method to find the solution.
  !> The input xerrbi influences the number of bisection iterations whilst
  !> xerrsec is used to set the stopping tolerance on the secant method.
  subroutine root(f, fval, a, b, xerrbi, xerrsec, ier, soln)
    implicit none
    interface
      real function f(x_)
        real, intent(in) :: x_
      end function f
    end interface
    real, intent(in) :: fval, a, b, xerrbi, xerrsec
    real, intent(out) :: soln
    integer, intent(out) :: ier
    real :: a1, b1, f1, f2, f3, trial, aold
    integer :: i, niter
    logical :: skip_bisection
    ier = 0 ; a1 = a ; b1 = b
    f1 = f(a1) - fval ;  f2 = f(b1) - fval
    skip_bisection = xerrbi <= 0 !Could be xerrbi <= abs(b1 - a1)?
    if (.not. skip_bisection .and. f1 * f2 > 0.) then
       write(*, *) 'f1 and f2 have same sign in bisection routine root'
       write(*, *) 'a1,f1,b1,f2=', a1, f1, b1, f2
       write(*, *) 'fval=', fval
       ier = 1
       skip_bisection = .true.
    end if

    if (.not. skip_bisection) then
       ! Presumably we require |b1 - a1| >= xerrbi here to avoid niter < 1?
       niter = 1 + int(log(abs(b1 - a1) / xerrbi) / log(2.))
       ! Repeatedly bisect between limits, swapping one limit
       ! for the bisection point each iteration. Helps to narrow
       ! down the bracketing region substantially.
       do i = 1, niter
          trial = 0.5 * (a1 + b1)
          f3 = f(trial) - fval
          if (f3 * f1 > 0.) then
             a1 = trial
             f1 = f3
          else
             b1 = trial
             f2 = f3
          end if
       end do
    end if

    !Swap so |f1| < |f2|
    if( abs(f1) > abs(f2) ) then
       f3 = f1 ; f1 = f2 ; f2 = f3
       aold = a1 ; a1 = b1 ; b1 = aold
    end if
    !     start secant method
    do i = 1, 10
       aold = a1 ; f3 = f2 - f1
       if (abs(f3) < 1.e-11)  f3 = 1.e-11 !Should we exit here? Saying f2~f1
       a1 = a1 - f1 * (b1 - a1) / f3 !Update guess for root.
       f2 = f1 ; b1 = aold ; f1 = f(a1) - fval
       if (abs(a1 - b1) < xerrsec) exit !Root found closely enough so exit
    end do
    soln = a1
    ! Detect instances where we have not converged to a solution
    if (abs(f1) > xerrsec) ier = ier + 2
  end subroutine root

  !> FIXME : Add documentation
  real function phi(rp)
    implicit none
    real, intent (in) :: rp
    integer, parameter :: nimax = 200
    real :: pbar, pb(nimax), dpb, qi, qnext
    integer :: i

    pbar = psifun(rp)
    ! Construct a uniformly space grid from 0 to pbar
    dpb = pbar / float(nimax - 1)
    do i = 1, nimax
       pb(i) = dpb * float(i - 1)
    end do

    ! Integrate area under q(pb) curve from pb = 0 to pbar
    phi = 0. ; qi = geom%qfun(pb(1))
    do i = 1, nimax - 1
       qnext = geom%qfun(pb(i+1))
       phi = phi + 0.5 * dpb * (qnext + qi)
       qi = qnext
    end do
    phi = phi * (rpmax - rpmin)
  end function phi

  !> Given an initial theta and gradpar (=dtheta/dl) compute a new
  !> theta grid, returned in arcl, which has gradpar = dg/dl  = constant.
  !> @note The input gradpar is modified on output.
  pure subroutine arclength (ntheta, nperiod, theta, gpar, arcl)
    use constants, only: twopi, pi
    implicit none
    integer, intent (in) :: ntheta, nperiod
    real, dimension(-ntgrid:), intent (out) :: arcl
    real, dimension(-ntgrid:), intent (in) :: theta
    real, dimension(-ntgrid:), intent (in out) :: gpar
    integer :: nth, j
    real :: arcfac
    nth = ntheta / 2
    arcl(-nth) = 0.
    do j = -nth, nth-1
       arcl(j+1) = arcl(j) + 0.5 * (theta(j+1) - theta(j)) * (1. / gpar(j) + 1. / gpar(j+1))
    end do
    arcfac = 1. / arcl(nth)

    do j = -nth, nth
       arcl(j) = arcl(j) * twopi * arcfac - pi
    end do

    call periodic_copy(arcl, twopi, nperiod)
    gpar = twopi * arcfac
  end subroutine arclength

  !> FIXME : Add documentation
  subroutine loftheta(rgrid, theta, ltheta)
    implicit none
    real, dimension (-ntgrid:), intent(in) :: rgrid, theta
    real, dimension (-ntgrid:), intent(out) :: ltheta
    integer :: i
    real :: rold, told, rpos_old, zpos_old, rpos_new, zpos_new
    real :: r, thet, rpos_orig, zpos_orig
    ltheta(0) = 0.
    rold = rgrid(0) ; told = theta(0)
    rpos_orig = geom%Rpos(rold, told) ; zpos_orig = geom%Zpos(rold, told)
    rpos_new = rpos_orig ; zpos_new = zpos_orig
    do i = 1, nth
       rpos_old = rpos_new ; zpos_old = zpos_new
       r = rgrid(i) ; thet = theta(i)
       rpos_new = geom%Rpos(r, thet) ; zpos_new = geom%Zpos(r, thet)
       ltheta(i) = ltheta(i-1) + hypot((Rpos_new - Rpos_old), (Zpos_new - Zpos_old))
    end do

    rpos_new = rpos_orig ; zpos_new = zpos_orig
    do i = -1, -nth, -1
       rpos_old = rpos_new ; zpos_old = zpos_new
       r = rgrid(i) ; thet = theta(i)
       rpos_new = geom%Rpos(r, thet) ; zpos_new = geom%Zpos(r, thet)
       ltheta(i) = ltheta(i+1) - hypot((Rpos_new - Rpos_old), (Zpos_new - Zpos_old))
    end do
  end subroutine loftheta

  !> Compute df / dltheta using a periodic central difference. The
  !> input `ext` can be used to control the effective jump in f
  !> between both ends of the "periodic" grid.
  pure subroutine gradl (ltheta, f, dfdl, ext, n)
    implicit none
    real, dimension (-ntgrid:), intent(in) :: ltheta, f
    real, dimension (-ntgrid:), intent(out) :: dfdl
    real, intent(in) :: ext
    integer, intent(in) :: n
    integer :: i
    dfdl(-n) = (f(n-1) + ext - f(-n+1)) / &
         (ltheta(n-1) - ltheta(-n+1) - ltheta(n) + ltheta(-n))
    do i = -n+1, n-1
       dfdl(i) = (f(i+1) - f(i-1)) / (ltheta(i+1) - ltheta(i-1))
    end do
    dfdl(n) = dfdl(-n)
  end subroutine gradl

  !> Finds the angle between the surface tangent vector in the poloidal plane
  !> and the direction of increasing major radius.
  pure subroutine th_bishop(rpgrad, th_bish, nth)
    use constants, only: twopi
    implicit none
    integer, intent (in) :: nth
    real, dimension (-ntgrid:, :), intent (in) :: rpgrad
    real, dimension (-ntgrid:), intent (out) :: th_bish
    real, dimension (-nth:nth) :: magrp
    real, dimension (-nth:nth, 2) :: tvec
    integer :: i
    real :: tvec_dot_gradR

    magrp = hypot(rpgrad(-nth:nth, 1), rpgrad(-nth:nth, 2))

    ! tvec is the tangent vector in the surface of the poloidal plane
    tvec(:, 1) = rpgrad(-nth:nth, 2) / magrp
    tvec(:, 2) = -rpgrad(-nth:nth, 1) / magrp

    ! need to find grad(R) and dot it with tvec.  but grad(R) is
    ! simple -- it is just (1,0) so simply return acos(tvec(:,1)) with
    ! corrections for left-hand quadrants
    do i = -nth, nth
       tvec_dot_gradR = min(1., max(-1., tvec(i, 1)))
       if (tvec(i, 2) > 0.) then
          th_bish(i) = twopi - acos(tvec_dot_gradR)
       else
          th_bish(i) = acos(tvec_dot_gradR)
       end if
    end do
  end subroutine th_bishop

  !> FIXME : Add documentation
  pure subroutine B_pol(invR, rpgrad, dpsidrp, bpol, nth)
    implicit none
    integer, intent (in) :: nth
    real, intent(in) :: dpsidrp
    real, dimension (-ntgrid:), intent (in) :: invR
    real, dimension (-ntgrid:, :), intent (in) :: rpgrad
    real, dimension (-ntgrid:)   , intent (out) :: bpol
    bpol(-nth:nth) = hypot(rpgrad(-nth:nth, 1), rpgrad(-nth:nth, 2)) * invR(-nth:nth) * dpsidrp
  end subroutine B_pol

  !> FIXME : Add documentation
  pure subroutine B_mod(invR, bpol, bmag, bi, nth)
    implicit none
    integer, intent (in) :: nth
    real, dimension (-ntgrid:), intent (in) :: invR, bpol
    real, dimension(-ntgrid:), intent (out) :: bmag
    real, intent(in) :: bi
    bmag(-nth:nth) = hypot((bi * invR(-nth:nth)), bpol(-nth:nth))
  end subroutine B_mod

  !> General calculation of the poloidal magnetic field radius of curvature
  !> (see section 2.3 of "GS2 analytic geometry specification") ! JRB
  !> @Note we actually calculate 1/Rpol as this is what we actually use.
  pure subroutine inv_Rpol(theta, th_bish, ltheta, invRpol, nth)
    use constants, only: twopi, pi
    implicit none
    integer, intent (in) :: nth
    real, dimension(-ntgrid:), intent (in) :: th_bish, theta, ltheta
    real, dimension(-ntgrid:), intent (out) :: invRpol
    real, dimension(-ntgrid:ntgrid) :: dthdl
    integer :: i
    real :: dth_bish, dth
    ! Calculate du/dtheta (where u is defined according to Miller PoP 5 973 (1998))
    ! Note that we must properly account for discontinuities in u. u is th_bish
    do i = -nth + 1, nth - 1
       dth_bish = th_bish(i + 1) - th_bish(i - 1) ; dth = theta(i + 1) - theta(i - 1)
       if ( dth_bish > pi) then
          invRpol(i) = (dth_bish - twopi) / dth
       else if(dth_bish < -pi) then
          invRpol(i) = (dth_bish + twopi) / dth
       else
          invRpol(i) = dth_bish / dth
       end if
    end do

    invRpol(nth) = (th_bish(-nth+1) - th_bish(nth-1)) / (theta(-nth+1) - theta(nth-1) + twopi)
    invRpol(-nth) = invRpol(nth)

    ! Calculate dtheta/dl (where l is the poloidal arc length)
    call gradl(ltheta, theta, dthdl, -twopi, nth)

    ! Scale by dthdl to complete calculation of  1/Rpol
    ! (the inverse radius of curvature of the poloidal magnetic field)
    invRpol(-nth:nth) = invRpol(-nth:nth) * dthdl(-nth:nth)
  end subroutine inv_Rpol

  !> FIXME : Add documentation
  subroutine finish_geometry
    implicit none
    if (allocated(geom)) then
       call geom%finalise
       deallocate(geom)
    end if
  end subroutine finish_geometry

  !> Initialse [[geometry:theta]] to a uniform grid with `2*((2*nperiod - 1)*(nt/2)) + 1`
  !> points. Also set [[geometry:nth]], [[geometry:ntheta]] and [[geometry::ntgrid]]
  !> consistently; importantly [[geometry:nth]] is forced to be even. If `ntheta` is odd, then
  !> `nth == ntheta - 1`.
  !>
  !> `ntheta` is set equal to [[geometry:ntheta]] on exit
  subroutine init_uniform_theta_grid(ntheta, theta, nperiod, ntgrid, nth)
    use constants, only: pi
    implicit none
    integer, intent(in out) :: ntheta
    real, dimension(:), allocatable, intent(in out) :: theta
    integer, intent(in) :: nperiod
    integer, intent(out) :: ntgrid, nth
    integer :: i

    nth = ntheta / 2
    ! Make sure ntheta is even, update argument and set ntgrid
    ntheta = nth * 2 ; ntgrid = (2 * nperiod - 1) * nth

    if (verb>3) write(6,*) "init_theta: allocated(theta),ntgrid,ntheta,nperiod=",&
      allocated(theta), ntgrid, ntheta, nperiod

    if (allocated(theta)) deallocate (theta)
    allocate(theta(-ntgrid:ntgrid))

    theta = [ (i * pi / real(nth), i=-ntgrid, ntgrid) ]
  end subroutine init_uniform_theta_grid

  !> Copy the `[-nth,nth] = [-pi,pi]` section of the passed array into
  !> the `theta > +/- pi` parts, adding on `k * ext` (with k an integer
  !> labelling the sections/periods). Here we expect `ext == array(nth) - array(-nth)`.
  pure subroutine periodic_copy(a, ext, nperiod)
    implicit none
    real, dimension(-ntgrid:), intent(in out) :: a
    real, intent(in) :: ext
    integer, intent(in) :: nperiod
    integer :: i, k, itot, ntheta
    if (nperiod <= 1) return
    ! Note we end up enforcing a(nth)-a(-nth) = ext, but we should probably check this is
    ! also true on input to the routine
    ntheta = 2 * nth
    do k = -nperiod + 1, -1
       do i = -nth, nth
          itot = i + k * ntheta
          a(itot) = a(i) + k * ext
       end do
    end do

    do k = 1, nperiod - 1
       do i = -nth, nth
          itot = i + k * ntheta
          a(itot) = a(i) + k * ext
       end do
    end do
  end subroutine periodic_copy

  !> FIXME : Add documentation
  pure subroutine bishop_gradB(bmag, Bpol, invRpol, th_bish, ltheta, rmajor, dp, bi, gradB)
    implicit none
    real, dimension(-ntgrid:), intent (in) :: bmag, bpol, invRpol
    real, dimension(-ntgrid:), intent (in) :: th_bish, ltheta, rmajor
    real, intent(in) :: dp, bi
    real, dimension(-ntgrid:,:), intent(out) :: gradB
    real :: bp
    integer :: i
    do i = -nth, nth
       bp = -bpol(i)
       gradB(i, 1) = (bp / bmag(i)) * (bp * invRpol(i) + dp * rmajor(i) &
            - bi**2 * sin(th_bish(i)) / (rmajor(i)**3 * bp)) !dBdrho
    end do
    call gradl(ltheta, bmag, gradB(:, 2), 0., nth) !dBdl
    gradB(:, 3) = 0.0
  end subroutine bishop_gradB

  !> Perform minor checks on consistency of flags
  subroutine check(geq, eeq, peq, leq, deq, ideq, ceq, gs2deq, teq, bishop)
    use mp, only: mp_abort
    implicit none
    logical, intent(in) :: geq, eeq, peq, leq, deq, ideq, ceq, gs2deq, teq
    integer, intent(in) :: bishop
    integer :: nflag_set
    if (.not.(geq .or. ideq) .and. bishop == 0) & !Possibly also need use_large_aspect if bishop == 0?
         call mp_abort('bishop = 0 only compatible with geq or ideq',.true.)
    nflag_set = count([geq, ceq, eeq, gs2deq, peq, teq, leq, deq, ideq])
    if (nflag_set == 0) then
       call mp_abort('Must set at least one equilibrium flag to true', .true.)
    else if (nflag_set > 1) then
       call mp_abort('More than one equilibrium flag to true, not allowed', .true.)
    end if
  end subroutine check

  !> Calculate the derivative of the flux label rho with respect to
  !> the internal flux label rp (the normalised poloidal flux)
  real function drho_drp(rp, dr)
    implicit none
    real, intent(in) :: rp, dr
    real :: rp1, rp2, rho1, rho2, phirpmax
    rp1 = rp * (1 - dr) ; rp2 = rp * (1 + dr)
    if (irho == 1) then
       phirpmax = phi(rpmax)
       rho1 = sqrt(abs(phi(rp1) / phirpmax)) ; rho2 = sqrt(abs(phi(rp2) / phirpmax))
    else if (irho == 2) then
       rho1 = 0.5 * diameter(rp1) ; rho2 = 0.5 * diameter(rp2)
    else if (irho == 3) then
       rho1 = psifun(rp1) ; rho2 = psifun(rp2)
    else if (irho == 4) then
       rho1 = rhofun(rp1) ; rho2 = rhofun(rp2)
    end if
    drho_drp = (rho2 - rho1) / (rp2 - rp1)
  end function drho_drp

  !> FIXME : Add documentation | Only used for writelots output and if irho=1
  real function drho_drhod(rp, drhodrp, dr)
    implicit none
    real, intent (in) :: drhodrp, rp, dr
    real :: rp1, rp2,  rhod1, rhod2
    !     compute d rho / d rb
    rp1 = rp * (1 - dr) ; rp2 = rp * (1 + dr)
    rhod1 = 0.5 * diameter(rp1) ; rhod2 = 0.5 * diameter(rp2)
    drho_drhod = drhodrp / ((rhod2 - rhod1) / (rp2 - rp1))
  end function drho_drhod

  !> FIXME : Add documentation
  pure real function fluxavg(f, theta, jacob)
    implicit none
    real, dimension(-nth:nth), intent(in) :: f, theta, jacob
    real, dimension(-nth:nth) :: delth
    delth(-nth+1:nth-1) = 0.5 * (theta(-nth+2:nth) - theta(-nth:nth-2))
    delth(-nth) = 0.5 * (theta(-nth+1) - theta(-nth))
    delth(nth) = 0.5 * (theta(nth) - theta(nth-1))

    ! We might want to write this to use trapezoidal_integration instead
    ! of this sum form. Currently only used in [[f_trap]], which itself
    ! is only used here for one piece of information written to screen.
    ! As it doesn't get into the simulated system we probably don't mind
    ! too much beyond readability concerns.
    fluxavg = sum(f * jacob * delth) / sum(jacob * delth)
  end function fluxavg

  !> FIXME : Add documentation
  pure real function f_trap(b_mag, theta, jacob)
    implicit none
    real, dimension(-nth:nth), intent (in) :: b_mag, theta, jacob
    real :: ftu, ftl, havg, h2avg, h(-nth:nth)
    ! use expressions from Y. R. Lin-Liu and Bob Miller, coding from Q. P. Liu
    h = min(b_mag / maxval(b_mag),1.0)
    havg = fluxavg(h, theta, jacob)
    h2avg = fluxavg(h**2, theta, jacob)
    ftu = 1.0 - (h2avg / havg**2) * (1.0 - sqrt(1.0 - havg) * (1.0 + 0.5 * havg))
    ftl = 1.0 - h2avg * fluxavg((1.0 - sqrt(1.0 - h) * (1.0 + 0.5 * h)) / h**2, theta, jacob)
    f_trap = 0.75 * ftu + 0.25 * ftl
  end function f_trap

  !> FIXME : Add documentation
  subroutine plotdata (rgrid, seik, grads, dpsidrho, dr, theta, &
       rplot, zplot, aplot, rprime, zprime, aprime)
    implicit none
    real, dimension (-ntgrid:), intent(in) :: rgrid, seik, theta
    real, dimension (-ntgrid:, :), intent(in) :: grads
    real, dimension (-ntgrid:), intent(out) :: rplot, zplot, aplot, rprime, zprime, aprime
    real, intent (in) :: dpsidrho, dr
    real :: rplus, rminus, dr_tot
    integer :: i
    ! At the present flux surface, calculate these three functions
    ! as a function of theta, together with their derivates with
    ! respect to rho:
    !
    ! R, Z, alpha-phi
    !
    ! Call these variables
    !
    ! Rplot, Zplot, aplot
    !
    ! and their derivatives:
    !
    ! Rprime, Zprime, aprime
    !
    ! BD: bug??
    ! Is this correct for EFIT data?  Need to check
    do i = -ntgrid, ntgrid
       Rplot(i) = geom%Rpos(rgrid(i), theta(i))
       Zplot(i) = geom%Zpos(rgrid(i), theta(i))
    end do
    aplot = seik

    do i = -ntgrid, ntgrid
       rplus  = rgrid(i) * (1 + dr) ; rminus = rgrid(i) * (1 - dr)
       dr_tot = 2 * dr * rgrid(i)
       Rprime(i) = (geom%Rpos(rplus, theta(i)) - geom%Rpos(rminus, theta(i))) / dr_tot
       Zprime(i) = (geom%Zpos(rplus, theta(i)) - geom%Zpos(rminus, theta(i))) / dr_tot
    end do
    aprime = grads(:, 1) * dpsidrho
  end subroutine plotdata
end module geometry