Notes on enthalpy flux implementation in CESM‐CAM‐CMEPS‐MOM6

git checkout enthalpy_revision_19sep2024

The overall goal is to move from CAM's current energy formula to a more rigorous energy formula using variable latent heats. The rigourous total energy formula is shown on this slide:

And the one used in CAM (constant latent heats and other assumptions) is:

Some observations on magnitude of latent heat terms

Plots below are from a 10 year FLT run (1995-2005)

The constant latent heat terms are much larger than the "ignored" T-dependent latent heat terms (lower left and right plots).

Rigorous implementation

After much back and forth this is the simplest implementation I could come up with (both in terms of coding and science). Any form of rigorous implementation will need the following information from each parameterization that performs phase changes and in any other way changes water in the column (even mixing water between levels):

• water changes due to phase changes that do not leave the grid-box
• water changes that leave the grid box due to falling precipitation and evaporation
• temperature at which phase changes occur + temperature of falling precipitation
• schemes that mix water species need to use variable latent heat formulas for energy conservation
• all the things I have not thought off …

Also, schemes that use energy fixer need to use new energy formula (e.g. clubb_intr.F90, gravity waves?, …)

As simple as possible implementation

The adjustments are done after CAM physics and using the full (dycore) total energy formula. Infrastructure has been put in to compute water fluxes and associated enthalpy fluxes after tphys_bc (in src/control/camsrfexch.F90) and after tphys_ac using physics buffers.

Algorithm:

• After CAM physics scale temperature tendencies for consistency with the dycore:

• Compute total energy after physics using dycore energy formula:

at this point the total energy budget is closed (assuming that pressure is constant) as shown in the Figure below:

• Do dry-mass adjustment

Associated total energy tendency in each column (10 year average)

• Update cp

Associated total energy tendency in each column (10 year average)

The mis-match in the energy budget at this point (dE/dt dme adjust + dE/dt cp update):

We will return to this imbalance in a moment.

Enthalpy fluxes (using SST for evaporation enthalpy flux and lowest atmosphere model level temperature for precipitation enthalpy fluxes) are:

(a) Enthalpy flux from ZM (b) Enthalpy flux from CLUBB-PUMAS (c) Enthakpy flux evaporation (d) Net enthalpy flux

Observation: if we subtract the total enthalpy flux from the "end of physics adjustments" then the energy budget is close to balanced in areas where evaporation dominates:

E.g. off the West coat of South America and South Africa.

However, over areas of heavy precipitation the energy budget is far from balanced.

Why?

• For evaporation the water vapor is added to lowest model levels and therefore the dme_adjust + cpdycore change closely matches the enthalpy flux.
• For precipitation it is much more complicated:
  • The temperature at which precipitation was formed is different than the lowest model level (we will show later that this effect is on the order of 3-4 W/m^2)
  • Precipitation in the parameterization may convert water vapor to falling precipitation so dme_adjust only sees a change in water vapor whereas the enthalpy flux sees a flux of liquid or snow with very different heat capacities. Hence there is a missing heating in the column associated with performing phase changes leading to falling precipitation under variable latent heat assumptions.

Fixer algorithm

In the following we ignore kinetic energy and the "PHIS" energy terms (they are very small compared to the enthalpy terms).

We attempt to fix this by adding heating in the column to balance the residual between the precipitation enthalpy flux and the "end of physics adjustments" shown in Figure (c) above. For that, we need to estimate where in the column precipitation is formed and how much is formed.

We split our fix into tphys_bc and tphys_ac:

• tphys_bc (ZM): for ZM we use rain production

• tphys_ac (CLUBB+PUMAS): we use cloud ice tendency from PUMAS and liquid water tendency from CLUBB

I also conducted experiments just using total heating from ZM and PUMAS and the F-case results were very similar:

https://webext.cgd.ucar.edu/FLTHIST/enthalpy/cam_enth_simple_ocnfrc_fix2/

Whatever the distribution function is we do the following: Compute the mass of precipitation production from tphys_bc and tphys_ac:

And thereafter we normalize the precipitation mass fluxes in each column:

We then compute the fraction of enthalpy precipitation flux from tphys_bc and tphys_ac:

And then we distribute the residual enthalpy flux as heating in the column according to our weightening functions:

The residual enthalpy flux is

shown in Figure (c) above.

In plots below the enthalpy flux in tphys_bc and tphys_ac are shown as a reference in (a) and (b) respectively and below the enthalpy flux fix:

And here is the residual that we are compensating for: (total enthalpy flux minus dEdt dme and dEdt cpdycore update) as well as the imbalance after applying the enthalpy heating fix (right plot):

This leaves a rather small imbalance to fix by the global energy fixer.

BUT ... only the ocean will see the enthalpyt fluxes so we need to fix enthalpy fluxes over land etc. with the global fixer in which case the imbalance is ~7 W/m^2 (plot shows the energy being fixed by the global energy fixer):

We have run a series of 10 year FLTHIST simulations 1995-2005 and ADF diagnostics can be found here:

https://webext.cgd.ucar.edu/FLTHIST/enthalpy/

The baseline run is cam6_4_015/

and the enthaly run (where we assume all surface components can receive enthalpy flux prescribed by the atmosphere is cam_enth_simple/

In all other runs we assume only the ocean can receive enthalpy fluxes using the suffix ocnfrc:

cam_enth_simple_ocnfrc/

We then performed a series of experiments varying various tuning parameters but none of them seem to remove the large precipitation biases!

In short, adding the enthalpy fix heating has a large effect on the circulation. The precipitation changes are large:

Baseline precipitation bias against OBS:

Enthalpy flux run precip biases (compared to OBS):

Implementation log - outdated (but maybe useful)

• Jim merged Thomas Toniazzo's CMEPS code into his branch: https://github.com/jedwards4b/CMEPS/pull/2/files
• Note: This file SourceMods/src.drv/esmFldsExchange_cesm_mod.F90 is still missing - Jim reviewing with Mariana
• Peter merged Thomas Toniazzo's code for adding hrain, hsnow, hevap to the coupler fields in CAM; note that fluxes are computed in control/camsrfexch.F90
• added diagnostics: HRAIN, HSNOW, HEVAP to cam_diagnositcs
• added debugging diagnostics:
  • HTOT_ATM_OCN: enthalpy flux + latent heat terms using ice reference state and only over ocean
  • HTOT_OCN_OCN: enthalpy flux + latent heat terms using liquid reference state and only over ocean

Code is in this branch (based on cam6_3_132 - the tag we are currently using for coupled runs):

git clone https://github.com/PeterHjortLauritzen/CAM cam-h
cd cam-h
git checkout cam-h
./manage_externals/checkout_externals

and remember to check out Jim's CMEPS branch:

cd components/cmeps
git remote add jimdev https://github.com/jedwards4b/CMEPS.git
git fetch jimdev
git checkout enthalpy_corrections

To turn on CMEPS enthalpy code:

compute_enthalpy=.true. in user_nl_drv

Fixer code here:

components/cmeps/mediator/med_enthalpy_mod.F90

The atmosphere uses an ice reference state:

where the reference temperature is 0K (and h_00^{(ice)} is zero).

The ocean uses a liquid reference state:

where the reference temperature is 273.15K (and h_00^{(liq)} is zero [note typo in equation above h_00^{(wv)}].

Equation (C13) can also be written as

Note that if we assume constant latent heats (L(T)=L_00) then this equation reduces to what MOM is currently using (cp(liq) for all the entalphy flux terms).

Note that if we compare the right-hand side terms (excluding the radiation and sensible heat flux terms) for the ice (below left - color scale from -100 to 800 W/m2) and liquid reference state (below right - color scale from -100 to 800 W/m2),

they should differ by the reference temperature term on the right-hand side of (C13) (color bar -20ish to 140ish W/m2):

(this is not entirely the case ... why?)

Plots above are instantaneous output at time-step 3 of an F-case (compute_enthalpy=.false.).

I still have some conceptual questions:

1* where should we account for the difference between cpsw and cpliq? 2* we still need to varify with MOM developers that it is consistent with MOM to use cpliq, cpice, cpwv for rain, snow and evaporation! The ocean obviously does not have phase changes so everything on the left-hand side of C13 is only sea water (liquid)! So are we assuming that there is a thin layer where everything is converted to sea water? (perhaps see https://journals.ametsoc.org/view/journals/clim/30/22/jcli-d-17-0137.1.xml) 3* note that since CAM computes hevap then the ocean receives hevap from the previous ocean time-step (instead of the current one) 4* The time-loop in CAM needs consideration: we call dme_adjust at the end of physics (after tphysac) and not at the end of tphysbc. hrain and hliq contains precipitation from tphysac from previous CAM time-step and tphysbc in the current time-step. So dme_adjust and coupler enthalpy flux are not in sync!

My (Thomas Toniazzo's; TT's) thoughts on these points:

1. we need to make sure the "translation" between atmosphere and ocean enthalpy reference states is correct: this depends on the ocean TD formulation. If the latter conforms to TEOS-10, where all TD coefficients are referred to a standard T,S,p for ocean water (which might be the reason why cpsw is used in MOM), and its potential temperature theta is computed accordingly, then I think this may also imply differences in the computation of latent heats as well as what we call material enthalpy. It seems to me we need more info on the ocean TD to work this out.
2. mostly cf 1, but on the last bit, the physical assumption in the current (and classic) air-sea flux computation is that there exists a thin layer of saturated air at the temperature of the ocean "surface" (whatever is choosen as TS; it may be a skin temperature as in shr_flux_atmocn_diurnal) 3.& 4. a time-step delay is inevitable, but note that hevap, hrain, hsnow are all at the same time-level at the end of tphysac, where (currently) dme_adjust is called. So are the (non-material) sensible heat fluxes, and the latent heat fluxes. So this way there's consistency between all air-sea fluxes in the atmospheric component. The ocean however will receive hevap, hrain, and hsnow from the previous time-step, plus non-penetrative and penetrative heat fluxes from atmocn at the current time-step. As long as the ocean time-step is a multiple of the atmo time-step, by virtue of averaginh (or accumulating) this causes a minor error compared to doing the converse.

Possible next steps

• Gustavo to check G-case
• Coupled: Current implementation adds hrain_a, hsnow_a, hevap_a, hrofl_a and hrofi_a to hcorr. Sanity check this setup!
• Coupled: If we want to remove hrain_a, hsnow_a, hevap_a from hcorr then we need to remove the entalphy flux over ocean in the atmosphere so that the atmosphere energy fixer does not fix that part! How do we do that? See comment on CAM time-loop above!

TT: I have devised a way to match hrain etc as computed or passed in camsrfexch to the ocean with the fluxes computed in the atmospheric column by physics_dme_bflux. Will test it soon.

updates (10/8/2023)

• In the testing of the new code we found that computing hevap_a, sending it to the coupler where it is mapped conervatively using the PCoM, and then subtracting the reference term (F^{(evap)*Tref*cpliq}) leads to noisy enthalpy fluxes because F^{(evap)} in the coupler is from the current time-step whereas F^{(evap)} in CAM is from the previous time-step.

Below is instantanous output (from G-case) of (F^{(evap)*SST*cpliq}) computed by atmosphere minus ``(F^{(evap)Trefcpliq})``` computed in coupler on ocean grid:

When computing (F^{(evap)*SST*cpliq}-F^{(evap)*Tref*cpliq}) in the atmosphere and sending it to the coupler the noise goes away

hevap still looks a little "blocky" because the atmosphere resolution is slightly less than MOM resolution. This could be an issue for configurations where the ocean resolution is much higher than the atmosphere resolution!

Next steps: hcorr in the coupler can no longer fix atmosphere enthalpy (since we only have enthalpy fluxes with ocean reference state in the coupler); in any case it does not seem to be an ideal solution to fix enthalpy fluxes with a sensible heat-flux from the ocean anyway!

=> need to go to next step where the energy fixer does not fix enthalpy flux over ocean (hence no need for hcorr). This can either be done by subtracting atmosphere enthalpy flux from state%te_cur(:,dyn_te_idx) before the call

    call pbuf_set_field(pbuf, teout_idx, state%te_cur(:,dyn_te_idx), (/1,itim_old/),(/pcols,1/))

or by using Thomas's dme_adjust.

Some newer implementation notes

test 1: change physics to using SE dycore energy formula (scaling in tphysac for SE/FV3 can go away, MPAS scaling changes and now FV will need scaling)

test 2: turn on enthalpy flux calculation:

• Do reference run so we have baseline numbers for dme_adjust and dynamical core dissipation errors (use se_ftype=1 so no physics-dynamics coupling errors)

• Accumulate enthalpy flux in current physics time-step and use that for adjustment at end of physics; T fixer for energy conservation encorporating dme_adjust; check how much the T fixer fixes column by column and global average!

Diagnostics needed:

1. enthalpy flux only
2. dme_adjust only
3. 1.+2.
4. temperature fixer (how does it compare to current global fixer?)

Changing physics_update to using variable cp

IMPORTANT: if we use generalized cp in physics_update then we can only pass energy check if we include enthalpy flux terms (even then the budget won't be entirely closed). The temperature dependence in L(T) makes the surface term temperature dependent! Similarly parameterizations like ZM need to be reformulated for variable cp. I started doing this work (for example one will need to update cp once q has been updated and then back out %s). This is not trivial. Equattions in ZM need to be re-derived with this assumption. Similarly the energy fix in CLUBB interface layer would need to be recoded. This is too complicated to tackle right now hence we keep treating CAM physics as a unit and scale for "effectively" latent heats at the end of physics.

test 3: move surface coupling to beginning or end of physics?

Some implementation notes

Comparing increment of total water change (instantaneous) with water actually leaving the column (computed from precipitation variables)

    fsnow(:) = 1000.0_r8*(precsc(:ncol)+precsl(:ncol))                           !snow flux
    frain(:) = 1000.0_r8*(precc (:ncol)-precsc(:ncol)+precl(:ncol)-precsl(:ncol))!rain flux