NEMURO Model

[arango]

ROMS has now the NEMURO lower trophic ecosystem model developed for the North Pacific. It has 11 functional compartments:

Code: Select all

     iSphy     small phytoplankton biomass, nanophytoplankton
     iLphy     large phytoplankton biomass, diatoms
     iSzoo     small zooplankton biomass, microzooplankton (ciliates)
     iLzoo     large zooplankton biomass, mesozooplankton (copepods)
     iPzoo     predator zooplankton biomass (euphausiids, etc)
     iNO3_     nitrate concentration, NO3
     iNH4_     ammonium concentration, NH4
     iPON_     particulate organic nitrogen
     iDON_     dissolved organic nitrogen
     iSiOH     silicate concentration, Si(OH)4 (silicic acid)
     iopal     particulate organic silica

The basic reference for this model is Kishi et al. (2007) , Ecological Modelling, 202, 12-25.

I coded this model from scratch using the same numerical scheme of other existing biological models in ROMS. The compartment being consumed is treated implicit to guarantee conservation, a positive definite and unconditionally stable algorithm. The model can be activated with the NEMURO cpp option.

Currently, it has three formulations for the grazing terms:

• Explit Ivlev grazing formulation, original code. Activated with option IVLEV_EXPLICIT.

• Holling-type s-shape grazing curve with uses a quadratic term on zooplankton. It is known to be numerically more stable and allows an implicit formulation. Activated with option HOLLING_GRAZING.

• Implicit Ivlev grazing formulation, default option. The implicit algorithm is possible after some algebraic transformation and approximation with Taylor series. Many thanks to Chris Edwards for his help in this formulation. This formulation does not work well yet. The model blows-up sometimes. It is recommended to use either IVLEV_EXPLICIT or HOLLING_GRAZING

The light formulation follows Platt et al. (1980) photoinhibition equation. As in all coupled models in ROMS, the sinking term is reconstructed with monotonic (PPM/WENO) parabolic segments and a semi-Lagrangian advective flux with no CFL constraints. This allows very fast sinking rates. The particulate fluxes reaching the seafloor (PON, Opal) are returned to the dissolved nutrient pool in the water column using the BIO_SEDIMENT option.

This model is available in the ROMS svn repository (revision 82 or higher). The input script nemuro.in can be found in ROMS/External. There is also a copy in User/External. As in all ecosystem models in ROMS, the compartment concentrations are in millimoles/m3. I added a simple initialization in ana_biology.h which is located in ROMS/Functionals with a copy in User/Functionals.

Warning: the model is not fully tested . It now goes to the beta-testers .

[hermann]

Hi group,

There appears to be a problem with units in the originally released nemuro.in, such that many of the values are off by a factor of 1.e3. I have confirmed this with Hernan. Here is a diff of the presently distributed file, and what I believe to be the correct values (that is, values corresponding to the Kishi et al. paper):

Code: Select all

% diff erroneous_nemuro.in corrected_nemuro.in
84,85c84,85
<        KNO3S == 1.0d-3                     ! small biomass
<        KNO3L == 3.00d-3                    ! large biomass
---
>        KNO3S == 1.0d0                     ! small biomass
>        KNO3L == 3.0d0                    ! large biomass
89,90c89,90
<        KNH4S == 0.1d-3                     ! small biomass
<        KNH4L == 0.30d-3                    ! large biomass
---
>        KNH4S == 0.1d0                     ! small biomass
>        KNH4L == 0.30d0                    ! large biomass
94c94
<         KSiL == 6.00d-3                    ! large biomass
---
>         KSiL == 6.00d0                    ! large biomass
98,99c98,99
<       PusaiS == 1.5d3                      ! small biomass
<       PusaiL == 1.5d3                      ! large biomass
---
>       PusaiS == 1.5d0                      ! small biomass
>       PusaiL == 1.5d0                      ! large biomass
124,125c124,125
<       MorPS0 == 58.5d0                     ! small biomass
<       MorPL0 == 29.0d0                     ! large biomass
---
>       MorPS0 == 58.5d-3                     ! small biomass
>       MorPL0 == 29.0d-3                    ! large biomass
150,152c150,152
<         LamS == 1.4d3                      ! small biomass
<         LamL == 1.4d3                      ! large biomass
<         LamP == 1.4d3                      ! predator biomass
---
>         LamS == 1.4d0                      ! small biomass
>         LamL == 1.4d0                      ! large biomass
>         LamP == 1.4d0                      ! predator biomass
167,173c167,173
<    PS2ZSstar == 4.3d-5                     ! small Zoo on small Phy
<    PS2ZLstar == 4.0d-5                     ! large Zoo on small Phy
<    PL2ZLstar == 4.0d-5                     ! large Zoo on large Phy
<    ZS2ZLstar == 4.0d-5                     ! large Zoo on small Zoo
<    PL2ZPstar == 4.0d-5                     ! predator Zoo on large Phy
<    ZS2ZPstar == 4.0d-5                     ! predator Zoo on small Zoo
<    ZL2ZPstar == 4.0d-5                     ! predator Zoo on large Zoo
---
>    PS2ZSstar == 4.3d-2                     ! small Zoo on small Phy
>    PS2ZLstar == 4.0d-2                     ! large Zoo on small Phy
>    PL2ZLstar == 4.0d-2                     ! large Zoo on large Phy
>    ZS2ZLstar == 4.0d-2                     ! large Zoo on small Zoo
>    PL2ZPstar == 4.0d-2                     ! predator Zoo on large Phy
>    ZS2ZPstar == 4.0d-2                     ! predator Zoo on small Zoo
>    ZL2ZPstar == 4.0d-2                     ! predator Zoo on large Zoo
177,178c177,178
<      PusaiPL == 4.605d3                    ! predator Zoo on large Phy
<      PusaiZS == 3.010d3                    ! predator Zoo on small Zoo
---
>      PusaiPL == 4.605d0                    ! predator Zoo on large Phy
>      PusaiZS == 3.010d0                    ! predator Zoo on small Zoo
182,184c182,184
<       MorZS0 == 58.5d0                     ! small biomass
<       MorZL0 == 58.5d0                     ! large biomass
<       MorZP0 == 58.5d0                     ! predator biomass
---
>       MorZS0 == 58.5d-3                    ! small biomass
>       MorZL0 == 58.5d-3                    ! large biomass
>       MorZP0 == 58.5d-3                    ! predator biomass

-Al

[arango]

Yes, thank you. I updated the repository with the correct values.

[hermann]

Several of us have recently been debugging the ROMS NEMURO code. One way of testing the code is to see whether all categories of nitrogen taken together ("total nitrogen") is conserved in a run. When I run my thus-far debugged version (using Holling grazing) in 1D (BIO_TOY), and set the sinking rates to zero, I get excellent (on the order of 0.1% per year) conservation of total nitrogen. But with sinking turned on at default rates (40m/d), there is a persistent leak of total nitrogen from the system, even when the BIO_SEDIMENT option is turned on (BIO_SEDIMENT is designed to add back in as NO3 whatever falls out the bottom of the domain as PON). I am wondering if the sinking algorithm is non-conserving in some way - it does involve spline interpolation; maybe this is part of the problem at high sinking rates. Does anyone have experience with this non-conservation problem in one of the other biological or sediment options?

Below is the diff of the the posted NEMURO code (nemuro-original.h), vs. what I am currently using (nemuro-modified.h):

---------------------------------------------------------------

/home/users/hermann/ROMS2007/ROMS/Nonlinear/ 7 % diff nemuro-original.h nemuro-modified.h
367c367
< RnewS=GppNPS/GppPS
---
> RnewS=GppNPS/(GppPS+MinVal)
404c404,405
< RnewL=cff4/(cff4+cff5)
---
> !AJH added MinVal in the next line
> RnewL=cff4/(cff4+cff5+MinVal)
406,407c407,411
< cff4=cff6*RnewL
< cff5=cff6*(1.0_r8-RnewL)
---
>
> !AJH added normalization in next two lines
> cff4=cff6*RnewL/(Bio(i,k,iNO3_)+MinVal)
> cff5=cff6*(1.0_r8-RnewL)/(Bio(i,k,iNH4_)+MinVal)
>
415c419
< Bio(i,k,iLphy)=Bio(i,k,iSphy)+GppPL
---
> Bio(i,k,iLphy)=Bio(i,k,iLphy)+GppPL
423c427
< RnewL=GppNPL/GppPL
---
> RnewL=GppNPL/(GppPL+MinVal)
535c539
< cff4=1.0_r8-EXP(LamS(ng)*(PS2ZSstar(ng)-Bio(i,k,iSphy))
---
> cff4=1.0_r8-EXP(LamS(ng)*(PS2ZSstar(ng)-Bio(i,k,iSphy)))
556c560
< cff4=1.0_r8-EXP(LamL(ng)*(PS2ZLstar(ng)-Bio(i,k,iSphy))
---
> cff4=1.0_r8-EXP(LamL(ng)*(PS2ZLstar(ng)-Bio(i,k,iSphy)))
611c615,616
< Bio(i,k,iSzoo)=Bio(i,k,iSzoo)*(1.0_r8+cff)
---
> !AJH fixed bug next line
> Bio(i,k,iSzoo)=Bio(i,k,iSzoo)/(1.0_r8+cff)
699,701c704,706
< ExcZS=(AlphaZS(ng)-BetaZS(ng))*ZS2PS
< ExcZL=(AlphaZL(ng)-BetaZL(ng))*(ZL2PS+ZL2PL+ZL2ZS)
< ExcZP=(AlphaZP(ng)-BetaZP(ng))*(ZP2PL+ZP2ZS+ZP2ZL)
---
> ExcZS=(AlphaZS(ng)-BetaZS(ng))*GraPS2ZS
> ExcZL=(AlphaZL(ng)-BetaZL(ng))*(GraPS2ZL+GraPL2ZL+GraZS2ZL)
> ExcZP=(AlphaZP(ng)-BetaZP(ng))*(GraPL2ZP+GraZS2ZP+GraZL2ZP)
738a744,746
> !AJH added the next line
> fac5=dtdays*VO2S0(ng)
>

------------------------------------------------

[wilkin]

Al,

I'm not familiar with the details of NEMURO but discussed several of the issues with Hernan when he was coding the option.

(We believe the implicit time stepping introduces accuracy in the time stepping that was wholly missing from the generic NEMURO code, so that in itself is a major advance.)

With regard to your question, I tested the bio_sediment option in the fasham option and it conserved perfectly. It is in fact rather more sophisticated, allowing for partition of the net sinking flux of nitrogen into ammonia, nitrate, and nitrogen gas.

You should possibly look at whether the currency of the Nemuro state variables are consistent with the bio_sediment parameterization. Namely, are the units of sinking opal actually in terms of nitrogen concentration, such that the sinking flux can be directly converted back to SiOH in the code lines:

cff1=FC(i,0)*Hz_inv(i,1)
Bio(i,1,iSiOH)=Bio(i,1,iSiOH)+cff1

Find the right biologist and its probably a straightforward question, but a total mystery to me.

When this is resolved, you might want to reconsider whether a parameterized microbial loop that instantly converts PON into nitrate and opal into SiOH is even rational (?!)

John.

[hermann]

John,

Thanks for your input; I certainly agree that the ROMS implementation is a big advance in terms of numerics (indeed, it conserves extremely well with sinking turned off).

The nitrogen and silicate parts of NEMURO are fundamentally separate, and one should not affect the other as regards conservation.

My question is actually not so much about bio_sediment (an admittedly crude way to close the nutrient cycle), but the semi-Lagrangian sinking algorithm throughout the water column, which I thought might be nonconserving. Your statement that the fasham option conserved perfectly with bio_sediment enabled suggests that the sinking algorithm should not be the source of my problem. I'll go back to the code and see what I can find......

I might try running with BioIter set to a number greater than one, as well, in case there is some overshooting going on with the sinking - is such a CFL-like error even possible, given what you know about the sinking algorithm?

-Al

[wilkin]

Al,

You prompted me to read a bit further and I see now (nicely illustrated in the diagram in Hernan's first post to this thread) that the silica cycle is separate and does carry any of the total nitrogen inventory - the currencies are separate.

To take sinking out of the list of suspects, you could try a very deep BIO_TOY set-up so that the BIO_SEDIMENT option never does anything, and then selectively disable various terms with input parameters that drive them to zero.

Does the imbalance scale with small perturbations to the fall velocity?

You are clearly doing the beta-tester duty that Hernan called for in his post.

Thanks, John.

[hermann]

Well, it turns out the BIO_SEDIMENT section was the source of the bug after all - the undefined variable "indx" was used where "ibio" should be! Things look MUCH better now, total nitrogen is conserving well. Here's the latest diff:

-----------------------------------------

/home/users/hermann/ROMS2007/ROMS/Nonlinear/ 5 % diff nemuro-original.h nemuro-modified.h
367c367
< RnewS=GppNPS/GppPS
---
> RnewS=GppNPS/(GppPS+MinVal)
404c404,405
< RnewL=cff4/(cff4+cff5)
---
> !AJH added MinVal in the next line
> RnewL=cff4/(cff4+cff5+MinVal)
406,407c407,411
< cff4=cff6*RnewL
< cff5=cff6*(1.0_r8-RnewL)
---
>
> !AJH added normalization in next two lines
> cff4=cff6*RnewL/(Bio(i,k,iNO3_)+MinVal)
> cff5=cff6*(1.0_r8-RnewL)/(Bio(i,k,iNH4_)+MinVal)
>
415c419
< Bio(i,k,iLphy)=Bio(i,k,iSphy)+GppPL
---
> Bio(i,k,iLphy)=Bio(i,k,iLphy)+GppPL
423c427
< RnewL=GppNPL/GppPL
---
> RnewL=GppNPL/(GppPL+MinVal)
535c539
< cff4=1.0_r8-EXP(LamS(ng)*(PS2ZSstar(ng)-Bio(i,k,iSphy))
---
> cff4=1.0_r8-EXP(LamS(ng)*(PS2ZSstar(ng)-Bio(i,k,iSphy)))
556c560
< cff4=1.0_r8-EXP(LamL(ng)*(PS2ZLstar(ng)-Bio(i,k,iSphy))
---
> cff4=1.0_r8-EXP(LamL(ng)*(PS2ZLstar(ng)-Bio(i,k,iSphy)))
611c615,616
< Bio(i,k,iSzoo)=Bio(i,k,iSzoo)*(1.0_r8+cff)
---
> !AJH fixed bug next line
> Bio(i,k,iSzoo)=Bio(i,k,iSzoo)/(1.0_r8+cff)
699,701c704,706
< ExcZS=(AlphaZS(ng)-BetaZS(ng))*ZS2PS
< ExcZL=(AlphaZL(ng)-BetaZL(ng))*(ZL2PS+ZL2PL+ZL2ZS)
< ExcZP=(AlphaZP(ng)-BetaZP(ng))*(ZP2PL+ZP2ZS+ZP2ZL)
---
> ExcZS=(AlphaZS(ng)-BetaZS(ng))*GraPS2ZS
> ExcZL=(AlphaZL(ng)-BetaZL(ng))*(GraPS2ZL+GraPL2ZL+GraZS2ZL)
> ExcZP=(AlphaZP(ng)-BetaZP(ng))*(GraPL2ZP+GraZS2ZP+GraZL2ZP)
738a744,746
> !AJH added the next line
> fac5=dtdays*VO2S0(ng)
>
970c978,980
< IF (indx.eq.iPON_) THEN
---
>
> !AJH fixed bug "indx" should be "ibio"
> IF (ibio.eq.iPON_) THEN
975c985
< ELSE IF (indx.eq.iopal) THEN
---
> ELSE IF (ibio.eq.iopal) THEN

[jlz]

I have run NEMURO with both IVLEV_EXPLICIT and HOLLING_GRAZING without problems, but I got lots of NaN with IVLEV_IMPLICIT. I wonder if it is because of terms like cff5=1.0_r8/(fac2*cff4). If it is, how to fix it?
Thanks,
Jinlun Zhang

[twainwright]

Hi,

I'm just starting to run ROMS NEMURO inside bio_toy to compare with my own NEMURO code. I found that I could run the HOLLING_GRAZING option OK, but IVLEV_EXPLICIT results in the following cpp error:

/usr/bin/cpp -P -traditional -DLINUX -DX86_64 -DGFORTRAN -D'ROOT_DIR="/home/twain/Code/ROMS"' -DBIO_TOY -D'HEADER="bio_toy.h"' -D'ROMS_HEADER="/home/twain/Code/MyROMS/test_Biotoy/Include/bio_toy.h"' -DNEMURO -DNestedGrids=1 -D'ANALYTICAL_DIR="/home/twain/Code/MyROMS/test_Biotoy/Functionals"' -D'MY_ANALYTICAL="on"' -D'SVN_REV="234M"' -IROMS/Include -IROMS/Nonlinear -IROMS/Utility -IROMS/Drivers -I/home/twain/Code/MyROMS/test_Biotoy/Functionals -IROMS/Functionals -I/home/twain/Code/MyROMS/test_Biotoy/Include -IMaster -ICompilers -D'HEADER_DIR="/home/twain/Code/MyROMS/test_Biotoy/Include"' -DNEMURO ROMS/Utility/checkvars.F > /home/twain/Code/MyROMS/test_Biotoy/Build/checkvars.f90
ROMS/Utility/checkdefs.F error: operator '&&' has no right operand

This seems to be a trivial typo in checkdefs.F. The following change fixed the problem--compiles and runs OK:

$ diff ROMS/Utility/checkdefs.F.orig ROMS/Utility/checkdefs.F
2891c2891
< #if defined NEMURO && defined HOLLING_GRAZING && IVLEV_EXPLICIT
---
> #if defined NEMURO && defined HOLLING_GRAZING && defined IVLEV_EXPLICIT

Nobody else seems to have hit this problem, maybe my cpp (GNU GCC 4.3.0) is pickier than others?

--Tom