Run all configure-time stability checks. Must run AFTER
configure_ocean_metrics + configure_ocean_land_mask (needs
the real filled ocean_metrics_t), before ocean_state_enter_data.
ierr present -> OCEAN_STATUS_ERR_SETUP on any hard-bound
violation (never error stop); warnings always just log,
whatever ierr does. Rank-0-only logging (mirrors every other
configure_ocean_* info/warning line); the ERROR path itself
always fires (every rank must agree the config is broken).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(config_t), | intent(in) | :: | cfg | |||
| type(ocean_metrics_t), | intent(in) | :: | metrics | |||
| type(hgrid_t), | intent(in) | :: | grid | |||
| integer, | intent(in) | :: | rank | |||
| integer, | intent(out), | optional | :: | ierr | ||
| real(kind=wp), | intent(in), | optional | :: | column(:,:) |
Reference column thickness (m) at T cells, ghosts included —
|
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| real(kind=wp), | private | :: | ah_max | ||||
| logical, | private | :: | auto_protected | ||||
| real(kind=wp), | private | :: | dt | ||||
| real(kind=wp), | private | :: | dt_therm | ||||
| real(kind=wp), | private | :: | dx_min | ||||
| integer, | private | :: | i_at | ||||
| logical, | private | :: | is_x | ||||
| integer, | private | :: | j_at | ||||
| real(kind=wp), | private | :: | kappa_h | ||||
| real(kind=wp), | private | :: | nu_h |
subroutine ocean_stability_audit(cfg, metrics, grid, rank, ierr, column) !! Run all configure-time stability checks. Must run AFTER !! `configure_ocean_metrics` + `configure_ocean_land_mask` (needs !! the real filled `ocean_metrics_t`), before `ocean_state_enter_data`. !! `ierr` present -> `OCEAN_STATUS_ERR_SETUP` on any hard-bound !! violation (never `error stop`); warnings always just log, !! whatever `ierr` does. Rank-0-only logging (mirrors every other !! `configure_ocean_*` info/warning line); the ERROR path itself !! always fires (every rank must agree the config is broken). type(config_t), intent(in) :: cfg type(ocean_metrics_t), intent(in) :: metrics type(hgrid_t), intent(in) :: grid integer, intent(in) :: rank integer, intent(out), optional :: ierr real(wp), intent(in), optional :: column(:, :) !! Reference column thickness (m) at T cells, ghosts included — !! `bt_work%bt_H_ref`, which is `b - z_draft` afloat and `0` !! where grounded. Present ⇒ the terrain-following stiffness !! check (Check 5) runs; absent ⇒ it is skipped, which is what a !! caller with no barotropic datum yet should do. real(wp) :: dx_min, dt, dt_therm, nu_h, kappa_h, ah_max integer :: i_at, j_at logical :: is_x logical :: auto_protected if (present(ierr)) ierr = OCEAN_STATUS_OK dt = cfg%dt_fixed nu_h = cfg%ocean%hvisc%nu_h kappa_h = cfg%ocean%hdiff%kappa_h ah_max = cfg%ocean%hvisc%ah_max ! dt_fixed <= 0 => adaptive-CFL dt, not known at configure time ! (matches the pre-existing kappa_h check's own gate). if (dt <= 0.0_wp) return call ocean_stability_min_cell(metrics, grid, dx_min, i_at, j_at, is_x) if (dx_min <= 0.0_wp .or. dx_min == huge(1.0_wp)) return ! ---- Check 1: viscous CFL (nu_h*dt/dx_min^2) ---- auto_protected = cfg%ocean%hvisc%bound_kh .or. cfg%ocean%hvisc%stress_tensor if (nu_h > 0.0_wp) then block real(wp) :: cfl, nu_h_max, dt_max cfl = ocean_viscous_cfl_number(nu_h, dt, dx_min) if (cfl > VISCOUS_CFL_LIMIT) then nu_h_max = ocean_viscous_cfl_max_nu_h(dt, dx_min, VISCOUS_CFL_LIMIT) dt_max = ocean_viscous_cfl_max_dt(nu_h, dx_min, VISCOUS_CFL_LIMIT) if (auto_protected) then if (rank == 0) then call logger%warning("&ocean_hvisc_nml nu_h = "//to_string(nu_h)// & " with dt = "//to_string(dt)//"s gives viscous CFL "// & to_string(cfl)//" at the smallest cell (dx = "// & to_string(dx_min/1000.0_wp)//" km, near i="//to_string(i_at)// & " j="//to_string(j_at)//"); limit is "// & to_string(VISCOUS_CFL_LIMIT)//" — but bound_kh/stress_tensor "// & "is enabled, so the runtime per-cell clamp will keep the "// & "EFFECTIVE viscosity within bound. Informational only.") end if else call fail("&ocean_hvisc_nml nu_h = "//to_string(nu_h)// & " with dt = "//to_string(dt)//"s gives viscous CFL "// & to_string(cfl)//" at the smallest cell (dx = "// & to_string(dx_min/1000.0_wp)//" km, near i="//to_string(i_at)// & " j="//to_string(j_at)//"); limit is "//to_string(VISCOUS_CFL_LIMIT)// & ". Set &ocean_hvisc_nml bound_kh = .true. for a per-cell clamp, "// & "or reduce nu_h below "//to_string(nu_h_max)// & ", or reduce dt below "//to_string(dt_max)//"s.", ierr, OCEAN_STATUS_ERR_SETUP) return end if end if end block end if ! ---- Check 2: tracer diffusive number (kappa_h) ---- ! Absorbed + fixed from the old rdb_config.F90 check (Cartesian-only, ! nominal dx/dy): now uses the REAL minimum cell, on ANY grid type. if (kappa_h > 0.0_wp) then block real(wp) :: dnum dt_therm = dt*real(cfg%ocean%vmix%dt_therm_ratio, wp) dnum = ocean_diffusive_number(kappa_h, dt_therm, dx_min) if (dnum > DIFFUSIVE_NUMBER_LIMIT) then call fail("&ocean_hdiff_nml kappa_h = "//to_string(kappa_h)// & " violates the explicit forward-Euler diffusive stability bound "// & "kappa_h*dt_therm*(1/dx_min^2+1/dy_min^2) <= "// & to_string(DIFFUSIVE_NUMBER_LIMIT)//" (computed "//to_string(dnum)// & " using dt_therm = "//to_string(dt_therm)//"s, dx_min = "// & to_string(dx_min)//"m near i="//to_string(i_at)//" j="//to_string(j_at)// & " — the ACTUAL smallest cell, not nominal dx/dy). Reduce kappa_h "// & "below "//to_string(DIFFUSIVE_NUMBER_LIMIT*dx_min**2/(2.0_wp*dt_therm))// & ", or reduce dt/dt_therm_ratio.", ierr, OCEAN_STATUS_ERR_SETUP) return end if end block end if ! ---- Check 3: Munk-layer resolution (WARNING) ---- if (nu_h > 0.0_wp) then block real(wp) :: ratio_min, beta_at, dx_at, nu_h_req integer :: j_munk call ocean_munk_worst_case(cfg, metrics, grid, nu_h, ratio_min, beta_at, dx_at, j_munk) if (ratio_min < MUNK_MIN_CELLS .and. ratio_min < huge(1.0_wp)) then nu_h_req = ocean_munk_required_nu_h(beta_at, dx_at, MUNK_MIN_CELLS) if (rank == 0) then call logger%warning("Munk sidewall boundary layer under-resolved: "// & "delta_M = (nu_h/beta)^(1/3) = "//to_string(ocean_munk_delta_m(nu_h, beta_at))// & "m spans only "//to_string(ratio_min)//" cells (dx = "// & to_string(dx_at)//"m, beta = "//to_string(beta_at)//" 1/(m*s) near j="// & to_string(j_munk)//"); need >= "//to_string(MUNK_MIN_CELLS)// & " cells or the wall carries grid-scale (2-delta) noise. Raise "// & "&ocean_hvisc_nml nu_h to at least "//to_string(nu_h_req)// & " m2/s (and raise ah_max to match — see the next check).") end if end if end block end if ! ---- Check 4: ah_max clamping nu_h inert (WARNING) ---- if (nu_h > 0.0_wp .and. ah_max > 0.0_wp .and. ah_max < nu_h) then if (rank == 0) then call logger%warning("&ocean_hvisc_nml ah_max = "//to_string(ah_max)// & " is BELOW nu_h = "//to_string(nu_h)//" — ah_max is a hard CAP on the "// & "per-face viscosity, so the configured nu_h is silently clamped down to "// & "ah_max and never actually applied. Set ah_max >= nu_h (e.g. ah_max = "// & to_string(nu_h)//" or higher).") end if end if ! ---- Check 5: terrain-following stiffness rx0 (WARNING) ---- ! Geometry only — no dt, no nz, no viscosity. A sigma-family ! coordinate divides the COLUMN into nz layers, so a big column ! contrast across one face offsets the two columns' K-th interfaces ! by Delta_e ~ rx0*(H_a+H_b), and the pressure-gradient truncation ! goes as Delta_e^3. See SIGMA_STIFFNESS_LIMIT for the citations and ! for the ISOMIP+ Ocean0 failure that motivated this. if (present(column)) then if (ocean_vcoord_is_terrain_following( & parse_vcoord_type(cfg%vcoord_type, VCOORD_SIGMA))) then block real(wp) :: rx0_max, h_thin, h_thick integer :: i_rx, j_rx, n_over, n_face logical :: rx_is_x call ocean_sigma_stiffness_worst( & size(column, 1), size(column, 2), & grid%nghost + 1, grid%nghost + grid%nx_phys, & grid%nghost + 1, grid%nghost + grid%ny_phys, & column, metrics%wet_T, rx0_max, i_rx, j_rx, rx_is_x, & h_thin, h_thick, n_over, n_face) if (rx0_max > SIGMA_STIFFNESS_LIMIT .and. rank == 0) then call logger%warning("Terrain-following stiffness rx0 = "// & to_string(rx0_max)//" exceeds "// & to_string(SIGMA_STIFFNESS_LIMIT)//" on the '"// & trim(adjustl(cfg%vcoord_type))// & "' vertical coordinate: a "//to_string(h_thin)// & " m column sits beside a "//to_string(h_thick)// & " m one across a single "// & merge("x", "y", rx_is_x)//" face near i="// & to_string(i_rx)//" j="//to_string(j_rx)//" ("// & to_string(n_over)//" of "//to_string(n_face)// & " wet-wet faces are over the bound). The sigma "// & "pressure-gradient truncation scales as the CUBE of the "// & "interface offset between neighbouring columns, so this "// & "is a LARGE spurious force, and a quiescent or long run "// & "over this geometry will measure it rather than the "// & "physics. A FORCED, viscous run is normally fine here — a "// & "constant lateral-viscosity floor arrests a steady spurious "// & "force rather than integrating it, which is what the shipped "// & "nu_h = 10000 double-gyre and neverworld2 cases rely on. A "// & "quiescent, weakly-damped or long spin-up run is NOT: smooth "// & "the topography to rx0 <= "//to_string(SIGMA_STIFFNESS_LIMIT)// & " (under an ice shelf, raising &ocean_cavity_dyn_nml "// & "h_min_cavity removes the thinnest columns), or use a "// & "coordinate whose interfaces do not follow the topography. "// & "Viscosity and a smaller dt DELAY a runaway, they do not "// & "remove the error.") end if end block end if end if if (present(ierr)) ierr = OCEAN_STATUS_OK end subroutine ocean_stability_audit