ocean_vcoord_build_zref_full Subroutine

   pure subroutine ocean_vcoord_build_zref_full(this, h_bed)
      !! Populate `z_ref(i, j, 0:nz_ml)` per column from the local
      !! bathymetry `h_bed(i, j)`.  Mirrors `zstar_full_build_column`
      !! from `src/ALE/rdb_vcoord.F90` but as a 2D loop owned by this
      !! slot — keeps the coastal helper untouched while letting the
      !! ocean path own its z_ref lifecycle.
      !!
      !! Top-down indexing inside the column: `z_ref(:, :, 0) = 0` is
      !! the surface, `z_ref(:, :, nz_ml) = h_bed(:, :)` is the bed.
      !! `compute_target_h(VCOORD_ZSTAR_FULL)` walks this table and
      !! emits ROMS-ordered `target_h(:, :, 1..nz_ml)`.
      !!
      !! Degenerate columns (`h_bed ≤ 0`) get a column of zeros — the
      !! ZSTAR_FULL branch then produces an all-h_min vanishing-layer
      !! result which the downstream dry-cell guards already handle.
      !!
      !! Call sites: once at setup, then any time bathymetry changes
      !! (which today is "never" — the ocean path doesn't move the bed).
      class(ocean_vcoord_t), intent(inout) :: this
      real(wp), intent(in) :: h_bed(:, :)
         !! Bed depth at cell centres (m, positive-down).
      integer  :: n_surf_eff, n_coarse, nz, i, j, k
      real(wp) :: h_surf_eff, h_fine, h_coarse, dz_uniform, r, base, w
      real(wp) :: r_lo, r_hi, r_mid, f_mid, target_ratio
      integer  :: it

      if (.not. this%is_init) return
      nz = this%nz_ml

      h_surf_eff = this%zstar_h_surf_target
      if (this%zstar_n_surf <= 0) then
         n_surf_eff = max(1, nz/3)
      else
         n_surf_eff = max(1, min(nz - 1, this%zstar_n_surf))
      end if
      n_coarse = nz - n_surf_eff

      ! The per-column work is column-local — different (i, j) cells
      ! don't read each other.  The do-concurrent `local()` clause keeps
      ! the scalar scratch private per thread; the bisection branch in
      ! the LOG stretching case requires a sequential inner loop so we
      ! keep it as a non-concurrent block.
      ! One-time z_ref setup runs on the HOST (plain do, not do concurrent):
      ! it writes this%z_ref, a vcoord allocatable component that is not yet
      ! device-mapped at seed time (build_zref_full runs before
      ! ocean_state_enter_data).  A device kernel writing an un-present
      ! derived-type component faults under OpenMP-target offload (the
      ! component can't be implicitly mapped); stdpar tolerates it but this
      ! is one-time init, so host is correct and costs nothing.
      do j = 1, this%ny_total
         do i = 1, this%nx_total
         if (h_bed(i, j) <= 0.0_wp) then
            ! Degenerate / dry column — leave z_ref at zero.
            do k = 0, nz
               this%z_ref(i, j, k) = 0.0_wp
            end do
         else if (h_surf_eff <= 0.0_wp .or. nz == 1) then
            ! Uniform spacing fallback.
            this%z_ref(i, j, 0) = 0.0_wp
            do k = 1, nz
               this%z_ref(i, j, k) = h_bed(i, j)*real(k, wp)/real(nz, wp)
            end do
         else
            this%z_ref(i, j, 0) = 0.0_wp
            if (this%zstar_stretching == STRETCH_LOG &
                .and. n_surf_eff >= 2 .and. n_coarse >= 1) then
               ! Geometric fine zone: layer k thickness = h_surf · r^(k-1).
               ! Bisect for r so total = h_bed.
               target_ratio = h_bed(i, j)/h_surf_eff
               r_lo = 1.000001_wp
               r_hi = 10.0_wp
               do it = 1, 60
                  r_mid = 0.5_wp*(r_lo + r_hi)
                  f_mid = (r_mid**n_surf_eff - 1.0_wp)/(r_mid - 1.0_wp) &
                          + r_mid**(n_surf_eff - 1)*real(n_coarse, wp)
                  if (f_mid > target_ratio) then
                     r_hi = r_mid
                  else
                     r_lo = r_mid
                  end if
                  if (r_hi - r_lo < 1.0e-9_wp) exit
               end do
               r = 0.5_wp*(r_lo + r_hi)
               w = 1.0_wp
               base = 0.0_wp
               do k = 1, n_surf_eff
                  base = base + h_surf_eff*w
                  this%z_ref(i, j, k) = base
                  w = w*r
               end do
               dz_uniform = h_surf_eff*r**(n_surf_eff - 1)
               do k = n_surf_eff + 1, nz
                  this%z_ref(i, j, k) = this%z_ref(i, j, k - 1) + dz_uniform
               end do
            else
               ! Uniform fine zone + uniform coarse fill.
               h_fine = h_surf_eff*real(n_surf_eff, wp)
               if (h_fine >= h_bed(i, j)) then
                  ! Column too shallow to honour h_surf for all fine layers.
                  ! Match MOM6 isopycnal NK=2: surface absorbs the water,
                  ! deeper layers vanish.  Reserve `zstar_h_min` for each
                  ! vanishing layer so target_h ≥ h_min downstream — kernels
                  ! that divide by h_layer can't see exact zero (would NaN).
                  ! Fine layers fill top-down at h_surf_eff each; bottom-most
                  ! fine grabs the residual; coarse vanishes to h_min.
                  base = 0.0_wp
                  do k = 1, n_surf_eff
                     base = base + h_surf_eff
                     ! Leave room for h_min of every layer below this one.
                     if (base > h_bed(i, j) - real(nz - k, wp)*this%zstar_h_min) then
                        base = h_bed(i, j) - real(nz - k, wp)*this%zstar_h_min
                     end if
                     this%z_ref(i, j, k) = base
                  end do
                  do k = n_surf_eff + 1, nz
                     this%z_ref(i, j, k) = this%z_ref(i, j, k - 1) + this%zstar_h_min
                  end do
               else
                  h_coarse = h_bed(i, j) - h_fine
                  do k = 1, n_surf_eff
                     this%z_ref(i, j, k) = h_surf_eff*real(k, wp)
                  end do
                  if (n_coarse > 0) then
                     dz_uniform = h_coarse/real(n_coarse, wp)
                     do k = n_surf_eff + 1, nz
                        this%z_ref(i, j, k) = h_fine + dz_uniform*real(k - n_surf_eff, wp)
                     end do
                  end if
               end if
            end if
            ! Pin the bed interface to h_bed exactly (round-off guard).
            this%z_ref(i, j, nz) = h_bed(i, j)
         end if
         end do
      end do
   end subroutine ocean_vcoord_build_zref_full