From 3df378d00243329631923733627d3f09f3d87209 Mon Sep 17 00:00:00 2001 From: Zhen Wu Date: Thu, 27 Jun 2024 23:41:25 +0800 Subject: [PATCH] double check nutrient tracers are non-negative --- src/Plankton/IronEnergyMode/growth_kernels.jl | 30 +++++++++---------- .../IronEnergyMode/plankton_growth.jl | 2 +- .../MacroMolecularMode/growth_kernels.jl | 28 +++++++++-------- .../MacroMolecularMode/plankton_growth.jl | 4 +-- src/Plankton/QuotaMode/growth_kernels.jl | 27 +++++++++-------- src/Plankton/QuotaMode/plankton_growth.jl | 4 +-- 6 files changed, 51 insertions(+), 44 deletions(-) diff --git a/src/Plankton/IronEnergyMode/growth_kernels.jl b/src/Plankton/IronEnergyMode/growth_kernels.jl index 0c1979d1..436317d3 100644 --- a/src/Plankton/IronEnergyMode/growth_kernels.jl +++ b/src/Plankton/IronEnergyMode/growth_kernels.jl @@ -90,52 +90,52 @@ function calc_carbon_fixation!(plank, nuts, p, ΔT, arch::Architecture) end ##### calculate ammonium uptake rate (mmolN/individual/second) -@inline function calc_NH4_uptake(NH4, temp, CH, qNH4, Bm, p, ac) +@inline function calc_NH4_uptake(NH4, temp, CH, qNH4, Bm, pop, p, ac, ΔT) Qn = qNH4/max(1.0f-30, Bm + CH) regQN = shape_func_dec(Qn, p.qNH4max, 1.0f-4) VNH4 = p.VNH4max * regQN * NH4/max(1.0f-30, NH4+p.KsatNH4) * tempFunc(temp, p) * Bm * ac - return VNH4 + return min(VNH4, NH4/ΔT/max(1.0f0,pop)) end ##### calculate nitrate uptake rate (mmolN/individual/second) -@inline function calc_NO3_uptake(NO3, temp, CH, qNO3, Bm, p, ac) +@inline function calc_NO3_uptake(NO3, temp, CH, qNO3, Bm, pop, p, ac, ΔT) Qn = qNO3/max(1.0f-30, Bm + CH) regQN = shape_func_dec(Qn, p.qNO3max, 1.0f-4) VNO3 = p.VNO3max * regQN * NO3/max(1.0f-30, NO3+p.KsatNO3) * tempFunc(temp, p) * Bm * ac - return VNO3 + return min(VNO3, NO3/ΔT/max(1.0f0,pop)) end ##### calculate phosphate uptake rate (mmolN/individual/second) -@inline function calc_P_uptake(PO4, temp, CH, qP, Bm, p, ac) +@inline function calc_P_uptake(PO4, temp, CH, qP, Bm, pop, p, ac, ΔT) Qp = (qP + Bm * p.R_PC)/max(1.0f-30, Bm + CH) regQP = shape_func_dec(Qp, p.qPmax, 1.0f-4) VPO4 = p.VPO4max * regQP * PO4/max(1.0f-30, PO4+p.KsatPO4) * tempFunc(temp, p) * Bm * ac - return VPO4 + return min(VPO4, PO4/ΔT/max(1.0f0,pop)) end ##### calculate iron uptake rate (mmolFe/individual/second) -@inline function calc_Fe_uptake(FeT, qFe, Bm, CH, Sz, p, ac) +@inline function calc_Fe_uptake(FeT, qFe, Bm, CH, Sz, pop, p, ac, ΔT) Qfe = qFe/max(1.0f-30, Bm +CH) regQFe = shape_func_dec(Qfe, p.qFemax, 1.0f-4) SA = p.SA * Sz^(2/3) VFe = p.KSAFe * SA * FeT * regQFe * ac - return VFe + return min(VFe, FeT/ΔT/max(1.0f0,pop)) end -@kernel function calc_nuts_uptake_kernel!(plank, nuts, p) +@kernel function calc_nuts_uptake_kernel!(plank, nuts, p, ΔT) i = @index(Global) @inbounds plank.VNH4[i] = calc_NH4_uptake(nuts.NH4[i], nuts.T[i], plank.CH[i], plank.qNH4[i], - plank.Bm[i], p, plank.ac[i]) + plank.Bm[i], nuts.pop[i], p, plank.ac[i], ΔT) @inbounds plank.VNO3[i] = calc_NO3_uptake(nuts.NO3[i], nuts.T[i], plank.CH[i], plank.qNO3[i], - plank.Bm[i], p, plank.ac[i]) + plank.Bm[i], nuts.pop[i], p, plank.ac[i], ΔT) @inbounds plank.VPO4[i] = calc_P_uptake(nuts.PO4[i], nuts.T[i], plank.CH[i], plank.qP[i], - plank.Bm[i], p, plank.ac[i]) + plank.Bm[i], nuts.pop[i], p, plank.ac[i], ΔT) @inbounds plank.VFe[i] = calc_Fe_uptake(nuts.FeT[i], plank.qFe[i], plank.Bm[i], plank.CH[i], - plank.Sz[i], p, plank.ac[i]) + plank.Sz[i], nuts.pop[i], p, plank.ac[i], ΔT) end -function calc_nuts_uptake!(plank, nuts, p, arch::Architecture) +function calc_nuts_uptake!(plank, nuts, p, ΔT, arch::Architecture) kernel! = calc_nuts_uptake_kernel!(device(arch), 256, (size(plank.ac,1))) - kernel!(plank, nuts, p) + kernel!(plank, nuts, p, ΔT) return nothing end diff --git a/src/Plankton/IronEnergyMode/plankton_growth.jl b/src/Plankton/IronEnergyMode/plankton_growth.jl index 220db75d..b1f601ea 100644 --- a/src/Plankton/IronEnergyMode/plankton_growth.jl +++ b/src/Plankton/IronEnergyMode/plankton_growth.jl @@ -4,7 +4,7 @@ function plankton_growth!(plank, nuts, rnd, p, ΔT, t, arch::Architecture) calc_repiration!(plank,nuts,p, ΔT, arch) update_state_1!(plank, ΔT, arch) calc_carbon_fixation!(plank, nuts, p, ΔT, arch) - calc_nuts_uptake!(plank, nuts, p, arch) + calc_nuts_uptake!(plank, nuts, p, ΔT, arch) update_state_2!(plank, ΔT, arch) calc_NO3_reduc!(plank, nuts, p, ΔT, arch) update_state_3!(plank, ΔT, arch) diff --git a/src/Plankton/MacroMolecularMode/growth_kernels.jl b/src/Plankton/MacroMolecularMode/growth_kernels.jl index 41b3f4b3..6f4f7261 100644 --- a/src/Plankton/MacroMolecularMode/growth_kernels.jl +++ b/src/Plankton/MacroMolecularMode/growth_kernels.jl @@ -23,7 +23,7 @@ end end ##### calculate nutrient uptake rate (mmolN/individual/second) -@inline function calc_NP_uptake(NH4, NO3, PO4, temp, NST, PST, PRO, DNA, RNA, Chl, p, ac) +@inline function calc_NP_uptake(NH4, NO3, PO4, temp, NST, PST, PRO, DNA, RNA, Chl, pop, p, ac, ΔT) N_tot = total_N_biomass(PRO, DNA, RNA, NST, Chl, p) P_tot = total_P_biomass(DNA, RNA, PST, p) R_NST = NST / max(1.0f-30, N_tot) @@ -33,21 +33,24 @@ end VNH4 = p.VNH4max * regQN * NH4/max(1.0f-30, NH4+p.KsatNH4) * tempFunc(temp, p) * PRO * ac VNO3 = p.VNO3max * regQN * NO3/max(1.0f-30, NO3+p.KsatNO3) * tempFunc(temp, p) * PRO * ac VPO4 = p.VPO4max * regQP * PO4/max(1.0f-30, PO4+p.KsatPO4) * tempFunc(temp, p) * PRO * ac - return VNH4, VNO3, VPO4 + return min(VNH4, NH4/ΔT/max(1.0f0,pop)), + min(VNO3, NO3/ΔT/max(1.0f0,pop)), + min(VPO4, PO4/ΔT/max(1.0f0,pop)) end -@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p) +@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p, ΔT) i = @index(Global) @inbounds plank.PS[i] = calc_PS(nuts.par[i], nuts.T[i], plank.Chl[i], plank.PRO[i], p) * plank.ac[i] @inbounds plank.VNH4[i], plank.VNO3[i], plank.VPO4[i] = calc_NP_uptake(nuts.NH4[i], nuts.NO3[i], nuts.PO4[i], nuts.T[i], plank.NST[i], plank.PST[i], plank.PRO[i], - plank.DNA[i], plank.RNA[i], plank.Chl[i], p, plank.ac[i]) + plank.DNA[i], plank.RNA[i], plank.Chl[i], + nuts.pop[i], p, plank.ac[i], ΔT) end -function calc_inorganic_uptake!(plank, nuts, p, arch::Architecture) +function calc_inorganic_uptake!(plank, nuts, p, ΔT, arch::Architecture) kernel! = calc_inorganic_uptake_kernel!(device(arch), 256, (size(plank.ac,1))) - kernel!(plank, nuts, p) + kernel!(plank, nuts, p, ΔT) return nothing end @@ -66,24 +69,25 @@ end ##### calculate DOC uptake rate (mmolC/individual/second) ##### DOC uptake needs support of photosynthesis for at least 5% of total C acquisition. -@inline function calc_DOC_uptake(DOC, temp, CH, PRO, DNA, RNA, Chl, p) +@inline function calc_DOC_uptake(DOC, temp, CH, PRO, DNA, RNA, Chl, pop, p, ΔT) C_tot = total_C_biomass(PRO, DNA, RNA, CH, Chl) R_CH = CH / max(1.0f-30, C_tot) regQ = shape_func_dec(R_CH, p.CHmax, 1.0f-4) VN = p.VDOCmax * regQ * DOC/max(1.0f-30, DOC+p.KsatDOC) * tempFunc(temp, p) * PRO - return VN + return min(VN, DOC/ΔT/max(1.0f0,pop)) end -@kernel function calc_organic_uptake_kernel!(plank, nuts, p) +@kernel function calc_organic_uptake_kernel!(plank, nuts, p, ΔT) i = @index(Global) @inbounds plank.VDOC[i] = calc_DOC_uptake(nuts.DOC[i], nuts.T[i], plank.CH[i], plank.PRO[i], plank.DNA[i], - plank.RNA[i], plank.Chl[i], p) * plank.ac[i] + plank.RNA[i], plank.Chl[i], + nuts.pop[i], p, ΔT) * plank.ac[i] @inbounds plank.VDOC[i] = plank.VDOC[i] * isless(0.05f0, plank.PS[i]/(plank.VDOC[i]+plank.PS[i])) end -function calc_organic_uptake!(plank, nuts, p, arch::Architecture) +function calc_organic_uptake!(plank, nuts, p, ΔT, arch::Architecture) kernel! = calc_organic_uptake_kernel!(device(arch), 256, (size(plank.ac,1))) - kernel!(plank, nuts, p) + kernel!(plank, nuts, p, ΔT) return nothing end diff --git a/src/Plankton/MacroMolecularMode/plankton_growth.jl b/src/Plankton/MacroMolecularMode/plankton_growth.jl index 58fbb7c7..43a009cd 100644 --- a/src/Plankton/MacroMolecularMode/plankton_growth.jl +++ b/src/Plankton/MacroMolecularMode/plankton_growth.jl @@ -1,10 +1,10 @@ ##### update physiological attributes of each individual function plankton_growth!(plank, nuts, rnd, p, ΔT, t, arch::Architecture) - calc_inorganic_uptake!(plank, nuts, p, arch) + calc_inorganic_uptake!(plank, nuts, p, ΔT, arch) update_quotas_1!(plank, ΔT, arch) - calc_organic_uptake!(plank, nuts, p, arch) + calc_organic_uptake!(plank, nuts, p, ΔT, arch) calc_ρChl!(plank, nuts.par, p, arch) diff --git a/src/Plankton/QuotaMode/growth_kernels.jl b/src/Plankton/QuotaMode/growth_kernels.jl index 2c40abeb..35f2e5ae 100644 --- a/src/Plankton/QuotaMode/growth_kernels.jl +++ b/src/Plankton/QuotaMode/growth_kernels.jl @@ -23,7 +23,7 @@ end end ##### calculate nutrient uptake rate (mmolN/individual/second) -@inline function calc_NP_uptake(NH4, NO3, PO4, temp, Cq, Nq, Pq, Bm, p, ac) +@inline function calc_NP_uptake(NH4, NO3, PO4, temp, Cq, Nq, Pq, Bm, pop, p, ac, ΔT) Qn = (Nq + Bm * p.R_NC)/max(1.0f-30, Bm + Cq) Qp = (Pq + Bm * p.R_PC)/max(1.0f-30, Bm + Cq) regQN = shape_func_dec(Qn, p.Nqmax, 1.0f-4) @@ -31,28 +31,30 @@ end VNH4 = p.VNH4max * regQN * NH4/max(1.0f-30, NH4+p.KsatNH4) * tempFunc(temp, p) * Bm * ac VNO3 = p.VNO3max * regQN * NO3/max(1.0f-30, NO3+p.KsatNO3) * tempFunc(temp, p) * Bm * ac VPO4 = p.VPO4max * regQP * PO4/max(1.0f-30, PO4+p.KsatPO4) * tempFunc(temp, p) * Bm * ac - return VNH4, VNO3, VPO4 + return min(VNH4, NH4/ΔT/max(1.0f0,pop)), + min(VNO3, NO3/ΔT/max(1.0f0,pop)), + min(VPO4, PO4/ΔT/max(1.0f0,pop)) end -@inline function calc_DOC_uptake(DOC, temp, Cq, Bm, p) +@inline function calc_DOC_uptake(DOC, temp, Cq, Bm, pop, p, ΔT) Qc = Cq/max(1.0f-30, Bm + Cq) regQ = shape_func_dec(Qc, p.Cqmax, 1.0f-4) VN = p.VDOCmax * regQ * DOC/max(1.0f-30, DOC+p.KsatDOC) * tempFunc(temp, p) * Bm - return VN + return min(VN, DOC/ΔT/max(1.0f0,pop)) end -@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p) +@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p, ΔT) i = @index(Global) @inbounds plank.PS[i] = calc_PS(nuts.par[i], nuts.T[i], plank.Chl[i], plank.Bm[i], p) * plank.ac[i] @inbounds plank.VNH4[i], plank.VNO3[i], plank.VPO4[i] = calc_NP_uptake(nuts.NH4[i], nuts.NO3[i], nuts.PO4[i], nuts.T[i], plank.Cq[i], plank.Nq[i], plank.Pq[i], - plank.Bm[i], p, plank.ac[i]) + plank.Bm[i], nuts.pop[i], p, plank.ac[i], ΔT) end -function calc_inorganic_uptake!(plank, nuts, p, arch::Architecture) +function calc_inorganic_uptake!(plank, nuts, p, ΔT, arch::Architecture) kernel! = calc_inorganic_uptake_kernel!(device(arch), 256, (size(plank.ac,1))) - kernel!(plank, nuts, p) + kernel!(plank, nuts, p, ΔT) return nothing end @@ -71,15 +73,16 @@ end ##### calculate DOC uptake rate (mmolC/individual/second) ##### DOC uptake needs support of photosynthesis for at least 1% of total C acquisition. -@kernel function calc_organic_uptake_kernel!(plank, nuts, p) +@kernel function calc_organic_uptake_kernel!(plank, nuts, p, ΔT) i = @index(Global) - @inbounds plank.VDOC[i] = calc_DOC_uptake(nuts.DOC[i], nuts.T[i], plank.Cq[i], plank.Bm[i], p) * plank.ac[i] + @inbounds plank.VDOC[i] = calc_DOC_uptake(nuts.DOC[i], nuts.T[i], plank.Cq[i], + plank.Bm[i], nuts.pop[i], p, ΔT) * plank.ac[i] @inbounds plank.VDOC[i] = plank.VDOC[i] * isless(1.0f-2, plank.PS[i]/(plank.VDOC[i]+plank.PS[i])) end -function calc_organic_uptake!(plank, nuts, p, arch::Architecture) +function calc_organic_uptake!(plank, nuts, p, ΔT, arch::Architecture) kernel! = calc_organic_uptake_kernel!(device(arch), 256, (size(plank.ac,1))) - kernel!(plank, nuts, p) + kernel!(plank, nuts, p, ΔT) return nothing end diff --git a/src/Plankton/QuotaMode/plankton_growth.jl b/src/Plankton/QuotaMode/plankton_growth.jl index e2f08f9a..7ab77d5e 100644 --- a/src/Plankton/QuotaMode/plankton_growth.jl +++ b/src/Plankton/QuotaMode/plankton_growth.jl @@ -1,10 +1,10 @@ ##### update physiological attributes of each individual function plankton_growth!(plank, nuts, rnd, p, ΔT, t, arch::Architecture) - calc_inorganic_uptake!(plank, nuts, p, arch) + calc_inorganic_uptake!(plank, nuts, p, ΔT, arch) update_quotas_1!(plank, ΔT, arch) - calc_organic_uptake!(plank, nuts, p, arch) + calc_organic_uptake!(plank, nuts, p, ΔT, arch) calc_ρChl!(plank, nuts.par, p, arch)