Host-side setup: decomposition -> grid -> god state -> IC seed
-> restart (optional) -> the 21 configure_ocean_*-family
stages -> ghost wraps -> halo init -> land mask -> wave
drag/porous -> sea-ice IC. Does NOT map the state onto the
device — call engine_enter_data next. Preserves
driver_run_ocean’s exact stage order; see the per-stage
comments below (carried over from the driver almost verbatim).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ocean_engine_t), | intent(inout) | :: | engine | |||
| type(config_t), | intent(inout) | :: | cfg | |||
| integer, | intent(out), | optional | :: | ierr |
Non-zero ( |
|
| integer, | intent(in), | optional | :: | compute_rank |
This rank’s 0-based compute-rank index. Default 0 (single-rank: the API/bench callers). |
|
| integer, | intent(in), | optional | :: | compute_size |
Total compute-rank count. Default 1. |
|
| integer, | intent(in), | optional | :: | mpi_rank |
Rank used for the restart filename convention
( |
|
| character(len=*), | intent(in), | optional | :: | restart_file |
Warm-restart source (file, or a run-directory — see
|
|
| real(kind=wp), | intent(out), | optional | :: | t_restart |
Restored simulation time (0 on a cold start). |
|
| integer, | intent(out), | optional | :: | step_restart |
Restored outer-step count (0 on a cold start). |
|
| logical, | intent(in), | optional | :: | validate_only |
|
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| character(len=:), | private, | allocatable | :: | bt_excl_reason | |||
| logical, | private | :: | bt_excluded | ||||
| integer, | private | :: | bt_halo_req | ||||
| integer, | private | :: | bt_halo_res | ||||
| integer, | private | :: | csize | ||||
| logical, | private | :: | did_restart | ||||
| logical, | private | :: | dist_fold | ||||
| integer, | private | :: | mrank | ||||
| logical, | private | :: | no_output | ||||
| real(kind=wp), | private, | allocatable | :: | q_heat_resume(:,:) | |||
| real(kind=wp), | private, | allocatable | :: | q_salt_resume(:,:) | |||
| integer, | private | :: | rank | ||||
| character(len=512), | private | :: | restart_filename | ||||
| logical, | private | :: | resume_q | ||||
| integer, | private | :: | step_restart_local | ||||
| real(kind=wp), | private | :: | t_restart_local |
subroutine engine_setup(engine, cfg, ierr, compute_rank, compute_size, mpi_rank, & restart_file, t_restart, step_restart, validate_only) !! Host-side setup: decomposition -> grid -> god state -> IC seed !! -> restart (optional) -> the 21 `configure_ocean_*`-family !! stages -> ghost wraps -> halo init -> land mask -> wave !! drag/porous -> sea-ice IC. Does NOT map the state onto the !! device — call `engine_enter_data` next. Preserves !! `driver_run_ocean`'s exact stage order; see the per-stage !! comments below (carried over from the driver almost verbatim). type(ocean_engine_t), intent(inout) :: engine type(config_t), intent(inout) :: cfg integer, intent(out), optional :: ierr !! Non-zero (`rdb_ocean_status` code) on a bad config/setup/IC !! failure when present; absent behaves as today — the !! offending `configure_ocean_*` stage (or this routine's own !! pre-flight checks) `error stop`s. integer, intent(in), optional :: compute_rank !! This rank's 0-based compute-rank index. Default 0 !! (single-rank: the API/bench callers). integer, intent(in), optional :: compute_size !! Total compute-rank count. Default 1. integer, intent(in), optional :: mpi_rank !! Rank used for the restart filename convention !! (`output_rank_filename`). Default = `compute_rank` (the !! ocean path has no separate I/O-server rank). character(len=*), intent(in), optional :: restart_file !! Warm-restart source (file, or a run-directory — see !! `driver_run_ocean`'s `.nc`-suffix convention). Absent/blank !! => cold start (matches the API/bench callers today). real(wp), intent(out), optional :: t_restart !! Restored simulation time (0 on a cold start). integer, intent(out), optional :: step_restart !! Restored outer-step count (0 on a cold start). logical, intent(in), optional :: validate_only !! `rdb --validate-only`: run every configure stage (and so every !! configure-time refusal) but create no output — the diag !! selection is still parsed and checked, the per-rank NetCDF !! stream is not opened and its directory is not created. !! Default `.false.`. integer :: rank, csize, mrank logical :: no_output character(len=512) :: restart_filename real(wp) :: t_restart_local integer :: step_restart_local logical :: did_restart logical :: resume_q real(wp), allocatable :: q_heat_resume(:, :), q_salt_resume(:, :) logical :: bt_excluded character(len=:), allocatable :: bt_excl_reason integer :: bt_halo_req, bt_halo_res logical :: dist_fold rank = 0 if (present(compute_rank)) rank = compute_rank csize = 1 if (present(compute_size)) csize = compute_size mrank = rank if (present(mpi_rank)) mrank = mpi_rank no_output = .false. if (present(validate_only)) no_output = validate_only if (present(ierr)) ierr = OCEAN_STATUS_OK t_restart_local = 0.0_wp step_restart_local = 0 ! ---- pre-flight (was comm_env_abort in the driver; a library ! cannot abort its host process, so these are now fail()-based) ---- if (cfg%dt_fixed <= 0.0_wp) then call fail("sim_type='ocean' requires dt_fixed > 0 (no adaptive CFL helper yet).", & ierr, OCEAN_STATUS_ERR_SETUP) return end if ! Process grid left unset (the `&mpi_nml` default px = py = 1) on more ! than one rank: choose it here, BEFORE the px*py check below (which ! used to reject it, leaving `decomp_init_from_config`'s own ! auto-factor unreachable on the ocean path). A folded grid is split ! north-south by default (px = 1: the fold stays local, no fold ! exchange; plan `tripolar_fold_px_gt_1` decision 9 — an east-west ! split is supported but explicit, `&mpi_nml px > 1`, until its cost ! is measured at scale); any other grid gets the perimeter-minimising ! factorisation. if (csize > 1 .and. cfg%px == 1 .and. cfg%py == 1) then if (ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_TRIPOLAR_FOLD) then cfg%px = 1 cfg%py = csize else call decomp_auto_factor(csize, cfg%nx, cfg%ny, cfg%px, cfg%py) end if if (rank == 0) call logger%info("Process grid auto: px x py = "// & to_string(cfg%px)//" x "//to_string(cfg%py)// & " (&mpi_nml px/py unset)") end if if (cfg%px*cfg%py /= csize) then call fail("Process grid px*py = "//to_string(cfg%px*cfg%py)// & " does not match number of compute ranks = "//to_string(csize)// & " (ocean path)", ierr, OCEAN_STATUS_ERR_SETUP) return end if if (csize > 1) then ! Tripolar north fold, east-west split (px > 1): the fold point ! (i, nj+d) mirrors (ni+1-i, nj+1-d), on another rank of the north ! row, so the fold runs through the owner-routed exchange of ! `rdb_ocean_fold_exchange` (plan `tripolar_fold_px_gt_1`). Tile ! WIDTH rule, symmetric with the height rule below (plan decision ! 11): refuse tiles narrower than nghost+1 columns. The exchange ! itself routes any width; the rule keeps every tile wider than its ! own ghost band, as the halo's seam stencils assume. The narrowest ! tile is nx/px (the remainder goes to the first columns), so the ! check is rank-invariant and every rank fails together. if (ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_TRIPOLAR_FOLD .and. & cfg%px > 1 .and. cfg%nx/cfg%px < cfg%nghost + 1) then call fail("Tripolar north fold with px = "//to_string(cfg%px)// & ": tiles of nx/px = "//to_string(cfg%nx/cfg%px)// & " columns are narrower than nghost+1 = "// & to_string(cfg%nghost + 1)//"; reduce px.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if ! The fold reads the nghost rows just below the fold line from the ! north tile itself, so that tile must be at least nghost+1 rows ! tall (the v fold row plus the nghost rows it mirrors). The ! smallest tile is ny/py (the remainder goes to the first rows), so ! the check is rank-invariant and every rank fails together. if (ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_TRIPOLAR_FOLD .and. & cfg%py > 1 .and. cfg%ny/cfg%py < cfg%nghost + 1) then call fail("Tripolar north fold with py = "//to_string(cfg%py)// & ": tiles of ny/py = "//to_string(cfg%ny/cfg%py)// & " rows are too short for the fold's mirror (need >= nghost+1 = "// & to_string(cfg%nghost + 1)//"); reduce py.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if ! E1: windowed tracer advect drain halo not yet wired for multi-rank. if (cfg%ocean%vmix%dt_tracer_advect_ratio > 1) then call fail("dt_tracer_advect_ratio > 1 multi-rank drain halo "// & "is deferred (E1); run single-rank or set ratio = 1.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if ! The single-rank features, keyed on the ACTUAL rank count. ! `validate_config` refuses them too, but on `px*py > 1` -- which ! an unset process grid (px = py = 1, auto-factored above) does ! not trip -- so this is the gate that always holds. Every rank ! evaluates the same condition and fails together. if (cfg%ocean%cavity_dyn%enable) then call fail("&ocean_cavity_dyn_nml enable = .true. is single-rank ("// & to_string(csize)//" ranks requested): the grounding statistics "// & "are global reductions the configure does not take. Run on 1 rank.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if if (cfg%ocean%wetdry%enable) then call fail("&ocean_wetdry_nml enable = .true. is single-rank ("// & to_string(csize)//" ranks requested): the wet-mask / outflow-"// & "limiter halo exchange is not implemented. Run on 1 rank.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if ! Sea ice runs on the ocean's decomposition: the category state, ! the EVP ice velocity, the transport's CAS state and the blended ! surface stress are halo-exchanged (`engine_step_ice`). The ! category state (part_size/m_ice/m_snow/enth_ice/sal_ice/ ! enth_snow/mca_ice/mca_snow) IS now folded across the tripolar ! north seam (`ocean_fold_wrap_centre_flat`, fold-seam fix) -- ! on `px = 1` this is exercised and tested ! (`test_ocean_ice_fold`). The refusal below stays for `px > 1` ! only: the ice fold rides the SAME px-aware dispatcher as the ! ocean prognostics (`ocean_fold_is_distributed`), but that path ! has (a) no MPI test coverage for ice (no decomp-bitid case), ! and (b) `ocean_fold_exchange_reserve` below is sized for the ! ml_state group only (`nz_layers*(3+ntr)`), not for the ice ! category group (up to `ncat*nk_ice`/`ncat+1` levels) -- an ! under-reservation `ocean_fold_begin` only warns and grows ! mid-run, not a correctness bug, but it does defeat the ! CUDA-aware-MPI IPC-handle-reuse this reserve exists for. EVP ! (`&ocean_ice_nml dynamics`) additionally has no fold of its ! own at any rank count and stays out of scope here. Lifting ! this for `px > 1` needs both the reserve sizing and a ! `test_ocean_decomp_bitid_mpi` ice+fold case; tracked as a ! follow-up, not done in this PR. if (cfg%ocean%ice%enable .and. & ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_TRIPOLAR_FOLD) then call fail("&ocean_ice_nml enable = .true. with north = 'tripolar_fold' on "// & to_string(csize)//" ranks: the ice category fold is untested under "// & "MPI (no decomp-bitid coverage) and its fold-exchange buffers are "// & "not reserved for the ice group. Run on 1 rank (px=1 is fully "// & "folded and tested).", & ierr, OCEAN_STATUS_ERR_SETUP) return end if ! In-memory geometry injection (the API's staged bathymetry / ! supergrid arrays) hands over WHOLE-grid arrays. if (allocated(engine%staged_bathymetry) .or. allocated(engine%staged_metrics_x)) then call fail("in-memory geometry injection (staged bathymetry / supergrid "// & "arrays) is single-rank ("//to_string(csize)//" ranks requested): "// & "the arrays describe the whole grid. Use the file readers "// & "(topo_config='file', grid_config='supergrid') or 1 rank.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if ! Chapman radiation keeps ONE edge-uniform eta target per edge, ! computed from the tile's own stretch of the edge (and a corner ! depth): split along a Chapman edge, every rank radiates toward a ! different target. if (any([ocean_bc_type_from_string(cfg%ocean%bc%west), & ocean_bc_type_from_string(cfg%ocean%bc%east), & ocean_bc_type_from_string(cfg%ocean%bc%south), & ocean_bc_type_from_string(cfg%ocean%bc%north)] == OBC_CHAPMAN)) then call fail("&ocean_bc_nml 'chapman' edges are single-rank ("// & to_string(csize)//" ranks requested): the edge-mean eta target "// & "is a per-rank partial mean. Use 'open' (Flather) or run on 1 rank.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if end if call decomp_init_from_config(engine%decomp, cfg, csize, rank) if (rank == 0 .and. csize > 1) call decomp_log_summary(engine%decomp, csize) ! Resolve the BT wide-halo march-in sentinel now that compute_size and ! the feature flags are both known (see rdb_config::resolve_bt_halo). bt_halo_req = cfg%ocean%bt%bt_halo call bt_halo_auto_exclusion(cfg, bt_excluded, bt_excl_reason) bt_halo_res = resolve_bt_halo(bt_halo_req, csize, bt_excluded) cfg%ocean%bt%bt_halo = bt_halo_res ! An EXPLICIT width must fit inside the smallest tile too (the wide ! clone's ghost ring is filled from the tile's own interior); refuse it ! here with a status instead of the enable step's `error stop`. if (bt_halo_res > 0 .and. & cfg%nghost + bt_halo_res - mod(bt_halo_res, 2) > min(cfg%nx/cfg%px, cfg%ny/cfg%py)) then call fail("&ocean_bt_nml bt_halo = "//to_string(bt_halo_res)//" does not fit: "// & "nghost + bt_halo = "//to_string(cfg%nghost + bt_halo_res - mod(bt_halo_res, 2))// & " exceeds the smallest tile ("//to_string(min(cfg%nx/cfg%px, cfg%ny/cfg%py))// & " cells). Reduce bt_halo (or leave it on auto) or use fewer ranks.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if if (rank == 0 .and. bt_halo_req == BT_HALO_AUTO_SENTINEL .and. csize > 1) then if (bt_excluded) then call logger%info("BT march-in: bt_halo auto -> 0 ("// & trim(bt_excl_reason)//" active)") else call logger%info("BT march-in: bt_halo auto -> 0 (opt-in: set "// & "&ocean_bt_nml bt_halo = 8 explicitly; it is not "// & "bit-reproducible across decompositions)") end if end if if (rank == 0 .and. bt_halo_res > 0 .and. csize > 1) then call logger%warning("BT march-in ON (bt_halo = "//to_string(bt_halo_res)// & "): answers are NOT bit-identical to the serial run over "// & "variable bathymetry or with open boundaries.") end if call engine%grid%init(engine%decomp%nx_local, engine%decomp%ny_local, & cfg%nghost, cfg%dx, cfg%dy) engine%grid%i_offset_global = engine%decomp%i_start - 1 engine%grid%j_offset_global = engine%decomp%j_start - 1 engine%grid%nx_global = engine%decomp%nx_global engine%grid%ny_global = engine%decomp%ny_global call engine%state%init_from_config(cfg, engine%grid) engine%state%multilayer%mass_out_efp_on = cfg%ocean%diag%reproducing_sums ! Wire vcoord parameters BEFORE the seed runs — the seed calls ! `vcoord%build_zref_full(b)` at its tail, which reads ! `zstar_h_surf_target` etc. if (engine%state%use_multilayer) then engine%state%vcoord%coord_type = parse_ocean_vcoord_type(cfg%vcoord_type) engine%state%vcoord%remap_method = parse_remap_method(cfg%remap_method) engine%state%vcoord%zstar_h_surf_target = cfg%zstar_h_surf_target engine%state%vcoord%zstar_h_min = cfg%zstar_h_min ! `z_fixed` nominal layering (uniform or a stretched profile): the ! cavity z_fixed seed below lays `h_layer` from the target builder. call configure_ocean_z_fixed_profile(cfg, engine%state, rank, log_it=.true.) end if ! Analytical IC from cfg scalars — or, when P2.5 geometry injection ! staged a bathymetry array (rdb_ocean_stage_bathymetry), that ! array overrides cfg%ocean%topo%topo_config entirely. ! ! The seed makes the static geometry seam-consistent (periodic wrap + ! north fold) the moment it exists, so it needs the grid topology NOW ! — before `configure_ocean_bc` below derives it. A staged topology ! overrides the tags axis-by-axis exactly as ! `ocean_bc_state_set_topology` will (an axis it marks periodic is ! periodic; one it does not keeps the namelist's edges). if (engine%has_staged_topology) then block logical :: per_x, per_y per_x = engine%staged_periodic_x .or. & (ocean_bc_type_from_string(cfg%ocean%bc%west) == OBC_PERIODIC .and. & ocean_bc_type_from_string(cfg%ocean%bc%east) == OBC_PERIODIC) per_y = engine%staged_periodic_y .or. & (ocean_bc_type_from_string(cfg%ocean%bc%south) == OBC_PERIODIC .and. & ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_PERIODIC) if (allocated(engine%staged_bathymetry)) then call ocean_state_seed_from_cfg(engine%state, engine%grid, cfg, ierr=ierr, & injected_b=engine%staged_bathymetry, & injected_b_convention=engine%staged_bathymetry_convention, & periodic_x=per_x, periodic_y=per_y) else call ocean_state_seed_from_cfg(engine%state, engine%grid, cfg, ierr=ierr, & periodic_x=per_x, periodic_y=per_y) end if end block else if (allocated(engine%staged_bathymetry)) then call ocean_state_seed_from_cfg(engine%state, engine%grid, cfg, ierr=ierr, & injected_b=engine%staged_bathymetry, & injected_b_convention=engine%staged_bathymetry_convention) else call ocean_state_seed_from_cfg(engine%state, engine%grid, cfg, ierr=ierr) end if if (setup_failed(ierr)) return ! PR-23 sponge: snapshot the reference state (target_source="ic") ! from the JUST-SEEDED initial condition, BEFORE any warm-restart ! read below (load-bearing ordering). call ocean_sponge_snapshot_reference(engine%state%sponge, engine%grid, engine%state%multilayer) ! Overlay tracer metadata + scalar IC values from cfg. call register_default_tracers( & engine%state%multilayer%tracers(engine%state%multilayer%idx_salinity), & engine%state%multilayer%tracers(engine%state%multilayer%idx_temperature), & cfg) ! Dynamic wetting/drying: MUST precede the restart read (persistent ! hysteresis mask must be allocated + registered before the restart ! registry walk). Default off => no-op, byte-identical. call configure_ocean_wetdry(cfg, engine%state, engine%grid, rank) ! The surface-flux component set is allocated BEFORE the restart read: ! the registry walk only sees the components (`sf_heat_cavity`/ ! `sf_salt_cavity`) and the carried assembly (`sf_Q_heat`/`sf_Q_salt`) ! once they exist -- allocated after the read they were written but ! never restored. call engine%state%surface_flux%set_components(engine%grid, & cfg%ocean%forcing%enable_components) ! c069 restart fix: the kappa-shear VERTEX corner carrier ! (`kshear%kd_corner`) MUST be allocated before the restart read too, ! same reason as wetdry/surface_flux above — `ocean_state_restart_read` ! builds its OWN registry (`ocean_state_build_restart_registry`), and ! that registry only sees a field once `allocated(...)` is true. The ! full `configure_ocean_kappa_shear` (which calls `init_vertex`) does ! not run until well after the restart read (it needs `cfg%ocean%vmix` ! validated first), so without this early, knob-only allocation a ! warm-restarted vertex run resumes `kd_corner` at the cold `0` ! `init_vertex` seeds it with — not the checkpointed value — and ! `visc_rem_precompute` (called at the START of the first resumed ! stage, BEFORE `configure_ocean_kappa_shear`'s later re-allocation ! no-ops over the real one) reads that cold `0` as the corner Kv ! source, perturbing `bt_visc_rem_u/v` at the `kappa_trunc` scale and, ! from there, every downstream field (compat row c069). `init_vertex` ! is idempotent (it only allocates when not already allocated), so ! this early call and `configure_ocean_kappa_shear`'s later one are ! harmless duplicates on a cold start. Knob-gated directly off `cfg` ! (no validation performed here) => inert whenever `at_vertex` is off. if (cfg%ocean%kshear%enable .and. cfg%ocean%kshear%at_vertex) then call engine%state%kshear%init_vertex(engine%grid, engine%state%multilayer%nz_ml) end if ! Warm restart (optional — absent/blank restart_file => cold start, ! matching the API/bench callers today). NetCDF-only: filename ! resolution needs `output_rank_filename`. did_restart = .false. if (present(restart_file)) then if (len_trim(restart_file) > 0) then did_restart = .true. engine%warm_restart = .true. #ifndef RDB_NO_NETCDF if (len_trim(restart_file) >= 3 .and. & restart_file(len_trim(restart_file) - 2:len_trim(restart_file)) == ".nc") then restart_filename = trim(restart_file) else restart_filename = output_rank_filename(trim(restart_file), "restart", mrank) end if call ocean_state_restart_read(engine%state, engine%grid, engine%decomp, & trim(restart_filename), t_restart_local, & step_restart_local, ierr=ierr) if (setup_failed(ierr)) return if (rank == 0) then call logger%info("Ocean warm restart from "//trim(restart_filename)// & " at t = "//to_string(t_restart_local)//" s (step "// & to_string(step_restart_local)//")") end if #else call fail("engine_setup: restart_file was given but this build has no "// & "NetCDF support (RDB_ENABLE_NETCDF=OFF)", ierr, OCEAN_STATUS_ERR_SETUP) return #endif end if end if if (present(t_restart)) t_restart = t_restart_local if (present(step_restart)) step_restart = step_restart_local ! Diag-manager: optional z-levels, default variable set, per-rank ! NetCDF stream. NetCDF-only (moved from driver_run_ocean's private ! helper of the same name). ! Under validate_only the stream is never opened, so teardown must ! not try to close it. engine%diag_enabled = cfg%ocean%diag%enabled .and. .not. no_output call engine_configure_diag(cfg, engine%state, mrank, engine%decomp, ierr=ierr, & open_stream_file=.not. no_output) if (setup_failed(ierr)) return ! Surface heat/salt/p_surf/sw-penetration/restore seeding — 2D ! fields device-mapped by ocean_state_enter_data; re-seeded on ! resume (configure-time static, not in the restart registry) -- ! except the ASSEMBLED Q_heat/Q_salt under the component set, which ! are carried state (`ocean_state_build_restart_registry`): when the ! checkpoint holds an assembly, the seed below must not replace it. if (rank == 0 .and. cfg%ocean%forcing%enable_components) then call logger%info("Forcing components: ON") end if resume_q = did_restart .and. engine%state%surface_flux%use_components .and. & engine%state%surface_flux%q_assembled > 0.5_wp if (resume_q) then q_heat_resume = engine%state%surface_flux%Q_heat q_salt_resume = engine%state%surface_flux%Q_salt end if call engine%state%surface_flux%set_surface_flux_const( & cfg%ocean%thermo%q_heat, cfg%ocean%thermo%q_salt) if (resume_q) then engine%state%surface_flux%Q_heat = q_heat_resume engine%state%surface_flux%Q_salt = q_salt_resume end if call engine%state%surface_flux%set_p_surf_const( & cfg%ocean%psurf%p_surf_const) ! E3 top-of-column IN-SITU EOS pressure: seed `ms%p_top` from the ! assembled `sf%p_surf` HERE, on the host and before `enter_data`, so ! the very first PGF of the run already sees the load. The ! outer-step driver refreshes it every step from the same source, ! which is what keeps it live once the ice mass-loading PR makes ! `p_surf` dynamic. Gated: with `in_eos = .false.` `p_top` stays the ! zero array it was allocated as and every EOS evaluation is ! bit-identical. ! ! P5.0 joins `&ocean_pgf_nml p_top_in_bc` — the load in the FV_MOM6 ! pressure-stack surface BC — to the SAME seed, so the unsplit ! driver (which has no per-step refresh) still gets the configure ! value rather than a zero array, and the split driver's step-1 PGF ! is already loaded. if ((cfg%ocean%psurf%in_eos .or. cfg%ocean%pgf%p_top_in_bc) .and. & allocated(engine%state%surface_flux%p_surf)) then engine%state%multilayer%p_top = engine%state%surface_flux%p_surf end if call engine%state%surface_flux%set_sw_penetration( & cfg%ocean%thermo%sw_pen_frac, cfg%ocean%thermo%sw_band_ratio, & cfg%ocean%thermo%sw_zeta1, cfg%ocean%thermo%sw_zeta2, & sw_source=cfg%ocean%thermo%sw_source) call engine%state%surface_flux%set_restore( & cfg%ocean%restore%enable_restore_temp, & cfg%ocean%restore%enable_restore_salt, & cfg%ocean%restore%piston_t, cfg%ocean%restore%piston_s, & cfg%ocean%restore%restore_sst, cfg%ocean%restore%restore_sss) if (rank == 0 .and. engine%state%surface_flux%has_restore_T) then call logger%info("Restore SST: piston = "// & to_string(cfg%ocean%restore%piston_t)//" m/day, target = "// & to_string(cfg%ocean%restore%restore_sst)//" degC") end if if (rank == 0 .and. engine%state%surface_flux%has_restore_S) then call logger%info("Restore SSS: piston = "// & to_string(cfg%ocean%restore%piston_s)//" m/day, target = "// & to_string(cfg%ocean%restore%restore_sss)//" PSU") end if if (rank == 0 .and. & (cfg%ocean%thermo%q_heat /= 0.0_wp .or. cfg%ocean%thermo%q_salt /= 0.0_wp)) then call logger%info("Surface flux: Q_heat = "// & to_string(cfg%ocean%thermo%q_heat)//" W/m2, Q_salt = "// & to_string(cfg%ocean%thermo%q_salt)//" kg/m2/s") end if ! Sea-ice: resume-fold the restart-carried brine/heat/shortwave ! contributions back into the just-reseeded Q_salt/Q_heat (cold ! start: exact +0.0). Not when a carried assembly was resumed above: ! it already holds them (summed in the assembler's order). if (engine%state%ice%enable) then if (.not. resume_q) then engine%state%surface_flux%Q_salt = engine%state%surface_flux%Q_salt & + engine%state%ice%salt_flux_diag engine%state%surface_flux%Q_heat = engine%state%surface_flux%Q_heat & + engine%state%ice%heat_flux_diag & + engine%state%ice%sw_thru_diag end if engine%state%surface_flux%has_salt = .true. engine%state%surface_flux%has_heat = .true. if (engine%state%surface_flux%use_components) then engine%state%surface_flux%has_q_sw = .true. end if end if ! Geothermal bottom heat flux (held on the engine, not a state slot; ! the split-driver's `geo` arg is optional). call engine%geo%init(engine%grid) engine%geo%enable = cfg%ocean%geothermal%enable engine%geo%q_geo_const = cfg%ocean%geothermal%q_geo if (rank == 0 .and. engine%geo%enable .and. cfg%ocean%geothermal%q_geo /= 0.0_wp) then call logger%info("Geothermal flux: Q_geo = "// & to_string(cfg%ocean%geothermal%q_geo)//" W/m2") end if ! P2.5 geometry injection: in-memory supergrid arrays ! (rdb_ocean_stage_metrics) bypass the cfg%ocean%grid%grid_config ! dispatch entirely — metrics_assemble_from_supergrid_arrays is the ! same battle-tested index-sum path metrics_fill_from_supergrid / ! the tripolar generator use, just fed in-memory arrays instead of a ! mosaic file. if (allocated(engine%staged_metrics_x)) then if (size(engine%staged_metrics_x, 1) /= 2*engine%grid%nx_phys + 1 .or. & size(engine%staged_metrics_x, 2) /= 2*engine%grid%ny_phys + 1) then call fail("engine_setup: staged metrics x/y shape mismatch — expected "// & "(2*nx_phys+1, 2*ny_phys+1) = ("//to_string(2*engine%grid%nx_phys + 1)// & ","//to_string(2*engine%grid%ny_phys + 1)//")", ierr, OCEAN_STATUS_ERR_BAD_SHAPE) return end if ! Same ghost-metric topology as the NetCDF mosaic reader: seam ! faces across a periodic seam, then the periodic / fold ghost ! images (`metrics_fold_periodic_ghosts`). Periodicity is the ! staged topology's or the tags', as for the seed above; the fold ! is the north tag's. No grid rotation is staged (`angle_dx` = 0). ! Injected supergrid arrays describe the WHOLE grid (the shape ! check above is against the tile), so they are single-rank. if (csize > 1) then call fail("engine_setup: injected (staged) supergrid metrics are "// & "single-rank — they describe the whole grid, not a tile.", & ierr, OCEAN_STATUS_ERR_SETUP) return end if block logical :: per_x, fold per_x = (ocean_bc_type_from_string(cfg%ocean%bc%west) == OBC_PERIODIC .and. & ocean_bc_type_from_string(cfg%ocean%bc%east) == OBC_PERIODIC) if (engine%has_staged_topology) per_x = per_x .or. engine%staged_periodic_x fold = ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_TRIPOLAR_FOLD call metrics_assemble_from_supergrid_arrays(engine%state%metrics, engine%grid, & engine%staged_metrics_x, engine%staged_metrics_y, & engine%staged_metrics_dx, engine%staged_metrics_dy, & engine%staged_metrics_area, periodic_x=per_x) if (per_x .or. fold) then call metrics_fold_periodic_ghosts(engine%state%metrics, engine%grid, & periodic_x=per_x, north_fold=fold) end if end block call metrics_finalize(engine%state%metrics) if (rank == 0) then call logger%info("Grid config: injected (in-memory supergrid arrays)") end if if (present(ierr)) ierr = OCEAN_STATUS_OK else call configure_ocean_metrics(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return end if ! Equilibrium body-force tide (C1): needs the filled geolatT/geolonT ! from configure_ocean_metrics. call configure_ocean_tides(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return ! Atmospheric surface-pressure loading / inverse barometer (PR-17). call configure_ocean_p_surf(cfg, engine%state, engine%grid, rank) call configure_ocean_forcing(cfg, engine%state, engine%grid, rank, decomp=engine%decomp) ! File-backed surface forcing (PR-15): must run after set_components ! (heat/salt destination) and after configure_ocean_forcing (whose ! formula wind it overrides), before enter_data (registration ! allocates the reader's bracket buffers). NetCDF-only. #ifndef RDB_NO_NETCDF call ocean_data_forcing_configure(cfg%ocean%dataovr, engine%state%data_input, engine%grid, & engine%state%surface_stress, engine%state%surface_flux, & engine%state%bc, engine%state%data_forcing, ierr=ierr) if (setup_failed(ierr)) return #else if (cfg%ocean%dataovr%enable) then call fail("engine_setup: &ocean_dataovr_nml enable = .true. but this build has no "// & "NetCDF support (RDB_ENABLE_NETCDF=OFF)", ierr, OCEAN_STATUS_ERR_SETUP) return end if #endif call configure_ocean_drag(cfg, engine%state, rank, ierr=ierr) if (setup_failed(ierr)) return call configure_ocean_hdiff(cfg, engine%state, rank) call configure_ocean_vmix(cfg, engine%state, rank, ierr=ierr) if (setup_failed(ierr)) return call configure_ocean_tracers(cfg, engine%state, rank) call configure_ocean_lateral(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return ! Sea-ice PR 5: EVP params + the atmospheric-stress snapshot ice ! feels (D7). MANDATED ORDER: the tau_a snapshot must read the ! pristine wind BEFORE the resume apply overwrites ! surface_stress%tau_x/y with the restart-carried blend. if (engine%state%ice%enable .and. engine%state%ice%dynamics) then engine%evp_params = ice_evp_params_from_config(cfg%ocean%ice%p0, cfg%ocean%ice%c0, & cfg%ocean%ice%ec, cfg%ocean%ice%cdw, & cfg%ocean%ice%rho_ocean, & cfg%ocean%ice%del_sh_min_scale, & cfg%ocean%ice%tdamp, cfg%ocean%ice%evp_sub_steps, & cfg%ocean%ice%a_face_stress, & cfg%ocean%ice%cfl_trunc, & cfg%ocean%ice%cfl_trunc_dyn_its, & cfg%ocean%ice%project_ci) engine%state%ice%tau_a_x = engine%state%surface_stress%tau_x engine%state%ice%tau_a_y = engine%state%surface_stress%tau_y call ice_ocean_stress_resume_apply(engine%state%surface_stress, engine%state%ice) end if call configure_ocean_pgf(cfg, engine%state, rank, ierr=ierr) if (setup_failed(ierr)) return ! The single configured ρ₀ (`&ocean_ic_nml rho_0` -> `eos%rho0`) out to ! the slots that keep their own copy: surface heat/salt flux, wind ! stress, KPP-vmix and the engine-held geothermal slot. Runs after ! every slot-specific configure (nothing downstream re-derives a ! reference density) and well before `ocean_state_enter_data`, which is ! what makes the on-device `vmix%rho0` reads correct without an ! explicit `!$acc update device` — see the routine's docstring. call configure_ocean_reference_density(engine%state, geo=engine%geo) call configure_ocean_bt(cfg, engine%state, engine%grid, rank) call configure_ocean_bt_split(cfg, engine%state, engine%grid, rank) ! may auto-set n_inner ! BC config: edge tags + tidal constituents + tracer inflow values. call configure_ocean_bc(cfg, engine%state, rank, ierr=ierr) if (setup_failed(ierr)) return ! Physical-domain-edge flags: at px*py=1 all four are .true. call ocean_bc_state_set_edges(engine%state%bc, engine%decomp%has_west, & engine%decomp%has_east, engine%decomp%has_south, & engine%decomp%has_north) ! P2.5 geometry injection: Oceananigans-style "the grid owns ! periodicity" (rdb_ocean_stage_topology) — overrides whatever ! configure_ocean_bc just derived from the namelist edge tags. Must ! run BEFORE the periodic ghost wrap / ocean_halo_init / ! configure_ocean_land_mask below, all of which read periodic_x/y. if (engine%has_staged_topology) then call ocean_bc_state_set_topology(engine%state%bc, engine%staged_periodic_x, & engine%staged_periodic_y, ierr=ierr) if (setup_failed(ierr)) return end if ! PR-23 sponge: build the per-cell idamp maps from the just-configured ! edge tags. Must run AFTER configure_ocean_bc, BEFORE enter_data. call configure_ocean_sponge(cfg, engine%state, engine%grid, rank) ! Seed the boundary data source from the same config so the per-step ! update is idempotent with configure_ocean_bc. engine%bc_source%u_west = cfg%ocean%bc%west_clamped_u engine%bc_source%u_east = cfg%ocean%bc%east_clamped_u engine%bc_source%v_south = cfg%ocean%bc%south_clamped_v engine%bc_source%v_north = cfg%ocean%bc%north_clamped_v engine%bc_source%eta_west = cfg%ocean%bc%west_clamped_eta engine%bc_source%eta_east = cfg%ocean%bc%east_clamped_eta engine%bc_source%eta_south = cfg%ocean%bc%south_clamped_eta engine%bc_source%eta_north = cfg%ocean%bc%north_clamped_eta ! Init-time periodic ghost wrap: after topo+IC are seeded and BC ! tags are configured, wrap bathymetry + multilayer state ghost ! cells so all kernels see correct periodic seam values on the ! first step. Run on the HOST here (before enter_data). ! ! Distributed fold (`px > 1`, this rank folds the north edge): the ! folds below need the routing plan, which `ocean_fold_exchange_init` ! builds further down, so on that path they are SKIPPED here and run ! after the host halo pass instead (same fields, same cold-start-only ! rule for the prognostics). `px = 1` is untouched. dist_fold = engine%state%bc%north_fold .and. engine%decomp%px > 1 if (engine%state%bc%periodic_x .or. engine%state%bc%periodic_y .or. & engine%state%bc%north_fold) then call ocean_periodic_wrap_centre_2d( & engine%state%barotropic%b, & engine%grid%nx_total, engine%grid%ny_total, & engine%grid%nx_phys, engine%grid%ny_phys, engine%grid%nghost, & engine%state%bc%periodic_x, engine%state%bc%periodic_y) ! The PROGNOSTIC state is wrapped on a cold start only: a ! checkpoint holds FULL local arrays, ghosts included, exactly as ! the run that wrote it carried them into its next step (see the ! ghost-cell policy in `ocean_state_build_restart_registry`). ! Those ghosts are not a pure function of the owned interior at a ! step boundary -- the duplicated seam faces of `u` hold each ! tile's own update -- so re-deriving them here resumed a ! DIFFERENT state from the one the writer stepped on (1/4-degree ! Southern Ocean, 4x1 periodic, 2026-10-01). Static geometry ! (`b`, the draft, `bt_H_ref`) is wrapped either way. if (.not. did_restart) then call ocean_periodic_wrap_state(engine%grid, engine%state%bc, engine%state%multilayer) end if if (.not. dist_fold) then call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, engine%state%barotropic%b) end if if (.not. did_restart .and. .not. dist_fold) then call ocean_fold_wrap_state(engine%grid, engine%state%bc, engine%state%multilayer) end if ! The ice draft is bathymetry-class static geometry, so it takes ! the bathymetry's ghost treatment VERBATIM: the analytic setter ! already filled the ghosts, and the wrap/fold then overwrites ! them with the seam-correct values on a periodic or folded edge. ! `bt_H_ref = b - z_draft` was latched from the UNWRAPPED pair, ! which is exactly why all three are re-wrapped here (and why ! they must be re-wrapped TOGETHER — a draft whose seam disagreed ! with the datum's would count the ice load twice at that face). if (engine%state%metrics%use_cavity) then call ocean_periodic_wrap_centre_2d( & engine%state%metrics%z_draft, & engine%grid%nx_total, engine%grid%ny_total, & engine%grid%nx_phys, engine%grid%ny_phys, engine%grid%nghost, & engine%state%bc%periodic_x, engine%state%bc%periodic_y) if (.not. dist_fold) call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, & engine%state%metrics%z_draft) call ocean_periodic_wrap_centre_2d( & engine%state%metrics%cover_frac, & engine%grid%nx_total, engine%grid%ny_total, & engine%grid%nx_phys, engine%grid%ny_phys, engine%grid%nghost, & engine%state%bc%periodic_x, engine%state%bc%periodic_y) if (.not. dist_fold) call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, & engine%state%metrics%cover_frac) end if ! bt_H_ref was snapshotted from the UNWRAPPED b inside ! configure_ocean_bt_split (above) — re-wrap it too. if (cfg%ocean%bt%n_inner >= 1) then call ocean_periodic_wrap_centre_2d( & engine%state%dyn%bt_work%bt_H_ref, & engine%grid%nx_total, engine%grid%ny_total, & engine%grid%nx_phys, engine%grid%ny_phys, engine%grid%nghost, & engine%state%bc%periodic_x, engine%state%bc%periodic_y) if (.not. dist_fold) call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, & engine%state%dyn%bt_work%bt_H_ref) end if end if ! Initialise the ocean halo module. Must run BEFORE any ocean_halo_* ! call, with the BC periodicity flags already set. call ocean_halo_init(engine%decomp, engine%grid%nghost, & engine%state%bc%periodic_x, engine%state%bc%periodic_y, ierr=ierr) if (setup_failed(ierr)) return ! The sea-ice category exchanges batch `ncat*nk_ice` (and ! `ncat+1`) levels per message: reserve for the widest so the ! buffers never regrow mid-run (a device reallocation breaks UCX IPC ! handle reuse). call ocean_halo_reserve(merge(max(cfg%nz_layers, cfg%ocean%ice%ncat*cfg%ocean%ice%nk_ice, & cfg%ocean%ice%ncat + 1), cfg%nz_layers, & cfg%ocean%ice%enable), & merge(engine%grid%nghost + cfg%ocean%bt%bt_halo, 0, & cfg%ocean%bt%bt_halo > 0), ierr=ierr) if (setup_failed(ierr)) return ! Distributed tripolar fold (px > 1): routing plan + buffers, built ! from the same decomposition. Inactive (no plan, no buffers) unless ! this rank folds the north edge of an east-west split, so a no-op on ! every px = 1 run. Reserve the largest group (the `ml_state` one: ! h, u, v and every tracer, at most ng+1 rows each). call ocean_fold_exchange_init(engine%decomp, engine%grid%nghost, & engine%state%bc%north_fold, ierr=ierr) if (setup_failed(ierr)) return block integer :: ntr ntr = 0 if (allocated(engine%state%multilayer%tracers)) ntr = size(engine%state%multilayer%tracers) call ocean_fold_exchange_reserve((engine%grid%nghost + 1)*cfg%nz_layers*(3 + ntr)) end block ! Host-side seam ghost fill (D0 init-halo, O2): single-rank ! non-periodic ⇒ no-op, periodic ⇒ local wrap. call ocean_halo_centre(engine%state%barotropic%b, device_resident=.false.) if (engine%state%metrics%use_cavity) then call ocean_halo_centre(engine%state%metrics%z_draft, device_resident=.false.) call ocean_halo_centre(engine%state%metrics%cover_frac, device_resident=.false.) end if ! Prognostic halo: cold start only, for the reason given at the wrap ! above -- a checkpoint already carries the writer's halo columns. if (.not. did_restart) then call ocean_halo_exchange_ml_state(engine%state%multilayer, device_resident=.false.) end if if (cfg%ocean%bt%n_inner >= 1) then call ocean_halo_centre(engine%state%dyn%bt_work%bt_H_ref, device_resident=.false.) end if ! px > 1 tripolar: refresh the static corner Coriolis array's x ghosts ! from their OWNERS before it is folded. The analytic generator ! evaluates `2*omega*sin(geolatBu)` over the whole local array, and ! gfortran -O3 -march=native vectorises that loop through libmvec: ! the scalar remainder lanes differ from the vector lanes by an ulp, ! and which columns are remainder depends on the TILE width. The tail ! lands in the ghost columns (owned f is unaffected), but those ghosts ! then differed from the serial run, whose ghosts are periodic copies ! of computed interior values. A corner array is face-type in x and, ! row-for-row, the same as a centre array in y (corner row j is the SW ! corner of T row j), so the face-x primitive over rows 1..ny_total is ! the correct exchange; the north-most row (ny_total+1) is a fold row ! on the north tile and beyond every stencil elsewhere. Collective ! over every rank (rank-uniform condition), host mode. if (engine%decomp%px > 1 .and. & ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_TRIPOLAR_FOLD) then block real(wp), allocatable :: f_rows(:, :) f_rows = engine%state%coriolis_adv%f_corner(:, 1:engine%grid%ny_total) call ocean_halo_face_x(f_rows, device_resident=.false.) engine%state%coriolis_adv%f_corner(:, 1:engine%grid%ny_total) = f_rows end block end if ! The deferred init-time folds of the distributed (`px > 1`) path — ! exchange → periodic wrap → fold, as at every step-time seam site; ! host mode, before the device map. Collective over the north rank ! row (every rank of it has `dist_fold` set). The fold is ! owner-routed, so it reads only owned points and writes every column ! of the north ghost rows: its place after the halo cannot change a ! bit. Prognostics on a cold start only (a checkpoint carries the ! writer's ghosts, see the wrap above); static geometry either way. if (dist_fold) then call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, engine%state%barotropic%b, & device_resident=.false.) if (.not. did_restart) then call ocean_fold_wrap_state(engine%grid, engine%state%bc, engine%state%multilayer, & device_resident=.false.) end if if (engine%state%metrics%use_cavity) then call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, & engine%state%metrics%z_draft, device_resident=.false.) call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, & engine%state%metrics%cover_frac, device_resident=.false.) end if if (cfg%ocean%bt%n_inner >= 1) then call ocean_fold_wrap_eta_2d(engine%grid, engine%state%bc, & engine%state%dyn%bt_work%bt_H_ref, device_resident=.false.) end if ! The static corner Coriolis array (plan site S3): on px = 1 ! `fill_f_corner_seam_ghosts` (configure_ocean_forcing) folds it — ! including the fold-line row projection, which is what makes the ! two copies of a fold-line vertex carry the same f bitwise even for ! planetary f. It returns early on px > 1 (its x ghosts come from ! the generator's global geography), so the fold — a scalar COPY, ! planetary or beta-plane alike — is done here, now that the plan ! exists. Nothing reads f_corner's north ghosts or fold row before ! this point (the BT CFL and the stability audit do not use it). call ocean_fold_north_corner(engine%state%coriolis_adv%f_corner, & engine%grid%nx_total + 1, engine%grid%ny_total + 1, & engine%grid%nx_phys, engine%grid%ny_phys, & engine%grid%nghost, .false., device_resident=.false.) end if ! The PGF keeps its OWN copy of the bathymetry (FV-MOM6 and gprime ! place the bottom interface from it at every face). Take it HERE, ! from the wrapped + halo-exchanged `b`, not in `configure_ocean_pgf`: ! a copy taken before this point froze whatever the seam ghosts held ! then — constant-extrapolated edge columns for a file/staged ! bathymetry, which on a periodic edge is the B2 seam jet. Before ! `enter_data` (the device copy is taken from the host values). call engine%state%pressure_force%set_bathymetry(engine%state%barotropic%b) ! The isopycnal-slopes slot keeps its own copy too: it is the bed ! datum of the geopotential interface heights whose across-face ! difference is the interface-tilt term, and a seam face reads the ! GHOST column — so it is taken from the same wrapped + halo- ! exchanged `b`, for the same reason. No-op when the slot is off. call engine%state%slopes%set_bathymetry(engine%state%barotropic%b) ! Same for the ZSTAR_FULL per-column reference table the seed built ! from `b`: rebuilt from the halo-exchanged field so a DECOMPOSED ! axis's seam ghost columns are right too (a pure function of `b` — ! identical wherever the seed already saw the right ghosts). call engine%state%vcoord%build_zref_full(engine%state%barotropic%b) ! Static land masking: derive the C-grid face/corner masks from the ! seeded wet_mask + zero the 6 face metrics at land faces. call configure_ocean_land_mask(cfg, engine%state, engine%grid, rank, & warm_restart=did_restart) ! Configure-time stability audit (viscous CFL / kappa_h diffusive ! number / Munk-layer resolution / ah_max-clamps-nu_h): MUST run ! after configure_ocean_metrics + configure_ocean_land_mask (needs ! the real per-cell metric arrays, not nominal &grid_nml dx/dy) and ! before enter_data. See rdb_ocean_stability_audit.F90 for the ! motivating failure (a global tripolar aquaplanet NaN, diagnosed ! only after the fact — this audit is the fix). ! `bt_H_ref` (b - z_draft afloat, 0 where grounded) is the COLUMN the ! vertical coordinate divides, so it — not the bathymetry — is what ! the terrain-following stiffness check must see under an ice shelf. ! It exists only once the barotropic datum has been built. if (cfg%ocean%bt%n_inner >= 1) then call ocean_stability_audit(cfg, engine%state%metrics, engine%grid, rank, & ierr=ierr, & column=engine%state%dyn%bt_work%bt_H_ref) else call ocean_stability_audit(cfg, engine%state%metrics, engine%grid, rank, ierr=ierr) end if if (setup_failed(ierr)) return ! Barotropic linear wave drag: host-side r_H map + h->face average. ! AFTER bathymetry + land masking, BEFORE enter_data. call configure_ocean_wave_drag(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return ! Porous barriers: grow + fill the static along-face subgrid ! statistics. Same ordering constraints as wave drag. call configure_ocean_porous(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return ! Ice-shelf cavity: build the static isostatic load from the draft ! the IC seed already laid down, and assert the counted-once datum ! invariant. AFTER configure_ocean_pgf (it needs the PGF reference ! density) and configure_ocean_bt_split (it checks that latch), ! BEFORE enter_data. call configure_ocean_cavity(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return ! Ice-shelf basal melt (P2b): copy the thermodynamic knobs onto the ! melt slot, resolve the gamma_s sentinel, build the per-column ! Coriolis array the hj99 law needs, and SEED ms%p_top from the ! isostatic load configure_ocean_cavity just built. Immediately ! after it (that is where p_ice_ref comes from), BEFORE enter_data. call configure_ocean_cavity_melt(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return ! Ice-shelf TOP drag (Phase 4a): coefficients + the static FACE ! cover masks projected from `metrics%cover_frac`. AFTER the ! cover_frac halo exchange and after configure_ocean_cavity_melt ! (it owns the one-C_d rule against the melt slot's cdrag_top), ! BEFORE enter_data (the face masks are host-filled and reach the ! device on the slot's `copyin` map). call configure_ocean_top_drag(cfg, engine%state, engine%grid, rank) ! Partial-step z-level face closure (&vcoord_nml ! zfixed_closed_faces). LAST of the static-geometry builders and ! still BEFORE enter_data: it needs `vcoord%z_top` (configure_ocean ! _cavity), `bt_work%bt_H_ref` (configure_ocean_bt_split), the ! land-masked `dy_cu`/`dx_cv` (configure_ocean_land_mask) and the ! periodic-wrap + halo pass above — the mask is built from the ! z_fixed target at eta = 0, and its ghost-band correctness IS the ! seam correctness of those two inputs. Knob off => literal no-op. call configure_ocean_closed_faces(cfg, engine%state, engine%grid, rank, ierr=ierr) if (setup_failed(ierr)) return ! The shared FIRST-LIVE-LAYER index `ms%k_top` (+ its two face ! twins) — what every top-side consumer reads instead of spelling ! `nz`, so that a `z_fixed` column whose top layers are inert ! fillers inside the ice draft is forced on the ice-adjacent LIVE ! layer and not on the filler. Same inputs and same ordering ! constraints as the closed-face mask above (it is built from the ! same `z_fixed` target at eta = 0), still before enter_data. ! Literal no-op on every coordinate but `z_fixed` under a cavity: ! the arrays already hold the `nz` fallback. call configure_ocean_k_top(cfg, engine%state, engine%grid, rank) ! Its bed-side mirror `ms%k_bot` (+ face twins, `max` rule): the ! first LIVE layer counting UP from the bed, read by every bed-side ! consumer (bottom drag, the vdiff bed row, geothermal, tidal-mixing ! bed anchor, MEKE bed speed, bed-reaching shortwave) instead of ! `k = 1`. Same inputs and ordering as `k_top`, but NOT gated on a ! cavity: every z_fixed column shallower than the nominal stack has ! bed fillers. Literal no-op off z_fixed (arrays hold the `1` ! fallback). call configure_ocean_k_bot(engine%state, engine%grid, rank) ! Sea-ice PR 24: analytic IC path. Host-side, run once, AFTER ! wet_mask/geolatT/wet_T are valid, BEFORE enter_data. Skips on a ! warm restart (the restart read already replaced the IC). if (engine%state%ice%enable) then if (.not. did_restart) then engine%ic_par = ice_ic_params_from_config(cfg%ocean%ice_ic%conc_config, & cfg%ocean%ice_ic%conc, & cfg%ocean%ice_ic%h_ice, cfg%ocean%ice_ic%h_snow, & cfg%ocean%ice_ic%t_ice, cfg%ocean%ice_ic%s_ice, & cfg%ocean%ice_ic%arctic_edge, & cfg%ocean%ice_ic%antarctic_edge) call ice_init_apply(engine%grid, engine%state%multilayer, engine%state%metrics, & engine%state%ice, engine%ic_par) ! X1 at cold start (host, before `enter_data`): the analytic ! IC is a pointwise formula, but its ghost band is only as ! right as the inputs' ghosts — make the category state the ! owner's by construction. Never on a warm restart: the ! checkpoint carries the writer's ghosts (8e1931f20). call ocean_halo_exchange_ice_state(engine%state%ice, engine%grid, engine%state%bc, & device_resident=.false.) if (rank == 0 .and. engine%ic_par%conc_config /= ICE_IC_CONC_ZERO) then call logger%info("Sea-ice IC: conc_config='"// & trim(cfg%ocean%ice_ic%conc_config)// & "', h_ice = "//to_string(cfg%ocean%ice_ic%h_ice)// & " m, conc = "//to_string(cfg%ocean%ice_ic%conc)) end if end if end if ! Resolve + store n_inner (does NOT fail loud here — see the ! n_inner docstring on ocean_engine_t; a caller wanting the C ! ABI's stricter "n_inner must be >= 1" contract checks it itself). engine%n_inner = cfg%ocean%bt%n_inner ! Phase-3 vanish_tol for console MaxCFL gating. engine%cfl_vtol = 0.0_wp if (engine%state%dyn%cfl_ignore_vanished) then if (engine%state%vcoord%coord_type == VCOORD_LAGRANGIAN) then engine%cfl_vtol = isopycnal_vanish_tol(engine%state%dyn%angstrom_h) end if end if engine%is_setup = .true. end subroutine engine_setup