C
C There are a total of 9 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 14 entries in the constant variable array.
C
C
C VOI is time in component environment (day).
C STATES(1) is B0 in component B0 (cells_per_GC).
C CONSTS(1) is pr in component kinetic_parameters (dimensionless).
C CONSTS(2) is mu in component kinetic_parameters (first_order_rate_constant).
C CONSTS(3) is rho in component kinetic_parameters (first_order_rate_constant).
C CONSTS(4) is delta_B in component kinetic_parameters (first_order_rate_constant).
C ALGBRC(1) is CT_star in component CT_star (cells_per_GC).
C STATES(2) is B1 in component B1 (cells_per_GC).
C ALGBRC(5) is alpha_B in component alpha_B (dimensionless).
C STATES(3) is B2 in component B2 (cells_per_GC).
C STATES(4) is B3 in component B3 (cells_per_GC).
C STATES(5) is B4 in component B4 (cells_per_GC).
C STATES(6) is B5 in component B5 (cells_per_GC).
C STATES(7) is B6 in component B6 (cells_per_GC).
C STATES(8) is B7 in component B7 (cells_per_GC).
C STATES(9) is B8 in component B8 (cells_per_GC).
C STATES(10) is B9 in component B9 (cells_per_GC).
C STATES(11) is B10 in component B10 (cells_per_GC).
C ALGBRC(2) is B_sum in component centroblasts_sum (cells_per_GC).
C STATES(12) is C in component C (cells_per_GC).
C CONSTS(5) is d in component C (first_order_rate_constant).
C STATES(13) is C_star in component C_star (cells_per_GC).
C ALGBRC(3) is CA in component CA (cells_per_GC).
C ALGBRC(6) is C_starsum in component centrocytes_sum (cells_per_GC).
C STATES(14) is M in component M (cells_per_GC).
C STATES(15) is A in component A (cells_per_GC).
C CONSTS(6) is z in component A (first_order_rate_constant).
C CONSTS(7) is u in component A (dimensionless).
C ALGBRC(4) is log_A in component A (dimensionless).
C STATES(16) is T in component T (cells_per_GC).
C CONSTS(8) is p in component T (first_order_rate_constant).
C CONSTS(9) is sigma in component T (first_order_rate_constant).
C CONSTS(10) is delta_T in component T (first_order_rate_constant).
C ALGBRC(7) is alpha_T in component alpha_T (dimensionless).
C CONSTS(11) is SA in component CA (dimensionless).
C CONSTS(12) is ST in component CT_star (dimensionless).
C CONSTS(13) is KB in component alpha_B (dimensionless).
C CONSTS(14) is KT in component alpha_T (dimensionless).
C ALGBRC(8) is total in component total (cells_per_GC).
C ALGBRC(9) is log_total in component total (dimensionless).
C RATES(1) is d/dt B0 in component B0 (cells_per_GC).
C RATES(2) is d/dt B1 in component B1 (cells_per_GC).
C RATES(3) is d/dt B2 in component B2 (cells_per_GC).
C RATES(4) is d/dt B3 in component B3 (cells_per_GC).
C RATES(5) is d/dt B4 in component B4 (cells_per_GC).
C RATES(6) is d/dt B5 in component B5 (cells_per_GC).
C RATES(7) is d/dt B6 in component B6 (cells_per_GC).
C RATES(8) is d/dt B7 in component B7 (cells_per_GC).
C RATES(9) is d/dt B8 in component B8 (cells_per_GC).
C RATES(10) is d/dt B9 in component B9 (cells_per_GC).
C RATES(11) is d/dt B10 in component B10 (cells_per_GC).
C RATES(12) is d/dt C in component C (cells_per_GC).
C RATES(13) is d/dt C_star in component C_star (cells_per_GC).
C RATES(14) is d/dt M in component M (cells_per_GC).
C RATES(15) is d/dt A in component A (cells_per_GC).
C RATES(16) is d/dt T in component T (cells_per_GC).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = 3
      CONSTS(1) = 0.15
      CONSTS(2) = 3
      CONSTS(3) = 4
      CONSTS(4) = 0.8
      STATES(2) = 0
      STATES(3) = 0
      STATES(4) = 0
      STATES(5) = 0
      STATES(6) = 0
      STATES(7) = 0
      STATES(8) = 0
      STATES(9) = 0
      STATES(10) = 0
      STATES(11) = 0
      STATES(12) = 0
      CONSTS(5) = 2
      STATES(13) = 0
      STATES(14) = 0
      STATES(15) = 500
      CONSTS(6) = 0.02
      CONSTS(7) = 0.15
      STATES(16) = 0
      CONSTS(8) = 2
      CONSTS(9) = 5
      CONSTS(10) = 0.8
      CONSTS(11) = 500
      CONSTS(12) = 50
      CONSTS(13) = 1e4
      CONSTS(14) = 100
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      RATES(12) =  CONSTS(5)*CONSTS(3)*STATES(11)*1.00000 -  CONSTS(2)*STATES(12)
      ALGBRC(1) = ( STATES(13)*STATES(16))/( CONSTS(12)*1.00000+STATES(13))
      RATES(1) =  CONSTS(1)*CONSTS(2)*ALGBRC(1) - ( CONSTS(3)*STATES(1)+ CONSTS(4)*STATES(1))
      ALGBRC(3) = ( STATES(12)*STATES(15))/( CONSTS(11)*1.00000+STATES(15))
      RATES(13) =  CONSTS(2)*ALGBRC(3) -  CONSTS(2)*STATES(13)
      RATES(14) =  (1.00000 - CONSTS(1))*CONSTS(2)*ALGBRC(1)
      RATES(15) =  - CONSTS(6)*STATES(15) -  CONSTS(7)*ALGBRC(3)*1.00000
      ALGBRC(2) = STATES(2)+STATES(2)+STATES(4)+STATES(5)+STATES(6)+STATES(7)+STATES(7)+STATES(8)+STATES(9)+STATES(10)+STATES(11)
      ALGBRC(5) = CONSTS(13)/(CONSTS(13)+ALGBRC(2)/1.00000)
      RATES(2) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(1) - ( CONSTS(3)*STATES(2)+ CONSTS(4)*STATES(2))
      RATES(3) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(2) - ( CONSTS(3)*STATES(3)+ CONSTS(4)*STATES(3))
      RATES(4) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(3) - ( CONSTS(3)*STATES(4)+ CONSTS(4)*STATES(4))
      RATES(5) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(4) - ( CONSTS(3)*STATES(5)+ CONSTS(4)*STATES(5))
      RATES(6) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(5) - ( CONSTS(3)*STATES(6)+ CONSTS(4)*STATES(6))
      RATES(7) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(6) - ( CONSTS(3)*STATES(7)+ CONSTS(4)*STATES(7))
      RATES(8) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(7) - ( CONSTS(3)*STATES(8)+ CONSTS(4)*STATES(8))
      RATES(9) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(8) - ( CONSTS(3)*STATES(9)+ CONSTS(4)*STATES(9))
      RATES(10) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(9) - ( CONSTS(3)*STATES(10)+ CONSTS(4)*STATES(10))
      RATES(11) =  CONSTS(3)*(1.00000+ALGBRC(5))*STATES(10) - ( CONSTS(3)*STATES(11)+ CONSTS(4)*STATES(11))
      ALGBRC(7) = CONSTS(14)/(CONSTS(14)+STATES(16)/1.00000)
      RATES(16) = ( CONSTS(9)*1.00000+ CONSTS(8)*ALGBRC(7)*ALGBRC(1)) -  CONSTS(10)*STATES(16)
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(1) = ( STATES(13)*STATES(16))/( CONSTS(12)*1.00000+STATES(13))
      ALGBRC(3) = ( STATES(12)*STATES(15))/( CONSTS(11)*1.00000+STATES(15))
      ALGBRC(2) = STATES(2)+STATES(2)+STATES(4)+STATES(5)+STATES(6)+STATES(7)+STATES(7)+STATES(8)+STATES(9)+STATES(10)+STATES(11)
      ALGBRC(5) = CONSTS(13)/(CONSTS(13)+ALGBRC(2)/1.00000)
      ALGBRC(7) = CONSTS(14)/(CONSTS(14)+STATES(16)/1.00000)
      ALGBRC(4) = arbitrary_log(STATES(15)/1.00000, 10)
      ALGBRC(6) = STATES(12)+STATES(13)
      ALGBRC(8) = ALGBRC(2)+ALGBRC(6)
      ALGBRC(9) = arbitrary_log(ALGBRC(8)/1.00000+1.00000e-12, 10)
      RETURN
      END