Fix dimensionally inconsistent mass regularisation in EVP dynamics (modes 1 & 2)

[patrickscholz]

Problem

In stress2rhs_m (called by both mEVP mode 1 and aEVP mode 2) and the inlined equivalent in EVPdynamics_m (mode 1 only), the column ice mass per unit area is regularised using:

mass = mass / (1.0_WP + mass*mass)

This expression is dimensionally inconsistent: the denominator adds a dimensionless 1 to mass² which has units of (kg/m²)². The formula has no physical meaning.

Effect

For typical Arctic sea ice (e.g., 1 m ice + 0.1 m snow):

mass ≈ 917 × 1.0 + 330 × 0.1 = 950 kg/m²
original:  mass / (1 + 950²) ≈ 1.1 × 10⁻⁶   (extremely over-damped)
correct:   1 / 950           ≈ 1.05 × 10⁻³   (true inverse mass)

The original formula suppresses ice acceleration by a factor of ~1000 more than physically correct for thick ice. For thin ice (mass ≈ 1 kg/m²) the expressions are closer, but the formula is still wrong.

Modes affected

• Mode 1 (mEVP, EVPdynamics_m): affected in two places — stress2rhs_m (shared RHS assembly) and the inlined fused loop inside EVPdynamics_m
• Mode 2 (aEVP, EVPdynamics_a): affected via stress2rhs_m (called from EVPdynamics_a)
• Mode 0 (legacy EVP, EVPdynamics): uses a different RHS path — not affected

Fix

Replace the dimensionally inconsistent expression with the physically correct inverse mass, floored at 9.0 kg/m² (≈ 1 cm of pure ice, ρ_ice ≈ 900 kg/m³) to prevent near-zero mass from causing numerical problems:

! stress2rhs_m  (modes 1 & 2)
mass = 1.0_WP / max(mass, 9.0_WP)   ! 9.0 kg/m² ≈ 1 cm ice; same floor as inv_thickness above

! EVPdynamics_m fused loop  (mode 1 only)
mass(i) = 1.0_WP / (max(mass(i), 9.0_WP) * area(1,i))   ! 9.0 kg/m² floor; matches inv_thickness floor above

The 9.0 kg/m² floor is already used in the same loop for inv_thickness (! Limit the mass comment, line 653), so this is consistent with the existing convention in the code.

The original expressions are kept as comments.

Files changed

• src/ice_maEVP.F90: lines 415–416 (stress2rhs_m) and 655–656 (EVPdynamics_m)

[patrickscholz]

i will make a test bias run for this one as well!

[patrickscholz]

done for mEVP (but aEVP looks the same)

• difference in Summer/Winter NH/SH sea ice concentration:

• maximum NH/SH MLD2

• temperature/salinity bias

Fazit: very minor differences!