!> 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, dimension(mxh_max_moments) :: c_mxh, s_mxh, dc_mxh_dr, ds_mxh_dr logical :: initialized = .false. real :: kp = -1. !> 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 !> Deprecataed input for zeroth cosine moment for the MXH local !> equilibrum. Now set as first element in [[c_mxh]], real :: c0_mxh = 0.0 !> Cosine moments for the MXH local equilibrum. Note the first !> element is the N=0 component. 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 !> 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 !> 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). 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 including the !> N=0 component. 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 !> 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 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 use warning_helpers, only: not_exactly_equal, is_not_zero implicit none type(theta_grid_parameters_config_type), intent(in), optional :: theta_grid_parameters_config_in integer, parameter :: verb=3 real :: c0_mxh 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 associate(self => theta_grid_parameters_config) #include "theta_grid_parameters_copy_out_auto_gen.inc" end associate if (kp > 0.) pk = 2.*kp ! Make sure ntheta is even ntheta = 2 * (ntheta / 2) ! 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 (not_exactly_equal(deltam, 1.0)) write (*,*) "WARNING: ignoring value of deltam, did you mean akappa?" if (not_exactly_equal(deltan, 1.0)) write (*,*) "WARNING: ignoring value of deltan, did you mean tri?" if (not_exactly_equal(delta2, 1.0)) write (*,*) "WARNING: ignoring value of delta2, did you mean akappa?" if (not_exactly_equal(delta3, 1.0)) write (*,*) "WARNING: ignoring value of delta3, did you mean tri?" if (is_not_zero(deltampri)) write (*,*) "WARNING: ignoring value of deltampri, did you mean akappri?" if (is_not_zero(deltanpri)) write (*,*) "WARNING: ignoring value of deltanpri, did you mean tripri?" if (is_not_zero(thetam)) write (*,*) "WARNING: ignoring value of thetam, did you mean thetak?" if (is_not_zero(thetan)) write (*,*) "WARNING: ignoring value of thetan, did you mean thetad?" if (is_not_zero(theta2)) write (*,*) "WARNING: ignoring value of theta2, did you mean thetak?" if (is_not_zero(theta3)) write (*,*) "WARNING: ignoring value of theta3, did you mean thetad?" if (not_exactly_equal(aSurf, 1.0)) write (*,*) "WARNING: ignoring value of aSurf, value not needed" if (not_exactly_equal(raxis, 3.0)) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?" if (is_not_zero(zaxis)) 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 (not_exactly_equal(deltam, 1.0)) write (*,*) "WARNING: ignoring value of deltam, did you mean delta2?" if (not_exactly_equal(deltan, 1.0)) write (*,*) "WARNING: ignoring value of deltan, did you mean delta3?" if (not_exactly_equal(akappa, 1.0)) write (*,*) "WARNING: ignoring value of akappa, did you mean delta2?" if (is_not_zero(tri)) write (*,*) "WARNING: ignoring value of tri, did you mean delta3?" if (is_not_zero(deltampri)) write (*,*) "WARNING: ignoring value of deltampri, value not needed" if (is_not_zero(deltanpri)) write (*,*) "WARNING: ignoring value of deltanpri, value not needed" if (is_not_zero(akappri)) write (*,*) "WARNING: ignoring value of akappri, value not needed" if (is_not_zero(tripri)) write (*,*) "WARNING: ignoring value of tripri, value not needed" if (is_not_zero(thetam)) write (*,*) "WARNING: ignoring value of thetam, did you mean theta2?" if (is_not_zero(thetan)) write (*,*) "WARNING: ignoring value of thetan, did you mean theta3?" if (is_not_zero(thetak)) write (*,*) "WARNING: ignoring value of thetak, did you mean theta2?" if (is_not_zero(thetad)) write (*,*) "WARNING: ignoring value of thetad, did you mean theta3?" if (is_not_zero(shift)) write (*,*) "WARNING: ignoring value of shift, did you mean Raxis?" if (is_not_zero(shiftVert)) 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 (not_exactly_equal(delta2, 1.0)) write (*,*) "WARNING: ignoring value of delta2, did you mean deltam?" if (not_exactly_equal(delta3, 1.0)) write (*,*) "WARNING: ignoring value of delta3, did you mean deltan?" if (not_exactly_equal(akappa, 1.0)) write (*,*) "WARNING: ignoring value of akappa, did you mean deltam?" if (is_not_zero(tri)) write (*,*) "WARNING: ignoring value of tri, did you mean deltan?" if (is_not_zero(akappri)) write (*,*) "WARNING: ignoring value of akappri, did you mean deltampri?" if (is_not_zero(tripri)) write (*,*) "WARNING: ignoring value of tripri, did you mean deltanpri?" if (is_not_zero(theta2)) write (*,*) "WARNING: ignoring value of theta2, did you mean thetam?" if (is_not_zero(theta3)) write (*,*) "WARNING: ignoring value of theta3, did you mean thetan?" if (is_not_zero(thetak)) write (*,*) "WARNING: ignoring value of thetak, did you mean thetam?" if (is_not_zero(thetad)) write (*,*) "WARNING: ignoring value of thetad, did you mean thetan?" if (not_exactly_equal(aSurf, 1.0)) write (*,*) "WARNING: ignoring value of aSurf, value not needed" if (not_exactly_equal(raxis, 3.0)) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?" if (is_not_zero(zaxis)) write (*,*) "WARNING: ignoring value of zaxis, did you mean shiftVert?" case (geo_type_fourier_series) if (not_exactly_equal(delta2, 1.0)) write (*,*) "WARNING: ignoring value of delta2, did you mean deltam?" if (not_exactly_equal(delta3, 1.0)) write (*,*) "WARNING: ignoring value of delta3, did you mean deltan?" if (not_exactly_equal(akappa, 1.0)) write (*,*) "WARNING: ignoring value of akappa, did you mean deltam?" if (is_not_zero(tri)) write (*,*) "WARNING: ignoring value of tri, did you mean deltan?" if (is_not_zero(akappri)) write (*,*) "WARNING: ignoring value of akappri, did you mean deltampri?" if (is_not_zero(tripri)) write (*,*) "WARNING: ignoring value of tripri, did you mean deltanpri?" if (is_not_zero(theta2)) write (*,*) "WARNING: ignoring value of theta2, did you mean thetam?" if (is_not_zero(theta3)) write (*,*) "WARNING: ignoring value of theta3, did you mean thetan?" if (is_not_zero(thetak)) write (*,*) "WARNING: ignoring value of thetak, did you mean thetam?" if (is_not_zero(thetad)) write (*,*) "WARNING: ignoring value of thetad, did you mean thetan?" if (not_exactly_equal(aSurf, 1.0)) write (*,*) "WARNING: ignoring value of aSurf, value not needed" if (not_exactly_equal(raxis, 3.0)) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?" if (is_not_zero(zaxis)) write (*,*) "WARNING: ignoring value of zaxis, did you mean shiftVert?" case (geo_type_miller_extended_harmonic) if (not_exactly_equal(deltam, 1.0)) write (*,*) "WARNING: ignoring value of deltam, did you mean akappa?" if (is_not_zero(deltampri)) 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 (not_exactly_equal(deltan, 1.0)) write (*,*) "WARNING: ignoring value of deltan, value not needed" if (not_exactly_equal(delta2, 1.0)) write (*,*) "WARNING: ignoring value of delta2, value not needed" if (not_exactly_equal(delta3, 1.0)) write (*,*) "WARNING: ignoring value of delta3, value not needed" if (is_not_zero(deltampri)) write (*,*) "WARNING: ignoring value of deltampri, value not needed" if (is_not_zero(deltanpri)) write (*,*) "WARNING: ignoring value of deltanpri, value not needed" if (is_not_zero(thetam)) write (*,*) "WARNING: ignoring value of thetam, value not needed" if (is_not_zero(thetan)) write (*,*) "WARNING: ignoring value of thetan, value not needed" if (is_not_zero(theta2)) write (*,*) "WARNING: ignoring value of theta2, value not needed" if (is_not_zero(theta3)) write (*,*) "WARNING: ignoring value of theta3, value not needed" if (not_exactly_equal(aSurf, 1.0)) write (*,*) "WARNING: ignoring value of aSurf, value not needed" if (not_exactly_equal(raxis, 3.0)) write (*,*) "WARNING: ignoring value of raxis, did you mean shift?" if (is_not_zero(zaxis)) 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 if (is_not_zero(theta_grid_parameters_config%c0_mxh)) then write(*,'("Warning: c0_mxh input is deprecated, set via c_mxh now.")') 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, 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 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 if (.not. mgeo_ov%is_initialised()) return 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 #include "theta_grid_parameters_auto_gen.inc" end module theta_grid_params