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_compress
	Calculates the compressibility of water from T, S, and P. More...
interface	calculate_density
	Calculates density of sea water from T, S and P. More...
interface	calculate_density_derivs
	Calculate the derivatives of density with temperature and salinity from T, S, and P. More...
interface	calculate_density_second_derivs
	Calculates the second derivatives of density with various combinations of temperature, salinity, and pressure from T, S and P. More...
interface	calculate_spec_vol
	Calculates specific volume 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, rho_ref, scale)
	Calls the appropriate subroutine to calculate density of sea water for scalar inputs. If rho_ref is present, the anomaly with respect to rho_ref is returned. More...
subroutine	calculate_density_array (T, S, pressure, rho, start, npts, EOS, rho_ref, scale)
	Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs. If rho_ref is present, the anomaly with respect to rho_ref is returned. More...
subroutine	calculate_spec_vol_scalar (T, S, pressure, specvol, EOS, spv_ref, scale)
	Calls the appropriate subroutine to calculate specific volume of sea water for scalar inputs. More...
subroutine	calculate_spec_vol_array (T, S, pressure, specvol, start, npts, EOS, spv_ref, scale)
	Calls the appropriate subroutine to calculate the specific volume of sea water for 1-D array inputs. 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	calculate_density_derivs_array (T, S, pressure, drho_dT, drho_dS, start, npts, EOS, scale)
	Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs. More...
subroutine	calculate_density_derivs_scalar (T, S, pressure, drho_dT, drho_dS, EOS, scale)
	Calls the appropriate subroutines to calculate density derivatives by promoting a scalar to a one-element array. More...
subroutine	calculate_density_second_derivs_array (T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, start, npts, EOS, scale)
	Calls the appropriate subroutine to calculate density second derivatives for 1-D array inputs. More...
subroutine	calculate_density_second_derivs_scalar (T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, EOS, scale)
	Calls the appropriate subroutine to calculate density second derivatives for scalar nputs. More...
subroutine, public	calculate_specific_vol_derivs (T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS, scale)
	Calls the appropriate subroutine to calculate specific volume derivatives for an array. More...
subroutine	calculate_compress_array (T, S, pressure, rho, drho_dp, start, npts, EOS)
	Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs. More...
subroutine	calculate_compress_scalar (T, S, pressure, rho, drho_dp, EOS)
	Calculate density and compressibility for a scalar. This just promotes the scalar to an array with a singleton dimension and calls calculate_compress_array. 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, bathyP, dP_tiny, useMassWghtInterp)
	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, bathyT, dz_neglect, useMassWghtInterp)
	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. 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_manual_init (EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
	Manually initialized an EOS type (intended for unit testing of routines which need a specific EOS) 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, public	int_density_dz_generic (T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp)
	This subroutine calculates (by numerical quadrature) integrals of pressure anomalies across layers, which are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. More...
subroutine, public	int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp)
	Compute pressure gradient force integrals by quadrature for the case where T and S are linear profiles. More...
subroutine, public	find_depth_of_pressure_in_cell (T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, rho_ref, G_e, EOS, P_b, z_out, z_tol)
	Find the depth at which the reconstructed pressure matches P_tgt. 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)
	Compute pressure gradient force integrals for the case where T and S are parabolic profiles. More...
subroutine	compute_integral_quadratic (x, y, f, integral)
	Compute the integral of the quadratic function. More...
subroutine	evaluate_shape_bilinear (xi, eta, phi, dphidxi, dphideta)
	Evaluation of the four bilinear shape fn and their gradients at (xi,eta) More...
subroutine	evaluate_shape_quadratic (xi, eta, phi, dphidxi, dphideta)
	Evaluation of the nine quadratic shape fn weights and their gradients at (xi,eta) More...
subroutine, public	int_spec_vol_dp_generic (T, S, p_t, p_b, alpha_ref, HI, EOS, dza, intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, useMassWghtInterp)
	This subroutine calculates integrals of specific volume anomalies in pressure across layers, 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 quadrature to do the integrals. More...
subroutine, public	int_spec_vol_dp_generic_plm (T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, dP_neglect, bathyP, HI, EOS, dza, intp_dza, intx_dza, inty_dza, useMassWghtInterp)
	This subroutine calculates integrals of specific volume anomalies in pressure across layers, 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 quadrature to do the integrals. More...
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...
subroutine, public	extract_member_eos (EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, Rho_T0_S0, drho_dT, dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp)
	Extractor routine for the EOS type if the members need to be accessed outside this module. More...

Variables
integer, parameter, public	eos_linear = 1
	A named integer specifying an equation of state. More...
integer, parameter, public	eos_unesco = 2
	A named integer specifying an equation of state. More...
integer, parameter, public	eos_wright = 3
	A named integer specifying an equation of state. More...
integer, parameter, public	eos_teos10 = 4
	A named integer specifying an equation of state. More...
integer, parameter, public	eos_nemo = 5
	A named integer specifying an equation of state. More...
character *(10), parameter	eos_linear_string = "LINEAR"
	A string for specifying the equation of state. More...
character *(10), parameter	eos_unesco_string = "UNESCO"
	A string for specifying the equation of state. More...
character *(10), parameter	eos_wright_string = "WRIGHT"
	A string for specifying the equation of state. More...
character *(10), parameter	eos_teos10_string = "TEOS10"
	A string for specifying the equation of state. More...
character *(10), parameter	eos_nemo_string = "NEMO"
	A string for specifying the equation of state. More...
character *(10), parameter	eos_default = EOS_WRIGHT_STRING
	The default equation of state. More...
integer, parameter	tfreeze_linear = 1
	A named integer specifying a freezing point expression. More...
integer, parameter	tfreeze_millero = 2
	A named integer specifying a freezing point expression. More...
integer, parameter	tfreeze_teos10 = 3
	A named integer specifying a freezing point expression. More...
character *(10), parameter	tfreeze_linear_string = "LINEAR"
	A string for specifying the freezing point expression. More...
character *(10), parameter	tfreeze_millero_string = "MILLERO_78"
	A string for specifying freezing point expression. More...
character *(10), parameter	tfreeze_teos10_string = "TEOS10"
	A string for specifying the freezing point expression. More...
character *(10), parameter	tfreeze_default = TFREEZE_LINEAR_STRING
	The default freezing point expression. More...

Function/Subroutine Documentation

calculate_compress_array()

subroutine mom_eos::calculate_compress_array ( 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  )

private

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

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [PSU]
[in]	pressure	Pressure [Pa]
[out]	rho	In situ density [kg m-3].
[out]	drho_dp	The partial derivative of density with pressure (also the inverse of the square of sound speed) [s2 m-2].
[in]	start	Starting index within the array
[in]	npts	The number of values to calculate
	eos	Equation of state structure

Definition at line 603 of file MOM_EOS.F90.

  real, dimension(:), intent(in)  :: T        !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S        !< Salinity [PSU]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: rho      !< In situ density [kg m-3].
  real, dimension(:), intent(out) :: drho_dp   !< The partial derivative of density with pressure
                                     !! (also the inverse of the square of sound speed) [s2 m-2].
  integer,            intent(in)  :: start    !< Starting index within the array
  integer,            intent(in)  :: npts     !< The number of values to calculate
  type(EOS_type),     pointer     :: EOS      !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_compress called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_compress_linear(t, s, pressure, rho, drho_dp, start, npts, &
                                     eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
    case (eos_unesco)
      call calculate_compress_unesco(t, s, pressure, rho, drho_dp, start, npts)
    case (eos_wright)
      call calculate_compress_wright(t, s, pressure, rho, drho_dp, start, npts)
    case (eos_teos10)
      call calculate_compress_teos10(t, s, pressure, rho, drho_dp, start, npts)
    case (eos_nemo)
      call calculate_compress_nemo(t, s, pressure, rho, drho_dp, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_compress: EOS%form_of_EOS is not valid.")
  end select
 

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

Referenced by calculate_compress_scalar().

Here is the call graph for this function:

Here is the caller graph for this function:

calculate_compress_scalar()

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

private

Calculate density and compressibility for a scalar. This just promotes the scalar to an array with a singleton dimension and calls calculate_compress_array.

Parameters

[in]	t	Potential temperature referenced to the surface (degC)
[in]	s	Salinity (PSU)
[in]	pressure	Pressure (Pa)
[out]	rho	In situ density in kg m-3.
[out]	drho_dp	The partial derivative of density with pressure (also the inverse of the square of sound speed) in s2 m-2.
	eos	Equation of state structure

Definition at line 638 of file MOM_EOS.F90.

  real, intent(in)        :: T        !< Potential temperature referenced to the surface (degC)
  real, intent(in)        :: S        !< Salinity (PSU)
  real, intent(in)        :: pressure !< Pressure (Pa)
  real, intent(out)       :: rho      !< In situ density in kg m-3.
  real, intent(out)       :: drho_dp  !< The partial derivative of density with pressure
                                      !! (also the inverse of the square of sound speed) in s2 m-2.
  type(EOS_type), pointer :: EOS      !< Equation of state structure
 
  real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_compress called with an unassociated EOS_type EOS.")
  ta(1) = t ; sa(1) = s; pa(1) = pressure
 
  call calculate_compress_array(ta, sa, pa, rhoa, drho_dpa, 1, 1, eos)
  rho = rhoa(1) ; drho_dp = drho_dpa(1)
 

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

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, real, intent(in), optional  rho_ref, real, intent(in), optional  scale  )

private

Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs. If rho_ref is present, the anomaly with respect to rho_ref is returned.

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	rho	Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3]
[in]	start	Start index for computation
[in]	npts	Number of point to compute
	eos	Equation of state structure
[in]	rho_ref	A reference density [kg m-3].
[in]	scale	A multiplicative factor by which to scale density from kg m-3 to the desired units [R m3 kg-1]

Definition at line 177 of file MOM_EOS.F90.

  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3]
  integer,            intent(in)  :: start !< Start index for computation
  integer,            intent(in)  :: npts !< Number of point to compute
  type(EOS_type),     pointer     :: EOS !< Equation of state structure
  real,     optional, intent(in)  :: rho_ref  !< A reference density [kg m-3].
  real,     optional, intent(in)  :: scale !< A multiplicative factor by which to scale density
                                           !! from kg m-3 to the desired units [R m3 kg-1]
  integer :: j
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_array called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_linear(t, s, pressure, rho, start, npts, &
                                      eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
    case (eos_unesco)
      call calculate_density_unesco(t, s, pressure, rho, start, npts, rho_ref)
    case (eos_wright)
      call calculate_density_wright(t, s, pressure, rho, start, npts, rho_ref)
    case (eos_teos10)
      call calculate_density_teos10(t, s, pressure, rho, start, npts, rho_ref)
    case (eos_nemo)
      call calculate_density_nemo  (t, s, pressure, rho, start, npts, rho_ref)
    case default
      call mom_error(fatal, &
           "calculate_density_array: EOS%form_of_EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then
    do j=start,start+npts-1 ; rho(j) = scale * rho(j) ; enddo
  endif ; endif
 

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().

Here is the call graph for this function:

Here is the caller graph for this function:

calculate_density_derivs_array()

subroutine mom_eos::calculate_density_derivs_array ( 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, real, intent(in), optional  scale  )

private

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

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	drho_dt	The partial derivative of density with potential temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1].
[out]	drho_ds	The partial derivative of density with salinity, in [kg m-3 ppt-1] or [R degC-1 ~> kg m-3 ppt-1].
[in]	start	Starting index within the array
[in]	npts	The number of values to calculate
	eos	Equation of state structure
[in]	scale	A multiplicative factor by which to scale density from kg m-3 to the desired units [R m3 kg-1]

Definition at line 367 of file MOM_EOS.F90.

  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: drho_dT !< The partial derivative of density with potential
                                             !! temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1].
  real, dimension(:), intent(out) :: drho_dS !< The partial derivative of density with salinity,
                                             !! in [kg m-3 ppt-1] or [R degC-1 ~> kg m-3 ppt-1].
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(EOS_type),     pointer     :: EOS !< Equation of state structure
  real,     optional, intent(in)  :: scale !< A multiplicative factor by which to scale density
                                           !! from kg m-3 to the desired units [R m3 kg-1]
  integer :: j
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_derivs_linear(t, s, pressure, drho_dt, drho_ds, eos%Rho_T0_S0, &
                                           eos%dRho_dT, eos%dRho_dS, start, npts)
    case (eos_unesco)
      call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
    case (eos_wright)
      call calculate_density_derivs_wright(t, s, pressure, drho_dt, drho_ds, start, npts)
    case (eos_teos10)
      call calculate_density_derivs_teos10(t, s, pressure, drho_dt, drho_ds, start, npts)
    case (eos_nemo)
      call calculate_density_derivs_nemo(t, s, pressure, drho_dt, drho_ds, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs_array: EOS%form_of_EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
    drho_dt(j) = scale * drho_dt(j)
    drho_ds(j) = scale * drho_ds(j)
  enddo ; endif ; endif
 

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

Here is the call graph for this function:

calculate_density_derivs_scalar()

subroutine mom_eos::calculate_density_derivs_scalar ( real, intent(in)  T, real, intent(in)  S, real, intent(in)  pressure, real, intent(out)  drho_dT, real, intent(out)  drho_dS, type(eos_type), pointer  EOS, real, intent(in), optional  scale  )

private

Calls the appropriate subroutines to calculate density derivatives by promoting a scalar to a one-element array.

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	drho_dt	The partial derivative of density with potential temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1]
[out]	drho_ds	The partial derivative of density with salinity, in [kg m-3 ppt-1] or [R ppt-1 ~> kg m-3 ppt-1].
	eos	Equation of state structure
[in]	scale	A multiplicative factor by which to scale density from kg m-3 to the desired units [R m3 kg-1]

Definition at line 411 of file MOM_EOS.F90.

  real,           intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: drho_dT !< The partial derivative of density with potential
                                         !! temperature [kg m-3 degC-1] or [R degC-1 ~> kg m-3 degC-1]
  real,           intent(out) :: drho_dS !< The partial derivative of density with salinity,
                                         !! in [kg m-3 ppt-1] or [R ppt-1 ~> kg m-3 ppt-1].
  type(EOS_type), pointer     :: EOS !< Equation of state structure
  real, optional, intent(in)  :: scale   !< A multiplicative factor by which to scale density
                                         !! from kg m-3 to the desired units [R m3 kg-1]
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_derivs_linear(t, s, pressure, drho_dt, drho_ds, &
                                           eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
    case (eos_wright)
      call calculate_density_derivs_wright(t, s, pressure, drho_dt, drho_ds)
    case (eos_teos10)
      call calculate_density_derivs_teos10(t, s, pressure, drho_dt, drho_ds)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs_scalar: EOS%form_of_EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then
    drho_dt = scale * drho_dt
    drho_ds = scale * drho_ds
  endif ; endif
 

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

Here is the call 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, real, intent(in), optional  rho_ref, real, intent(in), optional  scale  )

private

Calls the appropriate subroutine to calculate density of sea water for scalar inputs. If rho_ref is present, the anomaly with respect to rho_ref is returned.

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	rho	Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3]
	eos	Equation of state structure
[in]	rho_ref	A reference density [kg m-3].
[in]	scale	A multiplicative factor by which to scale density from kg m-3 to the desired units [R m3 kg-1]

Definition at line 139 of file MOM_EOS.F90.

  real,           intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3]
  type(EOS_type), pointer     :: EOS !< Equation of state structure
  real, optional, intent(in)  :: rho_ref  !< A reference density [kg m-3].
  real, optional, intent(in)  :: scale !< A multiplicative factor by which to scale density
                                       !! from kg m-3 to the desired units [R m3 kg-1]
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_linear(t, s, pressure, rho, &
                                      eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, rho_ref)
    case (eos_unesco)
      call calculate_density_unesco(t, s, pressure, rho, rho_ref)
    case (eos_wright)
      call calculate_density_wright(t, s, pressure, rho, rho_ref)
    case (eos_teos10)
      call calculate_density_teos10(t, s, pressure, rho, rho_ref)
    case (eos_nemo)
      call calculate_density_nemo(t, s, pressure, rho, rho_ref)
    case default
      call mom_error(fatal, &
           "calculate_density_scalar: EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then
    rho = scale * rho
  endif ; endif
 

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

Here is the call graph for this function:

calculate_density_second_derivs_array()

subroutine mom_eos::calculate_density_second_derivs_array ( real, dimension(:), intent(in)  T, real, dimension(:), intent(in)  S, real, dimension(:), intent(in)  pressure, real, dimension(:), intent(out)  drho_dS_dS, real, dimension(:), intent(out)  drho_dS_dT, real, dimension(:), intent(out)  drho_dT_dT, real, dimension(:), intent(out)  drho_dS_dP, real, dimension(:), intent(out)  drho_dT_dP, integer, intent(in)  start, integer, intent(in)  npts, type(eos_type), pointer  EOS, real, intent(in), optional  scale  )

private

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

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	drho_ds_ds	Partial derivative of beta with respect to S [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2]
[out]	drho_ds_dt	Partial derivative of beta with respect to T [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1]
[out]	drho_dt_dt	Partial derivative of alpha with respect to T [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2]
[out]	drho_ds_dp	Partial derivative of beta with respect to pressure [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1]
[out]	drho_dt_dp	Partial derivative of alpha with respect to pressure [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1]
[in]	start	Starting index within the array
[in]	npts	The number of values to calculate
	eos	Equation of state structure
[in]	scale	A multiplicative factor by which to scale density from kg m-3 to the desired units [R m3 kg-1]

Definition at line 448 of file MOM_EOS.F90.

  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure   !< Pressure [Pa]
  real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S
                                                !!  [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2]
  real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with respect to T
                                                !! [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1]
  real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T
                                                !! [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2]
  real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure
                                                !! [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1]
  real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure
                                                !! [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1]
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(EOS_type),     pointer     :: EOS !< Equation of state structure
  real,     optional, intent(in)  :: scale !< A multiplicative factor by which to scale density
                                           !! from kg m-3 to the desired units [R m3 kg-1]
  integer :: j
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_second_derivs_linear(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
    case (eos_wright)
      call calculate_density_second_derivs_wright(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
    case (eos_teos10)
      call calculate_density_second_derivs_teos10(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs: EOS%form_of_EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
    drho_ds_ds(j) = scale * drho_ds_ds(j)
    drho_ds_dt(j) = scale * drho_ds_dt(j)
    drho_dt_dt(j) = scale * drho_dt_dt(j)
    drho_ds_dp(j) = scale * drho_ds_dp(j)
    drho_dt_dp(j) = scale * drho_dt_dp(j)
  enddo ; endif ; endif
 

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

Here is the call graph for this function:

calculate_density_second_derivs_scalar()

subroutine mom_eos::calculate_density_second_derivs_scalar ( real, intent(in)  T, real, intent(in)  S, real, intent(in)  pressure, real, intent(out)  drho_dS_dS, real, intent(out)  drho_dS_dT, real, intent(out)  drho_dT_dT, real, intent(out)  drho_dS_dP, real, intent(out)  drho_dT_dP, type(eos_type), pointer  EOS, real, intent(in), optional  scale  )

private

Calls the appropriate subroutine to calculate density second derivatives for scalar nputs.

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	drho_ds_ds	Partial derivative of beta with respect to S [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2]
[out]	drho_ds_dt	Partial derivative of beta with respect to T [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1]
[out]	drho_dt_dt	Partial derivative of alpha with respect to T [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2]
[out]	drho_ds_dp	Partial derivative of beta with respect to pressure [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1]
[out]	drho_dt_dp	Partial derivative of alpha with respect to pressure [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1]
	eos	Equation of state structure
[in]	scale	A multiplicative factor by which to scale density from kg m-3 to the desired units [R m3 kg-1]

Definition at line 499 of file MOM_EOS.F90.

  real, intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, intent(in)  :: S !< Salinity [ppt]
  real, intent(in)  :: pressure   !< Pressure [Pa]
  real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S
                                  !! [kg m-3 ppt-2] or [R ppt-2 ~> kg m-3 ppt-2]
  real, intent(out) :: drho_dS_dT !< Partial derivative of beta with respect to T
                                  !! [kg m-3 ppt-1 degC-1] or [R ppt-1 degC-1 ~> kg m-3 ppt-1 degC-1]
  real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T
                                  !! [kg m-3 degC-2] or [R degC-2 ~> kg m-3 degC-2]
  real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure
                                  !! [kg m-3 ppt-1 Pa-1] or [R ppt-1 Pa-1 ~> kg m-3 ppt-1 Pa-1]
  real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure
                                  !! [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1]
  type(EOS_type), pointer    :: EOS !< Equation of state structure
  real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density
                                  !! from kg m-3 to the desired units [R m3 kg-1]
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_density_second_derivs_linear(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp)
    case (eos_wright)
      call calculate_density_second_derivs_wright(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp)
    case (eos_teos10)
      call calculate_density_second_derivs_teos10(t, s, pressure, drho_ds_ds, drho_ds_dt, &
                                                  drho_dt_dt, drho_ds_dp, drho_dt_dp)
    case default
      call mom_error(fatal, &
           "calculate_density_derivs: EOS%form_of_EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then
    drho_ds_ds = scale * drho_ds_ds
    drho_ds_dt = scale * drho_ds_dt
    drho_dt_dt = scale * drho_dt_dt
    drho_ds_dp = scale * drho_ds_dp
    drho_dt_dp = scale * drho_dt_dp
  endif ; endif
 

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

Here is the call graph for this function:

calculate_spec_vol_array()

subroutine mom_eos::calculate_spec_vol_array ( real, dimension(:), intent(in)  T, real, dimension(:), intent(in)  S, real, dimension(:), intent(in)  pressure, real, dimension(:), intent(out)  specvol, integer, intent(in)  start, integer, intent(in)  npts, type(eos_type), pointer  EOS, real, intent(in), optional  spv_ref, real, intent(in), optional  scale  )

private

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

Parameters

[in]	t	potential temperature relative to the surface [degC].
[in]	s	salinity [ppt].
[in]	pressure	pressure [Pa].
[out]	specvol	in situ specific volume [kg m-3] or [R-1 ~> m3 kg-1].
[in]	start	the starting point in the arrays.
[in]	npts	the number of values to calculate.
	eos	Equation of state structure
[in]	spv_ref	A reference specific volume [m3 kg-1].
[in]	scale	A multiplicative factor by which to scale specific volume from m3 kg-1 to the desired units [kg m-3 R-1]

Definition at line 264 of file MOM_EOS.F90.

  real, dimension(:), intent(in)  :: T        !< potential temperature relative to the surface
                                              !! [degC].
  real, dimension(:), intent(in)  :: S        !< salinity [ppt].
  real, dimension(:), intent(in)  :: pressure !< pressure [Pa].
  real, dimension(:), intent(out) :: specvol  !< in situ specific volume [kg m-3] or [R-1 ~> m3 kg-1].
  integer,            intent(in)  :: start    !< the starting point in the arrays.
  integer,            intent(in)  :: npts     !< the number of values to calculate.
  type(EOS_type),     pointer     :: EOS      !< Equation of state structure
  real,     optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
  real,     optional, intent(in)  :: scale    !< A multiplicative factor by which to scale specific volume
                                              !! from m3 kg-1 to the desired units [kg m-3 R-1]
 
  real, dimension(size(specvol)) :: rho
  integer :: j
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_spec_vol_array called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_spec_vol_linear(t, s, pressure, specvol, start, npts, &
               eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
    case (eos_unesco)
      call calculate_spec_vol_unesco(t, s, pressure, specvol, start, npts, spv_ref)
    case (eos_wright)
      call calculate_spec_vol_wright(t, s, pressure, specvol, start, npts, spv_ref)
    case (eos_teos10)
      call calculate_spec_vol_teos10(t, s, pressure, specvol, start, npts, spv_ref)
    case (eos_nemo)
      call calculate_density_nemo  (t, s, pressure, rho, start, npts)
      if (present(spv_ref)) then
        specvol(:) = 1.0 / rho(:) - spv_ref
      else
        specvol(:) = 1.0 / rho(:)
      endif
    case default
      call mom_error(fatal, &
           "calculate_spec_vol_array: EOS%form_of_EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
    specvol(j) = scale * specvol(j)
  enddo ;  endif ; endif
 

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

Here is the call graph for this function:

calculate_spec_vol_scalar()

subroutine mom_eos::calculate_spec_vol_scalar ( real, intent(in)  T, real, intent(in)  S, real, intent(in)  pressure, real, intent(out)  specvol, type(eos_type), pointer  EOS, real, intent(in), optional  spv_ref, real, intent(in), optional  scale  )

private

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

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	specvol	In situ? specific volume [m3 kg-1] or [R-1 ~> m3 kg-1]
	eos	Equation of state structure
[in]	spv_ref	A reference specific volume [m3 kg-1].
[in]	scale	A multiplicative factor by which to scale specific volume from m3 kg-1 to the desired units [kg m-3 R-1]

Definition at line 218 of file MOM_EOS.F90.

  real,           intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: specvol  !< In situ? specific volume [m3 kg-1] or [R-1 ~> m3 kg-1]
  type(EOS_type), pointer     :: EOS      !< Equation of state structure
  real, optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
  real, optional, intent(in)  :: scale    !< A multiplicative factor by which to scale specific volume
                                          !! from m3 kg-1 to the desired units [kg m-3 R-1]
 
  real :: rho
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_spec_vol_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_spec_vol_linear(t, s, pressure, specvol, &
               eos%rho_T0_S0, eos%drho_dT, eos%drho_dS, spv_ref)
    case (eos_unesco)
      call calculate_spec_vol_unesco(t, s, pressure, specvol, spv_ref)
    case (eos_wright)
      call calculate_spec_vol_wright(t, s, pressure, specvol, spv_ref)
    case (eos_teos10)
      call calculate_spec_vol_teos10(t, s, pressure, specvol, spv_ref)
    case (eos_nemo)
      call calculate_density_nemo(t, s, pressure, rho)
      if (present(spv_ref)) then
        specvol = 1.0 / rho - spv_ref
      else
        specvol = 1.0 / rho
      endif
    case default
      call mom_error(fatal, &
           "calculate_spec_vol_scalar: EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then
    specvol = scale * specvol
  endif ; endif
 

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

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,
		real, intent(in), optional	scale
	)

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

Parameters

[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	dsv_dt	The partial derivative of specific volume with potential temperature [m3 kg-1 degC-1] or [R-1 degC-1 ~> m3 kg-1 degC-1]
[out]	dsv_ds	The partial derivative of specific volume with salinity [m3 kg-1 ppt-1] or [R-1 ppt-1 ~> m3 kg-1 ppt-1].
[in]	start	Starting index within the array
[in]	npts	The number of values to calculate
	eos	Equation of state structure
[in]	scale	A multiplicative factor by which to scale specific volume from m3 kg-1 to the desired units [kg m-3 R-1]

Definition at line 546 of file MOM_EOS.F90.

  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: dSV_dT !< The partial derivative of specific volume with potential
                                            !! temperature [m3 kg-1 degC-1] or [R-1 degC-1 ~> m3 kg-1 degC-1]
  real, dimension(:), intent(out) :: dSV_dS !< The partial derivative of specific volume with salinity
                                            !! [m3 kg-1 ppt-1] or [R-1 ppt-1 ~> m3 kg-1 ppt-1].
  integer,            intent(in)  :: start  !< Starting index within the array
  integer,            intent(in)  :: npts   !< The number of values to calculate
  type(EOS_type),     pointer     :: EOS    !< Equation of state structure
  real, optional,     intent(in)  :: scale  !< A multiplicative factor by which to scale specific volume
                                            !! from m3 kg-1 to the desired units [kg m-3 R-1]
 
  ! Local variables
  real, dimension(size(T)) :: dRho_dT, dRho_dS, rho
  integer :: j
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_density_derivs called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_EOS)
    case (eos_linear)
      call calculate_specvol_derivs_linear(t, s, pressure, dsv_dt, dsv_ds, start, &
                                           npts, eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS)
    case (eos_unesco)
      call calculate_density_unesco(t, s, pressure, rho, start, npts)
      call calculate_density_derivs_unesco(t, s, pressure, drho_dt, drho_ds, start, npts)
      do j=start,start+npts-1
        dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
        dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
      enddo
    case (eos_wright)
      call calculate_specvol_derivs_wright(t, s, pressure, dsv_dt, dsv_ds, start, npts)
    case (eos_teos10)
      call calculate_specvol_derivs_teos10(t, s, pressure, dsv_dt, dsv_ds, start, npts)
    case (eos_nemo)
      call calculate_density_nemo(t, s, pressure, rho, start, npts)
      call calculate_density_derivs_nemo(t, s, pressure, drho_dt, drho_ds, start, npts)
      do j=start,start+npts-1
        dsv_dt(j) = -drho_dt(j)/(rho(j)**2)
        dsv_ds(j) = -drho_ds(j)/(rho(j)**2)
      enddo
    case default
      call mom_error(fatal, &
           "calculate_density_derivs: EOS%form_of_EOS is not valid.")
  end select
 
  if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1
    dsv_dt(j) = scale * dsv_dt(j)
    dsv_ds(j) = scale * dsv_ds(j)
  enddo ;  endif ; endif
 
 

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().

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]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	t_fr	Freezing point potential temperature referenced to the surface [degC]
[in]	start	Starting index within the array
[in]	npts	The number of values to calculate
	eos	Equation of state structure

Definition at line 339 of file MOM_EOS.F90.

  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
  real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
  real, dimension(:), intent(out) :: T_fr !< Freezing point potential temperature referenced
                                          !! to the surface [degC]
  integer,            intent(in)  :: start !< Starting index within the array
  integer,            intent(in)  :: npts !< The number of values to calculate
  type(EOS_type),     pointer     :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_TFreeze)
    case (tfreeze_linear)
      call calculate_tfreeze_linear(s, pressure, t_fr, start, npts, &
                                    eos%TFr_S0_P0, eos%dTFr_dS, eos%dTFr_dp)
    case (tfreeze_millero)
      call calculate_tfreeze_millero(s, pressure, t_fr, start, npts)
    case (tfreeze_teos10)
      call calculate_tfreeze_teos10(s, pressure, t_fr, start, npts)
    case default
      call mom_error(fatal, &
           "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
  end select
 

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

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]	s	Salinity [ppt]
[in]	pressure	Pressure [Pa]
[out]	t_fr	Freezing point potential temperature referenced to the surface [degC]
	eos	Equation of state structure

Definition at line 313 of file MOM_EOS.F90.

  real,           intent(in)  :: S !< Salinity [ppt]
  real,           intent(in)  :: pressure !< Pressure [Pa]
  real,           intent(out) :: T_fr !< Freezing point potential temperature referenced
                                      !! to the surface [degC]
  type(EOS_type), pointer     :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "calculate_TFreeze_scalar called with an unassociated EOS_type EOS.")
 
  select case (eos%form_of_TFreeze)
    case (tfreeze_linear)
      call calculate_tfreeze_linear(s, pressure, t_fr, eos%TFr_S0_P0, &
                                    eos%dTFr_dS, eos%dTFr_dp)
    case (tfreeze_millero)
      call calculate_tfreeze_millero(s, pressure, t_fr)
    case (tfreeze_teos10)
      call calculate_tfreeze_teos10(s, pressure, t_fr)
    case default
      call mom_error(fatal, &
           "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
  end select
 

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

Here is the call 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

Compute the integral of the quadratic function.

Parameters

[in]	x	The x-position of the corners
[in]	y	The y-position of the corners
[in]	f	The function at the quadrature points
[out]	integral	The returned integral

Definition at line 1834 of file MOM_EOS.F90.

  real, dimension(4), intent(in)  :: x  !< The x-position of the corners
  real, dimension(4), intent(in)  :: y  !< The y-position of the corners
  real, dimension(9), intent(in)  :: f  !< The function at the quadrature points
  real,               intent(out) :: integral !< The returned integral
 
  ! Local variables
  integer               :: i, k
  real, dimension(9)    :: weight, xi, eta          ! integration points
  real                  :: f_k
  real                  :: dxdxi, dxdeta
  real                  :: dydxi, dydeta
  real, dimension(4)    :: phiiso, dphiisodxi, dphiisodeta
  real, dimension(9)    :: phi, dphidxi, dphideta
  real                  :: jacobian_k
  real                  :: t
 
  ! Quadrature rule (4 points)
  !weight(:) = 1.0
  !xi(1) = - sqrt(3.0) / 3.0
  !xi(2) = sqrt(3.0) / 3.0
  !xi(3) = sqrt(3.0) / 3.0
  !xi(4) = - sqrt(3.0) / 3.0
  !eta(1) = - sqrt(3.0) / 3.0
  !eta(2) = - sqrt(3.0) / 3.0
  !eta(3) = sqrt(3.0) / 3.0
  !eta(4) = sqrt(3.0) / 3.0
 
  ! Quadrature rule (9 points)
  t = sqrt(3.0/5.0)
  weight(1) = 25.0/81.0 ; xi(1) = -t ; eta(1) = t
  weight(2) = 40.0/81.0 ; xi(2) = .0 ; eta(2) = t
  weight(3) = 25.0/81.0 ; xi(3) =  t ; eta(3) = t
  weight(4) = 40.0/81.0 ; xi(4) = -t ; eta(4) = .0
  weight(5) = 64.0/81.0 ; xi(5) = .0 ; eta(5) = .0
  weight(6) = 40.0/81.0 ; xi(6) =  t ; eta(6) = .0
  weight(7) = 25.0/81.0 ; xi(7) = -t ; eta(7) = -t
  weight(8) = 40.0/81.0 ; xi(8) = .0 ; eta(8) = -t
  weight(9) = 25.0/81.0 ; xi(9) =  t ; eta(9) = -t
 
  integral = 0.0
 
  ! Integration loop
  do k = 1,9
 
    ! Evaluate shape functions and gradients for isomorphism
    call evaluate_shape_bilinear( xi(k), eta(k), phiiso, &
                                  dphiisodxi, dphiisodeta )
 
    ! Determine gradient of global coordinate at integration point
    dxdxi  = 0.0
    dxdeta = 0.0
    dydxi  = 0.0
    dydeta = 0.0
 
    do i = 1,4
      dxdxi  = dxdxi  + x(i) * dphiisodxi(i)
      dxdeta = dxdeta + x(i) * dphiisodeta(i)
      dydxi  = dydxi  + y(i) * dphiisodxi(i)
      dydeta = dydeta + y(i) * dphiisodeta(i)
    enddo
 
    ! Evaluate Jacobian at integration point
    jacobian_k = dxdxi*dydeta - dydxi*dxdeta
 
    ! Evaluate shape functions for interpolation
    call evaluate_shape_quadratic( xi(k), eta(k), phi, dphidxi, dphideta )
 
    ! Evaluate function at integration point
    f_k = 0.0
    do i = 1,9
      f_k = f_k + f(i) * phi(i)
    enddo
 
    integral = integral + weight(k) * f_k * jacobian_k
 
  enddo ! end integration loop
 

References evaluate_shape_bilinear(), and evaluate_shape_quadratic().

Referenced by int_density_dz_generic_ppm().

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]	g	The ocean's grid structure
[in,out]	t	Potential temperature referenced to the surface [degC]
[in,out]	s	Salinity [ppt]
[in]	press	Pressure at the top of the layer [Pa].
	eos	Equation of state structure
[in]	mask_z	3d mask regulating which points to convert.
[in]	kd	The number of layers to work on

Definition at line 2438 of file MOM_EOS.F90.

  use mom_grid, only : ocean_grid_type
 
  type(ocean_grid_type), intent(in)    :: G   !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                         intent(inout) :: T   !< Potential temperature referenced to the surface [degC]
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                         intent(inout) :: S   !< Salinity [ppt]
  real, dimension(:),    intent(in)    :: press !< Pressure at the top of the layer [Pa].
  type(EOS_type),        pointer       :: EOS !< Equation of state structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                         intent(in)    :: mask_z !< 3d mask regulating which points to convert.
  integer,               intent(in)    :: kd  !< The number of layers to work on
 
  integer :: i,j,k
  real :: gsw_sr_from_sp, gsw_ct_from_pt, gsw_sa_from_sp
  real :: p
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "convert_temp_salt_to_TEOS10 called with an unassociated EOS_type EOS.")
 
  if ((eos%form_of_EOS /= eos_teos10) .and. (eos%form_of_EOS /= eos_nemo)) return
 
  do k=1,kd ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if (mask_z(i,j,k) >= 1.0) then
     s(i,j,k) = gsw_sr_from_sp(s(i,j,k))
!     p=press(k)/10000. !convert pascal to dbar
!     S(i,j,k) = gsw_sa_from_sp(S(i,j,k),p,G%geoLonT(i,j),G%geoLatT(i,j))
     t(i,j,k) = gsw_ct_from_pt(s(i,j,k),t(i,j,k))
    endif
  enddo ; enddo ; enddo

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

Here is the call graph for this function:

eos_allocate()

subroutine, public mom_eos::eos_allocate

(

type(eos_type), pointer

EOS

)

Allocates EOS_type.

Parameters

eos	Equation of state structure

Definition at line 940 of file MOM_EOS.F90.

  type(EOS_type), pointer :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) allocate(eos)

Referenced by eos_init().

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

eos	Equation of state structure

Definition at line 947 of file MOM_EOS.F90.

  type(EOS_type), pointer :: EOS !< Equation of state structure
 
  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_file	Parameter file structure
	eos	Equation of state structure

Definition at line 806 of file MOM_EOS.F90.

  type(param_file_type), intent(in) :: param_file !< Parameter file structure
  type(EOS_type),        pointer    :: EOS !< Equation of state structure
  ! Local variables
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_EOS" ! This module's name.
  character(len=40)  :: tmpstr
 
  if (.not.associated(eos)) call eos_allocate(eos)
 
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
 
  call get_param(param_file, mdl, "EQN_OF_STATE", tmpstr, &
                 "EQN_OF_STATE determines which ocean equation of state "//&
                 "should be used.  Currently, the valid choices are "//&
                 '"LINEAR", "UNESCO", "WRIGHT", "NEMO" and "TEOS10". '//&
                 "This is only used if USE_EOS is true.", default=eos_default)
  select case (uppercase(tmpstr))
    case (eos_linear_string)
      eos%form_of_EOS = eos_linear
    case (eos_unesco_string)
      eos%form_of_EOS = eos_unesco
    case (eos_wright_string)
      eos%form_of_EOS = eos_wright
    case (eos_teos10_string)
      eos%form_of_EOS = eos_teos10
    case (eos_nemo_string)
      eos%form_of_EOS = eos_nemo
    case default
      call mom_error(fatal, "interpret_eos_selection: EQN_OF_STATE "//&
                              trim(tmpstr) // "in input file is invalid.")
  end select
  call mom_mesg('interpret_eos_selection: equation of state set to "' // &
                trim(tmpstr)//'"', 5)
 
  if (eos%form_of_EOS == eos_linear) then
    eos%Compressible = .false.
    call get_param(param_file, mdl, "RHO_T0_S0", eos%Rho_T0_S0, &
                 "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
                 "this is the density at T=0, S=0.", units="kg m-3", &
                 default=1000.0)
    call get_param(param_file, mdl, "DRHO_DT", eos%dRho_dT, &
                 "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
                 "this is the partial derivative of density with "//&
                 "temperature.", units="kg m-3 K-1", default=-0.2)
    call get_param(param_file, mdl, "DRHO_DS", eos%dRho_dS, &
                 "When EQN_OF_STATE="//trim(eos_linear_string)//", "//&
                 "this is the partial derivative of density with "//&
                 "salinity.", units="kg m-3 PSU-1", default=0.8)
  endif
 
  call get_param(param_file, mdl, "EOS_QUADRATURE", eos%EOS_quadrature, &
                 "If true, always use the generic (quadrature) code "//&
                 "code for the integrals of density.", default=.false.)
 
  call get_param(param_file, mdl, "TFREEZE_FORM", tmpstr, &
                 "TFREEZE_FORM determines which expression should be "//&
                 "used for the freezing point.  Currently, the valid "//&
                 'choices are "LINEAR", "MILLERO_78", "TEOS10"', &
                 default=tfreeze_default)
  select case (uppercase(tmpstr))
    case (tfreeze_linear_string)
      eos%form_of_TFreeze = tfreeze_linear
    case (tfreeze_millero_string)
      eos%form_of_TFreeze = tfreeze_millero
    case (tfreeze_teos10_string)
      eos%form_of_TFreeze = tfreeze_teos10
    case default
      call mom_error(fatal, "interpret_eos_selection:  TFREEZE_FORM "//&
                              trim(tmpstr) // "in input file is invalid.")
  end select
 
  if (eos%form_of_TFreeze == tfreeze_linear) then
    call get_param(param_file, mdl, "TFREEZE_S0_P0",eos%TFr_S0_P0, &
                 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
                 "this is the freezing potential temperature at "//&
                 "S=0, P=0.", units="deg C", default=0.0)
    call get_param(param_file, mdl, "DTFREEZE_DS",eos%dTFr_dS, &
                 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
                 "this is the derivative of the freezing potential "//&
                 "temperature with salinity.", &
                 units="deg C PSU-1", default=-0.054)
    call get_param(param_file, mdl, "DTFREEZE_DP",eos%dTFr_dP, &
                 "When TFREEZE_FORM="//trim(tfreeze_linear_string)//", "//&
                 "this is the derivative of the freezing potential "//&
                 "temperature with pressure.", &
                 units="deg C Pa-1", default=0.0)
  endif
 
  if ((eos%form_of_EOS == eos_teos10 .OR. eos%form_of_EOS == eos_nemo) .AND. &
      eos%form_of_TFreeze /= tfreeze_teos10) then
      call mom_error(fatal, "interpret_eos_selection:  EOS_TEOS10 or EOS_NEMO \n" //&
      "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 .")
  endif
 
 

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().

Referenced by mom_ice_shelf::initialize_ice_shelf(), and mom::initialize_mom().

Here is the call graph for this function:

Here is the caller graph for this function:

eos_manual_init()

subroutine, public mom_eos::eos_manual_init	(	type(eos_type), pointer	EOS,
		integer, intent(in), optional	form_of_EOS,
		integer, intent(in), optional	form_of_TFreeze,
		logical, intent(in), optional	EOS_quadrature,
		logical, intent(in), optional	Compressible,
		real, intent(in), optional	Rho_T0_S0,
		real, intent(in), optional	drho_dT,
		real, intent(in), optional	dRho_dS,
		real, intent(in), optional	TFr_S0_P0,
		real, intent(in), optional	dTFr_dS,
		real, intent(in), optional	dTFr_dp
	)

Manually initialized an EOS type (intended for unit testing of routines which need a specific EOS)

Parameters

	eos	Equation of state structure
[in]	form_of_eos	A coded integer indicating the equation of state to use.
[in]	form_of_tfreeze	A coded integer indicating the expression for the potential temperature of the freezing point.
[in]	eos_quadrature	If true, always use the generic (quadrature) code for the integrals of density.
[in]	compressible	If true, in situ density is a function of pressure.
[in]	rho_t0_s0	Density at T=0 degC and S=0 ppt [kg m-3]
[in]	drho_dt	Partial derivative of density with temperature in [kg m-3 degC-1]
[in]	drho_ds	Partial derivative of density with salinity in [kg m-3 ppt-1]
[in]	tfr_s0_p0	The freezing potential temperature at S=0, P=0 [degC].
[in]	dtfr_ds	The derivative of freezing point with salinity in [degC ppt-1].
[in]	dtfr_dp	The derivative of freezing point with pressure in [degC Pa-1].

Definition at line 907 of file MOM_EOS.F90.

  type(EOS_type),    pointer    :: EOS !< Equation of state structure
  integer, optional, intent(in) :: form_of_EOS !< A coded integer indicating the equation of state to use.
  integer, optional, intent(in) :: form_of_TFreeze !< A coded integer indicating the expression for
                                       !! the potential temperature of the freezing point.
  logical, optional, intent(in) :: EOS_quadrature !< If true, always use the generic (quadrature)
                                       !! code for the integrals of density.
  logical, optional, intent(in) :: Compressible  !< If true, in situ density is a function of pressure.
  real   , optional, intent(in) :: Rho_T0_S0 !< Density at T=0 degC and S=0 ppt [kg m-3]
  real   , optional, intent(in) :: drho_dT   !< Partial derivative of density with temperature
                                             !! in [kg m-3 degC-1]
  real   , optional, intent(in) :: dRho_dS   !< Partial derivative of density with salinity
                                             !! in [kg m-3 ppt-1]
  real   , optional, intent(in) :: TFr_S0_P0 !< The freezing potential temperature at S=0, P=0 [degC].
  real   , optional, intent(in) :: dTFr_dS   !< The derivative of freezing point with salinity
                                             !! in [degC ppt-1].
  real   , optional, intent(in) :: dTFr_dp   !< The derivative of freezing point with pressure
                                             !! in [degC Pa-1].
 
  if (present(form_of_eos    ))  eos%form_of_EOS     = form_of_eos
  if (present(form_of_tfreeze))  eos%form_of_TFreeze = form_of_tfreeze
  if (present(eos_quadrature ))  eos%EOS_quadrature  = eos_quadrature
  if (present(compressible   ))  eos%Compressible    = compressible
  if (present(rho_t0_s0      ))  eos%Rho_T0_S0       = rho_t0_s0
  if (present(drho_dt        ))  eos%drho_dT         = drho_dt
  if (present(drho_ds        ))  eos%dRho_dS         = drho_ds
  if (present(tfr_s0_p0      ))  eos%TFr_S0_P0       = tfr_s0_p0
  if (present(dtfr_ds        ))  eos%dTFr_dS         = dtfr_ds
  if (present(dtfr_dp        ))  eos%dTFr_dp         = dtfr_dp
 

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_s0	Density at T=0 degC and S=0 ppt [kg m-3]
[in]	drho_dt	Partial derivative of density with temperature [kg m-3 degC-1]
[in]	drho_ds	Partial derivative of density with salinity [kg m-3 ppt-1]
[in]	use_quadrature	If true, always use the generic (quadrature) code for the integrals of density.
	eos	Equation of state structure

Definition at line 958 of file MOM_EOS.F90.

  real,              intent(in) :: Rho_T0_S0 !< Density at T=0 degC and S=0 ppt [kg m-3]
  real,              intent(in) :: dRho_dT   !< Partial derivative of density with temperature [kg m-3 degC-1]
  real,              intent(in) :: dRho_dS   !< Partial derivative of density with salinity [kg m-3 ppt-1]
  logical, optional, intent(in) :: use_quadrature !< If true, always use the generic (quadrature)
                                             !! code for the integrals of density.
  type(EOS_type),    pointer    :: EOS       !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "MOM_EOS.F90: EOS_use_linear() called with an unassociated EOS_type EOS.")
 
  eos%form_of_EOS = eos_linear
  eos%Compressible = .false.
  eos%Rho_T0_S0 = rho_t0_s0
  eos%dRho_dT = drho_dt
  eos%dRho_dS = drho_ds
  eos%EOS_quadrature = .false.
  if (present(use_quadrature)) eos%EOS_quadrature = use_quadrature
 

References eos_linear, and mom_error_handler::mom_error().

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

Evaluation of the four bilinear shape fn and their gradients at (xi,eta)

Parameters

[in]	xi	The x position to evaluate
[in]	eta	The z position to evaluate
[out]	phi	The weights of the four corners at this point
[out]	dphidxi	The x-gradient of the weights of the four corners at this point
[out]	dphideta	The z-gradient of the weights of the four corners at this point

Definition at line 1917 of file MOM_EOS.F90.

  real,               intent(in)  :: xi  !< The x position to evaluate
  real,               intent(in)  :: eta !< The z position to evaluate
  real, dimension(4), intent(out) :: phi !< The weights of the four corners at this point
  real, dimension(4), intent(out) :: dphidxi  !< The x-gradient of the weights of the four
                                         !! corners at this point
  real, dimension(4), intent(out) :: dphideta !< The z-gradient of the weights of the four
                                         !! corners at this point
 
  ! The shape functions within the parent element are defined as shown here:
  !
  !    (-1,1) 2 o------------o 1 (1,1)
  !             |            |
  !             |            |
  !             |            |
  !             |            |
  !   (-1,-1) 3 o------------o 4 (1,-1)
  !
 
  phi(1) = 0.25 * ( 1 + xi ) * ( 1 + eta )
  phi(2) = 0.25 * ( 1 - xi ) * ( 1 + eta )
  phi(3) = 0.25 * ( 1 - xi ) * ( 1 - eta )
  phi(4) = 0.25 * ( 1 + xi ) * ( 1 - eta )
 
  dphidxi(1) = 0.25 * ( 1 + eta )
  dphidxi(2) = - 0.25 * ( 1 + eta )
  dphidxi(3) = - 0.25 * ( 1 - eta )
  dphidxi(4) = 0.25 * ( 1 - eta )
 
  dphideta(1) = 0.25 * ( 1 + xi )
  dphideta(2) = 0.25 * ( 1 - xi )
  dphideta(3) = - 0.25 * ( 1 - xi )
  dphideta(4) = - 0.25 * ( 1 + xi )
 

Referenced by compute_integral_quadratic().

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

Evaluation of the nine quadratic shape fn weights and their gradients at (xi,eta)

Parameters

[in]	xi	The x position to evaluate
[in]	eta	The z position to evaluate
[out]	phi	The weights of the 9 bilinear quadrature points at this point
[out]	dphidxi	The x-gradient of the weights of the 9 bilinear quadrature points corners at this point
[out]	dphideta	The z-gradient of the weights of the 9 bilinear quadrature points corners at this point

Definition at line 1956 of file MOM_EOS.F90.

 
  ! Arguments
  real,               intent(in)  :: xi  !< The x position to evaluate
  real,               intent(in)  :: eta !< The z position to evaluate
  real, dimension(9), intent(out) :: phi !< The weights of the 9 bilinear quadrature points
                                         !! at this point
  real, dimension(9), intent(out) :: dphidxi  !< The x-gradient of the weights of the 9 bilinear
                                         !! quadrature points corners at this point
  real, dimension(9), intent(out) :: dphideta !< The z-gradient of the weights of the 9 bilinear
                                         !! quadrature points corners at this point
 
  ! The quadratic shape functions within the parent element are defined as shown here:
  !
  !                 5 (0,1)
  !    (-1,1) 2 o------o------o 1 (1,1)
  !             |             |
  !             |   9 (0,0)   |
  !    (-1,0) 6 o      o      o 8 (1,0)
  !             |             |
  !             |             |
  !   (-1,-1) 3 o------o------o 4 (1,-1)
  !                 7 (0,-1)
  !
 
  phi(:)   = 0.0
  dphidxi(:)  = 0.0
  dphideta(:) = 0.0
 
  phi(1) = 0.25 * xi * ( 1 + xi ) * eta * ( 1 + eta )
  phi(2) = - 0.25 * xi * ( 1 - xi ) * eta * ( 1 + eta )
  phi(3) = 0.25 * xi * ( 1 - xi ) * eta * ( 1 - eta )
  phi(4) = - 0.25 * xi * ( 1 + xi ) * eta * ( 1 - eta )
  phi(5) = 0.5 * ( 1 + xi ) * ( 1 - xi ) * eta * ( 1 + eta )
  phi(6) = - 0.5 * xi * ( 1 - xi ) * ( 1 - eta ) * ( 1 + eta )
  phi(7) = - 0.5 * ( 1 - xi ) * ( 1 + xi ) * eta * ( 1 - eta )
  phi(8) = 0.5 * xi * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
  phi(9) = ( 1 - xi ) * ( 1 + xi ) * ( 1 - eta ) * ( 1 + eta )
 
  !dphidxi(1) = 0.25 * ( 1 + 2*xi ) * eta * ( 1 + eta )
  !dphidxi(2) = - 0.25 * ( 1 - 2*xi ) * eta * ( 1 + eta )
  !dphidxi(3) = 0.25 * ( 1 - 2*xi ) * eta * ( 1 - eta )
  !dphidxi(4) = - 0.25 * ( 1 + 2*xi ) * eta * ( 1 - eta )
  !dphidxi(5) = - xi * eta * ( 1 + eta )
  !dphidxi(6) = - 0.5 * ( 1 - 2*xi ) * ( 1 - eta ) * ( 1 + eta )
  !dphidxi(7) = xi * eta * ( 1 - eta )
  !dphidxi(8) = 0.5 * ( 1 + 2*xi ) * ( 1 - eta ) * ( 1 + eta )
  !dphidxi(9) = - 2 * xi * ( 1 - eta ) * ( 1 + eta )
 
  !dphideta(1) = 0.25 * xi * ( 1 + xi ) * ( 1 + 2*eta )
  !dphideta(2) = - 0.25 * xi * ( 1 - xi ) * ( 1 + 2*eta )
  !dphideta(3) = 0.25 * xi * ( 1 - xi ) * ( 1 - 2*eta )
  !dphideta(4) = - 0.25 * xi * ( 1 + xi ) * ( 1 - 2*eta )
  !dphideta(5) = 0.5 * ( 1 + xi ) * ( 1 - xi ) * ( 1 + 2*eta )
  !dphideta(6) = xi * ( 1 - xi ) * eta
  !dphideta(7) = - 0.5 * ( 1 - xi ) * ( 1 + xi ) * ( 1 - 2*eta )
  !dphideta(8) = - xi * ( 1 + xi ) * eta
  !dphideta(9) = - 2 * ( 1 - xi ) * ( 1 + xi ) * eta
 

Referenced by compute_integral_quadratic().

Here is the caller graph for this function:

extract_member_eos()

subroutine, public mom_eos::extract_member_eos	(	type(eos_type), pointer	EOS,
		integer, intent(out), optional	form_of_EOS,
		integer, intent(out), optional	form_of_TFreeze,
		logical, intent(out), optional	EOS_quadrature,
		logical, intent(out), optional	Compressible,
		real, intent(out), optional	Rho_T0_S0,
		real, intent(out), optional	drho_dT,
		real, intent(out), optional	dRho_dS,
		real, intent(out), optional	TFr_S0_P0,
		real, intent(out), optional	dTFr_dS,
		real, intent(out), optional	dTFr_dp
	)

Extractor routine for the EOS type if the members need to be accessed outside this module.

Parameters

	eos	Equation of state structure
[out]	form_of_eos	A coded integer indicating the equation of state to use.
[out]	form_of_tfreeze	A coded integer indicating the expression for the potential temperature of the freezing point.
[out]	eos_quadrature	If true, always use the generic (quadrature) code for the integrals of density.
[out]	compressible	If true, in situ density is a function of pressure.
[out]	rho_t0_s0	Density at T=0 degC and S=0 ppt [kg m-3]
[out]	drho_dt	Partial derivative of density with temperature in [kg m-3 degC-1]
[out]	drho_ds	Partial derivative of density with salinity in [kg m-3 ppt-1]
[out]	tfr_s0_p0	The freezing potential temperature at S=0, P=0 [degC].
[out]	dtfr_ds	The derivative of freezing point with salinity [degC PSU-1].
[out]	dtfr_dp	The derivative of freezing point with pressure [degC Pa-1].

Definition at line 2473 of file MOM_EOS.F90.

  type(EOS_type),    pointer     :: EOS !< Equation of state structure
  integer, optional, intent(out) :: form_of_EOS !< A coded integer indicating the equation of state to use.
  integer, optional, intent(out) :: form_of_TFreeze !< A coded integer indicating the expression for
                                       !! the potential temperature of the freezing point.
  logical, optional, intent(out) :: EOS_quadrature !< If true, always use the generic (quadrature)
                                       !! code for the integrals of density.
  logical, optional, intent(out) :: Compressible !< If true, in situ density is a function of pressure.
  real   , optional, intent(out) :: Rho_T0_S0 !< Density at T=0 degC and S=0 ppt [kg m-3]
  real   , optional, intent(out) :: drho_dT   !< Partial derivative of density with temperature
                                              !! in [kg m-3 degC-1]
  real   , optional, intent(out) :: dRho_dS   !< Partial derivative of density with salinity
                                              !! in [kg m-3 ppt-1]
  real   , optional, intent(out) :: TFr_S0_P0 !< The freezing potential temperature at S=0, P=0 [degC].
  real   , optional, intent(out) :: dTFr_dS   !< The derivative of freezing point with salinity
                                              !! [degC PSU-1].
  real   , optional, intent(out) :: dTFr_dp   !< The derivative of freezing point with pressure
                                              !! [degC Pa-1].
 
  if (present(form_of_eos    ))  form_of_eos     = eos%form_of_EOS
  if (present(form_of_tfreeze))  form_of_tfreeze = eos%form_of_TFreeze
  if (present(eos_quadrature ))  eos_quadrature  = eos%EOS_quadrature
  if (present(compressible   ))  compressible    = eos%Compressible
  if (present(rho_t0_s0      ))  rho_t0_s0       = eos%Rho_T0_S0
  if (present(drho_dt        ))  drho_dt         = eos%drho_dT
  if (present(drho_ds        ))  drho_ds         = eos%dRho_dS
  if (present(tfr_s0_p0      ))  tfr_s0_p0       = eos%TFr_S0_P0
  if (present(dtfr_ds        ))  dtfr_ds         = eos%dTFr_dS
  if (present(dtfr_dp        ))  dtfr_dp         = eos%dTFr_dp
 

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,
		real, intent(in), optional	z_tol
	)

Find the depth at which the reconstructed pressure matches P_tgt.

Parameters

[in]	t_t	Potential temperatue at the cell top [degC]
[in]	t_b	Potential temperatue at the cell bottom [degC]
[in]	s_t	Salinity at the cell top [ppt]
[in]	s_b	Salinity at the cell bottom [ppt]
[in]	z_t	Absolute height of top of cell [Z ~> m]. (Boussinesq ????)
[in]	z_b	Absolute height of bottom of cell [Z ~> m].
[in]	p_t	Anomalous pressure of top of cell, relative to g*rho_ref*z_t [Pa]
[in]	p_tgt	Target pressure at height z_out, relative to g*rho_ref*z_out [Pa]
[in]	rho_ref	Reference density with which calculation are anomalous to
[in]	g_e	Gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
	eos	Equation of state structure
[out]	p_b	Pressure at the bottom of the cell [Pa]
[out]	z_out	Absolute depth at which anomalous pressure = p_tgt [Z ~> m].
[in]	z_tol	The tolerance in finding z_out [Z ~> m].

Definition at line 1470 of file MOM_EOS.F90.

  real,           intent(in)  :: T_t !< Potential temperatue at the cell top [degC]
  real,           intent(in)  :: T_b !< Potential temperatue at the cell bottom [degC]
  real,           intent(in)  :: S_t !< Salinity at the cell top [ppt]
  real,           intent(in)  :: S_b !< Salinity at the cell bottom [ppt]
  real,           intent(in)  :: z_t !< Absolute height of top of cell [Z ~> m].   (Boussinesq ????)
  real,           intent(in)  :: z_b !< Absolute height of bottom of cell [Z ~> m].
  real,           intent(in)  :: P_t !< Anomalous pressure of top of cell, relative to g*rho_ref*z_t [Pa]
  real,           intent(in)  :: P_tgt !< Target pressure at height z_out, relative to g*rho_ref*z_out [Pa]
  real,           intent(in)  :: rho_ref !< Reference density with which calculation are anomalous to
  real,           intent(in)  :: G_e !< Gravitational acceleration [m2 Z-1 s-2 ~> m s-2]
  type(EOS_type), pointer     :: EOS !< Equation of state structure
  real,           intent(out) :: P_b !< Pressure at the bottom of the cell [Pa]
  real,           intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m].
  real, optional, intent(in)  :: z_tol !< The tolerance in finding z_out [Z ~> m].
  ! Local variables
  real :: top_weight, bottom_weight, rho_anom, w_left, w_right, GxRho, dz, dp, F_guess, F_l, F_r
  real :: Pa, Pa_left, Pa_right, Pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz
 
  gxrho = g_e * rho_ref
 
  ! Anomalous pressure difference across whole cell
  dp = frac_dp_at_pos(t_t, t_b, s_t, s_b, z_t, z_b, rho_ref, g_e, 1.0, eos)
 
  p_b = p_t + dp ! Anomalous pressure at bottom of cell
 
  if (p_tgt <= p_t ) then
    z_out = z_t
    return
  endif
 
  if (p_tgt >= p_b) then
    z_out = z_b
    return
  endif
 
  f_l = 0.
  pa_left = p_t - p_tgt ! Pa_left < 0
  f_r = 1.
  pa_right = p_b - p_tgt ! Pa_right > 0
  pa_tol = gxrho * 1.e-5 ! 1e-5 has dimensions of m, but should be converted to the units of z.
  if (present(z_tol)) pa_tol = gxrho * z_tol
  f_guess = f_l - pa_left / ( pa_right -pa_left ) * ( f_r - f_l )
  pa = pa_right - pa_left ! To get into iterative loop
  do while ( abs(pa) > pa_tol )
 
    z_out = z_t + ( z_b - z_t ) * f_guess
    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 )
 
    if (pa<pa_left) then
      write(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
      stop 'Blurgh! Too negative'
    elseif (pa<0.) then
      pa_left = pa
      f_l = f_guess
    elseif (pa>pa_right) then
      write(0,*) pa_left,pa,pa_right,p_t-p_tgt,p_b-p_tgt
      stop 'Blurgh! Too positive'
    elseif (pa>0.) then
      pa_right = pa
      f_r = f_guess
    else ! Pa == 0
      return
    endif
    f_guess = f_l - pa_left / ( pa_right -pa_left ) * ( f_r - f_l )
 
  enddo
 

References frac_dp_at_pos().

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

[in]	t_t	Potential temperatue at the cell top [degC]
[in]	t_b	Potential temperatue at the cell bottom [degC]
[in]	s_t	Salinity at the cell top [ppt]
[in]	s_b	Salinity at the cell bottom [ppt]
[in]	z_t	The geometric height at the top of the layer [Z ~> m]
[in]	z_b	The geometric height at the bottom of the layer [Z ~> m]
[in]	rho_ref	A mean density [kg m-3], that is subtracted out to reduce the magnitude of each of the integrals.
[in]	g_e	The Earth's gravitational acceleration [m s-2]
[in]	pos	The fractional vertical position, 0 to 1 [nondim].
	eos	Equation of state structure

Definition at line 1542 of file MOM_EOS.F90.

  real,           intent(in)  :: T_t !< Potential temperatue at the cell top [degC]
  real,           intent(in)  :: T_b !< Potential temperatue at the cell bottom [degC]
  real,           intent(in)  :: S_t !< Salinity at the cell top [ppt]
  real,           intent(in)  :: S_b !< Salinity at the cell bottom [ppt]
  real,           intent(in)  :: z_t !< The geometric height at the top of the layer [Z ~> m]
  real,           intent(in)  :: z_b !< The geometric height at the bottom of the layer [Z ~> m]
  real,           intent(in)  :: rho_ref !< A mean density [kg m-3], that is subtracted out to
                                     !! reduce the magnitude of each of the integrals.
  real,           intent(in)  :: G_e !< The Earth's gravitational acceleration [m s-2]
  real,           intent(in)  :: pos !< The fractional vertical position, 0 to 1 [nondim].
  type(EOS_type), pointer     :: EOS !< Equation of state structure
  ! Local variables
  real, parameter :: C1_90 = 1.0/90.0  ! Rational constants.
  real :: dz, top_weight, bottom_weight, rho_ave
  real, dimension(5) :: T5, S5, p5, rho5
  integer :: n
 
  do n=1,5
    ! Evalute density at five quadrature points
    bottom_weight = 0.25*real(n-1) * pos
    top_weight = 1.0 - bottom_weight
    ! Salinity and temperature points are linearly interpolated
    s5(n) = top_weight * s_t + bottom_weight * s_b
    t5(n) = top_weight * t_t + bottom_weight * t_b
    p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( g_e * rho_ref )
  enddo
  call calculate_density_array(t5, s5, p5, rho5, 1, 5, eos)
  rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref
 
  ! Use Boole's rule to estimate the average density
  rho_ave = c1_90*(7.0*(rho5(1)+rho5(5)) + 32.0*(rho5(2)+rho5(4)) + 12.0*rho5(3))
 
  dz = ( z_t - z_b ) * pos
  frac_dp_at_pos = g_e * dz * rho_ave

References calculate_density_array().

Referenced by find_depth_of_pressure_in_cell().

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,
		real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in), optional	bathyT,
		real, intent(in), optional	dz_neglect,
		logical, intent(in), optional	useMassWghtInterp
	)

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.

Parameters

[in]	hii	Ocean horizontal index structures for the input arrays
[in]	hio	Ocean horizontal index structures for the output arrays
[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	z_t	Height at the top of the layer in depth units [Z ~> m].
[in]	z_b	Height at the bottom of the layer [Z ~> m].
[in]	rho_ref	A mean density [kg m-3], that is subtracted out to reduce the magnitude of each of the integrals.
[in]	rho_0	A density [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_e	The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
	eos	Equation of state structure
[out]	dpa	The change in the pressure anomaly across the layer [Pa].
[out]	intz_dpa	The integral through the thickness of the layer of
[out]	intx_dpa	The integral in x of the difference between the
[out]	inty_dpa	The integral in y of the difference between the
[in]	bathyt	The depth of the bathymetry [Z ~> m].
[in]	dz_neglect	A miniscule thickness change [Z ~> m].
[in]	usemasswghtinterp	If true, uses mass weighting to interpolate T/S for top and bottom integrals.

Definition at line 733 of file MOM_EOS.F90.

  type(hor_index_type), intent(in)  :: HII !< Ocean horizontal index structures for the input arrays
  type(hor_index_type), intent(in)  :: HIO !< Ocean horizontal index structures for the output arrays
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: T   !< Potential temperature referenced to the surface [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: S   !< Salinity [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< Height at the top of the layer in depth units [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is subtracted out to
                                           !! reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [kg m-3], that is used to calculate the
                                           !! pressure (as p~=-z*rho_0*G_e) used in the equation of state.
  real,                 intent(in)  :: G_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
  type(EOS_type),       pointer     :: EOS !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [Pa Z ~> Pa m].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the x grid spacing [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the y grid spacing [Pa].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
              optional, intent(in)  :: bathyT !< The depth of the bathymetry [Z ~> m].
  real,       optional, intent(in)  :: dz_neglect !< A miniscule thickness change [Z ~> m].
  logical,    optional, intent(in)  :: useMassWghtInterp !< If true, uses mass weighting to
                                           !! interpolate T/S for top and bottom integrals.
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "int_density_dz called with an unassociated EOS_type EOS.")
 
  if (eos%EOS_quadrature) then
    call int_density_dz_generic(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio, &
                                eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
                                bathyt, dz_neglect, usemasswghtinterp)
  else ; select case (eos%form_of_EOS)
    case (eos_linear)
      call int_density_dz_linear(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio,  &
                                 eos%Rho_T0_S0, eos%dRho_dT, eos%dRho_dS, &
                                 dpa, intz_dpa, intx_dpa, inty_dpa, &
                                 bathyt, dz_neglect, usemasswghtinterp)
    case (eos_wright)
      call int_density_dz_wright(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio,  &
                                 dpa, intz_dpa, intx_dpa, inty_dpa, &
                                 bathyt, dz_neglect, usemasswghtinterp)
    case default
      call int_density_dz_generic(t, s, z_t, z_b, rho_ref, rho_0, g_e, hii, hio,  &
                                  eos, dpa, intz_dpa, intx_dpa, inty_dpa, &
                                  bathyt, dz_neglect, usemasswghtinterp)
  end select ; endif
 

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

Here is the call graph for this function:

int_density_dz_generic()

subroutine, public 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,
		real, dimension(hii%isd:hii%ied,hii%jsd:hii%jed), intent(in), optional	bathyT,
		real, intent(in), optional	dz_neglect,
		logical, intent(in), optional	useMassWghtInterp
	)

This subroutine calculates (by numerical quadrature) integrals of pressure anomalies across layers, which are required for calculating the finite-volume form pressure accelerations in a Boussinesq model.

Parameters

[in]	hii	Horizontal index type for input variables.
[in]	hio	Horizontal index type for output variables.
[in]	t	Potential temperature of the layer [degC].
[in]	s	Salinity of the layer [ppt].
[in]	z_t	Height at the top of the layer in depth units [Z ~> m].
[in]	z_b	Height at the bottom of the layer [Z ~> m].
[in]	rho_ref	A mean density [kg m-3], that is subtracted out to reduce the magnitude of each of the integrals.
[in]	rho_0	A density [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_e	The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
	eos	Equation of state structure
[out]	dpa	The change in the pressure anomaly
[out]	intz_dpa	The integral through the thickness of the
[out]	intx_dpa	The integral in x of the difference between
[out]	inty_dpa	The integral in y of the difference between
[in]	bathyt	The depth of the bathymetry [Z ~> m].
[in]	dz_neglect	A miniscule thickness change [Z ~> m].
[in]	usemasswghtinterp	If true, uses mass weighting to interpolate T/S for top and bottom integrals.

Definition at line 984 of file MOM_EOS.F90.

  type(hor_index_type), intent(in)  :: HII !< Horizontal index type for input variables.
  type(hor_index_type), intent(in)  :: HIO !< Horizontal index type for output variables.
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: T  !< Potential temperature of the layer [degC].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: S  !< Salinity of the layer [ppt].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< Height at the top of the layer in depth units [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is
                                          !! subtracted out to reduce the magnitude
                                          !! of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [kg m-3], that is used
                                          !! to calculate the pressure (as p~=-z*rho_0*G_e)
                                          !! used in the equation of state.
  real,                 intent(in)  :: G_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
  type(EOS_type),       pointer     :: EOS !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly
                                          !! across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the
                                          !! layer of the pressure anomaly relative to the
                                          !! anomaly at the top of the layer [Pa Z ~> Pa m].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intx_dpa !< The integral in x of the difference between
                                          !! the pressure anomaly at the top and bottom of the
                                          !! layer divided by the x grid spacing [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: inty_dpa !< The integral in y of the difference between
                                          !! the pressure anomaly at the top and bottom of the
                                          !! layer divided by the y grid spacing [Pa].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
              optional, intent(in)  :: bathyT !< The depth of the bathymetry [Z ~> m].
  real,       optional, intent(in)  :: dz_neglect !< A miniscule thickness change [Z ~> m].
  logical,    optional, intent(in)  :: useMassWghtInterp !< If true, uses mass weighting to
                                          !! interpolate T/S for top and bottom integrals.
  real :: T5(5), S5(5), p5(5), r5(5)
  real :: rho_anom   ! The depth averaged density anomaly [kg m-3].
  real :: w_left, w_right
  real, parameter :: C1_90 = 1.0/90.0  ! Rational constants.
  real :: GxRho, I_Rho
  real :: dz         ! The layer thickness [Z ~> m].
  real :: hWght      ! A pressure-thickness below topography [Z ~> m].
  real :: hL, hR     ! Pressure-thicknesses of the columns to the left and right [Z ~> m].
  real :: iDenom     ! The inverse of the denominator in the weights [Z-2 ~> m-2].
  real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim].
  real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim].
  real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim].
  real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim].
  real :: intz(5)    ! The gravitational acceleration times the integrals of density
                     ! with height at the 5 sub-column locations [Pa].
  logical :: do_massWeight ! Indicates whether to do mass weighting.
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m, n, ioff, joff
 
  ioff = hio%idg_offset - hii%idg_offset
  joff = hio%jdg_offset - hii%jdg_offset
 
  ! These array bounds work for the indexing convention of the input arrays, but
  ! on the computational domain defined for the output arrays.
  isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
  jsq = hio%JscB + joff ; jeq = hio%JecB + joff
  is = hio%isc + ioff ; ie = hio%iec + ioff
  js = hio%jsc + joff ; je = hio%jec + joff
 
  gxrho = g_e * rho_0
  i_rho = 1.0 / rho_0
 
  do_massweight = .false.
  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
    do_massweight = .true.
    if (.not.present(bathyt)) call mom_error(fatal, "int_density_dz_generic: "//&
        "bathyT must be present if useMassWghtInterp is present and true.")
    if (.not.present(dz_neglect)) call mom_error(fatal, "int_density_dz_generic: "//&
        "dz_neglect must be present if useMassWghtInterp is present and true.")
  endif ; endif
 
  do j=jsq,jeq+1 ; do i=isq,ieq+1
    dz = z_t(i,j) - z_b(i,j)
    do n=1,5
      t5(n) = t(i,j) ; s5(n) = s(i,j)
      p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
    enddo
    call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
 
    ! Use Bode's rule to estimate the pressure anomaly change.
    rho_anom = c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))
    dpa(i-ioff,j-joff) = g_e*dz*rho_anom
    ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of
    ! the pressure anomaly.
    if (present(intz_dpa)) intz_dpa(i-ioff,j-joff) = 0.5*g_e*dz**2 * &
          (rho_anom - c1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
  enddo ; enddo
 
  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., -bathyt(i,j)-z_t(i+1,j), -bathyt(i+1,j)-z_t(i,j))
    if (hwght > 0.) then
      hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
      hr = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif
 
    intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
    do m=2,4
      ! T, S, and z are interpolated in the horizontal.  The z interpolation
      ! is linear, but for T and S it may be thickness weighted.
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
      dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i+1,j) - z_b(i+1,j))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
      p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i+1,j))
      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo
      call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
 
    ! Use Bode's rule to estimate the pressure anomaly change.
      intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
    enddo
    ! Use Bode's rule to integrate the bottom pressure anomaly values in x.
    intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                           12.0*intz(3))
  enddo ; enddo ; endif
 
  if (present(inty_dpa)) then ; do j=jsq,jeq ; do i=is,ie
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., -bathyt(i,j)-z_t(i,j+1), -bathyt(i,j+1)-z_t(i,j))
    if (hwght > 0.) then
      hl = (z_t(i,j) - z_b(i,j)) + dz_neglect
      hr = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif
 
    intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i-ioff,j-joff+1)
    do m=2,4
      ! T, S, and z are interpolated in the horizontal.  The z interpolation
      ! is linear, but for T and S it may be thickness weighted.
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
      dz = wt_l*(z_t(i,j) - z_b(i,j)) + wt_r*(z_t(i,j+1) - z_b(i,j+1))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
      p5(1) = -gxrho*(wt_l*z_t(i,j) + wt_r*z_t(i,j+1))
      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1)
        p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo
      call calculate_density(t5, s5, p5, r5, 1, 5, eos, rho_ref)
 
    ! Use Bode's rule to estimate the pressure anomaly change.
      intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)))
    enddo
    ! Use Bode's rule to integrate the values.
    inty_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                                     12.0*intz(3))
  enddo ; enddo ; endif

References mom_error_handler::mom_error().

Referenced by int_density_dz().

Here is the call graph for this function:

Here is the caller graph for this function:

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)	dz_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
	)

Compute pressure gradient force integrals by quadrature for the case where T and S are linear profiles.

Parameters

[in]	hii	Ocean horizontal index structures for the input arrays
[in]	hio	Ocean horizontal index structures for the output arrays
[in]	t_t	Potential temperatue at the cell top [degC]
[in]	t_b	Potential temperatue at the cell bottom [degC]
[in]	s_t	Salinity at the cell top [ppt]
[in]	s_b	Salinity at the cell bottom [ppt]
[in]	z_t	The geometric height at the top of the layer,
[in]	z_b	The geometric height at the bottom of the layer [Z ~> m].
[in]	rho_ref	A mean density [kg m-3], that is subtracted out to reduce the magnitude of each of the integrals.
[in]	rho_0	A density [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_e	The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
[in]	dz_subroundoff	A miniscule thickness change [Z ~> m].
[in]	bathyt	The depth of the bathymetry [Z ~> m].
	eos	Equation of state structure
[out]	dpa	The change in the pressure anomaly across the layer [Pa].
[out]	intz_dpa	The integral through the thickness of the layer of
[out]	intx_dpa	The integral in x of the difference between the
[out]	inty_dpa	The integral in y of the difference between the
[in]	usemasswghtinterp	If true, uses mass weighting to interpolate T/S for top and bottom integrals.

Definition at line 1171 of file MOM_EOS.F90.

  type(hor_index_type), intent(in)  :: HII !< Ocean horizontal index structures for the input arrays
  type(hor_index_type), intent(in)  :: HIO !< Ocean horizontal index structures for the output arrays
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: T_t !< Potential temperatue at the cell top [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: T_b !< Potential temperatue at the cell bottom [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: S_t !< Salinity at the cell top [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: S_b !< Salinity at the cell bottom [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< The geometric height at the top of the layer,
                                           !! in depth units [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< The geometric height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is subtracted out to
                                           !! reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [kg m-3], that is used to calculate the
                                           !! pressure (as p~=-z*rho_0*G_e) used in the equation of state.
  real,                 intent(in)  :: G_e !< The Earth's gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
  real,                 intent(in)  :: dz_subroundoff !< A miniscule thickness change [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: bathyT !< The depth of the bathymetry [Z ~> m].
  type(EOS_type),       pointer     :: EOS !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [Pa Z].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the x grid spacing [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the y grid spacing [Pa].
  logical,    optional, intent(in)  :: useMassWghtInterp !< If true, uses mass weighting to
                                           !! interpolate T/S for top and bottom integrals.
! This subroutine calculates (by numerical quadrature) 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.
!
! It is assumed that the salinity and temperature profiles are linear in the
! vertical. The top and bottom values within each layer are provided and
! a linear interpolation is used to compute intermediate values.
 
  ! Local variables
  real :: T5((5*HIO%iscB+1):(5*(HIO%iecB+2)))  ! Temperatures along a line of subgrid locations [degC].
  real :: S5((5*HIO%iscB+1):(5*(HIO%iecB+2)))  ! Salinities along a line of subgrid locations [ppt].
  real :: p5((5*HIO%iscB+1):(5*(HIO%iecB+2)))  ! Pressures along a line of subgrid locations [Pa].
  real :: r5((5*HIO%iscB+1):(5*(HIO%iecB+2)))  ! Densities along a line of subgrid locations [kg m-3].
  real :: T15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Temperatures at an array of subgrid locations [degC].
  real :: S15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Salinities at an array of subgrid locations [ppt].
  real :: p15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Pressures at an array of subgrid locations [Pa].
  real :: r15((15*HIO%iscB+1):(15*(HIO%iecB+1))) ! Densities at an array of subgrid locations [kg m-3].
  real :: wt_t(5), wt_b(5)          ! Top and bottom weights [nondim].
  real :: rho_anom                  ! A density anomaly [kg m-3].
  real :: w_left, w_right           ! Left and right weights [nondim].
  real :: intz(5)    ! The gravitational acceleration times the integrals of density
                     ! with height at the 5 sub-column locations [Pa].
  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant [nondim].
  real :: GxRho                     ! Gravitational acceleration times density [kg m-1 Z-1 s-2 ~> kg m-2 s-2].
  real :: I_Rho                     ! The inverse of the reference density [m3 kg-1].
  real :: dz(HIO%iscB:HIO%iecB+1)   ! Layer thicknesses at tracer points [Z ~> m].
  real :: dz_x(5,HIO%iscB:HIO%iecB) ! Layer thicknesses along an x-line of subrid locations [Z ~> m].
  real :: dz_y(5,HIO%isc:HIO%iec)   ! Layer thicknesses along a y-line of subrid locations [Z ~> m].
  real :: weight_t, weight_b        ! Nondimensional weights of the top and bottom.
  real :: massWeightToggle          ! A nondimensional toggle factor (0 or 1).
  real :: Ttl, Tbl, Ttr, Tbr        ! Temperatures at the velocity cell corners [degC].
  real :: Stl, Sbl, Str, Sbr        ! Salinities at the velocity cell corners [ppt].
  real :: hWght                     ! A topographically limited thicknes weight [Z ~> m].
  real :: hL, hR                    ! Thicknesses to the left and right [Z ~> m].
  real :: iDenom                    ! The denominator of the thickness weight expressions [Z-2 ~> m-2].
  integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n
  integer :: iin, jin, ioff, joff
  integer :: pos
 
  ioff = hio%idg_offset - hii%idg_offset
  joff = hio%jdg_offset - hii%jdg_offset
 
  isq = hio%IscB ; ieq = hio%IecB ; jsq = hio%JscB ; jeq = hio%JecB
 
  gxrho = g_e * rho_0
  i_rho = 1.0 / rho_0
  massweighttoggle = 0.
  if (present(usemasswghtinterp)) then
    if (usemasswghtinterp) massweighttoggle = 1.
  endif
 
  do n = 1, 5
    wt_t(n) = 0.25 * real(5-n)
    wt_b(n) = 1.0 - wt_t(n)
  enddo
 
  ! =============================
  ! 1. Compute vertical integrals
  ! =============================
  do j=jsq,jeq+1
    jin = j+joff
    do i = isq,ieq+1 ; iin = i+ioff
      dz(i) = z_t(iin,jin) - z_b(iin,jin)
      do n=1,5
        p5(i*5+n) = -gxrho*(z_t(iin,jin) - 0.25*real(n-1)*dz(i))
        ! Salinity and temperature points are linearly interpolated
        s5(i*5+n) = wt_t(n) * s_t(iin,jin) + wt_b(n) * s_b(iin,jin)
        t5(i*5+n) = wt_t(n) * t_t(iin,jin) + wt_b(n) * t_b(iin,jin)
      enddo
    enddo
    call calculate_density_array(t5, s5, p5, r5, 1, (ieq-isq+2)*5, eos, rho_ref )
 
    do i=isq,ieq+1 ; iin = i+ioff
    ! Use Bode's rule to estimate the pressure anomaly change.
      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))
      dpa(i,j) = g_e*dz(i)*rho_anom
      if (present(intz_dpa)) then
      ! Use a Bode's-rule-like fifth-order accurate estimate of
      ! the double integral of the pressure anomaly.
        intz_dpa(i,j) = 0.5*g_e*dz(i)**2 * &
                (rho_anom - c1_90*(16.0*(r5(i*5+4)-r5(i*5+2)) + 7.0*(r5(i*5+5)-r5(i*5+1))) )
      endif
    enddo
  enddo ! end loops on j
 
 
  ! ==================================================
  ! 2. Compute horizontal integrals in the x direction
  ! ==================================================
  if (present(intx_dpa)) then ; do j=hio%jsc,hio%jec ; jin = j+joff
    do i=isq,ieq ; iin = i+ioff
      ! Corner values of T and S
      ! hWght is the distance measure by which the cell is violation of
      ! hydrostatic consistency. For large hWght we bias the interpolation
      ! of T,S along the top and bottom integrals, almost like thickness
      ! weighting.
      ! Note: To work in terrain following coordinates we could offset
      ! this distance by the layer thickness to replicate other models.
      hwght = massweighttoggle * &
              max(0., -bathyt(iin,jin)-z_t(iin+1,jin), -bathyt(iin+1,jin)-z_t(iin,jin))
      if (hwght > 0.) then
        hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
        hr = (z_t(iin+1,jin) - z_b(iin+1,jin)) + dz_subroundoff
        hwght = hwght * ( (hl-hr)/(hl+hr) )**2
        idenom = 1./( hwght*(hr + hl) + hl*hr )
        ttl = ( (hwght*hr)*t_t(iin+1,jin) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
        ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin+1,jin) ) * idenom
        tbl = ( (hwght*hr)*t_b(iin+1,jin) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
        tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin+1,jin) ) * idenom
        stl = ( (hwght*hr)*s_t(iin+1,jin) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
        str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin+1,jin) ) * idenom
        sbl = ( (hwght*hr)*s_b(iin+1,jin) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
        sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin+1,jin) ) * idenom
      else
        ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin+1,jin); tbr = t_b(iin+1,jin)
        stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin+1,jin); sbr = s_b(iin+1,jin)
      endif
 
      do m=2,4
        w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
        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))
 
        ! Salinity and temperature points are linearly interpolated in
        ! the horizontal. The subscript (1) refers to the top value in
        ! the vertical profile while subscript (5) refers to the bottom
        ! value in the vertical profile.
        pos = i*15+(m-2)*5
        t15(pos+1) = w_left*ttl + w_right*ttr
        t15(pos+5) = w_left*tbl + w_right*tbr
 
        s15(pos+1) = w_left*stl + w_right*str
        s15(pos+5) = w_left*sbl + w_right*sbr
 
        p15(pos+1) = -gxrho*(w_left*z_t(iin,jin) + w_right*z_t(iin+1,jin))
 
        ! Pressure
        do n=2,5
          p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_x(m,i)
        enddo
 
        ! Salinity and temperature (linear interpolation in the vertical)
        do n=2,4
          weight_t = 0.25 * real(5-n)
          weight_b = 1.0 - weight_t
          s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
          t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
        enddo
      enddo
    enddo
 
    call calculate_density(t15, s15, p15, r15, 1, 15*(ieq-isq+1), eos, rho_ref)
 
    do i=isq,ieq ; iin = i+ioff
      intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
 
      ! Use Bode's rule to estimate the pressure anomaly change.
      do m = 2,4
        pos = i*15+(m-2)*5
        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)) + &
                          12.0*r15(pos+3)))
      enddo
      ! Use Bode's rule to integrate the bottom pressure anomaly values in x.
      intx_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                             12.0*intz(3))
    enddo
  enddo ; endif
 
  ! ==================================================
  ! 3. Compute horizontal integrals in the y direction
  ! ==================================================
  if (present(inty_dpa)) then ; do j=jsq,jeq ; jin = j+joff
    do i=hio%isc,hio%iec ; iin = i+ioff
    ! Corner values of T and S
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, almost like thickness
    ! weighting.
    ! Note: To work in terrain following coordinates we could offset
    ! this distance by the layer thickness to replicate other models.
      hwght = massweighttoggle * &
              max(0., -bathyt(i,j)-z_t(iin,jin+1), -bathyt(i,j+1)-z_t(iin,jin))
      if (hwght > 0.) then
        hl = (z_t(iin,jin) - z_b(iin,jin)) + dz_subroundoff
        hr = (z_t(iin,jin+1) - z_b(iin,jin+1)) + dz_subroundoff
        hwght = hwght * ( (hl-hr)/(hl+hr) )**2
        idenom = 1./( hwght*(hr + hl) + hl*hr )
        ttl = ( (hwght*hr)*t_t(iin,jin+1) + (hwght*hl + hr*hl)*t_t(iin,jin) ) * idenom
        ttr = ( (hwght*hl)*t_t(iin,jin) + (hwght*hr + hr*hl)*t_t(iin,jin+1) ) * idenom
        tbl = ( (hwght*hr)*t_b(iin,jin+1) + (hwght*hl + hr*hl)*t_b(iin,jin) ) * idenom
        tbr = ( (hwght*hl)*t_b(iin,jin) + (hwght*hr + hr*hl)*t_b(iin,jin+1) ) * idenom
        stl = ( (hwght*hr)*s_t(iin,jin+1) + (hwght*hl + hr*hl)*s_t(iin,jin) ) * idenom
        str = ( (hwght*hl)*s_t(iin,jin) + (hwght*hr + hr*hl)*s_t(iin,jin+1) ) * idenom
        sbl = ( (hwght*hr)*s_b(iin,jin+1) + (hwght*hl + hr*hl)*s_b(iin,jin) ) * idenom
        sbr = ( (hwght*hl)*s_b(iin,jin) + (hwght*hr + hr*hl)*s_b(iin,jin+1) ) * idenom
      else
        ttl = t_t(iin,jin); tbl = t_b(iin,jin); ttr = t_t(iin,jin+1); tbr = t_b(iin,jin+1)
        stl = s_t(iin,jin); sbl = s_b(iin,jin); str = s_t(iin,jin+1); sbr = s_b(iin,jin+1)
      endif
 
      do m=2,4
        w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
        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))
 
        ! Salinity and temperature points are linearly interpolated in
        ! the horizontal. The subscript (1) refers to the top value in
        ! the vertical profile while subscript (5) refers to the bottom
        ! value in the vertical profile.
        pos = i*15+(m-2)*5
        t15(pos+1) = w_left*ttl + w_right*ttr
        t15(pos+5) = w_left*tbl + w_right*tbr
 
        s15(pos+1) = w_left*stl + w_right*str
        s15(pos+5) = w_left*sbl + w_right*sbr
 
        p15(pos+1) = -gxrho*(w_left*z_t(iin,jin) + w_right*z_t(iin,jin+1))
 
        ! Pressure
        do n=2,5 ; p15(pos+n) = p15(pos+n-1) + gxrho*0.25*dz_y(m,i) ; enddo
 
        ! Salinity and temperature (linear interpolation in the vertical)
        do n=2,4
          weight_t = 0.25 * real(5-n)
          weight_b = 1.0 - weight_t
          s15(pos+n) = weight_t * s15(pos+1) + weight_b * s15(pos+5)
          t15(pos+n) = weight_t * t15(pos+1) + weight_b * t15(pos+5)
        enddo
      enddo
    enddo
 
    call calculate_density_array(t15(15*hio%isc+1:), s15(15*hio%isc+1:), p15(15*hio%isc+1:), &
                                 r15(15*hio%isc+1:), 1, 15*(hio%iec-hio%isc+1), eos, rho_ref)
    do i=hio%isc,hio%iec ; iin = i+ioff
      intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
 
      ! Use Bode's rule to estimate the pressure anomaly change.
      do m = 2,4
        pos = i*15+(m-2)*5
        intz(m) = g_e*dz_y(m,i)*( c1_90*(7.0*(r15(pos+1)+r15(pos+5)) + &
                                         32.0*(r15(pos+2)+r15(pos+4)) + &
                                         12.0*r15(pos+3)))
      enddo
      ! Use Bode's rule to integrate the values.
      inty_dpa(i,j) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                             12.0*intz(3))
    enddo
  enddo ; endif
 

References calculate_density_array().

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%jsd:hio%jed), intent(out), optional	intx_dpa,
		real, dimension(hio%isd:hio%ied,hio%jsdb:hio%jedb), intent(out), optional	inty_dpa
	)

Compute pressure gradient force integrals for the case where T and S are parabolic profiles.

Parameters

[in]	hii	Ocean horizontal index structures for the input arrays
[in]	hio	Ocean horizontal index structures for the output arrays
[in]	t	Potential temperature referenced to the surface [degC]
[in]	t_t	Potential temperatue at the cell top [degC]
[in]	t_b	Potential temperatue at the cell bottom [degC]
[in]	s	Salinity [ppt]
[in]	s_t	Salinity at the cell top [ppt]
[in]	s_b	Salinity at the cell bottom [ppt]
[in]	z_t	Height at the top of the layer [Z ~> m].
[in]	z_b	Height at the bottom of the layer [Z ~> m].
[in]	rho_ref	A mean density [kg m-3], that is subtracted out to reduce the magnitude of each of the integrals.
[in]	rho_0	A density [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_e	The Earth's gravitational acceleration [m s-2]
	eos	Equation of state structure
[out]	dpa	The change in the pressure anomaly across the layer [Pa].
[out]	intz_dpa	The integral through the thickness of the layer of
[out]	intx_dpa	The integral in x of the difference between the
[out]	inty_dpa	The integral in y of the difference between the

Definition at line 1585 of file MOM_EOS.F90.

 
  type(hor_index_type), intent(in)  :: HII !< Ocean horizontal index structures for the input arrays
  type(hor_index_type), intent(in)  :: HIO !< Ocean horizontal index structures for the output arrays
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: T   !< Potential temperature referenced to the surface [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: T_t !< Potential temperatue at the cell top [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: T_b !< Potential temperatue at the cell bottom [degC]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: S   !< Salinity [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: S_t !< Salinity at the cell top [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: S_b !< Salinity at the cell bottom [ppt]
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_t !< Height at the top of the layer [Z ~> m].
  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
                        intent(in)  :: z_b !< Height at the bottom of the layer [Z ~> m].
  real,                 intent(in)  :: rho_ref !< A mean density [kg m-3], that is subtracted out to
                                           !! reduce the magnitude of each of the integrals.
  real,                 intent(in)  :: rho_0 !< A density [kg m-3], that is used to calculate the
                                           !! pressure (as p~=-z*rho_0*G_e) used in the equation of state.
  real,                 intent(in)  :: G_e !< The Earth's gravitational acceleration [m s-2]
  type(EOS_type),       pointer     :: EOS !< Equation of state structure
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
                        intent(out) :: dpa !< The change in the pressure anomaly across the layer [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intz_dpa !< The integral through the thickness of the layer of
                                           !! the pressure anomaly relative to the anomaly at the
                                           !! top of the layer [Pa Z ~> Pa m].
  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
              optional, intent(out) :: intx_dpa !< The integral in x of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the x grid spacing [Pa].
  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
              optional, intent(out) :: inty_dpa !< The integral in y of the difference between the
                                           !! pressure anomaly at the top and bottom of the layer
                                           !! divided by the y grid spacing [Pa].
 
! This subroutine calculates (by numerical quadrature) 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.
!
! It is assumed that the salinity and temperature profiles are linear in the
! vertical. The top and bottom values within each layer are provided and
! a linear interpolation is used to compute intermediate values.
 
  ! Local variables
  real :: T5(5), S5(5), p5(5), r5(5)
  real :: rho_anom
  real :: w_left, w_right, intz(5)
  real, parameter :: C1_90 = 1.0/90.0  ! Rational constants.
  real :: GxRho, I_Rho
  real :: dz
  real :: weight_t, weight_b
  real :: s0, s1, s2                   ! parabola coefficients for S [ppt]
  real :: t0, t1, t2                   ! parabola coefficients for T [degC]
  real :: xi                           ! normalized coordinate
  real :: T_top, T_mid, T_bot
  real :: S_top, S_mid, S_bot
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m, n, ioff, joff
  real, dimension(4) :: x, y
  real, dimension(9) :: S_node, T_node, p_node, r_node
 
 
  call mom_error(fatal, &
    "int_density_dz_generic_ppm: the implementation is not done yet, contact developer")
 
  ioff = hio%idg_offset - hii%idg_offset
  joff = hio%jdg_offset - hii%jdg_offset
 
  ! These array bounds work for the indexing convention of the input arrays, but
  ! on the computational domain defined for the output arrays.
  isq = hio%IscB + ioff ; ieq = hio%IecB + ioff
  jsq = hio%JscB + joff ; jeq = hio%JecB + joff
  is = hio%isc + ioff ; ie = hio%iec + ioff
  js = hio%jsc + joff ; je = hio%jec + joff
 
  gxrho = g_e * rho_0
  i_rho = 1.0 / rho_0
 
  ! =============================
  ! 1. Compute vertical integrals
  ! =============================
  do j=jsq,jeq+1 ; do i=isq,ieq+1
    dz = z_t(i,j) - z_b(i,j)
 
    ! Coefficients of the parabola for S
    s0 = s_t(i,j)
    s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
    s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0*s(i,j) )
 
    ! Coefficients of the parabola for T
    t0 = t_t(i,j)
    t1 = 6.0 * t(i,j) - 4.0 * t_t(i,j) - 2.0 * t_b(i,j)
    t2 = 3.0 * ( t_t(i,j) + t_b(i,j) - 2.0*t(i,j) )
 
    do n=1,5
      p5(n) = -gxrho*(z_t(i,j) - 0.25*real(n-1)*dz)
 
      ! Parabolic reconstruction for T and S
      xi = 0.25 * ( n - 1 )
      s5(n) = s0 + s1 * xi + s2 * xi**2
      t5(n) = t0 + t1 * xi + t2 * xi**2
    enddo
 
    call calculate_density(t5, s5, p5, r5, 1, 5, eos)
 
    ! Use Bode's rule to estimate the pressure anomaly change.
    !rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) - &
    !       rho_ref
 
    rho_anom = 1000.0 + s(i,j) - rho_ref
    dpa(i-ioff,j-joff) = g_e*dz*rho_anom
 
    ! Use a Bode's-rule-like fifth-order accurate estimate of
    ! the double integral of the pressure anomaly.
    !r5 = r5 - rho_ref
    !if (present(intz_dpa)) intz_dpa(i,j) = 0.5*G_e*dz**2 * &
    !      (rho_anom - C1_90*(16.0*(r5(4)-r5(2)) + 7.0*(r5(5)-r5(1))) )
 
    intz_dpa(i-ioff,j-joff) = 0.5 * g_e * dz**2 * ( 1000.0 - rho_ref + s0 + s1/3.0 + &
                                    s2/6.0 )
  enddo ; enddo ! end loops on j and i
 
  ! ==================================================
  ! 2. Compute horizontal integrals in the x direction
  ! ==================================================
  if (present(intx_dpa)) then ; do j=js,je ; do i=isq,ieq
    intz(1) = dpa(i-ioff,j-joff) ; intz(5) = dpa(i+1-ioff,j-joff)
    do m=2,4
      w_left = 0.25*real(5-m) ; w_right = 1.0-w_left
      dz = w_left*(z_t(i,j) - z_b(i,j)) + w_right*(z_t(i+1,j) - z_b(i+1,j))
 
      ! Salinity and temperature points are linearly interpolated in
      ! the horizontal. The subscript (1) refers to the top value in
      ! the vertical profile while subscript (5) refers to the bottom
      ! value in the vertical profile.
      t_top = w_left*t_t(i,j) + w_right*t_t(i+1,j)
      t_mid = w_left*t(i,j)   + w_right*t(i+1,j)
      t_bot = w_left*t_b(i,j) + w_right*t_b(i+1,j)
 
      s_top = w_left*s_t(i,j) + w_right*s_t(i+1,j)
      s_mid = w_left*s(i,j)   + w_right*s(i+1,j)
      s_bot = w_left*s_b(i,j) + w_right*s_b(i+1,j)
 
      p5(1) = -gxrho*(w_left*z_t(i,j) + w_right*z_t(i+1,j))
 
      ! Pressure
      do n=2,5
        p5(n) = p5(n-1) + gxrho*0.25*dz
      enddo
 
      ! Coefficients of the parabola for S
      s0 = s_top
      s1 = 6.0 * s_mid - 4.0 * s_top - 2.0 * s_bot
      s2 = 3.0 * ( s_top + s_bot - 2.0*s_mid )
 
      ! Coefficients of the parabola for T
      t0 = t_top
      t1 = 6.0 * t_mid - 4.0 * t_top - 2.0 * t_bot
      t2 = 3.0 * ( t_top + t_bot - 2.0*t_mid )
 
      do n=1,5
        ! Parabolic reconstruction for T and S
        xi = 0.25 * ( n - 1 )
        s5(n) = s0 + s1 * xi + s2 * xi**2
        t5(n) = t0 + t1 * xi + t2 * xi**2
      enddo
 
      call calculate_density(t5, s5, p5, r5, 1, 5, eos)
 
    ! Use Bode's rule to estimate the pressure anomaly change.
      intz(m) = g_e*dz*( c1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + &
                            12.0*r5(3)) - rho_ref)
    enddo
    intx_dpa(i-ioff,j-joff) = c1_90*(7.0*(intz(1)+intz(5)) + 32.0*(intz(2)+intz(4)) + &
                           12.0*intz(3))
 
    ! Use Gauss quadrature rule to compute integral
 
    ! The following coordinates define the quadrilateral on which the integral
    ! is computed
    x(1) = 1.0
    x(2) = 0.0
    x(3) = 0.0
    x(4) = 1.0
    y(1) = z_t(i+1,j)
    y(2) = z_t(i,j)
    y(3) = z_b(i,j)
    y(4) = z_b(i+1,j)
 
    t_node = 0.0
    p_node = 0.0
 
    ! Nodal values for S
 
    ! Parabolic reconstruction on the left
    s0 = s_t(i,j)
    s1 = 6.0 * s(i,j) - 4.0 * s_t(i,j) - 2.0 * s_b(i,j)
    s2 = 3.0 * ( s_t(i,j) + s_b(i,j) - 2.0 * s(i,j) )
    s_node(2) = s0
    s_node(6) = s0 + 0.5 * s1 + 0.25 * s2
    s_node(3) = s0 + s1 + s2
 
    ! Parabolic reconstruction on the left
    s0 = s_t(i+1,j)
    s1 = 6.0 * s(i+1,j) - 4.0 * s_t(i+1,j) - 2.0 * s_b(i+1,j)
    s2 = 3.0 * ( s_t(i+1,j) + s_b(i+1,j) - 2.0 * s(i+1,j) )
    s_node(1) = s0
    s_node(8) = s0 + 0.5 * s1 + 0.25 * s2
    s_node(4) = s0 + s1 + s2
 
    s_node(5) = 0.5 * ( s_node(2) + s_node(1) )
    s_node(9) = 0.5 * ( s_node(6) + s_node(8) )
    s_node(7) = 0.5 * ( s_node(3) + s_node(4) )
 
    call calculate_density( t_node, s_node, p_node, r_node, 1, 9, eos )
    r_node = r_node - rho_ref
 
    call compute_integral_quadratic( x, y, r_node, intx_dpa(i-ioff,j-joff) )
 
    intx_dpa(i-ioff,j-joff) = intx_dpa(i-ioff,j-joff) * g_e
 
  enddo ; enddo ; endif
 
  ! ==================================================
  ! 3. Compute horizontal integrals in the y direction
  ! ==================================================
  if (present(inty_dpa)) then
    call mom_error(warning, "int_density_dz_generic_ppm still needs to be written for inty_dpa!")
    do j=jsq,jeq ; do i=is,ie
 
      inty_dpa(i-ioff,j-joff) = 0.0
 
    enddo ; enddo
  endif
 

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

Here is the call graph for this function:

int_spec_vol_dp_generic()

subroutine, public 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,
		real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in), optional	bathyP,
		real, intent(in), optional	dP_neglect,
		logical, intent(in), optional	useMassWghtInterp
	)

This subroutine calculates integrals of specific volume anomalies in pressure across layers, 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 quadrature to do the integrals.

Parameters

[in]	hi	A horizontal index type structure.
[in]	t	Potential temperature of the layer [degC].
[in]	s	Salinity of the layer [ppt].
[in]	p_t	Pressure atop the layer [Pa].
[in]	p_b	Pressure below the layer [Pa].
[in]	alpha_ref	A 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 alpha_ref alters the effects of roundoff, and answers do change.
	eos	Equation of state structure
[out]	dza	The change in the geopotential anomaly
[out]	intp_dza	The integral in pressure through the
[out]	intx_dza	The integral in x of the difference
[out]	inty_dza	The integral in y of the difference
[in]	halo_size	The width of halo points on which to calculate dza.
[in]	bathyp	The pressure at the bathymetry [Pa]
[in]	dp_neglect	A miniscule pressure change with the same units as p_t (Pa?)
[in]	usemasswghtinterp	If true, uses mass weighting to interpolate T/S for top and bottom integrals.

Definition at line 2024 of file MOM_EOS.F90.

  type(hor_index_type), intent(in)  :: HI !< A horizontal index type structure.
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: T  !< Potential temperature of the layer [degC].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: S  !< Salinity of the layer [ppt].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_t !< Pressure atop the layer [Pa].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_b !< Pressure below the layer [Pa].
  real,                 intent(in)  :: alpha_ref !< A 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 alpha_ref
                            !! alters the effects of roundoff, and answers do change.
  type(EOS_type),       pointer     :: EOS !< Equation of state structure
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(out) :: dza !< The change in the geopotential anomaly
                            !! across the layer [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(out) :: intp_dza !< The integral in pressure through the
                            !! layer of the geopotential anomaly relative to the anomaly
                            !! at the bottom of the layer [Pa m2 s-2].
  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
              optional, intent(out) :: intx_dza  !< The integral in x of the difference
                            !! between the geopotential anomaly at the top and bottom of
                            !! the layer divided by the x grid spacing [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
              optional, intent(out) :: inty_dza  !< The integral in y of the difference
                            !! between the geopotential anomaly at the top and bottom of
                            !! the layer divided by the y grid spacing [m2 s-2].
  integer,    optional, intent(in)  :: halo_size !< The width of halo points on which to calculate dza.
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(in)  :: bathyP !< The pressure at the bathymetry [Pa]
  real,       optional, intent(in)  :: dP_neglect !< A miniscule pressure change with
                                             !! the same units as p_t (Pa?)
  logical,    optional, intent(in)  :: useMassWghtInterp !< If true, uses mass weighting
                            !! to interpolate T/S for top and bottom integrals.
 
!   This subroutine calculates 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| < 0.34.
 
  real :: T5(5), S5(5), p5(5), a5(5)
  real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1].
  real :: dp         ! The pressure change through a layer [Pa].
!  real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa].
  real :: hWght      ! A pressure-thickness below topography [Pa].
  real :: hL, hR     ! Pressure-thicknesses of the columns to the left and right [Pa].
  real :: iDenom     ! The inverse of the denominator in the weights [Pa-2].
  real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim].
  real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim].
  real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim].
  real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim].
  real :: intp(5)    ! The integrals of specific volume with pressure at the
                     ! 5 sub-column locations [m2 s-2].
  logical :: do_massWeight ! Indicates whether to do mass weighting.
  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant.
  integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, n, halo
 
  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
  halo = 0 ; if (present(halo_size)) halo = max(halo_size,0)
  ish = hi%isc-halo ; ieh = hi%iec+halo ; jsh = hi%jsc-halo ; jeh = hi%jec+halo
  if (present(intx_dza)) then ; ish = min(isq,ish) ; ieh = max(ieq+1,ieh); endif
  if (present(inty_dza)) then ; jsh = min(jsq,jsh) ; jeh = max(jeq+1,jeh); endif
 
  do_massweight = .false.
  if (present(usemasswghtinterp)) then ; if (usemasswghtinterp) then
    do_massweight = .true.
    if (.not.present(bathyp)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
        "bathyP must be present if useMassWghtInterp is present and true.")
    if (.not.present(dp_neglect)) call mom_error(fatal, "int_spec_vol_dp_generic: "//&
        "dP_neglect must be present if useMassWghtInterp is present and true.")
  endif ; endif
 
  do j=jsh,jeh ; do i=ish,ieh
    dp = p_b(i,j) - p_t(i,j)
    do n=1,5
      t5(n) = t(i,j) ; s5(n) = s(i,j)
      p5(n) = p_b(i,j) - 0.25*real(n-1)*dp
    enddo
    call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode's rule to estimate the interface height anomaly change.
    alpha_anom = c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3))
    dza(i,j) = dp*alpha_anom
    ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of
    ! the interface height anomaly.
    if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
          (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
  enddo ; enddo
 
  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif
 
    intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
 
      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness wekghted.
      p5(1) = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
      dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i+1,j) - p_t(i+1,j))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i+1,j)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i+1,j)
 
      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - 0.25*dp
      enddo
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode's rule to estimate the interface height anomaly change.
      intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
                                12.0*a5(3)))
    enddo
    ! Use Bode's rule to integrate the interface height anomaly values in x.
    intx_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
                           12.0*intp(3))
  enddo ; enddo ; endif
 
  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation of
    ! T & S along the top and bottom integrals, akin to thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif
 
    intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
 
      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness wekghted.
      p5(1) = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
      dp = wt_l*(p_b(i,j) - p_t(i,j)) + wt_r*(p_b(i,j+1) - p_t(i,j+1))
      t5(1) = wtt_l*t(i,j) + wtt_r*t(i,j+1)
      s5(1) = wtt_l*s(i,j) + wtt_r*s(i,j+1)
      do n=2,5
        t5(n) = t5(1) ; s5(n) = s5(1) ; p5(n) = p5(n-1) - 0.25*dp
      enddo
      call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode's rule to estimate the interface height anomaly change.
      intp(m) = dp*( c1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + &
                                12.0*a5(3)))
    enddo
    ! Use Bode's rule to integrate the interface height anomaly values in y.
    inty_dza(i,j) = c1_90*(7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4)) + &
                           12.0*intp(3))
  enddo ; enddo ; endif
 

References mom_error_handler::mom_error().

Referenced by int_specific_vol_dp().

Here is the call graph for this function:

Here is the caller graph for this function:

int_spec_vol_dp_generic_plm()

subroutine, public mom_eos::int_spec_vol_dp_generic_plm	(	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)	p_t,
		real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)	p_b,
		real, intent(in)	alpha_ref,
		real, intent(in)	dP_neglect,
		real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in)	bathyP,
		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,
		logical, intent(in), optional	useMassWghtInterp
	)

This subroutine calculates integrals of specific volume anomalies in pressure across layers, 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 quadrature to do the integrals.

Parameters

[in]	hi	A horizontal index type structure.
[in]	t_t	Potential temperature at the top of the layer [degC].
[in]	t_b	Potential temperature at the bottom of the layer [degC].
[in]	s_t	Salinity at the top the layer [ppt].
[in]	s_b	Salinity at the bottom the layer [ppt].
[in]	p_t	Pressure atop the layer [Pa].
[in]	p_b	Pressure below the layer [Pa].
[in]	alpha_ref	A 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 alpha_ref alters the effects of roundoff, and answers do change.
[in]	dp_neglect	A miniscule pressure change with the same units as p_t (Pa?)
[in]	bathyp	The pressure at the bathymetry [Pa]
	eos	Equation of state structure
[out]	dza	The change in the geopotential anomaly
[out]	intp_dza	The integral in pressure through the
[out]	intx_dza	The integral in x of the difference
[out]	inty_dza	The integral in y of the difference
[in]	usemasswghtinterp	If true, uses mass weighting to interpolate T/S for top and bottom integrals.

Definition at line 2214 of file MOM_EOS.F90.

  type(hor_index_type), intent(in)  :: HI !< A horizontal index type structure.
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: T_t  !< Potential temperature at the top of the layer [degC].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: T_b  !< Potential temperature at the bottom of the layer [degC].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: S_t  !< Salinity at the top the layer [ppt].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: S_b  !< Salinity at the bottom the layer [ppt].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_t !< Pressure atop the layer [Pa].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_b !< Pressure below the layer [Pa].
  real,                 intent(in)  :: alpha_ref !< A 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 alpha_ref
                            !! alters the effects of roundoff, and answers do change.
  real,                 intent(in)  :: dP_neglect !< A miniscule pressure change with
                                             !! the same units as p_t (Pa?)
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: bathyP !< The pressure at the bathymetry [Pa]
  type(EOS_type),       pointer     :: EOS !< Equation of state structure
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(out) :: dza !< The change in the geopotential anomaly
                            !! across the layer [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(out) :: intp_dza !< The integral in pressure through the
                            !! layer of the geopotential anomaly relative to the anomaly
                            !! at the bottom of the layer [Pa m2 s-2].
  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
              optional, intent(out) :: intx_dza  !< The integral in x of the difference
                            !! between the geopotential anomaly at the top and bottom of
                            !! the layer divided by the x grid spacing [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
              optional, intent(out) :: inty_dza  !< The integral in y of the difference
                            !! between the geopotential anomaly at the top and bottom of
                            !! the layer divided by the y grid spacing [m2 s-2].
  logical,    optional, intent(in)  :: useMassWghtInterp !< If true, uses mass weighting
                            !! to interpolate T/S for top and bottom integrals.
 
!   This subroutine calculates 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| < 0.34.
 
  real, dimension(5) :: T5, S5, p5, a5
  real, dimension(15) :: T15, S15, p15, a15
  real :: wt_t(5), wt_b(5)
  real :: T_top, T_bot, S_top, S_bot, P_top, P_bot
 
  real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1].
  real :: dp         ! The pressure change through a layer [Pa].
  real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa].
  real :: hWght      ! A pressure-thickness below topography [Pa].
  real :: hL, hR     ! Pressure-thicknesses of the columns to the left and right [Pa].
  real :: iDenom     ! The inverse of the denominator in the weights [Pa-2].
  real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim].
  real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim].
  real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim].
  real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim].
  real :: intp(5)    ! The integrals of specific volume with pressure at the
                     ! 5 sub-column locations [m2 s-2].
  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant.
  logical :: do_massWeight ! Indicates whether to do mass weighting.
  integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos
 
  isq = hi%IscB ; ieq = hi%IecB ; jsq = hi%JscB ; jeq = hi%JecB
 
  do_massweight = .false.
  if (present(usemasswghtinterp)) do_massweight = usemasswghtinterp
 
  do n = 1, 5 ! Note that these are reversed from int_density_dz.
    wt_t(n) = 0.25 * real(n-1)
    wt_b(n) = 1.0 - wt_t(n)
  enddo
 
  ! =============================
  ! 1. Compute vertical integrals
  ! =============================
  do j=jsq,jeq+1; do i=isq,ieq+1
    dp = p_b(i,j) - p_t(i,j)
    do n=1,5 ! T, S and p are linearly interpolated in the vertical.
      p5(n) = wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j)
      s5(n) = wt_t(n) * s_t(i,j) + wt_b(n) * s_b(i,j)
      t5(n) = wt_t(n) * t_t(i,j) + wt_b(n) * t_b(i,j)
    enddo
    call calculate_spec_vol(t5, s5, p5, a5, 1, 5, eos, alpha_ref)
 
    ! Use Bode's rule to estimate the interface height anomaly change.
    alpha_anom = c1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3))
    dza(i,j) = dp*alpha_anom
    ! Use a Bode's-rule-like fifth-order accurate estimate of the double integral of
    ! the interface height anomaly.
    if (present(intp_dza)) intp_dza(i,j) = 0.5*dp**2 * &
          (alpha_anom - c1_90*(16.0*(a5(4)-a5(2)) + 7.0*(a5(5)-a5(1))) )
  enddo ; enddo
 
  ! ==================================================
  ! 2. Compute horizontal integrals in the x direction
  ! ==================================================
  if (present(intx_dza)) then ; do j=hi%jsc,hi%jec ; do i=isq,ieq
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, almost like thickness
    ! weighting. Note: To work in terrain following coordinates we could
    ! offset this distance by the layer thickness to replicate other models.
    hwght = 0.0
    if (do_massweight) &
      hwght =  max(0., bathyp(i,j)-p_t(i+1,j), bathyp(i+1,j)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i+1,j) - p_t(i+1,j)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif
 
    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
 
      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness wekghted.
      p_top = wt_l*p_t(i,j) + wt_r*p_t(i+1,j)
      p_bot = wt_l*p_b(i,j) + wt_r*p_b(i+1,j)
      t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i+1,j)
      t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i+1,j)
      s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i+1,j)
      s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i+1,j)
      dp_90(m) = c1_90*(p_bot - p_top)
 
      ! Salinity, temperature and pressure with linear interpolation in the vertical.
      pos = (m-2)*5
      do n=1,5
        p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot
        s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
        t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
      enddo
    enddo
 
    call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref)
 
    intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
    do m=2,4
      ! Use Bode's rule to estimate the interface height anomaly change.
      ! The integrals at the ends of the segment are already known.
      pos = (m-2)*5
      intp(m) = dp_90(m)*((7.0*(a15(pos+1)+a15(pos+5)) + &
                          32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
    enddo
    ! Use Bode's rule to integrate the interface height anomaly values in x.
    intx_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
                           12.0*intp(3))
  enddo ; enddo ; endif
 
  ! ==================================================
  ! 3. Compute horizontal integrals in the y direction
  ! ==================================================
  if (present(inty_dza)) then ; do j=jsq,jeq ; do i=hi%isc,hi%iec
    ! hWght is the distance measure by which the cell is violation of
    ! hydrostatic consistency. For large hWght we bias the interpolation
    ! of T,S along the top and bottom integrals, like thickness weighting.
    hwght = 0.0
    if (do_massweight) &
      hwght = max(0., bathyp(i,j)-p_t(i,j+1), bathyp(i,j+1)-p_t(i,j))
    if (hwght > 0.) then
      hl = (p_b(i,j) - p_t(i,j)) + dp_neglect
      hr = (p_b(i,j+1) - p_t(i,j+1)) + dp_neglect
      hwght = hwght * ( (hl-hr)/(hl+hr) )**2
      idenom = 1.0 / ( hwght*(hr + hl) + hl*hr )
      hwt_ll = (hwght*hl + hr*hl) * idenom ; hwt_lr = (hwght*hr) * idenom
      hwt_rr = (hwght*hr + hr*hl) * idenom ; hwt_rl = (hwght*hl) * idenom
    else
      hwt_ll = 1.0 ; hwt_lr = 0.0 ; hwt_rr = 1.0 ; hwt_rl = 0.0
    endif
 
    do m=2,4
      wt_l = 0.25*real(5-m) ; wt_r = 1.0-wt_l
      wtt_l = wt_l*hwt_ll + wt_r*hwt_rl ; wtt_r = wt_l*hwt_lr + wt_r*hwt_rr
 
      ! T, S, and p are interpolated in the horizontal.  The p interpolation
      ! is linear, but for T and S it may be thickness wekghted.
      p_top = wt_l*p_t(i,j) + wt_r*p_t(i,j+1)
      p_bot = wt_l*p_b(i,j) + wt_r*p_b(i,j+1)
      t_top = wtt_l*t_t(i,j) + wtt_r*t_t(i,j+1)
      t_bot = wtt_l*t_b(i,j) + wtt_r*t_b(i,j+1)
      s_top = wtt_l*s_t(i,j) + wtt_r*s_t(i,j+1)
      s_bot = wtt_l*s_b(i,j) + wtt_r*s_b(i,j+1)
      dp_90(m) = c1_90*(p_bot - p_top)
 
      ! Salinity, temperature and pressure with linear interpolation in the vertical.
      pos = (m-2)*5
      do n=1,5
        p15(pos+n) = wt_t(n) * p_top + wt_b(n) * p_bot
        s15(pos+n) = wt_t(n) * s_top + wt_b(n) * s_bot
        t15(pos+n) = wt_t(n) * t_top + wt_b(n) * t_bot
      enddo
    enddo
 
    call calculate_spec_vol(t15, s15, p15, a15, 1, 15, eos, alpha_ref)
 
    intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
    do m=2,4
      ! Use Bode's rule to estimate the interface height anomaly change.
      ! The integrals at the ends of the segment are already known.
      pos = (m-2)*5
      intp(m) = dp_90(m) * ((7.0*(a15(pos+1)+a15(pos+5)) + &
                            32.0*(a15(pos+2)+a15(pos+4))) + 12.0*a15(pos+3))
    enddo
    ! Use Bode's rule to integrate the interface height anomaly values in x.
    inty_dza(i,j) = c1_90*((7.0*(intp(1)+intp(5)) + 32.0*(intp(2)+intp(4))) + &
                           12.0*intp(3))
  enddo ; enddo ; endif
 

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,
		real, dimension(hi%isd:hi%ied,hi%jsd:hi%jed), intent(in), optional	bathyP,
		real, intent(in), optional	dP_tiny,
		logical, intent(in), optional	useMassWghtInterp
	)

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]	hi	The horizontal index structure
[in]	t	Potential temperature referenced to the surface [degC]
[in]	s	Salinity [ppt]
[in]	p_t	Pressure at the top of the layer [Pa].
[in]	p_b	Pressure at the bottom of the layer [Pa].
[in]	alpha_ref	A 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.
	eos	Equation of state structure
[out]	dza	The change in the geopotential anomaly across
[out]	intp_dza	The integral in pressure through the layer of the
[out]	intx_dza	The integral in x of the difference between the
[out]	inty_dza	The integral in y of the difference between the
[in]	halo_size	The width of halo points on which to calculate dza.
[in]	bathyp	The pressure at the bathymetry [Pa]
[in]	dp_tiny	A miniscule pressure change with the same units as p_t (Pa?)
[in]	usemasswghtinterp	If true, uses mass weighting to interpolate T/S for top and bottom integrals.

Definition at line 665 of file MOM_EOS.F90.

  type(hor_index_type), intent(in)  :: HI  !< The horizontal index structure
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: T   !< Potential temperature referenced to the surface [degC]
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: S   !< Salinity [ppt]
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_t !< Pressure at the top of the layer [Pa].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(in)  :: p_b !< Pressure at the bottom of the layer [Pa].
  real,                 intent(in)  :: alpha_ref !< A 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.
  type(EOS_type),       pointer     :: EOS !< Equation of state structure
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
                        intent(out) :: dza !< The change in the geopotential anomaly across
                            !! the layer [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(out) :: intp_dza !< The integral in pressure through the layer of the
                            !! geopotential anomaly relative to the anomaly at the bottom of the
                            !! layer [Pa m2 s-2].
  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
              optional, intent(out) :: intx_dza !< The integral in x of the difference between the
                            !! geopotential anomaly at the top and bottom of the layer divided by
                            !! the x grid spacing [m2 s-2].
  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
              optional, intent(out) :: inty_dza !< The integral in y of the difference between the
                            !! geopotential anomaly at the top and bottom of the layer divided by
                            !! the y grid spacing [m2 s-2].
  integer,    optional, intent(in)  :: halo_size !< The width of halo points on which to calculate dza.
  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
              optional, intent(in)  :: bathyP  !< The pressure at the bathymetry [Pa]
  real,       optional, intent(in)  :: dP_tiny !< A miniscule pressure change with
                                               !! the same units as p_t (Pa?)
  logical,    optional, intent(in)  :: useMassWghtInterp !< If true, uses mass weighting
                            !! to interpolate T/S for top and bottom integrals.
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "int_specific_vol_dp called with an unassociated EOS_type EOS.")
 
  if (eos%EOS_quadrature) then
    call int_spec_vol_dp_generic(t, s, p_t, p_b, alpha_ref, hi, eos, &
                                 dza, intp_dza, intx_dza, inty_dza, halo_size, &
                                 bathyp, dp_tiny, usemasswghtinterp)
  else ; select case (eos%form_of_EOS)
    case (eos_linear)
      call int_spec_vol_dp_linear(t, s, p_t, p_b, alpha_ref, hi, eos%Rho_T0_S0, &
                                  eos%dRho_dT, eos%dRho_dS, dza, intp_dza, &
                                  intx_dza, inty_dza, halo_size, &
                                  bathyp, dp_tiny, usemasswghtinterp)
    case (eos_wright)
      call int_spec_vol_dp_wright(t, s, p_t, p_b, alpha_ref, hi, dza, &
                                  intp_dza, intx_dza, inty_dza, halo_size, &
                                  bathyp, dp_tiny, usemasswghtinterp)
    case default
      call int_spec_vol_dp_generic(t, s, p_t, p_b, alpha_ref, hi, eos, &
                                   dza, intp_dza, intx_dza, inty_dza, halo_size, &
                                   bathyp, dp_tiny, usemasswghtinterp)
  end select ; endif
 

References eos_linear, eos_wright, int_spec_vol_dp_generic(), mom_eos_linear::int_spec_vol_dp_linear(), 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().

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

eos	Equation of state structure

Definition at line 796 of file MOM_EOS.F90.

  type(EOS_type), pointer :: EOS !< Equation of state structure
 
  if (.not.associated(eos)) call mom_error(fatal, &
    "query_compressible called with an unassociated EOS_type EOS.")
 
  query_compressible = eos%compressible

References mom_error_handler::mom_error().

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

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

The default equation of state.

Definition at line 123 of file MOM_EOS.F90.

character*(10), parameter :: EOS_DEFAULT = EOS_WRIGHT_STRING !< The default equation of state

Referenced by eos_init().

eos_linear

integer, parameter, public mom_eos::eos_linear = 1

A named integer specifying an equation of state.

Definition at line 112 of file MOM_EOS.F90.

integer, parameter, public :: EOS_LINEAR = 1 !< A named integer specifying an equation of state

Referenced by calculate_compress_array(), calculate_density_array(), calculate_density_derivs_array(), calculate_density_derivs_scalar(), calculate_density_scalar(), calculate_density_second_derivs_array(), calculate_density_second_derivs_scalar(), calculate_spec_vol_array(), calculate_spec_vol_scalar(), calculate_specific_vol_derivs(), eos_init(), eos_use_linear(), int_density_dz(), int_specific_vol_dp(), and mom_neutral_diffusion::ndiff_unit_tests_discontinuous().

eos_linear_string

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

private

A string for specifying the equation of state.

Definition at line 118 of file MOM_EOS.F90.

character*(10), parameter :: EOS_LINEAR_STRING = "LINEAR" !< A string for specifying the equation of state

Referenced by eos_init().

eos_nemo

integer, parameter, public mom_eos::eos_nemo = 5

A named integer specifying an equation of state.

Definition at line 116 of file MOM_EOS.F90.

integer, parameter, public :: EOS_NEMO   = 5 !< A named integer specifying an equation of state

Referenced by calculate_compress_array(), calculate_density_array(), calculate_density_derivs_array(), calculate_density_scalar(), calculate_spec_vol_array(), calculate_spec_vol_scalar(), calculate_specific_vol_derivs(), convert_temp_salt_for_teos10(), and eos_init().

eos_nemo_string

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

private

A string for specifying the equation of state.

Definition at line 122 of file MOM_EOS.F90.

character*(10), parameter :: EOS_NEMO_STRING   = "NEMO"   !< A string for specifying the equation of state

Referenced by eos_init().

eos_teos10

integer, parameter, public mom_eos::eos_teos10 = 4

A named integer specifying an equation of state.

Definition at line 115 of file MOM_EOS.F90.

integer, parameter, public :: EOS_TEOS10 = 4 !< A named integer specifying an equation of state

Referenced by calculate_compress_array(), calculate_density_array(), calculate_density_derivs_array(), calculate_density_derivs_scalar(), calculate_density_scalar(), calculate_density_second_derivs_array(), calculate_density_second_derivs_scalar(), calculate_spec_vol_array(), calculate_spec_vol_scalar(), calculate_specific_vol_derivs(), convert_temp_salt_for_teos10(), and eos_init().

eos_teos10_string

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

private

A string for specifying the equation of state.

Definition at line 121 of file MOM_EOS.F90.

character*(10), parameter :: EOS_TEOS10_STRING = "TEOS10" !< A string for specifying the equation of state

Referenced by eos_init().

eos_unesco

integer, parameter, public mom_eos::eos_unesco = 2

A named integer specifying an equation of state.

Definition at line 113 of file MOM_EOS.F90.

integer, parameter, public :: EOS_UNESCO = 2 !< A named integer specifying an equation of state

Referenced by calculate_compress_array(), calculate_density_array(), calculate_density_derivs_array(), calculate_density_scalar(), calculate_spec_vol_array(), calculate_spec_vol_scalar(), calculate_specific_vol_derivs(), and eos_init().

eos_unesco_string

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

private

A string for specifying the equation of state.

Definition at line 119 of file MOM_EOS.F90.

character*(10), parameter :: EOS_UNESCO_STRING = "UNESCO" !< A string for specifying the equation of state

Referenced by eos_init().

eos_wright

integer, parameter, public mom_eos::eos_wright = 3

A named integer specifying an equation of state.

Definition at line 114 of file MOM_EOS.F90.

integer, parameter, public :: EOS_WRIGHT = 3 !< A named integer specifying an equation of state

Referenced by calculate_compress_array(), calculate_density_array(), calculate_density_derivs_array(), calculate_density_derivs_scalar(), calculate_density_scalar(), calculate_density_second_derivs_array(), calculate_density_second_derivs_scalar(), calculate_spec_vol_array(), calculate_spec_vol_scalar(), calculate_specific_vol_derivs(), eos_init(), int_density_dz(), and int_specific_vol_dp().

eos_wright_string

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

private

A string for specifying the equation of state.

Definition at line 120 of file MOM_EOS.F90.

character*(10), parameter :: EOS_WRIGHT_STRING = "WRIGHT" !< A string for specifying the equation of state

Referenced by eos_init().

tfreeze_default

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

private

The default freezing point expression.

Definition at line 132 of file MOM_EOS.F90.

character*(10), parameter :: TFREEZE_DEFAULT = TFREEZE_LINEAR_STRING !< The default freezing point expression

Referenced by eos_init().

tfreeze_linear

integer, parameter mom_eos::tfreeze_linear = 1

private

A named integer specifying a freezing point expression.

Definition at line 125 of file MOM_EOS.F90.

integer, parameter :: TFREEZE_LINEAR = 1  !< A named integer specifying a freezing point expression

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

tfreeze_linear_string

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

private

A string for specifying the freezing point expression.

Definition at line 128 of file MOM_EOS.F90.

character*(10), parameter :: TFREEZE_LINEAR_STRING = "LINEAR" !< A string for specifying the freezing point expression

Referenced by eos_init().

tfreeze_millero

integer, parameter mom_eos::tfreeze_millero = 2

private

A named integer specifying a freezing point expression.

Definition at line 126 of file MOM_EOS.F90.

integer, parameter :: TFREEZE_MILLERO = 2 !< A named integer specifying a freezing point expression

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

tfreeze_millero_string

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

private

A string for specifying freezing point expression.

Definition at line 129 of file MOM_EOS.F90.

character*(10), parameter :: TFREEZE_MILLERO_STRING = "MILLERO_78" !< A string for specifying

Referenced by eos_init().

tfreeze_teos10

integer, parameter mom_eos::tfreeze_teos10 = 3

private

A named integer specifying a freezing point expression.

Definition at line 127 of file MOM_EOS.F90.

integer, parameter :: TFREEZE_TEOS10 = 3  !< A named integer specifying a freezing point expression

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

tfreeze_teos10_string

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

private

A string for specifying the freezing point expression.

Definition at line 131 of file MOM_EOS.F90.

character*(10), parameter :: TFREEZE_TEOS10_STRING = "TEOS10" !< A string for specifying the freezing point expression

Referenced by eos_init().