subroutine multilayer_state_init(this, grid, with_ideal_age)
!! Allocate per-layer C-grid arrays at the grid size and the
!! configured layer count (caller must set `this%nz_ml` first).
!! Registers salinity + temperature with default identity strings.
!! When `with_ideal_age` is present and true, also registers an
!! ideal-age tracer at index 3 (see `rdb_ocean_ideal_age`).
class(multilayer_state_t), intent(inout) :: this
type(hgrid_t), intent(in) :: grid
logical, intent(in), optional :: with_ideal_age
integer :: nx, ny, nz_ml, ntracers
logical :: age_on
age_on = .false.
if (present(with_ideal_age)) age_on = with_ideal_age
nx = grid%nx_total
ny = grid%ny_total
nz_ml = this%nz_ml
! Cell-centred per-layer arrays
allocate (this%h_layer(nx, ny, nz_ml), source=0.0_wp)
allocate (this%h_layer0(nx, ny, nz_ml), source=0.0_wp)
allocate (this%h_av_layer(nx, ny, nz_ml), source=0.0_wp)
! East-face per-layer arrays
allocate (this%u_face_x_layer(nx + 1, ny, nz_ml), source=0.0_wp)
allocate (this%hu_face_x_layer(nx + 1, ny, nz_ml), source=0.0_wp)
allocate (this%mass_flux_x_layer(nx + 1, ny, nz_ml), source=0.0_wp)
allocate (this%u_face_x_layer0(nx + 1, ny, nz_ml), source=0.0_wp)
allocate (this%u_av_layer(nx + 1, ny, nz_ml), source=0.0_wp)
! North-face per-layer arrays
allocate (this%v_face_y_layer(nx, ny + 1, nz_ml), source=0.0_wp)
allocate (this%hv_face_y_layer(nx, ny + 1, nz_ml), source=0.0_wp)
allocate (this%mass_flux_y_layer(nx, ny + 1, nz_ml), source=0.0_wp)
allocate (this%v_face_y_layer0(nx, ny + 1, nz_ml), source=0.0_wp)
allocate (this%v_av_layer(nx, ny + 1, nz_ml), source=0.0_wp)
! Cell-centred per-layer flux divergence (pure workspace)
allocate (this%flux_h_layer(nx, ny, nz_ml), source=0.0_wp)
! Budget contributor slots, drained at every budget eval cadence.
allocate (this%mass_budget_continuity(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_surface(nx, ny, nz_ml), source=0.0_wp)
allocate (this%salt_budget_surface(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_geothermal(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_sponge(nx, ny, nz_ml), source=0.0_wp)
allocate (this%salt_budget_sponge(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_vert_adv(nx, ny, nz_ml), source=0.0_wp)
allocate (this%salt_budget_vert_adv(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_vdiff(nx, ny, nz_ml), source=0.0_wp)
allocate (this%salt_budget_vdiff(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_hdiff(nx, ny, nz_ml), source=0.0_wp)
allocate (this%salt_budget_hdiff(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_horiz_adv(nx, ny, nz_ml), source=0.0_wp)
allocate (this%salt_budget_horiz_adv(nx, ny, nz_ml), source=0.0_wp)
allocate (this%mass_budget_remap(nx, ny, nz_ml), source=0.0_wp)
allocate (this%heat_budget_remap(nx, ny, nz_ml), source=0.0_wp)
allocate (this%salt_budget_remap(nx, ny, nz_ml), source=0.0_wp)
! Cell-centred density (filled by the EOS each step)
allocate (this%rho_layer(nx, ny, nz_ml), source=0.0_wp)
! Top-of-column in-situ EOS pressure. Unconditionally allocated +
! zeroed (nx*ny*8 B against the nx*ny*nz prognostics): the OFF path
! is then an unconditional `+ 0.0` inside the kernel rather than a
! branch or an optional dummy.
allocate (this%p_top(nx, ny), source=0.0_wp)
! Cross-layer vertical velocity at interfaces (k=1 bed .. k=nz+1 surface)
allocate (this%w_interface(nx, ny, nz_ml + 1), source=0.0_wp)
! Wet-cell mask: default to all-ocean (1.0) so analytical / flat-
! bottom tests don't see a behaviour change. Realistic-bathy runs
! overwrite this in `ocean_state_seed_from_cfg` after the bathy load.
allocate (this%wet_mask(nx, ny), source=1.0_wp)
! First-live-layer indices. Seeded at the fallback `nz_ml`
! everywhere, which IS the answer on every coordinate family that
! does not vanish a layer against the top — `configure_ocean_k_top`
! overwrites them only under `z_fixed` with a rigid top.
allocate (this%k_top(nx, ny), source=nz_ml)
allocate (this%k_top_u(nx + 1, ny), source=nz_ml)
allocate (this%k_top_v(nx, ny + 1), source=nz_ml)
! Bed-side mirror, seeded at its fallback `1` — the answer on every
! family without a static bed filler; `configure_ocean_k_bot`
! overwrites it only under `z_fixed`.
allocate (this%k_bot(nx, ny), source=1)
allocate (this%k_bot_u(nx + 1, ny), source=1)
allocate (this%k_bot_v(nx, ny + 1), source=1)
! Tracer registry: salinity at index 1, temperature at index 2,
! ideal-age at index 3 (optional).
ntracers = 2
if (age_on) ntracers = 3
allocate (this%tracers(ntracers))
this%idx_salinity = 1
this%idx_temperature = 2
call this%tracers(this%idx_salinity)%init(grid, nz_ml)
call this%tracers(this%idx_temperature)%init(grid, nz_ml)
this%tracers(this%idx_salinity)%name = "salinity"
this%tracers(this%idx_salinity)%long_name = "Sea water salinity"
this%tracers(this%idx_salinity)%units = "PSU"
this%tracers(this%idx_salinity)%standard_name = "sea_water_salinity"
this%tracers(this%idx_salinity)%budget_id = TRACER_BUDGET_SALT
this%tracers(this%idx_temperature)%name = "temperature"
this%tracers(this%idx_temperature)%long_name = "Sea water potential temperature"
this%tracers(this%idx_temperature)%units = "degC"
this%tracers(this%idx_temperature)%standard_name = "sea_water_potential_temperature"
this%tracers(this%idx_temperature)%budget_id = TRACER_BUDGET_HEAT
if (age_on) then
this%idx_age = 3
call this%tracers(this%idx_age)%init(grid, nz_ml)
this%tracers(this%idx_age)%name = "age"
this%tracers(this%idx_age)%long_name = "Ideal age of sea water"
this%tracers(this%idx_age)%units = "s"
this%tracers(this%idx_age)%standard_name = "age_of_sea_water"
end if
this%is_init = .true.
end subroutine multilayer_state_init