Errors with ModelingToolkit.MultithreadedForm

Here are the codes: Parameters and ODE function


#----------------------------------------------
# Number of substances
#----------------------------------------------
const nsolid = 11 #  # number of solid substances
const ndissolved = 13 #  # number of dissolved substances
const nsummed = 4 #  # number of summed substances
const nnoneq = 24 #  # number of solid + dissolved substances
const nspec = 28 #  # number of total substances
const Lwbdwth = 31 #  # lower bandwidth of jacobian matrix
const Upbdwth = 55 #  # upper bandwidth of jacobian matrix

#----------------------------------------------
# global parameters
#----------------------------------------------
const depth = 500.0 # m # water depth
const salinity = 35.0 # psu # bottom water salinity
const temp = 5.0 # Celsius # bottom water temperature
const ds_rho = 2.6 # g cm^-3 # dry sediment density
const O2BW = 8.5e-6 # mmol/cm3 # bottom water oxygen
const sw_dens = 1.0287324258804407 # g cm^-3 # seawater density

#----------------------------------------------
# grid parameters
#----------------------------------------------
const L = 50.0 # cm # model sediment section thickness
const Ngrid = 100 # integer # number of model grid
const pgrid = L / 5 # cm # constant in gridtran, attenuation scale
const ξ = range(0, step = L / (Ngrid), length = Ngrid + 1) # cm # uniform grid
const xᵥ = broadcast(x -> L * (exp(x * pgrid / L) - 1) / (exp(pgrid) - 1), ξ) # cm # no grid transformation
const x = (xᵥ[2:(Ngrid+1)] .+ xᵥ[1:Ngrid]) / 2 # cm # cell center
const dx = xᵥ[2:(Ngrid+1)] .- xᵥ[1:Ngrid] # cm # cell volume

#----------------------------------------------
# porosity parameters
#----------------------------------------------
const phi0 = 0.9472 # dimentionless # surface porosity
const phiL = 0.7884 # dimentionless # bottom porosity
const xphi = 107.3 # cm # porosity attenuation scale
const phi_Inf = 0.7884 # dimentionless # porosity at infinite sediment depth (normally where porosity stops changing). Needed to calculate burial velocities. If constant_porosity_profile = yes, then phi_Inf should be the same as the porosity constant. If constant_porosity_profile = no, then phi_Inf should be consistent with the depth dependent porosity function
const phif = broadcast(x -> phiL + (phi0 - phiL) * exp(-x / xphi), x) # dimentionless # fluid volume fraction
const phis = 1.0 .- phif # dimentionless # solid volume fraction
const pwtods = phif ./ phis # dimentionless # conversion from pore water to solid sediment volume unit
const dstopw = phis ./ phif # dimentionless # conversion from solid sediment to pore water volume unit

#----------------------------------------------
# burial parameters
#----------------------------------------------
const Fsed = 0.073 # g cm^-2 yr^-1 # total sediment flux
const w_Inf = Fsed / ds_rho / (1 - phi_Inf) # cm yr^-1 # solid sediment burial velocity at infinite depth
const uf = phi_Inf * w_Inf ./ phif # cm yr^-1 # pore water burial velocity
const us = Fsed / ds_rho ./ phis # cm yr^-1 # solid sediment burial velocity

#----------------------------------------------
# bioturbation parameters
#----------------------------------------------
const Dbt0 = 5.2 * 10^(0.76241122 - 0.00039724 * depth) # cm/yr^3 # bioturbation coefficient
const aO2bt = 1.0e-5 # mmol/cm3 # constant used in calculating oxygen dependence of bioturbation
const bO2bt = 5.0e-7 # mmol/cm3 # constant used in calculating oxygen dependence of bioturbation
const fO2bt = 0.5 + 0.5 * erf((O2BW - aO2bt) / bO2bt) # missing # oxygen dependence of bioturbation
const xbt = 3.0 # cm # attentuation scale of bioturbation
const Ds = broadcast(x -> Dbt0 * fO2bt * exp(-x / xbt) + 1e-8, x) # cm^2 yr^-1 # Bioturbation coefficient

#----------------------------------------------
# bioirrigation parameters
#----------------------------------------------
const Dbir0 = 465.0 # yr^-1 # bioirrigation constant
const aO2bir = 1.0e-5 # mmol/cm3 # constant used in calculating oxygen dependence of bioirrigation
const bO2bir = 5.0e-7 # mmol/cm3 # constant used in calculating oxygen dependence of bioirrigation
const fO2bir = 0.5 + 0.5 * erf((O2BW - aO2bir) / bO2bir) # missing # oxygen dependence of irrigation
const xbir = 2.0 # cm # attentuation scale of bioirrigation
const alpha = broadcast(x -> Dbir0 * fO2bir * exp(-x / xbir), x) # cm^2 yr^-1 # Bioirrigation coefficient

#----------------------------------------------
# adsorption parameters
#----------------------------------------------
const KNH4_ads = 1.6 * ds_rho # cm^3(porewater) cm^-3(dry sediment) # Adsorption constant
const KMn_ads = 28.0 * ds_rho # cm^3(porewater) cm^-3(dry sediment) # Adsorption constant
const KFe_ads = 268.0 * ds_rho # cm^3(porewater) cm^-3(dry sediment) # Adsorption constant

#----------------------------------------------
# diffusion parameters
#----------------------------------------------
const DO2 = 4.2700687755264715E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DNO3 = 3.4158080786693472E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DMn = 1.1809021810434676E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DFe = 1.2122535663809046E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DCH4 = 3.2443708884445152E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DNO2 = 3.5695791591339196E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DCa = 1.3421378770645731E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DAl = 1.6944677313329149E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DMo = 1.6349371584775281E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DNH4 = 3.4531311564520092E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DH3PO4 = 1.2376332592731156E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DSO4 = 1.8034511184582917E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DH4SiO4 = 1.7771493723450564E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DH = 1.8564498889096735E+03 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DOH = 9.3665996003371845E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DCO2 = 3.3679572156139625E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DHCO3 = 1.9213920442515075E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DCO3 = 1.5899631135414575E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DH2S = 3.0488407320129431E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DHS = 3.5128480809042725E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DH3BO3 = 1.9631421965620916E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DH4BO4 = 1.7177494219918299E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient

#----------------------------------------------
# Boundary Condition parameters
#----------------------------------------------
const delta = 0.05 # cm # thickness of the diffusive boundary layer
const FPOC0 = 0.31 # mmol cm^-2 yr^-1 # Flux of POC at the  TOP of sediment column
const FMnO20 = 0.005 # mmol cm^-2 yr^-1 # Flux of MnO2 at the  TOP of sediment column
const FFeOOH0 = 0.0235 # mmol cm^-2 yr^-1 # Flux of FeOOH at the  TOP of sediment column
const FFeS0 = 2.22044604925031e-16 # mmol cm^-2 yr^-1 # Flux of FeS at the  TOP of sediment column
const FFeS20 = 0.002 # mmol cm^-2 yr^-1 # Flux of FeS2 at the  TOP of sediment column
const FCaCO30 = 0.055 # mmol cm^-2 yr^-1 # Flux of CaCO3 at the  TOP of sediment column
const FMnCO30 = 2.22044604925031e-16 # mmol cm^-2 yr^-1 # Flux of MnCO3 at the  TOP of sediment column
const FFeCO30 = 2.22044604925031e-16 # mmol cm^-2 yr^-1 # Flux of FeCO3 at the  TOP of sediment column
const FAge0 = 0.0 # mmol cm^-2 yr^-1 # Flux of Age at the  TOP of sediment column
const FBSi0 = 0.2 # mmol cm^-2 yr^-1 # Flux of BSi at the  TOP of sediment column
const FAnnite0 = Fsed * 0.1 / 511.88 * 1000 # mmol cm^-2 yr^-1 # Flux of Annite at the  TOP of sediment column
const FMoS40 = 365 / 1e6 / 95.94 # mmol cm^-2 yr^-1 # Flux of MoS4 at the  TOP of sediment column
const FBasalt0 = Fsed * 0.01 / 87.383 * 1000 # mmol cm^-2 yr^-2 # Flux of Basalt at the  TOP of sediment column
const O2BW = 8.5e-6 # mmol cm^-3 # Bottom water concentration of O2
const NO3BW = 2.8e-5 # mmol cm^-3 # Bottom water concentration of NO3
const Mn0 = 6.0e-7 # mmol cm^-3 # Concentration of Mn at the TOP of sediment column
const Fe0 = 2.95e-8 # mmol cm^-3 # Concentration of Fe at the TOP of sediment column
const CH4BW = 2.22044604925031e-16 # mmol cm^-3 # Bottom water concentration of CH4
const NO2BW = 2.0e-7 # mmol cm^-3 # Bottom water concentration of NO2
const CaBW = 0.01033 # mmol cm^-3 # Bottom water concentration of Ca
const AlBW = 1.0e-9 # mmol cm^-3 # Bottom water concentration of Al
const MoBW = 1.2e-7 # mmol cm^-3 # Bottom water concentration of Mo
const NH40 = 1.0e-6 # mmol cm^-3 # Concentration of NH4 at the TOP of sediment column
const H3PO4BW = 3.7e-6 # mmol cm^-3 # Bottom water concentration of H3PO4
const SO4BW = 0.028 # mmol cm^-3 # Bottom water concentration of SO4
const H4SiO4BW = 0.00022 # mmol cm^-3 # Bottom water concentration of H4SiO4
const NdBW = 3.0e-11 # mmol cm^-3 # Bottom water concentration of Nd
const eNdBW = -2.0 # missing # Bottom water eNd
const Nd144BW = 7.1394e-12 # mmol cm^-3 # Bottom water concentration of Nd144
const Nd143BW = 3.65919575165256e-12 # mmol cm^-3 # Bottom water concentration of Nd143
const pHBW = 7.62 # free pH scale # Bottom water pH
const TCO2BW = 0.002345 # mmol cm^-3 # Bottom water concentration of TCO2
const TH2SBW = 2.22044604925031e-16 # mmol cm^-3 # Bottom water concentration of TH2S
const TH3BO3BW = 8.7062e-5 # mmol cm^-3 # Bottom water concentration of TH3BO3
const CH4L = 9.0e-5 # mmol cm^-3 # Concentration of CH4 at the BOTTOM of sediment column
const CaL = 0.008 # mmol cm^-3 # Concentration of Ca at the BOTTOM of sediment column
const NH4L = 0.0013 # mmol cm^-3 # Concentration of NH4 at the BOTTOM of sediment column
const H3PO4L = 0.0001 # mmol cm^-3 # Concentration of H3PO4 at the BOTTOM of sediment column
const SO4L = 0.015 # mmol cm^-3 # Concentration of SO4 at the BOTTOM of sediment column
const H4SiO4L = 0.0005 # mmol cm^-3 # Concentration of H4SiO4 at the BOTTOM of sediment column
const pHL = 7.76210992440118 # free pH scale # pH at the BOTTOM of sediment column
const TCO2L = 0.014 # mmol cm^-3 # Concentration of TCO2 at the BOTTOM of sediment column
const TH2SL = 0.002 # mmol cm^-3 # Concentration of TH2S at the BOTTOM of sediment column
const TH3BO3L = 8.7062e-5 # mmol cm^-3 # Concentration of TH3BO3 at the BOTTOM of sediment column
const HBW = 2.3988329190194897E-08 # mmol cm^-3 # Bottom water concentration of H
const OHBW = 3.1707849164901275E-07 # mmol cm^-3 # Bottom water concentration of OH
const HL = 1.7293785803045845E-08 # mmol cm^-3 # Concentration of H at the BOTTOM of sediment column
const OHL = 4.3982175582788539E-07 # mmol cm^-3 # Concentration of OH at the BOTTOM of sediment column
const CO2BW = 6.7843776627964569E-05 # mmol cm^-3 # Bottom water concentration of CO2
const HCO3BW = 2.2363282164579099E-03 # mmol cm^-3 # Bottom water concentration of HCO3
const CO3BW = 4.0828006914125728E-05 # mmol cm^-3 # Bottom water concentration of CO3
const CO2L = 2.9239174203595165E-04 # mmol cm^-3 # Concentration of CO2 at the BOTTOM of sediment column
const HCO3L = 1.3369050230563931E-02 # mmol cm^-3 # Concentration of HCO3 at the BOTTOM of sediment column
const CO3L = 3.3855802740011910E-04 # mmol cm^-3 # Concentration of CO3 at the BOTTOM of sediment column
const H2SBW = 8.0198948904725221E-18 # mmol cm^-3 # Bottom water concentration of H2S
const HSBW = 2.1402471003455848E-16 # mmol cm^-3 # Bottom water concentration of HS
const H2SL = 5.2607578304958858E-05 # mmol cm^-3 # Concentration of H2S at the BOTTOM of sediment column
const HSL = 1.9473924216950411E-03 # mmol cm^-3 # Concentration of HS at the BOTTOM of sediment column
const H3BO3BW = 8.2637561283801204E-05 # mmol cm^-3 # Bottom water concentration of H3BO3
const H4BO4BW = 4.4244387161987993E-06 # mmol cm^-3 # Bottom water concentration of H4BO4
const H3BO3L = 8.1043235029232773E-05 # mmol cm^-3 # Concentration of H3BO3 at the BOTTOM of sediment column
const H4BO4L = 6.0187649707672218E-06 # mmol cm^-3 # Concentration of H4BO4 at the BOTTOM of sediment column
const Mn_ads0 = KMn_ads * Mn0 # mmol cm^-3 # Concentration of adsorbed Mn at the TOP of sediment column
const Fe_ads0 = KFe_ads * Fe0 # mmol cm^-3 # Concentration of adsorbed Fe at the TOP of sediment column
const NH4_ads0 = KNH4_ads * NH40 # mmol cm^-3 # Concentration of adsorbed NH4 at the TOP of sediment column
const NH4_adsL = KNH4_ads * NH4L # mmol cm^-3 # Concentration of adsorbed NH4 at the BOTTOM of sediment column
const betaO2 = 8.5401375510529433E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaNO3 = 6.8316161573386944E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaCH4 = 6.4887417768890300E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaNO2 = 7.1391583182678387E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaCa = 2.6842757541291462E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaAl = 3.3889354626658296E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaMo = 3.2698743169550562E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaH3PO4 = 2.4752665185462311E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaSO4 = 3.6069022369165832E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaH4SiO4 = 3.5542987446901125E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaH = 3.7128997778193465E+04 # cm yr^-1 # solute mass transfer velocity across SWI
const betaOH = 1.8733199200674368E+04 # cm yr^-1 # solute mass transfer velocity across SWI
const betaCO2 = 6.7359144312279250E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaHCO3 = 3.8427840885030150E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaCO3 = 3.1799262270829149E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaH2S = 6.0976814640258863E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaHS = 7.0256961618085443E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaH3BO3 = 3.9262843931241832E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaH4BO4 = 3.4354988439836598E+03 # cm yr^-1 # solute mass transfer velocity across SWI

# finite volume complete flux
@inline Afun(zᵢ, zᵢ₊₁, Peᵢ, Peᵢ₊₁) =
    zᵢ * zᵢ₊₁ /
    ((-exp(-Peᵢ / 2 - Peᵢ₊₁ / 2) + exp(-Peᵢ / 2)) * zᵢ + (1 - exp(-Peᵢ / 2)) * zᵢ₊₁)

@inline Bfunₗ(hᵢ, uᵢ, φᵢ, uᵢ₊₁, φᵢ₊₁, Peᵢ, Peᵢ₊₁) =
    hᵢ * (exp(-Peᵢ / 2.0) / (uᵢ * Peᵢ) + 1.0 / (2.0*uᵢ) - 1 / (uᵢ * Peᵢ)) / (
        (1 - exp(-Peᵢ / 2)) / (uᵢ * φᵢ) +
        (exp(-Peᵢ / 2) - exp(-Peᵢ / 2 - Peᵢ₊₁ / 2)) / (uᵢ₊₁ * φᵢ₊₁)
    )

@inline Bfunᵤ(hᵢ₊₁, uᵢ, φᵢ, uᵢ₊₁, φᵢ₊₁, Peᵢ, Peᵢ₊₁) =
    hᵢ₊₁ * (
        exp(-Peᵢ / 2) / (uᵢ₊₁ * Peᵢ₊₁) - exp(-Peᵢ / 2 - Peᵢ₊₁ / 2) / (uᵢ₊₁ * Peᵢ₊₁) -
        1 / 2 / (uᵢ₊₁) * exp(-Peᵢ / 2 - Peᵢ₊₁ / 2)
    ) / (
        (1 - exp(-Peᵢ / 2)) / (uᵢ * φᵢ) +
        (exp(-Peᵢ / 2) - exp(-Peᵢ / 2 - Peᵢ₊₁ / 2)) / (uᵢ₊₁ * φᵢ₊₁)
    )

function fvcf(φ, D, u, dx, N)
    Pe = u .* dx ./ D
    Aₗ = zeros(N - 1)
    Aᵤ = zeros(N - 1)
    Aₘ = zeros(N)
    Aₘₗ = zeros(N)
    Aₘᵤ = zeros(N)

    Bₗ = zeros(N - 1)
    Bᵤ = zeros(N - 1)
    Bₘ = zeros(N)
    Bₘₗ = zeros(N)
    Bₘᵤ = zeros(N)

    for i = 1:(N-1)
        Aₗ[i] = Afun(u[i]φ[i], u[i+1]φ[i+1], Pe[i], Pe[i+1]) / (dx[i+1]φ[i+1])
        Aₘₗ[i+1] =
            -Afun(u[i]φ[i], u[i+1]φ[i+1], Pe[i], Pe[i+1]) * exp(-Pe[i] / 2 - Pe[i+1] / 2) /
            (dx[i+1]φ[i+1])
        Aₘᵤ[i] = -Afun(u[i]φ[i], u[i+1]φ[i+1], Pe[i], Pe[i+1]) / (dx[i]φ[i])
        Aᵤ[i] =
            Afun(u[i]φ[i], u[i+1]φ[i+1], Pe[i], Pe[i+1]) * exp(-Pe[i] / 2 - Pe[i+1] / 2) /
            (dx[i]φ[i])

        Bₗ[i] = Bfunₗ(dx[i], u[i], φ[i], u[i+1], φ[i+1], Pe[i], Pe[i+1]) / (dx[i+1]φ[i+1])
        Bₘₗ[i+1] =
            -Bfunᵤ(dx[i+1], u[i], φ[i], u[i+1], φ[i+1], Pe[i], Pe[i+1]) / (dx[i+1]φ[i+1])
        Bₘᵤ[i] = -Bfunₗ(dx[i], u[i], φ[i], u[i+1], φ[i+1], Pe[i], Pe[i+1]) / (dx[i]φ[i])
        Bᵤ[i] = Bfunᵤ(dx[i+1], u[i], φ[i], u[i+1], φ[i+1], Pe[i], Pe[i+1]) / (dx[i]φ[i])
    end

    Aₘ .= Aₘₗ .+ Aₘᵤ
    A = Tridiagonal(Aₗ, Aₘ, Aᵤ)
    Bₘ .= Bₘₗ .+ Bₘᵤ #.+ 1.0
    B = Tridiagonal(Bₗ, Bₘ, Bᵤ)

    return (A, B)
end

function fvcf_bc(φ, D, u, dx, BC::Tuple{Tuple,Tuple}, N, ads = false)
    A_bc = zeros(2)
    B_bc = zeros(2)
    b_bc = zeros(2)

    if ads
        A_bc[1] = u[1]φ[1] / (dx[1]φ[1])
        A_bc[2] = -u[N]φ[N] / (dx[N]φ[N])
        b_bc[1] = 0
        b_bc[2] = 0
    else
        α⁰, β⁰, γ⁰ = BC[1]
        αᴸ, βᴸ, γᴸ = BC[2]

        Pe = u .* dx ./ D

        A_bc[1] =
            u[1]φ[1] * exp(-Pe[1] / 2) * (α⁰ * D[1] + β⁰ * u[1]) /
            (α⁰ * D[1] * exp(-Pe[1] / 2) - α⁰ * D[1] + β⁰ * u[1] * exp(-Pe[1] / 2)) /
            (dx[1]φ[1])
        A_bc[2] =
            u[N]φ[N] * (αᴸ * D[N] + βᴸ * u[N]) /
            (αᴸ * D[N] * exp(-Pe[N] / 2) - αᴸ * D[N] - βᴸ * u[N]) / (dx[N]φ[N])

        B_bc[1] = 0.0
        B_bc[2] =
            (αᴸ * D[N] + βᴸ * u[N]) *
            (exp(-Pe[N] / 2) / Pe[N] - 1 / Pe[N] + 1 / 2) /
            (αᴸ * D[N] * exp(-Pe[N] / 2) - αᴸ * D[N] - βᴸ * u[N])

        b_bc[1] =
            φ[1]u[1]D[1] * γ⁰ /
            (-α⁰ * D[1] * exp(-Pe[1] / 2) + α⁰ * D[1] - β⁰ * u[1] * exp(-Pe[1] / 2)) /
            (dx[1]φ[1])
        b_bc[2] =
            φ[N]u[N]D[N] * γᴸ * exp(-Pe[N] / 2) /
            (-αᴸ * D[N] * exp(-Pe[N] / 2) + αᴸ * D[N] + βᴸ * u[N]) / (dx[N]φ[N])
    end
    return (A_bc, B_bc, b_bc)
end

const BcPOC = ((phis[1]us[1], -phis[1]Ds[1], FPOC0), (0.0, 1.0, 0.0)) #  # Boundary condition of POC
const BcMnO2 = ((phis[1]us[1], -phis[1]Ds[1], FMnO20), (0.0, 1.0, 0.0)) #  # Boundary condition of MnO2
const BcFeOOH = ((phis[1]us[1], -phis[1]Ds[1], FFeOOH0), (0.0, 1.0, 0.0)) #  # Boundary condition of FeOOH
const BcFeS = ((phis[1]us[1], -phis[1]Ds[1], FFeS0), (0.0, 1.0, 0.0)) #  # Boundary condition of FeS
const BcFeS2 = ((phis[1]us[1], -phis[1]Ds[1], FFeS20), (0.0, 1.0, 0.0)) #  # Boundary condition of FeS2
const BcCaCO3 = ((phis[1]us[1], -phis[1]Ds[1], FCaCO30), (0.0, 1.0, 0.0)) #  # Boundary condition of CaCO3
const BcMnCO3 = ((phis[1]us[1], -phis[1]Ds[1], FMnCO30), (0.0, 1.0, 0.0)) #  # Boundary condition of MnCO3
const BcFeCO3 = ((phis[1]us[1], -phis[1]Ds[1], FFeCO30), (0.0, 1.0, 0.0)) #  # Boundary condition of FeCO3
const BcAge =
    ((phis[1]us[1], -phis[1]Ds[1], FAge0), (0.0, 1.0, 1.0 / us[Ngrid])) #  # Boundary condition of Age
const BcBSi = ((phis[1]us[1], -phis[1]Ds[1], FBSi0), (0.0, 1.0, 0.0)) #  # Boundary condition of BSi
const BcMoS4 = ((phis[1]us[1], -phis[1]Ds[1], FMoS40), (0.0, 1.0, 0.0)) #  # Boundary condition of MoS4
const BcO2 =
    ((betaO2 + phif[1]uf[1], -phif[1]DO2[1], betaO2 * O2BW), (0.0, 1.0, 0.0)) #  # Boundary condition of O2
const BcNO3 = (
    (betaNO3 + phif[1]uf[1], -phif[1]DNO3[1], betaNO3 * NO3BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of NO3
const BcMn = ((1.0, 0.0, Mn0), (0.0, 1.0, 0.0)) #  # Boundary condition of Mn
const BcFe = ((1.0, 0.0, Fe0), (0.0, 1.0, 0.0)) #  # Boundary condition of Fe
const BcCH4 = (
    (betaCH4 + phif[1]uf[1], -phif[1]DCH4[1], betaCH4 * CH4BW),
    (1.0, 0.0, CH4L),
) #  # Boundary condition of CH4
const BcNO2 = (
    (betaNO2 + phif[1]uf[1], -phif[1]DNO2[1], betaNO2 * NO2BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of NO2
const BcCa =
    ((betaCa + phif[1]uf[1], -phif[1]DCa[1], betaCa * CaBW), (1.0, 0.0, CaL)) #  # Boundary condition of Ca
const BcAl =
    ((betaAl + phif[1]uf[1], -phif[1]DAl[1], betaAl * AlBW), (0.0, 1.0, 0.0)) #  # Boundary condition of Al
const BcMo =
    ((betaMo + phif[1]uf[1], -phif[1]DMo[1], betaMo * MoBW), (0.0, 1.0, 0.0)) #  # Boundary condition of Mo
const BcNH4 = ((1.0, 0.0, NH40), (1.0, 0.0, NH4L)) #  # Boundary condition of NH4
const BcH3PO4 = (
    (betaH3PO4 + phif[1]uf[1], -phif[1]DH3PO4[1], betaH3PO4 * H3PO4BW),
    (1.0, 0.0, H3PO4L),
) #  # Boundary condition of H3PO4
const BcSO4 = (
    (betaSO4 + phif[1]uf[1], -phif[1]DSO4[1], betaSO4 * SO4BW),
    (1.0, 0.0, SO4L),
) #  # Boundary condition of SO4
const BcH4SiO4 = (
    (betaH4SiO4 + phif[1]uf[1], -phif[1]DH4SiO4[1], betaH4SiO4 * H4SiO4BW),
    (1.0, 0.0, H4SiO4L),
) #  # Boundary condition of H4SiO4
const BcH = ((betaH + phif[1]uf[1], -phif[1]DH[1], betaH * HBW), (1.0, 0.0, HL)) #  # Boundary condition of H
const BcOH =
    ((betaOH + phif[1]uf[1], -phif[1]DOH[1], betaOH * OHBW), (1.0, 0.0, OHL)) #  # Boundary condition of H
const BcCO2 = (
    (betaCO2 + phif[1]uf[1], -phif[1]DCO2[1], betaCO2 * CO2BW),
    (1.0, 0.0, CO2L),
) #  # Boundary condition of TCO2
const BcHCO3 = (
    (betaHCO3 + phif[1]uf[1], -phif[1]DHCO3[1], betaHCO3 * HCO3BW),
    (1.0, 0.0, HCO3L),
) #  # Boundary condition of TCO2
const BcCO3 = (
    (betaCO3 + phif[1]uf[1], -phif[1]DCO3[1], betaCO3 * CO3BW),
    (1.0, 0.0, CO3L),
) #  # Boundary condition of TCO2
const BcH2S = (
    (betaH2S + phif[1]uf[1], -phif[1]DH2S[1], betaH2S * H2SBW),
    (1.0, 0.0, H2SL),
) #  # Boundary condition of TH2S
const BcHS =
    ((betaHS + phif[1]uf[1], -phif[1]DHS[1], betaHS * HSBW), (1.0, 0.0, HSL)) #  # Boundary condition of TH2S
const BcH3BO3 = (
    (betaH3BO3 + phif[1]uf[1], -phif[1]DH3BO3[1], betaH3BO3 * H3BO3BW),
    (1.0, 0.0, H3BO3L),
) #  # Boundary condition of TH3BO3
const BcH4BO4 = (
    (betaH4BO4 + phif[1]uf[1], -phif[1]DH4BO4[1], betaH4BO4 * H4BO4BW),
    (1.0, 0.0, H4BO4L),
) #  # Boundary condition of TH3BO3
const BcMn_ads = ((1.0, 0.0, Mn_ads0), (0.0, 1.0, 0.0)) #  # Boundary condition of Mn
const BcFe_ads = ((1.0, 0.0, Fe_ads0), (0.0, 1.0, 0.0)) #  # Boundary condition of Fe
const BcNH4_ads = ((1.0, 0.0, NH4_ads0), (1.0, 0.0, NH4_adsL)) #  # Boundary condition of NH4

#----------------------------------------------
# Boundary transport matrix
#----------------------------------------------
const BcAmPOC, BcBmPOC, BcCmPOC = fvcf_bc(phis, Ds, us, dx, BcPOC, Ngrid) #  # Boundary transport matrix of POC
const BcAmMnO2, BcBmMnO2, BcCmMnO2 = fvcf_bc(phis, Ds, us, dx, BcMnO2, Ngrid) #  # Boundary transport matrix of MnO2
const BcAmFeOOH, BcBmFeOOH, BcCmFeOOH =
    fvcf_bc(phis, Ds, us, dx, BcFeOOH, Ngrid) #  # Boundary transport matrix of FeOOH
const BcAmFeS, BcBmFeS, BcCmFeS = fvcf_bc(phis, Ds, us, dx, BcFeS, Ngrid) #  # Boundary transport matrix of FeS
const BcAmFeS2, BcBmFeS2, BcCmFeS2 = fvcf_bc(phis, Ds, us, dx, BcFeS2, Ngrid) #  # Boundary transport matrix of FeS2
const BcAmCaCO3, BcBmCaCO3, BcCmCaCO3 =
    fvcf_bc(phis, Ds, us, dx, BcCaCO3, Ngrid) #  # Boundary transport matrix of CaCO3
const BcAmMnCO3, BcBmMnCO3, BcCmMnCO3 =
    fvcf_bc(phis, Ds, us, dx, BcMnCO3, Ngrid) #  # Boundary transport matrix of MnCO3
const BcAmFeCO3, BcBmFeCO3, BcCmFeCO3 =
    fvcf_bc(phis, Ds, us, dx, BcFeCO3, Ngrid) #  # Boundary transport matrix of FeCO3
const BcAmAge, BcBmAge, BcCmAge = fvcf_bc(phis, Ds, us, dx, BcAge, Ngrid) #  # Boundary transport matrix of Age
const BcAmBSi, BcBmBSi, BcCmBSi = fvcf_bc(phis, Ds, us, dx, BcBSi, Ngrid) #  # Boundary transport matrix of BSi
const BcAmMoS4, BcBmMoS4, BcCmMoS4 = fvcf_bc(phis, Ds, us, dx, BcMoS4, Ngrid) #  # Boundary transport matrix of MoS4
const BcAmO2, BcBmO2, BcCmO2 = fvcf_bc(phif, DO2, uf, dx, BcO2, Ngrid) #  # Boundary transport matrix of O2
const BcAmNO3, BcBmNO3, BcCmNO3 = fvcf_bc(phif, DNO3, uf, dx, BcNO3, Ngrid) #  # Boundary transport matrix of NO3
const BcAmMn, BcBmMn, BcCmMn = fvcf_bc(phif, DMn, uf, dx, BcMn, Ngrid) #  # Boundary transport matrix of Mn
const BcAmFe, BcBmFe, BcCmFe = fvcf_bc(phif, DFe, uf, dx, BcFe, Ngrid) #  # Boundary transport matrix of Fe
const BcAmCH4, BcBmCH4, BcCmCH4 = fvcf_bc(phif, DCH4, uf, dx, BcCH4, Ngrid) #  # Boundary transport matrix of CH4
const BcAmNO2, BcBmNO2, BcCmNO2 = fvcf_bc(phif, DNO2, uf, dx, BcNO2, Ngrid) #  # Boundary transport matrix of NO2
const BcAmCa, BcBmCa, BcCmCa = fvcf_bc(phif, DCa, uf, dx, BcCa, Ngrid) #  # Boundary transport matrix of Ca
const BcAmAl, BcBmAl, BcCmAl = fvcf_bc(phif, DAl, uf, dx, BcAl, Ngrid) #  # Boundary transport matrix of Al
const BcAmMo, BcBmMo, BcCmMo = fvcf_bc(phif, DMo, uf, dx, BcMo, Ngrid) #  # Boundary transport matrix of Mo
const BcAmNH4, BcBmNH4, BcCmNH4 = fvcf_bc(phif, DNH4, uf, dx, BcNH4, Ngrid) #  # Boundary transport matrix of NH4
const BcAmH3PO4, BcBmH3PO4, BcCmH3PO4 =
    fvcf_bc(phif, DH3PO4, uf, dx, BcH3PO4, Ngrid) #  # Boundary transport matrix of H3PO4
const BcAmSO4, BcBmSO4, BcCmSO4 = fvcf_bc(phif, DSO4, uf, dx, BcSO4, Ngrid) #  # Boundary transport matrix of SO4
const BcAmH4SiO4, BcBmH4SiO4, BcCmH4SiO4 =
    fvcf_bc(phif, DH4SiO4, uf, dx, BcH4SiO4, Ngrid) #  # Boundary transport matrix of H4SiO4
const BcAmH, BcBmH, BcCmH = fvcf_bc(phif, DH, uf, dx, BcH, Ngrid) #  # Boundary transport matrix of H
const BcAmOH, BcBmOH, BcCmOH = fvcf_bc(phif, DOH, uf, dx, BcOH, Ngrid) #  # Boundary transport matrix of OH
const BcAmCO2, BcBmCO2, BcCmCO2 = fvcf_bc(phif, DCO2, uf, dx, BcCO2, Ngrid) #  # Boundary transport matrix of CO2
const BcAmHCO3, BcBmHCO3, BcCmHCO3 = fvcf_bc(phif, DHCO3, uf, dx, BcHCO3, Ngrid) #  # Boundary transport matrix of HCO3
const BcAmCO3, BcBmCO3, BcCmCO3 = fvcf_bc(phif, DCO3, uf, dx, BcCO3, Ngrid) #  # Boundary transport matrix of CO3
const BcAmH2S, BcBmH2S, BcCmH2S = fvcf_bc(phif, DH2S, uf, dx, BcH2S, Ngrid) #  # Boundary transport matrix of H2S
const BcAmHS, BcBmHS, BcCmHS = fvcf_bc(phif, DHS, uf, dx, BcHS, Ngrid) #  # Boundary transport matrix of HS
const BcAmH3BO3, BcBmH3BO3, BcCmH3BO3 =
    fvcf_bc(phif, DH3BO3, uf, dx, BcH3BO3, Ngrid) #  # Boundary transport matrix of H3BO3
const BcAmH4BO4, BcBmH4BO4, BcCmH4BO4 =
    fvcf_bc(phif, DH4BO4, uf, dx, BcH4BO4, Ngrid) #  # Boundary transport matrix of H4BO4
const BcAmMn_ads, BcBmMn_ads, BcCmMn_ads =
    fvcf_bc(phis, Ds, us, dx, BcMn_ads, Ngrid) #  # Boundary transport matrix of Mn_ads
const BcAmFe_ads, BcBmFe_ads, BcCmFe_ads =
    fvcf_bc(phis, Ds, us, dx, BcFe_ads, Ngrid) #  # Boundary transport matrix of Fe_ads
const BcAmNH4_ads, BcBmNH4_ads, BcCmNH4_ads =
    fvcf_bc(phis, Ds, us, dx, BcNH4_ads, Ngrid) #  # Boundary transport matrix of NH4_ads

#----------------------------------------------
# Interior transport matrix
#----------------------------------------------
const AmPOC, BmPOC = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of POC
const AmMnO2, BmMnO2 = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of MnO2
const AmFeOOH, BmFeOOH = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of FeOOH
const AmFeS, BmFeS = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of FeS
const AmFeS2, BmFeS2 = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of FeS2
const AmCaCO3, BmCaCO3 = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of CaCO3
const AmMnCO3, BmMnCO3 = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of MnCO3
const AmFeCO3, BmFeCO3 = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of FeCO3
const AmAge, BmAge = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of Age
const AmBSi, BmBSi = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of BSi
const AmMoS4, BmMoS4 = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of MoS4
const AmO2, BmO2 = fvcf(phif, DO2, uf, dx, Ngrid) #  # Interior transport matrix of O2
const AmNO3, BmNO3 = fvcf(phif, DNO3, uf, dx, Ngrid) #  # Interior transport matrix of NO3
const AmMn, BmMn = fvcf(phif, DMn, uf, dx, Ngrid) #  # Interior transport matrix of Mn
const AmFe, BmFe = fvcf(phif, DFe, uf, dx, Ngrid) #  # Interior transport matrix of Fe
const AmCH4, BmCH4 = fvcf(phif, DCH4, uf, dx, Ngrid) #  # Interior transport matrix of CH4
const AmNO2, BmNO2 = fvcf(phif, DNO2, uf, dx, Ngrid) #  # Interior transport matrix of NO2
const AmCa, BmCa = fvcf(phif, DCa, uf, dx, Ngrid) #  # Interior transport matrix of Ca
const AmAl, BmAl = fvcf(phif, DAl, uf, dx, Ngrid) #  # Interior transport matrix of Al
const AmMo, BmMo = fvcf(phif, DMo, uf, dx, Ngrid) #  # Interior transport matrix of Mo
const AmNH4, BmNH4 = fvcf(phif, DNH4, uf, dx, Ngrid) #  # Interior transport matrix of NH4
const AmH3PO4, BmH3PO4 = fvcf(phif, DH3PO4, uf, dx, Ngrid) #  # Interior transport matrix of H3PO4
const AmSO4, BmSO4 = fvcf(phif, DSO4, uf, dx, Ngrid) #  # Interior transport matrix of SO4
const AmH4SiO4, BmH4SiO4 = fvcf(phif, DH4SiO4, uf, dx, Ngrid) #  # Interior transport matrix of H4SiO4
const AmH, BmH = fvcf(phif, DH, uf, dx, Ngrid) #  # Interior transport matrix of H
const AmOH, BmOH = fvcf(phif, DOH, uf, dx, Ngrid) #  # Interior transport matrix of OH
const AmCO2, BmCO2 = fvcf(phif, DCO2, uf, dx, Ngrid) #  # Interior transport matrix of CO2
const AmHCO3, BmHCO3 = fvcf(phif, DHCO3, uf, dx, Ngrid) #  # Interior transport matrix of HCO3
const AmCO3, BmCO3 = fvcf(phif, DCO3, uf, dx, Ngrid) #  # Interior transport matrix of CO3
const AmH2S, BmH2S = fvcf(phif, DH2S, uf, dx, Ngrid) #  # Interior transport matrix of H2S
const AmHS, BmHS = fvcf(phif, DHS, uf, dx, Ngrid) #  # Interior transport matrix of HS
const AmH3BO3, BmH3BO3 = fvcf(phif, DH3BO3, uf, dx, Ngrid) #  # Interior transport matrix of H3BO3
const AmH4BO4, BmH4BO4 = fvcf(phif, DH4BO4, uf, dx, Ngrid) #  # Interior transport matrix of H4BO4
const AmMn_ads, BmMn_ads = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of Mn_ads
const AmFe_ads, BmFe_ads = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of Fe_ads
const AmNH4_ads, BmNH4_ads = fvcf(phis, Ds, us, dx, Ngrid) #  # Interior transport matrix of NH4_ads

#----------------------------------------------
# Acid dissociation constants
#----------------------------------------------
const KH2O = 7.6061832368069808E-15 # mmol cm^-3 # Water dissociation constant
const KCO2 = 7.9072510553018769E-07 # mmol cm^-3 # H2CO3 dissociation constant
const KHCO3 = 4.3794808956390707E-10 # mmol cm^-3 # HCO3 dissociation constant
const KH2S = 6.4016988617197496E-07 # mmol cm^-3 # H2S dissociation constant
const KH3BO3 = 1.2843420202288193E-09 # mmol cm^-3 # H3BO3 dissociation constant

#----------------------------------------------
# Reaction parameters
#----------------------------------------------
const KO2 = 1.0e-6 #  mmol cm-3 pw yr-1 # Monod constant
const nu = 0.12 # dimentionless # POC reactivity
const a = 3.2 # yr # initial POC age
const rNC = 0.125 # mol mol^-1 # N/C ratio Sediment trap 
const rPC = 0.005 # mol mol^-2 # P/C ratio Sediment trap
const KNO2 = 1.0e-5 # mmol cm-3 pw yr-1 # Monod constant
const KNO3 = 1.0e-5 # mmol cm-3 pw yr-1 # Monod constant
const KMnO2 = 16.0 * ds_rho / 1e3 # mmol cm-3 ds # Monod constant
const KFeOOH = 500.0 * ds_rho / 1e3 # mmol cm-3 ds # Monod constant
const KSO4 = 0.001 # mmol cm-3 pw yr-1 # Monod constant
const kO2NO2 = 1.0e7 # (mmol cm-3 pw)-1 yr-1 #  
const kO2NH4 = 1.0e7 # (mmol cm-3 pw)-1 yr-1 #  
const kO2Mn = 5.0e6 # (mmol cm-3 pw)-1 yr-1  #  
const kO2Mn_ads = 1.0e7 # (mmol cm-3 pw)-1 yr-1  #  
const kO2Fe = 5.0e8 # (mmol cm-3 pw)-1 yr-1 #  
const kO2Fe_ads = 5.0e6 # (mmol cm-3 pw)-1 yr-1  #  
const kO2H2S = 100000.0 # (mmol cm-3 pw)-1 yr-1  #  
const kO2FeS = 100000.0 # (mmol cm-3 pw)-1 yr-1  #  
const kO2CH4 = 1.0e10 # (mmol cm-3 pw)-1 yr-1  #  
const kNO2NH4 = 1.0e8 # (mmol cm-3 pw)-1 yr-1  #  
const kNO3Mn = 100000.0 # (mmol cm-3 pw)-1 yr-1  #  
const kNO3Fe = 150000.0 # (mmol cm-3 pw)-1 yr-1  #  
const kNO3H2S = 0.0 # (mmol cm-3 pw)-1 yr-1 #  
const kAOM = 0.04 # yr-1  #  
const KAOM = 0.001 # mmol cm-3 pw yr-1  #  
const kMnO2Fe = 1.0e7 # (mmol cm-3 pw)-1 yr-1  #  
const kMnO2H2S = 100000.0 # (mmol cm-3 pw)-1 yr-1  #  
const kFeOOHH2S = 4000.0 # (mmol cm-3 pw)-0.5 yr-1  #  
const kFeSH2S = 300.0 # (mmol cm-3 pw)-1 yr-1  #  
const kFeSdis = 0.001 # yr-1  #  
const KspFeS = 10^(-2.2) # (mmol cm-3 pw)^-1  # apparent solubility of FeS  
const kFeSpre = 200e-3 * ds_rho # mmol cm-3 ds yr-1  #  
const kCaCO3dis = 0.02 # yr^-1  # missing
const KspCaCO3_dis = 4.49609622396569E-07 * sw_dens^2 # (mmol cm^-3 pw)^2 # missing
const kCaCO3pre = 9.0e-5 # mmol cm^-3 ds yr^-1  # missing
const KspCaCO3_pre = 4.49609622396569E-07 * sw_dens^2 # (mmol cm^-3 pw)^2 # missing
const kMnCO3dis = 0.8 # yr-1  # missing
const KspMnCO3 = 10^(-9.5) * sw_dens^2 # (mmol cm-3 pw)^-2  # missing
const kMnCO3pre = 0.2 * ds_rho # mmol cm-3 ds yr-1  # missing
const kFeCO3dis = 0.5 # yr-1  # missing
const KspFeCO3 = 10^(-10.9 + 2.518 * 0.7^0.5 - 0.657 * 0.7) * sw_dens^2 # (mmol cm-3 pw)^-2   # missing
const kFeCO3pre = 0.0001 # mmol cm-3 ds yr-1  # missing
const kBSidis0 = 5e-6 * 365 * 24 # yr-1 # opal dissolution rate
const xBSi = 1.7 # cm # dissolution attenuation scale
const kBSidis = broadcast(x -> kBSidis0 * exp(-x / xBSi), x) # missing # depth dependent dissolution rate
const H4SiO4_dis_sat =
    2.754e6 * exp(
        1 / (temp + 273.15) * (
            -2229 - 3.688e-3 * (1500 - 22 * 60) + 0.20 * (depth / 10) -
            2.7e-4 * (depth / 10)^2 + 1.46e-7 * (depth / 10)^3
        ),
    ) # mmol cm-3 pw # solubility of opal
const kASipre = 5e-6 * 365 * 24 # yr-1  # authigenic silicate precipitation rate
const H4SiO4_pre_sat = 0.0002 # mmol cm-3 pw # authigenic silicate precipitation threshold
const SAnnite = 1.7 #  m2/g  # surface area
const MAnnite = 511.88 #  g/mol # missing
const EaAnnite = 49000.0 #  J/mol  # activation energy
const kAnnite_0 =
    1.9e-12 *
    365 *
    24 *
    3600 *
    exp(-EaAnnite / 8.314 * (1.0 / (273.15 + temp) - 1.0 / 298.15)) # mol m^-2 yr^-1 # Annite dissolution rate
const kAnnite = kAnnite_0 * SAnnite * MAnnite # yr^-1 # missing
const a_s0 = 1e6 # yr  # initial age of Annite
const pl = 0.6 #  exponent of age dependent silicate weathering rate R=R_0*t^p  # missing
const KBW = 0.010578038 # mmol cm-3 # missing
const KAnnite = 10^(39.35134272) #  apparen solubility # missing
const kMoS4_pre = 0.5e5 # (mmol cm-3)-1 yr-1 pw  # missing
const TH2S_Mo_pre = 0.1e-6 # mmol cm-3  # hreshold for MoS4 precipitation 
const Cl = 0.565772678 # mmol cm-3 # bottom water Cl concentration
const rNdSi = 8.37018234241649e-6 # dimentionless (mol/mol) # Nd:Si ratio in oceanic arc basalt
const rNdnrSi = 0.23798 * rNdSi # missing # missing
const eNd_Basalt = 10.0 # missing # missing
const rNdrSi = rNdnrSi * (eNd_Basalt / 1e4 + 1) * 0.512638 # missing # missing
const KspBasalt = 41.63251909642741 # missing # missing
const SABasalt = 250.0 # cm2 g^-1 # missing
const Mbasalt = 87.383 # g/mol # missing
const EaBasalt = 25500.0 # J mol^-1 # missing
const kBasalt_0 =
    10^(-5.6) *
    365 *
    24 *
    3600 *
    exp(-EaBasalt / 8.314 / (273.15 + temp)) *
    0.75 / 0.076^0.33 # mol cm^-2 yr^-1 # missing
const kBasalt = kBasalt_0 * SABasalt * Mbasalt # yr^-1 # missing
const rNdMn = 0.00038089295618413754 # missing # missing
const rNdnrMn = 0.23798 * rNdMn # missing # missing
const eNd_MnO2 = -2.5 # missing # missing
const rNdrMn = rNdnrMn * (eNd_MnO2 / 1e4 + 1) * 0.512638 # missing # missing
const rNdFe = 0.0003882418191902385 # missing # missing
const rNdnrFe = 0.23798 * rNdFe # missing # missing
const eNd_FeOOH = -2.5 # missing # missing
const rNdrFe = rNdnrFe * (eNd_FeOOH / 1e4 + 1) * 0.512638 # missing # missing

#----------------------------------------------
# Reaction parameters
#----------------------------------------------
const C_ini = [
    FPOC0 / (phis[1] * us[1]),
    FMnO20 / (phis[1] * us[1]),
    FFeOOH0 / (phis[1] * us[1]),
    FFeS0 / (phis[1] * us[1]),
    FFeS20 / (phis[1] * us[1]),
    FCaCO30 / (phis[1] * us[1]),
    FMnCO30 / (phis[1] * us[1]),
    FFeCO30 / (phis[1] * us[1]),
    FAge0 / (phis[1] * us[1]),
    FBSi0 / (phis[1] * us[1]),
    FMoS40 / (phis[1] * us[1]),
    O2BW,
    NO3BW,
    Mn0,
    Fe0,
    CH4BW,
    NO2BW,
    CaBW,
    AlBW,
    MoBW,
    NH40,
    H3PO4BW,
    SO4BW,
    H4SiO4BW,
    HBW,
    TCO2BW,
    TH2SBW,
    TH3BO3BW,
] #  # initial conditions
const C_uni = repeat(C_ini, Ngrid) #  # initial conditions

#----------------------------------------------
# Tolerance parameters
#----------------------------------------------

#----------------------------------------------
# Indices
#----------------------------------------------
const POCID = ((1:Ngrid) .- 1)nspec .+ 1 #  # index
const MnO2ID = ((1:Ngrid) .- 1)nspec .+ 2 #  # index
const FeOOHID = ((1:Ngrid) .- 1)nspec .+ 3 #  # index
const FeSID = ((1:Ngrid) .- 1)nspec .+ 4 #  # index
const FeS2ID = ((1:Ngrid) .- 1)nspec .+ 5 #  # index
const CaCO3ID = ((1:Ngrid) .- 1)nspec .+ 6 #  # index
const MnCO3ID = ((1:Ngrid) .- 1)nspec .+ 7 #  # index
const FeCO3ID = ((1:Ngrid) .- 1)nspec .+ 8 #  # index
const AgeID = ((1:Ngrid) .- 1)nspec .+ 9 #  # index
const BSiID = ((1:Ngrid) .- 1)nspec .+ 10 #  # index
const MoS4ID = ((1:Ngrid) .- 1)nspec .+ 11 #  # index
const O2ID = ((1:Ngrid) .- 1)nspec .+ 12 #  # index
const NO3ID = ((1:Ngrid) .- 1)nspec .+ 13 #  # index
const MnID = ((1:Ngrid) .- 1)nspec .+ 14 #  # index
const FeID = ((1:Ngrid) .- 1)nspec .+ 15 #  # index
const CH4ID = ((1:Ngrid) .- 1)nspec .+ 16 #  # index
const NO2ID = ((1:Ngrid) .- 1)nspec .+ 17 #  # index
const CaID = ((1:Ngrid) .- 1)nspec .+ 18 #  # index
const AlID = ((1:Ngrid) .- 1)nspec .+ 19 #  # index
const MoID = ((1:Ngrid) .- 1)nspec .+ 20 #  # index
const NH4ID = ((1:Ngrid) .- 1)nspec .+ 21 #  # index
const H3PO4ID = ((1:Ngrid) .- 1)nspec .+ 22 #  # index
const SO4ID = ((1:Ngrid) .- 1)nspec .+ 23 #  # index
const H4SiO4ID = ((1:Ngrid) .- 1)nspec .+ 24 #  # index
const HID = ((1:Ngrid) .- 1)nspec .+ 25 #  # index
const TCO2ID = ((1:Ngrid) .- 1)nspec .+ 26 #  # index
const TH2SID = ((1:Ngrid) .- 1)nspec .+ 27 #  # index
const TH3BO3ID = ((1:Ngrid) .- 1)nspec .+ 28 #  # index

function reactran_fvcf_auto(dC, C, parms, t)
    POC = @view C[POCID]
    MnO2 = @view C[MnO2ID]
    FeOOH = @view C[FeOOHID]
    FeS = @view C[FeSID]
    FeS2 = @view C[FeS2ID]
    CaCO3 = @view C[CaCO3ID]
    MnCO3 = @view C[MnCO3ID]
    FeCO3 = @view C[FeCO3ID]
    Age = @view C[AgeID]
    BSi = @view C[BSiID]
    MoS4 = @view C[MoS4ID]
    O2 = @view C[O2ID]
    NO3 = @view C[NO3ID]
    Mn = @view C[MnID]
    Fe = @view C[FeID]
    CH4 = @view C[CH4ID]
    NO2 = @view C[NO2ID]
    Ca = @view C[CaID]
    Al = @view C[AlID]
    Mo = @view C[MoID]
    NH4 = @view C[NH4ID]
    H3PO4 = @view C[H3PO4ID]
    SO4 = @view C[SO4ID]
    H4SiO4 = @view C[H4SiO4ID]
    H = @view C[HID]
    TCO2 = @view C[TCO2ID]
    TH2S = @view C[TH2SID]
    TH3BO3 = @view C[TH3BO3ID]

    dPOC = AmPOC * POC
    dMnO2 = AmMnO2 * MnO2
    dFeOOH = AmFeOOH * FeOOH
    dFeS = AmFeS * FeS
    dFeS2 = AmFeS2 * FeS2
    dCaCO3 = AmCaCO3 * CaCO3
    dMnCO3 = AmMnCO3 * MnCO3
    dFeCO3 = AmFeCO3 * FeCO3
    dAge = AmAge * Age
    dBSi = AmBSi * BSi
    dMoS4 = AmMoS4 * MoS4
    dO2 = AmO2 * O2
    dNO3 = AmNO3 * NO3
    dCH4 = AmCH4 * CH4
    dNO2 = AmNO2 * NO2
    dCa = AmCa * Ca
    dAl = AmAl * Al
    dMo = AmMo * Mo
    dH3PO4 = AmH3PO4 * H3PO4
    dSO4 = AmSO4 * SO4
    dH4SiO4 = AmH4SiO4 * H4SiO4

    dPOC[1] += BcAmPOC[1] * POC[1] + BcCmPOC[1]
    dMnO2[1] += BcAmMnO2[1] * MnO2[1] + BcCmMnO2[1]
    dFeOOH[1] += BcAmFeOOH[1] * FeOOH[1] + BcCmFeOOH[1]
    dFeS[1] += BcAmFeS[1] * FeS[1] + BcCmFeS[1]
    dFeS2[1] += BcAmFeS2[1] * FeS2[1] + BcCmFeS2[1]
    dCaCO3[1] += BcAmCaCO3[1] * CaCO3[1] + BcCmCaCO3[1]
    dMnCO3[1] += BcAmMnCO3[1] * MnCO3[1] + BcCmMnCO3[1]
    dFeCO3[1] += BcAmFeCO3[1] * FeCO3[1] + BcCmFeCO3[1]
    dAge[1] += BcAmAge[1] * Age[1] + BcCmAge[1]
    dBSi[1] += BcAmBSi[1] * BSi[1] + BcCmBSi[1]
    dMoS4[1] += BcAmMoS4[1] * MoS4[1] + BcCmMoS4[1]
    dO2[1] += BcAmO2[1] * O2[1] + BcCmO2[1]
    dNO3[1] += BcAmNO3[1] * NO3[1] + BcCmNO3[1]
    dCH4[1] += BcAmCH4[1] * CH4[1] + BcCmCH4[1]
    dNO2[1] += BcAmNO2[1] * NO2[1] + BcCmNO2[1]
    dCa[1] += BcAmCa[1] * Ca[1] + BcCmCa[1]
    dAl[1] += BcAmAl[1] * Al[1] + BcCmAl[1]
    dMo[1] += BcAmMo[1] * Mo[1] + BcCmMo[1]
    dH3PO4[1] += BcAmH3PO4[1] * H3PO4[1] + BcCmH3PO4[1]
    dSO4[1] += BcAmSO4[1] * SO4[1] + BcCmSO4[1]
    dH4SiO4[1] += BcAmH4SiO4[1] * H4SiO4[1] + BcCmH4SiO4[1]
    dPOC[Ngrid] += BcAmPOC[2] * POC[Ngrid] + BcCmPOC[2]
    dMnO2[Ngrid] += BcAmMnO2[2] * MnO2[Ngrid] + BcCmMnO2[2]
    dFeOOH[Ngrid] += BcAmFeOOH[2] * FeOOH[Ngrid] + BcCmFeOOH[2]
    dFeS[Ngrid] += BcAmFeS[2] * FeS[Ngrid] + BcCmFeS[2]
    dFeS2[Ngrid] += BcAmFeS2[2] * FeS2[Ngrid] + BcCmFeS2[2]
    dCaCO3[Ngrid] += BcAmCaCO3[2] * CaCO3[Ngrid] + BcCmCaCO3[2]
    dMnCO3[Ngrid] += BcAmMnCO3[2] * MnCO3[Ngrid] + BcCmMnCO3[2]
    dFeCO3[Ngrid] += BcAmFeCO3[2] * FeCO3[Ngrid] + BcCmFeCO3[2]
    dAge[Ngrid] += BcAmAge[2] * Age[Ngrid] + BcCmAge[2]
    dBSi[Ngrid] += BcAmBSi[2] * BSi[Ngrid] + BcCmBSi[2]
    dMoS4[Ngrid] += BcAmMoS4[2] * MoS4[Ngrid] + BcCmMoS4[2]
    dO2[Ngrid] += BcAmO2[2] * O2[Ngrid] + BcCmO2[2]
    dNO3[Ngrid] += BcAmNO3[2] * NO3[Ngrid] + BcCmNO3[2]
    dCH4[Ngrid] += BcAmCH4[2] * CH4[Ngrid] + BcCmCH4[2]
    dNO2[Ngrid] += BcAmNO2[2] * NO2[Ngrid] + BcCmNO2[2]
    dCa[Ngrid] += BcAmCa[2] * Ca[Ngrid] + BcCmCa[2]
    dAl[Ngrid] += BcAmAl[2] * Al[Ngrid] + BcCmAl[2]
    dMo[Ngrid] += BcAmMo[2] * Mo[Ngrid] + BcCmMo[2]
    dH3PO4[Ngrid] += BcAmH3PO4[2] * H3PO4[Ngrid] + BcCmH3PO4[2]
    dSO4[Ngrid] += BcAmSO4[2] * SO4[Ngrid] + BcCmSO4[2]
    dH4SiO4[Ngrid] += BcAmH4SiO4[2] * H4SiO4[Ngrid] + BcCmH4SiO4[2]

    @. dO2 += alpha * (O2BW - O2)
    @. dNO3 += alpha * (NO3BW - NO3)
    @. dCH4 += alpha * (CH4BW - CH4)
    @. dNO2 += alpha * (NO2BW - NO2)
    @. dCa += alpha * (CaBW - Ca)
    @. dAl += alpha * (AlBW - Al)
    @. dMo += alpha * (MoBW - Mo)
    @. dH3PO4 += alpha * (H3PO4BW - H3PO4)
    @. dSO4 += alpha * (SO4BW - SO4)
    @. dH4SiO4 += alpha * (H4SiO4BW - H4SiO4)

    Mn_ads = @. Mn * KMn_ads
    dMn = AmMn_ads * Mn_ads
    @. dMn *= dstopw
    dMn += AmMn * Mn
    dMn[1] += (BcAmMn_ads[1] * Mn_ads[1] + BcCmMn_ads[1]) * dstopw[1]
    dMn[1] += BcAmMn[1] * Mn[1] + BcCmMn[1]
    dMn[Ngrid] +=
        (BcAmMn_ads[2] * Mn_ads[Ngrid] + BcCmMn_ads[2]) * dstopw[Ngrid]
    dMn[Ngrid] += BcAmMn[2] * Mn[Ngrid] + BcCmMn[2]
    @. dMn += alpha * (Mn0 - Mn)
    @. dMn /= 1 + dstopw * KMn_ads

    Fe_ads = @. Fe * KFe_ads
    dFe = AmFe_ads * Fe_ads
    @. dFe *= dstopw
    dFe += AmFe * Fe
    dFe[1] += (BcAmFe_ads[1] * Fe_ads[1] + BcCmFe_ads[1]) * dstopw[1]
    dFe[1] += BcAmFe[1] * Fe[1] + BcCmFe[1]
    dFe[Ngrid] +=
        (BcAmFe_ads[2] * Fe_ads[Ngrid] + BcCmFe_ads[2]) * dstopw[Ngrid]
    dFe[Ngrid] += BcAmFe[2] * Fe[Ngrid] + BcCmFe[2]
    @. dFe += alpha * (Fe0 - Fe)
    @. dFe /= 1 + dstopw * KFe_ads

    NH4_ads = @. NH4 * KNH4_ads
    dNH4 = AmNH4_ads * NH4_ads
    @. dNH4 *= dstopw
    dNH4 += AmNH4 * NH4
    dNH4[1] += (BcAmNH4_ads[1] * NH4_ads[1] + BcCmNH4_ads[1]) * dstopw[1]
    dNH4[1] += BcAmNH4[1] * NH4[1] + BcCmNH4[1]
    dNH4[Ngrid] +=
        (BcAmNH4_ads[2] * NH4_ads[Ngrid] + BcCmNH4_ads[2]) * dstopw[Ngrid]
    dNH4[Ngrid] += BcAmNH4[2] * NH4[Ngrid] + BcCmNH4[2]
    @. dNH4 += alpha * (NH40 - NH4)
    @. dNH4 /= 1 + dstopw * KNH4_ads

    CO3 = @. KCO2 * KHCO3 * TCO2 / (H^2 + H * KCO2 + KCO2 * KHCO3)
    CO2 = @. H^2 * TCO2 / (H^2 + H * KCO2 + KCO2 * KHCO3)
    HCO3 = @. H * KCO2 * TCO2 / (H^2 + H * KCO2 + KCO2 * KHCO3)
    H2S = @. H * TH2S / (H + KH2S)
    HS = @. KH2S * TH2S / (H + KH2S)
    H4BO4 = @. KH3BO3 * TH3BO3 / (H + KH3BO3)
    H3BO3 = @. H * TH3BO3 / (H + KH3BO3)
    OH = @. KH2O / H

    dTA_dTCO2 = @. KCO2 * (H + 2 * KHCO3) / (H^2 + H * KCO2 + KCO2 * KHCO3)
    dTA_dTH2S = @. KH2S / (H + KH2S)
    dTA_dTH3BO3 = @. KH3BO3 / (H + KH3BO3)

    dTA_dH = @. -KCO2 * TCO2 * (H^2 + 4 * H * KHCO3 + KCO2 * KHCO3) /
       (H^2 + H * KCO2 + KCO2 * KHCO3)^2
    dTA_dH += @. -KH2S * TH2S / (H + KH2S)^2
    dTA_dH += @. -KH3BO3 * TH3BO3 / (H + KH3BO3)^2
    dTA_dH += @. -(H^2 + KH2O) / H^2

    HCO3_tran = AmHCO3 * HCO3
    HCO3_tran[1] += BcAmHCO3[1] * HCO3[1] + BcCmHCO3[1]
    HCO3_tran[Ngrid] += BcAmHCO3[2] * HCO3[Ngrid] + BcCmHCO3[2]
    @. HCO3_tran += alpha * (HCO3BW - HCO3)

    CO3_tran = AmCO3 * CO3
    CO3_tran[1] += BcAmCO3[1] * CO3[1] + BcCmCO3[1]
    CO3_tran[Ngrid] += BcAmCO3[2] * CO3[Ngrid] + BcCmCO3[2]
    @. CO3_tran += alpha * (CO3BW - CO3)

    CO2_tran = AmCO2 * CO2
    CO2_tran[1] += BcAmCO2[1] * CO2[1] + BcCmCO2[1]
    CO2_tran[Ngrid] += BcAmCO2[2] * CO2[Ngrid] + BcCmCO2[2]
    @. CO2_tran += alpha * (CO2BW - CO2)

    H2S_tran = AmH2S * H2S
    H2S_tran[1] += BcAmH2S[1] * H2S[1] + BcCmH2S[1]
    H2S_tran[Ngrid] += BcAmH2S[2] * H2S[Ngrid] + BcCmH2S[2]
    @. H2S_tran += alpha * (H2SBW - H2S)

    HS_tran = AmHS * HS
    HS_tran[1] += BcAmHS[1] * HS[1] + BcCmHS[1]
    HS_tran[Ngrid] += BcAmHS[2] * HS[Ngrid] + BcCmHS[2]
    @. HS_tran += alpha * (HSBW - HS)

    H3BO3_tran = AmH3BO3 * H3BO3
    H3BO3_tran[1] += BcAmH3BO3[1] * H3BO3[1] + BcCmH3BO3[1]
    H3BO3_tran[Ngrid] += BcAmH3BO3[2] * H3BO3[Ngrid] + BcCmH3BO3[2]
    @. H3BO3_tran += alpha * (H3BO3BW - H3BO3)

    H4BO4_tran = AmH4BO4 * H4BO4
    H4BO4_tran[1] += BcAmH4BO4[1] * H4BO4[1] + BcCmH4BO4[1]
    H4BO4_tran[Ngrid] += BcAmH4BO4[2] * H4BO4[Ngrid] + BcCmH4BO4[2]
    @. H4BO4_tran += alpha * (H4BO4BW - H4BO4)

    H_tran = AmH * H
    H_tran[1] += BcAmH[1] * H[1] + BcCmH[1]
    H_tran[Ngrid] += BcAmH[2] * H[Ngrid] + BcCmH[2]
    @. H_tran += alpha * (HBW - H)

    OH_tran = AmOH * OH
    OH_tran[1] += BcAmOH[1] * OH[1] + BcCmOH[1]
    OH_tran[Ngrid] += BcAmOH[2] * OH[Ngrid] + BcCmOH[2]
    @. OH_tran += alpha * (OHBW - OH)

    dTCO2 = @. HCO3_tran + CO3_tran + CO2_tran
    dTH2S = @. H2S_tran + HS_tran
    dTH3BO3 = @. H3BO3_tran + H4BO4_tran

    TA_tran = @. HCO3_tran + 2 * CO3_tran
    TA_tran += @. HS_tran
    TA_tran += @. H4BO4_tran
    TA_tran += @. OH_tran - H_tran

    dH = @. TA_tran
    dH -= @. dTCO2 * dTA_dTCO2
    dH -= @. dTH2S * dTA_dTH2S
    dH -= @. dTH3BO3 * dTA_dTH3BO3
    dH ./= @. dTA_dH
    # speciation
    Fe_free = @. Fe * H / (
        H * (
            794328234.724282 * CO3 * OH +
            4466.83592150963 * CO3 +
            0.758577575029184 * Cl +
            251188.643150958 * HS +
            9.1201083935591 * SO4
        ) +
        H +
        0.00630957344480193 * HS
    )
    Al_free =
        @. Al * H^4 / (H^4 + 9.31537441154163e-19 * H + 2.96802036656283e-24)

    # reaction rates
    RO2POC = @. O2 / (KO2 + O2) * nu / (a + Age) * POC
    RNO2POC = @. NO2 / (KNO2 + NO2) * KO2 / (KO2 + O2) * nu / (a + Age) * POC
    RNO3POC =
        @. NO3 / (KNO3 + NO3) * KNO2 / (KNO2 + NO2) * KO2 / (KO2 + O2) * nu /
           (a + Age) * POC
    RMnO2POC =
        @. MnO2 / (KMnO2 + MnO2) * KNO3 / (KNO3 + NO3) * KNO2 / (KNO2 + NO2) *
           KO2 / (KO2 + O2) * nu / (a + Age) * POC
    RFeOOHPOC = @. FeOOH / (KFeOOH + FeOOH) * KMnO2 / (KMnO2 + MnO2) * KNO3 /
       (KNO3 + NO3) * KNO2 / (KNO2 + NO2) * KO2 / (KO2 + O2) * nu /
       (a + Age) * POC
    RSO4POC = @. SO4 / (KSO4 + SO4) * KFeOOH / (KFeOOH + FeOOH) * KMnO2 /
       (KMnO2 + MnO2) * KNO3 / (KNO3 + NO3) * KNO2 / (KNO2 + NO2) * KO2 /
       (KO2 + O2) * nu / (a + Age) * POC
    RCH4POC = @. KSO4 / (KSO4 + SO4) * KFeOOH / (KFeOOH + FeOOH) * KMnO2 /
       (KMnO2 + MnO2) * KNO3 / (KNO3 + NO3) * KNO2 / (KNO2 + NO2) * KO2 /
       (KO2 + O2) * nu / (a + Age) * POC
    RO2NO2 = @. kO2NO2 * O2 * NO2
    RO2NH4 = @. kO2NH4 * O2 * NH4
    RO2Mn = @. kO2Mn * O2 * Mn
    RO2Mn_ads = @. kO2Mn_ads * O2 * Mn_ads
    RO2Fe = @. kO2Fe * O2 * Fe
    RO2Fe_ads = @. kO2Fe_ads * O2 * Fe_ads
    RO2H2S = @. kO2H2S * O2 * TH2S
    RO2FeS = @. kO2FeS * O2 * FeS
    RO2CH4 = @. kO2CH4 * O2 * CH4
    RNO2NH4 = @. kNO2NH4 * NO2 * NH4
    RNO3Mn = @. kNO3Mn * NO3 * Mn
    RNO3Fe = @. kNO3Fe * NO3 * Fe
    RNO3H2S = @. kNO3H2S * NO3 * TH2S
    RSO4CH4 = @. kAOM * CH4 * SO4 / (SO4 + KAOM)
    RMnO2Fe = @. kMnO2Fe * MnO2 * Fe
    RMnO2H2S = @. kMnO2H2S * MnO2 * TH2S
    RFeOOHH2S = @. kFeOOHH2S * FeOOH * TH2S
    RFeSH2S = @. kFeSH2S * FeS * TH2S
    Omega_RFeS_dis = @. Fe * HS / (H * KspFeS)
    RFeS_dis = @. (Omega_RFeS_dis < 1.0) * kFeSdis * FeS * (1 - Omega_RFeS_dis)
    Omega_RFeS_pre = @. Fe * HS / (H * KspFeS)
    RFeS_pre = @. (Omega_RFeS_pre > 1.0) * kFeSpre * (Omega_RFeS_pre - 1)
    Omega_RCaCO3_dis = @. Ca * CO3 / KspCaCO3_dis
    RCaCO3_dis =
        @. (Omega_RCaCO3_dis < 1.0) * kCaCO3dis * CaCO3 * (1 - Omega_RCaCO3_dis)
    Omega_RCaCO3_pre = @. Ca * CO3 / KspCaCO3_pre
    RCaCO3_pre =
        @. (Omega_RCaCO3_pre > 1.0) * kCaCO3pre * (Omega_RCaCO3_pre - 1)
    Omega_RMnCO3_dis = @. Mn * CO3 / KspMnCO3
    RMnCO3_dis =
        @. (Omega_RMnCO3_dis < 1.0) * kMnCO3dis * MnCO3 * (1 - Omega_RMnCO3_dis)
    Omega_RMnCO3_pre = @. Mn * CO3 / KspMnCO3
    RMnCO3_pre =
        @. (Omega_RMnCO3_pre > 1.0) * kMnCO3pre * (Omega_RMnCO3_pre - 1)
    Omega_RFeCO3_dis = @. Fe * CO3 / KspFeCO3
    RFeCO3_dis =
        @. (Omega_RFeCO3_dis < 1.0) * kFeCO3dis * FeCO3 * (1 - Omega_RFeCO3_dis)
    Omega_RFeCO3_pre = @. Fe * CO3 / KspFeCO3
    RFeCO3_pre =
        @. (Omega_RFeCO3_pre > 1.0) * kFeCO3pre * (Omega_RFeCO3_pre - 1)
    Omega_RBSi_dis = @. H4SiO4 / H4SiO4_dis_sat
    RBSi_dis = @. (Omega_RBSi_dis < 1.0) * kBSidis * (1 - Omega_RBSi_dis) * BSi
    Omega_RASi_pre = @. H4SiO4 / H4SiO4_pre_sat
    RASi_pre = @. (Omega_RASi_pre > 1.0) *
       kASipre *
       H4SiO4_pre_sat *
       (Omega_RASi_pre - 1)
    Omega_RMoS4_pre = @. TH2S / TH2S_Mo_pre
    RMoS4_pre = @. (Omega_RMoS4_pre > 1.0) * kMoS4_pre * Mo * TH2S

    # species rates
    S_POC = @. -1 * RO2POC +
       -1 * RNO2POC +
       -1 * RNO3POC +
       -1 * RMnO2POC +
       -1 * RFeOOHPOC +
       -1 * RSO4POC +
       -1 * RCH4POC
    S_O2 = @. -1 * RO2POC * dstopw +
       -1 / 2 * RO2NO2 +
       -3 / 2 * RO2NH4 +
       -1 / 2 * RO2Mn +
       -1 / 2 * RO2Mn_ads * dstopw +
       -1 / 4 * RO2Fe +
       -1 / 4 * RO2Fe_ads * dstopw +
       -2 * RO2H2S +
       -9 / 4 * RO2FeS * dstopw +
       -2 * RO2CH4
    S_TCO2 = @. 1 * RO2POC * dstopw +
       1 * RNO2POC * dstopw +
       1 * RNO3POC * dstopw +
       1 * RMnO2POC * dstopw +
       1 * RFeOOHPOC * dstopw +
       1 * RSO4POC * dstopw +
       1 / 2 * RCH4POC * dstopw +
       1 * RO2CH4 +
       1 * RSO4CH4 +
       1 * RCaCO3_dis * dstopw +
       -1 * RCaCO3_pre * dstopw +
       1 * RMnCO3_dis * dstopw +
       -1 * RMnCO3_pre * dstopw +
       1 * RFeCO3_dis * dstopw +
       -1 * RFeCO3_pre * dstopw
    S_NH4 = @. rNC * RO2POC / (pwtods + KNH4_ads) +
       rNC * RNO2POC / (pwtods + KNH4_ads) +
       rNC * RNO3POC / (pwtods + KNH4_ads) +
       rNC * RMnO2POC / (pwtods + KNH4_ads) +
       rNC * RFeOOHPOC / (pwtods + KNH4_ads) +
       rNC * RSO4POC / (pwtods + KNH4_ads) +
       rNC * RCH4POC / (pwtods + KNH4_ads) +
       -1 * RO2NH4 / (1 + dstopw * KNH4_ads) +
       -1 * RNO2NH4 / (1 + dstopw * KNH4_ads) +
       1 * RNO3H2S / (1 + dstopw * KNH4_ads)
    S_H3PO4 = @. rPC * RO2POC * dstopw +
       rPC * RNO2POC * dstopw +
       rPC * RNO3POC * dstopw +
       rPC * RMnO2POC * dstopw +
       rPC * RFeOOHPOC * dstopw +
       rPC * RSO4POC * dstopw +
       rPC * RCH4POC * dstopw
    S_NO2 = @. -4 / 3 * RNO2POC * dstopw +
       2 * RNO3POC * dstopw +
       -1 * RO2NO2 +
       1 * RO2NH4 +
       -1 * RNO2NH4
    S_NO3 = @. -2 * RNO3POC * dstopw +
       1 * RO2NO2 +
       -2 / 5 * RNO3Mn +
       -1 / 5 * RNO3Fe +
       -1 * RNO3H2S
    S_MnO2 = @. -2 * RMnO2POC +
       1 * RO2Mn * pwtods +
       1 * RO2Mn_ads +
       1 * RNO3Mn * pwtods +
       -1 / 2 * RMnO2Fe +
       -1 * RMnO2H2S
    S_Mn = @. 2 * RMnO2POC / (pwtods + KMn_ads) +
       -1 * RO2Mn / (1 + dstopw * KMn_ads) +
       -1 * RO2Mn_ads / (pwtods + KMn_ads) +
       -1 * RNO3Mn / (1 + dstopw * KMn_ads) +
       1 / 2 * RMnO2Fe / (pwtods + KMn_ads) +
       1 * RMnO2H2S / (pwtods + KMn_ads) +
       1 * RMnCO3_dis / (pwtods + KMn_ads) +
       -1 * RMnCO3_pre / (pwtods + KMn_ads)
    S_FeOOH = @. -4 * RFeOOHPOC +
       1 * RO2Fe * pwtods +
       1 * RO2Fe_ads +
       1 * RO2FeS +
       1 * RNO3Fe * pwtods +
       1 * RMnO2Fe +
       -2 * RFeOOHH2S
    S_Fe = @. 4 * RFeOOHPOC / (pwtods + KFe_ads) +
       -1 * RO2Fe / (1 + dstopw * KFe_ads) +
       -1 * RO2Fe_ads / (pwtods + KFe_ads) +
       -1 * RNO3Fe / (1 + dstopw * KFe_ads) +
       -1 * RMnO2Fe / (pwtods + KFe_ads) +
       2 * RFeOOHH2S / (pwtods + KFe_ads) +
       1 * RFeS_dis / (pwtods + KFe_ads) +
       -1 * RFeS_pre / (1 + dstopw * KFe_ads) +
       1 * RFeCO3_dis / (pwtods + KFe_ads) +
       -1 * RFeCO3_pre / (pwtods + KFe_ads)
    S_SO4 = @. -1 / 2 * RSO4POC * dstopw +
       1 * RO2H2S +
       1 * RO2FeS * dstopw +
       1 * RNO3H2S +
       -1 * RSO4CH4
    S_TH2S = @. 1 / 2 * RSO4POC * dstopw +
       -1 * RO2H2S +
       -1 * RNO3H2S +
       1 * RSO4CH4 +
       -1 * RMnO2H2S * dstopw +
       -1 * RFeOOHH2S * dstopw +
       -1 * RFeSH2S * dstopw +
       1 * RFeS_dis * dstopw +
       -1 * RFeS_pre +
       -4 * RMoS4_pre
    S_CH4 = @. 1 / 2 * RCH4POC * dstopw + -1 * RO2CH4 + -1 * RSO4CH4
    S_FeS =
        @. -1 * RO2FeS + -1 * RFeSH2S + -1 * RFeS_dis + 1 * RFeS_pre * pwtods
    S_FeS2 = @. 1 * RFeSH2S
    S_CaCO3 = @. -1 * RCaCO3_dis + 1 * RCaCO3_pre
    S_Ca = @. 1 * RCaCO3_dis * dstopw + -1 * RCaCO3_pre * dstopw
    S_MnCO3 = @. -1 * RMnCO3_dis + 1 * RMnCO3_pre
    S_FeCO3 = @. -1 * RFeCO3_dis + 1 * RFeCO3_pre
    S_BSi = @. -1 * RBSi_dis
    S_H4SiO4 = @. 1 * RBSi_dis * dstopw + -1 * RASi_pre
    S_Mo = @. -1 * RMoS4_pre
    S_MoS4 = @. 1 * RMoS4_pre * pwtods
    S_TA = @. 2 * RCaCO3_dis * dstopw +
       -2 * RCaCO3_pre * dstopw +
       2 * RMnCO3_dis * dstopw +
       -2 * RMnCO3_pre * dstopw +
       2 * RFeCO3_dis * dstopw +
       -2 * RFeCO3_pre * dstopw
    S_TA += @. 1 * RFeS_dis * dstopw + -1 * RFeS_pre
    S_TA += @. (rNC - rPC) * RO2POC * dstopw +
       (rNC - rPC + 4 / 3) * RNO2POC * dstopw +
       (rNC - rPC) * RNO3POC * dstopw +
       (rNC - rPC + 4) * RMnO2POC * dstopw +
       (rNC - rPC + 8) * RFeOOHPOC * dstopw +
       (rNC - rPC + 1) * RSO4POC * dstopw +
       (rNC - rPC) * RCH4POC * dstopw +
       -2 * RO2NH4 +
       -2 * RO2Mn +
       -1 * RO2Mn_ads * dstopw +
       -2 * RO2Fe +
       -1 * RO2Fe_ads * dstopw +
       -2 * RO2H2S +
       -2 * RO2FeS * dstopw +
       -8 / 5 * RNO3Mn +
       -9 / 5 * RNO3Fe +
       2 * RSO4CH4 +
       -1 * RMnO2Fe * dstopw +
       2 * RMnO2H2S * dstopw +
       4 * RFeOOHH2S * dstopw +
       1 * RFeS_dis * dstopw +
       -1 * RFeS_pre
    S_H = @. S_TA
    S_H -= @. S_TCO2 * dTA_dTCO2
    S_H -= @. S_TH2S * dTA_dTH2S
    S_H ./= @. dTA_dH
    S_Age = @. ones(Ngrid)

    @. dPOC += S_POC
    @. dMnO2 += S_MnO2
    @. dFeOOH += S_FeOOH
    @. dFeS += S_FeS
    @. dFeS2 += S_FeS2
    @. dCaCO3 += S_CaCO3
    @. dMnCO3 += S_MnCO3
    @. dFeCO3 += S_FeCO3
    @. dAge += S_Age
    @. dBSi += S_BSi
    @. dMoS4 += S_MoS4
    @. dO2 += S_O2
    @. dNO3 += S_NO3
    @. dMn += S_Mn
    @. dFe += S_Fe
    @. dCH4 += S_CH4
    @. dNO2 += S_NO2
    @. dCa += S_Ca
    @. dMo += S_Mo
    @. dNH4 += S_NH4
    @. dH3PO4 += S_H3PO4
    @. dSO4 += S_SO4
    @. dH4SiO4 += S_H4SiO4
    @. dH += S_H
    @. dTCO2 += S_TCO2
    @. dTH2S += S_TH2S

    dPOC += BmPOC * S_POC
    dMnO2 += BmMnO2 * S_MnO2
    dFeOOH += BmFeOOH * S_FeOOH
    dFeS += BmFeS * S_FeS
    dFeS2 += BmFeS2 * S_FeS2
    dCaCO3 += BmCaCO3 * S_CaCO3
    dMnCO3 += BmMnCO3 * S_MnCO3
    dFeCO3 += BmFeCO3 * S_FeCO3
    dAge += BmAge * S_Age
    dBSi += BmBSi * S_BSi
    dMoS4 += BmMoS4 * S_MoS4
    dPOC[1] += BcBmPOC[1] * S_POC[1]
    dMnO2[1] += BcBmMnO2[1] * S_MnO2[1]
    dFeOOH[1] += BcBmFeOOH[1] * S_FeOOH[1]
    dFeS[1] += BcBmFeS[1] * S_FeS[1]
    dFeS2[1] += BcBmFeS2[1] * S_FeS2[1]
    dCaCO3[1] += BcBmCaCO3[1] * S_CaCO3[1]
    dMnCO3[1] += BcBmMnCO3[1] * S_MnCO3[1]
    dFeCO3[1] += BcBmFeCO3[1] * S_FeCO3[1]
    dAge[1] += BcBmAge[1] * S_Age[1]
    dBSi[1] += BcBmBSi[1] * S_BSi[1]
    dMoS4[1] += BcBmMoS4[1] * S_MoS4[1]
    dPOC[Ngrid] += BcBmPOC[2] * S_POC[Ngrid]
    dMnO2[Ngrid] += BcBmMnO2[2] * S_MnO2[Ngrid]
    dFeOOH[Ngrid] += BcBmFeOOH[2] * S_FeOOH[Ngrid]
    dFeS[Ngrid] += BcBmFeS[2] * S_FeS[Ngrid]
    dFeS2[Ngrid] += BcBmFeS2[2] * S_FeS2[Ngrid]
    dCaCO3[Ngrid] += BcBmCaCO3[2] * S_CaCO3[Ngrid]
    dMnCO3[Ngrid] += BcBmMnCO3[2] * S_MnCO3[Ngrid]
    dFeCO3[Ngrid] += BcBmFeCO3[2] * S_FeCO3[Ngrid]
    dAge[Ngrid] += BcBmAge[2] * S_Age[Ngrid]
    dBSi[Ngrid] += BcBmBSi[2] * S_BSi[Ngrid]
    dMoS4[Ngrid] += BcBmMoS4[2] * S_MoS4[Ngrid]

    dC[POCID] .= dPOC
    dC[MnO2ID] .= dMnO2
    dC[FeOOHID] .= dFeOOH
    dC[FeSID] .= dFeS
    dC[FeS2ID] .= dFeS2
    dC[CaCO3ID] .= dCaCO3
    dC[MnCO3ID] .= dMnCO3
    dC[FeCO3ID] .= dFeCO3
    dC[AgeID] .= dAge
    dC[BSiID] .= dBSi
    dC[MoS4ID] .= dMoS4
    dC[O2ID] .= dO2
    dC[NO3ID] .= dNO3
    dC[MnID] .= dMn
    dC[FeID] .= dFe
    dC[CH4ID] .= dCH4
    dC[NO2ID] .= dNO2
    dC[CaID] .= dCa
    dC[AlID] .= dAl
    dC[MoID] .= dMo
    dC[NH4ID] .= dNH4
    dC[H3PO4ID] .= dH3PO4
    dC[SO4ID] .= dSO4
    dC[H4SiO4ID] .= dH4SiO4
    dC[HID] .= dH
    dC[TCO2ID] .= dTCO2
    dC[TH2SID] .= dTH2S
    dC[TH3BO3ID] .= dTH3BO3
    return nothing
end

ODE problem


using DifferentialEquations,ModelingToolkit

prob = ODEProblem(reactran_fvcf_auto, C_uni, (0.0, 1000.0), nothing);

@time sys = ModelingToolkit.modelingtoolkitize(prob)

@time f_MKL = ODEFunction(sys, C_uni, (0.0,1000.0), nothing; jac=true, sparse=true,simplify=true,parallel=ModelingToolkit.MultithreadedForm())

SciML / ModelingToolkit.jl

Errors with ModelingToolkit.MultithreadedForm #1057