vcoord_target_dz_column Subroutine

   pure subroutine vcoord_target_dz_column(coord_type, nz, H, dsig, z_ref, &
                                           depth_transition, blend_width, dz)
      !$acc routine seq
      !! Compute target layer thicknesses for a single water column.
      !! Pure, called from `do concurrent` (one thread per column).
      !!
      !! Output `dz` is bottom-up (ROMS): dz(1) = BOTTOM, dz(nz) = SURFACE;
      !! matches the solvers' `h_layer` indexing (consumers write
      !! h_layer(k,...) = dz(k) with no reversal). sum(dz) = H in all cases.
      !!   SIGMA:       dz(k) = dsig(k)*H (terrain-following).
      !!   ZSIGMA:      smooth blend sigma (shallow) → fixed z-levels (deep).
      !!   ZSTAR:       z*-lite, dz(k) = (z_ref(nz-k+1)-z_ref(nz-k))*H/z_ref(nz);
      !!                interfaces stay at fixed relative position as η changes.
      !!   ZSTAR_SIGMA: smoothstep blend of sigma (shallow) and z*-lite (deep);
      !!                each branch sums to H so no surface trim needed.
      integer, intent(in) :: coord_type
      integer, intent(in) :: nz
      real(wp), intent(in) :: H
         !! Total water depth at this column (m)
      real(wp), intent(in) :: dsig(nz)
         !! Reference sigma fractions (sum = 1), ROMS-ordered: dsig(1) bottom,
         !! dsig(nz) surface. Currently uniform 1/nz.
      real(wp), intent(in) :: z_ref(0:nz)
         !! Reference z-level interface depths (m, positive down).
         !! `z_ref(0) = 0` is the surface; `z_ref(nz)` is the deepest
         !! reference interface.  Used for VCOORD_ZSIGMA and VCOORD_ZSTAR.
      real(wp), intent(in) :: depth_transition
         !! Depth (m) below which blending begins
      real(wp), intent(in) :: blend_width
         !! Width of the blending zone (m)
      real(wp), intent(out) :: dz(nz)
         !! Output target layer thicknesses, ROMS-ordered (sum = H,
         !! `dz(1)` = bottom, `dz(nz)` = surface)

      real(wp) :: alpha, x, z_top_k, z_bot_k, dz_z, dz_sum, deficit
      integer :: k

      select case (coord_type)

      case (VCOORD_SIGMA)
         ! Pure terrain-following.  Uniform dsig means orientation is
         ! immaterial — the output is the same in either direction.
         do k = 1, nz
            dz(k) = dsig(k)*H
         end do

      case (VCOORD_ZSTAR)
         ! z*-lite: stretch the global reference pattern uniformly so
         ! sum(dz) = H (H = h_bed + η). ROMS-ordered: dz(1) deepest, dz(nz)
         ! surface. Reduces to uniform sigma for uniform z_ref.
         if (z_ref(nz) > 0.0_wp) then
            do k = 1, nz
               dz_z = z_ref(nz - k + 1) - z_ref(nz - k)
               dz(k) = dz_z*H/z_ref(nz)
            end do
         else
            ! Degenerate z_ref — fall back to uniform sigma
            do k = 1, nz
               dz(k) = dsig(k)*H
            end do
         end if

      case (VCOORD_ZSIGMA)
         ! Smooth blend sigma↔z-levels: alpha=0 shallow (H<=depth_transition,
         ! pure sigma) → alpha=1 deep (H>=depth_transition+blend_width, z-levels).
         ! ROMS order: layer k spans z_ref(nz-k)..z_ref(nz-k+1).
         if (H <= depth_transition) then
            ! Pure sigma — uniform dsig, orientation immaterial
            do k = 1, nz
               dz(k) = dsig(k)*H
            end do
         else
            ! Compute blending factor
            if (blend_width > 0.0_wp) then
               x = (H - depth_transition)/blend_width
               x = max(0.0_wp, min(1.0_wp, x))
               alpha = x*x*(3.0_wp - 2.0_wp*x)  ! smoothstep
            else
               alpha = 1.0_wp
            end if

            ! Compute z-level thicknesses (clip to column depth)
            dz_sum = 0.0_wp
            do k = 1, nz
               z_top_k = min(z_ref(nz - k), H)        ! shallower interface
               z_bot_k = min(z_ref(nz - k + 1), H)    ! deeper interface
               dz_z = max(z_bot_k - z_top_k, 0.0_wp)
               ! Blend: (1-alpha)*sigma + alpha*z-level
               dz(k) = (1.0_wp - alpha)*dsig(k)*H + alpha*dz_z
               dz_sum = dz_sum + dz(k)
            end do

            ! When H is shallower than the deepest z_ref, layers near the
            ! bottom (k=1 under ROMS) get clipped to zero.  Put the
            ! deficit back into the bottom layer so sum(dz) = H.
            deficit = H - dz_sum
            if (abs(deficit) > 0.0_wp) then
               dz(1) = dz(1) + deficit
            end if
         end if

      case (VCOORD_ZSTAR_SIGMA)
         ! z*/sigma smoothstep blend: alpha=0 shallow (pure sigma) → alpha=1 deep
         ! (pure z*-lite). Both branches sum to H so the blend does, no trim.
         if (H <= depth_transition .or. z_ref(nz) <= 0.0_wp) then
            ! Pure sigma — including the degenerate-z_ref fallback so
            ! mass conservation holds without an extra dz_sum fix-up.
            do k = 1, nz
               dz(k) = dsig(k)*H
            end do
         else
            if (blend_width > 0.0_wp) then
               x = (H - depth_transition)/blend_width
               x = max(0.0_wp, min(1.0_wp, x))
               alpha = x*x*(3.0_wp - 2.0_wp*x)  ! smoothstep
            else
               alpha = 1.0_wp
            end if
            ! z*-lite contribution: dz_z = (z_ref(nz-k+1)-z_ref(nz-k)) * H/z_ref(nz)
            do k = 1, nz
               dz_z = (z_ref(nz - k + 1) - z_ref(nz - k))*H/z_ref(nz)
               dz(k) = (1.0_wp - alpha)*dsig(k)*H + alpha*dz_z
            end do
         end if

      case default
         ! Fallback to sigma (uniform, orientation immaterial)
         do k = 1, nz
            dz(k) = dsig(k)*H
         end do

      end select
   end subroutine vcoord_target_dz_column