search

SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
535 stars 202 forks source link

Example of QNDF failing where CVODE_BDF succeeds #1496

Open ChrisRackauckas opened 2 years ago

ChrisRackauckas commented 2 years ago

This seems to be a problem at the start. I wonder if this needs fixed leading coefficient style.

using DifferentialEquations, Sundials, SpecialFunctions, PreallocationTools,MKL
using LinearAlgebra, SparseArrays, StaticArrays,ForwardDiff
using BenchmarkTools
using FastBroadcast
BLAS.set_num_threads(16)

const ⊕ = +
const ⊗ = *

# 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] =
            (α⁰ * D[1] + β⁰ * u[1]) *
            (-exp(-Pe[1] / 2) / Pe[1] + 1 / Pe[1] - 1 / 2 * exp(-Pe[1] / 2)) /
            (α⁰ * D[1] * exp(-Pe[1] / 2) - α⁰ * D[1] + β⁰ * u[1] * exp(-Pe[1] / 2))
        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

#----------------------------------------------
# Number of substances
#----------------------------------------------
const nsolid = 10 #  # number of solid substances
const ndissolved = 12 #  # number of dissolved substances
const nsummed = 6 #  # number of summed substances
const nnoneq = 22 #  # number of solid + dissolved substances
const nspec = 28 #  # number of total substances
const Lwbdwth = 33 #  # lower bandwidth of jacobian matrix
const Upbdwth = 55 #  # upper bandwidth of jacobian matrix

#----------------------------------------------
# global parameters
#----------------------------------------------
const depth = 3000.0 # m # water depth
const salinity = 35.0 # psu # bottom water salinity
const temp = 2.0 # Celsius # bottom water temperature
const ds_rho = 2.6 # g cm^-3 # dry sediment density
const O2BW = 8.0e-5 # mmol/cm3 # bottom water oxygen
const sw_dens = 1.0290834608199197 # 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 / 20 # 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.8 # dimentionless # surface porosity
const phiL = 0.7 # dimentionless # bottom porosity
const xphi = 100.0 # cm # porosity attenuation scale
const phi_Inf = 0.7 # 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.01014 # 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 = 0.37212703274831377 # cm2/yr # 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 = 2.0 # cm # attentuation scale of bioturbation
const Ds = broadcast(x -> Dbt0 * fO2bt * exp(-x / xbt), x) # cm^2 yr^-1 # Bioturbation coefficient

#----------------------------------------------
# bioirrigation parameters
#----------------------------------------------
const Dbir0 = 28.42212793142007 # 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 * 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 DNH4 = 3.4531311564520092E+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 DNdnr = 1.0763975632520105E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DNdr = 1.0763975632520105E+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 DH3SiO4 = 0.0000000000000000E+00 ./ (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
const DH2PO4 = 1.5332320353118095E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DHPO4 = 1.2376332592731156E+02 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DPO4 = 9.9577971524145752E+01 ./ (1.0 .- 2log.(phif)) .+ 15Ds # cm^2 yr^-1 # Sediment diffusion coefficient
const DH3PO4 = 1.4353131461911693E+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

#----------------------------------------------
# Boundary Condition parameters
#----------------------------------------------
const delta = 0.05 # cm # thickness of the diffusive boundary layer
const FPOC0 = 0.05416666666666667 # mmol cm^-2 yr^-1 # Flux of POC at the  TOP of sediment column
const FMnO20 = 0.00036912995995631593 # mmol cm^-2 yr^-1 # Flux of MnO2 at the  TOP of sediment column
const FFeOOH0 = 0.0005451612903225807 # mmol cm^-2 yr^-1 # Flux of FeOOH at the  TOP of sediment column
const FFeS0 = 0.0 # mmol cm^-2 yr^-1 # Flux of FeS at the  TOP of sediment column
const FFeS20 = 0.0 # mmol cm^-2 yr^-1 # Flux of FeS2 at the  TOP of sediment column
const FCaCO30 = 0.015 # mmol cm^-2 yr^-1 # Flux of CaCO3 at the  TOP of sediment column
const FMnCO30 = 0.0 # mmol cm^-2 yr^-1 # Flux of MnCO3 at the  TOP of sediment column
const FFeCO30 = 0.0 # mmol cm^-2 yr^-1 # Flux of FeCO3 at the  TOP of sediment column
const Age0 = 0.0 # year # missing
const FBSi0 = 0.10679957280170879 # 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 = 0.0 # mmol cm^-2 yr^-1 # Flux of MoS4 at the  TOP of sediment column
const FBasalt0 = Fsed * 0.02 / 87.383 * 1000 # mmol cm^-2 yr^-2 # Flux of Basalt at the  TOP of sediment column
const eNdPO4 = -2.3 # epsilon # missing
const FNdnrPO40 = 1.0e-20 # mmol cm^-2 yr^-2 # missing
const FNdrPO40 = 5.125200932600001e-21 # mmol cm^-2 yr^-2 # missing
const O2BW = 8.0e-5 # mmol cm^-3 # Bottom water concentration of O2
const NO3BW = 2.0e-5 # mmol cm^-3 # Bottom water concentration of NO3
const Mn0 = 5.0e-10 # mmol cm^-3 # Concentration of Mn at the TOP of sediment column
const Fe0 = 1.0e-9 # mmol cm^-3 # Concentration of Fe at the TOP of sediment column
const CH4BW = 0.0 # mmol cm^-3 # Bottom water concentration of CH4
const NO2BW = 2.0e-5 # 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 TH3PO4BW = 2.8e-6 # mmol cm^-3 # Bottom water concentration of H3PO4
const SO4BW = 0.028 # mmol cm^-3 # Bottom water concentration of SO4
const TH4SiO4BW = 0.00017 # mmol cm^-3 # Bottom water concentration of H4SiO4
const NdBW = 3.5e-11 # mmol cm^-3 # Bottom water concentration of Nd
const eNdBW = -2.3 # epsilon # Bottom water eNd
const NdnrBW = 2.314018845499301e-11 # mmol cm^-3 # Bottom water concentration of Nd144
const NdrBW = 1.1859811545006993e-11 # mmol cm^-3 # Bottom water concentration of Nd143
const pHBW = 7.62 # free pH scale # Bottom water pH
const TCO2BW = 0.00236 # mmol cm^-3 # Bottom water concentration of TCO2
const TH2SBW = 0.0 # mmol cm^-3 # Bottom water concentration of TH2S
const TH3BO3BW = 8.7062e-5 # mmol cm^-3 # Bottom water concentration of TH3BO3
const FSurfMn_Ndnr0 = 2.4404962393637638e-8 # missing # missing
const FSurfMn_Ndr0 = 1.2508033601993955e-8 # missing # missing
const FSurfFe_Ndnr0 = 3.604324284694764e-8 # missing # missing
const FSurfFe_Ndr0 = 1.8472886185310433e-8 # missing # missing
const HBW = 2.3988329190194897E-08 # mmol cm^-3 # Bottom water concentration of H
const OHBW = 3.0840219029483606E-07 # mmol cm^-3 # Bottom water concentration of OH
const CO2BW = 5.5211181337131811E-05 # mmol cm^-3 # Bottom water concentration of CO2
const HCO3BW = 2.2598877805390218E-03 # mmol cm^-3 # Bottom water concentration of HCO3
const CO3BW = 4.4901038123846659E-05 # mmol cm^-3 # Bottom water concentration of CO3
const H2SBW = 0.0000000000000000E+00 # mmol cm^-3 # Bottom water concentration of H2S
const HSBW = 0.0000000000000000E+00 # mmol cm^-3 # Bottom water concentration of HS
const H3BO3BW = 8.1450125183506411E-05 # mmol cm^-3 # Bottom water concentration of H3BO3
const H4BO4BW = 5.6118748164935929E-06 # mmol cm^-3 # Bottom water concentration of H4BO4
const H3PO4BW = 7.6017785673264449E-14 # mmol cm^-3 # Bottom water concentration of H3PO4
const H2PO4BW = 8.2097551403566648E-08 # mmol cm^-3 # Bottom water concentration of H2PO4
const HPO4BW = 2.6566241456605459E-06 # mmol cm^-3 # Bottom water concentration of HPO4
const PO4BW = 6.1278226918101791E-08 # mmol cm^-3 # Bottom water concentration of PO4
const H4SiO4BW = 1.6866695527305970E-04 # mmol cm^-3 # Bottom water concentration of H4SiO4
const H3SiO4BW = 1.3330447269403019E-06 # mmol cm^-3 # Bottom water concentration of H3SiO4
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 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 betaSO4 = 3.6069022369165832E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaNdnr = 2.1527951265040210E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaNdr = 2.1527951265040210E+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 betaH3SiO4 = 0.0000000000000000E+00 # 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
const betaH2PO4 = 3.0664640706236187E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaHPO4 = 2.4752665185462311E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaPO4 = 1.9915594304829149E+03 # cm yr^-1 # solute mass transfer velocity across SWI
const betaH3PO4 = 2.8706262923823383E+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 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 = ((1.0, 0.0, Age0), (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 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),
    (0.0, 1.0, 0.0),
) #  # 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), (0.0, 1.0, 0.0)) #  # 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 BcNH4 = ((1.0, 0.0, NH40), (0.0, 1.0, 0.0)) #  # Boundary condition of NH4
const BcSO4 = (
    (betaSO4 + phif[1]uf[1], -phif[1]DSO4[1], betaSO4 * SO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of SO4
const BcNdnr = (
    (betaNdnr + phif[1]uf[1], -phif[1]DNdnr[1], betaNdnr * NdnrBW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of Ndnr
const BcNdr = (
    (betaNdr + phif[1]uf[1], -phif[1]DNdr[1], betaNdr * NdrBW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of Ndr
const BcH4SiO4 = (
    (betaH4SiO4 + phif[1]uf[1], -phif[1]DH4SiO4[1], betaH4SiO4 * H4SiO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH4SiO4
const BcH3SiO4 = (
    (betaH3SiO4 + phif[1]uf[1], -phif[1]DH3SiO4[1], betaH3SiO4 * H3SiO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH4SiO4
const BcCO2 = (
    (betaCO2 + phif[1]uf[1], -phif[1]DCO2[1], betaCO2 * CO2BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TCO2
const BcHCO3 = (
    (betaHCO3 + phif[1]uf[1], -phif[1]DHCO3[1], betaHCO3 * HCO3BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TCO2
const BcCO3 = (
    (betaCO3 + phif[1]uf[1], -phif[1]DCO3[1], betaCO3 * CO3BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TCO2
const BcH2S = (
    (betaH2S + phif[1]uf[1], -phif[1]DH2S[1], betaH2S * H2SBW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH2S
const BcHS =
    ((betaHS + phif[1]uf[1], -phif[1]DHS[1], betaHS * HSBW), (0.0, 1.0, 0.0)) #  # Boundary condition of TH2S
const BcH3BO3 = (
    (betaH3BO3 + phif[1]uf[1], -phif[1]DH3BO3[1], betaH3BO3 * H3BO3BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH3BO3
const BcH4BO4 = (
    (betaH4BO4 + phif[1]uf[1], -phif[1]DH4BO4[1], betaH4BO4 * H4BO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH3BO3
const BcH2PO4 = (
    (betaH2PO4 + phif[1]uf[1], -phif[1]DH2PO4[1], betaH2PO4 * H2PO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH3PO4
const BcHPO4 = (
    (betaHPO4 + phif[1]uf[1], -phif[1]DHPO4[1], betaHPO4 * HPO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH3PO4
const BcPO4 = (
    (betaPO4 + phif[1]uf[1], -phif[1]DPO4[1], betaPO4 * PO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH3PO4
const BcH3PO4 = (
    (betaH3PO4 + phif[1]uf[1], -phif[1]DH3PO4[1], betaH3PO4 * H3PO4BW),
    (0.0, 1.0, 0.0),
) #  # Boundary condition of TH3PO4
const BcH =
    ((betaH + phif[1]uf[1], -phif[1]DH[1], betaH * HBW), (0.0, 1.0, 0.0)) #  # Boundary condition of H
const BcOH =
    ((betaOH + phif[1]uf[1], -phif[1]DOH[1], betaOH * OHBW), (0.0, 1.0, 0.0)) #  # Boundary condition of H
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), (0.0, 1.0, 0.0)) #  # 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 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 BcAmNH4, BcBmNH4, BcCmNH4 = fvcf_bc(phif, DNH4, uf, dx, BcNH4, Ngrid) #  # Boundary transport matrix of NH4
const BcAmSO4, BcBmSO4, BcCmSO4 = fvcf_bc(phif, DSO4, uf, dx, BcSO4, Ngrid) #  # Boundary transport matrix of SO4
const BcAmNdnr, BcBmNdnr, BcCmNdnr = fvcf_bc(phif, DNdnr, uf, dx, BcNdnr, Ngrid) #  # Boundary transport matrix of Ndnr
const BcAmNdr, BcBmNdr, BcCmNdr = fvcf_bc(phif, DNdr, uf, dx, BcNdr, Ngrid) #  # Boundary transport matrix of Ndr
const BcAmH4SiO4, BcBmH4SiO4, BcCmH4SiO4 =
    fvcf_bc(phif, DH4SiO4, uf, dx, BcH4SiO4, Ngrid) #  # Boundary transport matrix of H4SiO4
const BcAmH3SiO4, BcBmH3SiO4, BcCmH3SiO4 =
    fvcf_bc(phif, DH3SiO4, uf, dx, BcH3SiO4, Ngrid) #  # Boundary transport matrix of H3SiO4
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 BcAmH2PO4, BcBmH2PO4, BcCmH2PO4 =
    fvcf_bc(phif, DH2PO4, uf, dx, BcH2PO4, Ngrid) #  # Boundary transport matrix of H2PO4
const BcAmHPO4, BcBmHPO4, BcCmHPO4 = fvcf_bc(phif, DHPO4, uf, dx, BcHPO4, Ngrid) #  # Boundary transport matrix of HPO4
const BcAmPO4, BcBmPO4, BcCmPO4 = fvcf_bc(phif, DPO4, uf, dx, BcPO4, Ngrid) #  # Boundary transport matrix of PO4
const BcAmH3PO4, BcBmH3PO4, BcCmH3PO4 =
    fvcf_bc(phif, DH3PO4, uf, dx, BcH3PO4, Ngrid) #  # Boundary transport matrix of H3PO4
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 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 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 AmNH4, BmNH4 = fvcf(phif, DNH4, uf, dx, Ngrid) #  # Interior transport matrix of NH4
const AmSO4, BmSO4 = fvcf(phif, DSO4, uf, dx, Ngrid) #  # Interior transport matrix of SO4
const AmNdnr, BmNdnr = fvcf(phif, DNdnr, uf, dx, Ngrid) #  # Interior transport matrix of Ndnr
const AmNdr, BmNdr = fvcf(phif, DNdr, uf, dx, Ngrid) #  # Interior transport matrix of Ndr
const AmH4SiO4, BmH4SiO4 = fvcf(phif, DH4SiO4, uf, dx, Ngrid) #  # Interior transport matrix of H4SiO4
const AmH3SiO4, BmH3SiO4 = fvcf(phif, DH3SiO4, uf, dx, Ngrid) #  # Interior transport matrix of H3SiO4
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 AmH2PO4, BmH2PO4 = fvcf(phif, DH2PO4, uf, dx, Ngrid) #  # Interior transport matrix of H2PO4
const AmHPO4, BmHPO4 = fvcf(phif, DHPO4, uf, dx, Ngrid) #  # Interior transport matrix of HPO4
const AmPO4, BmPO4 = fvcf(phif, DPO4, uf, dx, Ngrid) #  # Interior transport matrix of PO4
const AmH3PO4, BmH3PO4 = fvcf(phif, DH3PO4, uf, dx, Ngrid) #  # Interior transport matrix of H3PO4
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 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.3980532637696576E-15 # mmol cm^-3 # Water dissociation constant
const KCO2 = 9.8188321096491150E-07 # mmol cm^-3 # H2CO3 dissociation constant
const KHCO3 = 4.7661697752063504E-10 # mmol cm^-3 # HCO3 dissociation constant
const KH2S = 7.7624695805430387E-07 # mmol cm^-3 # H2S dissociation constant
const KH3BO3 = 1.6527844514531606E-09 # mmol cm^-3 # H3BO3 dissociation constant
const KH3PO4 = 2.5906872600083355E-02 # mmol cm^-3 # H3PO4 dissociation constant
const KH2PO4 = 7.7624695805430387E-07 # mmol cm^-3 # H2PO4 dissociation constant
const KHPO4 = 5.5331962630242337E-10 # mmol cm^-3 # HPO4 dissociation constant
const KH4SiO4 = 1.8958968983182346E-10 # mmol cm^-3 # H4SiO4 dissociation constant

#----------------------------------------------
# Reaction parameters
#----------------------------------------------
const nu = 0.1 # dimentionless # POC reactivity
const a = 0.0001 # yr # initial POC age
const rNC = 0.13675213675213677 # mol mol^-1 # N/C ratio Sediment trap
const rPC = 0.008547008547008548 # mol mol^-2 # P/C ratio Sediment trap
const KO2 = 1.0e-6 #  mmol cm-3 pw yr-1 # Monod constant
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 = 2 / (86.93685 / ds_rho / 10) # mmol cm-3 ds # Monod constant
const KFeOOH = 5 / (88.85174 / ds_rho / 10) # mmol cm-3 ds # Monod constant
const KSO4 = 0.0005 # 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.0e7 # (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 = 100.0 # (mmol cm-3 pw)-0.5 yr-1  #
const kFeSH2S = 300.0 # (mmol cm-3 pw)-1 yr-1  #
const kFeSdis = 31604.559264 # yr-1  #
const KspFeS = 0.0013465481745351092 # (mmol cm-3 pw)^-1  # apparent solubility of FeS
const kFeSpre = 100000.0 # mmol cm-3 ds yr-1  #
const MCaCO3 = 100.09 # g/mol # Calcite molecular weight
const SACaCO3 = 0.35 # m2/g # Calcite specific surface area
const kCaCO3dis0 = 0.055368829236732454 # yr^-1  # close to equilibrium rate
const kCaCO3dis1 = 175.09161176499725 # yr^-1  # far from equilibrum rate
const nCaCO3dis = 5.0 # missing # far from equilibrium reaction order
const KspCaCO3_dis = 7.863574397229617e-7 # (mmol cm^-3 pw)^2 # missing
const kCaCO3pre = 1.0e-8 # mmol cm^-3 ds yr^-1  # missing
const KspCaCO3_pre = 7.863574397229617e-7 # (mmol cm^-3 pw)^2 # missing
const kMnCO3dis = 0.001 # yr-1  # missing
const KspMnCO3 = 4.665192547670257e-9 # (mmol cm-3 pw)^-2  # missing
const kMnCO3pre = 0.000005 * ds_rho # mmol cm-3 ds yr-1  # missing
const kFeCO3dis = 0.3 # yr-1  # missing
const KspFeCO3 = 1.066802767745367e-9 # (mmol cm-3 pw)^-2   # missing
const kFeCO3pre = 0.0001 # mmol cm-3 ds yr-1  # missing
const H4SiO4_dis_sat =
    2.754 * 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 nuBSi = 50e-6 * 365 * 24 # missing # missing
const aBSi = 1.0 # missing # missing
const kASipre = 20e-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 = 4.1 # missing # missing
const rNdrSi = rNdnrSi * (eNd_Basalt / 1e4 + 1) * 0.512638 # missing # missing
const KspBasalt = 4.323907704303024 # missing # missing
const SABasalt = 999.9999999999999 # 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.6473 / 0.0746^0.33 # mol cm^-2 yr^-1 # missing
const kBasalt = kBasalt_0 * SABasalt * Mbasalt # yr^-1 # missing
const aBasalt = 10000.0 # yr  # initial age of Annite
const pBasalt = 0.6 #  exponent of age dependent silicate weathering rate R=R_0*t^p  # missing
const rNdMn = 0.0001 # missing # missing
const eNd_MnO2 = -2.3 # missing # missing
const rNdFe = 0.0001 # missing # missing
const eNd_FeOOH = -2.3 # missing # missing
const KspNdPO4 = 4.6480703029633346e-20 # missing # missing
const kNdPO4_pre = 5.0e-9 # missing # missing
const rFeSi = 0.14120690979293393 # missing # Fe:Si ratio in oceanic arc basalt
const rMnSi = 0.0028532538704253542 # missing # missing
const kChlorite_pre = 1.0e-18 # missing # missing
const KspChlorite = 3.95640731525403e57 # 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]),
    Age0,
    FBSi0 / (phis[1] * us[1]),
    O2BW,
    NO3BW,
    Mn0,
    Fe0,
    CH4BW,
    NO2BW,
    CaBW,
    AlBW,
    NH40,
    SO4BW,
    NdnrBW,
    NdrBW,
    TH4SiO4BW,
    TCO2BW,
    TH2SBW,
    TH3BO3BW,
    TH3PO4BW,
    HBW,
] #  # initial conditions
const C_uni = repeat(C_ini, outer = 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 O2ID = ((1:Ngrid) .- 1)nspec .+ 11 #  # index
const NO3ID = ((1:Ngrid) .- 1)nspec .+ 12 #  # index
const MnID = ((1:Ngrid) .- 1)nspec .+ 13 #  # index
const FeID = ((1:Ngrid) .- 1)nspec .+ 14 #  # index
const CH4ID = ((1:Ngrid) .- 1)nspec .+ 15 #  # index
const NO2ID = ((1:Ngrid) .- 1)nspec .+ 16 #  # index
const CaID = ((1:Ngrid) .- 1)nspec .+ 17 #  # index
const AlID = ((1:Ngrid) .- 1)nspec .+ 18 #  # index
const NH4ID = ((1:Ngrid) .- 1)nspec .+ 19 #  # index
const SO4ID = ((1:Ngrid) .- 1)nspec .+ 20 #  # index
const NdnrID = ((1:Ngrid) .- 1)nspec .+ 21 #  # index
const NdrID = ((1:Ngrid) .- 1)nspec .+ 22 #  # index
const TH4SiO4ID = ((1:Ngrid) .- 1)nspec .+ 23 #  # index
const TCO2ID = ((1:Ngrid) .- 1)nspec .+ 24 #  # index
const TH2SID = ((1:Ngrid) .- 1)nspec .+ 25 #  # index
const TH3BO3ID = ((1:Ngrid) .- 1)nspec .+ 26 #  # index
const TH3PO4ID = ((1:Ngrid) .- 1)nspec .+ 27 #  # index
const HID = ((1:Ngrid) .- 1)nspec .+ 28 #  # index
nothing

mutable struct Reactran{T,chunk_size}
    Mn_ads::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Fe_ads::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    NH4_ads::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    HCO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    CO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    CO2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H3PO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    PO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H2PO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    HPO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    HS::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H4BO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H3BO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H4SiO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H3SiO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    OH::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    HCO3_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    CO3_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    CO2_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H3PO4_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H2PO4_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    HPO4_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    PO4_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H2S_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    HS_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H3BO3_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H4BO4_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H4SiO4_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H3SiO4_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    H_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    OH_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    TA_tran::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    dTA_dH::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    dTA_dTCO2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    dTA_dTH3PO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    dTA_dTH2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    dTA_dTH3BO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    dTA_dTH4SiO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Fe_free::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Al_free::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Ndnr_free::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Ndr_free::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RFeS_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RFeS_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RCaCO3_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RCaCO3_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RMnCO3_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RMnCO3_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RFeCO3_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RFeCO3_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    Omega_RBSi_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2POC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RNO2POC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RNO3POC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RMnO2POC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RFeOOHPOC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RSO4POC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RCH4POC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2NO2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2NH4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2Mn::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2Mn_ads::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2Fe::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2Fe_ads::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2H2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2FeS::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RO2CH4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RNO2NH4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RNO3Mn::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RNO3Fe::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RNO3H2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RSO4CH4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RMnO2Fe::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RMnO2H2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RFeOOHH2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RFeSH2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RFeS_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RFeS_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RCaCO3_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RCaCO3_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RMnCO3_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RMnCO3_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RFeCO3_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RFeCO3_pre::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    RBSi_dis::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_POC::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_O2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_TCO2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_NH4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_TH3PO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_NO2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_NO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_MnO2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_Mn::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_FeOOH::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_Fe::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_SO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_TH2S::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_CH4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_FeS::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_FeS2::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_CaCO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_Ca::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_MnCO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_FeCO3::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_BSi::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_TH4SiO4::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_TA::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_H::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
    S_Age::PreallocationTools.DiffCache{
        Array{T,1},
        Array{ForwardDiff.Dual{nothing,T,chunk_size},1},
    }
end

function init(
    u0::Array{T,1},
    Ngrid::Int,
    ::Val{chunk_size},
) where {T,chunk_size}
    Mn_ads = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Fe_ads = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    NH4_ads = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    HCO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    CO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    CO2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H3PO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    PO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H2PO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    HPO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    HS = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H4BO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H3BO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H4SiO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H3SiO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    OH = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    HCO3_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    CO3_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    CO2_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H3PO4_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H2PO4_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    HPO4_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    PO4_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H2S_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    HS_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H3BO3_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H4BO4_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H4SiO4_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H3SiO4_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    H_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    OH_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    TA_tran = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    dTA_dH = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    dTA_dTCO2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    dTA_dTH3PO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    dTA_dTH2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    dTA_dTH3BO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    dTA_dTH4SiO4 =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Fe_free = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Al_free = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Ndnr_free = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Ndr_free = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RFeS_dis =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RFeS_pre =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RCaCO3_dis =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RCaCO3_pre =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RMnCO3_dis =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RMnCO3_pre =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RFeCO3_dis =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RFeCO3_pre =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    Omega_RBSi_dis =
        PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2POC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RNO2POC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RNO3POC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RMnO2POC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RFeOOHPOC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RSO4POC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RCH4POC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2NO2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2NH4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2Mn = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2Mn_ads = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2Fe = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2Fe_ads = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2H2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2FeS = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RO2CH4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RNO2NH4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RNO3Mn = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RNO3Fe = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RNO3H2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RSO4CH4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RMnO2Fe = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RMnO2H2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RFeOOHH2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RFeSH2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RFeS_dis = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RFeS_pre = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RCaCO3_dis = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RCaCO3_pre = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RMnCO3_dis = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RMnCO3_pre = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RFeCO3_dis = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RFeCO3_pre = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    RBSi_dis = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_POC = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_O2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_TCO2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_NH4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_TH3PO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_NO2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_NO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_MnO2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_Mn = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_FeOOH = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_Fe = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_SO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_TH2S = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_CH4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_FeS = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_FeS2 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_CaCO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_Ca = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_MnCO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_FeCO3 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_BSi = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_TH4SiO4 = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_TA = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_H = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})
    S_Age = PreallocationTools.dualcache(zeros(T, Ngrid), Val{chunk_size})

    cache = Reactran(
        Mn_ads,
        Fe_ads,
        NH4_ads,
        HCO3,
        CO3,
        CO2,
        H3PO4,
        PO4,
        H2PO4,
        HPO4,
        HS,
        H2S,
        H4BO4,
        H3BO3,
        H4SiO4,
        H3SiO4,
        OH,
        HCO3_tran,
        CO3_tran,
        CO2_tran,
        H3PO4_tran,
        H2PO4_tran,
        HPO4_tran,
        PO4_tran,
        H2S_tran,
        HS_tran,
        H3BO3_tran,
        H4BO4_tran,
        H4SiO4_tran,
        H3SiO4_tran,
        H_tran,
        OH_tran,
        TA_tran,
        dTA_dH,
        dTA_dTCO2,
        dTA_dTH3PO4,
        dTA_dTH2S,
        dTA_dTH3BO3,
        dTA_dTH4SiO4,
        Fe_free,
        Al_free,
        Ndnr_free,
        Ndr_free,
        Omega_RFeS_dis,
        Omega_RFeS_pre,
        Omega_RCaCO3_dis,
        Omega_RCaCO3_pre,
        Omega_RMnCO3_dis,
        Omega_RMnCO3_pre,
        Omega_RFeCO3_dis,
        Omega_RFeCO3_pre,
        Omega_RBSi_dis,
        RO2POC,
        RNO2POC,
        RNO3POC,
        RMnO2POC,
        RFeOOHPOC,
        RSO4POC,
        RCH4POC,
        RO2NO2,
        RO2NH4,
        RO2Mn,
        RO2Mn_ads,
        RO2Fe,
        RO2Fe_ads,
        RO2H2S,
        RO2FeS,
        RO2CH4,
        RNO2NH4,
        RNO3Mn,
        RNO3Fe,
        RNO3H2S,
        RSO4CH4,
        RMnO2Fe,
        RMnO2H2S,
        RFeOOHH2S,
        RFeSH2S,
        RFeS_dis,
        RFeS_pre,
        RCaCO3_dis,
        RCaCO3_pre,
        RMnCO3_dis,
        RMnCO3_pre,
        RFeCO3_dis,
        RFeCO3_pre,
        RBSi_dis,
        S_POC,
        S_O2,
        S_TCO2,
        S_NH4,
        S_TH3PO4,
        S_NO2,
        S_NO3,
        S_MnO2,
        S_Mn,
        S_FeOOH,
        S_Fe,
        S_SO4,
        S_TH2S,
        S_CH4,
        S_FeS,
        S_FeS2,
        S_CaCO3,
        S_Ca,
        S_MnCO3,
        S_FeCO3,
        S_BSi,
        S_TH4SiO4,
        S_TA,
        S_H,
        S_Age,
    )
    return cache
end
nothing

function (f::Reactran)(dC, C, parms, t)
    @show t
    Mn_ads = PreallocationTools.get_tmp(f.Mn_ads, C)
    Fe_ads = PreallocationTools.get_tmp(f.Fe_ads, C)
    NH4_ads = PreallocationTools.get_tmp(f.NH4_ads, C)
    HCO3 = PreallocationTools.get_tmp(f.HCO3, C)
    CO3 = PreallocationTools.get_tmp(f.CO3, C)
    CO2 = PreallocationTools.get_tmp(f.CO2, C)
    H3PO4 = PreallocationTools.get_tmp(f.H3PO4, C)
    PO4 = PreallocationTools.get_tmp(f.PO4, C)
    H2PO4 = PreallocationTools.get_tmp(f.H2PO4, C)
    HPO4 = PreallocationTools.get_tmp(f.HPO4, C)
    HS = PreallocationTools.get_tmp(f.HS, C)
    H2S = PreallocationTools.get_tmp(f.H2S, C)
    H4BO4 = PreallocationTools.get_tmp(f.H4BO4, C)
    H3BO3 = PreallocationTools.get_tmp(f.H3BO3, C)
    H4SiO4 = PreallocationTools.get_tmp(f.H4SiO4, C)
    H3SiO4 = PreallocationTools.get_tmp(f.H3SiO4, C)
    OH = PreallocationTools.get_tmp(f.OH, C)
    HCO3_tran = PreallocationTools.get_tmp(f.HCO3_tran, C)
    CO3_tran = PreallocationTools.get_tmp(f.CO3_tran, C)
    CO2_tran = PreallocationTools.get_tmp(f.CO2_tran, C)
    H3PO4_tran = PreallocationTools.get_tmp(f.H3PO4_tran, C)
    H2PO4_tran = PreallocationTools.get_tmp(f.H2PO4_tran, C)
    HPO4_tran = PreallocationTools.get_tmp(f.HPO4_tran, C)
    PO4_tran = PreallocationTools.get_tmp(f.PO4_tran, C)
    H2S_tran = PreallocationTools.get_tmp(f.H2S_tran, C)
    HS_tran = PreallocationTools.get_tmp(f.HS_tran, C)
    H3BO3_tran = PreallocationTools.get_tmp(f.H3BO3_tran, C)
    H4BO4_tran = PreallocationTools.get_tmp(f.H4BO4_tran, C)
    H4SiO4_tran = PreallocationTools.get_tmp(f.H4SiO4_tran, C)
    H3SiO4_tran = PreallocationTools.get_tmp(f.H3SiO4_tran, C)
    H_tran = PreallocationTools.get_tmp(f.H_tran, C)
    OH_tran = PreallocationTools.get_tmp(f.OH_tran, C)
    TA_tran = PreallocationTools.get_tmp(f.TA_tran, C)
    dTA_dH = PreallocationTools.get_tmp(f.dTA_dH, C)
    dTA_dTCO2 = PreallocationTools.get_tmp(f.dTA_dTCO2, C)
    dTA_dTH3PO4 = PreallocationTools.get_tmp(f.dTA_dTH3PO4, C)
    dTA_dTH2S = PreallocationTools.get_tmp(f.dTA_dTH2S, C)
    dTA_dTH3BO3 = PreallocationTools.get_tmp(f.dTA_dTH3BO3, C)
    dTA_dTH4SiO4 = PreallocationTools.get_tmp(f.dTA_dTH4SiO4, C)
    Fe_free = PreallocationTools.get_tmp(f.Fe_free, C)
    Al_free = PreallocationTools.get_tmp(f.Al_free, C)
    Ndnr_free = PreallocationTools.get_tmp(f.Ndnr_free, C)
    Ndr_free = PreallocationTools.get_tmp(f.Ndr_free, C)
    Omega_RFeS_dis = PreallocationTools.get_tmp(f.Omega_RFeS_dis, C)
    Omega_RFeS_pre = PreallocationTools.get_tmp(f.Omega_RFeS_pre, C)
    Omega_RCaCO3_dis = PreallocationTools.get_tmp(f.Omega_RCaCO3_dis, C)
    Omega_RCaCO3_pre = PreallocationTools.get_tmp(f.Omega_RCaCO3_pre, C)
    Omega_RMnCO3_dis = PreallocationTools.get_tmp(f.Omega_RMnCO3_dis, C)
    Omega_RMnCO3_pre = PreallocationTools.get_tmp(f.Omega_RMnCO3_pre, C)
    Omega_RFeCO3_dis = PreallocationTools.get_tmp(f.Omega_RFeCO3_dis, C)
    Omega_RFeCO3_pre = PreallocationTools.get_tmp(f.Omega_RFeCO3_pre, C)
    Omega_RBSi_dis = PreallocationTools.get_tmp(f.Omega_RBSi_dis, C)
    RO2POC = PreallocationTools.get_tmp(f.RO2POC, C)
    RNO2POC = PreallocationTools.get_tmp(f.RNO2POC, C)
    RNO3POC = PreallocationTools.get_tmp(f.RNO3POC, C)
    RMnO2POC = PreallocationTools.get_tmp(f.RMnO2POC, C)
    RFeOOHPOC = PreallocationTools.get_tmp(f.RFeOOHPOC, C)
    RSO4POC = PreallocationTools.get_tmp(f.RSO4POC, C)
    RCH4POC = PreallocationTools.get_tmp(f.RCH4POC, C)
    RO2NO2 = PreallocationTools.get_tmp(f.RO2NO2, C)
    RO2NH4 = PreallocationTools.get_tmp(f.RO2NH4, C)
    RO2Mn = PreallocationTools.get_tmp(f.RO2Mn, C)
    RO2Mn_ads = PreallocationTools.get_tmp(f.RO2Mn_ads, C)
    RO2Fe = PreallocationTools.get_tmp(f.RO2Fe, C)
    RO2Fe_ads = PreallocationTools.get_tmp(f.RO2Fe_ads, C)
    RO2H2S = PreallocationTools.get_tmp(f.RO2H2S, C)
    RO2FeS = PreallocationTools.get_tmp(f.RO2FeS, C)
    RO2CH4 = PreallocationTools.get_tmp(f.RO2CH4, C)
    RNO2NH4 = PreallocationTools.get_tmp(f.RNO2NH4, C)
    RNO3Mn = PreallocationTools.get_tmp(f.RNO3Mn, C)
    RNO3Fe = PreallocationTools.get_tmp(f.RNO3Fe, C)
    RNO3H2S = PreallocationTools.get_tmp(f.RNO3H2S, C)
    RSO4CH4 = PreallocationTools.get_tmp(f.RSO4CH4, C)
    RMnO2Fe = PreallocationTools.get_tmp(f.RMnO2Fe, C)
    RMnO2H2S = PreallocationTools.get_tmp(f.RMnO2H2S, C)
    RFeOOHH2S = PreallocationTools.get_tmp(f.RFeOOHH2S, C)
    RFeSH2S = PreallocationTools.get_tmp(f.RFeSH2S, C)
    RFeS_dis = PreallocationTools.get_tmp(f.RFeS_dis, C)
    RFeS_pre = PreallocationTools.get_tmp(f.RFeS_pre, C)
    RCaCO3_dis = PreallocationTools.get_tmp(f.RCaCO3_dis, C)
    RCaCO3_pre = PreallocationTools.get_tmp(f.RCaCO3_pre, C)
    RMnCO3_dis = PreallocationTools.get_tmp(f.RMnCO3_dis, C)
    RMnCO3_pre = PreallocationTools.get_tmp(f.RMnCO3_pre, C)
    RFeCO3_dis = PreallocationTools.get_tmp(f.RFeCO3_dis, C)
    RFeCO3_pre = PreallocationTools.get_tmp(f.RFeCO3_pre, C)
    RBSi_dis = PreallocationTools.get_tmp(f.RBSi_dis, C)
    S_POC = PreallocationTools.get_tmp(f.S_POC, C)
    S_O2 = PreallocationTools.get_tmp(f.S_O2, C)
    S_TCO2 = PreallocationTools.get_tmp(f.S_TCO2, C)
    S_NH4 = PreallocationTools.get_tmp(f.S_NH4, C)
    S_TH3PO4 = PreallocationTools.get_tmp(f.S_TH3PO4, C)
    S_NO2 = PreallocationTools.get_tmp(f.S_NO2, C)
    S_NO3 = PreallocationTools.get_tmp(f.S_NO3, C)
    S_MnO2 = PreallocationTools.get_tmp(f.S_MnO2, C)
    S_Mn = PreallocationTools.get_tmp(f.S_Mn, C)
    S_FeOOH = PreallocationTools.get_tmp(f.S_FeOOH, C)
    S_Fe = PreallocationTools.get_tmp(f.S_Fe, C)
    S_SO4 = PreallocationTools.get_tmp(f.S_SO4, C)
    S_TH2S = PreallocationTools.get_tmp(f.S_TH2S, C)
    S_CH4 = PreallocationTools.get_tmp(f.S_CH4, C)
    S_FeS = PreallocationTools.get_tmp(f.S_FeS, C)
    S_FeS2 = PreallocationTools.get_tmp(f.S_FeS2, C)
    S_CaCO3 = PreallocationTools.get_tmp(f.S_CaCO3, C)
    S_Ca = PreallocationTools.get_tmp(f.S_Ca, C)
    S_MnCO3 = PreallocationTools.get_tmp(f.S_MnCO3, C)
    S_FeCO3 = PreallocationTools.get_tmp(f.S_FeCO3, C)
    S_BSi = PreallocationTools.get_tmp(f.S_BSi, C)
    S_TH4SiO4 = PreallocationTools.get_tmp(f.S_TH4SiO4, C)
    S_TA = PreallocationTools.get_tmp(f.S_TA, C)
    S_H = PreallocationTools.get_tmp(f.S_H, C)
    S_Age = PreallocationTools.get_tmp(f.S_Age, C)

    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]
    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]
    NH4 = @view C[NH4ID]
    SO4 = @view C[SO4ID]
    Ndnr = @view C[NdnrID]
    Ndr = @view C[NdrID]
    TH4SiO4 = @view C[TH4SiO4ID]
    TCO2 = @view C[TCO2ID]
    TH2S = @view C[TH2SID]
    TH3BO3 = @view C[TH3BO3ID]
    TH3PO4 = @view C[TH3PO4ID]
    H = @view C[HID]

    dPOC = @view dC[POCID]
    dMnO2 = @view dC[MnO2ID]
    dFeOOH = @view dC[FeOOHID]
    dFeS = @view dC[FeSID]
    dFeS2 = @view dC[FeS2ID]
    dCaCO3 = @view dC[CaCO3ID]
    dMnCO3 = @view dC[MnCO3ID]
    dFeCO3 = @view dC[FeCO3ID]
    dAge = @view dC[AgeID]
    dBSi = @view dC[BSiID]
    dO2 = @view dC[O2ID]
    dNO3 = @view dC[NO3ID]
    dMn = @view dC[MnID]
    dFe = @view dC[FeID]
    dCH4 = @view dC[CH4ID]
    dNO2 = @view dC[NO2ID]
    dCa = @view dC[CaID]
    dAl = @view dC[AlID]
    dNH4 = @view dC[NH4ID]
    dSO4 = @view dC[SO4ID]
    dNdnr = @view dC[NdnrID]
    dNdr = @view dC[NdrID]
    dTH4SiO4 = @view dC[TH4SiO4ID]
    dTCO2 = @view dC[TCO2ID]
    dTH2S = @view dC[TH2SID]
    dTH3BO3 = @view dC[TH3BO3ID]
    dTH3PO4 = @view dC[TH3PO4ID]
    dH = @view dC[HID]

    mul!(dPOC, AmPOC, POC)
    mul!(dMnO2, AmMnO2, MnO2)
    mul!(dFeOOH, AmFeOOH, FeOOH)
    mul!(dFeS, AmFeS, FeS)
    mul!(dFeS2, AmFeS2, FeS2)
    mul!(dCaCO3, AmCaCO3, CaCO3)
    mul!(dMnCO3, AmMnCO3, MnCO3)
    mul!(dFeCO3, AmFeCO3, FeCO3)
    mul!(dAge, AmAge, Age)
    mul!(dBSi, AmBSi, BSi)
    mul!(dO2, AmO2, O2)
    mul!(dNO3, AmNO3, NO3)
    mul!(dCH4, AmCH4, CH4)
    mul!(dNO2, AmNO2, NO2)
    mul!(dCa, AmCa, Ca)
    mul!(dAl, AmAl, Al)
    mul!(dSO4, AmSO4, SO4)
    mul!(dNdnr, AmNdnr, Ndnr)
    mul!(dNdr, AmNdr, Ndr)

    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]
    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]
    dSO4[1] += BcAmSO4[1] ⊗ SO4[1] ⊕ BcCmSO4[1]
    dNdnr[1] += BcAmNdnr[1] ⊗ Ndnr[1] ⊕ BcCmNdnr[1]
    dNdr[1] += BcAmNdr[1] ⊗ Ndr[1] ⊕ BcCmNdr[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]
    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]
    dSO4[Ngrid] += BcAmSO4[2] ⊗ SO4[Ngrid] ⊕ BcCmSO4[2]
    dNdnr[Ngrid] += BcAmNdnr[2] ⊗ Ndnr[Ngrid] ⊕ BcCmNdnr[2]
    dNdr[Ngrid] += BcAmNdr[2] ⊗ Ndr[Ngrid] ⊕ BcCmNdr[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)
    @.. dSO4 += alpha ⊗ (SO4BW - SO4)
    @.. dNdnr += alpha ⊗ (NdnrBW - Ndnr)
    @.. dNdr += alpha ⊗ (NdrBW - Ndr)

    @.. Mn_ads = Mn ⊗ KMn_ads
    mul!(dMn, AmMn_ads, Mn_ads)
    @.. dMn *= dstopw
    mul!(dMn, AmMn, Mn, 1.0, 1.0)
    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 / (1 ⊕ dstopw ⊗ KMn_ads)

    @.. Fe_ads = Fe ⊗ KFe_ads
    mul!(dFe, AmFe_ads, Fe_ads)
    @.. dFe *= dstopw
    mul!(dFe, AmFe, Fe, 1.0, 1.0)
    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 / (1 ⊕ dstopw ⊗ KFe_ads)

    @.. NH4_ads = NH4 ⊗ KNH4_ads
    mul!(dNH4, AmNH4_ads, NH4_ads)
    @.. dNH4 *= dstopw
    mul!(dNH4, AmNH4, NH4, 1.0, 1.0)
    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 / (1 ⊕ dstopw ⊗ KNH4_ads)

    @.. HCO3 = H ⊗ KCO2 ⊗ TCO2 / (H^2 ⊕ H ⊗ KCO2 ⊕ KCO2 ⊗ KHCO3)
    @.. CO3 = KCO2 ⊗ KHCO3 ⊗ TCO2 / (H^2 ⊕ H ⊗ KCO2 ⊕ KCO2 ⊗ KHCO3)
    @.. CO2 = H^2 ⊗ TCO2 / (H^2 ⊕ H ⊗ KCO2 ⊕ KCO2 ⊗ KHCO3)
    @.. H3PO4 =
        H^3 ⊗ TH3PO4 /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)
    @.. PO4 =
        KH2PO4 ⊗ KH3PO4 ⊗ KHPO4 ⊗ TH3PO4 /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)
    @.. H2PO4 =
        H^2 ⊗ KH3PO4 ⊗ TH3PO4 /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)
    @.. HPO4 =
        H ⊗ KH2PO4 ⊗ KH3PO4 ⊗ TH3PO4 /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)
    @.. HS = KH2S ⊗ TH2S / (H ⊕ KH2S)
    @.. H2S = H ⊗ TH2S / (H ⊕ KH2S)
    @.. H4BO4 = KH3BO3 ⊗ TH3BO3 / (H ⊕ KH3BO3)
    @.. H3BO3 = H ⊗ TH3BO3 / (H ⊕ KH3BO3)
    @.. H4SiO4 = H ⊗ TH4SiO4 / (H ⊕ KH4SiO4)
    @.. H3SiO4 = KH4SiO4 ⊗ TH4SiO4 / (H ⊕ KH4SiO4)
    @.. OH = KH2O / H

    @.. dTA_dTCO2 = KCO2 ⊗ (H ⊕ 2 ⊗ KHCO3) / (H^2 ⊕ H ⊗ KCO2 ⊕ KCO2 ⊗ KHCO3)
    @.. dTA_dTH3PO4 =
        H ⊗ KH2PO4 ⊗ KH3PO4 /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)
    @.. dTA_dTH3PO4 +=
        2 ⊗ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4 /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)
    @.. dTA_dTH3PO4 +=
        -H^3 /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)
    @.. dTA_dTH2S = KH2S / (H ⊕ KH2S)
    @.. dTA_dTH3BO3 = KH3BO3 / (H ⊕ KH3BO3)
    @.. dTA_dTH4SiO4 = KH4SiO4 / (H ⊕ KH4SiO4)

    @.. dTA_dH =
        -KCO2 ⊗ TCO2 ⊗ (H^2 ⊕ 4 ⊗ H ⊗ KHCO3 ⊕ KCO2 ⊗ KHCO3) /
        (H^2 ⊕ H ⊗ KCO2 ⊕ KCO2 ⊗ KHCO3)^2
    @.. dTA_dH +=
        -KH2PO4 ⊗ KH3PO4 ⊗ TH3PO4 ⊗
        (2 ⊗ H^3 ⊕ H^2 ⊗ KH3PO4 - KH2PO4 ⊗ KH3PO4 ⊗ KHPO4) /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)^2
    @.. dTA_dH +=
        -2 ⊗ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4 ⊗ TH3PO4 ⊗
        (3 ⊗ H^2 ⊕ 2 ⊗ H ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4) /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)^2
    @.. dTA_dH +=
        -H^2 ⊗ KH3PO4 ⊗ TH3PO4 ⊗ (H^2 ⊕ 2 ⊗ H ⊗ KH2PO4 ⊕ 3 ⊗ KH2PO4 ⊗ KHPO4) /
        (H^3 ⊕ H^2 ⊗ KH3PO4 ⊕ H ⊗ KH2PO4 ⊗ KH3PO4 ⊕ KH2PO4 ⊗ KH3PO4 ⊗ KHPO4)^2
    @.. dTA_dH += -KH2S ⊗ TH2S / (H ⊕ KH2S)^2
    @.. dTA_dH += -KH3BO3 ⊗ TH3BO3 / (H ⊕ KH3BO3)^2
    @.. dTA_dH += -KH4SiO4 ⊗ TH4SiO4 / (H ⊕ KH4SiO4)^2
    @.. dTA_dH += -(H^2 ⊕ KH2O) / H^2

    mul!(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)

    mul!(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)

    mul!(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)

    mul!(H3PO4_tran, AmH3PO4, H3PO4)
    H3PO4_tran[1] += BcAmH3PO4[1] ⊗ H3PO4[1] ⊕ BcCmH3PO4[1]
    H3PO4_tran[Ngrid] += BcAmH3PO4[2] ⊗ H3PO4[Ngrid] ⊕ BcCmH3PO4[2]
    @.. H3PO4_tran += alpha ⊗ (H3PO4BW - H3PO4)

    mul!(H2PO4_tran, AmH2PO4, H2PO4)
    H2PO4_tran[1] += BcAmH2PO4[1] ⊗ H2PO4[1] ⊕ BcCmH2PO4[1]
    H2PO4_tran[Ngrid] += BcAmH2PO4[2] ⊗ H2PO4[Ngrid] ⊕ BcCmH2PO4[2]
    @.. H2PO4_tran += alpha ⊗ (H2PO4BW - H2PO4)

    mul!(HPO4_tran, AmHPO4, HPO4)
    HPO4_tran[1] += BcAmHPO4[1] ⊗ HPO4[1] ⊕ BcCmHPO4[1]
    HPO4_tran[Ngrid] += BcAmHPO4[2] ⊗ HPO4[Ngrid] ⊕ BcCmHPO4[2]
    @.. HPO4_tran += alpha ⊗ (HPO4BW - HPO4)

    mul!(PO4_tran, AmPO4, PO4)
    PO4_tran[1] += BcAmPO4[1] ⊗ PO4[1] ⊕ BcCmPO4[1]
    PO4_tran[Ngrid] += BcAmPO4[2] ⊗ PO4[Ngrid] ⊕ BcCmPO4[2]
    @.. PO4_tran += alpha ⊗ (PO4BW - PO4)

    mul!(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)

    mul!(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)

    mul!(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)

    mul!(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)

    mul!(H4SiO4_tran, AmH4SiO4, H4SiO4)
    H4SiO4_tran[1] += BcAmH4SiO4[1] ⊗ H4SiO4[1] ⊕ BcCmH4SiO4[1]
    H4SiO4_tran[Ngrid] += BcAmH4SiO4[2] ⊗ H4SiO4[Ngrid] ⊕ BcCmH4SiO4[2]
    @.. H4SiO4_tran += alpha ⊗ (H4SiO4BW - H4SiO4)

    mul!(H3SiO4_tran, AmH3SiO4, H3SiO4)
    H3SiO4_tran[1] += BcAmH3SiO4[1] ⊗ H3SiO4[1] ⊕ BcCmH3SiO4[1]
    H3SiO4_tran[Ngrid] += BcAmH3SiO4[2] ⊗ H3SiO4[Ngrid] ⊕ BcCmH3SiO4[2]
    @.. H3SiO4_tran += alpha ⊗ (H3SiO4BW - H3SiO4)

    mul!(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)

    mul!(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
    @.. dTH3PO4 = H3PO4_tran ⊕ H2PO4_tran ⊕ HPO4_tran ⊕ PO4_tran
    @.. dTH2S = H2S_tran ⊕ HS_tran
    @.. dTH3BO3 = H3BO3_tran ⊕ H4BO4_tran
    @.. dTH4SiO4 = H4SiO4_tran ⊕ H3SiO4_tran

    @.. TA_tran = HCO3_tran ⊕ 2 ⊗ CO3_tran
    @.. TA_tran += HPO4_tran ⊕ 2 ⊗ PO4_tran - H3PO4_tran
    @.. TA_tran += HS_tran
    @.. TA_tran += H4BO4_tran
    @.. TA_tran += H3SiO4_tran
    @.. TA_tran += OH_tran - H_tran

    @.. dH = TA_tran
    @.. dH -= dTCO2 ⊗ dTA_dTCO2
    @.. dH -= dTH3PO4 ⊗ dTA_dTH3PO4
    @.. dH -= dTH2S ⊗ dTA_dTH2S
    @.. dH -= dTH3BO3 ⊗ dTA_dTH3BO3
    @.. dH -= dTH4SiO4 ⊗ dTA_dTH4SiO4
    @.. dH = dH / dTA_dH
    # speciation
    @.. Fe_free =
        Fe ⊗ H / (
            H ⊗ (
                38404.8078076036 ⊗ CO3^2 ⊕ 1141.13847737758 ⊗ CO3 ⊕
                0.541149470776101 ⊗ Cl ⊕ 3.75214432319968 ⊗ HCO3 ⊕
                409766.824665317 ⊗ HS ⊕ 390245.286156751 ⊗ OH^2 ⊕
                3971.47624308745 ⊗ OH ⊕ 6.626643393432 ⊗ SO4
            ) ⊕ H ⊕ 0.00150222599838656 ⊗ HS
        )
    @.. Al_free =
        Al / (
            1889348.44155961 ⊗ CO3 ⊕ 5.71374695230749e+39 ⊗ OH^4 ⊕
            1.2556441493638e+31 ⊗ OH^3 ⊕ 5.63992735079517e+20 ⊗ OH^2 ⊕
            7872134080.00216 ⊗ OH ⊕ 1
        )
    @.. Ndnr_free =
        H ⊗ Ndnr / (
            H ⊗ (
                3632465961.77468 ⊗ CO3^2 ⊕ 290431.690955537 ⊗ CO3 ⊕
                0.343269152592752 ⊗ Cl ⊕ 63095734448.0194 ⊗ H3SiO4^2 ⊕
                158489.319246111 ⊗ H3SiO4 ⊕ 7.38572662109916 ⊗ HCO3 ⊕
                56.0142660345257 ⊗ SO4
            ) ⊕ H ⊕ 9.62959782270152e-10
        )
    @.. Ndr_free =
        H ⊗ Ndr / (
            H ⊗ (
                3632465961.77468 ⊗ CO3^2 ⊕ 290431.690955537 ⊗ CO3 ⊕
                0.343269152592752 ⊗ Cl ⊕ 63095734448.0194 ⊗ H3SiO4^2 ⊕
                158489.319246111 ⊗ H3SiO4 ⊕ 7.38572662109917 ⊗ HCO3 ⊕
                56.0142660345257 ⊗ SO4
            ) ⊕ H ⊕ 9.62959782270152e-10
        )

    # reaction rates
    @.. Omega_RFeS_dis = Fe_free ⊗ HS / (H ⊗ KspFeS)
    @.. Omega_RFeS_pre = Fe_free ⊗ HS / (H ⊗ KspFeS)
    @.. Omega_RCaCO3_dis = Ca ⊗ CO3 / KspCaCO3_dis
    @.. Omega_RCaCO3_pre = Ca ⊗ CO3 / KspCaCO3_pre
    @.. Omega_RMnCO3_dis = Mn ⊗ CO3 / KspMnCO3
    @.. Omega_RMnCO3_pre = Mn ⊗ CO3 / KspMnCO3
    @.. Omega_RFeCO3_dis = Fe_free ⊗ CO3 / KspFeCO3
    @.. Omega_RFeCO3_pre = Fe_free ⊗ CO3 / KspFeCO3
    @.. Omega_RBSi_dis = H4SiO4 / H4SiO4_dis_sat
    @.. 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
    @.. RFeS_dis =
        (-tanh(100.0 ⊗ (Omega_RFeS_dis - 1.0)) / 2 ⊕ 0.5) ⊗
        (kFeSdis ⊗ FeS ⊗ (1 - Omega_RFeS_dis))
    @.. RFeS_pre =
        (tanh(100.0 ⊗ (Omega_RFeS_pre - 1.0)) / 2 ⊕ 0.5) ⊗
        (kFeSpre ⊗ Fe ⊗ TH2S ⊗ (Omega_RFeS_pre - 1))
    @.. RCaCO3_dis =
        (-tanh(100.0 ⊗ (Omega_RCaCO3_dis - 1.0)) / 2 ⊕ 0.5) ⊗ (
            kCaCO3dis0 ⊗ CaCO3 ⊕
            kCaCO3dis1 ⊗ CaCO3 ⊗ (1 - Omega_RCaCO3_dis)^nCaCO3dis
        )
    @.. RCaCO3_pre =
        (tanh(100.0 ⊗ (Omega_RCaCO3_pre - 1.0)) / 2 ⊕ 0.5) ⊗
        (kCaCO3pre ⊗ CaCO3 ⊗ (Omega_RCaCO3_pre - 1))
    @.. RMnCO3_dis =
        (-tanh(100.0 ⊗ (Omega_RMnCO3_dis - 1.0)) / 2 ⊕ 0.5) ⊗
        (kMnCO3dis ⊗ MnCO3 ⊗ (1 - Omega_RMnCO3_dis))
    @.. RMnCO3_pre =
        (tanh(100.0 ⊗ (Omega_RMnCO3_pre - 1.0)) / 2 ⊕ 0.5) ⊗
        (kMnCO3pre ⊗ (Omega_RMnCO3_pre - 1))
    @.. RFeCO3_dis =
        (-tanh(100.0 ⊗ (Omega_RFeCO3_dis - 1.0)) / 2 ⊕ 0.5) ⊗
        (kFeCO3dis ⊗ FeCO3 ⊗ (1 - Omega_RFeCO3_dis))
    @.. RFeCO3_pre =
        (tanh(100.0 ⊗ (Omega_RFeCO3_pre - 1.0)) / 2 ⊕ 0.5) ⊗
        (kFeCO3pre ⊗ (Omega_RFeCO3_pre - 1))
    @.. RBSi_dis =
        (-tanh(100.0 ⊗ (Omega_RBSi_dis - 1.0)) / 2 ⊕ 0.5) ⊗
        ((1 - Omega_RBSi_dis) ⊗ BSi ⊗ nuBSi / (aBSi ⊕ Age))

    # 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_TH3PO4 =
        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
    @.. 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_TH4SiO4 = 1 ⊗ RBSi_dis ⊗ dstopw
    @.. 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_TH3PO4 ⊗ dTA_dTH3PO4
    @.. S_H -= S_TH2S ⊗ dTA_dTH2S
    @.. S_H -= S_TH4SiO4 ⊗ dTA_dTH4SiO4
    @.. S_H = S_H / dTA_dH
    @.. S_Age = 1

    @.. 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
    @.. dO2 += S_O2
    @.. dNO3 += S_NO3
    @.. dMn += S_Mn
    @.. dFe += S_Fe
    @.. dCH4 += S_CH4
    @.. dNO2 += S_NO2
    @.. dCa += S_Ca
    @.. dNH4 += S_NH4
    @.. dSO4 += S_SO4
    @.. dTH4SiO4 += S_TH4SiO4
    @.. dTCO2 += S_TCO2
    @.. dTH2S += S_TH2S
    @.. dTH3PO4 += S_TH3PO4
    @.. dH += S_H

    mul!(dPOC, BmPOC, S_POC, 1.0, 1.0)
    mul!(dMnO2, BmMnO2, S_MnO2, 1.0, 1.0)
    mul!(dFeOOH, BmFeOOH, S_FeOOH, 1.0, 1.0)
    mul!(dFeS, BmFeS, S_FeS, 1.0, 1.0)
    mul!(dFeS2, BmFeS2, S_FeS2, 1.0, 1.0)
    mul!(dCaCO3, BmCaCO3, S_CaCO3, 1.0, 1.0)
    mul!(dMnCO3, BmMnCO3, S_MnCO3, 1.0, 1.0)
    mul!(dFeCO3, BmFeCO3, S_FeCO3, 1.0, 1.0)
    mul!(dAge, BmAge, S_Age, 1.0, 1.0)
    mul!(dBSi, BmBSi, S_BSi, 1.0, 1.0)
    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]
    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]

    return nothing
end

chunk_size = 10;
reactran_fvcf = init(zeros(Ngrid*nspec),Ngrid,Val(chunk_size));

tspan = (0.0, 15000.0);
prob = ODEProblem(reactran_fvcf, C_uni, tspan, nothing);
cb = TerminateSteadyState(1e-8, 1e-6, DiffEqCallbacks.allDerivPass);

@time sol = DifferentialEquations.solve(
    prob,
    CVODE_BDF(linear_solver = :LapackBand, jac_upper = Upbdwth, jac_lower = Lwbdwth),
    reltol = 1e-6,
    abstol = 1e-15,
    save_everystep = false,
    callback = cb,
    saveat = 50.0,
    tstops = 5.0,
    # dtmax = 0.01,
    maxiters = Int(1e6),
)

using BandedMatrices
jac_prototype = BandedMatrix(Ones(Ngrid*nspec,Ngrid*nspec),(Lwbdwth,Upbdwth));
reactran_band = ODEFunction(reactran_fvcf,jac_prototype=jac_prototype);

qr!(jac_prototype)

tspan = (0.0, 15.0);
prob_band = ODEProblem(reactran_band, sol[end], tspan, nothing);
cb = TerminateSteadyState(1e-8, 1e-6, DiffEqCallbacks.allDerivPass);

@time sol = DifferentialEquations.solve(
    prob_band,
    FBDF(autodiff=false),
#     QNDF(autodiff=true,chunk_size=chunk_size), using AD
    reltol = 1e-6,
    abstol = 1e-15,
    save_everystep = false,
    callback = cb,
    saveat = 50.0,
    tstops = 5.0,
    # dtmax = 0.01,
    maxiters = Int(1e6),
)
ChrisRackauckas commented 2 years ago

@JianghuiDu try the improved FBDF on the newest version with this and see if that works.

JianghuiDu commented 2 years ago

Where do I get the latest? Is it already registered?

ChrisRackauckas commented 2 years ago

It's registered.

JianghuiDu commented 2 years ago

Still doesn't work... Again, CVODE+LapackBand takes 3 mins to get to t=13000.

FBDF + BandedMatrices It seems to be working. But very slow. I waited for half an hour and still not reaching t=1.