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