double check nutrient tracers are non-negative

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
 

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)

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
 

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)
 

src/Plankton/QuotaMode/growth_kernels.jl

@@ -23,36 +23,38 @@ 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)
     regQP = shape_func_dec(Qp, p.Pqmax, 1.0f-4)
     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
 

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)