Make field manager use wet mixing ratios, physics procs use dry mixing ratios as needed

[PeterCaldwell]

Updated plan: we actually decided the field manager should always use wet mixing ratios because this makes more sense for weird variables like TKE that would be totally counterintuitive if they were per unit dry air. We want to have the field manager always do the same thing to keep code (and assumptions) simple. Most physics processes make more sense for dry mixing ratios, however. We will handle this by creating universal functions for converting between dry and wet and using these functions liberally in process interfaces.

OLD TEXT BELOW:

@bogensch , @brhillman , and @PeterCaldwell agree that it would be better to use mixing ratio (water mass divided by mass of dry air) rather than specific humidity (water mass divided by total mass). One reason for this is that a lot of empirical formulations used in P3 microphysics are based on mixing ratio rather than specific humidity. Another reason is that the variable q typically refers to a mixing ratio rather than specific humidity, so it is really confusing to use q as the core variable for moisture in our model and have it refer to specific humidity. Finally, mixing ratio is the more typical variable in the atm sci world.

Note, however, that this is different from what EAM does - it assumes q is a specific humidity.

P3 and SHOC are already formulated assuming mixing ratio. In particular, our qsat (saturation calculation) already uses mixing ratio. RRTMGP uses total cell mass (qv*pres) as its water variable and doesn't have the interface to qv set up yet - so will need to perform the correct conversion.

My understanding is that HOMME operates on total mass qv*pres and converts to specific humidity for passing to physics. That is probably the only concrete task for this issue.

[PeterCaldwell]

@mt5555 - I'm assigning you because I know you probably understand the issue. Feel free to delegate to someone else. @bartgol in particular maybe was the person who converted qdp to qv in hommexx so may be able to quickly fix this.

[PeterCaldwell]

Also - I talk about qv above just for concreteness and convenience. I suspect this issue will also affect qc, qi, nc, ni, etc.

[oksanaguba]

+1 for 'dry' water forms, would make relaxing const moist pressure approximation easier.

[mt5555]

You wrote SCREAM v0.1 above - but I think you mean SCREAM v1?

some initial comments:

For v0.1, using the EAM fortran code base, I dont think we want to make such a change, as we'd also have to upgrade the energy and diagnostics to make sure they can work with either wet or dry mixing ratios of water variables. I'm assuming right now they are all coded for wet mixing ratios. For v0.1, I think it would be better to have P3 and SHOCK interface layers convert to dry mixing ratios before calling the parameterization (if not being done already).

In SCREAM v1, I would propose the coupling layer communicates tracer mass - as everyone agrees on this independent of mixing ratios and amount of water loading. Also, for the np4/pg2 mapping, we want to conserve mass and be monotone w.r.t. wet mixing ratio, so I think it makes most sense to convert to dry mixing ratios at the end of the coupling.

[PeterCaldwell]

Great point @mt5555 . I was thinking we would make this change for v0.1 as well, but you're right that doing so requires digging into the spaghetti code we're trying to avoid. So we shouldn't. Converting from specific humidity to mixing ratio in the v0.1 interfaces is probably the way to go... though since we're doing things differently in C++, this will introduce non-bfb differences between the F90 and C++ versions of the interfaces.

In terms of a v1 plan: I'm not sure what you mean by "coupling" here. My assumption was that the dycore would do what it's always done (by handling qdp) and the dycore's interface layer would translate qdp into dry mixing ratio for storage in the field manager and communication to the rest of the model. Is this what you're advocating for, or did you want the field manager to store qdp and all the processes need to convert that to mixing ratio?

[mt5555]

By coupling layer - I think I should have said the "driver". I was assuming the SCREAM driver is similar to the EAM "dp coupling layer". All this code currently lives in the dp coupling layer, I think it would be most natural to port it into the new driver:

1. extract data from HOMME data structures (Qdp)
2. map np4->PG2. (conservative, monotone w.r.t. wet mixing ratio)
3. compute dry mixing ratios on PG2 grid. (depends on water loading choices)
4. store in field manager

Then after the physics, the driver would invert that procedure.

[PeterCaldwell]

Ok, thanks Mark. That makes sense.

Just FYI and to keep things confusing, we're referring to the dp coupling layer as the dynamics interface layer these days (since all processes have an interface layer and originally the mapping between grids was supposed to occur outside the interface layer).

[mt5555]

Sorry I'm not up on the status of the dynamics interface layer. Is this done already? Or should we start a confluence page to flush out a design? We should stick with mass in this layer, and then the conversion to dry mixing ratios should just be a subroutine call at the very end. The more complex code from d_p_coupling.F90 that will have to be ported is either the mapping from GLL duplicate degrees of freedom to the physics "uniquepoints" decomposition, or the GLL->PG2 mapping code (for when we do PG2).

[PeterCaldwell]

My understanding is that @bartgol is working on the dynamics interface layer right now and that the mapping stuff was done previously when we thought mapping should be independent of process. We aren't supporting PG2 for our first implementation of SCREAMv1, which simplifies things. I think converting to mixing ratio rather than specific humidity at the end of the dynamics interface layer should be trivial to add. @bartgol - I'm not sure if this makes sense. I can fill you in on slack if you haven't read this whole conversation.

[bartgol]

Sorry, I get a bit confused when you say stuff like 'convert to wet mixing ratio', cause I'm confused which of the two ends uses wet and which one uses dry.

So, for the sake of clarity, the only info I need is: what is stored inside the FM, what each process is expecting, and what is the conversion formula.

Digging in eam I found q_dry = q_wet*(dp_wet/dp_dry), and dp_dry = dp_wet(1-qv). The question now is: which version of dp are we storing in the FM?

Also, when converting dp from wet to dry, is qv supposed to be wet or dry? Assuming a 'dry' qv makes sense at all... If the total mass is dry, can qv_dry be larger than 1?

[mt5555]

Per @PeterCaldwell 's message, FM will be storing pseudo_density=dp_wet, and all tracers as dry mixing ratios.

The dp_dry calculation above should be a subroutine, since that formula depends on water loading (which we hope to add in the future).

formula, removing ambiguity about qv: dp_dry=dp_wet(1-qv_wet)=dp_wet/(1+qv_dry)

[bartgol]

And homme uses dry or wet tracers? Also, is the conversion wet<->dry to be applied to all tracers? Or is it just for some of them? E.g., I know ni and qi are different kind of ratios; is the conversion formula the same? What about other tracers (meaning other than q and n?

[mt5555]

Yes, the wet/dry conversion can be applied to any tracer.

Regarding HOMME: internally HOMME uses mass (not mixing ratios) as the prognostic variable, but does carry around (for convenience) the wet mixing ratios of all tracers.

[bartgol]

Right, Homme works with qdp. I guess my question was "how should I relate Homme's qdp with q and dp?". From what you said, I am guessing qdp=q_wet*dp_wet?

[mt5555]

the fortran dynamics interface layer communicates via tendencies for wet mixing ratios. Simplest would probably just have the dyn interface layer convert to wet mixing ratios, and then the rest of it can be a straight port of the fortran version. The fortran version will convert to a wet mixing ratio tendency, map to GLL data decomposition (possible including PG2->NP4), and then convert to a mass tendency, then call routines in HOMME to apply these tendencies in various ways based on "ftype".

[PeterCaldwell]

Yes, the wet/dry conversion can be applied to any tracer.

We need to be a bit careful here - the whole concept of wet vs dry pressure only applies to tracers with units of stuff per kg air. Thus we wouldn't want to apply this conversion to TKE, for example.

[PeterCaldwell]

Regarding HOMME: internally HOMME uses mass (not mixing ratios) as the prognostic variable, but does carry around (for convenience) the wet mixing ratios of all tracers.

This confuses me. I thought the q in qdp was the big Q = all tracers lumped into a single array. I'm surprised that wet mixing ratios are also carried around in HOMME. Does this require us to convert from dry to wet mixing ratios in the dynamics interface before passing variables into the actual dynamics code? Or can we delete this redundant carrying of both qdp and q?

[bartgol]

Regardless of whether Homme keeps q around or not, don't we need to convert dry->wet anyways, before calling Homme? I mean, the (backed out) forcing FQ is dry on the AD side (since Q is dry in the Field manager), but it has to be wet in Homme, so the dyn interface would need a dry->wet conversion (and back, after Homme returns) anyways, no?

[bartgol]

Yes, the wet/dry conversion can be applied to any tracer.

We need to be a bit careful here - the whole concept of wet vs dry pressure only applies to tracers with units of stuff per kg air. Thus we wouldn't want to apply this conversion to TKE, for example.

That was my question: I was wondering about stuff like radiation tracers, or even ni, nc, nblah (for instance).

Edit: nvm, I see ni, nc, nblah are stuff/kg.

[PeterCaldwell]

Regardless of whether Homme keeps q around or not, don't we need to convert dry->wet anyways, before calling Homme? I mean, the (backed out) forcing FQ is dry on the AD side (since Q is dry in the Field manager), but it has to be wet in Homme, so the dyn interface would need a dry->wet conversion (and back, after Homme returns) anyways, no?

Yes, I think "physics is dry, dynamics is wet => conversion is needed in the to and from dynamics directions" is the right way to think about it. I think you could multiply dry qv by dry pseudo_density to get qdp without taking both the dry qv to wet qv, then wet qv * wet pseudo_density = qdp if you wanted... but you'd end up making the same computations in the end anyways.

[PeterCaldwell]

I was wondering about stuff like radiation tracers, or even ni, nc, nblah (for instance).

Yeah, I think the correct guidance is that all tracers with units per kg air will be subject to this wet/dry conversion and all the other ones won't be.

[bartgol]

Ugh, this is a bit complicated. But I guess I can query the units of the individual tracers Field's, and figure out for which iq I need to convert Q(:,iq,:). I guess we found a real use for having fields units. :)

[PeterCaldwell]

Note that the fundamental equation relating "dry" mixing ratio and "wet" specific humidity is qv_wet = qv_dry/(1+qv_dry).

Fun fact: This can be derived by noting that qv_wet is the mass of vapor divided by the mass of dry air + the mass of vapor. Dividing the numerator and denominator of this eq by mass of dry air yields qv/(qv+1).

[mt5555]

this is a side issue, but I would claim that the concept of mixing ratio extends to any quantity that is transported in a conservative way, including things like number density. So any quantity that is transported by the atmospheric flow:

d/dt (Qdp) + div( U Qdp ) = 0

will also obey

D/Dt (Q) = 0

Qdp is the quantity which is conserved, while Q (the wet mixing ratio) is the quantity that is moved around in a Lagrangian way (and is the quantity for which we need to enforce monotonicity).

[PeterCaldwell]

I would claim that the concept of mixing ratio extends to any quantity that is transported in a conservative way, including things like number density

I think we're all in agreement that all variables in EAM's "Q" array need the dry/wet conversion. This includes qv, qc, qi, qr, qm, ni, nc, nr, nm... It does not include TKE or other variables that don't have units of per kg air though because it wouldn't make sense to multiply them by dp... right?

[PeterCaldwell]

Qdp is the quantity which is conserved, while Q (the wet mixing ratio) is the quantity that is moved around in a Lagrangian way (and is the quantity for which we need to enforce monotonicity).

I'm not sure if I understand what you're saying here. Are you arguing for why Qdp and Q are both needed in HOMME?

[bartgol]

this is a side issue, but I would claim that the concept of mixing ratio extends to any quantity that is transported in a conservative way, including things like number density. So any quantity that is transported by the atmospheric flow:

d/dt (Qdp) + div( U Qdp ) = 0

will also obey

D/Dt (Q) = 0

Qdp is the quantity which is conserved, while Q (the wet mixing ratio) is the quantity that is moved around in a Lagrangian way (and is the quantity for which we need to enforce monotonicity).

The two are equivalent cause dp (aka, the mass) is conserved, right? So this mean that transporting Q or Qdp is (at least analytically) the same? And therefore you can transport Qdp=Q*dp, regardless of Q's units?

[mt5555]

I was saying that any quantity that is transported by the flow (independent of what its units are), will have a natural "mass", Qdp, and a mixing ratio Q. So I dont think this depends on the units, but on how the constituent is transported in the atmosphere.

CLUBB for example I think has a mode where TKE is transported by the dycore, so it would have to be added to the Q array.

[mt5555]

So I think we all agree: everything in the Q array should be converted to a wet mixing ratio and than passed to the dynamics interface layer. (or this conversion could be done inside the dynamcis interface layer). The dynamics interface layer will need both the new state (after the physics) and the old state in order to back out tendencies. It will then do a few more complex steps outlined above, before populating the GLL-based HOMME data structures.

[bartgol]

I was saying that any quantity that is transported by the flow (independent of what its units are), will have a natural "mass", Qdp, and a mixing ratio Q. So I dont think this depends on the units, but on how the constituent is transported in the atmosphere.

CLUBB for example I think has a mode where TKE is transported by the dycore, so it would have to be added to the Q array.

Ok, this makes sense. But I still have one doubt: does the dry->wet conversion makes sense for stuff that has no blah/kg units? Shouldn't these tracers simply be multiplied by dp (wet) and stuffed in homme's qdp, rather than first getting mulitplied by dp_dry/dp_wet?

[mt5555]

as @PeterCaldwell said above: "I think we're all in agreement that all variables in EAM's "Q" array need the dry/wet conversion. "

[PeterCaldwell]

Way to use my words against me, @mt5555 ! My previous statement was incorrect - I was assuming that EAM's Q array didn't include stuff like TKE because I didn't understand that Q*dp made sense for TKE.

I totally agree with Luca that it doesn't make sense to convert wet<->dry for variables not in units of per kg air. Take TKE for example. It has units m2/s2 and doesn't have anything to do with mass of air - either wet or dry. What would it even mean to write TKE_wet = TKE_dry/(1+TKE_dry)?

Maybe this was why EAM did all its physics in "wet" values even though all the physics schemes are more naturally defined for dry values.

[mt5555]

how about, given our current plans for transported species, "I think we're all in agreement that all variables in EAM's "Q" array need the dry/wet conversion. "

[PeterCaldwell]

how about, given our current plans for transported species, "I think we're all in agreement that all variables in EAM's "Q" array need the dry/wet conversion. "

Isn't TKE a transported species for SHOC in SCREAM?

[bartgol]

On the other hand, if Q has arbitrary units (not stuff/kg), then Qdp is a "fictitius" mass. If we're ok with a fictitious mass, then what do we care if we transport Qdp_wet or Q(dp_dry/dp_wet)dp_wet = Qdp_dry "mass"?

[PeterCaldwell]

On the other hand, if Q has arbitrary units (not stuff/kg), then Qdp is a "fictitius" mass. If we're ok with a fictitious mass, then what do we care if we transport Qdp_wet or Q(dp_dry/dp_wet)dp_wet = Qdp_dry "mass"?

But dp_dry isn't conserved by the dycore, right?

[PeterCaldwell]

Maybe just for my own clarification but putting here in case it helps anyone else:

A nuance to my formula to convert between wet and dry mass for mixing ratios other than qv: qc_wet = qc_dry/(1+qv_dry). I was thinking you'd just replace qv with qc everywhere if you wanted to the dry value for qc (or any other variable), but that's not true.

With this detail hashed out, we see that qcdp = qc_wet pseudo_density_wet = qc_dry pseudo_density_dry since pseudo_density_wet = (1+qv_dry) pseudo_density_dry [because pseudo_density_wet = pseudo_density_dry + pseudo_density_vapor and pseudo_density_vapor=qv_drypseudo_density_dry].

[bartgol]

@PeterCaldwell , I think the way I see it (and the way I will implement it) is: compute conversion factor dp_wet/dp_dry first. If water_loading=q_v only, that is 1+qv. Then process all tracers in the same way, by doing q_wet=q_dry*conv_factor. I think this removes the apparent difference between qc and qv.

[PeterCaldwell]

@bartgol - I'm not opposed to your "scale by dp_wet/dp_dry" approach, but it takes a bit of grappling to understand what dp_wet and dp_dry even mean and the two cancel out to (1+qv) anyways. On the other hand, qX_wet = specific humidity of X = 1/(1+ qv_dry) * qX_dry can be found on wikipedia or in any atm thermodynamics text book. And it is simpler than calculating dp_wet and dp_dry separately.

Finally, we should make sure to put the conversion to and from wet and dry into a function since, as Mark said, we will eventually want to add the mass of liquid and ice.

And of course we're still trying to decide whether the field manager should store q*dp anyways.

[PeterCaldwell]

Interesting comments from @beharrop about forcing all vars to be "wet":

Yea, the trace gases and aerosols were “dry” tracers, but CLUBB and the GW parameterizations assumed everything was a wet tracer and this broke mass conservation. I suspect no matter what you choose someone will be unhappy and something will break. Even “wet” specific humidities could cause confusion with whether the total water mass or just the mass of vapor is used, and might be nasty to debug given how close those number would look to one another. One of the CO2 gotchas that NCAR helped us with is that simply switching CO2 to a wet mixing ratio before CLUBB meant we went from having a relatively constant, well-mixed profile of a dry tracer, to one that had a gradient owing to water vapor. So CLUBB started mixing an already well-mixed tracer and creating artificial maxima and minima.

I don’t have any particular opinion on which of those two would be better. I think regardless of which you end up choosing, it would be good to have some tracer tests for any new parameterization that gets hooked up to the model like conserving mass, not introducing new minima or maxima, preserving shape (where appropriate), etc. Out of curiosity, what is the vertical coordinate for the SCREAM physics? Is it still the hydrostatic hybrid sigma-pressure coordinate or is it something more akin to what one might find in a CRM?

[PeterCaldwell]

Before hearing Bryce's comment, I had settled on making everything in the field manager a "wet" specific humidity and having a handy function for converting wet stuff to dry. This has the advantage of being a more natural coordinate than "total mass" and means we don't have to make any changes to quantities like TKE which are naturally/implicitly defined as "per kg total air mass". As Bryce notes, though, this may introduce artificial gradients in things that would have been well mixed in dry mass coordinates. I"m not sure how big of a deal this is.

[mt5555]

The case mentioned above - CLUBB being sensitive to water vapor gradients, seems to be saying that if we store wet mixing ratios in the FM, then CLUBB has to use the wet mixing ratio in its calculations. But CLUBB could still continue to work with dry mixing rations and just convert before and after.

It might be interested to store mass in the FM - this will force all parameterizations to explicitly convert to the mixing ratio they think they should be working with.

[beharrop]

I will clarify, just to make sure we are all on the same page. CO2 is generally measured to be a well-mixed gas relative to "dry" air. In EAM, we keep CO2 as a "dry" mixing type, but CLUBB assumes all tracers are "wet" mixing types. This was breaking CO2 mass conservation in EAMv1, and the initial solution NCAR came up with (https://github.com/E3SM-Project/E3SM/pull/2883), which we adopted, was to convert CO2 to a wet mixing type right before CLUBB was called, and then convert back afterwards. Keith Lindsay found that this solution was not a complete fix because converting a well mixed dry tracer (zero vertical gradient) to a wet tracer gave it a vertical gradient because water vapor was included. CLUBB would then mix CO2, creating new maxima and minima on the CO2 profile. Thus, a new solution was implemented to conserve both the mass and the shape of the CO2 profile (https://github.com/E3SM-Project/E3SM/pull/3003).

This is one very specific example, and while it seems to suggest dry mixing ratios would be the ideal solution, there are other examples where that won't be the case (just looking back at the discussion on this thread). I think @mt5555's idea of storing the mass and forcing parameterizations to pay attention to what is coming in and going out it kind of nice. I will also reiterate what I sent @PeterCaldwell via email that it would be good to have some tests that could be run after each parameterization to check for things like conservation and shape.

[PeterCaldwell]

Thanks @beharrop ! Making sure there are good tests for whatever approach we take is a great idea. I also like @mt5555 's point that just because we store stuff "wet" doesn't mean CLUBB couldn't convert to dry before advecting. I'm currently seeing all options (wet, dry, total mass) as being essentially interchangeable as long as users are careful to keep track of what is needed.

Given that, I think we should stick with making everything "wet" for v1. In particular, I want to just choose something rather than getting mired down. Since we won't have prognostic CO2 in v1, worrying about supporting it isn't worthwhile. We can add support for dry or wet variables in v2 or whenever it becomes relevant.

[PeterCaldwell]

Replaced with the more action-oriented #1066