Questions about MARBL code time-stepping and vertical integration

mangelesg commented 3 years ago

Hello,

Im trying to implement MARBL in python with a similar code, however I have two questions about the implementation:

1- When MARBL is coupled to a GCM: which time step of the tracer should be passed on to MARBL? --- In the MARBL documentation it says that the current time step of a tracer "T" should be passed: T(current). --- However, looking at the BGC code from CESM2, the tracer is being passed at 0.5*(T(current) + T(previous)). Which one would be the correct form?

2- In the MARBL code: when the ecosystem equations are being solved for a given vertical level k=km; why do the particulates terms need to be loop over the vertical levels from k = 1, km?

Here is how it looks in the code:

some previous code ......

call compute_routing(km, zooplankton_derived_terms, autotroph_derived_terms)

call compute_dissolved_organic_matter (km, marbl_tracer_indices, num_PAR_subcols, &
     PAR, zooplankton_derived_terms, autotroph_derived_terms,   &
     delta_z1, tracer_local(:, :), dissolved_organic_matter)

here is where the sub-loop for particulate terms starts

do k = 1, km

   call compute_scavenging(k, km, marbl_tracer_indices, tracer_local(:,:), &
        POC, P_CaCO3, P_SiO2, dust, Fefree(:), Fe_scavenge_rate(:), &
        Fe_scavenge(:), Lig_scavenge(:), marbl_status_log)

........ some more code

end do ! k

..... more code

My question is: Wouldn't be enough to loop only one time at each z-level, solve the ecosystem equations and save Particulate_FLUX_OUT to use as Particulate_FLUX_IN for the next vertical level, without having to do an extra loop for particulates?

I would appreciate very much any help in these matters,

Angeles

mnlevy1981 commented 3 years ago

1- When MARBL is coupled to a GCM: which time step of the tracer should be passed on to MARBL? --- In the MARBL documentation it says that the current time step of a tracer "T" should be passed: T(current). --- However, looking at the BGC code from CESM2, the tracer is being passed at 0.5*(T(current) + T(previous)). Which one would be the correct form?

I believe this should actually depend on your GCM's time-stepping method. My understanding is that CESM used 0.5*(T(current) + T(previous)) because it had a leap-frog time-stepping scheme (though I don't know why that formulation remains in place with the move to the Robert filter).

2- In the MARBL code: when the ecosystem equations are being solved for a given vertical level k=km; why do the particulates terms need to be loop over the vertical levels from k = 1, km?

This is being addressed in #53 - the code MARBL was originally based on computed the interior tendency terms at all (lat,lon) grid points on the block for a given level before moving on to the next level. We did not want MARBL to need to know anything about the horizontal structure of the GCM grid, so we tranposed everything to be computed column-by-column instead of level-by-level. This included pushing the k loop into routines like compute_routing() and compute_dissolved_organic_matter(), but compute_scavenging(), compute_large_detritus_prod(), and compute_particulate_terms() are in the loop you highlighted because there is still some interdependency between those functions from level to level. I think your comment

Wouldn't be enough to loop only one time at each z-level, solve the ecosystem equations and save Particulate_FLUX_OUT to use as Particulate_FLUX_IN for the next vertical level, without having to do an extra loop for particulates?

is similar to the approach we plan to take when we finish addressing the issue - we need to rethink some of the data structures but should be able to remove that loop soon.

mangelesg commented 3 years ago

Once more, thanks a lot for your prompt and clear response.

mangelesg commented 3 years ago

Hello again,

Im having trouble implementing MARBL, and I would like to know if you have any advice on the following issue.

First, some context:

I have an advection-diffussion offline model for a small region, that borrows velocities from a parent model.
The advection-diffussion code works well for a single tracer.
Then, I have the MARBL code, which I coupled with the offline model: I advect and diffuse the tracers and apply the MARBL equations. I tested the coupled model with the file provided by MARBL "call_compute_subroutines.history.nc", to look at the initial biology-tendency terms, and I get very similar values.
The problem is when I run the model for several days. Then, in the first vertical level the nitrate starts going to zero (with some negative values) very fast, while phytoplankton grows very fast too (mostly small phyto.). I only run the coupled model for 120 days with a time step of 800 seconds.

I understand that it may be very hard to help me without having all the details and the code, but I thought maybe this issue is common and I could get some advice. In any case, here are my questions:

--- Could it be that the ecosystem has not reached equilibrium yet? Is there a reasonable amount of integration days where the system reaches an equilibrium?

--- When MARBL is implemented with the adv-diff model, should I apply a condition of minimum value for nutrients? or a saturation value for phytoplankton?

---Since I obtain the correct initial tendency terms from MARBL, does that mean that my issue has to be in the time stepping? Does MARBL require a filtering process?

I would very much appreciate any help

mangelesg commented 3 years ago

Actually, I set a minimum limit for all the nutrients of 1e-9, and that seemed to solve the problem!

marbl-ecosys / MARBL