C
C There are a total of 18 entries in the algebraic variable array.
C There are a total of 9 entries in each of the rate and state variable arrays.
C There are a total of 49 entries in the constant variable array.
C
C
C VOI is time in component Environment (second).
C CONSTS(1) is Eta_b in component Parameters (mmHg_sec).
C CONSTS(2) is Rho_b in component Parameters (gram_per_cubic_meter).
C CONSTS(3) is E_p in component Parameters (mmHg).
C CONSTS(4) is r_p in component Parameters (meter).
C CONSTS(5) is l_p in component Parameters (meter).
C CONSTS(6) is h_p in component Parameters (meter).
C CONSTS(7) is r_a0 in component Parameters (meter).
C CONSTS(8) is h_a0 in component Parameters (meter).
C CONSTS(9) is sigma_ae0 in component Parameters (mmHg).
C CONSTS(10) is K_a_sigma in component Parameters (dimensionless).
C CONSTS(11) is sigma_ac in component Parameters (mmHg).
C CONSTS(12) is T_a_max0 in component Parameters (mmHg_m).
C CONSTS(13) is r_am in component Parameters (meter).
C CONSTS(14) is r_at in component Parameters (meter).
C CONSTS(15) is n_am in component Parameters (dimensionless).
C CONSTS(16) is Eta_a in component Parameters (mmHg_sec).
C CONSTS(17) is Q_bl in component Parameters (cubic_m_per_sec).
C CONSTS(18) is P_ic in component Parameters (mmHg).
C CONSTS(19) is R_v in component Parameters (mmHg_sec_per_cubic_m).
C CONSTS(20) is C_v in component Parameters (cubic_m_per_mmHg).
C CONSTS(21) is P_ai in component Parameters (mmHg).
C CONSTS(22) is T_myo0 in component Parameters (mmHg_m).
C CONSTS(23) is T_myo_s in component Parameters (mmHg_m).
C CONSTS(24) is Tau_myo in component Parameters (second).
C CONSTS(25) is Tau_shear in component Parameters (second).
C CONSTS(26) is G_myo in component Parameters (dimensionless).
C CONSTS(27) is G_shear in component Parameters (dimensionless).
C CONSTS(28) is G_neuro in component Parameters (dimensionless).
C CONSTS(29) is G_meta in component Parameters (dimensionless).
C CONSTS(30) is x_ini in component Parameters (dimensionless).
C CONSTS(31) is Tau_neuro in component Parameters (second).
C CONSTS(32) is Tau_meta in component Parameters (second).
C CONSTS(33) is conc_CO2_a in component Parameters (mol_per_m3).
C CONSTS(34) is M_CO2_0 in component Parameters (mol_per_sec).
C CONSTS(35) is alpha_tv in component Parameters (dimensionless).
C CONSTS(36) is beta_tv in component Parameters (mol_per_m3).
C CONSTS(37) is vol_PVC in component Parameters (cubic_m).
C CONSTS(38) is K_shear in component Parameters (second).
C CONSTS(39) is K_Ra in component Parameters (mmHg_sec_m).
C CONSTS(40) is conc_CO2_t_sp in component Parameters (mol_per_m3).
C CONSTS(41) is f_Q in component Parameters (dimensionless).
C CONSTS(42) is n_QM in component Parameters (dimensionless).
C CONSTS(43) is K_Va in component Parameters (meter).
C CONSTS(44) is P_vo in component Parameters (mmHg).
C STATES(1) is Q_Lp in component PCA (cubic_m_per_sec).
C STATES(2) is P_p in component PCA (mmHg).
C CONSTS(47) is C_p in component PCA_Components (cubic_m_per_mmHg).
C CONSTS(48) is L_p in component PCA_Components (mmHg_sec2_per_m3).
C CONSTS(49) is R_p in component PCA_Components (mmHg_sec_per_cubic_m).
C ALGBRC(1) is R_a in component AC_Components (mmHg_sec_per_cubic_m).
C ALGBRC(12) is P_a in component Arteriolar_Circulation (mmHg).
C ALGBRC(3) is C_a in component AC_Components (cubic_m_per_mmHg).
C ALGBRC(14) is Q_a in component AC_Flow (cubic_m_per_sec).
C STATES(3) is P_v in component Venous_Circulation (mmHg).
C STATES(4) is r_a in component Arteriolar_Circulation (meter).
C ALGBRC(5) is T_ae in component AC_Tension (mmHg_m).
C ALGBRC(9) is T_am in component AC_Tension (mmHg_m).
C ALGBRC(10) is T_av in component AC_Tension (mmHg_m).
C ALGBRC(11) is T_a in component AC_Tension (mmHg_m).
C ALGBRC(4) is h_a in component AC_Tension_Variables (meter).
C ALGBRC(8) is T_a_max in component AC_Tension_Variables (mmHg_m).
C ALGBRC(7) is M_s in component CBF_Regulation (dimensionless).
C ALGBRC(6) is M_s1 in component CBF_Regulation (dimensionless).
C STATES(5) is x_myo in component Myo_Regulation (dimensionless).
C STATES(6) is x_shear in component Shear_Regulation (dimensionless).
C STATES(7) is x_neuro in component Neuro_Regulation (dimensionless).
C STATES(8) is x_meta in component Meta_Regulation (dimensionless).
C ALGBRC(13) is A_myo in component Myo_Regulation (dimensionless).
C ALGBRC(15) is A_shear in component Shear_Regulation (dimensionless).
C ALGBRC(16) is A_neuro in component Neuro_Regulation (dimensionless).
C ALGBRC(2) is A_Meta in component Meta_Regulation (dimensionless).
C CONSTS(45) is G_a_Meta in component Meta_Regulation (m3_per_mol).
C STATES(9) is conc_CO2_t in component Meta_CO2_Kinetics (mol_per_m3).
C ALGBRC(17) is conc_CO2_v in component Meta_CO2_Kinetics (mol_per_m3).
C ALGBRC(18) is M_CO2 in component VC_Meta_Stim (mol_per_sec).
C CONSTS(46) is unity in component VC_Meta_Stim (dimensionless).
C RATES(2) is d/dt P_p in component PCA (mmHg).
C RATES(1) is d/dt Q_Lp in component PCA (cubic_m_per_sec).
C RATES(4) is d/dt r_a in component Arteriolar_Circulation (meter).
C RATES(3) is d/dt P_v in component Venous_Circulation (mmHg).
C RATES(5) is d/dt x_myo in component Myo_Regulation (dimensionless).
C RATES(6) is d/dt x_shear in component Shear_Regulation (dimensionless).
C RATES(7) is d/dt x_neuro in component Neuro_Regulation (dimensionless).
C RATES(8) is d/dt x_meta in component Meta_Regulation (dimensionless).
C RATES(9) is d/dt conc_CO2_t in component Meta_CO2_Kinetics (mol_per_m3).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      CONSTS(1) = 3E-5
      CONSTS(2) = 1.05E6
      CONSTS(3) = 12E3
      CONSTS(4) = 1.05E-3
      CONSTS(5) = 8.6E-2
      CONSTS(6) = 2.6E-4
      CONSTS(7) = 7.5E-5
      CONSTS(8) = 2.5E-5
      CONSTS(9) = 11.19
      CONSTS(10) = 4.5
      CONSTS(11) = 41.32
      CONSTS(12) = 1.50E-2
      CONSTS(13) = 1.28E-4
      CONSTS(14) = 1.74E-4
      CONSTS(15) = 1.75
      CONSTS(16) = 47.8
      CONSTS(17) = 8.8E-7
      CONSTS(18) = 10
      CONSTS(19) = 1.6E7
      CONSTS(20) = 2.5E-8
      CONSTS(21) = 96.0
      CONSTS(22) = 3.6E-3
      CONSTS(23) = 3E-3
      CONSTS(24) = 7
      CONSTS(25) = 60
      CONSTS(26) = 4
      CONSTS(27) = 0
      CONSTS(28) = -0.7
      CONSTS(29) = 0
      CONSTS(30) = 0
      CONSTS(31) = 5
      CONSTS(32) = 15
      CONSTS(33) = 20.65
      CONSTS(34) = 1.35E-7
      CONSTS(35) = 0.96
      CONSTS(36) = 8.9
      CONSTS(37) = 5.99E-6
      CONSTS(38) = 5.2E-7
      CONSTS(39) = 2.5E-9
      CONSTS(40) = 12.41
      CONSTS(41) = 1.25
      CONSTS(42) = 2.2
      CONSTS(43) = 120
      CONSTS(44) = 14
      STATES(1) = 8.8E-7
      STATES(2) = 91
      STATES(3) = 21
      STATES(4) = 7.7E-5
      STATES(5) = 0
      STATES(6) = 0
      STATES(7) = 0
      STATES(8) = 0
      CONSTS(45) = 0.59
      STATES(9) = 12.41
      CONSTS(46) = 1
      CONSTS(47) =  (( 3.00000* 3.14159265358979*CONSTS(4) ** 2.00000*CONSTS(4)/CONSTS(6)+1.00000 ** 2.00000)/( CONSTS(3)*(( 2.00000*CONSTS(4))/CONSTS(6)+1.00000)))*CONSTS(5)
      CONSTS(48) =  (CONSTS(2)/(  3.14159265358979*CONSTS(4) ** 2.00000))*CONSTS(5)*7.50000e-06
      CONSTS(49) =  (( 8.00000*CONSTS(1))/(  3.14159265358979*CONSTS(4) ** 4.00000))*CONSTS(5)
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      RATES(1) = ((CONSTS(21) - STATES(2)) -  STATES(1)*CONSTS(49))/CONSTS(48)
      ALGBRC(2) =  CONSTS(45)*(STATES(9) - CONSTS(40))
      RATES(8) = (ALGBRC(2) - STATES(8))/CONSTS(32)
      ALGBRC(1) = CONSTS(39)/STATES(4) ** 4.00000
      ALGBRC(4) =  (STATES(4) ** 2.00000+ 2.00000*CONSTS(7)*CONSTS(8)+CONSTS(8) ** 2.00000) ** (1.0 / 2) - STATES(4)
      ALGBRC(5) =  ALGBRC(4)*( CONSTS(9)*(EXP(( CONSTS(10)*(STATES(4) - CONSTS(7)))/CONSTS(7)) - 1.00000) - CONSTS(11))
      ALGBRC(6) =  CONSTS(26)*STATES(5)+ CONSTS(27)*STATES(6)+ CONSTS(28)*STATES(7)+ CONSTS(29)*STATES(8)+CONSTS(30)
      ALGBRC(7) = (EXP( 2.00000*ALGBRC(6)) - 1.00000)/(EXP( 2.00000*ALGBRC(6))+1.00000)
      ALGBRC(8) =  CONSTS(12)*(1.00000+ALGBRC(7))
      ALGBRC(9) =  ALGBRC(8)*EXP(- ABS((STATES(4) - CONSTS(13))/(CONSTS(14) - CONSTS(13))) ** CONSTS(15))
      RATES(4) = ( CONSTS(7)*(( STATES(2)*(ALGBRC(1)+CONSTS(19))*STATES(4)+ STATES(3)*ALGBRC(1)*STATES(4)) -  (ALGBRC(5)+ALGBRC(9)+ CONSTS(18)*(STATES(4)+ALGBRC(4)))*( 2.00000*ALGBRC(1)+CONSTS(19))))/( CONSTS(7)*STATES(4) ** 2.00000*CONSTS(43)*ALGBRC(1)*(ALGBRC(1)+CONSTS(19))+ CONSTS(16)*ALGBRC(4)*( 2.00000*ALGBRC(1)+CONSTS(19)))
      ALGBRC(10) =  (CONSTS(16)/CONSTS(7))*RATES(4)*ALGBRC(4)
      ALGBRC(11) = ALGBRC(5)+ALGBRC(9)+ALGBRC(10)
      ALGBRC(12) = (ALGBRC(11)+ CONSTS(18)*(STATES(4)+ALGBRC(4)))/STATES(4)
      RATES(2) = ( 2.00000*(( STATES(1)*ALGBRC(1) -  2.00000*STATES(2))+ 2.00000*ALGBRC(12)))/( ALGBRC(1)*CONSTS(47))
      RATES(3) =  (2.00000/CONSTS(20))*((ALGBRC(12) - STATES(3))/(ALGBRC(1)+CONSTS(19)) - (STATES(3) - CONSTS(44))/CONSTS(19))
      ALGBRC(13) = (ALGBRC(11) - CONSTS(22))/CONSTS(23)
      RATES(5) = (ALGBRC(13) - STATES(5))/CONSTS(24)
      ALGBRC(14) = ( 2.00000*(ALGBRC(12) - STATES(3)))/(ALGBRC(1)+CONSTS(19))
      ALGBRC(15) = ( CONSTS(38)*ALGBRC(14))/STATES(4) ** 3.00000 - 1.00000
      RATES(6) = (ALGBRC(15) - STATES(6))/CONSTS(25)
      ALGBRC(16) = TERNRY(VOI.GE.1.00000, 1.00000, 0.00000)
      RATES(7) = (ALGBRC(16) - STATES(7))/CONSTS(31)
      ALGBRC(17) =  CONSTS(35)*STATES(9)+CONSTS(36)
      ALGBRC(18) =  CONSTS(34)*(CONSTS(46)+ ALGBRC(16)*((CONSTS(41) - CONSTS(46))/CONSTS(42)))
      RATES(9) =  (1.00000/CONSTS(37))*(ALGBRC(18) -  ALGBRC(14)*(ALGBRC(17) - CONSTS(33)))
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(2) =  CONSTS(45)*(STATES(9) - CONSTS(40))
      ALGBRC(1) = CONSTS(39)/STATES(4) ** 4.00000
      ALGBRC(4) =  (STATES(4) ** 2.00000+ 2.00000*CONSTS(7)*CONSTS(8)+CONSTS(8) ** 2.00000) ** (1.0 / 2) - STATES(4)
      ALGBRC(5) =  ALGBRC(4)*( CONSTS(9)*(EXP(( CONSTS(10)*(STATES(4) - CONSTS(7)))/CONSTS(7)) - 1.00000) - CONSTS(11))
      ALGBRC(6) =  CONSTS(26)*STATES(5)+ CONSTS(27)*STATES(6)+ CONSTS(28)*STATES(7)+ CONSTS(29)*STATES(8)+CONSTS(30)
      ALGBRC(7) = (EXP( 2.00000*ALGBRC(6)) - 1.00000)/(EXP( 2.00000*ALGBRC(6))+1.00000)
      ALGBRC(8) =  CONSTS(12)*(1.00000+ALGBRC(7))
      ALGBRC(9) =  ALGBRC(8)*EXP(- ABS((STATES(4) - CONSTS(13))/(CONSTS(14) - CONSTS(13))) ** CONSTS(15))
      ALGBRC(10) =  (CONSTS(16)/CONSTS(7))*RATES(4)*ALGBRC(4)
      ALGBRC(11) = ALGBRC(5)+ALGBRC(9)+ALGBRC(10)
      ALGBRC(12) = (ALGBRC(11)+ CONSTS(18)*(STATES(4)+ALGBRC(4)))/STATES(4)
      ALGBRC(13) = (ALGBRC(11) - CONSTS(22))/CONSTS(23)
      ALGBRC(14) = ( 2.00000*(ALGBRC(12) - STATES(3)))/(ALGBRC(1)+CONSTS(19))
      ALGBRC(15) = ( CONSTS(38)*ALGBRC(14))/STATES(4) ** 3.00000 - 1.00000
      ALGBRC(16) = TERNRY(VOI.GE.1.00000, 1.00000, 0.00000)
      ALGBRC(17) =  CONSTS(35)*STATES(9)+CONSTS(36)
      ALGBRC(18) =  CONSTS(34)*(CONSTS(46)+ ALGBRC(16)*((CONSTS(41) - CONSTS(46))/CONSTS(42)))
      ALGBRC(3) = 1.03120/ALGBRC(1)
      RETURN
      END
      REAL FUNCTION TERNRY(TEST, VALA, VALB)
      LOGICAL TEST
      REAL VALA, VALB
      IF (TEST) THEN
        TERNRY = VALA
      ELSE
        TERNRY = VALB
      ENDIF
      RETURN
      END