Am conservation added

[EdouardBarrier]

Adding angular momentum conservations capabilities as outlined in Barrier & Madhusudhan 2025a (+ Barrier & Madhusudhan in prep). Essentially, implementing the CAM6 corrections of Toniazzi+2020 for CAM4. Now has an option for level by level correction, instead of a global fixer.

Also included: some QOL additions to diff_exocam.py, a large set of initial condition files used for various sub-Neptune runs, an option to give an internal temperature, and some changes to the ccsm_utils I found necessary (these don't necessarily need to be approved in the PR).

[storyofthewolf]

Hi Edouard!

Apologies for the slow review. The PR lists as "ready to merge" but there are few targeted edits I was hoping you could address before we merge. Needless to say perhaps, this comment was aided by Claude Code and Codex agents.

Many of these I can do myself, except for item #1. And I want to make sure the PR comes through with provenance from you, so on the Github page you are credited.

1. wg_hiord is never assigned — cd_core.F90

line  329:  real(r8) :: wg_hiord                     ! declared
line 1295:  ... + wg_hiord*vf(1,j,k)*(dpn(1,j,k) - dpo(1,j,k))*0.5_r8
line 1301:  ... + wg_hiord*vf(i,j,k)*(dpn(i,j,k) - dpo(i,j,k))*0.5_r8

grep -n wg_hiord cd_core.F90 returns exactly these three lines, no
assignment anywhere. So lines 1295 and 1301 apply whatever value is in the memory slot. Clearly it didn't make your model crash, but it is a silent error.

**Are we missing some code in the PR that defines wg_hiord, or is this genuinely an undefined variable that has snuck into the calculation and pass silently?

2. Related, xakap (declared 306, assigned at 420) is never used for anything, and am_press_crrct (declared 328, never assigned or read) are orphaned.

I think these can be harmlessly removed, again, unless we are missing some code here that involves them later.

3. Missing deallocate — cd_core.F90

Seven arrays allocated under if (do_am_fixes) at 406–417 (help, kelp,
dpn, dpo, ddpu, dpns, ddus); grep -n deallocate cd_core.F90 returns
nothing. These allocated variables need a matching deallocate statement to avoid memory leaks.

4. tau4 halved outside the AM conditionals — cd_core.F90:641

-  tau4 = 0.01_r8  / (abs(dt))
+  tau4 = 0.005_r8 / (abs(dt))

Outside every do_am_fixes guard, so it halves 4th-order divergence damping for
all runs regardless of flag state. I'm guessing this is sub-Neptune stability tuning
that came along from your working tree. Please provide an explanation for this change. Consider, should we keep the new value for all atmospheres? Or is this something you think should be promoted to an parameter in exoplanet_mod.F90 for modifying between atmosphere types being run? I am fine either way, I just want a justifiation.

5. Slab ocean q-flux and internal heat — docn_comp_mod.F90:655,657

This is a great addition to the model, but it needs to be conditional and tied back to exoplanet_mod.F90 to facilitate different assumptions using the same code-base.

line 655:  !- avstrm%rAttr(kqbot ,n) &        ! prescribed q-flux commented out
line 657:  + shr_const_stebol*t_int**4 ) * &  ! internal heat added unconditionally

Both unconditional, so with t_int = 30.0_r8 every slab run gets σT_int⁴ ≈
46 W/m² added to every cell while prescribed ocean heat transport is removed.
For Earth-like aquaplanets that's a large spurious heat source and a silent
removal of ocean heat convergence.

Suggested:

! exoplanet_mod.F90
logical,  public, parameter :: do_internal_heat = .false.  ! sub-Neptune slab
real(r8), public, parameter :: t_int = 30.0_r8             ! internal temperature (K)

.false. restores qbot and omits the σT_int⁴ term. I'd want it permanently
default-off — different ocean configuration, not an improved one. (shr_const_stebol
itself is fine, use shr_const_mod is already in scope at line 8.)

6. Flag defaults for am_fix flags in exoplanet_mod.F90

Both am_fix flags default .true. (exoplanet_mod.F90:30,32). I'd like .false. on
merge, flipped later in a tagged release: work published against v1.0.0 would
otherwise stop reproducing on rebuild, and cam_aqua_fv would diverge from
cam_land_fv/cam_mixed_fv until I propagate the AM work.

I am glad that these flags are in place for the PR. Eventually I could see just deleting the flags and having your am fixes run as the default/only mode. But for now, they are good!

6. Two smaller things in the level-by-level path.

The AM GLOBAL FIXER log line is wrong when do_am_fix_lbl = .true.
dyn_comp.F90:2153 accumulates du_fix_g from the global increment am1/me0
regardless of which branch ran, and that's what prints at 2156. In
level-by-level mode the correction actually applied is dons(k)/dods(k) per
level, so the log shows a number that was computed but never used. The winds are
fine either way — u is correct in both branches, and du_fix_s is guarded
properly (2118 vs 2149). But this is the only AM diagnostic in atm.log, so it
would mislead anyone checking whether the fixer is working.

Also, the shr_reprosum_calc at 2133 runs unconditionally, and in
level-by-level mode its result only feeds that log line. Moving it inside the
global branch would fix the diagnostic and save a global sum every substep.

Level-by-level turns off pressure tapering. The else at 2084–2088 sets
taper(k) = 1 for all levels, so there's no pressure tapering in that mode at
all (apart from taper(1:kmtp) = 0 at 2091, which applies either way). Was that
intended? It's a second change riding along with the flag, and there's no way to
get level-by-level with tapering if someone wanted it.

Minor: sw_am_corr at 1184 turns off the sw_core half of the correction
if the transport orders aren't 4, while the cd_core fixer stays on — a startup
log message when that happens would save someone a confusing afternoon. And a
note on line 32 that do_am_fix_lbl does nothing when do_am_fixes = .false.
would be worth adding (the code already handles it correctly).

7. Machine configs — additive only, please

config_machines.xml is +38/−110 and config_compilers.xml +6/−61, deleting
hyak, summit, discover, computecanada and discover-mil. That removes
support for every machine documented in cesm1.2.1/instructions/, including
Discover, which all the namelist templates point at. You noted these don't
necessarily need approving — so please add ioa/ioa_2 and delete nothing.
mkbatch.ioa is good as-is.

8. NetCDF ICs — please drop from the PR

The 45 .nc files under cesm1.2.1/initial_files/cam_aqua_fv/ reverse the
git rm --cached cleanup from v1.0.0 and land in history permanently. Could you
keep them in your fork? I'm setting up proper data hosting on GoogleDrive, so not to overwhlem the size of the github repo. That is, I very much want all these 45 netcdf files, I just have to put them somewhere else (e.g. Google Drive).

9. Before and after regression tests.

Do you have any before and after tests to share with me, not for the PR here? In your papers you mention only high-temperature atmospheres are affected, T > 600 K. So presumably you have some climate simulation netcdfs somewhere for a variety of states, with and without the angular momentum fix.