src/equation_of_state/MOM_EOS.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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