eos_wright_impl Subroutine

private pure subroutine eos_wright_impl(h_layer, hS_layer, hT_layer, rho_layer, rho_0, p_ref, nx, ny, nz)

Wright (1997) rational EOS evaluated at a single reference pressure p_ref, which is a SCALAR by design — see the horizontal-uniformity contract in eos_compute_arrays. Same outer-shim signature as eos_linear_impl — bare 3D arrays, vanishing-layer fallback to rho_0.

α_0(T, S) = a0 + a1T + a2S p_0(T, S) = b0 + b1T + b2T^2 + b3T^3 + b4S + b5ST λ (T, S) = c0 + c1T + c2T^2 + c3T^3 + c4S + c5ST

ρ = (P + p_0) / (λ + α_0 * (P + p_0))

Reference behaviour (T=10, S=35, P=0): ρ ≈ 1027.3 kg/m^3, within 0.01 kg/m^3 of the surface seawater density used across the MOM6 test suite.

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) :: p_ref
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz

Calls

proc~~eos_wright_impl~~CallsGraph proc~eos_wright_impl eos_wright_impl local local proc~eos_wright_impl->local

Called by

proc~~eos_wright_impl~~CalledByGraph proc~eos_wright_impl eos_wright_impl proc~eos_compute_arrays eos_compute_arrays proc~eos_compute_arrays->proc~eos_wright_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_cu
real(kind=wp), private :: T_k
real(kind=wp), private :: T_sq
real(kind=wp), private :: alpha_0
real(kind=wp), private :: denom
integer, private :: i
real(kind=wp), private :: inv_h
integer, private :: j
integer, private :: k
real(kind=wp), private :: lambda
real(kind=wp), private :: p_0
real(kind=wp), private :: p_plus_p0

Source Code

   pure subroutine eos_wright_impl(h_layer, hS_layer, hT_layer, rho_layer, &
                                   rho_0, p_ref, nx, ny, nz)
      !! Wright (1997) rational EOS evaluated at a single reference
      !! pressure `p_ref`, which is a SCALAR by design — see the
      !! horizontal-uniformity contract in `eos_compute_arrays`.  Same
      !! outer-shim signature as `eos_linear_impl` — bare 3D arrays,
      !! vanishing-layer fallback to `rho_0`.
      !!
      !!   α_0(T, S) = a0 + a1*T + a2*S
      !!   p_0(T, S) = b0 + b1*T + b2*T^2 + b3*T^3 + b4*S + b5*S*T
      !!   λ  (T, S) = c0 + c1*T + c2*T^2 + c3*T^3 + c4*S + c5*S*T
      !!
      !!   ρ = (P + p_0) / (λ + α_0 * (P + p_0))
      !!
      !! Reference behaviour (T=10, S=35, P=0): ρ ≈ 1027.3 kg/m^3,
      !! within 0.01 kg/m^3 of the surface seawater density used
      !! across the MOM6 test suite.
      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, p_ref

      integer :: i, j, k
      real(wp) :: inv_h, S_k, T_k, T_sq, T_cu
      real(wp) :: alpha_0, p_0, lambda, p_plus_p0, denom

      do concurrent(k=1:nz, j=1:ny, i=1:nx) &
         local(inv_h, S_k, T_k, T_sq, T_cu, &
               alpha_0, p_0, lambda, p_plus_p0, denom)
         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
            T_sq = T_k*T_k
            T_cu = T_sq*T_k

            alpha_0 = WRIGHT_A0 + WRIGHT_A1*T_k + WRIGHT_A2*S_k
            p_0 = WRIGHT_B0 + WRIGHT_B1*T_k + WRIGHT_B2*T_sq + WRIGHT_B3*T_cu + &
                  WRIGHT_B4*S_k + WRIGHT_B5*S_k*T_k
            lambda = WRIGHT_C0 + WRIGHT_C1*T_k + WRIGHT_C2*T_sq + WRIGHT_C3*T_cu + &
                     WRIGHT_C4*S_k + WRIGHT_C5*S_k*T_k

            p_plus_p0 = p_ref + p_0
            denom = lambda + alpha_0*p_plus_p0
            rho_layer(i, j, k) = p_plus_p0/denom
         else
            rho_layer(i, j, k) = rho_0
         end if
      end do
   end subroutine eos_wright_impl