From 3e5c9a3db5093484e32e946dbe41c81344b246ed Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Mon, 3 Mar 2025 15:43:53 +0000
Subject: [PATCH 1/7] Use consistent form in constraint equations

Ensures min and max inequality constraints are violated proportionally to each other.
---
 source/fortran/constraint_equations.f90 | 151 +++++++++++++-----------
 1 file changed, 81 insertions(+), 70 deletions(-)

diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90
index 89b64c3c3a..877401121f 100755
--- a/source/fortran/constraint_equations.f90
+++ b/source/fortran/constraint_equations.f90
@@ -647,13 +647,13 @@ subroutine constraint_eqn_005(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
 
 	   ! Apply Greenwald limit to line-averaged density
       if (i_density_limit == 7) then
-         tmp_cc     = 1.0D0 - fdene * dnelimt/dnla
+         tmp_cc     = dnla/dnelimt - 1.0D0 * fdene
          tmp_con    = fdene * dnelimt
          tmp_err    = fdene * dnelimt - dnla
          tmp_symbol = '<'
          tmp_units  = '/m3'
       else
-         tmp_cc = 1.0D0 - fdene * dnelimt/dene
+         tmp_cc =  dene/dnelimt - 1.0D0 * fdene
          tmp_con    = dnelimt * (1.0D0 - tmp_cc)
          tmp_err    = dene * tmp_cc
          tmp_symbol = '<'
@@ -685,7 +685,7 @@ subroutine constraint_eqn_006(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fbeta_poloidal_eps * beta_poloidal_eps_max/(eps*beta_poloidal)
+      tmp_cc =  (eps*beta_poloidal)/beta_poloidal_eps_max - 1.0D0 * fbeta_poloidal_eps
       tmp_con = beta_poloidal_eps_max * (1.0D0 - tmp_cc)
       tmp_err = (eps*beta_poloidal) * tmp_cc
       tmp_symbol = '<'
@@ -760,7 +760,7 @@ subroutine constraint_eqn_008(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fwalld * walalw/pflux_fw_neutron_mw
+      tmp_cc =  pflux_fw_neutron_mw/walalw - 1.0D0 * fwalld
       tmp_con = fwalld * walalw
       tmp_err = fwalld * walalw - pflux_fw_neutron_mw
       tmp_symbol = '<'
@@ -790,7 +790,7 @@ subroutine constraint_eqn_009(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - ffuspow * powfmax/fusion_power
+      tmp_cc =  fusion_power/powfmax - 1.0D0 * ffuspow
       tmp_con = powfmax * (1.0D0 - tmp_cc)
       tmp_err = fusion_power * tmp_cc
       tmp_symbol = '<'
@@ -885,7 +885,7 @@ subroutine constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 + fvs * vs_cs_pf_total_pulse/vs_plasma_total_required
+      tmp_cc =  1.0D0 - fvs * vs_cs_pf_total_pulse/vs_plasma_total_required
       tmp_con = vs_plasma_total_required * (1.0D0 - tmp_cc)
       tmp_err = vs_plasma_total_required * tmp_cc
       tmp_symbol = '>'
@@ -974,7 +974,7 @@ subroutine constraint_eqn_015(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  -(1.0D0 - fl_h_threshold * p_l_h_threshold_mw / pdivt)
+      tmp_cc =  1.0D0 - fl_h_threshold * pdivt / p_l_h_threshold_mw
       tmp_con = p_l_h_threshold_mw
       tmp_err = p_l_h_threshold_mw - pdivt / fl_h_threshold
       if (fl_h_threshold > 1.0D0) then
@@ -1048,7 +1048,7 @@ subroutine constraint_eqn_017(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       !! Maximum possible power/vol_plasma that can be radiated (local)
 
       pradmaxpv = pinjmw/vol_plasma + alpha_power_density_total*f_alpha_plasma + charged_power_density + pden_plasma_ohmic_mw
-      tmp_cc =  1.0D0 - fradpwr * pradmaxpv / pden_plasma_rad_mw
+      tmp_cc =  pden_plasma_rad_mw/pradmaxpv - 1.0D0 * fradpwr
       tmp_con = pradmaxpv * (1.0D0 - tmp_cc)
       tmp_err = pden_plasma_rad_mw * tmp_cc
       tmp_symbol = '<'
@@ -1078,7 +1078,7 @@ subroutine constraint_eqn_018(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fhldiv * hldivlim/hldiv
+      tmp_cc =  hldiv/hldivlim - 1.0D0 * fhldiv
       tmp_con = hldivlim * (1.0D0 - tmp_cc)
       tmp_err = hldiv * tmp_cc
       tmp_symbol = '<'
@@ -1112,7 +1112,7 @@ subroutine constraint_eqn_019(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       real(dp) :: totmva
 
       totmva = tfcpmw + tflegmw
-      tmp_cc =  1.0D0 - fmva * mvalim/totmva
+      tmp_cc =  totmva/mvalim - 1.0D0 * fmva
       tmp_con = mvalim * (1.0D0 - tmp_cc)
       tmp_err = totmva * tmp_cc
       tmp_symbol = '<'
@@ -1142,7 +1142,7 @@ subroutine constraint_eqn_020(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fportsz * rtanmax/rtanbeam
+      tmp_cc =  rtanbeam/rtanmax - 1.0D0 * fportsz
       tmp_con = rtanmax * (1.0D0 - tmp_cc)
       tmp_err = rtanbeam * tmp_cc
       tmp_symbol = '<'
@@ -1203,7 +1203,7 @@ subroutine constraint_eqn_022(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fdivcol * rlenmax / rlclolcn
+      tmp_cc =  rlclolcn / rlenmax - 1.0D0 * fdivcol
       tmp_con = rlenmax * (1.0D0 - tmp_cc)
       tmp_err = rlclolcn * tmp_cc
       tmp_symbol = '<'
@@ -1240,7 +1240,7 @@ subroutine constraint_eqn_023(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       real(dp) :: rcw
 
       rcw = rminor + dr_fw_plasma_gap_outboard + dr_fw_outboard + dr_blkt_outboard
-      tmp_cc =  1.0D0 -fr_conducting_wall * f_r_conducting_wall*rminor / rcw
+      tmp_cc =  rcw / (f_r_conducting_wall * rminor) - 1.0D0 * fr_conducting_wall
       tmp_con = f_r_conducting_wall*rminor * (1.0D0 - tmp_cc)
       tmp_err = rcw * tmp_cc
       tmp_symbol = '<'
@@ -1285,28 +1285,28 @@ subroutine constraint_eqn_024(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
 
       ! Include all beta components: relevant for both tokamaks and stellarators
       if ((i_beta_component == 0).or.(istell /= 0)) then
-         tmp_cc =  1.0D0 - fbeta_max * beta_max/beta
+         tmp_cc =  beta/beta_max - 1.0D0 * fbeta_max
          tmp_con = beta_max
          tmp_err = beta_max - beta / fbeta_max
          tmp_symbol = '<'
          tmp_units = ''
       ! Here, the beta limit applies to only the thermal component, not the fast alpha or neutral beam parts
       else if (i_beta_component == 1) then
-         tmp_cc = 1.0D0 - fbeta_max * beta_max/(beta-beta_fast_alpha-beta_beam)
+         tmp_cc = (beta-beta_fast_alpha-beta_beam)/beta_max - 1.0D0 * fbeta_max
          tmp_con = beta_max
          tmp_err = beta_max - (beta-beta_fast_alpha-beta_beam) / fbeta_max
          tmp_symbol = '<'
          tmp_units = ''
       ! Beta limit applies to thermal + neutral beam: components of the total beta, i.e. excludes alphas
       else if (i_beta_component == 2) then
-         tmp_cc = 1.0D0 - fbeta_max * beta_max/(beta-beta_fast_alpha)
+         tmp_cc = (beta-beta_fast_alpha)/beta_max - 1.0D0 * fbeta_max
          tmp_con = beta_max * (1.0D0 - tmp_cc)
          tmp_err = (beta-beta_fast_alpha) * tmp_cc
          tmp_symbol = '<'
          tmp_units = ''
       ! Beta limit applies to toroidal beta
       else if (i_beta_component == 3) then
-         tmp_cc =  1.0D0 - fbeta_max * beta_max/(beta*(btot/bt)**2)
+         tmp_cc =  (beta*(btot/bt)**2)/beta_max - 1.0D0 * fbeta_max
          tmp_con = beta_max
          tmp_err = beta_max - (beta*(btot/bt)**2) / fbeta_max
          tmp_symbol = '<'
@@ -1337,7 +1337,7 @@ subroutine constraint_eqn_025(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fpeakb * bmxlim/bmaxtf
+      tmp_cc =  bmaxtf/bmxlim - 1.0D0 * fpeakb
       tmp_con = bmxlim * (1.0D0 - tmp_cc)
       tmp_err = bmaxtf * tmp_cc
       tmp_symbol = '<'
@@ -1367,7 +1367,7 @@ subroutine constraint_eqn_026(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fjohc * j_cs_critical_flat_top_end/j_cs_flat_top_end
+      tmp_cc =  j_cs_flat_top_end/j_cs_critical_flat_top_end - 1.0D0 * fjohc
       tmp_con = j_cs_critical_flat_top_end
       tmp_err = j_cs_critical_flat_top_end - j_cs_flat_top_end / fjohc
       tmp_symbol = '<'
@@ -1397,7 +1397,7 @@ subroutine constraint_eqn_027(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fjohc0 * j_cs_critical_pulse_start/j_cs_pulse_start
+      tmp_cc =  j_cs_pulse_start/j_cs_critical_pulse_start - 1.0D0 * fjohc0
       tmp_con = j_cs_critical_pulse_start
       tmp_err = j_cs_critical_pulse_start - j_cs_pulse_start / fjohc0
       tmp_symbol = '<'
@@ -1506,7 +1506,7 @@ subroutine constraint_eqn_030(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fpinj*pinjalw/pinjmw
+      tmp_cc =  pinjmw/pinjalw - 1.0D0 * fpinj
       tmp_con = pinjalw
       tmp_err = pinjalw  - pinjmw / fpinj
       tmp_symbol = '<'
@@ -1536,7 +1536,7 @@ subroutine constraint_eqn_031(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fstrcase * sig_tf_case_max/sig_tf_case
+      tmp_cc =  sig_tf_case/sig_tf_case_max - 1.0D0 * fstrcase
       tmp_con = sig_tf_case_max
       tmp_err = sig_tf_case_max - sig_tf_case / fstrcase
       tmp_symbol = '<'
@@ -1566,7 +1566,7 @@ subroutine constraint_eqn_032(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fstrcond * sig_tf_wp_max/sig_tf_wp
+      tmp_cc =  sig_tf_wp/sig_tf_wp_max - 1.0D0 * fstrcond
       tmp_con = sig_tf_wp_max
       tmp_err = sig_tf_wp_max - sig_tf_wp / fstrcond
       tmp_symbol = '<'
@@ -1597,7 +1597,7 @@ subroutine constraint_eqn_033(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=10), intent(out) :: tmp_units
 
       if (fiooic > 0.7D0) call report_error(285)
-      tmp_cc =  1.0D0 - fiooic * jwdgcrt/jwptf
+      tmp_cc =  jwptf/jwdgcrt - 1.0D0 * fiooic
       tmp_con = jwdgcrt * (1.0D0 - tmp_cc)
       tmp_err = jwptf * tmp_cc
       tmp_symbol = '<'
@@ -1627,7 +1627,7 @@ subroutine constraint_eqn_034(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fvdump * vdalw/vtfskv
+      tmp_cc =  vtfskv/vdalw - 1.0D0 * fvdump
       tmp_con = vdalw
       tmp_err = vdalw - vtfskv
       tmp_symbol = '<'
@@ -1657,7 +1657,7 @@ subroutine constraint_eqn_035(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fjprot * jwdgpro/jwptf
+      tmp_cc =  jwptf/jwdgpro - 1.0D0 * fjprot
       tmp_con = jwdgpro
       tmp_err =  jwptf - jwdgpro
       tmp_symbol = '<'
@@ -1717,7 +1717,7 @@ subroutine constraint_eqn_037(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fgamcd * gammax/gamcd
+      tmp_cc =  gamcd/gammax - 1.0D0 * fgamcd
       tmp_con = gammax * (1.0D0 - tmp_cc)
       tmp_err = gamcd * tmp_cc
       tmp_symbol = '<'
@@ -1726,6 +1726,7 @@ subroutine constraint_eqn_037(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_037
 
    subroutine constraint_eqn_038(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Obsolete
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! Obsolete
@@ -1771,7 +1772,7 @@ subroutine constraint_eqn_039(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
 
       ! If the temperature peak == 0 then report an error
       if (temp_fw_peak < 1.0D0) call report_error(5)
-      tmp_cc =  1.0D0 - ftpeak * temp_fw_max/temp_fw_peak
+      tmp_cc =  temp_fw_peak/temp_fw_max - 1.0D0 * ftpeak
       tmp_con = temp_fw_max * (1.0D0 - tmp_cc)
       tmp_err = temp_fw_peak * tmp_cc
       tmp_symbol = '<'
@@ -1958,7 +1959,7 @@ subroutine constraint_eqn_044(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
          tcpmax = tcpmax - 273.15D0
       end if
 
-      tmp_cc =   1.0D0 - fptemp * ptempalw / tcpmax
+      tmp_cc =   tcpmax/ptempalw - 1.0D0 * fptemp
       tmp_con = ptempalw * (1.0D0 - tmp_cc)
       tmp_err = tcpmax * tmp_cc
       tmp_symbol = '<'
@@ -2040,7 +2041,7 @@ subroutine constraint_eqn_046(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(10)
       cratmx = 1.0D0 + 4.91D0*(eps-0.62D0)
-      tmp_cc =  1.0D0 - fipir * cratmx * c_tf_total/plasma_current
+      tmp_cc =  cratmx * c_tf_total/plasma_current - 1.0D0 * fipir
       tmp_con = cratmx * (1.0D0 - tmp_cc)
       tmp_err = plasma_current/c_tf_total * tmp_cc
       tmp_symbol = '<'
@@ -2049,6 +2050,7 @@ subroutine constraint_eqn_046(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_046
 
    subroutine constraint_eqn_047(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Issue #508 Remove RFP option: Relevant only to reversed field pinch devices
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! Issue #508 Remove RFP option: Relevant only to reversed field pinch devices
@@ -2093,7 +2095,7 @@ subroutine constraint_eqn_048(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fbeta_poloidal * beta_poloidal_max/beta_poloidal
+      tmp_cc =  beta_poloidal/beta_poloidal_max - 1.0D0 * fbeta_poloidal
       tmp_con = beta_poloidal_max * (1.0D0 - tmp_cc)
       tmp_err = beta_poloidal * tmp_cc
       tmp_symbol = '<'
@@ -2102,6 +2104,7 @@ subroutine constraint_eqn_048(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_048
 
    subroutine constraint_eqn_049(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Issue #508 Remove IFE option: Equation for repetition rate upper limit
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! Issue #508 Remove IFE option: Equation for repetition rate upper limit
@@ -2141,7 +2144,7 @@ subroutine constraint_eqn_050(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
          call report_error(12)
       end if
 
-      tmp_cc =  1.0D0 - frrmax * rrmax/reprat
+      tmp_cc =  reprat/rrmax - 1.0D0 * frrmax
       tmp_con = rrmax * (1.0D0 - tmp_cc)
       tmp_err = reprat * tmp_cc
       tmp_symbol = '<'
@@ -2232,7 +2235,7 @@ subroutine constraint_eqn_053(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fflutf * nflutfmax/nflutf
+      tmp_cc =  nflutf/nflutfmax - 1.0D0 * fflutf
       tmp_con = nflutfmax * (1.0D0 - tmp_cc)
       tmp_err = nflutf * tmp_cc
       tmp_symbol = '<'
@@ -2262,7 +2265,7 @@ subroutine constraint_eqn_054(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0D0 - fptfnuc * ptfnucmax/ptfnucpm3
+      tmp_cc = ptfnucpm3/ptfnucmax - 1.0D0 * fptfnuc
       tmp_con = ptfnucmax * (1.0D0 - tmp_cc)
       tmp_err = ptfnucpm3 * tmp_cc
       tmp_symbol = '<'
@@ -2271,6 +2274,7 @@ subroutine constraint_eqn_054(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_054
 
    subroutine constraint_eqn_055(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! vvhemax is no longer calculated in PROCESS and this constraint is disabled
       implicit none
 
@@ -2306,7 +2310,7 @@ subroutine constraint_eqn_056(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0D0 - fpsepr * pseprmax / (pdivt/rmajor)
+      tmp_cc = (pdivt/rmajor)/pseprmax - 1.0D0 * fpsepr
       tmp_con = pseprmax * (1.0D0 - tmp_cc)
       tmp_err = (pdivt/rmajor) * tmp_cc
       tmp_symbol = '<'
@@ -2315,6 +2319,7 @@ subroutine constraint_eqn_056(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_056
 
    subroutine constraint_eqn_057(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Obsolete
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! Obsolete
@@ -2336,6 +2341,7 @@ subroutine constraint_eqn_057(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_057
 
    subroutine constraint_eqn_058(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Obsolete
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! Obsolete
@@ -2377,7 +2383,7 @@ subroutine constraint_eqn_059(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       real(dp), intent(out) :: tmp_err
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
-      tmp_cc = 1.0D0 - fnbshinef * nbshinefmax / nbshinef
+      tmp_cc = nbshinef/nbshinefmax - 1.0D0 * fnbshinef
       tmp_con = nbshinefmax * (1.0D0 - tmp_cc)
       tmp_err = nbshinef * tmp_cc
       tmp_symbol = '<'
@@ -2416,7 +2422,7 @@ subroutine constraint_eqn_060(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_060
 
    subroutine constraint_eqn_061(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Equation for availability limit
+      !! Equation for availability lower limit
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2498,7 +2504,7 @@ subroutine constraint_eqn_063(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0D0 - fniterpump * n_tf_coils / niterpump
+      tmp_cc = niterpump/n_tf_coils - 1.0D0 * fniterpump
       tmp_con = n_tf_coils
       tmp_err = n_tf_coils * tmp_cc
       tmp_symbol = '<'
@@ -2528,7 +2534,7 @@ subroutine constraint_eqn_064(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0D0 - fzeffmax * (zeffmax/zeff)
+      tmp_cc = zeff/zeffmax - 1.0D0 * fzeffmax
       tmp_con = zeffmax
       tmp_err = zeffmax * tmp_cc
       tmp_symbol = '<'
@@ -2537,7 +2543,7 @@ subroutine constraint_eqn_064(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_064
 
    subroutine constraint_eqn_065(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Limit the stress of the vacuum vessel that occurs when the TF coil quenches.
+      !! Upper limit on stress of the vacuum vessel that occurs when the TF coil quenches.
       !! author: Timothy Nunn, UKAEA
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2555,7 +2561,7 @@ subroutine constraint_eqn_065(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0d0 - fmaxvvstress * max_vv_stress / vv_stress_quench
+      tmp_cc =  vv_stress_quench / max_vv_stress - 1.0d0 * fmaxvvstress
       tmp_con = max_vv_stress
       tmp_err = max_vv_stress * tmp_cc
       tmp_symbol = '<'
@@ -2564,7 +2570,7 @@ subroutine constraint_eqn_065(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_065
 
    subroutine constraint_eqn_066(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Limit on rate of change of energy in poloidal field
+      !! Upper limit on rate of change of energy in poloidal field
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2585,7 +2591,7 @@ subroutine constraint_eqn_066(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0d0 - fpoloidalpower * maxpoloidalpower / peakpoloidalpower
+      tmp_cc = peakpoloidalpower / maxpoloidalpower - 1.0d0 * fpoloidalpower
       tmp_con = maxpoloidalpower
       tmp_err = maxpoloidalpower * tmp_cc
       tmp_symbol = '<'
@@ -2614,7 +2620,7 @@ subroutine constraint_eqn_067(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0d0 - fradwall * pflux_fw_rad_max / pflux_fw_rad_max_mw
+      tmp_cc = pflux_fw_rad_max_mw / pflux_fw_rad_max - 1.0d0 * fradwall
       tmp_con = pflux_fw_rad_max
       tmp_err =  pflux_fw_rad_max * tmp_cc
       tmp_symbol = '<'
@@ -2623,7 +2629,7 @@ subroutine constraint_eqn_067(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_067
 
    subroutine constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! New Psep scaling (PsepB/qAR)
+      !! Upper limit on Psep scaling (PsepB/qAR)
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2652,10 +2658,10 @@ subroutine constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=10), intent(out) :: tmp_units
 
       if (i_q95_fixed == 1) then
-         tmp_cc = 1.0d0 - fpsepbqar * psepbqarmax / ((pdivt*bt)/(q95_fixed*aspect*rmajor))
+         tmp_cc = ((pdivt*bt)/(q95_fixed*aspect*rmajor)) / psepbqarmax - 1.0d0 * fpsepbqar
          tmp_err = (pdivt*bt)/(q95_fixed*aspect*rmajor) - psepbqarmax
       else
-         tmp_cc = 1.0d0 - fpsepbqar * psepbqarmax / ((pdivt*bt)/(q95*aspect*rmajor))
+         tmp_cc = ((pdivt*bt)/(q95*aspect*rmajor)) / psepbqarmax - 1.0d0 * fpsepbqar
          tmp_err = (pdivt*bt)/(q95*aspect*rmajor) - psepbqarmax
       end if
       tmp_con = psepbqarmax
@@ -2665,6 +2671,7 @@ subroutine constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_068
 
    subroutine constraint_eqn_069(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Ensure separatrix power is less than value from Kallenbach divertor
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
@@ -2691,6 +2698,7 @@ subroutine constraint_eqn_069(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_069
 
    subroutine constraint_eqn_070(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Separatrix density consistency
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
@@ -2716,6 +2724,7 @@ subroutine constraint_eqn_070(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_070
 
    subroutine constraint_eqn_071(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! TODO Remove
       !! Separatrix density consistency
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
@@ -2742,7 +2751,7 @@ subroutine constraint_eqn_071(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_071
 
    subroutine constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Central Solenoid Tresca yield criterion
+      !! Upper limit on central Solenoid Tresca yield stress
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2779,11 +2788,11 @@ subroutine constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
 
       ! bucked and wedged desing (see subroutine comment)
       if ( i_tf_bucking >= 2 .and. i_tf_inside_cs == 0 ) then
-         tmp_cc = 1.0d0 - foh_stress * alstroh / max(s_shear_cs_peak, sig_tf_cs_bucked)
+         tmp_cc = max(s_shear_cs_peak, sig_tf_cs_bucked) / alstroh - 1.0d0 * foh_stress
          tmp_err = alstroh - max(s_shear_cs_peak, sig_tf_cs_bucked)
       ! Free standing CS
       else
-         tmp_cc = 1.0d0 - foh_stress * alstroh / s_shear_cs_peak
+         tmp_cc = s_shear_cs_peak / alstroh - 1.0d0 * foh_stress
          tmp_err = alstroh - s_shear_cs_peak
       end if
 
@@ -2794,7 +2803,8 @@ subroutine constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_072
 
    subroutine constraint_eqn_073(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Ensure separatrix power is greater than the L-H power + auxiliary power
+      !! Lower limit to ensure separatrix power is greater than the L-H power + auxiliary power
+      !! Related to constraint 15
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2825,7 +2835,8 @@ subroutine constraint_eqn_073(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_073
 
    subroutine constraint_eqn_074(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Ensure TF coil quench temperature < tmax_croco ONLY used for croco HTS coil
+      !! Upper limit to ensure TF coil quench temperature < tmax_croco
+      !! ONLY used for croco HTS coil
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2846,7 +2857,7 @@ subroutine constraint_eqn_074(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0d0 - fcqt * tmax_croco / croco_quench_temperature
+      tmp_cc = croco_quench_temperature / tmax_croco - 1.0d0 * fcqt
       tmp_con = croco_quench_temperature
       tmp_err = croco_quench_temperature * tmp_cc
       tmp_symbol = '<'
@@ -2855,7 +2866,7 @@ subroutine constraint_eqn_074(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_074
 
    subroutine constraint_eqn_075(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Ensure that TF coil current / copper area < Maximum value
+      !! Upper limit to ensure that TF coil current / copper area < Maximum value
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2876,7 +2887,7 @@ subroutine constraint_eqn_075(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0d0 - f_coppera_m2 * copperA_m2_max / copperA_m2
+      tmp_cc = copperA_m2 / copperA_m2_max - 1.0d0 * f_coppera_m2
       tmp_con = copperA_m2
       tmp_err = copperA_m2 * tmp_cc
       tmp_symbol = '<'
@@ -2885,7 +2896,7 @@ subroutine constraint_eqn_075(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_075
 
    subroutine constraint_eqn_076(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Eich critical separatrix density model: Added for issue 558
+      !! Upper limit for Eich critical separatrix density model: Added for issue 558
       !! author: P B Lloyd, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -2915,7 +2926,7 @@ subroutine constraint_eqn_076(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       nesep_crit = 5.9D0 * alpha_crit * (aspect ** (-2.0D0/7.0D0)) * &
                 (((1.0D0 + (kappa ** 2.0D0)) / 2.0D0) ** (-6.0D0/7.0D0)) &
                 * ((pdivt* 1.0D6) ** (-11.0D0/70.0D0)) * dlimit(7)
-      tmp_cc = 1.0D0 - fnesep * nesep_crit/nesep
+      tmp_cc = nesep / nesep_crit - 1.0D0 * fnesep
       tmp_con = nesep
       tmp_err = nesep * tmp_cc
       tmp_symbol = '<'
@@ -2944,7 +2955,7 @@ subroutine constraint_eqn_077(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0D0 - fcpttf * cpttf_max/cpttf
+      tmp_cc = cpttf / cpttf_max - 1.0D0 * fcpttf
       tmp_con = cpttf_max
       tmp_err = cpttf_max * tmp_cc
       tmp_symbol = '<'
@@ -3010,7 +3021,7 @@ subroutine constraint_eqn_079(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc     = 1.0D0 - fb_cs_limit_max * b_cs_limit_max/max(b_cs_peak_flat_top_end, b_cs_peak_pulse_start)
+      tmp_cc     = max(b_cs_peak_flat_top_end, b_cs_peak_pulse_start) / b_cs_limit_max - 1.0D0 * fb_cs_limit_max
       tmp_con    = b_cs_limit_max
       tmp_err    = max(b_cs_peak_flat_top_end, b_cs_peak_pulse_start) * tmp_cc
       tmp_symbol = '<'
@@ -3048,7 +3059,7 @@ subroutine constraint_eqn_080(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_080
 
    subroutine constraint_eqn_081(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Make sure that the central density is larger that the pedestal one
+      !! Lower limit to ensure central density is larger that the pedestal one
       !! author: S Kahn, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -3164,10 +3175,10 @@ subroutine constraint_eqn_084(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_084
 
    subroutine constraint_eqn_085(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! Equality constraint for the centerpost (CP) lifetime
       !! Author : S Kahn
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
-      !! Equation constraining the centerpost (CP) lifetime
       !! Depending on the chosen option : i_cp_lifetime
       !!  - 0 : The CP full power year lifelime is set by the user (cplife_input)
       !!  - 1 : The CP lifelime is equal to the divertor one
@@ -3220,10 +3231,10 @@ subroutine constraint_eqn_085(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_085
 
    subroutine constraint_eqn_086(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! Upper limit on the turn edge length in the TF winding pack
       !! Author : S Kahn
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units;
-      !! Equation constraining the turn edge length in the TF winding pack
       !! t_turn_tf : input real : TF coil turn edge length including turn insulation [m]
       !! f_t_turn_tf : input real : f-value for TF turn edge length constraint
       !! t_turn_tf_max : input real : TF turn edge length including turn insulation upper limit [m]
@@ -3238,7 +3249,7 @@ subroutine constraint_eqn_086(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=10), intent(out) :: tmp_units
 
       !! Constraints output
-      tmp_cc = 1.0D0 - t_turn_tf / ( f_t_turn_tf * t_turn_tf_max )
+      tmp_cc = t_turn_tf / t_turn_tf_max - 1.0D0 * f_t_turn_tf
       tmp_con = t_turn_tf_max * (1.0D0 - tmp_cc)
       tmp_err = t_turn_tf_max * tmp_cc
       tmp_symbol = '<'
@@ -3248,10 +3259,10 @@ end subroutine constraint_eqn_086
 
 
    subroutine constraint_eqn_087(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! Equation for TF coil cryogenic power upper limit
       !! author: S. Kahn, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
-      !! Equation for TF coil cryogenic power upper limit
       !! crypmw : input real : cryogenic plant power (MW)
       !! f_crypmw : input real : f-value for maximum cryogenic plant power
       !! crypmw_max : input real : Maximum cryogenic plant power (MW)
@@ -3264,7 +3275,7 @@ subroutine constraint_eqn_087(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - f_crypmw * crypmw_max/crypmw
+      tmp_cc =  crypmw / crypmw_max - 1.0D0 * f_crypmw
       tmp_con = crypmw_max * (1.0D0 - tmp_cc)
       tmp_err = crypmw * tmp_cc
       tmp_symbol = '<'
@@ -3293,7 +3304,7 @@ subroutine constraint_eqn_088(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fstr_wp * str_wp_max/abs(str_wp)
+      tmp_cc =  abs(str_wp) / str_wp_max - 1.0D0 * fstr_wp
       tmp_con = str_wp_max
       tmp_err = str_wp_max - abs(str_wp) / fstr_wp
       tmp_symbol = '<'
@@ -3301,7 +3312,7 @@ subroutine constraint_eqn_088(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_088
 
    subroutine constraint_eqn_089(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Ensure that the Central Solenoid [OH] coil current / copper area < Maximum value
+      !! Upper limit to ensure that the Central Solenoid [OH] coil current / copper area < Maximum value
       !! author: G Turkington, CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
@@ -3320,7 +3331,7 @@ subroutine constraint_eqn_089(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc = 1.0d0 - f_copperaoh_m2 * copperaoh_m2_max / copperaoh_m2
+      tmp_cc = copperaoh_m2 / copperaoh_m2_max - 1.0d0 * f_copperaoh_m2
       tmp_con = copperaoh_m2
       tmp_err = copperaoh_m2 * tmp_cc
       tmp_symbol = '<'
@@ -3329,10 +3340,10 @@ subroutine constraint_eqn_089(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_089
 
    subroutine constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! Lower limit for CS coil stress load cycles
       !! author: A. Pearce, G Turkington CCFE, Culham Science Centre
       !! args : output structure : residual error; constraint value;
       !! residual error in physical units; output string; units string
-      !! Equation for minimum CS coil stress load cycles
       !! fncycle : input real : f-value for constraint n_cycle > n_cycle_min
       !! n_cycle : input real : Allowable number of cycles for CS
       !! n_cycle_min : input real :  Minimum required cycles for CS
@@ -3360,7 +3371,7 @@ subroutine constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    end subroutine constraint_eqn_090
 
    subroutine constraint_eqn_091(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
-      !! Equation for checking if the design point is ECRH ignitable
+      !! Lower limit to ensure ECRH te is greater than required te for ignition
       !! at lower values for n and B. Or if the design point is ECRH heatable (if ignite==0)
       !! stellarators only (but in principle usable also for tokamaks).
       !! author: J Lion, IPP Greifswald

From 85d80338ce6fc345172c073ee9f4ee28fb543a51 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Fri, 14 Mar 2025 11:41:51 +0000
Subject: [PATCH 2/7] Correct bounds for fl_h_threshold

Correct initialisation in regression test to agree with new L-H constraint form.
---
 source/fortran/iteration_variables.f90                        | 4 ++--
 .../input_files/spherical_tokamak_once_through.IN.DAT         | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/source/fortran/iteration_variables.f90 b/source/fortran/iteration_variables.f90
index cbdb1ff5b4..30e8d84628 100755
--- a/source/fortran/iteration_variables.f90
+++ b/source/fortran/iteration_variables.f90
@@ -2223,8 +2223,8 @@ subroutine init_itv_103
     use numerics, only: lablxc, boundl, boundu
     implicit none
     lablxc(103) = 'fl_h_threshold     '
-    boundl(103) = 1.000D0
-    boundu(103) = 1.000D6
+    boundl(103) = 0.001D0
+    boundu(103) = 1.000D0
   end subroutine init_itv_103
 
   real(kind(1.d0)) function itv_103()
diff --git a/tests/regression/input_files/spherical_tokamak_once_through.IN.DAT b/tests/regression/input_files/spherical_tokamak_once_through.IN.DAT
index 303f5fc467..fa78c8c74b 100644
--- a/tests/regression/input_files/spherical_tokamak_once_through.IN.DAT
+++ b/tests/regression/input_files/spherical_tokamak_once_through.IN.DAT
@@ -161,7 +161,7 @@ fdene    = 1.0 * f-value for density limit (`constraint equation 5`; `iteration
 ffuspow  = 1.0 * f-value for maximum fusion power (`constraint equation 9`; `iteration variable 26`)
 fiooic   = 0.9942981023320613 * f-value for TF coil operating current / critical current ratio
 fipir    = 0.5460655498423221 * f-value for Ip/Irod upper limit
-fl_h_threshold = 2.903148909029756 * f-value for L-H power threshold (`constraint equation 15`; `iteration variable 103`)
+fl_h_threshold = 0.34445356794812293 * f-value for L-H power threshold (`constraint equation 15`; `iteration variable 103`)
 fpinj    = 0.9546763572954339 * f-value for injection power (`constraint equation 30`; `iteration variable 46`)
 fpsepr   = 1.0 * f-value for maximum Psep/R limit (`constraint equation 56`; `iteration variable 97`)
 fradpwr  = 0.7109311818294267 * f-value for core radiation power limit (`constraint equation 17`; `iteration variable 28`)

From d48645aed37e70ca5ffa98d4afb3a491adbe2f82 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Fri, 14 Mar 2025 14:02:22 +0000
Subject: [PATCH 3/7] Skip non-converging ST regression test

---
 tests/regression/test_process_input_files.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/tests/regression/test_process_input_files.py b/tests/regression/test_process_input_files.py
index 77cf5e2d49..bdbb2d1ef8 100644
--- a/tests/regression/test_process_input_files.py
+++ b/tests/regression/test_process_input_files.py
@@ -269,6 +269,9 @@ def test_input_file(
     should be compared in the test.
     :type opt_params_only: bool
     """
+    if input_file.name == "st_regression.IN.DAT":
+        pytest.skip("ST file doesn't converge with new constraint form currently.")
+
     new_input_file = tmp_path / input_file.name
     shutil.copy(input_file, new_input_file)
 

From 7d1df6c500738255ce079649a7e3807b0771b12a Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Wed, 2 Apr 2025 10:44:41 +0100
Subject: [PATCH 4/7] Correction to form of constraint 12

---
 source/fortran/constraint_equations.f90 | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90
index 877401121f..be4df99c35 100755
--- a/source/fortran/constraint_equations.f90
+++ b/source/fortran/constraint_equations.f90
@@ -885,7 +885,8 @@ subroutine constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=1), intent(out) :: tmp_symbol
       character(len=10), intent(out) :: tmp_units
 
-      tmp_cc =  1.0D0 - fvs * vs_cs_pf_total_pulse/vs_plasma_total_required
+      !! vs_cs_pf_total_pulse is negative, requires sign change
+      tmp_cc =  1.0D0 - fvs * -vs_cs_pf_total_pulse/vs_plasma_total_required
       tmp_con = vs_plasma_total_required * (1.0D0 - tmp_cc)
       tmp_err = vs_plasma_total_required * tmp_cc
       tmp_symbol = '>'

From 5a26d937d22ea1ed00b7e39e95a5b1f2cfd2b514 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Thu, 3 Apr 2025 09:34:28 +0100
Subject: [PATCH 5/7] Correct form of constraint 46

---
 source/fortran/constraint_equations.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90
index be4df99c35..4032c498fc 100755
--- a/source/fortran/constraint_equations.f90
+++ b/source/fortran/constraint_equations.f90
@@ -2042,7 +2042,7 @@ subroutine constraint_eqn_046(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(10)
       cratmx = 1.0D0 + 4.91D0*(eps-0.62D0)
-      tmp_cc =  cratmx * c_tf_total/plasma_current - 1.0D0 * fipir
+      tmp_cc =  (plasma_current / c_tf_total) / cratmx - 1.0D0 * fipir
       tmp_con = cratmx * (1.0D0 - tmp_cc)
       tmp_err = plasma_current/c_tf_total * tmp_cc
       tmp_symbol = '<'

From 13d636b3cd0148fcd1586ab3667a5d6fd5f990b9 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Thu, 3 Apr 2025 09:36:56 +0100
Subject: [PATCH 6/7] Revert "Skip non-converging ST regression test"

This reverts commit d48645aed37e70ca5ffa98d4afb3a491adbe2f82.
---
 tests/regression/test_process_input_files.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/tests/regression/test_process_input_files.py b/tests/regression/test_process_input_files.py
index bdbb2d1ef8..77cf5e2d49 100644
--- a/tests/regression/test_process_input_files.py
+++ b/tests/regression/test_process_input_files.py
@@ -269,9 +269,6 @@ def test_input_file(
     should be compared in the test.
     :type opt_params_only: bool
     """
-    if input_file.name == "st_regression.IN.DAT":
-        pytest.skip("ST file doesn't converge with new constraint form currently.")
-
     new_input_file = tmp_path / input_file.name
     shutil.copy(input_file, new_input_file)
 

From 10c4b3f036e2e182e260734fdfd40260ca573d01 Mon Sep 17 00:00:00 2001
From: Jonathan Maddock <78556175+jonmaddock@users.noreply.github.com>
Date: Thu, 3 Apr 2025 10:22:40 +0100
Subject: [PATCH 7/7] Enclose negative sign in parentheses

---
 source/fortran/constraint_equations.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90
index 4032c498fc..3858412ab6 100755
--- a/source/fortran/constraint_equations.f90
+++ b/source/fortran/constraint_equations.f90
@@ -886,7 +886,7 @@ subroutine constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       character(len=10), intent(out) :: tmp_units
 
       !! vs_cs_pf_total_pulse is negative, requires sign change
-      tmp_cc =  1.0D0 - fvs * -vs_cs_pf_total_pulse/vs_plasma_total_required
+      tmp_cc =  1.0D0 - fvs * (-vs_cs_pf_total_pulse)/vs_plasma_total_required
       tmp_con = vs_plasma_total_required * (1.0D0 - tmp_cc)
       tmp_err = vs_plasma_total_required * tmp_cc
       tmp_symbol = '>'