Nitrogen fixation update (#84)

* rename functions and arguments for better understanding

* fix respiration

* fix Fe uptake default parameter

* fix typos in boundary conditions

* update iron uptake, C fixation processes

* update compat and bump version

Author: zhenwu0728

Manifest.toml

@@ -1,19 +1,23 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.11.0"
+julia_version = "1.11.4"
 manifest_format = "2.0"
-project_hash = "a73a8d7764dc5026b1a2046f8ed34a77d5030cd7"
+project_hash = "e6539aca5cbaeb4dfe24d3f271a1eaf6e3cb52c3"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
+git-tree-sha1 = "cd8b948862abee8f3d3e9b73a102a9ca924debb0"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.4"
-weakdeps = ["StaticArrays"]
+version = "4.2.0"
 
     [deps.Adapt.extensions]
+    AdaptSparseArraysExt = "SparseArrays"
     AdaptStaticArraysExt = "StaticArrays"
 
+    [deps.Adapt.weakdeps]
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 version = "1.1.2"
@@ -24,19 +28,26 @@ version = "1.11.0"
 
 [[deps.Atomix]]
 deps = ["UnsafeAtomics"]
-git-tree-sha1 = "c06a868224ecba914baa6942988e2f2aade419be"
+git-tree-sha1 = "b5bb4dc6248fde467be2a863eb8452993e74d402"
 uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
-version = "0.1.0"
+version = "1.1.1"
+
+    [deps.Atomix.extensions]
+    AtomixCUDAExt = "CUDA"
+    AtomixMetalExt = "Metal"
+    AtomixOpenCLExt = "OpenCL"
+    AtomixoneAPIExt = "oneAPI"
+
+    [deps.Atomix.weakdeps]
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+    OpenCL = "08131aa3-fb12-5dee-8b74-c09406e224a2"
+    oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b"
 
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 version = "1.11.0"
 
-[[deps.CEnum]]
-git-tree-sha1 = "389ad5c84de1ae7cf0e28e381131c98ea87d54fc"
-uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
-version = "0.5.0"
-
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
@@ -79,9 +90,15 @@ version = "1.6.0"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "82d8afa92ecf4b52d78d869f038ebfb881267322"
+git-tree-sha1 = "b66970a70db13f45b7e57fbda1736e1cf72174ea"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.16.3"
+version = "1.17.0"
+
+    [deps.FileIO.extensions]
+    HTTPExt = "HTTP"
+
+    [deps.FileIO.weakdeps]
+    HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
@@ -99,21 +116,15 @@ version = "1.0.0"
 
 [[deps.JLD2]]
 deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "PrecompileTools", "Requires", "TranscodingStreams"]
-git-tree-sha1 = "aeab5c68eb2cf326619bf71235d8f4561c62fe22"
+git-tree-sha1 = "91d501cb908df6f134352ad73cde5efc50138279"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.5.5"
-
-[[deps.JLLWrappers]]
-deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "f389674c99bfcde17dc57454011aa44d5a260a40"
-uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.6.0"
+version = "0.5.11"
 
 [[deps.KernelAbstractions]]
-deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
-git-tree-sha1 = "04e52f596d0871fa3890170fa79cb15e481e4cd8"
+deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs"]
+git-tree-sha1 = "80d268b2f4e396edc5ea004d1e0f569231c71e9e"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.28"
+version = "0.9.34"
 
     [deps.KernelAbstractions.extensions]
     EnzymeExt = "EnzymeCore"
@@ -125,29 +136,6 @@ version = "0.9.28"
     LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-[[deps.LLVM]]
-deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"]
-git-tree-sha1 = "4ad43cb0a4bb5e5b1506e1d1f48646d7e0c80363"
-uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
-version = "9.1.2"
-
-    [deps.LLVM.extensions]
-    BFloat16sExt = "BFloat16s"
-
-    [deps.LLVM.weakdeps]
-    BFloat16s = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
-
-[[deps.LLVMExtra_jll]]
-deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
-git-tree-sha1 = "05a8bd5a42309a9ec82f700876903abce1017dd3"
-uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
-version = "0.0.34+0"
-
-[[deps.LazyArtifacts]]
-deps = ["Artifacts", "Pkg"]
-uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
-version = "1.11.0"
-
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
@@ -187,10 +175,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 version = "1.11.0"
 
 [[deps.MacroTools]]
-deps = ["Markdown", "Random"]
-git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df"
+git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472"
 uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.13"
+version = "0.5.15"
 
 [[deps.Markdown]]
 deps = ["Base64"]
@@ -220,9 +207,9 @@ uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
 version = "0.3.27+1"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
+git-tree-sha1 = "cc4054e898b852042d7b503313f7ad03de99c3dd"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.3"
+version = "1.8.0"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"]
@@ -259,9 +246,9 @@ version = "1.11.0"
 
 [[deps.Requires]]
 deps = ["UUIDs"]
-git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
+git-tree-sha1 = "62389eeff14780bfe55195b7204c0d8738436d64"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.3.0"
+version = "1.3.1"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -273,9 +260,9 @@ version = "1.11.0"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50"
+git-tree-sha1 = "0feb6b9031bd5c51f9072393eb5ab3efd31bf9e4"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.7"
+version = "1.9.13"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
@@ -304,19 +291,22 @@ version = "1.11.1"
 
 [[deps.StructArrays]]
 deps = ["ConstructionBase", "DataAPI", "Tables"]
-git-tree-sha1 = "f4dc295e983502292c4c3f951dbb4e985e35b3be"
+git-tree-sha1 = "5a3a31c41e15a1e042d60f2f4942adccba05d3c9"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-version = "0.6.18"
+version = "0.7.0"
 
     [deps.StructArrays.extensions]
     StructArraysAdaptExt = "Adapt"
-    StructArraysGPUArraysCoreExt = "GPUArraysCore"
+    StructArraysGPUArraysCoreExt = ["GPUArraysCore", "KernelAbstractions"]
+    StructArraysLinearAlgebraExt = "LinearAlgebra"
     StructArraysSparseArraysExt = "SparseArrays"
     StructArraysStaticArraysExt = "StaticArrays"
 
     [deps.StructArrays.weakdeps]
     Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
     StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
@@ -357,15 +347,15 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 version = "1.11.0"
 
 [[deps.UnsafeAtomics]]
-git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278"
+git-tree-sha1 = "b13c4edda90890e5b04ba24e20a310fbe6f249ff"
 uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f"
-version = "0.2.1"
+version = "0.3.0"
+
+    [deps.UnsafeAtomics.extensions]
+    UnsafeAtomicsLLVM = ["LLVM"]
 
-[[deps.UnsafeAtomicsLLVM]]
-deps = ["LLVM", "UnsafeAtomics"]
-git-tree-sha1 = "2d17fabcd17e67d7625ce9c531fb9f40b7c42ce4"
-uuid = "d80eeb9a-aca5-4d75-85e5-170c8b632249"
-version = "0.2.1"
+    [deps.UnsafeAtomics.weakdeps]
+    LLVM = "929cbde3-209d-540e-8aea-75f648917ca0"
 
 [[deps.Zlib_jll]]
 deps = ["Libdl"]

Project.toml

@@ -1,7 +1,7 @@
 name = "PlanktonIndividuals"
 uuid = "dfeb4ede-c75e-11e9-050d-41ab7a730e9e"
 authors = ["zhenwu0728 <zhenwu@mit.edu>", "gaelforget <gforget@mit.edu>"]
-version = "0.7.3"
+version = "0.7.4"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -31,9 +31,9 @@ Metal = "^1"
 GPUArrays = "^10,11"
 JLD2 = "^0.5"
 KernelAbstractions = "^0.9"
-StructArrays = "^0.6"
+StructArrays = "^0.7"
 Statistics = "^1"
-julia = "^1.7"
+julia = "^1.10"
 
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"

src/Diagnostics/diagnostics_kernels.jl

@@ -20,7 +20,7 @@ function diags_spcs!(diags_sp, plank, ac, x, y, z, mode::AbstractMode, arch::Arc
     elseif isa(mode, MacroMolecularMode)
         diags = (:PS, :VDOC, :VHN4, :VNO3, :VPO4, :S_PRO, :S_DNA, :S_RNA, :resp, :ρChl, :CH, :NST, :PST, :PRO, :DNA, :RNA, :Chl)
     elseif isa(mode, IronEnergyMode)
-        diags = (:PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :NF, :ENF, :BS, :VNH4, :VNO3, :VPO4, :VFe, :PS2ST, :ST2PS, :NR2ST, :ST2NR, :NF2ST, :ST2NF, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :qFePS, :qFeNR, :qFeNF, :Chl)
+        diags = (:PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :NF, :ENF, :BS, :VNH4, :VNO3, :VPO4, :VFe, :PS2ST, :ST2PS, :NR2ST, :ST2NR, :NF2ST, :ST2NF, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :qFePS, :qFeNR, :qFeNF, :Chl, :tdark)
     end
 
     for diag in keys(diags_sp)

src/Diagnostics/diagnostics_struct.jl

@@ -110,7 +110,7 @@ function plank_avail_diags(mode::AbstractMode)
     elseif isa(mode, MacroMolecularMode)
         plank_avail = (:num, :graz, :mort, :dvid, :PS, :VDOC, :VNH4, :VNO3, :VPO4, :resp, :ρChl, :S_PRO, :S_DNA, :S_RNA, :exu, :CH, :NST, :PST, :PRO, :DNA, :RNA, :Chl)
     elseif isa(mode, IronEnergyMode)
-        plank_avail = (:num, :graz, :mort, :dvid, :PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :NF, :ENF, :BS, :VNH4, :VNO3, :VPO4, :VFe, :PS2ST, :ST2PS, :NR2ST, :ST2NR, :NF2ST, :ST2NF, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :qFePS, :qFeNR, :qFeNF, :Chl)
+        plank_avail = (:num, :graz, :mort, :dvid, :PS, :CF, :ECF, :RS, :ERS, :NR, :ENR, :NF, :ENF, :BS, :VNH4, :VNO3, :VPO4, :VFe, :PS2ST, :ST2PS, :NR2ST, :ST2NR, :NF2ST, :ST2NF, :Bm, :En, :CH, :qNO3, :qNH4, :qP, :qFe, :qFePS, :qFeNR, :qFeNF, :Chl, :tdark)
     end
     return plank_avail
 end

src/Fields/apply_bcs.jl

@@ -2,8 +2,7 @@
 ##### apply flux boundary conditions to tracer tendency `Gc`
 ##### 
 
-# a positive west flux is associated with an *increase* in `Gc` near the west boundary.
-# same for south and bottom.
+# positive west/south/bottom bc values are associated with an *increase* in `Gc` near the boundary.
 @kernel function apply_west_bcs_kernel!(Gc, grid, west_bc, iter, ΔT)
     j, k = @index(Global, NTuple)
     jj = j + grid.Hy
@@ -25,8 +24,7 @@ end
     @inbounds Gc[ii, jj, grid.Nz+grid.Hz] += getbc(bottom_bc, i, j, iter) * ΔT * Az(ii, jj, grid.Nz+grid.Hz, grid) / volume(ii, jj, grid.Nz+grid.Hz, grid)
 end
 
-# a positive east flux is associated with an *decrease* in `Gc` near the east boundary.
-# same for north and top.
+# positive east/north/top bc values are associated with a *decrease* in `Gc` near the boundary.
 @kernel function apply_east_bcs_kernel!(Gc, grid, east_bc, iter, ΔT)
     j, k = @index(Global, NTuple)
     jj = j + grid.Hy

src/Fields/boundary_conditions.jl

@@ -57,14 +57,14 @@ function validate_bc(bc, bc_size, nΔT)
         if size(bc) == bc_size
             return nothing
         else
-            throw(ArgumentError("BC west: grid mismatch, size(bc) must equal to $(bc_size) for a constant flux boundary condition."))
+            throw(ArgumentError("BC: grid mismatch, size(bc) must equal to $(bc_size) for a constant flux boundary condition."))
         end
     end
     if typeof(bc) <: AbstractArray{eltype(bc),3}
         if size(bc) == (bc_size..., nΔT)
             return nothing
         else
-            throw(ArgumentError("BC west: grid mismatch, size(bc) must equal to $((bc_size..., nΔT)) for a time-dependent flux boundary condition."))
+            throw(ArgumentError("BC: grid mismatch, size(bc) must equal to $((bc_size..., nΔT)) for a time-dependent flux boundary condition."))
         end
     end
 end

src/Parameters/param_default.jl

@@ -148,7 +148,7 @@ Generate default phytoplankton parameter values based on `AbstractMode` and spec
 """
 function phyt_params_default(N::Int64, mode::IronEnergyMode)
     params=Dict(
-        "Nsuper"    => [1],       # Number of phyto cells each super individual represents
+        "Nsuper"    => [1],       # Number of phyto cells each super individual represents (cells)
         "Cquota"    => [1.8e-11], # C quota of phyto cells at size = 1.0 (mmolC/cell)
         "SA"        => [0.18],    # Surface area (μm²) Prochlorococcus
         "mean"      => [1.2],     # Mean of the normal distribution of initial phyto individuals
@@ -187,7 +187,7 @@ function phyt_params_default(N::Int64, mode::IronEnergyMode)
         "KsatNH4"   => [0.005],   # Half-saturation coeff (mmolN/m³)
         "KsatNO3"   => [0.010],   # Half-saturation coeff (mmolN/m³)
         "KsatPO4"   => [0.003],   # Half-saturation coeff (mmolP/m³)
-        "KSAFe"     => [1.16e-17],# Surface-area specific iron uptake rate (m³/μm²/second)
+        "KSAFe"     => [2.77e-19],# Surface-area specific iron uptake rate (m³/μm²/cell/second)
         "qNO3max"   => [0.25],    # Maximum NO3 quota in cell (mmolN/mmolC)
         "qNH4max"   => [0.25],    # Maximum NH4 quota in cell (mmolN/mmolC)
         "qPmax"     => [0.02],    # Maximum P quota in cell (mmolP/mmolC)

src/Plankton/CarbonMode/growth_kernels.jl

@@ -1,16 +1,16 @@
 ##### temperature function for photosynthesis
-@inline function tempFunc_PS(temp, p)
-    x = temp - p.Topt; xmax = p.Tmax - p.Topt
+@inline function tempFunc_PS(T, p)
+    x = T - p.Topt; xmax = p.Tmax - p.Topt
     regT = shape_func_dec(x, xmax, 1.0f-1)
-    k = exp(-p.Ea/(8.3145f0*(temp+273.15f0))) * regT
+    k = exp(-p.Ea/(8.3145f0*(T+273.15f0))) * regT
     k = max(0.0f0, k)
     OGT_rate = exp(-p.Ea/(8.3145f0*(p.Topt+273.15f0)))
     return min(1.0f0, k/OGT_rate)
 end
 
 ##### temperature function for metabolic rates
-@inline function tempFunc(temp, p)
-    k = exp(-p.Ea/(8.3145f0*(temp+273.15f0)))
+@inline function tempFunc(T, p)
+    k = exp(-p.Ea/(8.3145f0*(T+273.15f0)))
     k = max(0.0f0, k)
     OGT_rate = exp(-p.Ea/(8.3145f0*(p.Topt+273.15f0)))
     return min(1.0f0, k/OGT_rate)
@@ -20,19 +20,18 @@ end
 ##### only functional when damaged biomass is greater than 0.0
 ##### use Bd/Bm as the allocation, more damaged biomass means more allocation to repair
 @inline function gamma_alloc(Bm, Bd, p)
-    #γ = max(0.0, temp - p.Topt) / (p.Tmax - p.Topt) * isless(0.0, Bd)
     γ = Bd / max(1.0f-30, Bm) * isless(0.0f0, Bd)
     γ = min(1.0f0, γ) * isequal(1.0f0, p.thermal)
     return γ
 end
 
 ##### calculate photosynthesis rate (mmolC/individual/second)
-@inline function calc_photosynthesis(par, temp, Chl, Bm, Bd, p)
+@inline function calc_photosynthesis(par, T, Chl, Bm, Bd, p)
     αI  = par * p.α * p.Φ
     if p.thermal == 1.0f0
-        PCm = p.PCmax * tempFunc(temp, p)
+        PCm = p.PCmax * tempFunc(T, p)
     else
-        PCm = p.PCmax * tempFunc_PS(temp, p)
+        PCm = p.PCmax * tempFunc_PS(T, p)
     end
     PS  = PCm * (1.0f0 - exp(-αI / max(1.0f-30, PCm) * Chl / max(1.0f-30, Bm))) * max(0.0f0, Bm - Bd)
     return PS

src/Plankton/IronEnergyMode/division_death.jl

@@ -64,6 +64,7 @@ end
         @inbounds plank.dvid[idx[i]] = plank.dvid[i]
         @inbounds plank.graz[idx[i]] = plank.graz[i]
         @inbounds plank.mort[idx[i]] = plank.mort[i]
+        @inbounds plank.tdark[idx[i]]= plank.tdark[i]
     end
 end
 function copy_daughter_individuals!(plank, con, idx::AbstractArray{Int,1}, arch)