MOM6
MOM_EOS_linear.F90
1 !> A simple linear equation of state for sea water with constant coefficients
2 module mom_eos_linear
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_hor_index, only : hor_index_type
7 
8 implicit none ; private
9 
10 #include <MOM_memory.h>
11 
12 public calculate_compress_linear, calculate_density_linear, calculate_spec_vol_linear
13 public calculate_density_derivs_linear, calculate_density_derivs_scalar_linear
14 public calculate_specvol_derivs_linear
15 public calculate_density_scalar_linear, calculate_density_array_linear
16 public calculate_density_second_derivs_linear
17 public int_density_dz_linear, int_spec_vol_dp_linear
18 
19 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
20 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
21 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
22 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
23 
24 !> Compute the density of sea water (in kg/m^3), or its anomaly from a reference density,
25 !! using a simple linear equation of state from salinity (in psu), potential temperature (in deg C)
26 !! and pressure [Pa].
27 interface calculate_density_linear
28  module procedure calculate_density_scalar_linear, calculate_density_array_linear
29 end interface calculate_density_linear
30 
31 !> Compute the specific volume of sea water (in m^3/kg), or its anomaly from a reference value,
32 !! using a simple linear equation of state from salinity (in psu), potential temperature (in deg C)
33 !! and pressure [Pa].
34 interface calculate_spec_vol_linear
35  module procedure calculate_spec_vol_scalar_linear, calculate_spec_vol_array_linear
36 end interface calculate_spec_vol_linear
37 
38 !> For a given thermodynamic state, return the derivatives of density with temperature and
39 !! salinity using the simple linear equation of state
40 interface calculate_density_derivs_linear
41  module procedure calculate_density_derivs_scalar_linear, calculate_density_derivs_array_linear
42 end interface calculate_density_derivs_linear
43 
44 !> For a given thermodynamic state, return the second derivatives of density with various
45 !! combinations of temperature, salinity, and pressure. Note that with a simple linear
46 !! equation of state these second derivatives are all 0.
47 interface calculate_density_second_derivs_linear
48  module procedure calculate_density_second_derivs_scalar_linear, calculate_density_second_derivs_array_linear
49 end interface calculate_density_second_derivs_linear
50 
51 contains
52 
53 !> This subroutine computes the density of sea water with a trivial
54 !! linear equation of state (in [kg m-3]) from salinity (sal [PSU]),
55 !! potential temperature (T [degC]), and pressure [Pa].
56 subroutine calculate_density_scalar_linear(T, S, pressure, rho, &
57  Rho_T0_S0, dRho_dT, dRho_dS, rho_ref)
58  real, intent(in) :: t !< Potential temperature relative to the surface [degC].
59  real, intent(in) :: s !< Salinity [PSU].
60  real, intent(in) :: pressure !< pressure [Pa].
61  real, intent(out) :: rho !< In situ density [kg m-3].
62  real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
63  real, intent(in) :: drho_dt !< The derivatives of density with temperature
64  !! [kg m-3 degC-1].
65  real, intent(in) :: drho_ds !< The derivatives of density with salinity
66  !! in [kg m-3 ppt-1].
67  real, optional, intent(in) :: rho_ref !< A reference density [kg m-3].
68 
69  if (present(rho_ref)) then
70  rho = (rho_t0_s0 - rho_ref) + (drho_dt*t + drho_ds*s)
71  else
72  rho = rho_t0_s0 + drho_dt*t + drho_ds*s
73  endif
74 
75 end subroutine calculate_density_scalar_linear
76 
77 !> This subroutine computes the density of sea water with a trivial
78 !! linear equation of state (in kg/m^3) from salinity (sal in psu),
79 !! potential temperature (T [degC]), and pressure [Pa].
80 subroutine calculate_density_array_linear(T, S, pressure, rho, start, npts, &
81  Rho_T0_S0, dRho_dT, dRho_dS, rho_ref)
82  real, dimension(:), intent(in) :: t !< potential temperature relative to the surface [degC].
83  real, dimension(:), intent(in) :: s !< salinity [PSU].
84  real, dimension(:), intent(in) :: pressure !< pressure [Pa].
85  real, dimension(:), intent(out) :: rho !< in situ density [kg m-3].
86  integer, intent(in) :: start !< the starting point in the arrays.
87  integer, intent(in) :: npts !< the number of values to calculate.
88  real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
89  real, intent(in) :: drho_dt !< The derivatives of density with temperature
90  !! [kg m-3 degC-1].
91  real, intent(in) :: drho_ds !< The derivatives of density with salinity
92  !! in [kg m-3 ppt-1].
93  real, optional, intent(in) :: rho_ref !< A reference density [kg m-3].
94  ! Local variables
95  integer :: j
96 
97  if (present(rho_ref)) then ; do j=start,start+npts-1
98  rho(j) = (rho_t0_s0 - rho_ref) + (drho_dt*t(j) + drho_ds*s(j))
99  enddo ; else ; do j=start,start+npts-1
100  rho(j) = rho_t0_s0 + drho_dt*t(j) + drho_ds*s(j)
101  enddo ; endif
102 
103 end subroutine calculate_density_array_linear
104 
105 !> This subroutine computes the in situ specific volume of sea water (specvol in
106 !! [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC])
107 !! and pressure [Pa], using a trivial linear equation of state for density.
108 !! If spv_ref is present, specvol is an anomaly from spv_ref.
109 subroutine calculate_spec_vol_scalar_linear(T, S, pressure, specvol, &
110  Rho_T0_S0, dRho_dT, dRho_dS, spv_ref)
111  real, intent(in) :: T !< potential temperature relative to the surface
112  !! [degC].
113  real, intent(in) :: S !< Salinity [PSU].
114  real, intent(in) :: pressure !< Pressure [Pa].
115  real, intent(out) :: specvol !< In situ specific volume [m3 kg-1].
116  real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [kg m-3].
117  real, intent(in) :: dRho_dT !< The derivatives of density with temperature [kg m-3 degC-1].
118  real, intent(in) :: dRho_dS !< The derivatives of density with salinity [kg m-3 ppt-1].
119  real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
120  ! Local variables
121  integer :: j
122 
123  if (present(spv_ref)) then
124  specvol = ((1.0 - rho_t0_s0*spv_ref) + spv_ref*(drho_dt*t + drho_ds*s)) / &
125  ( rho_t0_s0 + (drho_dt*t + drho_ds*s))
126  else
127  specvol = 1.0 / ( rho_t0_s0 + (drho_dt*t + drho_ds*s))
128  endif
129 
130 end subroutine calculate_spec_vol_scalar_linear
131 
132 !> This subroutine computes the in situ specific volume of sea water (specvol in
133 !! [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC])
134 !! and pressure [Pa], using a trivial linear equation of state for density.
135 !! If spv_ref is present, specvol is an anomaly from spv_ref.
136 subroutine calculate_spec_vol_array_linear(T, S, pressure, specvol, start, npts, &
137  Rho_T0_S0, dRho_dT, dRho_dS, spv_ref)
138  real, dimension(:), intent(in) :: T !< potential temperature relative to the surface
139  !! [degC].
140  real, dimension(:), intent(in) :: S !< Salinity [PSU].
141  real, dimension(:), intent(in) :: pressure !< Pressure [Pa].
142  real, dimension(:), intent(out) :: specvol !< in situ specific volume [m3 kg-1].
143  integer, intent(in) :: start !< the starting point in the arrays.
144  integer, intent(in) :: npts !< the number of values to calculate.
145  real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [kg m-3].
146  real, intent(in) :: dRho_dT !< The derivatives of density with temperature [kg m-3 degC-1].
147  real, intent(in) :: dRho_dS !< The derivatives of density with salinity [kg m-3 ppt-1].
148  real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
149  ! Local variables
150  integer :: j
151 
152  if (present(spv_ref)) then ; do j=start,start+npts-1
153  specvol(j) = ((1.0 - rho_t0_s0*spv_ref) + spv_ref*(drho_dt*t(j) + drho_ds*s(j))) / &
154  ( rho_t0_s0 + (drho_dt*t(j) + drho_ds*s(j)))
155  enddo ; else ; do j=start,start+npts-1
156  specvol(j) = 1.0 / ( rho_t0_s0 + (drho_dt*t(j) + drho_ds*s(j)))
157  enddo ; endif
158 
159 end subroutine calculate_spec_vol_array_linear
160 
161 !> This subroutine calculates the partial derivatives of density *
162 !! with potential temperature and salinity.
163 subroutine calculate_density_derivs_array_linear(T, S, pressure, drho_dT_out, &
164  drho_dS_out, Rho_T0_S0, dRho_dT, dRho_dS, start, npts)
165  real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface
166  !! [degC].
167  real, intent(in), dimension(:) :: S !< Salinity [PSU].
168  real, intent(in), dimension(:) :: pressure !< Pressure [Pa].
169  real, intent(out), dimension(:) :: drho_dT_out !< The partial derivative of density with
170  !! potential temperature [kg m-3 degC-1].
171  real, intent(out), dimension(:) :: drho_dS_out !< The partial derivative of density with
172  !! salinity [kg m-3 ppt-1].
173  real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [kg m-3].
174  real, intent(in) :: dRho_dT !< The derivative of density with temperature [kg m-3 degC-1].
175  real, intent(in) :: dRho_dS !< The derivative of density with salinity [kg m-3 ppt-1].
176  integer, intent(in) :: start !< The starting point in the arrays.
177  integer, intent(in) :: npts !< The number of values to calculate.
178  ! Local variables
179  integer :: j
180 
181  do j=start,start+npts-1
182  drho_dt_out(j) = drho_dt
183  drho_ds_out(j) = drho_ds
184  enddo
185 
186 end subroutine calculate_density_derivs_array_linear
187 
188 !> This subroutine calculates the partial derivatives of density *
189 !! with potential temperature and salinity for a single point.
190 subroutine calculate_density_derivs_scalar_linear(T, S, pressure, drho_dT_out, &
191  drho_dS_out, Rho_T0_S0, dRho_dT, dRho_dS)
192  real, intent(in) :: t !< Potential temperature relative to the surface
193  !! [degC].
194  real, intent(in) :: s !< Salinity [PSU].
195  real, intent(in) :: pressure !< pressure [Pa].
196  real, intent(out) :: drho_dt_out !< The partial derivative of density with
197  !! potential temperature [kg m-3 degC-1].
198  real, intent(out) :: drho_ds_out !< The partial derivative of density with
199  !! salinity [kg m-3 ppt-1].
200  real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
201  real, intent(in) :: drho_dt !< The derivatives of density with temperature [kg m-3 degC-1].
202  real, intent(in) :: drho_ds !< The derivatives of density with salinity [kg m-3 ppt-1].
203  drho_dt_out = drho_dt
204  drho_ds_out = drho_ds
205 
206 end subroutine calculate_density_derivs_scalar_linear
207 
208 !> This subroutine calculates the five, partial second derivatives of density w.r.t.
209 !! potential temperature and salinity and pressure which for a linear equation of state should all be 0.
210 subroutine calculate_density_second_derivs_scalar_linear(T, S, pressure, drho_dS_dS, drho_dS_dT, &
211  drho_dT_dT, drho_dS_dP, drho_dT_dP)
212  real, intent(in) :: T !< Potential temperature relative to the surface [degC].
213  real, intent(in) :: S !< Salinity [PSU].
214  real, intent(in) :: pressure !< pressure [Pa].
215  real, intent(out) :: drho_dS_dS !< The second derivative of density with
216  !! salinity [kg m-3 PSU-2].
217  real, intent(out) :: drho_dS_dT !< The second derivative of density with
218  !! temperature and salinity [kg m-3 ppt-1 degC-1].
219  real, intent(out) :: drho_dT_dT !< The second derivative of density with
220  !! temperature [kg m-3 degC-2].
221  real, intent(out) :: drho_dS_dP !< The second derivative of density with
222  !! salinity and pressure [kg m-3 PSU-1 Pa-1].
223  real, intent(out) :: drho_dT_dP !< The second derivative of density with
224  !! temperature and pressure [kg m-3 degC-1 Pa-1].
225 
226  drho_ds_ds = 0.
227  drho_ds_dt = 0.
228  drho_dt_dt = 0.
229  drho_ds_dp = 0.
230  drho_dt_dp = 0.
231 
232 end subroutine calculate_density_second_derivs_scalar_linear
233 
234 !> This subroutine calculates the five, partial second derivatives of density w.r.t.
235 !! potential temperature and salinity and pressure which for a linear equation of state should all be 0.
236 subroutine calculate_density_second_derivs_array_linear(T, S,pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT,&
237  drho_dS_dP, drho_dT_dP, start, npts)
238  real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC].
239  real, dimension(:), intent(in) :: S !< Salinity [PSU].
240  real, dimension(:), intent(in) :: pressure !< pressure [Pa].
241  real, dimension(:), intent(out) :: drho_dS_dS !< The second derivative of density with
242  !! salinity [kg m-3 PSU-2].
243  real, dimension(:), intent(out) :: drho_dS_dT !< The second derivative of density with
244  !! temperature and salinity [kg m-3 ppt-1 degC-1].
245  real, dimension(:), intent(out) :: drho_dT_dT !< The second derivative of density with
246  !! temperature [kg m-3 degC-2].
247  real, dimension(:), intent(out) :: drho_dS_dP !< The second derivative of density with
248  !! salinity and pressure [kg m-3 PSU-1 Pa-1].
249  real, dimension(:), intent(out) :: drho_dT_dP !< The second derivative of density with
250  !! temperature and pressure [kg m-3 degC-1 Pa-1].
251  integer, intent(in) :: start !< The starting point in the arrays.
252  integer, intent(in) :: npts !< The number of values to calculate.
253  ! Local variables
254  integer :: j
255  do j=start,start+npts-1
256  drho_ds_ds(j) = 0.
257  drho_ds_dt(j) = 0.
258  drho_dt_dt(j) = 0.
259  drho_ds_dp(j) = 0.
260  drho_dt_dp(j) = 0.
261  enddo
262 
263 end subroutine calculate_density_second_derivs_array_linear
264 
265 !> Calculate the derivatives of specific volume with temperature and salinity
266 subroutine calculate_specvol_derivs_linear(T, S, pressure, dSV_dT, dSV_dS, &
267  start, npts, Rho_T0_S0, dRho_dT, dRho_dS)
268  real, intent(in), dimension(:) :: t !< Potential temperature relative to the surface
269  !! [degC].
270  real, intent(in), dimension(:) :: s !< Salinity [PSU].
271  real, intent(in), dimension(:) :: pressure !< pressure [Pa].
272  real, intent(out), dimension(:) :: dsv_ds !< The partial derivative of specific volume with
273  !! salinity [m3 kg-1 PSU-1].
274  real, intent(out), dimension(:) :: dsv_dt !< The partial derivative of specific volume with
275  !! potential temperature [m3 kg-1 degC-1].
276  integer, intent(in) :: start !< The starting point in the arrays.
277  integer, intent(in) :: npts !< The number of values to calculate.
278  real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
279  real, intent(in) :: drho_dt !< The derivative of density with
280  !! temperature, [kg m-3 degC-1].
281  real, intent(in) :: drho_ds !< The derivative of density with
282  !! salinity [kg m-3 ppt-1].
283  ! Local variables
284  real :: i_rho2
285  integer :: j
286 
287  do j=start,start+npts-1
288  ! Sv = 1.0 / (Rho_T0_S0 + dRho_dT*T(j) + dRho_dS*S(j))
289  i_rho2 = 1.0 / (rho_t0_s0 + (drho_dt*t(j) + drho_ds*s(j)))**2
290  dsv_dt(j) = -drho_dt * i_rho2
291  dsv_ds(j) = -drho_ds * i_rho2
292  enddo
293 
294 end subroutine calculate_specvol_derivs_linear
295 
296 !> This subroutine computes the in situ density of sea water (rho)
297 !! and the compressibility (drho/dp == C_sound^-2) at the given
298 !! salinity, potential temperature, and pressure.
299 subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts,&
300  Rho_T0_S0, dRho_dT, dRho_dS)
301  real, intent(in), dimension(:) :: t !< Potential temperature relative to the surface
302  !! [degC].
303  real, intent(in), dimension(:) :: s !< Salinity [PSU].
304  real, intent(in), dimension(:) :: pressure !< pressure [Pa].
305  real, intent(out), dimension(:) :: rho !< In situ density [kg m-3].
306  real, intent(out), dimension(:) :: drho_dp !< The partial derivative of density with pressure
307  !! (also the inverse of the square of sound speed)
308  !! [s2 m-2].
309  integer, intent(in) :: start !< The starting point in the arrays.
310  integer, intent(in) :: npts !< The number of values to calculate.
311  real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
312  real, intent(in) :: drho_dt !< The derivative of density with
313  !! temperature [kg m-3 degC-1].
314  real, intent(in) :: drho_ds !< The derivative of density with
315  !! salinity [kg m-3 ppt-1].
316  ! Local variables
317  integer :: j
318 
319  do j=start,start+npts-1
320  rho(j) = rho_t0_s0 + drho_dt*t(j) + drho_ds*s(j)
321  drho_dp(j) = 0.0
322  enddo
323 end subroutine calculate_compress_linear
324 
325 !> This subroutine calculates analytical and nearly-analytical integrals of
326 !! pressure anomalies across layers, which are required for calculating the
327 !! finite-volume form pressure accelerations in a Boussinesq model.
328 subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HII, HIO, &
329  Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa, &
330  bathyT, dz_neglect, useMassWghtInterp)
331  type(hor_index_type), intent(in) :: hii !< The horizontal index type for the input arrays.
332  type(hor_index_type), intent(in) :: hio !< The horizontal index type for the output arrays.
333  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
334  intent(in) :: t !< Potential temperature relative to the surface
335  !! [degC].
336  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
337  intent(in) :: s !< Salinity [PSU].
338  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
339  intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m].
340  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
341  intent(in) :: z_b !< Height at the top of the layer [Z ~> m].
342  real, intent(in) :: rho_ref !< A mean density [kg m-3], that is subtracted
343  !! out to reduce the magnitude of each of the
344  !! integrals.
345  real, intent(in) :: rho_0_pres !< A density [kg m-3], that is used to calculate
346  !! the pressure (as p~=-z*rho_0_pres*G_e) used in
347  !! the equation of state. rho_0_pres is not used
348  !! here.
349  real, intent(in) :: g_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
350  real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
351  real, intent(in) :: drho_dt !< The derivative of density with temperature,
352  !! [kg m-3 degC-1].
353  real, intent(in) :: drho_ds !< The derivative of density with salinity,
354  !! in [kg m-3 ppt-1].
355  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
356  intent(out) :: dpa !< The change in the pressure anomaly across the
357  !! layer [Pa].
358  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
359  optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer
360  !! of the pressure anomaly relative to the anomaly
361  !! at the top of the layer [Pa Z].
362  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
363  optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
364  !! pressure anomaly at the top and bottom of the
365  !! layer divided by the x grid spacing [Pa].
366  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
367  optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
368  !! pressure anomaly at the top and bottom of the
369  !! layer divided by the y grid spacing [Pa].
370  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
371  optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m].
372  real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m].
373  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
374  !! interpolate T/S for top and bottom integrals.
375  ! Local variables
376  real :: rho_anom ! The density anomaly from rho_ref [kg m-3].
377  real :: ral, rar ! rho_anom to the left and right [kg m-3].
378  real :: dz, dzl, dzr ! Layer thicknesses [Z ~> m].
379  real :: hwght ! A pressure-thickness below topography [Z ~> m].
380  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m].
381  real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2].
382  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
383  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
384  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
385  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
386  real :: intz(5) ! The integrals of density with height at the
387  ! 5 sub-column locations [Pa].
388  logical :: do_massweight ! Indicates whether to do mass weighting.
389  real, parameter :: c1_6 = 1.0/6.0, c1_90 = 1.0/90.0 ! Rational constants.
390  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, ioff, joff, m
391 
392  ioff = hio%idg_offset - hii%idg_offset
393  joff = hio%jdg_offset - hii%jdg_offset
394 
395  ! These array bounds work for the indexing convention of the input arrays, but
396  ! on the computational domain defined for the output arrays.
397  isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
398  jsq = hio%JscB + joff ; jeq = hio%JecB + joff
399  is = hio%isc + ioff ; ie = hio%iec + ioff
400  js = hio%jsc + joff ; je = hio%jec + joff
401 
402  do_massweight = .false.
403  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
404  do_massweight = .true.
405  ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//&
406  ! "bathyT must be present if useMassWghtInterp is present and true.")
407  ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//&
408  ! "dz_neglect must be present if useMassWghtInterp is present and true.")
409  endif ; endif
410 
411  do j=jsq,jeq+1 ; do i=isq,ieq+1
412  dz = z_t(i,j) - z_b(i,j)
413  rho_anom = (rho_t0_s0 - rho_ref) + drho_dt*t(i,j) + drho_ds*s(i,j)
414  dpa(i-ioff,j-joff) = g_e*rho_anom*dz
415  if (present(intz_dpa)) intz_dpa(i-ioff,j-joff) = 0.5*g_e*rho_anom*dz**2
416  enddo ; enddo
417 
418  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
419  ! hWght is the distance measure by which the cell is violation of
420  ! hydrostatic consistency. For large hWght we bias the interpolation of
421  ! T & S along the top and bottom integrals, akin to thickness weighting.
422  hwght = 0.0
423  if (do_massweight) &
424  hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
425 
426  if (hwght <= 0.0) then
427  dzl = z_t(i,j) - z_b(i,j) ; dzr = z_t(i+1,j) - z_b(i+1,j)
428  ral = (rho_t0_s0 - rho_ref) + (drho_dt*t(i,j) + drho_ds*s(i,j))
429  rar = (rho_t0_s0 - rho_ref) + (drho_dt*t(i+1,j) + drho_ds*s(i+1,j))
430 
431  intx_dpa(i-ioff,j-joff) = g_e*c1_6 * (dzl*(2.0*ral + rar) + dzr*(2.0*rar + ral))
432  else
433  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
434  hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
435  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
436  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
437  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
438  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
439 
440  intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
441  do m=2,4
442  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
443  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
444 
445  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
446  rho_anom = (rho_t0_s0 - rho_ref) + &
447  (drho_dt * (wtt_l*t(i,j) + wtt_r*t(i+1,j)) + &
448  drho_ds * (wtt_l*s(i,j) + wtt_r*s(i+1,j)))
449  intz(m) = g_e*rho_anom*dz
450  enddo
451  ! Use Bode's rule to integrate the values.
452  intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
453  12.0*intz(3))
454  endif
455  enddo ; enddo ; endif
456 
457  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
458  ! hWght is the distance measure by which the cell is violation of
459  ! hydrostatic consistency. For large hWght we bias the interpolation of
460  ! T & S along the top and bottom integrals, akin to thickness weighting.
461  hwght = 0.0
462  if (do_massweight) &
463  hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
464 
465  if (hwght <= 0.0) then
466  dzl = z_t(i,j) - z_b(i,j) ; dzr = z_t(i,j+1) - z_b(i,j+1)
467  ral = (rho_t0_s0 - rho_ref) + (drho_dt*t(i,j) + drho_ds*s(i,j))
468  rar = (rho_t0_s0 - rho_ref) + (drho_dt*t(i,j+1) + drho_ds*s(i,j+1))
469 
470  inty_dpa(i-ioff,j-joff) = g_e*c1_6 * (dzl*(2.0*ral + rar) + dzr*(2.0*rar + ral))
471  else
472  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
473  hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
474  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
475  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
476  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
477  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
478 
479  intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
480  do m=2,4
481  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
482  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
483 
484  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
485  rho_anom = (rho_t0_s0 - rho_ref) + &
486  (drho_dt * (wtt_l*t(i,j) + wtt_r*t(i,j+1)) + &
487  drho_ds * (wtt_l*s(i,j) + wtt_r*s(i,j+1)))
488  intz(m) = g_e*rho_anom*dz
489  enddo
490  ! Use Bode's rule to integrate the values.
491  inty_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
492  12.0*intz(3))
493  endif
494 
495  enddo ; enddo ; endif
496 end subroutine int_density_dz_linear
497 
498 !> Calculates analytical and nearly-analytical integrals in
499 !! pressure across layers of geopotential anomalies, which are required for
500 !! calculating the finite-volume form pressure accelerations in a non-Boussinesq
501 !! model. Specific volume is assumed to vary linearly between adjacent points.
502 subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, &
503  dRho_dT, dRho_dS, dza, intp_dza, intx_dza, inty_dza, halo_size, &
504  bathyP, dP_neglect, useMassWghtInterp)
505  type(hor_index_type), intent(in) :: hi !< The ocean's horizontal index type.
506  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
507  intent(in) :: t !< Potential temperature relative to the surface
508  !! [degC].
509  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
510  intent(in) :: s !< Salinity [PSU].
511  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
512  intent(in) :: p_t !< Pressure at the top of the layer [Pa].
513  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
514  intent(in) :: p_b !< Pressure at the top of the layer [Pa].
515  real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
516  !! to reduce the magnitude of each of the integrals, m3 kg-1. The calculation is
517  !! mathematically identical with different values of alpha_ref, but this reduces the
518  !! effects of roundoff.
519  real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
520  real, intent(in) :: drho_dt !< The derivative of density with temperature
521  !! [kg m-3 degC-1].
522  real, intent(in) :: drho_ds !< The derivative of density with salinity,
523  !! in [kg m-3 ppt-1].
524  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
525  intent(out) :: dza !< The change in the geopotential anomaly across
526  !! the layer [m2 s-2].
527  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
528  optional, intent(out) :: intp_dza !< The integral in pressure through the layer of
529  !! the geopotential anomaly relative to the anomaly
530  !! at the bottom of the layer [Pa m2 s-2].
531  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
532  optional, intent(out) :: intx_dza !< The integral in x of the difference between the
533  !! geopotential anomaly at the top and bottom of
534  !! the layer divided by the x grid spacing
535  !! [m2 s-2].
536  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
537  optional, intent(out) :: inty_dza !< The integral in y of the difference between the
538  !! geopotential anomaly at the top and bottom of
539  !! the layer divided by the y grid spacing
540  !! [m2 s-2].
541  integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
542  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
543  optional, intent(in) :: bathyp !< The pressure at the bathymetry [Pa]
544  real, optional, intent(in) :: dp_neglect !< A miniscule pressure change with
545  !! the same units as p_t [Pa]
546  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
547  !! to interpolate T/S for top and bottom integrals.
548  ! Local variables
549  real :: drho_ts ! The density anomaly due to T and S [kg m-3].
550  real :: alpha_anom ! The specific volume anomaly from 1/rho_ref [m3 kg-1].
551  real :: aal, aar ! rho_anom to the left and right [kg m-3].
552  real :: dp, dpl, dpr ! Layer pressure thicknesses [Pa].
553  real :: hwght ! A pressure-thickness below topography [Pa].
554  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Pa].
555  real :: idenom ! The inverse of the denominator in the weights [Pa-2].
556  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
557  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
558  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
559  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
560  real :: intp(5) ! The integrals of specific volume with pressure at the
561  ! 5 sub-column locations [m2 s-2].
562  logical :: do_massweight ! Indicates whether to do mass weighting.
563  real, parameter :: c1_6 = 1.0/6.0, c1_90 = 1.0/90.0 ! Rational constants.
564  integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, halo
565 
566  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
567  halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
568  ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
569  if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh); endif
570  if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh); endif
571 
572  do_massweight = .false.
573  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
574  do_massweight = .true.
575 ! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//&
576 ! "bathyP must be present if useMassWghtInterp is present and true.")
577 ! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//&
578 ! "dP_neglect must be present if useMassWghtInterp is present and true.")
579  endif ; endif
580 
581  do j=jsh,jeh ; do i=ish,ieh
582  dp = p_b(i,j) - p_t(i,j)
583  drho_ts = drho_dt*t(i,j) + drho_ds*s(i,j)
584  ! alpha_anom = 1.0/(Rho_T0_S0 + dRho_TS)) - alpha_ref
585  alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
586  dza(i,j) = alpha_anom*dp
587  if (present(intp_dza)) intp_dza(i,j) = 0.5*alpha_anom*dp**2
588  enddo ; enddo
589 
590  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
591  ! hWght is the distance measure by which the cell is violation of
592  ! hydrostatic consistency. For large hWght we bias the interpolation of
593  ! T & S along the top and bottom integrals, akin to thickness weighting.
594  hwght = 0.0
595  if (do_massweight) &
596  hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
597 
598  if (hwght <= 0.0) then
599  dpl = p_b(i,j) - p_t(i,j) ; dpr = p_b(i+1,j) - p_t(i+1,j)
600  drho_ts = drho_dt*t(i,j) + drho_ds*s(i,j)
601  aal = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
602  drho_ts = drho_dt*t(i+1,j) + drho_ds*s(i+1,j)
603  aar = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
604 
605  intx_dza(i,j) = c1_6 * (2.0*(dpl*aal + dpr*aar) + (dpl*aar + dpr*aal))
606  else
607  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
608  hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
609  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
610  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
611  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
612  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
613 
614  intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
615  do m=2,4
616  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
617  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
618 
619  ! T, S, and p are interpolated in the horizontal. The p interpolation
620  ! is linear, but for T and S it may be thickness wekghted.
621  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
622 
623  drho_ts = drho_dt*(wtt_l*t(i,j) + wtt_r*t(i+1,j)) + &
624  drho_ds*(wtt_l*s(i,j) + wtt_r*s(i+1,j))
625  ! alpha_anom = 1.0/(Rho_T0_S0 + dRho_TS)) - alpha_ref
626  alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
627  intp(m) = alpha_anom*dp
628  enddo
629  ! Use Bode's rule to integrate the interface height anomaly values in y.
630  intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
631  12.0*intp(3))
632  endif
633  enddo ; enddo ; endif
634 
635  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
636  ! hWght is the distance measure by which the cell is violation of
637  ! hydrostatic consistency. For large hWght we bias the interpolation of
638  ! T & S along the top and bottom integrals, akin to thickness weighting.
639  hwght = 0.0
640  if (do_massweight) &
641  hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
642 
643  if (hwght <= 0.0) then
644  dpl = p_b(i,j) - p_t(i,j) ; dpr = p_b(i,j+1) - p_t(i,j+1)
645  drho_ts = drho_dt*t(i,j) + drho_ds*s(i,j)
646  aal = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
647  drho_ts = drho_dt*t(i,j+1) + drho_ds*s(i,j+1)
648  aar = ((1.0 - rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
649 
650  inty_dza(i,j) = c1_6 * (2.0*(dpl*aal + dpr*aar) + (dpl*aar + dpr*aal))
651  else
652  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
653  hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
654  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
655  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
656  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
657  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
658 
659  intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
660  do m=2,4
661  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
662  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
663 
664  ! T, S, and p are interpolated in the horizontal. The p interpolation
665  ! is linear, but for T and S it may be thickness wekghted.
666  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
667 
668  drho_ts = drho_dt*(wtt_l*t(i,j) + wtt_r*t(i,j+1)) + &
669  drho_ds*(wtt_l*s(i,j) + wtt_r*s(i,j+1))
670  ! alpha_anom = 1.0/(Rho_T0_S0 + dRho_TS)) - alpha_ref
671  alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho_ts*alpha_ref) / (rho_t0_s0 + drho_ts)
672  intp(m) = alpha_anom*dp
673  enddo
674  ! Use Bode's rule to integrate the interface height anomaly values in y.
675  inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
676  12.0*intp(3))
677  endif
678  enddo ; enddo ; endif
679 end subroutine int_spec_vol_dp_linear
680 
681 end module mom_eos_linear
mom_eos_linear::calculate_density_second_derivs_linear
For a given thermodynamic state, return the second derivatives of density with various combinations o...
Definition: MOM_EOS_linear.F90:47
mom_eos_linear::calculate_density_derivs_linear
For a given thermodynamic state, return the derivatives of density with temperature and salinity usin...
Definition: MOM_EOS_linear.F90:40
mom_eos_linear
A simple linear equation of state for sea water with constant coefficients.
Definition: MOM_EOS_linear.F90:2
mom_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_eos_linear::calculate_spec_vol_linear
Compute the specific volume of sea water (in m^3/kg), or its anomaly from a reference value,...
Definition: MOM_EOS_linear.F90:34
mom_eos_linear::calculate_density_linear
Compute the density of sea water (in kg/m^3), or its anomaly from a reference density,...
Definition: MOM_EOS_linear.F90:27
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:15