seed_ts_linear_z Subroutine

   subroutine seed_ts_linear_z(ms, cfg, ierr, z_draft)
      !! Seed T and S from the ANALYTIC affine geopotential profiles
      !! `T(z) = lin_t_ref + lin_dt_dz*z`, `S(z) = lin_s_ref +
      !! lin_ds_dz*z` (`z` positive UP, zero at the `z = 0` datum),
      !! sampled at each layer centre's TRUE geopotential depth and
      !! written as `hTr = value * h_layer`.
      !!
      !! Opens no file, so — unlike the NetCDF path — it fills the FULL
      !! array INCLUDING the ghost rows and the dry columns: the formula
      !! is defined everywhere, and a ghost column left on some other
      !! profile is the wall-adjacent spurious-density-jump gotcha
      !! (CLAUDE.md, "formula bathymetry setters must fill ghost rows").
      !! There is no `land_fill_t/_s` on this path for the same reason.
      !!
      !! NOT `pure`: it logs one line reporting what it seeded, matching
      !! `seed_ts_from_zfile`.  The arithmetic itself is the `pure
      !! elemental` `linear_in_z`.
      type(multilayer_state_t), intent(inout) :: ms
      type(ocean_zinit_config_t), intent(in) :: cfg
      integer, intent(out), optional :: ierr
         !! Non-zero on a malformed `z_draft` when present; absent
         !! behaves like the rest of the seed path (`error stop`).
      real(wp), intent(in), optional :: z_draft(:, :)
         !! Ice-base depth (m, positive-down), FULL ghosted
         !! `(nx_total, ny_total)` shape.  Absent ⇒ open ocean.

      integer :: nxt, nyt, nz_ml, idx_S, idx_T, i, j, k
      real(wp) :: z_ctr(ms%nz_ml)
      real(wp) :: z_top
      logical :: have_draft

      have_draft = present(z_draft)
      if (have_draft) then
         if (.not. draft_shape_ok(z_draft, ms)) then
            call fail("ocean_zinit: z_draft must be the FULL ghosted "// &
                      "(nx_total, ny_total) array, matching h_layer", &
                      ierr, OCEAN_STATUS_ERR_IO)
            return
         end if
      end if

      nxt = size(ms%h_layer, 1)
      nyt = size(ms%h_layer, 2)
      nz_ml = ms%nz_ml
      idx_S = ms%idx_salinity
      idx_T = ms%idx_temperature

      ! Host-side, before enter_data: plain do loops, NOT do concurrent.
      do j = 1, nyt
         do i = 1, nxt
            z_top = 0.0_wp
            if (have_draft) z_top = z_draft(i, j)
            call build_z_ctr(ms%h_layer(i, j, :), nz_ml, z_top, z_ctr)
            if (idx_T > 0) then
               do k = 1, nz_ml
                  ms%tracers(idx_T)%hTr(i, j, k) = &
                     linear_in_z(cfg%lin_t_ref, cfg%lin_dt_dz, z_ctr(k))* &
                     ms%h_layer(i, j, k)
               end do
            end if
            if (idx_S > 0) then
               do k = 1, nz_ml
                  ms%tracers(idx_S)%hTr(i, j, k) = &
                     linear_in_z(cfg%lin_s_ref, cfg%lin_ds_dz, z_ctr(k))* &
                     ms%h_layer(i, j, k)
               end do
            end if
         end do
      end do

      if (have_draft) then
         call logger%info("ocean_zinit: seeded analytic linear-in-z T/S onto "// &
                          to_string(nz_ml)//" model layers, measuring the "// &
                          "geopotential depth from the ice base (cavity draft).")
      else
         call logger%info("ocean_zinit: seeded analytic linear-in-z T/S onto "// &
                          to_string(nz_ml)//" model layers.")
      end if
      call logger%info("ocean_zinit:   T = "//to_string(cfg%lin_t_ref)//" + "// &
                       to_string(cfg%lin_dt_dz)//"*z,  S = "// &
                       to_string(cfg%lin_s_ref)//" + "// &
                       to_string(cfg%lin_ds_dz)//"*z   (z positive UP from z = 0).")
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine seed_ts_linear_z