Fix inverted interpolation weights for Redi K_33 at vertical interfaces

[patrickscholz]

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.

[suvarchal]

is it another bug hunt of yours?

[patrickscholz]

its some older stuff from nikolay and sergey about some AI assisted code evaluation, that i remembered and work through now

[patrickscholz]

luckily the effect of this bug is not very large!!!

• Temperature bias

• Salinity bias

• MLD3 bias