search

E3SM-Project / E3SM

Energy Exascale Earth System Model source code. NOTE: use "maint" branches for your work. Head of master is not validated.
https://docs.e3sm.org/E3SM
Other
345 stars 352 forks source link

Non-Standard Exner Definition in dp_coupling.F90? #2680

Open PeterCaldwell opened 5 years ago

PeterCaldwell commented 5 years ago

@AaronDonahue noticed that the Exner function defined in the SE dycore's dp_coupling routine and passed into physics uses a non-standard (incorrect?) definition for the Exner function. The Exner function is the conversion factor (p/p_0)^{Rd/cp} needed to convert potential temperature into temperature. In this equation, p_0 is a reference pressure typically taken to be 1000mb (see e.g. http://glossary.ametsoc.org/wiki/Potential_temperature). In E3SM, p_0 is instead taken to be pint(i,pver+1), which is (at least approximately) the time-and-space-varying surface pressure. Here are the lines of code:

https://github.com/E3SM-Project/E3SM/blob/d0241119c45f6864d90fdee50dabccec3bdfa865/components/cam/src/dynamics/se/dp_coupling.F90#L511-L512

This weird definition caused conservation errors when Aaron implemented code which used the normal definition of theta in some places and this weird Exner in others... so at the very least, this current definition is dangerous.

Aaron used grep to identify locations using Exner and found it was used for the following:

To calculate the bottom potential temperature for export to the surface models. In CLUBB* In vertical diffusion Unicon UW-Shallow Convection

and points out:

CLUBB is an interesting case. CLUBB carries two versions of exner, there is the standard CAM value in the state structure and there is a local “exner_clubb” that is used for CLUBB specific calculations. I talked with Peter B about it and this was because the definition of TH in CLUBB theory was inconsistent with CAM and they were experiencing massive issues over topography. Defining exner_clubb fixed the problem.

I suspect the weird definition of Exner was intended for export to the surface models but is being improperly used in the other cases (except CLUBB, which noticed the problem and created exner-clubb as a workaround)?

Does anyone know the history of this "feature"?

oksanaguba commented 5 years ago

is name of the var misleading, too? Exner is (p/p0)^\kappa ~ p^{\kappa} but there Exner is (p_bottom/p)^\kappa, so, ~p^{-\kappa}?

PeterCaldwell commented 5 years ago

Yes, I was wondering about that as well while I wrote this issue, but thought cappa was probably defined as -Rd/cp. I just checked in control/physconst.F90 and you're right - cappa=Rd/cp so E3SM's definition is also the reciprocal of the normal definition.

AaronDonahue commented 5 years ago

To clarify the line for exporting values to the surface is components/cam/src/control/camsrfexch.F90:462: cam_out%thbot(i) = state%t(i,pver) * state%exner(i,pver)

PeterCaldwell commented 5 years ago

Ok, so the line exporting values to the surface uses the weird e3sm exner correctly in the sense that Texner = theta instead of thetaexner = T as typically defined. As a sanity check, I just looked through a bunch of websites and confirmed that they all define exner as T/theta rather than the way e3sm does.

mt5555 commented 5 years ago

The FV dycore is doing the same thing. So this looks like a quirk of CAM physics to define exner pressure in terms of surface pressure instead of p0. The variable name is very missleading - what about renaming it "rexner" , with "r" for reciprocal?

update: EUL dycore also using this nonstandard definition of 'exner'

mt5555 commented 5 years ago

do we still use the UW shallow scheme in v1?

PeterCaldwell commented 5 years ago

No, the UW shallow scheme is only used in v0. I like the idea of renaming it rexner, like the T Rex of Exners. Let's wait to do it until we understand whether using p0 = surface pressure is appropriate though.

AaronDonahue commented 5 years ago

Could it be that it was defined this way in the days of FV and EUL dycores to be consistent, but with the SE dycore we are no longer consistent?

mt5555 commented 5 years ago

it looks completely independent of the dycore - its a computation based on dycore state to be passed to the physics, and then (based on your comments above) the surface models. So the definition has to be consistent with what the surface models are expecting?

PeterCaldwell commented 5 years ago

Dave told me this weird definition is an ancient practice from nonhydrostatic modeling which was supposed to be more computationally efficient...

AaronDonahue commented 5 years ago

I did a little digging and it looks like the surface coupling value thbot that is used with this definition of exner is used to determine the stability of the solution in the ATM->OCN coupling. Since it is over the ocean there is probably not too much deviation in PS from P0 so this might be why the strange definition isn't too much of a problem.

That being said, I grabbed a 10 year ne30 climatology I had laying around and looked at PS over the ocean and found that in that case there can be about a 15 degree swing in potential temperature over the ocean depending on the definition you use (P0 vs. PS). Disclaimer the run was an F case so prescribed OCN.

Lastly, based on the definition of exner, the value of exner (in this context) at the surface will be a constant because pint at the surface is just ps and pmid at the bottom layer is just ps*hybm, so the ratio is always hybm(pver)**(-1).

AaronDonahue commented 5 years ago

UPDATE I did some more searching around and took a look at the other dycores and how they handle the "Exner" expression. The short summary is, they all handle it in basically the same way. The longer explanation:

The Eulerian (EUL), Semi-Lagrangian (SLD) and Spectral Element (SE) dycores all use the exact same expression for the definition of exner as Peter noted in the original issue statement. The Finite Volume (FV) dycore uses a slightly different formulation: phys_state(lchnk)%exner(i,k) = pic(i) / pkzxy(ic,jc,k) I am not a FV-dycore expert but I did search around for what pic and pkxzy represent and to the best of my knowledge, pic is the surface pressure raised to the cappa power and pkxzy is the pressure at mid-layer raised to the cappa power. So the same definition as the other dycores but in the FV syntax.

I want to note that the FV dycore also has exner defined in two other places, for gravity waves (fv/gravity_waves_sources.F90) and circulation diagnostics (fv/ctem.F90). In both of these cases exner is calculated locally to be used to determine the potential temperature. And in both of these cases exner is calculated in the way we would expect with a constant reference temperature instead of the surface temperature.

pexf = (ps0 / pm(i,k,j))**cappa

So it would seem that there was a initial definition of exner (which as we already noted is really the reciprocal of exner) in one of the dycores and this has been maintained from then on in the dp_coupling layer. And whenever an individual process actually needs to calculate the potential temperature it derives exner locally. The only exception is that exner is used to determine the bottom potential temperature for the atm->ocn coupling.

Lastly, on the atm-ocn coupling, the thbot value is used to calculate the potential temperature difference and the stability in the ocn/atm interface: delt = thbot(n) - ts(n) ! pot temp diff (K) !--- compute stability & evaluate all stability functions --- hol = loc_karman*loc_g*zbot(n)* & (tstar/thbot(n)+qstar/(1.0_R8/loc_zvir+qbot(n)))/ustar**2 in routines: SUBROUTINE shr_flux_atmOcn_diurnal and shr_flux_atmOcn

For ATM->ICE (which I'm not sure is called) it is also used to calculate the virtual potential temperature (see http://glossary.ametsoc.org/wiki/Virtual_potential_temperature) thvbot = thbot(n)*(1.0_R8 + loc_zvir * qbot(n)) ! virtual pot temp (K) in routine shr_flux_atmIce

mt5555 commented 5 years ago

Assigning this to @AaronDonahue since I cant tell if this is a bug or just a poor naming convention. At the very least we should rename it rexner_appox (reciprocal approximation) so that it is not accidently used as exner pressure in the future.

AaronDonahue commented 5 years ago

At the very least it seems like a poor naming convention, but it would be nice to hear from anyone in the Ocean group about whether or not it is also a bug since what we are calling the potential temperature thbot is really basically the temperature.

maltrud commented 5 years ago

at the ocean surface, potential temperature and in-situ temperature are the same. our prognostic field is potential temp. does this help?

AaronDonahue commented 5 years ago

@maltrud , yes, thank you. I think that answers the question about whether or not this definition is problematic.

I suppose the next question for the atmosphere folks is whether or not having something defined as "Exner" in the model could potentially cause bugs in the future if someone takes that to mean the true Exner formula when in fact it is not. For example, I only noticed this when I attempted to use exner in the phys_state structure to get potential temperature to use with P3 microphysics and got energy conservation errors. If we think this could be a problem we may want to fix the name to a) show that it is actually the inverse of the exner formula and b) maybe a name other than exner since it doesn't use a consistent reference pressure globally.