From 0452259aa68b8afd5db906f8f14a8d234662421a Mon Sep 17 00:00:00 2001 From: Patrick Scholz Date: Thu, 25 Jun 2026 14:59:04 +0200 Subject: [PATCH] Fix inverted interpolation weights for Redi K_33 at vertical interfaces MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Background ---------- The Redi (1982) isopycnal diffusion scheme replaces simple vertical diffusion with a full diffusion tensor that acts along isopycnal surfaces. Its vertical (33) component is the effective vertical diffusivity: K_33 = Kv + Ki * slope² where Kv is the background diapycnal diffusivity (defined at interfaces, nl levels), Ki is the isopycnal diffusivity and slope² is the squared isopycnal slope (both defined at cell centres, nl-1 levels). The implicit vertical diffusion solver assembles flux coefficients at layer interfaces. For the interface between layers nz and nz+1 this is: c(nz) = -(Kv(nz+1) + Ty1) * zinv2 * ... Kv is already at the interface. Ty1 = Ki*slope² must be linearly interpolated from the two bracketing cell centres Z(nz) and Z(nz+1) to the interface at zbar(nz+1). The bug ------- Standard linear interpolation of a quantity f from cell centres Z(nz) and Z(nz+1) to the interface at zbar(nz+1) is: f* = f(nz) * (zbar(nz+1) - Z(nz+1)) / (Z(nz) - Z(nz+1)) + f(nz+1) * (Z(nz) - zbar(nz+1)) / (Z(nz) - Z(nz+1)) i.e. each cell gets a weight equal to the distance to the *other* cell divided by the total distance between cell centres. The cell closer to the interface receives the larger weight. The code had the two distance factors swapped: Ty1 = (Z(nz) - zbar(nz+1)) * zinv2 * slope²(nz) * Ki(nz) ! WRONG + (zbar(nz+1) - Z(nz+1) ) * zinv2 * slope²(nz+1)* Ki(nz+1) ! WRONG giving Ki(nz) the weight of the *distance from the upper cell to the interface* instead of the distance from the interface to the lower cell. Concrete example (Z(nz)=-10m, zbar(nz+1)=-15m, Z(nz+1)=-30m, total=20m): correct weight on Ki(nz) : (zbar-Z(nz+1))/total = 15/20 = 0.75 code weight on Ki(nz) : (Z(nz)-zbar )/total = 5/20 = 0.25 <- wrong On a uniform vertical grid both distances are equal (hnode/2 each), so both weights evaluate to 0.5 and the error cancels exactly. The bug is therefore invisible in idealised test cases with uniform vertical spacing and only manifests on stretched grids. Impact ------ All production FESOM2 configurations use a stretched vertical grid (high resolution near the surface, coarser below). At each interface the interpolated K_33 is biased toward whichever adjacent cell is thicker, i.e. toward the coarser, deeper cell. This systematically overestimates K_33 at shallow interfaces and underestimates it at depth, effectively rotating the Redi diffusion tensor away from the isopycnal surface and introducing spurious diapycnal mixing — exactly the flux the Redi scheme is designed to suppress. The error is largest where the vertical grid changes most rapidly and where isopycnal slopes are steep: shelf breaks, the Antarctic Circumpolar Current, and the upper ocean seasonal thermocline. Note: the Redi scheme is not active by default. Impact is limited to runs with isopycnal (Redi/GM) mixing enabled. Fix --- Swap the two distance factors in all three instances of Ty/Ty1 in diff_ver_part_impl_ale (surface layer nzmin, interior loop nzmin+1 to nzmax-2, and bottom layer nzmax-1), so each cell centre is weighted by its distance to the opposite cell rather than to the interface: Ty1 = (zbar(nz+1) - Z(nz+1)) * zinv2 * slope²(nz) * Ki(nz) + (Z(nz) - zbar(nz+1)) * zinv2 * slope²(nz+1)* Ki(nz+1) No other logic is changed. On a uniform grid the fix is a no-op. --- src/oce_ale_tracer.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 9c4fb7a37..955c60887 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -952,8 +952,8 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) zinv=1.0_WP*dt ! no .../(zbar(1)-zbar(2)) because of ALE ! calculate isoneutral diffusivity : Kd*s^2 --> K_33 = Kv + Kd*s^2 - Ty1= (Z_n(nz) -zbar_n(nz+1))*zinv2 *slope_tapered(3,nz ,n)**2*Ki(nz ,n) + & - (zbar_n(nz+1)-Z_n( nz+1))*zinv2 *slope_tapered(3,nz+1,n)**2*Ki(nz+1,n) + Ty1= (zbar_n(nz+1)-Z_n( nz+1))*zinv2 *slope_tapered(3,nz ,n)**2*Ki(nz ,n) + & + (Z_n(nz) -zbar_n(nz+1))*zinv2 *slope_tapered(3,nz+1,n)**2*Ki(nz+1,n) Ty1=Ty1*isredi ! layer dependent coefficients for for solving dT(1)/dt+d/dz*K_33*d/dz*T(1) = ... @@ -985,10 +985,10 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) ! 1/dz(nz) zinv2=1.0_WP/(Z_n(nz)-Z_n(nz+1)) ! calculate isoneutral diffusivity : Kd*s^2 --> K_33 = Kv + Kd*s^2 - Ty = (Z_n(nz-1 )-zbar_n(nz ))*zinv1 *slope_tapered(3,nz-1,n)**2*Ki(nz-1,n)+ & - (zbar_n(nz )-Z_n(nz ))*zinv1 *slope_tapered(3,nz ,n)**2*Ki(nz ,n) - Ty1= (Z_n(nz )-zbar_n(nz+1))*zinv2 *slope_tapered(3,nz ,n)**2*Ki(nz ,n)+ & - (zbar_n(nz+1)-Z_n(nz+1 ))*zinv2 *slope_tapered(3,nz+1,n)**2*Ki(nz+1,n) + Ty = (zbar_n(nz )-Z_n(nz ))*zinv1 *slope_tapered(3,nz-1,n)**2*Ki(nz-1,n)+ & + (Z_n(nz-1 )-zbar_n(nz ))*zinv1 *slope_tapered(3,nz ,n)**2*Ki(nz ,n) + Ty1= (zbar_n(nz+1)-Z_n(nz+1 ))*zinv2 *slope_tapered(3,nz ,n)**2*Ki(nz ,n)+ & + (Z_n(nz )-zbar_n(nz+1))*zinv2 *slope_tapered(3,nz+1,n)**2*Ki(nz+1,n) Ty =Ty *isredi Ty1=Ty1*isredi @@ -1028,8 +1028,8 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) zinv=1.0_WP*dt ! no ... /(zbar(nzmax-1)-zbar(nzmax)) because of ale ! calculate isoneutral diffusivity : Kd*s^2 --> K_33 = Kv + Kd*s^2 - Ty= (Z_n(nz-1) -zbar_n(nz)) * zinv1 * slope_tapered(3,nz-1,n)**2 * Ki(nz-1,n) + & - (zbar_n(nz)-Z_n(nz) ) * zinv1 * slope_tapered(3,nz ,n)**2 * Ki(nz,n) + Ty= (zbar_n(nz)-Z_n(nz) ) * zinv1 * slope_tapered(3,nz-1,n)**2 * Ki(nz-1,n) + & + (Z_n(nz-1) -zbar_n(nz)) * zinv1 * slope_tapered(3,nz ,n)**2 * Ki(nz,n) Ty =Ty *isredi ! layer dependent coefficients for for solving dT(nz)/dt+d/dz*K_33*d/dz*T(nz) = ...