seed_eady_ic Subroutine

   subroutine seed_eady_ic(state, grid, cfg, ierr)
      !! Eady-front overlay IC.  Assumes flat bottom + uniform layer
      !! split already in place from `ocean_state_seed_from_cfg`.
      !!
      !! Sets:
      !!   T(i, j, k) = T_ref + dT_dz · z(k) + dT_dy · (y_phys(j) - y_mid)
      !!                                                       + ε(i, j, k)
      !!   u_face_x(i, j, k) = dU/dz · (z(k) - z_mid)
      !! with `dU/dz` from thermal-wind balance:
      !!     dU/dz = g · alpha_T / (rho_0 · f) · dT_dy
      !!
      !! Convention: k=1 is the bed, k=nz the surface.  z(k=1) =
      !! -H + dz/2, z(k=nz) = -dz/2.  y_phys uses cell centres.
      !!
      !! The perturbation ε is a uniform-random ±eady_pert_amp/2 noise
      !! seeded from `cfg%ocean%ic%eady_pert_seed`, applied to interior cells
      !! (j ∈ [ng+2, ng+ny_phys-1]) so wall-adjacent rows stay clean.
      !!
      !! Requires `coriolis_f /= 0`.  Errors if topo_config is not "flat".
      !!
      !! MPI: `y_mid` (the front centre) uses the GLOBAL meridional extent
      !! `grid%ny_global`, and the local row index is mapped to a global
      !! physical position with `grid%j_offset_global`.  Built from the
      !! LOCAL extents instead, every rank would put the front in the middle
      !! of its own tile.  On a single rank the offset is 0 and global ==
      !! local, so the seed is byte-identical.
      !!
      !! Not decomposition-invariant: the `eady_pert_amp` noise is drawn
      !! per-rank on the LOCAL array, so its realisation depends on the rank
      !! layout (the deterministic Eady profile underneath does not).
      type(ocean_state_t), intent(inout) :: state
      type(hgrid_t), intent(in) :: grid
      type(config_t), intent(in) :: cfg
      integer, intent(out), optional :: ierr
         !! Non-zero on an Eady-IC configuration conflict when present;
         !! absent behaves as today (`error stop`).

      integer :: i, j, k, nx_total, ny_total, nz_ml, ng
      integer :: ny_phys, j_phys, joff
      real(wp) :: H, dz, y_phys, y_mid, z_k, z_mid
      real(wp) :: dUdz, T_local, hk
      real(wp), allocatable :: noise(:, :, :)
      integer :: idx_T, n_seed, seed_val
      integer, allocatable :: seed_buf(:)

      if (trim(cfg%ocean%topo%topo_config) /= "flat") then
         call fail("seed_eady_ic: requires topo_config='flat'", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if
      if (abs(cfg%coriolis_f) < tiny(1.0_wp)) then
         call fail("seed_eady_ic: requires coriolis_f /= 0", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if

      nx_total = size(state%barotropic%b, 1)
      ny_total = size(state%barotropic%b, 2)
      nz_ml = state%multilayer%nz_ml
      ng = grid%nghost
      ! LOCAL ny_phys — used only to place the noise-free wall-adjacent
      ! rows, which is a per-rank array-bound question, not a coordinate.
      ny_phys = grid%ny_phys
      joff = grid%j_offset_global
      idx_T = state%multilayer%idx_temperature

      H = cfg%ocean%topo%max_depth
      dz = H/real(nz_ml, wp)
      y_mid = 0.5_wp*real(grid%ny_global, wp)*grid%dy
      z_mid = -0.5_wp*H

      ! Thermal-wind balance: ∂u_g/∂z = (g/(fρ₀))·∂ρ/∂y, with linear
      ! EOS ρ = ρ₀ − α(T−T_ref) giving ∂ρ/∂y = −α·∂T/∂y.  Therefore
      !     dU/dz = -g·α/(f·ρ₀) · dT/dy
      ! (cold-to-north dT/dy < 0 → eastward shear with z, surface jet
      !  westerly — the canonical mid-latitude convention).
      !
      ! Geostrophic balance for the y-independent part of -f·u would
      ! also demand a linearly tilted SSH (∂η/∂y = f·dU/dz·z_mid/g).
      ! Setting that η here AND refreshing h_layer accordingly produces
      ! ghost-row h_layer mismatch and tracer-mass corruption in the
      ! continuity-PPM path — the analytic balance doesn't exactly
      ! match the FV-lite PGF's discrete vertical integral.  Left for
      ! follow-up; for now the IC is hydrostatic-only and IG waves
      ! shed in the first ~24h as the system geostrophically adjusts.
      dUdz = -GRAVITY*cfg%ocean%ic%alpha_T/(cfg%ocean%ic%rho_0*cfg%coriolis_f) &
             *cfg%ocean%ic%eady_dT_dy

      ! Pre-generate noise on host (RNG isn't device-callable).  Fill
      ! the full array including ghosts; we zero ghost + bordering rows
      ! after to keep the perturbation cleanly inside the physical
      ! interior.
      allocate (noise(nx_total, ny_total, nz_ml))
      call random_seed(size=n_seed)
      allocate (seed_buf(n_seed))
      seed_val = cfg%ocean%ic%eady_pert_seed
      seed_buf = [(seed_val + 17*j, j=1, n_seed)]
      call random_seed(put=seed_buf)
      call random_number(noise)
      noise = cfg%ocean%ic%eady_pert_amp*(noise - 0.5_wp)
      ! Zero ghosts + first/last physical row to avoid wall-adjacent
      ! noise that would project onto the gravest mode.
      noise(:, 1:ng + 1, :) = 0.0_wp
      noise(:, ng + ny_phys:, :) = 0.0_wp

      ! Override T(i,j,k) with the Eady analytical profile + noise.
      if (idx_T > 0) then
         do k = 1, nz_ml
            z_k = -H + (real(k, wp) - 0.5_wp)*dz
            do j = 1, ny_total
               ! GLOBAL physical row index of this local row.
               j_phys = j - ng + joff
               y_phys = (real(j_phys, wp) - 0.5_wp)*grid%dy
               T_local = cfg%ocean%ic%eady_T_ref + cfg%ocean%ic%eady_dT_dz*z_k &
                         + cfg%ocean%ic%eady_dT_dy*(y_phys - y_mid)
               do i = 1, nx_total
                  hk = state%multilayer%h_layer(i, j, k)
                  state%multilayer%tracers(idx_T)%hTr(i, j, k) = &
                     (T_local + noise(i, j, k))*hk
               end do
            end do
         end do
      end if

      ! Thermal-wind-balanced u_face_x(z).  U depends on z only (no y, x).
      do k = 1, nz_ml
         z_k = -H + (real(k, wp) - 0.5_wp)*dz
         state%multilayer%u_face_x_layer(:, :, k) = dUdz*(z_k - z_mid)
         ! hu on x-faces: interpolate cell-centred h to the u-face (end faces
         ! take the adjacent cell), mirroring the production convention in
         ! face_depth_mean_u.  u_face/h_layer are staggered (nx+1 faces vs nx
         ! cells), so the old whole-array `u_face * h_layer` was non-conformant
         ! (silently OOB on non-checking compilers; LFortran flags it).
         associate (uf => state%multilayer%u_face_x_layer, &
                    hl => state%multilayer%h_layer, &
                    huf => state%multilayer%hu_face_x_layer)
            huf(1, :, k) = uf(1, :, k)*hl(1, :, k)
            huf(2:nx_total, :, k) = uf(2:nx_total, :, k) &
                                    *0.5_wp*(hl(1:nx_total - 1, :, k) + hl(2:nx_total, :, k))
            huf(nx_total + 1, :, k) = uf(nx_total + 1, :, k)*hl(nx_total, :, k)
         end associate
      end do

      deallocate (noise, seed_buf)
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine seed_eady_ic