set_bathymetry_neverworld2 Subroutine

public subroutine set_bathymetry_neverworld2(b, grid, max_depth, nl_continent_amp, nl_roughness_amp, min_depth)

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules.

Fill b(:,:) with the Neverworld2 idealized-basin bathymetry (Marques et al. 2022, GMD; MOM6-inspired). A single Pangaea-style basin spanning a 60°×140° spherical sector with a re-entrant (periodic-x) southern channel — the Drake-Passage analog. Depth is built in normalized coordinates

x = (i_phys − 0.5)/nx_phys ∈ [0,1],   y = (j_phys − 0.5)/ny_phys ∈ [0,1]

which equal MOM6’s (lon − west)/len_lon and (lat − south)/len_lat on a uniform grid, so no metrics access is needed. The fractional depth is

D_frac(x,y) = 1
   − 1.1·spike(y−1, 0.12)              ! great northern wall
   − 1.1·spike(y,   0.12)              ! Antarctica (south wall)
   − A_c·[ continents + ridges ]        ! A_c = nl_continent_amp
   − A_r·cos(14πx)·sin(14πy)            ! A_r = nl_roughness_amp
   − A_r·cos(20πx)·cos(20πy)

clamped D_frac = max(D_frac, 0), then D = D_frac·max_depth. The continent/ridge block carves two meridional barriers (S-America at the x-edges, Africa mid-basin) that wall the northern gyres, plus the Drake/Scotia ridge system that sets the channel sill; A_c=0 gives an aquaplanet with the southern channel only (useful for bring-up).

Sign convention: b is bottom depth positive-down (matching the spoon + flat branches). Land emerges where D_frac→0 (the wet/dry mask is seeded downstream from b ≥ LAND_DEPTH_THRESHOLD).

GHOST-ROW FOOTGUN (load-bearing): the loop runs over the FULL array incl. ghost rows. Leaving ghosts at the alloc-time zero gives h_layer = 0 there, sending the EOS into its vanishing-layer fallback (rho_layer = rho_0) → a spurious density jump at wall-adjacent faces → ~12-h e-fold blowup (the same trap the spoon/seamount setters guard). Host only — call enter_data afterwards.

MPI: the normalised coordinates divide by the GLOBAL extents grid%nx_global / grid%ny_global, and local cell indices are mapped to global physical indices with grid%i_offset_global / grid%j_offset_global, so x/y sweep [0,1] across the WHOLE basin. Normalised by the LOCAL extents instead, every rank would rebuild the entire Pangaea basin — walls, continents and all — inside its own tile. On a single rank the offsets are 0 and global == local, so the fill is byte-identical.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: b(:,:)
type(hgrid_t), intent(in) :: grid
real(kind=wp), intent(in) :: max_depth
real(kind=wp), intent(in) :: nl_continent_amp
real(kind=wp), intent(in) :: nl_roughness_amp
real(kind=wp), intent(in) :: min_depth

Floor on the resulting depth (m). The fractional-depth formula produces b in [0, ~1.1*max_depth]; the 1.1*spike walls and the continent terms drive b to 0 (land) at the basin edges. The ocean C-grid dyn-core does not yet carry a robust wet/dry path, and true zero-depth land + the resulting sub-metre surface layers blow up within ~12 h. Flooring every cell to min_depth (MOM6’s MINIMUM_DEPTH approach) turns the continents into shallow shelves that steer — rather than hard-block — the flow, keeping the basin all-wet and stable. Lower values give stronger topographic steering at the cost of thinner layers (min_depth/nz); true land barriers wait on the wet/dry-plumbing follow-up.

Calls

proc~~set_bathymetry_neverworld2~~CallsGraph proc~set_bathymetry_neverworld2 set_bathymetry_neverworld2 proc~nw2_cosbell nw2_cosbell proc~set_bathymetry_neverworld2->proc~nw2_cosbell proc~nw2_spike nw2_spike proc~set_bathymetry_neverworld2->proc~nw2_spike

Called by

proc~~set_bathymetry_neverworld2~~CalledByGraph proc~set_bathymetry_neverworld2 set_bathymetry_neverworld2 proc~ocean_state_seed_from_cfg ocean_state_seed_from_cfg proc~ocean_state_seed_from_cfg->proc~set_bathymetry_neverworld2 proc~engine_setup engine_setup proc~engine_setup->proc~ocean_state_seed_from_cfg proc~complete_ocean_create complete_ocean_create proc~complete_ocean_create->proc~engine_setup proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_setup proc~driver_validate driver_validate proc~driver_validate->proc~engine_setup proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean proc~rdb_ocean_create_finalize rdb_ocean_create_finalize proc~rdb_ocean_create_finalize->proc~complete_ocean_create proc~rdb_ocean_create_from_string rdb_ocean_create_from_string proc~rdb_ocean_create_from_string->proc~complete_ocean_create

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private, parameter :: PI = 4.0_wp*atan(1.0_wp)
real(kind=wp), private :: d_frac
integer, private :: i
integer, private :: i_phys
integer, private :: ioff
integer, private :: j
integer, private :: j_phys
integer, private :: joff
integer, private :: ng
real(kind=wp), private :: nxp
real(kind=wp), private :: nyp
real(kind=wp), private :: x
real(kind=wp), private :: y

Source Code

   subroutine set_bathymetry_neverworld2(b, grid, max_depth, nl_continent_amp, nl_roughness_amp, min_depth)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !!
      !! Fill `b(:,:)` with the **Neverworld2** idealized-basin bathymetry
      !! (Marques et al. 2022, GMD; MOM6-inspired).  A single Pangaea-style
      !! basin spanning a 60°×140° spherical sector with a re-entrant
      !! (periodic-x) southern channel — the Drake-Passage analog.  Depth is
      !! built in **normalized coordinates**
      !!
      !!     x = (i_phys − 0.5)/nx_phys ∈ [0,1],   y = (j_phys − 0.5)/ny_phys ∈ [0,1]
      !!
      !! which equal MOM6's `(lon − west)/len_lon` and `(lat − south)/len_lat`
      !! on a uniform grid, so no metrics access is needed.  The fractional
      !! depth is
      !!
      !!     D_frac(x,y) = 1
      !!        − 1.1·spike(y−1, 0.12)              ! great northern wall
      !!        − 1.1·spike(y,   0.12)              ! Antarctica (south wall)
      !!        − A_c·[ continents + ridges ]        ! A_c = nl_continent_amp
      !!        − A_r·cos(14πx)·sin(14πy)            ! A_r = nl_roughness_amp
      !!        − A_r·cos(20πx)·cos(20πy)
      !!
      !! clamped `D_frac = max(D_frac, 0)`, then `D = D_frac·max_depth`.  The
      !! continent/ridge block carves two meridional barriers (S-America at the
      !! x-edges, Africa mid-basin) that wall the northern gyres, plus the
      !! Drake/Scotia ridge system that sets the channel sill; `A_c=0` gives an
      !! aquaplanet with the southern channel only (useful for bring-up).
      !!
      !! Sign convention: `b` is bottom depth positive-down (matching the spoon
      !! + flat branches).  Land emerges where D_frac→0 (the wet/dry mask is
      !! seeded downstream from `b ≥ LAND_DEPTH_THRESHOLD`).
      !!
      !! GHOST-ROW FOOTGUN (load-bearing): the loop runs over the FULL array
      !! incl. ghost rows.  Leaving ghosts at the alloc-time zero gives
      !! `h_layer = 0` there, sending the EOS into its vanishing-layer fallback
      !! (`rho_layer = rho_0`) → a spurious density jump at wall-adjacent faces
      !! → ~12-h e-fold blowup (the same trap the spoon/seamount setters guard).
      !! Host only — call `enter_data` afterwards.
      !!
      !! MPI: the normalised coordinates divide by the GLOBAL extents
      !! `grid%nx_global` / `grid%ny_global`, and local cell indices are
      !! mapped to global physical indices with `grid%i_offset_global` /
      !! `grid%j_offset_global`, so `x`/`y` sweep [0,1] across the WHOLE
      !! basin.  Normalised by the LOCAL extents instead, every rank would
      !! rebuild the entire Pangaea basin — walls, continents and all —
      !! inside its own tile.  On a single rank the offsets are 0 and
      !! global == local, so the fill is byte-identical.
      real(wp), intent(inout) :: b(:, :)
      type(hgrid_t), intent(in) :: grid
      real(wp), intent(in) :: max_depth, nl_continent_amp, nl_roughness_amp
      real(wp), intent(in) :: min_depth
         !! Floor on the resulting depth (m).  The fractional-depth formula
         !! produces b in [0, ~1.1*max_depth]; the `1.1*spike` walls and the
         !! continent terms drive b to 0 (land) at the basin edges.  The ocean
         !! C-grid dyn-core does not yet carry a robust wet/dry path, and true
         !! zero-depth land + the resulting sub-metre surface layers blow up
         !! within ~12 h.  Flooring every cell to `min_depth` (MOM6's
         !! MINIMUM_DEPTH approach) turns the continents into shallow shelves
         !! that steer — rather than hard-block — the flow, keeping the basin
         !! all-wet and stable.  Lower values give stronger topographic
         !! steering at the cost of thinner layers (min_depth/nz); true land
         !! barriers wait on the wet/dry-plumbing follow-up.
      real(wp), parameter :: PI = 4.0_wp*atan(1.0_wp)
      real(wp) :: x, y, d_frac, nxp, nyp
      integer :: i, j, i_phys, j_phys, ng, ioff, joff

      ng = grid%nghost
      ioff = grid%i_offset_global
      joff = grid%j_offset_global
      nxp = real(grid%nx_global, wp)
      nyp = real(grid%ny_global, wp)

      do j = 1, size(b, 2)
         ! GLOBAL physical indices of this local cell.
         j_phys = j - ng + joff
         y = (real(j_phys, wp) - 0.5_wp)/nyp
         do i = 1, size(b, 1)
            i_phys = i - ng + ioff
            x = (real(i_phys, wp) - 0.5_wp)/nxp
            d_frac = 1.0_wp &
                     - 1.1_wp*nw2_spike(y - 1.0_wp, 0.12_wp) &   ! great northern wall
                     - 1.1_wp*nw2_spike(y, 0.12_wp) &            ! Antarctica (south wall)
                     - nl_continent_amp*( &
                     (1.2_wp*nw2_spike(x, 0.2_wp) + 1.2_wp*nw2_spike(x - 1.0_wp, 0.2_wp)) &
                     *nw2_spike(min(0.0_wp, y - 0.3_wp), 0.2_wp) &                          ! South America
                     + 1.2_wp*nw2_spike(x - 0.5_wp, 0.2_wp) &
                     *nw2_spike(min(0.0_wp, y - 0.55_wp), 0.2_wp) &                         ! Africa
                     + 1.2_wp*(nw2_spike(x, 0.12_wp) + nw2_spike(x - 1.0_wp, 0.12_wp)) &
                     *nw2_spike(max(0.0_wp, y - 0.06_wp), 0.12_wp) &                        ! Antarctic Peninsula
                     + 0.1_wp*(nw2_cosbell(x, 0.1_wp) + nw2_cosbell(x - 1.0_wp, 0.1_wp)) &  ! Drake Passage ridge
                     + 0.5_wp*nw2_cosbell(x - 0.16_wp, 0.05_wp)*(nw2_cosbell(y - 0.18_wp, 0.13_wp)**0.4_wp) &  ! Scotia Arc E
                     + 0.4_wp*(nw2_cosbell(x - 0.09_wp, 0.08_wp)**0.4_wp)*nw2_cosbell(y - 0.26_wp, 0.05_wp) &  ! Scotia Arc N
                     + 0.4_wp*(nw2_cosbell(x - 0.08_wp, 0.08_wp)**0.4_wp)*nw2_cosbell(y - 0.1_wp, 0.05_wp)) &  ! Scotia Arc S
                     - nl_roughness_amp*cos(14.0_wp*PI*x)*sin(14.0_wp*PI*y) &  ! roughness
                     - nl_roughness_amp*cos(20.0_wp*PI*x)*cos(20.0_wp*PI*y)
            if (d_frac < 0.0_wp) d_frac = 0.0_wp
            b(i, j) = max(d_frac*max_depth, min_depth)
         end do
      end do
   end subroutine set_bathymetry_neverworld2