zstar_full_build_column Subroutine

   pure subroutine zstar_full_build_column(nz, h_bed, h_surf_target, &
                                           n_surf, stretching, &
                                           z_ref_col)
      !! Build a per-column reference z-level pattern for VCOORD_ZSTAR_FULL.
      !! Output `z_ref_col(0:nz)` monotonically increasing (positive-down),
      !! z_ref_col(0)=0 surface, z_ref_col(nz)=h_bed. Built top-down; the caller
      !! handles ROMS ordering. Top `n_surf_use` layers use the stretching mode
      !! (log/uniform) toward h_surf_target; below that uniform to the bed;
      !! degenerate (h_bed≤0, nz≤0) ⇒ uniform.
      integer, intent(in) :: nz
      real(wp), intent(in) :: h_bed
         !! Local bed depth (m, positive down).  Must be > 0 in normal use.
      real(wp), intent(in) :: h_surf_target
         !! Target thickness of the surface layer (m).  ≤ 0 means "auto"
         !! and falls back to uniform spacing.
      integer, intent(in) :: n_surf
         !! Number of fine near-surface layers.  ≤ 0 means auto.
      integer, intent(in) :: stretching
         !! STRETCH_LOG or STRETCH_UNIFORM
      real(wp), intent(out) :: z_ref_col(0:nz)

      integer :: k, n_surf_use, n_coarse
      real(wp) :: h_fine, h_coarse, dz_uniform, r, base, w

      if (nz <= 0) return

      ! Degenerate / dry column: output zeros (never negative interfaces, which
      ! would propagate as negative dz_new and crash the solver via NaN CFL).
      if (h_bed <= 0.0_wp) then
         do k = 0, nz
            z_ref_col(k) = 0.0_wp
         end do
         return
      end if

      ! Auto n_surf: 1/3 of layers, at least 1, at most nz-1.
      if (n_surf <= 0) then
         n_surf_use = max(1, nz/3)
      else
         n_surf_use = max(1, min(nz - 1, n_surf))
      end if

      ! If no surface target requested, just uniform.
      if (h_surf_target <= 0.0_wp .or. nz == 1) then
         z_ref_col(0) = 0.0_wp
         do k = 1, nz
            z_ref_col(k) = h_bed*real(k, wp)/real(nz, wp)
         end do
         return
      end if

      n_coarse = nz - n_surf_use
      z_ref_col(0) = 0.0_wp

      select case (stretching)
      case (STRETCH_LOG)
         ! Geometric fine zone: layer k thickness = h_surf_target*r^(k-1). Pick r
         ! (bisection, monotonic) so the bottom fine layer matches the coarse
         ! thickness and fine+coarse sums = h_bed.
         if (n_surf_use == 1 .or. n_coarse == 0) then
            ! Single fine layer or no coarse — fall through to UNIFORM behaviour
            r = 1.0_wp
         else
            block
               real(wp) :: r_lo, r_hi, r_mid, f_mid, target_ratio
               integer :: it
               target_ratio = h_bed/h_surf_target
               r_lo = 1.000001_wp
               r_hi = 10.0_wp
               do it = 1, 60
                  r_mid = 0.5_wp*(r_lo + r_hi)
                  ! f(r) = (r^n - 1)/(r - 1) + r^(n-1) * n_coarse
                  f_mid = (r_mid**n_surf_use - 1.0_wp)/(r_mid - 1.0_wp) &
                          + r_mid**(n_surf_use - 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)
            end block
         end if
         w = 1.0_wp
         base = 0.0_wp
         do k = 1, n_surf_use
            base = base + h_surf_target*w
            z_ref_col(k) = base
            w = w*r
         end do
         ! Coarse zone: continue with dz = h_surf_target * r^(n_surf-1)
         dz_uniform = h_surf_target*r**(n_surf_use - 1)
         do k = n_surf_use + 1, nz
            z_ref_col(k) = z_ref_col(k - 1) + dz_uniform
         end do
      case default  ! STRETCH_UNIFORM
         ! Fine zone: n_surf_use layers each of thickness h_surf_target.
         ! Coarse zone: uniform fill of the remainder.
         h_fine = h_surf_target*real(n_surf_use, wp)
         if (h_fine >= h_bed) then
            ! Column too shallow to honour h_surf_target for all fine layers
            ! — fall back to uniform spacing across the whole column.
            do k = 1, nz
               z_ref_col(k) = h_bed*real(k, wp)/real(nz, wp)
            end do
            return
         end if
         h_coarse = h_bed - h_fine
         do k = 1, n_surf_use
            z_ref_col(k) = h_surf_target*real(k, wp)
         end do
         if (n_coarse > 0) then
            dz_uniform = h_coarse/real(n_coarse, wp)
            do k = n_surf_use + 1, nz
               z_ref_col(k) = h_fine + dz_uniform*real(k - n_surf_use, wp)
            end do
         end if
      end select

      ! Lock the bed interface exactly to h_bed (in case of round-off).
      z_ref_col(nz) = h_bed
   end subroutine zstar_full_build_column