Replace q with just q95

examples/data/csv_output_large_tokamak_MFILE.DAT

@@ -1425,7 +1425,7 @@ dr_cs = 0.5
 * q
 ixc = 18
 boundl(18) = 3.0
-q = 3.5
+q95 = 3.5
 
 * Machine dr_bore [m]
 ixc = 29

examples/data/large_tokamak_1_MFILE.DAT

@@ -1419,7 +1419,7 @@ dr_cs = 0.5
 * q
 ixc = 18
 boundl(18) = 3.0
-q = 3.5
+q95 = 3.5
 
 * Machine dr_bore [m]
 ixc = 29

examples/data/large_tokamak_2_MFILE.DAT

@@ -1419,7 +1419,7 @@ dr_cs = 0.5
 * q
 ixc = 18
 boundl(18) = 3.0
-q = 3.5
+q95 = 3.5
 
 * Machine dr_bore [m]
 ixc = 29

examples/data/large_tokamak_3_MFILE.DAT

@@ -1419,7 +1419,7 @@ dr_cs = 0.5
 * q
 ixc = 18
 boundl(18) = 3.0
-q = 3.5
+q95 = 3.5
 
 * Machine dr_bore [m]
 ixc = 29

examples/data/large_tokamak_4_MFILE.DAT

@@ -1419,7 +1419,7 @@ dr_cs = 0.5
 * q
 ixc = 18
 boundl(18) = 3.0
-q = 3.5
+q95 = 3.5
 
 * Machine dr_bore [m]
 ixc = 29

examples/data/large_tokamak_IN.DAT

@@ -226,10 +226,10 @@ ixc = 16
 boundl(16) = 0.3
 dr_cs = 0.5
 
-* q
+* q95
 ixc = 18
 boundl(18) = 3.0
-q = 3.5
+q95 = 3.5
 
 * Machine dr_bore [m]
 ixc = 29

examples/data/scan_MFILE.DAT

@@ -9259,7 +9259,7 @@ i_plasma_geometry   = 0 * Switch for plasma cross-sectional shape calculation: u
 *kappa = 1.7808
 kappa = 1.848
 triang   = 0.5 * Plasma separatrix triangularity (calculated if i_plasma_geometry=1; 3 or 4)
-q        = 3.247 * Safety factor 'near' plasma edge (iteration variable 18);
+q95        = 3.247 * Safety factor 'near' plasma edge (iteration variable 18);
 q0       = 1.0 * Safety factor on axis
 rmajor   = 9.072 * Plasma major radius (m) (iteration variable 3)
 i_single_null    = 1 * Switch for single null / double null plasma;
@@ -9323,7 +9323,7 @@ t_burn    = 1.0d4 * Burn time (s) (calculated if i_pulsed_plant=1)
                      dr_tf_inboard          =  1.2080E+00
                      fwalld         =  1.3100E-01
                      dr_cs          =  5.5242E-01
-                     q              =  3.5000E+00
+                     q95              =  3.5000E+00
                      dr_bore           =  2.3322E+00
                      fbeta_max       =  4.8251E-01
                      coheof         =  2.0726E+07

examples/data/scan_example_file_IN.DAT

@@ -226,10 +226,10 @@ ixc = 16
 boundl(16) = 0.3
 dr_cs = 0.5
 
-* Safety factor near plasma edge
+* Safety factor at 95% flux surface 
 ixc = 18
 boundl(18) = 3.0
-q = 3.5
+q95 = 3.5
 
 * Machine dr_bore [m]
 ixc = 29

process/divertor.py

@@ -116,15 +116,18 @@ def run(self, output: bool) -> None:
         #  (2.5 factor comes from normalization to ITER 1990)
 
         tconl = (
-            2.5e0 * pv.rmajor * pv.q * (1.0e0 + 1.0e0 / (pv.q * pv.aspect) ** 2) ** 0.5
+            2.5e0
+            * pv.rmajor
+            * pv.q95
+            * (1.0e0 + 1.0e0 / (pv.q95 * pv.aspect) ** 2) ** 0.5
         )
         dtheta = plsep / pv.rminor
         dconl = (
             2.5e0
             * bv.rspo
-            * pv.q
+            * pv.q95
             * dtheta
-            * (1.0e0 + 1.0e0 / (pv.q * pv.aspect) ** 2) ** 0.5
+            * (1.0e0 + 1.0e0 / (pv.q95 * pv.aspect) ** 2) ** 0.5
         )
         rconl = dconl / tconl
 

process/io/obsolete_vars.py

@@ -201,6 +201,7 @@
     "rli": "ind_plasma_internal_norm",
     "gamma": "ejima_coeff",
     "lpulse": "i_pulsed_plant",
+    "q": "q95",
 }
 
 OBS_VARS_HELP = {

process/physics.py

@@ -1545,8 +1545,6 @@ def physics(self):
           https://inis.iaea.org/search/search.aspx?orig_q=RN:45031642
         """
 
-        physics_variables.q95 = physics_variables.q
-
         # Calculate plasma composition
         # Issue #261 Remove old radiation model (imprad_model=0)
         self.plasma_composition()
@@ -1583,7 +1581,7 @@ def physics(self):
             physics_variables.p0,
             physics_variables.len_plasma_poloidal,
             physics_variables.q0,
-            physics_variables.q,
+            physics_variables.q95,
             physics_variables.ind_plasma_internal_norm,
             physics_variables.rmajor,
             physics_variables.rminor,
@@ -3215,7 +3213,7 @@ def calculate_plasma_current(
                 9 = FIESTA ST scaling
             iprofile (int): Switch for current profile consistency.
                 0: Use input values for alphaj, ind_plasma_internal_norm, beta_norm_max.
-                1: Make these consistent with input q, q_0 values.
+                1: Make these consistent with input q95, q_0 values.
                 2: Use input values for alphaj, ind_plasma_internal_norm. Scale beta_norm_max with aspect ratio (original scaling).
                 3: Use input values for alphaj, ind_plasma_internal_norm. Scale beta_norm_max with aspect ratio (Menard scaling).
                 4: Use input values for alphaj, beta_norm_max. Set ind_plasma_internal_norm from elongation (Menard scaling).
@@ -3679,12 +3677,12 @@ def outplas(self):
                 if physics_variables.iprofile == 1:
                     po.ocmmnt(
                         self.outfile,
-                        "Consistency between q0,q,alphaj,ind_plasma_internal_norm,beta_norm_max is enforced",
+                        "Consistency between q0,q95,alphaj,ind_plasma_internal_norm,beta_norm_max is enforced",
                     )
                 else:
                     po.ocmmnt(
                         self.outfile,
-                        "Consistency between q0,q,alphaj,ind_plasma_internal_norm,beta_norm_max is not enforced",
+                        "Consistency between q0,q95,alphaj,ind_plasma_internal_norm,beta_norm_max is not enforced",
                     )
 
                 po.oblnkl(self.outfile)
@@ -3769,7 +3767,10 @@ def outplas(self):
 
             if physics_variables.i_plasma_current == 2:
                 po.ovarrf(
-                    self.outfile, "Mean edge safety factor", "(q)", physics_variables.q
+                    self.outfile,
+                    "Mean edge safety factor",
+                    "(q95)",
+                    physics_variables.q95,
                 )
 
             po.ovarrf(
@@ -3790,7 +3791,7 @@ def outplas(self):
             if physics_variables.i_plasma_geometry == 1:
                 po.ovarrf(
                     self.outfile,
-                    "Lower limit for edge safety factor q",
+                    "Lower limit for edge safety factor q95",
                     "(q95_min)",
                     physics_variables.q95_min,
                     "OP ",
@@ -5891,7 +5892,7 @@ def output_confinement_comparison(self, istell: int) -> None:
                 physics_variables.rminor,
                 physics_variables.ten,
                 physics_variables.tin,
-                physics_variables.q,
+                physics_variables.q95,
                 physics_variables.qstar,
                 physics_variables.vol_plasma,
                 physics_variables.zeff,
@@ -6210,7 +6211,8 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float:
         # inverse_q = 1/safety factor
         # Parabolic q profile assumed
         inverse_q = 1 / (
-            physics_variables.q0 + (physics_variables.q - physics_variables.q0) * roa**2
+            physics_variables.q0
+            + (physics_variables.q95 - physics_variables.q0) * roa**2
         )
         # Create new array of average mass of fuel portion of ions
         amain = np.full_like(inverse_q, physics_variables.m_fuel_amu)
@@ -6690,7 +6692,7 @@ def fhz(hfact: float) -> float:
                 physics_variables.rminor,
                 physics_variables.ten,
                 physics_variables.tin,
-                physics_variables.q,
+                physics_variables.q95,
                 physics_variables.qstar,
                 physics_variables.vol_plasma,
                 physics_variables.zeff,
@@ -6745,7 +6747,7 @@ def calculate_confinement_time(
         rminor: float,
         ten: float,
         tin: float,
-        q: float,
+        q95: float,
         qstar: float,
         vol_plasma: float,
         zeff: float,
@@ -6770,7 +6772,7 @@ def calculate_confinement_time(
         :param pinjmw: Auxiliary power to ions and electrons (MW)
         :param plasma_current: Plasma current (A)
         :param pcoreradpv: Total core radiation power (MW/m3)
-        :param q: Edge safety factor (tokamaks), or rotational transform iotabar (stellarators)
+        :param q95: Edge safety factor (tokamaks), or rotational transform iotabar (stellarators)
         :param qstar: Equivalent cylindrical edge safety factor
         :param rmajor: Plasma major radius (m)
         :param rminor: Plasma minor radius (m)
@@ -7129,7 +7131,7 @@ def calculate_confinement_time(
                     dnla20,
                     bt,
                     p_plasma_loss_mw,
-                    q,
+                    q95,
                 )
             )
 
@@ -7324,7 +7326,8 @@ def calculate_confinement_time(
 
         # ISS95 stellarator scaling
         elif i_confinement_time == 37:
-            iotabar = q  # dummy argument q is actual argument iotabar for stellarators
+            # dummy argument q95 is actual argument iotabar for stellarators
+            iotabar = q95
             t_electron_confinement = confinement.iss95_stellarator_confinement_time(
                 rminor,
                 rmajor,
@@ -7338,7 +7341,8 @@ def calculate_confinement_time(
 
         # ISS04 stellarator scaling
         elif i_confinement_time == 38:
-            iotabar = q  # dummy argument q is actual argument iotabar for stellarators
+            # dummy argument q95 is actual argument iotabar for stellarators
+            iotabar = q95
             t_electron_confinement = confinement.iss04_stellarator_confinement_time(
                 rminor,
                 rmajor,
@@ -7401,7 +7405,7 @@ def calculate_confinement_time(
                 p_plasma_loss_mw,
                 rmajor,
                 rminor,
-                q,
+                q95,
                 qstar,
                 aspect,
                 m_fuel_amu,

process/stellarator.py

@@ -4121,8 +4121,6 @@ def stphys(self, output):
             * physics_variables.rmajor
         )
 
-        physics_variables.q95 = physics_variables.q
-
         #  Calculate poloidal field using rotation transform
         physics_variables.bp = (
             physics_variables.rminor

source/fortran/constraint_equations.f90

@@ -1981,14 +1981,14 @@ subroutine constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       !! and hence also optional here.
       !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
       !! fq : input real : f-value for edge safety factor
-      !! q : safety factor 'near' plasma edge: equal to q95
+      !! q95 : safety factor 'near' plasma edge
       !! (unless i_plasma_current = 2 (ST current scaling), in which case q = mean edge safety factor qbar)
       !! q95_min : input real :  lower limit for edge safety factor
       !! itart : input integer : switch for spherical tokamak (ST) models:<UL>
       !! <LI> = 0 use conventional aspect ratio models;
       !! <LI> = 1 use spherical tokamak models</UL>
       use constraint_variables, only: fq
-      use physics_variables, only: q, q95_min, itart
+      use physics_variables, only: q95, q95_min, itart
       implicit none
             real(dp), intent(out) :: tmp_cc
       real(dp), intent(out) :: tmp_con
@@ -1998,7 +1998,7 @@ subroutine constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
 
       ! if the machine isn't a ST then report error
       if (itart == 0) call report_error(9)
-      tmp_cc =   1.0D0 - fq * q/q95_min
+      tmp_cc =   1.0D0 - fq * q95/q95_min
       tmp_con = q95_min * (1.0D0 - tmp_cc)
       tmp_err = q95_min * tmp_cc
       tmp_symbol = '<'

source/fortran/input.f90

@@ -308,7 +308,7 @@ subroutine parse_input_file(in_file,out_file,show_changes)
       fpdivlim, beta_poloidal_eps_max, i_confinement_time, kappa95, aspect, f_r_conducting_wall, nesep, c_beta, csawth, dene, &
       ftar, plasma_res_factor, f_sync_reflect, f_nd_beam_electron, beta, neped, hfact, beta_norm_max, &
       fgwsep, rhopedn, tratio, q0, i_plasma_geometry, i_plasma_shape, fne0, ignite, f_tritium, &
-      i_beta_fast_alpha, tauee_in, alphaj, alphat, i_plasma_current, q, ti, tesep, ind_plasma_internal_norm, triang, &
+      i_beta_fast_alpha, tauee_in, alphaj, alphat, i_plasma_current, q95, ti, tesep, ind_plasma_internal_norm, triang, &
       itart, f_nd_alpha_electron, iprofile, triang95, rad_fraction_sol, betbm0, f_nd_protium_electrons, &
       teped, f_helium3, iwalld, ejima_coeff, f_alpha_plasma, fgwped, tbeta, i_bootstrap_current, &
       i_rad_loss, te, alphan, rmajor, plasma_square, kappa, fkzohm, beamfus0, &
@@ -700,9 +700,9 @@ subroutine parse_input_file(in_file,out_file,show_changes)
        case ('plasma_res_factor')
           call parse_real_variable('plasma_res_factor', plasma_res_factor, 0.0D0, 1.0D0, &
                'Plasma resistivity pre-factor')
-       case ('q')
-          call parse_real_variable('q', q, 1.00D0, 50.0D0, &
-               'Safety factor near plasma edge')
+      case ('q95')
+         call parse_real_variable('q95', q95, 1.00D0, 50.0D0, &
+               'Safety factor at 95% flux surface')
        case ('q0')
           call parse_real_variable('q0', q0, 0.01D0, 20.0D0, &
                'Safety factor on axis')

source/fortran/iteration_variables.f90

@@ -437,25 +437,25 @@ end subroutine set_itv_17
   !---------------------------------
 
   subroutine init_itv_18
-    !! <LI> (18) q
+    !! <LI> (18) q95
     use numerics, only: lablxc, boundl, boundu
     implicit none
-    lablxc(18) = 'q             '
+    lablxc(18) = 'q95             '
     boundl(18) = 2.000D0
     boundu(18) = 50.00D0
   end subroutine init_itv_18
 
   real(kind(1.d0)) function itv_18()
-    use physics_variables, only: q
+    use physics_variables, only: q95
     implicit none
-    itv_18 = q
+    itv_18 = q95
   end function itv_18
 
   subroutine set_itv_18(ratio)
-    use physics_variables, only: q
+    use physics_variables, only: q95
     implicit none
     real(kind(1.d0)) :: ratio
-    q = ratio
+    q95 = ratio
   end subroutine set_itv_18
 
   !---------------------------------

source/fortran/physics_variables.f90

@@ -400,7 +400,7 @@ module physics_variables
   !! switch for current profile consistency:
   !!
   !! - =0 use input values for alphaj, ind_plasma_internal_norm, beta_norm_max
-  !! - =1 make these consistent with input q, q_0 values (recommend `i_plasma_current=4` with this option)
+  !! - =1 make these consistent with input q95, q_0 values (recommend `i_plasma_current=4` with this option)
   !! - =2 use input values for alphaj, ind_plasma_internal_norm. Scale beta_norm_max with aspect ratio (original scaling)
   !! - =3 use input values for alphaj, ind_plasma_internal_norm. Scale beta_norm_max with aspect ratio (Menard scaling)
   !! - =4 use input values for alphaj, beta_norm_max. Set ind_plasma_internal_norm from elongation (Menard scaling)
@@ -728,15 +728,12 @@ module physics_variables
   real(dp) :: pden_ion_transport_loss_mw
   !! ion transport power per volume (MW/m3)
 
-  real(dp) :: q
-  !! Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95
-  !! (unless `i_plasma_current=2` (ST current scaling), in which case q = mean edge safety factor qbar)
-
   real(dp) :: q0
-  !! safety factor on axis
+  !! Safety factor on axis
 
   real(dp) :: q95
-  !! safety factor at 95% surface
+  !! Safety factor at 95% flux surface (iteration variable 18) (unless icurr=2 (ST current scaling),
+  !! in which case q95 = mean edge safety factor qbar)
 
   real(dp) :: qfuel
   !! plasma fuelling rate (nucleus-pairs/s)
@@ -1072,7 +1069,6 @@ subroutine init_physics_variables
     p_ion_transport_loss_mw = 0.0D0
     pscalingmw = 0.0D0
     pden_ion_transport_loss_mw = 0.0D0
-    q = 3.0D0
     q0 = 1.0D0
     q95 = 0.0D0
     qfuel = 0.0D0