diff --git a/documentation/docs/pages/Model_description/WOMBATlite_model_description.md b/documentation/docs/pages/Model_description/WOMBATlite_model_description.md
new file mode 100644
index 00000000..e76bc7f1
--- /dev/null
+++ b/documentation/docs/pages/Model_description/WOMBATlite_model_description.md
@@ -0,0 +1,984 @@
+# Description of the WOMBATlite ocean biogeochemical model
+
+```
+ (\___/) .-. .-. .--. .-..-..---. .--. .-----.
+ / o o \ : :.-.: :: ,. :: '' :: .; :: .; :'-. .-'
+ ( " ) : :: :: :: :: :: .. :: .': : : :
+ \__ __/ : '' '' ;: :; :: :; :: .; :: :: : : :
+ '.,'.,' '.__.':_;:_;:___.':_;:_; :_;
+
+ World Ocean Model of Biogeochemistry And Trophic-dynamics (WOMBAT)
+```
+
+_Contact Pearse J. Buchanan for any questions_
+
+_Pearse.Buchanan@csiro.au_
+
+---
+
+## Subroutine - "update_from_source"
+
+The subroutine `generic_WOMBATlite_update_from_source` is the heart of the World Ocean Model of Biogeochemistry And Trophic‑dynamics.
+Its purpose is to apply biological source–sink terms to ocean tracers (nutrients, phytoplankton, zooplankton, detritus, iron and carbon pools)
+at each tracer time‑step. The subroutine is documented internally by a list of numbered steps (see code comments). These steps are:
+
+1. Light attenuation through the water column.
+2. Nutrient limitation of phytoplankton.
+3. Temperature‑dependent autotrophy and heterotrophy.
+4. Light limitation of phytoplankton.
+5. Realized growth rate of phytoplankton.
+6. Growth of chlorophyll.
+7. Phytoplankton uptake of iron.
+8. Iron chemistry (precipitation, scavenging and coagulation).
+9. Mortality and remineralisation.
+10. Zooplankton grazing, growth, excretion and egestion.
+11. CaCO3 calculations.
+12. Tracer tendencies (update tracer concentrations).
+13. Check for conservation of mass.
+14. Additional operations on tracers.
+15. Sinking rate of particulates.
+16. Sedimentary processes.
+
+Below is a step‑by‑step explanation of each section together with the key equations. Variable names in grey follow the Fortran code, while
+variable names in math font are pointers to the equations; i,j,k refer to horizontal and vertical indices; square brackets denote units.
+If a variable is without i,j,k dimensions, this variable is held as a scalar and not an array.
+
+The model carries tracers in [mol kg-1]. Some calculations herein are performed by converting tracers to units of [mmol m-3] or in the case of dissolved iron [µmol m-3].
+However, we stress that all tracer tendency terms are converted back to [mol kg-1 s-1] when sources and sinks are applied.
+
+---
+
+
+### 1. Light attenuation through the water column.
+
+Photosynthetically available radiation (PAR) is split into blue, green and red wavelengths. The incoming visible (photosynthetically available) short wave radiation flux (PAR, [W m-2]) is received from the physical model, and is then split evenly into each of blue, green and red light bands.
+
+At the top (`par_bgr_top(k,b)`, $PAR^{top}$, [W m-2]) and mid‑point (`par_bgr_mid(k,b)`, $PAR^{mid}$, [W m-2]) of each layer `k` we calculate the downward irradiance by exponential decay of each band `b` through the layer thickness (`dzt(i,j,k)`, $\Delta z$, [m]) using band‑specific attenuation coefficients (`ek_bgr(k,b)`, $ex_{bgr}$, [m-1]). These attenuation coefficients are related to the concentration of chlorophyll (`chl`, [mg m-3]), organic detritus (`ndet`, [mg N m-3]) and calcium carbonate (`carb`, [kg m-3]) in the water column. For chlorophyll, attentuation coefficients for each of blue, green and red light are retrieved from the look-up table of [Morel & Maritorena (2001)](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2000JC000319) (their Table 2) that explicitly relates chlorophyll concentration to attenutation rates and accounts for the packaging effect of chlorophyll in larger cells. For organic detritus, attenutation coefficients for blue, green and red light are taken from [Dutkiewicz et al. (2015)](https://bg.copernicus.org/articles/12/4447/2015/bg-12-4447-2015.html) (their Fig. 1b), while for calcium carbonate we take the coefficients defined in [Soja-Wozniak et al. (2019)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JC014998). For both detritus and calcium carbonate, these studies provide concentration normalized attenuation coefficients, which must be multipled against concentrations to retrieve the correct units of [m-1].
+
+As an example, the PAR in the blue band (`b=1`) at the top of level k is computed as
+
+$PAR^{top}(k,1) = PAR^{top}(k-1,1) * e^{(-ex_{bgr}(k-1,1) * \Delta z(k-1))}$
+
+where the total attenutation rate of blue light in the grid cell above `k` is the sum of attenuation due to all particulates in that grid cell, which includes chlorophyll, detritus and calcium carbonate:
+
+$ex_{bgr}(k-1,1) = ex_{chl}(k-1,1) + ex_{det}(k-1,1) + ex_{CaCO_3}(k-1,1)$
+
+where
+
+- $ex_{chl}(k-1,1)$ is the attenuation rate of blue light (`b=1`) in the overlying grid cell (`k=k-1`) due to chlorophyll (`zbgr(2,ichl)`, [m-1])
+- $ex_{det}(k-1,1)$ is the attenuation rate of blue light (`b=1`) in the overlying grid cell (`k=k-1`) due to detritus (`ndet * dbgr(1)`, [m-1])
+- $ex_{CaCO_3}(k-1,1)$ is the attenuation rate of blue light (`b=1`) in the overlying grid cell (`k=k-1`) due to calcium carbonate (`carb * cbgr(1)`, [m-1])
+
+The irradiance in the red band (`b=3`) at the mid point of layer `k`, in contrast, is equal to
+
+$PAR^{mid}(k,3) = PAR^{mid}(k-1,3) * e^{(-0.5*(ex_{bgr}(k-1,3) * \Delta z(k-1) + ex_{bgr}(k,3) * \Delta z(k)))}$
+
+where
+
+- $PAR^{mid}(k-1,3)$ is the red light (`b=3`) at the mid-point of the overlying grid cell (`par_bgr_mid(k-1,3)`, [W m-2])
+- $ex_{bgr}(k-1,3)$ is the total attenuation of red light (`b=3`) in the overlying grid cell (`ek_bgr(k-1,3)`, [m-1])
+- $ex_{bgr}(k,3)$ is the total attenuation of red light (`b=3`) in the current grid cell (`ek_bgr(k,3)`, [m-1])
+- $\Delta z(k-1)$ and $\Delta z(k)$ are the grid cell thicknesses of the overlying and current grid cells (`dzt(i,j,k)`, [m])
+
+The total PAR available to phytoplantkon is assumed to be the sum of the blue, green and red bands. Because we assume that phytoplankton are
+homogenously distributed within a layer `k`, but we do not assume that light is homogenously distributed within that layer, we solve for the
+PAR that is seen by the average phytoplankton within that grid cell (`radbio`, $PAR$, [W m-2])
+
+$PAR(k) = \sum_{b=1}^3 \dfrac{(PAR^{top}(k,b) - PAR^{top}(k+1,b))}{ex_{bgr}(k,b) * \Delta z(k)}$
+
+where
+
+- $PAR^{top}(k,b)$ is the incoming photosynthetically active radiation at the top of grid cell `k` and light band `b` (`par_bgr_top(k,b)`, [W m-2])
+- $ex_{bgr}(k,b)$ is the attenuation rate of light band `b` in grid cell `k` (`ek_bgr(k,b)`, [m-1])
+- $\Delta z(k)$ is the grid cell thickness of grid cell `k` (`dzt(i,j,k)`, [m])
+
+This ensures phytoplankton growth in the model responds to the mean light they experience in the grid cell. See Eq. 19 from [Baird et al. (2020)](https://gmd.copernicus.org/articles/13/4503/2020/).
+
+The euphotic depth (`zeuphot(i,j)`, [m]) is defined as the depth where `radbio` falls below the 1% threshold of incidient shortwave radiation or below 0.01 W m-2, whichever is shallower.
+
+---
+
+
+### 2. Nutrient limitation of phytoplankton.
+
+At the start of each vertical loop the code computes biomass of phytoplankton (`biophy`, $B_{phy}^{C}$, [mmol C m-3]). Phytoplankton biomass is used to scale how both nitrate (`biono3`, $NO_{3}$, [mmol NO3 m-3]) and dissolved iron (`biofer`, $dFe$, [µmol Fe m-3]) affect the growth of phytoplankton. Using compilations of marine phytoplankton and zooplankton communities, [Wickman et al. (2024)](https://www.science.org/doi/10.1126/science.adk6901) show that the nutrient affinity, $\mathit{aff}$, of a phytoplankton cell is related to its volume, $V$, via
+
+$\mathit{aff} = V^{-0.57}$
+
+Additionally, the authors demonstrate that the volume of the average phytoplankton cell is related to the density (i.e., concentration) of phytoplankton via
+
+$V = (B_{phy}^{C})^{0.65}$
+
+when combining panels c and f of their Figure 1. This then relates the affinity of an average cell to the concentration of phytoplankton biomass as
+
+$\mathit{aff} = (B_{phy}^{C})^{-0.37}$
+
+With this information, we allow the half-saturation terms for nitrate (`phy_kni(i,j,k)`, $K_{phy}^{N}$, [mmol N m-3]) and dissolved iron (`phy_kfe(i,j,k)`, $K_{phy}^{Fe}$, [µmol Fe m-3]) uptake to vary as a function of phytoplankton biomass concentration. We set reference values for the half-saturation coefficient of nitrate (`phykn`, $K_{phy}^{N,0}$, [mmol N m-3]) and dissolved iron (`phykf`, $K_{phy}^{Fe,0}$, [µmol dFe m-3]) as input parameters to the model, and also set a threshold phytoplankton concentration (`phybiot`, $B_{phy}^{thresh}$, [mmol C m-3]) beneath which cell size cannot decrease and affinity can no longer increase. At this minimum, where affinity is maximised, the half-saturation coefficients are bounded to be 10% of their reference values.
+
+$K_{phy}^{N} = K_{phy}^{N,0} * \max(0.1, \max(0.0, (B_{phy}-B_{phy}^{thresh}))^{0.37} )$
+
+$K_{phy}^{Fe} = K_{phy}^{Fe,0} * \max(0.1, \max(0.0, (B_{phy}-B_{phy}^{thresh}))^{0.37} )$
+
+**Limitation of phytoplankton growth by nitrate** (`phy_lnit(i,j,k)`, $L_{phy}^{N}$), [dimensionless]) then follows the Monod equation:
+
+$L_{phy}^{N} = \dfrac{NO_3}{NO_3 + K_{phy}^{N}}$
+
+where
+
+- $NO_3$ is the ambient nitrate concentration (`biono3`, $NO_3$, [mmol N m-3])
+- $K_{phy}^{N}$ is the michaelis-menten half-saturation coefficient (`phy_kni(i,j,k)`, [mmol N m-3])
+
+**Limitation of phytoplankton growth by iron** follows an internal quota approach ([Droop, 1983](https://www.degruyterbrill.com/document/doi/10.1515/botm.1983.26.3.99/html)). Phytoplankton have a minimum iron quota (`phy_minqfe`, $Q_{phy}^{-Fe:C}$, [mol Fe (mol C)-1]) and an optimal quota for growth (`phy_optqfe`, $Q_{phy}^{*Fe:C}$, [mol Fe (mol C)-1]). The minimum iron quota, $Q_{phy}^{-Fe:C}$, is dependent on the chlorophyll content of the cell and the degree of nitrogen limitation according to
+
+$$
+\begin{align}
+Q_{phy}^{-Fe:C} = & 0.00167 / 55.85 * Q_{phy}^{Chl:C} * 12 \\
+& + 1.21 \times 10^{-5} * 14.0 / 55.85 / 7.625 * 0.5 * 1.5 * L_{phy}^{N} \\
+& + 1.15 \times 10^{-4} * 14.0 / 55.85 / 7.625 * 0.5 * L_{phy}^{N}
+\end{align}
+$$
+
+The first term reflects the amount of iron required for photosystems I and II. `0.00167/55.85` is equivalent to the grams of Fe per gram of chlorophyll divided by the grams of Fe per mol Fe, giving mol Fe per gram chlorophyll. This term is multipled by the chlorophyll to carbon ratio of the phytoplantkon cell (`phy_chlc`, $Q_{phy}^{Chl:C}$, [mol C (mol C)-1]) and grams of C per mol C, returning mol Fe per mol C. At a healthy chlorophyll:C ratio of 0.03, this term returns an Fe:C ratio of roughly 10 $\mu$mol:mol, which reproduces well known requirements of phytoplankton cells ([Morel, Rueter & Price, 1991](https://www.jstor.org/stable/43924569)). The second term, representing the respiratory iron requirement, is derived from [Flynn & Hipkin (1999)](https://onlinelibrary.wiley.com/doi/10.1046/j.1529-8817.1999.3561171.x) who estimated 1.21 $\times 10^{-5}$ grams Fe per gram N assimilated into the cell, which is converted to mol Fe per mol C with 14 g N per mol N divided by 55.85 g Fe per mol Fe $\times$ 7.625 mol C per mol N. This second term assumes that respiration is reduced as growth becomes more limited by available nitrogen (`phy_lnit(i,j,k)`, $L_{phy}^{N}$, [dimensionless]). Finally, the third term represents the iron required by nitrate/nitrite reduction. Nitrate assimilation requires roughly 1.8$\times$ more iron than ammonia assimilation ([Raven, 1988](https://nph.onlinelibrary.wiley.com/doi/abs/10.1111/j.1469-8137.1988.tb04196.x)). [Flynn & Hipkin (1999)](https://onlinelibrary.wiley.com/doi/10.1046/j.1529-8817.1999.3561171.x) estimated a demand of 1.15 $\times 10^{-4}$ g Fe per mol $NO_3$ reduced, which is accounted for by the nitrogen limitation term. Note that the 1.5 in the equation for $Q_{phy}^{-Fe:C}$ is designed to account for dark respiration (i.e., respiration when the cells are not growing) and the 0.5 refers to the fact that during cell division the cell must reinstate half of its Fe reserves.
+
+The Fe limitation factor (`phy_lfer(i,j,k)`, $L_{phy}^{Fe}$) is then computed from the present Fe:C quota of the phytoplankton cells (`phy_Fe2C`, $Q_{phy}^{Fe:C}$, [mol Fe (mol C)-1]) relative to the minimum and optimal quotas.
+
+$L_{phy}^{Fe} = \max\left(0.0, \min\left(1.0, \dfrac{ Q_{phy}^{Fe:C} - Q_{phy}^{-Fe:C} }{Q_{phy}^{*Fe:C}} \right)\right)$
+
+If the cell is Fe‑replete with a quota that exceeds the minimum quota by as much as the optimal quota, then Fe does not limit growth ($L_{phy}^{Fe}$ = 1). If the cell is Fe‑deplete with a quota equal to or less than the minimum quota, then the growth rate is reduced to zero. The optimal quota ($Q_{phy}^{*Fe:C}$) is therefore a measure of how much excess Fe is required to allow unrestricted growth.
+The minimum iron requirements of the cell for growth increases as the ratio of chlorophyll to carbon increases (`phy_chlc`, $Q_{phy}^{Chl:C}$, [mol C (mol C)-1]), as respiration increases and as the cell uses more nitrate. This formulation and the coefficients applied to chlorophyll content and nitrate use derive from [Flynn & Hipkin (1999)](https://onlinelibrary.wiley.com/doi/10.1046/j.1529-8817.1999.3561171.x).
+
+
+---
+
+
+### 3. Temperature‑dependent autotrophy and heterotrophy.
+
+**Autotrophy**
+
+The maximum potential growth rate for phytoplankton (`phy_mumax(i,j,k)`, $\mu_{phy}^{max}$, [s-1]) is prescribed by the temperature-dependent Eppley curve ([Eppley, 1972](https://spo.nmfs.noaa.gov/content/temperature-and-phytoplankton-growth-sea)). This formulation scales a reference growth rate at 0ºC via a power-law scaling with temperature (`Temp(i,j,k)`, $T$, [ºC]).
+
+$\mu_{phy}^{max} = \mu_{phy}^{0^{\circ}C} \cdot (β_{auto})^{T}$
+
+where
+
+- $\mu_{phy}^{0ºC}$ is the rate of phytoplankton growth at 0ºC (`abioa`, [s-1])
+- $β_{auto}$ is the base temperature-sensitivity coefficient for autotrophy (`bbioa`, [dimenionless])
+
+In the above, both $\mu_{phy}^{0ºC}$ and $β_{auto}$ are reference values input to the model and control how productive the ocean is and can therefore be modified to reflect known variations in the response of different phytoplankton types to temperature ([Anderson et al., 2021](https://www.nature.com/articles/s41467-021-26651-8)).
+
+**Heterotrophy**
+
+Heterotrophic processes include mortality of phytoplankton and zooplankton, grazing rates of zooplankton and the remineralisation rate of detritus in the water column and sediments. These processes are scaled similarly to autotrophy, where some reference rate at 0ºC ($\mu_{het}^{0^{\circ}C}$, [s-1]) is multiplied by a power-law with temperature ($β_{hete}$). Each heterotrophic process has a different $\mu_{het}^{0^{\circ}C}$ value and we expand on this later under mortality and grazing sections. However, the basic formulation for scaling heterotrophic metabolisms with temperature takes the form:
+
+$\mu_{het} = \mu_{het}^{0^{\circ}C} \cdot (β_{hete})^{T}$
+
+where
+
+- $\mu_{het}^{0ºC}$ is the rate of some heterotrophic metabolism at 0ºC ([s-1])
+- $β_{hete}$ is the base temperature-sensitivity coefficient for heterotrophy (`bbioh`, [dimenionless])
+
+See sections below for further details on heterotrophic metabolisms of mortality and grazing.
+
+---
+
+
+### 4. Light limitation of phytoplankton.
+
+Phytoplankton growth is limited by light through a photosynthesis–irradiance (P–I) relationship that links cellular chlorophyll content and available photosynthetically active radiation (`radbio`, $PAR$, [W m-2]).
+
+First, The initial slope of the P–I curve, (`phy_pisl`, $\alpha_{phy}$, [s-1 (W m-2)-1]), determines how efficiently phytoplankton convert light into carbon fixation. It is scaled by the cellular chlorophyll-to-carbon ratio.
+
+$\alpha_{phy} = \max(\alpha_{phy}^{Chl} \cdot Q_{phy}^{Chl:C} \ , \ \alpha_{phy}^{Chl} \cdot Q_{phy}^{-Chl:C}$ )
+
+where
+
+- $\alpha_{phy}^{Chl}$ is the photosynthetic efficiency per unit chlorophyll (`alphabio`, [s-1 (W m-2)-1 (mg/mg)-1])
+- $Q_{phy}^{Chl:C}$ is the in situ chlorophyll to carbon ratio (`phy_chlc`, [mol C (mol C)-1])
+- $Q_{phy}^{-Chl:C}$ is the minimum chlorophyll to carbon ratio of the cell (`phyminqc`, [mol C (mol C)-1]).
+
+This constraint prevents photosynthesis from collapsing unrealistically at low chlorophyll concentrations.
+
+Second, light limitation (`phy_lpar(i,j,k)`, $L_{phy}^{PAR}$), [dimensionless]) is calculated using an exponential P–I formulation.
+
+$L_{phy}^{PAR} = 1 - e^{- \alpha_{phy} PAR }$
+
+where
+
+- $PAR$ is the total photosynthetically available radiation (`radbio`, [W m-2])
+
+At low irradiance (PAR), growth increases approximately linearly with light, while at high irradiance photosynthesis asymptotically saturates. Photoinhibition is not included in this formulation.
+
+---
+
+
+### 5. Realized growth rate of phytoplankton.
+
+Realized growth of phytoplankton (`phy_mu(i,j,k)`, $\mu_{phy}$, [s-1]) is calculated as:
+
+$\mu_{phy} = \mu_{phy}^{max} L_{phy}^{PAR} \min(L_{phy}^{N}, L_{phy}^{Fe})$
+
+where:
+
+- $\mu_{phy}^{max}$ is the maximum potential rate of carbon fixation (`phy_mumax`, [s-1])
+- $L_{phy}^{PAR}$ is the growth limiter by light (`phy_lpar(i,j,k)`, [dimensionless])
+- $L_{phy}^{N}$ is the growth limiter by nitrogen (`phy_lnit(i,j,k)`, [dimensionless])
+- $L_{phy}^{Fe}$ is the growth limiter by iron (`phy_lfer(i,j,k)`, [dimensionless])
+
+We apply Liebig's law of the minimum ([Liebig, 1840](https://archive.org/details/organicchemistry00liebrich/mode/2up), [Blackman, 1905](https://doi.org/10.1093/oxfordjournals.aob.a089000)) to resources that are required for biomass synthesis (N and Fe), while light is considered multiplicative because it is an energy supply constraint that powers nutrient aquisition and biomass synthesis.
+
+The rate of carbon fixation by phytoplankton (`phygrow(i,j,k)`, $\mu_{phy}^{C}$, [mmol C m-3 s-1]) is
+
+$\mu_{phy}^{C} = \mu_{phy} B_{phy}^{C}$
+
+where
+
+- $B_{phy}^{C}$ is the carbon concentration of phytoplankton (`biophy`, [mmol C m-3])
+
+---
+
+
+### 6. Growth of chlorophyll.
+
+This step diagnoses the **rate of chlorophyll production** as a function of mixed-layer light, the phytoplankton growth rate and nutrient availability. The structure is consistent with the [Geider, MacIntyre & Kana (1997)](https://doi.org/10.3354/meps148187) formulation that relaxes the chlorophyll-to-carbon ratio towards an optimal value that supports photosynthetic growth under prevailing light and nutrient conditions.
+
+We first solve for the optimal chlorophyll-to-carbon ratio (`theta_opt`, $Q_{phy}^{*Chl:C}$, [mol C (mol C)-1]), which is diagnosed as the ratio required to support maximal photosynthetic carbon fixation under the ambient mean light level in the mixed layer, while accounting for nutrient limitation:
+
+$Q_{phy}^{*Chl:C} = \dfrac{Q_{phy}^{+Chl:C}}{1 + \dfrac{\alpha_{phy} PAR_{MLD} Q_{phy}^{+Chl:C}}{2 \mu_{phy}^{max} \min \left(L_{phy}^{N}, L_{phy}^{Fe} \right) }}$
+
+where:
+
+- $Q_{phy}^{+Chl:C}$ is the maximum allowable chlorophyll-to-carbon ratio (`phymaxqc`, [mol C (mol C)-1])
+- $\alpha_{phy}^{Chl}$ is the chlorophyll-specific initial slope of the P–I curve (`alphabio`, [s-1 (W m-2)-1 (mg/mg)-1])
+- $PAR_{MLD}$ is mean photosynthetically available radiation over the mixed layer (`radmld`, [W m-2>])
+- $\mu_{phy}^{max}$ is the temperature-dependent maximum phytoplankton growth rate (`phy_mumax`, [s-1])
+- $L_{phy}^{N}$ and $L_{phy}^{Fe}$ are the nutrient limitation factors for growth on N and Fe (`phy_lnit(i,j,k)`, [dimensionless]) (`phy_lfer(i,j,k)`, [dimensionless])
+
+We set a floor for the minimum chlorophyll-to-carbon ratio of phytoplankton via:
+
+$Q_{phy}^{*Chl:C} = \min \left( Q_{phy}^{*Chl:C}, Q_{phy}^{-Chl:C} \right)$
+
+where
+
+- $Q_{phy}^{-Chl:C}$ is the minimum allowable chlorophyll-to-carbon ratio (`phyminqc`, [mol C (mol C)-1])
+
+Growth of chlorophyll is then calculated as:
+
+$\dfrac{\delta Chl}{\delta t} = \mu_{phy} B_{phy}^{Chl} + \dfrac{ Q_{phy}^{*Chl:C} - Q_{phy}^{Chl:C} }{\tau_{phy}^{Chl}} \cdot B_{phy}^{C}$
+
+where:
+
+- $Q_{phy}^{Chl:C}$ is the in-situ chlorophyll-to-carbon ratio (`phy_chlc`, [mol C (mol C)-1])
+- $B_{phy}^{Chl}$ is the in-situ concentation of phytoplankton chlorophyll (`f_pchl(i,j,k)`, [mol C kg-1])
+- $\mu_{phy}$ is the realized growth rate of phytoplankton (`phy_mu(i,j,k)`, [s-1])
+- $\tau_{phy}^{Chl}$ is the rate of chlorophyll synthesis (`phytauqc`, [s-1])
+- $B_{phy}^{C}$ is the in-situ concentation of phytoplankton carbon (`f_phy(i,j,k)`, [mol C kg-1])
+
+This formulation elevates chlorophyll-to-carbon ratios in low light and supresses synthesis when nutrients are low. $\tau_{phy}^{Chl}$ is an input parameter at run time and should ideally be less than the doubling time of phytplankton given that phytoplankton can internally regulate their chlorophyll stores at rates greater than their overall growth.
+
+---
+
+
+### 7. Phytoplankton uptake of iron.
+
+Like chlorophyll, the iron content of phytoplankton is explicitly tracked as a tracer in WOMBAT-lite. First, a maximum Fe quota is found dependent on the maximum quota of Fe within the phytoplankton type and the phytoplankton concentration in the water column:
+
+$B_{phy}^{+Fe} = B_{phy}^{C} Q_{phy}^{+Fe:C}$
+
+where:
+
+- $B_{phy}^{+Fe}$ is the maximum Fe quota of the cell (`phy_maxqfe`, [mmol Fe m-3])
+- $Q_{phy}^{+Fe:C}$ is the maximum Fe:C ratio of the cell (`phymaxqf`, [mol Fe (mol C)-1])
+- $B_{phy}^{C}$ is the in-situ concentation of phytoplankton carbon (`f_phy(i,j,k)`, [mmol C m-3])
+
+Following [Aumont et al. (2015)](https://gmd.copernicus.org/articles/8/2465/2015/), this rate is scaled by three terms relating to (i) the michaelis-menten type affinity for dFe, (ii) up-regulation of dFe uptake representing investment in transporters when cell quotas are limiting to growth, and (iii) down regulation of dFe uptake associated with enriched cellular quotas:
+
+(i) $\dfrac{dFe}{dFe + K_{phy}^{Fe}}$
+
+(ii) $4 - \dfrac{4.5 L_{phy}^{Fe}}{0.5 + L_{phy}^{Fe}}$
+
+(iii) $\max\left(0, 1 - \dfrac{B_{phy}^{Fe} / B_{phy}^{+Fe}}{|1.05 - B_{phy}^{Fe} / B_{phy}^{+Fe}|} \right)$
+
+where:
+
+- $dFe$ is the in situ dissolved iron concentation (`biofer`, [µmol Fe m-3])
+- $K_{phy}^{Fe}$ is the half-saturation coefficient for dFe uptake by phytoplankton (`phy_kfe(i,j,k)`, [µmol Fe m-3])
+- $L_{phy}^{Fe}$ is the growth limiter by iron (`phy_lfer(i,j,k)`, [dimensionless])
+- $B_{phy}^{Fe}$ is the in situ Fe quota of the cell (`biophyfe`, [mmol Fe m-3])
+- $B_{phy}^{+Fe}$ is the maximum Fe quota of the cell (`phy_maxqfe`, [mmol Fe m-3])
+
+Note that we additionally include a fourth term that decreases the maximum dFe uptake of a cell under light limitation. This is informed by slower uptake of Fe by cells grown in darkness compared to those grown in light by roughly 10-fold ([Strzepek et al., 2025](https://doi.org/10.1093/ismejo/wraf015)), which may be due to physiological stimulation of Fe uptake machinery or photoreduction of ligand-bound iron complexes ([Kong et al., 2023](https://doi.org/10.1002/lno.12331); [Maldonado et al., 2005](https://doi.org/10.1029/2005GB002481)), or possibly a combination of both. To obtain a 10-fold relative increase in Fe uptake rates under light, we applied the following term:
+
+(iv) $\max\left(0.01, L_{phy}^{PAR}\right)^{0.5}$
+
+where:
+
+- $L_{phy}^{PAR}$ is the growth limiter by light (`phy_lpar(i,j,k)`, [dimensionless])
+
+Under very low light, this fourth term reduces maximum potential Fe uptake by 10-fold than what it otherwise would be. Dissolved iron uptake by phytoplankton is then calculated as:
+
+$\mu_{phy}^{Fe} = \mu_{phy}^{max} B_{phy}^{+Fe} \cdot (i) \cdot (ii) \cdot (iii) \cdot (iv)$
+
+where:
+
+- $\mu_{phy}^{Fe}$ is the realized uptake rate of dissolved iron by phytoplankton (`phy_dfeupt(i,j,k)`, [mmol Fe m-3])
+- $\mu_{phy}^{max}$ is the temperature-dependent maximum phytoplankton growth rate (`phy_mumax`, [s-1])
+
+
+---
+
+
+### 8. Iron chemistry (precipitation, scavenging and coagulation).
+
+Treatment of dissolved iron (`biofer`, $dFe$, [nmol Fe kg-1]) follows a combination of [Aumont et al. (2015)](https://gmd.copernicus.org/articles/8/2465/2015/) and [Tagliabue et al. (2023)](https://www.nature.com/articles/s41586-023-06210-5). Our calculations involve:
+1. Solving for the distinct pools of dissolved iron: free iron, ligand-bound iron and colloidal iron.
+2. Computing precipitation of free iron into nanoparticles that are permanently lost.
+3. Computing scavenging of free iron onto sinking organic particles.
+4. Computing coagulation of colloidal iron onto sinking organic particles.
+
+We first estimate the **solubility of free Fe from Fe3+** in solution using temperature, pH and salinity using the thermodynamic equilibrium equations of [Liu & Millero (2002)](https://www.sciencedirect.com/science/article/abs/pii/S0304420301000743).
+
+$T_K = \max(5.0, T) + 273.15$ \
+$T_K^{-1} = \dfrac{1}{T_K}$ \
+$I_{S} = \dfrac{19.924 S}{1000 - 1.005 S}$
+
+Solubility constants:
+
+$Fe_{sol1} = 10^{\left(-13.486 - 0.1856\sqrt{I_S} + 0.3073 I_S + 5254,T_K^{-1}\right)}$ \
+$Fe_{sol2} = 10^{\left(2.517 - 0.8885\sqrt{I_S} + 0.2139 I_S - 1320,T_K^{-1}\right)}$ \
+$Fe_{sol3} = 10^{\left(0.4511 - 0.3305\sqrt{I_S} - 1996,T_K^{-1}\right)}$ \
+$Fe_{sol4} = 10^{\left(-0.2965 - 0.7881\sqrt{I_S} - 4086,T_K^{-1}\right)}$ \
+$Fe_{sol5} = 10^{\left(4.4466 - 0.8505\sqrt{I_S} - 7980,T_K^{-1}\right)}$
+
+Final Fe(III) solubility:
+
+$Fe_{sol} = Fe_{sol1}\left([H^+]^3 + Fe_{sol2}[H^+]^2 + Fe_{sol3}[H^+] + Fe_{sol4} + \dfrac{Fe_{sol5}}{[H^+]}\right)\times10^{9}$
+
+where:
+
+- $T_K$ is in situ water temperature (`ztemk`, [ºK])
+- $I_{S}$ is a salinity coefficient (`zval`, [dimenionless])
+- $[H^+]$ is in situ hydrogen ion concentration (`hp`, [mol L-1])
+- $\times10^{9}$ converts [mol Fe kg-1] to [nmol Fe kg-1]
+- $Fe_{sol}$ is the final estimated solubility of dissolve iron in seawater (`fe3sol`, [nmol Fe kg-1]).
+
+Next we **estimate the concentration of colloidal iron** in solution following [Tagliabue et al. 2023](https://www.nature.com/articles/s41586-023-06210-5). Colloidal Fe (`fecol(i,j,k)`, $Fe_{col}$, [nmol Fe kg-1]) is whatever exceeds the inorganic solubility ceiling (`fe3sol`, $Fe_{sol}$, [nmol Fe kg-1]), but we enforce a hard minimum that colloids are at least 10% of total dissolved Fe (`biofer`, $dFe$, [nmol Fe kg-1]).
+
+$Fe_{col} = \max\left(0.1 dFe\, \ dFe - Fe_{sol} \right)$
+
+Following solving for colloidal Fe, we **partition the remaining dissolved Fe into ligand-bound and free iron**. To do so, we find the remaining dissolved iron not in colloidal form (`fe_sfe`, $Fe_{sFe}$, [nmol Fe kg-1]),
+
+$Fe_{sFe} = \max\left(0.0,\ dFe - Fe_{col} \right)$
+
+Partitioning is done using a standard quadratic form that drops out of mass balance + equilibrium for 1:1 complexation with a single ligand class. To do so, we need the temperature- and light-dependent conditional stability constant for Fe–ligand complexation (`fe_keq`, $Fe_{Keq}$, [nmol kg-1]). The temperature dependency comes from [Volker & Tagliabue (2015)](https://doi.org/10.1016/j.marchem.2014.11.008). Our light-dependency accounts for the photoreduction of photoreactive ligands, which was identified to reduce the conditional stability constant of aquachelin by 0.7 log10 units ([Barbeau et al., 2001](https://doi.org/10.1146/annurev.marine.010908.163712); [Vraspir & Butler, 2009](https://doi.org/10.1146/annurev.marine.010908.163712)):
+
+$Fe_{Keq} = \left( 10^{ \left(17.27 - 1565.7 T_K^{-1} \right) - 0.7\dfrac{PAR}{PAR + 10.0} } \right) \times 10^{-9}$,
+
+After finding $Fe_{Keq}$ we solve for the free Fe concentration (`feIII`, $Fe_{free}$, [nmol Fe kg-1]) using:
+
+$z = 1.0 + [Ligand] \cdot Fe_{Keq} - Fe_{sFe}\cdot Fe_{Keq}$ \
+$Fe_{free} = \dfrac{-z + \sqrt{z^2 + 4.0 Fe_{Keq} Fe_{sFe}}}{2 Fe_{Keq} + \varepsilon}$ \
+$Fe_{free} = \max\left(0,\ \min(Fe_{free},\ Fe_{sFe})\right)$
+
+where:
+
+- $[Ligand]$ is the concentration of ligand-bound dissolved iron (`felig`, [nmol Fe kg-1]).
+
+Whatever soluble Fe is not present as inorganic free iron is assigned to ligand-bound Fe:
+
+$Fe_{lig} = Fe_{sFe} - Fe_{free}$
+
+Now that we have separated the dissolved Fe pool into its subcomponents of free, ligand-bound and colloidal Fe all in units of [nmol Fe kg-1], we solve for precipiation of nanoparticles, scavenging and coagulation of dissolved Fe, all of which remove dFe from the water column. These are the major sinks outside of phytoplankton uptake.
+
+**Precipitation:**
+
+Precipitation of dissolved iron specifically affects free iron when the concentration is greater than that deemed soluble. This iron is permanently lost from the water column.
+
+$Fe_{nanop}^{→} = \max\left(0.0,\ Fe_{free} - Fe_{sol}\right)\, \gamma_{Fe}^{nano}$
+
+where
+
+- $Fe_{free}$ is the concentration of free dissolved iron (`feIII(i,j,k)`, [nmol Fe kg-1])
+- $Fe_{sol}$ is the final estimated solubility of dissolve iron in seawater (`fe3sol`, [nmol Fe kg-1])
+- $\gamma_{Fe}^{nano}$ is the rate constant of nanoparticle precipitation (`knano_dfe`, [s-1])
+- $Fe_{nanop}^{→}$ is the rate of loss of dFe via nanoparticle precipitation (`feprecip(i,j,k)`, [nmol Fe kg-1 s-1])
+
+**Scavenging:**
+
+Scavenging of dissolved iron specifically affects free iron, is accelerated by the presence of particles in the water column ([Ye et al., 2011](https://bg.copernicus.org/articles/8/2107/2011/); [Tagliabue et al., 2019](https://doi.org/10.1038/s41467-019-12775-5)) and we route this iron to sinking organic particles in WOMBAT-lite since we do not represent sinking authigenic particles.
+
+$Fe_{scav}^{→} = Fe_{free} \left(10^{-7} + \gamma_{Fe}^{scav} \cdot (B_{particles}^{M}) \right)$ \
+
+where
+
+- $Fe_{scav}^{→}$ is the rate of loss of dFe via scavenging onto all particulates (`fescaven(i,j,k)`, [nmol Fe kg-1 s-1])
+- $Fe_{free}$ is the concentration of free dissolved iron (`feIII(i,j,k)`, [nmol Fe kg-1])
+- $\gamma_{Fe}^{scav}$ is the rate constant of scavenging (`kscav_dfe`, [(mmol m-3)-1 s-1])
+- $B_{particles}^{M}$ is the in situ mass concentration of detrital particles in the water column (`partic`, [mmol m-3])
+
+The rate that scavenging adds to sinking organic detritus is (`fescadet(i,j,k)`, [nmol Fe kg-1 s-1]) is calculated as
+
+$Fe_{scav}^{→det} = Fe_{scav}^{→} \cdot \dfrac{ 2 \cdot B_{det}^{C} }{ B_{particles}^{M} }$\
+$B_{particles}^{M} = 2 \cdot B_{det}^{C} + 8.3 \cdot B_{CaCO_3}^{C}$
+
+where
+- $B_{det}^{C}$ is the concentration of particulate organic carbon (`biodet`, [mmol C m-3])
+- $B_{CaCO_3}^{C}$ is the concentration of particulate calcium carbonate in units of carbon (`biocaco3`, [mmol C m-3])
+
+Organic carbon-based particles ($B_{det}^{C}$) are multipled by 2 assuming that carbon represents half the mass of the particle, while inorganic carbon-based particles ($B_{CaCO_3}^{C}$) are multipled by 8.3 since the moleculate weight of calcium carbonate is 100 g mol-1.
+
+
+**Coagulation:**
+
+Similarly to scavenging of free iron, coagulation routes dissolved iron to sinking organic particles (since we do not represent sinking authigenic particles). However, coagulation acts on the colloidal fraction of dissolved iron ([Tagliabue et al., 2023)](https://www.nature.com/articles/s41586-023-06210-5)). Rates of coagulation of colloidal iron are of the form:
+
+$Fe_{coag}^{→ det} = Fe_{col} \gamma_{Fe}^{coag} \cdot S_{det}^{coag}$
+
+where
+- $Fe_{col}$ is the in situ concentration of colloidal iron (`fecol(i,j,k)`, [nmol Fe kg-1])
+- $\gamma_{Fe}^{coag}$ is the iron coagulation rate constant (`kcoag_dfe`, [(mmol m-3)-1 s-1])
+- $S_{det}^{coag}$ is the scaling coefficient to decelerate or accelerate coagulation (`zval`, [mmol C m-3])
+
+The coagulation scaling coefficient is dependent on the concentrations of dissolved organic carbon, particulate organic carbon, phytoplankton biomass and the rate of mixing via
+
+- $S_{det}^{coag} = H_{mix} (12 \cdot F_{coag} [DOC] + 9 \cdot B_{det}^{C}) + 2.5 \cdot B_{det}^{C} + 128 \cdot F_{coag} [DOC] + 725 \cdot B_{det}^{C} $
+
+where
+- $H_{mix}$ is a heavyside step function that is equal to 1 in the mixed layer and 0.01 beneath the mixed layer (`shear`, [dimensionless])
+- $F_{coag}$ is the phytoplankton-dependent coagulation factor
$F_{coag}= \max\left(\dfrac{1}{3},\ \dfrac{B_{phy}^{C}}{B_{phy}^{C} + 0.03}\right)$ (`biof`, [dimensionless])
+- $[DOC]$ is a proxy estimate of the concentration of dissolved organic carbon
$[DOC] = 10 + 40 \cdot\left(1 - \min(L_{phy}^{N},\ L_{phy}^{Fe})\right)$ (`biodoc`, [mmol C m-3])
+- $B_{phy}^{C}$ is the concentration of phytoplankton carbon biomass (`biophy`, [mmol C m-3])
+- $B_{det}^{C}$ is the concentration of detrital carbon biomass (`biodet`, [mmol C m-3])
+
+Together, these terms implement a biologically mediated coagulation pathway in which iron removal from the dissolved pool is tightly coupled to ecosystem state. The formulation reflects the central conclusion of [Tagliabue et al. (2023)](https://www.nature.com/articles/s41586-023-06210-5): that iron cycling is not governed solely by inorganic chemistry, but is strongly regulated by biological activity, organic matter dynamics, and particle ecology across the upper ocean. However, we note the absence of slowly sinking authigenic phases of iron in WOMBAT-lite that are represented by [Tagliabue et al. (2023)](https://www.nature.com/articles/s41586-023-06210-5). Instead, we note that colloidal coagulation passes dissolved iron directly to sinking organic detritus.
+
+---
+
+
+### 9. Mortality and remineralisation.
+
+Mortality of phytoplankton and zooplankton are affected by both linear ($\gamma$) and quadratic ($\Gamma$) terms. Linear terms are per-capita losses associated with the costs of basal metabolism. Quadratic, and thus density-dependent losses, are associated with disease, aggregation and coagulation, viruses, infection and canabalism. None of these processes are represented explicitly within the model, so we represent them implicitly.
+
+**Linear losses** for phytoplankton and zooplankton in [mmol C m-3 s-1] are modelled as
+
+$\gamma_{phy}^{C} = \gamma_{phy}^{0^{\circ}C} (β_{hete})^{T} B_{phy}^{C}$
+
+$\gamma_{zoo}^{C} = \gamma_{zoo}^{0^{\circ}C} (β_{hete})^{T} F_{zoo}^{\gamma} B_{zoo}^{C}$
+
+In the above, we scale down **linear zooplankton mortality** when zooplankton biomass is small, such that
+
+$F_{zoo}^{\gamma} = \dfrac{B_{zoo}^{C}}{B_{zoo}^{C} + K_{zoo}^{\gamma}}$,
+
+where
+
+- $\gamma_{phy}^{C}$ is the linear loss rate of phytoplankton biomass (`phymorl`, [mmol C m-3 s-1])
+- $\gamma_{zoo}^{C}$ is the linear loss rate of zooplankton biomass (`zoomorl`, [mmol C m-3 s-1])
+- $\gamma_{phy}^{0^{\circ}C}$ is the linear loss rate of phytoplankton biomass at 0ºC (`phylmor`, [s-1])
+- $\gamma_{zoo}^{0^{\circ}C}$ is the linear loss rate of zooplankton biomass at 0ºC (`zoolmor`, [s-1])
+- $β_{hete}$ is the base temperature-sensitivity coefficient for heterotrophy (`bbioh`, [dimenionless])
+- $T$ is the in situ temperature (`Temp(i,j,k)`, [ºC])
+- $B_{phy}^{C}$ is the concentration of phytoplankton carbon biomass (`biophy`, [mmol C m-3])
+- $B_{zoo}^{C}$ is the concentration of zooplankton carbon biomass (`biozoo`, [mmol C m-3])
+- $K_{zoo}^{\gamma}$ is the half-saturation coefficient for scaling down linear mortality losses (`zookz`, [mmolC m -3])
+
+
+**Quadratic losses** of phytoplankton and zooplankton in [mmol m-3 s-1] are modelled as
+
+$\Gamma_{phy}^{C} = \Gamma_{phy}^{0^{\circ}C} (β_{hete})^{T} \left(B_{phy}^{C}\right)^{2}$
+
+$\Gamma_{zoo}^{C} = \Gamma_{zoo}^{0^{\circ}C} (β_{hete})^{T} \left(B_{zoo}^{C}\right)^{2}$
+
+where
+
+- $\Gamma_{phy}^{C}$ is the quadratic (density-dependent) loss rate of phytoplankton biomass (`phymorq`, [mmol C m-3 s-1])
+- $\Gamma_{zoo}^{C}$ is the quadratic (density-dependent) loss rate of zooplankton biomass (`zoomorq`, [mmol C m-3 s-1])
+- $\Gamma_{phy}^{0^{\circ}C}$ is the quadratic (density-dependent) loss rate of phytoplankton biomass at 0ºC (`phyqmor`, [(mmol C m-3)-1 s-1])
+- $\Gamma_{zoo}^{0^{\circ}C}$ is the quadratic (density-dependent) loss rate of zooplankton biomass at 0ºC (`zooqmor`, [(mmol C m-3)-1 s-1])
+- $β_{hete}$ is the base temperature-sensitivity coefficient for heterotrophy (`bbioh`, [dimenionless])
+- $T$ is the in situ temperature (`Temp(i,j,k)`, [ºC])
+- $B_{phy}^{C}$ is the concentration of phytoplankton carbon biomass (`biophy`, [mmol C m-3])
+- $B_{zoo}^{C}$ is the concentration of zooplankton carbon biomass (`biozoo`, [mmol C m-3])
+
+
+**Remineralisation** of detritus is only affected by a quadratic, density-dependent loss term,
+
+$\Gamma_{det}^{C} = \Gamma_{det}^{0^{\circ}C} (β_{hete})^{T} \left(B_{det}^{C}\right)^{2}$
+
+- $\Gamma_{det}^{C}$ is the quadratic (density-dependent) loss rate of particulate organic detritus (`detremi`, [mmol C m-3 s-1])
+- $\Gamma_{det}^{0^{\circ}C}$ is the quadratic (density-dependent) loss rate of particulate organic detritus at 0ºC (`detlrem`, [(mmol C m-3)-1 s-1])
+- $β_{hete}$ is the base temperature-sensitivity coefficient for heterotrophy (`bbioh`, [dimenionless])
+- $T$ is the in situ temperature (`Temp(i,j,k)`, [ºC])
+- $B_{det}^{C}$ is the concentration of particulate organic detritus (`biodet`, [mmol C m-3])
+
+since hydrolyzation of organic detritus is performed by an heterotrophic bacterial population that is not explicitly resolved in the model and their acitivity is density-dependent.
+
+---
+
+
+### 10. Zooplankton grazing, egestion, excretion and assimilation.
+
+**Grazing by zooplankton** (`g_npz`, $g_{zoo}$, [s-1]) is computed using a Holling Type III functional response [Holling, 1959](https://doi.org/10.4039/Ent91385-7):
+
+$g_{zoo} = \dfrac{\mu_{zoo}^{max} (β_{hete})^{T} \varepsilon (B_{prey}^{C})^{2}}{\mu_{zoo}^{max} (β_{hete})^{T} + \varepsilon (B_{prey}^{C})^{2}}$
+
+where:
+
+- $\mu_{zoo}^{max}$ is the maximum rate of zooplankton grazing at 0ºC (`zoogmax`, [s-1])
+- $β_{hete}$ is the base temperature-sensitivity coefficient for heterotrophy (`bbioh`, [dimenionless])
+- $T$ is the in situ temperature (`Temp(i,j,k)`, [ºC])
+- $B_{prey}^{C}$ is the concentration of prey biomass (`zooprey`, [mmol C m-3])
+- $\varepsilon$ is the prey capture rate coefficient (`zooeps(i,j,k)`, [(mmol C m-3)-2])
+
+This formulation suppresses grazing at very low prey biomass ($B_{prey}^{C}$) due to reduced encounter and clearance rates, accelerates grazing at intermediate prey biomass as zooplankton effectively learn and switch to available prey, and saturates at high prey biomass due to handling-time limitation ([Gentleman and Neuheimer, 2008](https://doi.org/10.1093/plankt/fbn078); Rohr et al., [2022](https://doi.org/10.1016/j.pocean.2022.102878), [2024](https://doi.org/10.1029/2023GL107732)). This choice increases ecosystem stability and prolongs phytoplankton blooms relative to a Type II formulation.
+
+The application of $g_{zoo}^{max} (β_{hete})^{T}$ in both the numerator and denominator makes this grazing formula unique [(Rohr et al., 2023)](https://www.nature.com/articles/s43247-023-00871-w) and equivalent to a disk formulation, rather than a Michaelis–Menten formulation [(Rohr et al., 2022)](https://doi.org/10.1016/j.pocean.2022.102878). Practically, this amplifies grazing in warmer climes, but to a lesser extent than other formulations that apply the temperature amplification ($(β_{hete})^{T}$) only in the numerator [(Rohr et al., 2023)](https://www.nature.com/articles/s43247-023-00871-w). This dampens the effect that variations in temperature have on grazing activity, amplifying the effect of $\varepsilon$ and aligning with observations that the ratio of grazing to phytoplankton growth varies little between tropical and polar climes [(Calbet and Landry, 2004)](https://doi.org/10.4319/lo.2004.49.1.0051). Theoretically, this assumes some evolutionary adaptation to account for the physiological effects of temperature across environmental niches, such that the efficiency of prey capture and handling becomes more important to grazers than metabolic constraints due to temperature.
+
+The total prey biomass available to zooplankton is defined as a preference-weighted sum of phytoplankton and detritus that is normalized to reflect explicit dietary fractions ([Gentleman et al., 2003](https://doi.org/10.1016/j.dsr2.2003.07.001)):
+
+$B_{prey}^{C} = \phi_{zoo}^{phy} B_{phy}^{C} + \phi_{zoo}^{det} B_{det}^{C}$
+
+where:
+- $B_{phy}^{C}$ is the concentration of phytoplankton biomass (`biophy`, [mmol C m-3])
+- $B_{det}^{C}$ is the concentration of particulate organic detritus (`biodet`, [mmol C m-3])
+- $\phi_{zoo}^{phy}$ is the relative preference of zooplankton grazing on phytoplankton (`zooprefphy`, [dimensionless])
+- $\phi_{zoo}^{det}$ is the relative preference of zooplankton grazing on particulate detritus (`zooprefdet`, [dimensionless])
+
+and where:
+
+$\phi_{zoo}^{phy} + \phi_{zoo}^{det} = 1$
+
+The prey capture rate coefficient, $\varepsilon$ (`zooeps(i,j,k)`, $\varepsilon$, [(mmol C m-3)-2]), is allowed to vary as a function of prey biomass, following the prey-dependent behaviour described by [Rohr et al. (2024)](doi.org/10.1029/2023GL107732). This reflects a transition from microzooplankton-like feeding with higher prey capture rate coefficients at low prey biomass to mesozooplankton-like feeding with lower prey capture rate coefficients at high prey biomass.
+
+A prey-dependent scaling factor (`g_peffect`, $F_{prey}$, [dimensionless]) is defined as
+
+$F_{prey} = e^{\left(-B_{prey}^{C} \varepsilon_{shift} \right)}$
+
+and the effective capture rate coefficient, (`zooeps`, $\varepsilon$, [(mmol C m-3)-2]) is then computed as
+
+$\varepsilon = \varepsilon_{\min} + (\varepsilon_{\max} - \varepsilon_{\min}) F_{prey}$
+
+where:
+
+- $\varepsilon_{shift}$ is the rate at which the prey capture efficiency transitions from micro- to meso-zooplankton (`zooepsrat`, [(mmol C m-3)-1])
+- $\varepsilon_{\min}$ is the prey capture rate coefficient of a zooplankton community dominated by meso-zooplankton (`zooepsmin`, [(mmol C m-3)-2])
+- $\varepsilon_{\max}$ is the prey capture rate coefficient of a zooplankton community dominated by micro-zooplankton (`zooepsmax`, [(mmol C m-3)-2])
+
+At low prey biomass, $\varepsilon \rightarrow \varepsilon_{\max}$, enhancing grazing efficiency. At high prey biomass, $\varepsilon \rightarrow \varepsilon_{\min}$, reducing capture efficiency as handling time and feeding mode are more ineffective on average in a community with relatively more mesozooplankton.
+
+Zooplankton total grazing of biomass ([mmol C m-3 s-1]) is therefore
+
+$g_{zoo}^{C} = g_{zoo} B_{zoo}^{C} = \left(g_{zoo}^{→ phy} + g_{zoo}^{→ det}\right) B_{zoo}^{C}$
+
+where:
+
+- $g_{zoo}^{→ phy} = g_{zoo} \dfrac{\phi_{zoo}^{phy} B_{phy}^{C}}{B_{prey}^{C}}$ and is the proportion of zooplankton grazing of phytoplankton (`zoograzphy`, [mol C kg-1 s-1])
+- $g_{zoo}^{→ det} = g_{zoo} \dfrac{\phi_{zoo}^{det} B_{det}^{C}}{B_{prey}^{C}}$ and is the proportion of zooplankton grazing of detritus (`zoograzdet`, [mol C kg-1 s-1])
+
+
+**Zooplankton egestion, excretion and assimilation**, are then calculated assuming static assimilation coefficients. Grazed biomass is routed to either egestion or ingestion via an ingestion coefficient (`zooCingest`, $\lambda$, [mol C (mol C)-1]), with the egested fraction being equal to $1.0 - \lambda$. The biomass that is ingested is then split between assimilation and excretion based on an assimilation coefficient (`zooCassim`, $\eta$, [mol C (mol C)-1]) with the excreted fraction being equal to $1.0 - \eta$. Egestion ($E$), excretion ($X$) and assimilation ($A$) of organic carbon due to grazing of phytoplankton and detritus are:
+
+$E_{zoo}^{C,→ phy} = g_{zoo}^{→ phy} \left(1 - \lambda^{C}\right)$
+
+$E_{zoo}^{C,→ det} = g_{zoo}^{→ det} \left(1 - \lambda^{C}\right)$
+
+$X_{zoo}^{C,→ phy} = g_{zoo}^{→ phy} \lambda^{C} \left(1 - \eta^{C}\right)$
+
+$X_{zoo}^{C,→ det} = g_{zoo}^{→ det} \lambda^{C} \left(1 - \eta^{C}\right)$
+
+$A_{zoo}^{C,→ phy} = g_{zoo}^{→ phy} \lambda^{C} \eta^{C}$
+
+$A_{zoo}^{C,→ det} = g_{zoo}^{→ det} \lambda^{C} \eta^{C}$
+
+where:
+
+- $\lambda^{C}$ is the ingestion coefficient for carbon biomass by zooplankton (`zooCingest`, [mol C (mol C)-1])
+- $\eta^{C}$ is the assimilation coefficient for carbon biomass by zooplankton (`zooCassim`, [mol C (mol C)-1])
+- $g_{zoo}^{→ phy}$ is the zooplankton grazing rate of phytoplankton (`zoograzphy`, [mol C kg-1 s-1])
+- $g_{zoo}^{→ det}$ is the zooplankton grazing rate of detritus (`zoograzdet`, [mol C kg-1 s-1])
+- $E_{zoo}^{→ phy}$ is the egestion rate of carbon by zooplankton consuming phytoplankton (`zooegesphy(i,j,k)`, [mol C kg-1 s-1])
+- $E_{zoo}^{→ det}$ is the egestion rate of carbon by zooplankton consuming detritus (`zooegesdet(i,j,k)`, [mol C kg-1 s-1])
+- $X_{zoo}^{→ phy}$ is the excretion rate of carbon by zooplankton consuming phytoplankton (`zooexcrphy(i,j,k)`, [mol C kg-1 s-1])
+- $X_{zoo}^{→ det}$ is the excretion rate of carbon by zooplankton consuming detritus (`zooexcrdet(i,j,k)`, [mol C kg-1 s-1])
+- $A_{zoo}^{→ phy}$ is the assimilation rate of carbon by zooplankton consuming phytoplankton (`zooassiphy(i,j,k)`, [mol C kg-1 s-1])
+- $A_{zoo}^{→ det}$ is the assimilation rate of carbon by zooplankton consuming detritus (`zooassidet(i,j,k)`, [mol C kg-1 s-1])
+
+Because we track both carbon and iron through the ecosystem components, we assign unique ingestion and assimilation coefficients to carbon and iron. This separation of ingestion and assimilation coefficients for iron and carbon follows [Le Mézo & Galbraith (2021)](https://doi.org/10.1002/lno.11597). For iron, we apply unique ingestion (`zooFeingest`, $\lambda^{Fe}$, [mol Fe (mol Fe)-1]) and assimilation coefficients (`zooFeassim`, $\eta^{Fe}$, [mol Fe (mol Fe)-1]). [Le Mézo & Galbraith (2021)](https://doi.org/10.1002/lno.11597) show that if $\lambda^{Fe} << \lambda^{C}$ then egestion is enriched in Fe:C, and it follows that $\eta^{Fe} >> \eta^{C}$ so that zooplankton can absorb sufficient iron from their prey. Consequently:
+
+$E_{zoo}^{Fe,→ phy} = g_{zoo}^{→ phy} \left(1 - \lambda^{Fe}\right) \cdot \dfrac{B_{phy}^{Fe}}{B_{phy}^{C}}$
+
+$E_{zoo}^{Fe,→ det} = g_{zoo}^{→ det} \left(1 - \lambda^{Fe}\right) \cdot \dfrac{B_{det}^{Fe}}{B_{det}^{C}}$
+
+$X_{zoo}^{Fe,→ phy} = g_{zoo}^{→ phy} \lambda^{Fe} \left(1 - \eta^{Fe}\right) \cdot \dfrac{B_{phy}^{Fe}}{B_{phy}^{C}}$
+
+$X_{zoo}^{Fe,→ det} = g_{zoo}^{→ det} \lambda^{Fe} \left(1 - \eta^{Fe}\right) \cdot \dfrac{B_{det}^{Fe}}{B_{det}^{C}}$
+
+$A_{zoo}^{Fe,→ phy} = g_{zoo}^{→ phy} \lambda^{Fe} \eta^{Fe} \cdot \dfrac{B_{phy}^{Fe}}{B_{phy}^{C}}$
+
+$A_{zoo}^{Fe,→ det} = g_{zoo}^{→ det} \lambda^{Fe} \eta^{Fe} \cdot \dfrac{B_{det}^{Fe}}{B_{det}^{C}}$
+
+where:
+
+- $\lambda^{Fe}$ is the ingestion coefficient for iron biomass by zooplankton (`zooFeingest`, [mol Fe (mol Fe)-1])
+- $\eta^{Fe}$ is the assimilation coefficient for iron biomass by zooplankton (`zooFeassim`, [mol Fe (mol Fe)-1])
+- $g_{zoo}^{→ phy}$ is the zooplankton grazing rate of phytoplankton (`zoograzphy`, [mol C kg-1 s-1])
+- $g_{zoo}^{→ det}$ is the zooplankton grazing rate of detritus (`zoograzdet`, [mol C kg-1 s-1])
+- $\dfrac{B_{phy}^{Fe}}{B_{phy}^{C}}$ is the Fe:C ratio of in situ phytoplankton being grazed (`phy_Fe2C`, [mol Fe (mol C)-1])
+- $\dfrac{B_{det}^{Fe}}{B_{det}^{C}}$ is the Fe:C ratio of in situ detritus being grazed (`det_Fe2C`, [mol Fe (mol C)-1])
+- $E_{zoo}^{Fe,→ phy}$ is the egestion rate of iron by zooplankton consuming phytoplankton (`zooegesphyfe`, [mol Fe kg-1 s-1])
+- $E_{zoo}^{Fe,→ det}$ is the egestion rate of iron by zooplankton consuming detritus (`zooegesdetfe`, [mol Fe kg-1 s-1])
+- $X_{zoo}^{Fe,→ phy}$ is the excretion rate of iron by zooplankton consuming phytoplankton (`zooexcrphyfe`, [mol Fe kg-1 s-1])
+- $X_{zoo}^{Fe,→ det}$ is the excretion rate of iron by zooplankton consuming detritus (`zooexcrdetfe`, [mol Fe kg-1 s-1])
+- $A_{zoo}^{Fe,→ phy}$ is the assimilation rate of iron by zooplankton consuming phytoplankton (`zooassiphyfe`, [mol Fe kg-1 s-1])
+- $A_{zoo}^{Fe,→ det}$ is the assimilation rate of iron by zooplankton consuming detritus (`zooassidetfe`, [mol Fe kg-1 s-1])
+
+---
+
+
+### 11. CaCO3 calculations.
+
+**Dynamic $CaCO_3$ production and dissolution**
+
+When $CaCO_3$ dynamics are enabled (`do_caco3_dynamics = .true.`), the model computes both particulate inorganic carbon production (via the PIC:POC ratio) and $CaCO_3$ dissolution rates as functions of carbonate chemistry, temperature, and organic matter availability.
+
+
+**Production** of $CaCO_3$ in WOMBAT-lite comes from three sources: (1) density-dependent phytoplankton mortality ($P_{CaCO_3}^{\Gamma phy}$), (2) density-dependent zooplankton mortality ($P_{CaCO_3}^{\Gamma zoo}$) and (3) zooplankton egestion of grazed phytoplankton ($P_{CaCO_3}^{g_{zoo}^{→ phy}}$). Each term is multipled by the particulate inorganic to organic carbon production ratio (`pic2poc(i,j,k)`, $PIC:POC$, [mol C (mol C)-1]) to achieve a rate of $CaCO_3$ production [mol C kg-1 s-1]:
+
+(1) $P_{CaCO_3}^{\Gamma phy} = \Gamma_{phy}^{C} \cdot PIC:POC$
+
+(2) $P_{CaCO_3}^{\Gamma zoo} = \Gamma_{zoo}^{C} \cdot PIC:POC$
+
+(3) $P_{CaCO_3}^{g_{zoo}^{→ phy}} = g_{zoo}^{→ phy} \cdot PIC:POC \cdot \left(1 - F_{gut} \right)$
+
+where:
+
+- $\Gamma_{phy}^{C}$ is the quadratic (density-dependent) loss rate of phytoplankton biomass (`phymorq`, [mmol C m-3 s-1])
+- $\Gamma_{zoo}^{C}$ is the quadratic (density-dependent) loss rate of zooplankton biomass (`zoomorq`, [mmol C m-3 s-1])
+- $g_{zoo}^{→ phy}$ is the phytoplankton grazed by zooplankton (`zoograzphy`, [mol C kg-1 s-1])
+- $PIC:POC$ is the particulate inorganic (i.e., $CaCO_3$) to particulate organic ratio (`pic2poc(i,j,k)`, [mol C (mol C)-1])
+- $F_{gut}$ is the fraction of $CaCO_3$ that is dissolved within zooplankton guts (`fgutdiss`, [mol C (mol C)-1])
+
+In the above, the $PIC:POC$ ratio is formulated as:
+
+$PIC:POC = \min \left( 0.3, \left( f_{\text{inorg}} + 10^{-3 + 4.31 \times 10^{-6} \left( \dfrac{[HCO_3^-]}{[H^+]} \right)} \right) F_T \right)$
+
+where:
+
+- $f_{\text{inorg}}$ is the background PIC:POC ratio (`f_inorg`, [mol C (mol C)-1])
+- $[HCO_{3}^{-}]$ is the concentration of bicarbonate ions (`hco3`, [mol kg-1])
+- $[H^{+}]$ is concentration of free hydrogen ions (`htotal(i,j,k)`, [µmol kg-1])
+- $F_{T}$ is a temperature-dependent suppression term and if defined by $F_T = 0.55 + 0.45 \cdot \tanh \left(T - 4\right)$
+
+This formulation of $PIC:POC$ is therefore a function of the substrate–inhibitor ratio between bicarbonate and free hydrogen ions (`hco3 / htotal(i,j,k)`, $\dfrac{[HCO_3^-]}{[H^+]}$, [mol µmol-1]), following [Lehmann & Bach (2025)](https://www.nature.com/articles/s41561-025-01644-0). This reflects the sensitivity of calcification to carbonate system speciation, which is nonlinearly enhanced with increasing $[HCO_3^-]/[H^+]$. Moreover, the $F_{T}$ term strongly reduces $CaCO_3$ production in cold waters, enforcing near-zero calcification below approximately 3 °C consistent with observations of _Emiliania huxleyi_ growth limits in polar environments [(Fielding, 2013)](https://doi.org/10.4319/lo.2013.58.2.0663). Finally, we also cap the $PIC:POC$ ratio at an upper bound of 0.3 to prevent unrealistically high inorganic carbon production and accord with the highest measured ratios in the ocean.
+
+
+**Dissolution** of $CaCO_3$ is computed as the sum of four contributions:
+
+(i) the undersaturation-driven dissolution of calcite (`dissratcal`, $D_{\text{cal}}$, [s-1])
+
+(ii) the undersaturation-driven dissolution of aragonite (`dissratara`, $D_{\text{ara}}$, [s-1])
+
+(iii) the biologically mediated dissolution associated with degredation of detrital organic matter (`dissratpoc`, $D_{\text{det}}$, [s-1])
+
+(iv) the dissolution of $CaCO_3$ within zooplankton guts during consumption of detrital aggregates (`dissratzoo`, $D_{\text{gut}}$, [s-1])
+
+Total $CaCO_3$ dissolution is:
+
+$D_{\text{CaCO}3} = \left(D{\text{cal}} + D_{\text{ara}} + D_{\text{det}} \right) B_{CaCO_3}^{C} + D_{\text{gut}}$
+
+where:
+
+- $B_{CaCO_3}^{C}$ is the in situ concentration of $CaCO_3$ (`f_caco3(i,j,k)`, $B_{CaCO_3}$, [mol C kg-1]):
+
+The first three terms follow that suggested by [Kwon et al. (2024)](https://www.science.org/doi/full/10.1126/sciadv.adl0779).
+
+$D_{\text{cal}} = d_{\text{cal}} \max(0, 1 - \Omega_{\text{cal}})^{2.2}$
+
+$D_{\text{ara}} = d_{\text{ara}} \max(0,; 1 - \Omega_{\text{ara}})^{1.5}$
+
+$D_{\text{det}} = d_{\text{det}} \Gamma_{det}$
+
+where:
+
+- $\Omega_{\text{cal}}$ is the saturation state of calcite (`omega_cal(i,j,k)`, [dimenionless])
+- $\Omega_{\text{ara}}$ is the saturation state of aragonite (`omega_ara(i,j,k)`, [dimenionless])
+- $d_{\text{cal}}$ is the reference dissolution rate constant for calcite (`disscal`, [s-1])
+- $d_{\text{ara}}$ is the reference dissolution rate constant for aragonite (`dissara`, [s-1])
+- $\Gamma_{det}$ is the local remineralisation rate of organic matter (`reminr(i,j,k)`, [mmol C m-3 s-1])
+- $d_{\text{det}}$ is the reference dissolution rate constant per unit of organic detritus (`dissdet`, [(mmol C m-3)-1])
+
+For $D_{\text{cal}}$ and $D_{\text{ara}}$, dissolution is activated only under undersaturated conditions ($\Omega < 1$) and increases nonlinearly with increasing undersaturation. In contrast, $D_{\text{det}}$ represents shallow water $CaCO_3$ dissolution due to reducing microenvironments. In this scenario, $\Omega$ of calcite and aragonite in the water column tend to be > 1 [(Sulpis et al. (2021)](https://www.nature.com/articles/s41561-021-00743-y) but dissolution nonetheless occurs in microenvironments enriched in $CO_{2}^{*}$.
+
+The fourth term (`zoodiss`, $D_{\text{gut}}$, [mol C kg-1 s-1]), which represents dissolution of $CaCO_3$ during zooplankton consumption of particulates and is calculated as:
+
+$D_{\text{gut}} = g_{zoo}^{→ det} F_{gut} \dfrac{B_{CaCO_3}^{C}}{B_{det}^{C}}$
+
+where:
+
+- $g_{zoo}^{→ det}$ is the rate of zooplankton grazing of detritus (`zoograzdet`, [mol C kg-1 s-1])
+- $F_{gut}$ is the fraction of $CaCO_3$ that is dissolved within zooplankton guts (`fgutdiss`, [mol C (mol C)-1])
+- $\dfrac{B_{CaCO_3}^{C}}{B_{det}^{C}}$ is the in situ ratio of $CaCO_3$ to organic carbon detritus (`biocaco3/biodet`, [mol C (mol C)-1])
+
+Here we note that the processing of $CaCO_3$ by zooplankton grazing is treated differently to processing of organic carbon. For organic carbon, we route the biomass between zooplankton biomass (assimilation), inorganic nutrients (excretion) and particulate detritus (egestion). For $CaCO_3$ consumption by zooplankton the $CaCO_3$ is not assimilated since it does not contain nitrogen or other key elements for biosynthesis, and so is only routed between excretion to DIC and alkalinity or goes undissolved and remains $CaCO_3$ that sinks through the water column. This is supported by the fact that micro- and meso-zooplankton may dissolve 92±7% and 38-73% of coccolithophore calcite during feeding, respectively ([Smith et al., 2024](https://www.science.org/doi/10.1126/sciadv.adr5453); [White et al., 2018](https://www.nature.com/articles/s41598-018-28073-x); [Harris 1994](https://link.springer.com/article/10.1007/BF00347540)), and that the remainder is excreted and not assimilated ([Mayers et al., 2020](https://www.frontiersin.org/journals/marine-science/articles/10.3389/fmars.2020.569896/full)).
+
+
+**Static $CaCO_3$ production and dissolution**
+
+When $CaCO_3$ dynamics are disabled (`do_caco3_dynamics = .false.`), the model uses a static PIC:POC ratio (`f_inorg + 0.025`, [mol C (mol C)-1]) and $CaCO_3$ dissolution rate (`caco3lrem`, [s-1]). These are set as input parameters to the model.
+
+---
+
+
+### 12. Tracer tendencies.
+
+**Nitrate** (`f_no3(i,j,k)`, $NO_3$, [mol N kg-1])
+
+$\dfrac{\Delta NO_3}{\Delta t} = \left(\Gamma_{det}^{C} + \gamma_{zoo}^{C} + \gamma_{phy}^{C} + \eta_{zoo}^{C} - \mu_{phy}^{C} \right) \cdot R^{N:C}$
+
+**Oxygen** (`f_o2(i,j,k)`, $O_2$, [mol O2 kg-1])
+
+$\dfrac{\Delta O_2}{\Delta t} = \left( \mu_{phy}^{C} - \Gamma_{det}^{C} - \gamma_{zoo}^{C} - \gamma_{phy}^{C} - \eta_{zoo}^{C} \right) \cdot R^{O_2:C}$
+
+**Dissolved iron** (`f_fe(i,j,k)`, $dFe$, [mol Fe kg-1])
+
+$$
+\begin{align}
+\dfrac{\Delta dFe}{\Delta t} &= \Gamma_{det}^{C} Q_{det}^{Fe:C} + \gamma_{phy}^{C} Q_{phy}^{Fe:C} + \gamma_{zoo}^{C} Q_{zoo}^{Fe:C} +
+\left( g_{zoo}^{→ phy} \cdot Q_{phy}^{Fe:C} + g_{zoo}^{→ det} \cdot Q_{det}^{Fe:C} \right) \left(1 - \lambda \right) \cdot \eta \\
+&-
+\mu_{phy}^{Fe} - Fe_{nanop}^{→} - Fe_{scav}^{→} - Fe_{coag}^{→det}
+\end{align}
+$$
+
+**Phytoplankton** (`f_phy(i,j,k)`, $phy$, [mol C kg-1])
+
+$\dfrac{\Delta phy}{\Delta t} = \mu_{phy}^{C} - \Gamma_{phy}^{C} - \gamma_{phy}^{C} - g_{zoo}^{→ phy}$
+
+**Phytoplankton chlorophyll** (`f_pchl(i,j,k)`, $chl$, [mol Chl kg-1])
+
+$\dfrac{\Delta chl}{\Delta t} = \mu_{phy}^{chl} - \left( \Gamma_{phy}^{C} + \gamma_{phy}^{C} + g_{zoo}^{→ phy} \right) \cdot Q_{phy}^{Chl:C}$
+
+**Phytoplankton iron** (`f_phyfe(i,j,k)`, $phy^{Fe}$, [mol Fe kg-1])
+
+$\dfrac{\Delta phy^{Fe}}{\Delta t} = \mu_{phy}^{Fe} - \left( \Gamma_{phy}^{C} + \gamma_{phy}^{C} + g_{zoo}^{→ phy} \right) \cdot Q_{phy}^{Fe:C}$
+
+**Zoooplankton** (`f_zoo(i,j,k)`, $zoo$, [mol C kg-1])
+
+$\dfrac{\Delta zoo}{\Delta t} = \mu_{zoo}^{C} - \Gamma_{zoo}^{C} - \gamma_{zoo}^{C}$
+
+**Zoooplankton iron** (`f_zoofe(i,j,k)`, $zoo^{Fe}$, [mol Fe kg-1])
+
+$\dfrac{\Delta zoo^{Fe}}{\Delta t} = \left( g_{zoo}^{→ phy} \cdot Q_{phy}^{Fe:C} + g_{zoo}^{→ det} \cdot Q_{det}^{Fe:C} \right) \lambda - \left( \Gamma_{zoo}^{C} - \gamma_{zoo}^{C} \right) \cdot Q_{zoo}^{Fe:C}$
+
+**Detritus** (`f_det(i,j,k)`, $det$, [mol C kg-1])
+
+$\dfrac{\Delta det}{\Delta t} = \Gamma_{phy}^{C} + \Gamma_{zoo}^{C} + E_{zoo}^{C} - \left( \Gamma_{det}^{C} + g_{zoo}^{→ det} \right)$
+
+**Detritus iron** (`f_detfe(i,j,k)`, $det^{Fe}$, [mol Fe kg-1])
+
+$$
+\begin{align}
+\dfrac{\Delta det^{Fe}}{\Delta t} &= \Gamma_{phy}^{C} Q_{phy}^{Fe:C} + \Gamma_{zoo}^{C} Q_{zoo}^{Fe:C} +
+\left( g_{zoo}^{→ phy} \cdot Q_{phy}^{Fe:C} + g_{zoo}^{→ det} \cdot Q_{det}^{Fe:C} \right) \left(1 - \lambda \right) \left(1 - \eta \right) \\
+&+
+Fe_{scav}^{→det} + Fe_{coag}^{→det} -
+\left( \Gamma_{det}^{C} + g_{zoo}^{→ det} \right) \cdot Q_{det}^{Fe:C}
+\end{align}
+$$
+
+**Calcium Carbonate** (`f_caco3(i,j,k)`, $CaCO_3$, [mol C kg-1])
+
+$\dfrac{\Delta CaCO_3}{\Delta t} = \left( \Gamma_{zoo}^{C} + \Gamma_{phy}^{C} \right) PIC:POC + g_{zoo}^{→ phy} \left(1 - F_{zoo}^{diss} \right) PIC:POC - D_{\text{CaCO}3}$
+
+**Dissolved Inorganic Carbon** (`f_dic(i,j,k)`, $DIC$, [mol C kg-1])
+
+$\dfrac{\Delta DIC}{\Delta t} = \Gamma_{det}^{C} + \gamma_{zoo}^{C} + \gamma_{phy}^{C} + \eta_{zoo}^{C} - \mu_{phy}^{C} - \dfrac{\Delta CaCO_3}{\Delta t}$
+
+**Alkalinity** (`f_alk(i,j,k)`, $Alk$, [mol Eq kg-1])
+
+$\dfrac{\Delta Alk}{\Delta t} = - \dfrac{\Delta NO_3}{\Delta t} - 2 \cdot \dfrac{\Delta CaCO_3}{\Delta t}$
+
+
+---
+
+
+### 13. Check for conservation of mass.
+
+When checks for the conservation of mass is enabled (`do_check_n_conserve = .true.` or `do_check_c_conserve = .true.`), the model will calculate the budget of nitrogen or carbon both before and after the ecosystem equations have completed. This checks that the ecosystem equations detailed above have indeed conserved the mass of both nitrogen and carbon within the ocean. In WOMBATlite, both $NO_3$ and $DIC$ should be perfectly conserved during ecosystem cycling.
+
+---
+
+### 14. Additional operations on tracers.
+
+**First**, dissolved iron concentrations are set to equal 1 nM everywhere where the depth of the water column is less than 200 metres deep. WOMBAT-lite is not considered to be a model of the coastal ocean, but rather a model of the global pelagic ocean. Given that coastal waters are not limited in dissolved iron due to substantial interactions with sediments and exchange with the land, we universally set the dissolved iron concentration in these waters to 1 nM.
+
+**Second**, if dissolved iron concentrations dip below that measureable by operational detection limits, we reset these concentrations to this minimum (`zfermin`, $[dFe]^{min}$, [nmol Fe kg-1]):
+
+$[dFe]^{min} = \min\left( \max\left( 0.003 \cdot [NO_3]^{2}, 0.005 \right), 0.007 \right)$
+
+where:
+
+- $[NO_3]$ is the ambient nitrate concentration in units of [mmol m-3]
+
+---
+
+
+### 15. Sinking rate of particulates.
+
+WOMBAT-lite functions with a spatially variable sinking rate of both organic and inorganic (i.e., $CaCO_3$) particulate matter. The variable sinking rates of detritus and $CaCO_3$ are computed within the `generic_WOMBATlite_update_from_source` subroutine. Values are positive downward. Once computed, these 3D arrays are transfered to pointer arrays (`p_wdet(i,j,k)` and `p_wcaco3(i,j,k)`) that interface with a tridiagonal solver within the `g_tracer_vertdiff_G` subroutine.
+
+**Organic detritus**
+
+We consider that sinking rates of detritus are varied as a function of phytoplankton concentration (in a similar fashion to half-saturation coefficients described earlier), as well as the fraction of particulate matter that is $CaCO_3$. This approach is taken to emulate observations of varying sinking speeds [(Riley et al., 2012)](https://doi.org/10.1029/2011GB004085) and because such variations may be strongly dependent on phytoplankton community composition [(Bach et al., 2016)](https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016GB005372).
+
+In accordance with a more general Navier–Stokes drag equation and using a compilation of particle sinking speeds, [Cael et al. (2021)](https://doi.org/10.1029/2020GL091771) identified that the sinking velocity of detrital particles ($\omega_{det}$) in [m s-1]) is proportional to their diameter raised to the power of roughly 0.63, such that
+
+$\omega_{det} \propto d^{0.63}$.
+
+Knowing that
+
+$d = 6 V \pi \dfrac{1}{3}$
+
+and given that the average volume of phytoplankton cells can be approximated by $V = (B_{phy}^{C})^{0.65}$ [(Wickman et al., 2024)](https://doi.org/10.1126/science.adk6901), we can relate $\omega_{det}$ to the biomass concentration of phytoplankton multiplied by the scaler $\omega_{det}^{0}$
+
+$\omega_{det} = \omega_{det}^{0} \cdot \left( B_{phy}^{C,k=1} \right)^{0.21}$
+
+where:
+
+- $\omega_{det}^{0}$ is the sinking speed of organic detritus produced by a surface phytoplankton concentration of 1 mmol C m-3 (`wdetbio`, [m s-1])
+- $B_{phy}^{C}$ is the concentration of phytoplankton biomass (`biophy`, [mmol C m-3])
+
+This formula is identical to that presented by [Cael et al. (2021)](https://doi.org/10.1029/2020GL091771) in their Eq. (3), with the exception that we have related sinking rates to the biomass concentration of phytoplankton ($B_{phy}^{C}$) by assuming that $V = (B_{phy}^{C})^{0.65}$ based on marine phytoplankton data [(Wickman et al., 2024)](https://doi.org/10.1126/science.adk6901).
+
+As phytoplankton concentrations are negligible beneath the euphotic zone we use $B_{phy}^{C}$ only in the uppermost grid cell (k = 1). This assumes that the sinking velocities of marine aggregates can be related to phytoplankton community composition at the surface ([Bach et al., 2016](https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016GB005372); [Iversen and Lampitt, 2020](https://doi.org/10.1016/j.pocean.2020.102445)), which varies more horizontally across the ocean than vertically. Moreover, because we do not include dissolved/suspended organic matter as a tracer in WOMBAT-lite, we must also account for the large fraction of organics that are suspended and thus neutrally buoyant in the gyres. As such, we include a phytoplankton biomass threshold (`phybiot`, $B_{phy}^{thresh}$, [mmol C m-3]) above which sinking accelerates and beneath which any produced detritus emulates dissolved (neutrally buoyant) organic matter:
+
+$\omega_{det} = \omega_{det}^{0} \cdot \max \left( 0.0, B_{phy}^{C,k=1} - B_{phy}^{thresh} \right)^{0.21}$
+
+Sinking speeds are then accelerated depending on the fraction of $CaCO_3$ within the particulate mass by
+
+$\omega_{det} = \omega_{det} + \dfrac{10}{86400} \cdot \min \left( 1, \dfrac{B_{CaCO_3}^{C}}{B_{CaCO_3}^{C} + B_{det}^{C}} \right)$
+
+where:
+
+- $B_{phy}^{C,k=1}$ is the concentration of phytoplankton biomass at the surface (`biophy1`, [mmol C m-3])
+- $B_{phy}^{thresh}$ is the phytoplankton biomass threshold above which the community cell size begins to increase (`phybiot`, [mmol C m-3])
+- $B_{CaCO_3}^{C}$ is the in situ concentration of $CaCO_3$ biomass (`f_caco3(i,j,k)`, [mol C kg-1])
+- $B_{det}^{C}$ is the in situ concentration of organic detrital biomass (`f_det(i,j,k)`, [mol C kg-1])
+
+
+and where we chose a maximum increase of 10 metres per day based on the more modest effect observed in mesocosm experiments ([Bach et al., 2016](https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016GB005372).
+
+Finally, we apply a linear increase to sinking speeds with depth to ensure that the trend in the concentration of detritus with depth exhibits a power law behavior, which is widely observed ([Berelson, 2001](https://doi.org/10.1016/S0967-0645(01)00102-3); [Martin et al., 1987](https://doi.org/10.1016/0198-0149(87)90086-0)), thought to be associated with a greater attenuation of more slowly sinking particles, and shows better performance than a constant sinking rate in models ([Tjiputra et al., 2020](https://gmd.copernicus.org/articles/13/2393/2020/)). This is applied after the previous equation as:
+
+$\omega_{det} = \omega_{det} + \max\left(0.0, \dfrac{z}{5000} \cdot \left(\omega_{det}^{max} - \omega_{det} \right)\right)$
+
+where:
+
+- $z$ is the in situ depth (`zw(i,j,k)`, [m])
+- $\omega_{det}^{max}$ is the terminal sinking rate of organic detritus (`wdetmax`, [m s-1])
+
+
+**Calcium carbonate**
+
+The sinking rate of $CaCO_3$ is considered to always be a fraction of the sinking rate of organic detritus, such that
+
+$\omega_{CaCO_3} = \omega_{det} \cdot \dfrac{\omega_{CaCO_3}^{0}}{\omega_{det}^{0}}$,
+
+where:
+
+- $\omega_{det}^{0}$ is the sinking speed of organic detritus produced by a surface phytoplankton concentration of 1 mmol C m-3 (`wdetbio`, [m s-1])
+- $\omega_{CaCO_3}^{0}$ is the sinking speed of $CaCO_3$ produced by a surface phytoplankton concentration of 1 mmol C m-3 (`wcaco3`, [m s-1])
+
+Although more dense than organic matter, $CaCO_3$ particles tend to be smaller than organic aggregates and sink at a slower rate ([De La Rocha & Passow, 2007](https://doi.org/10.1016/j.dsr2.2007.01.004); [Zhang et al., 2018](https://doi.org/10.5194/bg-15-4759-2018)). Furthermore, the shedding of coccoliths by coccolithophores, which are near-neutrally bouyant, also contributes to a slower mean sinking speed of $CaCO_3$ ([Balch et al., 2009](https://doi.org/10.1029/2008JC004902)). Hence, $\omega_{CaCO_3}^{0} should be < \omega_{det}^{0}$.
+
+The degree to which $CaCO_3$ particles sink more slowly than organic detritus is controlled by the user through altering the ratio of $\omega_{CaCO_3}^{0}$ to $\omega_{det}^{0}$.
+
+
+**Detrital iron**
+
+WOMBAT-lite sinks organic detrital iron at the same rate as organic detrital carbon.
+
+---
+
+
+### 16. Sedimentary processes.
+
+WOMBAT-lite tracks the accumulation of organic detrital carbon (`p_det_sediment(i,j)`, $B_{det,sed}^{C}$, [mol m-2]), organic detrital iron (`p_detfe_sediment(i,j)`, $B_{det,sed}^{Fe}$, [mol m-2]) and $CaCO_3$ (`p_caco3_sediment(i,j)`, $B_{CaCO_3,sed}^{C}$, [mol m-2]) within sedimentary pools. The organic pools contribute to bottom fluxes of dissolved inorganic carbon (DIC), nitrate ($NO_3$), dissolved iron (dFe), oxygen ($O_2$) and alkalinity (Alk). Remineralisation of organic carbon ($\gamma_{det,sed}^{C}$) produces DIC and $NO_3$, but removes $O_2$ and Alk. Ratios of nitrogen to carbon and oxygen to carbon are static at 16:122 and -172:122. Remineralisation of organic iron produces dFe.
+
+$\gamma_{det,sed}^{C} = \gamma_{det,sed}^{0^{\circ}C} (β_{hete})^{T} B_{det,sed}^{C}$
+
+$\gamma_{det,sed}^{N} = \gamma_{det,sed}^{C} R^{N:C}$
+
+$\gamma_{det,sed}^{O_2} = \gamma_{det,sed}^{C} R^{O_2:C}$
+
+$\gamma_{det,sed}^{Alk} = -\gamma_{det,sed}^{N}$
+
+$\gamma_{det,sed}^{Fe} = \gamma_{det,sed}^{0^{\circ}C} (β_{hete})^{T} B_{det,sed}^{Fe}$
+
+where:
+
+- $\gamma_{det,sed}^{C}$ is the remineralisation rate of organic detrital carbon in the sediment pool (`det_sed_remin(i,j)`, [mol C m-2 s-1])
+- $\gamma_{det,sed}^{N}$ is the production of nitrate due to organic detrital remineralisation in the sediment pool ([mol N m-2 s-1])
+- $\gamma_{det,sed}^{O_2}$ is the consumption rate of oxygen during remineralisation of organics in the sediment pool ([mol $O_2$ m-2 s-1])
+- $\gamma_{det,sed}^{Alk}$ is the production of alkalinity due to organic detrital remineralisation in the sediment pool ([mol Alk m-2 s-1])
+- $\gamma_{det,sed}^{Fe}$ is the production of dissolved iron due to organic detrital remineralisation in the sediment pool ([mol Fe m-2 s-1])
+- $β_{hete}$ is the base temperature-sensitivity coefficient for heterotrophy (`bbioh`, [dimenionless])
+- $T$ is the in situ temperature at the sediment-water interface (`sedtemp(i,j)`, [ºC])
+- $B_{det,sed}^{C}$ is the total concentration of organic detrital carbon biomass (`p_det_sediment(i,j,1`), [mol C m-2])
+- $B_{det,sed}^{Fe}$ is the total concentration of organic detrital iron biomass (`p_detfe_sediment(i,j,1`), [mol Fe m-2])
+- $R^{N:C}$ is the static ratio of nitrogen to carbon in organic matter and is equal to 16:122
+- $R^{O_2:C}$ is the static ratio of dissolved oxygen utilisation to carbon content in organic matter and is equal to -172:122
+
+Alk and DIC are in turn affected by dissolution of the sedimenary $CaCO_3$ pool, which is considered as entirely calcite. Sedimentary dissolution is controlled by bottom water temperature and an estimate of the pore-water calcite saturation state ($\Omega_{cal}^{sed}$):
+
+$D_{cal,sed} = d_{cal,sed} (β_{hete})^{T} \max(1 - \Omega_{cal,+sed}, 1 - \Omega_{cal,sed})^{4.5}$
+
+where
+
+- $d_{cal,sed}$ is a base rate of dissolution (`caco3lrem_sed`, [s-1])
+- $\Omega_{cal,+sed}$ is the maximum possible saturation state within sediment pore water (`omegamax_sed`, [dimensionless])
+
+The $\Omega_{cal,sed}$ is calculated using the [`mocsy`](https://github.com/ACCESS-NRI/mocsy) package for solving carbonate chemistry of seawater ([Orr & Epitalon, 2015](https://doi.org/10.5194/gmd-8-485-2015)). These routines require Alk and DIC as inputs, along with nutrient concentrations and temperature and salinity of bottom waters. For DIC, we chose to sum the water column concentration of DIC and the organic carbon content of the sediment to approximate the interstitial (i.e., porewater) DIC concentration. We assume that the organic carbon content of the sediment (`p_det_sediment` [mol m-2]) is relevant over one meter, and therefore can be automatically converted to [mol m-3]. With this assumption these arrays can be added together and approximates the reducing conditions of organic-rich sediments, which have lower $\Omega_{cal,sed}$. This ensures a greater rate of $CaCO_3$ dissolution within the sediment as organic matter accumulates.
+
+Overall bottom fluxes of tracers are:
+
+$\dfrac{\Delta NO_3}{\Delta t} = \gamma_{det,sed}^{N}$
+
+$\dfrac{\Delta O_2}{\Delta t} = \gamma_{det,sed}^{O_2}$
+
+$\dfrac{\Delta dFe}{\Delta t} = \gamma_{det,sed}^{Fe}$
+
+$\dfrac{\Delta DIC}{\Delta t} = \gamma_{det,sed}^{C} + D_{cal,sed}$
+
+$\dfrac{\Delta Alk}{\Delta t} = \gamma_{det,sed}^{Alk} + 2 \cdot D_{cal,sed}$
+
+where:
+
+- $\dfrac{\Delta NO_3}{\Delta t}$ is the total flux of nitrate into the water column (`b_no3(i,j)`, [mol N m-2 s-1])
+- $\dfrac{\Delta O_2}{\Delta t}$ is the total flux of oxygen into the water column (`b_o2(i,j)`, [mol $O_2$ m-2 s-1])
+- $\dfrac{\Delta dFe}{\Delta t}$ is the total flux of dissolved iron into the water column (`b_fe(i,j)`, [mol Fe m-2 s-1])
+- $\dfrac{\Delta DIC}{\Delta t}$ is the total flux of dissolved inorganic carbon into the water column (`b_dic(i,j)`, [mol C m-2/s-1])
+- $\dfrac{\Delta Alk}{\Delta t}$ is the total flux of alkalinity into the water column (`b_alk(i,j)`, [mol Alk m-2 s-1])
+
+---
+
+## Subroutine - "update_from_bottom"
+
+The subroutine `generic_WOMBATlite_update_from_bottom` moves sinking organic material from the water column into the sediment pools.
+It is at this point that the model performs permanent burial of sinking organic matter if desired.
+
+---
+
+### Permanent burial of particulates.
+
+If `do_burial = .true.`, we compute the fraction of incident sinking organic matter and $CaCO_3$ that is permanently buried in the sediments. This permanently buried fraction is effectively removed from the model and therefore is not accumulated within the sedimentary pools.
+
+The fraction buried is calculated according to Equation 3 of [Dunne et al. (2007)](https://doi.org/10.1029/2006GB002907):
+
+$F_{bury} = 0.013 \cdot 0.53 \dfrac{(f_{org})^{2}}{\left(7 + f_{org}\right)^{2}}$
+
+where $f_{org}$ is the rain rate of organic carbon detritus on the seafloor in [mmol C m-2 s-1].
+
+As organic matter rains down at a more rapid rate, the fraction of incident organic carbon, organic iron and $CACO_3$ that is buried increases.
+
+If `do_conserve_tracers = .true.`, then we capture the total loss of both Alk and $NO_3$ via burial or organic detritus and $CaCO_3$ and redistribute the Alk and $NO_3$ amount back at the ocean surface. This amount of each tracer is redistributed uniformly to avoid strong gradients.
+
+---
diff --git a/documentation/docs/pages/index.md b/documentation/docs/pages/index.md
index d1d04561..76f8849b 100644
--- a/documentation/docs/pages/index.md
+++ b/documentation/docs/pages/index.md
@@ -1,3 +1,7 @@
-# WOMBAT documentation
+# Home
-Welcome to the documentation for the World Ocean Model of Biogeochemistry and Trophic dynamics (WOMBAT).
\ No newline at end of file
+Welcome to the documentation for the World Ocean Model of Biogeochemistry and Trophic dynamics (WOMBAT).
+
+You can navigate this documentation using the pane on the left. Some pages that may be of interest are:
+
+- [WOMBATlite model description](Model_description/WOMBATlite_model_description/)
diff --git a/documentation/mkdocs.yml b/documentation/mkdocs.yml
index 19480ac2..48f1fb86 100755
--- a/documentation/mkdocs.yml
+++ b/documentation/mkdocs.yml
@@ -95,9 +95,9 @@ markdown_extensions:
# Navigation
nav:
- Home: pages/index.md
- # List other pages here
- # Sub-pages can be listed as nested items
- # Documentation reference (https://www.mkdocs.org/user-guide/configuration/#nav)
+
+ - Model description:
+ - WOMBATlite: pages/Model_description/WOMBATlite_model_description.md
# Footer
extra:
diff --git a/generic_tracers/generic_WOMBATlite.F90 b/generic_tracers/generic_WOMBATlite.F90
index ebb5d74b..330f814c 100644
--- a/generic_tracers/generic_WOMBATlite.F90
+++ b/generic_tracers/generic_WOMBATlite.F90
@@ -160,8 +160,10 @@ module generic_WOMBATlite
phymaxqf, &
phylmor, &
phyqmor, &
- zooassi, &
- zooexcr, &
+ zooCingest, &
+ zooCassim, &
+ zooFeingest, &
+ zooFeassim, &
fgutdiss, &
zookz, &
zoogmax, &
@@ -187,7 +189,6 @@ module generic_WOMBATlite
dissara, &
dissdet, &
ligand, &
- fcolloid, &
knano_dfe, &
kscav_dfe, &
kcoag_dfe, &
@@ -301,23 +302,26 @@ module generic_WOMBATlite
phy_feupreg, &
phy_fedoreg, &
phygrow, &
- phyresp, &
- phymort, &
+ phymorl, &
+ phymorq, &
zooeps, &
zoograzphy, &
zoograzdet, &
- zooresp, &
- zoomort, &
+ zoomorl, &
+ zoomorq, &
zooexcrphy, &
zooexcrdet, &
- zooslopphy, &
- zooslopdet, &
+ zooassiphy, &
+ zooassidet, &
+ zooegesphy, &
+ zooegesdet, &
reminr, &
detremi, &
pic2poc, &
dissratcal, &
dissratara, &
dissratpoc, &
+ zoodiss, &
caldiss, &
aradiss, &
pocdiss, &
@@ -378,23 +382,26 @@ module generic_WOMBATlite
id_phy_feupreg = -1, &
id_phy_fedoreg = -1, &
id_phygrow = -1, &
- id_phyresp = -1, &
- id_phymort = -1, &
+ id_phymorl = -1, &
+ id_phymorq = -1, &
id_zooeps = -1, &
id_zoograzphy = -1, &
id_zoograzdet = -1, &
- id_zooresp = -1, &
- id_zoomort = -1, &
+ id_zoomorl = -1, &
+ id_zoomorq = -1, &
id_zooexcrphy = -1, &
id_zooexcrdet = -1, &
- id_zooslopphy = -1, &
- id_zooslopdet = -1, &
+ id_zooassiphy = -1, &
+ id_zooassidet = -1, &
+ id_zooegesphy = -1, &
+ id_zooegesdet = -1, &
id_reminr = -1, &
id_detremi = -1, &
id_pic2poc = -1, &
id_dissratcal = -1, &
id_dissratara = -1, &
id_dissratpoc = -1, &
+ id_zoodiss = -1, &
id_caldiss = -1, &
id_aradiss = -1, &
id_pocdiss = -1, &
@@ -865,13 +872,13 @@ subroutine generic_WOMBATlite_register_diag(diag_list)
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
- 'phyresp', 'Respiration of phytoplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
- wombat%id_phyresp = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ 'phymorl', 'Linear mortality of phytoplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_phymorl = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
- 'phymort', 'Mortality of phytoplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
- wombat%id_phymort = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ 'phymorq', 'Quadratic mortality of phytoplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_phymorq = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
@@ -890,13 +897,13 @@ subroutine generic_WOMBATlite_register_diag(diag_list)
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
- 'zooresp', 'Respiration of zooplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
- wombat%id_zooresp = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ 'zoomorl', 'Linear mortality of zooplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_zoomorl = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
- 'zoomort', 'Mortality of zooplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
- wombat%id_zoomort = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ 'zoomorq', 'Quadratic mortality of zooplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_zoomorq = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
@@ -910,13 +917,23 @@ subroutine generic_WOMBATlite_register_diag(diag_list)
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
- 'zooslopphy', 'Sloppy feeding of zooplankton on phytoplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
- wombat%id_zooslopphy = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ 'zooassiphy', 'Assimilation into biomass of zooplankton feeding on phytoplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_zooassiphy = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
- 'zooslopdet', 'Sloppy feeding of zooplankton on detritus', 'h', 'L', 's', 'molC/kg/s', 'f')
- wombat%id_zooslopdet = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ 'zooassidet', 'Assimilation into biomass of zooplankton feeding on detritus', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_zooassidet = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
+
+ vardesc_temp = vardesc( &
+ 'zooegesphy', 'Egestion of zooplankton feeding on phytoplankton', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_zooegesphy = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
+
+ vardesc_temp = vardesc( &
+ 'zooegesdet', 'Egestion of zooplankton feeding on detritus', 'h', 'L', 's', 'molC/kg/s', 'f')
+ wombat%id_zooegesdet = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
vardesc_temp = vardesc( &
@@ -949,6 +966,11 @@ subroutine generic_WOMBATlite_register_diag(diag_list)
wombat%id_dissratpoc = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
+ vardesc_temp = vardesc( &
+ 'zoodiss', 'Dissolution of CaCO3 due to zooplankton grazing', 'h', 'L', 's', 'molCaCO3/kg/s', 'f')
+ wombat%id_zoodiss = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
+ init_time, vardesc_temp%longname, vardesc_temp%units, missing_value=missing_value1)
+
vardesc_temp = vardesc( &
'caldiss', 'Dissolution of Calcite CaCO3', 'h', 'L', 's', 'molCaCO3/kg/s', 'f')
wombat%id_caldiss = register_diag_field(package_name, vardesc_temp%name, axes(1:3), &
@@ -1200,13 +1222,21 @@ subroutine user_add_params
!-----------------------------------------------------------------------
call g_tracer_add_param('phybiot', wombat%phybiot, 0.6)
- ! Zooplankton assimilation efficiency [1]
+ ! Zooplankton ingestion efficiency of prey carbon (the rest is egested) [1]
!-----------------------------------------------------------------------
- call g_tracer_add_param('zooassi', wombat%zooassi, 0.3)
+ call g_tracer_add_param('zooCingest', wombat%zooCingest, 0.8)
- ! Zooplankton excretion of unassimilated prey [0-1]
+ ! Zooplankton assimilation of ingested prey carbon (the rest is excreted) [0-1]
!-----------------------------------------------------------------------
- call g_tracer_add_param('zooexcr', wombat%zooexcr, 0.75)
+ call g_tracer_add_param('zooCassim', wombat%zooCassim, 0.3)
+
+ ! Zooplankton ingestion efficiency of prey carbon (the rest is egested) [1]
+ !-----------------------------------------------------------------------
+ call g_tracer_add_param('zooFeingest', wombat%zooFeingest, 0.2)
+
+ ! Zooplankton assimilation of ingested prey carbon (the rest is excreted) [0-1]
+ !-----------------------------------------------------------------------
+ call g_tracer_add_param('zooFeassim', wombat%zooFeassim, 0.9)
! Zooplankton dissolution efficiency of CaCO3 within guts [1]
!-----------------------------------------------------------------------
@@ -1307,17 +1337,19 @@ subroutine user_add_params
!-----------------------------------------------------------------------
call g_tracer_add_param('ligand', wombat%ligand, 0.5)
- ! Fraction of dissolved iron in colloidal form [0-1]
- !-----------------------------------------------------------------------
- call g_tracer_add_param('fcolloid', wombat%fcolloid, 0.5)
-
! Precipitation of Fe` as nanoparticles (in excess of solubility) [/d]
!-----------------------------------------------------------------------
call g_tracer_add_param('knano_dfe', wombat%knano_dfe, 0.1)
- ! Scavenging of Fe` onto biogenic particles [(mmolC/m3)-1 d-1]
+ ! Scavenging of Fe` onto biogenic particles [(mmol/m3)-1 d-1]
!-----------------------------------------------------------------------
- call g_tracer_add_param('kscav_dfe', wombat%kscav_dfe, 0.05)
+ ! Ye et al., 2011 (Biogeosciences) find scavenging rates of 30 - 750
+ ! (kg/m3)-1 day-1 in mesocosm experiments. Assuming that there are
+ ! 40,000 mmol C kg-1 (1 kg of pure carbon contains 83 mol and assuming
+ ! that half of marine organic particles are pure carbon by mass means
+ ! that roughly 40,000 mmol C kg-1), this translates to scavenging rates
+ ! of 0.001 to 0.02 (mmol mass particles / m3)-1 day-1.
+ call g_tracer_add_param('kscav_dfe', wombat%kscav_dfe, 0.01)
! Coagulation of dFe onto organic particles [(mmolC/m3)-1 d-1]
!-----------------------------------------------------------------------
@@ -1849,8 +1881,11 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
real, dimension(nbands) :: sw_pen
real :: swpar
real :: g_npz, g_peffect
+ real :: zooegesphyfe, zooegesdetfe
+ real :: zooassiphyfe, zooassidetfe
+ real :: zooexcrphyfe, zooexcrdetfe
real :: biophy, biozoo, biodet, biono3, biofer, biocaco3
- real :: biophyfe, biophy1, zooprey
+ real :: biophyfe, biophy1, zooprefphy, zooprefdet, zooprey
real :: fbc, zval
real, parameter :: epsi = 1.0e-30
integer :: ichl
@@ -1861,7 +1896,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
real, dimension(:), allocatable :: wsink, wsinkcal
real, dimension(4,61) :: zbgr
real, dimension(3) :: dbgr, cbgr
- real :: ztemk, I_ztemk, fe_keq, fe_par, fe_sfe, fe_tfe, partic
+ real :: ztemk, I_ztemk, fe_keq, fe_sfe, partic
real :: fesol1, fesol2, fesol3, fesol4, fesol5, hp, fe3sol
real :: biof, biodoc, zno3, zfermin
real :: phy_Fe2C, zoo_Fe2C, det_Fe2C
@@ -2078,23 +2113,26 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
wombat%phy_feupreg(:,:,:) = 0.0
wombat%phy_fedoreg(:,:,:) = 0.0
wombat%phygrow(:,:,:) = 0.0
- wombat%phyresp(:,:,:) = 0.0
- wombat%phymort(:,:,:) = 0.0
+ wombat%phymorl(:,:,:) = 0.0
+ wombat%phymorq(:,:,:) = 0.0
wombat%zooeps(:,:,:) = 0.0
wombat%zoograzphy(:,:,:) = 0.0
wombat%zoograzdet(:,:,:) = 0.0
- wombat%zooresp(:,:,:) = 0.0
- wombat%zoomort(:,:,:) = 0.0
+ wombat%zoomorl(:,:,:) = 0.0
+ wombat%zoomorq(:,:,:) = 0.0
wombat%zooexcrphy(:,:,:) = 0.0
wombat%zooexcrdet(:,:,:) = 0.0
- wombat%zooslopphy(:,:,:) = 0.0
- wombat%zooslopdet(:,:,:) = 0.0
+ wombat%zooassiphy(:,:,:) = 0.0
+ wombat%zooassidet(:,:,:) = 0.0
+ wombat%zooegesphy(:,:,:) = 0.0
+ wombat%zooegesdet(:,:,:) = 0.0
wombat%reminr(:,:,:) = 0.0
wombat%detremi(:,:,:) = 0.0
wombat%pic2poc(:,:,:) = 0.0
wombat%dissratcal(:,:,:) = 0.0
wombat%dissratara(:,:,:) = 0.0
wombat%dissratpoc(:,:,:) = 0.0
+ wombat%zoodiss(:,:,:) = 0.0
wombat%caldiss(:,:,:) = 0.0
wombat%aradiss(:,:,:) = 0.0
wombat%pocdiss(:,:,:) = 0.0
@@ -2187,17 +2225,18 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
! 2. Nutrient limitation of phytoplankton !
! 3. Temperature-dependence of heterotrophy !
! 4. Light limitation of phytoplankton !
- ! 5. Growth of chlorophyll !
- ! 6. Phytoplankton uptake of iron !
- ! 7. Iron chemistry !
- ! 8. Mortality and remineralisation !
- ! 9. Zooplankton grazing !
- ! 10. CaCO3 calculations !
- ! 11. Tracer tendencies !
- ! 12. Check for conservation by ecosystem component !
- ! 13. Additional operations on tracers !
- ! 14. Sinking rate of particulates !
- ! 15. Sedimentary processes !
+ ! 5. Realized growth of phytoplankton !
+ ! 6. Growth of chlorophyll !
+ ! 7. Phytoplankton uptake of iron !
+ ! 8. Iron chemistry !
+ ! 9. Mortality and remineralisation !
+ ! 10. Zooplankton grazing, egestion, excretion and assimilation !
+ ! 11. CaCO3 calculations !
+ ! 12. Tracer tendencies !
+ ! 13. Check for conservation by ecosystem component !
+ ! 14. Additional operations on tracers !
+ ! 15. Sinking rate of particulates !
+ ! 16. Sedimentary processes !
! !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -2395,10 +2434,21 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
! 1. Initial slope of Photosynthesis-Irradiance curve
! 2. Light limitation
- ! 3. Apply light and nutrient limitations to maximum growth rate
phy_pisl = max(wombat%alphabio * phy_chlc, wombat%alphabio * wombat%phyminqc)
wombat%phy_lpar(i,j,k) = (1. - exp(-phy_pisl * wombat%radbio(i,j,k)))
+
+
+ !-----------------------------------------------------------------------!
+ !-----------------------------------------------------------------------!
+ !-----------------------------------------------------------------------!
+ ! [Step 5] Realized growth rate of phytoplankton !
+ !-----------------------------------------------------------------------!
+ !-----------------------------------------------------------------------!
+ !-----------------------------------------------------------------------!
+
+ ! 1. Apply light and nutrient limitations to maximum growth rate
+
wombat%phy_mu(i,j,k) = wombat%phy_mumax(i,j,k) * wombat%phy_lpar(i,j,k) * &
min(wombat%phy_lnit(i,j,k), wombat%phy_lfer(i,j,k))
@@ -2412,7 +2462,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 5] Growth of chlorophyll !
+ ! [Step 6] Growth of chlorophyll !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -2439,36 +2489,37 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 6] Phytoplankton uptake of iron !
+ ! [Step 7] Phytoplankton uptake of iron !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
! 1. Maximum iron content of phytoplankton cell
! 2. Ensure that dFe uptake increases or decreases in response to cell quota
- ! 3. Iron uptake of phytoplankton (reduced to 20% at night and when N is limiting)
+ ! 3. Iron uptake of phytoplankton (reduced 10-fold in darkness)
phy_maxqfe = biophy * wombat%phymaxqf !mmol Fe / m3
wombat%phy_feupreg(i,j,k) = (4.0 - 4.5 * wombat%phy_lfer(i,j,k) / &
(wombat%phy_lfer(i,j,k) + 0.5) )
wombat%phy_fedoreg(i,j,k) = max(0.0, (1.0 - biophyfe/phy_maxqfe) / &
abs(1.05 - biophyfe/phy_maxqfe) )
- wombat%phy_dfeupt(i,j,k) = (wombat%phy_mumax(i,j,k) * wombat%phymaxqf * &
- max(0.2, wombat%phy_lpar(i,j,k) * wombat%phy_lnit(i,j,k)) * &
+ wombat%phy_dfeupt(i,j,k) = (wombat%phy_mumax(i,j,k) * phy_maxqfe * &
+ max(0.01, wombat%phy_lpar(i,j,k))**0.5 * &
biofer / (biofer + wombat%phy_kfe(i,j,k)) * &
wombat%phy_feupreg(i,j,k) * &
- wombat%phy_fedoreg(i,j,k) * biophy) * mmol_m3_to_mol_kg
+ wombat%phy_fedoreg(i,j,k)) * mmol_m3_to_mol_kg
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 7] Iron chemistry (Aumont et al., 2015; GMD) !
+ ! [Step 8] Iron chemistry (precipitation, scavenging & coagulation) !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! Estimate solubility of Fe3+ (free Fe) in solution using temperature, pH and salinity
+ ! Estimate solubility of Fe3+ (free Fe) in solution using temperature,
+ ! pH and salinity using the equations of Liu & Millero (2002)
ztemk = max(5.0, Temp(i,j,k)) + 273.15 ! temperature in kelvin
I_ztemk = 1.0 / ztemk
zval = 19.924 * Salt(i,j,k) / ( 1000. - 1.005 * Salt(i,j,k))
@@ -2485,28 +2536,34 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
endif
fe3sol = fesol1 * ( hp*hp*hp + fesol2*hp*hp + fesol3*hp + fesol4 + fesol5/hp ) *1e9
- ! Estimate total colloidal iron
- ! ... for now, we assume that 50% of all dFe is colloidal, and we separate this from the
- ! equilibrium fractionation between Fe' and Fe-L below
- wombat%fecol(i,j,k) = wombat%fcolloid * biofer
-
- ! Determine equilibriuim fractionation of the remain dFe (non-colloidal fraction) into Fe' and L-Fe
- fe_keq = 10.0**( 17.27 - 1565.7 * I_ztemk ) * 1e-9 ! Temperature reduces solubility
- fe_par = 0.47587 * wombat%radbio(i,j,k) ! Light increases solubility. dts: 0.47587=4.77e-7*0.5/10.0**(-6.3)
+ ! Estimate total colloidal iron following Tagliabue et al. (2023).
+ ! Colloidal dFe is considered to be whatever exceeds the inorganic solubility
+ ! ceiling, although there is always a hard lower limit of 10% of total dFe.
+ wombat%fecol(i,j,k) = max(0.1 * biofer, biofer - fe3sol)
+
+ ! Determine equilibriuim fractionation of the remaining dFe (non-colloidal)
+ ! between Fe' and ligand-bound iron (L-Fe). Below, temperature increases the
+ ! solubility constant (reducing free Fe) and light decreases the solubility
+ ! constant (increasing free Fe). The temperature-dependency comes from Volker
+ ! & Tagliabue (2015), while the light dependency is informed by Barbeau et al.
+ ! (2001) who saw a 0.7 log10 unit decrease in K in high light.
+ fe_keq = 1e-9 * 10.0**( (17.27 - 1565.7 * I_ztemk ) - 0.7 * &
+ wombat%radbio(i,j,k) / (wombat%radbio(i,j,k) + 10.0) )
fe_sfe = max(0.0, biofer - wombat%fecol(i,j,k))
- fe_tfe = (1.0 + fe_par) * fe_sfe
- wombat%feIII(i,j,k) = ( -( 1. + wombat%ligand * fe_keq + fe_par - fe_sfe * fe_keq ) &
- + SQRT( ( 1. + wombat%ligand * fe_keq + fe_par - fe_sfe * fe_keq )**2 &
- + 4. * fe_tfe * fe_keq) ) / ( 2. * fe_keq + epsi )
+ zval = 1.0 + wombat%ligand * fe_keq - fe_sfe * fe_keq
+ wombat%feIII(i,j,k) = ( -zval + SQRT( zval*zval + 4.0*fe_keq*fe_sfe ) ) &
+ / ( 2.*fe_keq + epsi )
+ wombat%feIII(i,j,k) = max(0.0, min(wombat%feIII(i,j,k), fe_sfe) )
wombat%felig(i,j,k) = max(0.0, fe_sfe - wombat%feIII(i,j,k))
+
! Precipitation of Fe' (creation of nanoparticles)
wombat%feprecip(i,j,k) = max(0.0, ( wombat%feIII(i,j,k) - fe3sol ) ) * wombat%knano_dfe/86400.0
! Scavenging of Fe` onto biogenic particles
- partic = (biodet + biocaco3)
+ partic = (biodet*2 + biocaco3*8.3)
wombat%fescaven(i,j,k) = wombat%feIII(i,j,k) * (1e-7 + wombat%kscav_dfe * partic) / 86400.0
- wombat%fescadet(i,j,k) = wombat%fescaven(i,j,k) * biodet / (partic+epsi)
+ wombat%fescadet(i,j,k) = wombat%fescaven(i,j,k) * biodet*2 / (partic+epsi)
! Coagulation of colloidal Fe (umol/m3) to form sinking particles (mmol/m3)
! Following Tagliabue et al. (2023), make coagulation rate dependent on DOC and Phytoplankton biomass
@@ -2532,27 +2589,27 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 8] Mortality and remineralisation !
+ ! [Step 9] Mortality and remineralisation !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
if (biophy>1e-3) then
- wombat%phyresp(i,j,k) = wombat%phylmor * fbc * wombat%f_phy(i,j,k) ! [molC/kg/s]
- wombat%phymort(i,j,k) = wombat%phyqmor / mmol_m3_to_mol_kg * wombat%f_phy(i,j,k) * wombat%f_phy(i,j,k) ! [molC/kg/s]
+ wombat%phymorl(i,j,k) = wombat%phylmor * fbc * wombat%f_phy(i,j,k) ! [molC/kg/s]
+ wombat%phymorq(i,j,k) = wombat%phyqmor / mmol_m3_to_mol_kg * wombat%f_phy(i,j,k) * wombat%f_phy(i,j,k) ! [molC/kg/s]
else
- wombat%phyresp(i,j,k) = 0.0
- wombat%phymort(i,j,k) = 0.0
+ wombat%phymorl(i,j,k) = 0.0
+ wombat%phymorq(i,j,k) = 0.0
endif
if (biozoo>1e-3) then
! reduce linear mortality (respiration losses) of zooplankton when there is low biomass
zoo_slmor = biozoo / (biozoo + wombat%zookz)
- wombat%zooresp(i,j,k) = wombat%zoolmor * fbc * wombat%f_zoo(i,j,k) * zoo_slmor ! [molC/kg/s]
- wombat%zoomort(i,j,k) = wombat%zooqmor / mmol_m3_to_mol_kg * wombat%f_zoo(i,j,k) * wombat%f_zoo(i,j,k) ! [molC/kg/s]
+ wombat%zoomorl(i,j,k) = wombat%zoolmor * fbc * wombat%f_zoo(i,j,k) * zoo_slmor ! [molC/kg/s]
+ wombat%zoomorq(i,j,k) = wombat%zooqmor / mmol_m3_to_mol_kg * wombat%f_zoo(i,j,k) * wombat%f_zoo(i,j,k) ! [molC/kg/s]
else
- wombat%zooresp(i,j,k) = 0.0
- wombat%zoomort(i,j,k) = 0.0
+ wombat%zoomorl(i,j,k) = 0.0
+ wombat%zoomorq(i,j,k) = 0.0
endif
if (wombat%f_det(i,j,k) > epsi) then
@@ -2565,37 +2622,52 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 9] Zooplankton grazing !
+ ! [Step 10] Zooplankton grazing, egestion, excretion and assimilation !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- zooprey = wombat%zprefphy * biophy + wombat%zprefdet * biodet
- ! Epsilon (prey capture rate coefficient) is made a function of
- ! phytoplankton biomass (Fig 2 of Rohr et al., 2024; GRL)
+ ! Calculate the prey biomass from dietary fractions (Gentleman et al., 2003)
+ zooprefphy = wombat%zprefphy / (wombat%zprefphy + wombat%zprefdet)
+ zooprefdet = wombat%zprefdet / (wombat%zprefphy + wombat%zprefdet)
+ zooprey = (zooprefphy * biophy + zooprefdet * biodet)
+ ! Epsilon (prey capture rate coefficient) is made a function of phytoplankton
+ ! biomass (Fig 2 of Rohr et al., 2024; GRL)
! - scales towards lower values (mesozooplankton) as prey biomass increases
g_peffect = exp(-zooprey * wombat%zooepsrat)
wombat%zooeps(i,j,k) = wombat%zooepsmin + (wombat%zooepsmax - wombat%zooepsmin) * g_peffect
g_npz = wombat%zoogmax * fbc * (wombat%zooeps(i,j,k) * zooprey*zooprey) / &
(wombat%zoogmax * fbc + (wombat%zooeps(i,j,k) * zooprey*zooprey))
+ ! We follow Le Mezo & Galbraith (2021) L&O - The fecal iron pump: ...
+ ! - egestion, assimilation and excretion of carbon and iron by zooplankton are calculated separately
+ ! - the idea is to enrich fecal pellets in iron compared to carbon
+ ! 1. zooplankton ingest C and Fe (the rest is egested)
+ ! 2. zooplankton assimilate the ingested C and Fe (the rest is excreted)
if (zooprey>1e-3) then
- wombat%zoograzphy(i,j,k) = g_npz * wombat%f_zoo(i,j,k) * (wombat%zprefphy*biophy)/zooprey ! [molC/kg/s]
- wombat%zoograzdet(i,j,k) = g_npz * wombat%f_zoo(i,j,k) * (wombat%zprefdet*biodet)/zooprey ! [molC/kg/s]
+ wombat%zoograzphy(i,j,k) = g_npz * wombat%f_zoo(i,j,k) * (zooprefphy*biophy)/zooprey ! [molC/kg/s]
+ wombat%zoograzdet(i,j,k) = g_npz * wombat%f_zoo(i,j,k) * (zooprefdet*biodet)/zooprey ! [molC/kg/s]
else
wombat%zoograzphy(i,j,k) = 0.0
wombat%zoograzdet(i,j,k) = 0.0
endif
- wombat%zooexcrphy(i,j,k) = wombat%zoograzphy(i,j,k) * (1.0 - wombat%zooassi)*wombat%zooexcr
- wombat%zooexcrdet(i,j,k) = wombat%zoograzdet(i,j,k) * (1.0 - wombat%zooassi)*wombat%zooexcr
- wombat%zooslopphy(i,j,k) = wombat%zoograzphy(i,j,k) * (1.0 - wombat%zooassi)*(1.0-wombat%zooexcr)
- wombat%zooslopdet(i,j,k) = wombat%zoograzdet(i,j,k) * (1.0 - wombat%zooassi)*(1.0-wombat%zooexcr)
-
+ wombat%zooegesphy(i,j,k) = wombat%zoograzphy(i,j,k) * (1.0-wombat%zooCingest)
+ wombat%zooegesdet(i,j,k) = wombat%zoograzdet(i,j,k) * (1.0-wombat%zooCingest)
+ wombat%zooassiphy(i,j,k) = wombat%zoograzphy(i,j,k) * wombat%zooCingest*wombat%zooCassim
+ wombat%zooassidet(i,j,k) = wombat%zoograzdet(i,j,k) * wombat%zooCingest*wombat%zooCassim
+ wombat%zooexcrphy(i,j,k) = wombat%zoograzphy(i,j,k) * wombat%zooCingest*(1.0-wombat%zooCassim)
+ wombat%zooexcrdet(i,j,k) = wombat%zoograzdet(i,j,k) * wombat%zooCingest*(1.0-wombat%zooCassim)
+ zooegesphyfe = wombat%zoograzphy(i,j,k) * phy_Fe2C * (1.0-wombat%zooFeingest)
+ zooegesdetfe = wombat%zoograzdet(i,j,k) * det_Fe2C * (1.0-wombat%zooFeingest)
+ zooassiphyfe = wombat%zoograzphy(i,j,k) * phy_Fe2C * wombat%zooFeingest*wombat%zooFeassim
+ zooassidetfe = wombat%zoograzdet(i,j,k) * det_Fe2C * wombat%zooFeingest*wombat%zooFeassim
+ zooexcrphyfe = wombat%zoograzphy(i,j,k) * phy_Fe2C * wombat%zooFeingest*(1.0-wombat%zooFeassim)
+ zooexcrdetfe = wombat%zoograzdet(i,j,k) * det_Fe2C * wombat%zooFeingest*(1.0-wombat%zooFeassim)
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 10] CaCO3 calculations !
+ ! [Step 11] CaCO3 calculations !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -2610,7 +2682,8 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
hco3 / wombat%htotal(i,j,k))) * &
(0.55 + 0.45 * tanh(Temp(i,j,k) - 4.0)) )
! The dissolution rate is a function of omegas for calcite and aragonite, as well the
- ! concentration of POC, following Kwon et al., 2024, Science Advances; Table S1
+ ! concentration of POC, following Kwon et al., 2024, Science Advances; Table S1, and
+ ! we account for the dissolution due to zooplankton grazing on particulates
wombat%dissratcal(i,j,k) = (wombat%disscal * max(0.0, 1.0 - wombat%omega_cal(i,j,k))**2.2) / 86400.0
wombat%dissratara(i,j,k) = (wombat%dissara * max(0.0, 1.0 - wombat%omega_ara(i,j,k))**1.5) / 86400.0
wombat%dissratpoc(i,j,k) = (wombat%dissdet * wombat%reminr(i,j,k) * biodet**2.0)
@@ -2622,10 +2695,12 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
endif
if (wombat%f_caco3(i,j,k) > epsi) then
+ wombat%zoodiss(i,j,k) = wombat%zoograzdet(i,j,k) * wombat%fgutdiss * biocaco3/biodet
wombat%caldiss(i,j,k) = wombat%dissratcal(i,j,k) * wombat%f_caco3(i,j,k) ! [mol/kg/s]
wombat%aradiss(i,j,k) = wombat%dissratara(i,j,k) * wombat%f_caco3(i,j,k) ! [mol/kg/s]
wombat%pocdiss(i,j,k) = wombat%dissratpoc(i,j,k) * wombat%f_caco3(i,j,k) ! [mol/kg/s]
else
+ wombat%zoodiss(i,j,k) = 0.0
wombat%caldiss(i,j,k) = 0.0
wombat%aradiss(i,j,k) = 0.0
wombat%pocdiss(i,j,k) = 0.0
@@ -2635,7 +2710,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 11] Tracer tendencies !
+ ! [Step 12] Tracer tendencies !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -2644,34 +2719,34 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!----------------------------------------------------------------------
wombat%f_no3(i,j,k) = wombat%f_no3(i,j,k) + dtsb * 16./122. * ( &
wombat%detremi(i,j,k) + &
- wombat%zooresp(i,j,k) + &
+ wombat%zoomorl(i,j,k) + &
wombat%zooexcrphy(i,j,k) + &
wombat%zooexcrdet(i,j,k) + &
- wombat%phyresp(i,j,k) - &
+ wombat%phymorl(i,j,k) - &
wombat%phygrow(i,j,k) )
! Phytoplankton equation ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_phy(i,j,k) = wombat%f_phy(i,j,k) + dtsb * ( &
wombat%phygrow(i,j,k) - &
- wombat%phyresp(i,j,k) - &
- wombat%phymort(i,j,k) - &
+ wombat%phymorl(i,j,k) - &
+ wombat%phymorq(i,j,k) - &
wombat%zoograzphy(i,j,k) )
! Phytoplankton chlorophyll equation ! [molChl/kg]
!-----------------------------------------------------------------------
wombat%f_pchl(i,j,k) = wombat%f_pchl(i,j,k) + dtsb * ( &
wombat%pchl_mu(i,j,k) - &
- wombat%phyresp(i,j,k) * phy_chlc - &
- wombat%phymort(i,j,k) * phy_chlc - &
+ wombat%phymorl(i,j,k) * phy_chlc - &
+ wombat%phymorq(i,j,k) * phy_chlc - &
wombat%zoograzphy(i,j,k) * phy_chlc )
! Phytoplankton iron equation ! [molFe/kg]
!-----------------------------------------------------------------------
wombat%f_phyfe(i,j,k) = wombat%f_phyfe(i,j,k) + dtsb * ( &
wombat%phy_dfeupt(i,j,k) - &
- wombat%phyresp(i,j,k) * phy_Fe2C - &
- wombat%phymort(i,j,k) * phy_Fe2C - &
+ wombat%phymorl(i,j,k) * phy_Fe2C - &
+ wombat%phymorq(i,j,k) * phy_Fe2C - &
wombat%zoograzphy(i,j,k) * phy_Fe2C )
! Net primary productivity ! [molC/kg/s]
@@ -2680,40 +2755,41 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
! Zooplankton equation ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_zoo(i,j,k) = wombat%f_zoo(i,j,k) + dtsb * ( &
- wombat%zooassi * wombat%zoograzphy(i,j,k) + &
- wombat%zooassi * wombat%zoograzdet(i,j,k) - &
- wombat%zooresp(i,j,k) - &
- wombat%zoomort(i,j,k) )
+ wombat%zooassiphy(i,j,k) + &
+ wombat%zooassidet(i,j,k) - &
+ wombat%zoomorl(i,j,k) - &
+ wombat%zoomorq(i,j,k) )
! Zooplankton iron equation ! [molFe/kg]
!-----------------------------------------------------------------------
wombat%f_zoofe(i,j,k) = wombat%f_zoofe(i,j,k) + dtsb * ( &
- wombat%zooassi * wombat%zoograzphy(i,j,k) * phy_Fe2C + &
- wombat%zooassi * wombat%zoograzdet(i,j,k) * det_Fe2C - &
- wombat%zooresp(i,j,k) * zoo_Fe2C - &
- wombat%zoomort(i,j,k) * zoo_Fe2C )
+ zooassiphyfe + &
+ zooassidetfe - &
+ wombat%zoomorl(i,j,k) * zoo_Fe2C - &
+ wombat%zoomorq(i,j,k) * zoo_Fe2C )
! Estimate secondary productivity from zooplankton growth ! [molC/kg/s]
wombat%zsp3d(i,j,k) = wombat%zsp3d(i,j,k) + dtsb * &
- wombat%zooassi * (wombat%zoograzphy(i,j,k) + wombat%zoograzdet(i,j,k))
+ wombat%zooCingest*wombat%zooCassim * &
+ (wombat%zoograzphy(i,j,k) + wombat%zoograzdet(i,j,k))
! Detritus equation ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_det(i,j,k) = wombat%f_det(i,j,k) + dtsb * ( &
- wombat%zooslopphy(i,j,k) + &
- wombat%zooslopdet(i,j,k) + &
- wombat%phymort(i,j,k) + &
- wombat%zoomort(i,j,k) - &
+ wombat%zooegesphy(i,j,k) + &
+ wombat%zooegesdet(i,j,k) + &
+ wombat%phymorq(i,j,k) + &
+ wombat%zoomorq(i,j,k) - &
wombat%zoograzdet(i,j,k) - &
wombat%detremi(i,j,k) )
! Detrital iron equation ! [molFe/kg]
!-----------------------------------------------------------------------
wombat%f_detfe(i,j,k) = wombat%f_detfe(i,j,k) + dtsb * ( &
- wombat%zooslopphy(i,j,k) * phy_Fe2C + &
- wombat%zooslopdet(i,j,k) * det_Fe2C + &
- wombat%phymort(i,j,k) * phy_Fe2C + &
- wombat%zoomort(i,j,k) * zoo_Fe2C - &
+ zooegesphyfe + &
+ zooegesdetfe + &
+ wombat%phymorq(i,j,k) * phy_Fe2C + &
+ wombat%zoomorq(i,j,k) * zoo_Fe2C - &
wombat%zoograzdet(i,j,k) * det_Fe2C - &
wombat%detremi(i,j,k) * det_Fe2C + &
wombat%fescadet(i,j,k) + &
@@ -2724,70 +2800,74 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
if (wombat%f_o2(i,j,k) > epsi) &
wombat%f_o2(i,j,k) = wombat%f_o2(i,j,k) - 172./122. * dtsb * ( &
wombat%detremi(i,j,k) + &
- wombat%zooresp(i,j,k) + &
+ wombat%zoomorl(i,j,k) + &
wombat%zooexcrphy(i,j,k) + &
wombat%zooexcrdet(i,j,k) + &
- wombat%phyresp(i,j,k) - &
+ wombat%phymorl(i,j,k) - &
wombat%phygrow(i,j,k) )
! Equation for CaCO3 ! [molCaCO3/kg]
!-----------------------------------------------------------------------
wombat%f_caco3(i,j,k) = wombat%f_caco3(i,j,k) + dtsb * ( &
- wombat%zoograzphy(i,j,k) * (1. - wombat%fgutdiss) * wombat%pic2poc(i,j,k) + &
- wombat%phymort(i,j,k) * wombat%pic2poc(i,j,k) + &
- wombat%zoomort(i,j,k) * wombat%pic2poc(i,j,k) - &
- wombat%caldiss(i,j,k) - wombat%aradiss(i,j,k) - wombat%pocdiss(i,j,k) )
+ wombat%phymorq(i,j,k) * wombat%pic2poc(i,j,k) + &
+ wombat%zoomorq(i,j,k) * wombat%pic2poc(i,j,k) + &
+ wombat%zoograzphy(i,j,k) * (1. - wombat%fgutdiss) * wombat%pic2poc(i,j,k) - &
+ wombat%zoodiss(i,j,k) - wombat%caldiss(i,j,k) - &
+ wombat%aradiss(i,j,k) - wombat%pocdiss(i,j,k) )
! Equation for DIC ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_dic(i,j,k) = wombat%f_dic(i,j,k) + dtsb * ( &
wombat%detremi(i,j,k) + &
- wombat%zooresp(i,j,k) + &
+ wombat%zoomorl(i,j,k) + &
wombat%zooexcrphy(i,j,k) + &
wombat%zooexcrdet(i,j,k) + &
- wombat%phyresp(i,j,k) - &
+ wombat%phymorl(i,j,k) - &
wombat%phygrow(i,j,k) - &
wombat%zoograzphy(i,j,k) * (1.0-wombat%fgutdiss) * wombat%pic2poc(i,j,k) - &
- wombat%phymort(i,j,k) * wombat%pic2poc(i,j,k) - &
- wombat%zoomort(i,j,k) * wombat%pic2poc(i,j,k) + &
- wombat%caldiss(i,j,k) + wombat%aradiss(i,j,k) + wombat%pocdiss(i,j,k) )
+ wombat%phymorq(i,j,k) * wombat%pic2poc(i,j,k) - &
+ wombat%zoomorq(i,j,k) * wombat%pic2poc(i,j,k) + &
+ wombat%zoodiss(i,j,k) + wombat%caldiss(i,j,k) + &
+ wombat%aradiss(i,j,k) + wombat%pocdiss(i,j,k) )
! Equation for DICr ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_dicr(i,j,k) = wombat%f_dicr(i,j,k) + dtsb * ( &
wombat%detremi(i,j,k) + &
- wombat%zooresp(i,j,k) + &
+ wombat%zoomorl(i,j,k) + &
wombat%zooexcrphy(i,j,k) + &
wombat%zooexcrdet(i,j,k) + &
- wombat%phyresp(i,j,k) - &
+ wombat%phymorl(i,j,k) - &
wombat%phygrow(i,j,k) - &
wombat%zoograzphy(i,j,k) * (1.0-wombat%fgutdiss) * wombat%pic2poc(i,j,k) - &
- wombat%phymort(i,j,k) * wombat%pic2poc(i,j,k) - &
- wombat%zoomort(i,j,k) * wombat%pic2poc(i,j,k) + &
- wombat%caldiss(i,j,k) + wombat%aradiss(i,j,k) + wombat%pocdiss(i,j,k) )
+ wombat%phymorq(i,j,k) * wombat%pic2poc(i,j,k) - &
+ wombat%zoomorq(i,j,k) * wombat%pic2poc(i,j,k) + &
+ wombat%zoodiss(i,j,k) + wombat%caldiss(i,j,k) + &
+ wombat%aradiss(i,j,k) + wombat%pocdiss(i,j,k) )
! Equation for ALK ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_alk(i,j,k) = wombat%f_alk(i,j,k) + dtsb * 16.0/122.0 * ( &
wombat%phygrow(i,j,k) - &
wombat%detremi(i,j,k) - &
- wombat%zooresp(i,j,k) - &
+ wombat%zoomorl(i,j,k) - &
wombat%zooexcrphy(i,j,k) - &
wombat%zooexcrdet(i,j,k) - &
- wombat%phyresp(i,j,k) ) + dtsb * 2.0 * ( &
- wombat%caldiss(i,j,k) + wombat%aradiss(i,j,k) + wombat%pocdiss(i,j,k) - &
+ wombat%phymorl(i,j,k) ) + dtsb * 2.0 * ( &
+ wombat%zoodiss(i,j,k) + wombat%caldiss(i,j,k) + &
+ wombat%aradiss(i,j,k) + wombat%pocdiss(i,j,k) + &
wombat%zoograzphy(i,j,k) * (1.0-wombat%fgutdiss) * wombat%pic2poc(i,j,k) - &
- wombat%phymort(i,j,k) * wombat%pic2poc(i,j,k) - &
- wombat%zoomort(i,j,k) * wombat%pic2poc(i,j,k) )
+ wombat%phymorq(i,j,k) * wombat%pic2poc(i,j,k) - &
+ wombat%zoomorq(i,j,k) * wombat%pic2poc(i,j,k) )
! Extra equation for iron ! [molFe/kg]
!-----------------------------------------------------------------------
wombat%f_fe(i,j,k) = wombat%f_fe(i,j,k) + dtsb * ( &
wombat%detremi(i,j,k) * det_Fe2C + &
- wombat%zooresp(i,j,k) * zoo_Fe2C + &
- wombat%zooexcrphy(i,j,k) * phy_Fe2C + &
- wombat%zooexcrdet(i,j,k) * det_Fe2C + &
- wombat%phyresp(i,j,k) * phy_Fe2C - &
+ wombat%zoomorl(i,j,k) * zoo_Fe2C + &
+ zooexcrphyfe + &
+ zooexcrdetfe + &
+ wombat%phymorl(i,j,k) * phy_Fe2C - &
wombat%phy_dfeupt(i,j,k) - &
wombat%feprecip(i,j,k) - &
wombat%fescaven(i,j,k) - &
@@ -2796,10 +2876,10 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
! Collect dFe sources and sinks for diagnostic output
wombat%fesources(i,j,k) = wombat%fesources(i,j,k) + dtsb * ( &
wombat%detremi(i,j,k) * det_Fe2C + &
- wombat%zooresp(i,j,k) * zoo_Fe2C + &
- wombat%zooexcrphy(i,j,k) * phy_Fe2C + &
- wombat%zooexcrdet(i,j,k) * det_Fe2C + &
- wombat%phyresp(i,j,k) * phy_Fe2C)
+ wombat%zoomorl(i,j,k) * zoo_Fe2C + &
+ zooexcrphyfe + &
+ zooexcrdetfe + &
+ wombat%phymorl(i,j,k) * phy_Fe2C)
wombat%fesinks(i,j,k) = wombat%fesinks(i,j,k) + dtsb * ( &
wombat%phy_dfeupt(i,j,k) + &
wombat%feprecip(i,j,k) + &
@@ -2810,7 +2890,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 12] Check for conservation of mass by ecosystem component !
+ ! [Step 13] Check for conservation of mass by ecosystem component !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -2890,7 +2970,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 13] Additional operations on tracers !
+ ! [Step 14] Additional operations on tracers !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -2931,7 +3011,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 14] Sinking of particulates !
+ ! [Step 15] Sinking of particulates !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -2970,7 +3050,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
- ! [Step 14] Sedimentary processes !
+ ! [Step 16] Sedimentary processes !
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
@@ -3229,12 +3309,12 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
used = g_send_data(wombat%id_phygrow, wombat%phygrow, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
- if (wombat%id_phyresp > 0) &
- used = g_send_data(wombat%id_phyresp, wombat%phyresp, model_time, &
+ if (wombat%id_phymorl > 0) &
+ used = g_send_data(wombat%id_phymorl, wombat%phymorl, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
- if (wombat%id_phymort > 0) &
- used = g_send_data(wombat%id_phymort, wombat%phymort, model_time, &
+ if (wombat%id_phymorq > 0) &
+ used = g_send_data(wombat%id_phymorq, wombat%phymorq, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (wombat%id_zooeps > 0) &
@@ -3249,12 +3329,12 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
used = g_send_data(wombat%id_zoograzdet, wombat%zoograzdet, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
- if (wombat%id_zooresp > 0) &
- used = g_send_data(wombat%id_zooresp, wombat%zooresp, model_time, &
+ if (wombat%id_zoomorl > 0) &
+ used = g_send_data(wombat%id_zoomorl, wombat%zoomorl, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
- if (wombat%id_zoomort > 0) &
- used = g_send_data(wombat%id_zoomort, wombat%zoomort, model_time, &
+ if (wombat%id_zoomorq > 0) &
+ used = g_send_data(wombat%id_zoomorq, wombat%zoomorq, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (wombat%id_zooexcrphy > 0) &
@@ -3265,12 +3345,20 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
used = g_send_data(wombat%id_zooexcrdet, wombat%zooexcrdet, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
- if (wombat%id_zooslopphy > 0) &
- used = g_send_data(wombat%id_zooslopphy, wombat%zooslopphy, model_time, &
+ if (wombat%id_zooassiphy > 0) &
+ used = g_send_data(wombat%id_zooassiphy, wombat%zooassiphy, model_time, &
+ rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
+
+ if (wombat%id_zooassidet > 0) &
+ used = g_send_data(wombat%id_zooassidet, wombat%zooassidet, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
- if (wombat%id_zooslopdet > 0) &
- used = g_send_data(wombat%id_zooslopdet, wombat%zooslopdet, model_time, &
+ if (wombat%id_zooegesphy > 0) &
+ used = g_send_data(wombat%id_zooegesphy, wombat%zooegesphy, model_time, &
+ rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
+
+ if (wombat%id_zooegesdet > 0) &
+ used = g_send_data(wombat%id_zooegesdet, wombat%zooegesdet, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (wombat%id_reminr > 0) &
@@ -3297,6 +3385,10 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
used = g_send_data(wombat%id_dissratpoc, wombat%dissratpoc, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
+ if (wombat%id_zoodiss > 0) &
+ used = g_send_data(wombat%id_zoodiss, wombat%zoodiss, model_time, &
+ rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
+
if (wombat%id_caldiss > 0) &
used = g_send_data(wombat%id_caldiss, wombat%caldiss, model_time, &
rmask=grid_tmask, is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
@@ -3680,23 +3772,26 @@ subroutine user_allocate_arrays
allocate(wombat%phy_feupreg(isd:ied, jsd:jed, 1:nk)); wombat%phy_feupreg(:,:,:)=0.0
allocate(wombat%phy_fedoreg(isd:ied, jsd:jed, 1:nk)); wombat%phy_fedoreg(:,:,:)=0.0
allocate(wombat%phygrow(isd:ied, jsd:jed, 1:nk)); wombat%phygrow(:,:,:)=0.0
- allocate(wombat%phyresp(isd:ied, jsd:jed, 1:nk)); wombat%phyresp(:,:,:)=0.0
- allocate(wombat%phymort(isd:ied, jsd:jed, 1:nk)); wombat%phymort(:,:,:)=0.0
+ allocate(wombat%phymorl(isd:ied, jsd:jed, 1:nk)); wombat%phymorl(:,:,:)=0.0
+ allocate(wombat%phymorq(isd:ied, jsd:jed, 1:nk)); wombat%phymorq(:,:,:)=0.0
allocate(wombat%zooeps(isd:ied, jsd:jed, 1:nk)); wombat%zooeps(:,:,:)=0.0
allocate(wombat%zoograzphy(isd:ied, jsd:jed, 1:nk)); wombat%zoograzphy(:,:,:)=0.0
allocate(wombat%zoograzdet(isd:ied, jsd:jed, 1:nk)); wombat%zoograzdet(:,:,:)=0.0
- allocate(wombat%zooresp(isd:ied, jsd:jed, 1:nk)); wombat%zooresp(:,:,:)=0.0
- allocate(wombat%zoomort(isd:ied, jsd:jed, 1:nk)); wombat%zoomort(:,:,:)=0.0
+ allocate(wombat%zoomorl(isd:ied, jsd:jed, 1:nk)); wombat%zoomorl(:,:,:)=0.0
+ allocate(wombat%zoomorq(isd:ied, jsd:jed, 1:nk)); wombat%zoomorq(:,:,:)=0.0
allocate(wombat%zooexcrphy(isd:ied, jsd:jed, 1:nk)); wombat%zooexcrphy(:,:,:)=0.0
allocate(wombat%zooexcrdet(isd:ied, jsd:jed, 1:nk)); wombat%zooexcrdet(:,:,:)=0.0
- allocate(wombat%zooslopphy(isd:ied, jsd:jed, 1:nk)); wombat%zooslopphy(:,:,:)=0.0
- allocate(wombat%zooslopdet(isd:ied, jsd:jed, 1:nk)); wombat%zooslopdet(:,:,:)=0.0
+ allocate(wombat%zooassiphy(isd:ied, jsd:jed, 1:nk)); wombat%zooassiphy(:,:,:)=0.0
+ allocate(wombat%zooassidet(isd:ied, jsd:jed, 1:nk)); wombat%zooassidet(:,:,:)=0.0
+ allocate(wombat%zooegesphy(isd:ied, jsd:jed, 1:nk)); wombat%zooegesphy(:,:,:)=0.0
+ allocate(wombat%zooegesdet(isd:ied, jsd:jed, 1:nk)); wombat%zooegesdet(:,:,:)=0.0
allocate(wombat%reminr(isd:ied, jsd:jed, 1:nk)); wombat%reminr(:,:,:)=0.0
allocate(wombat%detremi(isd:ied, jsd:jed, 1:nk)); wombat%detremi(:,:,:)=0.0
allocate(wombat%pic2poc(isd:ied, jsd:jed, 1:nk)); wombat%pic2poc(:,:,:)=0.0
allocate(wombat%dissratcal(isd:ied, jsd:jed, 1:nk)); wombat%dissratcal(:,:,:)=0.0
allocate(wombat%dissratara(isd:ied, jsd:jed, 1:nk)); wombat%dissratara(:,:,:)=0.0
allocate(wombat%dissratpoc(isd:ied, jsd:jed, 1:nk)); wombat%dissratpoc(:,:,:)=0.0
+ allocate(wombat%zoodiss(isd:ied, jsd:jed, 1:nk)); wombat%zoodiss(:,:,:)=0.0
allocate(wombat%caldiss(isd:ied, jsd:jed, 1:nk)); wombat%caldiss(:,:,:)=0.0
allocate(wombat%aradiss(isd:ied, jsd:jed, 1:nk)); wombat%aradiss(:,:,:)=0.0
allocate(wombat%pocdiss(isd:ied, jsd:jed, 1:nk)); wombat%pocdiss(:,:,:)=0.0
@@ -3801,22 +3896,25 @@ subroutine user_deallocate_arrays
wombat%phy_feupreg, &
wombat%phy_fedoreg, &
wombat%phygrow, &
- wombat%phyresp, &
- wombat%phymort, &
+ wombat%phymorl, &
+ wombat%phymorq, &
wombat%zoograzphy, &
wombat%zoograzdet, &
- wombat%zooresp, &
- wombat%zoomort, &
+ wombat%zoomorl, &
+ wombat%zoomorq, &
wombat%zooexcrphy, &
wombat%zooexcrdet, &
- wombat%zooslopphy, &
- wombat%zooslopdet, &
+ wombat%zooassiphy, &
+ wombat%zooassidet, &
+ wombat%zooegesphy, &
+ wombat%zooegesdet, &
wombat%reminr, &
wombat%detremi, &
wombat%pic2poc, &
wombat%dissratcal, &
wombat%dissratara, &
wombat%dissratpoc, &
+ wombat%zoodiss, &
wombat%caldiss, &
wombat%aradiss, &
wombat%pocdiss, &
diff --git a/generic_tracers/generic_WOMBATmid.F90 b/generic_tracers/generic_WOMBATmid.F90
index 359b71a8..03dd14f5 100644
--- a/generic_tracers/generic_WOMBATmid.F90
+++ b/generic_tracers/generic_WOMBATmid.F90
@@ -4056,6 +4056,7 @@ subroutine generic_WOMBATmid_update_from_source(tracer_list, Temp, Salt, &
real :: biof, zno3, zfermin, shear
real :: phy_Fe2C, dia_Fe2C, zoo_Fe2C, mes_Fe2C, det_Fe2C, bdet_Fe2C
real :: dom_N2C, dia_Si2C, bdet_Si2C
+ real :: theta_opt
real :: phy_minqfe, phy_maxqfe
real :: dia_minqfe, dia_maxqfe, dia_maxqsi
real :: zoo_slmor, mes_slmor