seed_ts_from_zfile Subroutine

   subroutine seed_ts_from_zfile(ms, grid, cfg, ierr, z_draft)
      !! Read T/S from the configured z-level NetCDF and overwrite the
      !! tracer slots with the depth-interpolated profiles.  Wet columns
      !! interpolate; dry columns get the namelist land-fill constants.
      !! Must be called AFTER bathymetry, the uniform `h_layer` seed, and
      !! `wet_mask` are in place.  Takes the multilayer slot directly (not
      !! the full ocean state) to avoid a circular dependency.
      type(multilayer_state_t), intent(inout) :: ms
      type(hgrid_t), intent(in) :: grid
      type(ocean_zinit_config_t), intent(in) :: cfg
      integer, intent(out), optional :: ierr
         !! Non-zero on a malformed/mismatched z-level IC file when
         !! present; absent behaves as today (`error stop`).
      real(wp), intent(in), optional :: z_draft(:, :)
         !! Ice-base depth (m, positive-down, `>= 0`), FULL ghosted
         !! `(nx_total, ny_total)` shape — `metrics%z_draft`.  The
         !! geopotential depth of every layer centre is measured from
         !! `z = 0`, so the column top sits at `z_draft` under a shelf.
         !! Absent ⇒ `z_top = 0` (open ocean), bit-identical to the
         !! pre-cavity behaviour.

      integer :: ncid, ng, nx, ny, nz_ml, nz_src, local_ierr
      integer :: idx_S, idx_T, i, j, k
      integer :: t_varid, s_varid, z_varid
      real(wp), allocatable :: t_src(:, :, :), s_src(:, :, :), z_src(:)
      real(wp) :: z_ctr(ms%nz_ml)
      real(wp) :: t_col(ms%nz_ml), s_col(ms%nz_ml)
      real(wp) :: z_top
      logical :: needs_transpose, 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
      ng = grid%nghost
      nx = grid%nx_phys
      ny = grid%ny_phys
      nz_ml = ms%nz_ml
      idx_S = ms%idx_salinity
      idx_T = ms%idx_temperature

      call logger%info("Loading z-level T/S IC from: "//trim(cfg%file))

      call nc_open_read(trim(cfg%file), ncid, ierr=local_ierr)
      if (.not. zinit_io_ok(local_ierr, ierr)) return

      ! Locate the three variables: namelist override then name-fallbacks.
      ! ierr threaded down ONLY when THIS routine's own ierr is present:
      ! otherwise find_var/read_dims must keep reaching their own
      ! `error stop` (specific text) rather than the generic wrapper
      ! message below (P0.1 review F2).
      if (present(ierr)) then
         call find_var(ncid, cfg%t_var, &
                       [character(len=16) :: "temp", "T", "temperature"], "temperature", t_varid, &
                       ierr=local_ierr)
         if (local_ierr /= 0) then
            ierr = local_ierr
            return
         end if
      else
         call find_var(ncid, cfg%t_var, &
                       [character(len=16) :: "temp", "T", "temperature"], "temperature", t_varid)
      end if
      if (present(ierr)) then
         call find_var(ncid, cfg%s_var, &
                       [character(len=16) :: "salt", "S", "salinity"], "salinity", s_varid, &
                       ierr=local_ierr)
         if (local_ierr /= 0) then
            ierr = local_ierr
            return
         end if
      else
         call find_var(ncid, cfg%s_var, &
                       [character(len=16) :: "salt", "S", "salinity"], "salinity", s_varid)
      end if
      if (present(ierr)) then
         call find_var(ncid, cfg%z_var, &
                       [character(len=16) :: "z_src", "z", "depth", "lev"], "source axis", z_varid, &
                       ierr=local_ierr)
         if (local_ierr /= 0) then
            ierr = local_ierr
            return
         end if
      else
         call find_var(ncid, cfg%z_var, &
                       [character(len=16) :: "z_src", "z", "depth", "lev"], "source axis", z_varid)
      end if

      ! Validate dims + determine the storage-order permutation against
      ! the model grid; error-stops on mismatch.
      if (present(ierr)) then
         call read_dims(ncid, t_varid, grid, nz_src, needs_transpose, ierr=local_ierr)
         if (local_ierr /= 0) then
            ierr = local_ierr
            return
         end if
      else
         call read_dims(ncid, t_varid, grid, nz_src, needs_transpose)
      end if

      ! Read the source axis + assert monotonic increase (positive-down).
      allocate (z_src(nz_src))
      call nc_get_var_1d(ncid, z_varid, z_src, ierr=local_ierr)
      if (.not. zinit_io_ok(local_ierr, ierr, ncid)) return
      do k = 2, nz_src
         if (z_src(k) <= z_src(k - 1)) then
            call nc_close(ncid)
            call fail("ocean_zinit: z_src not monotonically increasing at level "// &
                      to_string(k)//" ("//to_string(z_src(k))//" <= "// &
                      to_string(z_src(k - 1))//")", ierr, OCEAN_STATUS_ERR_IO)
            return
         end if
      end do

      ! Read T/S into model-grid (x, y, z) interior arrays, undoing the
      ! C/Fortran dimension reversal when the file is C-ordered (z, y, x).
      allocate (t_src(nx, ny, nz_src), s_src(nx, ny, nz_src))
      call read_field_xyz(ncid, t_varid, t_src, nx, ny, nz_src, needs_transpose, &
                          grid%i_offset_global, grid%j_offset_global, ierr=local_ierr)
      if (.not. zinit_io_ok(local_ierr, ierr, ncid)) return
      call read_field_xyz(ncid, s_varid, s_src, nx, ny, nz_src, needs_transpose, &
                          grid%i_offset_global, grid%j_offset_global, ierr=local_ierr)
      if (.not. zinit_io_ok(local_ierr, ierr, ncid)) return
      call nc_close(ncid)

      ! Host-side per-column interpolation.  Plain do loops: this runs
      ! before enter_data, so the arrays are host-resident.
      do j = 1, ny
         do i = 1, nx
            ! Layer-centre GEOPOTENTIAL depths (positive-down from z = 0)
            ! from h_layer, offset by the depth of the column top.
            ! Bottom-up: k=1 bed (deepest), k=nz_ml surface (shallowest).
            z_top = 0.0_wp
            if (have_draft) z_top = z_draft(ng + i, ng + j)
            call build_z_ctr(ms%h_layer(ng + i, ng + j, :), nz_ml, z_top, z_ctr)

            if (ms%wet_mask(ng + i, ng + j) <= 0.0_wp) then
               ! Dry column — fill with the namelist land-fill constants.
               t_col = cfg%land_fill_t
               s_col = cfg%land_fill_s
            else
               ! Interp in DEPTH SPACE so the source's ascending
               ! positive-down order and the model's bottom-up order
               ! never need reconciling.
               call interp_column_linear_z(z_src, t_src(i, j, :), nz_src, z_ctr, nz_ml, t_col)
               call interp_column_linear_z(z_src, s_src(i, j, :), nz_src, z_ctr, nz_ml, s_col)
            end if

            ! Roundabout tracer convention: hTr = value * h_layer.
            if (idx_T > 0) then
               do k = 1, nz_ml
                  ms%tracers(idx_T)%hTr(ng + i, ng + j, k) = &
                     t_col(k)*ms%h_layer(ng + i, ng + j, k)
               end do
            end if
            if (idx_S > 0) then
               do k = 1, nz_ml
                  ms%tracers(idx_S)%hTr(ng + i, ng + j, k) = &
                     s_col(k)*ms%h_layer(ng + i, ng + j, k)
               end do
            end if
         end do
      end do

      deallocate (t_src, s_src, z_src)

      call logger%info("ocean_zinit: seeded T/S from "//to_string(nz_src)// &
                       " source z-levels onto "//to_string(nz_ml)//" model layers.")
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine seed_ts_from_zfile