Skip to content

workbench_idemix2: IDEMIX2 - Internal Wave Energy Parameterisation for FESOM2#927

Draft
patrickscholz wants to merge 73 commits into
mainfrom
workbench_idemix2
Draft

workbench_idemix2: IDEMIX2 - Internal Wave Energy Parameterisation for FESOM2#927
patrickscholz wants to merge 73 commits into
mainfrom
workbench_idemix2

Conversation

@patrickscholz

Copy link
Copy Markdown
Contributor

IDEMIX2: Internal Wave Energy Parameterisation for FESOM2

IDEMIX2 (Internal wave Dissipation, Energy and MIXing, version 2) is a prognostic parameterisation of internal gravity wave energy and its contribution to diapycnal mixing, based on Olbers & Eden (2013) and Eden et al. (2014). It is implemented here as an extension to the existing IDEMIX1/TKE mixing scheme in FESOM2.

Concept

Unlike IDEMIX1, which tracks a single vertically integrated internal wave energy field, IDEMIX2 resolves the horizontal propagation direction of internal waves through a spectral angular dimension (nfbin directional bins, default 50 interior + 2 ghost bins covering 0–2π). This captures anisotropic wave fields near generation sites such as rough topography and coastlines, and allows wave energy to be tracked as it propagates across ocean basins.

Wave energy compartments

Two spectral wave energy fields are evolved independently:

  Surface wind forcing                 M2 tidal generation
  (NIW, near-inertial)                 (rough topography)
          │                                    │
          ▼                                    ▼
  ┌───────────────┐                  ┌─────────────────┐
  │  E_niw(φ, x)  │                  │   E_M2(φ, x)    │
  │  Near-inertial│                  │   M2 tidal IW   │
  │  waves        │                  │   (~12.4 h)     │
  └───────┬───────┘                  └────────┬────────┘
          │  ω ≈ f (lat-dependent)            │  ω = ω_M2 (fixed)
          └──────────────┬────────────────────┘
                         │  wave-wave interaction
                         ▼
               ┌──────────────────┐
               │   E_iw(z, x)     │   background IW continuum
               │   (vertically    │   (GM spectrum)
               │    resolved)     │
               └────────┬─────────┘
                        │  dissipation → diapycnal mixing
                        ▼
               ┌──────────────────┐
               │   TKE scheme     │   iw_diss added as
               │   Kv, Av         │   energy source term
               └──────────────────┘

Timestep structure

Each call to calc_cvmix_idemix2 follows this sequence:

  ┌─────────────────────────────────────────────────────────────────┐
  │ 1. Node parameters                                              │
  │    cn (baroclinic wave speed), α_c, v0, c0, τ_M2, τ_niw,        │
  │    vertical structure functions  [loop: node, OMP parallel]     │
  ├─────────────────────────────────────────────────────────────────┤
  │ 2. Element group velocities                                     │
  │    u_M2(φ), v_M2(φ), w_M2(φ)  — spectral propagation speeds     │
  │    scatter elem → node (area-weighted)                          │
  ├─────────────────────────────────────────────────────────────────┤
  │ 3. Horizontal spectral integration  [for M2, then NIW]          │
  │                                                                 │
  │   E^n  ──▶  superbee horizontal advection  ──▶  flx_uv(φ,edge)  │
  │              │                                                  │
  │              ├──▶  reflective coastal BC                        │
  │              │     (redirect toward-coast flux into             │
  │              │      mirror spectral bin at interior node)       │
  │              ▼                                                  │
  │         superbee cross-spectral advection ──▶  flx_cs(φ,node)   │
  │              │                                                  │
  │              ▼                                                  │
  │         flux → divergence  +  reflect_src injection             │
  │              │                                                  │
  │              ▼                                                  │
  │         E^(n+1) = (E^n + dt·(Edivh + Edivs + forc))             │
  │                   / (1 + dt·τ)     [AB2 or Euler]               │
  ├─────────────────────────────────────────────────────────────────┤
  │ 4. Eiw vertical diffusion + implicit dissipation                │
  │    tridiagonal solve per node  [OMP parallel]                   │
  ├─────────────────────────────────────────────────────────────────┤
  │ 5. Eiw horizontal diffusion  (optional)                         │
  ├─────────────────────────────────────────────────────────────────┤
  │ 6. Wave-wave interaction  E_M2/niw ↔ E_iw                       │
  │    [loop: node, OMP parallel]                                   │
  ├─────────────────────────────────────────────────────────────────┤
  │ 7. iw_diss → TKE (Kv, Av)                                       │
  └─────────────────────────────────────────────────────────────────┘

Reflective coastal boundary condition

Wave energy reaching a coastline is redirected into the spectrally mirrored propagation direction at the last interior node, rather than leaking into land:

  LAND │ n_coast ──── edge ──── n_int (interior node)
       │
       │  toward-coast flux in bin φ_i  →  zeroed
       │                                       │
       │                                       ▼
       │                            energy re-injected at n_int
       │                            in mirror bin φ_k = edge normal
       │                            (net energy change = 0,
       │                             spectral redistribution only)

The mirror bin φ_k for each coast edge is pre-computed once at init from the edge geometry and stored in iwe2_refl_bin(edge).

Forcing inputs

Field Source Role
M2 tidal forcing Observed anisotropic tidal generation (modes 1+2) E_M2 surface source
NIW / surface forcing Wind-derived near-inertial flux E_niw surface source or E_iw surface BC
Bottom tidal forcing Near-field tidal dissipation dataset E_iw bottom source
Topographic HRMS / λ Seafloor roughness fields Scattering timescale τ

patrickscholz and others added 30 commits January 28, 2026 16:23
…sages if something is wrong in the forcing files
…ested yet, it is discarded for the moment for the idemix2 implementation
…ter, fix some bugs in the horizontal advection of E_niw
patrickscholz and others added 24 commits June 4, 2026 12:54
Remove annual_event debug prints (fixes #930)
…ng option Eiw_diss, equivalent to the elem-->node interpolation present in idemix1
…een by tke. It turns out that the idemix2 only mode (no EM2, Eniw) suffers from occasionally winter blowup hickups that require the smoothing of these fields to get the blowups away
…een by tke. It turns out that the idemix2 only mode (no EM2, Eniw) suffers from occasionally winter blowup hickups that require the smoothing of these fields to get the blowups away
…erface copies

Problem
-------
The idemix2_smooth_Eiw / _Eiw_diss / _alpha_c options were introduced to
suppress a recurring cflz blowup at ~72-75°N in February caused by the
following chain:

  deep winter mixed layer → Eiw builds to cap (Eiw_maxthresh = 0.1 m²/s²)
  → alpha_c spikes at ML base (quadratic dissipation: alpha_c × Eiw_old × Eiw_new)
  → Eiw_diss ≈ Eiw_maxthresh / dt ≈ 3.7e-5 m²/s³ at a single node
  → tke_in3d_iwdis → TKE source spike → Kv → cflz

The three smooth_nod calls were applied inside gen_modules_cvmix_idemix2.F90
directly to the prognostic state iwe2_E_iw(:,:,tip1) and to the diagnostic
fields iwe2_E_iw_diss and iwe2_alpha_c.  This suppressed the blowup but had
a catastrophic side effect: smooth_nod applied every timestep to a prognostic
field is mathematically equivalent to adding a horizontal diffusion term with
effective diffusivity

  D_eff ≈ Δx² / dt  ≈  (50 km)² / 2700 s  ≈  930 m²/s

The physical horizontal Eiw diffusion in the model is

  D_h = v0² × tau_h  ≈  (0.009 m/s)² × 1.296×10⁶ s  ≈  105 m²/s   (ML floor)
                      ≈  3×10⁴ m²/s                                  (pycnocline)

The artificial smoothing is therefore 10–10000× stronger than the physical
horizontal diffusion.  Over the years of spin-up integration, it diffuses Eiw
laterally over hundreds of kilometres, homogenising the field and reducing the
global Eiw integral by roughly an order of magnitude — well below IDEMIX1 —
and destroying all horizontal heterogeneity that IDEMIX2 is supposed to
produce from locally variable N² and bottom forcing patterns.

Root cause
----------
Applying smooth_nod to iwe2_E_iw(:,:,tip1) every timestep is cumulative:
the smoothed state becomes the initial condition for the next step, which is
smoothed again, and so on.  This accumulates over thousands of timesteps.

By contrast, iwe2_E_iw_diss and iwe2_alpha_c are diagnostic fields recomputed
from scratch from N² at every timestep.  Smoothing them in-place before TKE
reads them has no memory and does not accumulate — but smoothing them ALSO had
no effect on preventing blowup, because TKE reads these fields only after the
IDEMIX2 step is complete, at which point the spike in Eiw_diss originates from
the unsmoothed iwe2_E_iw.

Fix
---
Move all three smooth_nod calls from gen_modules_cvmix_idemix2.F90 into
gen_modules_cvmix_tke.F90, where they are applied to the TKE-interface copies
(tke_in3d_iwe, tke_in3d_iwdis, tke_in3d_iwealphac) immediately after these
arrays are filled from the IDEMIX2 module variables:

  tke_in3d_iwe       = iwe2_E_iw(:,:,iwe2_tip1)
  tke_in3d_iwdis     = -iwe2_E_iw_diss
  tke_in3d_iwealphac = iwe2_alpha_c
  if (idemix2_smooth_Eiw)      smooth_nod(tke_in3d_iwe,       ...)
  if (idemix2_smooth_Eiw_diss) smooth_nod(tke_in3d_iwdis,     ...)
  if (idemix2_smooth_alpha_c)  smooth_nod(tke_in3d_iwealphac, ...)

These three tke_in3d_* arrays are temporary diagnostic copies refilled from
scratch every timestep.  Smoothing them damps the Eiw_diss spike seen by TKE
without touching the prognostic iwe2_E_iw field.  There is no accumulation,
no effective diffusion of Eiw, and no memory across timesteps.

The prognostic Eiw field is never modified by smooth_nod under this scheme.

Result
------
- February 1959 cflz blowup at 72-75°N fully suppressed with
  idemix2_smooth_Eiw = .true., idemix2_smooth_Eiw_diss = .true. (niter=1)
- Global Eiw field recovers the horizontal heterogeneity of the unsmoothed
  IDEMIX2 run: spatial structure at ridge, boundary current and tidal
  hotspot scales is preserved
- IDEMIX2 Eiw consistently higher than IDEMIX1 over most of the ocean
  (idemix1 - idemix2 difference predominantly positive), as expected
  physically from the additional spectral structure in IDEMIX2
- All three smoothing flags (idemix2_smooth_Eiw, idemix2_smooth_Eiw_diss,
  idemix2_smooth_alpha_c) and idemix2_smooth_niter remain fully functional
  from the namelist; their semantics are unchanged
… mode

Bug
---
iwe2_E_iw_hdif (horizontal Eiw diffusion tendency) was always zero when
idemix2_enable_hor_diff_impl_iter = .true., even with idemix2_diag_Eiw = .true.

Root cause: Fortran aliasing in compute_hdiff_Eiw.  The iterative implicit
call passes the same array for both the Eiw_old (intent IN) and Eiw (intent
INOUT) arguments:

  call compute_hdiff_Eiw(iwe2_E_iw(:,:,iwe2_tip1),   ! Eiw_old
                          iwe2_E_iw(:,:,iwe2_tip1),   ! Eiw  — same array
                          ...)

The edge loop updates Eiw in-place.  The end-of-subroutine diagnostic then
computed Eiw - Eiw_old, but since both point to the same memory, all
differences were identically zero.  The explicit mode (Eiw_old = iwe2_ti,
Eiw = iwe2_tip1) would not have had this problem, but idemix2_enable_hor_-
diff_expl is always .false. in practice.

Fix
---
Removed the optional Eiw_dt and Eiw_hdif arguments from compute_hdiff_Eiw
entirely and deleted the aliasing-prone end-of-subroutine diagnostic block.
Diagnostics are now computed at the driver level using a before/after
difference on iwe2_E_iw(:,:,iwe2_tip1):

  ! before both diffusion modes:
  iwe2_E_iw_hdif = -iwe2_E_iw(:,:,iwe2_tip1)   ! save negated pre-state

  ! ... explicit and/or iterative compute_hdiff_Eiw calls ...

  ! after both diffusion modes:
  iwe2_E_iw_hdif = (iwe2_E_iw_hdif + iwe2_E_iw(:,:,iwe2_tip1)) / dt
  iwe2_E_iw_dt  +=  iwe2_E_iw_hdif

This gives iwe2_E_iw_hdif = (Eiw_after - Eiw_before) / dt in m²/s³,
correct for both explicit (single-pass) and iterative implicit modes, and
without any aliasing issue.

Refactor
--------
The previously separate sections 5 (explicit) and 6 (iterative) and their
duplicated if (idemix2_diag_Eiw) / else branches are consolidated into a
single section 5+6 with one shared before/after diagnostic wrapper.
compute_hdiff_Eiw is now a pure physics routine with no optional arguments.
Bug 1: compute_hdiff_Eiw always produced zero increments
---------------------------------------------------------
vol_wcelli (inverse nodal cell volume) was allocated and zeroed in
init_cvmix_idemix2 but compute_vol_wcell() was never called.  Every
flux increment in compute_hdiff_Eiw multiplies vol_wcelli:

  Eiw(nz,node) += dt * tau_h / niter * vol_wcelli(nz,node) * vflux

With vol_wcelli = 0 everywhere, horizontal diffusion was a silent
no-op regardless of Eiw gradients or group velocity.  The previous
commit (2ca4457) introduced the correct before/after diagnostic, but
the diagnostic faithfully reported zero because hdiff genuinely changed
nothing.

Fix: call compute_vol_wcell(vol_wcelli, mesh, partit) in calc_idemix2
at the start of each timestep.  The call must be per-timestep because
compute_vol_wcell uses Z_3d_n and hnode, which evolve under ALE.

Bug 2: wrong units in compute_Eiw_waveinteract diagnostics
-----------------------------------------------------------
M2_diss is a dissipation rate [m²/s³] — erroneously multiplied by dt,
giving units of [m²/s²].  IW_diss is an energy increment [m²/s²]
accumulated inside the timestep loop — erroneously stored as-is,
giving units of [m²/s²] instead of the required tendency [m²/s³].

Fix:
  E_M2_dt      -= M2_diss          (was: -= dt*M2_diss)
  E_M2_diss_wwi = -M2_diss          (was: = -dt*M2_diss)
  E_iw_dt      += IW_diss/dt        (was: += IW_diss)
  E_iw_diss_M2  = IW_diss/dt        (was: = IW_diss)
  E_iw_diss_niw = IW_diss/dt        (was: = IW_diss)  [NIW branch]

Bug 3: Edt tendency missing /dt in hsintegrate_Ecompart
--------------------------------------------------------
(E_tip1 - E_ti) is an energy difference [m²/s²]; dividing by dt gives
the tendency [m²/s³] that the Edt diagnostic variable represents.

  Edt = (E_tip1 - E_ti) / dt       (was: Edt = E_tip1 - E_ti)
…ostic units

Bug 1: E_M2_struct / E_niw_struct wrong level slice at wave-wave call sites
---------------------------------------------------------------------------
(gen_modules_cvmix_idemix2.F90, section 7, 8 call sites)

iwe2_E_M2_struct and iwe2_E_niw_struct are filled at levels (uln:nln, node) by
compute_vert_struct_fct, but were passed as (:, node) to compute_Eiw_waveinteract.
For nodes with uln > 1 (ice-shelf cavities, steep shelf), the subroutine read
level 1 (zero-initialised) as the first active level and missed the deepest active
level entirely, silently zeroing the wave-wave energy injection at those nodes.
Fixed by passing (uln:nln, node) at all 8 call sites.

Bug 2: tau_compart uninitialised before floor application
---------------------------------------------------------
(cvmix_idemix2.F90, compute_compart_interact_tscale)

tau_compart is declared intent(out); max(1/(50 days), tau_compart) on line 512 was
reached without ever assigning tau_compart when omega_compart <= |f| (M2 poleward
of ~74.5 N).  Fortran standard leaves intent(out) undefined on entry.  Fixed by
adding tau_compart = 0.0_cvmix_r8 before the conditional block.

Bug 3: wrong unit labels for IDEMIX1 tendency diagnostics
---------------------------------------------------------
(io_meandata.F90)

iwe_Ttot, iwe_Tsur, iwe_Tbot and iwe_Thdi were labeled 'm^2/s^2' (energy per unit
mass) instead of 'm^2/s^3' (tendency/production rate).  iwe_Tdif and iwe_Tdis were
already correct.  Fixed unit strings to 'm^2/s^3' for all four diagnostics.
Bug 1: NIW spectral refraction kdot inflated x1000 near equator
---------------------------------------------------------------
(cvmix_idemix2.F90, compute_compart_groupvel)

The phi_dot formula has two structurally different terms:
  1st: -cn / sqrt(omega^2-f^2) / omega * f * grad_f   <- singular as f->omega, needs floor
  2nd: -sqrt(omega^2-f^2) / omega * grad_cn            <- non-singular, goes to zero naturally

A single shared fxa = sqrt(max(1d-10, omega^2-f^2)) was used for both. The 1d-10
floor is correct for the 1st term (prevents 1/fxa divergence), but applying it to
the 2nd term inflates sqrt(fxa) near the equator: for NIW with omega_niw = 1.05|f|,
omega_niw^2-f^2 = 0.1025*f^2 -> 0 at equator, but the floor gives fxa = 1e-5 instead
of the physical ~0. With omega_niw = 1e-8 at the equator this means fxa/omega = 1000
instead of ~0 -- the cn-gradient kdot term was x1000 too large in the equatorial band.

Fix: split into fxa (floored, for the singular 1/fxa term) and fxa2 (physical zero,
max(0d0,...), for the sqrt(fxa) terms). Matches pyOM2's idemix_group_vel.f90 exactly,
which uses two separate fxa computations with different floors.

Bug 2: hard +/-5 deg equatorial mask on NIW velocities causing accumulation band
--------------------------------------------------------------------------------
(gen_modules_cvmix_idemix2.F90, section 2 group velocity loop)

NIW spatial and spectral velocities were zeroed for all elements with |lat| < 5 deg
to avoid the equatorial singularity. This created an energy flux wall at +/-5 deg:
NIW energy advecting equatorward from higher latitudes could not cross the boundary,
accumulating at the outer edge of the masked region. The accumulated E_niw drove
anomalously large iwe2_Eiw_dissniw = tau_niw * sint * E_niw_struct at the +/-5 deg
band, visible as sharp stripes in diagnostics at all depth levels.

The fix from Bug 1 (fxa2 -> 0 at equator) makes the mask unnecessary. pyOM2 has no
such mask; it relies solely on the 1d-10 floor in the singular 1/fxa term. Removed.

Bug 3: IDEMIX1 iwe_v0 diagnostic stream pointing to iwe_c0 array
-----------------------------------------------------------------
(io_meandata.F90, CASE 'IDEMIX    ')

The def_stream call for 'iwe_v0' (horizontal group velocity) passed iwe_c0(:,:)
as the data array instead of iwe_v0(:,:) -- a copy-paste error from the preceding
iwe_c0 registration. Since c0 << v0 (gofx2 << hofx2 for typical N/|f| >> 1), the
iwe_v0 output file was ~10-100x smaller than the physically correct values, making
IDEMIX1 v0 appear negligible compared to IDEMIX2 iwe2_v0. Fixed to iwe_v0(:,:).
@patrickscholz patrickscholz force-pushed the workbench_idemix2 branch 3 times, most recently from 25152b2 to dbfd056 Compare June 19, 2026 10:39
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants