MOM6
MOM_EOS.F90
Go to the documentation of this file.
1 !> Provides subroutines for quantities specific to the equation of state
2 module mom_eos
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
27 use mom_eos_teos10, only : gsw_sp_from_sr, gsw_pt_from_ct
30 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg
33 use mom_hor_index, only : hor_index_type
34 
35 implicit none ; private
36 
37 #include <MOM_memory.h>
38 
46 public int_spec_vol_dp_generic_plm !, int_spec_vol_dz_generic_ppm
49 public calculate_tfreeze
51 public gsw_sp_from_sr, gsw_pt_from_ct
52 public extract_member_eos
53 
54 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
55 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
56 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
57 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
58 
59 !> Calculates density of sea water from T, S and P
62 end interface calculate_density
63 
64 !> Calculates specific volume of sea water from T, S and P
67 end interface calculate_spec_vol
68 
69 !> Calculate the derivatives of density with temperature and salinity from T, S, and P
72 end interface calculate_density_derivs
73 
74 !> Calculates the second derivatives of density with various combinations of temperature,
75 !! salinity, and pressure from T, S and P
79 
80 !> Calculates the freezing point of sea water from T, S and P
83 end interface calculate_tfreeze
84 
85 !> Calculates the compressibility of water from T, S, and P
88 end interface calculate_compress
89 
90 !> A control structure for the equation of state
91 type, public :: eos_type ; private
92  integer :: form_of_eos = 0 !< The equation of state to use.
93  integer :: form_of_tfreeze = 0 !< The expression for the potential temperature
94  !! of the freezing point.
95  logical :: eos_quadrature !< If true, always use the generic (quadrature)
96  !! code for the integrals of density.
97  logical :: compressible = .true. !< If true, in situ density is a function of pressure.
98 ! The following parameters are used with the linear equation of state only.
99  real :: rho_t0_s0 !< The density at T=0, S=0 [kg m-3].
100  real :: drho_dt !< The partial derivative of density with temperature [kg m-3 degC-1]
101  real :: drho_ds !< The partial derivative of density with salinity [kg m-3 ppt-1].
102 ! The following parameters are use with the linear expression for the freezing
103 ! point only.
104  real :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC].
105  real :: dtfr_ds !< The derivative of freezing point with salinity [degC ppt-1].
106  real :: dtfr_dp !< The derivative of freezing point with pressure [degC Pa-1].
107 
108 ! logical :: test_EOS = .true. ! If true, test the equation of state
109 end type eos_type
110 
111 ! The named integers that might be stored in eqn_of_state_type%form_of_EOS.
112 integer, parameter, public :: eos_linear = 1 !< A named integer specifying an equation of state
113 integer, parameter, public :: eos_unesco = 2 !< A named integer specifying an equation of state
114 integer, parameter, public :: eos_wright = 3 !< A named integer specifying an equation of state
115 integer, parameter, public :: eos_teos10 = 4 !< A named integer specifying an equation of state
116 integer, parameter, public :: eos_nemo = 5 !< A named integer specifying an equation of state
117 
118 character*(10), parameter :: eos_linear_string = "LINEAR" !< A string for specifying the equation of state
119 character*(10), parameter :: eos_unesco_string = "UNESCO" !< A string for specifying the equation of state
120 character*(10), parameter :: eos_wright_string = "WRIGHT" !< A string for specifying the equation of state
121 character*(10), parameter :: eos_teos10_string = "TEOS10" !< A string for specifying the equation of state
122 character*(10), parameter :: eos_nemo_string = "NEMO" !< A string for specifying the equation of state
123 character*(10), parameter :: eos_default = eos_wright_string !< The default equation of state
124 
125 integer, parameter :: tfreeze_linear = 1 !< A named integer specifying a freezing point expression
126 integer, parameter :: tfreeze_millero = 2 !< A named integer specifying a freezing point expression
127 integer, parameter :: tfreeze_teos10 = 3 !< A named integer specifying a freezing point expression
128 character*(10), parameter :: tfreeze_linear_string = "LINEAR" !< A string for specifying the freezing point expression
129 character*(10), parameter :: tfreeze_millero_string = "MILLERO_78" !< A string for specifying
130  !! freezing point expression
131 character*(10), parameter :: tfreeze_teos10_string = "TEOS10" !< A string for specifying the freezing point expression
132 character*(10), parameter :: tfreeze_default = tfreeze_linear_string !< The default freezing point expression
133 
134 contains
135 
136 !> Calls the appropriate subroutine to calculate density of sea water for scalar inputs.
137 !! If rho_ref is present, the anomaly with respect to rho_ref is returned.
138 subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, scale)
139  real, intent(in) :: T !< Potential temperature referenced to the surface [degC]
140  real, intent(in) :: S !< Salinity [ppt]
141  real, intent(in) :: pressure !< Pressure [Pa]
142  real, intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3]
143  type(eos_type), pointer :: EOS !< Equation of state structure
144  real, optional, intent(in) :: rho_ref !< A reference density [kg m-3].
145  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
146  !! from kg m-3 to the desired units [R m3 kg-1]
147 
148  if (.not.associated(eos)) call mom_error(fatal, &
149  "calculate_density_scalar called with an unassociated EOS_type EOS.")
150 
151  select case (eos%form_of_EOS)
152  case (eos_linear)
153  call calculate_density_linear(t, s, pressure, rho, &
154  eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
155  case (eos_unesco)
156  call calculate_density_unesco(t, s, pressure, rho, rho_ref)
157  case (eos_wright)
158  call calculate_density_wright(t, s, pressure, rho, rho_ref)
159  case (eos_teos10)
160  call calculate_density_teos10(t, s, pressure, rho, rho_ref)
161  case (eos_nemo)
162  call calculate_density_nemo(t, s, pressure, rho, rho_ref)
163  case default
164  call mom_error(fatal, &
165  "calculate_density_scalar: EOS is not valid.")
166  end select
167 
168  if (present(scale)) then ; if (scale /= 1.0) then
169  rho = scale * rho
170  endif ; endif
171 
172 end subroutine calculate_density_scalar
173 
174 !> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs.
175 !! If rho_ref is present, the anomaly with respect to rho_ref is returned.
176 subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_ref, scale)
177  real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC]
178  real, dimension(:), intent(in) :: S !< Salinity [ppt]
179  real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
180  real, dimension(:), intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3]
181  integer, intent(in) :: start !< Start index for computation
182  integer, intent(in) :: npts !< Number of point to compute
183  type(eos_type), pointer :: EOS !< Equation of state structure
184  real, optional, intent(in) :: rho_ref !< A reference density [kg m-3].
185  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
186  !! from kg m-3 to the desired units [R m3 kg-1]
187  integer :: j
188 
189  if (.not.associated(eos)) call mom_error(fatal, &
190  "calculate_density_array called with an unassociated EOS_type EOS.")
191 
192  select case (eos%form_of_EOS)
193  case (eos_linear)
194  call calculate_density_linear(t, s, pressure, rho, start, npts, &
195  eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
196  case (eos_unesco)
197  call calculate_density_unesco(t, s, pressure, rho, start, npts, rho_ref)
198  case (eos_wright)
199  call calculate_density_wright(t, s, pressure, rho, start, npts, rho_ref)
200  case (eos_teos10)
201  call calculate_density_teos10(t, s, pressure, rho, start, npts, rho_ref)
202  case (eos_nemo)
203  call calculate_density_nemo (t, s, pressure, rho, start, npts, rho_ref)
204  case default
205  call mom_error(fatal, &
206  "calculate_density_array: EOS%form_of_EOS is not valid.")
207  end select
208 
209  if (present(scale)) then ; if (scale /= 1.0) then
210  do j=start,start+npts-1 ; rho(j) = scale * rho(j) ; enddo
211  endif ; endif
212 
213 end subroutine calculate_density_array
214 
215 !> Calls the appropriate subroutine to calculate specific volume of sea water
216 !! for scalar inputs.
217 subroutine calculate_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, scale)
218  real, intent(in) :: T !< Potential temperature referenced to the surface [degC]
219  real, intent(in) :: S !< Salinity [ppt]
220  real, intent(in) :: pressure !< Pressure [Pa]
221  real, intent(out) :: specvol !< In situ? specific volume [m3 kg-1] or [R-1 ~> m3 kg-1]
222  type(eos_type), pointer :: EOS !< Equation of state structure
223  real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
224  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific volume
225  !! from m3 kg-1 to the desired units [kg m-3 R-1]
226 
227  real :: rho
228 
229  if (.not.associated(eos)) call mom_error(fatal, &
230  "calculate_spec_vol_scalar called with an unassociated EOS_type EOS.")
231 
232  select case (eos%form_of_EOS)
233  case (eos_linear)
234  call calculate_spec_vol_linear(t, s, pressure, specvol, &
235  eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
236  case (eos_unesco)
237  call calculate_spec_vol_unesco(t, s, pressure, specvol, spv_ref)
238  case (eos_wright)
239  call calculate_spec_vol_wright(t, s, pressure, specvol, spv_ref)
240  case (eos_teos10)
241  call calculate_spec_vol_teos10(t, s, pressure, specvol, spv_ref)
242  case (eos_nemo)
243  call calculate_density_nemo(t, s, pressure, rho)
244  if (present(spv_ref)) then
245  specvol = 1.0 / rho - spv_ref
246  else
247  specvol = 1.0 / rho
248  endif
249  case default
250  call mom_error(fatal, &
251  "calculate_spec_vol_scalar: EOS is not valid.")
252  end select
253 
254  if (present(scale)) then ; if (scale /= 1.0) then
255  specvol = scale * specvol
256  endif ; endif
257 
258 end subroutine calculate_spec_vol_scalar
259 
260 
261 !> Calls the appropriate subroutine to calculate the specific volume of sea water
262 !! for 1-D array inputs.
263 subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref, scale)
264  real, dimension(:), intent(in) :: T !< potential temperature relative to the surface
265  !! [degC].
266  real, dimension(:), intent(in) :: S !< salinity [ppt].
267  real, dimension(:), intent(in) :: pressure !< pressure [Pa].
268  real, dimension(:), intent(out) :: specvol !< in situ specific volume [kg m-3] or [R-1 ~> m3 kg-1].
269  integer, intent(in) :: start !< the starting point in the arrays.
270  integer, intent(in) :: npts !< the number of values to calculate.
271  type(eos_type), pointer :: EOS !< Equation of state structure
272  real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1].
273  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific volume
274  !! from m3 kg-1 to the desired units [kg m-3 R-1]
275 
276  real, dimension(size(specvol)) :: rho
277  integer :: j
278 
279  if (.not.associated(eos)) call mom_error(fatal, &
280  "calculate_spec_vol_array called with an unassociated EOS_type EOS.")
281 
282  select case (eos%form_of_EOS)
283  case (eos_linear)
284  call calculate_spec_vol_linear(t, s, pressure, specvol, start, npts, &
285  eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
286  case (eos_unesco)
287  call calculate_spec_vol_unesco(t, s, pressure, specvol, start, npts, spv_ref)
288  case (eos_wright)
289  call calculate_spec_vol_wright(t, s, pressure, specvol, start, npts, spv_ref)
290  case (eos_teos10)
291  call calculate_spec_vol_teos10(t, s, pressure, specvol, start, npts, spv_ref)
292  case (eos_nemo)
293  call calculate_density_nemo (t, s, pressure, rho, start, npts)
294  if (present(spv_ref)) then
295  specvol(:) = 1.0 / rho(:) - spv_ref
296  else
297  specvol(:) = 1.0 / rho(:)
298  endif
299  case default
300  call mom_error(fatal, &
301  "calculate_spec_vol_array: EOS%form_of_EOS is not valid.")
302  end select
303 
304  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
305  specvol(j) = scale * specvol(j)
306  enddo ; endif ; endif
307 
308 end subroutine calculate_spec_vol_array
309 
310 
311 !> Calls the appropriate subroutine to calculate the freezing point for scalar inputs.
312 subroutine calculate_tfreeze_scalar(S, pressure, T_fr, EOS)
313  real, intent(in) :: S !< Salinity [ppt]
314  real, intent(in) :: pressure !< Pressure [Pa]
315  real, intent(out) :: T_fr !< Freezing point potential temperature referenced
316  !! to the surface [degC]
317  type(eos_type), pointer :: EOS !< Equation of state structure
318 
319  if (.not.associated(eos)) call mom_error(fatal, &
320  "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
321 
322  select case (eos%form_of_TFreeze)
323  case (tfreeze_linear)
324  call calculate_tfreeze_linear(s, pressure, t_fr, eos%TFr_S0_P0, &
325  eos%dTFr_dS, eos%dTFr_dp)
326  case (tfreeze_millero)
327  call calculate_tfreeze_millero(s, pressure, t_fr)
328  case (tfreeze_teos10)
329  call calculate_tfreeze_teos10(s, pressure, t_fr)
330  case default
331  call mom_error(fatal, &
332  "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
333  end select
334 
335 end subroutine calculate_tfreeze_scalar
336 
337 !> Calls the appropriate subroutine to calculate the freezing point for a 1-D array.
338 subroutine calculate_tfreeze_array(S, pressure, T_fr, start, npts, EOS)
339  real, dimension(:), intent(in) :: S !< Salinity [ppt]
340  real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
341  real, dimension(:), intent(out) :: T_fr !< Freezing point potential temperature referenced
342  !! to the surface [degC]
343  integer, intent(in) :: start !< Starting index within the array
344  integer, intent(in) :: npts !< The number of values to calculate
345  type(eos_type), pointer :: EOS !< Equation of state structure
346 
347  if (.not.associated(eos)) call mom_error(fatal, &
348  "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
349 
350  select case (eos%form_of_TFreeze)
351  case (tfreeze_linear)
352  call calculate_tfreeze_linear(s, pressure, t_fr, start, npts, &
353  eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
354  case (tfreeze_millero)
355  call calculate_tfreeze_millero(s, pressure, t_fr, start, npts)
356  case (tfreeze_teos10)
357  call calculate_tfreeze_teos10(s, pressure, t_fr, start, npts)
358  case default
359  call mom_error(fatal, &
360  "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
361  end select
362 
363 end subroutine calculate_tfreeze_array
364 
365 !> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.
366 subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS, scale)
367  real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC]
368  real, dimension(:), intent(in) :: S !< Salinity [ppt]
369  real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
370  real, dimension(:), intent(out) :: drho_dT !< The partial derivative of density with potential
371  !! temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1].
372  real, dimension(:), intent(out) :: drho_dS !< The partial derivative of density with salinity,
373  !! in [kg m-3 ppt-1] or [R degC-1 ~> kg m-3 ppt-1].
374  integer, intent(in) :: start !< Starting index within the array
375  integer, intent(in) :: npts !< The number of values to calculate
376  type(eos_type), pointer :: EOS !< Equation of state structure
377  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
378  !! from kg m-3 to the desired units [R m3 kg-1]
379  integer :: j
380 
381  if (.not.associated(eos)) call mom_error(fatal, &
382  "calculate_density_derivs called with an unassociated EOS_type EOS.")
383 
384  select case (eos%form_of_EOS)
385  case (eos_linear)
386  call calculate_density_derivs_linear(t, s, pressure, drho_dt, drho_ds, eos%Rho_T0_S0, &
387  eos%dRho_dT, eos%dRho_dS, start, npts)
388  case (eos_unesco)
389  call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
390  case (eos_wright)
391  call calculate_density_derivs_wright(t, s, pressure, drho_dt, drho_ds, start, npts)
392  case (eos_teos10)
393  call calculate_density_derivs_teos10(t, s, pressure, drho_dt, drho_ds, start, npts)
394  case (eos_nemo)
395  call calculate_density_derivs_nemo(t, s, pressure, drho_dt, drho_ds, start, npts)
396  case default
397  call mom_error(fatal, &
398  "calculate_density_derivs_array: EOS%form_of_EOS is not valid.")
399  end select
400 
401  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
402  drho_dt(j) = scale * drho_dt(j)
403  drho_ds(j) = scale * drho_ds(j)
404  enddo ; endif ; endif
405 
406 end subroutine calculate_density_derivs_array
407 
408 !> Calls the appropriate subroutines to calculate density derivatives by promoting a scalar
409 !! to a one-element array
410 subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS, scale)
411  real, intent(in) :: T !< Potential temperature referenced to the surface [degC]
412  real, intent(in) :: S !< Salinity [ppt]
413  real, intent(in) :: pressure !< Pressure [Pa]
414  real, intent(out) :: drho_dT !< The partial derivative of density with potential
415  !! temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1]
416  real, intent(out) :: drho_dS !< The partial derivative of density with salinity,
417  !! in [kg m-3 ppt-1] or [R ppt-1 ~> kg m-3 ppt-1].
418  type(eos_type), pointer :: EOS !< Equation of state structure
419  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
420  !! from kg m-3 to the desired units [R m3 kg-1]
421 
422  if (.not.associated(eos)) call mom_error(fatal, &
423  "calculate_density_derivs called with an unassociated EOS_type EOS.")
424 
425  select case (eos%form_of_EOS)
426  case (eos_linear)
427  call calculate_density_derivs_linear(t, s, pressure, drho_dt, drho_ds, &
428  eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
429  case (eos_wright)
430  call calculate_density_derivs_wright(t, s, pressure, drho_dt, drho_ds)
431  case (eos_teos10)
432  call calculate_density_derivs_teos10(t, s, pressure, drho_dt, drho_ds)
433  case default
434  call mom_error(fatal, &
435  "calculate_density_derivs_scalar: EOS%form_of_EOS is not valid.")
436  end select
437 
438  if (present(scale)) then ; if (scale /= 1.0) then
439  drho_dt = scale * drho_dt
440  drho_ds = scale * drho_ds
441  endif ; endif
442 
443 end subroutine calculate_density_derivs_scalar
444 
445 !> Calls the appropriate subroutine to calculate density second derivatives for 1-D array inputs.
446 subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, &
447  drho_dS_dP, drho_dT_dP, start, npts, EOS, scale)
448  real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC]
449  real, dimension(:), intent(in) :: S !< Salinity [ppt]
450  real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
451  real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S
452  !! [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2]
453  real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with respect to T
454  !! [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1]
455  real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T
456  !! [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2]
457  real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure
458  !! [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1]
459  real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure
460  !! [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1]
461  integer, intent(in) :: start !< Starting index within the array
462  integer, intent(in) :: npts !< The number of values to calculate
463  type(eos_type), pointer :: EOS !< Equation of state structure
464  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
465  !! from kg m-3 to the desired units [R m3 kg-1]
466  integer :: j
467 
468  if (.not.associated(eos)) call mom_error(fatal, &
469  "calculate_density_derivs called with an unassociated EOS_type EOS.")
470 
471  select case (eos%form_of_EOS)
472  case (eos_linear)
473  call calculate_density_second_derivs_linear(t, s, pressure, drho_ds_ds, drho_ds_dt, &
474  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
475  case (eos_wright)
476  call calculate_density_second_derivs_wright(t, s, pressure, drho_ds_ds, drho_ds_dt, &
477  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
478  case (eos_teos10)
479  call calculate_density_second_derivs_teos10(t, s, pressure, drho_ds_ds, drho_ds_dt, &
480  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
481  case default
482  call mom_error(fatal, &
483  "calculate_density_derivs: EOS%form_of_EOS is not valid.")
484  end select
485 
486  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
487  drho_ds_ds(j) = scale * drho_ds_ds(j)
488  drho_ds_dt(j) = scale * drho_ds_dt(j)
489  drho_dt_dt(j) = scale * drho_dt_dt(j)
490  drho_ds_dp(j) = scale * drho_ds_dp(j)
491  drho_dt_dp(j) = scale * drho_dt_dp(j)
492  enddo ; endif ; endif
493 
495 
496 !> Calls the appropriate subroutine to calculate density second derivatives for scalar nputs.
497 subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, &
498  drho_dS_dP, drho_dT_dP, EOS, scale)
499  real, intent(in) :: T !< Potential temperature referenced to the surface [degC]
500  real, intent(in) :: S !< Salinity [ppt]
501  real, intent(in) :: pressure !< Pressure [Pa]
502  real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S
503  !! [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2]
504  real, intent(out) :: drho_dS_dT !< Partial derivative of beta with respect to T
505  !! [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1]
506  real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T
507  !! [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2]
508  real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure
509  !! [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1]
510  real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure
511  !! [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1]
512  type(eos_type), pointer :: EOS !< Equation of state structure
513  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
514  !! from kg m-3 to the desired units [R m3 kg-1]
515 
516  if (.not.associated(eos)) call mom_error(fatal, &
517  "calculate_density_derivs called with an unassociated EOS_type EOS.")
518 
519  select case (eos%form_of_EOS)
520  case (eos_linear)
521  call calculate_density_second_derivs_linear(t, s, pressure, drho_ds_ds, drho_ds_dt, &
522  drho_dt_dt, drho_ds_dp, drho_dt_dp)
523  case (eos_wright)
524  call calculate_density_second_derivs_wright(t, s, pressure, drho_ds_ds, drho_ds_dt, &
525  drho_dt_dt, drho_ds_dp, drho_dt_dp)
526  case (eos_teos10)
527  call calculate_density_second_derivs_teos10(t, s, pressure, drho_ds_ds, drho_ds_dt, &
528  drho_dt_dt, drho_ds_dp, drho_dt_dp)
529  case default
530  call mom_error(fatal, &
531  "calculate_density_derivs: EOS%form_of_EOS is not valid.")
532  end select
533 
534  if (present(scale)) then ; if (scale /= 1.0) then
535  drho_ds_ds = scale * drho_ds_ds
536  drho_ds_dt = scale * drho_ds_dt
537  drho_dt_dt = scale * drho_dt_dt
538  drho_ds_dp = scale * drho_ds_dp
539  drho_dt_dp = scale * drho_dt_dp
540  endif ; endif
541 
543 
544 !> Calls the appropriate subroutine to calculate specific volume derivatives for an array.
545 subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS, scale)
546  real, dimension(:), intent(in) :: t !< Potential temperature referenced to the surface [degC]
547  real, dimension(:), intent(in) :: s !< Salinity [ppt]
548  real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
549  real, dimension(:), intent(out) :: dsv_dt !< The partial derivative of specific volume with potential
550  !! temperature [m3 kg-1 degC-1] or [R-1 degC-1 ~> m3 kg-1 degC-1]
551  real, dimension(:), intent(out) :: dsv_ds !< The partial derivative of specific volume with salinity
552  !! [m3 kg-1 ppt-1] or [R-1 ppt-1 ~> m3 kg-1 ppt-1].
553  integer, intent(in) :: start !< Starting index within the array
554  integer, intent(in) :: npts !< The number of values to calculate
555  type(eos_type), pointer :: eos !< Equation of state structure
556  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific volume
557  !! from m3 kg-1 to the desired units [kg m-3 R-1]
558 
559  ! Local variables
560  real, dimension(size(T)) :: drho_dt, drho_ds, rho
561  integer :: j
562 
563  if (.not.associated(eos)) call mom_error(fatal, &
564  "calculate_density_derivs called with an unassociated EOS_type EOS.")
565 
566  select case (eos%form_of_EOS)
567  case (eos_linear)
568  call calculate_specvol_derivs_linear(t, s, pressure, dsv_dt, dsv_ds, start, &
569  npts, eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
570  case (eos_unesco)
571  call calculate_density_unesco(t, s, pressure, rho, start, npts)
572  call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
573  do j=start,start+npts-1
574  dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
575  dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
576  enddo
577  case (eos_wright)
578  call calculate_specvol_derivs_wright(t, s, pressure, dsv_dt, dsv_ds, start, npts)
579  case (eos_teos10)
580  call calculate_specvol_derivs_teos10(t, s, pressure, dsv_dt, dsv_ds, start, npts)
581  case (eos_nemo)
582  call calculate_density_nemo(t, s, pressure, rho, start, npts)
583  call calculate_density_derivs_nemo(t, s, pressure, drho_dt, drho_ds, start, npts)
584  do j=start,start+npts-1
585  dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
586  dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
587  enddo
588  case default
589  call mom_error(fatal, &
590  "calculate_density_derivs: EOS%form_of_EOS is not valid.")
591  end select
592 
593  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
594  dsv_dt(j) = scale * dsv_dt(j)
595  dsv_ds(j) = scale * dsv_ds(j)
596  enddo ; endif ; endif
597 
598 
599 end subroutine calculate_specific_vol_derivs
600 
601 !> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs.
602 subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, EOS)
603  real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC]
604  real, dimension(:), intent(in) :: S !< Salinity [PSU]
605  real, dimension(:), intent(in) :: pressure !< Pressure [Pa]
606  real, dimension(:), intent(out) :: rho !< In situ density [kg m-3].
607  real, dimension(:), intent(out) :: drho_dp !< The partial derivative of density with pressure
608  !! (also the inverse of the square of sound speed) [s2 m-2].
609  integer, intent(in) :: start !< Starting index within the array
610  integer, intent(in) :: npts !< The number of values to calculate
611  type(eos_type), pointer :: EOS !< Equation of state structure
612 
613  if (.not.associated(eos)) call mom_error(fatal, &
614  "calculate_compress called with an unassociated EOS_type EOS.")
615 
616  select case (eos%form_of_EOS)
617  case (eos_linear)
618  call calculate_compress_linear(t, s, pressure, rho, drho_dp, start, npts, &
619  eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
620  case (eos_unesco)
621  call calculate_compress_unesco(t, s, pressure, rho, drho_dp, start, npts)
622  case (eos_wright)
623  call calculate_compress_wright(t, s, pressure, rho, drho_dp, start, npts)
624  case (eos_teos10)
625  call calculate_compress_teos10(t, s, pressure, rho, drho_dp, start, npts)
626  case (eos_nemo)
627  call calculate_compress_nemo(t, s, pressure, rho, drho_dp, start, npts)
628  case default
629  call mom_error(fatal, &
630  "calculate_compress: EOS%form_of_EOS is not valid.")
631  end select
632 
633 end subroutine calculate_compress_array
634 
635 !> Calculate density and compressibility for a scalar. This just promotes the scalar to an array with a singleton
636 !! dimension and calls calculate_compress_array
637 subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS)
638  real, intent(in) :: T !< Potential temperature referenced to the surface (degC)
639  real, intent(in) :: S !< Salinity (PSU)
640  real, intent(in) :: pressure !< Pressure (Pa)
641  real, intent(out) :: rho !< In situ density in kg m-3.
642  real, intent(out) :: drho_dp !< The partial derivative of density with pressure
643  !! (also the inverse of the square of sound speed) in s2 m-2.
644  type(eos_type), pointer :: EOS !< Equation of state structure
645 
646  real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa
647 
648  if (.not.associated(eos)) call mom_error(fatal, &
649  "calculate_compress called with an unassociated EOS_type EOS.")
650  ta(1) = t ; sa(1) = s; pa(1) = pressure
651 
652  call calculate_compress_array(ta, sa, pa, rhoa, drho_dpa, 1, 1, eos)
653  rho = rhoa(1) ; drho_dp = drho_dpa(1)
654 
655 end subroutine calculate_compress_scalar
656 !> Calls the appropriate subroutine to alculate analytical and nearly-analytical
657 !! integrals in pressure across layers of geopotential anomalies, which are
658 !! required for calculating the finite-volume form pressure accelerations in a
659 !! non-Boussinesq model. There are essentially no free assumptions, apart from the
660 !! use of Bode's rule to do the horizontal integrals, and from a truncation in the
661 !! series for log(1-eps/1+eps) that assumes that |eps| < .
662 subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, &
663  dza, intp_dza, intx_dza, inty_dza, halo_size, &
664  bathyP, dP_tiny, useMassWghtInterp)
665  type(hor_index_type), intent(in) :: hi !< The horizontal index structure
666  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
667  intent(in) :: t !< Potential temperature referenced to the surface [degC]
668  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
669  intent(in) :: s !< Salinity [ppt]
670  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
671  intent(in) :: p_t !< Pressure at the top of the layer [Pa].
672  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
673  intent(in) :: p_b !< Pressure at the bottom of the layer [Pa].
674  real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out
675  !! to reduce the magnitude of each of the integrals, m3 kg-1. The
676  !! calculation is mathematically identical with different values of
677  !! alpha_ref, but this reduces the effects of roundoff.
678  type(eos_type), pointer :: eos !< Equation of state structure
679  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
680  intent(out) :: dza !< The change in the geopotential anomaly across
681  !! the layer [m2 s-2].
682  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
683  optional, intent(out) :: intp_dza !< The integral in pressure through the layer of the
684  !! geopotential anomaly relative to the anomaly at the bottom of the
685  !! layer [Pa m2 s-2].
686  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
687  optional, intent(out) :: intx_dza !< The integral in x of the difference between the
688  !! geopotential anomaly at the top and bottom of the layer divided by
689  !! the x grid spacing [m2 s-2].
690  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
691  optional, intent(out) :: inty_dza !< The integral in y of the difference between the
692  !! geopotential anomaly at the top and bottom of the layer divided by
693  !! the y grid spacing [m2 s-2].
694  integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
695  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
696  optional, intent(in) :: bathyp !< The pressure at the bathymetry [Pa]
697  real, optional, intent(in) :: dp_tiny !< A miniscule pressure change with
698  !! the same units as p_t (Pa?)
699  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
700  !! to interpolate T/S for top and bottom integrals.
701 
702  if (.not.associated(eos)) call mom_error(fatal, &
703  "int_specific_vol_dp called with an unassociated EOS_type EOS.")
704 
705  if (eos%EOS_quadrature) then
706  call int_spec_vol_dp_generic(t, s, p_t, p_b, alpha_ref, hi, eos, &
707  dza, intp_dza, intx_dza, inty_dza, halo_size, &
708  bathyp, dp_tiny, usemasswghtinterp)
709  else ; select case (eos%form_of_EOS)
710  case (eos_linear)
711  call int_spec_vol_dp_linear(t, s, p_t, p_b, alpha_ref, hi, eos%Rho_T0_S0, &
712  eos%dRho_dT, eos%dRho_dS, dza, intp_dza, &
713  intx_dza, inty_dza, halo_size, &
714  bathyp, dp_tiny, usemasswghtinterp)
715  case (eos_wright)
716  call int_spec_vol_dp_wright(t, s, p_t, p_b, alpha_ref, hi, dza, &
717  intp_dza, intx_dza, inty_dza, halo_size, &
718  bathyp, dp_tiny, usemasswghtinterp)
719  case default
720  call int_spec_vol_dp_generic(t, s, p_t, p_b, alpha_ref, hi, eos, &
721  dza, intp_dza, intx_dza, inty_dza, halo_size, &
722  bathyp, dp_tiny, usemasswghtinterp)
723  end select ; endif
724 
725 end subroutine int_specific_vol_dp
726 
727 !> This subroutine calculates analytical and nearly-analytical integrals of
728 !! pressure anomalies across layers, which are required for calculating the
729 !! finite-volume form pressure accelerations in a Boussinesq model.
730 subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, &
731  dpa, intz_dpa, intx_dpa, inty_dpa, &
732  bathyT, dz_neglect, useMassWghtInterp)
733  type(hor_index_type), intent(in) :: hii !< Ocean horizontal index structures for the input arrays
734  type(hor_index_type), intent(in) :: hio !< Ocean horizontal index structures for the output arrays
735  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
736  intent(in) :: t !< Potential temperature referenced to the surface [degC]
737  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
738  intent(in) :: s !< Salinity [ppt]
739  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
740  intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m].
741  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
742  intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m].
743  real, intent(in) :: rho_ref !< A mean density [kg m-3], that is subtracted out to
744  !! reduce the magnitude of each of the integrals.
745  real, intent(in) :: rho_0 !< A density [kg m-3], that is used to calculate the
746  !! pressure (as p~=-z*rho_0*G_e) used in the equation of state.
747  real, intent(in) :: g_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
748  type(eos_type), pointer :: eos !< Equation of state structure
749  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
750  intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
751  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
752  optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
753  !! the pressure anomaly relative to the anomaly at the
754  !! top of the layer [Pa Z ~> Pa m].
755  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
756  optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
757  !! pressure anomaly at the top and bottom of the layer
758  !! divided by the x grid spacing [Pa].
759  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
760  optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
761  !! pressure anomaly at the top and bottom of the layer
762  !! divided by the y grid spacing [Pa].
763  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
764  optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m].
765  real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m].
766  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
767  !! interpolate T/S for top and bottom integrals.
768 
769  if (.not.associated(eos)) call mom_error(fatal, &
770  "int_density_dz called with an unassociated EOS_type EOS.")
771 
772  if (eos%EOS_quadrature) then
773  call int_density_dz_generic(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
774  eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
775  bathyt, dz_neglect, usemasswghtinterp)
776  else ; select case (eos%form_of_EOS)
777  case (eos_linear)
778  call int_density_dz_linear(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
779  eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, &
780  dpa, intz_dpa, intx_dpa, inty_dpa, &
781  bathyt, dz_neglect, usemasswghtinterp)
782  case (eos_wright)
783  call int_density_dz_wright(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
784  dpa, intz_dpa, intx_dpa, inty_dpa, &
785  bathyt, dz_neglect, usemasswghtinterp)
786  case default
787  call int_density_dz_generic(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
788  eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
789  bathyt, dz_neglect, usemasswghtinterp)
790  end select ; endif
791 
792 end subroutine int_density_dz
793 
794 !> Returns true if the equation of state is compressible (i.e. has pressure dependence)
795 logical function query_compressible(EOS)
796  type(eos_type), pointer :: eos !< Equation of state structure
797 
798  if (.not.associated(eos)) call mom_error(fatal, &
799  "query_compressible called with an unassociated EOS_type EOS.")
800 
801  query_compressible = eos%compressible
802 end function query_compressible
803 
804 !> Initializes EOS_type by allocating and reading parameters
805 subroutine eos_init(param_file, EOS)
806  type(param_file_type), intent(in) :: param_file !< Parameter file structure
807  type(eos_type), pointer :: eos !< Equation of state structure
808  ! Local variables
809 #include "version_variable.h"
810  character(len=40) :: mdl = "MOM_EOS" ! This module's name.
811  character(len=40) :: tmpstr
812 
813  if (.not.associated(eos)) call eos_allocate(eos)
814 
815  ! Read all relevant parameters and write them to the model log.
816  call log_version(param_file, mdl, version, "")
817 
818  call get_param(param_file, mdl, "EQN_OF_STATE", tmpstr, &
819  "EQN_OF_STATE determines which ocean equation of state "//&
820  "should be used. Currently, the valid choices are "//&
821  '"LINEAR", "UNESCO", "WRIGHT", "NEMO" and "TEOS10". '//&
822  "This is only used if USE_EOS is true.", default=eos_default)
823  select case (uppercase(tmpstr))
824  case (eos_linear_string)
825  eos%form_of_EOS = eos_linear
826  case (eos_unesco_string)
827  eos%form_of_EOS = eos_unesco
828  case (eos_wright_string)
829  eos%form_of_EOS = eos_wright
830  case (eos_teos10_string)
831  eos%form_of_EOS = eos_teos10
832  case (eos_nemo_string)
833  eos%form_of_EOS = eos_nemo
834  case default
835  call mom_error(fatal, "interpret_eos_selection: EQN_OF_STATE "//&
836  trim(tmpstr) // "in input file is invalid.")
837  end select
838  call mom_mesg('interpret_eos_selection: equation of state set to "' // &
839  trim(tmpstr)//'"', 5)
840 
841  if (eos%form_of_EOS == eos_linear) then
842  eos%Compressible = .false.
843  call get_param(param_file, mdl, "RHO_T0_S0", eos%Rho_T0_S0, &
844  "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
845  "this is the density at T=0, S=0.", units="kg m-3", &
846  default=1000.0)
847  call get_param(param_file, mdl, "DRHO_DT", eos%dRho_dT, &
848  "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
849  "this is the partial derivative of density with "//&
850  "temperature.", units="kg m-3 K-1", default=-0.2)
851  call get_param(param_file, mdl, "DRHO_DS", eos%dRho_dS, &
852  "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
853  "this is the partial derivative of density with "//&
854  "salinity.", units="kg m-3 PSU-1", default=0.8)
855  endif
856 
857  call get_param(param_file, mdl, "EOS_QUADRATURE", eos%EOS_quadrature, &
858  "If true, always use the generic (quadrature) code "//&
859  "code for the integrals of density.", default=.false.)
860 
861  call get_param(param_file, mdl, "TFREEZE_FORM", tmpstr, &
862  "TFREEZE_FORM determines which expression should be "//&
863  "used for the freezing point. Currently, the valid "//&
864  'choices are "LINEAR", "MILLERO_78", "TEOS10"', &
865  default=tfreeze_default)
866  select case (uppercase(tmpstr))
867  case (tfreeze_linear_string)
868  eos%form_of_TFreeze = tfreeze_linear
870  eos%form_of_TFreeze = tfreeze_millero
871  case (tfreeze_teos10_string)
872  eos%form_of_TFreeze = tfreeze_teos10
873  case default
874  call mom_error(fatal, "interpret_eos_selection: TFREEZE_FORM "//&
875  trim(tmpstr) // "in input file is invalid.")
876  end select
877 
878  if (eos%form_of_TFreeze == tfreeze_linear) then
879  call get_param(param_file, mdl, "TFREEZE_S0_P0",eos%TFr_S0_P0, &
880  "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
881  "this is the freezing potential temperature at "//&
882  "S=0, P=0.", units="deg C", default=0.0)
883  call get_param(param_file, mdl, "DTFREEZE_DS",eos%dTFr_dS, &
884  "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
885  "this is the derivative of the freezing potential "//&
886  "temperature with salinity.", &
887  units="deg C PSU-1", default=-0.054)
888  call get_param(param_file, mdl, "DTFREEZE_DP",eos%dTFr_dP, &
889  "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
890  "this is the derivative of the freezing potential "//&
891  "temperature with pressure.", &
892  units="deg C Pa-1", default=0.0)
893  endif
894 
895  if ((eos%form_of_EOS == eos_teos10 .OR. eos%form_of_EOS == eos_nemo) .AND. &
896  eos%form_of_TFreeze /= tfreeze_teos10) then
897  call mom_error(fatal, "interpret_eos_selection: EOS_TEOS10 or EOS_NEMO \n" //&
898  "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 .")
899  endif
900 
901 
902 end subroutine eos_init
903 
904 !> Manually initialized an EOS type (intended for unit testing of routines which need a specific EOS)
905 subroutine eos_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, &
906  Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
907  type(eos_type), pointer :: eos !< Equation of state structure
908  integer, optional, intent(in) :: form_of_eos !< A coded integer indicating the equation of state to use.
909  integer, optional, intent(in) :: form_of_tfreeze !< A coded integer indicating the expression for
910  !! the potential temperature of the freezing point.
911  logical, optional, intent(in) :: eos_quadrature !< If true, always use the generic (quadrature)
912  !! code for the integrals of density.
913  logical, optional, intent(in) :: compressible !< If true, in situ density is a function of pressure.
914  real , optional, intent(in) :: rho_t0_s0 !< Density at T=0 degC and S=0 ppt [kg m-3]
915  real , optional, intent(in) :: drho_dt !< Partial derivative of density with temperature
916  !! in [kg m-3 degC-1]
917  real , optional, intent(in) :: drho_ds !< Partial derivative of density with salinity
918  !! in [kg m-3 ppt-1]
919  real , optional, intent(in) :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC].
920  real , optional, intent(in) :: dtfr_ds !< The derivative of freezing point with salinity
921  !! in [degC ppt-1].
922  real , optional, intent(in) :: dtfr_dp !< The derivative of freezing point with pressure
923  !! in [degC Pa-1].
924 
925  if (present(form_of_eos )) eos%form_of_EOS = form_of_eos
926  if (present(form_of_tfreeze)) eos%form_of_TFreeze = form_of_tfreeze
927  if (present(eos_quadrature )) eos%EOS_quadrature = eos_quadrature
928  if (present(compressible )) eos%Compressible = compressible
929  if (present(rho_t0_s0 )) eos%Rho_T0_S0 = rho_t0_s0
930  if (present(drho_dt )) eos%drho_dT = drho_dt
931  if (present(drho_ds )) eos%dRho_dS = drho_ds
932  if (present(tfr_s0_p0 )) eos%TFr_S0_P0 = tfr_s0_p0
933  if (present(dtfr_ds )) eos%dTFr_dS = dtfr_ds
934  if (present(dtfr_dp )) eos%dTFr_dp = dtfr_dp
935 
936 end subroutine eos_manual_init
937 
938 !> Allocates EOS_type
939 subroutine eos_allocate(EOS)
940  type(eos_type), pointer :: eos !< Equation of state structure
941 
942  if (.not.associated(eos)) allocate(eos)
943 end subroutine eos_allocate
944 
945 !> Deallocates EOS_type
946 subroutine eos_end(EOS)
947  type(eos_type), pointer :: eos !< Equation of state structure
948 
949  if (associated(eos)) deallocate(eos)
950 end subroutine eos_end
951 
952 !> Set equation of state structure (EOS) to linear with given coefficients
953 !!
954 !! \note This routine is primarily for testing and allows a local copy of the
955 !! EOS_type (EOS argument) to be set to use the linear equation of state
956 !! independent from the rest of the model.
957 subroutine eos_use_linear(Rho_T0_S0, dRho_dT, dRho_dS, EOS, use_quadrature)
958  real, intent(in) :: rho_t0_s0 !< Density at T=0 degC and S=0 ppt [kg m-3]
959  real, intent(in) :: drho_dt !< Partial derivative of density with temperature [kg m-3 degC-1]
960  real, intent(in) :: drho_ds !< Partial derivative of density with salinity [kg m-3 ppt-1]
961  logical, optional, intent(in) :: use_quadrature !< If true, always use the generic (quadrature)
962  !! code for the integrals of density.
963  type(eos_type), pointer :: eos !< Equation of state structure
964 
965  if (.not.associated(eos)) call mom_error(fatal, &
966  "MOM_EOS.F90: EOS_use_linear() called with an unassociated EOS_type EOS.")
967 
968  eos%form_of_EOS = eos_linear
969  eos%Compressible = .false.
970  eos%Rho_T0_S0 = rho_t0_s0
971  eos%dRho_dT = drho_dt
972  eos%dRho_dS = drho_ds
973  eos%EOS_quadrature = .false.
974  if (present(use_quadrature)) eos%EOS_quadrature = use_quadrature
975 
976 end subroutine eos_use_linear
977 
978 !> This subroutine calculates (by numerical quadrature) integrals of
979 !! pressure anomalies across layers, which are required for calculating the
980 !! finite-volume form pressure accelerations in a Boussinesq model.
981 subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
982  EOS, dpa, intz_dpa, intx_dpa, inty_dpa, &
983  bathyT, dz_neglect, useMassWghtInterp)
984  type(hor_index_type), intent(in) :: hii !< Horizontal index type for input variables.
985  type(hor_index_type), intent(in) :: hio !< Horizontal index type for output variables.
986  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
987  intent(in) :: t !< Potential temperature of the layer [degC].
988  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
989  intent(in) :: s !< Salinity of the layer [ppt].
990  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
991  intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m].
992  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
993  intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m].
994  real, intent(in) :: rho_ref !< A mean density [kg m-3], that is
995  !! subtracted out to reduce the magnitude
996  !! of each of the integrals.
997  real, intent(in) :: rho_0 !< A density [kg m-3], that is used
998  !! to calculate the pressure (as p~=-z*rho_0*G_e)
999  !! used in the equation of state.
1000  real, intent(in) :: g_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
1001  type(eos_type), pointer :: eos !< Equation of state structure
1002  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1003  intent(out) :: dpa !< The change in the pressure anomaly
1004  !! across the layer [Pa].
1005  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1006  optional, intent(out) :: intz_dpa !< The integral through the thickness of the
1007  !! layer of the pressure anomaly relative to the
1008  !! anomaly at the top of the layer [Pa Z ~> Pa m].
1009  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
1010  optional, intent(out) :: intx_dpa !< The integral in x of the difference between
1011  !! the pressure anomaly at the top and bottom of the
1012  !! layer divided by the x grid spacing [Pa].
1013  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
1014  optional, intent(out) :: inty_dpa !< The integral in y of the difference between
1015  !! the pressure anomaly at the top and bottom of the
1016  !! layer divided by the y grid spacing [Pa].
1017  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1018  optional, intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m].
1019  real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m].
1020  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
1021  !! interpolate T/S for top and bottom integrals.
1022  real :: t5(5), s5(5), p5(5), r5(5)
1023  real :: rho_anom ! The depth averaged density anomaly [kg m-3].
1024  real :: w_left, w_right
1025  real, parameter :: c1_90 = 1.0/90.0 ! Rational constants.
1026  real :: gxrho, i_rho
1027  real :: dz ! The layer thickness [Z ~> m].
1028  real :: hwght ! A pressure-thickness below topography [Z ~> m].
1029  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Z ~> m].
1030  real :: idenom ! The inverse of the denominator in the weights [Z-2 ~> m-2].
1031  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
1032  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
1033  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
1034  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
1035  real :: intz(5) ! The gravitational acceleration times the integrals of density
1036  ! with height at the 5 sub-column locations [Pa].
1037  logical :: do_massweight ! Indicates whether to do mass weighting.
1038  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n, ioff, joff
1039 
1040  ioff = hio%idg_offset - hii%idg_offset
1041  joff = hio%jdg_offset - hii%jdg_offset
1042 
1043  ! These array bounds work for the indexing convention of the input arrays, but
1044  ! on the computational domain defined for the output arrays.
1045  isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
1046  jsq = hio%JscB + joff ; jeq = hio%JecB + joff
1047  is = hio%isc + ioff ; ie = hio%iec + ioff
1048  js = hio%jsc + joff ; je = hio%jec + joff
1049 
1050  gxrho = g_e * rho_0
1051  i_rho = 1.0 / rho_0
1052 
1053  do_massweight = .false.
1054  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
1055  do_massweight = .true.
1056  if (.not.present(bathyt)) call mom_error(fatal, "int_density_dz_generic: "//&
1057  "bathyT must be present if useMassWghtInterp is present and true.")
1058  if (.not.present(dz_neglect)) call mom_error(fatal, "int_density_dz_generic: "//&
1059  "dz_neglect must be present if useMassWghtInterp is present and true.")
1060  endif ; endif
1061 
1062  do j=jsq,jeq+1 ; do i=isq,ieq+1
1063  dz = z_t(i,j) - z_b(i,j)
1064  do n=1,5
1065  t5(n) = t(i,j) ; s5(n) = s(i,j)
1066  p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
1067  enddo
1068  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
1069 
1070  ! Use Bode's rule to estimate the pressure anomaly change.
1071  rho_anom = c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))
1072  dpa(i-ioff,j-joff) = g_e*dz*rho_anom
1073  ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of
1074  ! the pressure anomaly.
1075  if (present(intz_dpa)) intz_dpa(i-ioff,j-joff) = 0.5*g_e*dz**2 * &
1076  (rho_anom - c1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
1077  enddo ; enddo
1078 
1079  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
1080  ! hWght is the distance measure by which the cell is violation of
1081  ! hydrostatic consistency. For large hWght we bias the interpolation of
1082  ! T & S along the top and bottom integrals, akin to thickness weighting.
1083  hwght = 0.0
1084  if (do_massweight) &
1085  hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
1086  if (hwght > 0.) then
1087  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
1088  hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
1089  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1090  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1091  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1092  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1093  else
1094  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1095  endif
1096 
1097  intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
1098  do m=2,4
1099  ! T, S, and z are interpolated in the horizontal. The z interpolation
1100  ! is linear, but for T and S it may be thickness weighted.
1101  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1102  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1103  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
1104  t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
1105  s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
1106  p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i+1,j))
1107  do n=2,5
1108  t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) + gxrho*0.25*dz
1109  enddo
1110  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
1111 
1112  ! Use Bode's rule to estimate the pressure anomaly change.
1113  intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
1114  enddo
1115  ! Use Bode's rule to integrate the bottom pressure anomaly values in x.
1116  intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1117  12.0*intz(3))
1118  enddo ; enddo ; endif
1119 
1120  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
1121  ! hWght is the distance measure by which the cell is violation of
1122  ! hydrostatic consistency. For large hWght we bias the interpolation of
1123  ! T & S along the top and bottom integrals, akin to thickness weighting.
1124  hwght = 0.0
1125  if (do_massweight) &
1126  hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
1127  if (hwght > 0.) then
1128  hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
1129  hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
1130  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1131  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
1132  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
1133  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
1134  else
1135  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
1136  endif
1137 
1138  intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i-ioff,j-joff+1)
1139  do m=2,4
1140  ! T, S, and z are interpolated in the horizontal. The z interpolation
1141  ! is linear, but for T and S it may be thickness weighted.
1142  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
1143  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
1144  dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
1145  t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
1146  s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
1147  p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i,j+1))
1148  do n=2,5
1149  t5(n) = t5(1) ; s5(n) = s5(1)
1150  p5(n) = p5(n-1) + gxrho*0.25*dz
1151  enddo
1152  call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
1153 
1154  ! Use Bode's rule to estimate the pressure anomaly change.
1155  intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
1156  enddo
1157  ! Use Bode's rule to integrate the values.
1158  inty_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1159  12.0*intz(3))
1160  enddo ; enddo ; endif
1161 end subroutine int_density_dz_generic
1162 
1163 
1164 ! ==========================================================================
1165 !> Compute pressure gradient force integrals by quadrature for the case where
1166 !! T and S are linear profiles.
1167 subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, &
1168  rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, &
1169  intz_dpa, intx_dpa, inty_dpa, &
1170  useMassWghtInterp)
1171  type(hor_index_type), intent(in) :: hii !< Ocean horizontal index structures for the input arrays
1172  type(hor_index_type), intent(in) :: hio !< Ocean horizontal index structures for the output arrays
1173  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1174  intent(in) :: t_t !< Potential temperatue at the cell top [degC]
1175  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1176  intent(in) :: t_b !< Potential temperatue at the cell bottom [degC]
1177  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1178  intent(in) :: s_t !< Salinity at the cell top [ppt]
1179  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1180  intent(in) :: s_b !< Salinity at the cell bottom [ppt]
1181  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1182  intent(in) :: z_t !< The geometric height at the top of the layer,
1183  !! in depth units [Z ~> m].
1184  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1185  intent(in) :: z_b !< The geometric height at the bottom of the layer [Z ~> m].
1186  real, intent(in) :: rho_ref !< A mean density [kg m-3], that is subtracted out to
1187  !! reduce the magnitude of each of the integrals.
1188  real, intent(in) :: rho_0 !< A density [kg m-3], that is used to calculate the
1189  !! pressure (as p~=-z*rho_0*G_e) used in the equation of state.
1190  real, intent(in) :: g_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
1191  real, intent(in) :: dz_subroundoff !< A miniscule thickness change [Z ~> m].
1192  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1193  intent(in) :: bathyt !< The depth of the bathymetry [Z ~> m].
1194  type(eos_type), pointer :: eos !< Equation of state structure
1195  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1196  intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
1197  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1198  optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
1199  !! the pressure anomaly relative to the anomaly at the
1200  !! top of the layer [Pa Z].
1201  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
1202  optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
1203  !! pressure anomaly at the top and bottom of the layer
1204  !! divided by the x grid spacing [Pa].
1205  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
1206  optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
1207  !! pressure anomaly at the top and bottom of the layer
1208  !! divided by the y grid spacing [Pa].
1209  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting to
1210  !! interpolate T/S for top and bottom integrals.
1211 ! This subroutine calculates (by numerical quadrature) integrals of
1212 ! pressure anomalies across layers, which are required for calculating the
1213 ! finite-volume form pressure accelerations in a Boussinesq model. The one
1214 ! potentially dodgy assumtion here is that rho_0 is used both in the denominator
1215 ! of the accelerations, and in the pressure used to calculated density (the
1216 ! latter being -z*rho_0*G_e). These two uses could be separated if need be.
1217 !
1218 ! It is assumed that the salinity and temperature profiles are linear in the
1219 ! vertical. The top and bottom values within each layer are provided and
1220 ! a linear interpolation is used to compute intermediate values.
1221 
1222  ! Local variables
1223  real :: t5((5*hio%iscb+1):(5*(hio%iecb+2))) ! Temperatures along a line of subgrid locations [degC].
1224  real :: s5((5*hio%iscb+1):(5*(hio%iecb+2))) ! Salinities along a line of subgrid locations [ppt].
1225  real :: p5((5*hio%iscb+1):(5*(hio%iecb+2))) ! Pressures along a line of subgrid locations [Pa].
1226  real :: r5((5*hio%iscb+1):(5*(hio%iecb+2))) ! Densities along a line of subgrid locations [kg m-3].
1227  real :: t15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Temperatures at an array of subgrid locations [degC].
1228  real :: s15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Salinities at an array of subgrid locations [ppt].
1229  real :: p15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Pressures at an array of subgrid locations [Pa].
1230  real :: r15((15*hio%iscb+1):(15*(hio%iecb+1))) ! Densities at an array of subgrid locations [kg m-3].
1231  real :: wt_t(5), wt_b(5) ! Top and bottom weights [nondim].
1232  real :: rho_anom ! A density anomaly [kg m-3].
1233  real :: w_left, w_right ! Left and right weights [nondim].
1234  real :: intz(5) ! The gravitational acceleration times the integrals of density
1235  ! with height at the 5 sub-column locations [Pa].
1236  real, parameter :: c1_90 = 1.0/90.0 ! A rational constant [nondim].
1237  real :: gxrho ! Gravitational acceleration times density [kg m-1 Z-1 s-2 ~> kg m-2 s-2].
1238  real :: i_rho ! The inverse of the reference density [m3 kg-1].
1239  real :: dz(hio%iscb:hio%iecb+1) ! Layer thicknesses at tracer points [Z ~> m].
1240  real :: dz_x(5,hio%iscb:hio%iecb) ! Layer thicknesses along an x-line of subrid locations [Z ~> m].
1241  real :: dz_y(5,hio%isc:hio%iec) ! Layer thicknesses along a y-line of subrid locations [Z ~> m].
1242  real :: weight_t, weight_b ! Nondimensional weights of the top and bottom.
1243  real :: massweighttoggle ! A nondimensional toggle factor (0 or 1).
1244  real :: ttl, tbl, ttr, tbr ! Temperatures at the velocity cell corners [degC].
1245  real :: stl, sbl, str, sbr ! Salinities at the velocity cell corners [ppt].
1246  real :: hwght ! A topographically limited thicknes weight [Z ~> m].
1247  real :: hl, hr ! Thicknesses to the left and right [Z ~> m].
1248  real :: idenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2].
1249  integer :: isq, ieq, jsq, jeq, i, j, m, n
1250  integer :: iin, jin, ioff, joff
1251  integer :: pos
1252 
1253  ioff = hio%idg_offset - hii%idg_offset
1254  joff = hio%jdg_offset - hii%jdg_offset
1255 
1256  isq = hio%IscB ; ieq = hio%IecB ; jsq = hio%JscB ; jeq = hio%JecB
1257 
1258  gxrho = g_e * rho_0
1259  i_rho = 1.0 / rho_0
1260  massweighttoggle = 0.
1261  if (present(usemasswghtinterp)) then
1262  if (usemasswghtinterp) massweighttoggle = 1.
1263  endif
1264 
1265  do n = 1, 5
1266  wt_t(n) = 0.25 * real(5-n)
1267  wt_b(n) = 1.0 - wt_t(n)
1268  enddo
1269 
1270  ! =============================
1271  ! 1. Compute vertical integrals
1272  ! =============================
1273  do j=jsq,jeq+1
1274  jin = j+joff
1275  do i = isq,ieq+1 ; iin = i+ioff
1276  dz(i) = z_t(iin,jin) - z_b(iin,jin)
1277  do n=1,5
1278  p5(i*5+n) = -gxrho*(z_t(iin,jin) - 0.25*real(n-1)*dz(i))
1279  ! Salinity and temperature points are linearly interpolated
1280  s5(i*5+n) = wt_t(n) * s_t(iin,jin) + wt_b(n) * s_b(iin,jin)
1281  t5(i*5+n) = wt_t(n) * t_t(iin,jin) + wt_b(n) * t_b(iin,jin)
1282  enddo
1283  enddo
1284  call calculate_density_array(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho_ref )
1285 
1286  do i=isq,ieq+1 ; iin = i+ioff
1287  ! Use Bode's rule to estimate the pressure anomaly change.
1288  rho_anom = c1_90*(7.0*(r5(i*5+1)+r5(i*5+5)) + 32.0*(r5(i*5+2)+r5(i*5+4)) + 12.0*r5(i*5+3))
1289  dpa(i,j) = g_e*dz(i)*rho_anom
1290  if (present(intz_dpa)) then
1291  ! Use a Bode's-rule-like fifth-order accurate estimate of
1292  ! the double integral of the pressure anomaly.
1293  intz_dpa(i,j) = 0.5*g_e*dz(i)**2 * &
1294  (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )
1295  endif
1296  enddo
1297  enddo ! end loops on j
1298 
1299 
1300  ! ==================================================
1301  ! 2. Compute horizontal integrals in the x direction
1302  ! ==================================================
1303  if (present(intx_dpa)) then ; do j=hio%jsc,hio%jec ; jin = j+joff
1304  do i=isq,ieq ; iin = i+ioff
1305  ! Corner values of T and S
1306  ! hWght is the distance measure by which the cell is violation of
1307  ! hydrostatic consistency. For large hWght we bias the interpolation
1308  ! of T,S along the top and bottom integrals, almost like thickness
1309  ! weighting.
1310  ! Note: To work in terrain following coordinates we could offset
1311  ! this distance by the layer thickness to replicate other models.
1312  hwght = massweighttoggle * &
1313  max(0., -bathyt(iin,jin)-z_t(iin+1,jin), -bathyt(iin+1,jin)-z_t(iin,jin))
1314  if (hwght > 0.) then
1315  hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
1316  hr = (z_t(iin+1,jin) - z_b(iin+1,jin)) + dz_subroundoff
1317  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1318  idenom = 1./( hwght*(hr + hl) + hl*hr )
1319  ttl = ( (hwght*hr)*t_t(iin+1,jin) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
1320  ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin+1,jin) ) * idenom
1321  tbl = ( (hwght*hr)*t_b(iin+1,jin) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
1322  tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin+1,jin) ) * idenom
1323  stl = ( (hwght*hr)*s_t(iin+1,jin) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
1324  str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin+1,jin) ) * idenom
1325  sbl = ( (hwght*hr)*s_b(iin+1,jin) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
1326  sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin+1,jin) ) * idenom
1327  else
1328  ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin+1,jin); tbr = t_b(iin+1,jin)
1329  stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin+1,jin); sbr = s_b(iin+1,jin)
1330  endif
1331 
1332  do m=2,4
1333  w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
1334  dz_x(m,i) = w_left*(z_t(iin,jin) - z_b(iin,jin)) + w_right*(z_t(iin+1,jin) - z_b(iin+1,jin))
1335 
1336  ! Salinity and temperature points are linearly interpolated in
1337  ! the horizontal. The subscript (1) refers to the top value in
1338  ! the vertical profile while subscript (5) refers to the bottom
1339  ! value in the vertical profile.
1340  pos = i*15+(m-2)*5
1341  t15(pos+1) = w_left*ttl + w_right*ttr
1342  t15(pos+5) = w_left*tbl + w_right*tbr
1343 
1344  s15(pos+1) = w_left*stl + w_right*str
1345  s15(pos+5) = w_left*sbl + w_right*sbr
1346 
1347  p15(pos+1) = -gxrho*(w_left*z_t(iin,jin) + w_right*z_t(iin+1,jin))
1348 
1349  ! Pressure
1350  do n=2,5
1351  p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)
1352  enddo
1353 
1354  ! Salinity and temperature (linear interpolation in the vertical)
1355  do n=2,4
1356  weight_t = 0.25 * real(5-n)
1357  weight_b = 1.0 - weight_t
1358  s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
1359  t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
1360  enddo
1361  enddo
1362  enddo
1363 
1364  call calculate_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho_ref)
1365 
1366  do i=isq,ieq ; iin = i+ioff
1367  intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
1368 
1369  ! Use Bode's rule to estimate the pressure anomaly change.
1370  do m = 2,4
1371  pos = i*15+(m-2)*5
1372  intz(m) = g_e*dz_x(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + 32.0*(r15(pos+2)+r15(pos+4)) + &
1373  12.0*r15(pos+3)))
1374  enddo
1375  ! Use Bode's rule to integrate the bottom pressure anomaly values in x.
1376  intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1377  12.0*intz(3))
1378  enddo
1379  enddo ; endif
1380 
1381  ! ==================================================
1382  ! 3. Compute horizontal integrals in the y direction
1383  ! ==================================================
1384  if (present(inty_dpa)) then ; do j=jsq,jeq ; jin = j+joff
1385  do i=hio%isc,hio%iec ; iin = i+ioff
1386  ! Corner values of T and S
1387  ! hWght is the distance measure by which the cell is violation of
1388  ! hydrostatic consistency. For large hWght we bias the interpolation
1389  ! of T,S along the top and bottom integrals, almost like thickness
1390  ! weighting.
1391  ! Note: To work in terrain following coordinates we could offset
1392  ! this distance by the layer thickness to replicate other models.
1393  hwght = massweighttoggle * &
1394  max(0., -bathyt(i,j)-z_t(iin,jin+1), -bathyt(i,j+1)-z_t(iin,jin))
1395  if (hwght > 0.) then
1396  hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
1397  hr = (z_t(iin,jin+1) - z_b(iin,jin+1)) + dz_subroundoff
1398  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
1399  idenom = 1./( hwght*(hr + hl) + hl*hr )
1400  ttl = ( (hwght*hr)*t_t(iin,jin+1) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
1401  ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin,jin+1) ) * idenom
1402  tbl = ( (hwght*hr)*t_b(iin,jin+1) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
1403  tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin,jin+1) ) * idenom
1404  stl = ( (hwght*hr)*s_t(iin,jin+1) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
1405  str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin,jin+1) ) * idenom
1406  sbl = ( (hwght*hr)*s_b(iin,jin+1) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
1407  sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin,jin+1) ) * idenom
1408  else
1409  ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin,jin+1); tbr = t_b(iin,jin+1)
1410  stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin,jin+1); sbr = s_b(iin,jin+1)
1411  endif
1412 
1413  do m=2,4
1414  w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
1415  dz_y(m,i) = w_left*(z_t(iin,jin) - z_b(iin,jin)) + w_right*(z_t(iin,jin+1) - z_b(iin,jin+1))
1416 
1417  ! Salinity and temperature points are linearly interpolated in
1418  ! the horizontal. The subscript (1) refers to the top value in
1419  ! the vertical profile while subscript (5) refers to the bottom
1420  ! value in the vertical profile.
1421  pos = i*15+(m-2)*5
1422  t15(pos+1) = w_left*ttl + w_right*ttr
1423  t15(pos+5) = w_left*tbl + w_right*tbr
1424 
1425  s15(pos+1) = w_left*stl + w_right*str
1426  s15(pos+5) = w_left*sbl + w_right*sbr
1427 
1428  p15(pos+1) = -gxrho*(w_left*z_t(iin,jin) + w_right*z_t(iin,jin+1))
1429 
1430  ! Pressure
1431  do n=2,5 ; p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i) ; enddo
1432 
1433  ! Salinity and temperature (linear interpolation in the vertical)
1434  do n=2,4
1435  weight_t = 0.25 * real(5-n)
1436  weight_b = 1.0 - weight_t
1437  s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
1438  t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
1439  enddo
1440  enddo
1441  enddo
1442 
1443  call calculate_density_array(t15(15*hio%isc+1:), s15(15*hio%isc+1:), p15(15*hio%isc+1:), &
1444  r15(15*hio%isc+1:), 1, 15*(hio%iec-hio%isc+1), eos, rho_ref)
1445  do i=hio%isc,hio%iec ; iin = i+ioff
1446  intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
1447 
1448  ! Use Bode's rule to estimate the pressure anomaly change.
1449  do m = 2,4
1450  pos = i*15+(m-2)*5
1451  intz(m) = g_e*dz_y(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &
1452  32.0*(r15(pos+2)+r15(pos+4)) + &
1453  12.0*r15(pos+3)))
1454  enddo
1455  ! Use Bode's rule to integrate the values.
1456  inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1457  12.0*intz(3))
1458  enddo
1459  enddo ; endif
1460 
1461 end subroutine int_density_dz_generic_plm
1462 ! ==========================================================================
1463 ! Above is the routine where only the S and T profiles are modified
1464 ! The real topography is still used
1465 ! ==========================================================================
1466 
1467 !> Find the depth at which the reconstructed pressure matches P_tgt
1468 subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, &
1469  rho_ref, G_e, EOS, P_b, z_out, z_tol)
1470  real, intent(in) :: t_t !< Potential temperatue at the cell top [degC]
1471  real, intent(in) :: t_b !< Potential temperatue at the cell bottom [degC]
1472  real, intent(in) :: s_t !< Salinity at the cell top [ppt]
1473  real, intent(in) :: s_b !< Salinity at the cell bottom [ppt]
1474  real, intent(in) :: z_t !< Absolute height of top of cell [Z ~> m]. (Boussinesq ????)
1475  real, intent(in) :: z_b !< Absolute height of bottom of cell [Z ~> m].
1476  real, intent(in) :: p_t !< Anomalous pressure of top of cell, relative to g*rho_ref*z_t [Pa]
1477  real, intent(in) :: p_tgt !< Target pressure at height z_out, relative to g*rho_ref*z_out [Pa]
1478  real, intent(in) :: rho_ref !< Reference density with which calculation are anomalous to
1479  real, intent(in) :: g_e !< Gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
1480  type(eos_type), pointer :: eos !< Equation of state structure
1481  real, intent(out) :: p_b !< Pressure at the bottom of the cell [Pa]
1482  real, intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m].
1483  real, optional, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m].
1484  ! Local variables
1485  real :: top_weight, bottom_weight, rho_anom, w_left, w_right, gxrho, dz, dp, f_guess, f_l, f_r
1486  real :: pa, pa_left, pa_right, pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz
1487 
1488  gxrho = g_e * rho_ref
1489 
1490  ! Anomalous pressure difference across whole cell
1491  dp = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, 1.0, eos)
1492 
1493  p_b = p_t + dp ! Anomalous pressure at bottom of cell
1494 
1495  if (p_tgt <= p_t ) then
1496  z_out = z_t
1497  return
1498  endif
1499 
1500  if (p_tgt >= p_b) then
1501  z_out = z_b
1502  return
1503  endif
1504 
1505  f_l = 0.
1506  pa_left = p_t - p_tgt ! Pa_left < 0
1507  f_r = 1.
1508  pa_right = p_b - p_tgt ! Pa_right > 0
1509  pa_tol = gxrho * 1.e-5 ! 1e-5 has dimensions of m, but should be converted to the units of z.
1510  if (present(z_tol)) pa_tol = gxrho * z_tol
1511  f_guess = f_l - pa_left / ( pa_right -pa_left ) * ( f_r - f_l )
1512  pa = pa_right - pa_left ! To get into iterative loop
1513  do while ( abs(pa) > pa_tol )
1514 
1515  z_out = z_t + ( z_b - z_t ) * f_guess
1516  pa = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, f_guess, eos) - ( p_tgt - p_t )
1517 
1518  if (pa<pa_left) then
1519  write(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
1520  stop 'Blurgh! Too negative'
1521  elseif (pa<0.) then
1522  pa_left = pa
1523  f_l = f_guess
1524  elseif (pa>pa_right) then
1525  write(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
1526  stop 'Blurgh! Too positive'
1527  elseif (pa>0.) then
1528  pa_right = pa
1529  f_r = f_guess
1530  else ! Pa == 0
1531  return
1532  endif
1533  f_guess = f_l - pa_left / ( pa_right -pa_left ) * ( f_r - f_l )
1534 
1535  enddo
1536 
1537 end subroutine find_depth_of_pressure_in_cell
1538 
1539 !> Returns change in anomalous pressure change from top to non-dimensional
1540 !! position pos between z_t and z_b
1541 real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS)
1542  real, intent(in) :: t_t !< Potential temperatue at the cell top [degC]
1543  real, intent(in) :: t_b !< Potential temperatue at the cell bottom [degC]
1544  real, intent(in) :: s_t !< Salinity at the cell top [ppt]
1545  real, intent(in) :: s_b !< Salinity at the cell bottom [ppt]
1546  real, intent(in) :: z_t !< The geometric height at the top of the layer [Z ~> m]
1547  real, intent(in) :: z_b !< The geometric height at the bottom of the layer [Z ~> m]
1548  real, intent(in) :: rho_ref !< A mean density [kg m-3], that is subtracted out to
1549  !! reduce the magnitude of each of the integrals.
1550  real, intent(in) :: g_e !< The Earth's gravitational acceleration [m s-2]
1551  real, intent(in) :: pos !< The fractional vertical position, 0 to 1 [nondim].
1552  type(eos_type), pointer :: eos !< Equation of state structure
1553  ! Local variables
1554  real, parameter :: c1_90 = 1.0/90.0 ! Rational constants.
1555  real :: dz, top_weight, bottom_weight, rho_ave
1556  real, dimension(5) :: t5, s5, p5, rho5
1557  integer :: n
1558 
1559  do n=1,5
1560  ! Evalute density at five quadrature points
1561  bottom_weight = 0.25*real(n-1) * pos
1562  top_weight = 1.0 - bottom_weight
1563  ! Salinity and temperature points are linearly interpolated
1564  s5(n) = top_weight * s_t + bottom_weight * s_b
1565  t5(n) = top_weight * t_t + bottom_weight * t_b
1566  p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( g_e * rho_ref )
1567  enddo
1568  call calculate_density_array(t5, s5, p5, rho5, 1, 5, eos)
1569  rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref
1570 
1571  ! Use Boole's rule to estimate the average density
1572  rho_ave = c1_90*(7.0*(rho5(1)+rho5(5)) + 32.0*(rho5(2)+rho5(4)) + 12.0*rho5(3))
1573 
1574  dz = ( z_t - z_b ) * pos
1575  frac_dp_at_pos = g_e * dz * rho_ave
1576 end function frac_dp_at_pos
1577 
1578 
1579 ! ==========================================================================
1580 !> Compute pressure gradient force integrals for the case where T and S
1581 !! are parabolic profiles
1582 subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, &
1583  z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
1584  EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
1586  type(hor_index_type), intent(in) :: hii !< Ocean horizontal index structures for the input arrays
1587  type(hor_index_type), intent(in) :: hio !< Ocean horizontal index structures for the output arrays
1588  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1589  intent(in) :: t !< Potential temperature referenced to the surface [degC]
1590  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1591  intent(in) :: t_t !< Potential temperatue at the cell top [degC]
1592  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1593  intent(in) :: t_b !< Potential temperatue at the cell bottom [degC]
1594  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1595  intent(in) :: s !< Salinity [ppt]
1596  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1597  intent(in) :: s_t !< Salinity at the cell top [ppt]
1598  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1599  intent(in) :: s_b !< Salinity at the cell bottom [ppt]
1600  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1601  intent(in) :: z_t !< Height at the top of the layer [Z ~> m].
1602  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
1603  intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m].
1604  real, intent(in) :: rho_ref !< A mean density [kg m-3], that is subtracted out to
1605  !! reduce the magnitude of each of the integrals.
1606  real, intent(in) :: rho_0 !< A density [kg m-3], that is used to calculate the
1607  !! pressure (as p~=-z*rho_0*G_e) used in the equation of state.
1608  real, intent(in) :: g_e !< The Earth's gravitational acceleration [m s-2]
1609  type(eos_type), pointer :: eos !< Equation of state structure
1610  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1611  intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
1612  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
1613  optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
1614  !! the pressure anomaly relative to the anomaly at the
1615  !! top of the layer [Pa Z ~> Pa m].
1616  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
1617  optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
1618  !! pressure anomaly at the top and bottom of the layer
1619  !! divided by the x grid spacing [Pa].
1620  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
1621  optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
1622  !! pressure anomaly at the top and bottom of the layer
1623  !! divided by the y grid spacing [Pa].
1624 
1625 ! This subroutine calculates (by numerical quadrature) integrals of
1626 ! pressure anomalies across layers, which are required for calculating the
1627 ! finite-volume form pressure accelerations in a Boussinesq model. The one
1628 ! potentially dodgy assumtion here is that rho_0 is used both in the denominator
1629 ! of the accelerations, and in the pressure used to calculated density (the
1630 ! latter being -z*rho_0*G_e). These two uses could be separated if need be.
1631 !
1632 ! It is assumed that the salinity and temperature profiles are linear in the
1633 ! vertical. The top and bottom values within each layer are provided and
1634 ! a linear interpolation is used to compute intermediate values.
1635 
1636  ! Local variables
1637  real :: t5(5), s5(5), p5(5), r5(5)
1638  real :: rho_anom
1639  real :: w_left, w_right, intz(5)
1640  real, parameter :: c1_90 = 1.0/90.0 ! Rational constants.
1641  real :: gxrho, i_rho
1642  real :: dz
1643  real :: weight_t, weight_b
1644  real :: s0, s1, s2 ! parabola coefficients for S [ppt]
1645  real :: t0, t1, t2 ! parabola coefficients for T [degC]
1646  real :: xi ! normalized coordinate
1647  real :: t_top, t_mid, t_bot
1648  real :: s_top, s_mid, s_bot
1649  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j, m, n, ioff, joff
1650  real, dimension(4) :: x, y
1651  real, dimension(9) :: s_node, t_node, p_node, r_node
1652 
1653 
1654  call mom_error(fatal, &
1655  "int_density_dz_generic_ppm: the implementation is not done yet, contact developer")
1656 
1657  ioff = hio%idg_offset - hii%idg_offset
1658  joff = hio%jdg_offset - hii%jdg_offset
1659 
1660  ! These array bounds work for the indexing convention of the input arrays, but
1661  ! on the computational domain defined for the output arrays.
1662  isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
1663  jsq = hio%JscB + joff ; jeq = hio%JecB + joff
1664  is = hio%isc + ioff ; ie = hio%iec + ioff
1665  js = hio%jsc + joff ; je = hio%jec + joff
1666 
1667  gxrho = g_e * rho_0
1668  i_rho = 1.0 / rho_0
1669 
1670  ! =============================
1671  ! 1. Compute vertical integrals
1672  ! =============================
1673  do j=jsq,jeq+1 ; do i=isq,ieq+1
1674  dz = z_t(i,j) - z_b(i,j)
1675 
1676  ! Coefficients of the parabola for S
1677  s0 = s_t(i,j)
1678  s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
1679  s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0*s(i,j) )
1680 
1681  ! Coefficients of the parabola for T
1682  t0 = t_t(i,j)
1683  t1 = 6.0 * t(i,j) - 4.0 * t_t(i,j) - 2.0 * t_b(i,j)
1684  t2 = 3.0 * ( t_t(i,j) + t_b(i,j) - 2.0*t(i,j) )
1685 
1686  do n=1,5
1687  p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
1688 
1689  ! Parabolic reconstruction for T and S
1690  xi = 0.25 * ( n - 1 )
1691  s5(n) = s0 + s1 * xi + s2 * xi**2
1692  t5(n) = t0 + t1 * xi + t2 * xi**2
1693  enddo
1694 
1695  call calculate_density(t5, s5, p5, r5, 1, 5, eos)
1696 
1697  ! Use Bode's rule to estimate the pressure anomaly change.
1698  !rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) - &
1699  ! rho_ref
1700 
1701  rho_anom = 1000.0 + s(i,j) - rho_ref
1702  dpa(i-ioff,j-joff) = g_e*dz*rho_anom
1703 
1704  ! Use a Bode's-rule-like fifth-order accurate estimate of
1705  ! the double integral of the pressure anomaly.
1706  !r5 = r5 - rho_ref
1707  !if (present(intz_dpa)) intz_dpa(i,j) = 0.5*G_e*dz**2 * &
1708  ! (rho_anom - C1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
1709 
1710  intz_dpa(i-ioff,j-joff) = 0.5 * g_e * dz**2 * ( 1000.0 - rho_ref + s0 + s1/3.0 + &
1711  s2/6.0 )
1712  enddo ; enddo ! end loops on j and i
1713 
1714  ! ==================================================
1715  ! 2. Compute horizontal integrals in the x direction
1716  ! ==================================================
1717  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
1718  intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
1719  do m=2,4
1720  w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
1721  dz = w_left*(z_t(i,j) - z_b(i,j)) + w_right*(z_t(i+1,j) - z_b(i+1,j))
1722 
1723  ! Salinity and temperature points are linearly interpolated in
1724  ! the horizontal. The subscript (1) refers to the top value in
1725  ! the vertical profile while subscript (5) refers to the bottom
1726  ! value in the vertical profile.
1727  t_top = w_left*t_t(i,j) + w_right*t_t(i+1,j)
1728  t_mid = w_left*t(i,j) + w_right*t(i+1,j)
1729  t_bot = w_left*t_b(i,j) + w_right*t_b(i+1,j)
1730 
1731  s_top = w_left*s_t(i,j) + w_right*s_t(i+1,j)
1732  s_mid = w_left*s(i,j) + w_right*s(i+1,j)
1733  s_bot = w_left*s_b(i,j) + w_right*s_b(i+1,j)
1734 
1735  p5(1) = -gxrho*(w_left*z_t(i,j) + w_right*z_t(i+1,j))
1736 
1737  ! Pressure
1738  do n=2,5
1739  p5(n) = p5(n-1) + gxrho*0.25*dz
1740  enddo
1741 
1742  ! Coefficients of the parabola for S
1743  s0 = s_top
1744  s1 = 6.0 * s_mid - 4.0 * s_top - 2.0 * s_bot
1745  s2 = 3.0 * ( s_top + s_bot - 2.0*s_mid )
1746 
1747  ! Coefficients of the parabola for T
1748  t0 = t_top
1749  t1 = 6.0 * t_mid - 4.0 * t_top - 2.0 * t_bot
1750  t2 = 3.0 * ( t_top + t_bot - 2.0*t_mid )
1751 
1752  do n=1,5
1753  ! Parabolic reconstruction for T and S
1754  xi = 0.25 * ( n - 1 )
1755  s5(n) = s0 + s1 * xi + s2 * xi**2
1756  t5(n) = t0 + t1 * xi + t2 * xi**2
1757  enddo
1758 
1759  call calculate_density(t5, s5, p5, r5, 1, 5, eos)
1760 
1761  ! Use Bode's rule to estimate the pressure anomaly change.
1762  intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + &
1763  12.0*r5(3)) - rho_ref)
1764  enddo
1765  intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
1766  12.0*intz(3))
1767 
1768  ! Use Gauss quadrature rule to compute integral
1769 
1770  ! The following coordinates define the quadrilateral on which the integral
1771  ! is computed
1772  x(1) = 1.0
1773  x(2) = 0.0
1774  x(3) = 0.0
1775  x(4) = 1.0
1776  y(1) = z_t(i+1,j)
1777  y(2) = z_t(i,j)
1778  y(3) = z_b(i,j)
1779  y(4) = z_b(i+1,j)
1780 
1781  t_node = 0.0
1782  p_node = 0.0
1783 
1784  ! Nodal values for S
1785 
1786  ! Parabolic reconstruction on the left
1787  s0 = s_t(i,j)
1788  s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
1789  s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0 * s(i,j) )
1790  s_node(2) = s0
1791  s_node(6) = s0 + 0.5 * s1 + 0.25 * s2
1792  s_node(3) = s0 + s1 + s2
1793 
1794  ! Parabolic reconstruction on the left
1795  s0 = s_t(i+1,j)
1796  s1 = 6.0 * s(i+1,j) - 4.0 * s_t(i+1,j) - 2.0 * s_b(i+1,j)
1797  s2 = 3.0 * ( s_t(i+1,j) + s_b(i+1,j) - 2.0 * s(i+1,j) )
1798  s_node(1) = s0
1799  s_node(8) = s0 + 0.5 * s1 + 0.25 * s2
1800  s_node(4) = s0 + s1 + s2
1801 
1802  s_node(5) = 0.5 * ( s_node(2) + s_node(1) )
1803  s_node(9) = 0.5 * ( s_node(6) + s_node(8) )
1804  s_node(7) = 0.5 * ( s_node(3) + s_node(4) )
1805 
1806  call calculate_density( t_node, s_node, p_node, r_node, 1, 9, eos )
1807  r_node = r_node - rho_ref
1808 
1809  call compute_integral_quadratic( x, y, r_node, intx_dpa(i-ioff,j-joff) )
1810 
1811  intx_dpa(i-ioff,j-joff) = intx_dpa(i-ioff,j-joff) * g_e
1812 
1813  enddo ; enddo ; endif
1814 
1815  ! ==================================================
1816  ! 3. Compute horizontal integrals in the y direction
1817  ! ==================================================
1818  if (present(inty_dpa)) then
1819  call mom_error(warning, "int_density_dz_generic_ppm still needs to be written for inty_dpa!")
1820  do j=jsq,jeq ; do i=is,ie
1821 
1822  inty_dpa(i-ioff,j-joff) = 0.0
1823 
1824  enddo ; enddo
1825  endif
1826 
1827 end subroutine int_density_dz_generic_ppm
1828 
1829 
1830 
1831 ! =============================================================================
1832 !> Compute the integral of the quadratic function
1833 subroutine compute_integral_quadratic( x, y, f, integral )
1834  real, dimension(4), intent(in) :: x !< The x-position of the corners
1835  real, dimension(4), intent(in) :: y !< The y-position of the corners
1836  real, dimension(9), intent(in) :: f !< The function at the quadrature points
1837  real, intent(out) :: integral !< The returned integral
1838 
1839  ! Local variables
1840  integer :: i, k
1841  real, dimension(9) :: weight, xi, eta ! integration points
1842  real :: f_k
1843  real :: dxdxi, dxdeta
1844  real :: dydxi, dydeta
1845  real, dimension(4) :: phiiso, dphiisodxi, dphiisodeta
1846  real, dimension(9) :: phi, dphidxi, dphideta
1847  real :: jacobian_k
1848  real :: t
1849 
1850  ! Quadrature rule (4 points)
1851  !weight(:) = 1.0
1852  !xi(1) = - sqrt(3.0) / 3.0
1853  !xi(2) = sqrt(3.0) / 3.0
1854  !xi(3) = sqrt(3.0) / 3.0
1855  !xi(4) = - sqrt(3.0) / 3.0
1856  !eta(1) = - sqrt(3.0) / 3.0
1857  !eta(2) = - sqrt(3.0) / 3.0
1858  !eta(3) = sqrt(3.0) / 3.0
1859  !eta(4) = sqrt(3.0) / 3.0
1860 
1861  ! Quadrature rule (9 points)
1862  t = sqrt(3.0/5.0)
1863  weight(1) = 25.0/81.0 ; xi(1) = -t ; eta(1) = t
1864  weight(2) = 40.0/81.0 ; xi(2) = .0 ; eta(2) = t
1865  weight(3) = 25.0/81.0 ; xi(3) = t ; eta(3) = t
1866  weight(4) = 40.0/81.0 ; xi(4) = -t ; eta(4) = .0
1867  weight(5) = 64.0/81.0 ; xi(5) = .0 ; eta(5) = .0
1868  weight(6) = 40.0/81.0 ; xi(6) = t ; eta(6) = .0
1869  weight(7) = 25.0/81.0 ; xi(7) = -t ; eta(7) = -t
1870  weight(8) = 40.0/81.0 ; xi(8) = .0 ; eta(8) = -t
1871  weight(9) = 25.0/81.0 ; xi(9) = t ; eta(9) = -t
1872 
1873  integral = 0.0
1874 
1875  ! Integration loop
1876  do k = 1,9
1877 
1878  ! Evaluate shape functions and gradients for isomorphism
1879  call evaluate_shape_bilinear( xi(k), eta(k), phiiso, &
1880  dphiisodxi, dphiisodeta )
1881 
1882  ! Determine gradient of global coordinate at integration point
1883  dxdxi = 0.0
1884  dxdeta = 0.0
1885  dydxi = 0.0
1886  dydeta = 0.0
1887 
1888  do i = 1,4
1889  dxdxi = dxdxi + x(i) * dphiisodxi(i)
1890  dxdeta = dxdeta + x(i) * dphiisodeta(i)
1891  dydxi = dydxi + y(i) * dphiisodxi(i)
1892  dydeta = dydeta + y(i) * dphiisodeta(i)
1893  enddo
1894 
1895  ! Evaluate Jacobian at integration point
1896  jacobian_k = dxdxi*dydeta - dydxi*dxdeta
1897 
1898  ! Evaluate shape functions for interpolation
1899  call evaluate_shape_quadratic( xi(k), eta(k), phi, dphidxi, dphideta )
1900 
1901  ! Evaluate function at integration point
1902  f_k = 0.0
1903  do i = 1,9
1904  f_k = f_k + f(i) * phi(i)
1905  enddo
1906 
1907  integral = integral + weight(k) * f_k * jacobian_k
1908 
1909  enddo ! end integration loop
1910 
1911 end subroutine compute_integral_quadratic
1912 
1913 
1914 ! =============================================================================
1915 !> Evaluation of the four bilinear shape fn and their gradients at (xi,eta)
1916 subroutine evaluate_shape_bilinear( xi, eta, phi, dphidxi, dphideta )
1917  real, intent(in) :: xi !< The x position to evaluate
1918  real, intent(in) :: eta !< The z position to evaluate
1919  real, dimension(4), intent(out) :: phi !< The weights of the four corners at this point
1920  real, dimension(4), intent(out) :: dphidxi !< The x-gradient of the weights of the four
1921  !! corners at this point
1922  real, dimension(4), intent(out) :: dphideta !< The z-gradient of the weights of the four
1923  !! corners at this point
1924 
1925  ! The shape functions within the parent element are defined as shown here:
1926  !
1927  ! (-1,1) 2 o------------o 1 (1,1)
1928  ! | |
1929  ! | |
1930  ! | |
1931  ! | |
1932  ! (-1,-1) 3 o------------o 4 (1,-1)
1933  !
1934 
1935  phi(1) = 0.25 * ( 1 + xi ) * ( 1 + eta )
1936  phi(2) = 0.25 * ( 1 - xi ) * ( 1 + eta )
1937  phi(3) = 0.25 * ( 1 - xi ) * ( 1 - eta )
1938  phi(4) = 0.25 * ( 1 + xi ) * ( 1 - eta )
1939 
1940  dphidxi(1) = 0.25 * ( 1 + eta )
1941  dphidxi(2) = - 0.25 * ( 1 + eta )
1942  dphidxi(3) = - 0.25 * ( 1 - eta )
1943  dphidxi(4) = 0.25 * ( 1 - eta )
1944 
1945  dphideta(1) = 0.25 * ( 1 + xi )
1946  dphideta(2) = 0.25 * ( 1 - xi )
1947  dphideta(3) = - 0.25 * ( 1 - xi )
1948  dphideta(4) = - 0.25 * ( 1 + xi )
1949 
1950 end subroutine evaluate_shape_bilinear
1951 
1952 
1953 ! =============================================================================
1954 !> Evaluation of the nine quadratic shape fn weights and their gradients at (xi,eta)
1955 subroutine evaluate_shape_quadratic ( xi, eta, phi, dphidxi, dphideta )
1957  ! Arguments
1958  real, intent(in) :: xi !< The x position to evaluate
1959  real, intent(in) :: eta !< The z position to evaluate
1960  real, dimension(9), intent(out) :: phi !< The weights of the 9 bilinear quadrature points
1961  !! at this point
1962  real, dimension(9), intent(out) :: dphidxi !< The x-gradient of the weights of the 9 bilinear
1963  !! quadrature points corners at this point
1964  real, dimension(9), intent(out) :: dphideta !< The z-gradient of the weights of the 9 bilinear
1965  !! quadrature points corners at this point
1966 
1967  ! The quadratic shape functions within the parent element are defined as shown here:
1968  !
1969  ! 5 (0,1)
1970  ! (-1,1) 2 o------o------o 1 (1,1)
1971  ! | |
1972  ! | 9 (0,0) |
1973  ! (-1,0) 6 o o o 8 (1,0)
1974  ! | |
1975  ! | |
1976  ! (-1,-1) 3 o------o------o 4 (1,-1)
1977  ! 7 (0,-1)
1978  !
1979 
1980  phi(:) = 0.0
1981  dphidxi(:) = 0.0
1982  dphideta(:) = 0.0
1983 
1984  phi(1) = 0.25 * xi * ( 1 + xi ) * eta * ( 1 + eta )
1985  phi(2) = - 0.25 * xi * ( 1 - xi ) * eta * ( 1 + eta )
1986  phi(3) = 0.25 * xi * ( 1 - xi ) * eta * ( 1 - eta )
1987  phi(4) = - 0.25 * xi * ( 1 + xi ) * eta * ( 1 - eta )
1988  phi(5) = 0.5 * ( 1 + xi ) * ( 1 - xi ) * eta * ( 1 + eta )
1989  phi(6) = - 0.5 * xi * ( 1 - xi ) * ( 1 - eta ) * ( 1 + eta )
1990  phi(7) = - 0.5 * ( 1 - xi ) * ( 1 + xi ) * eta * ( 1 - eta )
1991  phi(8) = 0.5 * xi * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
1992  phi(9) = ( 1 - xi ) * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
1993 
1994  !dphidxi(1) = 0.25 * ( 1 + 2*xi ) * eta * ( 1 + eta )
1995  !dphidxi(2) = - 0.25 * ( 1 - 2*xi ) * eta * ( 1 + eta )
1996  !dphidxi(3) = 0.25 * ( 1 - 2*xi ) * eta * ( 1 - eta )
1997  !dphidxi(4) = - 0.25 * ( 1 + 2*xi ) * eta * ( 1 - eta )
1998  !dphidxi(5) = - xi * eta * ( 1 + eta )
1999  !dphidxi(6) = - 0.5 * ( 1 - 2*xi ) * ( 1 - eta ) * ( 1 + eta )
2000  !dphidxi(7) = xi * eta * ( 1 - eta )
2001  !dphidxi(8) = 0.5 * ( 1 + 2*xi ) * ( 1 - eta ) * ( 1 + eta )
2002  !dphidxi(9) = - 2 * xi * ( 1 - eta ) * ( 1 + eta )
2003 
2004  !dphideta(1) = 0.25 * xi * ( 1 + xi ) * ( 1 + 2*eta )
2005  !dphideta(2) = - 0.25 * xi * ( 1 - xi ) * ( 1 + 2*eta )
2006  !dphideta(3) = 0.25 * xi * ( 1 - xi ) * ( 1 - 2*eta )
2007  !dphideta(4) = - 0.25 * xi * ( 1 + xi ) * ( 1 - 2*eta )
2008  !dphideta(5) = 0.5 * ( 1 + xi ) * ( 1 - xi ) * ( 1 + 2*eta )
2009  !dphideta(6) = xi * ( 1 - xi ) * eta
2010  !dphideta(7) = - 0.5 * ( 1 - xi ) * ( 1 + xi ) * ( 1 - 2*eta )
2011  !dphideta(8) = - xi * ( 1 + xi ) * eta
2012  !dphideta(9) = - 2 * ( 1 - xi ) * ( 1 + xi ) * eta
2013 
2014 end subroutine evaluate_shape_quadratic
2015 ! ==============================================================================
2016 
2017 !> This subroutine calculates integrals of specific volume anomalies in
2018 !! pressure across layers, which are required for calculating the finite-volume
2019 !! form pressure accelerations in a non-Boussinesq model. There are essentially
2020 !! no free assumptions, apart from the use of Bode's rule quadrature to do the integrals.
2021 subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, &
2022  dza, intp_dza, intx_dza, inty_dza, halo_size, &
2023  bathyP, dP_neglect, useMassWghtInterp)
2024  type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
2025  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2026  intent(in) :: t !< Potential temperature of the layer [degC].
2027  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2028  intent(in) :: s !< Salinity of the layer [ppt].
2029  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2030  intent(in) :: p_t !< Pressure atop the layer [Pa].
2031  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2032  intent(in) :: p_b !< Pressure below the layer [Pa].
2033  real, intent(in) :: alpha_ref !< A mean specific volume that is
2034  !! subtracted out to reduce the magnitude of each of the
2035  !! integrals [m3 kg-1]. The calculation is mathematically
2036  !! identical with different values of alpha_ref, but alpha_ref
2037  !! alters the effects of roundoff, and answers do change.
2038  type(eos_type), pointer :: eos !< Equation of state structure
2039  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2040  intent(out) :: dza !< The change in the geopotential anomaly
2041  !! across the layer [m2 s-2].
2042  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2043  optional, intent(out) :: intp_dza !< The integral in pressure through the
2044  !! layer of the geopotential anomaly relative to the anomaly
2045  !! at the bottom of the layer [Pa m2 s-2].
2046  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
2047  optional, intent(out) :: intx_dza !< The integral in x of the difference
2048  !! between the geopotential anomaly at the top and bottom of
2049  !! the layer divided by the x grid spacing [m2 s-2].
2050  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
2051  optional, intent(out) :: inty_dza !< The integral in y of the difference
2052  !! between the geopotential anomaly at the top and bottom of
2053  !! the layer divided by the y grid spacing [m2 s-2].
2054  integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
2055  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2056  optional, intent(in) :: bathyp !< The pressure at the bathymetry [Pa]
2057  real, optional, intent(in) :: dp_neglect !< A miniscule pressure change with
2058  !! the same units as p_t (Pa?)
2059  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
2060  !! to interpolate T/S for top and bottom integrals.
2061 
2062 ! This subroutine calculates analytical and nearly-analytical integrals in
2063 ! pressure across layers of geopotential anomalies, which are required for
2064 ! calculating the finite-volume form pressure accelerations in a non-Boussinesq
2065 ! model. There are essentially no free assumptions, apart from the use of
2066 ! Bode's rule to do the horizontal integrals, and from a truncation in the
2067 ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
2068 
2069  real :: t5(5), s5(5), p5(5), a5(5)
2070  real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1].
2071  real :: dp ! The pressure change through a layer [Pa].
2072 ! real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa].
2073  real :: hwght ! A pressure-thickness below topography [Pa].
2074  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Pa].
2075  real :: idenom ! The inverse of the denominator in the weights [Pa-2].
2076  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
2077  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
2078  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
2079  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
2080  real :: intp(5) ! The integrals of specific volume with pressure at the
2081  ! 5 sub-column locations [m2 s-2].
2082  logical :: do_massweight ! Indicates whether to do mass weighting.
2083  real, parameter :: c1_90 = 1.0/90.0 ! A rational constant.
2084  integer :: isq, ieq, jsq, jeq, ish, ieh, jsh, jeh, i, j, m, n, halo
2085 
2086  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
2087  halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
2088  ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
2089  if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh); endif
2090  if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh); endif
2091 
2092  do_massweight = .false.
2093  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
2094  do_massweight = .true.
2095  if (.not.present(bathyp)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
2096  "bathyP must be present if useMassWghtInterp is present and true.")
2097  if (.not.present(dp_neglect)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
2098  "dP_neglect must be present if useMassWghtInterp is present and true.")
2099  endif ; endif
2100 
2101  do j=jsh,jeh ; do i=ish,ieh
2102  dp = p_b(i,j) - p_t(i,j)
2103  do n=1,5
2104  t5(n) = t(i,j) ; s5(n) = s(i,j)
2105  p5(n) = p_b(i,j) - 0.25*real(n-1)*dp
2106  enddo
2107  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
2108 
2109  ! Use Bode's rule to estimate the interface height anomaly change.
2110  alpha_anom = c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3))
2111  dza(i,j) = dp*alpha_anom
2112  ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of
2113  ! the interface height anomaly.
2114  if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
2115  (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
2116  enddo ; enddo
2117 
2118  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
2119  ! hWght is the distance measure by which the cell is violation of
2120  ! hydrostatic consistency. For large hWght we bias the interpolation of
2121  ! T & S along the top and bottom integrals, akin to thickness weighting.
2122  hwght = 0.0
2123  if (do_massweight) &
2124  hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
2125  if (hwght > 0.) then
2126  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2127  hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
2128  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2129  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2130  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2131  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2132  else
2133  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2134  endif
2135 
2136  intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
2137  do m=2,4
2138  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2139  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2140 
2141  ! T, S, and p are interpolated in the horizontal. The p interpolation
2142  ! is linear, but for T and S it may be thickness wekghted.
2143  p5(1) = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
2144  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
2145  t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
2146  s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
2147 
2148  do n=2,5
2149  t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - 0.25*dp
2150  enddo
2151  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
2152 
2153  ! Use Bode's rule to estimate the interface height anomaly change.
2154  intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
2155  12.0*a5(3)))
2156  enddo
2157  ! Use Bode's rule to integrate the interface height anomaly values in x.
2158  intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
2159  12.0*intp(3))
2160  enddo ; enddo ; endif
2161 
2162  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
2163  ! hWght is the distance measure by which the cell is violation of
2164  ! hydrostatic consistency. For large hWght we bias the interpolation of
2165  ! T & S along the top and bottom integrals, akin to thickness weighting.
2166  hwght = 0.0
2167  if (do_massweight) &
2168  hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
2169  if (hwght > 0.) then
2170  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2171  hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
2172  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2173  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2174  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2175  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2176  else
2177  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2178  endif
2179 
2180  intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
2181  do m=2,4
2182  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2183  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2184 
2185  ! T, S, and p are interpolated in the horizontal. The p interpolation
2186  ! is linear, but for T and S it may be thickness wekghted.
2187  p5(1) = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
2188  dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
2189  t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
2190  s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
2191  do n=2,5
2192  t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - 0.25*dp
2193  enddo
2194  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
2195 
2196  ! Use Bode's rule to estimate the interface height anomaly change.
2197  intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
2198  12.0*a5(3)))
2199  enddo
2200  ! Use Bode's rule to integrate the interface height anomaly values in y.
2201  inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
2202  12.0*intp(3))
2203  enddo ; enddo ; endif
2204 
2205 end subroutine int_spec_vol_dp_generic
2206 
2207 !> This subroutine calculates integrals of specific volume anomalies in
2208 !! pressure across layers, which are required for calculating the finite-volume
2209 !! form pressure accelerations in a non-Boussinesq model. There are essentially
2210 !! no free assumptions, apart from the use of Bode's rule quadrature to do the integrals.
2211 subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, &
2212  dP_neglect, bathyP, HI, EOS, dza, &
2213  intp_dza, intx_dza, inty_dza, useMassWghtInterp)
2214  type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
2215  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2216  intent(in) :: t_t !< Potential temperature at the top of the layer [degC].
2217  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2218  intent(in) :: t_b !< Potential temperature at the bottom of the layer [degC].
2219  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2220  intent(in) :: s_t !< Salinity at the top the layer [ppt].
2221  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2222  intent(in) :: s_b !< Salinity at the bottom the layer [ppt].
2223  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2224  intent(in) :: p_t !< Pressure atop the layer [Pa].
2225  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2226  intent(in) :: p_b !< Pressure below the layer [Pa].
2227  real, intent(in) :: alpha_ref !< A mean specific volume that is
2228  !! subtracted out to reduce the magnitude of each of the
2229  !! integrals [m3 kg-1]. The calculation is mathematically
2230  !! identical with different values of alpha_ref, but alpha_ref
2231  !! alters the effects of roundoff, and answers do change.
2232  real, intent(in) :: dp_neglect !< A miniscule pressure change with
2233  !! the same units as p_t (Pa?)
2234  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2235  intent(in) :: bathyp !< The pressure at the bathymetry [Pa]
2236  type(eos_type), pointer :: eos !< Equation of state structure
2237  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2238  intent(out) :: dza !< The change in the geopotential anomaly
2239  !! across the layer [m2 s-2].
2240  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
2241  optional, intent(out) :: intp_dza !< The integral in pressure through the
2242  !! layer of the geopotential anomaly relative to the anomaly
2243  !! at the bottom of the layer [Pa m2 s-2].
2244  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
2245  optional, intent(out) :: intx_dza !< The integral in x of the difference
2246  !! between the geopotential anomaly at the top and bottom of
2247  !! the layer divided by the x grid spacing [m2 s-2].
2248  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
2249  optional, intent(out) :: inty_dza !< The integral in y of the difference
2250  !! between the geopotential anomaly at the top and bottom of
2251  !! the layer divided by the y grid spacing [m2 s-2].
2252  logical, optional, intent(in) :: usemasswghtinterp !< If true, uses mass weighting
2253  !! to interpolate T/S for top and bottom integrals.
2254 
2255 ! This subroutine calculates analytical and nearly-analytical integrals in
2256 ! pressure across layers of geopotential anomalies, which are required for
2257 ! calculating the finite-volume form pressure accelerations in a non-Boussinesq
2258 ! model. There are essentially no free assumptions, apart from the use of
2259 ! Bode's rule to do the horizontal integrals, and from a truncation in the
2260 ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
2261 
2262  real, dimension(5) :: t5, s5, p5, a5
2263  real, dimension(15) :: t15, s15, p15, a15
2264  real :: wt_t(5), wt_b(5)
2265  real :: t_top, t_bot, s_top, s_bot, p_top, p_bot
2266 
2267  real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1].
2268  real :: dp ! The pressure change through a layer [Pa].
2269  real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa].
2270  real :: hwght ! A pressure-thickness below topography [Pa].
2271  real :: hl, hr ! Pressure-thicknesses of the columns to the left and right [Pa].
2272  real :: idenom ! The inverse of the denominator in the weights [Pa-2].
2273  real :: hwt_ll, hwt_lr ! hWt_LA is the weighted influence of A on the left column [nondim].
2274  real :: hwt_rl, hwt_rr ! hWt_RA is the weighted influence of A on the right column [nondim].
2275  real :: wt_l, wt_r ! The linear weights of the left and right columns [nondim].
2276  real :: wtt_l, wtt_r ! The weights for tracers from the left and right columns [nondim].
2277  real :: intp(5) ! The integrals of specific volume with pressure at the
2278  ! 5 sub-column locations [m2 s-2].
2279  real, parameter :: c1_90 = 1.0/90.0 ! A rational constant.
2280  logical :: do_massweight ! Indicates whether to do mass weighting.
2281  integer :: isq, ieq, jsq, jeq, i, j, m, n, pos
2282 
2283  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
2284 
2285  do_massweight = .false.
2286  if (present(usemasswghtinterp)) do_massweight = usemasswghtinterp
2287 
2288  do n = 1, 5 ! Note that these are reversed from int_density_dz.
2289  wt_t(n) = 0.25 * real(n-1)
2290  wt_b(n) = 1.0 - wt_t(n)
2291  enddo
2292 
2293  ! =============================
2294  ! 1. Compute vertical integrals
2295  ! =============================
2296  do j=jsq,jeq+1; do i=isq,ieq+1
2297  dp = p_b(i,j) - p_t(i,j)
2298  do n=1,5 ! T, S and p are linearly interpolated in the vertical.
2299  p5(n) = wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j)
2300  s5(n) = wt_t(n) * s_t(i,j) + wt_b(n) * s_b(i,j)
2301  t5(n) = wt_t(n) * t_t(i,j) + wt_b(n) * t_b(i,j)
2302  enddo
2303  call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
2304 
2305  ! Use Bode's rule to estimate the interface height anomaly change.
2306  alpha_anom = c1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3))
2307  dza(i,j) = dp*alpha_anom
2308  ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of
2309  ! the interface height anomaly.
2310  if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
2311  (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
2312  enddo ; enddo
2313 
2314  ! ==================================================
2315  ! 2. Compute horizontal integrals in the x direction
2316  ! ==================================================
2317  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
2318  ! hWght is the distance measure by which the cell is violation of
2319  ! hydrostatic consistency. For large hWght we bias the interpolation
2320  ! of T,S along the top and bottom integrals, almost like thickness
2321  ! weighting. Note: To work in terrain following coordinates we could
2322  ! offset this distance by the layer thickness to replicate other models.
2323  hwght = 0.0
2324  if (do_massweight) &
2325  hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
2326  if (hwght > 0.) then
2327  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2328  hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
2329  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2330  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2331  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2332  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2333  else
2334  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2335  endif
2336 
2337  do m=2,4
2338  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2339  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2340 
2341  ! T, S, and p are interpolated in the horizontal. The p interpolation
2342  ! is linear, but for T and S it may be thickness wekghted.
2343  p_top = wt_l*p_t(i,j) + wt_r*p_t(i+1,j)
2344  p_bot = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
2345  t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i+1,j)
2346  t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i+1,j)
2347  s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i+1,j)
2348  s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i+1,j)
2349  dp_90(m) = c1_90*(p_bot - p_top)
2350 
2351  ! Salinity, temperature and pressure with linear interpolation in the vertical.
2352  pos = (m-2)*5
2353  do n=1,5
2354  p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot
2355  s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
2356  t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
2357  enddo
2358  enddo
2359 
2360  call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref)
2361 
2362  intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
2363  do m=2,4
2364  ! Use Bode's rule to estimate the interface height anomaly change.
2365  ! The integrals at the ends of the segment are already known.
2366  pos = (m-2)*5
2367  intp(m) = dp_90(m)*((7.0*(a15(pos+1)+a15(pos+5)) + &
2368  32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
2369  enddo
2370  ! Use Bode's rule to integrate the interface height anomaly values in x.
2371  intx_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
2372  12.0*intp(3))
2373  enddo ; enddo ; endif
2374 
2375  ! ==================================================
2376  ! 3. Compute horizontal integrals in the y direction
2377  ! ==================================================
2378  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
2379  ! hWght is the distance measure by which the cell is violation of
2380  ! hydrostatic consistency. For large hWght we bias the interpolation
2381  ! of T,S along the top and bottom integrals, like thickness weighting.
2382  hwght = 0.0
2383  if (do_massweight) &
2384  hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
2385  if (hwght > 0.) then
2386  hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
2387  hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
2388  hwght = hwght * ( (hl-hr)/(hl+hr) )**2
2389  idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
2390  hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
2391  hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
2392  else
2393  hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
2394  endif
2395 
2396  do m=2,4
2397  wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
2398  wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
2399 
2400  ! T, S, and p are interpolated in the horizontal. The p interpolation
2401  ! is linear, but for T and S it may be thickness wekghted.
2402  p_top = wt_l*p_t(i,j) + wt_r*p_t(i,j+1)
2403  p_bot = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
2404  t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i,j+1)
2405  t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i,j+1)
2406  s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i,j+1)
2407  s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i,j+1)
2408  dp_90(m) = c1_90*(p_bot - p_top)
2409 
2410  ! Salinity, temperature and pressure with linear interpolation in the vertical.
2411  pos = (m-2)*5
2412  do n=1,5
2413  p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot
2414  s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
2415  t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
2416  enddo
2417  enddo
2418 
2419  call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref)
2420 
2421  intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
2422  do m=2,4
2423  ! Use Bode's rule to estimate the interface height anomaly change.
2424  ! The integrals at the ends of the segment are already known.
2425  pos = (m-2)*5
2426  intp(m) = dp_90(m) * ((7.0*(a15(pos+1)+a15(pos+5)) + &
2427  32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
2428  enddo
2429  ! Use Bode's rule to integrate the interface height anomaly values in x.
2430  inty_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
2431  12.0*intp(3))
2432  enddo ; enddo ; endif
2433 
2434 end subroutine int_spec_vol_dp_generic_plm
2435 
2436 !> Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10
2437 subroutine convert_temp_salt_for_teos10(T, S, press, G, kd, mask_z, EOS)
2439 
2440  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
2441  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
2442  intent(inout) :: t !< Potential temperature referenced to the surface [degC]
2443  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
2444  intent(inout) :: s !< Salinity [ppt]
2445  real, dimension(:), intent(in) :: press !< Pressure at the top of the layer [Pa].
2446  type(eos_type), pointer :: eos !< Equation of state structure
2447  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
2448  intent(in) :: mask_z !< 3d mask regulating which points to convert.
2449  integer, intent(in) :: kd !< The number of layers to work on
2450 
2451  integer :: i,j,k
2452  real :: gsw_sr_from_sp, gsw_ct_from_pt, gsw_sa_from_sp
2453  real :: p
2454 
2455  if (.not.associated(eos)) call mom_error(fatal, &
2456  "convert_temp_salt_to_TEOS10 called with an unassociated EOS_type EOS.")
2457 
2458  if ((eos%form_of_EOS /= eos_teos10) .and. (eos%form_of_EOS /= eos_nemo)) return
2459 
2460  do k=1,kd ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
2461  if (mask_z(i,j,k) >= 1.0) then
2462  s(i,j,k) = gsw_sr_from_sp(s(i,j,k))
2463 ! p=press(k)/10000. !convert pascal to dbar
2464 ! S(i,j,k) = gsw_sa_from_sp(S(i,j,k),p,G%geoLonT(i,j),G%geoLatT(i,j))
2465  t(i,j,k) = gsw_ct_from_pt(s(i,j,k),t(i,j,k))
2466  endif
2467  enddo ; enddo ; enddo
2468 end subroutine convert_temp_salt_for_teos10
2469 
2470 !> Extractor routine for the EOS type if the members need to be accessed outside this module
2471 subroutine extract_member_eos(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, &
2472  Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
2473  type(eos_type), pointer :: eos !< Equation of state structure
2474  integer, optional, intent(out) :: form_of_eos !< A coded integer indicating the equation of state to use.
2475  integer, optional, intent(out) :: form_of_tfreeze !< A coded integer indicating the expression for
2476  !! the potential temperature of the freezing point.
2477  logical, optional, intent(out) :: eos_quadrature !< If true, always use the generic (quadrature)
2478  !! code for the integrals of density.
2479  logical, optional, intent(out) :: compressible !< If true, in situ density is a function of pressure.
2480  real , optional, intent(out) :: rho_t0_s0 !< Density at T=0 degC and S=0 ppt [kg m-3]
2481  real , optional, intent(out) :: drho_dt !< Partial derivative of density with temperature
2482  !! in [kg m-3 degC-1]
2483  real , optional, intent(out) :: drho_ds !< Partial derivative of density with salinity
2484  !! in [kg m-3 ppt-1]
2485  real , optional, intent(out) :: tfr_s0_p0 !< The freezing potential temperature at S=0, P=0 [degC].
2486  real , optional, intent(out) :: dtfr_ds !< The derivative of freezing point with salinity
2487  !! [degC PSU-1].
2488  real , optional, intent(out) :: dtfr_dp !< The derivative of freezing point with pressure
2489  !! [degC Pa-1].
2490 
2491  if (present(form_of_eos )) form_of_eos = eos%form_of_EOS
2492  if (present(form_of_tfreeze)) form_of_tfreeze = eos%form_of_TFreeze
2493  if (present(eos_quadrature )) eos_quadrature = eos%EOS_quadrature
2494  if (present(compressible )) compressible = eos%Compressible
2495  if (present(rho_t0_s0 )) rho_t0_s0 = eos%Rho_T0_S0
2496  if (present(drho_dt )) drho_dt = eos%drho_dT
2497  if (present(drho_ds )) drho_ds = eos%dRho_dS
2498  if (present(tfr_s0_p0 )) tfr_s0_p0 = eos%TFr_S0_P0
2499  if (present(dtfr_ds )) dtfr_ds = eos%dTFr_dS
2500  if (present(dtfr_dp )) dtfr_dp = eos%dTFr_dp
2501 
2502 end subroutine extract_member_eos
2503 
2504 end module mom_eos
2505 
2506 !> \namespace mom_eos
2507 !!
2508 !! The MOM_EOS module is a wrapper for various equations of state (e.g. Linear,
2509 !! Wright, UNESCO) and provides a uniform interface to the rest of the model
2510 !! independent of which equation of state is being used.
mom_eos::int_density_dz_generic_plm
subroutine, public int_density_dz_generic_plm(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp)
Compute pressure gradient force integrals by quadrature for the case where T and S are linear profile...
Definition: MOM_EOS.F90:1171
mom_eos_wright::calculate_specvol_derivs_wright
subroutine, public calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts)
For a given thermodynamic state, return the partial derivatives of specific volume with temperature a...
Definition: MOM_EOS_Wright.F90:343
mom_eos_linear::int_density_dz_linear
subroutine, public int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HII, HIO, Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
This subroutine calculates analytical and nearly-analytical integrals of pressure anomalies across la...
Definition: MOM_EOS_linear.F90:331
mom_eos_unesco::calculate_spec_vol_unesco
Compute the in situ specific volume of sea water (in [m3 kg-1]), or an anomaly with respect to a refe...
Definition: MOM_EOS_UNESCO.F90:29
mom_eos_linear::calculate_density_second_derivs_linear
For a given thermodynamic state, return the second derivatives of density with various combinations o...
Definition: MOM_EOS_linear.F90:47
mom_eos::evaluate_shape_bilinear
subroutine evaluate_shape_bilinear(xi, eta, phi, dphidxi, dphideta)
Evaluation of the four bilinear shape fn and their gradients at (xi,eta)
Definition: MOM_EOS.F90:1917
mom_eos::convert_temp_salt_for_teos10
subroutine, public convert_temp_salt_for_teos10(T, S, press, G, kd, mask_z, EOS)
Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10.
Definition: MOM_EOS.F90:2438
mom_eos::eos_end
subroutine, public eos_end(EOS)
Deallocates EOS_type.
Definition: MOM_EOS.F90:947
mom_eos_unesco
The equation of state using the Jackett and McDougall fits to the UNESCO EOS.
Definition: MOM_EOS_UNESCO.F90:2
mom_eos::int_specific_vol_dp
subroutine, public int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp)
Calls the appropriate subroutine to alculate analytical and nearly-analytical integrals in pressure a...
Definition: MOM_EOS.F90:665
mom_eos::calculate_density_derivs_array
subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS, scale)
Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.
Definition: MOM_EOS.F90:367
mom_eos::eos_wright_string
character *(10), parameter eos_wright_string
A string for specifying the equation of state.
Definition: MOM_EOS.F90:120
mom_eos_linear::calculate_compress_linear
subroutine, public calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts, Rho_T0_S0, dRho_dT, dRho_dS)
This subroutine computes the in situ density of sea water (rho) and the compressibility (drho/dp == C...
Definition: MOM_EOS_linear.F90:301
mom_eos::query_compressible
logical function, public query_compressible(EOS)
Returns true if the equation of state is compressible (i.e. has pressure dependence)
Definition: MOM_EOS.F90:796
mom_eos::eos_linear_string
character *(10), parameter eos_linear_string
A string for specifying the equation of state.
Definition: MOM_EOS.F90:118
mom_eos::calculate_spec_vol_scalar
subroutine calculate_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, scale)
Calls the appropriate subroutine to calculate specific volume of sea water for scalar inputs.
Definition: MOM_EOS.F90:218
mom_eos_wright::int_density_dz_wright
subroutine, public int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
This subroutine calculates analytical and nearly-analytical integrals of pressure anomalies across la...
Definition: MOM_EOS_Wright.F90:412
mom_eos::int_spec_vol_dp_generic
subroutine, public int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, useMassWghtInterp)
This subroutine calculates integrals of specific volume anomalies in pressure across layers,...
Definition: MOM_EOS.F90:2024
mom_eos::eos_use_linear
subroutine, public eos_use_linear(Rho_T0_S0, dRho_dT, dRho_dS, EOS, use_quadrature)
Set equation of state structure (EOS) to linear with given coefficients.
Definition: MOM_EOS.F90:958
mom_eos::calculate_compress
Calculates the compressibility of water from T, S, and P.
Definition: MOM_EOS.F90:86
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_eos::int_density_dz_generic
subroutine, public int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
This subroutine calculates (by numerical quadrature) integrals of pressure anomalies across layers,...
Definition: MOM_EOS.F90:984
mom_eos::calculate_tfreeze_scalar
subroutine calculate_tfreeze_scalar(S, pressure, T_fr, EOS)
Calls the appropriate subroutine to calculate the freezing point for scalar inputs.
Definition: MOM_EOS.F90:313
mom_eos_linear::calculate_density_derivs_linear
For a given thermodynamic state, return the derivatives of density with temperature and salinity usin...
Definition: MOM_EOS_linear.F90:40
mom_eos::tfreeze_millero_string
character *(10), parameter tfreeze_millero_string
A string for specifying freezing point expression.
Definition: MOM_EOS.F90:129
mom_eos_linear::calculate_specvol_derivs_linear
subroutine, public calculate_specvol_derivs_linear(T, S, pressure, dSV_dT, dSV_dS, start, npts, Rho_T0_S0, dRho_dT, dRho_dS)
Calculate the derivatives of specific volume with temperature and salinity.
Definition: MOM_EOS_linear.F90:268
mom_eos_linear
A simple linear equation of state for sea water with constant coefficients.
Definition: MOM_EOS_linear.F90:2
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_eos::int_density_dz
subroutine, public int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
This subroutine calculates analytical and nearly-analytical integrals of pressure anomalies across la...
Definition: MOM_EOS.F90:733
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_eos::eos_nemo
integer, parameter, public eos_nemo
A named integer specifying an equation of state.
Definition: MOM_EOS.F90:116
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_eos_teos10::calculate_spec_vol_teos10
Compute the in situ specific volume of sea water (in [m3 kg-1]), or an anomaly with respect to a refe...
Definition: MOM_EOS_TEOS10.F90:35
mom_eos::int_density_dz_generic_ppm
subroutine, public int_density_dz_generic_ppm(T, T_t, T_b, S, S_t, S_b, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
Compute pressure gradient force integrals for the case where T and S are parabolic profiles.
Definition: MOM_EOS.F90:1585
mom_eos::calculate_density_scalar
subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, scale)
Calls the appropriate subroutine to calculate density of sea water for scalar inputs....
Definition: MOM_EOS.F90:139
mom_string_functions::uppercase
character(len=len(input_string)) function, public uppercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:42
mom_eos::eos_teos10_string
character *(10), parameter eos_teos10_string
A string for specifying the equation of state.
Definition: MOM_EOS.F90:121
mom_eos::calculate_compress_scalar
subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS)
Calculate density and compressibility for a scalar. This just promotes the scalar to an array with a ...
Definition: MOM_EOS.F90:638
mom_eos::calculate_tfreeze
Calculates the freezing point of sea water from T, S and P.
Definition: MOM_EOS.F90:81
mom_eos_linear::int_spec_vol_dp_linear
subroutine, public int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, dRho_dT, dRho_dS, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, useMassWghtInterp)
Calculates analytical and nearly-analytical integrals in pressure across layers of geopotential anoma...
Definition: MOM_EOS_linear.F90:505
mom_eos_teos10::calculate_density_teos10
Compute the in situ density of sea water ([kg m-3]), or its anomaly with respect to a reference densi...
Definition: MOM_EOS_TEOS10.F90:28
mom_eos_wright
The equation of state using the Wright 1997 expressions.
Definition: MOM_EOS_Wright.F90:2
mom_eos::eos_teos10
integer, parameter, public eos_teos10
A named integer specifying an equation of state.
Definition: MOM_EOS.F90:115
mom_eos_wright::calculate_compress_wright
subroutine, public calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts)
This subroutine computes the in situ density of sea water (rho in [kg m-3]) and the compressibility (...
Definition: MOM_EOS_Wright.F90:380
mom_eos::int_spec_vol_dp_generic_plm
subroutine, public int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, dP_neglect, bathyP, HI, EOS, dza, intp_dza, intx_dza, inty_dza, useMassWghtInterp)
This subroutine calculates integrals of specific volume anomalies in pressure across layers,...
Definition: MOM_EOS.F90:2214
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
mom_eos_teos10
The equation of state using the TEOS10 expressions.
Definition: MOM_EOS_TEOS10.F90:2
mom_eos::tfreeze_linear_string
character *(10), parameter tfreeze_linear_string
A string for specifying the freezing point expression.
Definition: MOM_EOS.F90:128
mom_eos::tfreeze_millero
integer, parameter tfreeze_millero
A named integer specifying a freezing point expression.
Definition: MOM_EOS.F90:126
mom_eos_nemo
The equation of state using the expressions of Roquet et al. that are used in NEMO.
Definition: MOM_EOS_NEMO.F90:2
mom_eos_wright::calculate_density_wright
Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to a reference de...
Definition: MOM_EOS_Wright.F90:32
mom_tfreeze
Freezing point expressions.
Definition: MOM_TFreeze.F90:2
mom_eos::calculate_spec_vol_array
subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref, scale)
Calls the appropriate subroutine to calculate the specific volume of sea water for 1-D array inputs.
Definition: MOM_EOS.F90:264
mom_tfreeze::calculate_tfreeze_millero
Compute the freezing point potential temperature [degC] from salinity [PSU] and pressure [Pa] using t...
Definition: MOM_TFreeze.F90:27
mom_eos_nemo::calculate_density_derivs_nemo
For a given thermodynamic state, return the derivatives of density with conservative temperature and ...
Definition: MOM_EOS_NEMO.F90:34
mom_eos_nemo::calculate_density_nemo
Compute the in situ density of sea water ([kg m-3]), or its anomaly with respect to a reference densi...
Definition: MOM_EOS_NEMO.F90:28
mom_eos::eos_linear
integer, parameter, public eos_linear
A named integer specifying an equation of state.
Definition: MOM_EOS.F90:112
mom_tfreeze::calculate_tfreeze_linear
Compute the freezing point potential temperature [degC] from salinity [ppt] and pressure [Pa] using a...
Definition: MOM_TFreeze.F90:18
mom_eos::calculate_density_second_derivs
Calculates the second derivatives of density with various combinations of temperature,...
Definition: MOM_EOS.F90:76
mom_eos::tfreeze_teos10_string
character *(10), parameter tfreeze_teos10_string
A string for specifying the freezing point expression.
Definition: MOM_EOS.F90:131
mom_eos::eos_unesco
integer, parameter, public eos_unesco
A named integer specifying an equation of state.
Definition: MOM_EOS.F90:113
mom_eos_teos10::calculate_density_second_derivs_teos10
For a given thermodynamic state, return the second derivatives of density with various combinations o...
Definition: MOM_EOS_TEOS10.F90:47
mom_eos::calculate_density_second_derivs_scalar
subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, EOS, scale)
Calls the appropriate subroutine to calculate density second derivatives for scalar nputs.
Definition: MOM_EOS.F90:499
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_eos_linear::calculate_spec_vol_linear
Compute the specific volume of sea water (in m^3/kg), or its anomaly from a reference value,...
Definition: MOM_EOS_linear.F90:34
mom_eos_teos10::calculate_compress_teos10
subroutine, public calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts)
This subroutine computes the in situ density of sea water (rho in [kg m-3]) and the compressibility (...
Definition: MOM_EOS_TEOS10.F90:314
mom_eos_unesco::calculate_compress_unesco
subroutine, public calculate_compress_unesco(T, S, pressure, rho, drho_dp, start, npts)
This subroutine computes the in situ density of sea water (rho) and the compressibility (drho/dp == C...
Definition: MOM_EOS_UNESCO.F90:284
mom_eos_teos10::calculate_density_derivs_teos10
For a given thermodynamic state, return the derivatives of density with conservative temperature and ...
Definition: MOM_EOS_TEOS10.F90:41
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_eos_linear::calculate_density_linear
Compute the density of sea water (in kg/m^3), or its anomaly from a reference density,...
Definition: MOM_EOS_linear.F90:27
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:15
mom_eos_wright::calculate_density_derivs_wright
For a given thermodynamic state, return the derivatives of density with temperature and salinity.
Definition: MOM_EOS_Wright.F90:44
mom_eos::eos_unesco_string
character *(10), parameter eos_unesco_string
A string for specifying the equation of state.
Definition: MOM_EOS.F90:119
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_eos::evaluate_shape_quadratic
subroutine evaluate_shape_quadratic(xi, eta, phi, dphidxi, dphideta)
Evaluation of the nine quadratic shape fn weights and their gradients at (xi,eta)
Definition: MOM_EOS.F90:1956
mom_eos::eos_allocate
subroutine, public eos_allocate(EOS)
Allocates EOS_type.
Definition: MOM_EOS.F90:940
mom_eos::eos_default
character *(10), parameter eos_default
The default equation of state.
Definition: MOM_EOS.F90:123
mom_eos::eos_manual_init
subroutine, public eos_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
Manually initialized an EOS type (intended for unit testing of routines which need a specific EOS)
Definition: MOM_EOS.F90:907
mom_eos::compute_integral_quadratic
subroutine compute_integral_quadratic(x, y, f, integral)
Compute the integral of the quadratic function.
Definition: MOM_EOS.F90:1834
mom_eos::calculate_density_derivs_scalar
subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS, scale)
Calls the appropriate subroutines to calculate density derivatives by promoting a scalar to a one-ele...
Definition: MOM_EOS.F90:411
mom_eos::tfreeze_linear
integer, parameter tfreeze_linear
A named integer specifying a freezing point expression.
Definition: MOM_EOS.F90:125
mom_eos_nemo::calculate_compress_nemo
subroutine, public calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts)
Compute the in situ density of sea water (rho in [kg m-3]) and the compressibility (drho/dp = C_sound...
Definition: MOM_EOS_NEMO.F90:369
mom_eos::calculate_specific_vol_derivs
subroutine, public calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS, scale)
Calls the appropriate subroutine to calculate specific volume derivatives for an array.
Definition: MOM_EOS.F90:546
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_eos::calculate_tfreeze_array
subroutine calculate_tfreeze_array(S, pressure, T_fr, start, npts, EOS)
Calls the appropriate subroutine to calculate the freezing point for a 1-D array.
Definition: MOM_EOS.F90:339
mom_eos::frac_dp_at_pos
real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS)
Returns change in anomalous pressure change from top to non-dimensional position pos between z_t and ...
Definition: MOM_EOS.F90:1542
mom_eos::find_depth_of_pressure_in_cell
subroutine, public find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, rho_ref, G_e, EOS, P_b, z_out, z_tol)
Find the depth at which the reconstructed pressure matches P_tgt.
Definition: MOM_EOS.F90:1470
mom_eos::calculate_density_second_derivs_array
subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, start, npts, EOS, scale)
Calls the appropriate subroutine to calculate density second derivatives for 1-D array inputs.
Definition: MOM_EOS.F90:448
mom_eos_wright::calculate_spec_vol_wright
Compute the in situ specific volume of sea water (in [m3 kg-1]), or an anomaly with respect to a refe...
Definition: MOM_EOS_Wright.F90:39
mom_eos::eos_init
subroutine, public eos_init(param_file, EOS)
Initializes EOS_type by allocating and reading parameters.
Definition: MOM_EOS.F90:806
mom_eos_wright::calculate_density_second_derivs_wright
For a given thermodynamic state, return the second derivatives of density with various combinations o...
Definition: MOM_EOS_Wright.F90:50
mom_eos::extract_member_eos
subroutine, public extract_member_eos(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
Extractor routine for the EOS type if the members need to be accessed outside this module.
Definition: MOM_EOS.F90:2473
mom_eos::tfreeze_default
character *(10), parameter tfreeze_default
The default freezing point expression.
Definition: MOM_EOS.F90:132
mom_eos_unesco::calculate_density_unesco
Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to a reference de...
Definition: MOM_EOS_UNESCO.F90:22
mom_eos::tfreeze_teos10
integer, parameter tfreeze_teos10
A named integer specifying a freezing point expression.
Definition: MOM_EOS.F90:127
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_eos_teos10::calculate_specvol_derivs_teos10
subroutine, public calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts)
For a given thermodynamic state, calculate the derivatives of specific volume with conservative tempe...
Definition: MOM_EOS_TEOS10.F90:221
mom_tfreeze::calculate_tfreeze_teos10
Compute the freezing point conservative temperature [degC] from absolute salinity [g/kg] and pressure...
Definition: MOM_TFreeze.F90:33
mom_eos::calculate_spec_vol
Calculates specific volume of sea water from T, S and P.
Definition: MOM_EOS.F90:65
mom_eos_wright::int_spec_vol_dp_wright
subroutine, public int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, useMassWghtInterp)
This subroutine calculates analytical and nearly-analytical integrals in pressure across layers of ge...
Definition: MOM_EOS_Wright.F90:617
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_eos::eos_nemo_string
character *(10), parameter eos_nemo_string
A string for specifying the equation of state.
Definition: MOM_EOS.F90:122
mom_eos::calculate_compress_array
subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, EOS)
Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs.
Definition: MOM_EOS.F90:603
mom_eos::calculate_density_array
subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_ref, scale)
Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs....
Definition: MOM_EOS.F90:177
mom_eos_unesco::calculate_density_derivs_unesco
subroutine, public calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts)
This subroutine calculates the partial derivatives of density with potential temperature and salinity...
Definition: MOM_EOS_UNESCO.F90:213
mom_eos::eos_wright
integer, parameter, public eos_wright
A named integer specifying an equation of state.
Definition: MOM_EOS.F90:114
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60