MOM6
Data Types | Functions/Subroutines | Variables
mom_eos Module Reference

Detailed Description

Provides subroutines for quantities specific to the equation of state.

The MOM_EOS module is a wrapper for various equations of state (e.g. Linear, Wright, UNESCO) and provides a uniform interface to the rest of the model independent of which equation of state is being used.

Data Types

interface  calculate_density
 Calculates density of sea water from T, S and P. More...
 
interface  calculate_tfreeze
 Calculates the freezing point of sea water from T, S and P. More...
 
type  eos_type
 A control structure for the equation of state. More...
 

Functions/Subroutines

subroutine calculate_density_scalar (T, S, pressure, rho, EOS)
 Calls the appropriate subroutine to calculate density of sea water for scalar inputs. More...
 
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. More...
 
subroutine calculate_tfreeze_scalar (S, pressure, T_fr, EOS)
 Calls the appropriate subroutine to calculate the freezing point for scalar inputs. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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 across layers of geopotential anomalies, which are required for calculating the finite-volume form pressure accelerations in a non-Boussinesq model. There are essentially no free assumptions, apart from the use of Bode's rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) that assumes that |eps| < . More...
 
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 layers, which are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. The one potentially dodgy assumtion here is that rho_0 is used both in the denominator of the accelerations, and in the pressure used to calculated density (the latter being -z*rho_0*G_e). These two uses could be separated if need be. More...
 
logical function, public query_compressible (EOS)
 Returns true if the equation of state is compressible (i.e. has pressure dependence) More...
 
subroutine, public eos_init (param_file, EOS)
 Initializes EOS_type by allocating and reading parameters. More...
 
subroutine, public eos_allocate (EOS)
 Allocates EOS_type. More...
 
subroutine, public eos_end (EOS)
 Deallocates EOS_type. More...
 
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. More...
 
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)
 
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)
 
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)
 
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. More...
 
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 z_b. More...
 
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)
 
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)
 
subroutine compute_integral_bilinear (x, y, f, integral)
 
subroutine compute_integral_quadratic (x, y, f, integral)
 
subroutine evaluate_shape_bilinear (xi, eta, phi, dphidxi, dphideta)
 
subroutine evaluate_shape_quadratic (xi, eta, phi, dphidxi, dphideta)
 
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)
 
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. More...
 

Variables

integer, parameter eos_linear = 1
 
integer, parameter eos_unesco = 2
 
integer, parameter eos_wright = 3
 
integer, parameter eos_teos10 = 4
 
integer, parameter eos_nemo = 5
 
character *(10), parameter eos_linear_string = "LINEAR"
 
character *(10), parameter eos_unesco_string = "UNESCO"
 
character *(10), parameter eos_wright_string = "WRIGHT"
 
character *(10), parameter eos_teos10_string = "TEOS10"
 
character *(10), parameter eos_nemo_string = "NEMO"
 
character *(10), parameter eos_default = EOS_WRIGHT_STRING
 
integer, parameter tfreeze_linear = 1
 
integer, parameter tfreeze_millero = 2
 
integer, parameter tfreeze_teos10 = 3
 
character *(10), parameter tfreeze_linear_string = "LINEAR"
 
character *(10), parameter tfreeze_millero_string = "MILLERO_78"
 
character *(10), parameter tfreeze_teos10_string = "TEOS10"
 
character *(10), parameter tfreeze_default = TFREEZE_LINEAR_STRING
 

Function/Subroutine Documentation

◆ calculate_compress()

subroutine, public mom_eos::calculate_compress ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(out)  rho,
real, dimension(:), intent(out)  drho_dp,
integer, intent(in)  start,
integer, intent(in)  npts,
type(eos_type), pointer  EOS 
)

Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs.

Parameters
[in]tPotential temperature referenced to the surface (degC)
[in]sSalinity (PSU)
[in]pressurePressure (Pa)
[out]rhoIn situ density in kg m-3.
[out]drho_dpThe partial derivative of density with pressure (also the inverse of the square of sound speed) in s2 m-2.
[in]startStarting index within the array
[in]nptsThe number of values to calculate
eosEquation of state structure

Definition at line 293 of file MOM_EOS.F90.

References mom_eos_linear::calculate_compress_linear(), mom_eos_nemo::calculate_compress_nemo(), mom_eos_unesco::calculate_compress_unesco(), mom_eos_wright::calculate_compress_wright(), eos_linear, eos_nemo, eos_teos10, eos_unesco, eos_wright, and mom_error_handler::mom_error().

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 
Here is the call graph for this function:

◆ calculate_density_array()

subroutine mom_eos::calculate_density_array ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(out)  rho,
integer, intent(in)  start,
integer, intent(in)  npts,
type(eos_type), pointer  EOS 
)
private

Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs.

Parameters
[in]tPotential temperature referenced to the surface (degC)
[in]sSalinity (PSU)
[in]pressurePressure (Pa)
[out]rhoDensity (in-situ if pressure is local) (kg m-3)
[in]startStart index for computation
[in]nptsNumber of point to compute
eosEquation of state structure

Definition at line 130 of file MOM_EOS.F90.

References eos_linear, eos_nemo, eos_teos10, eos_unesco, eos_wright, and mom_error_handler::mom_error().

Referenced by frac_dp_at_pos(), and int_density_dz_generic_plm().

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 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ calculate_density_derivs()

subroutine, public mom_eos::calculate_density_derivs ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(out)  drho_dT,
real, dimension(:), intent(out)  drho_dS,
integer, intent(in)  start,
integer, intent(in)  npts,
type(eos_type), pointer  EOS 
)

Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.

Parameters
[in]tPotential temperature referenced to the surface (degC)
[in]sSalinity (PSU)
[in]pressurePressure (Pa)
[out]drho_dtThe partial derivative of density with potential tempetature, in kg m-3 K-1.
[out]drho_dsThe partial derivative of density with salinity, in kg m-3 psu-1.
[in]startStarting index within the array
[in]nptsThe number of values to calculate
eosEquation of state structure

Definition at line 214 of file MOM_EOS.F90.

References eos_linear, eos_nemo, eos_teos10, eos_unesco, eos_wright, and mom_error_handler::mom_error().

Referenced by mom_diabatic_aux::applyboundaryfluxesinout(), benchmark_initialization::benchmark_init_temperature_salinity(), benchmark_initialization::benchmark_initialize_thickness(), coord_adapt::build_adapt_column(), coord_slight::build_slight_column(), mom_bulk_mixed_layer::bulkmixedlayer(), mom_isopycnal_slopes::calc_isoneutral_slopes(), mom_kappa_shear::calculate_kappa_shear(), mom_forcing_type::calculatebuoyancyflux1d(), midas_vertmap::determine_temperature(), dome_initialization::dome_set_obc_data(), mom_entrain_diffusive::entrainment_diffusive(), mom_int_tide_input::find_n2_bottom(), mom_geothermal::geothermal(), isomip_initialization::isomip_initialize_temperature_salinity(), mom_neutral_diffusion::neutral_diffusion_calc_coeffs(), mom_set_visc::set_viscous_bbl(), mom_set_visc::set_viscous_ml(), mom_thickness_diffuse::thickness_diffuse_full(), mom_wave_speed::wave_speed(), and mom_wave_structure::wave_structure().

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 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ calculate_density_scalar()

subroutine mom_eos::calculate_density_scalar ( real, intent(in)  T,
real, intent(in)  S,
real, intent(in)  pressure,
real, intent(out)  rho,
type(eos_type), pointer  EOS 
)
private

Calls the appropriate subroutine to calculate density of sea water for scalar inputs.

Parameters
[in]tPotential temperature referenced to the surface (degC)
[in]sSalinity (PSU)
[in]pressurePressure (Pa)
[out]rhoDensity (in-situ if pressure is local) (kg m-3)
eosEquation of state structure

Definition at line 100 of file MOM_EOS.F90.

References eos_linear, eos_nemo, eos_teos10, eos_unesco, eos_wright, and mom_error_handler::mom_error().

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 
Here is the call graph for this function:

◆ calculate_specific_vol_derivs()

subroutine, public mom_eos::calculate_specific_vol_derivs ( real, dimension(:), intent(in)  T,
real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(out)  dSV_dT,
real, dimension(:), intent(out)  dSV_dS,
integer, intent(in)  start,
integer, intent(in)  npts,
type(eos_type), pointer  EOS 
)

Calls the appropriate subroutine to calculate specific volume derivatives for an array.

Parameters
[in]tPotential temperature referenced to the surface (degC)
[in]sSalinity (PSU)
[in]pressurePressure (Pa)
[out]dsv_dtThe partial derivative of specific volume with potential temperature, in m3 kg-1 K-1.
[out]dsv_dsThe partial derivative of specific volume with salinity, in m3 kg-1 / (g/kg).
[in]startStarting index within the array
[in]nptsThe number of values to calculate
eosEquation of state structure

Definition at line 247 of file MOM_EOS.F90.

References eos_linear, eos_nemo, eos_teos10, eos_unesco, eos_wright, and mom_error_handler::mom_error().

Referenced by mom_diabatic_aux::applyboundaryfluxesinout(), and mom_diapyc_energy_req::diapyc_energy_req_calc().

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 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ calculate_tfreeze_array()

subroutine mom_eos::calculate_tfreeze_array ( real, dimension(:), intent(in)  S,
real, dimension(:), intent(in)  pressure,
real, dimension(:), intent(out)  T_fr,
integer, intent(in)  start,
integer, intent(in)  npts,
type(eos_type), pointer  EOS 
)
private

Calls the appropriate subroutine to calculate the freezing point for a 1-D array.

Parameters
[in]sSalinity (PSU)
[in]pressurePressure (Pa)
[out]t_frFreezing point potential temperature referenced to the surface (degC)
[in]startStarting index within the array
[in]nptsThe number of values to calculate
eosEquation of state structure

Definition at line 187 of file MOM_EOS.F90.

References mom_error_handler::mom_error(), tfreeze_linear, tfreeze_millero, and tfreeze_teos10.

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 
Here is the call graph for this function:

◆ calculate_tfreeze_scalar()

subroutine mom_eos::calculate_tfreeze_scalar ( real, intent(in)  S,
real, intent(in)  pressure,
real, intent(out)  T_fr,
type(eos_type), pointer  EOS 
)
private

Calls the appropriate subroutine to calculate the freezing point for scalar inputs.

Parameters
[in]sSalinity (PSU)
[in]pressurePressure (Pa)
[out]t_frFreezing point potential temperature referenced to the surface (degC)
eosEquation of state structure

Definition at line 162 of file MOM_EOS.F90.

References mom_error_handler::mom_error(), tfreeze_linear, tfreeze_millero, and tfreeze_teos10.

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 
Here is the call graph for this function:

◆ compute_integral_bilinear()

subroutine mom_eos::compute_integral_bilinear ( real, dimension(4), intent(in)  x,
real, dimension(4), intent(in)  y,
real, dimension(4), intent(in)  f,
real, intent(out)  integral 
)
private

Definition at line 1756 of file MOM_EOS.F90.

References evaluate_shape_bilinear().

Referenced by int_density_dz_generic_plm_analytic().

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 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ compute_integral_quadratic()

subroutine mom_eos::compute_integral_quadratic ( real, dimension(4), intent(in)  x,
real, dimension(4), intent(in)  y,
real, dimension(9), intent(in)  f,
real, intent(out)  integral 
)
private

Definition at line 1822 of file MOM_EOS.F90.

References evaluate_shape_bilinear(), and evaluate_shape_quadratic().

Referenced by int_density_dz_generic_ppm().

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 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ convert_temp_salt_for_teos10()

subroutine, public mom_eos::convert_temp_salt_for_teos10 ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  S,
real, dimension(:), intent(in)  press,
type(ocean_grid_type), intent(in)  G,
integer, intent(in)  kd,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  mask_z,
type(eos_type), pointer  EOS 
)

Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10.

Parameters
[in]gThe horizontal index structure
[in]gThe ocean's grid structure
[in,out]tPotential temperature referenced to the surface (degC)
[in,out]sSalinity (PSU)
[in]pressPressure at the top of the layer in Pa.
eosEquation of state structure
[in]mask_z3d mask

Definition at line 2129 of file MOM_EOS.F90.

References eos_nemo, eos_teos10, and mom_error_handler::mom_error().

Referenced by mom_state_initialization::mom_temp_salt_initialize_from_z().

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
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
Here is the call graph for this function:
Here is the caller graph for this function:

◆ eos_allocate()

subroutine, public mom_eos::eos_allocate ( type(eos_type), pointer  EOS)

Allocates EOS_type.

Parameters
eosEquation of state structure

Definition at line 558 of file MOM_EOS.F90.

Referenced by eos_init().

558  type(eos_type), pointer :: eos !< Equation of state structure
559 
560  if (.not.associated(eos)) allocate(eos)
Here is the caller graph for this function:

◆ eos_end()

subroutine, public mom_eos::eos_end ( type(eos_type), pointer  EOS)

Deallocates EOS_type.

Parameters
eosEquation of state structure

Definition at line 565 of file MOM_EOS.F90.

565  type(eos_type), pointer :: eos !< Equation of state structure
566 
567  if (associated(eos)) deallocate(eos)

◆ eos_init()

subroutine, public mom_eos::eos_init ( type(param_file_type), intent(in)  param_file,
type(eos_type), pointer  EOS 
)

Initializes EOS_type by allocating and reading parameters.

Parameters
[in]param_fileParameter file structure
eosEquation of state structure

Definition at line 459 of file MOM_EOS.F90.

References eos_allocate(), eos_default, eos_linear, eos_linear_string, eos_nemo, eos_nemo_string, eos_teos10, eos_teos10_string, eos_unesco, eos_unesco_string, eos_wright, eos_wright_string, mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), tfreeze_default, tfreeze_linear, tfreeze_linear_string, tfreeze_millero, tfreeze_millero_string, tfreeze_teos10, tfreeze_teos10_string, and mom_string_functions::uppercase().

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 
Here is the call graph for this function:

◆ eos_use_linear()

subroutine, public mom_eos::eos_use_linear ( real, intent(in)  Rho_T0_S0,
real, intent(in)  dRho_dT,
real, intent(in)  dRho_dS,
type(eos_type), pointer  EOS,
logical, intent(in), optional  use_quadrature 
)

Set equation of state structure (EOS) to linear with given coefficients.

Note
This routine is primarily for testing and allows a local copy of the EOS_type (EOS argument) to be set to use the linear equation of state independent from the rest of the model.
Parameters
[in]rho_t0_s0Density at T=0 degC and S=0 ppt (kg m-3)
[in]drho_dtPartial derivative of density with temperature (kg m-3 degC-1)
[in]drho_dsPartial derivative of density with salinity (kg m-3 ppt-1)
[in]use_quadraturePartial derivative of density with salinity (kg m-3 ppt-1)
eosEquation of state structure

Definition at line 576 of file MOM_EOS.F90.

References eos_linear, and mom_error_handler::mom_error().

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 
Here is the call graph for this function:

◆ evaluate_shape_bilinear()

subroutine mom_eos::evaluate_shape_bilinear ( real, intent(in)  xi,
real, intent(in)  eta,
real, dimension(4), intent(out)  phi,
real, dimension(4), intent(out)  dphidxi,
real, dimension(4), intent(out)  dphideta 
)
private

Definition at line 1907 of file MOM_EOS.F90.

Referenced by compute_integral_bilinear(), and compute_integral_quadratic().

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 
Here is the caller graph for this function:

◆ evaluate_shape_quadratic()

subroutine mom_eos::evaluate_shape_quadratic ( real, intent(in)  xi,
real, intent(in)  eta,
real, dimension(9), intent(out)  phi,
real, dimension(9), intent(out)  dphidxi,
real, dimension(9), intent(out)  dphideta 
)
private

Definition at line 1945 of file MOM_EOS.F90.

Referenced by compute_integral_quadratic().

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 
Here is the caller graph for this function:

◆ find_depth_of_pressure_in_cell()

subroutine, public mom_eos::find_depth_of_pressure_in_cell ( real, intent(in)  T_t,
real, intent(in)  T_b,
real, intent(in)  S_t,
real, intent(in)  S_b,
real, intent(in)  z_t,
real, intent(in)  z_b,
real, intent(in)  P_t,
real, intent(in)  P_tgt,
real, intent(in)  rho_ref,
real, intent(in)  G_e,
type(eos_type), pointer  EOS,
real, intent(out)  P_b,
real, intent(out)  z_out 
)

Find the depth at which the reconstructed pressure matches P_tgt.

Parameters
[in]t_tPotential temperatue at the cell top (degC)
[in]t_bPotential temperatue at the cell bottom (degC)
[in]s_tSalinity at the cell top (ppt)
[in]s_bSalinity at the cell bottom (ppt)
[in]z_tAbsolute height of top of cell (m) (Boussinesq ????)
[in]z_bAbsolute height of bottom of cell (m)
[in]p_tAnomalous pressure of top of cell, relative to g*rho_ref*z_t (Pa)
[in]p_tgtTarget pressure at height z_out, relative to g*rho_ref*z_out (Pa)
[in]rho_refReference density with which calculation are anomalous to
[in]g_eGravitational acceleration (m/s2)
eosEquation of state structure
[out]p_bPressure at the bottom of the cell (Pa)
[out]z_outAbsolute depth at which anomalous pressure = p_tgt (m)

Definition at line 1235 of file MOM_EOS.F90.

References frac_dp_at_pos().

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 
Here is the call graph for this function:

◆ frac_dp_at_pos()

real function mom_eos::frac_dp_at_pos ( real, intent(in)  T_t,
real, intent(in)  T_b,
real, intent(in)  S_t,
real, intent(in)  S_b,
real, intent(in)  z_t,
real, intent(in)  z_b,
real, intent(in)  rho_ref,
real, intent(in)  G_e,
real, intent(in)  pos,
type(eos_type), pointer  EOS 
)
private

Returns change in anomalous pressure change from top to non-dimensional position pos between z_t and z_b.

Parameters
eosEquation of state structure

Definition at line 1304 of file MOM_EOS.F90.

References calculate_density_array().

Referenced by find_depth_of_pressure_in_cell().

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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ int_density_dz()

subroutine, public mom_eos::int_density_dz ( real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  T,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  S,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_b,
real, intent(in)  rho_ref,
real, intent(in)  rho_0,
real, intent(in)  G_e,
type(hor_index_type), intent(in)  HII,
type(hor_index_type), intent(in)  HIO,
type(eos_type), pointer  EOS,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out)  dpa,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out), optional  intz_dpa,
real, dimension(hio%isdb:hio%iedb,hio%jsd:hio%jed), intent(out), optional  intx_dpa,
real, dimension(hio%isd:hio%ied,hio%jsdb:hio%jedb), intent(out), optional  inty_dpa 
)

This subroutine calculates analytical and nearly-analytical integrals of pressure anomalies across layers, which are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. The one potentially dodgy assumtion here is that rho_0 is used both in the denominator of the accelerations, and in the pressure used to calculated density (the latter being -z*rho_0*G_e). These two uses could be separated if need be.

Parameters
[in]hiiOcean horizontal index structures for the input arrays
[in]hioOcean horizontal index structures for the output arrays
[in]tPotential temperature referenced to the surface (degC)
[in]sSalinity (PSU)
[in]z_tHeight at the top of the layer in m.
[in]z_bHeight at the bottom of the layer in m.
[in]rho_refA mean density, in kg m-3, that is subtracted out to reduce the magnitude of each of the integrals. (The pressure is calculated as p~=-z*rho_0*G_e.)
[in]rho_0A density, in kg m-3, that is used to calculate the pressure (as p~=-z*rho_0*G_e) used in the equation of state.
[in]g_eThe Earth's gravitational acceleration, in m s-2.
eosEquation of state structure
[out]dpaThe change in the pressure anomaly across the layer, in Pa.
[out]intz_dpaThe integral through the thickness of the layer of the pressure anomaly relative to the anomaly at the top of the layer, in Pa m.
[out]intx_dpaThe integral in x of the difference between the pressure anomaly at the top and bottom of the layer divided by the x grid spacing, in Pa.
[out]inty_dpaThe integral in y of the difference between the pressure anomaly at the top and bottom of the layer divided by the y grid spacing, in Pa.

Definition at line 392 of file MOM_EOS.F90.

References eos_linear, eos_wright, int_density_dz_generic(), and mom_error_handler::mom_error().

Referenced by mom_diagnostics::calculate_vertical_integrals().

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 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ int_density_dz_generic()

subroutine mom_eos::int_density_dz_generic ( real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  T,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  S,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_b,
real, intent(in)  rho_ref,
real, intent(in)  rho_0,
real, intent(in)  G_e,
type(hor_index_type), intent(in)  HII,
type(hor_index_type), intent(in)  HIO,
type(eos_type), pointer  EOS,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out)  dpa,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out), optional  intz_dpa,
real, dimension(hio%isdb:hio%iedb,hio%jsd:hio%jed), intent(out), optional  intx_dpa,
real, dimension(hio%isd:hio%ied,hio%jsdb:hio%jedb), intent(out), optional  inty_dpa 
)
private
Parameters
eosEquation of state structure

Definition at line 597 of file MOM_EOS.F90.

Referenced by int_density_dz().

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
Here is the caller graph for this function:

◆ int_density_dz_generic_cell()

subroutine mom_eos::int_density_dz_generic_cell ( real, dimension(2), intent(in)  T_t_arg,
real, dimension(2), intent(in)  T_b_arg,
real, dimension(2), intent(in)  S_t_arg,
real, dimension(2), intent(in)  S_b_arg,
real, dimension(2), intent(inout)  z_t_arg,
real, dimension(2), intent(inout)  z_b_arg,
real, dimension(2), intent(in)  depth,
real, intent(in)  rho_ref,
real, intent(in)  rho_0,
real, intent(in)  G_e,
type(eos_type), pointer  EOS,
real, dimension(2), intent(out)  dpa,
real, dimension(2), intent(out)  intz_dpa,
real, intent(out)  intx_dpa,
real, intent(out)  inty_dpa 
)
private
Parameters
eosEquation of state structure

Definition at line 733 of file MOM_EOS.F90.

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 

◆ int_density_dz_generic_plm()

subroutine, public mom_eos::int_density_dz_generic_plm ( real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  T_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  T_b,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  S_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  S_b,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_b,
real, intent(in)  rho_ref,
real, intent(in)  rho_0,
real, intent(in)  G_e,
real, intent(in)  H_subroundoff,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  bathyT,
type(hor_index_type), intent(in)  HII,
type(hor_index_type), intent(in)  HIO,
type(eos_type), pointer  EOS,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out)  dpa,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out), optional  intz_dpa,
real, dimension(hio%isdb:hio%iedb,hio%jsd:hio%jed), intent(out), optional  intx_dpa,
real, dimension(hio%isd:hio%ied,hio%jsdb:hio%jedb), intent(out), optional  inty_dpa,
logical, intent(in), optional  useMassWghtInterp 
)
Parameters
eosEquation of state structure

Definition at line 938 of file MOM_EOS.F90.

References calculate_density_array().

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 
Here is the call graph for this function:

◆ int_density_dz_generic_plm_analytic()

subroutine, public mom_eos::int_density_dz_generic_plm_analytic ( real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  T_t,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  T_b,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  S_t,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  S_b,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  z_t,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  z_b,
real, intent(in)  rho_ref,
real, intent(in)  rho_0,
real, intent(in)  G_e,
type(hor_index_type), intent(in)  HI,
type(eos_type), pointer  EOS,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(out)  dpa,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(out), optional  intz_dpa,
real, dimension(hi%isdb:hi%iedb,hi%jsd:hi%jed), intent(out), optional  intx_dpa,
real, dimension(hi%isd:hi%ied,hi%jsdb:hi%jedb), intent(out), optional  inty_dpa 
)
Parameters
eosEquation of state structure

Definition at line 1591 of file MOM_EOS.F90.

References compute_integral_bilinear().

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 
Here is the call graph for this function:

◆ int_density_dz_generic_ppm()

subroutine, public mom_eos::int_density_dz_generic_ppm ( real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  T,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  T_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  T_b,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  S,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  S_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  S_b,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_t,
real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in)  z_b,
real, intent(in)  rho_ref,
real, intent(in)  rho_0,
real, intent(in)  G_e,
type(hor_index_type), intent(in)  HII,
type(hor_index_type), intent(in)  HIO,
type(eos_type), pointer  EOS,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out)  dpa,
real, dimension(hio%isd:hio%ied,hio%jsd:hio%jed), intent(out), optional  intz_dpa,
real, dimension(hio%isdb:hio%iedb,hio%jsdb:hio%jedb), intent(out), optional  intx_dpa,
real, dimension(hio%isd:hio%ied,hio%jsdb:hio%jedb), intent(out), optional  inty_dpa 
)
Parameters
eosEquation of state structure

Definition at line 1339 of file MOM_EOS.F90.

References compute_integral_quadratic(), and mom_error_handler::mom_error().

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 
Here is the call graph for this function:

◆ int_spec_vol_dp_generic()

subroutine mom_eos::int_spec_vol_dp_generic ( real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  T,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  S,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  p_t,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  p_b,
real, intent(in)  alpha_ref,
type(hor_index_type), intent(in)  HI,
type(eos_type), pointer  EOS,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(out)  dza,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(out), optional  intp_dza,
real, dimension(hi%isdb:hi%iedb,hi%jsd:hi%jed), intent(out), optional  intx_dza,
real, dimension(hi%isd:hi%ied,hi%jsdb:hi%jedb), intent(out), optional  inty_dza,
integer, intent(in), optional  halo_size 
)
private
Parameters
eosEquation of state structure

Definition at line 2003 of file MOM_EOS.F90.

Referenced by int_specific_vol_dp().

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 
Here is the caller graph for this function:

◆ int_specific_vol_dp()

subroutine, public mom_eos::int_specific_vol_dp ( real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  T,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  S,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  p_t,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)  p_b,
real, intent(in)  alpha_ref,
type(hor_index_type), intent(in)  HI,
type(eos_type), pointer  EOS,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(out)  dza,
real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(out), optional  intp_dza,
real, dimension(hi%isdb:hi%iedb,hi%jsd:hi%jed), intent(out), optional  intx_dza,
real, dimension(hi%isd:hi%ied,hi%jsdb:hi%jedb), intent(out), optional  inty_dza,
integer, intent(in), optional  halo_size 
)

Calls the appropriate subroutine to alculate analytical and nearly-analytical integrals in pressure across layers of geopotential anomalies, which are required for calculating the finite-volume form pressure accelerations in a non-Boussinesq model. There are essentially no free assumptions, apart from the use of Bode's rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) that assumes that |eps| < .

Parameters
[in]hiThe horizontal index structure
[in]tPotential temperature referenced to the surface (degC)
[in]sSalinity (PSU)
[in]p_tPressure at the top of the layer in Pa.
[in]p_bPressure at the bottom of the layer in Pa.
[in]alpha_refA mean specific volume that is subtracted out to reduce the magnitude of each of the integrals, m3 kg-1. The calculation is mathematically identical with different values of alpha_ref, but this reduces the effects of roundoff.
eosEquation of state structure
[out]dzaThe change in the geopotential anomaly across the layer, in m2 s-2.
[out]intp_dzaThe integral in pressure through the layer of the geopotential anomaly relative to the anomaly at the bottom of the layer, in Pa m2 s-2.
[out]intx_dzaThe integral in x of the difference between the geopotential anomaly at the top and bottom of the layer divided by the x grid spacing, in m2 s-2.
[out]inty_dzaThe integral in y of the difference between the geopotential anomaly at the top and bottom of the layer divided by the y grid spacing, in m2 s-2.
[in]halo_sizeThe width of halo points on which to calculate dza.

Definition at line 333 of file MOM_EOS.F90.

References eos_linear, eos_wright, int_spec_vol_dp_generic(), mom_eos_linear::int_spec_vol_dp_linear(), mom_eos_wright::int_spec_vol_dp_wright(), and mom_error_handler::mom_error().

Referenced by mom_state_initialization::convert_thickness(), mom_interface_heights::find_eta_2d(), mom_interface_heights::find_eta_3d(), and mom_pressureforce_mont::pressureforce_mont_nonbouss().

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 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ query_compressible()

logical function, public mom_eos::query_compressible ( type(eos_type), pointer  EOS)

Returns true if the equation of state is compressible (i.e. has pressure dependence)

Parameters
eosEquation of state structure

Definition at line 449 of file MOM_EOS.F90.

References mom_error_handler::mom_error().

Referenced by mom_pressureforce_mont::pressureforce_mont_bouss(), and mom_pressureforce_mont::pressureforce_mont_nonbouss().

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
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ eos_default

character*(10), parameter mom_eos::eos_default = EOS_WRIGHT_STRING
private

Definition at line 86 of file MOM_EOS.F90.

Referenced by eos_init().

86 character*(10), parameter :: EOS_DEFAULT = EOS_WRIGHT_STRING

◆ eos_linear

integer, parameter mom_eos::eos_linear = 1

Definition at line 75 of file MOM_EOS.F90.

Referenced by calculate_compress(), calculate_density_array(), calculate_density_derivs(), calculate_density_scalar(), calculate_specific_vol_derivs(), eos_init(), eos_use_linear(), int_density_dz(), and int_specific_vol_dp().

75 integer, parameter :: eos_linear = 1

◆ eos_linear_string

character*(10), parameter mom_eos::eos_linear_string = "LINEAR"
private

Definition at line 81 of file MOM_EOS.F90.

Referenced by eos_init().

81 character*(10), parameter :: EOS_LINEAR_STRING = "LINEAR"

◆ eos_nemo

integer, parameter mom_eos::eos_nemo = 5
private

Definition at line 79 of file MOM_EOS.F90.

Referenced by calculate_compress(), calculate_density_array(), calculate_density_derivs(), calculate_density_scalar(), calculate_specific_vol_derivs(), convert_temp_salt_for_teos10(), and eos_init().

79 integer, parameter :: eos_nemo = 5

◆ eos_nemo_string

character*(10), parameter mom_eos::eos_nemo_string = "NEMO"
private

Definition at line 85 of file MOM_EOS.F90.

Referenced by eos_init().

85 character*(10), parameter :: EOS_NEMO_STRING = "NEMO"

◆ eos_teos10

integer, parameter mom_eos::eos_teos10 = 4
private

Definition at line 78 of file MOM_EOS.F90.

Referenced by calculate_compress(), calculate_density_array(), calculate_density_derivs(), calculate_density_scalar(), calculate_specific_vol_derivs(), convert_temp_salt_for_teos10(), and eos_init().

78 integer, parameter :: eos_teos10 = 4

◆ eos_teos10_string

character*(10), parameter mom_eos::eos_teos10_string = "TEOS10"
private

Definition at line 84 of file MOM_EOS.F90.

Referenced by eos_init().

84 character*(10), parameter :: EOS_TEOS10_STRING = "TEOS10"

◆ eos_unesco

integer, parameter mom_eos::eos_unesco = 2
private

Definition at line 76 of file MOM_EOS.F90.

Referenced by calculate_compress(), calculate_density_array(), calculate_density_derivs(), calculate_density_scalar(), calculate_specific_vol_derivs(), and eos_init().

76 integer, parameter :: eos_unesco = 2

◆ eos_unesco_string

character*(10), parameter mom_eos::eos_unesco_string = "UNESCO"
private

Definition at line 82 of file MOM_EOS.F90.

Referenced by eos_init().

82 character*(10), parameter :: EOS_UNESCO_STRING = "UNESCO"

◆ eos_wright

integer, parameter mom_eos::eos_wright = 3
private

Definition at line 77 of file MOM_EOS.F90.

Referenced by calculate_compress(), calculate_density_array(), calculate_density_derivs(), calculate_density_scalar(), calculate_specific_vol_derivs(), eos_init(), int_density_dz(), and int_specific_vol_dp().

77 integer, parameter :: eos_wright = 3

◆ eos_wright_string

character*(10), parameter mom_eos::eos_wright_string = "WRIGHT"
private

Definition at line 83 of file MOM_EOS.F90.

Referenced by eos_init().

83 character*(10), parameter :: EOS_WRIGHT_STRING = "WRIGHT"

◆ tfreeze_default

character*(10), parameter mom_eos::tfreeze_default = TFREEZE_LINEAR_STRING
private

Definition at line 94 of file MOM_EOS.F90.

Referenced by eos_init().

94 character*(10), parameter :: TFREEZE_DEFAULT = TFREEZE_LINEAR_STRING

◆ tfreeze_linear

integer, parameter mom_eos::tfreeze_linear = 1
private

Definition at line 88 of file MOM_EOS.F90.

Referenced by calculate_tfreeze_array(), calculate_tfreeze_scalar(), and eos_init().

88 integer, parameter :: tfreeze_linear = 1

◆ tfreeze_linear_string

character*(10), parameter mom_eos::tfreeze_linear_string = "LINEAR"
private

Definition at line 91 of file MOM_EOS.F90.

Referenced by eos_init().

91 character*(10), parameter :: TFREEZE_LINEAR_STRING = "LINEAR"

◆ tfreeze_millero

integer, parameter mom_eos::tfreeze_millero = 2
private

Definition at line 89 of file MOM_EOS.F90.

Referenced by calculate_tfreeze_array(), calculate_tfreeze_scalar(), and eos_init().

89 integer, parameter :: tfreeze_millero = 2

◆ tfreeze_millero_string

character*(10), parameter mom_eos::tfreeze_millero_string = "MILLERO_78"
private

Definition at line 92 of file MOM_EOS.F90.

Referenced by eos_init().

92 character*(10), parameter :: TFREEZE_MILLERO_STRING = "MILLERO_78"

◆ tfreeze_teos10

integer, parameter mom_eos::tfreeze_teos10 = 3
private

Definition at line 90 of file MOM_EOS.F90.

Referenced by calculate_tfreeze_array(), calculate_tfreeze_scalar(), and eos_init().

90 integer, parameter :: tfreeze_teos10 = 3

◆ tfreeze_teos10_string

character*(10), parameter mom_eos::tfreeze_teos10_string = "TEOS10"
private

Definition at line 93 of file MOM_EOS.F90.

Referenced by eos_init().

93 character*(10), parameter :: TFREEZE_TEOS10_STRING = "TEOS10"