MOM6
Data Types | Functions/Subroutines
mom_eos_wright Module Reference

Detailed Description

The equation of state using the Wright 1997 expressions.

Data Types

interface  calculate_density_derivs_wright
 For a given thermodynamic state, return the derivatives of density with temperature and salinity. More...
 
interface  calculate_density_second_derivs_wright
 For a given thermodynamic state, return the second derivatives of density with various combinations of temperature, salinity, and pressure. More...
 
interface  calculate_density_wright
 Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to a reference density, from salinity (in psu), potential temperature (in deg C), and pressure [Pa], using the expressions from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. More...
 
interface  calculate_spec_vol_wright
 Compute the in situ specific volume of sea water (in [m3 kg-1]), or an anomaly with respect to a reference specific volume, from salinity (in psu), potential temperature (in deg C), and pressure [Pa], using the expressions from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. More...
 

Functions/Subroutines

subroutine calculate_density_scalar_wright (T, S, pressure, rho, rho_ref)
 This subroutine computes the in situ density of sea water (rho in [kg m-3]) from salinity (S [PSU]), potential temperature (T [degC]), and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. More...
 
subroutine calculate_density_array_wright (T, S, pressure, rho, start, npts, rho_ref)
 This subroutine computes the in situ density of sea water (rho in [kg m-3]) from salinity (S [PSU]), potential temperature (T [degC]), and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. More...
 
subroutine calculate_spec_vol_scalar_wright (T, S, pressure, specvol, spv_ref)
 This subroutine computes the in situ specific volume of sea water (specvol in [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC]) and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. If spv_ref is present, specvol is an anomaly from spv_ref. More...
 
subroutine calculate_spec_vol_array_wright (T, S, pressure, specvol, start, npts, spv_ref)
 This subroutine computes the in situ specific volume of sea water (specvol in [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC]) and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. If spv_ref is present, specvol is an anomaly from spv_ref. More...
 
subroutine calculate_density_derivs_array_wright (T, S, pressure, drho_dT, drho_dS, start, npts)
 For a given thermodynamic state, return the thermal/haline expansion coefficients. More...
 
subroutine calculate_density_derivs_scalar_wright (T, S, pressure, drho_dT, drho_dS)
 The scalar version of calculate_density_derivs which promotes scalar inputs to a 1-element array and then demotes the output back to a scalar. More...
 
subroutine calculate_density_second_derivs_array_wright (T, S, P, drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
 Second derivatives of density with respect to temperature, salinity, and pressure. More...
 
subroutine calculate_density_second_derivs_scalar_wright (T, S, P, drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp)
 Second derivatives of density with respect to temperature, salinity, and pressure for scalar inputs. Inputs promoted to 1-element array and output demoted to scalar. More...
 
subroutine, public calculate_specvol_derivs_wright (T, S, pressure, dSV_dT, dSV_dS, start, npts)
 For a given thermodynamic state, return the partial derivatives of specific volume with temperature and salinity. More...
 
subroutine, public calculate_compress_wright (T, S, pressure, rho, drho_dp, start, npts)
 This subroutine computes the in situ density of sea water (rho in [kg m-3]) and the compressibility (drho/dp = C_sound^-2) (drho_dp [s2 m-2]) from salinity (sal in psu), potential temperature (T [degC]), and pressure [Pa]. It uses the expressions from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. Coded by R. Hallberg, 1/01. More...
 
subroutine, public int_density_dz_wright (T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale)
 This subroutine calculates analytical and nearly-analytical integrals of pressure anomalies across layers, which are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. More...
 
subroutine, public int_spec_vol_dp_wright (T, S, p_t, p_b, spv_ref, HI, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, useMassWghtInterp, SV_scale, pres_scale)
 This subroutine calculates analytical and nearly-analytical integrals in pressure across layers of geopotential anomalies, which are required for calculating the finite-volume form pressure accelerations in a non-Boussinesq model. There are essentially no free assumptions, apart from the use of Boole's rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34. More...
 

Variables

real, parameter a0 = 7.057924e-4
 Parameters in the Wright equation of state.
 
real, parameter a1 = 3.480336e-7
 Parameters in the Wright equation of state.
 
real, parameter a2 = -1.112733e-7
 Parameters in the Wright equation of state.
 
real, parameter b0 = 5.790749e8
 Parameters in the Wright equation of state.
 
real, parameter b1 = 3.516535e6
 Parameters in the Wright equation of state.
 
real, parameter b2 = -4.002714e4
 Parameters in the Wright equation of state.
 
real, parameter b3 = 2.084372e2
 Parameters in the Wright equation of state.
 
real, parameter b4 = 5.944068e5
 Parameters in the Wright equation of state.
 
real, parameter b5 = -9.643486e3
 Parameters in the Wright equation of state.
 
real, parameter c0 = 1.704853e5
 Parameters in the Wright equation of state.
 
real, parameter c1 = 7.904722e2
 Parameters in the Wright equation of state.
 
real, parameter c2 = -7.984422
 Parameters in the Wright equation of state.
 
real, parameter c3 = 5.140652e-2
 Parameters in the Wright equation of state.
 
real, parameter c4 = -2.302158e2
 Parameters in the Wright equation of state.
 
real, parameter c5 = -3.079464
 Parameters in the Wright equation of state.
 

Function/Subroutine Documentation

◆ calculate_compress_wright()

subroutine, public mom_eos_wright::calculate_compress_wright ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(inout)  rho,
real, dimension(:), intent(inout)  drho_dp,
integer, intent(in)  start,
integer, intent(in)  npts 
)

This subroutine computes the in situ density of sea water (rho in [kg m-3]) and the compressibility (drho/dp = C_sound^-2) (drho_dp [s2 m-2]) from salinity (sal in psu), potential temperature (T [degC]), and pressure [Pa]. It uses the expressions from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. Coded by R. Hallberg, 1/01.

Parameters
[in]tPotential temperature relative to the surface [degC].
[in]sSalinity [PSU].
[in]pressurepressure [Pa].
[in,out]rhoIn situ density [kg m-3].
[in,out]drho_dpThe partial derivative of density with pressure (also the inverse of the square of sound speed) [s2 m-2].
[in]startThe starting point in the arrays.
[in]nptsThe number of values to calculate.

Definition at line 380 of file MOM_EOS_Wright.F90.

380  real, intent(in), dimension(:) :: t !< Potential temperature relative to the surface [degC].
381  real, intent(in), dimension(:) :: s !< Salinity [PSU].
382  real, intent(in), dimension(:) :: pressure !< pressure [Pa].
383  real, intent(inout), dimension(:) :: rho !< In situ density [kg m-3].
384  real, intent(inout), dimension(:) :: drho_dp !< The partial derivative of density with pressure
385  !! (also the inverse of the square of sound speed)
386  !! [s2 m-2].
387  integer, intent(in) :: start !< The starting point in the arrays.
388  integer, intent(in) :: npts !< The number of values to calculate.
389 
390  ! Coded by R. Hallberg, 1/01
391  ! Local variables
392  real :: al0, p0, lambda, i_denom
393  integer :: j
394 
395  do j=start,start+npts-1
396  al0 = (a0 + a1*t(j)) +a2*s(j)
397  p0 = (b0 + b4*s(j)) + t(j) * (b1 + t(j)*((b2 + b3*t(j))) + b5*s(j))
398  lambda = (c0 +c4*s(j)) + t(j) * (c1 + t(j)*((c2 + c3*t(j))) + c5*s(j))
399 
400  i_denom = 1.0 / (lambda + al0*(pressure(j) + p0))
401  rho(j) = (pressure(j) + p0) * i_denom
402  drho_dp(j) = lambda * i_denom * i_denom
403  enddo

◆ calculate_density_array_wright()

subroutine mom_eos_wright::calculate_density_array_wright ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(inout)  rho,
integer, intent(in)  start,
integer, intent(in)  npts,
real, intent(in), optional  rho_ref 
)
private

This subroutine computes the in situ density of sea water (rho in [kg m-3]) from salinity (S [PSU]), potential temperature (T [degC]), and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740.

Parameters
[in]tpotential temperature relative to the surface [degC].
[in]ssalinity [PSU].
[in]pressurepressure [Pa].
[in,out]rhoin situ density [kg m-3].
[in]startthe starting point in the arrays.
[in]nptsthe number of values to calculate.
[in]rho_refA reference density [kg m-3].

Definition at line 111 of file MOM_EOS_Wright.F90.

111  real, dimension(:), intent(in) :: t !< potential temperature relative to the surface [degC].
112  real, dimension(:), intent(in) :: s !< salinity [PSU].
113  real, dimension(:), intent(in) :: pressure !< pressure [Pa].
114  real, dimension(:), intent(inout) :: rho !< in situ density [kg m-3].
115  integer, intent(in) :: start !< the starting point in the arrays.
116  integer, intent(in) :: npts !< the number of values to calculate.
117  real, optional, intent(in) :: rho_ref !< A reference density [kg m-3].
118 
119  ! Original coded by R. Hallberg, 7/00, anomaly coded in 3/18.
120  ! Local variables
121  real :: al0, p0, lambda
122  real :: al_ts, p_tsp, lam_ts, pa_000
123  integer :: j
124 
125  if (present(rho_ref)) pa_000 = (b0*(1.0 - a0*rho_ref) - rho_ref*c0)
126  if (present(rho_ref)) then ; do j=start,start+npts-1
127  al_ts = a1*t(j) +a2*s(j)
128  al0 = a0 + al_ts
129  p_tsp = pressure(j) + (b4*s(j) + t(j) * (b1 + (t(j)*(b2 + b3*t(j)) + b5*s(j))))
130  lam_ts = c4*s(j) + t(j) * (c1 + (t(j)*(c2 + c3*t(j)) + c5*s(j)))
131 
132  ! The following two expressions are mathematically equivalent.
133  ! rho(j) = (b0 + p0_TSp) / ((c0 + lam_TS) + al0*(b0 + p0_TSp)) - rho_ref
134  rho(j) = (pa_000 + (p_tsp - rho_ref*(p_tsp*al0 + (b0*al_ts + lam_ts)))) / &
135  ( (c0 + lam_ts) + al0*(b0 + p_tsp) )
136  enddo ; else ; do j=start,start+npts-1
137  al0 = (a0 + a1*t(j)) +a2*s(j)
138  p0 = (b0 + b4*s(j)) + t(j) * (b1 + t(j)*(b2 + b3*t(j)) + b5*s(j))
139  lambda = (c0 +c4*s(j)) + t(j) * (c1 + t(j)*(c2 + c3*t(j)) + c5*s(j))
140  rho(j) = (pressure(j) + p0) / (lambda + al0*(pressure(j) + p0))
141  enddo ; endif
142 

◆ calculate_density_derivs_array_wright()

subroutine mom_eos_wright::calculate_density_derivs_array_wright ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(inout)  drho_dT,
real, dimension(:), intent(inout)  drho_dS,
integer, intent(in)  start,
integer, intent(in)  npts 
)
private

For a given thermodynamic state, return the thermal/haline expansion coefficients.

Parameters
[in]tPotential temperature relative to the surface [degC].
[in]sSalinity [PSU].
[in]pressurepressure [Pa].
[in,out]drho_dtThe partial derivative of density with potential temperature [kg m-3 degC-1].
[in,out]drho_dsThe partial derivative of density with salinity, in [kg m-3 PSU-1].
[in]startThe starting point in the arrays.
[in]nptsThe number of values to calculate.

Definition at line 200 of file MOM_EOS_Wright.F90.

200  real, intent(in), dimension(:) :: t !< Potential temperature relative to the
201  !! surface [degC].
202  real, intent(in), dimension(:) :: s !< Salinity [PSU].
203  real, intent(in), dimension(:) :: pressure !< pressure [Pa].
204  real, intent(inout), dimension(:) :: drho_dt !< The partial derivative of density with potential
205  !! temperature [kg m-3 degC-1].
206  real, intent(inout), dimension(:) :: drho_ds !< The partial derivative of density with salinity,
207  !! in [kg m-3 PSU-1].
208  integer, intent(in) :: start !< The starting point in the arrays.
209  integer, intent(in) :: npts !< The number of values to calculate.
210 
211  ! Local variables
212  real :: al0, p0, lambda, i_denom2
213  integer :: j
214 
215  do j=start,start+npts-1
216  al0 = (a0 + a1*t(j)) + a2*s(j)
217  p0 = (b0 + b4*s(j)) + t(j) * (b1 + t(j)*((b2 + b3*t(j))) + b5*s(j))
218  lambda = (c0 +c4*s(j)) + t(j) * (c1 + t(j)*((c2 + c3*t(j))) + c5*s(j))
219 
220  i_denom2 = 1.0 / (lambda + al0*(pressure(j) + p0))
221  i_denom2 = i_denom2 *i_denom2
222  drho_dt(j) = i_denom2 * &
223  (lambda* (b1 + t(j)*(2.0*b2 + 3.0*b3*t(j)) + b5*s(j)) - &
224  (pressure(j)+p0) * ( (pressure(j)+p0)*a1 + &
225  (c1 + t(j)*(c2*2.0 + c3*3.0*t(j)) + c5*s(j)) ))
226  drho_ds(j) = i_denom2 * (lambda* (b4 + b5*t(j)) - &
227  (pressure(j)+p0) * ( (pressure(j)+p0)*a2 + (c4 + c5*t(j)) ))
228  enddo
229 

◆ calculate_density_derivs_scalar_wright()

subroutine mom_eos_wright::calculate_density_derivs_scalar_wright ( real, intent(in)  T,
real, intent(in)  S,
real, intent(in)  pressure,
real, intent(out)  drho_dT,
real, intent(out)  drho_dS 
)
private

The scalar version of calculate_density_derivs which promotes scalar inputs to a 1-element array and then demotes the output back to a scalar.

Parameters
[in]tPotential temperature relative to the surface [degC].
[in]sSalinity [PSU].
[in]pressurepressure [Pa].
[out]drho_dtThe partial derivative of density with potential temperature [kg m-3 degC-1].
[out]drho_dsThe partial derivative of density with salinity, in [kg m-3 PSU-1].

Definition at line 235 of file MOM_EOS_Wright.F90.

235  real, intent(in) :: t !< Potential temperature relative to the surface [degC].
236  real, intent(in) :: s !< Salinity [PSU].
237  real, intent(in) :: pressure !< pressure [Pa].
238  real, intent(out) :: drho_dt !< The partial derivative of density with potential
239  !! temperature [kg m-3 degC-1].
240  real, intent(out) :: drho_ds !< The partial derivative of density with salinity,
241  !! in [kg m-3 PSU-1].
242 
243  ! Local variables needed to promote the input/output scalars to 1-element arrays
244  real, dimension(1) :: t0, s0, p0
245  real, dimension(1) :: drdt0, drds0
246 
247  t0(1) = t
248  s0(1) = s
249  p0(1) = pressure
250  call calculate_density_derivs_array_wright(t0, s0, p0, drdt0, drds0, 1, 1)
251  drho_dt = drdt0(1)
252  drho_ds = drds0(1)
253 

◆ calculate_density_scalar_wright()

subroutine mom_eos_wright::calculate_density_scalar_wright ( real, intent(in)  T,
real, intent(in)  S,
real, intent(in)  pressure,
real, intent(out)  rho,
real, intent(in), optional  rho_ref 
)
private

This subroutine computes the in situ density of sea water (rho in [kg m-3]) from salinity (S [PSU]), potential temperature (T [degC]), and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740.

Parameters
[in]tPotential temperature relative to the surface [degC].
[in]sSalinity [PSU].
[in]pressurepressure [Pa].
[out]rhoIn situ density [kg m-3].
[in]rho_refA reference density [kg m-3].

Definition at line 81 of file MOM_EOS_Wright.F90.

81  real, intent(in) :: t !< Potential temperature relative to the surface [degC].
82  real, intent(in) :: s !< Salinity [PSU].
83  real, intent(in) :: pressure !< pressure [Pa].
84  real, intent(out) :: rho !< In situ density [kg m-3].
85  real, optional, intent(in) :: rho_ref !< A reference density [kg m-3].
86 
87 ! *====================================================================*
88 ! * This subroutine computes the in situ density of sea water (rho in *
89 ! * [kg m-3]) from salinity (S [PSU]), potential temperature *
90 ! * (T [degC]), and pressure [Pa]. It uses the expression from *
91 ! * Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. *
92 ! * Coded by R. Hallberg, 7/00 *
93 ! *====================================================================*
94 
95  real, dimension(1) :: t0, s0, pressure0, rho0
96 
97  t0(1) = t
98  s0(1) = s
99  pressure0(1) = pressure
100 
101  call calculate_density_array_wright(t0, s0, pressure0, rho0, 1, 1, rho_ref)
102  rho = rho0(1)
103 

◆ calculate_density_second_derivs_array_wright()

subroutine mom_eos_wright::calculate_density_second_derivs_array_wright ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  P,
real, dimension(:), intent(inout)  drho_ds_ds,
real, dimension(:), intent(inout)  drho_ds_dt,
real, dimension(:), intent(inout)  drho_dt_dt,
real, dimension(:), intent(inout)  drho_ds_dp,
real, dimension(:), intent(inout)  drho_dt_dp,
integer, intent(in)  start,
integer, intent(in)  npts 
)
private

Second derivatives of density with respect to temperature, salinity, and pressure.

Parameters
[in]tPotential temperature referenced to 0 dbar [degC]
[in]sSalinity [PSU]
[in]pPressure [Pa]
[in,out]drho_ds_dsPartial derivative of beta with respect to S [kg m-3 PSU-2]
[in,out]drho_ds_dtPartial derivative of beta with respcct to T [kg m-3 PSU-1 degC-1]
[in,out]drho_dt_dtPartial derivative of alpha with respect to T [kg m-3 degC-2]
[in,out]drho_ds_dpPartial derivative of beta with respect to pressure [kg m-3 PSU-1 Pa-1]
[in,out]drho_dt_dpPartial derivative of alpha with respect to pressure [kg m-3 degC-1 Pa-1]
[in]startStarting index in T,S,P
[in]nptsNumber of points to loop over

Definition at line 259 of file MOM_EOS_Wright.F90.

259  real, dimension(:), intent(in ) :: t !< Potential temperature referenced to 0 dbar [degC]
260  real, dimension(:), intent(in ) :: s !< Salinity [PSU]
261  real, dimension(:), intent(in ) :: p !< Pressure [Pa]
262  real, dimension(:), intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect
263  !! to S [kg m-3 PSU-2]
264  real, dimension(:), intent(inout) :: drho_ds_dt !< Partial derivative of beta with respcct
265  !! to T [kg m-3 PSU-1 degC-1]
266  real, dimension(:), intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect
267  !! to T [kg m-3 degC-2]
268  real, dimension(:), intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect
269  !! to pressure [kg m-3 PSU-1 Pa-1]
270  real, dimension(:), intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect
271  !! to pressure [kg m-3 degC-1 Pa-1]
272  integer, intent(in ) :: start !< Starting index in T,S,P
273  integer, intent(in ) :: npts !< Number of points to loop over
274 
275  ! Local variables
276  real :: z0, z1, z2, z3, z4, z5, z6 ,z7, z8, z9, z10, z11, z2_2, z2_3
277  integer :: j
278  ! Based on the above expression with common terms factored, there probably exists a more numerically stable
279  ! and/or efficient expression
280 
281  do j = start,start+npts-1
282  z0 = t(j)*(b1 + b5*s(j) + t(j)*(b2 + b3*t(j)))
283  z1 = (b0 + p(j) + b4*s(j) + z0)
284  z3 = (b1 + b5*s(j) + t(j)*(2.*b2 + 2.*b3*t(j)))
285  z4 = (c0 + c4*s(j) + t(j)*(c1 + c5*s(j) + t(j)*(c2 + c3*t(j))))
286  z5 = (b1 + b5*s(j) + t(j)*(b2 + b3*t(j)) + t(j)*(b2 + 2.*b3*t(j)))
287  z6 = c1 + c5*s(j) + t(j)*(c2 + c3*t(j)) + t(j)*(c2 + 2.*c3*t(j))
288  z7 = (c4 + c5*t(j) + a2*z1)
289  z8 = (c1 + c5*s(j) + t(j)*(2.*c2 + 3.*c3*t(j)) + a1*z1)
290  z9 = (a0 + a2*s(j) + a1*t(j))
291  z10 = (b4 + b5*t(j))
292  z11 = (z10*z4 - z1*z7)
293  z2 = (c0 + c4*s(j) + t(j)*(c1 + c5*s(j) + t(j)*(c2 + c3*t(j))) + z9*z1)
294  z2_2 = z2*z2
295  z2_3 = z2_2*z2
296 
297  drho_ds_ds(j) = (z10*(c4 + c5*t(j)) - a2*z10*z1 - z10*z7)/z2_2 - (2.*(c4 + c5*t(j) + z9*z10 + a2*z1)*z11)/z2_3
298  drho_ds_dt(j) = (z10*z6 - z1*(c5 + a2*z5) + b5*z4 - z5*z7)/z2_2 - (2.*(z6 + z9*z5 + a1*z1)*z11)/z2_3
299  drho_dt_dt(j) = (z3*z6 - z1*(2.*c2 + 6.*c3*t(j) + a1*z5) + (2.*b2 + 4.*b3*t(j))*z4 - z5*z8)/z2_2 - &
300  (2.*(z6 + z9*z5 + a1*z1)*(z3*z4 - z1*z8))/z2_3
301  drho_ds_dp(j) = (-c4 - c5*t(j) - 2.*a2*z1)/z2_2 - (2.*z9*z11)/z2_3
302  drho_dt_dp(j) = (-c1 - c5*s(j) - t(j)*(2.*c2 + 3.*c3*t(j)) - 2.*a1*z1)/z2_2 - (2.*z9*(z3*z4 - z1*z8))/z2_3
303  enddo
304 

◆ calculate_density_second_derivs_scalar_wright()

subroutine mom_eos_wright::calculate_density_second_derivs_scalar_wright ( real, intent(in)  T,
real, intent(in)  S,
real, intent(in)  P,
real, intent(out)  drho_ds_ds,
real, intent(out)  drho_ds_dt,
real, intent(out)  drho_dt_dt,
real, intent(out)  drho_ds_dp,
real, intent(out)  drho_dt_dp 
)
private

Second derivatives of density with respect to temperature, salinity, and pressure for scalar inputs. Inputs promoted to 1-element array and output demoted to scalar.

Parameters
[in]tPotential temperature referenced to 0 dbar
[in]sSalinity [PSU]
[in]ppressure [Pa]
[out]drho_ds_dsPartial derivative of beta with respect to S [kg m-3 PSU-2]
[out]drho_ds_dtPartial derivative of beta with respcct to T [kg m-3 PSU-1 degC-1]
[out]drho_dt_dtPartial derivative of alpha with respect to T [kg m-3 degC-2]
[out]drho_ds_dpPartial derivative of beta with respect to pressure [kg m-3 PSU-1 Pa-1]
[out]drho_dt_dpPartial derivative of alpha with respect to pressure [kg m-3 degC-1 Pa-1]

Definition at line 311 of file MOM_EOS_Wright.F90.

311  real, intent(in ) :: t !< Potential temperature referenced to 0 dbar
312  real, intent(in ) :: s !< Salinity [PSU]
313  real, intent(in ) :: p !< pressure [Pa]
314  real, intent( out) :: drho_ds_ds !< Partial derivative of beta with respect
315  !! to S [kg m-3 PSU-2]
316  real, intent( out) :: drho_ds_dt !< Partial derivative of beta with respcct
317  !! to T [kg m-3 PSU-1 degC-1]
318  real, intent( out) :: drho_dt_dt !< Partial derivative of alpha with respect
319  !! to T [kg m-3 degC-2]
320  real, intent( out) :: drho_ds_dp !< Partial derivative of beta with respect
321  !! to pressure [kg m-3 PSU-1 Pa-1]
322  real, intent( out) :: drho_dt_dp !< Partial derivative of alpha with respect
323  !! to pressure [kg m-3 degC-1 Pa-1]
324  ! Local variables
325  real, dimension(1) :: t0, s0, p0
326  real, dimension(1) :: drdsds, drdsdt, drdtdt, drdsdp, drdtdp
327 
328  t0(1) = t
329  s0(1) = s
330  p0(1) = p
331  call calculate_density_second_derivs_array_wright(t0, s0, p0, drdsds, drdsdt, drdtdt, drdsdp, drdtdp, 1, 1)
332  drho_ds_ds = drdsds(1)
333  drho_ds_dt = drdsdt(1)
334  drho_dt_dt = drdtdt(1)
335  drho_ds_dp = drdsdp(1)
336  drho_dt_dp = drdtdp(1)
337 

◆ calculate_spec_vol_array_wright()

subroutine mom_eos_wright::calculate_spec_vol_array_wright ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(inout)  specvol,
integer, intent(in)  start,
integer, intent(in)  npts,
real, intent(in), optional  spv_ref 
)
private

This subroutine computes the in situ specific volume of sea water (specvol in [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC]) and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. If spv_ref is present, specvol is an anomaly from spv_ref.

Parameters
[in]tpotential temperature relative to the surface [degC].
[in]ssalinity [PSU].
[in]pressurepressure [Pa].
[in,out]specvolin situ specific volume [m3 kg-1].
[in]startthe starting point in the arrays.
[in]nptsthe number of values to calculate.
[in]spv_refA reference specific volume [m3 kg-1].

Definition at line 172 of file MOM_EOS_Wright.F90.

172  real, dimension(:), intent(in) :: t !< potential temperature relative to the
173  !! surface [degC].
174  real, dimension(:), intent(in) :: s !< salinity [PSU].
175  real, dimension(:), intent(in) :: pressure !< pressure [Pa].
176  real, dimension(:), intent(inout) :: specvol !< in situ specific volume [m3 kg-1].
177  integer, intent(in) :: start !< the starting point in the arrays.
178  integer, intent(in) :: npts !< the number of values to calculate.
179  real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
180 
181  ! Local variables
182  real :: al0, p0, lambda
183  integer :: j
184 
185  do j=start,start+npts-1
186  al0 = (a0 + a1*t(j)) +a2*s(j)
187  p0 = (b0 + b4*s(j)) + t(j) * (b1 + t(j)*((b2 + b3*t(j))) + b5*s(j))
188  lambda = (c0 +c4*s(j)) + t(j) * (c1 + t(j)*((c2 + c3*t(j))) + c5*s(j))
189 
190  if (present(spv_ref)) then
191  specvol(j) = (lambda + (al0 - spv_ref)*(pressure(j) + p0)) / (pressure(j) + p0)
192  else
193  specvol(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0)
194  endif
195  enddo

◆ calculate_spec_vol_scalar_wright()

subroutine mom_eos_wright::calculate_spec_vol_scalar_wright ( real, intent(in)  T,
real, intent(in)  S,
real, intent(in)  pressure,
real, intent(out)  specvol,
real, intent(in), optional  spv_ref 
)
private

This subroutine computes the in situ specific volume of sea water (specvol in [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC]) and pressure [Pa]. It uses the expression from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. If spv_ref is present, specvol is an anomaly from spv_ref.

Parameters
[in]tpotential temperature relative to the surface [degC].
[in]ssalinity [PSU].
[in]pressurepressure [Pa].
[out]specvolin situ specific volume [m3 kg-1].
[in]spv_refA reference specific volume [m3 kg-1].

Definition at line 151 of file MOM_EOS_Wright.F90.

151  real, intent(in) :: t !< potential temperature relative to the surface [degC].
152  real, intent(in) :: s !< salinity [PSU].
153  real, intent(in) :: pressure !< pressure [Pa].
154  real, intent(out) :: specvol !< in situ specific volume [m3 kg-1].
155  real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
156 
157  ! Local variables
158  real, dimension(1) :: t0, s0, pressure0, spv0
159 
160  t0(1) = t ; s0(1) = s ; pressure0(1) = pressure
161 
162  call calculate_spec_vol_array_wright(t0, s0, pressure0, spv0, 1, 1, spv_ref)
163  specvol = spv0(1)

◆ calculate_specvol_derivs_wright()

subroutine, public mom_eos_wright::calculate_specvol_derivs_wright ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(inout)  dSV_dT,
real, dimension(:), intent(inout)  dSV_dS,
integer, intent(in)  start,
integer, intent(in)  npts 
)

For a given thermodynamic state, return the partial derivatives of specific volume with temperature and salinity.

Parameters
[in]tPotential temperature relative to the surface [degC].
[in]sSalinity [PSU].
[in]pressurepressure [Pa].
[in,out]dsv_dtThe partial derivative of specific volume with potential temperature [m3 kg-1 degC-1].
[in,out]dsv_dsThe partial derivative of specific volume with salinity [m3 kg-1 / Pa].
[in]startThe starting point in the arrays.
[in]nptsThe number of values to calculate.

Definition at line 343 of file MOM_EOS_Wright.F90.

343  real, intent(in), dimension(:) :: t !< Potential temperature relative to the surface [degC].
344  real, intent(in), dimension(:) :: s !< Salinity [PSU].
345  real, intent(in), dimension(:) :: pressure !< pressure [Pa].
346  real, intent(inout), dimension(:) :: dsv_dt !< The partial derivative of specific volume with
347  !! potential temperature [m3 kg-1 degC-1].
348  real, intent(inout), dimension(:) :: dsv_ds !< The partial derivative of specific volume with
349  !! salinity [m3 kg-1 / Pa].
350  integer, intent(in) :: start !< The starting point in the arrays.
351  integer, intent(in) :: npts !< The number of values to calculate.
352 
353  ! Local variables
354  real :: al0, p0, lambda, i_denom
355  integer :: j
356 
357  do j=start,start+npts-1
358 ! al0 = (a0 + a1*T(j)) + a2*S(j)
359  p0 = (b0 + b4*s(j)) + t(j) * (b1 + t(j)*((b2 + b3*t(j))) + b5*s(j))
360  lambda = (c0 +c4*s(j)) + t(j) * (c1 + t(j)*((c2 + c3*t(j))) + c5*s(j))
361 
362  ! SV = al0 + lambda / (pressure(j) + p0)
363 
364  i_denom = 1.0 / (pressure(j) + p0)
365  dsv_dt(j) = (a1 + i_denom * (c1 + t(j)*((2.0*c2 + 3.0*c3*t(j))) + c5*s(j))) - &
366  (i_denom**2 * lambda) * (b1 + t(j)*((2.0*b2 + 3.0*b3*t(j))) + b5*s(j))
367  dsv_ds(j) = (a2 + i_denom * (c4 + c5*t(j))) - &
368  (i_denom**2 * lambda) * (b4 + b5*t(j))
369  enddo
370 

◆ int_density_dz_wright()

subroutine, public mom_eos_wright::int_density_dz_wright ( real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  T,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  S,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  z_t,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  z_b,
real, intent(in)  rho_ref,
real, intent(in)  rho_0,
real, intent(in)  G_e,
type(hor_index_type), intent(in)  HI,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(inout)  dpa,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(inout), optional  intz_dpa,
real, dimension(hi%isdb:hi%iedb,hi%jsd:hi%jed), intent(inout), optional  intx_dpa,
real, dimension(hi%isd:hi%ied,hi%jsdb:hi%jedb), intent(inout), optional  inty_dpa,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in), optional  bathyT,
real, intent(in), optional  dz_neglect,
logical, intent(in), optional  useMassWghtInterp,
real, intent(in), optional  rho_scale,
real, intent(in), optional  pres_scale 
)

This subroutine calculates analytical and nearly-analytical integrals of pressure anomalies across layers, which are required for calculating the finite-volume form pressure accelerations in a Boussinesq model.

Parameters
[in]hiThe horizontal index type for the arrays.
[in]tPotential temperature relative to the surface
[in]sSalinity [PSU].
[in]z_tHeight at the top of the layer in depth units [Z ~> m].
[in]z_bHeight at the top of the layer [Z ~> m].
[in]rho_refA mean density [R ~> kg m-3] or [kg m-3], that is subtracted out to reduce the magnitude of each of the integrals. (The pressure is calucated as p~=-z*rho_0*G_e.)
[in]rho_0Density [R ~> kg m-3] or [kg m-3], that is used to calculate the pressure (as p~=-z*rho_0*G_e) used in the equation of state.
[in]g_eThe Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2].
[in,out]dpaThe change in the pressure anomaly across the
[in,out]intz_dpaThe integral through the thickness of the layer
[in,out]intx_dpaThe integral in x of the difference between the
[in,out]inty_dpaThe integral in y of the difference between the
[in]bathytThe depth of the bathymetry [Z ~> m].
[in]dz_neglectA miniscule thickness change [Z ~> m].
[in]usemasswghtinterpIf true, uses mass weighting to interpolate T/S for top and bottom integrals.
[in]rho_scaleA multiplicative factor by which to scale density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
[in]pres_scaleA multiplicative factor to convert pressure into Pa [Pa T2 R-1 L-2 ~> 1].

Definition at line 412 of file MOM_EOS_Wright.F90.

412  type(hor_index_type), intent(in) :: hi !< The horizontal index type for the arrays.
413  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
414  intent(in) :: t !< Potential temperature relative to the surface
415  !! [degC].
416  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
417  intent(in) :: s !< Salinity [PSU].
418  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
419  intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m].
420  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
421  intent(in) :: z_b !< Height at the top of the layer [Z ~> m].
422  real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3] or [kg m-3], that is subtracted
423  !! out to reduce the magnitude of each of the integrals.
424  !! (The pressure is calucated as p~=-z*rho_0*G_e.)
425  real, intent(in) :: rho_0 !< Density [R ~> kg m-3] or [kg m-3], that is used
426  !! to calculate the pressure (as p~=-z*rho_0*G_e)
427  !! used in the equation of state.
428  real, intent(in) :: g_e !< The Earth's gravitational acceleration
429  !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2].
430  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
431  intent(inout) :: dpa !< The change in the pressure anomaly across the
432  !! layer [R L2 T-2 ~> Pa] or [Pa].
433  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
434  optional, intent(inout) :: intz_dpa !< The integral through the thickness of the layer
435  !! of the pressure anomaly relative to the anomaly
436  !! at the top of the layer [R Z L2 T-2 ~> Pa m].
437  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
438  optional, intent(inout) :: intx_dpa !< The integral in x of the difference between the
439  !! pressure anomaly at the top and bottom of the
440  !! layer divided by the x grid spacing [R L2 T-2 ~> Pa].
441  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
442  optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the
443  !! pressure anomaly at the top and bottom of the
444  !! layer divided by the y grid spacing [R L2 T-2 ~> Pa].
445  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
446  optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m].
447  real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m].
448  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
449  !! interpolate T/S for top and bottom integrals.
450  real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density
451  !! from kg m-3 to the desired units [R m3 kg-1 ~> 1]
452  real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
453  !! into Pa [Pa T2 R-1 L-2 ~> 1].
454 
455  ! Local variables
456  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d
457  real :: al0, p0, lambda
458  real :: rho_anom ! The density anomaly from rho_ref [kg m-3].
459  real :: eps, eps2, rem
460  real :: gxrho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2]
461  real :: g_earth ! The gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
462  real :: i_rho ! The inverse of the Boussinesq density [m3 kg-1]
463  real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3]
464  real :: p_ave, i_al0, i_lzz
465  real :: dz ! The layer thickness [Z ~> m].
466  real :: hwght ! A pressure-thickness below topography [Z ~> m].
467  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m].
468  real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2].
469  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
470  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
471  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
472  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
473  real :: intz(5) ! The integrals of density with height at the
474  ! 5 sub-column locations [R L2 T-2 ~> Pa].
475  real :: pa_to_rl2_t2 ! A conversion factor of pressures from Pa to the output units indicated by
476  ! pres_scale [R L2 T-2 Pa-1 ~> 1] or [1].
477  logical :: do_massweight ! Indicates whether to do mass weighting.
478  real, parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0 ! Rational constants.
479  real, parameter :: c1_9 = 1.0/9.0, c1_90 = 1.0/90.0 ! Rational constants.
480  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m
481 
482  ! These array bounds work for the indexing convention of the input arrays, but
483  ! on the computational domain defined for the output arrays.
484  isq = hi%IscB ; ieq = hi%IecB
485  jsq = hi%JscB ; jeq = hi%JecB
486  is = hi%isc ; ie = hi%iec
487  js = hi%jsc ; je = hi%jec
488 
489  if (present(pres_scale)) then
490  gxrho = pres_scale * g_e * rho_0 ; g_earth = pres_scale * g_e
491  pa_to_rl2_t2 = 1.0 / pres_scale
492  else
493  gxrho = g_e * rho_0 ; g_earth = g_e
494  pa_to_rl2_t2 = 1.0
495  endif
496  if (present(rho_scale)) then
497  g_earth = g_earth * rho_scale
498  rho_ref_mks = rho_ref / rho_scale ; i_rho = rho_scale / rho_0
499  else
500  rho_ref_mks = rho_ref ; i_rho = 1.0 / rho_0
501  endif
502 
503  do_massweight = .false.
504  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
505  do_massweight = .true.
506  ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//&
507  ! "bathyT must be present if useMassWghtInterp is present and true.")
508  ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//&
509  ! "dz_neglect must be present if useMassWghtInterp is present and true.")
510  endif ; endif
511 
512  do j=jsq,jeq+1 ; do i=isq,ieq+1
513  al0_2d(i,j) = (a0 + a1*t(i,j)) + a2*s(i,j)
514  p0_2d(i,j) = (b0 + b4*s(i,j)) + t(i,j) * (b1 + t(i,j)*((b2 + b3*t(i,j))) + b5*s(i,j))
515  lambda_2d(i,j) = (c0 +c4*s(i,j)) + t(i,j) * (c1 + t(i,j)*((c2 + c3*t(i,j))) + c5*s(i,j))
516 
517  al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
518 
519  dz = z_t(i,j) - z_b(i,j)
520  p_ave = -0.5*gxrho*(z_t(i,j)+z_b(i,j))
521 
522  i_al0 = 1.0 / al0
523  i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
524  eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
525 
526 ! rho(j) = (pressure(j) + p0) / (lambda + al0*(pressure(j) + p0))
527 
528  rho_anom = (p0 + p_ave)*(i_lzz*i_al0) - rho_ref_mks
529  rem = i_rho * (lambda * i_al0**2) * eps2 * &
530  (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
531  dpa(i,j) = pa_to_rl2_t2 * (g_earth*rho_anom*dz - 2.0*eps*rem)
532  if (present(intz_dpa)) &
533  intz_dpa(i,j) = pa_to_rl2_t2 * (0.5*g_earth*rho_anom*dz**2 - dz*(1.0+eps)*rem)
534  enddo ; enddo
535 
536  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
537  ! hWght is the distance measure by which the cell is violation of
538  ! hydrostatic consistency. For large hWght we bias the interpolation of
539  ! T & S along the top and bottom integrals, akin to thickness weighting.
540  hwght = 0.0
541  if (do_massweight) &
542  hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
543  if (hwght > 0.) then
544  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
545  hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
546  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
547  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
548  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
549  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
550  else
551  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
552  endif
553 
554  intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
555  do m=2,4
556  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
557  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
558 
559  al0 = wtt_l*al0_2d(i,j) + wtt_r*al0_2d(i+1,j)
560  p0 = wtt_l*p0_2d(i,j) + wtt_r*p0_2d(i+1,j)
561  lambda = wtt_l*lambda_2d(i,j) + wtt_r*lambda_2d(i+1,j)
562 
563  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
564  p_ave = -0.5*gxrho*(wt_l*(z_t(i,j)+z_b(i,j)) + &
565  wt_r*(z_t(i+1,j)+z_b(i+1,j)))
566 
567  i_al0 = 1.0 / al0
568  i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
569  eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
570 
571  intz(m) = pa_to_rl2_t2 * ( g_earth*dz*((p0 + p_ave)*(i_lzz*i_al0) - rho_ref_mks) - 2.0*eps * &
572  i_rho * (lambda * i_al0**2) * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2))) )
573  enddo
574  ! Use Boole's rule to integrate the values.
575  intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))
576  enddo ; enddo ; endif
577 
578  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
579  ! hWght is the distance measure by which the cell is violation of
580  ! hydrostatic consistency. For large hWght we bias the interpolation of
581  ! T & S along the top and bottom integrals, akin to thickness weighting.
582  hwght = 0.0
583  if (do_massweight) &
584  hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
585  if (hwght > 0.) then
586  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
587  hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
588  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
589  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
590  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
591  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
592  else
593  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
594  endif
595 
596  intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
597  do m=2,4
598  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
599  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
600 
601  al0 = wtt_l*al0_2d(i,j) + wtt_r*al0_2d(i,j+1)
602  p0 = wtt_l*p0_2d(i,j) + wtt_r*p0_2d(i,j+1)
603  lambda = wtt_l*lambda_2d(i,j) + wtt_r*lambda_2d(i,j+1)
604 
605  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
606  p_ave = -0.5*gxrho*(wt_l*(z_t(i,j)+z_b(i,j)) + &
607  wt_r*(z_t(i,j+1)+z_b(i,j+1)))
608 
609  i_al0 = 1.0 / al0
610  i_lzz = 1.0 / (p0 + (lambda * i_al0) + p_ave)
611  eps = 0.5*gxrho*dz*i_lzz ; eps2 = eps*eps
612 
613  intz(m) = pa_to_rl2_t2 * ( g_earth*dz*((p0 + p_ave)*(i_lzz*i_al0) - rho_ref_mks) - 2.0*eps * &
614  i_rho * (lambda * i_al0**2) * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2))) )
615  enddo
616  ! Use Boole's rule to integrate the values.
617  inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + 12.0*intz(3))
618  enddo ; enddo ; endif
619 

◆ int_spec_vol_dp_wright()

subroutine, public mom_eos_wright::int_spec_vol_dp_wright ( real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  T,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  S,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  p_t,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  p_b,
real, intent(in)  spv_ref,
type(hor_index_type), intent(in)  HI,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(inout)  dza,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(inout), optional  intp_dza,
real, dimension(hi%isdb:hi%iedb,hi%jsd:hi%jed), intent(inout), optional  intx_dza,
real, dimension(hi%isd:hi%ied,hi%jsdb:hi%jedb), intent(inout), optional  inty_dza,
integer, intent(in), optional  halo_size,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in), optional  bathyP,
real, intent(in), optional  dP_neglect,
logical, intent(in), optional  useMassWghtInterp,
real, intent(in), optional  SV_scale,
real, intent(in), optional  pres_scale 
)

This subroutine calculates analytical and nearly-analytical integrals in pressure across layers of geopotential anomalies, which are required for calculating the finite-volume form pressure accelerations in a non-Boussinesq model. There are essentially no free assumptions, apart from the use of Boole's rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34.

Parameters
[in]hiThe ocean's horizontal index type.
[in]tPotential temperature relative to the surface
[in]sSalinity [PSU].
[in]p_tPressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa].
[in]p_bPressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa].
[in]spv_refA mean specific volume that is subtracted out to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]. The calculation is mathematically identical with different values of spv_ref, but this reduces the effects of roundoff.
[in,out]dzaThe change in the geopotential anomaly across
[in,out]intp_dzaThe integral in pressure through the layer of
[in,out]intx_dzaThe integral in x of the difference between the
[in,out]inty_dzaThe integral in y of the difference between the
[in]halo_sizeThe width of halo points on which to calculate dza.
[in]bathypThe pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
[in]dp_neglectA miniscule pressure change with the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
[in]usemasswghtinterpIf true, uses mass weighting to interpolate T/S for top and bottom integrals.
[in]sv_scaleA multiplicative factor by which to scale specific volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
[in]pres_scaleA multiplicative factor to convert pressure into Pa [Pa T2 R-1 L-2 ~> 1].

Definition at line 631 of file MOM_EOS_Wright.F90.

631  type(hor_index_type), intent(in) :: hi !< The ocean's horizontal index type.
632  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
633  intent(in) :: t !< Potential temperature relative to the surface
634  !! [degC].
635  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
636  intent(in) :: s !< Salinity [PSU].
637  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
638  intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa].
639  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
640  intent(in) :: p_b !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa].
641  real, intent(in) :: spv_ref !< A mean specific volume that is subtracted out
642  !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1].
643  !! The calculation is mathematically identical with different values of
644  !! spv_ref, but this reduces the effects of roundoff.
645  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
646  intent(inout) :: dza !< The change in the geopotential anomaly across
647  !! the layer [T-2 ~> m2 s-2] or [m2 s-2].
648  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
649  optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of
650  !! the geopotential anomaly relative to the anomaly
651  !! at the bottom of the layer [R L4 T-4 ~> Pa m2 s-2]
652  !! or [Pa m2 s-2].
653  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
654  optional, intent(inout) :: intx_dza !< The integral in x of the difference between the
655  !! geopotential anomaly at the top and bottom of
656  !! the layer divided by the x grid spacing
657  !! [L2 T-2 ~> m2 s-2] or [m2 s-2].
658  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
659  optional, intent(inout) :: inty_dza !< The integral in y of the difference between the
660  !! geopotential anomaly at the top and bottom of
661  !! the layer divided by the y grid spacing
662  !! [L2 T-2 ~> m2 s-2] or [m2 s-2].
663  integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate
664  !! dza.
665  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
666  optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa]
667  real, optional, intent(in) :: dp_neglect !< A miniscule pressure change with
668  !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa]
669  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
670  !! to interpolate T/S for top and bottom integrals.
671  real, optional, intent(in) :: sv_scale !< A multiplicative factor by which to scale specific
672  !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
673  real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
674  !! into Pa [Pa T2 R-1 L-2 ~> 1].
675 
676  ! Local variables
677  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d
678  real :: al0 ! A term in the Wright EOS [R-1 ~> m3 kg-1]
679  real :: p0 ! A term in the Wright EOS [R L2 T-2 ~> Pa]
680  real :: lambda ! A term in the Wright EOS [L2 T-2 ~> m2 s-2]
681  real :: al0_scale ! Scaling factor to convert al0 from MKS units [R-1 kg m-3 ~> 1]
682  real :: p0_scale ! Scaling factor to convert p0 from MKS units [R L2 T-2 Pa-1 ~> 1]
683  real :: lam_scale ! Scaling factor to convert lambda from MKS units [L2 s2 T-2 m-2 ~> 1]
684  real :: p_ave ! The layer average pressure [R L2 T-2 ~> Pa]
685  real :: rem ! [L2 T-2 ~> m2 s-2]
686  real :: eps, eps2 ! A nondimensional ratio and its square [nondim]
687  real :: alpha_anom ! The depth averaged specific volume anomaly [R-1 ~> m3 kg-1].
688  real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa].
689  real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa].
690  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa].
691  real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2].
692  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
693  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
694  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
695  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
696  real :: intp(5) ! The integrals of specific volume with pressure at the
697  ! 5 sub-column locations [L2 T-2 ~> m2 s-2].
698  logical :: do_massweight ! Indicates whether to do mass weighting.
699  real, parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0 ! Rational constants.
700  real, parameter :: c1_9 = 1.0/9.0, c1_90 = 1.0/90.0 ! Rational constants.
701  integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, halo
702 
703  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
704  halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
705  ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
706  if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh); endif
707  if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh); endif
708 
709 
710  al0_scale = 1.0 ; if (present(sv_scale)) al0_scale = sv_scale
711  p0_scale = 1.0
712  if (present(pres_scale)) then ; if (pres_scale /= 1.0) then
713  p0_scale = 1.0 / pres_scale
714  endif ; endif
715  lam_scale = al0_scale * p0_scale
716 
717  do_massweight = .false.
718  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
719  do_massweight = .true.
720 ! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//&
721 ! "bathyP must be present if useMassWghtInterp is present and true.")
722 ! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//&
723 ! "dP_neglect must be present if useMassWghtInterp is present and true.")
724  endif ; endif
725 
726  ! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0)
727  do j=jsh,jeh ; do i=ish,ieh
728  al0_2d(i,j) = al0_scale * ( (a0 + a1*t(i,j)) + a2*s(i,j) )
729  p0_2d(i,j) = p0_scale * ( (b0 + b4*s(i,j)) + t(i,j) * (b1 + t(i,j)*((b2 + b3*t(i,j))) + b5*s(i,j)) )
730  lambda_2d(i,j) = lam_scale * ( (c0 + c4*s(i,j)) + t(i,j) * (c1 + t(i,j)*((c2 + c3*t(i,j))) + c5*s(i,j)) )
731 
732  al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j)
733  dp = p_b(i,j) - p_t(i,j)
734  p_ave = 0.5*(p_t(i,j)+p_b(i,j))
735 
736  eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
737  alpha_anom = al0 + lambda / (p0 + p_ave) - spv_ref
738  rem = lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
739  dza(i,j) = alpha_anom*dp + 2.0*eps*rem
740  if (present(intp_dza)) &
741  intp_dza(i,j) = 0.5*alpha_anom*dp**2 - dp*(1.0-eps)*rem
742  enddo ; enddo
743 
744  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
745  ! hWght is the distance measure by which the cell is violation of
746  ! hydrostatic consistency. For large hWght we bias the interpolation of
747  ! T & S along the top and bottom integrals, akin to thickness weighting.
748  hwght = 0.0
749  if (do_massweight) &
750  hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
751  if (hwght > 0.) then
752  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
753  hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
754  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
755  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
756  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
757  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
758  else
759  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
760  endif
761 
762  intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
763  do m=2,4
764  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
765  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
766 
767  ! T, S, and p are interpolated in the horizontal. The p interpolation
768  ! is linear, but for T and S it may be thickness wekghted.
769  al0 = wtt_l*al0_2d(i,j) + wtt_r*al0_2d(i+1,j)
770  p0 = wtt_l*p0_2d(i,j) + wtt_r*p0_2d(i+1,j)
771  lambda = wtt_l*lambda_2d(i,j) + wtt_r*lambda_2d(i+1,j)
772 
773  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
774  p_ave = 0.5*(wt_l*(p_t(i,j)+p_b(i,j)) + wt_r*(p_t(i+1,j)+p_b(i+1,j)))
775 
776  eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
777  intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
778  lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
779  enddo
780  ! Use Boole's rule to integrate the values.
781  intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
782  12.0*intp(3))
783  enddo ; enddo ; endif
784 
785  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
786  ! hWght is the distance measure by which the cell is violation of
787  ! hydrostatic consistency. For large hWght we bias the interpolation of
788  ! T & S along the top and bottom integrals, akin to thickness weighting.
789  hwght = 0.0
790  if (do_massweight) &
791  hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
792  if (hwght > 0.) then
793  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
794  hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
795  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
796  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
797  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
798  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
799  else
800  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
801  endif
802 
803  intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
804  do m=2,4
805  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
806  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
807 
808  ! T, S, and p are interpolated in the horizontal. The p interpolation
809  ! is linear, but for T and S it may be thickness wekghted.
810  al0 = wt_l*al0_2d(i,j) + wt_r*al0_2d(i,j+1)
811  p0 = wt_l*p0_2d(i,j) + wt_r*p0_2d(i,j+1)
812  lambda = wt_l*lambda_2d(i,j) + wt_r*lambda_2d(i,j+1)
813 
814  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
815  p_ave = 0.5*(wt_l*(p_t(i,j)+p_b(i,j)) + wt_r*(p_t(i,j+1)+p_b(i,j+1)))
816 
817  eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
818  intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
819  lambda * eps2 * (c1_3 + eps2*(0.2 + eps2*(c1_7 + c1_9*eps2)))
820  enddo
821  ! Use Boole's rule to integrate the values.
822  inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
823  12.0*intp(3))
824  enddo ; enddo ; endif