configure_ocean_forcing Subroutine

   subroutine configure_ocean_forcing(cfg, ocean_state, grid, compute_rank, decomp)
      !! Wire surface wind stress, horizontal-viscosity coefficients, and the
      !! Coriolis beta-plane + PV-scheme variant from cfg into the ocean slots.
      !!
      !! When `decomp` is present the wind-stress setters receive the rank's
      !! global meridional offset and global domain extent so each subdomain
      !! seeds the correct portion of the global wind profile.  Absent ⇒
      !! single-rank byte-identical path.
      type(config_t), intent(in) :: cfg
      type(ocean_state_t), intent(inout) :: ocean_state
      type(hgrid_t), intent(in) :: grid
      integer, intent(in) :: compute_rank
      type(decomp_t), intent(in), optional :: decomp

      integer :: joff, nyg

      joff = 0
      nyg = grid%ny_phys
      if (present(decomp)) then
         joff = decomp%j_start - 1
         nyg = decomp%ny_global
      end if

      ! Wind stress: dispatch on cfg%ocean%topo%wind_config.
      select case (trim(cfg%ocean%topo%wind_config))
      case ("2gyre")
         call ocean_state%surface_stress%set_wind_stress_2gyre( &
            grid, cfg%ocean%topo%taux_magnitude, &
            j_offset=joff, ny_global=nyg)
      case ("neverworld2")
         call ocean_state%surface_stress%set_wind_stress_neverworld2( &
            grid, cfg%ocean%topo%taux_magnitude, &
            j_offset=joff, ny_global=nyg)
      case default  ! "constant"
         ! Warn loudly when a wind was configured on the knob this path does
         ! NOT read. `taux_magnitude` is only consumed by the "2gyre" and
         ! "neverworld2" profiles above; the constant path reads
         ! `&physics_nml wind_stress_x/y`. Setting the former under
         ! `wind_config="constant"` is silently inert — which is exactly how
         ! `coriolis_coast.nml` shipped a case whose header promises a
         ! "wind-spun-up gyre" that it has never generated: it ran its full
         ! 10 days at En = 0.000E+00.
         if (abs(cfg%ocean%topo%taux_magnitude) > 0.0_wp .and. &
             abs(cfg%wind_stress_x) <= 0.0_wp .and. &
             abs(cfg%wind_stress_y) <= 0.0_wp) then
            call logger%warning( &
               "&ocean_topo_nml taux_magnitude = "// &
               to_string(cfg%ocean%topo%taux_magnitude)//" is INERT under "// &
               'wind_config="constant" (it is only read by "2gyre" / '// &
               '"neverworld2"), and &physics_nml wind_stress_x/y are both '// &
               "zero, so this run has NO wind forcing at all. Set "// &
               "&physics_nml wind_stress_x to apply a constant wind.")
         end if
         call ocean_state%surface_stress%set_wind_stress_const( &
            cfg%wind_stress_x, cfg%wind_stress_y)
      end select

      ! Make the freshly-seeded wind valid in the ghost bands, then derive
      ! stress_mag from it.  Host-side (device_resident=.false.): this runs
      ! in the driver's configure phase, ahead of ocean_state_enter_data.
      !
      ! The exchange matters even though the wind is configure-static.
      ! `set_wind_stress_2gyre` / `_neverworld2` fill PHYSICAL rows only and
      ! leave ghost rows at zero ("so wall faces see no spurious stress"),
      ! which is right at a true wall and wrong at an MPI seam — there the
      ! ghost belongs to the neighbour rank and must carry its wind, not 0.
      ! Nothing else ever repairs it, and both `stress_mag` (KPP/EPBL u_*)
      ! and the MLE corner average read one cell beyond their own.  At a
      ! true wall the exchange is a no-op, so the zero-ghost intent of the
      ! analytic setters is preserved exactly (and `test_wind_2gyre_ghosts`
      ! still holds).
      call ocean_seam_refresh_surface_stress(ocean_state%surface_stress, grid, &
                                             ocean_state%bc, device_resident=.false.)

      ! Horizontal momentum viscosity (default 0 keeps analytical tests bit-
      ! identical; realistic wind-driven runs need a non-zero value).
      ocean_state%hvisc%nu_h = cfg%ocean%hvisc%nu_h
      ocean_state%hvisc%nu_4 = cfg%ocean%hvisc%nu_4
      ocean_state%hvisc%stress_tensor = cfg%ocean%hvisc%stress_tensor
      ocean_state%hvisc%bound_coef = cfg%ocean%hvisc%bound_coef
      ocean_state%hvisc%bound_kh = cfg%ocean%hvisc%bound_kh
      ocean_state%hvisc%no_slip = cfg%ocean%hvisc%no_slip
      ! Fail-loud guard (2026-07-28 dt=800 forensics): bound_kh clamps the
      ! effective per-cell viscosity to
      ! kh_max = bound_coef·0.125/(dt·(1/dx² + 1/dy²)), so a configured
      ! nu_h far above kh_max is silently unreachable and the run is
      ! under-damped relative to its nml intent — the dt=800 blow-up
      ! family was exactly this (bound_coef=0.15 capped the effective nu
      ! at ~177 m²/s against nu_h=10000).  Warning, not error: flow-aware
      ! closures legitimately over-provision constant floors/ceilings.
      ! CARTESIAN only: `grid%dx`/`grid%dy` are the cell size there, but a
      ! placeholder on a curvilinear grid (1 m on a supergrid, degrees on a
      ! spherical sector) — which made the estimate ~1e-5 m2/s and the
      ! warning fire, falsely, on every global run.  Curvilinear grids are
      ! covered per wet cell by the viscous-CFL stability audit.
      if (cfg%ocean%hvisc%bound_kh .and. cfg%dt_fixed > 0.0_wp &
          .and. parse_grid_config(cfg%ocean%grid%grid_config) == GRID_CONFIG_CARTESIAN &
          .and. grid%dx > 0.0_wp .and. grid%dy > 0.0_wp) then
         block
            real(wp) :: kh_max_est
            kh_max_est = cfg%ocean%hvisc%bound_coef*0.125_wp/ &
                         (cfg%dt_fixed*(1.0_wp/grid%dx**2 + 1.0_wp/grid%dy**2))
            if (cfg%ocean%hvisc%nu_h > 10.0_wp*kh_max_est .and. compute_rank == 0) then
               call logger%warning("hvisc: nu_h = "//to_string(cfg%ocean%hvisc%nu_h)// &
                                   " m2/s is unreachable — bound_kh clamps the effective "// &
                                   "viscosity to ~"//to_string(kh_max_est)// &
                                   " m2/s at this dx/dt (bound_coef = "// &
                                   to_string(cfg%ocean%hvisc%bound_coef)// &
                                   ").  The run will be under-damped relative to the nml "// &
                                   "intent; raise bound_coef (MOM6 default 0.8), reduce dt, "// &
                                   "or lower nu_h to what the clamp admits.")
            end if
         end block
      end if
      ! Anisotropic viscosity (Smith & McWilliams 2003) — stress-tensor
      ! path only.  Default kh_aniso=0 ⇒ isotropic, bit-identical.  The
      ! constant direction tensor is precomputed once from aniso_dir
      ! (mode 0 = grid-relative).  `aniso_mode /= 0` is rejected at
      ! configure (validate_config); the else below is a defensive guard.
      ocean_state%hvisc%kh_aniso = cfg%ocean%hvisc%kh_aniso
      if (cfg%ocean%hvisc%aniso_mode == 0) then
         call ocean_hvisc_set_aniso_direction(ocean_state%hvisc, &
                                              cfg%ocean%hvisc%aniso_dir(1), &
                                              cfg%ocean%hvisc%aniso_dir(2))
      else
         call ocean_hvisc_set_aniso_direction(ocean_state%hvisc, 1.0_wp, 0.0_wp)
      end if

      ! Coriolis: beta = 0 → uniform f_0 = cfg%coriolis_f; non-zero beta engages
      ! f(y) = f_0 + beta*(y - y_ref).  Filled by the single generator-driven
      ! routine (D7) — beta_plane is bit-identical to the legacy setter; the
      ! `f_0`/`beta` diagnostics are recorded here (the generator only fills
      ! the array).
      ocean_state%coriolis_adv%f_0 = cfg%coriolis_f
      ocean_state%coriolis_adv%beta = cfg%ocean%topo%coriolis_beta
      call fill_coriolis_corner(cfg, ocean_state%metrics, grid, &
                                ocean_state%coriolis_adv%f_corner)
      ! Tripolar: fold the STATIC f_corner once at configure (it is read
      ! straight from f_corner in the BT substep + Coriolis-adv, never
      ! re-wrapped per step).  f is a SCALAR under the fold — the reflected
      ! corner sits at the SAME latitude, so it COPIES (negate=.false.), no
      ! sign flip.  Periodic-x of the corner columns first (Appendix A
      ! ordering), then the north fold.  No-op for non-tripolar runs.
      if (trim(cfg%ocean%bc%north) == "tripolar_fold") then
         call fill_f_corner_seam_ghosts(grid, ocean_state%coriolis_adv%f_corner)
      end if
      ocean_state%coriolis_adv%pv_variant = parse_pv_variant(cfg%ocean%coriolis%form)
      ocean_state%coriolis_adv%use_hk_correction = &
         (ocean_state%coriolis_adv%pv_variant == PV_VARIANT_SADOURNY_HK)
      ocean_state%coriolis_adv%pv_adv_scheme = &
         parse_pv_adv_scheme(cfg%ocean%coriolis%pv_adv_scheme)
      ! HK pair floor on the coordinates with STATIC bed fillers, closed
      ! faces or not (see `coriolis_adv_t%hk_pair_floor`).
      ocean_state%coriolis_adv%hk_pair_floor = &
         any(parse_ocean_vcoord_type(cfg%vcoord_type) == &
             [VCOORD_Z_FIXED, VCOORD_ZSTAR, VCOORD_ZSTAR_FULL])
      ! Coastal lateral BC (land mask, C1): free-slip default, shared knob.
      ocean_state%coriolis_adv%no_slip = cfg%ocean%hvisc%no_slip
      ! Mass-consistent CorAdCalc (MOM6 parity); config fail-loud
      ! guarantees form="sadourny_energy" + split_scheme="pred_corr".
      ocean_state%coriolis_adv%state_fluxes = cfg%ocean%coriolis%use_state_fluxes
      ! BOUND_CORIOLIS velocity-form clamp (energy scheme only; config
      ! fail-loud guarantees form="sadourny_energy").
      ocean_state%coriolis_adv%bound_coriolis = cfg%ocean%coriolis%bound_coriolis
      ! PV corner-thickness construction (energy scheme only, config fail-loud).
      if (trim(cfg%ocean%coriolis%corner_h) == "mom6_area") then
         ocean_state%coriolis_adv%corner_h_variant = CORNER_H_MOM6_AREA
      else
         ocean_state%coriolis_adv%corner_h_variant = CORNER_H_CELL_MEAN
      end if
      if (compute_rank == 0) then
         call logger%info("Coriolis form:    "//trim(cfg%ocean%coriolis%form))
         if (cfg%ocean%coriolis%use_state_fluxes) then
            call logger%info("Coriolis CorAdv:  mass-consistent state fluxes "// &
                             "in the mom6 corrector (MOM6 uh/vh parity)")
         end if
         if (cfg%ocean%coriolis%bound_coriolis) then
            call logger%info("Coriolis bound:   BOUND_CORIOLIS ON — energy-scheme accel "// &
                             "clamped to the (f+zeta)*v velocity-form range (MOM6)")
         end if
         if (trim(cfg%ocean%coriolis%corner_h) == "mom6_area") then
            call logger%info("Coriolis corner:  PV corner-h = mom6_area (MOM6 "// &
                             "Area_q/(hArea_q+vol_neglect); round-off-identical to "// &
                             "cell_mean above H_MIN_PV)")
         end if
      end if
   end subroutine configure_ocean_forcing