seed_baroclinic_jet_ic Subroutine

   subroutine seed_baroclinic_jet_ic(state, grid, cfg, ierr)
      !! Two-layer reduced-gravity baroclinic-instability IC — a
      !! geostrophically-balanced tanh jet in the upper (surface) layer of
      !! a spherical re-entrant channel, plus a front-localised sech²
      !! meander seed that the instability grows.  Reproduces the `bc_inst`
      !! spec (SIM_DETAILS.md §5), mapped to Roundabout's bottom-up layer
      !! convention (k=1 bed, k=nz surface).  The spec's layers FLIP:
      !!
      !!   * spec upper layer 1 (the jet) → Roundabout k=2 (surface)
      !!   * spec lower layer 2 (at rest) → Roundabout k=1 (bed)
      !!
      !! With `y = R·(φ − φ₀)` (m) about the jet centre and
      !! `xm = R·cosφ₀·(λ − λ_west)` (m):
      !!
      !!   ξ(x,y) = Δξ·tanh(y/L) + A_pert·sech²(y/L)·cos(k_x·xm)
      !!   η(y)   = −(g'/g_FS)·ξ                       (free-surface signature)
      !!   h_layer(:,:,1) = H_bed  + ξ                 (lower, k=1)
      !!   h_layer(:,:,2) = H_surf + (η − ξ)           (upper, k=2)
      !!   u_face_x_layer(:,:,2) = (g'/f)·(Δξ/L)·sech²(y/L)  (base tanh only)
      !!   u_face_x_layer(:,:,1) = 0 ;  v_face_y_layer = 0
      !!
      !! `g'` and `g_FS` are read from the SAME reduced-gravity parameters
      !! the gprime PGF consumes (`cfg%ocean%pgf%gprime_gint` /
      !! `gprime_gfs`), so the IC balance and the PGF that maintains it use
      !! one source of truth.  `f = 2Ω·sin(φ)` is the full variable
      !! Coriolis at the u-face (T-row) latitude — NOT the scalar
      !! `coriolis_f`.  The geostrophic jet uses the BASE tanh only; the
      !! perturbation is an unbalanced interface displacement the
      !! instability feeds on.
      !!
      !! Requires `nz_layers == 2`, `&ocean_pgf_nml form="gprime"`, and a
      !! spherical grid.  Fills the FULL arrays incl. ghost rows: φ/λ are
      !! continued analytically past the physical cells, which for a
      !! periodic-x channel wraps exactly because `k_x·L_x = 2π·pert_nx`
      !! (integer `pert_nx`), so the x-ghost columns match the physical
      !! ones.  Wall-y ghosts continue the smooth tanh/sech² profile.
      !!
      !! Host-compute (deterministic, RNG-free): plain host loops writing
      !! `state%multilayer` BEFORE `enter_data`, per the setup-code
      !! convention (a `do concurrent` here would round-trip unmapped
      !! arrays through the device).
      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 baroclinic-jet IC configuration conflict when
         !! present; absent behaves as today (`error stop`).
      !!
      !! MPI: the global physical-index offsets and the global domain
      !! extents come off `grid` (`i_offset_global` / `j_offset_global`,
      !! `nx_global` / `ny_global`).  All four are load-bearing: the
      !! jet-centre latitude, the perturbation wavenumber AND both
      !! index→coordinate maps must describe the WHOLE domain, or every
      !! rank reproduces the entire jet inside its own tile (with
      !! `pert_nx = 3`, three wavelengths per TILE instead of three per
      !! DOMAIN).  On a single rank the offsets are 0 and global == local,
      !! so the seed is byte-identical.

      integer :: i, j, k, nx_total, ny_total, nz_ml, ng
      integer :: idx_S, idx_T
      integer :: ioff, joff, nxg, nyg
      real(wp) :: gprime, g_fs, H_total, H_surf, H_bed
      real(wp) :: jet_L, dxi, a_pert, kx, phi0, cos_phi0, Lx
      real(wp) :: rad_earth, omega, lon_west, lat_south, dlon, dlat
      real(wp) :: phi, y_arg, tanh_y, sech2, f_u, u_jet
      real(wp) :: lam_deg, xm, xi_full, eta_local

      nz_ml = state%multilayer%nz_ml
      if (nz_ml /= 2) then
         call fail("seed_baroclinic_jet_ic: requires nz_layers == 2 "// &
                   "(two-layer reduced gravity); got "//to_string(nz_ml), ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if
      if (trim(cfg%ocean%pgf%form) /= "gprime") then
         call fail("seed_baroclinic_jet_ic: requires &ocean_pgf_nml "// &
                   "form='gprime'; got '"//trim(cfg%ocean%pgf%form)//"'", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if
      if (trim(cfg%ocean%grid%grid_config) /= "spherical") then
         call fail("seed_baroclinic_jet_ic: requires &ocean_grid_nml "// &
                   "grid_config='spherical'; got '"// &
                   trim(cfg%ocean%grid%grid_config)//"'", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if
      if (cfg%ocean%grid%omega <= 0.0_wp) then
         call fail("seed_baroclinic_jet_ic: requires &ocean_grid_nml "// &
                   "omega > 0 (planetary Coriolis f = 2*Omega*sin(phi))", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if

      nx_total = grid%nx_total
      ny_total = grid%ny_total
      ng = grid%nghost

      ! Reduced-gravity parameters — read from the gprime PGF's own config
      ! (single source of truth for g' and the free-surface gravity g_FS).
      gprime = cfg%ocean%pgf%gprime_gint
      g_fs = cfg%ocean%pgf%gprime_gfs

      ! Rest layer split: upper (surface, k=2) H_surf; lower (bed, k=1) H_bed.
      H_total = cfg%ocean%topo%max_depth
      H_surf = cfg%ocean%ic%upper_layer_rest
      H_bed = H_total - H_surf
      if (H_surf <= 0.0_wp .or. H_bed <= 0.0_wp) then
         call fail("seed_baroclinic_jet_ic: upper_layer_rest ("// &
                   to_string(H_surf)//" m) must lie in (0, max_depth="// &
                   to_string(H_total)//" m)", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if

      ioff = grid%i_offset_global
      joff = grid%j_offset_global
      nxg = grid%nx_global
      nyg = grid%ny_global

      jet_L = cfg%ocean%ic%jet_half_width
      dxi = cfg%ocean%ic%interface_amp
      a_pert = cfg%ocean%ic%pert_amp_frac*dxi

      ! Spherical geometry.  dx/dy are dlon/dlat in DEGREES for spherical.
      rad_earth = cfg%ocean%grid%rad_earth
      omega = cfg%ocean%grid%omega
      lon_west = cfg%ocean%grid%lon_west
      lat_south = cfg%ocean%grid%lat_south
      dlon = grid%dx
      dlat = grid%dy
      ! Jet-centre latitude = domain centre; φ = (lat_south + ny/2·dlat).
      phi0 = (lat_south + 0.5_wp*real(nyg, wp)*dlat)*DEG2RAD
      cos_phi0 = cos(phi0)

      ! Zonal perturbation wavenumber.  L_x = R·cosφ₀·(lon span, rad);
      ! integer pert_nx ⇒ k_x·L_x = 2π·pert_nx makes the x-continuation
      ! periodic (ghost columns match the physical wrap).
      Lx = rad_earth*cos_phi0*(real(nxg, wp)*dlon*DEG2RAD)
      kx = TWO_PI*real(cfg%ocean%ic%pert_nx, wp)/Lx

      ! ---- h_layer (cell centres), full array incl. ghosts ----
      do j = 1, ny_total
         phi = (lat_south + (real(j - ng + joff, wp) - 0.5_wp)*dlat)*DEG2RAD
         y_arg = rad_earth*(phi - phi0)/jet_L
         tanh_y = tanh(y_arg)
         sech2 = 1.0_wp/cosh(y_arg)**2
         do i = 1, nx_total
            lam_deg = lon_west + (real(i - ng + ioff, wp) - 0.5_wp)*dlon
            xm = rad_earth*cos_phi0*((lam_deg - lon_west)*DEG2RAD)
            xi_full = dxi*tanh_y + a_pert*sech2*cos(kx*xm)
            eta_local = -(gprime/g_fs)*xi_full
            state%multilayer%h_layer(i, j, 1) = H_bed + xi_full            ! lower (bed)
            state%multilayer%h_layer(i, j, 2) = H_surf + (eta_local - xi_full)  ! upper (surface)
         end do
      end do

      ! ---- balanced zonal jet on the upper-layer u-faces (row-only) ----
      ! u depends on the T-row latitude only (base tanh; NOT the pert).  The
      ! u-face shares its T-row's latitude, so f = 2Ω·sin(φ_row).  Lower
      ! layer + all v-faces stay at rest.
      state%multilayer%u_face_x_layer = 0.0_wp
      state%multilayer%v_face_y_layer = 0.0_wp
      do j = 1, ny_total
         phi = (lat_south + (real(j - ng + joff, wp) - 0.5_wp)*dlat)*DEG2RAD
         y_arg = rad_earth*(phi - phi0)/jet_L
         sech2 = 1.0_wp/cosh(y_arg)**2
         f_u = 2.0_wp*omega*sin(phi)
         u_jet = (gprime/f_u)*(dxi/jet_L)*sech2
         do i = 1, nx_total + 1
            state%multilayer%u_face_x_layer(i, j, 2) = u_jet
         end do
      end do

      ! ---- hu on x-faces: u · face-interpolated h (mirror seed_eady_ic) ----
      ! u_face/h_layer are staggered (nx+1 faces vs nx cells); end faces take
      ! the adjacent cell.  Lower layer u=0 ⇒ hu=0.  hv stays zero (v=0).
      state%multilayer%hv_face_y_layer = 0.0_wp
      do k = 1, nz_ml
         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

      ! ---- re-seed S/T uniform against the new h_layer ----
      ! The default seed set hTr = const·h_layer against the OLD uniform
      ! split; refresh so hTr/h_layer stays at the configured scalars.
      idx_S = state%multilayer%idx_salinity
      idx_T = state%multilayer%idx_temperature
      if (idx_S > 0) then
         call seed_tracer_uniform_impl( &
            state%multilayer%tracers(idx_S)%hTr, &
            state%multilayer%h_layer, cfg%initial_salinity, nz_ml)
      end if
      if (idx_T > 0) then
         call seed_tracer_uniform_impl( &
            state%multilayer%tracers(idx_T)%hTr, &
            state%multilayer%h_layer, cfg%initial_temperature, nz_ml)
      end if

      if (comm_env_rank() == 0) then
         call logger%info("seed_baroclinic_jet_ic: 2-layer reduced-gravity jet — "// &
                          "L="//to_string(jet_L)//" m, dxi="//to_string(dxi)// &
                          " m, g'="//to_string(gprime)//" m/s^2, H_surf="// &
                          to_string(H_surf)//" m, H_bed="//to_string(H_bed)// &
                          " m, n_x="//to_string(cfg%ocean%ic%pert_nx))
      end if
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine seed_baroclinic_jet_ic