qice007 (last ice-enthalpy layer) corrupted on ~99% of ice-covered cells every Soca2Cice cycle, even with a zero aice/hice increment

[guillaumevernieres]

Summary

In PostProcessIce::applyThermoStage (src/soca/PostProcess/PostProcessIce.cc),
the "surface ice layer" index used to cap ice enthalpy against the surface
temperature is computed as the last layer instead of the first:

// PostProcessIce.cc:807
const std::size_t lSurf = iceLev - 1;

CICE's layer convention (confirmed via siceLayerCice4's salinity profile
in IcePhysics.cc: layer 1 → salinity ≈ 0/fresh/surface-adjacent, layer
nlyr → salinity ≈ saltmax/salty/ocean-adjacent) is layer 1 = top
(surface), layer nlyr = bottom (ocean interface). So lSurf = iceLev - 1
selects the bottom layer, not the surface one.

At L852–859, whenever a category has ice (aice > 0) and Tsfcn < 0 (true
almost everywhere, always):

// Cap the surface ice layer enthalpy by iceEnthalpyBL99(Tsfcn, sice).
const double sice = siceCatLev[k][lSurf](jnode, 0);
const double sBL  = (sice > 0.0) ? sice : sLayer[lSurf];
const double qCap = icephysics::iceEnthalpyBL99(T, sBL);
double & qIce = qiceCatLev[k][lSurf](jnode, 0);
qIce = std::min(qIce, qCap);

this clamps the bottom, ocean-adjacent layer's enthalpy (qice00N,
N = iceLev) to whatever the cold atmosphere-adjacent surface
temperature (Tsfcn, e.g. -20 to -30 °C in winter) implies — instead of the
much warmer, physically correct value near the ocean freezing point
(~-1.8 °C). This is a large, spurious, one-directional cooling of the
ice-ocean interface layer.

Impact

• Fires unconditionally on essentially every ice-covered cell, every cycle
  the postprocessor runs, independent of whether there is any real DA
  increment — confirmed in a cycle with a true zero aicen/vicen
  increment (0.00000 fractional change, full-array diff), where the last
  layer's enthalpy was still altered on 99.98% of ice-covered cells.
• Under ktherm=2 (mushy-layer thermodynamics), which is sensitive to the
  ice's internal enthalpy profile, this drives anomalous basal/conductive
  (congelation) ice growth: observed as unrealistic, monotonic sea-ice
  volume growth (Arctic ~3.4x, Antarctic ~8.8x over ~9 days of cycling)
  relative to a baseline run without this code path, compounding every
  cycle, present in both hemispheres, and independent of ice-observation
  coverage (reproduced with AMSR2/ice obs excluded).

Evidence

Diffing the ice restart's per-layer enthalpy, background vs. postprocessed
analysis, at a divergent cycle and independently at a cycle with a
confirmed-zero ice increment:

layer	cells changed (real DA cycle)	cells changed (zero aice/hice increment cycle)
layer 1–6 (qice001–qice006)	0.135% (legitimate ITD rebins)	0%
layer 7 (qice007)	99.4%	99.98%

Where changed, the magnitude is comparable to the entire quantity itself
(background O(-2×10⁸) to O(-3×10⁸) J/m³, deltas up to ~3.3×10⁸ J/m³).
sice, qsno, and Tsfcn are unaffected — sice[lSurf] is only read here
to inform the salinity term, never written.

Suggested fix

const std::size_t lSurf = 0;

Suggested verification

• TestIcePhysics.cc / soca2cice ctests at low resolution.
• Re-run a cycle with a zero ice increment and confirm qice00N (all
  layers) stays bit-identical to background.
• Re-run a cycling experiment and confirm congel_h/frzmlt_h stay flat
  and matching a baseline without this code path.

[shlyaeva]

Note: this is from the experiments that use PostProcessIce, not Fortran soca2cice, so the code references above are wrong.

[shlyaeva]

Could it be the surface level enthalpy that's updated in

soca/src/soca/PostProcess/PostProcessIce.cc

Lines 855 to 862 in 3b4b23a

	// Cap the surface ice layer enthalpy by iceEnthalpyBL99(Tsfcn, sice).
	if (T < 0.0) {
	const double sice = siceCatLev[k][lSurf](jnode, 0);
	const double sBL = (sice > 0.0) ? sice : sLayer[lSurf];
	const double qCap = icephysics::iceEnthalpyBL99(T, sBL);
	double & qIce = qiceCatLev[k][lSurf](jnode, 0);
	qIce = std::min(qIce, qCap);
	}

?

[guillaumevernieres]

Note: this is from the experiments that use PostProcessIce, not Fortran soca2cice, so the code references above are wrong.

Oops, I'll update.