eos_linear_impl Subroutine

private pure subroutine eos_linear_impl(h_layer, hS_layer, hT_layer, rho_layer, rho_0, beta_S, S_ref, alpha_T, T_ref, nx, ny, nz)

Linear two-tracer EOS, flat-impl form:

rho_k = rho_0 + beta_S * (S_k - S_ref) - alpha_T * (T_k - T_ref)

where S_k = hS_k / h_k, T_k = hT_k / h_k. beta_S and alpha_T are pre-multiplied sensitivities in kg/m³ per unit S / T — standard seawater values are 0.78 and 0.17.

For vanishing layers (h_layer <= H_VANISHED) the cell falls back to rho_0 — same defensive branch as the coastal kernel; keeps the EOS finite under ZSTAR_FULL when bed-side layers can pinch out. The gate is > H_VANISHED (not > 0): during an active drain the PPM positivity limiter guarantees h >= 0 but NOT h >= H_VANISHED, so a layer at e.g. h = 1e-8 with hS ≈ 35·1e-8 would pass a > 0 gate and give S = hS/h ≈ 5e2 PSU ⇒ corrupted ρ ⇒ garbage PGF. > H_VANISHED (the D4 vanished-layer role) returns rho_0 for any layer in (0, H_VANISHED]. Bit-identical for any config whose layers all exceed H_VANISHED. CAVEAT: ZSTAR_FULL floors vanishing bed layers to zstar_h_min (type default 1.0e-4; the shipped namelists set 1.5e-4 == H_VANISHED exactly, and validate_config warns on anything above it for that family — see rdb_vcoord :: vcoord_h_min_role), so such a bed layer takes the rho_0 fallback instead of the computed density. That is the INTENT, not a casualty: those layers are below the bed and hold no water. A dynamically negligible change on a 0.15 mm layer (PGF contribution ~1e-4 of a normal layer); the shipped anchors (ocean_analytical 8/8, dyn_split, baroclinic_longrun, double-gyre helpers) pass unchanged on both toolchains.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: hS_layer(nx,ny,nz)
real(kind=wp), intent(in) :: hT_layer(nx,ny,nz)
real(kind=wp), intent(out) :: rho_layer(nx,ny,nz)
real(kind=wp), intent(in) :: rho_0
real(kind=wp), intent(in) :: beta_S
real(kind=wp), intent(in) :: S_ref
real(kind=wp), intent(in) :: alpha_T
real(kind=wp), intent(in) :: T_ref
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz

Calls

proc~~eos_linear_impl~~CallsGraph proc~eos_linear_impl eos_linear_impl local local proc~eos_linear_impl->local

Called by

proc~~eos_linear_impl~~CalledByGraph proc~eos_linear_impl eos_linear_impl proc~eos_compute_arrays eos_compute_arrays proc~eos_compute_arrays->proc~eos_linear_impl proc~ocean_eos_compute ocean_eos_compute proc~ocean_eos_compute->proc~eos_compute_arrays proc~rdb_ocean_set_h rdb_ocean_set_h proc~rdb_ocean_set_h->proc~ocean_eos_compute proc~rdb_ocean_set_tracer rdb_ocean_set_tracer proc~rdb_ocean_set_tracer->proc~ocean_eos_compute proc~run_stage run_stage proc~run_stage->proc~ocean_eos_compute proc~run_stage_split run_stage_split proc~run_stage_split->proc~ocean_eos_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 proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_split

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: S_k
real(kind=wp), private :: T_k
integer, private :: i
real(kind=wp), private :: inv_h
integer, private :: j
integer, private :: k

Source Code

   pure subroutine eos_linear_impl(h_layer, hS_layer, hT_layer, rho_layer, &
                                   rho_0, beta_S, S_ref, &
                                   alpha_T, T_ref, &
                                   nx, ny, nz)
      !! Linear two-tracer EOS, flat-impl form:
      !!
      !!   rho_k = rho_0 + beta_S * (S_k - S_ref) - alpha_T * (T_k - T_ref)
      !!
      !! where S_k = hS_k / h_k, T_k = hT_k / h_k.  `beta_S` and
      !! `alpha_T` are pre-multiplied sensitivities in kg/m³ per unit
      !! S / T — standard seawater values are 0.78 and 0.17.
      !!
      !! For vanishing layers (`h_layer <= H_VANISHED`) the cell falls
      !! back to rho_0 — same defensive branch as the coastal kernel; keeps
      !! the EOS finite under ZSTAR_FULL when bed-side layers can pinch out.
      !! The gate is `> H_VANISHED` (not `> 0`): during an active drain the
      !! PPM positivity limiter guarantees `h >= 0` but NOT `h >= H_VANISHED`,
      !! so a layer at e.g. `h = 1e-8` with `hS ≈ 35·1e-8` would pass a
      !! `> 0` gate and give `S = hS/h ≈ 5e2 PSU` ⇒ corrupted ρ ⇒ garbage
      !! PGF.  `> H_VANISHED` (the D4 vanished-layer role) returns rho_0 for
      !! any layer in `(0, H_VANISHED]`.  Bit-identical for any config whose
      !! layers all exceed H_VANISHED.  CAVEAT: ZSTAR_FULL floors vanishing
      !! bed layers to `zstar_h_min` (type default 1.0e-4; the shipped
      !! namelists set 1.5e-4 == H_VANISHED exactly, and `validate_config`
      !! warns on anything above it for that family — see
      !! `rdb_vcoord :: vcoord_h_min_role`), so such a bed layer takes the
      !! rho_0 fallback instead of the computed density.  That is the
      !! INTENT, not a casualty: those layers are below the bed and hold no
      !! water.  A dynamically negligible change on a 0.15 mm
      !! layer (PGF contribution ~1e-4 of a normal layer); the shipped
      !! anchors (ocean_analytical 8/8, dyn_split, baroclinic_longrun,
      !! double-gyre helpers) pass unchanged on both toolchains.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: hS_layer(nx, ny, nz)
      real(wp), intent(in) :: hT_layer(nx, ny, nz)
      real(wp), intent(out) :: rho_layer(nx, ny, nz)
      real(wp), intent(in) :: rho_0, beta_S, S_ref, alpha_T, T_ref

      integer :: i, j, k
      real(wp) :: inv_h, S_k, T_k

      do concurrent(k=1:nz, j=1:ny, i=1:nx) local(inv_h, S_k, T_k)
         if (h_layer(i, j, k) > H_VANISHED) then
            inv_h = 1.0_wp/h_layer(i, j, k)
            S_k = hS_layer(i, j, k)*inv_h
            T_k = hT_layer(i, j, k)*inv_h
         else
            S_k = S_ref
            T_k = T_ref
         end if
         rho_layer(i, j, k) = rho_0 + beta_S*(S_k - S_ref) - alpha_T*(T_k - T_ref)
      end do
   end subroutine eos_linear_impl