MOM_EOS.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!> Provides subroutines for quantities specific to the equation of state
6module mom_eos
7
18use mom_eos_unesco, only : unesco_eos
21use mom_eos_teos10, only : teos10_eos
22use mom_eos_teos10, only : gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp, gsw_ct_from_pt
24use mom_tfreeze, only : calculate_tfreeze_linear, calculate_tfreeze_millero
25use mom_tfreeze, only : calculate_tfreeze_teos10, calculate_tfreeze_teos_poly
26use mom_error_handler, only : mom_error, fatal, warning, mom_mesg
27use mom_file_parser, only : get_param, log_version, param_file_type
29use mom_io, only : stdout, stderr
32
33implicit none ; private
34
35public eos_domain
36public eos_init
37public eos_manual_init
38public eos_quadrature
39! public EOS_use_linear
40public eos_fit_range
41public eos_unit_tests
45public calculate_compress
47public calculate_density
48public calculate_density_derivs
49public calculate_density_second_derivs
50public calculate_spec_vol
51public calculate_specific_vol_derivs
52public calculate_tfreeze
58public gsw_sp_from_sr
59public gsw_sr_from_sp
60public gsw_pt_from_ct
62public get_eos_name
63
64! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
65! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
66! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
67! vary with the Boussinesq approximation, the Boussinesq variant is given first.
68
69!> Calculates density of sea water from T, S and P
70interface calculate_density
71 module procedure calculate_density_scalar
72 module procedure calculate_density_1d
74 module procedure calculate_stanley_density_1d
75end interface calculate_density
76
77!> Calculates specific volume of sea water from T, S and P
78interface calculate_spec_vol
79 module procedure calc_spec_vol_scalar
80 module procedure calc_spec_vol_1d
81end interface calculate_spec_vol
82
83!> Calculate the derivatives of density with temperature and salinity from T, S, and P
84interface calculate_density_derivs
86 module procedure calculate_density_derivs_1d
87end interface calculate_density_derivs
88
89!> Calculate the derivatives of specific volume with temperature and salinity from T, S, and P
90interface calculate_specific_vol_derivs
91 module procedure calc_spec_vol_derivs_1d
92end interface calculate_specific_vol_derivs
93
94!> Calculates the second derivatives of density with various combinations of temperature,
95!! salinity, and pressure from T, S and P
96interface calculate_density_second_derivs
98end interface calculate_density_second_derivs
99
100!> Calculates the freezing point of sea water from T, S and P
101interface calculate_tfreeze
103end interface calculate_tfreeze
104
105!> Calculates the compressibility of water from T, S, and P
106interface calculate_compress
108end interface calculate_compress
109
110!> A control structure for the equation of state
111type, public :: eos_type ; private
112 integer :: form_of_eos = 0 !< The equation of state to use.
113 integer :: form_of_tfreeze = 0 !< The expression for the potential temperature
114 !! of the freezing point.
115 logical :: eos_quadrature !< If true, always use the generic (quadrature)
116 !! code for the integrals of density.
117 logical :: compressible = .true. !< If true, in situ density is a function of pressure.
118! The following parameters are used with the linear equation of state only.
119 real :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3]
120 real :: drho_dt !< The partial derivative of density with temperature [kg m-3 degC-1]
121 real :: drho_ds !< The partial derivative of density with salinity [kg m-3 ppt-1]
122 real :: drho_dp !< The partial derivative of density with pressure [s2 m-2]
123! The following parameters are use with the linear expression for the freezing
124! point only.
125 real :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC]
126 real :: dtfr_ds !< The derivative of freezing point with salinity [degC ppt-1]
127 real :: dtfr_dp !< The derivative of freezing point with pressure [degC Pa-1]
128! The following are logicals pertaining to definitions of the thermodynamic state variables
129 logical :: use_cont_abss =.false. !< True if the model internal temperature is the conservative temperature and
130 !! the salinity is absolute salinity. These could be separated into two flags,
131 !! but right now it is controlled by one input parameter and there is no known
132 !! need to have one True and one False.
133 logical :: tfreeze_s_is_pracs =.true. !< True if the freezing point expression is formulated from practical salinity
134 logical :: tfreeze_t_is_pott = .true. !< True if the freezing point expression yields a potential temperature
135
136 logical :: use_wright_2nd_deriv_bug = .false. !< If true, use a separate subroutine that
137 !! retains a buggy version of the calculations of the second
138 !! derivative of density with temperature and with temperature and
139 !! pressure. This bug is corrected in the default version.
140
141! Unit conversion factors (normally used for dimensional testing but could also allow for
142! change of units of arguments to functions)
143 real :: m_to_z = 1. !< A constant that translates distances in meters to the units of depth [Z m-1 ~> 1]
144 real :: kg_m3_to_r = 1. !< A constant that translates kilograms per meter cubed to the
145 !! units of density [R m3 kg-1 ~> 1]
146 real :: r_to_kg_m3 = 1. !< A constant that translates the units of density to
147 !! kilograms per meter cubed [kg m-3 R-1 ~> 1]
148 real :: rl2_t2_to_pa = 1.!< Convert pressures from R L2 T-2 to Pa [Pa T2 R-1 L-2 ~> 1]
149 real :: l_t_to_m_s = 1. !< Convert lateral velocities from L T-1 to m s-1 [m T s-1 L-1 ~> 1]
150 real :: degc_to_c = 1. !< A constant that translates degrees Celsius to the units of temperature [C degC-1 ~> 1]
151 real :: c_to_degc = 1. !< A constant that translates the units of temperature to degrees Celsius [degC C-1 ~> 1]
152 real :: ppt_to_s = 1. !< A constant that translates parts per thousand to the units of salinity [S ppt-1 ~> 1]
153 real :: s_to_ppt = 1. !< A constant that translates the units of salinity to parts per thousand [ppt S-1 ~> 1]
154
155 !> The instance of the actual equation of state
156 class(eos_base), allocatable :: type
157
158end type eos_type
159
160! The named integers that might be stored in eqn_of_state_type%form_of_EOS.
161integer, parameter, public :: eos_linear = 1 !< A named integer specifying an equation of state
162integer, parameter, public :: eos_unesco = 2 !< A named integer specifying an equation of state
163integer, parameter, public :: eos_wright = 3 !< A named integer specifying an equation of state
164integer, parameter, public :: eos_wright_full = 4 !< A named integer specifying an equation of state
165integer, parameter, public :: eos_wright_reduced = 5 !< A named integer specifying an equation of state
166integer, parameter, public :: eos_teos10 = 6 !< A named integer specifying an equation of state
167integer, parameter, public :: eos_roquet_rho = 7 !< A named integer specifying an equation of state
168integer, parameter, public :: eos_roquet_spv = 8 !< A named integer specifying an equation of state
169integer, parameter, public :: eos_jackett06 = 9 !< A named integer specifying an equation of state
170!> A list of all the available EOS
171integer, dimension(9), public :: list_of_eos = (/ eos_linear, eos_unesco, &
172 eos_wright, eos_wright_full, eos_wright_reduced, &
173 eos_teos10, eos_roquet_rho, eos_roquet_spv, eos_jackett06 /)
174
175character*(12), parameter :: eos_linear_string = "LINEAR" !< A string for specifying the equation of state
176character*(12), parameter :: eos_unesco_string = "UNESCO" !< A string for specifying the equation of state
177character*(12), parameter :: eos_jackett_string = "JACKETT_MCD" !< A string for specifying the equation of state
178character*(12), parameter :: eos_wright_string = "WRIGHT" !< A string for specifying the equation of state
179character*(16), parameter :: eos_wright_red_string = "WRIGHT_REDUCED" !< A string for specifying the equation of state
180character*(12), parameter :: eos_wright_full_string = "WRIGHT_FULL" !< A string for specifying the equation of state
181character*(12), parameter :: eos_teos10_string = "TEOS10" !< A string for specifying the equation of state
182character*(12), parameter :: eos_nemo_string = "NEMO" !< A string for specifying the equation of state
183character*(12), parameter :: eos_roquet_rho_string = "ROQUET_RHO" !< A string for specifying the equation of state
184character*(12), parameter :: eos_roquet_spv_string = "ROQUET_SPV" !< A string for specifying the equation of state
185character*(12), parameter :: eos_jackett06_string = "JACKETT_06" !< A string for specifying the equation of state
186character*(12), parameter :: eos_default = eos_wright_full_string !< The default equation of state
187
188integer, parameter :: tfreeze_linear = 1 !< A named integer specifying a freezing point expression
189integer, parameter :: tfreeze_millero = 2 !< A named integer specifying a freezing point expression
190integer, parameter :: tfreeze_teos10 = 3 !< A named integer specifying a freezing point expression
191integer, parameter :: tfreeze_teospoly = 4 !< A named integer specifying a freezing point expression
192character*(10), parameter :: tfreeze_linear_string = "LINEAR" !< A string for specifying the freezing point expression
193character*(10), parameter :: tfreeze_millero_string = "MILLERO_78" !< A string for specifying the
194 !! freezing point expression
195character*(10), parameter :: tfreeze_teospoly_string = "TEOS_POLY" !< A string for specifying the
196 !! freezing point expression
197character*(10), parameter :: tfreeze_teos10_string = "TEOS10" !< A string for specifying the freezing point expression
198
199contains
200
201!> Density of sea water (in-situ if pressure is local) [R ~> kg m-3]
202!!
203!! If rho_ref is present, the anomaly with respect to rho_ref is returned. The pressure and
204!! density can be rescaled with the values stored in EOS. If the scale argument is present the density
205!! scaling uses the product of the two scaling factors.
206real elemental function calculate_density_elem(eos, t, s, pressure, rho_ref, scale)
207 type(eos_type), intent(in) :: eos !< Equation of state structure
208 real, intent(in) :: t !< Potential temperature referenced to the surface [C ~> degC]
209 real, intent(in) :: s !< Salinity [S ~> ppt]
210 real, intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
211 real, optional, intent(in) :: rho_ref !< A reference density [R ~> kg m-3]
212 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale output density in
213 !! combination with scaling stored in EOS [various]
214 real :: ta ! An array of temperatures [degC]
215 real :: sa ! An array of salinities [ppt]
216 real :: pres ! An mks version of the pressure to use [Pa]
217 real :: rho_mks ! An mks version of the density to be returned [kg m-3]
218 real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
219
220 pres = eos%RL2_T2_to_Pa * pressure
221 ta = eos%C_to_degC * t
222 sa = eos%S_to_ppt * s
223
224 if (present(rho_ref)) then
225 rho_mks = eos%type%density_anomaly_elem(ta, sa, pres, eos%R_to_kg_m3*rho_ref)
226 else
227 rho_mks = eos%type%density_elem(ta, sa, pres)
228 endif
229
230 ! Rescale the output density to the desired units.
231 rho_scale = eos%kg_m3_to_R
232 if (present(scale)) rho_scale = rho_scale * scale
233 calculate_density_elem = rho_scale * rho_mks
234
235end function calculate_density_elem
236
237!> Calls the appropriate subroutine to calculate density of sea water for scalar inputs.
238!! If rho_ref is present, the anomaly with respect to rho_ref is returned. The pressure and
239!! density can be rescaled with the values stored in EOS. If the scale argument is present the density
240!! scaling uses the product of the two scaling factors.
241subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, scale)
242 real, intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
243 real, intent(in) :: S !< Salinity [S ~> ppt]
244 real, intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
245 real, intent(out) :: rho !< Density (in-situ if pressure is local) [R ~> kg m-3]
246 type(eos_type), intent(in) :: EOS !< Equation of state structure
247 real, optional, intent(in) :: rho_ref !< A reference density [R ~> kg m-3]
248 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale output density in
249 !! combination with scaling stored in EOS [various]
250
251 real :: Ta ! An array of temperatures [degC]
252 real :: Sa ! An array of salinities [ppt]
253 real :: pres ! An mks version of the pressure to use [Pa]
254 real :: rho_mks ! An mks version of the density to be returned [kg m-3]
255 real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
256
257 pres = eos%RL2_T2_to_Pa * pressure
258 ta = eos%C_to_degC * t
259 sa = eos%S_to_ppt * s
260
261 if (present(rho_ref)) then
262 rho_mks = eos%type%density_anomaly_elem(ta, sa, pres, eos%R_to_kg_m3*rho_ref)
263 else
264 rho_mks = eos%type%density_elem(ta, sa, pres)
265 endif
266
267 ! Rescale the output density to the desired units.
268 rho_scale = eos%kg_m3_to_R
269 if (present(scale)) rho_scale = rho_scale * scale
270 rho = rho_scale * rho_mks
271
272end subroutine calculate_density_scalar
273
274!> Calls the appropriate subroutine to calculate density of sea water for scalar inputs
275!! including the variance of T, S and covariance of T-S.
276!! The calculation uses only the second order correction in a series as discussed
277!! in Stanley et al., 2020.
278!! If rho_ref is present, the anomaly with respect to rho_ref is returned. The
279!! density can be rescaled using rho_ref.
280subroutine calculate_stanley_density_scalar(T, S, pressure, Tvar, TScov, Svar, rho, EOS, rho_ref, scale)
281 real, intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
282 real, intent(in) :: S !< Salinity [S ~> ppt]
283 real, intent(in) :: Tvar !< Variance of potential temperature referenced to the surface [C2 ~> degC2]
284 real, intent(in) :: TScov !< Covariance of potential temperature and salinity [C S ~> degC ppt]
285 real, intent(in) :: Svar !< Variance of salinity [S2 ~> ppt2]
286 real, intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
287 real, intent(out) :: rho !< Density (in-situ if pressure is local) [R ~> kg m-3]
288 type(eos_type), intent(in) :: EOS !< Equation of state structure
289 real, optional, intent(in) :: rho_ref !< A reference density [R ~> kg m-3].
290 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale output density in
291 !! combination with scaling stored in EOS [various]
292 ! Local variables
293 real :: d2RdTT ! Second derivative of density with temperature [R C-2 ~> kg m-3 degC-2]
294 real :: d2RdST ! Second derivative of density with temperature and salinity [R S-1 C-1 ~> kg m-3 degC-1 ppt-1]
295 real :: d2RdSS ! Second derivative of density with salinity [R S-2 ~> kg m-3 ppt-2]
296 real :: d2RdSp ! Second derivative of density with salinity and pressure [T2 S-1 L-2 ~> kg m-3 ppt-1 Pa-1]
297 real :: d2RdTp ! Second derivative of density with temperature and pressure [T2 C-1 L-2 ~> kg m-3 degC-1 Pa-1]
298
299 call calculate_density_scalar(t, s, pressure, rho, eos, rho_ref)
300 call calculate_density_second_derivs_scalar(t, s, pressure, d2rdss, d2rdst, d2rdtt, d2rdsp, d2rdtp, eos)
301
302 ! Equation 25 of Stanley et al., 2020.
303 rho = rho + ( 0.5 * d2rdtt * tvar + ( d2rdst * tscov + 0.5 * d2rdss * svar ) )
304
305 if (present(scale)) rho = rho * scale
306
308
309!> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs,
310!! potentially limiting the domain of indices that are worked on.
311!! If rho_ref is present, the anomaly with respect to rho_ref is returned.
312subroutine calculate_density_1d(T, S, pressure, rho, EOS, dom, rho_ref, scale)
313 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
314 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
315 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
316 real, dimension(:), intent(inout) :: rho !< Density (in-situ if pressure is local) [R ~> kg m-3]
317 type(eos_type), intent(in) :: EOS !< Equation of state structure
318 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
319 !! into account that arrays start at 1.
320 real, optional, intent(in) :: rho_ref !< A reference density [R ~> kg m-3]
321 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
322 !! in combination with scaling stored in EOS [various]
323 ! Local variables
324 real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
325 real, dimension(size(rho)) :: pres ! Pressure converted to [Pa]
326 real, dimension(size(rho)) :: Ta ! Temperature converted to [degC]
327 real, dimension(size(rho)) :: Sa ! Salinity converted to [ppt]
328 integer :: i, is, ie, npts
329
330 if (present(dom)) then
331 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
332 else
333 is = 1 ; ie = size(rho) ; npts = 1 + ie - is
334 endif
335
336 if ((eos%RL2_T2_to_Pa == 1.0) .and. (eos%R_to_kg_m3 == 1.0) .and. &
337 (eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
338 call eos%type%calculate_density_array(t, s, pressure, rho, is, npts, rho_ref=rho_ref)
339 else ! This is the same as above, but with some extra work to rescale variables.
340 do i=is,ie
341 pres(i) = eos%RL2_T2_to_Pa * pressure(i)
342 ta(i) = eos%C_to_degC * t(i)
343 sa(i) = eos%S_to_ppt * s(i)
344 enddo
345 if (present(rho_ref)) then
346 call eos%type%calculate_density_array(ta, sa, pres, rho, is, npts, rho_ref=eos%R_to_kg_m3*rho_ref)
347 else
348 call eos%type%calculate_density_array(ta, sa, pres, rho, is, npts)
349 endif
350 endif
351
352 rho_scale = eos%kg_m3_to_R
353 if (present(scale)) rho_scale = rho_scale * scale
354 if (rho_scale /= 1.0) then ; do i=is,ie
355 rho(i) = rho_scale * rho(i)
356 enddo ; endif
357
358end subroutine calculate_density_1d
359
360!> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs
361!! including the variance of T, S and covariance of T-S,
362!! potentially limiting the domain of indices that are worked on.
363!! The calculation uses only the second order correction in a series as discussed
364!! in Stanley et al., 2020.
365!! If rho_ref is present, the anomaly with respect to rho_ref is returned.
366subroutine calculate_stanley_density_1d(T, S, pressure, Tvar, TScov, Svar, rho, EOS, dom, rho_ref, scale)
367 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
368 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
369 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
370 real, dimension(:), intent(in) :: Tvar !< Variance of potential temperature [C2 ~> degC2]
371 real, dimension(:), intent(in) :: TScov !< Covariance of potential temperature and salinity [C S ~> degC ppt]
372 real, dimension(:), intent(in) :: Svar !< Variance of salinity [S2 ~> ppt2]
373 real, dimension(:), intent(inout) :: rho !< Density (in-situ if pressure is local) [R ~> kg m-3]
374 type(eos_type), intent(in) :: EOS !< Equation of state structure
375 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
376 !! into account that arrays start at 1.
377 real, optional, intent(in) :: rho_ref !< A reference density [R ~> kg m-3]
378 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
379 !! in combination with scaling stored in EOS [various]
380 ! Local variables
381 real, dimension(size(T)) :: &
382 d2RdTT, & ! Second derivative of density with temperature [R C-2 ~> kg m-3 degC-2]
383 d2RdST, & ! Second derivative of density with temperature and salinity [R S-1 C-1 ~> kg m-3 degC-1 ppt-1]
384 d2RdSS, & ! Second derivative of density with salinity [R S-2 ~> kg m-3 ppt-2]
385 d2RdSp, & ! Second derivative of density with salinity and pressure [T2 S-1 L-2 ~> kg m-3 ppt-1 Pa-1]
386 d2RdTp ! Second derivative of density with temperature and pressure [T2 C-1 L-2 ~> kg m-3 degC-1 Pa-1]
387 integer :: i, is, ie, npts
388
389 if (present(dom)) then
390 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
391 else
392 is = 1 ; ie = size(rho) ; npts = 1 + ie - is
393 endif
394
395 call calculate_density_1d(t, s, pressure, rho, eos, dom, rho_ref)
396 call calculate_density_second_derivs_1d(t, s, pressure, d2rdss, d2rdst, d2rdtt, d2rdsp, d2rdtp, eos, dom)
397
398 ! Equation 25 of Stanley et al., 2020.
399 do i=is,ie
400 rho(i) = rho(i) + ( 0.5 * d2rdtt(i) * tvar(i) + ( d2rdst(i) * tscov(i) + 0.5 * d2rdss(i) * svar(i) ) )
401 enddo
402
403 if (present(scale)) then ; if (scale /= 1.0) then ; do i=is,ie
404 rho(i) = scale * rho(i)
405 enddo ; endif ; endif
406
407end subroutine calculate_stanley_density_1d
408
409!> Calls the appropriate subroutine to calculate the specific volume of sea water
410!! for 1-D array inputs.
411subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref, scale)
412 real, dimension(:), intent(in) :: T !< potential temperature relative to the surface [degC]
413 real, dimension(:), intent(in) :: S !< salinity [ppt]
414 real, dimension(:), intent(in) :: pressure !< pressure [Pa]
415 real, dimension(:), intent(inout) :: specvol !< in situ specific volume [kg m-3]
416 integer, intent(in) :: start !< the starting point in the arrays.
417 integer, intent(in) :: npts !< the number of values to calculate.
418 type(eos_type), intent(in) :: EOS !< Equation of state structure
419 real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]
420 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific
421 !! volume in combination with scaling stored in EOS [various]
422
423 integer :: j
424
425 if (.not. allocated(eos%type)) call mom_error(fatal, &
426 "calculate_spec_vol_array: EOS%form_of_EOS is not valid.")
427
428 call eos%type%calculate_spec_vol_array(t, s, pressure, specvol, start, npts, spv_ref)
429
430 if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
431 specvol(j) = scale * specvol(j)
432 enddo ; endif ; endif
433
434end subroutine calculate_spec_vol_array
435
436!> Calls the appropriate subroutine to calculate specific volume of sea water
437!! for scalar inputs.
438subroutine calc_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, scale)
439 real, intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
440 real, intent(in) :: S !< Salinity [S ~> ppt]
441 real, intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
442 real, intent(out) :: specvol !< In situ or potential specific volume [R-1 ~> m3 kg-1]
443 !! or other units determined by the scale argument
444 type(eos_type), intent(in) :: EOS !< Equation of state structure
445 real, optional, intent(in) :: spv_ref !< A reference specific volume [R-1 ~> m3 kg-1]
446 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific
447 !! volume in combination with scaling stored in EOS [various]
448
449 real, dimension(1) :: Ta ! Rescaled single element array version of temperature [degC]
450 real, dimension(1) :: Sa ! Rescaled single element array version of salinity [ppt]
451 real, dimension(1) :: pres ! Rescaled single element array version of pressure [Pa]
452 real, dimension(1) :: spv ! Rescaled single element array version of specific volume [m3 kg-1]
453 real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg R-1 m-3 ~> 1]
454
455 pres(1) = eos%RL2_T2_to_Pa * pressure
456 ta(1) = eos%C_to_degC * t ; sa(1) = eos%S_to_ppt * s
457
458 if (present(spv_ref)) then
459 call calculate_spec_vol_array(ta, sa, pres, spv, 1, 1, eos, eos%kg_m3_to_R*spv_ref)
460 else
461 call calculate_spec_vol_array(ta, sa, pres, spv, 1, 1, eos)
462 endif
463 specvol = spv(1)
464
465 spv_scale = eos%R_to_kg_m3
466 if (present(scale)) spv_scale = spv_scale * scale
467 if (spv_scale /= 1.0) then
468 specvol = spv_scale * specvol
469 endif
470
471end subroutine calc_spec_vol_scalar
472
473!> Calls the appropriate subroutine to calculate the specific volume of sea water for 1-D array
474!! inputs, potentially limiting the domain of indices that are worked on.
475subroutine calc_spec_vol_1d(T, S, pressure, specvol, EOS, dom, spv_ref, scale)
476 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
477 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
478 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
479 real, dimension(:), intent(inout) :: specvol !< In situ specific volume [R-1 ~> m3 kg-1]
480 type(eos_type), intent(in) :: EOS !< Equation of state structure
481 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
482 !! into account that arrays start at 1.
483 real, optional, intent(in) :: spv_ref !< A reference specific volume [R-1 ~> m3 kg-1]
484 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale
485 !! output specific volume in combination with
486 !! scaling stored in EOS [various]
487 ! Local variables
488 real, dimension(size(T)) :: pres ! Pressure converted to [Pa]
489 real, dimension(size(T)) :: Ta ! Temperature converted to [degC]
490 real, dimension(size(T)) :: Sa ! Salinity converted to [ppt]
491 real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
492 integer :: i, is, ie, npts
493
494 if (present(dom)) then
495 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
496 else
497 is = 1 ; ie = size(specvol) ; npts = 1 + ie - is
498 endif
499
500 if ((eos%RL2_T2_to_Pa == 1.0) .and. (eos%kg_m3_to_R == 1.0) .and. &
501 (eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
502 call calculate_spec_vol_array(t, s, pressure, specvol, is, npts, eos, spv_ref)
503 else ! This is the same as above, but with some extra work to rescale variables.
504 do i=is,ie
505 pres(i) = eos%RL2_T2_to_Pa * pressure(i)
506 ta(i) = eos%C_to_degC * t(i)
507 sa(i) = eos%S_to_ppt * s(i)
508 enddo
509 if (present(spv_ref)) then
510 call calculate_spec_vol_array(ta, sa, pres, specvol, is, npts, eos, eos%kg_m3_to_R*spv_ref)
511 else
512 ! There is rescaling of variables, but spv_ref is not present. Passing a 0 value of spv_ref
513 ! changes answers at roundoff for some equations of state, like Wright and UNESCO.
514 call calculate_spec_vol_array(ta, sa, pres, specvol, is, npts, eos)
515 endif
516 endif
517
518 spv_scale = eos%R_to_kg_m3
519 if (present(scale)) spv_scale = spv_scale * scale
520 if (spv_scale /= 1.0) then ; do i=is,ie
521 specvol(i) = spv_scale * specvol(i)
522 enddo ; endif
523
524end subroutine calc_spec_vol_1d
525
526
527!> Calls the appropriate subroutine to calculate the freezing point for scalar inputs.
528subroutine calculate_tfreeze_scalar(S, pressure, T_fr, EOS, pres_scale, scale_from_EOS)
529 real, intent(in) :: S !< Salinity, [ppt] or [S ~> ppt] depending on scale_from_EOS
530 real, intent(in) :: pressure !< Pressure, in [Pa] or [R L2 T-2 ~> Pa] depending on
531 !! pres_scale or scale_from_EOS
532 real, intent(out) :: T_fr !< Freezing point potential temperature referenced to the
533 !! surface [degC] or [C ~> degC] depending on scale_from_EOS
534 type(eos_type), intent(in) :: EOS !< Equation of state structure
535 real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
536 !! into Pa [Pa T2 R-1 L-2 ~> 1].
537 logical, optional, intent(in) :: scale_from_EOS !< If present true use the dimensional scaling
538 !! factors stored in EOS. Omission is the same .false.
539
540 ! Local variables
541 real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1]
542 real :: S_scale ! A factor to convert salinity to units of ppt [ppt S-1 ~> 1]
543 real :: iS_scale! A factor to convert salinity to units of S [S ppt-1 ~> 1]
544 real :: absS ! A salinity converted to absolute salinity, only used in specific scenarios [ppt]
545 real :: TFreeze_S ! The salinity for the freezing equation in model units [S ~> PSU or ppt]
546
547 p_scale = 1.0 ; s_scale = 1.0 ; is_scale = 1.0
548 if (present(pres_scale)) p_scale = pres_scale
549 if (present(scale_from_eos)) then ; if (scale_from_eos) then
550 p_scale = eos%RL2_T2_to_Pa
551 s_scale = eos%S_to_ppt
552 is_scale = eos%ppt_to_S
553 endif ; endif
554
555 if (eos%use_conT_absS) then
556 ! Otherwise absS is unneeded and therefore unset
557 abss = s*s_scale
558 if (eos%TFreeze_S_is_pracS) then
559 tfreeze_s = gsw_sp_from_sr(abss)*is_scale
560 else
561 tfreeze_s = s
562 endif
563 else
564 tfreeze_s = s
565 endif
566
567 select case (eos%form_of_TFreeze)
568 case (tfreeze_linear)
569 call calculate_tfreeze_linear(s_scale*tfreeze_s, p_scale*pressure, t_fr, eos%TFr_S0_P0, &
570 eos%dTFr_dS, eos%dTFr_dp)
571 case (tfreeze_millero)
572 call calculate_tfreeze_millero(s_scale*tfreeze_s, p_scale*pressure, t_fr)
573 case (tfreeze_teospoly)
574 call calculate_tfreeze_teos_poly(s_scale*tfreeze_s, p_scale*pressure, t_fr)
575 case (tfreeze_teos10)
576 call calculate_tfreeze_teos10(s_scale*tfreeze_s, p_scale*pressure, t_fr)
577 case default
578 call mom_error(fatal, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
579 end select
580
581 if (eos%use_conT_absS .and. eos%TFreeze_T_is_potT) then
582 ! absS is set only if EOS%use_conT_absS is True
583 ! absS and T_fr have physical units here and don't need converted
584 t_fr = gsw_ct_from_pt(abss,t_fr)
585 endif
586
587 if (present(scale_from_eos)) then ; if (scale_from_eos) then
588 t_fr = eos%degC_to_C * t_fr
589 endif ; endif
590
591end subroutine calculate_tfreeze_scalar
592
593!> Calls the appropriate subroutine to calculate the freezing point for a 1-D array.
594subroutine calculate_tfreeze_array(S, pressure, T_fr, start, npts, EOS, pres_scale)
595 real, dimension(:), intent(in) :: S !< Salinity [ppt]
596 real, dimension(:), intent(in) :: pressure !< Pressure, in [Pa] or [R L2 T-2 ~> Pa] depending on pres_scale
597 real, dimension(:), intent(inout) :: T_fr !< Freezing point, either potential temperature referenced to the
598 !! surface or conservative temperature depending on settings [degC]
599 integer, intent(in) :: start !< Starting index within the array
600 integer, intent(in) :: npts !< The number of values to calculate
601 type(eos_type), intent(in) :: EOS !< Equation of state structure
602 real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure
603 !! into Pa [Pa T2 R-1 L-2 ~> 1].
604
605 ! Local variables
606 real, dimension(size(pressure)) :: pres ! Pressure converted to [Pa]
607 real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1]
608 real, dimension(size(S)) :: absS ! A salinity converted to absolute salinity, only used in specific scenarios [ppt]
609 real, dimension(size(S)) :: TFreeze_S ! The salinity for the freezing equation in model units [S ~> PSU or ppt]
610 integer :: j
611
612 p_scale = 1.0 ; if (present(pres_scale)) p_scale = pres_scale
613
614 if (eos%use_conT_absS) then
615 ! Otherwise absS is unneeded and therefore unset
616 abss(:) = s(:)
617 if (eos%TFreeze_S_is_pracS) then
618 tfreeze_s(:) = gsw_sp_from_sr(abss(:))
619 else
620 tfreeze_s(:) = s(:)
621 endif
622 else
623 tfreeze_s(:) = s(:)
624 endif
625
626 if (p_scale == 1.0) then
627 select case (eos%form_of_TFreeze)
628 case (tfreeze_linear)
629 call calculate_tfreeze_linear(tfreeze_s, pressure, t_fr, start, npts, &
630 eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
631 case (tfreeze_millero)
632 call calculate_tfreeze_millero(tfreeze_s, pressure, t_fr, start, npts)
633 case (tfreeze_teospoly)
634 call calculate_tfreeze_teos_poly(tfreeze_s, pressure, t_fr, start, npts)
635 case (tfreeze_teos10)
636 call calculate_tfreeze_teos10(tfreeze_s, pressure, t_fr, start, npts)
637 case default
638 call mom_error(fatal, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
639 end select
640 else
641 do j=start,start+npts-1 ; pres(j) = p_scale * pressure(j) ; enddo
642 select case (eos%form_of_TFreeze)
643 case (tfreeze_linear)
644 call calculate_tfreeze_linear(tfreeze_s, pres, t_fr, start, npts, &
645 eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
646 case (tfreeze_millero)
647 call calculate_tfreeze_millero(tfreeze_s, pres, t_fr, start, npts)
648 case (tfreeze_teos10)
649 call calculate_tfreeze_teos10(tfreeze_s, pres, t_fr, start, npts)
650 case (tfreeze_teospoly)
651 call calculate_tfreeze_teos_poly(tfreeze_s, pres, t_fr, start, npts)
652 case default
653 call mom_error(fatal, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
654 end select
655 endif
656
657 if (eos%use_conT_absS .and. eos%TFreeze_T_is_potT) then
658 ! absS is set only if EOS%use_conT_absS is True!
659 t_fr(:) = gsw_ct_from_pt(abss(:),t_fr(:))
660 endif
661
662
663end subroutine calculate_tfreeze_array
664
665!> Calls the appropriate subroutine to calculate the freezing point for a 1-D array, taking
666!! dimensionally rescaled arguments with factors stored in EOS.
667subroutine calculate_tfreeze_1d(S, pressure, T_fr, EOS, dom)
668 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
669 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
670 real, dimension(:), intent(inout) :: T_fr !< Freezing point, either potential temperature referenced to the
671 !! surface or conservative temperature depending on settings
672 !! [C ~> degC]
673 type(eos_type), intent(in) :: EOS !< Equation of state structure
674 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
675 !! into account that arrays start at 1.
676
677 ! Local variables
678 real, dimension(size(T_fr)) :: pres ! Pressure converted to [Pa]
679 real, dimension(size(T_fr)) :: Sa ! Salinity converted to [ppt]
680 real, dimension(size(T_fr)) :: absS ! Salinity converted to absoluate salinity [ppt]
681 real, dimension(size(T_fr)) :: TFreeze_S ! The salinity for the freezing equation in model units [S ~> PSU or ppt]
682 integer :: i, is, ie, npts
683
684 if (present(dom)) then
685 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
686 else
687 is = 1 ; ie = size(t_fr) ; npts = 1 + ie - is
688 endif
689
690 if (eos%use_conT_absS) then
691 ! Otherwise absS is unneeded and therefore unset
692 abss(:) = s(:)*eos%S_to_ppt
693 if (eos%TFreeze_S_is_pracS) then
694 tfreeze_s(:) = gsw_sp_from_sr(abss(:))*eos%ppt_to_S
695 else
696 tfreeze_s(:) = s(:)
697 endif
698 else
699 tfreeze_s(:) = s(:)
700 endif
701
702 if ((eos%RL2_T2_to_Pa == 1.0) .and. (eos%S_to_ppt == 1.0)) then
703 select case (eos%form_of_TFreeze)
704 case (tfreeze_linear)
705 call calculate_tfreeze_linear(tfreeze_s, pressure, t_fr, is, npts, &
706 eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
707 case (tfreeze_millero)
708 call calculate_tfreeze_millero(tfreeze_s, pressure, t_fr, is, npts)
709 case (tfreeze_teospoly)
710 call calculate_tfreeze_teos_poly(tfreeze_s, pressure, t_fr, is, npts)
711 case (tfreeze_teos10)
712 call calculate_tfreeze_teos10(tfreeze_s, pressure, t_fr, is, npts)
713 case default
714 call mom_error(fatal, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
715 end select
716 else
717 do i=is,ie
718 pres(i) = eos%RL2_T2_to_Pa * pressure(i)
719 sa(i) = eos%S_to_ppt * tfreeze_s(i)
720 enddo
721 select case (eos%form_of_TFreeze)
722 case (tfreeze_linear)
723 call calculate_tfreeze_linear(sa, pres, t_fr, is, npts, &
724 eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
725 case (tfreeze_millero)
726 call calculate_tfreeze_millero(sa, pres, t_fr, is, npts)
727 case (tfreeze_teospoly)
728 call calculate_tfreeze_teos_poly(sa, pres, t_fr, is, npts)
729 case (tfreeze_teos10)
730 call calculate_tfreeze_teos10(sa, pres, t_fr, is, npts)
731 case default
732 call mom_error(fatal, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
733 end select
734 endif
735
736 if (eos%use_conT_absS .and. eos%TFreeze_T_is_potT) then
737 ! absS is set only if EOS%use_conT_absS is True!
738 ! absS is in ppt and T_fr is in degC at this point.
739 t_fr(:) = gsw_ct_from_pt(abss(:),t_fr(:))
740 endif
741
742
743 if (eos%degC_to_C /= 1.0) then
744 do i=is,ie ; t_fr(i) = eos%degC_to_C * t_fr(i) ; enddo
745 endif
746
747end subroutine calculate_tfreeze_1d
748
749
750!> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.
751subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS, scale)
752 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC]
753 real, dimension(:), intent(in) :: S !< Salinity [ppt]
754 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
755 real, dimension(:), intent(inout) :: drho_dT !< The partial derivative of density with potential
756 !! temperature [kg m-3 degC-1] or other units determined
757 !! by the optional scale argument
758 real, dimension(:), intent(inout) :: drho_dS !< The partial derivative of density with salinity,
759 !! in [kg m-3 ppt-1] or other units determined
760 !! by the optional scale argument
761 integer, intent(in) :: start !< Starting index within the array
762 integer, intent(in) :: npts !< The number of values to calculate
763 type(eos_type), intent(in) :: EOS !< Equation of state structure
764 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
765 !! in combination with scaling stored in EOS [various]
766
767 ! Local variables
768 integer :: j
769
770 if (.not. allocated(eos%type)) call mom_error(fatal, &
771 "calculate_density_derivs_array: EOS%form_of_EOS is not valid.")
772
773 call eos%type%calculate_density_derivs_array(t, s, pressure, drho_dt, drho_ds, start, npts)
774
775 if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
776 drho_dt(j) = scale * drho_dt(j)
777 drho_ds(j) = scale * drho_ds(j)
778 enddo ; endif ; endif
779
781
782
783!> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.
784subroutine calculate_density_derivs_1d(T, S, pressure, drho_dT, drho_dS, EOS, dom, scale)
785 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
786 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
787 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
788 real, dimension(:), intent(inout) :: drho_dT !< The partial derivative of density with potential
789 !! temperature [R C-1 ~> kg m-3 degC-1]
790 real, dimension(:), intent(inout) :: drho_dS !< The partial derivative of density with salinity
791 !! [R S-1 ~> kg m-3 ppt-1]
792 type(eos_type), intent(in) :: EOS !< Equation of state structure
793 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
794 !! into account that arrays start at 1.
795 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
796 !! in combination with scaling stored in EOS [various]
797 ! Local variables
798 real, dimension(size(drho_dT)) :: pres ! Pressure converted to [Pa]
799 real, dimension(size(drho_dT)) :: Ta ! Temperature converted to [degC]
800 real, dimension(size(drho_dT)) :: Sa ! Salinity converted to [ppt]
801 real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
802 real :: dRdT_scale ! A factor to convert drho_dT to the desired units [R degC m3 C-1 kg-1 ~> 1]
803 real :: dRdS_scale ! A factor to convert drho_dS to the desired units [R ppt m3 S-1 kg-1 ~> 1]
804 integer :: i, is, ie, npts
805
806 if (present(dom)) then
807 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
808 else
809 is = 1 ; ie = size(drho_dt) ; npts = 1 + ie - is
810 endif
811
812 if ((eos%RL2_T2_to_Pa == 1.0) .and. (eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
813 call calculate_density_derivs_array(t, s, pressure, drho_dt, drho_ds, is, npts, eos)
814 else
815 do i=is,ie
816 pres(i) = eos%RL2_T2_to_Pa * pressure(i)
817 ta(i) = eos%C_to_degC * t(i)
818 sa(i) = eos%S_to_ppt * s(i)
819 enddo
820 call calculate_density_derivs_array(ta, sa, pres, drho_dt, drho_ds, is, npts, eos)
821 endif
822
823 rho_scale = eos%kg_m3_to_R
824 if (present(scale)) rho_scale = rho_scale * scale
825 drdt_scale = rho_scale * eos%C_to_degC
826 drds_scale = rho_scale * eos%S_to_ppt
827 if ((drdt_scale /= 1.0) .or. (drds_scale /= 1.0)) then ; do i=is,ie
828 drho_dt(i) = drdt_scale * drho_dt(i)
829 drho_ds(i) = drds_scale * drho_ds(i)
830 enddo ; endif
831
832end subroutine calculate_density_derivs_1d
833
834
835!> Calls the appropriate subroutines to calculate density derivatives by promoting a scalar
836!! to a one-element array
837subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS, scale)
838 real, intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
839 real, intent(in) :: S !< Salinity [S ~> ppt]
840 real, intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
841 real, intent(out) :: drho_dT !< The partial derivative of density with potential
842 !! temperature [R C-1 ~> kg m-3 degC-1] or other
843 !! units determined by the optional scale argument
844 real, intent(out) :: drho_dS !< The partial derivative of density with salinity,
845 !! in [R S-1 ~> kg m-3 ppt-1] or other units
846 !! determined by the optional scale argument
847 type(eos_type), intent(in) :: EOS !< Equation of state structure
848 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
849 !! in combination with scaling stored in EOS [various]
850 ! Local variables
851 real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
852 real :: dRdT_scale ! A factor to convert drho_dT to the desired units [R degC m3 C-1 kg-1 ~> 1]
853 real :: dRdS_scale ! A factor to convert drho_dS to the desired units [R ppt m3 S-1 kg-1 ~> 1]
854 real :: pres(1) ! Pressure converted to [Pa]
855 real :: Ta(1) ! Temperature converted to [degC]
856 real :: Sa(1) ! Salinity converted to [ppt]
857
858 pres(1) = eos%RL2_T2_to_Pa*pressure
859 ta(1) = eos%C_to_degC * t
860 sa(1) = eos%S_to_ppt * s
861
862 call eos%type%calculate_density_derivs_scalar(ta(1), sa(1), pres(1), drho_dt, drho_ds)
863
864 rho_scale = eos%kg_m3_to_R
865 if (present(scale)) rho_scale = rho_scale * scale
866 drdt_scale = rho_scale * eos%C_to_degC
867 drds_scale = rho_scale * eos%S_to_ppt
868 if ((drdt_scale /= 1.0) .or. (drds_scale /= 1.0)) then
869 drho_dt = drdt_scale * drho_dt
870 drho_ds = drds_scale * drho_ds
871 endif
872
874
875!> Calls the appropriate subroutine to calculate density second derivatives for 1-D array inputs.
876subroutine calculate_density_second_derivs_1d(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, &
877 drho_dS_dP, drho_dT_dP, EOS, dom, scale)
878 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
879 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
880 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
881 real, dimension(:), intent(inout) :: drho_dS_dS !< Partial derivative of beta with respect to S
882 !! [R S-2 ~> kg m-3 ppt-2]
883 real, dimension(:), intent(inout) :: drho_dS_dT !< Partial derivative of beta with respect to T
884 !! [R S-1 C-1 ~> kg m-3 ppt-1 degC-1]
885 real, dimension(:), intent(inout) :: drho_dT_dT !< Partial derivative of alpha with respect to T
886 !! [R C-2 ~> kg m-3 degC-2]
887 real, dimension(:), intent(inout) :: drho_dS_dP !< Partial derivative of beta with respect to pressure
888 !! [T2 S-1 L-2 ~> kg m-3 ppt-1 Pa-1]
889 real, dimension(:), intent(inout) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure
890 !! [T2 C-1 L-2 ~> kg m-3 degC-1 Pa-1]
891 type(eos_type), intent(in) :: EOS !< Equation of state structure
892 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
893 !! into account that arrays start at 1.
894 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
895 !! in combination with scaling stored in EOS [various]
896 ! Local variables
897 real, dimension(size(T)) :: pres ! Pressure converted to [Pa]
898 real, dimension(size(T)) :: Ta ! Temperature converted to [degC]
899 real, dimension(size(T)) :: Sa ! Salinity converted to [ppt]
900 real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
901 integer :: i, is, ie, npts
902
903 if (.not. allocated(eos%type)) call mom_error(fatal, &
904 "calculate_density_second_derivs: EOS%form_of_EOS is not valid.")
905
906 if (present(dom)) then
907 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
908 else
909 is = 1 ; ie = size(t) ; npts = 1 + ie - is
910 endif
911
912 if ((eos%RL2_T2_to_Pa == 1.0) .and. (eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
913 call eos%type%calculate_density_second_derivs_array(t, s, pressure, &
914 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp, is, npts)
915 else
916 do i=is,ie
917 pres(i) = eos%RL2_T2_to_Pa * pressure(i)
918 ta(i) = eos%C_to_degC * t(i)
919 sa(i) = eos%S_to_ppt * s(i)
920 enddo
921 call eos%type%calculate_density_second_derivs_array(ta, sa, pres, drho_ds_ds, drho_ds_dt, &
922 drho_dt_dt, drho_ds_dp, drho_dt_dp, is, npts)
923 endif
924
925 rho_scale = eos%kg_m3_to_R
926 if (present(scale)) rho_scale = rho_scale * scale
927 if (rho_scale /= 1.0) then ; do i=is,ie
928 drho_ds_ds(i) = rho_scale * drho_ds_ds(i)
929 drho_ds_dt(i) = rho_scale * drho_ds_dt(i)
930 drho_dt_dt(i) = rho_scale * drho_dt_dt(i)
931 drho_ds_dp(i) = rho_scale * drho_ds_dp(i)
932 drho_dt_dp(i) = rho_scale * drho_dt_dp(i)
933 enddo ; endif
934
935 if (eos%RL2_T2_to_Pa /= 1.0) then ; do i=is,ie
936 drho_ds_dp(i) = eos%RL2_T2_to_Pa * drho_ds_dp(i)
937 drho_dt_dp(i) = eos%RL2_T2_to_Pa * drho_dt_dp(i)
938 enddo ; endif
939
940 if (eos%C_to_degC /= 1.0) then ; do i=is,ie
941 drho_ds_dt(i) = eos%C_to_degC * drho_ds_dt(i)
942 drho_dt_dt(i) = eos%C_to_degC**2 * drho_dt_dt(i)
943 drho_dt_dp(i) = eos%C_to_degC * drho_dt_dp(i)
944 enddo ; endif
945
946 if (eos%S_to_ppt /= 1.0) then ; do i=is,ie
947 drho_ds_ds(i) = eos%S_to_ppt**2 * drho_ds_ds(i)
948 drho_ds_dt(i) = eos%S_to_ppt * drho_ds_dt(i)
949 drho_ds_dp(i) = eos%S_to_ppt * drho_ds_dp(i)
950 enddo ; endif
951
953
954!> Calls the appropriate subroutine to calculate density second derivatives for scalar inputs.
955subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, &
956 drho_dS_dP, drho_dT_dP, EOS, scale)
957 real, intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
958 real, intent(in) :: S !< Salinity [S ~> ppt]
959 real, intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
960 real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S
961 !! [R S-2 ~> kg m-3 ppt-2]
962 real, intent(out) :: drho_dS_dT !< Partial derivative of beta with respect to T
963 !! [R S-1 C-1 ~> kg m-3 ppt-1 degC-1]
964 real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T
965 !! [R C-2 ~> kg m-3 degC-2]
966 real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure
967 !! [T2 S-1 L-2 ~> kg m-3 ppt-1 Pa-1]
968 real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure
969 !! [T2 C-1 L-2 ~> kg m-3 degC-1 Pa-1]
970 type(eos_type), intent(in) :: EOS !< Equation of state structure
971 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
972 !! in combination with scaling stored in EOS [various]
973 ! Local variables
974 real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
975 real :: pres ! Pressure converted to [Pa]
976 real :: Ta ! Temperature converted to [degC]
977 real :: Sa ! Salinity converted to [ppt]
978
979 if (.not. allocated(eos%type)) call mom_error(fatal, &
980 "calculate_density_second_derivs: EOS%form_of_EOS is not valid.")
981
982 pres = eos%RL2_T2_to_Pa*pressure
983 ta = eos%C_to_degC * t
984 sa = eos%S_to_ppt * s
985
986 call eos%type%calculate_density_second_derivs_scalar(ta, sa, pres, drho_ds_ds, drho_ds_dt, &
987 drho_dt_dt, drho_ds_dp, drho_dt_dp)
988
989 rho_scale = eos%kg_m3_to_R
990 if (present(scale)) rho_scale = rho_scale * scale
991 if (rho_scale /= 1.0) then
992 drho_ds_ds = rho_scale * drho_ds_ds
993 drho_ds_dt = rho_scale * drho_ds_dt
994 drho_dt_dt = rho_scale * drho_dt_dt
995 drho_ds_dp = rho_scale * drho_ds_dp
996 drho_dt_dp = rho_scale * drho_dt_dp
997 endif
998
999 if (eos%RL2_T2_to_Pa /= 1.0) then
1000 drho_ds_dp = eos%RL2_T2_to_Pa * drho_ds_dp
1001 drho_dt_dp = eos%RL2_T2_to_Pa * drho_dt_dp
1002 endif
1003
1004 if (eos%C_to_degC /= 1.0) then
1005 drho_ds_dt = eos%C_to_degC * drho_ds_dt
1006 drho_dt_dt = eos%C_to_degC**2 * drho_dt_dt
1007 drho_dt_dp = eos%C_to_degC * drho_dt_dp
1008 endif
1009
1010 if (eos%S_to_ppt /= 1.0) then
1011 drho_ds_ds = eos%S_to_ppt**2 * drho_ds_ds
1012 drho_ds_dt = eos%S_to_ppt * drho_ds_dt
1013 drho_ds_dp = eos%S_to_ppt * drho_ds_dp
1014 endif
1015
1017
1018!> Calls the appropriate subroutine to calculate specific volume derivatives for an array.
1019subroutine calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS)
1020 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC]
1021 real, dimension(:), intent(in) :: S !< Salinity [ppt]
1022 real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
1023 real, dimension(:), intent(inout) :: dSV_dT !< The partial derivative of specific volume with potential
1024 !! temperature [m3 kg-1 degC-1]
1025 real, dimension(:), intent(inout) :: dSV_dS !< The partial derivative of specific volume with salinity
1026 !! [m3 kg-1 ppt-1]
1027 integer, intent(in) :: start !< Starting index within the array
1028 integer, intent(in) :: npts !< The number of values to calculate
1029 type(eos_type), intent(in) :: EOS !< Equation of state structure
1030
1031 if (.not. allocated(eos%type)) call mom_error(fatal, &
1032 "calculate_spec_vol_derivs_array: EOS%form_of_EOS is not valid.")
1033
1034 call eos%type%calculate_specvol_derivs_array(t, s, pressure, dsv_dt, dsv_ds, start, npts)
1035
1037
1038!> Calls the appropriate subroutine to calculate specific volume derivatives for 1-d array inputs,
1039!! potentially limiting the domain of indices that are worked on.
1040subroutine calc_spec_vol_derivs_1d(T, S, pressure, dSV_dT, dSV_dS, EOS, dom, scale)
1041 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
1042 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
1043 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
1044 real, dimension(:), intent(inout) :: dSV_dT !< The partial derivative of specific volume with potential
1045 !! temperature [R-1 C-1 ~> m3 kg-1 degC-1]
1046 real, dimension(:), intent(inout) :: dSV_dS !< The partial derivative of specific volume with salinity
1047 !! [R-1 S-1 ~> m3 kg-1 ppt-1]
1048 type(eos_type), intent(in) :: EOS !< Equation of state structure
1049 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
1050 !! into account that arrays start at 1.
1051 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific
1052 !! volume in combination with scaling stored in EOS [various]
1053
1054 ! Local variables
1055 real, dimension(size(T)) :: pres ! Pressure converted to [Pa]
1056 real, dimension(size(T)) :: Ta ! Temperature converted to [degC]
1057 real, dimension(size(T)) :: Sa ! Salinity converted to [ppt]
1058 real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg R-1 m-3 ~> 1]
1059 real :: dSVdT_scale ! A factor to convert dSV_dT to the desired units [kg degC R-1 C-1 m-3 ~> 1]
1060 real :: dSVdS_scale ! A factor to convert dSV_dS to the desired units [kg ppt R-1 S-1 m-3 ~> 1]
1061 integer :: i, is, ie, npts
1062
1063 if (present(dom)) then
1064 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
1065 else
1066 is = 1 ; ie = size(dsv_dt) ; npts = 1 + ie - is
1067 endif
1068
1069 if ((eos%RL2_T2_to_Pa == 1.0) .and. (eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
1070 call calculate_spec_vol_derivs_array(t, s, pressure, dsv_dt, dsv_ds, is, npts, eos)
1071 else
1072 do i=is,ie
1073 pres(i) = eos%RL2_T2_to_Pa * pressure(i)
1074 ta(i) = eos%C_to_degC * t(i)
1075 sa(i) = eos%S_to_ppt * s(i)
1076 enddo
1077 call calculate_spec_vol_derivs_array(ta, sa, pres, dsv_dt, dsv_ds, is, npts, eos)
1078 endif
1079
1080 spv_scale = eos%R_to_kg_m3
1081 if (present(scale)) spv_scale = spv_scale * scale
1082 dsvdt_scale = spv_scale * eos%C_to_degC
1083 dsvds_scale = spv_scale * eos%S_to_ppt
1084 if ((dsvdt_scale /= 1.0) .or. (dsvds_scale /= 1.0)) then ; do i=is,ie
1085 dsv_dt(i) = dsvdt_scale * dsv_dt(i)
1086 dsv_ds(i) = dsvds_scale * dsv_ds(i)
1087 enddo ; endif
1088
1089end subroutine calc_spec_vol_derivs_1d
1090
1091
1092!> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array
1093!! inputs. The inputs and outputs use dimensionally rescaled units.
1094subroutine calculate_compress_1d(T, S, pressure, rho, drho_dp, EOS, dom)
1095 real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
1096 real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
1097 real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
1098 real, dimension(:), intent(inout) :: rho !< In situ density [R ~> kg m-3]
1099 real, dimension(:), intent(inout) :: drho_dp !< The partial derivative of density with pressure
1100 !! (also the inverse of the square of sound speed)
1101 !! [T2 L-2 ~> s2 m-2]
1102 type(eos_type), intent(in) :: EOS !< Equation of state structure
1103 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
1104 !! into account that arrays start at 1.
1105
1106 ! Local variables
1107 real, dimension(size(T)) :: pres ! Pressure converted to [Pa]
1108 real, dimension(size(T)) :: Ta ! Temperature converted to [degC]
1109 real, dimension(size(T)) :: Sa ! Salinity converted to [ppt]
1110 integer :: i, is, ie, npts
1111
1112 if (.not. allocated(eos%type)) call mom_error(fatal, &
1113 "calculate_compress_1d: EOS%form_of_EOS is not valid.")
1114
1115 if (present(dom)) then
1116 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
1117 else
1118 is = 1 ; ie = size(rho) ; npts = 1 + ie - is
1119 endif
1120
1121 do i=is,ie
1122 pres(i) = eos%RL2_T2_to_Pa * pressure(i)
1123 ta(i) = eos%C_to_degC * t(i)
1124 sa(i) = eos%S_to_ppt * s(i)
1125 enddo
1126
1127 call eos%type%calculate_compress_array(ta, sa, pres, rho, drho_dp, is, npts)
1128
1129 if (eos%kg_m3_to_R /= 1.0) then ; do i=is,ie
1130 rho(i) = eos%kg_m3_to_R * rho(i)
1131 enddo ; endif
1132 if (eos%L_T_to_m_s /= 1.0) then ; do i=is,ie
1133 drho_dp(i) = eos%L_T_to_m_s**2 * drho_dp(i)
1134 enddo ; endif
1135
1136end subroutine calculate_compress_1d
1137
1138!> Calculate density and compressibility for a scalar. This just promotes the scalar to an array
1139!! with a singleton dimension and calls calculate_compress_1d. The inputs and outputs use
1140!! dimensionally rescaled units.
1141subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS)
1142 real, intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
1143 real, intent(in) :: S !< Salinity [S ~> ppt]
1144 real, intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
1145 real, intent(out) :: rho !< In situ density [R ~> kg m-3]
1146 real, intent(out) :: drho_dp !< The partial derivative of density with pressure (also the
1147 !! inverse of the square of sound speed) [T2 L-2 ~> s2 m-2]
1148 type(eos_type), intent(in) :: EOS !< Equation of state structure
1149
1150 ! Local variables
1151 ! These arrays use the same units as their counterparts in calculate_compress_1d.
1152 real, dimension(1) :: pa ! Pressure in a size-1 1d array [R L2 T-2 ~> Pa]
1153 real, dimension(1) :: Ta ! Temperature in a size-1 1d array [C ~> degC]
1154 real, dimension(1) :: Sa ! Salinity in a size-1 1d array [S ~> ppt]
1155 real, dimension(1) :: rhoa ! In situ density in a size-1 1d array [R ~> kg m-3]
1156 real, dimension(1) :: drho_dpa ! The partial derivative of density with pressure (also the
1157 ! inverse of the square of sound speed) in a 1d array [T2 L-2 ~> s2 m-2]
1158
1159 ta(1) = t ; sa(1) = s ; pa(1) = pressure
1160
1161 call calculate_compress_1d(ta, sa, pa, rhoa, drho_dpa, eos)
1162 rho = rhoa(1) ; drho_dp = drho_dpa(1)
1163
1164end subroutine calculate_compress_scalar
1165
1166!> Calls the appropriate subroutine to calculate the layer averaged specific volume either using
1167!! Boole's rule quadrature or analytical and nearly-analytical averages in pressure.
1168subroutine average_specific_vol(T, S, p_t, dp, SpV_avg, EOS, dom, scale)
1169 real, dimension(:), intent(in) :: t !< Potential temperature referenced to the surface [C ~> degC]
1170 real, dimension(:), intent(in) :: s !< Salinity [S ~> ppt]
1171 real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
1172 real, dimension(:), intent(in) :: dp !< Pressure change in the layer [R L2 T-2 ~> Pa]
1173 real, dimension(:), intent(inout) :: spv_avg !< The vertical average specific volume
1174 !! in the layer [R-1 ~> m3 kg-1]
1175 type(eos_type), intent(in) :: eos !< Equation of state structure
1176 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
1177 !! into account that arrays start at 1.
1178 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale
1179 !! output specific volume in combination with
1180 !! scaling stored in EOS [various]
1181
1182 ! Local variables
1183 real, dimension(size(T)) :: pres ! Layer-top pressure converted to [Pa]
1184 real, dimension(size(T)) :: dpres ! Pressure change converted to [Pa]
1185 real, dimension(size(T)) :: ta ! Temperature converted to [degC]
1186 real, dimension(size(T)) :: sa ! Salinity converted to [ppt]
1187 real :: t5(5) ! Temperatures at five quadrature points [C ~> degC]
1188 real :: s5(5) ! Salinities at five quadrature points [S ~> ppt]
1189 real :: p5(5) ! Pressures at five quadrature points [R L2 T-2 ~> Pa]
1190 real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1]
1191 real, parameter :: c1_90 = 1.0/90.0 ! A rational constant [nondim]
1192 real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1]
1193 integer :: i, n, is, ie, npts
1194
1195 if (present(dom)) then
1196 is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is
1197 else
1198 is = 1 ; ie = size(t) ; npts = 1 + ie - is
1199 endif
1200
1201 if (eos%EOS_quadrature) then
1202 do i=is,ie
1203 do n=1,5
1204 t5(n) = t(i) ; s5(n) = s(i)
1205 p5(n) = p_t(i) + 0.25*real(5-n)*dp(i)
1206 enddo
1207 call calculate_spec_vol(t5, s5, p5, a5, eos)
1208
1209 ! Use Boole's rule to estimate the average specific volume.
1210 spv_avg(i) = c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3))
1211 enddo
1212 elseif ((eos%RL2_T2_to_Pa == 1.0) .and. (eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
1213 select case (eos%form_of_EOS)
1214 case (eos_linear)
1215 call avg_spec_vol_linear(t, s, p_t, dp, spv_avg, is, npts, eos%Rho_T0_S0, &
1216 eos%dRho_dT, eos%dRho_dS, eos%dRho_dp)
1217 case (eos_wright)
1218 call avg_spec_vol_buggy_wright(t, s, p_t, dp, spv_avg, is, npts)
1219 case (eos_wright_full)
1220 call avg_spec_vol_wright_full(t, s, p_t, dp, spv_avg, is, npts)
1221 case (eos_wright_reduced)
1222 call avg_spec_vol_wright_red(t, s, p_t, dp, spv_avg, is, npts)
1223 case default
1224 call mom_error(fatal, "No analytic average specific volume option is available with this EOS!")
1225 end select
1226 else
1227 do i=is,ie
1228 pres(i) = eos%RL2_T2_to_Pa * p_t(i)
1229 dpres(i) = eos%RL2_T2_to_Pa * dp(i)
1230 ta(i) = eos%C_to_degC * t(i)
1231 sa(i) = eos%S_to_ppt * s(i)
1232 enddo
1233 select case (eos%form_of_EOS)
1234 case (eos_linear)
1235 call avg_spec_vol_linear(ta, sa, pres, dpres, spv_avg, is, npts, eos%Rho_T0_S0, &
1236 eos%dRho_dT, eos%dRho_dS, eos%dRho_dp)
1237 case (eos_wright)
1238 call avg_spec_vol_buggy_wright(ta, sa, pres, dpres, spv_avg, is, npts)
1239 case (eos_wright_full)
1240 call avg_spec_vol_wright_full(ta, sa, pres, dpres, spv_avg, is, npts)
1241 case (eos_wright_reduced)
1242 call avg_spec_vol_wright_red(ta, sa, pres, dpres, spv_avg, is, npts)
1243 case default
1244 call mom_error(fatal, "No analytic average specific volume option is available with this EOS!")
1245 end select
1246 endif
1247
1248 spv_scale = eos%R_to_kg_m3
1249 if (eos%EOS_quadrature) spv_scale = 1.0
1250 if (present(scale)) spv_scale = spv_scale * scale
1251 if (spv_scale /= 1.0) then ; do i=is,ie
1252 spv_avg(i) = spv_scale * spv_avg(i)
1253 enddo ; endif
1254
1255end subroutine average_specific_vol
1256
1257!> Return the range of temperatures, salinities and pressures for which the equation of state that
1258!! is being used has been fitted to observations. Care should be taken when applying
1259!! this equation of state outside of its fit range.
1260subroutine eos_fit_range(EOS, T_min, T_max, S_min, S_max, p_min, p_max)
1261 type(eos_type), intent(in) :: eos !< Equation of state structure
1262 real, optional, intent(out) :: t_min !< The minimum temperature over which this EoS is fitted [degC]
1263 real, optional, intent(out) :: t_max !< The maximum temperature over which this EoS is fitted [degC]
1264 real, optional, intent(out) :: s_min !< The minimum salinity over which this EoS is fitted [ppt]
1265 real, optional, intent(out) :: s_max !< The maximum salinity over which this EoS is fitted [ppt]
1266 real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
1267 real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
1268
1269 if (.not. allocated(eos%type)) call mom_error(fatal, &
1270 "calculate_compress: EOS%form_of_EOS is not valid.")
1271
1272 call eos%type%EoS_fit_range(t_min, t_max, s_min, s_max, p_min, p_max)
1273
1274end subroutine eos_fit_range
1275
1276
1277!> This subroutine returns a two point integer array indicating the domain of i-indices
1278!! to work on in EOS calls based on information from a hor_index type
1279function eos_domain(HI, halo) result(EOSdom)
1280 type(hor_index_type), intent(in) :: hi !< The horizontal index structure
1281 integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0.
1282 integer, dimension(2) :: eosdom !< The index domain that the EOS will work on, taking into account
1283 !! that the arrays inside the EOS routines will start at 1.
1284
1285 ! Local variables
1286 integer :: halo_sz
1287
1288 halo_sz = 0 ; if (present(halo)) halo_sz = halo
1289
1290 eosdom(1) = hi%isc - (hi%isd-1) - halo_sz
1291 eosdom(2) = hi%iec - (hi%isd-1) + halo_sz
1292
1293end function eos_domain
1294
1295!> Calls the appropriate subroutine to calculate analytical and nearly-analytical
1296!! integrals in pressure across layers of geopotential anomalies, which are
1297!! required for calculating the finite-volume form pressure accelerations in a
1298!! non-Boussinesq model. There are essentially no free assumptions, apart from the
1299!! use of Boole's rule to do the horizontal integrals, and from a truncation in the
1300!! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
1301subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, &
1302 dza, intp_dza, intx_dza, inty_dza, halo_size, &
1303 bathyP, P_surf, dP_tiny, MassWghtInterp)
1304 type(hor_index_type), intent(in) :: hi !< The horizontal index structure
1305 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1306 intent(in) :: t !< Potential temperature referenced to the surface [C ~> degC]
1307 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1308 intent(in) :: s !< Salinity [S ~> ppt]
1309 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1310 intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa]
1311 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1312 intent(in) :: p_b !< Pressure at the bottom of the layer [R L2 T-2 ~> Pa]
1313 real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
1314 !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]
1315 !! The calculation is mathematically identical with different values of
1316 !! alpha_ref, but this reduces the effects of roundoff.
1317 type(eos_type), intent(in) :: eos !< Equation of state structure
1318 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1319 intent(inout) :: dza !< The change in the geopotential anomaly across
1320 !! the layer [L2 T-2 ~> m2 s-2]
1321 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1322 optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the
1323 !! geopotential anomaly relative to the anomaly at the bottom of the
1324 !! layer [R L4 T-4 ~> Pa m2 s-2]
1325 real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
1326 optional, intent(inout) :: intx_dza !< The integral in x of the difference between the
1327 !! geopotential anomaly at the top and bottom of the layer divided by
1328 !! the x grid spacing [L2 T-2 ~> m2 s-2]
1329 real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
1330 optional, intent(inout) :: inty_dza !< The integral in y of the difference between the
1331 !! geopotential anomaly at the top and bottom of the layer divided by
1332 !! the y grid spacing [L2 T-2 ~> m2 s-2]
1333 integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
1334 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1335 optional, intent(in) :: bathyp !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
1336 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1337 optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
1338 real, optional, intent(in) :: dp_tiny !< A miniscule pressure change with
1339 !! the same units as p_t [R L2 T-2 ~> Pa]
1340 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use
1341 !! mass weighting to interpolate T/S in integrals
1342
1343 ! Local variables
1344 real :: drdt_scale ! A factor to convert drho_dT to the desired units [R degC m3 C-1 kg-1 ~> 1]
1345 real :: drds_scale ! A factor to convert drho_dS to the desired units [R ppt m3 S-1 kg-1 ~> 1]
1346 real :: drdp_scale ! A factor to convert drho_dp to the desired units [T-2 L2 s2 m-2 ~> 1]
1347
1348 ! We should never reach this point with quadrature. EOS_quadrature indicates that numerical
1349 ! integration be used instead of analytic. This is a safety check.
1350 if (eos%EOS_quadrature) call mom_error(fatal, "EOS_quadrature is set!")
1351
1352 select case (eos%form_of_EOS)
1353 case (eos_linear)
1354 drdt_scale = eos%kg_m3_to_R * eos%C_to_degC
1355 drds_scale = eos%kg_m3_to_R * eos%S_to_ppt
1356 drdp_scale = eos%kg_m3_to_R * eos%RL2_T2_to_Pa
1357 call int_spec_vol_dp_linear(t, s, p_t, p_b, alpha_ref, hi, eos%kg_m3_to_R*eos%Rho_T0_S0, &
1358 drdt_scale*eos%dRho_dT, drds_scale*eos%dRho_dS, drdp_scale*eos%dRho_dp, &
1359 dza, intp_dza, intx_dza, inty_dza, halo_size, &
1360 bathyp, p_surf, dp_tiny, masswghtinterp)
1361 case (eos_wright)
1362 call int_spec_vol_dp_wright(t, s, p_t, p_b, alpha_ref, hi, dza, intp_dza, intx_dza, &
1363 inty_dza, halo_size, bathyp, p_surf, dp_tiny, masswghtinterp, &
1364 sv_scale=eos%R_to_kg_m3, pres_scale=eos%RL2_T2_to_Pa, &
1365 temp_scale=eos%C_to_degC, saln_scale=eos%S_to_ppt)
1366 case (eos_wright_full)
1367 call int_spec_vol_dp_wright_full(t, s, p_t, p_b, alpha_ref, hi, dza, intp_dza, intx_dza, &
1368 inty_dza, halo_size, bathyp, p_surf, dp_tiny, masswghtinterp, &
1369 sv_scale=eos%R_to_kg_m3, pres_scale=eos%RL2_T2_to_Pa, &
1370 temp_scale=eos%C_to_degC, saln_scale=eos%S_to_ppt)
1371 case (eos_wright_reduced)
1372 call int_spec_vol_dp_wright_red(t, s, p_t, p_b, alpha_ref, hi, dza, intp_dza, intx_dza, &
1373 inty_dza, halo_size, bathyp, p_surf, dp_tiny, masswghtinterp, &
1374 sv_scale=eos%R_to_kg_m3, pres_scale=eos%RL2_T2_to_Pa, &
1375 temp_scale=eos%C_to_degC, saln_scale=eos%S_to_ppt)
1376 case default
1377 call mom_error(fatal, "No analytic integration option is available with this EOS!")
1378 end select
1379
1380end subroutine analytic_int_specific_vol_dp
1381
1382!> This subroutine calculates analytical and nearly-analytical integrals of
1383!! pressure anomalies across layers, which are required for calculating the
1384!! finite-volume form pressure accelerations in a Boussinesq model.
1385subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, &
1386 intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p)
1387 type(hor_index_type), intent(in) :: hi !< Ocean horizontal index structure
1388 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1389 intent(in) :: t !< Potential temperature referenced to the surface [C ~> degC]
1390 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1391 intent(in) :: s !< Salinity [S ~> ppt]
1392 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1393 intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m]
1394 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1395 intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m]
1396 real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is
1397 !! subtracted out to reduce the magnitude of each of the
1398 !! integrals.
1399 real, intent(in) :: rho_0 !< A density [R ~> kg m-3], that is used
1400 !! to calculate the pressure (as p~=-z*rho_0*G_e)
1401 !! used in the equation of state.
1402 real, intent(in) :: g_e !< The Earth's gravitational acceleration
1403 !! [L2 Z-1 T-2 ~> m s-2]
1404 type(eos_type), intent(in) :: eos !< Equation of state structure
1405 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1406 intent(inout) :: dpa !< The change in the pressure anomaly
1407 !! across the layer [R L2 T-2 ~> Pa]
1408 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1409 optional, intent(inout) :: intz_dpa !< The integral through the thickness of the
1410 !! layer of the pressure anomaly relative to the
1411 !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]
1412 real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
1413 optional, intent(inout) :: intx_dpa !< The integral in x of the difference between
1414 !! the pressure anomaly at the top and bottom of the
1415 !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]
1416 real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
1417 optional, intent(inout) :: inty_dpa !< The integral in y of the difference between
1418 !! the pressure anomaly at the top and bottom of the
1419 !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
1420 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1421 optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m]
1422 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1423 optional, intent(in) :: ssh !< The sea surface height [Z ~> m]
1424 real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]
1425 integer, optional, intent(in) :: masswghtinterp !< A flag indicating whether and how to use
1426 !! mass weighting to interpolate T/S in integrals
1427 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
1428 optional, intent(in) :: z_0p !< The height at which the pressure is 0 [Z ~> m]
1429
1430 ! Local variables
1431 real :: rho_scale ! A multiplicative factor by which to scale density from kg m-3 to the
1432 ! desired units [R m3 kg-1 ~> 1]
1433 real :: drdt_scale ! A factor to convert drho_dT to the desired units [R degC m3 C-1 kg-1 ~> 1]
1434 real :: drds_scale ! A factor to convert drho_dS to the desired units [R ppt m3 S-1 kg-1 ~> 1]
1435 real :: drdp_scale ! A factor to convert drho_dp to the desired units [T-2 L2 s2 m-2 ~> 1]
1436 real :: pres_scale ! A multiplicative factor to convert pressure into Pa [Pa T2 R-1 L-2 ~> 1]
1437
1438 ! We should never reach this point with quadrature. EOS_quadrature indicates that numerical
1439 ! integration be used instead of analytic. This is a safety check.
1440 if (eos%EOS_quadrature) call mom_error(fatal, "EOS_quadrature is set!")
1441
1442 select case (eos%form_of_EOS)
1443 case (eos_linear)
1444 rho_scale = eos%kg_m3_to_R
1445 drdt_scale = eos%kg_m3_to_R * eos%C_to_degC
1446 drds_scale = eos%kg_m3_to_R * eos%S_to_ppt
1447 drdp_scale = eos%kg_m3_to_R * eos%RL2_T2_to_Pa
1448 if ((rho_scale /= 1.0) .or. (drdt_scale /= 1.0) .or. (drds_scale /= 1.0) .or. (drdp_scale /= 1.0)) then
1449 call int_density_dz_linear(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, rho_scale*eos%Rho_T0_S0, &
1450 drdt_scale*eos%dRho_dT, drds_scale*eos%dRho_dS, drdp_scale*eos%dRho_dp, &
1451 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, dz_neglect, masswghtinterp, z_0p=z_0p)
1452 else
1453 call int_density_dz_linear(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, &
1454 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, eos%dRho_dp, &
1455 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, dz_neglect, masswghtinterp, z_0p=z_0p)
1456 endif
1457 case (eos_wright)
1458 rho_scale = eos%kg_m3_to_R
1459 pres_scale = eos%RL2_T2_to_Pa
1460 if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (eos%C_to_degC /= 1.0) .or. (eos%S_to_ppt /= 1.0)) then
1461 call int_density_dz_wright(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, &
1462 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, &
1463 dz_neglect, masswghtinterp, rho_scale, pres_scale, &
1464 temp_scale=eos%C_to_degC, saln_scale=eos%S_to_ppt, z_0p=z_0p)
1465 else
1466 call int_density_dz_wright(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, &
1467 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, &
1468 dz_neglect, masswghtinterp, z_0p=z_0p)
1469 endif
1470 case (eos_wright_full)
1471 rho_scale = eos%kg_m3_to_R
1472 pres_scale = eos%RL2_T2_to_Pa
1473 if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (eos%C_to_degC /= 1.0) .or. (eos%S_to_ppt /= 1.0)) then
1474 call int_density_dz_wright_full(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, &
1475 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, &
1476 dz_neglect, masswghtinterp, rho_scale, pres_scale, &
1477 temp_scale=eos%C_to_degC, saln_scale=eos%S_to_ppt, z_0p=z_0p)
1478 else
1479 call int_density_dz_wright_full(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, &
1480 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, &
1481 dz_neglect, masswghtinterp, z_0p=z_0p)
1482 endif
1483 case (eos_wright_reduced)
1484 rho_scale = eos%kg_m3_to_R
1485 pres_scale = eos%RL2_T2_to_Pa
1486 if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (eos%C_to_degC /= 1.0) .or. (eos%S_to_ppt /= 1.0)) then
1487 call int_density_dz_wright_red(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, &
1488 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, &
1489 dz_neglect, masswghtinterp, rho_scale, pres_scale, &
1490 temp_scale=eos%C_to_degC, saln_scale=eos%S_to_ppt, z_0p=z_0p)
1491 else
1492 call int_density_dz_wright_red(t, s, z_t, z_b, rho_ref, rho_0, g_e, hi, &
1493 dpa, intz_dpa, intx_dpa, inty_dpa, bathyt, ssh, &
1494 dz_neglect, masswghtinterp, z_0p=z_0p)
1495 endif
1496 case default
1497 call mom_error(fatal, "No analytic integration option is available with this EOS!")
1498 end select
1499
1500end subroutine analytic_int_density_dz
1501
1502!> Returns true if the equation of state is compressible (i.e. has pressure dependence)
1503logical function query_compressible(EOS)
1504 type(eos_type), intent(in) :: eos !< Equation of state structure
1505
1506 query_compressible = eos%compressible
1507end function query_compressible
1508
1509!> Returns the string identifying the equation of state with enumeration "id"
1510function get_eos_name(id) result (eos_name)
1511 integer, optional, intent(in) :: id !< Enumerated ID
1512 character(:), allocatable :: eos_name !< The name of the EOS
1513
1514 select case (id)
1515 case (eos_linear)
1516 eos_name = eos_linear_string
1517 case (eos_unesco)
1518 eos_name = eos_unesco_string
1519 case (eos_wright)
1520 eos_name = eos_wright_string
1521 case (eos_wright_reduced)
1522 eos_name = eos_wright_red_string
1523 case (eos_wright_full)
1524 eos_name = eos_wright_full_string
1525 case (eos_teos10)
1526 eos_name = eos_teos10_string
1527 case (eos_roquet_rho)
1528 eos_name = eos_roquet_rho_string
1529 case (eos_roquet_spv)
1530 eos_name = eos_roquet_spv_string
1531 case (eos_jackett06)
1532 eos_name = eos_jackett06_string
1533 case default
1534 call mom_error(fatal, "get_EOS_name: something went wrong internally - enumeration is not valid.")
1535 end select
1536
1537end function get_eos_name
1538
1539!> Initializes EOS_type by allocating and reading parameters. The scaling factors in
1540!! US are stored in EOS for later use.
1541subroutine eos_init(param_file, EOS, US, use_conT_absS)
1542 type(param_file_type), intent(in) :: param_file !< Parameter file structure
1543 type(eos_type), intent(inout) :: eos !< Equation of state structure
1544 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1545 logical, intent(in), optional :: use_cont_abss !< True if the model is formulated for
1546 !! conservative temp and absolute salinity
1547 optional :: us
1548 ! Local variables
1549# include "version_variable.h"
1550 character(len=40) :: mdl = "MOM_EOS" ! This module's name.
1551 character(len=12) :: tfreeze_default ! The default freezing point expression
1552 character(len=40) :: tmpstr
1553 logical :: eos_quad_default, eos_ts_default
1554 real :: rho_tref_sref ! Density at Tref degC and Sref ppt [kg m-3]
1555 real :: tref ! Reference temperature [degC]
1556 real :: sref ! Reference salinity [psu]
1557 real :: pref ! Reference pressure [Pa]
1558 real :: rho0 ! Density at T=0, S=0 and p=0 [kg m-3]
1559
1560 ! Read all relevant parameters and write them to the model log.
1561 call log_version(param_file, mdl, version, "")
1562
1563 call get_param(param_file, mdl, "EQN_OF_STATE", tmpstr, &
1564 "EQN_OF_STATE determines which ocean equation of state should be used. "//&
1565 'Currently, the valid choices are "LINEAR", "UNESCO", "JACKETT_MCD", '//&
1566 '"WRIGHT", "WRIGHT_REDUCED", "WRIGHT_FULL", "NEMO", "ROQUET_RHO", "ROQUET_SPV" '//&
1567 'and "TEOS10". This is only used if USE_EOS is true.', default=eos_default)
1568 select case (uppercase(tmpstr))
1569 case (eos_linear_string)
1570 call eos_manual_init(eos, form_of_eos=eos_linear)
1571 case (eos_unesco_string)
1572 call eos_manual_init(eos, form_of_eos=eos_unesco)
1573 case (eos_jackett_string)
1574 call eos_manual_init(eos, form_of_eos=eos_unesco)
1575 case (eos_wright_string)
1576 call eos_manual_init(eos, form_of_eos=eos_wright)
1577 case (eos_wright_red_string)
1578 call eos_manual_init(eos, form_of_eos=eos_wright_reduced)
1579 case (eos_wright_full_string)
1580 call eos_manual_init(eos, form_of_eos=eos_wright_full)
1581 case (eos_teos10_string)
1582 call eos_manual_init(eos, form_of_eos=eos_teos10)
1583 case (eos_nemo_string)
1584 call eos_manual_init(eos, form_of_eos=eos_roquet_rho)
1585 case (eos_roquet_rho_string)
1586 call eos_manual_init(eos, form_of_eos=eos_roquet_rho)
1587 case (eos_roquet_spv_string)
1588 call eos_manual_init(eos, form_of_eos=eos_roquet_spv)
1589 case (eos_jackett06_string)
1590 call eos_manual_init(eos, form_of_eos=eos_jackett06)
1591 case default
1592 call mom_error(fatal, "interpret_eos_selection: EQN_OF_STATE "//&
1593 trim(tmpstr) // " in input file is invalid.")
1594 end select
1595 call mom_mesg('interpret_eos_selection: equation of state set to "' // &
1596 trim(tmpstr)//'"', 5)
1597
1598 if (eos%form_of_EOS == eos_linear) then
1599 ! RHO(T,S) = RHO_REF + DRHO_DT*(T-T_REF) + DRHO_DS*(S-S_REF) + DRHO_DP*(P-P_REF)
1600 ! = RHO_REF - (DRHO_DT*T_REF + DRHO_DS*SREF + DRHO_DP*PREF) + (DRHO_DT*T + DRHO_DS*S + DRHO_DP*P)
1601 ! = RHO_T0_S0 + (DRHO_DT*T + DRHO_DS*S + DRHO_DP*P)
1602 call get_param(param_file, mdl, "RHO_REF_LINEAR_EOS", rho_tref_sref, &
1603 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is the density "//&
1604 "at T=T_REF_LINEAR_EOS, S=S_REF_LINEAR_EOS and p=P_REF_LINEAR_EOS", &
1605 units="kg m-3", default=1000.0, old_name="RHO_TREF_SREF")
1606 call get_param(param_file, mdl, "T_REF_LINEAR_EOS", tref, &
1607 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is the reference "//&
1608 "temperature.", units="degC", default=0.0, old_name="TREF")
1609 call get_param(param_file, mdl, "S_REF_LINEAR_EOS", sref, &
1610 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is the reference "//&
1611 "salinity.", units="psu", default=0.0, old_name="SREF")
1612 call get_param(param_file, mdl, "P_REF_LINEAR_EOS", pref, &
1613 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is the reference "//&
1614 "pressure.", units="Pa", default=0.0)
1615 call get_param(param_file, mdl, "DRHO_DT", eos%dRho_dT, &
1616 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is "//&
1617 "the partial derivative of density with temperature.", &
1618 units="kg m-3 K-1", default=-0.2)
1619 call get_param(param_file, mdl, "DRHO_DS", eos%dRho_dS, &
1620 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is "//&
1621 "the partial derivative of density with salinity.", &
1622 units="kg m-3 ppt-1", default=0.8)
1623 call get_param(param_file, mdl, "DRHO_DP", eos%dRho_dp, &
1624 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is "//&
1625 "the partial derivative of density with pressure (the inverse of "//&
1626 "sound speed squared).", units="s2 m-2", default=0.0)
1627 rho0 = rho_tref_sref - ((eos%dRho_dT * tref + eos%dRho_dS * sref) + eos%dRho_dp * pref)
1628 call get_param(param_file, mdl, "RHO_T0_S0", eos%Rho_T0_S0, &
1629 "When EQN_OF_STATE="//trim(eos_linear_string)//", this is the density "//&
1630 "at T=0, S=0 and p=0. If RHO_TO_SO is specified, RHO_REF_LINEAR_EOS, "//&
1631 "T_REF_LINEAR_EOS, S_REF_LINEAR_EOS and P_REF_LINEAR_EOS are not used.", &
1632 units="kg m-3", default=rho0)
1633 eos%Compressible = (eos%dRho_dp/=0.0)
1634 call eos_manual_init(eos, form_of_eos=eos_linear, rho_t0_s0=eos%Rho_T0_S0, &
1635 drho_dt=eos%dRho_dT, drho_ds=eos%dRho_dS, drho_dp=eos%dRho_dp)
1636 endif
1637 if (eos%form_of_EOS == eos_wright) then
1638 call get_param(param_file, mdl, "USE_WRIGHT_2ND_DERIV_BUG", eos%use_Wright_2nd_deriv_bug, &
1639 "If true, use a bug in the calculation of the second derivatives of density "//&
1640 "with temperature and with temperature and pressure that causes some terms "//&
1641 "to be only 2/3 of what they should be.", default=.false.)
1642 call eos_manual_init(eos, form_of_eos=eos_wright, use_wright_2nd_deriv_bug=eos%use_Wright_2nd_deriv_bug)
1643 endif
1644
1645 if (present(use_cont_abss)) then
1646 eos%use_conT_absS = use_cont_abss
1647 else
1648 eos%use_conT_absS = .false. ! Assuming it is not needed, it is set to false
1649 endif
1650
1651 eos_quad_default = .not.((eos%form_of_EOS == eos_linear) .or. &
1652 (eos%form_of_EOS == eos_wright) .or. &
1653 (eos%form_of_EOS == eos_wright_reduced) .or. &
1654 (eos%form_of_EOS == eos_wright_full))
1655 call get_param(param_file, mdl, "EOS_QUADRATURE", eos%EOS_quadrature, &
1656 "If true, always use the generic (quadrature) code "//&
1657 "code for the integrals of density.", default=eos_quad_default)
1658
1659 tfreeze_default = tfreeze_linear_string
1660 if ((eos%form_of_EOS == eos_teos10 .or. eos%form_of_EOS == eos_roquet_rho .or. &
1661 eos%form_of_EOS == eos_roquet_spv)) &
1662 tfreeze_default = tfreeze_teos10_string
1663 call get_param(param_file, mdl, "TFREEZE_FORM", tmpstr, &
1664 "TFREEZE_FORM determines which expression should be "//&
1665 "used for the freezing point. Currently, the valid "//&
1666 'choices are "LINEAR", "MILLERO_78", "TEOS_POLY", "TEOS10"', &
1667 default=tfreeze_default)
1668 select case (uppercase(tmpstr))
1669 case (tfreeze_linear_string)
1670 eos%form_of_TFreeze = tfreeze_linear
1671 case (tfreeze_millero_string)
1672 eos%form_of_TFreeze = tfreeze_millero
1673 case (tfreeze_teospoly_string)
1674 eos%form_of_TFreeze = tfreeze_teospoly
1675 case (tfreeze_teos10_string)
1676 eos%form_of_TFreeze = tfreeze_teos10
1677 case default
1678 call mom_error(fatal, "interpret_eos_selection: TFREEZE_FORM "//&
1679 trim(tmpstr) // "in input file is invalid.")
1680 end select
1681
1682 if (eos%form_of_TFreeze == tfreeze_linear) then
1683 call get_param(param_file, mdl, "TFREEZE_S0_P0",eos%TFr_S0_P0, &
1684 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
1685 "this is the freezing potential temperature at "//&
1686 "S=0, P=0.", units="degC", default=0.0)
1687 call get_param(param_file, mdl, "DTFREEZE_DS",eos%dTFr_dS, &
1688 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
1689 "this is the derivative of the freezing potential "//&
1690 "temperature with salinity.", &
1691 units="degC ppt-1", default=-0.054)
1692 call get_param(param_file, mdl, "DTFREEZE_DP",eos%dTFr_dP, &
1693 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
1694 "this is the derivative of the freezing potential "//&
1695 "temperature with pressure.", &
1696 units="degC Pa-1", default=0.0)
1697 endif
1698
1699 if ((eos%form_of_TFreeze==tfreeze_teospoly) .or. (eos%form_of_TFreeze==tfreeze_teos10)) then
1700 ! Which default is appropriate for Millero?
1701 eos_ts_default = .false.
1702 else
1703 eos_ts_default = .true.
1704 endif
1705 call get_param(param_file, mdl, "TFREEZE_S_IS_PRACS", eos%TFreeze_S_is_pracS, &
1706 "When True, the model will check if the model internal salinity is "//&
1707 "practical salinity. If the model uses absolute salinity, a "//&
1708 "conversion will be applied.", default=eos_ts_default)
1709 call get_param(param_file, mdl, "TFREEZE_T_IS_POTT", eos%TFreeze_T_is_potT, &
1710 "When True, the model will check if the model internal temperature is "//&
1711 "potential temperature. If the model uses conservative temperature, a "//&
1712 "conversion will be applied.", default=eos_ts_default)
1713
1714 if ((eos%form_of_EOS == eos_teos10 .or. eos%form_of_EOS == eos_roquet_rho .or. &
1715 eos%form_of_EOS == eos_roquet_spv) .and. &
1716 .not.((eos%form_of_TFreeze == tfreeze_teos10) .or. (eos%form_of_TFreeze == tfreeze_teospoly)) ) then
1717 call mom_error(warning, "interpret_eos_selection: EOS_TEOS10 or EOS_ROQUET_RHO or EOS_ROQUET_SPV "//&
1718 "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 or TFREEZE_TEOSPOLY.")
1719 endif
1720
1721 ! Unit conversions
1722 eos%m_to_Z = 1. ; if (present(us)) eos%m_to_Z = us%m_to_Z
1723 eos%kg_m3_to_R = 1. ; if (present(us)) eos%kg_m3_to_R = us%kg_m3_to_R
1724 eos%R_to_kg_m3 = 1. ; if (present(us)) eos%R_to_kg_m3 = us%R_to_kg_m3
1725 eos%RL2_T2_to_Pa = 1. ; if (present(us)) eos%RL2_T2_to_Pa = us%RL2_T2_to_Pa
1726 eos%L_T_to_m_s = 1. ; if (present(us)) eos%L_T_to_m_s = us%L_T_to_m_s
1727 eos%degC_to_C = 1. ; if (present(us)) eos%degC_to_C = us%degC_to_C
1728 eos%C_to_degC = 1. ; if (present(us)) eos%C_to_degC = us%C_to_degC
1729 eos%ppt_to_S = 1. ; if (present(us)) eos%ppt_to_S = us%ppt_to_S
1730 eos%S_to_ppt = 1. ; if (present(us)) eos%S_to_ppt = us%S_to_ppt
1731
1732end subroutine eos_init
1733
1734!> Manually initialized an EOS type (intended for unit testing of routines which need a specific EOS)
1735subroutine eos_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, &
1736 Rho_T0_S0, drho_dT, dRho_dS, dRho_dp, TFr_S0_P0, dTFr_dS, dTFr_dp, &
1737 use_Wright_2nd_deriv_bug)
1738 type(eos_type), intent(inout) :: eos !< Equation of state structure
1739 integer, optional, intent(in) :: form_of_eos !< A coded integer indicating the equation of state to use.
1740 integer, optional, intent(in) :: form_of_tfreeze !< A coded integer indicating the expression for
1741 !! the potential temperature of the freezing point.
1742 logical, optional, intent(in) :: eos_quadrature !< If true, always use the generic (quadrature)
1743 !! code for the integrals of density.
1744 logical, optional, intent(in) :: compressible !< If true, in situ density is a function of pressure.
1745 real , optional, intent(in) :: rho_t0_s0 !< Density at T=0 degC and S=0 ppt [kg m-3]
1746 real , optional, intent(in) :: drho_dt !< Partial derivative of density with temperature
1747 !! in [kg m-3 degC-1]
1748 real , optional, intent(in) :: drho_ds !< Partial derivative of density with salinity
1749 !! in [kg m-3 ppt-1]
1750 real , optional, intent(in) :: drho_dp !< Partial derivative of density with pressure
1751 !! in [s2 m-2]
1752 real , optional, intent(in) :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC]
1753 real , optional, intent(in) :: dtfr_ds !< The derivative of freezing point with salinity
1754 !! in [degC ppt-1]
1755 real , optional, intent(in) :: dtfr_dp !< The derivative of freezing point with pressure
1756 !! in [degC Pa-1]
1757 logical, optional, intent(in) :: use_wright_2nd_deriv_bug !< Allow the Wright 2nd deriv bug
1758
1759 if (present(form_of_eos)) then
1760 eos%form_of_EOS = form_of_eos
1761 if (allocated(eos%type)) deallocate(eos%type) ! Needed during testing which re-initializes
1762 select case (eos%form_of_EOS)
1763 case (eos_linear)
1764 allocate(linear_eos :: eos%type)
1765 case (eos_unesco)
1766 allocate(unesco_eos :: eos%type)
1767 case (eos_wright)
1768 allocate(buggy_wright_eos :: eos%type)
1769 case (eos_wright_full)
1770 allocate(wright_full_eos :: eos%type)
1771 case (eos_wright_reduced)
1772 allocate(wright_red_eos :: eos%type)
1773 case (eos_jackett06)
1774 allocate(jackett06_eos :: eos%type)
1775 case (eos_teos10)
1776 allocate(teos10_eos :: eos%type)
1777 case (eos_roquet_rho)
1778 allocate(roquet_rho_eos :: eos%type)
1779 case (eos_roquet_spv)
1780 allocate(roquet_spv_eos :: eos%type)
1781 end select
1782 select type (t => eos%type)
1783 type is (linear_eos)
1784 call t%set_params_linear(rho_t0_s0, drho_dt, drho_ds, drho_dp)
1785 type is (buggy_wright_eos)
1786 call t%set_params_buggy_Wright(use_wright_2nd_deriv_bug)
1787 end select
1788 endif
1789 if (present(form_of_tfreeze)) eos%form_of_TFreeze = form_of_tfreeze
1790 if (present(eos_quadrature )) eos%EOS_quadrature = eos_quadrature
1791 if (present(compressible )) eos%Compressible = compressible
1792 if (present(rho_t0_s0 )) eos%Rho_T0_S0 = rho_t0_s0
1793 if (present(drho_dt )) eos%drho_dT = drho_dt
1794 if (present(drho_ds )) eos%dRho_dS = drho_ds
1795 if (present(drho_dp )) eos%dRho_dp = drho_dp
1796 if (present(tfr_s0_p0 )) eos%TFr_S0_P0 = tfr_s0_p0
1797 if (present(dtfr_ds )) eos%dTFr_dS = dtfr_ds
1798 if (present(dtfr_dp )) eos%dTFr_dp = dtfr_dp
1799 if (present(use_wright_2nd_deriv_bug)) eos%use_Wright_2nd_deriv_bug = use_wright_2nd_deriv_bug
1800
1801end subroutine eos_manual_init
1802
1803! !> Set equation of state structure (EOS) to linear with given coefficients
1804! !!
1805! !! \note This routine is primarily for testing and allows a local copy of the
1806! !! EOS_type (EOS argument) to be set to use the linear equation of state
1807! !! independent from the rest of the model.
1808! subroutine EOS_use_linear(Rho_T0_S0, dRho_dT, dRho_dS, EOS, use_quadrature)
1809! real, intent(in) :: Rho_T0_S0 !< Density at T=0 degC and S=0 ppt [kg m-3]
1810! real, intent(in) :: dRho_dT !< Partial derivative of density with temperature [kg m-3 degC-1]
1811! real, intent(in) :: dRho_dS !< Partial derivative of density with salinity [kg m-3 ppt-1]
1812! logical, optional, intent(in) :: use_quadrature !< If true, always use the generic (quadrature)
1813! !! code for the integrals of density.
1814! type(EOS_type), intent(inout) :: EOS !< Equation of state structure
1815
1816! call EOS_manual_init(EOS, form_of_EOS=EOS_LINEAR, Rho_T0_S0=Rho_T0_S0, dRho_dT=dRho_dT, dRho_dS=dRho_dS)
1817! EOS%Compressible = .false.
1818! EOS%EOS_quadrature = .false.
1819! if (present(use_quadrature)) EOS%EOS_quadrature = use_quadrature
1820
1821! end subroutine EOS_use_linear
1822
1823!> Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10
1824subroutine convert_temp_salt_for_teos10(T, S, HI, kd, mask_z, EOS)
1825 integer, intent(in) :: kd !< The number of layers to work on
1826 type(hor_index_type), intent(in) :: hi !< The horizontal index structure
1827 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
1828 intent(inout) :: t !< Potential temperature referenced to the surface [C ~> degC]
1829 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
1830 intent(inout) :: s !< Salinity [S ~> ppt]
1831 real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
1832 intent(in) :: mask_z !< 3d mask regulating which points to convert [nondim]
1833 type(eos_type), intent(in) :: eos !< Equation of state structure
1834
1835 real, parameter :: sref_sprac = (35.16504/35.0) ! The TEOS 10 conversion factor to go from
1836 ! practical salinity to reference salinity [PSU ppt-1]
1837 integer :: i, j, k
1838
1839 if ((eos%form_of_EOS /= eos_teos10) .and. (eos%form_of_EOS /= eos_roquet_rho) .and. &
1840 (eos%form_of_EOS /= eos_roquet_spv)) return
1841
1842 do k=1,kd ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
1843 if (mask_z(i,j,k) >= 1.0) then
1844 s(i,j,k) = sref_sprac * s(i,j,k)
1845 t(i,j,k) = eos%degC_to_C*potemp_to_constemp(eos%C_to_degC*t(i,j,k), eos%S_to_ppt*s(i,j,k))
1846 endif
1847 enddo ; enddo ; enddo
1848end subroutine convert_temp_salt_for_teos10
1849
1850
1851!> Converts an array of conservative temperatures to potential temperatures. The input arguments
1852!! use the dimensionally rescaling as specified within the EOS type. The output potential
1853!! temperature uses this same scaling, but this can be replaced by the factor given by scale.
1854subroutine cons_temp_to_pot_temp(T, S, poTemp, EOS, dom, scale)
1855 real, dimension(:), intent(in) :: t !< Conservative temperature [C ~> degC]
1856 real, dimension(:), intent(in) :: s !< Absolute salinity [S ~> ppt]
1857 real, dimension(:), intent(inout) :: potemp !< The potential temperature with a reference pressure
1858 !! of 0 Pa, [C ~> degC]
1859 type(eos_type), intent(in) :: eos !< Equation of state structure
1860 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
1861 !! into account that arrays start at 1.
1862 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale the output
1863 !! potential temperature in place of with scaling stored
1864 !! in EOS. A value of 1.0 returns temperatures in [degC],
1865 !! while the default is equivalent to EOS%degC_to_C.
1866
1867 ! Local variables
1868 real, dimension(size(T)) :: ta ! Temperature converted to [degC]
1869 real, dimension(size(S)) :: sa ! Salinity converted to [ppt]
1870 real :: t_scale ! A factor to convert potential temperature from degC to the desired units [C degC-1 ~> 1]
1871 integer :: i, is, ie
1872
1873 if (present(dom)) then
1874 is = dom(1) ; ie = dom(2)
1875 else
1876 is = 1 ; ie = size(t)
1877 endif
1878
1879 if ((eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
1880 potemp(is:ie) = constemp_to_potemp(t(is:ie), s(is:ie))
1881 else
1882 do i=is,ie
1883 ta(i) = eos%C_to_degC * t(i)
1884 sa(i) = eos%S_to_ppt * s(i)
1885 enddo
1886 potemp(is:ie) = constemp_to_potemp(ta(is:ie), sa(is:ie))
1887 endif
1888
1889 t_scale = eos%degC_to_C
1890 if (present(scale)) t_scale = scale
1891 if (t_scale /= 1.0) then ; do i=is,ie
1892 potemp(i) = t_scale * potemp(i)
1893 enddo ; endif
1894
1895end subroutine cons_temp_to_pot_temp
1896
1897
1898!> Converts an array of potential temperatures to conservative temperatures. The input arguments
1899!! use the dimensionally rescaling as specified within the EOS type. The output potential
1900!! temperature uses this same scaling, but this can be replaced by the factor given by scale.
1901subroutine pot_temp_to_cons_temp(T, S, consTemp, EOS, dom, scale)
1902 real, dimension(:), intent(in) :: t !< Potential temperature [C ~> degC]
1903 real, dimension(:), intent(in) :: s !< Absolute salinity [S ~> ppt]
1904 real, dimension(:), intent(inout) :: constemp !< The conservative temperature [C ~> degC]
1905 type(eos_type), intent(in) :: eos !< Equation of state structure
1906 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
1907 !! into account that arrays start at 1.
1908 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale the output
1909 !! potential temperature in place of with scaling stored
1910 !! in EOS. A value of 1.0 returns temperatures in [degC],
1911 !! while the default is equivalent to EOS%degC_to_C.
1912
1913 ! Local variables
1914 real, dimension(size(T)) :: tp ! Potential temperature converted to [degC]
1915 real, dimension(size(S)) :: sa ! Absolute salinity converted to [ppt]
1916 real :: t_scale ! A factor to convert potential temperature from degC to the desired units [C degC-1 ~> 1]
1917 integer :: i, is, ie
1918
1919 if (present(dom)) then
1920 is = dom(1) ; ie = dom(2)
1921 else
1922 is = 1 ; ie = size(t)
1923 endif
1924
1925
1926 if ((eos%C_to_degC == 1.0) .and. (eos%S_to_ppt == 1.0)) then
1927 constemp(is:ie) = potemp_to_constemp(t(is:ie), s(is:ie))
1928 else
1929 do i=is,ie
1930 tp(i) = eos%C_to_degC * t(i)
1931 sa(i) = eos%S_to_ppt * s(i)
1932 enddo
1933 constemp(is:ie) = potemp_to_constemp(tp(is:ie), sa(is:ie))
1934 endif
1935
1936 t_scale = eos%degC_to_C
1937 if (present(scale)) t_scale = scale
1938 if (t_scale /= 1.0) then ; do i=is,ie
1939 constemp(i) = t_scale * constemp(i)
1940 enddo ; endif
1941
1942end subroutine pot_temp_to_cons_temp
1943
1944
1945!> Converts an array of absolute salinity to practical salinity. The input arguments
1946!! use the dimensionally rescaling as specified within the EOS type. The output potential
1947!! temperature uses this same scaling, but this can be replaced by the factor given by scale.
1948subroutine abs_saln_to_prac_saln(S, prSaln, EOS, dom, scale)
1949 real, dimension(:), intent(in) :: s !< Absolute salinity [S ~> ppt]
1950 real, dimension(:), intent(inout) :: prsaln !< Practical salinity [S ~> PSU]
1951 type(eos_type), intent(in) :: eos !< Equation of state structure
1952 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
1953 !! into account that arrays start at 1.
1954 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale the output
1955 !! practical salinities in place of with scaling stored
1956 !! in EOS. A value of 1.0 returns salinities in [PSU],
1957 !! while the default is equivalent to EOS%ppt_to_S.
1958
1959 ! Local variables
1960 real :: s_scale ! A factor to convert practical salinity from ppt to the desired units [S PSU-1 ~> 1]
1961 real, parameter :: sprac_sref = (35.0/35.16504) ! The TEOS 10 conversion factor to go from
1962 ! reference salinity to practical salinity [PSU ppt-1]
1963 integer :: i, is, ie
1964
1965 if (present(dom)) then
1966 is = dom(1) ; ie = dom(2)
1967 else
1968 is = 1 ; ie = size(s)
1969 endif
1970
1971 if (present(scale)) then
1972 s_scale = sprac_sref * scale
1973 do i=is,ie
1974 prsaln(i) = s_scale * s(i)
1975 enddo
1976 else
1977 do i=is,ie
1978 prsaln(i) = sprac_sref * s(i)
1979 enddo
1980 endif
1981
1982end subroutine abs_saln_to_prac_saln
1983
1984
1985!> Converts an array of absolute salinity to practical salinity. The input arguments
1986!! use the dimensionally rescaling as specified within the EOS type. The output potential
1987!! temperature uses this same scaling, but this can be replaced by the factor given by scale.
1988subroutine prac_saln_to_abs_saln(S, absSaln, EOS, dom, scale)
1989 real, dimension(:), intent(in) :: s !< Practical salinity [S ~> PSU]
1990 real, dimension(:), intent(inout) :: abssaln !< Absolute salinity [S ~> ppt]
1991 type(eos_type), intent(in) :: eos !< Equation of state structure
1992 integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
1993 !! into account that arrays start at 1.
1994 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale the output
1995 !! absolute salnities in place of with scaling stored
1996 !! in EOS. A value of 1.0 returns salinities in [ppt],
1997 !! while the default is equivalent to EOS%ppt_to_S.
1998
1999 ! Local variables
2000 real :: s_scale ! A factor to convert absolute salinity from ppt to the desired units [S ppt-1 ~> 1]
2001 real, parameter :: sref_sprac = (35.16504/35.0) ! The TEOS 10 conversion factor to go from
2002 ! practical salinity to reference salinity [PSU ppt-1]
2003 integer :: i, is, ie
2004
2005 if (present(dom)) then
2006 is = dom(1) ; ie = dom(2)
2007 else
2008 is = 1 ; ie = size(s)
2009 endif
2010
2011 if (present(scale)) then
2012 s_scale = sref_sprac * scale
2013 do i=is,ie
2014 abssaln(i) = s_scale * s(i)
2015 enddo
2016 else
2017 do i=is,ie
2018 abssaln(i) = sref_sprac * s(i)
2019 enddo
2020 endif
2021
2022end subroutine prac_saln_to_abs_saln
2023
2024
2025!> Return value of EOS_quadrature
2026logical function eos_quadrature(EOS)
2027 type(eos_type), intent(in) :: eos !< Equation of state structure
2028
2029 eos_quadrature = eos%EOS_quadrature
2030
2031end function eos_quadrature
2032
2033!> Runs unit tests for consistency on the equations of state.
2034!! This should only be called from a single/root thread.
2035!! It returns True if any test fails, otherwise it returns False.
2036logical function eos_unit_tests(verbose)
2037 logical, intent(in) :: verbose !< If true, write results to stdout
2038 ! Local variables
2039 type(eos_type) :: eos_tmp
2040 logical :: fail
2041
2042 if (verbose) write(stdout,*) '==== MOM_EOS: EOS_unit_tests ===='
2043 eos_unit_tests = .false. ! Normally return false
2044
2045 call eos_manual_init(eos_tmp, form_of_eos=eos_teos10)
2046 fail = test_ts_conversion_consistency(t_cons=9.989811727177308, s_abs=35.16504, &
2047 t_pot=10.0, s_prac=35.0, eos=eos_tmp, verbose=verbose)
2048 if (verbose .and. fail) call mom_error(warning, "Some EOS variable conversions tests have failed.")
2049 eos_unit_tests = eos_unit_tests .or. fail
2050
2051 call eos_manual_init(eos_tmp, form_of_eos=eos_unesco)
2052 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "UNESCO", &
2053 rho_check=1027.54345796120*eos_tmp%kg_m3_to_R)
2054 if (verbose .and. fail) call mom_error(warning, "UNESCO EOS has failed some self-consistency tests.")
2055 eos_unit_tests = eos_unit_tests .or. fail
2056
2057 call eos_manual_init(eos_tmp, form_of_eos=eos_wright_full)
2058 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "WRIGHT_FULL", &
2059 rho_check=1027.55177447616*eos_tmp%kg_m3_to_R, avg_sv_check=.true.)
2060 if (verbose .and. fail) call mom_error(warning, "WRIGHT_FULL EOS has failed some self-consistency tests.")
2061 eos_unit_tests = eos_unit_tests .or. fail
2062
2063 call eos_manual_init(eos_tmp, form_of_eos=eos_wright_reduced)
2064 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "WRIGHT_REDUCED", &
2065 rho_check=1027.54303596346*eos_tmp%kg_m3_to_R, avg_sv_check=.true.)
2066 if (verbose .and. fail) call mom_error(warning, "WRIGHT_REDUCED EOS has failed some self-consistency tests.")
2067 eos_unit_tests = eos_unit_tests .or. fail
2068
2069 ! This test is deliberately outside of the fit range for WRIGHT_REDUCED, and it results in the expected warnings.
2070 ! call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT_REDUCED)
2071 ! fail = test_EOS_consistency(25.0, 15.0, 1.0e7, EOS_tmp, verbose, "WRIGHT_REDUCED", &
2072 ! rho_check=1012.625699301455*EOS_tmp%kg_m3_to_R)
2073 ! if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT_REDUCED EOS has failed some self-consistency tests.")
2074 ! EOS_unit_tests = EOS_unit_tests .or. fail
2075
2076 call eos_manual_init(eos_tmp, form_of_eos=eos_wright, use_wright_2nd_deriv_bug=.true.)
2077 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "WRIGHT", &
2078 rho_check=1027.54303596346*eos_tmp%kg_m3_to_R, avg_sv_check=.true.)
2079 ! These last test is a known failure and since MPI is not necessarily initializaed when running these tests
2080 ! we need to avoid flagging the fails.
2081 !if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT EOS has failed some self-consistency tests.")
2082 !EOS_unit_tests = EOS_unit_tests .or. fail
2083
2084 call eos_manual_init(eos_tmp, form_of_eos=eos_roquet_rho)
2085 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "ROQUET_RHO", &
2086 rho_check=1027.42385663668*eos_tmp%kg_m3_to_R)
2087 if (verbose .and. fail) call mom_error(warning, "ROQUET_RHO EOS has failed some self-consistency tests.")
2088 eos_unit_tests = eos_unit_tests .or. fail
2089
2090 call eos_manual_init(eos_tmp, form_of_eos=eos_roquet_spv)
2091 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "ROQUET_SPV", &
2092 rho_check=1027.42387475199*eos_tmp%kg_m3_to_R)
2093 if (verbose .and. fail) call mom_error(warning, "ROQUET_SPV EOS has failed some self-consistency tests.")
2094 eos_unit_tests = eos_unit_tests .or. fail
2095
2096 call eos_manual_init(eos_tmp, form_of_eos=eos_jackett06)
2097 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "JACKETT06", &
2098 rho_check=1027.539690758425*eos_tmp%kg_m3_to_R)
2099 if (verbose .and. fail) call mom_error(warning, "JACKETT06 EOS has failed some self-consistency tests.")
2100 eos_unit_tests = eos_unit_tests .or. fail
2101
2102 ! The TEOS10 equation of state is not passing the self consistency tests for dho_dS_dp due
2103 ! to a bug (a missing division by the square root of offset-salinity) on line 111 of
2104 ! pkg/GSW-Fortan/toolbox/gsw_specvol_second_derivatives.f90. This bug has been highlighted in an
2105 ! issue posted to the TEOS-10/GSW-Fortran page at github.com/TEOS-10/GSW-Fortran/issues/26, and
2106 ! it will be corrected by github.com/mom-ocean/GSW-Fortran/pull/1 .
2107 call eos_manual_init(eos_tmp, form_of_eos=eos_teos10)
2108 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "TEOS10", skip_2nd=.true., &
2109 rho_check=1027.42355961492*eos_tmp%kg_m3_to_R)
2110 if (verbose .and. fail) call mom_error(warning, "TEOS10 EOS has failed some self-consistency tests.")
2111 eos_unit_tests = eos_unit_tests .or. fail
2112
2113 call eos_manual_init(eos_tmp, form_of_eos=eos_roquet_rho)
2114 fail = test_eos_consistency(10.0, 30.0, 1.0e7, eos_tmp, verbose, "ROQUET_RHO", &
2115 rho_check=1027.45140117152*eos_tmp%kg_m3_to_R)
2116 ! The corresponding check value published by Roquet et al. (2015) is 1027.45140 [kg m-3].
2117 if (verbose .and. fail) call mom_error(warning, "Roquet_rho EOS has failed some self-consistency tests.")
2118 eos_unit_tests = eos_unit_tests .or. fail
2119
2120 call eos_manual_init(eos_tmp, form_of_eos=eos_roquet_spv)
2121 fail = test_eos_consistency(10.0, 30.0, 1.0e7, eos_tmp, verbose, "ROQUET_SPV", &
2122 spv_check=9.73282046614623e-04*eos_tmp%R_to_kg_m3)
2123 ! The corresponding check value here published by Roquet et al. (2015) is 9.732819628e-04 [m3 kg-1],
2124 ! but the order of arithmetic there was not completely specified with parentheses.
2125 if (verbose .and. fail) call mom_error(warning, "ROQUET_SPV EOS has failed some self-consistency tests.")
2126 eos_unit_tests = eos_unit_tests .or. fail
2127
2128 call eos_manual_init(eos_tmp, form_of_eos=eos_linear, rho_t0_s0=1000.0, drho_dt=-0.2, drho_ds=0.8, drho_dp=5.0e-7)
2129 fail = test_eos_consistency(25.0, 35.0, 1.0e7, eos_tmp, verbose, "LINEAR", &
2130 rho_check=1028.0*eos_tmp%kg_m3_to_R, avg_sv_check=.true.)
2131 if (verbose .and. fail) call mom_error(warning, "LINEAR EOS has failed some self-consistency tests.")
2132 eos_unit_tests = eos_unit_tests .or. fail
2133
2134 ! Test the freezing point calculations
2135
2136 call eos_manual_init(eos_tmp, form_of_tfreeze=tfreeze_linear, tfr_s0_p0=0.0, dtfr_ds=-0.054, &
2137 dtfr_dp=-7.6e-8)
2138 fail = test_tfr_consistency(35.0, 1.0e7, eos_tmp, verbose, "LINEAR", tfr_check=-2.65*eos_tmp%degC_to_C)
2139 if (verbose .and. fail) call mom_error(warning, "LINEAR TFr has failed some self-consistency tests.")
2140 eos_unit_tests = eos_unit_tests .or. fail
2141
2142 call eos_manual_init(eos_tmp, form_of_tfreeze=tfreeze_millero)
2143 fail = test_tfr_consistency(35.0, 1.0e7, eos_tmp, verbose, "MILLERO_78", &
2144 tfr_check=-2.69730134114106*eos_tmp%degC_to_C)
2145 if (verbose .and. fail) call mom_error(warning, "MILLERO_78 TFr has failed some self-consistency tests.")
2146 eos_unit_tests = eos_unit_tests .or. fail
2147
2148 call eos_manual_init(eos_tmp, form_of_tfreeze=tfreeze_teos10)
2149 fail = test_tfr_consistency(35.0, 1.0e7, eos_tmp, verbose, "TEOS10", &
2150 tfr_check=-2.69099996992861*eos_tmp%degC_to_C)
2151 if (verbose .and. fail) call mom_error(warning, "TEOS10 TFr has failed some self-consistency tests.")
2152 eos_unit_tests = eos_unit_tests .or. fail
2153
2154 call eos_manual_init(eos_tmp, form_of_tfreeze=tfreeze_teospoly)
2155 fail = test_tfr_consistency(35.0, 1.0e7, eos_tmp, verbose, "TEOS_POLY", &
2156 tfr_check=-2.691165259327735*eos_tmp%degC_to_C)
2157 if (verbose .and. fail) call mom_error(warning, "TEOS_POLY TFr has failed some self-consistency tests.")
2158 eos_unit_tests = eos_unit_tests .or. fail
2159
2160 if (eos_unit_tests) then
2161 call mom_error(warning, "EOS_unit_tests: One or more EOS tests have failed!")
2162 else
2163 if (verbose) call mom_mesg("EOS_unit_tests: All EOS consistency tests have passed.")
2164 endif
2165
2166end function eos_unit_tests
2167
2168logical function test_ts_conversion_consistency(T_cons, S_abs, T_pot, S_prac, EOS, verbose) &
2169 result(inconsistent)
2170 real, intent(in) :: t_cons !< Conservative temperature [degC]
2171 real, intent(in) :: s_abs !< Absolute salinity [g kg-1]
2172 real, intent(in) :: t_pot !< Potential temperature [degC]
2173 real, intent(in) :: s_prac !< Practical salinity [PSU]
2174 type(eos_type), intent(in) :: eos !< Equation of state structure
2175 logical, intent(in) :: verbose !< If true, write results to stdout
2176
2177 ! Local variables
2178 real :: sabs(1) ! Absolute or reference salinity [g kg-1]
2179 real :: sprac(1) ! Practical salinity [PSU]
2180 real :: stest(1) ! A converted salinity [ppt]
2181 real :: tcons(1) ! Conservative temperature [degC]
2182 real :: tpot(1) ! Potential temperature [degC]
2183 real :: ttest(1) ! A converted temperature [degC]
2184 real :: stol ! Roundoff error on a typical value of salinities [ppt]
2185 real :: ttol ! Roundoff error on a typical value of temperatures [degC]
2186 logical :: test_ok ! True if a particular test is consistent.
2187 logical :: ok ! True if all checks so far are consistent.
2188
2189 ok = .true.
2190
2191 ! Copy scalar input values into the corresponding arrays
2192 sabs(1) = s_abs ; sprac(1) = s_prac ; tcons(1) = t_cons ; tpot(1) = t_pot
2193
2194 ! Set tolerances for the conversions.
2195 ttol = 2.0 * 400.0*epsilon(ttol)
2196 stol = 35.0 * 400.0*epsilon(stol)
2197
2198 ! Check that the converted salinities agree
2199 call abs_saln_to_prac_saln(sabs, stest, eos)
2200 test_ok = (abs(stest(1) - sprac(1)) <= stol)
2201 if (verbose) call write_check_msg("MOM6 Sprac", stest(1), sprac(1), stol, test_ok)
2202 ok = ok .and. test_ok
2203
2204 call prac_saln_to_abs_saln(sprac, stest, eos)
2205 test_ok = (abs(stest(1) - sabs(1)) <= stol)
2206 if (verbose) call write_check_msg("MOM6 Sabs", stest(1), sabs(1), stol, test_ok)
2207 ok = ok .and. test_ok
2208
2209 call cons_temp_to_pot_temp(tcons, sabs, ttest, eos)
2210 test_ok = (abs(ttest(1) - tpot(1)) <= ttol)
2211 if (verbose) call write_check_msg("MOM6 Tpot", ttest(1), tpot(1), ttol, test_ok)
2212 ok = ok .and. test_ok
2213
2214 call pot_temp_to_cons_temp(tpot, sabs, ttest, eos)
2215 test_ok = (abs(ttest(1) - tcons(1)) <= ttol)
2216 if (verbose) call write_check_msg("MOM6 Tcons", ttest(1), tcons(1), ttol, test_ok)
2217 ok = ok .and. test_ok
2218
2219 inconsistent = .not.ok
2221
2222logical function test_tfr_consistency(S_test, p_test, EOS, verbose, EOS_name, TFr_check) &
2223 result(inconsistent)
2224 real, intent(in) :: s_test !< Salinity or absolute salinity [S ~> ppt]
2225 real, intent(in) :: p_test !< Pressure [R L2 T-2 ~> Pa]
2226 type(eos_type), intent(in) :: eos !< Equation of state structure
2227 logical, intent(in) :: verbose !< If true, write results to stdout
2228 character(len=*), intent(in) :: eos_name !< A name used in error messages to describe the EoS
2229 real, optional, intent(in) :: tfr_check !< A check value for the Freezing point [C ~> degC]
2230
2231 ! Local variables
2232 real, dimension(-3:3,-3:3) :: s ! Salinities at the test value and perturbed points [S ~> ppt]
2233 real, dimension(-3:3,-3:3) :: p ! Pressures at the test value and perturbed points [R L2 T-2 ~> Pa]
2234 real, dimension(-3:3,-3:3,2) :: tfr ! Freezing point at the test value and perturbed points [C ~> degC]
2235 real :: ds ! Magnitude of salinity perturbations [S ~> ppt]
2236 real :: dp ! Magnitude of pressure perturbations [R L2 T-2 ~> Pa]
2237 ! real :: tol ! The nondimensional tolerance from roundoff [nondim]
2238 real :: tfr_tol ! Roundoff error on a typical value of TFreeze [C ~> degC]
2239 logical :: test_ok ! True if a particular test is consistent.
2240 logical :: ok ! True if all checks so far are consistent.
2241 integer :: i, j, n
2242
2243 ok = .true.
2244
2245 ds = 0.5*eos%ppt_to_S ! Salinity perturbations [S ~> ppt]
2246 dp = 10.0e4 / eos%RL2_T2_to_Pa ! Pressure perturbations [R L2 T-2 ~> Pa]
2247
2248 ! TEOS 10 requires a tolerance that is ~20 times larger than other freezing point
2249 ! expressions because it lacks parentheses.
2250 tfr_tol = 2.0*eos%degC_to_C * 400.0*epsilon(tfr_tol)
2251
2252 do n=1,2
2253 ! Calculate density values with a wide enough stencil to estimate first and second derivatives
2254 ! with up to 6th order accuracy. Doing this twice with different sizes of perturbations allows
2255 ! the evaluation of whether the finite differences are converging to the calculated values at a
2256 ! rate that is consistent with the order of accuracy of the finite difference forms, and hence
2257 ! the consistency of the calculated values.
2258 do j=-3,3 ; do i=-3,3
2259 s(i,j) = max(s_test + n*ds*i, 0.0)
2260 p(i,j) = max(p_test + n*dp*j, 0.0)
2261 enddo ; enddo
2262 do j=-3,3
2263 call calculate_tfreeze(s(:,j), p(:,j), tfr(:,j,n), eos)
2264 enddo
2265 enddo
2266
2267 ! Check that the freezing point agrees with the provided check value
2268 if (present(tfr_check)) then
2269 test_ok = (abs(tfr_check - tfr(0,0,1)) <= tfr_tol)
2270 ok = ok .and. test_ok
2271 if (verbose) call write_check_msg(trim(eos_name)//" TFr", tfr(0,0,1), tfr_check, tfr_tol, test_ok)
2272 endif
2273
2274 inconsistent = .not.ok
2275end function test_tfr_consistency
2276
2277!> Write a message indicating how well a value matches its check value.
2278subroutine write_check_msg(var_name, val, val_chk, val_tol, test_OK)
2279 character(len=*), intent(in) :: var_name !< The name of the variable being tested.
2280 real, intent(in) :: val !< The value being checked [various]
2281 real, intent(in) :: val_chk !< The value being checked [various]
2282 real, intent(in) :: val_tol !< The value being checked [various]
2283 logical, intent(in) :: test_OK !< True if the values are within their tolerance
2284
2285 character(len=200) :: mesg
2286
2287 write(mesg, '(ES24.16," vs. ",ES24.16,", diff=",ES12.4,", tol=",ES12.4)') &
2288 val, val_chk, val-val_chk, val_tol
2289 if (test_ok) then
2290 write(stdout,*) trim(var_name)//" agrees with its check value :"//trim(mesg)
2291 else
2292 write(stderr,*) trim(var_name)//" disagrees with its check value :"//trim(mesg)
2293 endif
2294end subroutine write_check_msg
2295
2296!> Test an equation of state for self-consistency and consistency with check values, returning false
2297!! if it is consistent by all tests, and true if it fails any test.
2298logical function test_eos_consistency(T_test, S_test, p_test, EOS, verbose, &
2299 EOS_name, rho_check, spv_check, skip_2nd, avg_Sv_check) result(inconsistent)
2300 real, intent(in) :: t_test !< Potential temperature or conservative temperature [C ~> degC]
2301 real, intent(in) :: s_test !< Salinity or absolute salinity [S ~> ppt]
2302 real, intent(in) :: p_test !< Pressure [R L2 T-2 ~> Pa]
2303 type(eos_type), intent(in) :: eos !< Equation of state structure
2304 logical, intent(in) :: verbose !< If true, write results to stdout
2305 character(len=*), intent(in) :: eos_name !< A name used in error messages to describe the EoS
2306 real, optional, intent(in) :: rho_check !< A check value for the density [R ~> kg m-3]
2307 real, optional, intent(in) :: spv_check !< A check value for the specific volume [R-1 ~> m3 kg-1]
2308 logical, optional, intent(in) :: skip_2nd !< If present and true, do not check the 2nd derivatives.
2309 logical, optional, intent(in) :: avg_sv_check !< If present and true, compare analytical and numerical
2310 !! quadrature estimates of the layer-averaged specific volume.
2311
2312 ! Local variables
2313 real, dimension(-3:3,-3:3,-3:3) :: t ! Temperatures at the test value and perturbed points [C ~> degC]
2314 real, dimension(-3:3,-3:3,-3:3) :: s ! Salinities at the test value and perturbed points [S ~> ppt]
2315 real, dimension(-3:3,-3:3,-3:3) :: p ! Pressures at the test value and perturbed points [R L2 T-2 ~> Pa]
2316 real, dimension(-3:3,-3:3,-3:3,2) :: rho ! Densities relative to rho_ref at the test value and
2317 ! perturbed points [R ~> kg m-3]
2318 real, dimension(-3:3,-3:3,-3:3,2) :: spv ! Specific volumes relative to spv_ref at the test value and
2319 ! perturbed points [R-1 ~> m3 kg-1]
2320 real :: dt ! Magnitude of temperature perturbations [C ~> degC]
2321 real :: ds ! Magnitude of salinity perturbations [S ~> ppt]
2322 real :: dp ! Magnitude of pressure perturbations [R L2 T-2 ~> Pa]
2323 real :: rho_ref ! A reference density that is extracted for greater accuracy [R ~> kg m-3]
2324 real :: spv_ref ! A reference specific volume that is extracted for greater accuracy [R-1 ~> m3 kg-1]
2325 real :: rho_nooff ! Density with no reference offset [R ~> kg m-3]
2326 real :: spv_nooff ! Specific volume with no reference offset [R-1 ~> m3 kg-1]
2327 real :: drho_dt ! The partial derivative of density with potential
2328 ! temperature [R C-1 ~> kg m-3 degC-1]
2329 real :: drho_ds ! The partial derivative of density with salinity
2330 ! in [R S-1 ~> kg m-3 ppt-1]
2331 real :: drho_dp ! The partial derivative of density with pressure (also the
2332 ! inverse of the square of sound speed) [T2 L-2 ~> s2 m-2]
2333 real :: dsv_dt(1) ! The partial derivative of specific volume with potential
2334 ! temperature [R-1 C-1 ~> m3 kg-1 degC-1]
2335 real :: dsv_ds(1) ! The partial derivative of specific volume with salinity
2336 ! [R-1 S-1 ~> m3 kg-1 ppt-1]
2337 real :: spv_avg_a(1) ! The pressure-averaged specific volume determined analytically [R-1 ~> m3 kg-1]
2338 real :: spv_avg_q(1) ! The pressure-averaged specific volume determined via quadrature [R-1 ~> m3 kg-1]
2339 real :: drho_ds_ds ! Second derivative of density with respect to S [R S-2 ~> kg m-3 ppt-2]
2340 real :: drho_ds_dt ! Second derivative of density with respect to T and S [R S-1 C-1 ~> kg m-3 ppt-1 degC-1]
2341 real :: drho_dt_dt ! Second derivative of density with respect to T [R C-2 ~> kg m-3 degC-2]
2342 real :: drho_ds_dp ! Second derivative of density with respect to salinity and pressure
2343 ! [T2 S-1 L-2 ~> kg m-3 ppt-1 Pa-1]
2344 real :: drho_dt_dp ! Second derivative of density with respect to temperature and pressure
2345 ! [T2 C-1 L-2 ~> kg m-3 degC-1 Pa-1]
2346
2347 real :: drho_dt_fd(2) ! Two 6th order finite difference estimates of the partial derivative of density
2348 ! with potential temperature [R C-1 ~> kg m-3 degC-1]
2349 real :: drho_ds_fd(2) ! Two 6th order finite difference estimates of the partial derivative of density
2350 ! with salinity [R S-1 ~> kg m-3 ppt-1]
2351 real :: drho_dp_fd(2) ! Two 6th order finite difference estimates of the partial derivative of density
2352 ! with pressure (also the inverse of the square of sound speed) [T2 L-2 ~> s2 m-2]
2353 real :: dsv_dt_fd(2) ! Two 6th order finite difference estimates of the partial derivative of
2354 ! specific volume with potential temperature [R-1 C-1 ~> m3 kg-1 degC-1]
2355 real :: dsv_ds_fd(2) ! Two 6th order finite difference estimates of the partial derivative of
2356 ! specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
2357 real :: drho_ds_ds_fd(2) ! Two 6th order finite difference estimates of the second derivative of
2358 ! density with respect to salinity [R S-2 ~> kg m-3 ppt-2]
2359 real :: drho_ds_dt_fd(2) ! Two 6th order finite difference estimates of the second derivative of density
2360 ! with respect to temperature and salinity [R S-1 C-1 ~> kg m-3 ppt-1 degC-1]
2361 real :: drho_dt_dt_fd(2) ! Two 6th order finite difference estimates of the second derivative of
2362 ! density with respect to temperature [R C-2 ~> kg m-3 degC-2]
2363 real :: drho_ds_dp_fd(2) ! Two 6th order finite difference estimates of the second derivative of density
2364 ! with respect to salinity and pressure [T2 S-1 L-2 ~> kg m-3 ppt-1 Pa-1]
2365 real :: drho_dt_dp_fd(2) ! Two 6th order finite difference estimates of the second derivative of density
2366 ! with respect to temperature and pressure [T2 C-1 L-2 ~> kg m-3 degC-1 Pa-1]
2367 real :: rho_tmp ! A temporary copy of the situ density [R ~> kg m-3]
2368 real :: tol ! The nondimensional tolerance from roundoff [nondim]
2369 real :: r_tol ! Roundoff error on a typical value of density anomaly [R ~> kg m-3]
2370 real :: sv_tol ! Roundoff error on a typical value of specific volume anomaly [R-1 ~> m3 kg-1]
2371 real :: tol_here ! The tolerance for each check, in various units [various]
2372 real :: t_min, t_max ! The minimum and maximum temperature over which this EoS is fitted [degC]
2373 real :: s_min, s_max ! The minimum and maximum temperature over which this EoS is fitted [ppt]
2374 real :: p_min, p_max ! The minimum and maximum temperature over which this EoS is fitted [Pa]
2375 real :: count_fac ! A factor in the roundoff estimates based on the factors in the numerator and
2376 ! denominator in the finite difference derivative expression [nondim]
2377 real :: count_fac2 ! A factor in the roundoff estimates based on the factors in the numerator and
2378 ! denominator in the finite difference second derivative expression [nondim]
2379 character(len=200) :: mesg
2380 type(eos_type) :: eos_tmp
2381 logical :: test_ok ! True if a particular test is consistent.
2382 logical :: ok ! True if all checks so far are consistent.
2383 logical :: test_2nd ! If true, do tests on the 2nd derivative calculations
2384 logical :: test_avg_sv ! If true, compare numerical and analytical estimates of the vertically
2385 ! averaged specific volume
2386 integer :: order ! The order of accuracy of the centered finite difference estimates (2, 4 or 6).
2387 integer :: i, j, k, n
2388
2389 test_2nd = .true. ; if (present(skip_2nd)) test_2nd = .not.skip_2nd
2390 test_avg_sv = .false. ; if (present(avg_sv_check)) test_avg_sv = avg_sv_check
2391
2392 dt = 0.1*eos%degC_to_C ! Temperature perturbations [C ~> degC]
2393 ds = 0.5*eos%ppt_to_S ! Salinity perturbations [S ~> ppt]
2394 dp = 10.0e4 / eos%RL2_T2_to_Pa ! Pressure perturbations [R L2 T-2 ~> Pa]
2395
2396 r_tol = 50.0*eos%kg_m3_to_R * 10.*epsilon(r_tol)
2397 sv_tol = 5.0e-5*eos%R_to_kg_m3 * 10.*epsilon(sv_tol)
2398 rho_ref = 1000.0*eos%kg_m3_to_R
2399 spv_ref = 1.0 / rho_ref
2400
2401 order = 4 ! This should be 2, 4 or 6.
2402
2403 ! Check whether the consistency test is being applied outside of the value range of this EoS.
2404 call eos_fit_range(eos, t_min, t_max, s_min, s_max, p_min, p_max)
2405 if ((t_test < t_min) .or. (t_test > t_max)) then
2406 write(mesg, '(ES12.4," [degC] which is outside of the fit range of ",ES12.4," to ",ES12.4)') t_test, t_min, t_max
2407 call mom_error(warning, trim(eos_name)//" is being evaluated at a temperature of "//trim(mesg))
2408 endif
2409 if ((s_test < s_min) .or. (s_test > s_max)) then
2410 write(mesg, '(ES12.4," [ppt] which is outside of the fit range of ",ES12.4," to ",ES12.4)') s_test, s_min, s_max
2411 call mom_error(warning, trim(eos_name)//" is being evaluated at a salinity of "//trim(mesg))
2412 endif
2413 if ((p_test < p_min) .or. (p_test > p_max)) then
2414 write(mesg, '(ES12.4," [Pa] which is outside of the fit range of ",ES12.4," to ",ES12.4)') p_test, p_min, p_max
2415 call mom_error(warning, trim(eos_name)//" is being evaluated at a pressure of "//trim(mesg))
2416 endif
2417
2418 do n=1,2
2419 ! Calculate density values with a wide enough stencil to estimate first and second derivatives
2420 ! with up to 6th order accuracy. Doing this twice with different sizes of perturbations allows
2421 ! the evaluation of whether the finite differences are converging to the calculated values at a
2422 ! rate that is consistent with the order of accuracy of the finite difference forms, and hence
2423 ! the consistency of the calculated values.
2424 do k=-3,3 ; do j=-3,3 ; do i=-3,3
2425 t(i,j,k) = t_test + n*dt*i
2426 s(i,j,k) = s_test + n*ds*j
2427 p(i,j,k) = p_test + n*dp*k
2428 enddo ; enddo ; enddo
2429 do k=-3,3 ; do j=-3,3
2430 call calculate_density(t(:,j,k), s(:,j,k), p(:,j,k), rho(:,j,k,n), eos, rho_ref=rho_ref)
2431 call calculate_spec_vol(t(:,j,k), s(:,j,k), p(:,j,k), spv(:,j,k,n), eos, spv_ref=spv_ref)
2432 enddo ; enddo
2433
2434 drho_dt_fd(n) = first_deriv(rho(:,0,0,n), n*dt, order)
2435 drho_ds_fd(n) = first_deriv(rho(0,:,0,n), n*ds, order)
2436 drho_dp_fd(n) = first_deriv(rho(0,0,:,n), n*dp, order)
2437 dsv_dt_fd(n) = first_deriv(spv(:,0,0,n), n*dt, order)
2438 dsv_ds_fd(n) = first_deriv(spv(0,:,0,n), n*ds, order)
2439 if (test_2nd) then
2440 drho_dt_dt_fd(n) = second_deriv(rho(:,0,0,n), n*dt, order)
2441 drho_ds_ds_fd(n) = second_deriv(rho(0,:,0,n), n*ds, order)
2442 drho_ds_dt_fd(n) = derivs_2d(rho(:,:,0,n), n**2*dt*ds, order)
2443 drho_dt_dp_fd(n) = derivs_2d(rho(:,0,:,n), n**2*dt*dp, order)
2444 drho_ds_dp_fd(n) = derivs_2d(rho(0,:,:,n), n**2*ds*dp, order)
2445 endif
2446 enddo
2447
2448 call calculate_density_derivs(t(0,0,0), s(0,0,0), p(0,0,0), drho_dt, drho_ds, eos)
2449 ! The first indices here are "0:0" because there is no scalar form of calculate_specific_vol_derivs.
2450 call calculate_specific_vol_derivs(t(0:0,0,0), s(0:0,0,0), p(0:0,0,0), dsv_dt, dsv_ds, eos)
2451 if (test_2nd) &
2452 call calculate_density_second_derivs(t(0,0,0), s(0,0,0), p(0,0,0), &
2453 drho_ds_ds, drho_ds_dt, drho_dt_dt, drho_ds_dp, drho_dt_dp, eos)
2454 call calculate_compress(t(0,0,0), s(0,0,0), p(0,0,0), rho_tmp, drho_dp, eos)
2455
2456 if (test_avg_sv) then
2457 eos_tmp = eos
2458 call eos_manual_init(eos_tmp, eos_quadrature=.false.)
2459 call average_specific_vol(t(0:0,0,0), s(0:0,0,0), p(0:0,0,0), p(0:0,0,0), spv_avg_a, eos_tmp)
2460 call eos_manual_init(eos_tmp, eos_quadrature=.true.)
2461 call average_specific_vol(t(0:0,0,0), s(0:0,0,0), p(0:0,0,0), p(0:0,0,0), spv_avg_q, eos_tmp)
2462 endif
2463
2464 ok = .true.
2465
2466 tol = 1000.0*epsilon(tol)
2467
2468 ! Check that the density agrees with the provided check value
2469 if (present(rho_check)) then
2470 test_ok = (abs(rho_check - (rho_ref + rho(0,0,0,1))) < tol*(rho_ref + rho(0,0,0,1)))
2471 ok = ok .and. test_ok
2472 if (verbose) &
2473 call write_check_msg(trim(eos_name)//" rho", rho_ref+rho(0,0,0,1), rho_check, tol*rho(0,0,0,1), test_ok)
2474 endif
2475
2476 ! Check that the specific volume agrees with the provided check value or the inverse of density
2477 if (present(spv_check)) then
2478 test_ok = (abs(spv_check - (spv_ref + spv(0,0,0,1))) < tol*abs(spv_ref + spv(0,0,0,1)))
2479 if (verbose) &
2480 call write_check_msg(trim(eos_name)//" spv", spv_ref+spv(0,0,0,1), spv_check, tol*spv(0,0,0,1), test_ok)
2481 ok = ok .and. test_ok
2482 else
2483 test_ok = (abs((rho_ref+rho(0,0,0,1)) * (spv_ref + spv(0,0,0,1)) - 1.0) < tol)
2484 ok = ok .and. test_ok
2485 if (verbose) then
2486 write(mesg, '(ES16.8," and ",ES16.8,", ratio - 1 = ",ES16.8)') &
2487 rho_ref+rho(0,0,0,1), 1.0/(spv_ref + spv(0,0,0,1)), &
2488 (rho_ref+rho(0,0,0,1)) * (spv_ref + spv(0,0,0,1)) - 1.0
2489 if (test_ok) then
2490 write(stdout,*) "The values of "//trim(eos_name)//" rho and 1/spv agree. "//trim(mesg)
2491 else
2492 write(stderr,*) "The values of "//trim(eos_name)//" rho and 1/spv disagree. "//trim(mesg)
2493 endif
2494 endif
2495 endif
2496
2497 ! Check that the densities are consistent when the reference value is extracted
2498 call calculate_density(t(0,0,0), s(0,0,0), p(0,0,0), rho_nooff, eos)
2499 test_ok = (abs(rho_nooff - (rho_ref + rho(0,0,0,1))) < tol*rho_nooff)
2500 ok = ok .and. test_ok
2501 if (verbose .and. .not.test_ok) then
2502 write(mesg, '(ES24.16," vs. ",ES24.16," with tolerance ",ES12.4)') &
2503 rho_ref+rho(0,0,0,1), rho_nooff, tol*rho_nooff
2504 write(stderr,*) "For "//trim(eos_name)//&
2505 " rho with and without a reference value disagree: "//trim(mesg)
2506 endif
2507
2508 ! Check that the specific volumes are consistent when the reference value is extracted
2509 call calculate_spec_vol(t(0,0,0), s(0,0,0), p(0,0,0), spv_nooff, eos)
2510 test_ok = (abs(spv_nooff - (spv_ref + spv(0,0,0,1))) < tol*rho_nooff)
2511 ok = ok .and. test_ok
2512 if (verbose .and. .not.test_ok) then
2513 write(mesg, '(ES24.16," vs. ",ES24.16," with tolerance ",ES12.4)') &
2514 spv_ref + spv(0,0,0,1), spv_nooff, tol*spv_nooff
2515 write(stderr,*) "For "//trim(eos_name)//&
2516 " spv with and without a reference value disagree: "//trim(mesg)
2517 endif
2518
2519 ! Account for the factors of terms in the numerator and denominator when estimating roundoff
2520 if (order == 6) then
2521 count_fac = 110.0/60.0 ; count_fac2 = 1088.0/180.0
2522 elseif (order == 4) then ! Use values appropriate for 4th order schemes.
2523 count_fac = 18.0/12.0 ; count_fac2 = 64.0/12.0
2524 else ! Use values appropriate for 2nd order schemes.
2525 count_fac = 2.0/2.0 ; count_fac2 = 4.0
2526 endif
2527
2528 ! Check for the rate of convergence expected with a 4th or 6th order accurate discretization
2529 ! with a 20% margin of error and a tolerance for contributions from roundoff.
2530 tol_here = tol*abs(drho_dt) + count_fac*r_tol/dt
2531 ok = ok .and. check_fd(drho_dt, drho_dt_fd, tol_here, verbose, trim(eos_name)//" drho_dT", order)
2532 tol_here = tol*abs(drho_ds) + count_fac*r_tol/ds
2533 ok = ok .and. check_fd(drho_ds, drho_ds_fd, tol_here, verbose, trim(eos_name)//" drho_dS", order)
2534 tol_here = tol*abs(drho_dp) + count_fac*r_tol/dp
2535 ok = ok .and. check_fd(drho_dp, drho_dp_fd, tol_here, verbose, trim(eos_name)//" drho_dp", order)
2536 tol_here = tol*abs(dsv_dt(1)) + count_fac*sv_tol/dt
2537 ok = ok .and. check_fd(dsv_dt(1), dsv_dt_fd, tol_here, verbose, trim(eos_name)//" dSV_dT", order)
2538 tol_here = tol*abs(dsv_ds(1)) + count_fac*sv_tol/ds
2539 ok = ok .and. check_fd(dsv_ds(1), dsv_ds_fd, tol_here, verbose, trim(eos_name)//" dSV_dS", order)
2540 if (test_2nd) then
2541 tol_here = tol*abs(drho_dt_dt) + count_fac2*r_tol/dt**2
2542 ok = ok .and. check_fd(drho_dt_dt, drho_dt_dt_fd, tol_here, verbose, trim(eos_name)//" drho_dT_dT", order)
2543 ! The curvature in salinity is relatively weak, so looser tolerances are needed for some forms of EOS?
2544 tol_here = 10.0*(tol*abs(drho_ds_ds) + count_fac2*r_tol/ds**2)
2545 ok = ok .and. check_fd(drho_ds_ds, drho_ds_ds_fd, tol_here, verbose, trim(eos_name)//" drho_dS_dS", order)
2546 tol_here = tol*abs(drho_ds_dt) + count_fac**2*r_tol/(ds*dt)
2547 ok = ok .and. check_fd(drho_ds_dt, drho_ds_dt_fd, tol_here, verbose, trim(eos_name)//" drho_dS_dT", order)
2548 tol_here = tol*abs(drho_dt_dp) + count_fac**2*r_tol/(dt*dp)
2549 ok = ok .and. check_fd(drho_dt_dp, drho_dt_dp_fd, tol_here, verbose, trim(eos_name)//" drho_dT_dP", order)
2550 tol_here = tol*abs(drho_ds_dp) + count_fac**2*r_tol/(ds*dp)
2551 ok = ok .and. check_fd(drho_ds_dp, drho_ds_dp_fd, tol_here, verbose, trim(eos_name)//" drho_dS_dP", order)
2552 endif
2553
2554 if (test_avg_sv) then
2555 tol_here = 0.5*tol*(abs(spv_avg_a(1)) + abs(spv_avg_q(1)))
2556 test_ok = (abs(spv_avg_a(1) - spv_avg_q(1)) < tol_here)
2557 if (verbose) then
2558 write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2,"), tol=",ES16.8)') &
2559 spv_avg_a(1), spv_avg_q(1), spv_avg_a(1) - spv_avg_q(1), &
2560 2.0*(spv_avg_a(1) - spv_avg_q(1)) / (abs(spv_avg_a(1)) + abs(spv_avg_q(1)) + tiny(spv_avg_a(1))), &
2561 tol_here
2562 if (verbose .and. .not.test_ok) then
2563 write(stderr,*) "The values of "//trim(eos_name)//" SpV_avg disagree. "//trim(mesg)
2564 elseif (verbose) then
2565 write(stdout,*) "The values of "//trim(eos_name)//" SpV_avg agree: "//trim(mesg)
2566 endif
2567 endif
2568 ok = ok .and. test_ok
2569 endif
2570
2571 inconsistent = .not.ok
2572
2573 contains
2574
2575 !> Return a finite difference estimate of the first derivative of a field in arbitrary units [A B-1]
2576 real function first_deriv(r, dx, order)
2577 real, intent(in) :: r(-3:3) !< The field whose derivative is being taken, in arbitrary units [A]
2578 real, intent(in) :: dx !< The spacing in parameter space, in different arbitrary units [B]
2579 integer, intent(in) :: order !< The order of accuracy of the centered finite difference estimates (2, 4 or 6)
2580
2581 if (order == 6) then ! Find a 6th order accurate first derivative on a regular grid.
2582 first_deriv = (45.0*(r(1)-r(-1)) + (-9.0*(r(2)-r(-2)) + (r(3)-r(-3))) ) / (60.0 * dx)
2583 elseif (order == 4) then ! Find a 4th order accurate first derivative on a regular grid.
2584 first_deriv = (8.0*(r(1)-r(-1)) - (r(2)-r(-2)) ) / (12.0 * dx)
2585 else ! Find a 2nd order accurate first derivative on a regular grid.
2586 first_deriv = (r(1)-r(-1)) / (2.0 * dx)
2587 endif
2588 end function first_deriv
2589
2590 !> Return a finite difference estimate of the second derivative of a field in arbitrary units [A B-2]
2591 real function second_deriv(r, dx, order)
2592 real, intent(in) :: r(-3:3) !< The field whose derivative is being taken, in arbitrary units [A]
2593 real, intent(in) :: dx !< The spacing in parameter space, in different arbitrary units [B]
2594 integer, intent(in) :: order !< The order of accuracy of the centered finite difference estimates (2, 4 or 6)
2595
2596 if (order == 6) then ! Find a 6th order accurate second derivative on a regular grid.
2597 second_deriv = ( -490.0*r(0) + (270.0*(r(1)+r(-1)) + (-27.0*(r(2)+r(-2)) + 2.0*(r(3)+r(-3))) )) / (180.0 * dx**2)
2598 elseif (order == 4) then ! Find a 4th order accurate second derivative on a regular grid.
2599 second_deriv = ( -30.0*r(0) + (16.0*(r(1)+r(-1)) - (r(2)+r(-2))) ) / (12.0 * dx**2)
2600 else ! Find a 2nd order accurate second derivative on a regular grid.
2601 second_deriv = ( -2.0*r(0) + (r(1)+r(-1)) ) / dx**2
2602 endif
2603 end function second_deriv
2604
2605 !> Return a finite difference estimate of the second derivative with respect to two different
2606 !! parameters of a field in arbitrary units [A B-1 C-1]
2607 real function derivs_2d(r, dxdy, order)
2608 real, intent(in) :: r(-3:3,-3:3) !< The field whose derivative is being taken in arbitrary units [A]
2609 real, intent(in) :: dxdy !< The spacing in two directions in parameter space in different arbitrary units [B C]
2610 integer, intent(in) :: order !< The order of accuracy of the centered finite difference estimates (2, 4 or 6)
2611
2612 real :: drdx(-3:3) ! The first derivative in one direction times the grid spacing in that direction [A]
2613 integer :: i
2614
2615 do i=-3,3
2616 drdx(i) = first_deriv(r(:,i), 1.0, order)
2617 enddo
2618 derivs_2d = first_deriv(drdx, dxdy, order)
2619
2620 end function derivs_2d
2621
2622 !> Check for the rate of convergence expected with a finite difference discretization
2623 !! with a 20% margin of error and a tolerance for contributions from roundoff.
2624 logical function check_fd(val, val_fd, tol, verbose, field_name, order)
2625 real, intent(in) :: val !< The derivative being checked, in arbitrary units [arbitrary]
2626 real, intent(in) :: val_fd(2) !< Two finite difference estimates of val taken with a spacing
2627 !! in parameter space and twice this spacing, in the same
2628 !! arbitrary units as val [arbitrary]
2629 real, intent(in) :: tol !< An estimated fractional tolerance due to roundoff [arbitrary]
2630 logical, intent(in) :: verbose !< If true, write results to stdout
2631 character(len=*), intent(in) :: field_name !< A name used to describe the field in error messages
2632 integer, intent(in) :: order !< The order of accuracy of the centered finite difference estimates (2, 4 or 6)
2633
2634 character(len=200) :: mesg
2635
2636 check_fd = ( abs(val_fd(1) - val) < (1.2*abs(val_fd(2) - val)/2**order + abs(tol)) )
2637
2638 write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2,"), tol=",ES16.8)') &
2639 val, val_fd(1), val - val_fd(1), &
2640 2.0*(val - val_fd(1)) / (abs(val) + abs(val_fd(1)) + tiny(val)), &
2641 (1.2*abs(val_fd(2) - val)/2**order + abs(tol))
2642 ! This message is useful for debugging the two estimates:
2643 ! write(mesg, '(ES16.8," and ",ES16.8," or ",ES16.8," differ by ",2ES16.8," (",2ES10.2"), tol=",ES16.8)') &
2644 ! val, val_fd(1), val_fd(2), val - val_fd(1), val - val_fd(2), &
2645 ! 2.0*(val - val_fd(1)) / (abs(val) + abs(val_fd(1)) + tiny(val)), &
2646 ! 2.0*(val - val_fd(2)) / (abs(val) + abs(val_fd(2)) + tiny(val)), &
2647 ! (1.2*abs(val_fd(2) - val)/2**order + abs(tol))
2648 if (verbose .and. .not.check_fd) then
2649 write(stderr,*) "The values of "//trim(field_name)//" disagree. "//trim(mesg)
2650 elseif (verbose) then
2651 write(stdout,*) "The values of "//trim(field_name)//" agree: "//trim(mesg)
2652 endif
2653 end function check_fd
2654
2655end function test_eos_consistency
2656
2657end module mom_eos
2658
2659!> \namespace mom_eos
2660!!
2661!! The MOM_EOS module is a wrapper for various equations of state (i.e. Linear, Wright,
2662!! Wright_full, Wright_red, UNESCO, TEOS10, Roquet_SpV or Roquet_rho) and provides a uniform
2663!! interface to the rest of the model independent of which equation of state is being used.