seed_geostrophic_adjustment_ic Subroutine

   subroutine seed_geostrophic_adjustment_ic(state, grid, cfg, ierr)
      !! Rossby's classic geostrophic-adjustment problem.  Overlays a
      !! Gaussian SSH bump on a flat-bottom, single-layer (barotropic)
      !! state at rest:
      !!     η(x, y) = A · exp(-r² / L²)
      !!     h(i, j) = b(i, j) + η(i, j)
      !!     u = v = 0
      !! where r is the radial distance from the bump centre.
      !!
      !! Expected evolution (analytical):
      !!   * Length-scale ratio L / Rd with Rd = √(gH)/f governs the
      !!     split between radiated IG-wave energy and retained
      !!     balanced geostrophic ring.
      !!     L >> Rd → most energy retained, surface stays bumped with
      !!                 a balanced anticyclonic flow around it.
      !!     L << Rd → most energy radiates, the bump flattens out.
      !!     L ~ Rd  → partial, ~50/50.
      !!   * Time-series of η at the centre oscillates with period
      !!     2π/f, damping toward the residual value.
      !!   * IG-wave fronts propagate outward at c = √(gH).
      !!
      !! Diagnostic interpretation:
      !!   * Bump amplitude growing in time → numerical instability
      !!     (wrong PGF sign, bad continuity).
      !!   * Oscillation period != 2π/f → wrong Coriolis coupling.
      !!   * Energy not conserved (after damping out the IG transient)
      !!     → spurious sources / sinks in the bt step or Coriolis-adv.
      !!
      !! Requires `topo_config = "flat"` so the bump sits on a
      !! uniform reference depth.  Single-layer is recommended (nz=1)
      !! to isolate the barotropic adjustment; multi-layer works but
      !! the bump distributes uniformly across layers.
      !!
      !! MPI: the default bump centre uses the GLOBAL extents
      !! `grid%nx_global` / `grid%ny_global`, and local cell indices are
      !! mapped to global physical positions with `grid%i_offset_global` /
      !! `grid%j_offset_global`.  Built from the LOCAL extents instead,
      !! every rank would drop a full-amplitude bump in the middle of its
      !! own tile.  On a single rank the offsets are 0 and global == local,
      !! so the seed is byte-identical.
      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 a geostrophic-adjustment IC configuration conflict
         !! when present; absent behaves as today (`error stop`).

      integer :: i, j, nx_total, ny_total, nz_ml, ng
      integer :: i_phys, j_phys, ioff, joff
      real(wp) :: x_phys, y_phys, x_c, y_c, r2, eta_local, inv_L2

      if (trim(cfg%ocean%topo%topo_config) /= "flat") then
         call fail("seed_geostrophic_adjustment_ic: requires topo_config='flat'", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if
      if (cfg%ocean%ic%ga_length_scale <= 0.0_wp) then
         call fail("seed_geostrophic_adjustment_ic: ga_length_scale must be > 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
      ioff = grid%i_offset_global
      joff = grid%j_offset_global

      ! Default centre = GLOBAL basin midpoint (negative cfg values).
      x_c = cfg%ocean%ic%ga_x_center
      y_c = cfg%ocean%ic%ga_y_center
      if (x_c < 0.0_wp) x_c = 0.5_wp*real(grid%nx_global, wp)*grid%dx
      if (y_c < 0.0_wp) y_c = 0.5_wp*real(grid%ny_global, wp)*grid%dy

      inv_L2 = 1.0_wp/(cfg%ocean%ic%ga_length_scale*cfg%ocean%ic%ga_length_scale)

      ! Override h_total = b + η(x, y) and redistribute the per-layer
      ! thickness evenly.  Velocities + face fluxes left at zero from
      ! the default seed.  Velocities are unchanged at zero, so hu/hv
      ! stay at zero too.
      do j = 1, ny_total
         ! GLOBAL physical indices of this local cell.
         j_phys = j - ng + joff
         y_phys = (real(j_phys, wp) - 0.5_wp)*grid%dy
         do i = 1, nx_total
            i_phys = i - ng + ioff
            x_phys = (real(i_phys, wp) - 0.5_wp)*grid%dx
            r2 = (x_phys - x_c)**2 + (y_phys - y_c)**2
            eta_local = cfg%ocean%ic%ga_eta_amp*exp(-r2*inv_L2)
            state%barotropic%h(i, j) = state%barotropic%b(i, j) + eta_local
         end do
      end do

      call seed_h_layer_uniform_impl(state%multilayer%h_layer, &
                                     state%barotropic%h, nz_ml)
      ! Re-seed S, T so hTr/h_layer stays at the configured uniform
      ! values regardless of the η bump.
      if (state%multilayer%idx_salinity > 0) then
         call seed_tracer_uniform_impl( &
            state%multilayer%tracers(state%multilayer%idx_salinity)%hTr, &
            state%multilayer%h_layer, cfg%initial_salinity, nz_ml)
      end if
      if (state%multilayer%idx_temperature > 0) then
         call seed_tracer_uniform_impl( &
            state%multilayer%tracers(state%multilayer%idx_temperature)%hTr, &
            state%multilayer%h_layer, cfg%initial_temperature, nz_ml)
      end if
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine seed_geostrophic_adjustment_ic