Merge pull request #201 from ciadht/hmw-ct-jsw/PISCES

Hmw ct jsw/pisces

src/Models/AdvectedPopulations/PISCES/DOC.jl

@@ -69,9 +69,9 @@ end
     K_Feᵇᵃᶜᵗ = bgc.Fe_half_saturation_const_for_DOC_remin
     K_DOC = bgc.half_saturation_const_for_DOC_remin
 
-    L_DOCᵇᵃᶜᵗ = K_mondo(DOC, K_DOC) #34b
-    L_Feᵇᵃᶜᵗ = K_mondo(bFe, K_Feᵇᵃᶜᵗ) #34d
-    Lₚₒ₄ᵇᵃᶜᵗ = K_mondo(PO₄, Kₚₒ₄ᵇᵃᶜᵗ) #34e
+    L_DOCᵇᵃᶜᵗ = concentration_limitation(DOC, K_DOC) #34b
+    L_Feᵇᵃᶜᵗ = concentration_limitation(bFe, K_Feᵇᵃᶜᵗ) #34d
+    Lₚₒ₄ᵇᵃᶜᵗ = concentration_limitation(PO₄, Kₚₒ₄ᵇᵃᶜᵗ) #34e
 
     Lₙₕ₄ᵇᵃᶜᵗ = L_NH₄(NO₃, NH₄, Kₙₒ₃ᵇᵃᶜᵗ, Kₙₕ₄ᵇᵃᶜᵗ) #34g
     Lₙₒ₃ᵇᵃᶜᵗ = L_NO₃(NO₃, NH₄, Kₙₒ₃ᵇᵃᶜᵗ, Kₙₕ₄ᵇᵃᶜᵗ) #34h
@@ -141,5 +141,6 @@ end
 
     Φ₁ᴰᴼᶜ, Φ₂ᴰᴼᶜ, Φ₃ᴰᴼᶜ = Φᴰᴼᶜ(DOC, POC, GOC, sh, bgc)
 
-    return (1 - γᶻ)*(1 - eᶻ - σᶻ)*∑ᵢgᵢᶻ*Z + (1 - γᴹ)*(1 - eᴹ - σᴹ)*(∑ᵢgᵢᴹ + g_GOC_FFᴹ + gₚₒ_FFᴹ)*M + δᴰ*μᴰ*D + δᴾ*μᴾ*P + λₚₒ¹*POC + (1 - γᴹ)*Rᵤₚᴹ - Remin - Denit - Φ₁ᴰᴼᶜ - Φ₂ᴰᴼᶜ - Φ₃ᴰᴼᶜ #32
+    return ((1 - γᶻ)*(1 - eᶻ - σᶻ)*∑ᵢgᵢᶻ*Z + (1 - γᴹ)*(1 - eᴹ - σᴹ)*(∑ᵢgᵢᴹ + g_GOC_FFᴹ + gₚₒ_FFᴹ)*M + 
+           δᴰ*μᴰ*D + δᴾ*μᴾ*P + λₚₒ¹*POC + (1 - γᴹ)*Rᵤₚᴹ - Remin - Denit - Φ₁ᴰᴼᶜ - Φ₂ᴰᴼᶜ - Φ₃ᴰᴼᶜ) #32
 end #changed this to include gₚₒ_FF

src/Models/AdvectedPopulations/PISCES/PISCES.jl

@@ -176,7 +176,7 @@ struct PISCES{FT, PD, ZM, OT, W, CF, ZF} <: AbstractContinuousFormBiogeochemistr
     coefficient_of_bacterial_uptake_of_iron_in_POC :: FT
     coefficient_of_bacterial_uptake_of_iron_in_GOC :: FT
     max_FeC_ratio_of_bacteria :: FT
-    Fe_half_saturation_const_for_PLACEHOLDER :: FT    #not sure what this should be called
+    Fe_half_saturation_const_for_Bacteria :: FT    #not sure what this should be called
     proportion_of_sinking_grazed_shells :: ZM
 
     mixed_layer_depth :: CF
@@ -298,7 +298,7 @@ struct PISCES{FT, PD, ZM, OT, W, CF, ZF} <: AbstractContinuousFormBiogeochemistr
                     coefficient_of_bacterial_uptake_of_iron_in_POC :: FT,
                     coefficient_of_bacterial_uptake_of_iron_in_GOC :: FT,
                     max_FeC_ratio_of_bacteria :: FT,
-                    Fe_half_saturation_const_for_PLACEHOLDER :: FT,    #not sure what this should be called
+                    Fe_half_saturation_const_for_Bacteria :: FT,    #not sure what this should be called
                     proportion_of_sinking_grazed_shells :: ZM,
 
                     mixed_layer_depth :: CF,
@@ -420,7 +420,7 @@ struct PISCES{FT, PD, ZM, OT, W, CF, ZF} <: AbstractContinuousFormBiogeochemistr
                             coefficient_of_bacterial_uptake_of_iron_in_POC,
                             coefficient_of_bacterial_uptake_of_iron_in_GOC,
                             max_FeC_ratio_of_bacteria,
-                            Fe_half_saturation_const_for_PLACEHOLDER,    #not sure what this should be called
+                            Fe_half_saturation_const_for_Bacteria,    #not sure what this should be called
                             proportion_of_sinking_grazed_shells,
 
                             mixed_layer_depth,
@@ -547,7 +547,7 @@ end
                    coefficient_of_bacterial_uptake_of_iron_in_POC :: FT = 0.5,
                    coefficient_of_bacterial_uptake_of_iron_in_GOC :: FT = 0.5,
                    max_FeC_ratio_of_bacteria :: FT = 10.0e-6,     #or 6
-                   Fe_half_saturation_const_for_PLACEHOLDER :: FT = 2.5e-10, #or 2.5e-10    #not sure what this should be called
+                   Fe_half_saturation_const_for_Bacteria :: FT = 2.5e-10, #or 2.5e-10    #not sure what this should be called
                    proportion_of_sinking_grazed_shells :: ZM = (Z = 0.3, M = 0.3),  # 0.3 for both? not sure
 
                    mixed_layer_depth :: CF = ConstantField(100),
@@ -716,14 +716,14 @@ function PISCES(; grid, # finally the function
                    coefficient_of_bacterial_uptake_of_iron_in_POC :: FT = 0.5,
                    coefficient_of_bacterial_uptake_of_iron_in_GOC :: FT = 0.5,
                    max_FeC_ratio_of_bacteria :: FT = 10.0e-3,     #or 6
-                   Fe_half_saturation_const_for_PLACEHOLDER :: FT = 2.5e-10, #or 2.5e-10    #not sure what this should be called
+                   Fe_half_saturation_const_for_Bacteria :: FT = 0.03, #or 2.5e-10    #not sure what this should be called
                    proportion_of_sinking_grazed_shells :: ZM = (Z = 0.3, M = 0.3),  # 0.3 for both? not sure
 
                    mixed_layer_depth :: CF = ConstantField(100),
                    euphotic_layer_depth :: CF = ConstantField(50),
                    vertical_diffusivity :: CF  = ConstantField(1),
                    yearly_maximum_silicate :: FT = 1.0,
-                   dust_deposition :: FT = 1.0,
+                   dust_deposition :: FT = 0.0,
 
                   surface_photosynthetically_active_radiation = default_surface_PAR,
 
@@ -860,7 +860,7 @@ function PISCES(; grid, # finally the function
                                         coefficient_of_bacterial_uptake_of_iron_in_POC,
                                         coefficient_of_bacterial_uptake_of_iron_in_GOC,
                                         max_FeC_ratio_of_bacteria,
-                                        Fe_half_saturation_const_for_PLACEHOLDER,    #not sure what this should be called
+                                        Fe_half_saturation_const_for_Bacteria,    #not sure what this should be called
                                         proportion_of_sinking_grazed_shells,
 
                                         mixed_layer_depth,

src/Models/AdvectedPopulations/PISCES/POC_and_GOC.jl

@@ -37,7 +37,12 @@ end
     gₚₒ_FFᴹ = g_FF*(bₘ^T)*wₚₒ*POC #29a
     Φ = get_Φ(POC, GOC, sh, bgc)
 
-    return σᶻ*∑gᶻ*Z + 0.5*mᴰ*K_mondo(D, Kₘ)*D + rᶻ*(b_Z^T)*(K_mondo(Z, Kₘ) + 3*ΔO₂(O₂, bgc))*Z + mᶻ*(b_Z^T)*Z^2 + (1 - 0.5*R_CaCO₃)*(mᴾ*K_mondo(P, Kₘ)*P + sh*wᴾ*P^2) + λₚₒ¹*GOC + Φ₁ᴰᴼᶜ + Φ₃ᴰᴼᶜ - (gₚₒᴹ + gₚₒ_FFᴹ)*M - gₚₒᶻ*Z - λₚₒ¹*POC - Φ  #37
+    #println("Total change in POC is ", (σᶻ*∑gᶻ*Z + 0.5*mᴰ*concentration_limitation(D, Kₘ)*D + rᶻ*(b_Z^T)*(concentration_limitation(Z, Kₘ) + 3*ΔO₂(O₂, bgc))*Z + mᶻ*(b_Z^T)*Z^2 + (1 - 0.5*R_CaCO₃)*(mᴾ*concentration_limitation(P, Kₘ)*P + sh*wᴾ*P^2) + λₚₒ¹*GOC + Φ₁ᴰᴼᶜ + Φ₃ᴰᴼᶜ - (gₚₒᴹ + gₚₒ_FFᴹ)*M - gₚₒᶻ*Z - λₚₒ¹*POC - Φ))
+    #println("Positive terms: a = $(σᶻ*∑gᶻ*Z), b = $(0.5*mᴰ*concentration_limitation(D, Kₘ)*D), c = $(rᶻ*(b_Z^T)*(concentration_limitation(Z, Kₘ) + 3*ΔO₂(O₂, bgc))*Z), " )
+
+    return (σᶻ*∑gᶻ*Z + 0.5*mᴰ*concentration_limitation(D, Kₘ)*D + rᶻ*(b_Z^T)*(concentration_limitation(Z, Kₘ) + 3*ΔO₂(O₂, bgc))*Z 
+            + mᶻ*(b_Z^T)*Z^2 + (1 - 0.5*R_CaCO₃)*(mᴾ*concentration_limitation(P, Kₘ)*P + sh*wᴾ*P^2) 
+            + λₚₒ¹*GOC + Φ₁ᴰᴼᶜ + Φ₃ᴰᴼᶜ - (gₚₒᴹ + gₚₒ_FFᴹ)*M - gₚₒᶻ*Z - λₚₒ¹*POC - Φ)  #37
 end
 
 @inline function (bgc::PISCES)(::Val{:GOC}, x, y, z, t, P, D, Z, M, Pᶜʰˡ, Dᶜʰˡ, Pᶠᵉ, Dᶠᵉ, Dˢⁱ, DOC, POC, GOC, SFe, BFe, PSi, NO₃, NH₄, PO₄, Fe, Si, CaCO₃, DIC, Alk, O₂, T, zₘₓₗ, zₑᵤ, Si̅, D_dust, Ω, PAR, PAR¹, PAR², PAR³)
@@ -69,5 +74,7 @@ end
     g_GOC_FFᴹ = g_FF*bₘ^T*w_GOC*GOC #29b
     λₚₒ¹ = λ¹(T, O₂, bgc)
 
-    return σᴹ*(∑gᴹ + ∑g_FFᴹ)*M + rᴹ*bₘ^T*(K_mondo(M, Kₘ) + 3*ΔO₂(O₂, bgc))*M + Pᵤₚᴹ + 0.5*R_CaCO₃*(mᴾ*K_mondo(P, Kₘ)*P + sh*wᴾ*P^2) + 0.5*mᴰ*K_mondo(D, Kₘ)*D + sh*D^2*wᴰ + Φ + Φ₂ᴰᴼᶜ - g_GOC_FFᴹ*M - λₚₒ¹*GOC     #40   assumed that there is a typo in the D^2 instead of D^3
+    return (σᴹ*(∑gᴹ + ∑g_FFᴹ)*M + rᴹ*bₘ^T*(concentration_limitation(M, Kₘ) + 3*ΔO₂(O₂, bgc))*M 
+            + Pᵤₚᴹ + 0.5*R_CaCO₃*(mᴾ*concentration_limitation(P, Kₘ)*P + sh*wᴾ*P^2) + 0.5*mᴰ*concentration_limitation(D, Kₘ)*D
+            + sh*D^2*wᴰ + Φ + Φ₂ᴰᴼᶜ - g_GOC_FFᴹ*M - λₚₒ¹*GOC)     #40 
 end

src/Models/AdvectedPopulations/PISCES/calcite.jl

@@ -23,7 +23,7 @@ end
     P₂ = I₂(P, Pₘₐₓ)
     Lₙᴾ = Lᴾ(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, bgc)[5]
     Kₙₕ₄ᴾ = Kᵢᴶ(Kₙₕ₄ᴾᵐⁱⁿ, P₁, P₂, Sᵣₐₜᴾ)
-    Lₗᵢₘᶜᵃᶜᵒ³ = min(Lₙᴾ, K_mondo(Fe, 6e-11), K_mondo(PO₄, Kₙₕ₄ᴾ))
+    Lₗᵢₘᶜᵃᶜᵒ³ = min(Lₙᴾ, concentration_limitation(Fe, 6e-11), concentration_limitation(PO₄, Kₙₕ₄ᴾ))
     return r_CaCO₃*Lₗᵢₘᶜᵃᶜᵒ³*T*max(1, P/2)*max(0, PAR - 1)*30*(1 + exp((-(T-10)^2)/25))*min(1, 50/(zₘₓₗ + eps(0.0)))/((0.1 + T)*(4 + PAR)*(30 + PAR)) #eq77
 end
 
@@ -35,7 +35,7 @@ end
     ηᴹ = bgc.proportion_of_sinking_grazed_shells.M
     sh = get_sh(z, zₘₓₗ)
 
-    return get_R_CaCO₃(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, T, PAR, zₘₓₗ, bgc)*(ηᶻ*get_grazingᶻ(P, D, POC, T, bgc)[2]*Z+ηᴹ*get_grazingᴹ(P, D, Z, POC, T, bgc)[2]*M + 0.5*(mᴾ*K_mondo(P, Kₘ)*P + sh*wᴾ*P^2)) #eq76
+    return get_R_CaCO₃(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, T, PAR, zₘₓₗ, bgc)*(ηᶻ*get_grazingᶻ(P, D, POC, T, bgc)[2]*Z+ηᴹ*get_grazingᴹ(P, D, Z, POC, T, bgc)[2]*M + 0.5*(mᴾ*concentration_limitation(P, Kₘ)*P + sh*wᴾ*P^2)) #eq76
 end
 
 @inline function (bgc::PISCES)(::Val{:CaCO₃}, x, y, z, t, P, D, Z, M, Pᶜʰˡ, Dᶜʰˡ, Pᶠᵉ, Dᶠᵉ, Dˢⁱ, DOC, POC, GOC, SFe, BFe, PSi, NO₃, NH₄, PO₄, Fe, Si, CaCO₃, DIC, Alk, O₂, T, zₘₓₗ, zₑᵤ, Si̅, D_dust, Ω, PAR, PAR¹, PAR², PAR³)

src/Models/AdvectedPopulations/PISCES/carbonate_system.jl

@@ -45,7 +45,9 @@
     μᴾ = μᴵ(P, Pᶜʰˡ, PARᴾ, L_day, T, αᴾ, Lₗᵢₘᴾ, zₘₓₗ, zₑᵤ, t_darkᴾ, bgc)
     μᴰ = μᴵ(D, Dᶜʰˡ, PARᴰ, L_day, T, αᴰ, Lₗᵢₘᴰ, zₘₓₗ, zₑᵤ, t_darkᴰ, bgc)
 
-    return γᶻ*(1 - eᶻ - σᶻ)*∑gᶻ*Z + γᴹ*(1 - eᴹ - σᴹ)*(∑gᴹ + ∑g_FFᴹ)*M + γᴹ*Rᵤₚ(M, T, bgc) + get_Remin(O₂, NO₃, PO₄, NH₄, DOC, T, bFe, Bact, bgc) + get_Denit(NO₃, PO₄, NH₄, DOC, O₂, T, bFe, Bact, bgc) + λ_CaCO₃¹(CaCO₃, bgc, Ω)*CaCO₃ - P_CaCO₃(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, D, Z, M, POC, T, PAR, zₘₓₗ, z, bgc) - μᴰ*D - μᴾ*P #eq59
+    return (γᶻ*(1 - eᶻ - σᶻ)*∑gᶻ*Z + γᴹ*(1 - eᴹ - σᴹ)*(∑gᴹ + ∑g_FFᴹ)*M + γᴹ*Rᵤₚ(M, T, bgc) + 
+    get_Remin(O₂, NO₃, PO₄, NH₄, DOC, T, bFe, Bact, bgc) + get_Denit(NO₃, PO₄, NH₄, DOC, O₂, T, bFe, Bact, bgc) + 
+    λ_CaCO₃¹(CaCO₃, bgc, Ω)*CaCO₃ - P_CaCO₃(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, D, Z, M, POC, T, PAR, zₘₓₗ, z, bgc) - μᴰ*D - μᴾ*P) #eq59
 end
 
 @inline function (bgc::PISCES)(::Val{:Alk}, x, y, z, t, P, D, Z, M, Pᶜʰˡ, Dᶜʰˡ, Pᶠᵉ, Dᶠᵉ, Dˢⁱ, DOC, POC, GOC, SFe, BFe, PSi, NO₃, NH₄, PO₄, Fe, Si, CaCO₃, DIC, Alk, O₂, T, zₘₓₗ, zₑᵤ, Si̅, D_dust, Ω, PAR, PAR¹, PAR², PAR³) # eq59
@@ -92,5 +94,8 @@ end
     μₙₕ₄ᴾ = get_μₙₕ₄ᴾ(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, T, zₘₓₗ, zₑᵤ, L_day, PARᴾ, t_darkᴾ, Si̅, bgc)
     μₙₕ₄ᴰ = get_μₙₕ₄ᴰ(D, PO₄, NO₃, NH₄, Si, Dᶜʰˡ, Dᶠᵉ, T, zₘₓₗ, zₑᵤ, L_day, PARᴰ, t_darkᴰ, Si̅, bgc)
 
-    return θᴺᶜ*get_Remin(O₂, NO₃, PO₄, NH₄, DOC, T, bFe, Bact, bgc) + θᴺᶜ*(rₙₒ₃¹ + 1)*get_Denit(NO₃, PO₄, NH₄, DOC, O₂, T, bFe, Bact, bgc) + θᴺᶜ*γᶻ*(1 - eᶻ - σᶻ)*∑gᶻ*Z + θᴺᶜ*γᴹ*(1 - eᴹ - σᴹ)*(∑gᴹ + ∑g_FFᴹ + θᴺᶜ*γᴹ*Rᵤₚ(M, T, bgc))*M + θᴺᶜ*μₙₒ₃ᴾ*P + θᴺᶜ*μₙₒ₃ᴰ*D + θᴺᶜ*N_fix(bFe, PO₄, T, P, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, PAR, bgc) + 2*λ_CaCO₃¹(CaCO₃, bgc, Ω)*CaCO₃ + θᴺᶜ*ΔO₂(O₂, bgc)*(rₙₕ₄¹ - 1)*λₙₕ₄*NH₄ - θᴺᶜ*μₙₕ₄ᴾ*P - θᴺᶜ*μₙₕ₄ᴰ*D- 2*θᴺᶜ*Nitrif(NH₄, O₂, λₙₕ₄, PAR, bgc) - 2*P_CaCO₃(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, D, Z, M, POC, T, PAR, zₘₓₗ, z, bgc)
+    return (θᴺᶜ*get_Remin(O₂, NO₃, PO₄, NH₄, DOC, T, bFe, Bact, bgc) + θᴺᶜ*(rₙₒ₃¹ + 1)*get_Denit(NO₃, PO₄, NH₄, DOC, O₂, T, bFe, Bact, bgc) + 
+    θᴺᶜ*γᶻ*(1 - eᶻ - σᶻ)*∑gᶻ*Z + θᴺᶜ*γᴹ*(1 - eᴹ - σᴹ)*(∑gᴹ + ∑g_FFᴹ + θᴺᶜ*γᴹ*Rᵤₚ(M, T, bgc))*M + θᴺᶜ*μₙₒ₃ᴾ*P + θᴺᶜ*μₙₒ₃ᴰ*D + 
+    θᴺᶜ*N_fix(bFe, PO₄, T, P, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, PAR, bgc) + 2*λ_CaCO₃¹(CaCO₃, bgc, Ω)*CaCO₃ + θᴺᶜ*ΔO₂(O₂, bgc)*(rₙₕ₄¹ - 1)*λₙₕ₄*NH₄ - 
+    θᴺᶜ*μₙₕ₄ᴾ*P - θᴺᶜ*μₙₕ₄ᴰ*D- 2*θᴺᶜ*Nitrif(NH₄, O₂, λₙₕ₄, PAR, bgc) - 2*P_CaCO₃(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, D, Z, M, POC, T, PAR, zₘₓₗ, z, bgc))
 end

src/Models/AdvectedPopulations/PISCES/iron.jl

@@ -7,7 +7,7 @@
 
 @inline function get_Fe¹(Fe, DOC, T)
     Lₜ = max(0.09*(DOC + 40) - 3, 0.6) # bgc.total_concentration_of_iron_ligands
-    K_eqᶠᵉ = 10^(16.27 - 1565.7/max(T + 273.15, 5)) #check this value
+    K_eqᶠᵉ = exp(16.27 - 1565.7/max(T + 273.15, 5)) #check this value
     Δ = 1 +  K_eqᶠᵉ*Lₜ -  K_eqᶠᵉ*Fe
 
     return (-Δ + sqrt(Δ^2 + 4*K_eqᶠᵉ*Fe))/(2*K_eqᶠᵉ + eps(0.0)) #eq65
@@ -34,11 +34,11 @@ end
 @inline function Aggfe(Fe, DOC, T, bgc)
     λᶠᵉ = 1e-3 * bgc.slope_of_scavenging_rate_of_iron #parameter not defined in parameter list. Assumed scaled version λ_Fe to fit dimensions of Fe¹.
     Lₜ = max(0.09*(DOC + 40) - 3, 0.6)
-    return 1000*λᶠᵉ*max(0, Fe - Lₜ)*get_Fe¹(Fe, DOC, T)
+    return λᶠᵉ*max(0, Fe - Lₜ)*get_Fe¹(Fe, DOC, T)
 end
 
 @inline function get_Bactfe(μₘₐₓ⁰, z, Z, M, Fe, DOC, PO₄, NO₃, NH₄, bFe, T, zₘₐₓ, bgc)
-    K_Feᴮ¹ = bgc.Fe_half_saturation_const_for_PLACEHOLDER
+    K_Feᴮ¹ = bgc.Fe_half_saturation_const_for_Bacteria
     θₘₐₓᶠᵉᵇᵃᶜᵗ = bgc.max_FeC_ratio_of_bacteria
     Bact = get_Bact(zₘₐₓ, z, Z, M) 
     Lₗᵢₘᵇᵃᶜᵗ = Lᵇᵃᶜᵗ(DOC, PO₄, NO₃, NH₄, bFe, bgc)[2]
@@ -49,6 +49,7 @@ end
     σᶻ = bgc.non_assimilated_fraction.Z
     γᴹ = bgc.excretion_as_DOM.M
     σᴹ = bgc.non_assimilated_fraction.M
+    eₘₐₓᶻ, eₘₐₓᴹ = bgc.max_growth_efficiency_of_zooplankton
     δᴾ = bgc.exudation_of_DOC.P
     δᴰ = bgc.exudation_of_DOC.D
     θᶠᵉᶻ = bgc.FeC_ratio_of_zooplankton
@@ -98,16 +99,14 @@ end
     g_GOC_FFᴹ = g_FF*bₘ^T*w_GOC*GOC 
 
     #Gross growth efficiency
-    eₙᶻ = get_eₙᴶ(gₚᶻ, g_Dᶻ, gₚₒᶻ, 0, Pᶠᵉ, Dᶠᵉ, SFe, P, D, POC, bgc)
-    eₙᴹ = get_eₙᴶ(gₚᴹ, g_Dᴹ, gₚₒᴹ, g_Zᴹ, Pᶠᵉ, Dᶠᵉ, SFe, P, D, POC, bgc)
-
-    #println("Sum of positive terms in Iron is ",  max(0, (1-σᶻ)*(∑θᶠᵉⁱgᵢᶻ/(∑gᶻ + eps(0.0)) - eₙᶻ*θᶠᵉᶻ))*∑gᶻ*Z + max(0, (1-σᴹ)*(∑θᶠᵉⁱgᵢᴹ + θᶠᵉᴾᴼᶜ*gₚₒ_FF + θᶠᵉᴳᴼᶜ*g_GOC_FFᴹ )/(∑gᴹ+∑g_FFᴹ + eps(0.0)) - eₙᴹ*θᶠᵉᶻ)*(∑gᴹ+∑g_FFᴹ)*M + γᴹ*θᶠᵉᶻ*Rᵤₚ(M, T, bgc) + λₚₒ¹*SFe)
-    #println("term a1 = $(max(0, (1-σᶻ)*(∑θᶠᵉⁱgᵢᶻ/(∑gᶻ + eps(0.0)) - eₙᶻ*θᶠᵉᶻ))*∑gᶻ*Z), term b = $(max(0, (1-σᴹ)*(∑θᶠᵉⁱgᵢᴹ + θᶠᵉᴾᴼᶜ*gₚₒ_FF + θᶠᵉᴳᴼᶜ*g_GOC_FFᴹ )/(∑gᴹ+∑g_FFᴹ + eps(0.0)) - eₙᴹ*θᶠᵉᶻ)*(∑gᴹ+∑g_FFᴹ)*M ), term c = $( γᴹ*θᶠᵉᶻ*Rᵤₚ(M, T, bgc)), term d = $(λₚₒ¹*SFe)")
-    #println("Sum of negative terms is ", (1 - δᴾ)*μᴾᶠᵉ*P + (1 - δᴰ)*μᴰᶠᵉ*D + Scav(POC, GOC, CaCO₃, PSi, D_dust, DOC, T, Fe, bgc) + Cgfe1(sh, Fe, POC, DOC, T, bgc) + Cgfe2(sh, Fe, T, DOC, GOC, bgc) + Aggfe(Fe, DOC, T, bgc) + Bactfe)
-    #println("term 1 = $((1 - δᴾ)*μᴾᶠᵉ*P ), term 2 = $((1 - δᴰ)*μᴰᶠᵉ*D), term 3 = $(Scav(POC, GOC, CaCO₃, PSi, D_dust, DOC, T, Fe, bgc)), term 4 = $(Cgfe1(sh, Fe, POC, DOC, T, bgc)), term 5 = $(Cgfe2(sh, Fe, T, DOC, GOC, bgc)), term 6 = $(Aggfe(Fe, DOC, T, bgc)), term 6 = $(Bactfe)")
-    #println("--")
-    #println("Total change is ",max(0, (1-σᶻ)*(∑θᶠᵉⁱgᵢᶻ/(∑gᶻ + eps(0.0)) - eₙᶻ*θᶠᵉᶻ))*∑gᶻ*Z + max(0, (1-σᴹ)*(∑θᶠᵉⁱgᵢᴹ + θᶠᵉᴾᴼᶜ*gₚₒ_FF + θᶠᵉᴳᴼᶜ*g_GOC_FFᴹ )/(∑gᴹ+∑g_FFᴹ + eps(0.0)) - eₙᴹ*θᶠᵉᶻ)*(∑gᴹ+∑g_FFᴹ)*M + γᴹ*θᶠᵉᶻ*Rᵤₚ(M, T, bgc) + λₚₒ¹*SFe - (1 - δᴾ)*μᴾᶠᵉ*P - (1 - δᴰ)*μᴰᶠᵉ*D - Scav(POC, GOC, CaCO₃, PSi, D_dust, DOC, T, Fe, bgc) - Cgfe1(sh, Fe, POC, DOC, T, bgc) - Cgfe2(sh, Fe, T, DOC, GOC, bgc) - Aggfe(Fe, DOC, T, bgc) - Bactfe)
-    #println("-------------------------------------")
-
-    return max(0, (1-σᶻ)*(∑θᶠᵉⁱgᵢᶻ/(∑gᶻ + eps(0.0)) - eₙᶻ*θᶠᵉᶻ))*∑gᶻ*Z + max(0, (1-σᴹ)*(∑θᶠᵉⁱgᵢᴹ + θᶠᵉᴾᴼᶜ*gₚₒ_FF + θᶠᵉᴳᴼᶜ*g_GOC_FFᴹ )/(∑gᴹ+∑g_FFᴹ + eps(0.0)) - eₙᴹ*θᶠᵉᶻ)*(∑gᴹ+∑g_FFᴹ)*M + γᴹ*θᶠᵉᶻ*Rᵤₚ(M, T, bgc) + λₚₒ¹*SFe - (1 - δᴾ)*μᴾᶠᵉ*P - (1 - δᴰ)*μᴰᶠᵉ*D - Scav(POC, GOC, CaCO₃, PSi, D_dust, DOC, T, Fe, bgc) - Cgfe1(sh, Fe, POC, DOC, T, bgc) - Cgfe2(sh, Fe, T, DOC, GOC, bgc) - Aggfe(Fe, DOC, T, bgc) - Bactfe
+    eᶻ = eᴶ(eₘₐₓᶻ, σᶻ, gₚᶻ, g_Dᶻ, gₚₒᶻ, 0, Pᶠᵉ, Dᶠᵉ, SFe, P, D, POC, bgc) #eₘₐₓᶻ used in paper but changed here to be consistent with eqs 24, 28
+    eᴹ =  eᴶ(eₘₐₓᴹ, σᴹ, gₚᴹ, g_Dᴹ, gₚₒᴹ, g_Zᴹ,Pᶠᵉ, Dᶠᵉ, SFe, P, D, POC, bgc)
+
+   # println("Scav = $(Scav(POC, GOC, CaCO₃, PSi, D_dust, DOC, T, Fe, bgc)), Aggfe = $(Aggfe(Fe, DOC, T, bgc) )")
+
+    return (max(0, (1-σᶻ)*(∑θᶠᵉⁱgᵢᶻ/(∑gᶻ + eps(0.0)) - eᶻ*θᶠᵉᶻ))*∑gᶻ*Z 
+            + max(0, (1-σᴹ)*(∑θᶠᵉⁱgᵢᴹ + θᶠᵉᴾᴼᶜ*gₚₒ_FF + θᶠᵉᴳᴼᶜ*g_GOC_FFᴹ )/(∑gᴹ+∑g_FFᴹ + eps(0.0)) - eᴹ*θᶠᵉᶻ)*(∑gᴹ+∑g_FFᴹ)*M 
+            + γᴹ*θᶠᵉᶻ*Rᵤₚ(M, T, bgc) + λₚₒ¹*SFe 
+            - (1 - δᴾ)*μᴾᶠᵉ*P - (1 - δᴰ)*μᴰᶠᵉ*D 
+            - Scav(POC, GOC, CaCO₃, PSi, D_dust, DOC, T, Fe, bgc) - Cgfe1(sh, Fe, POC, DOC, T, bgc) - Cgfe2(sh, Fe, T, DOC, GOC, bgc) - Aggfe(Fe, DOC, T, bgc) - Bactfe)
 end