boundary_edges_linear Subroutine

private pure subroutine boundary_edges_linear(h_self, h_nbr, q_self, dq_up, q_t, q_b)

Linear-exact one-sided edge pair for a BOUNDARY layer (k=1 or k=nz), where a centred slope has no second neighbour.

dq_up is the layer-mean increment toward the SURFACE across the two cell centres (q(2)-q(1) at the bed, q(nz)-q(nz-1) at the surface). The centres are (h_self + h_nbr)/2 apart, so the per-metre slope is dq_up/((h_self+h_nbr)/2) and the half-jump across this layer is

d = dq_up * h_self / (h_self + h_nbr)

giving q_t = q + d (shallower edge) and q_b = q - d. For a profile that is linear in z this reproduces the true edge values EXACTLY, for any thickness pair — which is the property the FV pressure-gradient quadrature needs (Adcroft, Hallberg & Harrison 2008; White, Adcroft & Hallberg 2009 §2): a PCM flatten here leaves the full terrain-following truncation error in the layers next to the tilted boundary.

Limiter: |d| <= |dq_up|, i.e. the edge never leaves the interval the two cell means span on the other side. Since h_self/(h_self+h_nbr) < 1 it never bites on a real thickness pair — it is armour against a degenerate h_nbr <= 0, and it keeps the extrapolation from manufacturing a density inversion.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_self

Thickness of the boundary layer itself (m).
real(kind=wp), intent(in) :: h_nbr

Thickness of its single interior neighbour (m).
real(kind=wp), intent(in) :: q_self

Layer mean of the boundary layer.
real(kind=wp), intent(in) :: dq_up

Layer-mean increment toward the surface, neighbour -> self at the surface layer, self -> neighbour at the bed layer.
real(kind=wp), intent(out) :: q_t

Top (shallower) edge value.
real(kind=wp), intent(out) :: q_b

Bottom (deeper) edge value.

Called by

proc~~boundary_edges_linear~~CalledByGraph proc~boundary_edges_linear boundary_edges_linear proc~plm_edges_layer plm_edges_layer proc~plm_edges_layer->proc~boundary_edges_linear proc~ppm_edges_layer ppm_edges_layer proc~ppm_edges_layer->proc~boundary_edges_linear proc~compute_fv_mom6_reconstruct_impl compute_fv_mom6_reconstruct_impl proc~compute_fv_mom6_reconstruct_impl->proc~plm_edges_layer proc~compute_fv_mom6_reconstruct_impl->proc~ppm_edges_layer proc~ocean_pressure_force_compute ocean_pressure_force_compute proc~ocean_pressure_force_compute->proc~compute_fv_mom6_reconstruct_impl proc~run_stage run_stage proc~run_stage->proc~ocean_pressure_force_compute proc~run_stage_split run_stage_split proc~run_stage_split->proc~ocean_pressure_force_compute proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~run_stage proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private, parameter :: H_TINY = 1.0e-30_wp
real(kind=wp), private :: d

Source Code

   pure subroutine boundary_edges_linear(h_self, h_nbr, q_self, dq_up, q_t, q_b)
      !$acc routine seq
      !! Linear-exact one-sided edge pair for a BOUNDARY layer (k=1 or
      !! k=nz), where a centred slope has no second neighbour.
      !!
      !! `dq_up` is the layer-mean increment toward the SURFACE across the
      !! two cell centres (`q(2)-q(1)` at the bed, `q(nz)-q(nz-1)` at the
      !! surface).  The centres are `(h_self + h_nbr)/2` apart, so the
      !! per-metre slope is `dq_up/((h_self+h_nbr)/2)` and the half-jump
      !! across this layer is
      !!
      !!     d = dq_up * h_self / (h_self + h_nbr)
      !!
      !! giving `q_t = q + d` (shallower edge) and `q_b = q - d`.  For a
      !! profile that is linear in z this reproduces the true edge values
      !! EXACTLY, for any thickness pair — which is the property the FV
      !! pressure-gradient quadrature needs (Adcroft, Hallberg & Harrison
      !! 2008; White, Adcroft & Hallberg 2009 §2): a PCM flatten here
      !! leaves the full terrain-following truncation error in the layers
      !! next to the tilted boundary.
      !!
      !! Limiter: `|d| <= |dq_up|`, i.e. the edge never leaves the
      !! interval the two cell means span on the other side.  Since
      !! `h_self/(h_self+h_nbr) < 1` it never bites on a real thickness
      !! pair — it is armour against a degenerate `h_nbr <= 0`, and it
      !! keeps the extrapolation from manufacturing a density inversion.
      real(wp), intent(in)  :: h_self
         !! Thickness of the boundary layer itself (m).
      real(wp), intent(in)  :: h_nbr
         !! Thickness of its single interior neighbour (m).
      real(wp), intent(in)  :: q_self
         !! Layer mean of the boundary layer.
      real(wp), intent(in)  :: dq_up
         !! Layer-mean increment toward the surface, neighbour -> self at
         !! the surface layer, self -> neighbour at the bed layer.
      real(wp), intent(out) :: q_t
         !! Top (shallower) edge value.
      real(wp), intent(out) :: q_b
         !! Bottom (deeper) edge value.

      real(wp), parameter :: H_TINY = 1.0e-30_wp
      real(wp) :: d

      d = dq_up*h_self/max(h_self + h_nbr, H_TINY)
      d = sign(min(abs(d), abs(dq_up)), d)
      q_t = q_self + d
      q_b = q_self - d
   end subroutine boundary_edges_linear