C
C There are a total of 9 entries in the algebraic variable array.
C There are a total of 7 entries in each of the rate and state variable arrays.
C There are a total of 35 entries in the constant variable array.
C
C
C VOI is time in component environment (second).
C STATES(1) is IP3 in component IP3_dynamics (micromolar).
C ALGBRC(1) is j_IP3 in component IP3_dynamics (micromolar_micrometre_per_second).
C CONSTS(1) is J_IP3 in component IP3_dynamics (micromolar_micrometre_per_second).
C CONSTS(2) is k_0 in component IP3_dynamics (first_order_rate_constant).
C CONSTS(3) is k_degr in component IP3_dynamics (first_order_rate_constant).
C CONSTS(4) is IP3_0 in component IP3_dynamics (micromolar).
C CONSTS(5) is Ca_ER in component ER (micromolar).
C STATES(2) is Ca in component Calcium_dynamics (micromolar).
C CONSTS(6) is alpha in component Calcium_dynamics (dimensionless).
C ALGBRC(2) is J_channel in component Channel_kinetics (flux).
C ALGBRC(8) is J_pump in component SERCA_pump_kinetics (flux).
C ALGBRC(9) is J_leak in component Leak (flux).
C CONSTS(28) is R_buffering in component Calcium_buffering (flux).
C CONSTS(7) is J_max in component Channel_kinetics (flux).
C STATES(3) is h in component Channel_kinetics (dimensionless).
C CONSTS(8) is K_act in component Channel_kinetics (micromolar).
C CONSTS(9) is K_IP3 in component Channel_kinetics (micromolar).
C CONSTS(10) is K_inh in component Channel_kinetics (micromolar).
C CONSTS(11) is k_on in component Channel_kinetics (second_order_rate_constant).
C CONSTS(12) is V_max in component SERCA_pump_kinetics (flux).
C CONSTS(13) is K_p in component SERCA_pump_kinetics (micromolar).
C CONSTS(14) is L in component Leak (flux).
C CONSTS(15) is R1 in component Calcium_buffering (flux).
C CONSTS(16) is R2 in component Calcium_buffering (flux).
C STATES(4) is B1 in component Calcium_buffering (micromolar).
C STATES(5) is B2 in component Calcium_buffering (micromolar).
C STATES(6) is CaB1 in component Calcium_buffering (micromolar).
C STATES(7) is CaB2 in component Calcium_buffering (micromolar).
C ALGBRC(3) is k1_on in component Calcium_buffering (second_order_rate_constant).
C ALGBRC(4) is k1_off in component Calcium_buffering (first_order_rate_constant).
C ALGBRC(5) is k2_on in component Calcium_buffering (second_order_rate_constant).
C ALGBRC(6) is k2_off in component Calcium_buffering (first_order_rate_constant).
C CONSTS(17) is K1 in component Calcium_buffering (micromolar).
C CONSTS(18) is K2 in component Calcium_buffering (micromolar).
C CONSTS(19) is soma_or_neurite in component Plasma_membrane_extrusion_mechanisms (dimensionless).
C ALGBRC(7) is j_Ca in component Plasma_membrane_extrusion_mechanisms (micromolar_micrometre_per_second).
C CONSTS(20) is gamma_0 in component Plasma_membrane_extrusion_mechanisms (micrometre_per_second).
C CONSTS(31) is gamma in component Plasma_membrane_extrusion_mechanisms (micrometre_per_second).
C CONSTS(29) is gamma_s in component Plasma_membrane_extrusion_mechanisms (micrometre_per_second).
C CONSTS(30) is gamma_n in component Plasma_membrane_extrusion_mechanisms (micrometre_per_second).
C CONSTS(21) is delta in component Plasma_membrane_extrusion_mechanisms (dimensionless).
C CONSTS(22) is sigma in component Plasma_membrane_extrusion_mechanisms (per_micrometre).
C CONSTS(23) is w_n in component Plasma_membrane_extrusion_mechanisms (dimensionless).
C CONSTS(24) is w_s in component Plasma_membrane_extrusion_mechanisms (dimensionless).
C CONSTS(25) is sigma_soma_2D in component Plasma_membrane_extrusion_mechanisms (per_micrometre).
C CONSTS(26) is sigma_neurite_2D in component Plasma_membrane_extrusion_mechanisms (per_micrometre).
C CONSTS(27) is Ca_c in component Plasma_membrane_extrusion_mechanisms (micromolar).
C RATES(1) is d/dt IP3 in component IP3_dynamics (micromolar).
C RATES(2) is d/dt Ca in component Calcium_dynamics (micromolar).
C RATES(3) is d/dt h in component Channel_kinetics (dimensionless).
C RATES(4) is d/dt B1 in component Calcium_buffering (micromolar).
C RATES(6) is d/dt CaB1 in component Calcium_buffering (micromolar).
C RATES(5) is d/dt B2 in component Calcium_buffering (micromolar).
C RATES(7) is d/dt CaB2 in component Calcium_buffering (micromolar).
C
      SUBROUTINE initConsts(CONSTS, RATES, STATES)
      REAL CONSTS(*), RATES(*), STATES(*)
      STATES(1) = 3.0
      CONSTS(1) = 20.86
      CONSTS(2) = 1.188
      CONSTS(3) = 0.14
      CONSTS(4) = 0.16
      CONSTS(5) = 400.0
      STATES(2) = 0.05
      CONSTS(6) = 0.0
      CONSTS(7) = 3500.0
      STATES(3) = 0.8
      CONSTS(8) = 0.3
      CONSTS(9) = 0.8
      CONSTS(10) = 0.2
      CONSTS(11) = 2.7
      CONSTS(12) = 3.75
      CONSTS(13) = 0.27
      CONSTS(14) = 0.1
      CONSTS(15) = 0.1
      CONSTS(16) = 0.1
      STATES(4) = 450.0
      STATES(5) = 75.0
      STATES(6) = 0
      STATES(7) = 0
      CONSTS(17) = 10.0
      CONSTS(18) = 0.24
      CONSTS(19) = -1
      CONSTS(20) = 8.0
      CONSTS(21) = 1.45
      CONSTS(22) = 0.263
      CONSTS(23) = 0.377
      CONSTS(24) = 0.623
      CONSTS(25) = 0.132
      CONSTS(26) = 0.479
      CONSTS(27) = 0.2
      CONSTS(28) = CONSTS(15)+CONSTS(16)
      CONSTS(29) = ( CONSTS(20)*CONSTS(22))/( CONSTS(21)*CONSTS(26)*CONSTS(23)+ CONSTS(25)*CONSTS(24))
      CONSTS(31) = CONSTS(15)
      CONSTS(32) = - CONSTS(15)
      CONSTS(33) = CONSTS(16)
      CONSTS(34) = - CONSTS(16)
      CONSTS(30) = ( CONSTS(20)*CONSTS(22)*CONSTS(21))/( CONSTS(21)*CONSTS(26)*CONSTS(23)+ CONSTS(25)*CONSTS(24))
      CONSTS(31) = TERNRY(CONSTS(19).LE.0.00000, CONSTS(29), CONSTS(30))
      RETURN
      END
      SUBROUTINE computeRates(VOI, CONSTS,  RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      RATES(4) = CONSTS(31)
      RATES(6) = CONSTS(32)
      RATES(5) = CONSTS(33)
      RATES(7) = CONSTS(34)
      RATES(1) = - ( CONSTS(3)*(STATES(1) - CONSTS(4)))
      RATES(3) =  CONSTS(11)*(CONSTS(10) -  STATES(3)*(STATES(2)+CONSTS(10)))
      ALGBRC(2) =  CONSTS(7)* (STATES(1)/(STATES(1)+CONSTS(9)))*(STATES(2)/(STATES(2)+CONSTS(8)))*STATES(3) ** 3.00000*(1.00000 - STATES(2)/CONSTS(5))
      ALGBRC(8) =  CONSTS(12)*(STATES(2) ** 2.00000/(STATES(2) ** 2.00000+CONSTS(13) ** 2.00000))
      ALGBRC(9) =  CONSTS(14)*(1.00000 - STATES(2)/CONSTS(5))
      RATES(2) =  CONSTS(6)*(ALGBRC(2)+- ALGBRC(8)+ALGBRC(9))+CONSTS(28)
      RETURN
      END
      SUBROUTINE computeVariables(VOI, CONSTS, RATES, STATES, ALGBRC)
      REAL VOI, CONSTS(*), RATES(*), STATES(*), ALGBRC(*)
      ALGBRC(2) =  CONSTS(7)* (STATES(1)/(STATES(1)+CONSTS(9)))*(STATES(2)/(STATES(2)+CONSTS(8)))*STATES(3) ** 3.00000*(1.00000 - STATES(2)/CONSTS(5))
      ALGBRC(8) =  CONSTS(12)*(STATES(2) ** 2.00000/(STATES(2) ** 2.00000+CONSTS(13) ** 2.00000))
      ALGBRC(9) =  CONSTS(14)*(1.00000 - STATES(2)/CONSTS(5))
      ALGBRC(1) =  CONSTS(1)*EXP( - CONSTS(2)*VOI)
CALL minimize(minfunc_0, CONSTS, VARIABLES, <VAR>)
CALL minimize(minfunc_1, CONSTS, VARIABLES, <VAR>)
      ALGBRC(7) = TERNRY(STATES(2).GT.CONSTS(27),  CONSTS(31)*(STATES(2) - CONSTS(27)), 0.00000)
      RETURN
      END
REAL FUNCTION minfunc_0(CONSTS, VARIABLES)
            REAL CONSTS(*), VARIABLES(*) 
      minfunc_0[1] = abs(CONSTS(15) - (- ( ALGBRC(3)*STATES(2)*STATES(4))+ ALGBRC(4)*STATES(6)))
      minfunc_0[2] = abs(CONSTS(17) - ALGBRC(4)/ALGBRC(3))
      RETURN
      END
REAL FUNCTION minfunc_1(CONSTS, VARIABLES)
            REAL CONSTS(*), VARIABLES(*) 
      minfunc_1[1] = abs(CONSTS(16) - (- ( ALGBRC(5)*STATES(2)*STATES(5))+ ALGBRC(6)*STATES(7)))
      minfunc_1[2] = abs(CONSTS(18) - ALGBRC(6)/ALGBRC(5))
      RETURN
      END
      REAL FUNCTION TERNRY(TEST, VALA, VALB)
      LOGICAL TEST
      REAL VALA, VALB
      IF (TEST) THEN
        TERNRY = VALA
      ELSE
        TERNRY = VALB
      ENDIF
      RETURN
      END