theta_grid_params.f90 Source File


Contents

Source Code


Source Code

!> FIXME : Add documentation
module theta_grid_params
  use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN
  use geo_utils, only: mxh_max_moments

  implicit none

  private

  public :: init_theta_grid_params, finish_theta_grid_params, wnml_theta_grid_params
  public :: rhoc, rmaj, r_geo, eps, epsl, qinp, shat, alpmhd
  public :: pk, geoType, aSurf, shift, shiftVert, mMode, nMode, deltam, deltan
  public :: deltampri, deltanpri, thetam, thetan
  public :: btor_slab, betaprim, ntheta, nperiod
  public :: set_overrides, get_parameters_as_surf

  public :: theta_grid_parameters_config_type
  public :: set_theta_grid_parameters_config
  public :: get_theta_grid_parameters_config

  real :: rhoc, rmaj, r_geo, eps, epsl
  real :: qinp, shat, alpmhd, pk
  integer :: geoType, mMode, nMode
  real :: aSurf, shift, shiftVert, raxis, zaxis, deltam, deltan
  real :: akappa, tri, delta2, delta3, deltampri, deltanpri
  real :: akappri, tripri, thetam, thetan, thetak, thetad, theta2, theta3
  real :: btor_slab, betaprim
  integer :: ntheta, nperiod, n_mxh
  real :: c0_mxh
  real, dimension(mxh_max_moments) :: c_mxh, s_mxh, dc_mxh_dr, ds_mxh_dr

  logical :: initialized = .false.
  real :: kp = -1.
  logical :: exist

  !> Used to represent the input configuration of theta_grid. Sets a
  !> number of parameters used by the different theta grid
  !> implementations. Not all parameters are active for a given theta
  !> grid type.
  type, extends(abstract_config_type) :: theta_grid_parameters_config_type
     ! namelist : theta_grid_parameters
     ! indexed : false
     !> The flux surface elongation (only used when `geoType=0`, which
     !> selects the Generalised Miller analytic geometry
     !> specification).
     real :: akappa = 1.0
     !> The radial gradient flux surface elongation (only used when
     !> `geoType=0`, which selects the Generalised Miller analytic
     !> geometry specification).  `akappri` = \(d\kappa/d\rho\)
     real :: akappri = 0.0
     !> Used in conjunction with [[theta_grid_salpha_knobs:alpmhdfac]]
     !> to override `shift`, set as `shift=-alpmhd*alpmhdfac`. Do not
     !> use unless you know what you are doing.
     real :: alpmhd = 0.0
     !> Minor radius of the flux surface that receives the specified
     !> shaping (only used when `geoType=1`, which selects the Global
     !> MHD analytic geometry specification). See Section 2.1.2 of the
     !> [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: asurf = 1.0
     !> In the slab limit, determines the angle between the field and
     !> the background flow (which flows in an imaginary toroidal
     !> direction). It is effectively equal to
     !> \(\frac{B_t}{B}=\frac{u_{\parallel}}{u}\).  `btor_slab =
     !> btor/bpol` defines direction of a flow relative to B in slab
     !> geometry, where flow is by definition in the toroidal
     !> direction.
     real :: btor_slab = 0.0
     !> Zeroth cosine moment for the MXH local equilibrum. See [[geoType]] and [[local_eq]]
     real :: c0_mxh = 0.0
     !> Cosine moments for the MXH local equilibrum. See [[geoType]] and [[local_eq]]
     real, dimension(mxh_max_moments) :: c_mxh = 0.0
     !> Radial derivatives of cosine moments for the MXH local equilibrum. See [[geoType]] and [[local_eq]]
     real, dimension(mxh_max_moments) :: dc_mxh_dr = 0.0
     !> The elongation of the flux surface labeled by aSurf (only used
     !> when `geoType`=1, which selects the Global MHD analytic
     !> geometry specifications). See Section 2.1.2 of the [Analytic
     !> Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: delta2 = 1.0
     !> The triangularity of the flux surface labeled by aSurf (only
     !> used when geoType=1, which selects the Global MHD analytic
     !> geometry specifications). See Section 2.1.2 of the [Analytic
     !> Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: delta3 = 1.0
     !> The magnitude of the mMode shaping effect (only used when
     !> `geoType`=2 or 3, which selects the Generalised Elongation or
     !> Fourier analytic geometry specifications). See Sections 2.1.3
     !> and 2.1.4 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: deltam = 1.0
     !> Radial derivative of the magnitude of the mMode shaping effect
     !> (only used when `geoType`=2 or 3, which selects the
     !> Generalised Elongation or Fourier analytic geometry
     !> specifications). See Sections 2.1.3 and 2.1.4 of the [Analytic
     !> Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: deltampri = 0.0
     !> The magnitude of the nMode shaping effect (only used when
     !> `geoType`=2 or 3, which selects the Generalised Elongation or
     !> Fourier analytic geometry specifications). See Sections 2.1.3
     !> and 2.1.4 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: deltan = 1.0
     !> Radial derivative of the magnitude of the nMode shaping effect
     !> (only used when `geoType`=2 or 3, which selects the
     !> Generalised Elongation or Fourier analytic geometry
     !> specifications). See Sections 2.1.3 and 2.1.4 of the [Analytic
     !> Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: deltanpri = 0.0
     !> Controls particle trapping (among other things) in simple
     !> geometric models. \(\epsilon = r/R\)
     real :: eps = 0.3
     !>  \(\epsilon_\ell = \frac{2 L_\textrm{ref}}{ R}\) Sets
     !>  curvature drift in s-alpha model, where \(L_\textrm{ref}\) is
     !>  the GS2 equilibrium reference normalisation length and \(R\)
     !>  is the major radius at the centre of the flux surface.
     real :: epsl = 0.3
     !> Selects between different analytic geometry specifications
     !> (only used when `local_eq = T`):
     !>
     !> - `geoType = 0`: Generalised Miller
     !> - `geoType = 1`: Global MHD
     !> - `geoType = 2`: Generalised Elongation
     !> - `geoType = 3`: Fourier specification.
     !> - `geoType = 4`: Miller Extended Harmonic (MXH)
     !>
     !> See Section 2.1 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details on options 0-3, and [PPCF 63
     !> (2021) 012001 (5pp)](https://doi.org/10.1088/1361-6587/abc63b)
     !> for details of 4 (MXH).
     !>
     !> @note The default for this *MUST* be zero otherwise the
     !> Trinity interface will break.
     integer :: geotype = 0
     !> `kp` sets parallel wavenumber and box length in slab
     !> geometry. \(k_p = \frac{2 \pi L_\textrm{ref}}{L_z}\).
     !>
     !> - always use `kp` instead of `pk` in slab geometry (if `kp > 0`
     !> then gs2 sets `pk = 2*kp`)
     !> - in slab limit \(\textrm{shat} = \frac{L_z}{2 \pi L_s} =
     !> \frac{L_\textrm{ref}}{k_p L_s}\) : nb if `kp = 1`, \(\textrm{shat} =
     !> \frac{L_\textrm{ref}}{L_s}\), where \(L_s\) is the magnetic shear
     !> scale length.
     !>
     real :: kp = -1.0
     !> First flux surface shaping mode number (only used when
     !> `geoType`=2 or 3, which selects the Generalised Elongation or
     !> Fourier analytic geometry specifications). See Sections 2.1.3
     !> and 2.1.4 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     integer :: mmode = 2
     !> Second flux surface shaping mode number (only used when
     !> `geoType`=2 or 3, which selects the Generalised Elongation or
     !> Fourier analytic geometry specifications). See Sections 2.1.3
     !> and 2.1.4 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     integer :: nmode = 3
     !> Number of moments for the MXH local equilibrum. Maximum value
     !> is [[mxh_max_moments]]. See [[geoType]] and [[local_eq]] for
     !> more details.
     integer :: n_mxh = 0
     !> Sets the number of \(2\pi\) segments along equilibrium
     !> magnetic field to include in simulation
     !> domain. \(N_\textrm{segments} = (2n_\textrm{period} - 1)\).
     !> Ignored in some cases
     !>
     !> @todo Document when this is ignored
     integer :: nperiod = 2
     !> Rough number of grid points along equilibrium magnetic field between \(\theta=[-\pi,\pi]\).
     !> Actual number of grid points determined as follows:
     !>
     !> - number of points in GS2 theta grid always odd because grid
     !>   must contain both bounce points of trapped particles and one
     !>   grid point at \(\theta=0\). For even values, an extra point
     !>   is added
     !> - theta grid in code runs from `-ntgrid:ntgrid`, with
     !>   `ntgrid=int(ntheta/2)` for `nperiod=1`.
     !> - `ntheta` is used for local equilibria, s-alpha, and EFIT,
     !>   and ignored for all other numerical equilibria
     integer :: ntheta = 24
     !> \(p_k = \frac{\textrm{epsl}}{q} = \frac{2 L_\textrm{ref}}{ q R}\)
     !> Sets \(q\), the magnetic safety factor, and therefore also
     !> sets the connection length, i.e. the length of the box in the
     !> parallel direction, in terms of \(L_\textrm{ref}\). Used only in high
     !> aspect ratio \(\hat{s}-\alpha\) equilibrium model.
     real :: pk = 0.3
     !> Sets value of the safety factor when using local toroidal
     !> equilibrium model.
     real :: qinp = 1.5
     !> When not local_eq: Centerpoint of LCFS (normalized to
     !> \(L_\textrm{ref}\)) - When local_eq: Major radius of
     !> magnetic field reference point (normalized to
     !> \(L_\textrm{ref}\)). Specifically, the reference magnetic
     !> field is defined to be the value of the toroidal magnetic
     !> field at `R=r_geo` on the flux surface labeled by `rhoc`.
     real :: r_geo = 3.0
     !> Major radial location of the magnetic axis (only used when
     !> `geoType`=1, which selects the Global MHD analytic geometry
     !> specification). See Section 2.1.2 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: raxis = 3.0
     !> `rhoc` is the flux surface label (0< `rhoc`< 1). Its exact
     !> meaning depends on `irho`. Usually rho = midplane
     !> diameter/midplane diameter of LCFS (here, "LCFS" refers to the
     !> flux surface at a radius of \(L_\textrm{ref}\). If using
     !> Miller geometry with \(L_\textrm{ref} \neq a\), this is
     !> **not** the physical last closed flux surface.
     !>
     !> - When `irho` = 1, `rhoc` = sqrt(toroidal flux)/sqrt(toroidal flux of LCFS)
     !> - When `irho` = 2, `rhoc` = midplane diameter/(midplane diameter of LCFS). recommended
     !> - When `irho` = 3, `rhoc` = poloidal flux/(poloidal flux of LCFS)
     !>
     real :: rhoc = 0.5
     !> When not local_eq: Position of magnetic axis (normalized to \(L_\textrm{ref}\)).
     !> When local_eq: Major radius of the centre of the
     !> flux surface of interest (normalized to \(L_\textrm{ref}\))
     real :: rmaj = 3.0
     !> Sine moments for the MXH local equilibrum. See [[geoType]] and [[local_eq]]
     real, dimension(mxh_max_moments) :: s_mxh = 0.0
     !> Radial derivatives of sine moments for the MXH local equilibrum. See [[geoType]] and [[local_eq]]
     real, dimension(mxh_max_moments) :: ds_mxh_dr = 0.0
     !> Sets value of magnetic shear in simple geometric models.
     !> Over-ridden by `s_hat_input` in [[theta_grid_eik_knobs]] for most values of `bishop`.
     real :: shat = 0.75
     !> shift is related to minor radial derivatives of the major
     !>radial location of the flux surface centers (i.e. the Shafranov
     !>shift), but this input variable has **different physical
     !>definitions** in s-alpha and other analytic equilbrium models:
     !>
     !> - In s-alpha (i.e. equilibrium_option='s-alpha'), shift
     !> \(\propto p^\prime\) is a parameter for local \(J_{\phi}\)
     !> (and not \(B_{\theta}\) which is constant): \(\textrm{shift} =
     !> \frac{2\textrm{epsl}}{\textrm{pk}^2}\frac{d\beta}{d\rho} =
     !> -\frac{q^2R}{L_\textrm{ref}}\frac{d\beta}{d\rho} > 0\)
     !>- In other analytic specifications
     !> (i.e. equilibrium_option='eik' and `geoType`=0, 2, or 3),
     !> shift is a parameter for local \(B_{\theta}\) (and NOT for
     !> \(J_{\phi}\)): \(\textrm{shift} = \frac{1}{L_\textrm{ref}} \frac{dR}{d\rho} <
     !> 0\)
     !>
     !> NB in Miller shift contains the *1st* radial derivative of the
     !> Shafranov shift, BUT in s-alpha shift is related to a *2nd*
     !> radial derivative of the Shafranov shift.
     !>
     !> - in s-alpha(i.e. equilibrium_option='s-alpha'), no additional parameter
     !> is required as the piece of \(J_{\phi} \propto\) p' is
     !> specified by shift.
     !> - in other analytic specifications (i.e. equilibrium_option='eik'),
     !> an additional parameter(beta_prime) is required to specify the piece of
     !> \(J_{\phi} \propto\) p'
     !>
     real :: shift = 0.0
     !> Minor radial derivative of the axial location of the flux
     !> surface centers (i.e. the vertical Shafranov shift). It is
     !> only used when equilibrium_option='eik' and `geoType`=0, 2, or
     !> 3 (which selects the Generalised Miller, Generalised
     !> Elongation, or Fourier analytic geometry specifications). See
     !> Sections 2.1.1, 2.1.3, and 2.1.4 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: shiftvert = 0.0
     !> The tilt angle of the elongation of the flux surface labeled
     !> by `aSurf` (only used when `geoType`=1, which selects the
     !> Global MHD analytic geometry specification). See Section 2.1.2
     !> of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: theta2 = 0.0
     !> The tilt angle of the triangularity of the flux surface
     !> labeled by `aSurf` (only used when `geoType`=1, which selects
     !> the Global MHD analytic geometry specification). See Section
     !> 2.1.2 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: theta3 = 0.0
     !> The tilt angle of the triangularity (only used when
     !> `geoType`=0, which selects the Generalised Miller
     !> specification). See Section 2.1.1 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: thetad = 0.0
     !> The tilt angle of the elongation (only used when `geoType`=0,
     !> which selects the Generalised Miller specification). See
     !> Section 2.1.1 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: thetak = 0.0
     !> The tilt angle of the mMode shaping effect (only used when
     !> `geoType`=2 or 3, which selects the Generalised Elongation or
     !> Fourier analytic geometry specifications). See Sections 2.1.3
     !> and 2.1.4 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: thetam = 0.0
     !> The tilt angle of the nMode shaping effect (only used when
     !> `geoType`=2 or 3, which selects the Generalised Elongation or
     !> Fourier analytic geometry specifications). See Sections 2.1.3
     !> and 2.1.4 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: thetan = 0.0
     !>The flux surface triangularity (only used when `geoType`=0,
     !>which selects the Generalised Miller analytic geometry
     !>specification).  triangularity is `tri =
     !>arcsin[(R(max(Z))-R_major)/r_mid]`
     real :: tri = 0.0
     !> The radial gradient of the flux surface triangularity (only
     !> used when `geoType`=0, which selects the Generalised Miller
     !> analytic geometry specification). `tripri = dtri/drho`
     real :: tripri = 0.0
     !> Axial location of the magnetic axis (only used when geoType=1,
     !> which selects the Global MHD analytic geometry
     !> specification). See Section 2.1.2 of the [Analytic Geometry
     !> Specification](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
     !> documentation for more details.
     real :: zaxis = 0.0
   contains
     procedure, public :: read => read_theta_grid_parameters_config
     procedure, public :: write => write_theta_grid_parameters_config
     procedure, public :: reset => reset_theta_grid_parameters_config
     procedure, public :: broadcast => broadcast_theta_grid_parameters_config
     procedure, public, nopass :: get_default_name => get_default_name_theta_grid_parameters_config
     procedure, public, nopass :: get_default_requires_index => get_default_requires_index_theta_grid_parameters_config
  end type theta_grid_parameters_config_type

  type(theta_grid_parameters_config_type) :: theta_grid_parameters_config

contains

  !> FIXME : Add documentation
  subroutine init_theta_grid_params(theta_grid_parameters_config_in)
    use unit_tests, only: debug_message
    implicit none
    type(theta_grid_parameters_config_type), intent(in), optional :: theta_grid_parameters_config_in
    integer, parameter :: verb=3

    call debug_message(verb, "theta_grid_params::init_theta_grid_params start")
    if (initialized) return
    initialized = .true.

    call debug_message(verb, &
      "theta_grid_params::init_theta_grid_params call read_parameters")
    call read_parameters(theta_grid_parameters_config_in)

    betaprim = 0.0

    call debug_message(verb, "theta_grid_params::init_theta_grid_params end")
  end subroutine init_theta_grid_params

  !> FIXME : Add documentation
  subroutine finish_theta_grid_params
    implicit none
    initialized = .false.
    call theta_grid_parameters_config%reset()
  end subroutine finish_theta_grid_params

  !> FIXME : Add documentation
  subroutine read_parameters(theta_grid_parameters_config_in)
    use file_utils, only: input_unit, input_unit_exist
    use unit_tests, only: debug_message
    use geo_utils, only: geo_type_miller, geo_type_global, &
         geo_type_generalized_elongation, geo_type_fourier_series, &
         geo_type_miller_extended_harmonic, mxh_max_moments
    implicit none
    type(theta_grid_parameters_config_type), intent(in), optional :: theta_grid_parameters_config_in
    integer, parameter :: verb=3

    call debug_message(verb, "theta_grid_params::read_parameters start")

    if (present(theta_grid_parameters_config_in)) theta_grid_parameters_config = theta_grid_parameters_config_in

    call theta_grid_parameters_config%init(name = 'theta_grid_parameters', requires_index = .false.)

    ! Copy out internal values into module level parameters
    akappa = theta_grid_parameters_config%akappa
    akappri = theta_grid_parameters_config%akappri
    alpmhd = theta_grid_parameters_config%alpmhd
    asurf = theta_grid_parameters_config%asurf
    btor_slab = theta_grid_parameters_config%btor_slab
    c0_mxh = theta_grid_parameters_config%c0_mxh
    c_mxh = theta_grid_parameters_config%c_mxh
    dc_mxh_dr = theta_grid_parameters_config%dc_mxh_dr
    delta2 = theta_grid_parameters_config%delta2
    delta3 = theta_grid_parameters_config%delta3
    deltam = theta_grid_parameters_config%deltam
    deltampri = theta_grid_parameters_config%deltampri
    deltan = theta_grid_parameters_config%deltan
    deltanpri = theta_grid_parameters_config%deltanpri
    ds_mxh_dr = theta_grid_parameters_config%ds_mxh_dr
    eps = theta_grid_parameters_config%eps
    epsl = theta_grid_parameters_config%epsl
    geotype = theta_grid_parameters_config%geotype
    kp = theta_grid_parameters_config%kp
    mmode = theta_grid_parameters_config%mmode
    n_mxh = theta_grid_parameters_config%n_mxh
    nmode = theta_grid_parameters_config%nmode
    nperiod = theta_grid_parameters_config%nperiod
    ntheta = theta_grid_parameters_config%ntheta
    pk = theta_grid_parameters_config%pk
    qinp = theta_grid_parameters_config%qinp
    r_geo = theta_grid_parameters_config%r_geo
    raxis = theta_grid_parameters_config%raxis
    rhoc = theta_grid_parameters_config%rhoc
    rmaj = theta_grid_parameters_config%rmaj
    s_mxh = theta_grid_parameters_config%s_mxh
    shat = theta_grid_parameters_config%shat
    shift = theta_grid_parameters_config%shift
    shiftvert = theta_grid_parameters_config%shiftvert
    theta2 = theta_grid_parameters_config%theta2
    theta3 = theta_grid_parameters_config%theta3
    thetad = theta_grid_parameters_config%thetad
    thetak = theta_grid_parameters_config%thetak
    thetam = theta_grid_parameters_config%thetam
    thetan = theta_grid_parameters_config%thetan
    tri = theta_grid_parameters_config%tri
    tripri = theta_grid_parameters_config%tripri
    zaxis = theta_grid_parameters_config%zaxis

    exist = theta_grid_parameters_config%exist

    if (kp > 0.) pk = 2.*kp

    ! Print warning if the user has specified non-default values for the
    ! geometrical input parameters of other geometry types ! JRB
    select case (geoType)
       case (geo_type_miller)
          if (mMode/=2) write (*,*) "WARNING: ignoring value of mMode, value not needed"
          if (nMode/=3) write (*,*) "WARNING: ignoring value of nMode, value not needed"
          if (deltam/=1.0) write (*,*) "WARNING: ignoring value of deltam, did you mean akappa?"
          if (deltan/=1.0) write (*,*) "WARNING: ignoring value of deltan, did you mean tri?"
          if (delta2/=1.0) write (*,*) "WARNING: ignoring value of delta2, did you mean akappa?"
          if (delta3/=1.0) write (*,*) "WARNING: ignoring value of delta3, did you mean tri?"
          if (deltampri/=0.0) write (*,*) "WARNING: ignoring value of deltampri, did you mean akappri?"
          if (deltanpri/=0.0) write (*,*) "WARNING: ignoring value of deltanpri, did you mean tripri?"
          if (thetam/=0.0) write (*,*) "WARNING: ignoring value of thetam, did you mean thetak?"
          if (thetan/=0.0) write (*,*) "WARNING: ignoring value of thetan, did you mean thetad?"
          if (theta2/=0.0) write (*,*) "WARNING: ignoring value of theta2, did you mean thetak?"
          if (theta3/=0.0) write (*,*) "WARNING: ignoring value of theta3, did you mean thetad?"
          if (aSurf/=1.0) write (*,*) "WARNING: ignoring value of aSurf, value not needed"
          if (raxis/=3.0) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?"
          if (zaxis/=0.0) write (*,*) "WARNING: ignoring value of zaxis, did you mean shiftVert?"
          deltam = akappa
          deltan = tri
          deltampri = akappri
          deltanpri = tripri
          thetam = thetak
          thetan = thetad
       case (geo_type_global)
          if (mMode/=2) write (*,*) "WARNING: ignoring value of mMode, value not needed"
          if (nMode/=3) write (*,*) "WARNING: ignoring value of nMode, value not needed"
          if (deltam/=1.0) write (*,*) "WARNING: ignoring value of deltam, did you mean delta2?"
          if (deltan/=1.0) write (*,*) "WARNING: ignoring value of deltan, did you mean delta3?"
          if (akappa/=1.0) write (*,*) "WARNING: ignoring value of akappa, did you mean delta2?"
          if (tri/=0.0) write (*,*) "WARNING: ignoring value of tri, did you mean delta3?"
          if (deltampri/=0.0) write (*,*) "WARNING: ignoring value of deltampri, value not needed"
          if (deltanpri/=0.0) write (*,*) "WARNING: ignoring value of deltanpri, value not needed"
          if (akappri/=0.0) write (*,*) "WARNING: ignoring value of akappri, value not needed"
          if (tripri/=0.0) write (*,*) "WARNING: ignoring value of tripri, value not needed"
          if (thetam/=0.0) write (*,*) "WARNING: ignoring value of thetam, did you mean theta2?"
          if (thetan/=0.0) write (*,*) "WARNING: ignoring value of thetan, did you mean theta3?"
          if (thetak/=0.0) write (*,*) "WARNING: ignoring value of thetak, did you mean theta2?"
          if (thetad/=0.0) write (*,*) "WARNING: ignoring value of thetad, did you mean theta3?"
          if (shift/=0.0) write (*,*) "WARNING: ignoring value of shift, did you mean Raxis?"
          if (shiftVert/=0.0) write (*,*) "WARNING: ignoring value of shiftVert, did you mean Zaxis?"

          if (raxis>=Rmaj+aSurf) write (*,*) "WARNING: Raxis>=Rmag+aSurf, may specify non-nested flux surfaces"
          if (raxis<=Rmaj-aSurf) write (*,*) "WARNING: Raxis<=Rmag-aSurf, may specify non-nested flux surfaces"
          if (zaxis>=aSurf) write (*,*) "WARNING: Zaxis>=aSurf, may specify non-nested flux surfaces"
          if (zaxis<=-aSurf) write (*,*) "WARNING: Zaxis<=-aSurf, may specify non-nested flux surfaces"

          shift = raxis
          shiftVert = zaxis
          deltam = delta2
          deltan = delta3
          thetam = theta2
          thetan = theta3
       case (geo_type_generalized_elongation)
          if (delta2/=1.0) write (*,*) "WARNING: ignoring value of delta2, did you mean deltam?"
          if (delta3/=1.0) write (*,*) "WARNING: ignoring value of delta3, did you mean deltan?"
          if (akappa/=1.0) write (*,*) "WARNING: ignoring value of akappa, did you mean deltam?"
          if (tri/=0.0) write (*,*) "WARNING: ignoring value of tri, did you mean deltan?"
          if (akappri/=0.0) write (*,*) "WARNING: ignoring value of akappri, did you mean deltampri?"
          if (tripri/=0.0) write (*,*) "WARNING: ignoring value of tripri, did you mean deltanpri?"
          if (theta2/=0.0) write (*,*) "WARNING: ignoring value of theta2, did you mean thetam?"
          if (theta3/=0.0) write (*,*) "WARNING: ignoring value of theta3, did you mean thetan?"
          if (thetak/=0.0) write (*,*) "WARNING: ignoring value of thetak, did you mean thetam?"
          if (thetad/=0.0) write (*,*) "WARNING: ignoring value of thetad, did you mean thetan?"
          if (aSurf/=1.0) write (*,*) "WARNING: ignoring value of aSurf, value not needed"
          if (raxis/=3.0) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?"
          if (zaxis/=0.0) write (*,*) "WARNING: ignoring value of zaxis, did you mean shiftVert?"

       case (geo_type_fourier_series)
          if (delta2/=1.0) write (*,*) "WARNING: ignoring value of delta2, did you mean deltam?"
          if (delta3/=1.0) write (*,*) "WARNING: ignoring value of delta3, did you mean deltan?"
          if (akappa/=1.0) write (*,*) "WARNING: ignoring value of akappa, did you mean deltam?"
          if (tri/=0.0) write (*,*) "WARNING: ignoring value of tri, did you mean deltan?"
          if (akappri/=0.0) write (*,*) "WARNING: ignoring value of akappri, did you mean deltampri?"
          if (tripri/=0.0) write (*,*) "WARNING: ignoring value of tripri, did you mean deltanpri?"
          if (theta2/=0.0) write (*,*) "WARNING: ignoring value of theta2, did you mean thetam?"
          if (theta3/=0.0) write (*,*) "WARNING: ignoring value of theta3, did you mean thetan?"
          if (thetak/=0.0) write (*,*) "WARNING: ignoring value of thetak, did you mean thetam?"
          if (thetad/=0.0) write (*,*) "WARNING: ignoring value of thetad, did you mean thetan?"
          if (aSurf/=1.0) write (*,*) "WARNING: ignoring value of aSurf, value not needed"
          if (raxis/=3.0) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?"
          if (zaxis/=0.0) write (*,*) "WARNING: ignoring value of zaxis, did you mean shiftVert?"

       case (geo_type_miller_extended_harmonic)
          if (deltam/=1.0) write (*,*) "WARNING: ignoring value of deltam, did you mean akappa?"
          if (deltampri/=0.0) write (*,*) "WARNING: ignoring value of deltampri, did you mean akappri?"

          if (mMode/=2) write (*,*) "WARNING: ignoring value of mMode, value not needed"
          if (nMode/=3) write (*,*) "WARNING: ignoring value of nMode, value not needed"
          if (deltan/=1.0) write (*,*) "WARNING: ignoring value of deltan, value not needed"
          if (delta2/=1.0) write (*,*) "WARNING: ignoring value of delta2, value not needed"
          if (delta3/=1.0) write (*,*) "WARNING: ignoring value of delta3, value not needed"
          if (deltampri/=0.0) write (*,*) "WARNING: ignoring value of deltampri, value not needed"
          if (deltanpri/=0.0) write (*,*) "WARNING: ignoring value of deltanpri, value not needed"
          if (thetam/=0.0) write (*,*) "WARNING: ignoring value of thetam, value not needed"
          if (thetan/=0.0) write (*,*) "WARNING: ignoring value of thetan, value not needed"
          if (theta2/=0.0) write (*,*) "WARNING: ignoring value of theta2, value not needed"
          if (theta3/=0.0) write (*,*) "WARNING: ignoring value of theta3, value not needed"
          if (aSurf/=1.0) write (*,*) "WARNING: ignoring value of aSurf, value not needed"
          if (raxis/=3.0) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?"
          if (zaxis/=0.0) write (*,*) "WARNING: ignoring value of zaxis, did you mean shiftVert?"

          deltam = akappa
          deltampri = akappri

          if (theta_grid_parameters_config%n_mxh > mxh_max_moments) then
             write(*, '(A, I0, A, I0, A, I0)') &
                  "WARNING: ''n_mxh'' (", &
                  theta_grid_parameters_config%n_mxh, &
                  ") is larger than maximum number of moments for MXH geometry (", &
                  mxh_max_moments, "), truncating to ", mxh_max_moments
          end if

       case default ! We should abort here, but we don't know if local_eq is true so maybe this is fine.
          write (*,*) "ERROR: invalid analytic geometry specification"
    end select

  end subroutine read_parameters

  !> Packs all relevant parameters into a flux_surface_type instance.
  !> Note we can't fully populate the surface instance as we don't hold delrho
  function get_parameters_as_surf() result(surf)
    use geo_utils, only: flux_surface_type
    implicit none
    type(flux_surface_type) :: surf
    surf = flux_surface_type(geotype, &
         mMode=mMode, nMode=nMode, &
         rmaj=rmaj, r_geo=R_geo, r=rhoc, &
         aSurf=aSurf, sHorz=shift, &
         sVert=shiftVert, delm=deltam,&
         deln=deltan, delmp=deltampri, &
         delnp=deltanpri, thm=thetam, &
         thn=thetan, q=qinp, shat=shat, &
         nt=ntheta,& ! Might not want to pack this here
         n_mxh=n_mxh, c0_mxh=c0_mxh,&
         c_mxh=c_mxh, s_mxh=s_mxh,&
         dc_mxh_dr=dc_mxh_dr, ds_mxh_dr=ds_mxh_dr)
  end function get_parameters_as_surf

  !> FIXME : Add documentation
  subroutine wnml_theta_grid_params(unit)
    implicit none
    integer, intent(in) :: unit
    if (.not.exist) return
    write (unit, *)
    write (unit, fmt="(' &',a)") "theta_grid_parameters"
    write (unit, fmt="(' ntheta =   ',i4)") ntheta
    write (unit, fmt="(' nperiod =  ',i4)") nperiod
    write (unit, fmt="(' rhoc =     ',f7.4)") rhoc
    write (unit, fmt="(' Rmaj =     ',f7.4)") rmaj
    write (unit, fmt="(' R_geo =    ',f7.4)") r_geo
    write (unit, fmt="(' eps =      ',f7.4)") eps
    write (unit, fmt="(' epsl =     ',f7.4)") epsl
    write (unit, fmt="(' qinp =     ',f7.4)") qinp
    write (unit, fmt="(' shat =     ',f7.4)") shat
    write (unit, fmt="(' alpmhd =   ',f7.4)") alpmhd
    write (unit, fmt="(' pk =       ',f7.4)") pk
    write (unit, fmt="(' kp =       ',f7.4)") kp
    write (unit, fmt="(' geoType =  ',i4)") geoType
    write (unit, fmt="(' aSurf =    ',f7.4)") aSurf
    write (unit, fmt="(' shift =    ',f7.4)") shift
    write (unit, fmt="(' shiftVert =',f7.4)") shiftVert
    write (unit, fmt="(' mMode =    ',i4)") mMode
    write (unit, fmt="(' nMode =    ',i4)") nMode
    write (unit, fmt="(' deltam =   ',f7.4)") deltam
    write (unit, fmt="(' deltan =   ',f7.4)") deltan
    write (unit, fmt="(' deltampri =',f7.4)") deltampri
    write (unit, fmt="(' deltanpri =',f7.4)") deltanpri
    write (unit, fmt="(' thetam =   ',f7.4)") thetam
    write (unit, fmt="(' thetan =   ',f7.4)") thetan
    write (unit, fmt="(' btor_slab =',f7.4)") btor_slab
    write (unit, fmt="(' /')")
  end subroutine wnml_theta_grid_params

  !> FIXME : Add documentation
  subroutine set_overrides(mgeo_ov)
    use overrides, only: miller_geometry_overrides_type
    type(miller_geometry_overrides_type), intent(in) :: mgeo_ov
          !write (*,*) 'Calling tgpso'
    if (mgeo_ov%override_rhoc) rhoc = mgeo_ov%rhoc
    if (mgeo_ov%override_qinp) qinp = mgeo_ov%qinp
    if (mgeo_ov%override_shat) shat = mgeo_ov%shat
    if (mgeo_ov%override_rgeo_lcfs) r_geo = mgeo_ov%rgeo_lcfs
    if (mgeo_ov%override_rgeo_local) rmaj = mgeo_ov%rgeo_local
    if (mgeo_ov%override_geoType) geoType = mgeo_ov%geoType
    if (mgeo_ov%override_aSurf) aSurf = mgeo_ov%aSurf
    if (mgeo_ov%override_shift) shift = mgeo_ov%shift
    if (mgeo_ov%override_shiftVert) shiftVert = mgeo_ov%shiftVert
    if (mgeo_ov%override_mMode) mMode = mgeo_ov%mMode
    if (mgeo_ov%override_nMode) nMode = mgeo_ov%nMode
    if (mgeo_ov%override_deltam) then
      deltam = mgeo_ov%deltam
    else if (mgeo_ov%override_akappa) then
      deltam = mgeo_ov%akappa
    end if
    if (mgeo_ov%override_deltan) then
      deltan = mgeo_ov%deltan
    else if (mgeo_ov%override_tri) then
      deltan = mgeo_ov%tri
    end if
    if (mgeo_ov%override_deltampri) then
      deltampri = mgeo_ov%deltampri
    else if (mgeo_ov%override_akappri) then
      deltampri = mgeo_ov%akappri
    end if
    if (mgeo_ov%override_deltanpri) then
      deltanpri = mgeo_ov%deltanpri
    else if (mgeo_ov%override_tripri) then
      deltanpri = mgeo_ov%tripri
    end if
    if (mgeo_ov%override_thetam) thetam = mgeo_ov%thetam
    if (mgeo_ov%override_thetan) thetan = mgeo_ov%thetan
    if (mgeo_ov%override_betaprim) betaprim = mgeo_ov%betaprim
  end subroutine set_overrides

  !> Set the module level config type
  !> Will abort if the module has already been initialised to avoid
  !> inconsistencies.
  subroutine set_theta_grid_parameters_config(theta_grid_parameters_config_in)
    use mp, only: mp_abort
    type(theta_grid_parameters_config_type), intent(in), optional :: theta_grid_parameters_config_in
    if (initialized) then
       call mp_abort("Trying to set theta_grid_parameters_config when already initialized.", to_screen = .true.)
    end if
    if (present(theta_grid_parameters_config_in)) then
       theta_grid_parameters_config = theta_grid_parameters_config_in
    end if
  end subroutine set_theta_grid_parameters_config

#include "theta_grid_parameters_auto_gen.inc"
end module theta_grid_params