C
C There are a total of 47 entries in the algebraic variable array.
C There are a total of 13 entries in each of the rate and state variable arrays.
C There are a total of 55 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 R in component membrane (joule_per_mole_kelvin).
C CONSTS(2) is T in component membrane (kelvin).
C CONSTS(3) is F in component membrane (coulomb_per_millimole).
C ALGBRC(21) is i_Na in component fast_sodium_current (microA_per_microF).
C ALGBRC(40) is i_Ca in component L_type_Ca_current (microA_per_microF).
C ALGBRC(41) is i_CaK in component L_type_Ca_current (microA_per_microF).
C ALGBRC(27) is i_Kr in component rapid_activating_delayed_rectifiyer_K_current (microA_per_microF).
C ALGBRC(28) is i_Ks in component slow_activating_delayed_rectifiyer_K_current (microA_per_microF).
C ALGBRC(29) is i_to in component transient_outward_potassium_current (microA_per_microF).
C ALGBRC(25) is i_K1 in component time_independent_potassium_current (microA_per_microF).
C ALGBRC(31) is i_Kp in component plateau_potassium_current (microA_per_microF).
C ALGBRC(34) is i_NaCa in component Na_Ca_exchanger (microA_per_microF).
C ALGBRC(33) is i_NaK in component sodium_potassium_pump (microA_per_microF).
C ALGBRC(35) is i_p_Ca in component sarcolemmal_calcium_pump (microA_per_microF).
C ALGBRC(37) is i_Ca_b in component calcium_background_current (microA_per_microF).
C ALGBRC(38) is i_Na_b in component sodium_background_current (microA_per_microF).
C ALGBRC(11) is i_Stim in component membrane (microA_per_microF).
C CONSTS(4) is stim_start in component membrane (millisecond).
C CONSTS(5) is stim_end in component membrane (millisecond).
C CONSTS(6) is stim_period in component membrane (millisecond).
C CONSTS(7) is stim_duration in component membrane (millisecond).
C CONSTS(8) is stim_amplitude in component membrane (microA_per_microF).
C CONSTS(51) is E_Na in component fast_sodium_current (millivolt).
C CONSTS(9) is g_Na in component fast_sodium_current (milliS_per_microF).
C CONSTS(10) is Na_o in component standard_ionic_concentrations (millimolar).
C CONSTS(11) is Na_i in component standard_ionic_concentrations (millimolar).
C STATES(2) is m in component fast_sodium_current_m_gate (dimensionless).
C STATES(3) is h in component fast_sodium_current_h_gate (dimensionless).
C STATES(4) is j in component fast_sodium_current_j_gate (dimensionless).
C ALGBRC(12) is alpha_m in component fast_sodium_current_m_gate (per_millisecond).
C ALGBRC(22) is beta_m in component fast_sodium_current_m_gate (per_millisecond).
C ALGBRC(1) is E0_m in component fast_sodium_current_m_gate (millivolt).
C ALGBRC(2) is alpha_h in component fast_sodium_current_h_gate (per_millisecond).
C ALGBRC(13) is beta_h in component fast_sodium_current_h_gate (per_millisecond).
C CONSTS(12) is shift_h in component fast_sodium_current_h_gate (millivolt).
C ALGBRC(3) is alpha_j in component fast_sodium_current_j_gate (per_millisecond).
C ALGBRC(14) is beta_j in component fast_sodium_current_j_gate (per_millisecond).
C CONSTS(13) is shift_j in component fast_sodium_current_j_gate (millivolt).
C CONSTS(14) is g_K1 in component time_independent_potassium_current (milliS_per_microF).
C CONSTS(15) is K_mK1 in component time_independent_potassium_current (millimolar).
C CONSTS(52) is E_K in component rapid_activating_delayed_rectifiyer_K_current (millivolt).
C CONSTS(16) is K_o in component standard_ionic_concentrations (millimolar).
C ALGBRC(24) is K1_infinity in component time_independent_potassium_current_K1_gate (dimensionless).
C CONSTS(17) is g_Kr in component rapid_activating_delayed_rectifiyer_K_current (milliS_per_microF).
C ALGBRC(26) is R_V in component rapid_activating_delayed_rectifiyer_K_current (dimensionless).
C CONSTS(18) is K_i in component standard_ionic_concentrations (millimolar).
C STATES(5) is X_kr in component rapid_activating_delayed_rectifiyer_K_current_X_kr_gate (dimensionless).
C ALGBRC(4) is X_kr_inf in component rapid_activating_delayed_rectifiyer_K_current_X_kr_gate (dimensionless).
C ALGBRC(15) is tau_X_kr in component rapid_activating_delayed_rectifiyer_K_current_X_kr_gate (millisecond).
C CONSTS(19) is g_Ks in component slow_activating_delayed_rectifiyer_K_current (milliS_per_microF).
C CONSTS(53) is E_Ks in component slow_activating_delayed_rectifiyer_K_current (millivolt).
C STATES(6) is X_ks in component slow_activating_delayed_rectifiyer_K_current_X_ks_gate (dimensionless).
C ALGBRC(16) is tau_X_ks in component slow_activating_delayed_rectifiyer_K_current_X_ks_gate (millisecond).
C ALGBRC(5) is X_ks_infinity in component slow_activating_delayed_rectifiyer_K_current_X_ks_gate (dimensionless).
C CONSTS(20) is g_to in component transient_outward_potassium_current (milliS_per_microF).
C STATES(7) is X_to in component transient_outward_potassium_current_X_to_gate (dimensionless).
C STATES(8) is Y_to in component transient_outward_potassium_current_Y_to_gate (dimensionless).
C ALGBRC(6) is alpha_X_to in component transient_outward_potassium_current_X_to_gate (per_millisecond).
C ALGBRC(17) is beta_X_to in component transient_outward_potassium_current_X_to_gate (per_millisecond).
C ALGBRC(7) is alpha_Y_to in component transient_outward_potassium_current_Y_to_gate (per_millisecond).
C ALGBRC(18) is beta_Y_to in component transient_outward_potassium_current_Y_to_gate (per_millisecond).
C CONSTS(21) is g_Kp in component plateau_potassium_current (milliS_per_microF).
C ALGBRC(30) is Kp_V in component plateau_potassium_current_Kp_gate (dimensionless).
C CONSTS(22) is i_NaK_max in component sodium_potassium_pump (microA_per_microF).
C ALGBRC(32) is f_NaK in component sodium_potassium_pump (dimensionless).
C CONSTS(23) is K_mNai in component sodium_potassium_pump (millimolar).
C CONSTS(24) is K_mKo in component sodium_potassium_pump (millimolar).
C CONSTS(54) is sigma in component sodium_potassium_pump (dimensionless).
C CONSTS(25) is K_mCa in component Na_Ca_exchanger (micromolar).
C CONSTS(26) is K_mNa in component Na_Ca_exchanger (millimolar).
C CONSTS(27) is K_NaCa in component Na_Ca_exchanger (microA_per_microF).
C CONSTS(28) is K_sat in component Na_Ca_exchanger (dimensionless).
C CONSTS(29) is eta in component Na_Ca_exchanger (dimensionless).
C STATES(9) is Ca_i in component calcium_dynamics (micromolar).
C CONSTS(30) is Ca_o in component standard_ionic_concentrations (micromolar).
C CONSTS(31) is K_mpCa in component sarcolemmal_calcium_pump (micromolar).
C CONSTS(32) is i_pCa_max in component sarcolemmal_calcium_pump (microA_per_microF).
C CONSTS(33) is g_Cab in component calcium_background_current (milliS_per_microF).
C ALGBRC(36) is E_Ca in component calcium_background_current (millivolt).
C CONSTS(34) is g_Nab in component sodium_background_current (milliS_per_microF).
C CONSTS(35) is P_Ca in component L_type_Ca_current (cm_per_millisecond).
C CONSTS(36) is P_CaK in component L_type_Ca_current (cm_per_millisecond).
C CONSTS(37) is i_Ca_half in component L_type_Ca_current (microA_per_microF).
C ALGBRC(39) is i_Ca_max in component L_type_Ca_current (microA_per_microF).
C CONSTS(38) is C_sc in component L_type_Ca_current (microF_per_cm2).
C STATES(10) is f in component L_type_Ca_current_f_gate (dimensionless).
C STATES(11) is d in component L_type_Ca_current_d_gate (dimensionless).
C STATES(12) is f_Ca in component L_type_Ca_current_f_Ca_gate (dimensionless).
C ALGBRC(8) is f_infinity in component L_type_Ca_current_f_gate (dimensionless).
C ALGBRC(19) is tau_f in component L_type_Ca_current_f_gate (millisecond).
C ALGBRC(9) is d_infinity in component L_type_Ca_current_d_gate (dimensionless).
C ALGBRC(23) is tau_d in component L_type_Ca_current_d_gate (millisecond).
C ALGBRC(20) is E0_m in component L_type_Ca_current_d_gate (millivolt).
C CONSTS(55) is tau_f_Ca in component L_type_Ca_current_f_Ca_gate (millisecond).
C ALGBRC(10) is f_Ca_infinity in component L_type_Ca_current_f_Ca_gate (dimensionless).
C CONSTS(39) is K_mfCa in component L_type_Ca_current_f_Ca_gate (micromolar).
C ALGBRC(46) is beta_i in component calcium_dynamics (dimensionless).
C CONSTS(40) is K_mCMDN in component calcium_dynamics (micromolar).
C CONSTS(41) is CMDN_tot in component calcium_dynamics (micromolar).
C CONSTS(42) is V_myo in component calcium_dynamics (microlitre).
C CONSTS(43) is A_Cap in component calcium_dynamics (cm2).
C ALGBRC(44) is J_rel in component calcium_dynamics (micromolar_per_millisecond).
C ALGBRC(45) is J_leak in component calcium_dynamics (micromolar_per_millisecond).
C ALGBRC(42) is J_up in component calcium_dynamics (micromolar_per_millisecond).
C STATES(13) is Ca_SR in component calcium_dynamics (micromolar).
C CONSTS(44) is P_rel in component calcium_dynamics (per_millisecond).
C CONSTS(45) is P_leak in component calcium_dynamics (per_millisecond).
C CONSTS(46) is K_mCSQN in component calcium_dynamics (micromolar).
C CONSTS(47) is CSQN_tot in component calcium_dynamics (micromolar).
C CONSTS(48) is V_SR in component calcium_dynamics (microlitre).
C CONSTS(49) is V_up in component calcium_dynamics (micromolar_per_millisecond).
C CONSTS(50) is K_mup in component calcium_dynamics (micromolar).
C ALGBRC(43) is gamma in component calcium_dynamics (dimensionless).
C ALGBRC(47) is beta_SR in component calcium_dynamics (dimensionless).
C RATES(1) is d/dt V in component membrane (millivolt).
C RATES(2) is d/dt m in component fast_sodium_current_m_gate (dimensionless).
C RATES(3) is d/dt h in component fast_sodium_current_h_gate (dimensionless).
C RATES(4) is d/dt j in component fast_sodium_current_j_gate (dimensionless).
C RATES(5) is d/dt X_kr in component rapid_activating_delayed_rectifiyer_K_current_X_kr_gate (dimensionless).
C RATES(6) is d/dt X_ks in component slow_activating_delayed_rectifiyer_K_current_X_ks_gate (dimensionless).
C RATES(7) is d/dt X_to in component transient_outward_potassium_current_X_to_gate (dimensionless).
C RATES(8) is d/dt Y_to in component transient_outward_potassium_current_Y_to_gate (dimensionless).
C RATES(10) is d/dt f in component L_type_Ca_current_f_gate (dimensionless).
C RATES(11) is d/dt d in component L_type_Ca_current_d_gate (dimensionless).
C RATES(12) is d/dt f_Ca in component L_type_Ca_current_f_Ca_gate (dimensionless).
C RATES(13) is d/dt Ca_SR in component calcium_dynamics (micromolar).
C RATES(9) is d/dt Ca_i in component calcium_dynamics (micromolar).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = -94.7
      CONSTS(1) = 8.314
      CONSTS(2) = 310
      CONSTS(3) = 96.5
      CONSTS(4) = 50
      CONSTS(5) = 9000
      CONSTS(6) = 1000
      CONSTS(7) = 1
      CONSTS(8) = -80
      CONSTS(9) = 12.8
      CONSTS(10) = 138
      CONSTS(11) = 10
      STATES(2) = 0.00024676
      STATES(3) = 0.99869
      STATES(4) = 0.99887
      CONSTS(12) = 0
      CONSTS(13) = 0
      CONSTS(14) = 2.8
      CONSTS(15) = 13
      CONSTS(16) = 4
      CONSTS(17) = 0.0136
      CONSTS(18) = 149.4
      STATES(5) = 0.229
      CONSTS(19) = 0.0245
      STATES(6) = 0.0001
      CONSTS(20) = 0.23815
      STATES(7) = 0.00003742
      STATES(8) = 1
      CONSTS(21) = 0.002216
      CONSTS(22) = 0.693
      CONSTS(23) = 10
      CONSTS(24) = 1.5
      CONSTS(25) = 1380
      CONSTS(26) = 87.5
      CONSTS(27) = 1500
      CONSTS(28) = 0.2
      CONSTS(29) = 0.35
      STATES(9) = 0.0472
      CONSTS(30) = 2000
      CONSTS(31) = 0.05
      CONSTS(32) = 0.05
      CONSTS(33) = 0.0003842
      CONSTS(34) = 0.0031
      CONSTS(35) = 0.0000226
      CONSTS(36) = 0.000000579
      CONSTS(37) = -0.265
      CONSTS(38) = 1
      STATES(10) = 0.983
      STATES(11) = 0.0001
      STATES(12) = 0.942
      CONSTS(39) = 0.18
      CONSTS(40) = 2
      CONSTS(41) = 10
      CONSTS(42) = 0.00002584
      CONSTS(43) = 0.0001534
      STATES(13) = 320
      CONSTS(44) = 6
      CONSTS(45) = 0.000001
      CONSTS(46) = 600
      CONSTS(47) = 10000
      CONSTS(48) = 0.000002
      CONSTS(49) = 0.1
      CONSTS(50) = 0.32
      CONSTS(51) =  (( CONSTS(1)*CONSTS(2))/CONSTS(3))*log(CONSTS(10)/CONSTS(11))
      CONSTS(52) =  (( CONSTS(1)*CONSTS(2))/CONSTS(3))*log(CONSTS(16)/CONSTS(18))
      CONSTS(53) =  (( CONSTS(1)*CONSTS(2))/CONSTS(3))*log((CONSTS(16)+ 0.0183300*CONSTS(10))/(CONSTS(18)+ 0.0183300*CONSTS(11)))
      CONSTS(54) =  (1.00000/7.00000)*(EXP(CONSTS(10)/67.3000) - 1.00000)
      CONSTS(55) = 30.0000
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(10) = 1.00000/(1.00000+STATES(9)/CONSTS(39) ** 3.00000)
      RATES(12) = (ALGBRC(10) - STATES(12))/CONSTS(55)
      ALGBRC(2) =  0.135000*EXP(((STATES(1)+80.0000) - CONSTS(12))/- 6.80000)
      ALGBRC(13) = 7.50000/(1.00000+EXP( - 0.100000*((STATES(1)+11.0000) - CONSTS(12))))
      RATES(3) =  ALGBRC(2)*(1.00000 - STATES(3)) -  ALGBRC(13)*STATES(3)
      ALGBRC(3) = ( 0.175000*EXP(((STATES(1)+100.000) - CONSTS(13))/- 23.0000))/(1.00000+EXP( 0.150000*((STATES(1)+79.0000) - CONSTS(13))))
      ALGBRC(14) = 0.300000/(1.00000+EXP( - 0.100000*((STATES(1)+32.0000) - CONSTS(13))))
      RATES(4) =  ALGBRC(3)*(1.00000 - STATES(4)) -  ALGBRC(14)*STATES(4)
      ALGBRC(4) = 1.00000/(1.00000+EXP(- 2.18200 -  0.181900*STATES(1)))
      ALGBRC(15) = 43.0000+1.00000/(EXP(- 5.49500+ 0.169100*STATES(1))+EXP(- 7.67700 -  0.0128000*STATES(1)))
      RATES(5) = (ALGBRC(4) - STATES(5))/ALGBRC(15)
      ALGBRC(16) = 1.00000/(( 7.19000e-05*(STATES(1) - 10.0000))/(1.00000 - EXP( - 0.148000*(STATES(1) - 10.0000)))+( 0.000131000*(STATES(1) - 10.0000))/(EXP( 0.0687000*(STATES(1) - 10.0000)) - 1.00000))
      ALGBRC(5) = 1.00000/(1.00000+EXP((STATES(1) - 16.0000)/- 13.6000))
      RATES(6) = (ALGBRC(5) - STATES(6))/ALGBRC(16)
      ALGBRC(6) =  0.0451600*EXP( 0.0357700*STATES(1))
      ALGBRC(17) =  0.0989000*EXP( - 0.0623700*STATES(1))
      RATES(7) =  ALGBRC(6)*(1.00000 - STATES(7)) -  ALGBRC(17)*STATES(7)
      ALGBRC(7) = ( 0.00541500*EXP((STATES(1)+33.5000)/- 5.00000))/(1.00000+ 0.0513350*EXP((STATES(1)+33.5000)/- 5.00000))
      ALGBRC(18) = ( 0.00541500*EXP((STATES(1)+33.5000)/5.00000))/(1.00000+ 0.0513350*EXP((STATES(1)+33.5000)/5.00000))
      RATES(8) =  ALGBRC(7)*(1.00000 - STATES(8)) -  ALGBRC(18)*STATES(8)
      ALGBRC(8) = 1.00000/(1.00000+EXP((STATES(1)+12.5000)/5.00000))
      ALGBRC(19) = 30.0000+200.000/(1.00000+EXP((STATES(1)+20.0000)/9.50000))
      RATES(10) = (ALGBRC(8) - STATES(10))/ALGBRC(19)
      ALGBRC(1) = STATES(1)+47.1300
      ALGBRC(12) = ( 0.320000*ALGBRC(1))/(1.00000 - EXP( - 0.100000*ALGBRC(1)))
      ALGBRC(22) =  0.0800000*EXP(- STATES(1)/11.0000)
      RATES(2) =  ALGBRC(12)*(1.00000 - STATES(2)) -  ALGBRC(22)*STATES(2)
      ALGBRC(9) = 1.00000/(1.00000+EXP((STATES(1)+10.0000)/- 6.24000))
      ALGBRC(20) = STATES(1)+40.0000
      ALGBRC(23) = 1.00000/(( 0.250000*EXP( - 0.0100000*STATES(1)))/(1.00000+EXP( - 0.0700000*STATES(1)))+( 0.0700000*EXP( - 0.0500000*ALGBRC(20)))/(1.00000+EXP( 0.0500000*ALGBRC(20))))
      RATES(11) = (ALGBRC(9) - STATES(11))/ALGBRC(23)
      ALGBRC(21) =  CONSTS(9)*STATES(2) ** 3.00000*STATES(3)*STATES(4)*(STATES(1) - CONSTS(51))
      ALGBRC(39) = ( (( (CONSTS(35)/CONSTS(38))*4.00000*STATES(1)*CONSTS(3) ** 2.00000)/( CONSTS(1)*CONSTS(2)))*( STATES(9)*EXP(( 2.00000*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) -  0.341000*CONSTS(30)))/(EXP(( 2.00000*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) - 1.00000)
      ALGBRC(40) =  ALGBRC(39)*STATES(10)*STATES(11)*STATES(12)
      ALGBRC(41) = ( (( (( (CONSTS(36)/CONSTS(38))*STATES(10)*STATES(11)*STATES(12))/(1.00000+ALGBRC(39)/CONSTS(37)))*1000.00*STATES(1)*CONSTS(3) ** 2.00000)/( CONSTS(1)*CONSTS(2)))*( CONSTS(18)*EXP(( STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) - CONSTS(16)))/(EXP(( STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) - 1.00000)
      ALGBRC(26) = 1.00000/(1.00000+ 2.50000*EXP( 0.100000*(STATES(1)+28.0000)))
      ALGBRC(27) =  CONSTS(17)*ALGBRC(26)*STATES(5)* (CONSTS(16)/4.00000) ** (1.0 / 2)*(STATES(1) - CONSTS(52))
      ALGBRC(28) =  CONSTS(19)*STATES(6) ** 2.00000*(STATES(1) - CONSTS(53))
      ALGBRC(29) =  CONSTS(20)*STATES(7)*STATES(8)*(STATES(1) - CONSTS(52))
      ALGBRC(24) = 1.00000/(2.00000+EXP( (( 1.62000*CONSTS(3))/( CONSTS(1)*CONSTS(2)))*(STATES(1) - CONSTS(52))))
      ALGBRC(25) =  (( CONSTS(14)*ALGBRC(24)*CONSTS(16))/(CONSTS(16)+CONSTS(15)))*(STATES(1) - CONSTS(52))
      ALGBRC(30) = 1.00000/(1.00000+EXP((7.48800 - STATES(1))/5.98000))
      ALGBRC(31) =  CONSTS(21)*ALGBRC(30)*(STATES(1) - CONSTS(52))
      ALGBRC(34) =  (CONSTS(27)/( (CONSTS(26) ** 3.00000+CONSTS(10) ** 3.00000)*(CONSTS(25)+CONSTS(30))*(1.00000+ CONSTS(28)*EXP(( (CONSTS(29) - 1.00000)*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))))))*( EXP(( CONSTS(29)*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2)))*CONSTS(11) ** 3.00000*CONSTS(30) -  EXP(( (CONSTS(29) - 1.00000)*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2)))*CONSTS(10) ** 3.00000*STATES(9))
      ALGBRC(32) = 1.00000/(1.00000+ 0.124500*EXP(( - 0.100000*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2)))+ 0.0365000*CONSTS(54)*EXP(( - STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))))
      ALGBRC(33) = ( (( CONSTS(22)*ALGBRC(32))/(1.00000+CONSTS(23)/CONSTS(11) ** 1.50000))*CONSTS(16))/(CONSTS(16)+CONSTS(24))
      ALGBRC(35) = ( CONSTS(32)*STATES(9))/(CONSTS(31)+STATES(9))
      ALGBRC(36) =  (( CONSTS(1)*CONSTS(2))/( 2.00000*CONSTS(3)))*log(CONSTS(30)/STATES(9))
      ALGBRC(37) =  CONSTS(33)*(STATES(1) - ALGBRC(36))
      ALGBRC(38) =  CONSTS(34)*(STATES(1) - CONSTS(51))
      ALGBRC(11) = TERNRY(VOI.GE.CONSTS(4).AND.VOI.LE.CONSTS(5).AND.(VOI - CONSTS(4)) -  INT((VOI - CONSTS(4))/CONSTS(6))*CONSTS(6).LE.CONSTS(7), CONSTS(8), 0.00000)
      RATES(1) = - (ALGBRC(21)+ALGBRC(40)+ALGBRC(41)+ALGBRC(27)+ALGBRC(28)+ALGBRC(29)+ALGBRC(25)+ALGBRC(31)+ALGBRC(34)+ALGBRC(33)+ALGBRC(35)+ALGBRC(38)+ALGBRC(37)+ALGBRC(11))
      ALGBRC(43) = 1.00000/(1.00000+2000.00/STATES(13) ** 3.00000)
      ALGBRC(44) = ( CONSTS(44)*STATES(10)*STATES(11)*STATES(12)*( ALGBRC(43)*STATES(13) - STATES(9)))/(1.00000+ 1.65000*EXP(STATES(1)/20.0000))
      ALGBRC(45) =  CONSTS(45)*(STATES(13) - STATES(9))
      ALGBRC(42) = CONSTS(49)/(1.00000+CONSTS(50)/STATES(9) ** 2.00000)
      ALGBRC(47) = 1.00000/(1.00000+( CONSTS(47)*CONSTS(46))/CONSTS(46)+STATES(13) ** 2.00000)
      RATES(13) = ( ALGBRC(47)*((ALGBRC(42) - ALGBRC(45)) - ALGBRC(44))*CONSTS(42))/CONSTS(48)
      ALGBRC(46) = 1.00000/(1.00000+( CONSTS(41)*CONSTS(40))/CONSTS(40)+STATES(9) ** 2.00000)
      RATES(9) =  ALGBRC(46)*(((ALGBRC(44)+ALGBRC(45)) - ALGBRC(42)) -  (( CONSTS(43)*CONSTS(38))/( 2.00000*CONSTS(3)*CONSTS(42)))*((ALGBRC(40)+ALGBRC(37)+ALGBRC(35)) -  2.00000*ALGBRC(34)))
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(10) = 1.00000/(1.00000+STATES(9)/CONSTS(39) ** 3.00000)
      ALGBRC(2) =  0.135000*EXP(((STATES(1)+80.0000) - CONSTS(12))/- 6.80000)
      ALGBRC(13) = 7.50000/(1.00000+EXP( - 0.100000*((STATES(1)+11.0000) - CONSTS(12))))
      ALGBRC(3) = ( 0.175000*EXP(((STATES(1)+100.000) - CONSTS(13))/- 23.0000))/(1.00000+EXP( 0.150000*((STATES(1)+79.0000) - CONSTS(13))))
      ALGBRC(14) = 0.300000/(1.00000+EXP( - 0.100000*((STATES(1)+32.0000) - CONSTS(13))))
      ALGBRC(4) = 1.00000/(1.00000+EXP(- 2.18200 -  0.181900*STATES(1)))
      ALGBRC(15) = 43.0000+1.00000/(EXP(- 5.49500+ 0.169100*STATES(1))+EXP(- 7.67700 -  0.0128000*STATES(1)))
      ALGBRC(16) = 1.00000/(( 7.19000e-05*(STATES(1) - 10.0000))/(1.00000 - EXP( - 0.148000*(STATES(1) - 10.0000)))+( 0.000131000*(STATES(1) - 10.0000))/(EXP( 0.0687000*(STATES(1) - 10.0000)) - 1.00000))
      ALGBRC(5) = 1.00000/(1.00000+EXP((STATES(1) - 16.0000)/- 13.6000))
      ALGBRC(6) =  0.0451600*EXP( 0.0357700*STATES(1))
      ALGBRC(17) =  0.0989000*EXP( - 0.0623700*STATES(1))
      ALGBRC(7) = ( 0.00541500*EXP((STATES(1)+33.5000)/- 5.00000))/(1.00000+ 0.0513350*EXP((STATES(1)+33.5000)/- 5.00000))
      ALGBRC(18) = ( 0.00541500*EXP((STATES(1)+33.5000)/5.00000))/(1.00000+ 0.0513350*EXP((STATES(1)+33.5000)/5.00000))
      ALGBRC(8) = 1.00000/(1.00000+EXP((STATES(1)+12.5000)/5.00000))
      ALGBRC(19) = 30.0000+200.000/(1.00000+EXP((STATES(1)+20.0000)/9.50000))
      ALGBRC(1) = STATES(1)+47.1300
      ALGBRC(12) = ( 0.320000*ALGBRC(1))/(1.00000 - EXP( - 0.100000*ALGBRC(1)))
      ALGBRC(22) =  0.0800000*EXP(- STATES(1)/11.0000)
      ALGBRC(9) = 1.00000/(1.00000+EXP((STATES(1)+10.0000)/- 6.24000))
      ALGBRC(20) = STATES(1)+40.0000
      ALGBRC(23) = 1.00000/(( 0.250000*EXP( - 0.0100000*STATES(1)))/(1.00000+EXP( - 0.0700000*STATES(1)))+( 0.0700000*EXP( - 0.0500000*ALGBRC(20)))/(1.00000+EXP( 0.0500000*ALGBRC(20))))
      ALGBRC(21) =  CONSTS(9)*STATES(2) ** 3.00000*STATES(3)*STATES(4)*(STATES(1) - CONSTS(51))
      ALGBRC(39) = ( (( (CONSTS(35)/CONSTS(38))*4.00000*STATES(1)*CONSTS(3) ** 2.00000)/( CONSTS(1)*CONSTS(2)))*( STATES(9)*EXP(( 2.00000*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) -  0.341000*CONSTS(30)))/(EXP(( 2.00000*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) - 1.00000)
      ALGBRC(40) =  ALGBRC(39)*STATES(10)*STATES(11)*STATES(12)
      ALGBRC(41) = ( (( (( (CONSTS(36)/CONSTS(38))*STATES(10)*STATES(11)*STATES(12))/(1.00000+ALGBRC(39)/CONSTS(37)))*1000.00*STATES(1)*CONSTS(3) ** 2.00000)/( CONSTS(1)*CONSTS(2)))*( CONSTS(18)*EXP(( STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) - CONSTS(16)))/(EXP(( STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))) - 1.00000)
      ALGBRC(26) = 1.00000/(1.00000+ 2.50000*EXP( 0.100000*(STATES(1)+28.0000)))
      ALGBRC(27) =  CONSTS(17)*ALGBRC(26)*STATES(5)* (CONSTS(16)/4.00000) ** (1.0 / 2)*(STATES(1) - CONSTS(52))
      ALGBRC(28) =  CONSTS(19)*STATES(6) ** 2.00000*(STATES(1) - CONSTS(53))
      ALGBRC(29) =  CONSTS(20)*STATES(7)*STATES(8)*(STATES(1) - CONSTS(52))
      ALGBRC(24) = 1.00000/(2.00000+EXP( (( 1.62000*CONSTS(3))/( CONSTS(1)*CONSTS(2)))*(STATES(1) - CONSTS(52))))
      ALGBRC(25) =  (( CONSTS(14)*ALGBRC(24)*CONSTS(16))/(CONSTS(16)+CONSTS(15)))*(STATES(1) - CONSTS(52))
      ALGBRC(30) = 1.00000/(1.00000+EXP((7.48800 - STATES(1))/5.98000))
      ALGBRC(31) =  CONSTS(21)*ALGBRC(30)*(STATES(1) - CONSTS(52))
      ALGBRC(34) =  (CONSTS(27)/( (CONSTS(26) ** 3.00000+CONSTS(10) ** 3.00000)*(CONSTS(25)+CONSTS(30))*(1.00000+ CONSTS(28)*EXP(( (CONSTS(29) - 1.00000)*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))))))*( EXP(( CONSTS(29)*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2)))*CONSTS(11) ** 3.00000*CONSTS(30) -  EXP(( (CONSTS(29) - 1.00000)*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2)))*CONSTS(10) ** 3.00000*STATES(9))
      ALGBRC(32) = 1.00000/(1.00000+ 0.124500*EXP(( - 0.100000*STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2)))+ 0.0365000*CONSTS(54)*EXP(( - STATES(1)*CONSTS(3))/( CONSTS(1)*CONSTS(2))))
      ALGBRC(33) = ( (( CONSTS(22)*ALGBRC(32))/(1.00000+CONSTS(23)/CONSTS(11) ** 1.50000))*CONSTS(16))/(CONSTS(16)+CONSTS(24))
      ALGBRC(35) = ( CONSTS(32)*STATES(9))/(CONSTS(31)+STATES(9))
      ALGBRC(36) =  (( CONSTS(1)*CONSTS(2))/( 2.00000*CONSTS(3)))*log(CONSTS(30)/STATES(9))
      ALGBRC(37) =  CONSTS(33)*(STATES(1) - ALGBRC(36))
      ALGBRC(38) =  CONSTS(34)*(STATES(1) - CONSTS(51))
      ALGBRC(11) = TERNRY(VOI.GE.CONSTS(4).AND.VOI.LE.CONSTS(5).AND.(VOI - CONSTS(4)) -  INT((VOI - CONSTS(4))/CONSTS(6))*CONSTS(6).LE.CONSTS(7), CONSTS(8), 0.00000)
      ALGBRC(43) = 1.00000/(1.00000+2000.00/STATES(13) ** 3.00000)
      ALGBRC(44) = ( CONSTS(44)*STATES(10)*STATES(11)*STATES(12)*( ALGBRC(43)*STATES(13) - STATES(9)))/(1.00000+ 1.65000*EXP(STATES(1)/20.0000))
      ALGBRC(45) =  CONSTS(45)*(STATES(13) - STATES(9))
      ALGBRC(42) = CONSTS(49)/(1.00000+CONSTS(50)/STATES(9) ** 2.00000)
      ALGBRC(47) = 1.00000/(1.00000+( CONSTS(47)*CONSTS(46))/CONSTS(46)+STATES(13) ** 2.00000)
      ALGBRC(46) = 1.00000/(1.00000+( CONSTS(41)*CONSTS(40))/CONSTS(40)+STATES(9) ** 2.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