seed_cavity_draft Subroutine

   subroutine seed_cavity_draft(state, grid, cfg, ierr)
      !! Fill `metrics%z_draft` (and its `cover_frac` companion) from
      !! `&ocean_cavity_dyn_nml`, then cross-validate the geometry against
      !! the seeded bathymetry.  Runs from `ocean_state_seed_from_cfg`
      !! IMMEDIATELY after the bathymetry and BEFORE the layer split and
      !! the wet-mask seed, which both read `b − z_draft`.
      !!
      !! Non-`pure` on purpose (the only cavity routine that is): it
      !! reports the grounded / over-land column counts through the
      !! logger and fails loud through the error ring.  The arithmetic it
      !! drives is in `rdb_ocean_cavity`, where every routine IS `pure`.
      !!
      !! UNITS.  The namelist carries metres; the setters work in GRID
      !! coordinate units (metres on Cartesian, DEGREES on
      !! spherical/curvilinear), so the box corners are converted with
      !! `topo_length_to_grid_units` and the dimensionless slope is
      !! converted the inverse way — the same trap `&ocean_topo_nml
      !! slope_scale` documents, where a metres length against a degrees
      !! position collapsed a seamount to a flat basin.
      type(ocean_state_t), intent(inout) :: state
      type(hgrid_t), intent(in) :: grid
      type(config_t), intent(in) :: cfg
      integer, intent(out) :: ierr

      integer :: draft_code, source_code, nx, ny, ng
      integer :: n_over_land, n_grounded, n_interior
      integer :: sign_code, local_ierr
      real(wp) :: per_metre, x0_g, x1_g, y0_g, y1_g, slope_g, amp
      real(wp) :: grounded_frac

      ierr = OCEAN_STATUS_OK
      nx = size(state%metrics%z_draft, 1)
      ny = size(state%metrics%z_draft, 2)
      ng = grid%nghost
      if (nx /= size(state%barotropic%b, 1) .or. ny /= size(state%barotropic%b, 2)) then
         call fail("&ocean_cavity_dyn_nml: z_draft is at its placeholder size — "// &
                   "metrics%use_cavity must be latched BEFORE metrics%init (it is "// &
                   "latched in ocean_state_init_from_config)", &
                   ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if

      draft_code = parse_cavity_draft_config(cfg%ocean%cavity_dyn%draft_config)
      source_code = parse_cavity_draft_source(cfg%ocean%cavity_dyn%draft_source)
      ! `validate_config` already refuses every spelling outside the v1
      ! envelope with a full explanation; this is the same gate one level
      ! down, for direct (test / API) callers that bypass it.
      if (draft_code /= CAVITY_DRAFT_NONE .and. draft_code /= CAVITY_DRAFT_FLAT &
          .and. draft_code /= CAVITY_DRAFT_LINEAR .and. draft_code /= CAVITY_DRAFT_FILE) then
         call fail("&ocean_cavity_dyn_nml draft_config='"// &
                   trim(cfg%ocean%cavity_dyn%draft_config)//"' is not available "// &
                   "(none|flat|linear|file)", ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if
      if (source_code /= CAVITY_SOURCE_DRAFT .and. source_code /= CAVITY_SOURCE_THICKNESS) then
         call fail("&ocean_cavity_dyn_nml draft_source='"// &
                   trim(cfg%ocean%cavity_dyn%draft_source)//"' is not available "// &
                   "(draft|thickness; 'in_situ' isostasy is deferred)", &
                   ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if

      ! GRID units per metre (1 on Cartesian; degrees-per-metre otherwise).
      ! The "no limit" sentinels pass through UNCONVERTED so they stay
      ! sentinels on every grid.
      per_metre = topo_length_to_grid_units(1.0_wp, cfg%ocean%grid%grid_config, &
                                            cfg%ocean%grid%rad_earth)
      x0_g = cavity_bound_to_grid(cfg%ocean%cavity_dyn%draft_x0, per_metre)
      x1_g = cavity_bound_to_grid(cfg%ocean%cavity_dyn%draft_x1, per_metre)
      y0_g = cavity_bound_to_grid(cfg%ocean%cavity_dyn%draft_y0, per_metre)
      y1_g = cavity_bound_to_grid(cfg%ocean%cavity_dyn%draft_y1, per_metre)
      ! `draft_slope` is d(draft [m]) / d(x [m]); the setter wants
      ! d(draft [m]) / d(x [grid units]) = slope / (grid units per metre).
      slope_g = cfg%ocean%cavity_dyn%draft_slope/per_metre

      ! `draft_source = "thickness"`: the formula amplitude is an ice
      ! THICKNESS, converted by the Boussinesq-isostatic (flotation)
      ! relation `z_draft = rho_ice*h_ice/rho_0`.  Scaling the amplitude
      ! (and the slope with it) is exact because both setters are LINEAR
      ! in the amplitude — and it keeps one draft field, so everything
      ! downstream stays source-agnostic.
      amp = cfg%ocean%cavity_dyn%draft_depth
      if (source_code == CAVITY_SOURCE_THICKNESS) then
         amp = amp*cfg%ocean%cavity_dyn%rho_ice/state%eos%rho0
         slope_g = slope_g*cfg%ocean%cavity_dyn%rho_ice/state%eos%rho0
      end if

      select case (draft_code)
      case (CAVITY_DRAFT_FLAT)
         call set_draft_flat(state%metrics%z_draft, grid, amp, x0_g, x1_g, y0_g, y1_g)
      case (CAVITY_DRAFT_LINEAR)
         call set_draft_linear(state%metrics%z_draft, grid, amp, slope_g, &
                               x0_g, x1_g, y0_g, y1_g)
      case (CAVITY_DRAFT_FILE)
#ifndef RDB_NO_NETCDF
         ! Static 2-D NetCDF draft, through the PR-14 reader.  SINGLE
         ! RANK: the loader itself applies the global offset correctly,
         ! but the grounding statistics a few lines below are single-rank
         ! reductions and the whole cavity is fenced that way, so the
         ! restriction is asserted here rather than left implicit.
         if (grid%nx_phys /= grid%nx_global .or. grid%ny_phys /= grid%ny_global) then
            call fail("&ocean_cavity_dyn_nml draft_config='file' is single-rank "// &
                      "only (the cavity's grounding statistics are single-rank "// &
                      "reductions).  Run on one rank or use an analytic draft.", &
                      ierr, OCEAN_STATUS_ERR_IC_SEED)
            return
         end if
         sign_code = parse_cavity_draft_sign(cfg%ocean%cavity_dyn%draft_sign)
         if (sign_code == CAVITY_SIGN_INVALID) then
            call fail("&ocean_cavity_dyn_nml draft_sign='"// &
                      trim(adjustl(cfg%ocean%cavity_dyn%draft_sign))// &
                      "' is not recognised (depth|positive_down|elevation|"// &
                      "positive_up)", ierr, OCEAN_STATUS_ERR_IC_SEED)
            return
         end if
         ! Interior first (the reader writes the physical window only)...
         state%metrics%z_draft = 0.0_wp
         call ocean_data_input_load_static_2d( &
            trim(cfg%ocean%cavity_dyn%draft_file), &
            trim(cfg%ocean%cavity_dyn%draft_var), grid, nx, ny, &
            ng + 1, ng + 1, state%metrics%z_draft, ierr=local_ierr)
         if (local_ierr /= OCEAN_STATUS_OK) then
            ierr = local_ierr
            return
         end if
         ! ...sign-normalise onto DEPTH positive down...
         call cavity_draft_apply_sign(state%metrics%z_draft, nx, ny, sign_code)
         ! ...then fill the ghost band by constant extrapolation, the
         ! SAME routine and the same order the file bathymetry uses
         ! (`load_bathymetry_into_array` -> `fill_bathymetry_ghosts_array`).
         ! The periodic/fold re-wrap and the halo exchange that
         ! `metrics%z_draft` gets in `rdb_ocean_engine` run later and are
         ! shared with the formula path, so a file draft and a formula
         ! draft see an identical boundary treatment.
         call bathymetry_fill_ghosts_array(state%metrics%z_draft, grid)
#else
         call fail("&ocean_cavity_dyn_nml draft_config='file' requires "// &
                   "RDB_ENABLE_NETCDF=ON at build time (the static-2-D reader "// &
                   "lives in the NetCDF-gated rdb_ocean_data_input).", &
                   ierr, OCEAN_STATUS_ERR_IO)
         return
#endif
      case default  ! CAVITY_DRAFT_NONE
         state%metrics%z_draft = 0.0_wp
      end select

      ! Guard the field itself before anything derives geometry from it.
      ! Written as `.not. (z >= 0)` inside the helper so a NaN FAILS
      ! rather than sliding through a `z < 0` test that is false for NaN.
      if (.not. cavity_draft_is_finite_nonneg(state%metrics%z_draft, nx, ny)) then
         call fail("&ocean_cavity_dyn_nml: z_draft must be finite and >= 0 "// &
                   "everywhere (it is a DEPTH below z = 0, positive down)", &
                   ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if

      ! NO ICE OVER LAND: zero the draft wherever the bathymetry already
      ! says land, so land columns keep the datum they always had
      ! (`bt_H_ref = b`) and the counted-once invariant stays exact on
      ! every column.
      call cavity_apply_land_exclusion(state%metrics%z_draft, state%barotropic%b, &
                                       nx, ny, n_over_land)

      call cavity_count_grounded(state%barotropic%b, state%metrics%z_draft, &
                                 cfg%ocean%cavity_dyn%h_min_cavity, ng, &
                                 grid%nx_phys, grid%ny_phys, nx, ny, &
                                 n_grounded, n_interior)
      grounded_frac = real(n_grounded, wp)/real(max(n_interior, 1), wp)
      if (grounded_frac > cfg%ocean%cavity_dyn%grounded_max_frac) then
         call fail("&ocean_cavity_dyn_nml: the prescribed draft grounds "// &
                   to_string(n_grounded)//" of "//to_string(n_interior)// &
                   " interior columns ("//to_string(grounded_frac)//"), above "// &
                   "grounded_max_frac = "// &
                   to_string(cfg%ocean%cavity_dyn%grounded_max_frac)// &
                   ".  Either the draft is too deep for this bathymetry or the "// &
                   "shelf box is in the wrong place (check the UNITS of "// &
                   "draft_x0/x1 — they are metres, converted to grid units).", &
                   ierr, OCEAN_STATUS_ERR_IC_SEED)
         return
      end if

      call cavity_fill_cover_frac(state%metrics%cover_frac, state%metrics%z_draft, nx, ny)

      if (comm_env_rank() == 0) then
         call logger%info("Ice-shelf cavity: ON  draft_config='"// &
                          trim(cfg%ocean%cavity_dyn%draft_config)//"' source='"// &
                          trim(cfg%ocean%cavity_dyn%draft_source)//"' max draft = "// &
                          to_string(maxval(state%metrics%z_draft))//" m, "// &
                          "h_min_cavity = "// &
                          to_string(cfg%ocean%cavity_dyn%h_min_cavity)//" m")
         call logger%info("                  datum bt_H_ref = b - z_draft afloat, "// &
                          "0 where grounded; "// &
                          to_string(n_grounded)//" of "//to_string(n_interior)// &
                          " interior columns grounded (-> LAND via the wet mask)")
         if (n_over_land > 0) then
            call logger%warning("&ocean_cavity_dyn_nml: draft zeroed on "// &
                                to_string(n_over_land)//" column(s) whose bed is "// &
                                "already land (b < LAND_DEPTH_THRESHOLD) — no ice "// &
                                "over land.  Check the shelf box if that is a surprise.")
         end if
      end if
   end subroutine seed_cavity_draft