I am coding a module to do individual-based modeling in ROMS. But I ran into a problem that I can't figure out. I need help from the ROMS' Gurus!!!!
THE SETUP:
As a simple example, I added an oyster in the middle of a 3 by 3 by 3 grid with walls on all sides. Each cell is roughly 100m x 100m x 1m. Atmospheric forcing is provided in a NetCDF forcing file, and I have no tides.
THE PROBLEM:
Nitrogen is not conserved in the system.
The oyster grows... while phytoplankton is eaten (in more sophisticated versions, ammonia is excreted, oyster babies are produced ... and so on). It works pretty cool EXCEPT for the not-so-minor problem that when I add all the nitrogen I find it decreases over time.
Insight in the solution:
I think the problem resides in that all biological tracers continue into step3d_t.F where they are modified using the new zeta and Hz... while the oyster biomass doesn't enter step3d_t.F (I don't want the oyster to get diffused or advected). Perhaps this means that there is a unit disagreement in my code where all biological tracers are in, say, transport units... while the oyster is in biomass per individual. But I am just guessing...
I was able to "fix" the nitrogen leakage by dividing oyster biomass by 4 BEFORE the reaction equations in my code ...followed by multiplying oyster biomass by 4 AFTER the reaction equations (see example below).
CASE 1: This case results in nitrogen leakage.
! -----------------------------------------------------------------------
! Extract phytoplankton from tracer array -----
Phyto = t(2,2,2,3,iPhyto)
! Extract Individual Biomass and Parameters from
! variables already stored in ROMS
Oyster = MyOyster
Rate = MyFiltrationRate
! React biological tracers with individual-based biology
! Note: dtdays = time-step (in days) and V is the cell volume
! calculated using Hz and grdscl
Phyto = Phyto - (Rate * Phyto * dtdays)
Oyster = Oyster + (Rate * Phyto * (1/V) * dtdays)
! Update tracer array and Oyster biomass
MyOyster = Oyster
t(2,2,2,3,iPhyto) = Phyto
! -----------------------------------------------------------------------
CASE 2: This case is FIXED (i.e. nitrogen is conserved).
! -----------------------------------------------------------------------
! Extract phytoplankton from tracer array -----
Phyto = t(2,2,2,3,iPhyto)
! Extract Individual Biomass and Parameters from
! variables already stored in ROMS
Oyster = MyOyster / 4.0_r8
Rate = MyFiltrationRate
! React biological tracers with individual-based biology
! Note: dtdays = time-step (in days) and V is the cell volume
! calculated using Hz and grdscl
Phyto = Phyto - (Rate * Phyto * dtdays)
Oyster = Oyster + (Rate * Phyto * (1/V) * dtdays)
! Update tracer array and Oyster biomass
MyOyster = Oyster * 4.0_r8
t(2,2,2,3,iPhyto) = Phyto
! -----------------------------------------------------------------------
Why 4???...
... I have no idea!
Here are some properties of the problem...
1. Setting MyFiltrationRate = 0 results in conserved nitrogen over time (i.e. the nitrogen leak stops)... however (obviously) the oyster doesn't grow and phytoplankton is not eaten.
2. Changing MyFiltrationRate or MyOyster does not affect nitrogen conservation.
3. Changing grid configuration DOES affect nitrogen conservation. The magic /4 and *4 solution ONLY works for the 3-by-3-by-3 grid. Changing the grid size or the number of cells in the domain... requires dividing/multiplying by a factor different than 4 (sometimes with decimals).
I will keep doing tests... and will post here any clues I find.
Thank you all for your help!!!
I know this is a hard one.
Cheers,
Diego