!> 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