constraint_equations.f90 Source File


Contents


Source Code

! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module constraints
  !! author: J Morris
  !!
  !! Module defining the constraint equations and the routine that evaluates the
  !! constraint equations.

  ! Import modules
#ifndef dp
  use, intrinsic :: iso_fortran_env, only: dp=>real64
#endif
  use error_handling, only: report_error, idiags, fdiags

  implicit none

  public :: constraint_eqns

!   type constraint_args_type
!     real(dp) :: cc
!     !! Residual error in constraint equation
!     real(dp) :: con
!     !! constraint value for constraint equation in physical units
!     real(dp) :: err
!     !! residual error in constraint equation in physical units
!     character(len=1)  :: symbol
!     !! `=<`, `>`, `<` symbol for constraint equation denoting its type
!     character(len=10) :: units
!     !! constraint units in constraint equation
!   end type

contains

   subroutine constraint_eqns(m,ieqn,cc,con,err,symbol,units)
    !! Routine that formulates the constraint equations
    !!
    !! **author: P J Knight** (UKAEA)
    !!
    !! **author: J Morris** (UKAEA)
    !!
    !! if `ieqn` is zero or negative, evaluate all the constraint equations, otherwise
    !! evaluate only the `ieqn`th equation. The code attempts to make \( cc(i) = 0 \) for all
    !! \( i \in \{1,\cdots,m\} \) equations. All relevant consistency equations should be active in
    !! order to make a self-consistent machine.
    !!
    !! **References**
    !!
    !! 1.
    use numerics, only: icc
    use maths_library, only: variable_error

    implicit none

    integer, intent(in) :: m
    !! Number of constraint equations to solve

    integer, intent(in) :: ieqn
    !! Switch for constraint equations to evaluate;

    real(dp), dimension(m), intent(out) :: cc
    !! Residual error in equation i

    real(dp), optional, dimension(m), intent(out) :: con
    !! constraint value for equation i in physical units

    real(dp), optional, dimension(m), intent(out) :: err
    !! residual error in equation i in physical units

    character(len=1),  optional, dimension(m), intent(out) :: symbol
    !! `=<`, `>`, `<` symbol for equation i denoting its type

    character*10, optional, dimension(m), intent(out) :: units
    !! constraint units in equation i

    ! Local variables
    integer :: i, i1, i2

    real(dp) :: tmp_cc = 0
    !! Residual error in constraint equation
    real(dp) :: tmp_con = 0
    !! constraint value for constraint equation in physical units
    real(dp) :: tmp_err = 0
    !! residual error in constraint equation in physical units
    character(len=1)  :: tmp_symbol = ' '
    !! `=<`, `>`, `<` symbol for constraint equation denoting its type
    character(len=10) :: tmp_units = ' '
    !! constraint units in constraint equation


    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! If ieqn is positive, only evaluate the 'ieqn'th constraint residue,
    ! otherwise evaluate all m constraint residues
    if (ieqn > 0) then
      i1 = ieqn ; i2 = ieqn
    else
      i1 = 1 ; i2 = m
    end if

    ! Consistency (equality) constraints should converge to zero.
    do i = i1,i2

      ! The constraint value in physical units is
      ! a) for consistency equations, the quantity to be equated, or
      ! b) for limit equations, the limiting value.
      ! The symbol is = for a consistency equation, < for an upper limit
      ! or > for a lower limit.
      select case (icc(i))

	      ! Relationship between beta, temperature (keV) and density
        case (1); call constraint_eqn_001(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Global plasma power balance equation
        case (2); call constraint_eqn_002(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Global power balance equation for ions
        case (3); call constraint_eqn_003(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Global power balance equation for electrons
        case (4); call constraint_eqn_004(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for density upper limit
        case (5); call constraint_eqn_005(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for epsilon beta-poloidal upper limit
        case (6); call constraint_eqn_006(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for hot beam ion density
        case (7); call constraint_eqn_007(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for neutron wall load upper limit
        case (8); call constraint_eqn_008(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for fusion power upper limit
        case (9); call constraint_eqn_009(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Obsolete
        case (10); call constraint_eqn_010(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for radial build
        case (11); call constraint_eqn_011(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for volt-second capability lower limit
        case (12); call constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for burn time lower limit
        case (13); call constraint_eqn_013(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation to fix number of NBI decay lengths to plasma centre
        case (14); call constraint_eqn_014(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for L-H power threshold limit
        case (15); call constraint_eqn_015(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for net electric power lower limit
        case (16); call constraint_eqn_016(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for radiation power upper limit
        case (17); call constraint_eqn_017(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for divertor heat load upper limit
        case (18); call constraint_eqn_018(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for MVA upper limit
        case (19); call constraint_eqn_019(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for neutral beam tangency radius upper limit
        case (20); call constraint_eqn_020(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for minor radius lower limit
        case (21); call constraint_eqn_021(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for divertor collision/connection length ratio upper limit
        case (22); call constraint_eqn_022(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for conducting shell radius / rminor upper limit
        case (23); call constraint_eqn_023(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for beta upper limit
        case (24); call constraint_eqn_024(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for peak toroidal field upper limit
        case (25); call constraint_eqn_025(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for Central Solenoid current density upper limit at EOF
        case (26); call constraint_eqn_026(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for Central Solenoid current density upper limit at BOP
        case (27); call constraint_eqn_027(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for fusion gain (big Q) lower limit
        case (28); call constraint_eqn_028(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for inboard major radius
        case (29); call constraint_eqn_029(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for injection power upper limit
        case (30); call constraint_eqn_030(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for TF coil case stress upper limit (SCTF)
        case (31); call constraint_eqn_031(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for TF coil conduit stress upper limit (SCTF)
        case (32); call constraint_eqn_032(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for TF coil operating/critical J upper limit (SCTF)
        case (33); call constraint_eqn_033(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for TF coil dump voltage upper limit (SCTF)
        case (34); call constraint_eqn_034(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for TF coil J_wp/J_prot upper limit (SCTF)
        case (35); call constraint_eqn_035(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for TF coil s/c temperature margin lower limit (SCTF)
        case (36); call constraint_eqn_036(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for current drive gamma upper limit
        case (37); call constraint_eqn_037(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for first wall temperature upper limit
        case (39); call constraint_eqn_039(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for auxiliary power lower limit
        case (40); call constraint_eqn_040(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for plasma current ramp-up time lower limit
        case (41); call constraint_eqn_041(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for cycle time lower limit
        case (42); call constraint_eqn_042(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for average centrepost temperature
        case (43); call constraint_eqn_043(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for centrepost temperature upper limit (TART)
        case (44); call constraint_eqn_044(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for edge safety factor lower limit (TART)
        case (45); call constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for Ip/Irod upper limit (TART)
        case (46); call constraint_eqn_046(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for TF coil toroidal thickness upper limit
        case (47); call constraint_eqn_047(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for poloidal beta upper limit
        case (48); call constraint_eqn_048(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Issue #508 Remove RFP option Equation to ensure reversal parameter F is negative
        case (49); call constraint_eqn_049(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! IFE option: Equation for repetition rate upper limit
        case (50); call constraint_eqn_050(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation to enforce startup flux = available startup flux
        case (51); call constraint_eqn_051(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for tritium breeding ratio lower limit
        case (52); call constraint_eqn_052(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for fast neutron fluence on TF coil upper limit
        case (53); call constraint_eqn_053(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for peak TF coil nuclear heating upper limit
        case (54); call constraint_eqn_054(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for helium concentration in vacuum vessel upper limit
        case (55); call constraint_eqn_055(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for power through separatrix / major radius upper limit
        case (56); call constraint_eqn_056(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Obsolete
        case (57); call constraint_eqn_057(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Obsolete
        case (58); call constraint_eqn_058(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for neutral beam shine-through fraction upper limit
        case (59); call constraint_eqn_059(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for Central Solenoid s/c temperature margin lower limit
        case (60); call constraint_eqn_060(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Equation for availability limit
        case (61); call constraint_eqn_061(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Lower limit on taup/taueff the ratio of alpha particle to energy confinement times
        case (62); call constraint_eqn_062(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Upper limit on niterpump (vacuum_model = simple)
        case (63); call constraint_eqn_063(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Upper limit on Zeff
        case (64); call constraint_eqn_064(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Limit TF dump time to calculated quench time
        case (65); call constraint_eqn_065(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Limit on rate of change of energy in poloidal field
        case (66); call constraint_eqn_066(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Simple upper limit on radiation wall load
        case (67); call constraint_eqn_067(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! New Psep scaling (PsepB/qAR)
        case (68); call constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Ensure separatrix power is less than value from Kallenbach divertor
        case (69); call constraint_eqn_069(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Separatrix temperature consistency
        case (70); call constraint_eqn_070(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Separatrix density consistency
        case (71); call constraint_eqn_071(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	      ! Central Solenoid Tresca yield criterion
        case (72); call constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! ensure separatrix power is greater than the L-H power + auxiliary power
        case (73); call constraint_eqn_073(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! ensure TF coil quench temperature < tmax_croco
        case (74); call constraint_eqn_074(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	     ! ensure that TF coil current / copper area < Maximum value ONLY used for croco HTS coil
        case (75); call constraint_eqn_075(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Eich critical separatrix density model
        case (76); call constraint_eqn_076(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Equation for maximum TF current per turn upper limit
        case (77); call constraint_eqn_077(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
	  	  ! Equation for Reinke criterion, divertor impurity fraction lower limit
        case (78); call constraint_eqn_078(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Equation for maximum CS field
        case (79); call constraint_eqn_079(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Lower limit pdivt
        case (80); call constraint_eqn_080(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Constraint equation making sure that ne(0) > ne(ped)
        case (81); call constraint_eqn_081(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Constraint equation making sure that stellarator coils dont touch in toroidal direction
        case (82); call constraint_eqn_082(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Constraint ensuring radial build consistency for stellarators
        case (83); call constraint_eqn_083(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        ! Constraint for lower limit of beta
        case (84); call constraint_eqn_084(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
         ! Constraint for CP lifetime
        case (85); call constraint_eqn_085(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
         ! Constraint for turn dimension
        case (86); call constraint_eqn_086(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
         ! Constraint for cryogenic power
        case (87); call constraint_eqn_087(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
         ! Constraint for TF coil strain
        case (88); call constraint_eqn_088(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
         ! ensure that OH coil current / copper area < Maximum value ONLY used for croco HTS coil
        case (89); call constraint_eqn_089(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
         ! Constraint for minimum CS stress load cycles
        case (90); call constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
         ! Constraint for indication of ECRH ignitability
        case (91); call constraint_eqn_091(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       case default

        idiags(1) = icc(i)
        call report_error(13)
        tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ' '
        tmp_units = ' '

      end select

      cc(i) = tmp_cc
      if (present(con)) then
        con(i) = tmp_con
        if (present(err))    err(i)    = tmp_err
        if (present(symbol)) symbol(i) = tmp_symbol
        if (present(units))  units(i)  = tmp_units
      end if

      ! Crude method of catching NaN errors
      !if ((abs(cc(i)) > 9.99D99).or.(cc(i) /= cc(i))) then
      if (variable_error(cc(i))) then

        ! Add debugging lines as appropriate...
        select case (icc(i))

          ! Relationship between beta, temperature (keV) and density (consistency equation)
          case (1); call constraint_err_001()
          ! Equation for net electric power lower limit
          case (16); call constraint_err_016()
          ! Equation for injection power upper limit
          case (30); call constraint_err_030()
          ! Limit on rate of change of energy in poloidal field
          case (66); call constraint_err_066()

        end select

        idiags(1) = icc(i) ; fdiags(1) = cc(i)
        call report_error(14)

      end if

    end do
    ! Issue 505 Reverse the sign so it works as an inequality constraint (cc(i) > 0)
    ! This will have no effect if it is used as an equality constraint because it will be squared.
    cc = -cc

   end subroutine constraint_eqns

   !--- Error-handling routines

   subroutine constraint_err_001()
    !! Error in: Relationship between beta, temperature (keV) and density (consistency equation)
    !! author: P B Lloyd, CCFE, Culham Science Centre
    use physics_variables, only: betaft, beta_beam, dene, ten, dnitot, tin, btot, beta
    write(*,*) 'betaft = ', betaft
    write(*,*) 'beta_beam = ', beta_beam
    write(*,*) 'dene = ', dene
    write(*,*) 'ten = ', ten
    write(*,*) 'dnitot = ', dnitot
    write(*,*) 'tin = ', tin
    write(*,*) 'btot = ',btot
    write(*,*) 'beta = ', beta
   end subroutine

   subroutine constraint_err_016()
    !! Error in: Equation for net electric power lower limit
    !! author: P B Lloyd, CCFE, Culham Science Centre
    use constraint_variables, only: fpnetel, pnetelin
    use heat_transport_variables, only: pnetelmw
    implicit none
    write(*,*) 'fpnetel = ', fpnetel
    write(*,*) 'pnetelmw = ', pnetelmw
    write(*,*) 'pnetelin = ', pnetelin
   end subroutine

   subroutine constraint_err_030()
    !! Error in: Equation for injection power upper limit
    !! author: P B Lloyd, CCFE, Culham Science Centre
    use current_drive_variables, only: pinjmw, pinjalw
    use constraint_variables, only: fpinj
    implicit none
    write(*,*) 'fpinj = ', fpinj
    write(*,*) 'pinjalw = ', pinjalw
    write(*,*) 'pinjmw = ', pinjmw
   end subroutine

   subroutine constraint_err_066()
    !! Error in: Limit on rate of change of energy in poloidal field
    !! author: P B Lloyd, CCFE, Culham Science Centre
    use constraint_variables, only: fpoloidalpower
    use pf_power_variables, only: maxpoloidalpower, peakpoloidalpower
    implicit none
    write(*,*) 'fpoloidalpower = ', fpoloidalpower
    write(*,*) 'maxpoloidalpower = ', maxpoloidalpower
    write(*,*) 'peakpoloidalpower = ', peakpoloidalpower
   end subroutine constraint_err_066

   !---

   subroutine constraint_eqn_001(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    !! author: J Morris
    !! category: equality constraint
    !!
    !! Relationship between beta, temperature (keV) and density
    !!
    !! \begin{equation}
    !! c_i = 1 - \frac{1}{\beta}\left( \beta_{ft} + \beta_{NBI} + 2 \times 10^3 \mu_0 e
    !! \left( \frac{n_e T_e + n_i T_i}{B_{tot}^2} \right) \right)
    !! \end{equation}
    !!
    !! - \( \beta \) -- total plasma beta
    !! - \( \beta_{ft} \) -- fast alpha beta component
    !! - \( \beta_{NBI} \) -- neutral beam beta component
    !! - \( n_e \) -- electron density [m\(^3\)]
    !! - \( n_i \) -- total ion density [m\(^3\)]
    !! - \( T_e \) -- density weighted average electron temperature [keV]
    !! - \( T_i \) -- density weighted average ion temperature [keV]
    !! - \( B_{tot} \) -- total toroidal + poloidal field [T]

    use physics_variables, only: betaft, beta_beam, dene, ten, dnitot, tin, btot, beta
    use constants, only: electron_charge,rmu0

    implicit none

   !  type(constraint_args_type), intent(out) :: args
      real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

    !! constraint derived type

      tmp_cc = 1.0D0 - (betaft + beta_beam + &
        2.0D3*rmu0*electron_charge * (dene*ten + dnitot*tin)/btot**2 )/beta
      tmp_con = beta * (1.0D0 - tmp_cc)
      tmp_err = beta * tmp_cc
      tmp_symbol = '='
      tmp_units  = ''

   end subroutine constraint_eqn_001

   subroutine constraint_eqn_002(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
    !! author: J. Morris
    !! category: equality constraint
    !!
    !! Global plasma power balance equation
    !!
    !! \begin{equation} c_i =
    !! \end{equation}
    !!
    !! iradloss : input integer : switch for radiation loss term usage in power balance (see User Guide):<UL>
    !! <LI> = 0 total power lost is scaling power plus radiation (needed for ipedestal=2,3)
    !! <LI> = 1 total power lost is scaling power plus core radiation only
    !! <LI> = 2 total power lost is scaling power only, with no additional
    !! allowance for radiation. This is not recommended for power plant models.</UL>
    !! ignite : input integer : switch for ignition assumption:<UL>
    !! <LI> = 0 do not assume plasma ignition;
    !! <LI> = 1 assume ignited (but include auxiliary power in costs)</UL>
    !! ptrepv : input real : electron transport power per volume (MW/m3)
    !! ptripv : input real :  ion transport power per volume (MW/m3)
    !! pradpv : input real : total radiation power per volume (MW/m3)
    !! pcoreradpv : input real : total core radiation power per volume (MW/m3)
    !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma
    !! alpha_power_density_total : input real : alpha power per volume (MW/m3)
    !! charged_power_density : input real : non-alpha charged particle fusion power per volume (MW/m3)
    !! pohmpv : input real : ohmic heating power per volume (MW/m3)
    !! pinjmw : input real : total auxiliary injected power (MW)
    !! plasma_volume : input real : plasma volume (m3)

    use physics_variables, only: iradloss, ignite, ptrepv, ptripv, pradpv, &
                                  pcoreradpv, f_alpha_plasma, alpha_power_density_total, charged_power_density, &
                                  pohmpv, plasma_volume
    use current_drive_variables, only: pinjmw

    implicit none

          real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
    !! constraint derived type

    ! pscaling : Local real : total transport power per volume (MW/m3)
    real(dp) :: pscaling
    real(dp) :: pnumerator, pdenom
    pscaling = ptrepv + ptripv
    ! Total power lost is scaling power plus radiation:
    if (iradloss == 0) then
        pnumerator = pscaling + pradpv
    else if (iradloss == 1) then
        pnumerator = pscaling + pcoreradpv
    else
        pnumerator = pscaling
    end if

    ! if plasma not ignited include injected power
    if (ignite == 0) then
      pdenom = f_alpha_plasma*alpha_power_density_total + charged_power_density + pohmpv + pinjmw/plasma_volume
    else
      ! if plasma ignited
      pdenom = f_alpha_plasma*alpha_power_density_total + charged_power_density + pohmpv
    end if

    tmp_cc = 1.0D0 - pnumerator / pdenom
    tmp_con = pdenom * (1.0D0 - tmp_cc)
    tmp_err = pdenom * tmp_cc
    tmp_symbol = '='
    tmp_units = 'MW/m3'

   end subroutine constraint_eqn_002

   subroutine constraint_eqn_003(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Global power balance equation for ions
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Global power balance equation for ions
      !! This is a consistency equation (NBI)
      !! #=# physics
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ignite : input integer : switch for ignition assumption:<UL>
      !! <LI> = 0 do not assume plasma ignition;
      !! <LI> = 1 assume ignited (but include auxiliary power in costs)</UL>
      !! ptripv : input real :  ion transport power per volume (MW/m3)
      !! piepv : input real : ion/electron equilibration power per volume (MW/m3)
      !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma
      !! alpha_power_ions_density : input real : alpha power per volume to ions (MW/m3)
      !! pinjimw : input real : auxiliary injected power to ions (MW)
      !! plasma_volume : input real : plasma volume (m3)
      use physics_variables, only: ignite, ptripv, piepv, f_alpha_plasma, alpha_power_ions_density, plasma_volume
      use current_drive_variables, only: pinjimw
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

	   ! No assume plasma ignition:
      if (ignite == 0) then
         tmp_cc     = 1.0D0 - (ptripv + piepv) / (f_alpha_plasma*alpha_power_ions_density + pinjimw/plasma_volume)
         tmp_con    = (f_alpha_plasma*alpha_power_ions_density + pinjimw/plasma_volume) * (1.0D0 - tmp_cc)
         tmp_err    = (f_alpha_plasma*alpha_power_ions_density + pinjimw/plasma_volume) * tmp_cc
         tmp_symbol = '='
         tmp_units  = 'MW/m3'
	   ! Plasma ignited:
      else
         tmp_cc     = 1.0D0 - (ptripv+piepv) / (f_alpha_plasma*alpha_power_ions_density)
         tmp_con    = (f_alpha_plasma*alpha_power_ions_density) * (1.0D0 - tmp_cc)
         tmp_err    = (f_alpha_plasma*alpha_power_ions_density) * tmp_cc
         tmp_symbol = '='
         tmp_units  = 'MW/m3'
      end if

   end subroutine constraint_eqn_003

   subroutine constraint_eqn_004(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Global power balance equation for electrons
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Global power balance equation for electrons
      !! This is a consistency equation
      !! N.B. This constraint is currently NOT RECOMMENDED for use.
      !! #=# physics
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! iradloss : input integer : switch for radiation loss term usage in power balance (see User Guide):<UL>
      !! <LI> = 0 total power lost is scaling power plus radiation (needed for ipedestal=2,3)
      !! <LI> = 1 total power lost is scaling power plus core radiation only
      !! <LI> = 2 total power lost is scaling power only, with no additional
      !! allowance for radiation. This is not recommended for power plant models.</UL>
      !! ignite : input integer : switch for ignition assumption:<UL>
      !! <LI> = 0 do not assume plasma ignition;
      !! <LI> = 1 assume ignited (but include auxiliary power in costs)</UL>
      !! ptrepv : input real : electron transport power per volume (MW/m3)
      !! pradpv : input real : total radiation power per volume (MW/m3)
      !! pcoreradpv : input real : total core radiation power per volume (MW/m3)
      !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma
      !! alpha_power_electron_density : input real : alpha power per volume to electrons (MW/m3)
      !! piepv : input real : ion/electron equilibration power per volume (MW/m3)
      !! pinjemw : input real : auxiliary injected power to electrons (MW)
      !! plasma_volume : input real : plasma volume (m3)
      use physics_variables, only: iradloss, ignite, ptrepv, pcoreradpv, f_alpha_plasma, &
                                 alpha_power_electron_density, piepv, plasma_volume, pradpv
      use current_drive_variables, only: pinjemw
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! pscaling : Local real : total transport power per volume (MW/m3)
      real(dp) :: pscaling
      real(dp) :: pnumerator, pdenom
      pscaling = ptrepv
	   ! Total power lost is scaling power plus radiation:
      if (iradloss == 0) then
         pnumerator = pscaling + pradpv
      else if (iradloss == 1) then
         pnumerator = pscaling + pcoreradpv
      else
         pnumerator = pscaling
      end if

      ! if plasma not ignited include injected power
      if (ignite == 0) then
         pdenom = f_alpha_plasma*alpha_power_electron_density + piepv + pinjemw/plasma_volume
      else
      ! if plasma ignited
         pdenom = f_alpha_plasma*alpha_power_electron_density + piepv
      end if

      tmp_cc     = 1.0D0 - pnumerator / pdenom
      tmp_con    = pdenom * (1.0D0 - tmp_cc)
      tmp_err    = pdenom * tmp_cc
      tmp_symbol = '='
      tmp_units  = 'MW/m3'

   end subroutine constraint_eqn_004

   subroutine constraint_eqn_005(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for density upper 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
      !! Equation for density upper limit
      !! #=# physics
      !! #=#=# fdene, dnelimt
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! idensl : input integer : switch for density limit to enforce (constraint equation 5):<UL>
      !! <LI> = 1 old ASDEX;
      !! <LI> = 2 Borrass model for ITER (I);
      !! <LI> = 3 Borrass model for ITER (II);
      !! <LI> = 4 JET edge radiation;
      !! <LI> = 5 JET simplified;
      !! <LI> = 6 Hugill-Murakami Mq limit;
      !! <LI> = 7 Greenwald limit</UL>
      !! fdene : input real : f-value for density limit
      !! dene : input real : electron density (/m3)
      !! dnelimt : input real : density limit (/m3)
      !! dnla : input real : line averaged electron density (m-3)
      use physics_variables, only: idensl, dnelimt, dnla, dene
      use constraint_variables, only: fdene
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

	   ! Apply Greenwald limit to line-averaged density
      if (idensl == 7) then
         tmp_cc     = 1.0D0 - fdene * dnelimt/dnla
         tmp_con    = fdene * dnelimt
         tmp_err    = fdene * dnelimt - dnla
         tmp_symbol = '<'
         tmp_units  = '/m3'
      else
         tmp_cc = 1.0D0 - fdene * dnelimt/dene
         tmp_con    = dnelimt * (1.0D0 - tmp_cc)
         tmp_err    = dene * tmp_cc
         tmp_symbol = '<'
         tmp_units  = '/m3'
      end if

   end subroutine constraint_eqn_005

   subroutine constraint_eqn_006(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for epsilon beta-poloidal upper 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
      !! Equation for epsilon beta-poloidal upper limit
      !! #=# physics
      !! #=#=# fbeta, epbetmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fbeta : input real : f-value for epsilon beta-poloidal
      !! epbetmax : input real : maximum (eps*beta_poloidal)
      !! eps : input real : inverse aspect ratio
      !! betap : input real : poloidal beta
      use physics_variables, only: epbetmax, eps, betap
      use constraint_variables, only: fbeta, fbeta
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fbeta * epbetmax/(eps*betap)
      tmp_con = epbetmax * (1.0D0 - tmp_cc)
      tmp_err = (eps*betap) * tmp_cc
      tmp_symbol = '<'
      tmp_units = ''

   end subroutine constraint_eqn_006

   subroutine constraint_eqn_007(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for hot beam ion density
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for hot beam ion density
      !! This is a consistency equation (NBI)
      !! #=# physics
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ignite : input integer : switch for ignition assumption:<UL>
      !! <LI> = 0 do not assume plasma ignition;
      !! <LI> = 1 assume ignited (but include auxiliary power in costs)</UL>
      !! Obviously, ignite must be zero if current drive is required.
      !! If ignite=1, any auxiliary power is assumed to be used only
      !! during plasma start-up, and is excluded from all steady-state
      !! power balance calculations.
      !! beam_density_out : input real :  hot beam ion density from calculation (/m3)
      !! dnbeam : input real : hot beam ion density, variable (/m3)
      use physics_variables, only: ignite, beam_density_out, dnbeam
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

	   ! Do not assume plasma ignition:
      if (ignite == 0) then
         tmp_cc     = 1.0D0 - beam_density_out/dnbeam
         tmp_con    = dnbeam * (1.0D0 - tmp_cc)
         tmp_err    = dnbeam * tmp_cc
         tmp_symbol = '='
         tmp_units  = '/m3'
      else
         tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ''
        tmp_units = ''
         call report_error(1)
      end if

   end subroutine constraint_eqn_007

   subroutine constraint_eqn_008(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for neutron wall load upper 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
      !! Equation for neutron wall load upper limit
      !! #=# physics
      !! #=#=# fwalld, walalw
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fwalld : input real : f-value for maximum wall load
      !! walalw : input real : allowable wall-load (MW/m2)
      !! wallmw : input real : average neutron wall load (MW/m2)
      use constraint_variables, only: fwalld, walalw
      use physics_variables, only: wallmw
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fwalld * walalw/wallmw
      tmp_con = fwalld * walalw
      tmp_err = fwalld * walalw - wallmw
      tmp_symbol = '<'
      tmp_units = 'MW/m2'

   end subroutine constraint_eqn_008

   subroutine constraint_eqn_009(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for fusion power upper 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
      !! Equation for fusion power upper limit
      !! #=# physics
      !! #=#=# ffuspow, powfmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ffuspow : input real : f-value for maximum fusion power
      !! powfmax : input real : maximum fusion power (MW)
      !! fusion_power : input real : fusion power (MW)
      use constraint_variables, only: ffuspow, powfmax
      use physics_variables, only: fusion_power
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - ffuspow * powfmax/fusion_power
      tmp_con = powfmax * (1.0D0 - tmp_cc)
      tmp_err = fusion_power * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW'

   end subroutine constraint_eqn_009

   subroutine constraint_eqn_010(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! Equation for field at TF coil
      !! This is a consistency equation
      !! (do not use for stellarators)
      !! #=# tfcoil
      !! #=#=# consistency
      !! rmajor |  plasma major radius (m)
      !! bt     |  toroidal field on axis (T)
      !! rbmax  |  radius of maximum toroidal field (m)
      !! bmaxtf |  peak field at toroidal field coil (T)

      !! This constraint is depreciated

      implicit none

            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
      !! Constraints output

      ! This constraint is depreciated
      call report_error(236)

      tmp_con = 1.0D0
      tmp_err = 0.0D0
      tmp_symbol = '='
      tmp_units = ''

   end subroutine constraint_eqn_010

   subroutine constraint_eqn_011(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for radial build
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for radial build
      !! (This is a consistency equation.)
      !! #=# build
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! rbld : input real : sum of thicknesses to the major radius (m)
      !! rmajor : input real : plasma major radius (m)
      use build_variables, only: rbld
      use physics_variables, only: rmajor
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - rbld/rmajor
      tmp_con = rmajor * (1.0D0 - tmp_cc)
      tmp_err = rmajor * tmp_cc
      tmp_symbol = '='
      tmp_units = 'm'

   end subroutine constraint_eqn_011

   subroutine constraint_eqn_012(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for volt-second capability 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
      !! Equation for volt-second capability lower limit
      !! #=# pfcoil
      !! #=#=# fvs, vsstt
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! vsstt : input real : total V-s needed (Wb)
      !! vsstt (lower limit) is positive; vstot (available) is negative
      !! fvs : input real : f-value for flux-swing (V-s) requirement (STEADY STATE)
      !! vstot : input real :   total flux swing for pulse (Wb)
      use physics_variables, only: vsstt
      use constraint_variables, only: fvs
      use pfcoil_variables, only: vstot
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 + fvs * vstot/vsstt
      tmp_con = vsstt * (1.0D0 - tmp_cc)
      tmp_err = vsstt * tmp_cc
      tmp_symbol = '>'
      tmp_units = 'V.sec'

   end subroutine constraint_eqn_012

   subroutine constraint_eqn_013(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for burn time 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
      !! Equation for burn time lower limit
      !! #=# times
      !! #=#=# ftburn, tbrnmn
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftburn : input real : f-value for minimum burn time
      !! tburn : input real : burn time (s) (calculated if lpulse=1)
      !! tbrnmn : input real :  minimum burn time (s)
      use constraint_variables, only: ftburn,tbrnmn
      use times_variables, only: tburn
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - ftburn * tburn/tbrnmn
      tmp_con = tbrnmn / ftburn
      tmp_err = tbrnmn / ftburn  - tburn
      tmp_symbol = '>'
      tmp_units = 'sec'

   end subroutine constraint_eqn_013

   subroutine constraint_eqn_014(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation to fix number of NBI decay lengths to plasma centre
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation to fix number of NBI decay lengths to plasma centre
      !! This is a consistency equation
      !! #=# current_drive
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! taubeam : input real : neutral beam e-decay lengths to plasma centre
      !! tbeamin : input real : permitted neutral beam e-decay lengths to plasma centre
      use current_drive_variables, only: taubeam, tbeamin
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - taubeam/tbeamin
      tmp_con = tbeamin * (1.0D0 - tmp_cc)
      tmp_err = tbeamin * tmp_cc
      tmp_symbol = '='
      tmp_units = ''

   end subroutine constraint_eqn_014

   subroutine constraint_eqn_015(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for L-H power threshold 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
      !! Equation for L-H power threshold limit
      !! #=# physics
      !! #=#=# flhthresh, plhthresh
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! flhthresh : input real : f-value for L-H power threshold
      !! plhthresh : input real : L-H mode power threshold (MW)
      !! pdivt : input real : power to conducted to the divertor region (MW)
      use constraint_variables, only: flhthresh
      use physics_variables, only: plhthresh, pdivt
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  -(1.0D0 - flhthresh * plhthresh / pdivt)
      tmp_con = plhthresh
      tmp_err = plhthresh - pdivt / flhthresh
      if (flhthresh > 1.0D0) then
         tmp_symbol = '>'
      else
         tmp_symbol = '<'
      end if
      tmp_units = 'MW'

   end subroutine constraint_eqn_015

   subroutine constraint_eqn_016(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for net electric power 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
      !! Equation for net electric power lower limit
      !! #=# heat_transport
      !! #=#=# fpnetel, pnetelin
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fpnetel : input real : f-value for net electric power
      !! pnetelmw : input real : net electric power (MW)
      !! pnetelin : input real : required net electric power (MW)
      use constraint_variables, only: fpnetel, pnetelin
      use heat_transport_variables, only: pnetelmw
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fpnetel * pnetelmw / pnetelin
      tmp_con = pnetelin
      tmp_err = pnetelmw - pnetelin / fpnetel
      tmp_symbol = '>'
      tmp_units = 'MW'

   end subroutine constraint_eqn_016

   subroutine constraint_eqn_017(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for radiation power upper 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
      !! Equation for radiation power upper limit
      !! #=# physics
      !! #=#=# fradpwr, pradmaxpv
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! f_alpha_plasma : input real : fraction of alpha power deposited in plasma
      !! pinjmw : input real : total auxiliary injected power (MW)
      !! plasma_volume : input real : plasma volume (m3)
      !! alpha_power_density_total : input real : alpha power per volume (MW/m3)
      !! charged_power_density :  input real : non-alpha charged particle fusion power per volume (MW/m3)
      !! pohmpv : input real : ohmic heating power per volume (MW/m3)
      !! fradpwr : input real : f-value for core radiation power limit
      !! pradpv : input real : total radiation power per volume (MW/m3)
      use physics_variables, only: f_alpha_plasma, plasma_volume, alpha_power_density_total, charged_power_density, pohmpv, pradpv
      use current_drive_variables, only: pinjmw
      use constraint_variables, only: fradpwr
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      real(dp) :: pradmaxpv
      !! Maximum possible power/plasma_volume that can be radiated (local)

      pradmaxpv = pinjmw/plasma_volume + alpha_power_density_total*f_alpha_plasma + charged_power_density + pohmpv
      tmp_cc =  1.0D0 - fradpwr * pradmaxpv / pradpv
      tmp_con = pradmaxpv * (1.0D0 - tmp_cc)
      tmp_err = pradpv * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW/m3'

   end subroutine constraint_eqn_017

   subroutine constraint_eqn_018(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for divertor heat load upper 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
      !! Equation for divertor heat load upper limit
      !! #=# divertor
      !! #=#=# fhldiv, hldivlim
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fhldiv : input real : f-value for divertor heat load
      !! hldivlim : input real : heat load limit (MW/m2)
      !! hldiv : input real : divertor heat load (MW/m2)
      use constraint_variables, only: fhldiv
      use divertor_variables, only: hldivlim, hldiv
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fhldiv * hldivlim/hldiv
      tmp_con = hldivlim * (1.0D0 - tmp_cc)
      tmp_err = hldiv * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW/m2'

   end subroutine constraint_eqn_018

   subroutine constraint_eqn_019(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for MVA (power) upper limit: resistive TF coil set
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for MVA upper limit
      !! #=# tfcoil
      !! #=#=# fmva, mvalim
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! tfcpmw : input real : peak resistive TF coil inboard leg power (total) (MW)
      !! tflegmw : input real : TF coil outboard leg resistive power (total) (MW)
      !! fmva : input real : f-value for maximum MVA
      !! mvalim : input real : MVA limit for resistive TF coil set (total) (MW)
      use tfcoil_variables, only: tfcpmw, tflegmw
      use constraint_variables, only: fmva, mvalim
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
      ! totmva : local real : total MVA in TF coil (MW)
      real(dp) :: totmva

      totmva = tfcpmw + tflegmw
      tmp_cc =  1.0D0 - fmva * mvalim/totmva
      tmp_con = mvalim * (1.0D0 - tmp_cc)
      tmp_err = totmva * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MVA'

   end subroutine constraint_eqn_019

   subroutine constraint_eqn_020(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for neutral beam tangency radius upper 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
      !! Equation for neutral beam tangency radius upper limit
      !! #=# current_drive
      !! #=#=# fportsz, rtanmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fportsz : input real : f-value for neutral beam tangency radius limit
      !! rtanmax : input real : maximum tangency radius for centreline of beam (m)
      !! rtanbeam : input real : neutral beam centreline tangency radius (m)
      use constraint_variables, only: fportsz
      use current_drive_variables, only: rtanmax, rtanbeam
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fportsz * rtanmax/rtanbeam
      tmp_con = rtanmax * (1.0D0 - tmp_cc)
      tmp_err = rtanbeam * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'm'

   end subroutine constraint_eqn_020

   subroutine constraint_eqn_021(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for minor radius 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
      !! Equation for minor radius lower limit
      !! #=# physics
      !! #=#=# frminor, aplasmin
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! frminor : input real : f-value for minor radius limit
      !! rminor : input real : plasma minor radius (m)
      !! aplasmin : input real : minimum minor radius (m)
      use constraint_variables, only: frminor
      use physics_variables, only: rminor
      use build_variables, only: aplasmin
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - frminor * rminor/aplasmin
      tmp_con = aplasmin * (1.0D0 - tmp_cc)
      tmp_err = aplasmin * tmp_cc
      tmp_symbol = '>'
      tmp_units = ''

   end subroutine constraint_eqn_021

   subroutine constraint_eqn_022(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for divertor collision/connection length ratio upper 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
      !! Equation for divertor collision/connection length ratio upper limit
      !! #=# divertor
      !! #=#=# fdivcol, rlenmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fdivcol : input real : f-value for divertor collisionality
      !! rlenmax : input real : maximum value for length ratio (rlclolcn)
      !! rlclolcn : input real : ratio of collision length / connection length
      use constraint_variables, only: fdivcol
      use divertor_variables, only: rlenmax, rlclolcn
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fdivcol * rlenmax / rlclolcn
      tmp_con = rlenmax * (1.0D0 - tmp_cc)
      tmp_err = rlclolcn * tmp_cc
      tmp_symbol = '<'
      tmp_units = ''

   end subroutine constraint_eqn_022

   subroutine constraint_eqn_023(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for conducting shell radius / rminor upper 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
      !! Equation for conducting shell radius / rminor upper limit
      !! #=# physics
      !! #=#=# fcwr, cwrmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! rminor : input real : plasma minor radius (m)
      !! scraplo : input real : gap between plasma and first wall, outboard side (m)
      !! fwoth : input real : outboard first wall thickness, initial estimate (m)
      !! blnkoth : input real : outboard blanket thickness (m)
      !! fcwr : input real : f-value for conducting wall radius / rminor limit
      !! cwrmax : input real : maximum ratio of conducting wall distance to plasma minor radius for vertical stability
      use physics_variables, only: rminor, cwrmax
      use build_variables, only: scraplo, fwoth, blnkoth
      use constraint_variables, only: fcwr
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
      ! rcw : local real : conducting shell radius (m)
      real(dp) :: rcw

      rcw = rminor + scraplo + fwoth + blnkoth
      tmp_cc =  1.0D0 - fcwr * cwrmax*rminor / rcw
      tmp_con = cwrmax*rminor * (1.0D0 - tmp_cc)
      tmp_err = rcw * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'm'

   end subroutine constraint_eqn_023

   subroutine constraint_eqn_024(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for beta upper 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
      !! Equation for beta upper limit
      !! #=# physics
      !! #=#=# fbetatry, betalim
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! iculbl : input integer : switch for beta limit scaling (constraint equation  24):<UL>
      !! <LI> = 0 apply limit to total beta;
      !! <LI> = 1 apply limit to thermal beta;
      !! <LI> = 2 apply limit to thermal + neutral beam beta
      !! <LI> = 3 apply limit to toroidal beta </UL>
      !! istell : input integer : switch for stellarator option (set via <CODE>device.dat</CODE>):<UL>
      !! <LI> = 0 use tokamak model;
      !! <LI> = 1 use stellarator model</UL>
      !! fbetatry : input real : f-value for beta limit
      !! betalim : input real : allowable beta
      !! beta : input real : total plasma beta (calculated if ipedestal =3)
      !! betaft : input real : fast alpha beta component
      !! beta_beam : input real : neutral beam beta component
      !! bt : input real : toroidal field
      !! btot : input real : total field
      use physics_variables, only: iculbl, betalim, beta, beta_beam, betaft, bt, btot
      use stellarator_variables, only: istell
      use constraint_variables, only: fbetatry
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! Include all beta components: relevant for both tokamaks and stellarators
      if ((iculbl == 0).or.(istell /= 0)) then
         tmp_cc =  1.0D0 - fbetatry * betalim/beta
         tmp_con = betalim
         tmp_err = betalim - beta / fbetatry
         tmp_symbol = '<'
         tmp_units = ''
      ! Here, the beta limit applies to only the thermal component, not the fast alpha or neutral beam parts
      else if (iculbl == 1) then
         tmp_cc = 1.0D0 - fbetatry * betalim/(beta-betaft-beta_beam)
         tmp_con = betalim
         tmp_err = betalim - (beta-betaft-beta_beam) / fbetatry
         tmp_symbol = '<'
         tmp_units = ''
      ! Beta limit applies to thermal + neutral beam: components of the total beta, i.e. excludes alphas
      else if (iculbl == 2) then
         tmp_cc = 1.0D0 - fbetatry * betalim/(beta-betaft)
         tmp_con = betalim * (1.0D0 - tmp_cc)
         tmp_err = (beta-betaft) * tmp_cc
         tmp_symbol = '<'
         tmp_units = ''
      ! Beta limit applies to toroidal beta
      else if (iculbl == 3) then
         tmp_cc =  1.0D0 - fbetatry * betalim/(beta*(btot/bt)**2)
         tmp_con = betalim
         tmp_err = betalim - (beta*(btot/bt)**2) / fbetatry
         tmp_symbol = '<'
         tmp_units = ''
      end if

   end subroutine constraint_eqn_024

   subroutine constraint_eqn_025(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for peak toroidal field upper 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
      !! Equation for peak toroidal field upper limit
      !! #=# tfcoil
      !! #=#=# fpeakb, bmxlim
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fpeakb : input real : f-value for maximum toroidal field
      !! bmxlim : input real : maximum peak toroidal field (T)
      !! bmaxtf : input real : mean peak field at TF coil (T)
      use constraint_variables, only: fpeakb, bmxlim
      use tfcoil_variables, only: bmaxtf
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fpeakb * bmxlim/bmaxtf
      tmp_con = bmxlim * (1.0D0 - tmp_cc)
      tmp_err = bmaxtf * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'T'

   end subroutine constraint_eqn_025

   subroutine constraint_eqn_026(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for Central Solenoid current density upper limit at EOF
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for Central Solenoid current density upper limit at EOF
      !! #=# pfcoil
      !! #=#=# fjohc, rjohc
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fjohc : input real : f-value for central solenoid current at end-of-flattop
      !! rjohc : input real : allowable central solenoid current density at end of flat-top (A/m2)
      !! coheof : input real : central solenoid overall current density at end of flat-top (A/m2)
      use constraint_variables, only: fjohc
      use pfcoil_variables, only: rjohc, coheof
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fjohc * rjohc/coheof
      tmp_con = rjohc
      tmp_err = rjohc - coheof / fjohc
      tmp_symbol = '<'
      tmp_units = 'A/m2'

   end subroutine constraint_eqn_026

   subroutine constraint_eqn_027(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for Central Solenoid current density upper limit at BOP
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for Central Solenoid current density upper limit at BOP
      !! #=# pfcoil
      !! #=#=# fjohc0, rjohc0
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fjohc0 : input real : f-value for central solenoid current at beginning of pulse
      !! rjohc0 : input real : allowable central solenoid current density at beginning of pulse (A/m2)
      !! cohbop : input real : central solenoid overall current density at beginning of pulse (A/m2)
      use constraint_variables, only: fjohc0
      use pfcoil_variables, only: rjohc0, cohbop
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fjohc0 * rjohc0/cohbop
      tmp_con = rjohc0
      tmp_err = rjohc0 - cohbop / fjohc0
      tmp_symbol = '<'
      tmp_units = 'A/m2'

   end subroutine constraint_eqn_027

   subroutine constraint_eqn_028(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for fusion gain (big Q) 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
      !! Equation for fusion gain (big Q) lower limit
      !! #=# physics
      !! #=#=# fqval, bigqmin
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fqval : input real : pf-value for Q
      !! bigq : input real : Fusion gain; P_fusion / (P_injection + P_ohmic)
      !! bigqmin : input real : minimum fusion gain Q
      !! ignite : input integer : switch for ignition assumption:<UL>
      !! <LI> = 0 do not assume plasma ignition;
      !! <LI> = 1 assume ignited (but include auxiliary power in costs)</UL>
      !! Obviously, ignite must be zero if current drive is required.
      !! If ignite=1, any auxiliary power is assumed to be used only
      !! during plasma start-up, and is excluded from all steady-state
      !! power balance calculations.
      use constraint_variables, only: fqval, bigqmin
      use current_drive_variables, only: bigq
      use physics_variables, only: ignite
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! if plasma is not ignited ...
      if (ignite == 0) then
         tmp_cc =  1.0D0 - fqval * bigq/bigqmin
         tmp_con = bigqmin * (1.0D0 - tmp_cc)
         tmp_err = bigqmin * tmp_cc
         tmp_symbol = '>'
         tmp_units = ''
      ! if plasma is ignited report error
      else
         tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ''
        tmp_units = ''
         call report_error(4)
      end if

   end subroutine constraint_eqn_028

   subroutine constraint_eqn_029(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for inboard major radius: This is a consistency equation
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for inboard major radius: This is a consistency equation
      !! #=# build
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! rmajor : input real : plasma major radius (m) (iteration variable 3)
      !! rminor : input real : plasma minor radius (m)
      !! rinboard : input real : plasma inboard radius (m)
      use physics_variables, only: rmajor, rminor
      use build_variables, only: rinboard
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - (rmajor - rminor) / rinboard
      tmp_con = rinboard * (1.0D0 - tmp_cc)
      tmp_err = rinboard * tmp_cc
      tmp_symbol = '='
      tmp_units = 'm'

   end subroutine constraint_eqn_029

   subroutine constraint_eqn_030(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for injection power upper 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
      !! Equation for injection power upper limit
      !! #=# current_drive
      !! #=#=# fpinj, pinjalw
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! pinjmw : input real : total auxiliary injected power (MW)
      !! fpinj : input real : f-value for injection power
      !! pinjalw : input real : Maximum allowable value for injected power (MW)
      use current_drive_variables, only: pinjmw, pinjalw
      use constraint_variables, only: fpinj
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fpinj*pinjalw/pinjmw
      tmp_con = pinjalw
      tmp_err = pinjalw  - pinjmw / fpinj
      tmp_symbol = '<'
      tmp_units = 'MW'

   end subroutine constraint_eqn_030

   subroutine constraint_eqn_031(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for TF coil case stress upper limit (SCTF)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for TF coil case stress upper limit (SCTF)
      !! #=# tfcoil
      !! #=#=# fstrcase, sig_tf_case_max
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fstrcase : input real : f-value for TF coil case stress
      !! sig_tf_case_max : input real : Allowable maximum shear stress in TF coil case (Tresca criterion) (Pa)
      !! sig_tf_case : input real : Constrained stress in TF coil case (Pa)
      use constraint_variables, only: fstrcase
      use tfcoil_variables, only: sig_tf_case_max, sig_tf_case
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      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_con = sig_tf_case_max
      tmp_err = sig_tf_case_max - sig_tf_case / fstrcase
      tmp_symbol = '<'
      tmp_units = 'Pa'

   end subroutine constraint_eqn_031

   subroutine constraint_eqn_032(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for TF coil conduit stress upper limit (SCTF)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for TF coil conduit stress upper limit (SCTF)
      !! #=# tfcoil
      !! #=#=# fstrcond, sig_tf_wp_max
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fstrcond : input real : f-value for TF coil conduit stress
      !! sig_tf_wp_max : input real : Allowable maximum shear stress in TF coil conduit (Tresca criterion) (Pa)
      !! sig_tf_wp : input real : Constrained stress in TF conductor conduit (Pa)
      use constraint_variables, only: fstrcond
      use tfcoil_variables, only: sig_tf_wp_max, sig_tf_wp
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      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_con = sig_tf_wp_max
      tmp_err = sig_tf_wp_max - sig_tf_wp / fstrcond
      tmp_symbol = '<'
      tmp_units = 'Pa'

   end subroutine constraint_eqn_032

   subroutine constraint_eqn_033(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for TF coil operating/critical J upper limit (SCTF)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for TF coil operating/critical J upper limit (SCTF)
      !! #=# tfcoil
      !! #=#=# fiooic, jwdgcrt
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fiooic : input real : f-value for TF coil operating current / critical
      !! jwdgcrt : input real : critical current density for winding pack (A/m2)
      !! jwptf : input real : winding pack current density (A/m2)
      use constraint_variables, only: fiooic
      use tfcoil_variables, only: jwdgcrt, jwptf
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      if (fiooic > 0.7D0) call report_error(285)
      tmp_cc =  1.0D0 - fiooic * jwdgcrt/jwptf
      tmp_con = jwdgcrt * (1.0D0 - tmp_cc)
      tmp_err = jwptf * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'A/m2'

   end subroutine constraint_eqn_033

   subroutine constraint_eqn_034(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for TF coil dump voltage upper limit (SCTF)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for TF coil dump voltage upper limit (SCTF)
      !! #=# tfcoil
      !! #=#=# fvdump, vdalw
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fvdump : input real : f-value for dump voltage
      !! vdalw : input real : max voltage across TF coil during quench (kV)
      !! vtfskv : input real : voltage across a TF coil during quench (kV)
      use constraint_variables, only: fvdump
      use tfcoil_variables, only: vdalw, vtfskv
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fvdump * vdalw/vtfskv
      tmp_con = vdalw
      tmp_err = vdalw - vtfskv
      tmp_symbol = '<'
      tmp_units = 'V'

   end subroutine constraint_eqn_034

   subroutine constraint_eqn_035(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for TF coil J_wp/J_prot upper limit (SCTF)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for TF coil J_wp/J_prot upper limit (SCTF)
      !! #=# tfcoil
      !! #=#=# fjprot, jwdgpro
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fjprot : input real : f-value for TF coil winding pack current density
      !! jwdgpro : input real : allowable TF coil winding pack current density, for dump temperature rise protection (A/m2)
      !! jwptf : input real : winding pack current density (A/m2)
      use constraint_variables, only: fjprot
      use tfcoil_variables, only: jwdgpro, jwptf
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fjprot * jwdgpro/jwptf
      tmp_con = jwdgpro
      tmp_err =  jwptf - jwdgpro
      tmp_symbol = '<'
      tmp_units = 'A/m2'

   end subroutine constraint_eqn_035

   subroutine constraint_eqn_036(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for TF coil s/c temperature margin lower limit (SCTF)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for TF coil s/c temperature margin lower limit (SCTF)
      !! #=# tfcoil
      !! #=#=# ftmargtf, tmargmin_tf
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftmargtf : input real : f-value for TF coil temperature margin
      !! tmargtf : input real : TF coil temperature margin (K)
      !! tmargmin_tf : input real : minimum allowable temperature margin : TF coils (K)
      use constraint_variables, only: ftmargtf
      use tfcoil_variables, only: tmargtf, tmargmin_tf
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - ftmargtf * tmargtf/tmargmin_tf
      tmp_con = tmargmin_tf
      tmp_err = tmargmin_tf - tmargtf
      tmp_symbol = '>'
      tmp_units = 'K'

   end subroutine constraint_eqn_036

   subroutine constraint_eqn_037(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for current drive gamma upper 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
      !! Equation for current drive gamma upper limit
      !! #=# current_drive
      !! #=#=# fgamcd, gammax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fgamcd : input real : f-value for current drive gamma
      !! gammax : input real : maximum current drive gamma
      !! gamcd : input real : normalised current drive efficiency (1.0e20 A/W-m2)
      use constraint_variables, only: fgamcd, gammax
      use current_drive_variables, only: gamcd
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fgamcd * gammax/gamcd
      tmp_con = gammax * (1.0D0 - tmp_cc)
      tmp_err = gamcd * tmp_cc
      tmp_symbol = '<'
      tmp_units = '1E20 A/Wm2'

   end subroutine constraint_eqn_037

   subroutine constraint_eqn_038(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Obsolete
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! Obsolete
      !! #=# empty
      !! #=#=# empty
      implicit none
	   ! Dummy formal arguments, for compliance with interface
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ''
        tmp_units = ''

   end subroutine constraint_eqn_038

   subroutine constraint_eqn_039(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for first wall temperature upper 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
      !! Equation for first wall temperature upper limit
      !! #=# fwbs
      !! #=#=# ftpeak, tfwmatmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftpeak : input real : f-value for first wall peak temperature
      !! tfwmatmax : input real : maximum temperature of first wall material (K) (secondary_cycle>1)
      !! tpeak : input real : peak first wall temperature (K)
      use constraint_variables, only: ftpeak
      use fwbs_variables, only: tfwmatmax, tpeak
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! If the temperature peak == 0 then report an error
      if (tpeak < 1.0D0) call report_error(5)
      tmp_cc =  1.0D0 - ftpeak * tfwmatmax/tpeak
      tmp_con = tfwmatmax * (1.0D0 - tmp_cc)
      tmp_err = tpeak * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'K'

   end subroutine constraint_eqn_039

   subroutine constraint_eqn_040(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for auxiliary power 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
      !! Equation for auxiliary power lower limit
      !! #=# current_drive
      !! #=#=# fauxmn, auxmin
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fauxmn : input real : f-value for minimum auxiliary power
      !! pinjmw : input real : total auxiliary injected power (MW)
      !! auxmin : input real : minimum auxiliary power (MW)
      use constraint_variables, only: fauxmn, auxmin
      use current_drive_variables, only: pinjmw
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fauxmn * pinjmw/auxmin
      tmp_con = auxmin * (1.0D0 - tmp_cc)
      tmp_err = auxmin * tmp_cc
      tmp_symbol = '>'
      tmp_units = 'MW'

   end subroutine constraint_eqn_040

   subroutine constraint_eqn_041(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for plasma current ramp-up time 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
      !! Equation for plasma current ramp-up time lower limit
      !! #=# times
      !! #=#=# ftohs, tohsmn
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftohs : input real : f-value for plasma current ramp-up time
      !! tohs : input real : plasma current ramp-up time for current initiation (s)
      !! tohsmn : input real : minimum plasma current ramp-up time (s)
      use constraint_variables, only: ftohs, tohsmn
      use times_variables, only: tohs
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - ftohs * tohs/tohsmn
      tmp_con = tohsmn * (1.0D0 - tmp_cc)
      tmp_err = tohsmn * tmp_cc
      tmp_symbol = '>'
      tmp_units = 'sec'

   end subroutine constraint_eqn_041

   subroutine constraint_eqn_042(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for cycle time 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
      !! Equation for cycle time lower limit
      !! #=# times
      !! #=#=# ftcycl, tcycmn
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftcycl : input real : f-value for cycle time
      !! tcycle : input real : full cycle time (s)
      !! tcycmn : input real : minimum cycle time (s)
      use constraint_variables, only: ftcycl, tcycmn
      use times_variables, only: tcycle
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! if the minimum cycle time == 0 report an error
      if (tcycmn < 1.0D0) call report_error(6)
      tmp_cc =  1.0D0 - ftcycl * tcycle/tcycmn
      tmp_con = tcycmn * (1.0D0 - tmp_cc)
      tmp_err = tcycmn * tmp_cc
      tmp_symbol = '>'
      tmp_units = 'sec'

   end subroutine constraint_eqn_042

   subroutine constraint_eqn_043(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for average centrepost temperature: This is a consistency equation (TART)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for average centrepost temperature: This is a consistency equation (TART)
      !! #=# tfcoil
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! tcpav : input real : average temp of TF coil inboard leg conductor (C)e
      !! tcpav2 : input real : centrepost average temperature (C) (for consistency)
      !! 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 tfcoil_variables, only: tcpav, tcpav2
      use physics_variables, only: itart
      use tfcoil_variables, only:  i_tf_sup

      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! if the machine isn't a ST then report error
      if (itart == 0) call report_error(7)

      ! For some reasons these lines are needed to make VMCON CONVERGE ....
      if ( i_tf_sup == 0 ) then ! Copper case
         tcpav = tcpav - 273.15D0
         tcpav2 = tcpav2 - 273.15D0
      end if

      tmp_cc =   1.0D0 - tcpav/tcpav2
      tmp_con = tcpav2 * (1.0D0 - tmp_cc)
      tmp_err = tcpav2 * tmp_cc
      tmp_symbol = '='
      tmp_units = 'deg C'

      ! For some reasons these lines are needed to make VMCON CONVERGE ....
      if ( i_tf_sup == 0 ) then ! Copper case
         tcpav = tcpav + 273.15D0
         tcpav2 = tcpav2 + 273.15D0
      end if



   end subroutine constraint_eqn_043

   subroutine constraint_eqn_044(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for centrepost temperature upper limit (TART)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for centrepost temperature upper limit (TART)
      !! #=# tfcoil
      !! #=#=# fptemp, ptempalw
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fptemp : input real : f-value for peak centrepost temperature
      !! ptempalw : input real : maximum peak centrepost temperature (K)
      !! tcpmax : input real :  peak centrepost temperature (K)
      !! 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: fptemp
      use tfcoil_variables, only: ptempalw, tcpmax
      use physics_variables, only: itart
      use tfcoil_variables, only:  i_tf_sup
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! if the machine isn't a ST then report error
      if (itart == 0) call report_error(8)

      ! For some reasons these lines are needed to make VMCON CONVERGE ....
      if ( i_tf_sup == 0 ) then ! Copper case
         ptempalw = ptempalw - 273.15D0
         tcpmax = tcpmax - 273.15D0
      end if

      tmp_cc =   1.0D0 - fptemp * ptempalw / tcpmax
      tmp_con = ptempalw * (1.0D0 - tmp_cc)
      tmp_err = tcpmax * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'deg C'

      ! For some reasons these lines are needed to make VMCON CONVERGE ....
      if ( i_tf_sup == 0 ) then ! Copper case
         ptempalw = ptempalw + 273.15D0
         tcpmax = tcpmax + 273.15D0
      end if

   end subroutine constraint_eqn_044

   subroutine constraint_eqn_045(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for edge safety factor lower limit (TART)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for edge safety factor lower limit (TART)
      !! #=# tfcoil
      !! #=#=# fq, qlim
      !! 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
      !! (unless i_plasma_current = 2 (ST current scaling), in which case q = mean edge safety factor qbar)
      !! qlim : 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, qlim, itart
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: 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/qlim
      tmp_con = qlim * (1.0D0 - tmp_cc)
      tmp_err = qlim * tmp_cc
      tmp_symbol = '<'
      tmp_units = ''

   end subroutine constraint_eqn_045

   subroutine constraint_eqn_046(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for Ip/Irod upper limit (TART)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for Ip/Irod upper limit (TART)
      !! #=# tfcoil
      !! #=#=# fipir, cratmx
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! eps : input real :  inverse aspect ratio
      !! fipir : input real : f-value for Ip/Irod upper limit
      !! ritfc : input real : total (summed) current in TF coils (A)
      !! plasma_current : input real :  plasma current (A)
      !! 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 physics_variables, only: eps, plasma_current, itart
      use constraint_variables, only: fipir
      use tfcoil_variables, only: ritfc
      implicit none
      ! cratmx : local real : maximum ratio of plasma current to centrepost current
      real(dp) :: cratmx
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: 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 * ritfc/plasma_current
      tmp_con = cratmx * (1.0D0 - tmp_cc)
      tmp_err = plasma_current/ritfc * tmp_cc
      tmp_symbol = '<'
      tmp_units = ''

   end subroutine constraint_eqn_046

   subroutine constraint_eqn_047(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! 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
      !! Equation for TF coil toroidal thickness upper limit
      !! #=# empty
      !! #=#=# empty
      implicit none
      ! Dummy formal arguments, just to comply with the subroutine interface
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ''
        tmp_units = ''

   end subroutine constraint_eqn_047

   subroutine constraint_eqn_048(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for poloidal beta upper 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
      !! Equation for poloidal beta upper limit
      !! #=# physics
      !! #=#=# fbetap, betpmx
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fbetap : input real : rf-value for poloidal beta
      !! betpmx : input real :  maximum poloidal beta
      !! betap : input real :  poloidal beta
      use constraint_variables, only: fbetap, betpmx
      use physics_variables, only: betap
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fbetap * betpmx/betap
      tmp_con = betpmx * (1.0D0 - tmp_cc)
      tmp_err = betap * tmp_cc
      tmp_symbol = '<'
      tmp_units = ''

   end subroutine constraint_eqn_048

   subroutine constraint_eqn_049(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! 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
      !! #=# empty
      !! #=#=# empty
      ! Dummy formal arguments, just to comply with the subroutine interface
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ''
        tmp_units = ''

   end subroutine constraint_eqn_049

   subroutine constraint_eqn_050(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! IFE option: Equation for repetition rate upper limit
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! author: S I Muldrew, CCFE, Culham Science Centre
      !! IFE option: Equation for repetition rate upper limit
      !! #=# IFE
      !! #=#=# frrmax, rrmax
      use ife_variables, only: frrmax, ife, rrmax, reprat
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      if (ife /= 1) then
         call report_error(12)
      end if

      tmp_cc =  1.0D0 - frrmax * rrmax/reprat
      tmp_con = rrmax * (1.0D0 - tmp_cc)
      tmp_err = reprat * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'Hz'

   end subroutine constraint_eqn_050

   subroutine constraint_eqn_051(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation to enforce startup flux = available startup flux
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation to enforce startup flux = available startup flux
      !! #=# pfcoil
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! vsres : input real : resistive losses in startup V-s (Wb)
      !! vsind : input real :  internal and external plasma inductance V-s (Wb))
      !! vssu : input real :  total flux swing for startup (Wb)
      use physics_variables, only: vsres, vsind
      use pfcoil_variables, only: vssu, fvssu
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fvssu * abs((vsres+vsind) / vssu)
      tmp_con = vssu * (1.0D0 - tmp_cc)
      tmp_err = vssu * tmp_cc
      tmp_symbol = '='
      tmp_units = 'V.s'

   end subroutine constraint_eqn_051

   subroutine constraint_eqn_052(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for tritium breeding ratio 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
      !! Equation for tritium breeding ratio lower limit
      !! #=# fwbs
      !! #=#=# ftbr, tbrmin
      !! ? TODO should this only be for certain blanket models ?
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftbr : input real : f-value for minimum tritium breeding ratio
      !! tbr : input real :  tritium breeding ratio (iblanket=2,3 (KIT HCPB/HCLL))
      !! tbrmin : input real :  minimum tritium breeding ratio (If iblanket=1, tbrmin=minimum 5-year time-averaged tritium breeding ratio)
      use constraint_variables, only: ftbr, tbrmin
      use fwbs_variables, only: tbr
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - ftbr * tbr/tbrmin
      tmp_con = tbrmin * (1.0D0 - tmp_cc)
      tmp_err = tbrmin * tmp_cc
      tmp_symbol = '>'
      tmp_units = ''

   end subroutine constraint_eqn_052

   subroutine constraint_eqn_053(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for fast neutron fluence on TF coil upper 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
      !! Equation for fast neutron fluence on TF coil upper limit
      !! #=# fwbs
      !! #=#=# fflutf, nflutfmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fflutf : input real : f-value for maximum TF coil nuclear heating
      !! nflutfmax : input real :  max fast neutron fluence on TF coil (n/m2)
      !! nflutf : input real :  peak fast neutron fluence on TF coil superconductor (n/m2)
      use constraint_variables, only: fflutf, nflutfmax
      use fwbs_variables, only: nflutf
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - fflutf * nflutfmax/nflutf
      tmp_con = nflutfmax * (1.0D0 - tmp_cc)
      tmp_err = nflutf * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'neutron/m2'

   end subroutine constraint_eqn_053

   subroutine constraint_eqn_054(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for peak TF coil nuclear heating upper 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
      !! Equation for peak TF coil nuclear heating upper limit
      !! #=# fwbs
      !! #=#=# fptfnuc, ptfnucmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fptfnuc : input real : f-value for maximum TF coil nuclear heating
      !! ptfnucmax : input real :  maximum nuclear heating in TF coil (MW/m3)
      !! ptfnucpm3 : input real :  nuclear heating in the TF coil (MW/m3) (blktmodel>0)
      use constraint_variables, only: fptfnuc, ptfnucmax
      use fwbs_variables, only: ptfnucpm3
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - fptfnuc * ptfnucmax/ptfnucpm3
      tmp_con = ptfnucmax * (1.0D0 - tmp_cc)
      tmp_err = ptfnucpm3 * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW/m3'

   end subroutine constraint_eqn_054

   subroutine constraint_eqn_055(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! vvhemax is no longer calculated in PROCESS and this constraint is disabled
      implicit none

      real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      call report_error(173)
   end subroutine constraint_eqn_055

   subroutine constraint_eqn_056(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for power through separatrix / major radius upper 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
      !! Equation for power through separatrix / major radius upper limit
      !! #=# current_drive
      !! #=#=# fnbshinef, nbshinefmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fpsepr : input real : f-value for maximum Psep/R limit
      !! pseprmax : input real :  maximum ratio of power crossing the separatrix to plasma major radius (Psep/R) (MW/m)
      !! pdivt : input real :  power to be conducted to the divertor region (MW)
      !! rmajor : input real :  plasma major radius (m)
      use constraint_variables, only: fpsepr, pseprmax
      use physics_variables, only: pdivt, rmajor
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - fpsepr * pseprmax / (pdivt/rmajor)
      tmp_con = pseprmax * (1.0D0 - tmp_cc)
      tmp_err = (pdivt/rmajor) * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW/m'

   end subroutine constraint_eqn_056

   subroutine constraint_eqn_057(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Obsolete
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! Obsolete
      !! #=# empty
      !! #=#=# empty
      ! Dummy formal arguments, just to comply with the subroutine interface
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ''
        tmp_units = ''

   end subroutine constraint_eqn_057

   subroutine constraint_eqn_058(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Obsolete
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! Obsolete
      !! #=# empty
      !! #=#=# empty
      ! Dummy formal arguments, just to comply with the subroutine interface
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 0
        tmp_con = 0
        tmp_err = 0
        tmp_symbol = ''
        tmp_units = ''

   end subroutine constraint_eqn_058

   subroutine constraint_eqn_059(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for neutral beam shine-through fraction upper 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
      !! Equation for neutral beam shine-through fraction upper limit
      !! #=# current_drive
      !! #=#=# fnbshinef, nbshinefmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fnbshinef : input real : f-value for maximum neutral beam shine-through fraction
      !! nbshinefmax : input real :  maximum neutral beam shine-through fraction
      !! nbshinef : input real :  neutral beam shine-through fraction
      use constraint_variables, only: fnbshinef, nbshinefmax
      use current_drive_variables, only: nbshinef
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      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_con = nbshinefmax * (1.0D0 - tmp_cc)
      tmp_err = nbshinef * tmp_cc
      tmp_symbol = '<'
      tmp_units = ''
   end subroutine constraint_eqn_059

   subroutine constraint_eqn_060(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for Central Solenoid s/c temperature margin 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
      !! Equation for Central Solenoid s/c temperature margin lower limit
      !! #=# tfcoil
      !! #=#=# ftmargoh, tmargmin_cs
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftmargoh : input real :  f-value for central solenoid temperature margin
      !! tmargoh : input real :  Central solenoid temperature margin (K)
      !! tmargmin_cs : input real :  Minimum allowable temperature margin : CS (K)
      use constraint_variables, only: ftmargoh
      use pfcoil_variables, only: tmargoh
      use tfcoil_variables, only: tmargmin_cs
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - ftmargoh * tmargoh/tmargmin_cs
      tmp_con = tmargmin_cs
      tmp_err = tmargmin_cs - tmargoh
      tmp_symbol = '>'
      tmp_units = 'K'

   end subroutine constraint_eqn_060

   subroutine constraint_eqn_061(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for availability 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
      !! Equation for availability limit
      !! #=# cost
      !! #=#=# favail, avail_min
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! favail : input real : F-value for minimum availability
      !! cfactr : input real : Total plant availability fraction
      !! avail_min : input real : Minimum availability
      use cost_variables, only: favail, cfactr, avail_min
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - favail * cfactr / avail_min
      tmp_con = avail_min * (1.0D0 - tmp_cc)
      tmp_err = cfactr * tmp_cc
      tmp_symbol = '>'
      tmp_units = ''

   end subroutine constraint_eqn_061

   subroutine constraint_eqn_062(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Lower limit on taup/taueff the ratio of alpha particle to energy confinement times
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Lower limit on taup/taueff the ratio of alpha particle to energy confinement times
      !! #=# physics
      !! #=#=# ftaulimit, taulimit
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftaulimit : input real : f-value for lower limit on taup/taueff the ratio of alpha particle to energy confinement
      !! taup : input real : alpha particle confinement time (s)
      !! taueff : input real : global thermal energy confinement time (sec)
      !! taulimit : input real : Lower limit on taup/taueff the ratio of alpha particle to energy confinement times
      use constraint_variables, only: ftaulimit, taulimit
      use physics_variables, only: taup, taueff
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - ftaulimit * (taup / taueff) / taulimit
      tmp_con = taulimit
      tmp_err = (taup / taueff) * tmp_cc
      tmp_symbol = '>'
      tmp_units = ''

   end subroutine constraint_eqn_062

   subroutine constraint_eqn_063(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Upper limit on niterpump (vacuum_model = simple)
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Upper limit on niterpump (vacuum_model = simple)
      !! #=# vacuum
      !! #=#=# fniterpump, tfno
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fniterpump : input real : f-value for constraint that number of pumps < tfno
      !! tfno : input real : number of TF coils (default = 50 for stellarators)
      !! niterpump : input real : number of high vacuum pumps (real number), each with the throughput
      use constraint_variables, only: fniterpump
      use tfcoil_variables, only: n_tf
      use vacuum_variables, only: niterpump
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - fniterpump * n_tf / niterpump
      tmp_con = n_tf
      tmp_err = n_tf * tmp_cc
      tmp_symbol = '<'
      tmp_units = ''

   end subroutine constraint_eqn_063

   subroutine constraint_eqn_064(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Upper limit on Zeff
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Upper limit on Zeff
      !! #=# physics
      !! #=#=# fzeffmax, zeffmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fzeffmax : input real : f-value for maximum zeff
      !! zeffmax : input real : maximum value for Zeff
      !! zeff : input real : plasma effective charge
      use constraint_variables, only: fzeffmax, zeffmax
      use physics_variables, only: zeff
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - fzeffmax * (zeffmax/zeff)
      tmp_con = zeffmax
      tmp_err = zeffmax * tmp_cc
      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.
      !! author: Timothy Nunn, UKAEA
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fmaxvvstress : input real : f-value for constraint on maximum VV stress
      !! max_vv_stress : input real : Maximum permitted stress of the VV (Pa)
      !! vv_stress_quench : input real : Stress of the VV (Pa)
      use constraint_variables, only: fmaxvvstress
      use tfcoil_variables, only: max_vv_stress
      use sctfcoil_module, only: vv_stress_quench
      implicit none
      real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      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_con = max_vv_stress
      tmp_err = max_vv_stress * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'Pa'

   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
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Limit on rate of change of energy in poloidal field
      !! #=# pfcoil
      !! #=#=# fpoloidalpower, maxpoloidalpower
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fpoloidalpower : input real : f-value for constraint on rate of change of energy in poloidal field
      !! maxpoloidalpower : input real : Maximum permitted absolute rate of change of stored energy in poloidal field (MW)
      !! peakpoloidalpower : input real : Peak absolute rate of change of stored energy in poloidal field (MW) (11/01/16)
      use constraint_variables, only: fpoloidalpower
      use pf_power_variables, only: maxpoloidalpower, peakpoloidalpower
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0d0 - fpoloidalpower * maxpoloidalpower / peakpoloidalpower
      tmp_con = maxpoloidalpower
      tmp_err = maxpoloidalpower * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW'

   end subroutine constraint_eqn_066

   subroutine constraint_eqn_067(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Simple upper limit on radiation wall load
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Simple upper limit on radiation wall load
      !! #=# physics
      !! #=#=# fradwall, maxradwallload
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fradwall : input real : f-value for upper limit on radiation wall load
      !! maxradwallload : input real : Maximum permitted radiation wall load (MW/m^2)
      !! peakradwallload : input real : Peak radiation wall load (MW/m^2)
      use constraint_variables, only: fradwall, maxradwallload, peakradwallload
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0d0 - fradwall * maxradwallload / peakradwallload
      tmp_con = maxradwallload
      tmp_err =  maxradwallload * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW/m^2'

   end subroutine constraint_eqn_067

   subroutine constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! New 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
      !! New Psep scaling (PsepB/qAR)
      !! Issue #464
      !! #=# physics
      !! #=#=# fpsepbqar, psepbqarmax
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fpsepbqar : input real : f-value for upper limit on psepbqar, maximum Psep*Bt/qAR limit
      !! psepbqarmax : input real : maximum permitted value of ratio of Psep*Bt/qAR (MWT/m)
      !! pdivt : input real : Power to conducted to the divertor region (MW)
      !! bt : input real : toroidal field on axis (T) (iteration variable 2)
      !! q95 : input real : safety factor q at 95% flux surface
      !! aspect : input real : aspect ratio (iteration variable 1)
      !! rmajor : input real : plasma major radius (m) (iteration variable 3)
      use constraint_variables, only: fpsepbqar, psepbqarmax
      use physics_variables, only: pdivt, bt, q95, aspect, rmajor
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0d0 - fpsepbqar * psepbqarmax / ((pdivt*bt)/(q95*aspect*rmajor))
      tmp_con = psepbqarmax
      tmp_err = (pdivt*bt)/(q95*aspect*rmajor) - psepbqarmax
      tmp_symbol = '<'
      tmp_units = 'MWT/m'

   end subroutine constraint_eqn_068

   subroutine constraint_eqn_069(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! 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;
      !! residual error in physical units; output string; units string
      !! Ensure separatrix power is less than value from Kallenbach divertor
      !! #=# divertor_kallenbach
      !! #=#=# consistency, psep_kallenbach
      !! fpsep has been removed from the equation.
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! psep_kallenbach : input real : Power conducted through the separatrix, as calculated by the divertor model [W]
      !! pdivt : input real :  power to conducted to the divertor region (MW)
      ! use div_kal_vars, only: psep_kallenbach
      use physics_variables, only: pdivt
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! From Kallenbach model, should be reserved if the model is going to be added back

   end subroutine constraint_eqn_069

   subroutine constraint_eqn_070(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Separatrix density consistency
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Separatrix density consistency
      !! #=# divertor_kallenbach
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! teomp : input real : Separatrix temperature calculated by the Kallenbach divertor model [eV]
      !! tesep : input real : Electron temperature at separatrix [keV]
      ! use div_kal_vars, only: teomp
      use  physics_variables, only: tesep
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! From Kallenbach model, should be reserved if the model is going to be added back

   end subroutine constraint_eqn_070

   subroutine constraint_eqn_071(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Separatrix density consistency
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Separatrix density consistency
      !! #=# divertor_kallenbach
      !! #=#=# consistency
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! neomp : input real : Mean SOL density at OMP calculated by the Kallenbach divertor model [m-3]
      !! nesep : input real :  electron density at separatrix [m-3] (ipedestal=1,2, calculated if 3)
      !! neratio : input real : Ratio of mean SOL density at OMP to separatrix density at OMP (iteration variable 121)
      ! use div_kal_vars, only: neomp, neratio
      use physics_variables, only: nesep
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! From Kallenbach model, should be reserved if the model is going to be added back

   end subroutine constraint_eqn_071

   subroutine constraint_eqn_072(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Central Solenoid Tresca yield criterion
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Central Solenoid Tresca yield criterion
      !! #=# pfcoil
      !! #=#=# foh_stress, alstroh
      !! In the case if the bucked and wedged option ( i_tf_bucking >= 2 ) the constrained
      !! stress is the largest the largest stress of the
      !!  - CS stress at maximum current (conservative as the TF inward pressure is not taken
      !!    into account)
      !!  - CS stress at flux swing (no current in CS) from the TF inward pressure
      !! This allow to cover the 2 worst stress scenario in the bucked and wedged design
      !! Otherwise (free standing TF), the stress limits are only set by the CS stress at max current
      !! Reverse the sign so it works as an inequality constraint (tmp_cc > 0)
      !! This will have no effect if it is used as an equality constraint because it will be squared.
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! foh_stress : input real : f-value for Tresca yield criterion in Central Solenoid
      !! alstroh : input real :  allowable hoop stress in Central Solenoid structural material (Pa)
      !! s_tresca_oh : input real : Maximum shear stress coils/central solenoid (Pa)
      !! sig_tf_cs_bucked : input real : Maximum shear stress in CS case at flux swing (no current in CS)
      !!                       can be significant for the bucked and weged design
      !! i_tf_bucking : input integer : switch for TF structure design
      use constraint_variables, only: foh_stress
      use pfcoil_variables, only: alstroh, s_tresca_oh
      use tfcoil_variables, only: sig_tf_cs_bucked, i_tf_bucking
      use build_variables, only: tf_in_cs
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! bucked and wedged desing (see subroutine comment)
      if ( i_tf_bucking >= 2 .and. tf_in_cs == 0 ) then
         tmp_cc = 1.0d0 - foh_stress * alstroh / max(s_tresca_oh, sig_tf_cs_bucked)
         tmp_err = alstroh - max(s_tresca_oh, sig_tf_cs_bucked)
      ! Free standing CS
      else
         tmp_cc = 1.0d0 - foh_stress * alstroh / s_tresca_oh
         tmp_err = alstroh - s_tresca_oh
      end if

      tmp_con = alstroh
      tmp_symbol = '<'
      tmp_units = 'Pa'

   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
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Ensure separatrix power is greater than the L-H power + auxiliary power
      !! #=# physics
      !! #=#=# fplhsep, pdivt
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fplhsep : input real : F-value for Psep >= Plh + Paux : for consistency of two values of separatrix power
      !! plhthresh : input real : L-H mode power threshold (MW)
      !! pdivt : input real : power to be conducted to the divertor region (MW)
      !! pinjmw : inout real : total auxiliary injected power (MW)
      use physics_variables, only: fplhsep, plhthresh, pdivt
      use current_drive_variables, only: pinjmw
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0d0 - fplhsep * pdivt / (plhthresh+pinjmw)
      tmp_con = pdivt
      tmp_err = pdivt * tmp_cc
      tmp_symbol = '>'
      tmp_units = 'MW'

   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
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Ensure TF coil quench temperature < tmax_croco ONLY used for croco HTS coil
      !! #=# physics
      !! #=#=# fcqt, tmax_croco
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fcqt : input real : f-value: TF coil quench temparature remains below tmax_croco
      !! croco_quench_temperature : input real : CroCo strand: Actual temp reached during a quench (K)
      !! tmax_croco : input real : CroCo strand: maximum permitted temp during a quench (K)
      use constraint_variables, only: fcqt
      use tfcoil_variables, only: croco_quench_temperature, tmax_croco
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      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_con = croco_quench_temperature
      tmp_err = croco_quench_temperature * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'K'

   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
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Ensure that TF coil current / copper area < Maximum value
      !! ONLY used for croco HTS coil
      !! #=# physics
      !! #=#=# f_coppera_m2, copperA_m2_max
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! copperA_m2 : input real : TF coil current / copper area (A/m2)
      !! copperA_m2_max : input real : Maximum TF coil current / copper area (A/m2)
      !! f_coppera_m2 : input real : f-value for TF coil current / copper area < copperA_m2_max
      use rebco_variables, only: copperA_m2, copperA_m2_max, f_coppera_m2
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      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_con = copperA_m2
      tmp_err = copperA_m2 * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'A/m2'

   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
      !! author: P B Lloyd, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Eich critical separatrix density model
      !! Added for issue 558 with ref to http://iopscience.iop.org/article/10.1088/1741-4326/aaa340/pdf
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! alpha_crit : output real : critical ballooning parameter value
      !! nesep_crit : output real : critical electron density at separatrix [m-3]
      !! kappa : input real : plasma separatrix elongation (calculated if ishape = 1-5, 7 or 9)
      !! triang : input real : plasma separatrix triangularity (calculated if ishape = 1, 3-5 or 7)
      !! aspect : input real : aspect ratio (iteration variable 1)
      !! pdivt : input real : power to conducted to the divertor region (MW)
      !! dlimit(7) : input real array : density limit (/m3) as calculated using various models
      !! fnesep : input real : f-value for Eich critical separatrix density
      use physics_variables, only: alpha_crit, nesep_crit, kappa, triang, &
                                   aspect, pdivt, dlimit, nesep
      use constraint_variables, only: fnesep
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      alpha_crit = (kappa ** 1.2D0) * (1.0D0 + 1.5D0 * triang)
      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_con = nesep
      tmp_err = nesep * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'm-3'

   end subroutine constraint_eqn_076

   subroutine constraint_eqn_077(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for maximum TF current per turn upper 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
      !! Equation for maximum TF current per turn upper limit
      !! #=# tfcoil
      !! #=#=# fcpttf, cpttf, cpttf_max
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fcpttf : input : f-value for TF coil current per turn
      !! cpttf_max  : input : allowable TF coil current per turn [A/turn]
      !! cpttf  : input : TF coil current per turn [A/turn]
      use constraint_variables, only: fcpttf
      use tfcoil_variables, only: cpttf_max, cpttf
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc = 1.0D0 - fcpttf * cpttf_max/cpttf
      tmp_con = cpttf_max
      tmp_err = cpttf_max * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'A/turn'

   end subroutine constraint_eqn_077

   subroutine constraint_eqn_078(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for Reinke criterion, divertor impurity fraction 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
      !! Equation for Reinke criterion, divertor impurity fraction lower limit
      !! #=# divertor
      !! #=#=# freinke, fzactual, fzmin
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present;
      !! and con will be printed out only if present. Thesw conditions were missing.
      !! freinke : input : f-value for Reinke criterion (itv 147)
      !! fzmin : input : minimum impurity fraction from Reinke model
      !! fzactual : input : actual impurity fraction
      use constraint_variables, only: freinke
      use reinke_variables, only: fzactual, fzmin
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! write(*,*) 'freinke, fzact, fzmin = ', freinke, ', ', fzactual, ', ', fzmin
      !            1.0,    0.0,   value
      tmp_cc = 1.0D0 - freinke *  fzactual/fzmin
      !The following two pre-existing lines are not understood:
      !KE note - cc is always 1, code never enters IF statement...
      tmp_con = fzmin * (1.0D0 - tmp_cc)
      tmp_err = fzmin * tmp_cc
      tmp_symbol = '>'
      tmp_units  = ''
      ! write(*,*) 'cc, con = ', tmp_cc, ', ', tmp_con
      ! write(*,*) 'freinke, fzactual, fzmin = ', freinke, ', ', fzactual, ', ', fzmin

   end subroutine constraint_eqn_078

   subroutine constraint_eqn_079(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for maximum CS 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
      !! Equation for maximum CS field
      !! #=# pfcoil
      !! #=#=# fbmaxcs, bmaxoh, bmaxoh0, bmaxcs_lim
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fbmaxcs : input : F-value for CS mmax field (cons. 79, itvar 149)
      !! bmaxcs_lim : input : Central solenoid max field limit [T]
      !! bmaxoh0 : input : maximum field in central solenoid at beginning of pulse (T)
      !! bmaxoh : input real : maximum field in central solenoid at end of flat-top (EoF) (T)
      !! (Note: original code has "bmaxoh/bmaxoh0 |  peak CS field [T]".)
      use pfcoil_variables, only: fbmaxcs, bmaxcs_lim, bmaxoh0, bmaxoh
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc     = 1.0D0 - fbmaxcs * bmaxcs_lim/max(bmaxoh, bmaxoh0)
      tmp_con    = bmaxcs_lim
      tmp_err    = max(bmaxoh, bmaxoh0) * tmp_cc
      tmp_symbol = '<'
      tmp_units  = 'A/turn'

   end subroutine constraint_eqn_079

   subroutine constraint_eqn_080(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for pdivt lower limit
      !! author: J Morris, Culham Science Centre
      !! args : output structure : residual error; constraint value; residual error in physical units;
      !! output string; units string
      !! Lower limit pdivt
      !! #=# physics
      !! #=#=# fpdivlim, pdivt
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fpdivlim : input : F-value for lower limit on pdivt (cons. 80, itvar 153)
      !! pdivtlim : input : Minimum power crossing separatrix pdivt [MW]
      !! pdivt : input : Power crossing separatrix [MW]
      use physics_variables, only: fpdivlim, pdivt
      use constraint_variables, only : pdivtlim
      implicit none

            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
      tmp_cc     = 1.0D0 - fpdivlim * pdivt / pdivtlim
      tmp_con    = pdivtlim
      tmp_err    = pdivt * tmp_cc
      tmp_symbol = '>'
      tmp_units  = 'MW'

   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
      !! author: S Kahn, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Lower limit ne0 > neped
      !! !#=# physics
      !! !#=#=# ne0, neped
      !! Logic change during pre-factoring: err, symbol, units will be
      !! assigned only if present.
      !! fne0  : input : F-value for constraint on ne0 > neped
      !! ne0   : input : Central electron density [m-3]
      !! neped : input : Electron density at pedestal [m-3]
      use physics_variables, only: ne0, fne0, neped
      implicit none

            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
      tmp_cc     = 1.0D0 - fne0 * ne0/neped
      tmp_con    = fne0
      tmp_err    = fne0 * tmp_cc
      tmp_symbol = '>'
      tmp_units  = '/m3'

   end subroutine constraint_eqn_081

   subroutine constraint_eqn_082(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for toroidal consistency of stellarator build
      !! author: J Lion, IPP Greifswald
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! toroidalgap > tftort
      !! #=# tfcoil
      !! #=#=# tftort, ftoroidalgap
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! ftoroidalgap : input real : f-value for constraint toroidalgap > tftort
      !! toroidalgap : input real :  minimal gap between two stellarator coils
      !! tftort : input real :  total toroidal width of a tf coil
      use tfcoil_variables, only: tftort,ftoroidalgap,toroidalgap
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - ftoroidalgap * toroidalgap/tftort
      tmp_con = toroidalgap
      tmp_err = toroidalgap - tftort/ftoroidalgap
      tmp_symbol = '<'
      tmp_units = 'm'

   end subroutine constraint_eqn_082

   subroutine constraint_eqn_083(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for radial consistency of stellarator build
      !! author: J Lion, IPP Greifswald
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! available_radial_space > required_radial_space
      !! #=# build
      !! #=#=# required_radial_space, f_avspace
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! f_avspace : input real : f-value for constraint available_radial_space > required_radial_space
      !! available_radial_space : input real :  avaible space in radial direction as given by each s.-configuration
      !! required_radial_space : input real :  required space in radial direction
      use build_variables, only: available_radial_space, required_radial_space, f_avspace
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - f_avspace  * available_radial_space/required_radial_space
      tmp_con = available_radial_space * (1.0D0 - tmp_cc)
      tmp_err = required_radial_space * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'm'
   end subroutine constraint_eqn_083

   subroutine constraint_eqn_084(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for the lower limit of beta
      !! author: J Lion, IPP Greifswald
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !!  (beta-betaft) > betalim_lower
      !! #=# physics
      !! #=#=# betaft, beta, fbetatry_lower
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fbetatry_lower : input real : f-value for constraint beta-betaft > betalim_lower
      !! betalim_lower : input real :  Lower limit for beta
      !! beta : input real :  plasma beta
      !! betaft : input real : Alpha particle beta

      use physics_variables, only: betalim_lower, beta, betaft
      use constraint_variables, only: fbetatry_lower
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units


      tmp_cc = 1.0D0 - fbetatry_lower * (beta-betaft)/betalim_lower
      tmp_con = betalim_lower * (1.0D0 - tmp_cc)
      tmp_err = (beta-betaft) * tmp_cc
      tmp_symbol = '>'
      tmp_units = ''


   end subroutine constraint_eqn_084

   subroutine constraint_eqn_085(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! 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
      !!  - 2 : The CP lifetime is equal to the breeding blankets one
      !!  - 3 : The CP lifetime is equal to the plant one
      !! #=# availability
      !! #=#=# consistency
      !! Logic change during pre-factoring: err, symbol, units will be assigned
      !! only if present.
      !! cplife : input real : calculated CP full power year lifetime (years)
      !! bktlife : input real : calculated first wall/blanket power year lifetime (years)
      !! divlife : input real : calculated divertor  power year lifetime (years)
      !! i_cp_lifetime : input integer : switch chosing which plant element the CP
      !!                                 the CP lifetime must equate
      use cost_variables, only : cplife, divlife, cplife_input, &
         tlife, i_cp_lifetime
      use fwbs_variables, only : bktlife

      implicit none

            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
      !! Constraints output


      ! The CP lifetime is equal to the the divertor one
      if  ( i_cp_lifetime == 0 ) then
         tmp_cc = 1.0D0 - cplife/cplife_input

      else if ( i_cp_lifetime == 1 ) then
         tmp_cc = 1.0D0 - cplife/divlife

      ! The CP lifetime is equal to the tritium breeding blankets / FW one
      else if ( i_cp_lifetime == 2 ) then
         tmp_cc = 1.0D0 - cplife/bktlife

      ! The CP lifetime is equal to the
      else if ( i_cp_lifetime == 3 ) then
         tmp_cc = 1.0D0 - cplife/tlife
      end if

      tmp_con = divlife * (1.0D0 - tmp_cc)
      tmp_err = divlife * tmp_cc
      tmp_symbol = '='
      tmp_units = 'years'

   end subroutine constraint_eqn_085

   subroutine constraint_eqn_086(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! 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]
      use tfcoil_variables, only : t_turn_tf, f_t_turn_tf, t_turn_tf_max

      implicit none

            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      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_con = t_turn_tf_max * (1.0D0 - tmp_cc)
      tmp_err = t_turn_tf_max * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'm'

   end subroutine constraint_eqn_086


   subroutine constraint_eqn_087(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! 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)

      use heat_transport_variables, only: crypmw, crypmw_max, f_crypmw
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      tmp_cc =  1.0D0 - f_crypmw * crypmw_max/crypmw
      tmp_con = crypmw_max * (1.0D0 - tmp_cc)
      tmp_err = crypmw * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW'
   end subroutine constraint_eqn_087

   subroutine constraint_eqn_088(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! Equation for TF coil vertical strain upper limit (absolute value)
      !! author: CPS Swanson, PPPL, USA
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! Equation for TF coil vertical strain upper limit (absolute value)
      !! #=# tfcoil
      !! #=#=# fstr_wp, str_wp_max
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! fstr_wp : input real : f-value for TF coil strain
      !! str_wp_max : input real : Allowable maximum TF coil vertical strain
      !! str_wp : input real : Constrained TF coil vertical strain
      use constraint_variables, only: fstr_wp
      use tfcoil_variables, only: str_wp_max, str_wp
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      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_con = str_wp_max
      tmp_err = str_wp_max - abs(str_wp) / fstr_wp
      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
      !! author: G Turkington, CCFE, Culham Science Centre
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !! #=# physics
      !! #=#=# f_copperaoh_m2, copperaoh_m2_max
      !! and hence also optional here.
      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
      !! copperaoh_m2 : input real : CS coil current at EOF / copper area [A/m2]
      !! copperaoh_m2_max : input real : maximum coil current / copper area [A/m2]
      !! f_copperaoh_m2 : input real : f-value for CS coil current / copper area
      use rebco_variables, only: copperaoh_m2, copperaoh_m2_max, f_copperaoh_m2
      implicit none
                  real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      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_con = copperaoh_m2
      tmp_err = copperaoh_m2 * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'A/m2'

   end subroutine constraint_eqn_089

   subroutine constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
      !! 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
      use CS_fatigue_variables, only: n_cycle, n_cycle_min, bkt_life_csf
      use constraint_variables, only: fncycle
      use cost_variables, only: ibkt_life, bktcycles
      implicit none
            real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units
      !! Switch to relay the calculated fw/blanket lifetime cycles as the minimum required CS stress cycles.
      !! bkt_life_cycle = 1 turns on the relay. Otherwise the models run independently.
      if (ibkt_life == 1 .and. bkt_life_csf == 1 ) then
         n_cycle_min = bktcycles
      end if

      tmp_cc =  1.0D0 - fncycle * n_cycle / n_cycle_min
      tmp_con = n_cycle_min * (1.0D0 - tmp_cc)
      tmp_err = n_cycle * tmp_cc
      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
      !! 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
      !! args : output structure : residual error; constraint value;
      !! residual error in physical units; output string; units string
      !!  powerht_local > powerscaling
      !! #=# physics
      !! #=#=# fecrh_ignition, powerht_local, powerscaling
      !! fecrh_ignition : input real : f-value for constraint powerht_local > powerscaling
      !! max_gyrotron_frequency : input real :  Max. av. gyrotron frequency
      !! te0_ecrh_achievable : input real : Max. achievable electron temperature at ignition point
      use constraint_variables, only: fecrh_ignition
      use stellarator_variables, only: max_gyrotron_frequency, te0_ecrh_achievable, powerscaling_constraint, powerht_constraint
      use physics_variables, only: ignite
      use current_drive_variables, only: pheat
      implicit none
      real(dp), intent(out) :: tmp_cc
      real(dp), intent(out) :: tmp_con
      real(dp), intent(out) :: tmp_err
      character(len=1), intent(out) :: tmp_symbol
      character(len=10), intent(out) :: tmp_units

      ! Achievable ECRH te needs to be larger than needed te for igntion
      if(ignite==0) then
         tmp_cc = 1.0D0 - fecrh_ignition* (powerht_constraint+pheat)/powerscaling_constraint
      else
         tmp_cc = 1.0D0 - fecrh_ignition* powerht_constraint/powerscaling_constraint
      endif

      tmp_con = powerscaling_constraint * (1.0D0 - tmp_cc)
      tmp_err = powerht_constraint * tmp_cc
      tmp_symbol = '<'
      tmp_units = 'MW'
   end subroutine constraint_eqn_091



end module constraints