The shared time-varying NetCDF input reader (PR-14).
!! The shared time-varying NetCDF input reader (PR-14). module rdb_ocean_data_input !! Generic, decomposition-aware, device-resident reader for !! time-varying NetCDF fields `field(x, y[, z], t)`. Owns the file !! handles, the time axis, the bracket bookkeeping, the H->D motion !! (only when the bracket advances), and the per-step linear-in-time !! blend kernel. Owns NO field names and NO state-slot mapping: a !! consumer registers its own `(file, variable, destination)` triple !! via `ocean_data_input_register_2d/_3d` (or the edge-segment !! variants for OBC files) and gets back an opaque `id`. !! !! **The binding seam (see `PLAN_PR14_netcdf_input_reader.md` §13):** !! !! * The reader never touches state. It blends on-device into the !! caller's own, caller-MAPPED, explicit-shape, WHOLE (never a !! section) array, at a registration-time offset !! `(dest_i0, dest_j0)`. Passing a section into an explicit-shape !! dummy would trigger a host copy-in/out and silently detach !! from the device map — never do that. !! * The consumer owns the halo exchange (physical cells only are !! filled) and the vertical (k) orientation — a source file's z !! axis is read AS STORED; nothing here flips it. A 3-D field !! read from a surface-first source file is handed back !! bottom-last; if that's wrong for your consumer, flip it there. !! * `&ocean_data_nml` carries no per-field entries — only !! `max_fields`/`verbose`. Each consumer defines its own file/ !! variable keys in its own namelist group. !! * Out-of-range time is an ABORT by default (`DATA_OOR_ERROR`); !! `DATA_OOR_CLAMP` is opt-in per field. !! * No in-core horizontal interpolation (files are pre-regridded !! to the model's global horizontal extent), no vertical !! remapping of a source z axis, no calendar. Time interpolation !! is linear between the two bracketing records; `DATA_TIME_CYCLIC` !! wraps a climatology through an explicit period; `DATA_TIME_STATIC` !! reads record 1 once and never touches the file again. !! * Time interpolation is a PURE function of `t` — no reader-side !! restart state. A resumed run reproduces every field !! bit-for-bit from `t_current` alone. !! * `ocean_data_input_update_all(this, t)` refreshes every !! registered field's bracket (+ pushes any new slab to the !! device) and is the driver's one-line per-step hook — a no-op !! when `nfields == 0` (every shipped namelist today). It does !! NOT write into any consumer array (the reader stores the !! registration-time `dest` shape as metadata only, never a !! pointer). Each consumer then calls `update_2d/_3d` with its !! OWN array at its own point in the step, and that call is !! checked against the `t` `update_all` most recently refreshed !! with — calling `update_2d/_3d` before `update_all` has run !! this step's `t` is a fail-loud ordering bug, not a silent !! stale read. !! * `ocean_data_input_fill_static_host` is a one-shot HOST-side !! fill for `DATA_TIME_STATIC` fields, legal at SETUP time before !! the destination is device-mapped (PR-24 ice IC, PR-30 tidal !! maps). Touches no device memory; fails loud on a non-static !! field. !! !! Storage order: files must be FORTRAN-ordered `(x, y[, z], t)` — !! v1 detects a horizontal-dimension mismatch and aborts rather than !! silently transposing (a transposed field is exactly the class of !! bug the analytical tests are built to catch); regrid/reorder the !! file if this trips. use, intrinsic :: iso_fortran_env, only: int64 use rdb_constants, only: wp use rdb_grid, only: hgrid_t use rdb_config, only: ocean_data_config_t use rdb_io_netcdf, only: nc_check, nc_open_read, nc_close, nc_get_var_1d, & nc_get_var_slab_3d, nc_get_att_text use netcdf, only: nf90_noerr, nf90_inq_varid, nf90_inquire_variable, & nf90_inquire_dimension use rdb_mem_report, only: arr_bytes use pic_logger, only: logger => global_logger use pic_strings, only: to_string use pic_ascii, only: to_lower use rdb_ocean_status, only: OCEAN_STATUS_OK, OCEAN_STATUS_ERR_SETUP, OCEAN_STATUS_ERR_IO use rdb_error_ring, only: fail implicit none private public :: ocean_data_input_t public :: data_input_field_t ! --- time modes --- public :: DATA_TIME_LINEAR, DATA_TIME_CYCLIC, DATA_TIME_STATIC ! --- out-of-range policy --- public :: DATA_OOR_ERROR, DATA_OOR_CLAMP ! --- edge tags (segment registration, PR-22) --- public :: DATA_EDGE_WEST, DATA_EDGE_EAST, DATA_EDGE_SOUTH, DATA_EDGE_NORTH public :: ocean_data_input_register_2d public :: ocean_data_input_register_3d public :: ocean_data_input_register_segment_2d public :: ocean_data_input_register_segment_3d public :: ocean_data_input_update_2d public :: ocean_data_input_update_3d public :: ocean_data_input_fill_static_host public :: ocean_data_input_fill_static_host_3d public :: ocean_data_input_load_static_2d public :: ocean_data_input_update_all public :: data_input_locate public :: data_input_time_scale_from_units public :: data_input_time_mode_is_implemented public :: data_input_time_mode_from_string public :: data_input_dims_ok integer, parameter :: DATA_TIME_LINEAR = 1 integer, parameter :: DATA_TIME_CYCLIC = 2 integer, parameter :: DATA_TIME_STATIC = 3 integer, parameter :: DATA_OOR_ERROR = 1 integer, parameter :: DATA_OOR_CLAMP = 2 integer, parameter :: DATA_EDGE_WEST = 1 integer, parameter :: DATA_EDGE_EAST = 2 integer, parameter :: DATA_EDGE_SOUTH = 3 integer, parameter :: DATA_EDGE_NORTH = 4 type :: data_input_field_t !! One registered `(file, variable, destination-shape)` triple. !! Array-of-derived-type element — NEVER dereference this from !! inside a `do concurrent`/device-kernel body (outer-shim + !! flat-impl rule); every place this type is touched below is a !! host-side registry walk that hands plain explicit-shape arrays !! to a `pure` flat-impl kernel. logical :: active = .false. integer :: ncid = -1 integer :: varid = -1 logical :: is_3d = .false. integer :: time_mode = DATA_TIME_LINEAR integer :: oor = DATA_OOR_ERROR logical :: oor_warned = .false. ! Source-slab extents (this rank's own subdomain slab, source-k ! order, source-z count for a 3-D field, else 1). integer :: nx = 0, ny = 0, nz = 1 ! File-side 1-based start offsets for THIS rank's horizontal slab. integer :: i0 = 1, j0 = 1 ! Where the slab lands in the consumer's `dest` array. integer :: dest_i0 = 1, dest_j0 = 1 ! Shape metadata of the consumer's `dest` — validated, not stored ! as a pointer (see module docstring). integer :: dest_n1 = 0, dest_n2 = 0, dest_n3 = 1 real(wp) :: scale = 1.0_wp real(wp) :: add_offset = 0.0_wp real(wp) :: t_offset = 0.0_wp real(wp) :: cycle_period = 0.0_wp integer :: nt = 0 real(wp), allocatable :: t_axis(:) !! File time axis, converted to seconds (scale + t_offset NOT !! applied here — t_offset is applied to the QUERY time, not !! the axis; see `data_input_refresh_brackets`). real(wp), allocatable :: f0(:, :, :), f1(:, :, :) !! Bracketing records, shape `(nx, ny, nz)`. Device-mapped by !! `ocean_data_input_t%enter_data`. integer :: rec0 = -1, rec1 = -1 !! Currently-loaded record indices (1-based file record !! numbers); -1 = nothing loaded yet. real(wp) :: w = 0.0_wp !! Current blend weight: `f = (1-w)*f0 + w*f1`. real(wp) :: t_last = -huge(1.0_wp) !! Model time `update_all` most recently refreshed this field !! at. `update_2d/_3d` asserts the caller's `t` matches (fail !! loud on an ordering bug — see module docstring). Ignored !! for `DATA_TIME_STATIC` (never refreshed after registration). integer :: nreads = 0 !! Slab-read counter (test hook — T4 asserts exactly one read !! for a STATIC field's whole run). end type data_input_field_t type :: ocean_data_input_t logical :: is_init = .false. logical :: verbose = .false. integer :: nfields = 0 integer :: nfields_max = 0 type(data_input_field_t), allocatable :: fields(:) contains procedure, non_overridable :: init => ocean_data_input_init procedure, non_overridable :: destroy => ocean_data_input_destroy procedure, non_overridable :: enter_data => ocean_data_input_enter_data procedure, non_overridable :: exit_data => ocean_data_input_exit_data procedure, non_overridable :: bytes => ocean_data_input_bytes end type ocean_data_input_t ! Module-level host read-scratch (registration-time-only allocation ! per CLAUDE.md rule 5 — no per-step allocate on `-stdpar=gpu`; this ! buffer never touches the device, so it carries no `!$acc` directives ! at all). Copy of the `remap_workspace_ensure` idiom ! (`rdb_ml_dynamics.F90`). real(wp), allocatable :: data_input_ws(:, :, :) contains ! ====================================================================== ! Pure, independently-testable helpers ! ====================================================================== pure subroutine data_input_locate(t_axis, nt, mode, cycle_period, t, n0, n1, w, out_of_range) !! Bracket search + blend weight for a query time `t` (already in !! file-time units/offset — the caller applies `t_offset`/`t_scale` !! before calling). `nt == 1` is degenerate: always returns !! `n0 = n1 = 1`, `w = 0`, never out of range. integer, intent(in) :: nt, mode real(wp), intent(in) :: t_axis(nt) real(wp), intent(in) :: cycle_period, t integer, intent(out) :: n0, n1 real(wp), intent(out) :: w logical, intent(out) :: out_of_range integer :: k real(wp) :: t_eff, gap out_of_range = .false. if (nt <= 1) then n0 = 1 n1 = 1 w = 0.0_wp return end if if (mode == DATA_TIME_CYCLIC) then t_eff = t_axis(1) + modulo(t - t_axis(1), cycle_period) if (t_eff >= t_axis(nt)) then ! The seam bracket (nt, 1). Gap is NOT the mean record ! spacing — it is the wrap distance from the last record to ! the first record of the NEXT cycle. n0 = nt n1 = 1 gap = (t_axis(1) + cycle_period) - t_axis(nt) if (gap > 0.0_wp) then w = (t_eff - t_axis(nt))/gap else w = 0.0_wp end if else n0 = 1 do k = 1, nt - 1 if (t_eff >= t_axis(k) .and. t_eff < t_axis(k + 1)) then n0 = k exit end if end do n1 = n0 + 1 w = (t_eff - t_axis(n0))/(t_axis(n1) - t_axis(n0)) end if else ! DATA_TIME_LINEAR (STATIC never reaches here — the caller ! returns before calling this for a static field). if (t < t_axis(1)) then n0 = 1 n1 = 1 w = 0.0_wp out_of_range = .true. else if (t > t_axis(nt)) then n0 = nt n1 = nt w = 0.0_wp out_of_range = .true. else if (t >= t_axis(nt)) then n0 = nt n1 = nt w = 0.0_wp else n0 = 1 do k = 1, nt - 1 if (t >= t_axis(k) .and. t < t_axis(k + 1)) then n0 = k exit end if end do n1 = n0 + 1 w = (t - t_axis(n0))/(t_axis(n1) - t_axis(n0)) end if end if end subroutine data_input_locate pure subroutine data_input_time_scale_from_units(units, scale, ok) !! CF `units` attribute ("seconds since ...", "hours since ...", ...) !! -> a multiplier converting the raw file time axis to seconds. !! Only the leading unit word matters (no calendar, no reference !! date — see module docstring). `ok = .false.` for an !! unrecognised leading word; caller decides the fallback. character(len=*), intent(in) :: units real(wp), intent(out) :: scale logical, intent(out) :: ok character(len=:), allocatable :: u u = to_lower(adjustl(units)) ok = .true. if (index(u, "second") == 1 .or. index(u, "sec") == 1) then scale = 1.0_wp else if (index(u, "minute") == 1 .or. index(u, "min") == 1) then scale = 60.0_wp else if (index(u, "hour") == 1) then scale = 3600.0_wp else if (index(u, "day") == 1) then scale = 86400.0_wp else scale = 1.0_wp ok = .false. end if end subroutine data_input_time_scale_from_units pure logical function data_input_time_mode_is_implemented(tag) result(ok) !! `.true.` iff `tag` (case-insensitive) names a shipped time mode. !! Drives the fail-loud dispatch in `data_input_time_mode_from_string` !! — house idiom, see `lateral_closure_is_implemented`. character(len=*), intent(in) :: tag select case (to_lower(trim(tag))) case ("linear", "cyclic", "static") ok = .true. case default ok = .false. end select end function data_input_time_mode_is_implemented function data_input_time_mode_from_string(tag) result(mode) !! Translate a namelist/registration-time string into a !! `DATA_TIME_*` code. Fail-loud on anything not covered by !! `data_input_time_mode_is_implemented` — an unrecognised tag !! must never silently fall back to a default mode. character(len=*), intent(in) :: tag integer :: mode select case (to_lower(trim(tag))) case ("linear") mode = DATA_TIME_LINEAR case ("cyclic") mode = DATA_TIME_CYCLIC case ("static") mode = DATA_TIME_STATIC case default call logger%error("ocean_data_input: unimplemented time_mode '"//trim(tag)// & "' (implemented: linear, cyclic, static)") error stop "ocean_data_input: unimplemented time_mode" end select end function data_input_time_mode_from_string pure logical function dims_geometry_ok(nd, expect_nd, d_horiz1, d_horiz2, & expect1, expect2) result(ok) !! Low-level geometry check used inline by `register_common`: !! `.true.` iff the variable has the expected rank and its two !! validated-length dims are each at least as long as required !! (decomposition-safe: a subdomain slab only needs `offset+extent` !! to fit, not to equal the file's global length). Never aborts. integer, intent(in) :: nd, expect_nd, d_horiz1, d_horiz2, expect1, expect2 ok = (nd == expect_nd) .and. (d_horiz1 >= expect1) .and. (d_horiz2 >= expect2) end function dims_geometry_ok function data_input_dims_ok(filename, var, nx_phys, ny_phys, is_3d, nz_src) result(ok) !! Self-contained, non-erroring dimension validator — the !! single-rank (`i_offset_global = j_offset_global = 0`) testable !! twin of `register_common`'s inline checks; mirrors !! `zinit_dims_ok`. Opens `filename`, inspects `var`'s rank and !! lengths, and returns `.false.` on any mismatch, missing !! variable, or I/O error — NEVER aborts, so test code can probe !! the false branch without a subprocess/death-test harness (the !! house convention — see `zinit_dims_ok`). use netcdf, only: nf90_open, nf90_close, nf90_nowrite character(len=*), intent(in) :: filename, var integer, intent(in) :: nx_phys, ny_phys logical, intent(in) :: is_3d integer, intent(in), optional :: nz_src logical :: ok integer :: ncid, varid, ierr, var_ndims, expect_nd integer :: var_dimids(4) integer :: d1, d2, d3 ok = .false. ierr = nf90_open(trim(filename), nf90_nowrite, ncid) if (ierr /= nf90_noerr) return expect_nd = merge(4, 3, is_3d) ierr = nf90_inq_varid(ncid, trim(var), varid) if (ierr == nf90_noerr) then ierr = nf90_inquire_variable(ncid, varid, ndims=var_ndims, & dimids=var_dimids(1:expect_nd)) end if if (ierr == nf90_noerr .and. var_ndims /= expect_nd) ierr = -1 if (ierr == nf90_noerr) ierr = nf90_inquire_dimension(ncid, var_dimids(1), len=d1) if (ierr == nf90_noerr) ierr = nf90_inquire_dimension(ncid, var_dimids(2), len=d2) if (is_3d .and. ierr == nf90_noerr) then ierr = nf90_inquire_dimension(ncid, var_dimids(3), len=d3) end if if (ierr == nf90_noerr) then ok = (d1 >= nx_phys) .and. (d2 >= ny_phys) if (is_3d .and. present(nz_src)) ok = ok .and. (d3 == nz_src) end if ierr = nf90_close(ncid) end function data_input_dims_ok ! ====================================================================== ! Blend kernels — pure, explicit-shape, integer dims declared first ! (decl-order), a single `do concurrent` each. Called ONLY from the ! outer-shim `update_2d/_3d` below with plain arrays pulled out of the ! registry — never pass `this%fields(id)%f0` etc. into a kernel body. ! ====================================================================== pure subroutine data_input_blend_2d_impl(nxs, nys, dn1, dn2, i0, j0, f0, f1, w, dest) integer, intent(in) :: nxs, nys, dn1, dn2, i0, j0 real(wp), intent(in) :: f0(nxs, nys), f1(nxs, nys) real(wp), intent(in) :: w real(wp), intent(inout) :: dest(dn1, dn2) integer :: i, j do concurrent(j=1:nys, i=1:nxs) dest(i0 + i - 1, j0 + j - 1) = (1.0_wp - w)*f0(i, j) + w*f1(i, j) end do end subroutine data_input_blend_2d_impl pure subroutine data_input_blend_3d_impl(nxs, nys, nzs, dn1, dn2, dn3, i0, j0, f0, f1, w, dest) integer, intent(in) :: nxs, nys, nzs, dn1, dn2, dn3, i0, j0 real(wp), intent(in) :: f0(nxs, nys, nzs), f1(nxs, nys, nzs) real(wp), intent(in) :: w real(wp), intent(inout) :: dest(dn1, dn2, dn3) integer :: i, j, k do concurrent(k=1:nzs, j=1:nys, i=1:nxs) dest(i0 + i - 1, j0 + j - 1, k) = (1.0_wp - w)*f0(i, j, k) + w*f1(i, j, k) end do end subroutine data_input_blend_3d_impl ! ====================================================================== ! Lifecycle ! ====================================================================== subroutine ocean_data_input_init(this, cfg) !! Allocate the field registry. `cfg` optional so `ocean_state_init` !! (which many tests call directly with no `config_t` in hand) can !! construct a default-sized (16 slots, quiet) reader identically to !! the scaffold it replaces; `ocean_state_init_from_config` re-calls !! this with the real `cfg%ocean%data` once `cfg` is available — !! safe because nothing is registered between the two calls. class(ocean_data_input_t), intent(inout) :: this type(ocean_data_config_t), intent(in), optional :: cfg integer :: max_fields logical :: verbose max_fields = 16 verbose = .false. if (present(cfg)) then max_fields = cfg%max_fields verbose = cfg%verbose end if call this%destroy() allocate (this%fields(max(max_fields, 1))) this%nfields = 0 this%nfields_max = max_fields this%verbose = verbose this%is_init = .true. end subroutine ocean_data_input_init subroutine ocean_data_input_destroy(this) class(ocean_data_input_t), intent(inout) :: this integer :: i if (allocated(this%fields)) then do i = 1, this%nfields if (this%fields(i)%ncid >= 0) call nc_close(this%fields(i)%ncid) end do deallocate (this%fields) end if this%nfields = 0 this%nfields_max = 0 this%is_init = .false. call data_input_workspace_cleanup() end subroutine ocean_data_input_destroy subroutine ocean_data_input_enter_data(this) !! Type-bound wrapper — delegates to the non-polymorphic impl (the !! AMD libomptarget cross-slot-overlap fix; see !! `rdb_ocean_surface_stress.F90`). class(ocean_data_input_t), intent(inout) :: this select type (this) type is (ocean_data_input_t) call ocean_data_input_enter_data_impl(this) end select end subroutine ocean_data_input_enter_data subroutine ocean_data_input_enter_data_impl(this) type(ocean_data_input_t), intent(inout) :: this integer :: i if (.not. this%is_init) return do i = 1, this%nfields if (.not. allocated(this%fields(i)%f0)) cycle !$acc enter data copyin(this%fields(i)%f0, this%fields(i)%f1) ! Rule (2): create-mapped arrays do NOT carry the host values ! read at registration time — push them explicitly. !$acc update device(this%fields(i)%f0, this%fields(i)%f1) end do end subroutine ocean_data_input_enter_data_impl subroutine ocean_data_input_exit_data(this) class(ocean_data_input_t), intent(inout) :: this select type (this) type is (ocean_data_input_t) call ocean_data_input_exit_data_impl(this) end select end subroutine ocean_data_input_exit_data subroutine ocean_data_input_exit_data_impl(this) type(ocean_data_input_t), intent(inout) :: this integer :: i if (.not. this%is_init) return do i = 1, this%nfields if (.not. allocated(this%fields(i)%f0)) cycle !$acc exit data delete(this%fields(i)%f0, this%fields(i)%f1) end do end subroutine ocean_data_input_exit_data_impl pure function ocean_data_input_bytes(this) result(nbytes) !! Counted allocatable footprint (0 when unallocated). Summed over !! every registered field's `f0`/`f1` — the `t_axis` is small !! (`nt` reals) and intentionally excluded, matching the house !! convention of counting device-resident footprint only. class(ocean_data_input_t), intent(in) :: this integer(int64) :: nbytes integer :: i nbytes = 0_int64 if (.not. allocated(this%fields)) return do i = 1, this%nfields nbytes = nbytes + arr_bytes(this%fields(i)%f0) + arr_bytes(this%fields(i)%f1) end do end function ocean_data_input_bytes ! ====================================================================== ! Registration ! ====================================================================== subroutine ocean_data_input_register_2d(this, file, var, grid, dest_n1, dest_n2, & dest_i0, dest_j0, id, time_mode, cycle_period, & t_offset, t_scale, scale, add_offset, oor, & nx_extra, ny_extra, ierr) !! Register a 2-D time-varying field `f(x, y, t)`. Opens the file !! now and keeps the handle for the run. See the module docstring !! for the destination-offset contract. !! !! `nx_extra`/`ny_extra` (default 0) widen this rank's slab by that !! many values in x/y. They exist for C-grid FACE fields: an !! x-face array spans `nx_phys + 1` faces, not `nx_phys`, and the !! trailing face is OWNED by this rank (the west/south-owns-the-seam !! rule), so no exchange can supply it — it has to come from the !! file. A cell-centred field leaves both at 0. class(ocean_data_input_t), intent(inout) :: this character(len=*), intent(in) :: file, var type(hgrid_t), intent(in) :: grid integer, intent(in) :: dest_n1, dest_n2, dest_i0, dest_j0 integer, intent(out) :: id character(len=*), intent(in), optional :: time_mode real(wp), intent(in), optional :: cycle_period, t_offset, t_scale, scale, add_offset integer, intent(in), optional :: oor integer, intent(in), optional :: nx_extra, ny_extra integer, intent(out), optional :: ierr !! Non-zero (`OCEAN_STATUS_ERR_SETUP`/`OCEAN_STATUS_ERR_IO`) on a !! registry/file/dimension failure when present; absent behaves !! as today (`error stop`). integer :: nxe, nye nxe = 0 if (present(nx_extra)) nxe = nx_extra nye = 0 if (present(ny_extra)) nye = ny_extra call register_common(this, file, var, .false., & grid%i_offset_global + 1, grid%j_offset_global + 1, & grid%nx_phys + nxe, grid%ny_phys + nye, 1, & dest_i0, dest_j0, dest_n1, dest_n2, 1, & time_mode, cycle_period, t_offset, t_scale, scale, add_offset, oor, id, & ierr=ierr) end subroutine ocean_data_input_register_2d subroutine ocean_data_input_register_3d(this, file, var, grid, nz_src, dest_n1, dest_n2, & dest_n3, dest_i0, dest_j0, id, time_mode, & cycle_period, t_offset, t_scale, scale, add_offset, oor, & ierr) !! Register a 3-D time-varying field `f(x, y, z, t)`. `nz_src` is !! the source z-level count (the consumer's own concern — the !! reader validates it against the file's z dim and does NOT flip !! or remap k; see the module docstring). class(ocean_data_input_t), intent(inout) :: this character(len=*), intent(in) :: file, var type(hgrid_t), intent(in) :: grid integer, intent(in) :: nz_src, dest_n1, dest_n2, dest_n3, dest_i0, dest_j0 integer, intent(out) :: id character(len=*), intent(in), optional :: time_mode real(wp), intent(in), optional :: cycle_period, t_offset, t_scale, scale, add_offset integer, intent(in), optional :: oor integer, intent(out), optional :: ierr call register_common(this, file, var, .true., & grid%i_offset_global + 1, grid%j_offset_global + 1, & grid%nx_phys, grid%ny_phys, nz_src, & dest_i0, dest_j0, dest_n1, dest_n2, dest_n3, & time_mode, cycle_period, t_offset, t_scale, scale, add_offset, oor, id, & ierr=ierr) end subroutine ocean_data_input_register_3d subroutine ocean_data_input_register_segment_2d(this, file, var, grid, edge, & dest_n1, dest_n2, id, time_mode, & cycle_period, t_offset, t_scale, & scale, add_offset, oor, ierr) !! Register a 2-D OBC-segment field. The degenerate horizontal !! axis (x for west/east, y for south/north) reads as `start=1, !! count=1`; the along-edge axis slices from the global index !! range exactly as `register_2d` does. `dest_i0`/`dest_j0` are !! implied by `edge` (degenerate axis -> index 1; along-edge axis !! -> `grid%nghost + 1`) — not arguments (see module docstring). class(ocean_data_input_t), intent(inout) :: this character(len=*), intent(in) :: file, var type(hgrid_t), intent(in) :: grid integer, intent(in) :: edge, dest_n1, dest_n2 integer, intent(out) :: id character(len=*), intent(in), optional :: time_mode real(wp), intent(in), optional :: cycle_period, t_offset, t_scale, scale, add_offset integer, intent(in), optional :: oor integer, intent(out), optional :: ierr integer :: i0, j0, nx, ny, di0, dj0 call segment_geometry(grid, edge, i0, j0, nx, ny, di0, dj0, ierr=ierr) if (present(ierr)) then if (ierr /= OCEAN_STATUS_OK) return end if call register_common(this, file, var, .false., i0, j0, nx, ny, 1, & di0, dj0, dest_n1, dest_n2, 1, & time_mode, cycle_period, t_offset, t_scale, scale, add_offset, oor, id, & edge=edge, ierr=ierr) end subroutine ocean_data_input_register_segment_2d subroutine ocean_data_input_register_segment_3d(this, file, var, grid, edge, nz_src, & dest_n1, dest_n2, dest_n3, id, time_mode, & cycle_period, t_offset, t_scale, & scale, add_offset, oor, ierr) !! 3-D twin of `register_segment_2d` (PR-22 OBC segment files). class(ocean_data_input_t), intent(inout) :: this character(len=*), intent(in) :: file, var type(hgrid_t), intent(in) :: grid integer, intent(in) :: edge, nz_src, dest_n1, dest_n2, dest_n3 integer, intent(out) :: id character(len=*), intent(in), optional :: time_mode real(wp), intent(in), optional :: cycle_period, t_offset, t_scale, scale, add_offset integer, intent(in), optional :: oor integer, intent(out), optional :: ierr integer :: i0, j0, nx, ny, di0, dj0 call segment_geometry(grid, edge, i0, j0, nx, ny, di0, dj0, ierr=ierr) if (present(ierr)) then if (ierr /= OCEAN_STATUS_OK) return end if call register_common(this, file, var, .true., i0, j0, nx, ny, nz_src, & di0, dj0, dest_n1, dest_n2, dest_n3, & time_mode, cycle_period, t_offset, t_scale, scale, add_offset, oor, id, & edge=edge, ierr=ierr) end subroutine ocean_data_input_register_segment_3d subroutine segment_geometry(grid, edge, i0, j0, nx, ny, dest_i0, dest_j0, ierr) !! Degenerate-axis + along-edge slab geometry for an OBC-segment !! registration. `i0`/`j0`/`nx`/`ny` are FILE-side (start, count); !! `dest_i0`/`dest_j0` are where the slab lands in the consumer's !! staging array. type(hgrid_t), intent(in) :: grid integer, intent(in) :: edge integer, intent(out) :: i0, j0, nx, ny, dest_i0, dest_j0 integer, intent(out), optional :: ierr !! Non-zero (`OCEAN_STATUS_ERR_SETUP`) on an unrecognised `edge` !! tag when present; absent behaves as today (`error stop`). if (present(ierr)) ierr = OCEAN_STATUS_OK select case (edge) case (DATA_EDGE_WEST, DATA_EDGE_EAST) i0 = 1 nx = 1 j0 = grid%j_offset_global + 1 ny = grid%ny_phys dest_i0 = 1 dest_j0 = grid%nghost + 1 case (DATA_EDGE_SOUTH, DATA_EDGE_NORTH) i0 = grid%i_offset_global + 1 nx = grid%nx_phys j0 = 1 ny = 1 dest_i0 = grid%nghost + 1 dest_j0 = 1 case default call fail("ocean_data_input: unknown edge tag "//to_string(edge)// & " (expected DATA_EDGE_WEST/EAST/SOUTH/NORTH)", ierr, OCEAN_STATUS_ERR_SETUP) return end select end subroutine segment_geometry subroutine register_common(this, file, var, is_3d, i0, j0, nx, ny, nz, & dest_i0, dest_j0, dest_n1, dest_n2, dest_n3, & time_mode, cycle_period, t_offset, t_scale, scale, add_offset, & oor, id, edge, ierr) !! Shared registration body for register_2d/_3d/_segment_2d/_segment_3d. !! `(i0, j0, nx, ny, nz)` are FILE-side start/count (already resolved !! by the caller — plain global-offset slicing for the base !! variants, degenerate-axis geometry for the segment variants). class(ocean_data_input_t), intent(inout) :: this character(len=*), intent(in) :: file, var logical, intent(in) :: is_3d integer, intent(in) :: i0, j0, nx, ny, nz integer, intent(in) :: dest_i0, dest_j0, dest_n1, dest_n2, dest_n3 character(len=*), intent(in), optional :: time_mode real(wp), intent(in), optional :: cycle_period, t_offset, t_scale, scale, add_offset integer, intent(in), optional :: oor integer, intent(out) :: id integer, intent(in), optional :: edge integer, intent(out), optional :: ierr !! Non-zero (`OCEAN_STATUS_ERR_SETUP`/`OCEAN_STATUS_ERR_IO`) on a !! registry/file/dimension failure when present; absent behaves !! as today (`error stop`). integer :: ncid, varid, tvarid, tdimid, status integer :: var_ndims integer :: var_dimids(4) integer :: expect_nd, d1, d2, dz, dt, nt character(len=:), allocatable :: units_str logical :: units_ok real(wp) :: t_scale_eff integer :: k integer :: local_ierr if (present(ierr)) ierr = OCEAN_STATUS_OK if (this%nfields >= this%nfields_max) then call fail("ocean_data_input: registry full (max_fields = "// & to_string(this%nfields_max)//"); raise &ocean_data_nml max_fields", & ierr, OCEAN_STATUS_ERR_SETUP) return end if call nc_open_read(trim(file), ncid, ierr=local_ierr) if (.not. reg_io_ok(local_ierr, ierr)) return call nc_check(nf90_inq_varid(ncid, trim(var), varid), & "ocean_data_input: finding variable '"//trim(var)//"' in "//trim(file), local_ierr) if (.not. reg_io_ok(local_ierr, ierr, ncid)) return call nc_check(nf90_inquire_variable(ncid, varid, ndims=var_ndims, & dimids=var_dimids(1:merge(4, 3, is_3d))), & "ocean_data_input: querying variable '"//trim(var)//"'", local_ierr) if (.not. reg_io_ok(local_ierr, ierr, ncid)) return expect_nd = merge(4, 3, is_3d) if (var_ndims /= expect_nd) then call nc_close(ncid) call fail("ocean_data_input: '"//trim(var)//"' in "//trim(file)// & " has "//to_string(var_ndims)//" dims; expected "// & to_string(expect_nd)//" (x,y[,z],t) — Fortran storage order required, "// & "no in-core transpose (regrid/reorder the file)", ierr, OCEAN_STATUS_ERR_IO) return end if call nc_check(nf90_inquire_dimension(ncid, var_dimids(1), len=d1), & "ocean_data_input: dim 1", local_ierr) if (.not. reg_io_ok(local_ierr, ierr, ncid)) return call nc_check(nf90_inquire_dimension(ncid, var_dimids(2), len=d2), & "ocean_data_input: dim 2", local_ierr) if (.not. reg_io_ok(local_ierr, ierr, ncid)) return if (is_3d) then call nc_check(nf90_inquire_dimension(ncid, var_dimids(3), len=dz), & "ocean_data_input: dim 3 (z)", local_ierr) if (.not. reg_io_ok(local_ierr, ierr, ncid)) return call nc_check(nf90_inquire_dimension(ncid, var_dimids(4), len=dt), & "ocean_data_input: dim 4 (t)", local_ierr) if (.not. reg_io_ok(local_ierr, ierr, ncid)) return else dz = 1 call nc_check(nf90_inquire_dimension(ncid, var_dimids(3), len=dt), & "ocean_data_input: dim 3 (t)", local_ierr) if (.not. reg_io_ok(local_ierr, ierr, ncid)) return end if if (present(edge)) then ! Degenerate-axis fail-loud check (§5.7.1): the horizontal ! extent this registration claims is 1 (a segment file) must ! actually BE 1 in the file. if (nx == 1 .and. d1 /= 1) then call nc_close(ncid) call fail("ocean_data_input: segment '"//trim(var)//"' edge "// & to_string(edge)//" expects a degenerate x axis (extent 1); "// & "file has extent "//to_string(d1), ierr, OCEAN_STATUS_ERR_IO) return end if if (ny == 1 .and. d2 /= 1) then call nc_close(ncid) call fail("ocean_data_input: segment '"//trim(var)//"' edge "// & to_string(edge)//" expects a degenerate y axis (extent 1); "// & "file has extent "//to_string(d2), ierr, OCEAN_STATUS_ERR_IO) return end if end if if (.not. dims_geometry_ok(var_ndims, expect_nd, d1, d2, i0 + nx - 1, j0 + ny - 1)) then call nc_close(ncid) call fail("ocean_data_input: '"//trim(var)//"' in "//trim(file)// & " horizontal dims too small for this rank's slab: file ("// & to_string(d1)//","//to_string(d2)//"), need offset+extent ("// & to_string(i0 + nx - 1)//","//to_string(j0 + ny - 1)//")", ierr, OCEAN_STATUS_ERR_IO) return end if if (is_3d .and. dz /= nz) then call nc_close(ncid) call fail("ocean_data_input: '"//trim(var)//"' in "//trim(file)// & " has "//to_string(dz)//" z-levels; caller expects nz_src = "// & to_string(nz), ierr, OCEAN_STATUS_ERR_IO) return end if nt = dt if (nt < 1) then call nc_close(ncid) call fail("ocean_data_input: '"//trim(var)//"' in "//trim(file)//" has no time records", & ierr, OCEAN_STATUS_ERR_IO) return end if ! Time axis: the time DIMENSION and its coordinate VARIABLE share ! the CF name by convention; try that, then fall back to "time". tdimid = var_dimids(expect_nd) call resolve_time_var(ncid, tdimid, tvarid, status) if (status /= nf90_noerr) then call nc_close(ncid) call fail("ocean_data_input: no time coordinate variable found for '"// & trim(var)//"' in "//trim(file), ierr, OCEAN_STATUS_ERR_IO) return end if id = this%nfields + 1 this%nfields = id associate (fld => this%fields(id)) fld%active = .true. fld%ncid = ncid fld%varid = varid fld%is_3d = is_3d fld%nx = nx fld%ny = ny fld%nz = merge(nz, 1, is_3d) fld%i0 = i0 fld%j0 = j0 fld%dest_i0 = dest_i0 fld%dest_j0 = dest_j0 fld%dest_n1 = dest_n1 fld%dest_n2 = dest_n2 fld%dest_n3 = dest_n3 fld%nt = nt fld%scale = 1.0_wp if (present(scale)) fld%scale = scale fld%add_offset = 0.0_wp if (present(add_offset)) fld%add_offset = add_offset fld%t_offset = 0.0_wp if (present(t_offset)) fld%t_offset = t_offset fld%cycle_period = 0.0_wp if (present(cycle_period)) fld%cycle_period = cycle_period fld%oor = DATA_OOR_ERROR if (present(oor)) fld%oor = oor fld%time_mode = DATA_TIME_LINEAR if (present(time_mode)) fld%time_mode = data_input_time_mode_from_string(time_mode) if (fld%time_mode == DATA_TIME_CYCLIC .and. fld%cycle_period <= 0.0_wp) then fld%active = .false. this%nfields = this%nfields - 1 call nc_close(ncid) call fail("ocean_data_input: '"//trim(var)//"' time_mode='cyclic' requires "// & "cycle_period > 0", ierr, OCEAN_STATUS_ERR_SETUP) return end if allocate (fld%t_axis(nt)) call nc_get_var_1d(ncid, tvarid, fld%t_axis) t_scale_eff = 1.0_wp if (present(t_scale)) then t_scale_eff = t_scale else call nc_get_att_text(ncid, tvarid, "units", units_str, units_ok) if (units_ok) then call data_input_time_scale_from_units(trim(units_str), t_scale_eff, units_ok) end if ! units_ok = .false. (absent/unrecognised attribute) silently ! keeps the seconds-default — the common case for a ! model-native file with no CF units attribute. end if fld%t_axis = fld%t_axis*t_scale_eff do k = 2, nt if (fld%t_axis(k) <= fld%t_axis(k - 1)) then call logger%error("ocean_data_input: '"//trim(var)//"' time axis not "// & "monotonically increasing at record "//to_string(k)//" ("// & to_string(fld%t_axis(k))//" <= "//to_string(fld%t_axis(k - 1))//")") error stop "ocean_data_input: non-monotonic time axis" end if end do allocate (fld%f0(fld%nx, fld%ny, fld%nz)) allocate (fld%f1(fld%nx, fld%ny, fld%nz)) if (fld%time_mode == DATA_TIME_STATIC) then call data_input_read_slab_impl(fld, 1) fld%f1 = fld%f0 fld%rec0 = 1 fld%rec1 = 1 fld%w = 0.0_wp end if end associate end subroutine register_common function reg_io_ok(local_ierr, ierr, ncid) result(ok) !! Translate a raw `nc_check`-style status (0 = ok) from one of the !! `nc_*` reader calls in `register_common` into the caller's `ierr` !! contract: `.true.` on success; on failure, `.false.` with !! `ierr = OCEAN_STATUS_ERR_IO` when `ierr` is present (closing !! `ncid` first, when given, so a mid-registration failure does not !! leak the file handle), or `error stop`s with a generic message !! when `ierr` is absent — the SPECIFIC reason was already logged by !! `nc_check`/`fail` before this returns, so the legacy (no `ierr`) !! log output is unchanged; only the raw `error stop` text is !! generic (same idiom as `rdb_bathymetry::bathy_io_ok`, P0.1 F1/F2). integer, intent(in) :: local_ierr integer, intent(out), optional :: ierr integer, intent(in), optional :: ncid logical :: ok integer :: discard_ierr ok = (local_ierr == 0) if (ok) return if (present(ierr)) then if (present(ncid)) call nc_close(ncid, ierr=discard_ierr) ierr = OCEAN_STATUS_ERR_IO return end if error stop "ocean_data_input: NetCDF registration failure" end function reg_io_ok subroutine resolve_time_var(ncid, tdimid, tvarid, status) !! Look up the time coordinate variable by the CF convention (dim !! name == var name); fall back to a variable literally named !! "time". `status` is `nf90_noerr` iff found. integer, intent(in) :: ncid, tdimid integer, intent(out) :: tvarid, status character(len=64) :: dname status = nf90_inquire_dimension(ncid, tdimid, name=dname) if (status == nf90_noerr) status = nf90_inq_varid(ncid, trim(dname), tvarid) if (status /= nf90_noerr) status = nf90_inq_varid(ncid, "time", tvarid) end subroutine resolve_time_var ! ====================================================================== ! Per-step path ! ====================================================================== subroutine ocean_data_input_update_all(this, t) !! Driver hook: refresh every registered field's bracket for model !! time `t` (reading + pushing a new slab to the device only when !! the bracket actually advances). No-op when `nfields == 0` — !! every shipped namelist today. Does NOT write into any !! consumer array (see module docstring) — call `update_2d/_3d` !! afterwards for that. class(ocean_data_input_t), intent(inout) :: this real(wp), intent(in) :: t integer :: i if (this%nfields == 0) return do i = 1, this%nfields call data_input_refresh_brackets(this, i, t) end do end subroutine ocean_data_input_update_all subroutine data_input_refresh_brackets(this, id, t) !! Host-side registry walk (outer shim) for field `id`. class(ocean_data_input_t), intent(inout) :: this integer, intent(in) :: id real(wp), intent(in) :: t integer :: n0, n1 real(wp) :: w, t_query logical :: oor associate (fld => this%fields(id)) if (fld%time_mode == DATA_TIME_STATIC) then fld%t_last = t return end if t_query = t + fld%t_offset call data_input_locate(fld%t_axis, fld%nt, fld%time_mode, fld%cycle_period, & t_query, n0, n1, w, oor) if (oor) then if (fld%oor == DATA_OOR_ERROR) then call logger%error("ocean_data_input: query time "//to_string(t_query)// & "s is outside the file time axis ["// & to_string(fld%t_axis(1))//", "// & to_string(fld%t_axis(fld%nt))//"]s (field id "// & to_string(id)//"); pass oor=DATA_OOR_CLAMP to clamp instead") error stop "ocean_data_input: query time out of range" end if if (.not. fld%oor_warned) then call logger%warning("ocean_data_input: field id "//to_string(id)// & "clamped at the time-axis boundary (query time "// & to_string(t_query)//"s out of range)") fld%oor_warned = .true. end if end if ! `!$acc update device(...)` is only correct HERE (the runtime ! per-step path, always called after `ocean_state_enter_data` ! has mapped `fld%f0`/`fld%f1`) — never inside ! `data_input_read_slab_impl` itself, which is also called from ! `register_common` at SETUP time for a DATA_TIME_STATIC field, ! before anything is device-mapped (CLAUDE.md gotcha: an ! `update device` on an unmapped array is undefined). if (n0 /= fld%rec0) then call data_input_read_slab_impl(fld, n0, into_f1=.false.) !$acc update device(fld%f0) end if if (n1 /= fld%rec1) then if (n1 == n0) then fld%f1 = fld%f0 !$acc update device(fld%f1) else call data_input_read_slab_impl(fld, n1, into_f1=.true.) !$acc update device(fld%f1) end if end if if (n0 /= fld%rec0 .or. n1 /= fld%rec1) then if (this%verbose) then call logger%info("ocean_data_input: field id "//to_string(id)// & " bracket advanced to records ("//to_string(n0)//", "// & to_string(n1)//"), w = "//to_string(w)// & ", t = "//to_string(t)) end if end if fld%rec0 = n0 fld%rec1 = n1 fld%w = w fld%t_last = t end associate end subroutine data_input_refresh_brackets subroutine data_input_read_slab_impl(fld, rec, into_f1) !! Host-ONLY NetCDF slab read for file record `rec` -> `fld%f0` !! (default) or `fld%f1` (`into_f1 = .true.`), applying !! `scale`/`add_offset` once at read. Uses the module-level host !! workspace (registration/bracket-advance-time only — never a !! per-step allocation). Deliberately carries NO `!$acc` directive: !! called both before `enter_data` (STATIC field, at registration) !! and after it (LINEAR/CYCLIC bracket advance) — the caller pushes !! to the device itself, only when that is actually correct !! (`data_input_refresh_brackets`). type(data_input_field_t), intent(inout) :: fld integer, intent(in) :: rec logical, intent(in), optional :: into_f1 integer :: start4(4), count4(4) logical :: to_f1 to_f1 = .false. if (present(into_f1)) to_f1 = into_f1 call data_input_workspace_ensure(fld%nx, fld%ny, fld%nz) if (fld%is_3d) then start4 = [fld%i0, fld%j0, 1, rec] count4 = [fld%nx, fld%ny, fld%nz, 1] call nc_get_var_slab_3d(fld%ncid, fld%varid, start4, count4, data_input_ws) else start4(1:3) = [fld%i0, fld%j0, rec] count4(1:3) = [fld%nx, fld%ny, 1] call nc_get_var_slab_3d(fld%ncid, fld%varid, start4(1:3), count4(1:3), data_input_ws) end if data_input_ws = fld%scale*data_input_ws + fld%add_offset if (to_f1) then fld%f1 = data_input_ws(1:fld%nx, 1:fld%ny, 1:fld%nz) else fld%f0 = data_input_ws(1:fld%nx, 1:fld%ny, 1:fld%nz) end if fld%nreads = fld%nreads + 1 end subroutine data_input_read_slab_impl subroutine data_input_workspace_ensure(nx, ny, nz) integer, intent(in) :: nx, ny, nz logical :: need_alloc need_alloc = .false. if (.not. allocated(data_input_ws)) then need_alloc = .true. else if (size(data_input_ws, 1) /= nx .or. size(data_input_ws, 2) /= ny .or. & size(data_input_ws, 3) /= nz) then deallocate (data_input_ws) need_alloc = .true. end if if (need_alloc) allocate (data_input_ws(nx, ny, nz)) end subroutine data_input_workspace_ensure subroutine data_input_workspace_cleanup() if (allocated(data_input_ws)) deallocate (data_input_ws) end subroutine data_input_workspace_cleanup subroutine ocean_data_input_update_2d(this, id, t, n1, n2, dest) !! Blend field `id`'s current bracket into the caller's WHOLE, !! device-mapped `dest(n1, n2)` array at the registration-time !! offset. `t` must equal the value `update_all` most recently !! refreshed this field with (fail-loud ordering check — a !! consumer calling this before `update_all` has run for the !! current step is a real bug, not a silent stale read). Exempt !! for `DATA_TIME_STATIC` fields, whose bracket never changes. class(ocean_data_input_t), intent(in) :: this integer, intent(in) :: id, n1, n2 real(wp), intent(in) :: t real(wp), intent(inout) :: dest(n1, n2) call check_registered(this, id, is_3d=.false., n1=n1, n2=n2) call check_fresh(this, id, t) call data_input_blend_2d_impl(this%fields(id)%nx, this%fields(id)%ny, n1, n2, & this%fields(id)%dest_i0, this%fields(id)%dest_j0, & this%fields(id)%f0, this%fields(id)%f1, & this%fields(id)%w, dest) end subroutine ocean_data_input_update_2d subroutine ocean_data_input_update_3d(this, id, t, n1, n2, n3, dest) !! 3-D twin of `update_2d`. class(ocean_data_input_t), intent(in) :: this integer, intent(in) :: id, n1, n2, n3 real(wp), intent(in) :: t real(wp), intent(inout) :: dest(n1, n2, n3) call check_registered(this, id, is_3d=.true., n1=n1, n2=n2, n3=n3) call check_fresh(this, id, t) call data_input_blend_3d_impl(this%fields(id)%nx, this%fields(id)%ny, this%fields(id)%nz, & n1, n2, n3, & this%fields(id)%dest_i0, this%fields(id)%dest_j0, & this%fields(id)%f0, this%fields(id)%f1, & this%fields(id)%w, dest) end subroutine ocean_data_input_update_3d subroutine ocean_data_input_fill_static_host(this, id, n1, n2, dest) !! One-shot HOST-side fill of a `DATA_TIME_STATIC` field. Copies !! the already-read record-1 slab into `dest` at the !! registration-time offsets. Touches NO device memory — !! callable at setup, before `dest` is mapped. Plain host `do` !! loops (NOT `do concurrent` — on `-stdpar=gpu` a bare `do !! concurrent` is unconditionally offloaded regardless of whether !! `dest` is mapped, which is exactly the silent-stale-write bug !! this routine exists to avoid). Fails loud for any non-static !! field. class(ocean_data_input_t), intent(in) :: this integer, intent(in) :: id, n1, n2 real(wp), intent(inout) :: dest(n1, n2) integer :: i, j call check_registered(this, id, is_3d=.false., n1=n1, n2=n2) if (this%fields(id)%time_mode /= DATA_TIME_STATIC) then call logger%error("ocean_data_input: fill_static_host called on field id "// & to_string(id)//", which is not DATA_TIME_STATIC") error stop "ocean_data_input: fill_static_host on a non-static field" end if do j = 1, this%fields(id)%ny do i = 1, this%fields(id)%nx dest(this%fields(id)%dest_i0 + i - 1, this%fields(id)%dest_j0 + j - 1) = & this%fields(id)%f0(i, j, 1) end do end do end subroutine ocean_data_input_fill_static_host subroutine ocean_data_input_load_static_2d(file, var, grid, n1, n2, dest_i0, dest_j0, & dest, ierr) !! One-shot STATIC 2-D read: register, fill, close. The setup-time !! convenience wrapper over `register_2d(time_mode="static")` + !! `fill_static_host` for a field that is read ONCE and never !! updated — a prescribed, time-constant geometry rather than a !! forcing. `&ocean_cavity_dyn_nml draft_config="file"` is the !! first consumer. !! !! The reader is LOCAL and destroyed on the way out: nothing is left !! on `ocean_state_t%data_input`, which exists for the per-step !! forcing tags and whose registry would otherwise be permanently !! one slot smaller for a field never read again. It also keeps the !! static geometry out of `enter_data` — the destination array is !! the state's own, and IT is what gets mapped. !! !! FILE CONTRACT (inherited from `register_common`, not weakened !! here): `var` must be rank 3 in FORTRAN storage order `(x, y, t)` !! — i.e. `(nTime, ny, nx)` as a C/Python writer or `ncdump` spells !! it — with a time COORDINATE VARIABLE, and its horizontal extents !! must cover this rank's slab. Record 1 is read. There is NO !! horizontal interpolation: the file must already be on the model !! grid, exactly as `bathymetry_file` and `&ocean_zinit_nml file` !! require. !! !! HOST-ONLY: runs at setup, before `dest` is device-mapped. character(len=*), intent(in) :: file, var type(hgrid_t), intent(in) :: grid integer, intent(in) :: n1, n2 !! Shape of `dest` — the FULL ghosted `(nx_total, ny_total)`. integer, intent(in) :: dest_i0, dest_j0 !! Destination index of the first PHYSICAL cell (`nghost + 1`). real(wp), intent(inout) :: dest(n1, n2) integer, intent(out), optional :: ierr !! Non-zero on a registry / file / dimension failure when !! present; absent behaves as the rest of the reader !! (`error stop`). type(ocean_data_input_t) :: reader integer :: id, local_ierr if (present(ierr)) ierr = OCEAN_STATUS_OK call reader%init() if (present(ierr)) then call ocean_data_input_register_2d(reader, file, var, grid, & dest_n1=n1, dest_n2=n2, & dest_i0=dest_i0, dest_j0=dest_j0, & id=id, time_mode="static", ierr=local_ierr) if (local_ierr /= OCEAN_STATUS_OK) then call reader%destroy() ierr = local_ierr return end if else call ocean_data_input_register_2d(reader, file, var, grid, & dest_n1=n1, dest_n2=n2, & dest_i0=dest_i0, dest_j0=dest_j0, & id=id, time_mode="static") end if call ocean_data_input_fill_static_host(reader, id, n1, n2, dest) call reader%destroy() end subroutine ocean_data_input_load_static_2d subroutine ocean_data_input_fill_static_host_3d(this, id, n1, n2, n3, dest) !! 3-D twin of `fill_static_host`. class(ocean_data_input_t), intent(in) :: this integer, intent(in) :: id, n1, n2, n3 real(wp), intent(inout) :: dest(n1, n2, n3) integer :: i, j, k call check_registered(this, id, is_3d=.true., n1=n1, n2=n2, n3=n3) if (this%fields(id)%time_mode /= DATA_TIME_STATIC) then call logger%error("ocean_data_input: fill_static_host_3d called on field id "// & to_string(id)//", which is not DATA_TIME_STATIC") error stop "ocean_data_input: fill_static_host_3d on a non-static field" end if do k = 1, this%fields(id)%nz do j = 1, this%fields(id)%ny do i = 1, this%fields(id)%nx dest(this%fields(id)%dest_i0 + i - 1, this%fields(id)%dest_j0 + j - 1, k) = & this%fields(id)%f0(i, j, k) end do end do end do end subroutine ocean_data_input_fill_static_host_3d subroutine check_registered(this, id, is_3d, n1, n2, n3) !! Fail-loud guard shared by every per-step/fill accessor: `id` !! must name an active field of the right rank, and the caller's !! `dest` shape must match what was declared at registration. class(ocean_data_input_t), intent(in) :: this integer, intent(in) :: id, n1, n2 logical, intent(in) :: is_3d integer, intent(in), optional :: n3 integer :: dn3 if (id < 1 .or. id > this%nfields) then call logger%error("ocean_data_input: invalid field id "//to_string(id)) error stop "ocean_data_input: invalid field id" end if if (.not. this%fields(id)%active) then call logger%error("ocean_data_input: field id "//to_string(id)//" is not active") error stop "ocean_data_input: field id not active" end if if (this%fields(id)%is_3d .neqv. is_3d) then call logger%error("ocean_data_input: field id "//to_string(id)// & " rank mismatch (registered "// & merge("3-D", "2-D", this%fields(id)%is_3d)//", called as "// & merge("3-D", "2-D", is_3d)//")") error stop "ocean_data_input: field rank mismatch" end if dn3 = 1 if (present(n3)) dn3 = n3 if (n1 /= this%fields(id)%dest_n1 .or. n2 /= this%fields(id)%dest_n2 .or. & dn3 /= this%fields(id)%dest_n3) then call logger%error("ocean_data_input: field id "//to_string(id)// & " dest shape ("//to_string(n1)//","//to_string(n2)//","// & to_string(dn3)//") does not match the registered shape ("// & to_string(this%fields(id)%dest_n1)//","// & to_string(this%fields(id)%dest_n2)//","// & to_string(this%fields(id)%dest_n3)//")") error stop "ocean_data_input: dest shape mismatch" end if end subroutine check_registered subroutine check_fresh(this, id, t) !! Ordering guard: `update_2d/_3d` must be called with the same !! `t` `update_all` most recently refreshed this field's bracket !! with (skipped for STATIC — its bracket never changes, so !! "freshness" is meaningless). Catches a consumer calling !! `update_2d/_3d` before `update_all` has run this step, which !! would otherwise silently blend a stale bracket. class(ocean_data_input_t), intent(in) :: this integer, intent(in) :: id real(wp), intent(in) :: t if (this%fields(id)%time_mode == DATA_TIME_STATIC) return if (t /= this%fields(id)%t_last) then call logger%error("ocean_data_input: field id "//to_string(id)// & " queried at t = "//to_string(t)// & "s but its bracket was last refreshed at t = "// & to_string(this%fields(id)%t_last)// & "s — call ocean_data_input_update_all(this, t) first") error stop "ocean_data_input: stale bracket (update_all not called for this t)" end if end subroutine check_fresh end module rdb_ocean_data_input