C
C There are a total of 6 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 51 entries in the constant variable array.
C
C
C VOI is time in component environment (day).
C ALGBRC(1) is A_B in component A_B (mm2).
C STATES(1) is r_B in component r_B (mm).
C CONSTS(1) is X_c in component r_B (pM).
C CONSTS(2) is X_b in component r_B (pM).
C CONSTS(3) is X_y in component r_B (pM).
C CONSTS(4) is k_rB in component r_B (first_order_rate_constant).
C CONSTS(5) is k_for in component r_B (mm_per_day).
C CONSTS(6) is k_res in component r_B (mm_per_day).
C STATES(2) is xc in component xc (pM).
C STATES(3) is xb in component xb (pM).
C STATES(4) is xy in component xy (pM).
C CONSTS(7) is X_bss in component xy (pM).
C CONSTS(8) is X_yss in component xy (pM).
C CONSTS(9) is k_byp in component xy (first_order_rate_constant).
C CONSTS(10) is k_yd in component xy (first_order_rate_constant).
C ALGBRC(6) is F_sti in component F_sti (pM_N_per_mm2).
C CONSTS(11) is k_Fs in component F_sti (mm2_per_N).
C CONSTS(12) is k_y in component F_sti (per_pM).
C ALGBRC(5) is F_s in component F_s (N_per_mm2).
C ALGBRC(3) is F_a in component F_s (newton).
C STATES(5) is x_no in component x_no (pM).
C CONSTS(13) is X_noe in component x_no (flux).
C CONSTS(14) is k_yno in component x_no (mm2_day_per_dyn).
C CONSTS(15) is k_nod in component x_no (first_order_rate_constant).
C STATES(6) is x_pge in component x_pge (pM).
C CONSTS(16) is X_pgex in component x_pge (flux).
C CONSTS(17) is k_pged in component x_pge (first_order_rate_constant).
C CONSTS(18) is k_nopge in component x_pge (first_order_rate_constant).
C CONSTS(19) is k_ypge in component x_pge (mm2_day_per_dyn).
C STATES(7) is x_opg in component x_opg (pM).
C CONSTS(20) is k_opgd in component x_opg (first_order_rate_constant).
C CONSTS(21) is k_nopg in component x_opg (first_order_rate_constant).
C CONSTS(22) is Io in component model_parameters (flux).
C CONSTS(23) is K_o_p in component model_parameters (picomole_per_day_picomole_cells).
C STATES(8) is xr in component xr (pM).
C ALGBRC(4) is pi_c in component model_parameters (dimensionless).
C STATES(9) is x_kl in component x_kl (pM).
C ALGBRC(2) is pi_L in component x_kl (dimensionless).
C CONSTS(24) is K_l_p in component x_kl (dimensionless).
C CONSTS(25) is K in component x_kl (pM).
C CONSTS(26) is k_nokl in component x_kl (first_order_rate_constant).
C CONSTS(27) is Il in component x_kl (flux).
C CONSTS(28) is rl in component x_kl (flux).
C CONSTS(29) is k1 in component x_kl (second_order_rate_constant).
C CONSTS(30) is k2 in component x_kl (first_order_rate_constant).
C CONSTS(31) is k3 in component x_kl (second_order_rate_constant).
C CONSTS(32) is k4 in component x_kl (first_order_rate_constant).
C CONSTS(33) is ko in component x_kl (first_order_rate_constant).
C CONSTS(51) is pi_p in component model_parameters (dimensionless).
C CONSTS(34) is D_R in component xr (flux).
C CONSTS(35) is k_pger in component xr (first_order_rate_constant).
C CONSTS(47) is D_B in component model_parameters (first_order_rate_constant).
C CONSTS(36) is k_B in component xb (first_order_rate_constant).
C CONSTS(37) is D_C in component xc (flux).
C CONSTS(38) is D_A in component xc (first_order_rate_constant).
C CONSTS(39) is f0 in component model_parameters (dimensionless).
C CONSTS(40) is dB in component model_parameters (first_order_rate_constant).
C CONSTS(41) is C_s in component model_parameters (pM).
C CONSTS(48) is P in component model_parameters (pM).
C CONSTS(49) is P_0 in component model_parameters (pM).
C CONSTS(50) is P_s in component model_parameters (pM).
C CONSTS(42) is SP in component model_parameters (flux).
C CONSTS(43) is k5 in component model_parameters (second_order_rate_constant).
C CONSTS(44) is k6 in component model_parameters (first_order_rate_constant).
C CONSTS(45) is IP in component model_parameters (flux).
C CONSTS(46) is kP in component model_parameters (first_order_rate_constant).
C RATES(1) is d/dt r_B in component r_B (mm).
C RATES(4) is d/dt xy in component xy (pM).
C RATES(5) is d/dt x_no in component x_no (pM).
C RATES(6) is d/dt x_pge in component x_pge (pM).
C RATES(7) is d/dt x_opg in component x_opg (pM).
C RATES(9) is d/dt x_kl in component x_kl (pM).
C RATES(8) is d/dt xr in component xr (pM).
C RATES(3) is d/dt xb in component xb (pM).
C RATES(2) is d/dt xc in component xc (pM).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = 0.9912
      CONSTS(1) = 9.127e-4
      CONSTS(2) = 7.282e-4
      CONSTS(3) = 7.300e-3
      CONSTS(4) = 1.0e0
      CONSTS(5) = 1.0e-3
      CONSTS(6) = 10.0e-3
      STATES(2) = 9.127e-4
      STATES(3) = 7.282e-4
      STATES(4) = 7.300e-3
      CONSTS(7) = 7.282e-4
      CONSTS(8) = 7.300e-3
      CONSTS(9) = 1.00e-1
      CONSTS(10) = 1.00e0
      CONSTS(11) = 1.00e0
      CONSTS(12) = 1.00e0
      STATES(5) = 0.02
      CONSTS(13) = 0e0
      CONSTS(14) = 2e4
      CONSTS(15) = 1e3
      STATES(6) = 0.01
      CONSTS(16) = 0e0
      CONSTS(17) = 1e2
      CONSTS(18) = 1e0
      CONSTS(19) = 1e2
      STATES(7) = 0.01
      CONSTS(20) = 3.5e-1
      CONSTS(21) = 1e1
      CONSTS(22) = 0.0
      CONSTS(23) = 2e5
      STATES(8) = 7.734e-4
      STATES(9) = 0.01
      CONSTS(24) = 3e6
      CONSTS(25) = 10
      CONSTS(26) = 1e2
      CONSTS(27) = 0
      CONSTS(28) = 1e3
      CONSTS(29) = 1e-2
      CONSTS(30) = 10
      CONSTS(31) = 5.8e-4
      CONSTS(32) = 1.7e-2
      CONSTS(33) = 0.35
      CONSTS(34) = 7e-4
      CONSTS(35) = 1e-4
      CONSTS(36) = 0.189
      CONSTS(37) = 0.189
      CONSTS(38) = 0.7
      CONSTS(39) = 0.05
      CONSTS(40) = 0.7
      CONSTS(41) = 5e-3
      CONSTS(42) = 250
      CONSTS(43) = 0.02
      CONSTS(44) = 3
      CONSTS(45) = 0
      CONSTS(46) = 86
      CONSTS(47) =  CONSTS(39)*CONSTS(40)
      CONSTS(48) = CONSTS(45)/CONSTS(46)
      CONSTS(49) = CONSTS(42)/CONSTS(46)
      CONSTS(50) = CONSTS(44)/CONSTS(43)
      CONSTS(51) = (CONSTS(48)+CONSTS(49))/(CONSTS(48)+CONSTS(50))
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      RATES(1) = ( (CONSTS(5)/CONSTS(2))*STATES(3)+ (1.00000/CONSTS(3))*STATES(4)) - ( (CONSTS(6)/CONSTS(1))*STATES(2)+ CONSTS(4)*STATES(1))
      RATES(4) =  CONSTS(9)*(STATES(3) - CONSTS(7)) -  CONSTS(10)*(STATES(4) - CONSTS(8))
      RATES(9) = (CONSTS(28)+CONSTS(27)) - ( CONSTS(26)*STATES(5)+ CONSTS(28)*((1.00000+ (CONSTS(29)/CONSTS(30))*STATES(7)+ (CONSTS(31)/CONSTS(32))*CONSTS(25))/( CONSTS(24)*CONSTS(51)*STATES(3)))*STATES(9))
      ALGBRC(4) = (STATES(2)+ CONSTS(39)*CONSTS(41))/(STATES(2)+CONSTS(41))
      RATES(7) = ( (CONSTS(23)/1.00000)*ALGBRC(4)*STATES(8)+CONSTS(22)+ CONSTS(21)*STATES(5)) -  CONSTS(20)*STATES(7)
      RATES(8) = ( CONSTS(34)*ALGBRC(4)+ CONSTS(35)*STATES(6)) -  (CONSTS(47)/ALGBRC(4))*STATES(8)
      RATES(3) =  (CONSTS(47)/ALGBRC(4))*STATES(8) -  CONSTS(36)*STATES(3)
      ALGBRC(2) =  (( ( (CONSTS(31)/CONSTS(32))*CONSTS(24))*CONSTS(51)*STATES(3))/(1.00000+( CONSTS(31)*CONSTS(25))/CONSTS(32)+ (CONSTS(29)/( CONSTS(30)*CONSTS(33)))*( ((CONSTS(23)/1.00000)/CONSTS(51))*STATES(8)+CONSTS(22))))*(1.00000+CONSTS(27)/CONSTS(28))
      RATES(2) =  CONSTS(37)*ALGBRC(2) -  CONSTS(38)*ALGBRC(4)*STATES(2)
      ALGBRC(1) =  STATES(1) ** 2.00000* 3.14159265358979
      ALGBRC(3) = TERNRY(VOI.GE.100.000.AND.VOI.LT.130.000, 10000.0, 100.000)
      ALGBRC(5) = ALGBRC(3)/ALGBRC(1)
      ALGBRC(6) = ( ALGBRC(5)*STATES(4))/(1.00000+EXP(- ( CONSTS(11)*ALGBRC(5)+ CONSTS(12)*STATES(4))))
      RATES(5) = ( CONSTS(14)*ALGBRC(6)+CONSTS(13)) -  CONSTS(15)*STATES(5)
      RATES(6) = ( CONSTS(19)*ALGBRC(6)+ CONSTS(18)*STATES(5)+CONSTS(16)) -  CONSTS(17)*STATES(6)
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(4) = (STATES(2)+ CONSTS(39)*CONSTS(41))/(STATES(2)+CONSTS(41))
      ALGBRC(2) =  (( ( (CONSTS(31)/CONSTS(32))*CONSTS(24))*CONSTS(51)*STATES(3))/(1.00000+( CONSTS(31)*CONSTS(25))/CONSTS(32)+ (CONSTS(29)/( CONSTS(30)*CONSTS(33)))*( ((CONSTS(23)/1.00000)/CONSTS(51))*STATES(8)+CONSTS(22))))*(1.00000+CONSTS(27)/CONSTS(28))
      ALGBRC(1) =  STATES(1) ** 2.00000* 3.14159265358979
      ALGBRC(3) = TERNRY(VOI.GE.100.000.AND.VOI.LT.130.000, 10000.0, 100.000)
      ALGBRC(5) = ALGBRC(3)/ALGBRC(1)
      ALGBRC(6) = ( ALGBRC(5)*STATES(4))/(1.00000+EXP(- ( CONSTS(11)*ALGBRC(5)+ CONSTS(12)*STATES(4))))
      RETURN
      END
      REAL FUNCTION TERNRY(TEST, VALA, VALB)
      LOGICAL TEST
      REAL VALA, VALB
      IF (TEST) THEN
        TERNRY = VALA
      ELSE
        TERNRY = VALB
      ENDIF
      RETURN
      END