Example of QNDF failing where CVODE_BDF succeeds #1496
Description
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),
)