load_bathymetry_into_array Subroutine

   subroutine load_bathymetry_into_array(filename, b, grid, ierr)
      !! Load bathymetry from a NetCDF file into a target 2D array
      !! shaped `(grid%nx_total, grid%ny_total)`.
      !!
      !! The file must contain:
      !!   - Dimensions: x (nx_phys), y (ny_phys)
      !!   - Variable: b(x, y) or elevation(x, y) or depth(x, y)
      !!
      !! Handles NetCDF C/Fortran dimension reversal: files written by
      !! Python/C store b(x,y) in C order, which Fortran reads as b(y,x).
      !! If the direct dimensions don't match but the transposed ones do,
      !! we read transposed and copy correctly.
      !!
      !! Bathymetry lands in the interior; ghost cells are filled by
      !! constant extrapolation from the nearest interior cell via
      !! `fill_bathymetry_ghosts_array`.
      !!
      !! The file holds the WHOLE grid (`grid%nx_global x grid%ny_global`).
      !! A decomposed tile reads only the full-width band of rows its
      !! storage covers (physical rows +- nghost, clipped to the grid) with a
      !! start/count window, extrapolates that band's ghosts exactly as the
      !! whole grid would (at a GLOBAL edge the band edge is the grid edge;
      !! elsewhere the band's rows ARE the tile's ghost rows, read from the
      !! file), and keeps its own columns: every tile is bit-identical to
      !! its slice of the single-rank array, seam ghosts included.
      character(len=*), intent(in) :: filename
      real(wp), intent(inout) :: b(:, :)
      type(hgrid_t), intent(in) :: grid
      integer, intent(out), optional :: ierr
         !! Non-zero on a missing/unreadable file or a dimension mismatch
         !! when present; absent behaves as today (`error stop`).

      integer :: ncid, varid, file_nx, file_ny
      integer :: ng, i, j, ni, nj, jb0, jb1, nb, jo_b
      type(hgrid_t) :: gband
      real(wp), allocatable :: bband(:, :)
      integer :: var_dimids(2)
      integer :: var_ndims, dim1_len, dim2_len
      integer :: local_ierr
      character(len=64) :: dim1_name, dim2_name
      logical :: needs_transpose
      real(wp), allocatable :: b_interior(:, :)

      call logger%info("Loading bathymetry from: "//trim(filename))

      call nc_open_read(filename, ncid, ierr=local_ierr)
      if (.not. bathy_io_ok(local_ierr, ierr)) return

      ! Read grid dimensions (by name — these are the logical sizes)
      call nc_get_dim_len(ncid, "x", file_nx, ierr=local_ierr)
      if (.not. bathy_io_ok(local_ierr, ierr, ncid)) return
      call nc_get_dim_len(ncid, "y", file_ny, ierr=local_ierr)
      if (.not. bathy_io_ok(local_ierr, ierr, ncid)) return

      ! Try variable names in order: b, elevation, depth
      call try_get_bathymetry_var(ncid, varid, ierr=local_ierr)
      if (.not. bathy_io_ok(local_ierr, ierr, ncid)) return

      ! Query the variable's actual dimension order in the file.
      ! Files written by Python/C store b(x,y) in C-order. Fortran's
      ! NetCDF library reverses this, so the variable appears as b(y,x)
      ! in Fortran. We detect this by checking the first dimension name.
      call nc_check(nf90_inquire_variable(ncid, varid, ndims=var_ndims, &
                                          dimids=var_dimids), &
                    "querying bathymetry variable", local_ierr)
      if (.not. bathy_io_ok(local_ierr, ierr, ncid)) return
      call nc_check(nf90_inquire_dimension(ncid, var_dimids(1), &
                                           name=dim1_name, len=dim1_len), &
                    "querying dim 1", local_ierr)
      if (.not. bathy_io_ok(local_ierr, ierr, ncid)) return
      call nc_check(nf90_inquire_dimension(ncid, var_dimids(2), &
                                           name=dim2_name, len=dim2_len), &
                    "querying dim 2", local_ierr)
      if (.not. bathy_io_ok(local_ierr, ierr, ncid)) return

      call logger%info("Bathymetry variable: "//trim(dim1_name)//"="// &
                       to_string(dim1_len)//" x "//trim(dim2_name)//"="// &
                       to_string(dim2_len)//" (Fortran order)")

      ! In Fortran, dim1 is the first array index (varies fastest).
      ! If dim1 is "x" and matches nx, no transpose needed.
      ! If dim1 is "y" (C-order file reversed), we need to transpose.
      needs_transpose = (trim(dim1_name) == "y")

      ! The file describes the WHOLE grid.
      ni = grid%nx_global
      nj = grid%ny_global

      ! Validate dimensions
      if (needs_transpose) then
         if (dim1_len /= nj .or. dim2_len /= ni) then
            call logger%error("Bathymetry grid mismatch: file has "// &
                              to_string(dim2_len)//" x "//to_string(dim1_len)// &
                              " but simulation expects "// &
                              to_string(ni)//" x "// &
                              to_string(nj))
            call nc_close(ncid)
            if (present(ierr)) then
               ierr = OCEAN_STATUS_ERR_IO
               return
            end if
            error stop "Bathymetry grid mismatch"
         end if
      else
         if (dim1_len /= ni .or. dim2_len /= nj) then
            call logger%error("Bathymetry grid mismatch: file has "// &
                              to_string(dim1_len)//" x "//to_string(dim2_len)// &
                              " but simulation expects "// &
                              to_string(ni)//" x "// &
                              to_string(nj))
            call nc_close(ncid)
            if (present(ierr)) then
               ierr = OCEAN_STATUS_ERR_IO
               return
            end if
            error stop "Bathymetry grid mismatch"
         end if
      end if

      ng = grid%nghost
      ! Rows of the whole grid this tile's storage covers (all of them on an
      ! undecomposed grid), full width.
      jb0 = max(1, grid%j_offset_global + 1 - ng)
      jb1 = min(nj, grid%j_offset_global + grid%ny_phys + ng)
      nb = jb1 - jb0 + 1

      ! Read the band into a temporary matching the Fortran storage order.
      if (needs_transpose) then
         allocate (b_interior(nb, ni))
         call nc_get_var_slab_2d(ncid, varid, [jb0, 1], [nb, ni], b_interior, ierr=local_ierr)
      else
         allocate (b_interior(ni, nb))
         call nc_get_var_slab_2d(ncid, varid, [1, jb0], [ni, nb], b_interior, ierr=local_ierr)
      end if
      if (.not. bathy_io_ok(local_ierr, ierr, ncid)) return
      call nc_close(ncid)

      ! Place the band in a whole-width, band-height array with ghosts and
      ! extrapolate ITS ghosts -- exactly the undecomposed construction when
      ! the band is the whole grid.
      call gband%init(ni, nb, ng, grid%dx, grid%dy)
      allocate (bband(gband%nx_total, gband%ny_total))
      bband = 0.0_wp
      if (needs_transpose) then
         ! b_interior is (y, x) in Fortran — transpose to (x, y)
         do j = 1, nb
            do i = 1, ni
               bband(ng + i, ng + j) = b_interior(j, i)
            end do
         end do
      else
         bband(ng + 1:ng + ni, ng + 1:ng + nb) = b_interior
      end if
      deallocate (b_interior)

      ! Fill ghost cells by constant extrapolation from nearest interior cell.
      ! Without this, ghost cells retain b=0 which creates artificial cliffs
      ! against real bathymetry (e.g. b=-50m interior vs b=0 ghost), generating
      ! extreme velocities and tiny CFL timesteps.
      call fill_bathymetry_ghosts_array(bband, gband)

      ! Cut this tile's storage window (ghosts included) out of the band.
      jo_b = grid%j_offset_global - (jb0 - 1)
      do j = 1, grid%ny_total
         do i = 1, grid%nx_total
            b(i, j) = bband(i + grid%i_offset_global, j + jo_b)
         end do
      end do
      deallocate (bband)

      call logger%info("Bathymetry loaded: min = "// &
                       to_string(minval(b(ng + 1:ng + grid%nx_phys, &
                                          ng + 1:ng + grid%ny_phys)))// &
                       " m, max = "// &
                       to_string(maxval(b(ng + 1:ng + grid%nx_phys, &
                                          ng + 1:ng + grid%ny_phys)))// &
                       " m")

      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine load_bathymetry_into_array