configure_ocean_sponge Subroutine

   subroutine configure_ocean_sponge(cfg, ocean_state, grid, compute_rank)
      !! Populate `ocean_state%sponge`'s per-cell `idamp_h`/`idamp_u`/
      !! `idamp_v` maps from `cfg%ocean%sponge` (`&ocean_sponge_nml`) + the
      !! already-configured `ocean_state%bc` edge tags. Run AFTER
      !! `configure_ocean_bc` (reads the edge tags + `has_*` flags) and
      !! BEFORE `ocean_state_enter_data`.
      !!
      !! `damp_source = "band"` (the only value implemented in v1): cosine
      !! ramp from every `OBC_SPONGE`-tagged edge, at the exact cell/u-face/
      !! v-face offsets the legacy `rdb_ocean_sponge::ocean_sponge_apply{,
      !! _tracers}` kernels use (§3.2 of the plan — reproduces today's band
      !! so `enable=.true., damp_source="band"` is *physically* the same
      !! sponge as today). Overlapping edges (corners) SUM their rates — the
      !! exact continuation of the legacy kernel's SEQUENTIAL per-edge decay
      !! (`exp(-a*dt)*exp(-b*dt) = exp(-(a+b)*dt)`).
      !!
      !! Per-edge `west_width`/`west_strength` etc. (sentinel `< 0` ⇒
      !! inherit the single global `&ocean_bc_nml sponge_width`/
      !! `sponge_strength` already fanned out onto `bc%<edge>%sponge_*` by
      !! `configure_ocean_bc`) close the audit gap that the per-edge tag
      !! fields were previously unreachable (no config path could ever make
      !! them differ).
      !!
      !! Does NOT snapshot the reference state — `ocean_sponge_snapshot_reference`
      !! does that, called MUCH earlier from the driver (right after the IC
      !! seed, before any restart read — see that subroutine's docstring for
      !! why the ordering is load-bearing).
      !!
      !! No-op when `.not. ocean_state%sponge%enable` (default) — the maps
      !! stay at their `init`-time zero fill (or unallocated, if `enable`
      !! was false at `init` time too).
      type(config_t), intent(in) :: cfg
      type(ocean_state_t), intent(inout) :: ocean_state
      type(hgrid_t), intent(in) :: grid
      integer, intent(in) :: compute_rank

      integer :: i0, i1, j0, j1, band, ramp
      real(wp) :: strength
      integer :: nonzero_h, nonzero_u, nonzero_v

      associate (sp => ocean_state%sponge, bc => ocean_state%bc, s_cfg => cfg%ocean%sponge)

         if (.not. sp%enable) return

         sp%relax_uv = s_cfg%relax_uv
         sp%relax_tracers = s_cfg%relax_tracers
         sp%relax_h = s_cfg%relax_h
         sp%damp_source = s_cfg%damp_source
         sp%target_source = s_cfg%target_source
         sp%lin_t_ref = s_cfg%lin_t_ref
         sp%lin_dt_dz = s_cfg%lin_dt_dz
         sp%lin_s_ref = s_cfg%lin_s_ref
         sp%lin_ds_dz = s_cfg%lin_ds_dz
         ramp = merge(SPONGE_RAMP_LINEAR, SPONGE_RAMP_COSINE, &
                      trim(s_cfg%ramp) == "linear")

         ! Latch the two registry indices the analytic `linear_z` refresh
         ! needs, so the per-step kernel never walks the tracer registry
         ! (the array-of-derived-types device-indirection rule).
         sp%idx_t = ocean_state%multilayer%idx_temperature
         sp%idx_s = ocean_state%multilayer%idx_salinity

         ! Latch the geopotential depth of the column TOP.  The draft is
         ! static by design, so this is a configure-time copy and the
         ! sponge slot stays self-contained at run time.  No cavity ⇒
         ! `z_draft` is the (1,1) placeholder and `z_top` keeps its
         ! init-time zero (the open-ocean free-surface datum).
         if (ocean_state%metrics%use_cavity .and. &
             size(ocean_state%metrics%z_draft, 1) == size(sp%z_top, 1) .and. &
             size(ocean_state%metrics%z_draft, 2) == size(sp%z_top, 2)) then
            sp%z_top = ocean_state%metrics%z_draft
         else
            sp%z_top = 0.0_wp
         end if

         i0 = grid%nghost + 1
         i1 = grid%nghost + grid%nx_phys
         j0 = grid%nghost + 1
         j1 = grid%nghost + grid%ny_phys

         if (trim(sp%damp_source) == "band") then
            ! ---- West edge ----
            if (bc%west%bc_type == OBC_SPONGE .and. bc%has_west) then
               band = merge(s_cfg%west_width, bc%west%sponge_width, s_cfg%west_width >= 0)
               strength = merge(s_cfg%west_strength, bc%west%sponge_strength, &
                                s_cfg%west_strength >= 0.0_wp)
               if (band > 0) then
                  call sponge_add_band_x(sp%idamp_h, sp%idamp_u, sp%idamp_v, &
                                         wall_face=grid%nghost + 1, band=band, &
                                         strength=strength, side=+1, j0=j0, j1=j1, &
                                         ramp=ramp)
               end if
            end if
            ! ---- East edge ----
            if (bc%east%bc_type == OBC_SPONGE .and. bc%has_east) then
               band = merge(s_cfg%east_width, bc%east%sponge_width, s_cfg%east_width >= 0)
               strength = merge(s_cfg%east_strength, bc%east%sponge_strength, &
                                s_cfg%east_strength >= 0.0_wp)
               if (band > 0) then
                  call sponge_add_band_x(sp%idamp_h, sp%idamp_u, sp%idamp_v, &
                                         wall_face=grid%nghost + grid%nx_phys + 1, band=band, &
                                         strength=strength, side=-1, j0=j0, j1=j1, &
                                         ramp=ramp)
               end if
            end if
            ! ---- South edge ----
            if (bc%south%bc_type == OBC_SPONGE .and. bc%has_south) then
               band = merge(s_cfg%south_width, bc%south%sponge_width, s_cfg%south_width >= 0)
               strength = merge(s_cfg%south_strength, bc%south%sponge_strength, &
                                s_cfg%south_strength >= 0.0_wp)
               if (band > 0) then
                  call sponge_add_band_y(sp%idamp_h, sp%idamp_u, sp%idamp_v, &
                                         wall_face=grid%nghost + 1, band=band, &
                                         strength=strength, side=+1, i0=i0, i1=i1, &
                                         ramp=ramp)
               end if
            end if
            ! ---- North edge ----
            if (bc%north%bc_type == OBC_SPONGE .and. bc%has_north) then
               band = merge(s_cfg%north_width, bc%north%sponge_width, s_cfg%north_width >= 0)
               strength = merge(s_cfg%north_strength, bc%north%sponge_strength, &
                                s_cfg%north_strength >= 0.0_wp)
               if (band > 0) then
                  call sponge_add_band_y(sp%idamp_h, sp%idamp_u, sp%idamp_v, &
                                         wall_face=grid%nghost + grid%ny_phys + 1, band=band, &
                                         strength=strength, side=-1, i0=i0, i1=i1, &
                                         ramp=ramp)
               end if
            end if
         end if

         if (compute_rank == 0) then
            nonzero_h = count(sp%idamp_h > 0.0_wp)
            nonzero_u = count(sp%idamp_u > 0.0_wp)
            nonzero_v = count(sp%idamp_v > 0.0_wp)
            call logger%info("Sponge (map-driven): damp_source="//trim(sp%damp_source)// &
                             " target_source="//trim(sp%target_source)// &
                             " idamp_h max="//to_string(maxval(sp%idamp_h))//" 1/s"// &
                             " nonzero(h/u/v)="//to_string(nonzero_h)//"/"// &
                             to_string(nonzero_u)//"/"//to_string(nonzero_v))
         end if
      end associate
   end subroutine configure_ocean_sponge