From f5a0914cf95ae6d2b14149f1e8e2f65478004e62 Mon Sep 17 00:00:00 2001 From: Patrick Scholz Date: Fri, 26 Jun 2026 15:36:07 +0200 Subject: [PATCH] Fix dimensionally inconsistent mass regularisation in EVP stress-to-RHS (modes 1 and 2) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The original expression `mass/(1+mass²)` mixes dimensionless 1 with [kg/m²]², making the formula physically meaningless. Replace with the correct `1/max(mass, 9.0)`, where 9.0 kg/m² is the same floor already used for inv_thickness in the same loop (≈1 cm of pure ice, ρ_ice ≈ 900 kg/m²). The fix is applied in stress2rhs_m (shared by mEVP mode 1 and aEVP mode 2) and the inlined equivalent in EVPdynamics_m (mode 1). --- src/ice_maEVP.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/ice_maEVP.F90 b/src/ice_maEVP.F90 index aec61801c..3de24a770 100644 --- a/src/ice_maEVP.F90 +++ b/src/ice_maEVP.F90 @@ -413,7 +413,8 @@ subroutine stress2rhs_m(ice, partit, mesh) if ( ulevels_nod2d(row)>1 ) cycle mass=(m_ice(row)*rhoice+m_snow(row)*rhosno) - mass=mass/(1.0_WP+mass*mass) + mass=1.0_WP/max(mass, 9.0_WP) ! 9.0 kg/m² per GRID area — numerical floor to prevent near-zero inertia + !mass=mass/(1.0_WP+mass*mass) ! original: dimensionally inconsistent (1 + [kg/m²]²) u_rhs_ice(row)=(u_rhs_ice(row)*mass + rhs_a(row))/area(1,row) v_rhs_ice(row)=(v_rhs_ice(row)*mass + rhs_m(row))/area(1,row) end do @@ -649,10 +650,11 @@ subroutine EVPdynamics_m(ice, partit, mesh) if (a_ice(i) >= 0.01_WP) then inv_thickness(i) = (rhoice*m_ice(i)+rhosno*m_snow(i))/a_ice(i) - inv_thickness(i) = 1.0_WP/max(inv_thickness(i), 9.0_WP) ! Limit the mass + inv_thickness(i) = 1.0_WP/max(inv_thickness(i), 9.0_WP) ! 9.0 kg/m² per ICE area (≈1 cm ice per unit ice area) mass(i) = (m_ice(i)*rhoice+m_snow(i)*rhosno) - mass(i) = mass(i)/((1.0_WP+mass(i)*mass(i))*area(1,i)) + mass(i) = 1.0_WP/(max(mass(i), 9.0_WP)*area(1,i)) ! 9.0 kg/m² per GRID area (different quantity; same floor for numerical safety) + !mass(i) = mass(i)/((1.0_WP+mass(i)*mass(i))*area(1,i)) ! original: dimensionally inconsistent ! scale rhs_a, rhs_m, too. rhs_a(i) = rhs_a(i)/area(1,i)