C
C There are a total of 31 entries in the algebraic variable array.
C There are a total of 16 entries in each of the rate and state variable arrays.
C There are a total of 72 entries in the constant variable array.
C
C
C VOI is time in component environment (millisecond).
C STATES(1) is V in component membrane (millivolt).
C CONSTS(1) is C in component membrane (uF).
C ALGBRC(1) is i_Na in component sodium_current (nanoA).
C ALGBRC(16) is i_NaP in component persistent_sodium_current (nanoA).
C ALGBRC(23) is i_K in component delayed_rectifier_current (nanoA).
C ALGBRC(24) is i_leak in component leak_current (nanoA).
C ALGBRC(25) is i_T in component LVA_calcium_current (nanoA).
C ALGBRC(26) is i_N in component N_HVA_calcium_current (nanoA).
C ALGBRC(27) is i_P in component P_HVA_calcium_current (nanoA).
C ALGBRC(28) is i_SK in component calcium_dependent_potassium_current (nanoA).
C ALGBRC(29) is i_A in component fast_transient_potassium_current (nanoA).
C ALGBRC(30) is i_H in component hyperpolarization_activated_current (nanoA).
C ALGBRC(31) is i_app in component stimulus_protocol (nanoA).
C CONSTS(2) is g_Na in component sodium_current (uS).
C CONSTS(3) is E_Na in component sodium_current (millivolt).
C STATES(2) is m in component sodium_current_m_gate (dimensionless).
C STATES(3) is h in component sodium_current_h_gate (dimensionless).
C CONSTS(4) is theta_h in component sodium_current_h_gate (millivolt).
C CONSTS(5) is sigma_h in component sodium_current_h_gate (millivolt).
C CONSTS(6) is theta_1 in component sodium_current_h_gate (millivolt).
C CONSTS(7) is sigma_1 in component sodium_current_h_gate (millivolt).
C CONSTS(8) is sigma_2 in component sodium_current_h_gate (millivolt).
C ALGBRC(2) is tau_h in component sodium_current_h_gate (millisecond).
C ALGBRC(17) is h_infinity in component sodium_current_h_gate (dimensionless).
C CONSTS(9) is theta_m in component sodium_current_m_gate (millivolt).
C CONSTS(10) is sigma_m in component sodium_current_m_gate (millivolt).
C CONSTS(11) is tau_m in component sodium_current_m_gate (millisecond).
C ALGBRC(3) is m_infinity in component sodium_current_m_gate (dimensionless).
C CONSTS(12) is g_NaP in component persistent_sodium_current (uS).
C STATES(4) is m in component persistent_sodium_current_m_gate (dimensionless).
C STATES(5) is h in component persistent_sodium_current_h_gate (dimensionless).
C CONSTS(13) is theta_h in component persistent_sodium_current_h_gate (millivolt).
C CONSTS(14) is sigma_h in component persistent_sodium_current_h_gate (millivolt).
C CONSTS(15) is tau_h in component persistent_sodium_current_h_gate (millisecond).
C ALGBRC(4) is h_infinity in component persistent_sodium_current_h_gate (dimensionless).
C CONSTS(16) is theta_m in component persistent_sodium_current_m_gate (millivolt).
C CONSTS(17) is sigma_m in component persistent_sodium_current_m_gate (millivolt).
C CONSTS(18) is tau_m in component persistent_sodium_current_m_gate (millisecond).
C ALGBRC(5) is m_infinity in component persistent_sodium_current_m_gate (dimensionless).
C CONSTS(19) is g_K in component delayed_rectifier_current (uS).
C CONSTS(20) is E_K in component delayed_rectifier_current (millivolt).
C STATES(6) is n in component delayed_rectifier_current_n_gate (dimensionless).
C CONSTS(21) is theta_n in component delayed_rectifier_current_n_gate (millivolt).
C CONSTS(22) is sigma_n in component delayed_rectifier_current_n_gate (millivolt).
C CONSTS(23) is theta_1 in component delayed_rectifier_current_n_gate (millivolt).
C CONSTS(24) is sigma_1 in component delayed_rectifier_current_n_gate (millivolt).
C CONSTS(25) is sigma_2 in component delayed_rectifier_current_n_gate (millivolt).
C ALGBRC(6) is tau_n in component delayed_rectifier_current_n_gate (millisecond).
C ALGBRC(18) is n_infinity in component delayed_rectifier_current_n_gate (dimensionless).
C CONSTS(26) is g_leak in component leak_current (uS).
C CONSTS(27) is E_leak in component leak_current (millivolt).
C CONSTS(28) is g_T in component LVA_calcium_current (uS).
C CONSTS(29) is E_Ca in component LVA_calcium_current (millivolt).
C STATES(7) is m in component LVA_calcium_current_m_gate (dimensionless).
C STATES(8) is h in component LVA_calcium_current_h_gate (dimensionless).
C CONSTS(30) is theta_m in component LVA_calcium_current_m_gate (millivolt).
C CONSTS(31) is sigma_m in component LVA_calcium_current_m_gate (millivolt).
C CONSTS(32) is theta_1 in component LVA_calcium_current_m_gate (millivolt).
C CONSTS(33) is sigma_1 in component LVA_calcium_current_m_gate (millivolt).
C CONSTS(34) is sigma_2 in component LVA_calcium_current_m_gate (millivolt).
C ALGBRC(7) is tau_m in component LVA_calcium_current_m_gate (millisecond).
C ALGBRC(19) is m_infinity in component LVA_calcium_current_m_gate (dimensionless).
C CONSTS(35) is theta_h in component LVA_calcium_current_h_gate (millivolt).
C CONSTS(36) is sigma_h in component LVA_calcium_current_h_gate (millivolt).
C CONSTS(37) is theta_1 in component LVA_calcium_current_h_gate (millivolt).
C CONSTS(38) is sigma_1 in component LVA_calcium_current_h_gate (millivolt).
C ALGBRC(8) is tau_h in component LVA_calcium_current_h_gate (millisecond).
C ALGBRC(20) is h_infinity in component LVA_calcium_current_h_gate (dimensionless).
C CONSTS(39) is g_N in component N_HVA_calcium_current (uS).
C STATES(9) is m in component N_HVA_calcium_current_m_gate (dimensionless).
C STATES(10) is h in component N_HVA_calcium_current_h_gate (dimensionless).
C CONSTS(40) is theta_m in component N_HVA_calcium_current_m_gate (millivolt).
C CONSTS(41) is sigma_m in component N_HVA_calcium_current_m_gate (millivolt).
C CONSTS(42) is tau_m in component N_HVA_calcium_current_m_gate (millisecond).
C ALGBRC(9) is m_infinity in component N_HVA_calcium_current_m_gate (dimensionless).
C CONSTS(43) is theta_h in component N_HVA_calcium_current_h_gate (millivolt).
C CONSTS(44) is sigma_h in component N_HVA_calcium_current_h_gate (millivolt).
C CONSTS(45) is tau_h in component N_HVA_calcium_current_h_gate (millisecond).
C ALGBRC(10) is h_infinity in component N_HVA_calcium_current_h_gate (dimensionless).
C CONSTS(46) is g_P in component P_HVA_calcium_current (uS).
C STATES(11) is m in component P_HVA_calcium_current_m_gate (dimensionless).
C CONSTS(47) is theta_m in component P_HVA_calcium_current_m_gate (millivolt).
C CONSTS(48) is sigma_m in component P_HVA_calcium_current_m_gate (millivolt).
C CONSTS(49) is tau_m in component P_HVA_calcium_current_m_gate (millisecond).
C ALGBRC(11) is m_infinity in component P_HVA_calcium_current_m_gate (dimensionless).
C CONSTS(50) is g_SK in component calcium_dependent_potassium_current (uS).
C STATES(12) is z in component calcium_dependent_potassium_current_z_gate (dimensionless).
C CONSTS(51) is K1 in component calcium_dependent_potassium_current_z_gate (uM_per_nanocoulomb).
C CONSTS(52) is K2 in component calcium_dependent_potassium_current_z_gate (per_ms).
C STATES(13) is Ca_conc in component calcium_dependent_potassium_current_z_gate (uM).
C CONSTS(53) is tau_z in component calcium_dependent_potassium_current_z_gate (millisecond).
C ALGBRC(12) is z_infinity in component calcium_dependent_potassium_current_z_gate (dimensionless).
C CONSTS(54) is g_A in component fast_transient_potassium_current (uS).
C STATES(14) is m in component fast_transient_potassium_current_m_gate (dimensionless).
C STATES(15) is h in component fast_transient_potassium_current_h_gate (dimensionless).
C CONSTS(55) is theta_m in component fast_transient_potassium_current_m_gate (millivolt).
C CONSTS(56) is sigma_m in component fast_transient_potassium_current_m_gate (millivolt).
C CONSTS(57) is theta_1 in component fast_transient_potassium_current_m_gate (millivolt).
C CONSTS(58) is theta_2 in component fast_transient_potassium_current_m_gate (millivolt).
C CONSTS(59) is sigma_1 in component fast_transient_potassium_current_m_gate (millivolt).
C CONSTS(60) is sigma_2 in component fast_transient_potassium_current_m_gate (millivolt).
C ALGBRC(13) is tau_m in component fast_transient_potassium_current_m_gate (millisecond).
C ALGBRC(21) is m_infinity in component fast_transient_potassium_current_m_gate (dimensionless).
C CONSTS(61) is theta_h in component fast_transient_potassium_current_h_gate (millivolt).
C CONSTS(62) is sigma_h in component fast_transient_potassium_current_h_gate (millivolt).
C CONSTS(63) is tau_h in component fast_transient_potassium_current_h_gate (millisecond).
C ALGBRC(14) is h_infinity in component fast_transient_potassium_current_h_gate (dimensionless).
C CONSTS(64) is g_H in component hyperpolarization_activated_current (uS).
C CONSTS(65) is E_H in component hyperpolarization_activated_current (millivolt).
C STATES(16) is m in component hyperpolarization_activated_current_m_gate (dimensionless).
C CONSTS(66) is theta_m in component hyperpolarization_activated_current_m_gate (millivolt).
C CONSTS(67) is sigma_m in component hyperpolarization_activated_current_m_gate (millivolt).
C CONSTS(68) is theta_1 in component hyperpolarization_activated_current_m_gate (millivolt).
C CONSTS(69) is sigma_1 in component hyperpolarization_activated_current_m_gate (millivolt).
C ALGBRC(15) is tau_m in component hyperpolarization_activated_current_m_gate (millisecond).
C ALGBRC(22) is m_infinity in component hyperpolarization_activated_current_m_gate (dimensionless).
C CONSTS(70) is i_stimStart in component stimulus_protocol (millisecond).
C CONSTS(71) is i_stimEnd in component stimulus_protocol (millisecond).
C CONSTS(72) is i_stimAmplitude in component stimulus_protocol (nanoA).
C RATES(1) is d/dt V in component membrane (millivolt).
C RATES(3) is d/dt h in component sodium_current_h_gate (dimensionless).
C RATES(2) is d/dt m in component sodium_current_m_gate (dimensionless).
C RATES(5) is d/dt h in component persistent_sodium_current_h_gate (dimensionless).
C RATES(4) is d/dt m in component persistent_sodium_current_m_gate (dimensionless).
C RATES(6) is d/dt n in component delayed_rectifier_current_n_gate (dimensionless).
C RATES(7) is d/dt m in component LVA_calcium_current_m_gate (dimensionless).
C RATES(8) is d/dt h in component LVA_calcium_current_h_gate (dimensionless).
C RATES(9) is d/dt m in component N_HVA_calcium_current_m_gate (dimensionless).
C RATES(10) is d/dt h in component N_HVA_calcium_current_h_gate (dimensionless).
C RATES(11) is d/dt m in component P_HVA_calcium_current_m_gate (dimensionless).
C RATES(13) is d/dt Ca_conc in component calcium_dependent_potassium_current_z_gate (uM).
C RATES(12) is d/dt z in component calcium_dependent_potassium_current_z_gate (dimensionless).
C RATES(14) is d/dt m in component fast_transient_potassium_current_m_gate (dimensionless).
C RATES(15) is d/dt h in component fast_transient_potassium_current_h_gate (dimensionless).
C RATES(16) is d/dt m in component hyperpolarization_activated_current_m_gate (dimensionless).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = -71.847
      CONSTS(1) = 0.04
      CONSTS(2) = 0.7
      CONSTS(3) = 60
      STATES(2) = 0.015
      STATES(3) = 0.981
      CONSTS(4) = 44.1
      CONSTS(5) = 7
      CONSTS(6) = 35
      CONSTS(7) = 4
      CONSTS(8) = 25
      CONSTS(9) = 36
      CONSTS(10) = 8.5
      CONSTS(11) = 0.1
      CONSTS(12) = 0.05
      STATES(4) = 0.002
      STATES(5) = 0.797
      CONSTS(13) = 65
      CONSTS(14) = 5
      CONSTS(15) = 150
      CONSTS(16) = 47.1
      CONSTS(17) = 4.1
      CONSTS(18) = 0.1
      CONSTS(19) = 1.3
      CONSTS(20) = -80
      STATES(6) = 0.158
      CONSTS(21) = 30
      CONSTS(22) = 25
      CONSTS(23) = 30
      CONSTS(24) = 40
      CONSTS(25) = 50
      CONSTS(26) = 0.005
      CONSTS(27) = -50
      CONSTS(28) = 0.1
      CONSTS(29) = 40
      STATES(7) = 0.001
      STATES(8) = 0.562
      CONSTS(30) = 38
      CONSTS(31) = 5
      CONSTS(32) = 28
      CONSTS(33) = 25
      CONSTS(34) = 70
      CONSTS(35) = 70.1
      CONSTS(36) = 7
      CONSTS(37) = 70
      CONSTS(38) = 65
      CONSTS(39) = 0.05
      STATES(9) = 0.001
      STATES(10) = 0.649
      CONSTS(40) = 30
      CONSTS(41) = 6
      CONSTS(42) = 5
      CONSTS(43) = 70
      CONSTS(44) = 3
      CONSTS(45) = 25
      CONSTS(46) = 0.05
      STATES(11) = 0
      CONSTS(47) = 17
      CONSTS(48) = 3
      CONSTS(49) = 10
      CONSTS(50) = 0.3
      STATES(12) = 0
      CONSTS(51) = -500
      CONSTS(52) = 0.04
      STATES(13) = 0.0604
      CONSTS(53) = 1
      CONSTS(54) = 1
      STATES(14) = 0.057
      STATES(15) = 0.287
      CONSTS(55) = 27
      CONSTS(56) = 16
      CONSTS(57) = 40
      CONSTS(58) = 74
      CONSTS(59) = 5
      CONSTS(60) = 7.5
      CONSTS(61) = 80
      CONSTS(62) = 11
      CONSTS(63) = 20
      CONSTS(64) = 0.005
      CONSTS(65) = -38.8
      STATES(16) = 0.182
      CONSTS(66) = 79.8
      CONSTS(67) = 5.3
      CONSTS(68) = 70
      CONSTS(69) = 11
      CONSTS(70) = 10
      CONSTS(71) = 11
      CONSTS(72) = 10
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(3) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(9))/- CONSTS(10)))
      RATES(2) = (ALGBRC(3) - STATES(2))/CONSTS(11)
      ALGBRC(4) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(13))/CONSTS(14)))
      RATES(5) = (ALGBRC(4) - STATES(5))/CONSTS(15)
      ALGBRC(5) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(16))/CONSTS(17)))
      RATES(4) = (ALGBRC(5) - STATES(4))/CONSTS(18)
      ALGBRC(9) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(40))/CONSTS(41)))
      RATES(9) = (ALGBRC(9) - STATES(9))/CONSTS(42)
      ALGBRC(10) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(43))/CONSTS(44)))
      RATES(10) = (ALGBRC(10) - STATES(10))/CONSTS(45)
      ALGBRC(11) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(47))/CONSTS(48)))
      RATES(11) = (ALGBRC(11) - STATES(11))/CONSTS(49)
      ALGBRC(12) = 1.00000/(1.00000+0.00300000/STATES(13) ** 2.00000)
      RATES(12) = (ALGBRC(12) - STATES(12))/CONSTS(53)
      ALGBRC(14) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(61))/CONSTS(62)))
      RATES(15) = (ALGBRC(14) - STATES(15))/CONSTS(63)
      ALGBRC(2) = 3.50000/(EXP((STATES(1)+CONSTS(6))/CONSTS(7))+EXP(- (STATES(1)+CONSTS(6))/CONSTS(8)))+1.00000
      ALGBRC(17) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(4))/CONSTS(5)))
      RATES(3) = (ALGBRC(17) - STATES(3))/ALGBRC(2)
      ALGBRC(6) = 2.50000/(EXP((STATES(1)+CONSTS(23))/CONSTS(24))+EXP(- (STATES(1)+CONSTS(23))/CONSTS(25)))+0.0100000
      ALGBRC(18) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(21))/CONSTS(22)))
      RATES(6) = (ALGBRC(18) - STATES(6))/ALGBRC(6)
      ALGBRC(7) = 5.00000/(EXP((STATES(1)+CONSTS(32))/CONSTS(33))+EXP(- (STATES(1)+CONSTS(32))/CONSTS(34)))+2.00000
      ALGBRC(19) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(30))/CONSTS(31)))
      RATES(7) = (ALGBRC(19) - STATES(7))/ALGBRC(7)
      ALGBRC(8) = 20.0000/(EXP((STATES(1)+CONSTS(37))/CONSTS(38))+EXP(- (STATES(1)+CONSTS(37))/CONSTS(38)))+1.00000
      ALGBRC(20) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(35))/CONSTS(36)))
      RATES(8) = (ALGBRC(20) - STATES(8))/ALGBRC(8)
      ALGBRC(13) = 1.00000/(EXP((STATES(1)+CONSTS(57))/CONSTS(59))+EXP(- (STATES(1)+CONSTS(58))/CONSTS(60)))+0.370000
      ALGBRC(21) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(55))/CONSTS(56)))
      RATES(14) = (ALGBRC(21) - STATES(14))/ALGBRC(13)
      ALGBRC(15) = 1.00000/(EXP((STATES(1)+CONSTS(68))/CONSTS(69))+EXP(- (STATES(1)+CONSTS(68))/CONSTS(69)))+50.0000
      ALGBRC(22) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(66))/CONSTS(67)))
      RATES(16) = (ALGBRC(22) - STATES(16))/ALGBRC(15)
      ALGBRC(25) =  CONSTS(28)*STATES(7)*STATES(8)*(STATES(1) - CONSTS(29))
      ALGBRC(26) =  CONSTS(39)*STATES(9)*STATES(10)*(STATES(1) - CONSTS(29))
      ALGBRC(27) =  CONSTS(46)*STATES(11)*(STATES(1) - CONSTS(29))
      RATES(13) =  (1.00000/1000.00)*CONSTS(51)*(ALGBRC(25)+ALGBRC(26)+ALGBRC(27)) -  CONSTS(52)*STATES(13)
      ALGBRC(1) =  CONSTS(2)*STATES(2) ** 3.00000*STATES(3)*(STATES(1) - CONSTS(3))
      ALGBRC(16) =  CONSTS(12)*STATES(4)*STATES(5)*(STATES(1) - CONSTS(3))
      ALGBRC(23) =  CONSTS(19)*STATES(6) ** 4.00000*(STATES(1) - CONSTS(20))
      ALGBRC(24) =  CONSTS(26)*(STATES(1) - CONSTS(27))
      ALGBRC(28) =  CONSTS(50)*STATES(12) ** 2.00000*(STATES(1) - CONSTS(20))
      ALGBRC(29) =  CONSTS(54)*STATES(14)*STATES(15)*(STATES(1) - CONSTS(20))
      ALGBRC(30) =  CONSTS(64)*STATES(16)*(STATES(1) - CONSTS(65))
      ALGBRC(31) = TERNRY(VOI.GE.CONSTS(70).AND.VOI.LE.CONSTS(71), CONSTS(72), 0.00000)
      RATES(1) = (- (ALGBRC(1)+ALGBRC(16)+ALGBRC(23)+ALGBRC(24)+ALGBRC(25)+ALGBRC(26)+ALGBRC(27)+ALGBRC(28)+ALGBRC(29)+ALGBRC(30))+ALGBRC(31))/CONSTS(1)
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(3) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(9))/- CONSTS(10)))
      ALGBRC(4) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(13))/CONSTS(14)))
      ALGBRC(5) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(16))/CONSTS(17)))
      ALGBRC(9) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(40))/CONSTS(41)))
      ALGBRC(10) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(43))/CONSTS(44)))
      ALGBRC(11) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(47))/CONSTS(48)))
      ALGBRC(12) = 1.00000/(1.00000+0.00300000/STATES(13) ** 2.00000)
      ALGBRC(14) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(61))/CONSTS(62)))
      ALGBRC(2) = 3.50000/(EXP((STATES(1)+CONSTS(6))/CONSTS(7))+EXP(- (STATES(1)+CONSTS(6))/CONSTS(8)))+1.00000
      ALGBRC(17) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(4))/CONSTS(5)))
      ALGBRC(6) = 2.50000/(EXP((STATES(1)+CONSTS(23))/CONSTS(24))+EXP(- (STATES(1)+CONSTS(23))/CONSTS(25)))+0.0100000
      ALGBRC(18) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(21))/CONSTS(22)))
      ALGBRC(7) = 5.00000/(EXP((STATES(1)+CONSTS(32))/CONSTS(33))+EXP(- (STATES(1)+CONSTS(32))/CONSTS(34)))+2.00000
      ALGBRC(19) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(30))/CONSTS(31)))
      ALGBRC(8) = 20.0000/(EXP((STATES(1)+CONSTS(37))/CONSTS(38))+EXP(- (STATES(1)+CONSTS(37))/CONSTS(38)))+1.00000
      ALGBRC(20) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(35))/CONSTS(36)))
      ALGBRC(13) = 1.00000/(EXP((STATES(1)+CONSTS(57))/CONSTS(59))+EXP(- (STATES(1)+CONSTS(58))/CONSTS(60)))+0.370000
      ALGBRC(21) = 1.00000/(1.00000+EXP(- (STATES(1)+CONSTS(55))/CONSTS(56)))
      ALGBRC(15) = 1.00000/(EXP((STATES(1)+CONSTS(68))/CONSTS(69))+EXP(- (STATES(1)+CONSTS(68))/CONSTS(69)))+50.0000
      ALGBRC(22) = 1.00000/(1.00000+EXP((STATES(1)+CONSTS(66))/CONSTS(67)))
      ALGBRC(25) =  CONSTS(28)*STATES(7)*STATES(8)*(STATES(1) - CONSTS(29))
      ALGBRC(26) =  CONSTS(39)*STATES(9)*STATES(10)*(STATES(1) - CONSTS(29))
      ALGBRC(27) =  CONSTS(46)*STATES(11)*(STATES(1) - CONSTS(29))
      ALGBRC(1) =  CONSTS(2)*STATES(2) ** 3.00000*STATES(3)*(STATES(1) - CONSTS(3))
      ALGBRC(16) =  CONSTS(12)*STATES(4)*STATES(5)*(STATES(1) - CONSTS(3))
      ALGBRC(23) =  CONSTS(19)*STATES(6) ** 4.00000*(STATES(1) - CONSTS(20))
      ALGBRC(24) =  CONSTS(26)*(STATES(1) - CONSTS(27))
      ALGBRC(28) =  CONSTS(50)*STATES(12) ** 2.00000*(STATES(1) - CONSTS(20))
      ALGBRC(29) =  CONSTS(54)*STATES(14)*STATES(15)*(STATES(1) - CONSTS(20))
      ALGBRC(30) =  CONSTS(64)*STATES(16)*(STATES(1) - CONSTS(65))
      ALGBRC(31) = TERNRY(VOI.GE.CONSTS(70).AND.VOI.LE.CONSTS(71), CONSTS(72), 0.00000)
      RETURN
      END
      REAL FUNCTION TERNRY(TEST, VALA, VALB)
      LOGICAL TEST
      REAL VALA, VALB
      IF (TEST) THEN
        TERNRY = VALA
      ELSE
        TERNRY = VALB
      ENDIF
      RETURN
      END