ocean_vcoord_geometric_target Subroutine

private pure subroutine ocean_vcoord_geometric_target(coord_type, nx, ny, nz, target_h, total_h, eta, dsig, z_ref_global, z_ref, zsigma_depth_transition, zsigma_blend_width, zstar_h_min)

Geometric target-grid kernels (EULERIAN_Z, SIGMA, ZSIGMA, ZSTAR_SIGMA, ZSTAR_FULL; formulae documented on ocean_vcoord_compute_target_h_impl). Flat on purpose — every array an explicit-shape dummy, every knob a scalar dummy — so no derived-type component and no associate-name reaches a do concurrent (see ocean_vcoord_rho_target for the GPU fault that shape caused there).

Arguments

Type IntentOptional Attributes Name
integer, intent(in), value :: coord_type

VCOORD_* family (not LAGRANGIAN / Z_FIXED / ZSTAR: the dispatcher owns those).
integer, intent(in), value :: nx

i-extent of every horizontal array (total, incl. halos).
integer, intent(in), value :: ny

j-extent of every horizontal array (total, incl. halos).
integer, intent(in), value :: nz

Number of layers; k = 1 is the bed, k = nz the surface.
real(kind=wp), intent(inout) :: target_h(nx,ny,nz)

Target layer thickness (m), bottom-up.
real(kind=wp), intent(in) :: total_h(nx,ny)

Column-total depth H (m).
real(kind=wp), intent(in) :: eta(nx,ny)

Free-surface anomaly η (m).
real(kind=wp), intent(in) :: dsig(nz)

Nominal layer fractions, bottom-up.
real(kind=wp), intent(in) :: z_ref_global(0:nz)

Global reference interface depths (ZSIGMA / ZSTAR_SIGMA).
real(kind=wp), intent(in) :: z_ref(nx,ny,0:nz)

Per-column reference interface depths (ZSTAR_FULL).
real(kind=wp), intent(in), value :: zsigma_depth_transition

Sigma → z* transition depth (m).
real(kind=wp), intent(in), value :: zsigma_blend_width

Smoothstep blend width (m).
real(kind=wp), intent(in), value :: zstar_h_min

Vanished-layer thickness (m).

Calls

proc~~ocean_vcoord_geometric_target~~CallsGraph proc~ocean_vcoord_geometric_target ocean_vcoord_geometric_target local local proc~ocean_vcoord_geometric_target->local

Called by

proc~~ocean_vcoord_geometric_target~~CalledByGraph proc~ocean_vcoord_geometric_target ocean_vcoord_geometric_target proc~ocean_vcoord_compute_target_h_impl ocean_vcoord_compute_target_h_impl proc~ocean_vcoord_compute_target_h_impl->proc~ocean_vcoord_geometric_target proc~ocean_vcoord_eta0_target ocean_vcoord_eta0_target proc~ocean_vcoord_eta0_target->proc~ocean_vcoord_geometric_target proc~configure_ocean_closed_faces configure_ocean_closed_faces proc~configure_ocean_closed_faces->proc~ocean_vcoord_eta0_target proc~ocean_state_seed_from_cfg ocean_state_seed_from_cfg proc~ocean_state_seed_from_cfg->proc~ocean_vcoord_eta0_target proc~ocean_vcoord_compute_target_h ocean_vcoord_t%ocean_vcoord_compute_target_h proc~ocean_vcoord_compute_target_h->proc~ocean_vcoord_compute_target_h_impl proc~engine_setup engine_setup proc~engine_setup->proc~configure_ocean_closed_faces proc~engine_setup->proc~ocean_state_seed_from_cfg proc~ocean_apply_ale_remap_centres ocean_apply_ale_remap_centres proc~ocean_apply_ale_remap_centres->proc~ocean_vcoord_compute_target_h proc~ocean_apply_ale_remap_step ocean_apply_ale_remap_step proc~ocean_apply_ale_remap_step->proc~ocean_vcoord_compute_target_h proc~complete_ocean_create complete_ocean_create proc~complete_ocean_create->proc~engine_setup proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_setup proc~driver_validate driver_validate proc~driver_validate->proc~engine_setup proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~ocean_apply_ale_remap_step proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step_split proc~rdb_ocean_create_finalize rdb_ocean_create_finalize proc~rdb_ocean_create_finalize->proc~complete_ocean_create proc~rdb_ocean_create_from_string rdb_ocean_create_from_string proc~rdb_ocean_create_from_string->proc~complete_ocean_create

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: H_eff
real(kind=wp), private :: alpha
real(kind=wp), private :: column_total
real(kind=wp), private :: deficit
real(kind=wp), private :: dz_sum
real(kind=wp), private :: dz_z
real(kind=wp), private :: eta_loc
real(kind=wp), private :: h_bed_ref
integer, private :: i
integer, private :: j
integer, private :: k
real(kind=wp), private :: sum_dz
real(kind=wp), private :: x
real(kind=wp), private :: z_bot_k
real(kind=wp), private :: z_lower
real(kind=wp), private :: z_ref_nz_inv
real(kind=wp), private :: z_top_k
real(kind=wp), private :: z_upper

Source Code

   pure subroutine ocean_vcoord_geometric_target(coord_type, nx, ny, nz, target_h, &
                                                 total_h, eta, dsig, z_ref_global, &
                                                 z_ref, zsigma_depth_transition, &
                                                 zsigma_blend_width, zstar_h_min)
      !! Geometric target-grid kernels (EULERIAN_Z, SIGMA, ZSIGMA,
      !! ZSTAR_SIGMA, ZSTAR_FULL; formulae documented on
      !! `ocean_vcoord_compute_target_h_impl`).  Flat on purpose — every
      !! array an explicit-shape dummy, every knob a scalar dummy — so no
      !! derived-type component and no `associate`-name reaches a
      !! `do concurrent` (see `ocean_vcoord_rho_target` for the GPU fault
      !! that shape caused there).
      integer, intent(in), value :: coord_type
         !! `VCOORD_*` family (not LAGRANGIAN / Z_FIXED / ZSTAR: the
         !! dispatcher owns those).
      integer, intent(in), value :: nx
         !! i-extent of every horizontal array (total, incl. halos).
      integer, intent(in), value :: ny
         !! j-extent of every horizontal array (total, incl. halos).
      integer, intent(in), value :: nz
         !! Number of layers; `k = 1` is the bed, `k = nz` the surface.
      real(wp), intent(inout) :: target_h(nx, ny, nz)
         !! Target layer thickness (m), bottom-up.
      real(wp), intent(in) :: total_h(nx, ny)
         !! Column-total depth H (m).
      real(wp), intent(in) :: eta(nx, ny)
         !! Free-surface anomaly η (m).
      real(wp), intent(in) :: dsig(nz)
         !! Nominal layer fractions, bottom-up.
      real(wp), intent(in) :: z_ref_global(0:nz)
         !! Global reference interface depths (ZSIGMA / ZSTAR_SIGMA).
      real(wp), intent(in) :: z_ref(nx, ny, 0:nz)
         !! Per-column reference interface depths (ZSTAR_FULL).
      real(wp), intent(in), value :: zsigma_depth_transition
         !! Sigma → z* transition depth (m).
      real(wp), intent(in), value :: zsigma_blend_width
         !! Smoothstep blend width (m).
      real(wp), intent(in), value :: zstar_h_min
         !! Vanished-layer thickness (m).
      integer :: i, j, k
      real(wp) :: column_total, alpha, x, z_top_k, z_bot_k, dz_z, dz_sum, deficit
      real(wp) :: z_ref_nz_inv
      real(wp) :: h_bed_ref, eta_loc, H_eff, z_upper, z_lower, sum_dz

      select case (coord_type)

      case (VCOORD_EULERIAN_Z)
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            target_h(i, j, k) = total_h(i, j)*dsig(k)
         end do

      case (VCOORD_SIGMA)
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            column_total = total_h(i, j) + eta(i, j)
            target_h(i, j, k) = column_total*dsig(k)
         end do

      case (VCOORD_ZSIGMA)
         ! Smoothstep blend: sigma in shallow, fixed z-levels in deep.
         ! Mirrors `vcoord_target_dz_column` in `src/ALE/rdb_vcoord.F90`
         ! but built directly on the 2D (H, η) fields.  The deep-branch
         ! z-level intervals are clipped to the local column total so
         ! sum_k target_h = H + η exactly even when the column is shallower
         ! than the deepest reference interface; any residual deficit is
         ! deposited in the bed-side layer (k=1) to preserve the sum.
         do concurrent(j=1:ny, i=1:nx) &
            local(column_total, alpha, x, k, z_top_k, z_bot_k, dz_z, dz_sum, deficit)
            column_total = total_h(i, j) + eta(i, j)
            if (column_total <= zsigma_depth_transition) then
               do k = 1, nz
                  target_h(i, j, k) = dsig(k)*column_total
               end do
            else
               if (zsigma_blend_width > 0.0_wp) then
                  x = (column_total - zsigma_depth_transition)/zsigma_blend_width
                  x = max(0.0_wp, min(1.0_wp, x))
                  alpha = x*x*(3.0_wp - 2.0_wp*x)
               else
                  alpha = 1.0_wp
               end if
               dz_sum = 0.0_wp
               do k = 1, nz
                  z_top_k = min(z_ref_global(nz - k), column_total)
                  z_bot_k = min(z_ref_global(nz - k + 1), column_total)
                  dz_z = max(z_bot_k - z_top_k, 0.0_wp)
                  target_h(i, j, k) = (1.0_wp - alpha)*dsig(k)*column_total &
                                      + alpha*dz_z
                  dz_sum = dz_sum + target_h(i, j, k)
               end do
               deficit = column_total - dz_sum
               target_h(i, j, 1) = target_h(i, j, 1) + deficit
            end if
         end do

      case (VCOORD_ZSTAR_SIGMA)
         z_ref_nz_inv = 0.0_wp
         if (z_ref_global(nz) > 0.0_wp) then
            z_ref_nz_inv = 1.0_wp/z_ref_global(nz)
         end if
         do concurrent(j=1:ny, i=1:nx) &
            local(column_total, alpha, x, k, dz_z)
            column_total = total_h(i, j) + eta(i, j)
            if (column_total <= zsigma_depth_transition .or. z_ref_nz_inv == 0.0_wp) then
               do k = 1, nz
                  target_h(i, j, k) = dsig(k)*column_total
               end do
            else
               if (zsigma_blend_width > 0.0_wp) then
                  x = (column_total - zsigma_depth_transition)/zsigma_blend_width
                  x = max(0.0_wp, min(1.0_wp, x))
                  alpha = x*x*(3.0_wp - 2.0_wp*x)
               else
                  alpha = 1.0_wp
               end if
               do k = 1, nz
                  dz_z = (z_ref_global(nz - k + 1) - z_ref_global(nz - k)) &
                         *column_total*z_ref_nz_inv
                  target_h(i, j, k) = (1.0_wp - alpha)*dsig(k)*column_total &
                                      + alpha*dz_z
               end do
            end if
         end do

      case (VCOORD_ZSTAR_FULL)
         ! Per-column z*-full: walk the cached `z_ref(i, j, 0:nz)` from
         ! `build_zref_full`.  Surface layer absorbs η when η ≥ 0;
         ! bed-side layers vanish to `zstar_h_min` and the surface gets
         ! trimmed for exact conservation when η < 0.  Mirrors
         ! `vcoord_target_dz_column_zstar_full` (coastal) but emits
         ! ROMS order (k=1 bed, k=nz surface) directly.
         do concurrent(j=1:ny, i=1:nx) &
            local(k, h_bed_ref, eta_loc, H_eff, z_upper, z_lower, sum_dz, deficit)
            h_bed_ref = z_ref(i, j, nz)
            eta_loc = (total_h(i, j) + eta(i, j)) - h_bed_ref
            H_eff = max(total_h(i, j) + eta(i, j), 0.0_wp)
            if (h_bed_ref <= 0.0_wp) then
               ! Degenerate column: emit a single vanishing-layer stack.
               do k = 1, nz
                  target_h(i, j, k) = zstar_h_min
               end do
            else if (eta_loc >= 0.0_wp) then
               ! Column at or above reference: subsurface = z_ref intervals,
               ! surface (k=nz) gets the +η.
               do k = 1, nz
                  ! k_top = nz - k + 1 in the top-down z_ref convention.
                  target_h(i, j, k) = max( &
                                      z_ref(i, j, nz - k + 1) - z_ref(i, j, nz - k), &
                                      0.0_wp)
               end do
               target_h(i, j, nz) = target_h(i, j, nz) + eta_loc
            else
               ! Column shallower than reference (η < 0).  Walk top-down,
               ! clip layers to H_eff, vanish below.
               do k = 1, nz
                  z_upper = z_ref(i, j, nz - k)
                  z_lower = z_ref(i, j, nz - k + 1)
                  if (z_lower <= H_eff) then
                     target_h(i, j, k) = z_lower - z_upper
                  else if (z_upper < H_eff) then
                     target_h(i, j, k) = H_eff - z_upper
                  else
                     target_h(i, j, k) = zstar_h_min
                  end if
               end do
               ! Surface trim: drop the vanishing-layer overhead from the
               ! surface to make sum = H exactly.  If the surface would
               ! itself fall below h_min, leave it at h_min and let the
               ! downstream dry-cell guards handle the deficit.
               sum_dz = 0.0_wp
               do k = 1, nz
                  sum_dz = sum_dz + target_h(i, j, k)
               end do
               deficit = sum_dz - H_eff
               if (deficit > 0.0_wp) then
                  if (target_h(i, j, nz) - deficit >= zstar_h_min) then
                     target_h(i, j, nz) = target_h(i, j, nz) - deficit
                  else
                     target_h(i, j, nz) = zstar_h_min
                  end if
               end if
            end if
         end do

      case default
         error stop "ocean_vcoord_geometric_target: unknown coord_type."
      end select

   end subroutine ocean_vcoord_geometric_target