MOM_EOS_linear.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> A simple linear equation of state for sea water with constant coefficients
7
10
11implicit none ; private
12
13public linear_eos
17
18!> The EOS_base implementation of a linear equation of state
19type, extends (eos_base) :: linear_eos
20
21 real :: rho_t0_s0 !< The density at T=0, S=0 and p=0 [kg m-3].
22 real :: drho_dt !< The derivative of density with temperature [kg m-3 degC-1].
23 real :: drho_ds !< The derivative of density with salinity [kg m-3 ppt-1].
24 real :: drho_dp !< The derivative of density with pressure [s2 m-2].
25
26contains
27 !> Implementation of the in-situ density as an elemental function [kg m-3]
28 procedure :: density_elem => density_elem_linear
29 !> Implementation of the in-situ density anomaly as an elemental function [kg m-3]
30 procedure :: density_anomaly_elem => density_anomaly_elem_linear
31 !> Implementation of the in-situ specific volume as an elemental function [m3 kg-1]
32 procedure :: spec_vol_elem => spec_vol_elem_linear
33 !> Implementation of the in-situ specific volume anomaly as an elemental function [m3 kg-1]
34 procedure :: spec_vol_anomaly_elem => spec_vol_anomaly_elem_linear
35 !> Implementation of the calculation of derivatives of density
36 procedure :: calculate_density_derivs_elem => calculate_density_derivs_elem_linear
37 !> Implementation of the calculation of second derivatives of density
38 procedure :: calculate_density_second_derivs_elem => calculate_density_second_derivs_elem_linear
39 !> Implementation of the calculation of derivatives of specific volume
40 procedure :: calculate_specvol_derivs_elem => calculate_specvol_derivs_elem_linear
41 !> Implementation of the calculation of compressibility
42 procedure :: calculate_compress_elem => calculate_compress_elem_linear
43 !> Implementation of the range query function
44 procedure :: eos_fit_range => eos_fit_range_linear
45
46 !> Instance specific function to set internal parameters
48
49 !> Local implementation of generic calculate_density_array for efficiency
50 procedure :: calculate_density_array => calculate_density_array_linear
51 !> Local implementation of generic calculate_spec_vol_array for efficiency
52 procedure :: calculate_spec_vol_array => calculate_spec_vol_array_linear
53
54end type linear_eos
55
56contains
57
58!> Density computed as a linear function of T and S [kg m-3]
59!!
60!! This is an elemental function that can be applied to any combination of
61!! scalar and array inputs.
62real elemental function density_elem_linear(this, t, s, pressure)
63 class(linear_eos), intent(in) :: this !< This EOS
64 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
65 real, intent(in) :: s !< Salinity [ppt]
66 real, intent(in) :: pressure !< Pressure [Pa]
67
68 density_elem_linear = this%Rho_T0_S0 + this%dRho_dT*t + this%dRho_dS*s + this%dRho_dp*pressure
69
70end function density_elem_linear
71
72!> Density anomaly computed as a linear function of T and S [kg m-3]
73!!
74!! This is an elemental function that can be applied to any combination of
75!! scalar and array inputs.
76real elemental function density_anomaly_elem_linear(this, t, s, pressure, rho_ref)
77 class(linear_eos), intent(in) :: this !< This EOS
78 real, intent(in) :: t !< Potential temperature relative to the surface [degC]
79 real, intent(in) :: s !< Salinity [ppt]
80 real, intent(in) :: pressure !< Pressure [Pa]
81 real, intent(in) :: rho_ref !< A reference density [kg m-3]
82
84 (this%Rho_T0_S0 - rho_ref) + ((this%dRho_dT*t + this%dRho_dS*s) + this%dRho_dp*pressure)
85
87
88!> Specific volume using a linear equation of state for density [m3 kg-1]
89!!
90!! This is an elemental function that can be applied to any combination of
91!! scalar and array inputs.
92real elemental function spec_vol_elem_linear(this, t, s, pressure)
93 class(linear_eos), intent(in) :: this !< This EOS
94 real, intent(in) :: t !< Potential temperature relative to the surface [degC].
95 real, intent(in) :: s !< Salinity [ppt].
96 real, intent(in) :: pressure !< Pressure [Pa].
97
99 1.0 / ( this%Rho_T0_S0 + ((this%dRho_dT*t + this%dRho_dS*s) + this%dRho_dp*pressure) )
100
101end function spec_vol_elem_linear
102
103!> Specific volume anomaly using a linear equation of state for density [m3 kg-1]
104!!
105!! This is an elemental function that can be applied to any combination of
106!! scalar and array inputs.
107real elemental function spec_vol_anomaly_elem_linear(this, t, s, pressure, spv_ref)
108 class(linear_eos), intent(in) :: this !< This EOS
109 real, intent(in) :: t !< Potential temperature relative to the surface [degC].
110 real, intent(in) :: s !< Salinity [ppt].
111 real, intent(in) :: pressure !< Pressure [Pa].
112 real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
113
115 ((1.0 - this%Rho_T0_S0*spv_ref) - &
116 spv_ref*((this%dRho_dT*t + this%dRho_dS*s) + this%dRho_dp*pressure)) / &
117 ( this%Rho_T0_S0 + ((this%dRho_dT*t + this%dRho_dS*s) + this%dRho_dp*pressure) )
118
120
121!> This subroutine calculates the partial derivatives of density
122!! with potential temperature and salinity.
123elemental subroutine calculate_density_derivs_elem_linear(this, T, S, pressure, dRho_dT, dRho_dS)
124 class(linear_eos), intent(in) :: this !< This EOS
125 real, intent(in) :: t !< Potential temperature relative to the surface [degC].
126 real, intent(in) :: s !< Salinity [ppt].
127 real, intent(in) :: pressure !< Pressure [Pa].
128 real, intent(out) :: drho_dt !< The partial derivative of density with
129 !! potential temperature [kg m-3 degC-1].
130 real, intent(out) :: drho_ds !< The partial derivative of density with
131 !! salinity [kg m-3 ppt-1].
132
133 drho_dt = this%dRho_dT
134 drho_ds = this%dRho_dS
135
137
138!> This subroutine calculates the five, partial second derivatives of density w.r.t.
139!! potential temperature and salinity and pressure which for a linear equation of state should all be 0.
140elemental subroutine calculate_density_second_derivs_elem_linear(this, T, S, pressure, &
141 drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP)
142 class(linear_eos), intent(in) :: this !< This EOS
143 real, intent(in) :: t !< Potential temperature relative to the surface [degC].
144 real, intent(in) :: s !< Salinity [ppt].
145 real, intent(in) :: pressure !< pressure [Pa].
146 real, intent(inout) :: drho_ds_ds !< The second derivative of density with
147 !! salinity [kg m-3 ppt-2].
148 real, intent(inout) :: drho_ds_dt !< The second derivative of density with
149 !! temperature and salinity [kg m-3 ppt-1 degC-1].
150 real, intent(inout) :: drho_dt_dt !< The second derivative of density with
151 !! temperature [kg m-3 degC-2].
152 real, intent(inout) :: drho_ds_dp !< The second derivative of density with
153 !! salinity and pressure [kg m-3 ppt-1 Pa-1].
154 real, intent(inout) :: drho_dt_dp !< The second derivative of density with
155 !! temperature and pressure [kg m-3 degC-1 Pa-1].
156
157 drho_ds_ds = 0.
158 drho_ds_dt = 0.
159 drho_dt_dt = 0.
160 drho_ds_dp = 0.
161 drho_dt_dp = 0.
162
164
165!> Calculate the derivatives of specific volume with temperature and salinity
166elemental subroutine calculate_specvol_derivs_elem_linear(this, T, S, pressure, dSV_dT, dSV_dS)
167 class(linear_eos), intent(in) :: this !< This EOS
168 real, intent(in) :: t !< Potential temperature [degC]
169 real, intent(in) :: s !< Salinity [ppt]
170 real, intent(in) :: pressure !< pressure [Pa]
171 real, intent(inout) :: dsv_ds !< The partial derivative of specific volume with
172 !! salinity [m3 kg-1 ppt-1]
173 real, intent(inout) :: dsv_dt !< The partial derivative of specific volume with
174 !! potential temperature [m3 kg-1 degC-1]
175 ! Local variables
176 real :: i_rho2 ! The inverse of density squared [m6 kg-2]
177
178 ! Sv = 1.0 / (Rho_T0_S0 + dRho_dT*T + dRho_dS*S)
179 i_rho2 = 1.0 / (this%Rho_T0_S0 + ((this%dRho_dT*t + this%dRho_dS*s) + this%dRho_dp*pressure))**2
180 dsv_dt = -this%dRho_dT * i_rho2
181 dsv_ds = -this%dRho_dS * i_rho2
182
184
185!> This subroutine computes the in situ density of sea water (rho)
186!! and the compressibility (drho/dp == C_sound^-2) at the given
187!! salinity, potential temperature, and pressure.
188elemental subroutine calculate_compress_elem_linear(this, T, S, pressure, rho, drho_dp)
189 class(linear_eos), intent(in) :: this !< This EOS
190 real, intent(in) :: t !< Potential temperature relative to the surface [degC].
191 real, intent(in) :: s !< Salinity [ppt].
192 real, intent(in) :: pressure !< pressure [Pa].
193 real, intent(out) :: rho !< In situ density [kg m-3].
194 real, intent(out) :: drho_dp !< The partial derivative of density with pressure
195 !! (also the inverse of the square of sound speed)
196 !! [s2 m-2].
197
198 rho = this%Rho_T0_S0 + this%dRho_dT*t + this%dRho_dS*s + this%dRho_dp*pressure
199 drho_dp = this%dRho_dp
200
202
203!> Calculates the layer average specific volumes. The analytical solution is
204!! SpV_avg = 1 / (drho_dp*dp) * ln[(1+eps)/(1-eps)] and the expression here is the first five terms of its
205!! Taylor series with a trunction error of O(eps**10). |eps|<0.02 for real ocean parameters.
206subroutine avg_spec_vol_linear(T, S, p_t, dp, SpV_avg, start, npts, Rho_T0_S0, dRho_dT, dRho_dS, dRho_dp)
207 real, dimension(:), intent(in) :: t !< Potential temperature [degC]
208 real, dimension(:), intent(in) :: s !< Salinity [ppt]
209 real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa]
210 real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa]
211 real, dimension(:), intent(inout) :: spv_avg !< The vertical average specific volume
212 !! in the layer [m3 kg-1]
213 integer, intent(in) :: start !< the starting point in the arrays.
214 integer, intent(in) :: npts !< the number of values to calculate.
215 real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3]
216 real, intent(in) :: drho_dt !< The derivative of density with temperature
217 !! [kg m-3 degC-1]
218 real, intent(in) :: drho_ds !< The derivative of density with salinity
219 !! [kg m-3 ppt-1]
220 real, intent(in) :: drho_dp !< The derivative of density with pressure
221 !! [s2 m-2]
222 ! Local variables
223 real :: eps2 ! The square of a nondimensional ratio [nondim]
224 real :: alpha_p_ave ! The specific volume at pressure mid-point [R-1 ~> m3 kg-1]
225 real, parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0, c1_9 = 1.0/9.0 ! Rational constants [nondim]
226 integer :: j
227
228 do j=start,start+npts-1
229 alpha_p_ave = &
230 1.0 / (rho_t0_s0 + ((drho_dt*t(j) + drho_ds*s(j)) + drho_dp*(p_t(j) + 0.5 * dp(j))))
231 eps2 = (0.5 * (drho_dp * dp(j)) * alpha_p_ave)**2
232 spv_avg(j) = alpha_p_ave * (1.0 + eps2 * (c1_3 + eps2 * (0.2 + eps2 * (c1_7 + c1_9 * eps2))))
233 enddo
234end subroutine avg_spec_vol_linear
235
236!> Return the range of temperatures, salinities and pressures permitted for linear equation of state.
237!! Care should be taken when applying this equation of state outside of its fit range.
238subroutine eos_fit_range_linear(this, T_min, T_max, S_min, S_max, p_min, p_max)
239 class(linear_eos), intent(in) :: this !< This EOS
240 real, optional, intent(out) :: T_min !< The minimum potential temperature over which this EoS is fitted [degC]
241 real, optional, intent(out) :: T_max !< The maximum potential temperature over which this EoS is fitted [degC]
242 real, optional, intent(out) :: S_min !< The minimum salinity over which this EoS is fitted [ppt]
243 real, optional, intent(out) :: S_max !< The maximum salinity over which this EoS is fitted [ppt]
244 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
245 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
246
247 if (present(t_min)) t_min = -273.0
248 if (present(t_max)) t_max = 100.0
249 if (present(s_min)) s_min = 0.0
250 if (present(s_max)) s_max = 1000.0
251 if (present(p_min)) p_min = 0.0
252 if (present(p_max)) p_max = 1.0e9
253
254end subroutine eos_fit_range_linear
255
256!> Set coefficients for the linear equation of state
257subroutine set_params_linear(this, Rho_T0_S0, dRho_dT, dRho_dS, dRho_dp)
258 class(linear_eos), intent(inout) :: this !< This EOS
259 real, optional, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [kg m-3]
260 real, optional, intent(in) :: dRho_dT !< The derivative of density with temperature,
261 !! [kg m-3 degC-1]
262 real, optional, intent(in) :: dRho_dS !< The derivative of density with salinity,
263 !! in [kg m-3 ppt-1]
264 real, optional, intent(in) :: dRho_dp !< The derivative of density with pressure,
265 !! in [s2 m-2]
266
267 if (present(rho_t0_s0)) this%Rho_T0_S0 = rho_t0_s0
268 if (present(drho_dt)) this%dRho_dT = drho_dt
269 if (present(drho_ds)) this%dRho_dS = drho_ds
270 if (present(drho_dp)) this%dRho_dp = drho_dp
271
272end subroutine set_params_linear
273
274!> This subroutine calculates analytical and nearly-analytical integrals of
275!! pressure anomalies across layers, which are required for calculating the
276!! finite-volume form pressure accelerations in a Boussinesq model.
277subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
278 Rho_T0_S0, dRho_dT, dRho_dS, dRho_dp, dpa, intz_dpa, intx_dpa, inty_dpa, &
279 bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p)
280 type(hor_index_type), intent(in) :: hi !< The horizontal index type for the arrays.
281 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
282 intent(in) :: t !< Potential temperature relative to the surface
283 !! [C ~> degC].
284 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
285 intent(in) :: s !< Salinity [S ~> ppt].
286 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
287 intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m].
288 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
289 intent(in) :: z_b !< Height at the top of the layer [Z ~> m].
290 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that
291 !! is subtracted out to reduce the magnitude of
292 !! each of the integrals.
293 real, intent(in) :: rho_0 !< A density [R ~> kg m-3], used to calculate
294 !! the pressure (as p~=-z*rho_0*G_e) used in
295 !! the equation of state.
296 real, intent(in) :: g_e !< The Earth's gravitational acceleration
297 !! [L2 Z-1 T-2 ~> m s-2]
298 real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [R ~> kg m-3]
299 real, intent(in) :: drho_dt !< The derivative of density with temperature,
300 !! [R C-1 ~> kg m-3 degC-1]
301 real, intent(in) :: drho_ds !< The derivative of density with salinity,
302 !! in [R S-1 ~> kg m-3 ppt-1]
303 real, intent(in) :: drho_dp !< The derivative of density with pressure,
304 !! in [L-2 T2 ~> m-2 s2]
305 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
306 intent(out) :: dpa !< The change in the pressure anomaly across the
307 !! layer [R L2 T-2 ~> Pa]
308 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
309 optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer
310 !! of the pressure anomaly relative to the anomaly
311 !! at the top of the layer [R L2 Z T-2 ~> Pa m]
312 real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
313 optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
314 !! pressure anomaly at the top and bottom of the
315 !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
316 real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
317 optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
318 !! pressure anomaly at the top and bottom of the
319 !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
320 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
321 optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m].
322 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
323 optional, intent(in) :: ssh !< The sea surface height [Z ~> m]
324 real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m].
325 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use
326 !! mass weighting to interpolate T/S in integrals
327 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
328 optional, intent(in) :: z_0p !< The height at which the pressure is 0 [Z ~> m]
329
330 ! Local variables
331 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: z0pres ! The height at which the pressure is zero [Z ~> m]
332 real :: rho_anom ! The density anomaly from rho_ref [R ~> kg m-3].
333 real :: ral, rar ! rho_anom to the left and right [R ~> kg m-3].
334 real :: dz, dzl, dzr ! Layer thicknesses [Z ~> m].
335 real :: gxrho ! The gravitational acceleration times mean ocean density [R L2 Z-1 T-2 ~> Pa m-1]
336 real :: p_ave ! The layer averaged pressure [R L2 T-2 ~> Pa]
337 real :: hwght ! A pressure-thickness below topography [Z ~> m].
338 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m].
339 real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2].
340 real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
341 real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
342 real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
343 real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
344 real :: intz(5) ! The integrals of density with height at the
345 ! 5 sub-column locations [R L2 T-2 ~> Pa]
346 logical :: do_massweight ! Indicates whether to do mass weighting.
347 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface
348 real, parameter :: c1_6 = 1.0/6.0, c1_90 = 1.0/90.0 ! Rational constants [nondim].
349 integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m
350
351 ! These array bounds work for the indexing convention of the input arrays, but
352 ! on the computational domain defined for the output arrays.
353 isq = hi%IscB ; ieq = hi%IecB
354 jsq = hi%JscB ; jeq = hi%JecB
355 is = hi%isc ; ie = hi%iec
356 js = hi%jsc ; je = hi%jec
357
358 gxrho = g_e * rho_0
359
360 if (present(z_0p)) then
361 do j=jsq,jeq+1 ; do i=isq,ieq+1
362 z0pres(i,j) = z_0p(i,j)
363 enddo ; enddo
364 else
365 z0pres(:,:) = 0.0
366 endif
367
368 do_massweight = .false. ; top_massweight = .false.
369 if (present(masswghtinterp)) then
370 do_massweight = btest(masswghtinterp, 0) ! True for odd values
371 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set
372 endif
373
374 do j=jsq,jeq+1 ; do i=isq,ieq+1
375 dz = z_t(i,j) - z_b(i,j)
376 p_ave = -gxrho * (0.5 * (z_t(i,j) + z_b(i,j)) - z0pres(i,j))
377 rho_anom = (rho_t0_s0 - rho_ref) + drho_dt * t(i,j) + drho_ds * s(i,j) + drho_dp * p_ave
378 dpa(i,j) = g_e * rho_anom * dz
379 if (present(intz_dpa)) &
380 intz_dpa(i,j) = 0.5 * g_e * (rho_anom - c1_6 * drho_dp * (gxrho * dz)) * dz**2
381 enddo ; enddo
382
383 if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
384 ! hWght is the distance measure by which the cell is violation of
385 ! hydrostatic consistency. For large hWght we bias the interpolation of
386 ! T & S along the top and bottom integrals, akin to thickness weighting.
387 hwght = 0.0
388 if (do_massweight) &
389 hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
390 if (top_massweight) &
391 hwght = max(hwght, z_b(i+1,j)-ssh(i,j), z_b(i,j)-ssh(i+1,j))
392
393 if (hwght <= 0.0) then
394 dzl = z_t(i,j) - z_b(i,j) ; dzr = z_t(i+1,j) - z_b(i+1,j)
395
396 p_ave = -gxrho * (0.5 * (z_t(i,j) + z_b(i,j)) - z0pres(i,j))
397 ral = (rho_t0_s0 - rho_ref) + ((drho_dt*t(i,j) + drho_ds*s(i,j)) + drho_dp*p_ave)
398
399 p_ave = -gxrho * (0.5 * (z_t(i+1,j) + z_b(i+1,j)) - z0pres(i+1,j))
400 rar = (rho_t0_s0 - rho_ref) + ((drho_dt*t(i+1,j) + drho_ds*s(i+1,j)) + drho_dp*p_ave)
401
402 intx_dpa(i,j) = g_e*c1_6 * ((dzl*(2.0*ral + rar)) + (dzr*(2.0*rar + ral)))
403 else
404 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
405 hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
406 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
407 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
408 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
409 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
410
411 intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
412 do m=2,4
413 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
414 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
415
416 dz = (wt_l*(z_t(i,j) - z_b(i,j))) + (wt_r*(z_t(i+1,j) - z_b(i+1,j)))
417 p_ave = -gxrho * ((wt_l * (0.5 * (z_t(i,j) + z_b(i,j)) - z0pres(i,j))) + &
418 (wt_r * (0.5 * (z_t(i+1,j) + z_b(i+1,j)) - z0pres(i+1,j))))
419 rho_anom = (rho_t0_s0 - rho_ref) + &
420 ((drho_dt * ((wtt_l*t(i,j)) + (wtt_r*t(i+1,j))) + &
421 drho_ds * ((wtt_l*s(i,j)) + (wtt_r*s(i+1,j)))) + drho_dp * p_ave)
422 intz(m) = g_e*rho_anom*dz
423 enddo
424 ! Use Boole's rule to integrate the values.
425 intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
426 12.0*intz(3))
427 endif
428 enddo ; enddo ; endif
429
430 if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
431 ! hWght is the distance measure by which the cell is violation of
432 ! hydrostatic consistency. For large hWght we bias the interpolation of
433 ! T & S along the top and bottom integrals, akin to thickness weighting.
434 hwght = 0.0
435 if (do_massweight) &
436 hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
437 if (top_massweight) &
438 hwght = max(hwght, z_b(i,j+1)-ssh(i,j), z_b(i,j)-ssh(i,j+1))
439
440 if (hwght <= 0.0) then
441 dzl = z_t(i,j) - z_b(i,j) ; dzr = z_t(i,j+1) - z_b(i,j+1)
442
443 p_ave = -gxrho * (0.5 * (z_t(i,j) + z_b(i,j)) - z0pres(i,j))
444 ral = (rho_t0_s0 - rho_ref) + ((drho_dt*t(i,j) + drho_ds*s(i,j)) + drho_dp*p_ave)
445
446 p_ave = -gxrho * (0.5 * (z_t(i,j+1) + z_b(i,j+1)) - z0pres(i,j+1))
447 rar = (rho_t0_s0 - rho_ref) + ((drho_dt*t(i,j+1) + drho_ds*s(i,j+1)) + drho_dp*p_ave)
448
449 inty_dpa(i,j) = g_e*c1_6 * ((dzl*(2.0*ral + rar)) + (dzr*(2.0*rar + ral)))
450 else
451 hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
452 hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
453 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
454 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
455 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
456 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
457
458 intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
459 do m=2,4
460 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
461 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
462
463 dz = (wt_l*(z_t(i,j) - z_b(i,j))) + (wt_r*(z_t(i,j+1) - z_b(i,j+1)))
464 p_ave = -gxrho * ((wt_l * (0.5 * (z_t(i,j) + z_b(i,j)) - z0pres(i,j))) + &
465 (wt_r * (0.5 * (z_t(i,j+1) + z_b(i,j+1)) - z0pres(i,j+1))))
466 rho_anom = (rho_t0_s0 - rho_ref) + &
467 ((drho_dt * ((wtt_l*t(i,j)) + (wtt_r*t(i,j+1))) + &
468 drho_ds * ((wtt_l*s(i,j)) + (wtt_r*s(i,j+1)))) + drho_dp * p_ave)
469 intz(m) = g_e*rho_anom*dz
470 enddo
471 ! Use Boole's rule to integrate the values.
472 inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
473 12.0*intz(3))
474 endif
475
476 enddo ; enddo ; endif
477end subroutine int_density_dz_linear
478
479!> Calculates analytical and nearly-analytical integrals in
480!! pressure across layers of geopotential anomalies, which are required for
481!! calculating the finite-volume form pressure accelerations in a non-Boussinesq
482!! model. Specific volume is assumed to vary linearly between adjacent points.
483subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, &
484 dRho_dT, dRho_dS, dRho_dp, dza, intp_dza, intx_dza, inty_dza, halo_size, &
485 bathyP, P_surf, dP_neglect, MassWghtInterp)
486 type(hor_index_type), intent(in) :: hi !< The ocean's horizontal index type.
487 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
488 intent(in) :: t !< Potential temperature relative to the surface
489 !! [C ~> degC].
490 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
491 intent(in) :: s !< Salinity [S ~> ppt].
492 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
493 intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
494 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
495 intent(in) :: p_b !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
496 real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
497 !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1].
498 !! The calculation is mathematically identical with different values of
499 !! alpha_ref, but this reduces the effects of roundoff.
500 real, intent(in) :: rho_t0_s0 !< The density at T=0, S=0 [R ~> kg m-3]
501 real, intent(in) :: drho_dt !< The derivative of density with temperature
502 !! [R C-1 ~> kg m-3 degC-1]
503 real, intent(in) :: drho_ds !< The derivative of density with salinity,
504 !! in [R S-1 ~> kg m-3 ppt-1]
505 real, intent(in) :: drho_dp !< The derivative of density with pressure,
506 !! in [L-2 T2 ~> m-2 s2]
507 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
508 intent(out) :: dza !< The change in the geopotential anomaly across
509 !! the layer [L2 T-2 ~> m2 s-2]
510 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
511 optional, intent(out) :: intp_dza !< The integral in pressure through the layer of the
512 !! geopotential anomaly relative to the anomaly at the
513 !! bottom of the layer [R L4 T-4 ~> Pa m2 s-2]
514 real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
515 optional, intent(out) :: intx_dza !< The integral in x of the difference between the
516 !! geopotential anomaly at the top and bottom of
517 !! the layer divided by the x grid spacing
518 !! [L2 T-2 ~> m2 s-2]
519 real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
520 optional, intent(out) :: inty_dza !< The integral in y of the difference between the
521 !! geopotential anomaly at the top and bottom of
522 !! the layer divided by the y grid spacing
523 !! [L2 T-2 ~> m2 s-2]
524 integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
525 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
526 optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
527 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
528 optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
529 real, optional, intent(in) :: dp_neglect !< A miniscule pressure change with
530 !! the same units as p_t [R L2 T-2 ~> Pa]
531 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use
532 !! mass weighting to interpolate T/S in integrals
533 ! Local variables
534 real :: drho ! The density anomaly due to T, S and p [R ~> kg m-3]
535 real :: lambda ! The sound speed squared [L2 T-2 ~> m2 s-2]
536 real :: eps, eps2 ! A nondimensional ratio and its square [nondim]
537 real :: rem ! [L2 T-2 ~> m2 s-2]
538 real :: p_ave ! The layer averaged pressure [R L2 T-2 ~> Pa]
539 real :: alpha_p_ave ! The specific volume at p_ave [R-1 ~> m3 kg-1]
540 real :: alpha_anom ! The specific volume anomaly from 1/rho_ref [R-1 ~> m3 kg-1]
541 real :: aal, aar ! The specific volume anomaly to the left and right [R-1 ~> m3 kg-1]
542 real :: dp, dpl, dpr ! Layer pressure thicknesses [R L2 T-2 ~> Pa]
543 real :: hwght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]
544 real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]
545 real :: idenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]
546 real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
547 real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
548 real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
549 real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
550 real :: intp(5) ! The integrals of specific volume with pressure at the
551 ! 5 sub-column locations [L2 T-2 ~> m2 s-2]
552 logical :: do_massweight ! Indicates whether to do mass weighting.
553 logical :: top_massweight ! Indicates whether to do mass weighting the sea surface
554 logical :: massweight_bug ! If true, use an incorrect expression to determine where to apply mass weighting
555 real, parameter :: c1_3 = 1.0/3.0, c1_7 = 1.0/7.0, c1_9 = 1.0/9.0 ! Rational constants [nondim]
556 real, parameter :: c1_6 = 1.0/6.0, c1_90 = 1.0/90.0 ! Rational constants [nondim].
557 integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, halo
558
559 isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
560 halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
561 ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
562 if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh) ; endif
563 if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh) ; endif
564
565 do_massweight = .false. ; massweight_bug = .false. ; top_massweight = .false.
566 if (present(masswghtinterp)) then
567 do_massweight = btest(masswghtinterp, 0) ! True for odd values
568 top_massweight = btest(masswghtinterp, 1) ! True if the 2 bit is set
569 massweight_bug = btest(masswghtinterp, 3) ! True if the 8 bit is set
570 endif
571
572 lambda = 0.0 ; if (drho_dp/=0.0) lambda = 1.0 / drho_dp
573 do j=jsh,jeh ; do i=ish,ieh
574 dp = p_b(i,j) - p_t(i,j)
575 p_ave = 0.5 * (p_t(i,j) + p_b(i,j))
576
577 drho = (drho_dt * t(i,j) + drho_ds * s(i,j)) + drho_dp * p_ave
578 alpha_p_ave = 1.0 / (rho_t0_s0 + drho)
579
580 ! A realistic upbound of eps is ~0.02, using dRho_dp ~ (1500 m/s)**(-2), alpha_p_ave ~ 1/(1030 kg/m3)
581 ! and dp ~ 1e8 Pa [~dz=10000m]. And if we use dp ~ 1e6 [~dz=100m], eps ~ 2e-4.
582 ! Analytically dza = 1/dRho_dp * ln[(1+eps)/(1-eps)] - alpha_ref * dp, and the expression here gives the first
583 ! five terms from its Taylor series with a truncation error of O(eps**11), which is beyond double floating
584 ! point precision.
585 eps = 0.5 * (drho_dp * dp) * alpha_p_ave ; eps2 = eps * eps
586 ! alpha_anom = 1.0/(Rho_T0_S0 + dRho) - alpha_ref
587 alpha_anom = ((1.0 - rho_t0_s0 * alpha_ref) - drho * alpha_ref) / (rho_t0_s0 + drho)
588 ! The following expression would be more efficient but I suspect it changes answer.
589 ! alpha_anom = ((1.0 - Rho_T0_S0 * alpha_ref) - drho * alpha_ref) * alpha_p_ave
590 rem = (lambda * eps2) * (c1_3 + eps2 * (0.2 + eps2 * (c1_7 + c1_9 * eps2)))
591 dza(i,j) = alpha_anom * dp + 2.0 * eps * rem
592 if (present(intp_dza)) &
593 intp_dza(i,j) = 0.5 * alpha_anom * dp**2 - dp * ((1.0 - eps) * rem)
594 enddo ; enddo
595
596 if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
597 ! hWght is the distance measure by which the cell is violation of
598 ! hydrostatic consistency. For large hWght we bias the interpolation of
599 ! T & S along the top and bottom integrals, akin to thickness weighting.
600 hwght = 0.0
601 if (do_massweight .and. massweight_bug) then
602 hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
603 elseif (do_massweight) then
604 hwght = max(0., p_t(i+1,j)-bathyp(i,j), p_t(i,j)-bathyp(i+1,j))
605 endif
606 if (top_massweight) &
607 hwght = max(hwght, p_surf(i,j)-p_b(i+1,j), p_surf(i+1,j)-p_b(i,j))
608
609 if (hwght <= 0.0) then
610 dpl = p_b(i,j) - p_t(i,j) ; dpr = p_b(i+1,j) - p_t(i+1,j)
611
612 p_ave = 0.5 * (p_b(i,j) + p_t(i,j))
613 drho = (drho_dt*t(i,j) + drho_ds*s(i,j)) + drho_dp * p_ave
614 aal = ((1.0 - rho_t0_s0*alpha_ref) - drho*alpha_ref) / (rho_t0_s0 + drho)
615
616 p_ave = 0.5 * (p_b(i+1,j) + p_t(i+1,j))
617 drho = (drho_dt*t(i+1,j) + drho_ds*s(i+1,j)) + drho_dp * p_ave
618 aar = ((1.0 - rho_t0_s0*alpha_ref) - drho*alpha_ref) / (rho_t0_s0 + drho)
619
620 intx_dza(i,j) = c1_6 * (2.0*((dpl*aal) + (dpr*aar)) + ((dpl*aar) + (dpr*aal)))
621 else
622 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
623 hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
624 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
625 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
626 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
627 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
628
629 intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
630 do m=2,4
631 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
632 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
633
634 ! T, S, and p are interpolated in the horizontal. The p interpolation
635 ! is linear, but for T and S it may be thickness weighted.
636 dp = (wt_l*(p_b(i,j) - p_t(i,j))) + (wt_r*(p_b(i+1,j) - p_t(i+1,j)))
637 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))))
638
639 drho = (drho_dt*((wtt_l*t(i,j)) + (wtt_r*t(i+1,j))) + &
640 drho_ds*((wtt_l*s(i,j)) + (wtt_r*s(i+1,j)))) + drho_dp * p_ave
641 ! alpha_anom = 1.0/(Rho_T0_S0 + drho)) - alpha_ref
642 alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho*alpha_ref) / (rho_t0_s0 + drho)
643 intp(m) = alpha_anom*dp
644 enddo
645 ! Use Boole's rule to integrate the interface height anomaly values in y.
646 intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
647 12.0*intp(3))
648 endif
649 enddo ; enddo ; endif
650
651 if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
652 ! hWght is the distance measure by which the cell is violation of
653 ! hydrostatic consistency. For large hWght we bias the interpolation of
654 ! T & S along the top and bottom integrals, akin to thickness weighting.
655 hwght = 0.0
656 if (do_massweight .and. massweight_bug) then
657 hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
658 elseif (do_massweight) then
659 hwght = max(0., p_t(i,j+1)-bathyp(i,j), p_t(i,j)-bathyp(i,j+1))
660 endif
661 if (top_massweight) &
662 hwght = max(hwght, p_surf(i,j)-p_b(i,j+1), p_surf(i,j+1)-p_b(i,j))
663
664 if (hwght <= 0.0) then
665 dpl = p_b(i,j) - p_t(i,j) ; dpr = p_b(i,j+1) - p_t(i,j+1)
666
667 p_ave = 0.5 * (p_b(i,j) + p_t(i,j)) + drho_dp * p_ave
668 drho = (drho_dt*t(i,j) + drho_ds*s(i,j)) + drho_dp * p_ave
669 aal = ((1.0 - rho_t0_s0*alpha_ref) - drho*alpha_ref) / (rho_t0_s0 + drho)
670
671 p_ave = 0.5 * (p_b(i,j+1) + p_t(i,j+1)) + drho_dp * p_ave
672 drho = (drho_dt*t(i,j+1) + drho_ds*s(i,j+1)) + drho_dp * p_ave
673 aar = ((1.0 - rho_t0_s0*alpha_ref) - drho*alpha_ref) / (rho_t0_s0 + drho)
674
675 inty_dza(i,j) = c1_6 * (2.0*((dpl*aal) + (dpr*aar)) + ((dpl*aar) + (dpr*aal)))
676 else
677 hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
678 hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
679 hwght = hwght * ( (hl-hr)/(hl+hr) )**2
680 idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
681 hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
682 hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
683
684 intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
685 do m=2,4
686 wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
687 wtt_l = (wt_l*hwt_ll) + (wt_r*hwt_rl) ; wtt_r = (wt_l*hwt_lr) + (wt_r*hwt_rr)
688
689 ! T, S, and p are interpolated in the horizontal. The p interpolation
690 ! is linear, but for T and S it may be thickness weighted.
691 dp = (wt_l*(p_b(i,j) - p_t(i,j))) + (wt_r*(p_b(i,j+1) - p_t(i,j+1)))
692 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))))
693
694 drho = (drho_dt*((wtt_l*t(i,j)) + (wtt_r*t(i,j+1))) + &
695 drho_ds*((wtt_l*s(i,j)) + (wtt_r*s(i,j+1)))) + drho_dp * p_ave
696 ! alpha_anom = 1.0/(Rho_T0_S0 + drho)) - alpha_ref
697 alpha_anom = ((1.0-rho_t0_s0*alpha_ref) - drho*alpha_ref) / (rho_t0_s0 + drho)
698 intp(m) = alpha_anom*dp
699 enddo
700 ! Use Boole's rule to integrate the interface height anomaly values in y.
701 inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
702 12.0*intp(3))
703 endif
704 enddo ; enddo ; endif
705end subroutine int_spec_vol_dp_linear
706
707!> Calculate the in-situ density for 1D arraya inputs and outputs.
708subroutine calculate_density_array_linear(this, T, S, pressure, rho, start, npts, rho_ref)
709 class(linear_eos), intent(in) :: this !< This EOS
710 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
711 real, dimension(:), intent(in) :: S !< Salinity [ppt]
712 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
713 real, dimension(:), intent(out) :: rho !< In situ density [kg m-3]
714 integer, intent(in) :: start !< The starting index for calculations
715 integer, intent(in) :: npts !< The number of values to calculate
716 real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]
717
718 ! Local variables
719 integer :: j
720
721 if (present(rho_ref)) then
722 do j = start, start+npts-1
723 rho(j) = density_anomaly_elem_linear(this, t(j), s(j), pressure(j), rho_ref)
724 enddo
725 else
726 do j = start, start+npts-1
727 rho(j) = density_elem_linear(this, t(j), s(j), pressure(j))
728 enddo
729 endif
730
732
733!> Calculate the in-situ specific volume for 1D array inputs and outputs.
734subroutine calculate_spec_vol_array_linear(this, T, S, pressure, specvol, start, npts, spv_ref)
735 class(linear_eos), intent(in) :: this !< This EOS
736 real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface [degC]
737 real, dimension(:), intent(in) :: S !< Salinity [ppt]
738 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
739 real, dimension(:), intent(out) :: specvol !< In situ specific volume [m3 kg-1]
740 integer, intent(in) :: start !< The starting index for calculations
741 integer, intent(in) :: npts !< The number of values to calculate
742 real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
743
744 ! Local variables
745 integer :: j
746
747 if (present(spv_ref)) then
748 do j = start, start+npts-1
749 specvol(j) = spec_vol_anomaly_elem_linear(this, t(j), s(j), pressure(j), spv_ref)
750 enddo
751 else
752 do j = start, start+npts-1
753 specvol(j) = spec_vol_elem_linear(this, t(j), s(j), pressure(j) )
754 enddo
755 endif
756
758
759end module mom_eos_linear