EAMxx: calculate vertical derivatives in diagnostics

[mahf708]

[BFB]

--

this diagnostics is desired for the CESS2 campaign (cc HB+CT). I am requesting reviews from dycore experts to double check this conforms with how the 'divergence' is calculated there. I opted to generalize it just in case someone wants a value other than d(omega)/d(p)... also, DRAFT for now until one of us validates this is actually what we want (maybe in dpxx)

[github-actions]

PR Preview Action v1.6.2
🚀 View preview at https://E3SM-Project.github.io/E3SM/pr-preview/pr-7491/
Built to branch gh-pages at 2025-08-14 22:52 UTC. Preview will be ready when the GitHub Pages deployment is complete.

[bartgol]

I am requesting reviews from dycore experts to double check this conforms with how the 'divergence' is calculated there.

We don't have d/dz in the dycore, only d/dp, IIRC. And for that case, the formula seems to be the same as what you have here.

[mahf708]

@bartgol, I haven't thought deeply about this entire block. I just wanted to put in a skeleton for the dz counterpart of dp derivative (it will be needed eventually, but the main focus is dp for now). I think this should work, though I am not sure about recycling the same zm_icol and doing away with team barriers (in the surface coupling there are both fields and barriers). Let me know what you think!

[bartgol]

You do need the barriers, since you cannot compute z_mid until z_ int is done and you can't compute z_int until dz has been computed (there may be a race condition). OTOH, you can use the same view for both dz and zm, since the barriers ensure that dz has been used completely before you start computing zm.

Edit: I noticed you use the same vector for zmid and zint. You cannot do that. You need two separate fields.

[mahf708]

Update: I'd like to go back to usin dp and dz directly. I will need to double check the alignment. Then, I will have this ready for a merge

[mahf708]

@bartgol @whannah1 I'd like your opinion on renaming these to:

      # in this example, we use T_mid in units of K
      - T_mid_pvert_derivative  # K / Pa
      - T_mid_zvert_derivative  # K / m

Basically, getting rid of the "d" (redundant with "derivative") and attaching the p/z to vert since it is fundamentally part of that and to better distinguish it from other things.

Questions:

1. What do you think?
2. Should I shorten pvert to pver? (leaning yes)
3. Should I shorten derivative to deriv? (leaning no)

[bartgol]

Seems reasonable. I don't think there is a beautiful string that would make everyone satisfied anyways, so I think your suggestion is as good as it gets. Perhaps we could shorten derivative->deriv, since piping diags can make for veeeery long names already, so we may as well save 5 chars :)

[whannah1]

I think I prefer the verbose version - I also prefer "vert" to "ver"

[mahf708]

This is ready for final review and merge.

[bartgol]

I think it would be more efficient to make this diag depend on the dz diag (for "z" deriv method). And in both cases (z and p), make the denom alias the dz/dp input, and avoid update calls.

I also have a thought on testing, but I don't have much momentum on that topic, considering we've been doing "code repetition testing" for a while...

[bartgol]

I am not sure if this way of testing our diags (which, tbh, is a widespread issue, so it's not a critique to this PR per se) is great. I think it is only giving us a false sense of security. We are basically replicating the same piece of code twice, and verifying it gives the right answer. So if by chance we use a bad formula, all we are testing is that we used the same formula twice...

I would like more property based tests or analytical solution tests. Some examples:

• f(k+1)>f(k) then f'>0 (and viceversa)
• f(k+1) = a*f(k), dp(k+1)=b*dp(k) with a,b>0, then f'>1 if a>b, f'<1 otherwise (I think we should get something like f' ~ a/b)
• set f_mid(k) = g(p_mid(k)), and verify that f'(k) ~ g'(p_mid(k)). This of course holds if dp_mid(k)<<1. So for instance, with p(k) =k/1000, f(k)=e^p_mid, f' ~ e^p_mid (the computed value should be e^p_mid * (exp(0.001)-1)/0.001, but the extra factor is ~1. Or maybe something with g(p)=sin(p), or a polynomial.

[mahf708]

@bartgol addressed all your comments except the testing ... I will be adding these to a cime test soon, and yes, we should think about better testing diags. There's some utility in doing the tests we do, because they do find common coding mistakes, so the code "works" (but ofc we don't know if it works like intended)

[bartgol]

Up to you, but you could also consider recycling some existing functions.

parallel_for ("", policy, KOKKOS_LAMBDA(const MT& team) {
   // get cols subview
   
   // Compute f at interfaces
   ColumnOps::compute_interface_values_linear(team,nlevs,f_icol,dpicol,f_icol(0),f_icol(nlevs-1),temp_int);
   team.team_barrier();
   
   // Compute deriv numerator
   ColumnOps::compute_midpoint_deltas(team,ncol,temp_int,o_icol);
});
get_diagnostic().update<CombineOps::Divide>(denominator,1,1);

This would require a temporary field (could also use workspace, but that's not that important imho). You could also do the division in the kernel, but we don't have a ColumnOps::combine util (we may add it though).

Another important diff would be the boundary treatment. Your impl uses fwd/bwd euler at boundary, but the denominator is not really the right one (if you use fwd/bwd euler, you should use dp_int(1) and dp_int(nlevs-1), as dp(0) does not measure the thickness between 1st/2nd midpoints (and similarly for dp_int(nlevs-1)). The impl I show here is what rrtmgp uses to interpolate T_mid->T_int. At the boundaries, we use f_int(0)=f_mid(0) and f_int(nlevs)=f_mid(nlevs-1).

[bartgol]

I would encourage using existing tools/routines for interpolation, but besides that, I think the BC you have is not quite right. The reason is that the denominator is not really measuring the distance between the 1st and 2nd midpoints. You'd prob need to divide by dp_int(1), which needs to be computed from pmid(1)-pmid(0) or something like that.

Alternatively, I think one can use a constant extrap (f_int(0)=f_mid(0)) or maybe linear (which can be done as f_int(0) = 2f_mid(0) - f_int(1)).

[mahf708]

For the record, @hassanbeydoun sent me this validating the scientific impl of this derivative:

• div = omega_pvert_derivative_where_qi_gt_0.0001_horiz_avg * count_where_qi_gt_0.0001_horiz_avg
• csc = omega_pvert_derivative_where_qi_le_0.0001_horiz_avg*(1-count_where_qi_gt_0.0001_horiz_avg)

the divergence theorem stipulates the intergals of those two must match, which is evident in the plot