vcoord_target_dz_column_zstar_full Subroutine

   pure subroutine vcoord_target_dz_column_zstar_full(nz, H, z_ref_col, &
                                                      h_min, dz)
      !$acc routine seq
      !! Compute target layer thicknesses for VCOORD_ZSTAR_FULL.
      !! Inputs: z_ref_col(0:nz) local reference (top-down, 0=surface, nz=h_bed);
      !! H current total depth (m) = h_bed + η; h_min vanishing-layer floor (m).
      !! Output: dz(1:nz) ROMS-ordered (dz(1) bottom, dz(nz) surface),
      !! sum(dz) = H exactly, vanishing rows set to h_min.
      !! Surface layer absorbs η; if H < h_bed the deepest layers clip to h_min
      !! and the surface is trimmed to keep sum = H.
      integer, intent(in) :: nz
      real(wp), intent(in) :: H
      real(wp), intent(in) :: z_ref_col(0:nz)
      real(wp), intent(in) :: h_min
      real(wp), intent(out) :: dz(nz)

      real(wp) :: dz_top_kt, z_upper, z_lower, H_eff, h_bed_ref, eta
      real(wp) :: sum_dz, deficit
      integer :: kt, k

      h_bed_ref = z_ref_col(nz)
      eta = H - h_bed_ref
      H_eff = max(H, 0.0_wp)

      if (eta >= 0.0_wp) then
         ! Column at or above reference depth.  Subsurface layers keep
         ! their reference thicknesses; the surface layer absorbs the
         ! SSH offset.  sum(dz) = h_bed_ref + eta = H exactly.
         do kt = 1, nz
            k = nz - kt + 1
            dz_top_kt = max(z_ref_col(kt) - z_ref_col(kt - 1), 0.0_wp)
            if (kt == 1) then
               dz(k) = dz_top_kt + eta     ! ROMS surface gets +η
            else
               dz(k) = dz_top_kt
            end if
         end do
      else
         ! Column shallower than reference (H < h_bed_ref). Walk top-down:
         ! layer above bed keeps full thickness; straddling layer gets
         ! H_eff - z_upper; below-bed layer is vanishing (h_min). Then trim the
         ! surface so sum = H exactly (remap mass conservation needs this).
         do kt = 1, nz
            k = nz - kt + 1
            z_upper = z_ref_col(kt - 1)
            z_lower = z_ref_col(kt)
            if (z_lower <= H_eff) then
               dz(k) = z_lower - z_upper
            else if (z_upper < H_eff) then
               dz(k) = H_eff - z_upper
            else
               dz(k) = h_min
            end if
         end do
         ! Re-balance via the surface layer.
         sum_dz = 0.0_wp
         do k = 1, nz
            sum_dz = sum_dz + dz(k)
         end do
         deficit = sum_dz - H_eff
         if (deficit > 0.0_wp) then
            if (dz(nz) - deficit >= h_min) then
               dz(nz) = dz(nz) - deficit
            else
               ! Too thin for the floor — set surface to h_min, let downstream
               ! physics guards (DRY_TOLERANCE) handle it; sum = H not enforced.
               dz(nz) = h_min
            end if
         end if
      end if
   end subroutine vcoord_target_dz_column_zstar_full