clarification on prt_mass diagnostic function

[aschuh]

I posed this question to the UFS-fv3atm github repos but was pointed to this one as possibly being more appropriate.

#########

I was hoping to get clarification on the prt_mass and z_sum subroutines in the FV3/atmos_cubed_sphere/tools/fv_diagnostics.F90 file. I'm dumping emissions of a passive tracer into the q tracer array and hoping to sum up the atmos burden to verify that everything is working.

First, I'm confused by the is, ie, js, and je. I'm assuming they are tied to domain decomp and MPI and hence should try to "leave them alone" as much as possible and mimic existing code. Nevertheless, I'm also confused by the argument "n_g" and use the back and forth use of different subarrays of the q tracer array. Some sort of buffer or offset?

For example, q and psq are allocated early as:

real, intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km, nq)
real psq(is:ie,js:je,nwat)

But these are often fed into "z_sum" subroutine as:

if (liq_wat > 0) &
call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,liq_wat), psq(is,js,liq_wat))

In this case, q() basically becomes a scalar input. However, it is also assigned from is-n_g,js-n_g position to is,js position in psq. Furthermore in the z_sum function, the expected input is a non-trivial array (see q() below):

subroutine z_sum(is, ie, js, je, km, n_g, delp, q, sum2)
integer, intent(in):: is, ie, js, je, n_g, km
real, intent(in):: delp(is-n_g:ie+n_g, js-n_g:je+n_g, km)
real, intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)

This back and forth between q(is-n_g:ie+n_g, js-n_g:je+n_g, km) and q(is-n_g, js-n_g, km) and inconsistency in expected args makes me concerned that it isn't summing my tracers right (which it doesn't look to be).

Furthermore, it often appears that the tracer array is often only being summed over the 1st vertical level (unless I'm missing something?). It seems like it should be q(is-n_g,js-n_g,:,liq_wat) instead of q(is-n_g,js-n_g,1,liq_wat) if you want ALL liq water summed.

if (liq_wat > 0) &
call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,liq_wat), psq(is,js,liq_wat))

Basically, a lot of confusing stuff here and I can't follow the domain decomp/MPI stuff well enough to understand the indexing that is going on.

Any help or advice would be appreciate, thanks.
andrew

[lharris4]

Hi @aschuh . This is an excellent question that gets down to how FMS (used by FV3 and MOM6) divides a large domain amongst processor cores, and also some of the subtleties of Fortran.

The whole horizontal domain, which on the cubed sphere is (npx-1) x (npy-1) x 6 grid cells, is decomposed into a number of subdomains, each of which has a "compute" extent of isc:iec,jsc:jec on which the solution is computed upon. Often the 'c' is dropped from these names.

Each processor's domain also has a halo of "ghost" grid cells copied from the adjacent subdomain acting as lateral boundary conditions; in FV3 the halo has width 3, which is called n_g. The compute domain plus halo is called the "data" domain, which has extent isd:ied,jsd:jed. (Note that this differs a bit for staggered variables, like u, v, divg, etc.) Usually is-n_g is the same as isd, ie-n_g the same as ied, etc, which is also used in parts of the code.

In FV3, an entire column is kept on the same processor; we do not decompose in the vertical, although we do use multi-threading in the vertical. Also note that the number of vertical levels may be referred to as either km or npz; these are always equal.

As for the call to z_sum(), Fortran will allow you to reference an array in a subroutine call by specifying just its first element in each dimension, which is what is done when referring to q(is-n_g,js-n_g,1). Since the other arguments specify the size of the full array the subroutine needs, Fortran can read off all of the other elements in memory without making a new copy of the data. This is then like C's "pass-by-reference" technique.

This admittedly tricky syntax allows for a more flexible routine and also faster memory access, at the cost of a bit of readability; it appears a lot in FV3. You could call it with explicit array bounds instead (q(is-n_g:ie+n_g,js-n_g:je+n_g,1:km)), but it could be slightly slower. You can also just call with q (no subscripts) if you mean to pass the whole array, but do check the array bounds to make sure you intend to send the whole array.

There are a couple of tutorials that explain FMS and its domain decompositions more completely. @bensonr and @laurenchilutti may be able to point you to these.

Thanks,
Lucas

[aschuh]

Thanks Lucas. I've spoken to a few other NOAA folks who use the routine and I imagine I just need to move forward with the understanding that it "works" despite some of the coding nuances. I'm having some trouble verifying emission inputs vs atmos concentrations. I've been told that the NEXUS/ESMF regrid to cubed sphere is mass conserving and the prt_mass should correctly calculate atmos burden relative to wet mixing ratios (delp), meaning I'll have to look somewhere else for an explanation ,thanks again.