roquet_pcm_dpa_intz Subroutine

public pure subroutine roquet_pcm_dpa_intz(t, s, e_top, dz, rho0, rho_ref, dpa, intz_dpa)

Vertical integral of the Roquet et al. (2015) SpV in-situ density anomaly over one constant-T/S (PCM) layer: the 5-point Boole rule of boole_dpa_intz_layer (MOM6 int_density_dz_generic_pcm, which is also what MOM6 runs for ROQUET_SPV under a Boussinesq PGF), with the EOS FACTORED.

WHY NOT A CLOSED FORM. The Boussinesq FV PGF integrates DENSITY in height, int rho dz with p = -g*rho0*z. Roquet’s polynomial is in SPECIFIC VOLUME, SV(T, S, p) – a degree-6 polynomial in p – so the integrand is 1/SV(p), a rational function of z with no useful antiderivative (partial fractions over six complex roots). The integral that IS exact in closed form, int SV dp, is the NON-Boussinesq one (MOM6 int_spec_vol_dp); roundabout does not integrate it. MOM6 sets EOS_QUADRATURE = True by default for every Roquet/TEOS-10 form for the same reason.

WHAT IS FACTORED. In a PCM layer T and S are the same at all five Boole points, so everything expensive in the EOS – two sqrt, the PT->CT polynomial, the ~50-term (T, S) sums – is the same at all five. rdb_roq_ts_coeffs computes it ONCE; each point then costs the degree-6 pressure Horner and one division. The five densities are the numbers boole_dpa_intz_layer evaluates (up to its wt_t*t + wt_b*t blend of equal edge values, which can move T by an ulp), and the Boole weights and moment are its own, so dpa / intz_dpa agree with the generic path to round-off. Cost per layer: 1 T/S polynomial + 5 pressure Horners, against 5 full roquet_spv_point calls (each also computing the derivatives the density throws away) through the generic eos_t dispatch.

ACCURACY. The rule’s truncation error is set by the pressure curvature of 1/SV over the layer: at most 1.1e-13 of g*rho0*dz for any layer up to 1000 m thick and 8e-10 for the whole 6000 m column as one layer (test_ocean_pgf_eos_fast, against 64-panel Gauss-Legendre of the generic EOS).

Arguments

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

Layer potential temperature (degC), constant through the layer.
real(kind=wp), intent(in) :: s

Layer practical salinity (PSU), constant through the layer.
real(kind=wp), intent(in) :: e_top

Height of the SHALLOWER interface (m, geopotential).
real(kind=wp), intent(in) :: dz

Layer thickness (m, >= 0); the layer spans [e_top - dz, e_top].
real(kind=wp), intent(in) :: rho0

Boussinesq reference density of the pressure estimate (kg/m^3).
real(kind=wp), intent(in) :: rho_ref

Anomaly reference subtracted from the in-situ density (kg/m^3).
real(kind=wp), intent(out) :: dpa

g * int (rho - rho_ref) dz over the layer (Pa).
real(kind=wp), intent(out) :: intz_dpa

First moment from the top (Pa*m), as boole_dpa_intz_layer.

Calls

proc~~roquet_pcm_dpa_intz~~CallsGraph proc~roquet_pcm_dpa_intz roquet_pcm_dpa_intz rdb_roq_spv_p rdb_roq_spv_p proc~roquet_pcm_dpa_intz->rdb_roq_spv_p rdb_roq_ts_coeffs rdb_roq_ts_coeffs proc~roquet_pcm_dpa_intz->rdb_roq_ts_coeffs

Called by

proc~~roquet_pcm_dpa_intz~~CalledByGraph proc~roquet_pcm_dpa_intz roquet_pcm_dpa_intz proc~compute_fv_mom6_insitu_pcm_impl compute_fv_mom6_insitu_pcm_impl proc~compute_fv_mom6_insitu_pcm_impl->proc~roquet_pcm_dpa_intz proc~roquet_pcm_dpa_face roquet_pcm_dpa_face proc~compute_fv_mom6_insitu_pcm_impl->proc~roquet_pcm_dpa_face proc~roquet_pcm_dpa_face->proc~roquet_pcm_dpa_intz proc~ocean_pressure_force_compute ocean_pressure_force_compute proc~ocean_pressure_force_compute->proc~compute_fv_mom6_insitu_pcm_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 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 :: gxrho
real(kind=wp), private :: r1
real(kind=wp), private :: r2
real(kind=wp), private :: r3
real(kind=wp), private :: r4
real(kind=wp), private :: r5
real(kind=wp), private :: rho_anom
real(kind=wp), private :: sv0
real(kind=wp), private :: sv1
real(kind=wp), private :: sv2
real(kind=wp), private :: sv3

Source Code

   pure subroutine roquet_pcm_dpa_intz(t, s, e_top, dz, rho0, rho_ref, dpa, intz_dpa)
      !$acc routine seq
      !! Vertical integral of the Roquet et al. (2015) SpV in-situ density
      !! anomaly over one constant-T/S (PCM) layer: the 5-point Boole rule
      !! of `boole_dpa_intz_layer` (MOM6 `int_density_dz_generic_pcm`, which
      !! is also what MOM6 runs for `ROQUET_SPV` under a Boussinesq PGF),
      !! with the EOS FACTORED.
      !!
      !! WHY NOT A CLOSED FORM.  The Boussinesq FV PGF integrates DENSITY in
      !! height, `int rho dz` with `p = -g*rho0*z`.  Roquet's polynomial is
      !! in SPECIFIC VOLUME, `SV(T, S, p)` -- a degree-6 polynomial in `p`
      !! -- so the integrand is `1/SV(p)`, a rational function of `z` with
      !! no useful antiderivative (partial fractions over six complex
      !! roots).  The integral that IS exact in closed form, `int SV dp`, is
      !! the NON-Boussinesq one (MOM6 `int_spec_vol_dp`); roundabout does not
      !! integrate it.  MOM6 sets `EOS_QUADRATURE = True` by default for
      !! every Roquet/TEOS-10 form for the same reason.
      !!
      !! WHAT IS FACTORED.  In a PCM layer T and S are the same at all five
      !! Boole points, so everything expensive in the EOS -- two sqrt, the
      !! PT->CT polynomial, the ~50-term (T, S) sums -- is the same at all
      !! five.  `rdb_roq_ts_coeffs` computes it ONCE; each point then
      !! costs the degree-6 pressure Horner and one division.  The five
      !! densities are the numbers `boole_dpa_intz_layer` evaluates (up to
      !! its `wt_t*t + wt_b*t` blend of equal edge values, which can move
      !! T by an ulp), and the Boole weights and moment are its own, so
      !! `dpa` / `intz_dpa` agree with the generic path to round-off.
      !! Cost per layer: 1 T/S polynomial + 5 pressure Horners, against 5
      !! full `roquet_spv_point` calls (each also computing the derivatives
      !! the density throws away) through the generic `eos_t` dispatch.
      !!
      !! ACCURACY.  The rule's truncation error is set by the pressure
      !! curvature of `1/SV` over the layer: at most 1.1e-13 of `g*rho0*dz`
      !! for any layer up to 1000 m thick and 8e-10 for the whole 6000 m
      !! column as one layer (`test_ocean_pgf_eos_fast`, against 64-panel
      !! Gauss-Legendre of the generic EOS).
      real(wp), intent(in)  :: t
         !! Layer potential temperature (degC), constant through the layer.
      real(wp), intent(in)  :: s
         !! Layer practical salinity (PSU), constant through the layer.
      real(wp), intent(in)  :: e_top
         !! Height of the SHALLOWER interface (m, geopotential).
      real(wp), intent(in)  :: dz
         !! Layer thickness (m, >= 0); the layer spans `[e_top - dz, e_top]`.
      real(wp), intent(in)  :: rho0
         !! Boussinesq reference density of the pressure estimate (kg/m^3).
      real(wp), intent(in)  :: rho_ref
         !! Anomaly reference subtracted from the in-situ density (kg/m^3).
      real(wp), intent(out) :: dpa
         !! `g * int (rho - rho_ref) dz` over the layer (Pa).
      real(wp), intent(out) :: intz_dpa
         !! First moment from the top (Pa*m), as `boole_dpa_intz_layer`.

      real(wp) :: sv0, sv1, sv2, sv3, gxrho, rho_anom
      real(wp) :: r1, r2, r3, r4, r5

      call rdb_roq_ts_coeffs(t, s, sv0, sv1, sv2, sv3)
      gxrho = GRAVITY*rho0
      ! The five Boole points, top (n = 1) to bottom (n = 5), at the same
      ! `z5 = e_top - 0.25*(n-1)*dz` as `boole_dpa_intz_layer`.
      r1 = 1.0_wp/rdb_roq_spv_p(sv0, sv1, sv2, sv3, -gxrho*e_top) - rho_ref
      r2 = 1.0_wp/rdb_roq_spv_p(sv0, sv1, sv2, sv3, -gxrho*(e_top - 0.25_wp*dz)) - rho_ref
      r3 = 1.0_wp/rdb_roq_spv_p(sv0, sv1, sv2, sv3, -gxrho*(e_top - 0.5_wp*dz)) - rho_ref
      r4 = 1.0_wp/rdb_roq_spv_p(sv0, sv1, sv2, sv3, -gxrho*(e_top - 0.75_wp*dz)) - rho_ref
      r5 = 1.0_wp/rdb_roq_spv_p(sv0, sv1, sv2, sv3, -gxrho*(e_top - dz)) - rho_ref

      rho_anom = (1.0_wp/90.0_wp)*(7.0_wp*(r1 + r5) + 32.0_wp*(r2 + r4) + 12.0_wp*r3)
      dpa = GRAVITY*dz*rho_anom
      intz_dpa = 0.5_wp*GRAVITY*dz*dz*(rho_anom &
                                       - (1.0_wp/90.0_wp)*(16.0_wp*(r4 - r2) + 7.0_wp*(r5 - r1)))
   end subroutine roquet_pcm_dpa_intz