EAMxx: PBL entrainment budget diags

[mahf708]

Adds planetary boundary layer entrainment budget diagnostics as in mülmenstädt ... mahfouz ... et al (2024), ACP, doi: 10.5194/acp-24-13633-2024.

[BFB]

[mahf708]

@bartgol @tcclevenger: There are two elements that I'd like to address before merging (or soon thereafter):

• the one marked as TODO (returning multiple accumulations from the parallel_reduce); I wrote this code before we had kokkos 4.2 in, so I am just reviving it as-is
• A lot of the manual labor with the precalculated fields could be automated:
  • the atm_backtend could be used instead of doing stuff here
  • I also added theta_l (liquid potential temperature) that could also be used here instead of doing a manual calculation

Because the above require some thinking, I am punting at least for a bit. I wanted to get this on your radar for any potential feedback before digging deeper into this. Obviously, another big one that BIG TODOs on my list is to allow these classes to return individual fields rather than one big composited field...

[bartgol]

@mahf708 how urgent is it to have this diag in master? It seems to me that this is yet another case (along with water/number path, aodvis, precip_surf_mass_flux, etc) where supporting multiple-output diagnostics would be quite beneficial. Given the impl, I would say it may be even more beneficial, since it's actually computing all diags and stuffing them in a single field, which may be confusing.

So, if having this diag is not SUPER urgent, I would encourage you to work on the multi-diags support first (assuming you're still interested in working on it, as you previously said).

[mahf708]

@mahf708 how urgent is it to have this diag in master? It seems to me that this is yet another case (along with water/number path, aodvis, precip_surf_mass_flux, etc) where supporting multiple-output diagnostics would be quite beneficial. Given the impl, I would say it may be even more beneficial, since it's actually computing all diags and stuffing them in a single field, which may be confusing.

So, if having this diag is not SUPER urgent, I would encourage you to work on the multi-diags support first (assuming you're still interested in working on it, as you previously said).

Not urgent. I am just starting to (scientifically) validate this in the sense of ensuring I trust that these (relatively new) diagnostics are outputting meaningful results that I can trust, so I wanted to make sure the code isn't messed up.

And yeah, I am actively working on the multi-out diag too btw; the most urgent diags to convert are aerocom_cld and pbl_entrainment_budget (this one) becaues they are hacky.

What I would appreciate from you and @tcclevenger is if you have a moment to skim through it, please do, and let me know if you spot anything sus. Otherwise, no urgency at all. We have a lot of time to get this through ...

[rljacob]

BFB or non-BFB for eamxx?

[bartgol]

I have some questions, mostly b/c I don't know what the PBL is and had to google most of this stuff :)

[bartgol]

It doesn't change much, but I like to rely on the grid for the layouts, since it hides some implementation detail about how they are created:

auto scalar2d_layout = grid->get_2d_scalar_layout();
auto vector2d_layout = grid->get_2d_vector_layout(m_ndiag);

[bartgol]

You can just do

m_prev_qt = qv.clone("qtot_prev");

Similarly for m_prev_tl.

[bartgol]

This allocates memory at every evaluation. You could move the creation of the wsm to initialize_impl.

[bartgol]

do you mean "temperature gradient inversion"? Cause a negative T would be a bit concerning...

[bartgol]

Looking at the code below, you are looking for the minimum gradient. From your comment, I thought you wanted to find the place where the inversion starts. Is the comment misleading? Or is the code wrong?

[mahf708]

This should read temperature change/gradient goes positive/negative, you're right.

        // (for now). First, we can simply find the first level from the surface
-       // that has a a temperature "inversion" (temperature goes positive
+       // that has a temperature "inversion" (temperature gradient goes positive
        // instead of negative). Second, we can find the level which has the

I think the definition is about the sign of the gradient (so like the corrected comment) but in practice, it likely doesn not matter (because the minimum gradient will likely happen in the first layer this inversion occurs)...

I need to double check this. I don't actually know much more than you about PBL tbh... PBL stands for planetary boundary layer

[bartgol]

This can also happen if the ||reduce above finds NO suitable levels. In that case, minloc will remain with the initial values, which is an identify for the min operation. As one would expect, minloc.val would contain a large number (DBL_MAX). Luckily, minloc.loc is also initialized to a large number (INT_MAX), so the check would indeed work.

Still, I would consider adding a comment specifying we hit this branch if NO level is >70kPa.

[bartgol]

You never check that tm_grad(k)<0. I guess this usually happens at some point in the v direction, but as a software dev with little atm modeling knowlege I wonder "what if it never occurs?"

[mahf708]

I am confused. Why wouldn't it check tm_grad(k)<0? The other conditional is on pm_icol, not tm_grad. I definitely want this to check for <0; am I messing up?

[bartgol]

When using minloc_value_t in a parallel reduce, Kokkos initializes the val member to "something really large" (it uses DBL_MAX). That way, at the first iteration, whatever value you get, it WILL be smaller than what's already stored. However, nothing in youre code requires that tm_grad<0. If you had tm_grad>0 for all k, the code would still work. But it would not find a level where the grad is negative.

[bartgol]

If you are only interested in the "last level (from the top) where grad(T) is neg", you don't need MinLoc. You can just do

        int switch_lev = -1;
        Kokkos::parallel_reduce(
            Kokkos::TeamVectorRange(team, 1, num_levs),
            [&](const int &k, int &result) {
              if(tm_grad(k)<0)
                result = max(result,k);
            },Kokkos::Max<int>(switch_lev));

[mahf708]

Ah, I see what you mean!

So, I decided to go this route for generality. The other two methods (qt, theta_l) explicitly look for the "biggest jump" (and these two jumps have different signs):

        // We first want to find the PBL inversion. There are three methods to
        // do so. All our methods here rely on the *gradient* of state fields
        // (for now). First, we can simply find the first level from the surface
        // that has a temperature "inversion" (temperature gradient goes positive
        // instead of negative). Second, we can find the level which has the
        // biggest positive jump in theta_l. Third, we can find the level which
        // has the biggest negative jump in qt.

I misinterpreted what you meant as it never goes below <0; now, rereading it, I see that you meant, "you don't ensure it is negative" --- yes, I don't do that because I am looking for the biggest jump (and the sign of the jump depends on which method we are using).

Amazing hawk eye catching all these minor details; I appreciate it :) this will help me improve the PR once validated!

[bartgol]

Given what discussed above, the integral now starts from the point where tm_grad is the most negative. Is that the idea? Or was it supposed to start from the point where tm_grad changes sign?

[bartgol]

The hardcoded indices are a bit obscure. Perhaps consider using some variables and add comments, so we know why:

// In the input iop dataset the inversion is just below level 123 (recall that lev=0 is at the top)
constexpr int expected_inv_idx = 123;
REQUIRE (out_hv(0,0) == Real(iop_pm[expected_inv_idx]));
...
for (int ilev = expected_inv_idx+1; ilev<nlevs; ++ilev) {
  ...

[mahf708]

Consensus: Let's work on merging #6935 first. The other PR will take a bit of careful work, so it may be a while. I will either keep this up or close it depending on how fast progress goes on the other PR.