diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 27484aa536..12b87ebc64 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -428,7 +428,7 @@ subroutine calculate_stanley_density_1d(T, S, pressure, Tvar, TScov, Svar, rho, real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1] real :: T2_scale ! A factor to convert temperature variance to units of degC2 [degC2 C-2 ~> 1] real :: S2_scale ! A factor to convert salinity variance to units of ppt2 [ppt2 S-2 ~> 1] - real :: TS_scale ! A factor to convert temperture-salinity covariance to units of + real :: TS_scale ! A factor to convert temperature-salinity covariance to units of ! degC ppt [degC ppt C-1 S-1 ~> 1] real :: rho_reference ! rho_ref converted to [kg m-3] real, dimension(size(rho)) :: pres ! Pressure converted to [Pa] @@ -1023,7 +1023,7 @@ subroutine calculate_density_second_derivs_1d(T, S, pressure, drho_dS_dS, drho_d end subroutine calculate_density_second_derivs_1d -!> Calls the appropriate subroutine to calculate density second derivatives for scalar nputs. +!> Calls the appropriate subroutine to calculate density second derivatives for scalar inputs. 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) real, intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -1266,7 +1266,7 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) type(EOS_type), intent(in) :: EOS !< Equation of state structure ! Local variables - ! These arrays use the same units as their counterparts in calcluate_compress_1d. + ! These arrays use the same units as their counterparts in calculate_compress_1d. real, dimension(1) :: pa ! Pressure in a size-1 1d array [R L2 T-2 ~> Pa] real, dimension(1) :: Ta ! Temperature in a size-1 1d array [C ~> degC] real, dimension(1) :: Sa ! Salinity in a size-1 1d array [S ~> ppt] @@ -1734,7 +1734,7 @@ subroutine abs_saln_to_prac_saln(S, prSaln, EOS, dom, scale) ! Local variables real, dimension(size(S)) :: Sa ! Salinity converted to [ppt] - real :: S_scale ! A factor to convert practical salnity from ppt to the desired units [S ppt-1 ~> 1] + real :: S_scale ! A factor to convert practical salinity from ppt to the desired units [S ppt-1 ~> 1] integer :: i, is, ie if (present(dom)) then diff --git a/src/equation_of_state/MOM_EOS_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 index 476fda6b70..de4a715489 100644 --- a/src/equation_of_state/MOM_EOS_NEMO.F90 +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -35,7 +35,7 @@ module MOM_EOS_NEMO module procedure calculate_density_derivs_scalar_nemo, calculate_density_derivs_array_nemo end interface calculate_density_derivs_nemo -real, parameter :: Pa2db = 1.e-4 !< Conversion factor between Pa and dbar +real, parameter :: Pa2db = 1.e-4 !< Conversion factor between Pa and dbar [Pa dbar-1] !>@{ Parameters in the NEMO equation of state real, parameter :: rdeltaS = 32. real, parameter :: r1_S0 = 0.875/35.16504 @@ -184,8 +184,10 @@ subroutine calculate_density_scalar_nemo(T, S, pressure, rho, rho_ref) real, intent(out) :: rho !< In situ density [kg m-3]. real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. - real, dimension(1) :: T0, S0, pressure0 - real, dimension(1) :: rho0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the conservative temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the absolute salinity [g kg-1] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: rho0 ! A 1-d array with a copy of the density [kg m-3] T0(1) = T S0(1) = S @@ -345,8 +347,13 @@ subroutine calculate_density_derivs_scalar_nemo(T, S, pressure, drho_dt, drho_ds real, intent(out) :: drho_dS !< The partial derivative of density with salinity, !! in [kg m-3 ppt-1]. ! Local variables - real, dimension(1) :: T0, S0, pressure0 - real, dimension(1) :: drdt0, drds0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the conservative temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the absolute salinity [g kg-1] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: drdt0 ! A 1-d array with a copy of the derivative of density + ! with potential temperature [kg m-3 degC-1] + real, dimension(1) :: drds0 ! A 1-d array with a copy of the derivative of density + ! with salinity [kg m-3 ppt-1] T0(1) = T S0(1) = S @@ -358,12 +365,12 @@ subroutine calculate_density_derivs_scalar_nemo(T, S, pressure, drho_dt, drho_ds end subroutine calculate_density_derivs_scalar_nemo !> Compute the in situ density of sea water (rho in [kg m-3]) and the compressibility -!! (drho/dp = C_sound^-2, stored as drho_dp [s2 m-2]) from absolute salinity -!! (sal in g/kg), conservative temperature (T [degC]), and pressure [Pa], using the expressions +!! (drho/dp = C_sound^-2, stored as drho_dp [s2 m-2]) from absolute salinity (sal [g kg-1]), +!! conservative temperature (T [degC]), and pressure [Pa], using the expressions !! derived for use with NEMO. subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) real, intent(in), dimension(:) :: T !< Conservative temperature [degC]. - real, intent(in), dimension(:) :: S !< Absolute salinity [g/kg]. + real, intent(in), dimension(:) :: S !< Absolute salinity [g kg-1]. real, intent(in), dimension(:) :: pressure !< pressure [Pa]. real, intent(out), dimension(:) :: rho !< In situ density [kg m-3]. real, intent(out), dimension(:) :: drho_dp !< The partial derivative of density with pressure @@ -373,7 +380,9 @@ subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: zs,zt,zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] integer :: j call calculate_density_array_nemo(T, S, pressure, rho, start, npts) @@ -384,7 +393,7 @@ subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) do j=start,start+npts-1 !Conversions zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potential temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j)) enddo diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 index bbe9982b6f..4c7483c068 100644 --- a/src/equation_of_state/MOM_EOS_TEOS10.F90 +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -48,14 +48,13 @@ module MOM_EOS_TEOS10 module procedure calculate_density_second_derivs_scalar_teos10, calculate_density_second_derivs_array_teos10 end interface calculate_density_second_derivs_teos10 -real, parameter :: Pa2db = 1.e-4 !< The conversion factor from Pa to dbar. +real, parameter :: Pa2db = 1.e-4 !< The conversion factor from Pa to dbar [dbar Pa-1] contains -!> This subroutine computes the in situ density of sea water (rho in -!! [kg m-3]) from absolute salinity (S [g kg-1]), conservative temperature -!! (T [degC]), and pressure [Pa]. It uses the expression from the -!! TEOS10 website. +!> This subroutine computes the in situ density of sea water (rho in [kg m-3]) +!! from absolute salinity (S [g kg-1]), conservative temperature (T [degC]), +!! and pressure [Pa]. It uses the expression from the TEOS10 website. subroutine calculate_density_scalar_teos10(T, S, pressure, rho, rho_ref) real, intent(in) :: T !< Conservative temperature [degC]. real, intent(in) :: S !< Absolute salinity [g kg-1]. @@ -64,8 +63,10 @@ subroutine calculate_density_scalar_teos10(T, S, pressure, rho, rho_ref) real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. ! Local variables - real, dimension(1) :: T0, S0, pressure0 - real, dimension(1) :: rho0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the conservative temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the absolute salinity [g kg-1] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: rho0 ! A 1-d array with a copy of the density [kg m-3] T0(1) = T S0(1) = S @@ -76,9 +77,9 @@ subroutine calculate_density_scalar_teos10(T, S, pressure, rho, rho_ref) end subroutine calculate_density_scalar_teos10 -!> This subroutine computes the in situ density of sea water (rho in -!! [kg m-3]) from absolute salinity (S [g kg-1]), conservative temperature -!! (T [degC]), and pressure [Pa]. It uses the expression from the +!> This subroutine computes the in situ density of sea water (rho in [kg m-3]) +!! from absolute salinity (S [g kg-1]), conservative temperature (T [degC]), +!! and pressure [Pa]. It uses the expression from the !! TEOS10 website. subroutine calculate_density_array_teos10(T, S, pressure, rho, start, npts, rho_ref) real, dimension(:), intent(in) :: T !< Conservative temperature [degC]. @@ -90,13 +91,15 @@ subroutine calculate_density_array_teos10(T, S, pressure, rho, start, npts, rho_ real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. ! Local variables - real :: zs, zt, zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] integer :: j do j=start,start+npts-1 !Conversions zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potential temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar if (S(j) < -1.0e-10) then !Can we assume safely that this is a missing value? @@ -120,7 +123,10 @@ subroutine calculate_spec_vol_scalar_teos10(T, S, pressure, specvol, spv_ref) real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. ! Local variables - real, dimension(1) :: T0, S0, pressure0, spv0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the conservative temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the absolute salinity [g kg-1] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: spv0 ! A 1-d array with a copy of the specific volume [m3 kg-1] T0(1) = T ; S0(1) = S ; pressure0(1) = pressure @@ -134,8 +140,7 @@ end subroutine calculate_spec_vol_scalar_teos10 !! and pressure [Pa], using the TEOS10 equation of state. !! If spv_ref is present, specvol is an anomaly from spv_ref. subroutine calculate_spec_vol_array_teos10(T, S, pressure, specvol, start, npts, spv_ref) - real, dimension(:), intent(in) :: T !< Conservative temperature relative to the surface - !! [degC]. + real, dimension(:), intent(in) :: T !< Conservative temperature [degC]. real, dimension(:), intent(in) :: S !< salinity [g kg-1]. real, dimension(:), intent(in) :: pressure !< pressure [Pa]. real, dimension(:), intent(out) :: specvol !< in situ specific volume [m3 kg-1]. @@ -144,13 +149,15 @@ subroutine calculate_spec_vol_array_teos10(T, S, pressure, specvol, start, npts, real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. ! Local variables - real :: zs, zt, zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] integer :: j do j=start,start+npts-1 !Conversions zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potential temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar if (S(j) < -1.0e-10) then @@ -177,15 +184,17 @@ subroutine calculate_density_derivs_array_teos10(T, S, pressure, drho_dT, drho_d integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: zs, zt, zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] integer :: j do j=start,start+npts-1 !Conversions zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potential temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar - if (S(j) < -1.0e-10) then ; !Can we assume safely that this is a missing value? + if (S(j) < -1.0e-10) then !Can we assume safely that this is a missing value? drho_dT(j) = 0.0 ; drho_dS(j) = 0.0 else call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS(j), drho_dct=drho_dT(j)) @@ -206,10 +215,13 @@ subroutine calculate_density_derivs_scalar_teos10(T, S, pressure, drho_dT, drho_ !! [kg m-3 (g/kg)-1]. ! Local variables - real :: zs, zt, zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] + !Conversions zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity - zt = T !gsw_ct_from_pt(S,T) !Convert potantial temp to conservative temp + zt = T !gsw_ct_from_pt(S,T) !Convert potential temp to conservative temp zp = pressure* Pa2db !Convert pressure from Pascal to decibar if (S < -1.0e-10) return !Can we assume safely that this is a missing value? call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS, drho_dct=drho_dT) @@ -229,15 +241,17 @@ subroutine calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: zs, zt, zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] integer :: j do j=start,start+npts-1 !Conversions zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potential temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar - if (S(j) < -1.0e-10) then ; !Can we assume safely that this is a missing value? + if (S(j) < -1.0e-10) then !Can we assume safely that this is a missing value? dSV_dT(j) = 0.0 ; dSV_dS(j) = 0.0 else call gsw_specvol_first_derivatives(zs,zt,zp, v_sa=dSV_dS(j), v_ct=dSV_dT(j)) @@ -252,18 +266,25 @@ subroutine calculate_density_second_derivs_scalar_teos10(T, S, pressure, drho_dS real, intent(in) :: T !< Conservative temperature [degC] real, intent(in) :: S !< Absolute Salinity [g kg-1] real, intent(in) :: pressure !< pressure [Pa]. - real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S - real, intent(out) :: drho_dS_dT !< Partial derivative of beta with resepct to T - real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T - real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure - real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure + real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect + !! to S [kg m-3 (g/kg)-2] + real, intent(out) :: drho_dS_dT !< Partial derivative of beta with respect + !! to T [kg m-3 (g/kg)-1 degC-1] + real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect + !! to T [kg m-3 degC-2] + real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect + !! to pressure [kg m-3 (g/kg)-1 Pa-1] = [s2 m-2 (g/kg)-1] + real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect + !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] ! Local variables - real :: zs, zt, zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] !Conversions zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity - zt = T !gsw_ct_from_pt(S,T) !Convert potantial temp to conservative temp + zt = T !gsw_ct_from_pt(S,T) !Convert potential temp to conservative temp zp = pressure* Pa2db !Convert pressure from Pascal to decibar if (S < -1.0e-10) return !Can we assume safely that this is a missing value? call gsw_rho_second_derivatives(zs, zt, zp, rho_sa_sa=drho_dS_dS, rho_sa_ct=drho_dS_dT, & @@ -277,24 +298,31 @@ subroutine calculate_density_second_derivs_array_teos10(T, S, pressure, drho_dS_ real, dimension(:), intent(in) :: T !< Conservative temperature [degC] real, dimension(:), intent(in) :: S !< Absolute Salinity [g kg-1] real, dimension(:), intent(in) :: pressure !< pressure [Pa]. - real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S - real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with resepct to T - real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T - real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure - real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure + real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect + !! to S [kg m-3 (g/kg)-2] + real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with respect + !! to T [kg m-3 (g/kg)-1 degC-1] + real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect + !! to T [kg m-3 degC-2] + real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect + !! to pressure [kg m-3 (g/kg)-1 Pa-1] = [s2 m-2 (g/kg)-1] + real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect + !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] integer, intent(in) :: start !< The starting point in the arrays. integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: zs, zt, zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] integer :: j do j=start,start+npts-1 !Conversions zs = S(j) !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity - zt = T(j) !gsw_ct_from_pt(S,T) !Convert potantial temp to conservative temp + zt = T(j) !gsw_ct_from_pt(S,T) !Convert potential temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar - if (S(j) < -1.0e-10) then ; !Can we assume safely that this is a missing value? + if (S(j) < -1.0e-10) then !Can we assume safely that this is a missing value? drho_dS_dS(j) = 0.0 ; drho_dS_dT(j) = 0.0 ; drho_dT_dT(j) = 0.0 drho_dS_dP(j) = 0.0 ; drho_dT_dP(j) = 0.0 else @@ -307,7 +335,7 @@ end subroutine calculate_density_second_derivs_array_teos10 !> This subroutine computes the in situ density of sea water (rho in !! [kg m-3]) and the compressibility (drho/dp = C_sound^-2) -!! (drho_dp [s2 m-2]) from absolute salinity (sal in g/kg), +!! (drho_dp [s2 m-2]) from absolute salinity (sal [g kg-1]), !! conservative temperature (T [degC]), and pressure [Pa]. It uses the !! subroutines from TEOS10 website subroutine calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts) @@ -322,15 +350,17 @@ subroutine calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts) integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: zs,zt,zp + real :: zs ! Absolute salinity [g kg-1] + real :: zt ! Conservative temperature [degC] + real :: zp ! Pressure converted to decibars [dbar] integer :: j do j=start,start+npts-1 !Conversions zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potential temp to conservative temp zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar - if (S(j) < -1.0e-10) then ; !Can we assume safely that this is a missing value? + if (S(j) < -1.0e-10) then !Can we assume safely that this is a missing value? rho(j) = 1000.0 ; drho_dp(j) = 0.0 else rho(j) = gsw_rho(zs,zt,zp) diff --git a/src/equation_of_state/MOM_EOS_UNESCO.F90 b/src/equation_of_state/MOM_EOS_UNESCO.F90 index a296cfc382..59ebb92c7a 100644 --- a/src/equation_of_state/MOM_EOS_UNESCO.F90 +++ b/src/equation_of_state/MOM_EOS_UNESCO.F90 @@ -17,45 +17,80 @@ module MOM_EOS_UNESCO public calculate_density_scalar_UNESCO, calculate_density_array_UNESCO !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to -!! a reference density, from salinity [PSU], potential temperature [degC], and pressure [Pa], -!! using the UNESCO (1981) equation of state. +!! a reference density, from salinity [PSU], potential temperature [degC] and pressure [Pa], +!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995). interface calculate_density_UNESCO module procedure calculate_density_scalar_UNESCO, calculate_density_array_UNESCO end interface calculate_density_UNESCO !> Compute the in situ specific volume of sea water (in [m3 kg-1]), or an anomaly with respect !! to a reference specific volume, from salinity [PSU], potential temperature [degC], and -!! pressure [Pa], using the UNESCO (1981) equation of state. +!! pressure [Pa], using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995). interface calculate_spec_vol_UNESCO module procedure calculate_spec_vol_scalar_UNESCO, calculate_spec_vol_array_UNESCO end interface calculate_spec_vol_UNESCO -!>@{ Parameters in the UNESCO equation of state -! The following constants are used to calculate rho0. The notation -! is Rab for the contribution to rho0 from T^aS^b. -real, parameter :: R00 = 999.842594, R10 = 6.793952e-2, R20 = -9.095290e-3, & - R30 = 1.001685e-4, R40 = -1.120083e-6, R50 = 6.536332e-9, R01 = 0.824493, & - R11 = -4.0899e-3, R21 = 7.6438e-5, R31 = -8.2467e-7, R41 = 5.3875e-9, & - R032 = -5.72466e-3, R132 = 1.0227e-4, R232 = -1.6546e-6, R02 = 4.8314e-4 - -! The following constants are used to calculate the secant bulk mod- -! ulus. The notation here is Sab for terms proportional to T^a*S^b, -! Spab for terms proportional to p*T^a*S^b, and SPab for terms +!>@{ Parameters in the UNESCO equation of state, as published in appendix A3 of Gill, 1982. +! The following constants are used to calculate rho0, the density of seawater at 1 +! atmosphere pressure. The notation is Rab for the contribution to rho0 from T^a*S^b. +real, parameter :: R00 = 999.842594 ! A coefficient in the fit for rho0 [kg m-3] +real, parameter :: R10 = 6.793952e-2 ! A coefficient in the fit for rho0 [kg m-3 degC-1] +real, parameter :: R20 = -9.095290e-3 ! A coefficient in the fit for rho0 [kg m-3 degC-2] +real, parameter :: R30 = 1.001685e-4 ! A coefficient in the fit for rho0 [kg m-3 degC-3] +real, parameter :: R40 = -1.120083e-6 ! A coefficient in the fit for rho0 [kg m-3 degC-4] +real, parameter :: R50 = 6.536332e-9 ! A coefficient in the fit for rho0 [kg m-3 degC-5] +real, parameter :: R01 = 0.824493 ! A coefficient in the fit for rho0 [kg m-3 PSU-1] +real, parameter :: R11 = -4.0899e-3 ! A coefficient in the fit for rho0 [kg m-3 degC-1 PSU-1] +real, parameter :: R21 = 7.6438e-5 ! A coefficient in the fit for rho0 [kg m-3 degC-2 PSU-1] +real, parameter :: R31 = -8.2467e-7 ! A coefficient in the fit for rho0 [kg m-3 degC-3 PSU-1] +real, parameter :: R41 = 5.3875e-9 ! A coefficient in the fit for rho0 [kg m-3 degC-4 PSU-1] +real, parameter :: R032 = -5.72466e-3 ! A coefficient in the fit for rho0 [kg m-3 PSU-3/2] +real, parameter :: R132 = 1.0227e-4 ! A coefficient in the fit for rho0 [kg m-3 PSU-3/2] +real, parameter :: R232 = -1.6546e-6 ! A coefficient in the fit for rho0 [kg m-3 PSU-3/2] +real, parameter :: R02 = 4.8314e-4 ! A coefficient in the fit for rho0 [kg m-3 PSU-2] + +! The following constants are used to calculate the secant bulk modulus. +! The notation here is Sab for terms proportional to T^a*S^b, +! Spab for terms proportional to p*T^a*S^b, and SP0ab for terms ! proportional to p^2*T^a*S^b. -real, parameter :: S00 = 1.965933e4, S10 = 1.444304e2, S20 = -1.706103, & - S30 = 9.648704e-3, S40 = -4.190253e-5, S01 = 52.84855, S11 = -3.101089e-1, & - S21 = 6.283263e-3, S31 = -5.084188e-5, S032 = 3.886640e-1, S132 = 9.085835e-3, & - S232 = -4.619924e-4, Sp00 = 3.186519, Sp10 = 2.212276e-2, Sp20 = -2.984642e-4, & - Sp30 = 1.956415e-6, Sp01 = 6.704388e-3, Sp11 = -1.847318e-4, Sp21 = 2.059331e-7, & - Sp032 = 1.480266e-4, SP000 = 2.102898e-4, SP010 = -1.202016e-5, SP020 = 1.394680e-7, & - SP001 = -2.040237e-6, SP011 = 6.128773e-8, SP021 = 6.207323e-10 +! Note that these values differ from those in Appendix A of Gill (1982) because the expressions +! from Jackett and MacDougall (1995) use potential temperature, rather than in situ temperature. +real, parameter :: S00 = 1.965933e4 ! A coefficient in the secant bulk modulus fit [bar] +real, parameter :: S10 = 1.444304e2 ! A coefficient in the secant bulk modulus fit [bar degC-1] +real, parameter :: S20 = -1.706103 ! A coefficient in the secant bulk modulus fit [bar degC-2] +real, parameter :: S30 = 9.648704e-3 ! A coefficient in the secant bulk modulus fit [bar degC-3] +real, parameter :: S40 = -4.190253e-5 ! A coefficient in the secant bulk modulus fit [bar degC-4] +real, parameter :: S01 = 52.84855 ! A coefficient in the secant bulk modulus fit [bar PSU-1] +real, parameter :: S11 = -3.101089e-1 ! A coefficient in the secant bulk modulus fit [bar degC-1 PSU-1] +real, parameter :: S21 = 6.283263e-3 ! A coefficient in the secant bulk modulus fit [bar degC-2 PSU-1] +real, parameter :: S31 = -5.084188e-5 ! A coefficient in the secant bulk modulus fit [bar degC-3 PSU-1] +real, parameter :: S032 = 3.886640e-1 ! A coefficient in the secant bulk modulus fit [bar PSU-3/2] +real, parameter :: S132 = 9.085835e-3 ! A coefficient in the secant bulk modulus fit [bar degC-1 PSU-3/2] +real, parameter :: S232 = -4.619924e-4 ! A coefficient in the secant bulk modulus fit [bar degC-2 PSU-3/2] + +real, parameter :: Sp00 = 3.186519 ! A coefficient in the secant bulk modulus fit [nondim] +real, parameter :: Sp10 = 2.212276e-2 ! A coefficient in the secant bulk modulus fit [degC-1] +real, parameter :: Sp20 = -2.984642e-4 ! A coefficient in the secant bulk modulus fit [degC-2] +real, parameter :: Sp30 = 1.956415e-6 ! A coefficient in the secant bulk modulus fit [degC-3] +real, parameter :: Sp01 = 6.704388e-3 ! A coefficient in the secant bulk modulus fit [PSU-1] +real, parameter :: Sp11 = -1.847318e-4 ! A coefficient in the secant bulk modulus fit [degC-1 PSU-1] +real, parameter :: Sp21 = 2.059331e-7 ! A coefficient in the secant bulk modulus fit [degC-2 PSU-1] +real, parameter :: Sp032 = 1.480266e-4 ! A coefficient in the secant bulk modulus fit [PSU-3/2] + +real, parameter :: SP000 = 2.102898e-4 ! A coefficient in the secant bulk modulus fit [bar-1] +real, parameter :: SP010 = -1.202016e-5 ! A coefficient in the secant bulk modulus fit [bar-1 degC-1] +real, parameter :: SP020 = 1.394680e-7 ! A coefficient in the secant bulk modulus fit [bar-1 degC-2] +real, parameter :: SP001 = -2.040237e-6 ! A coefficient in the secant bulk modulus fit [bar-1 PSU-1] +real, parameter :: SP011 = 6.128773e-8 ! A coefficient in the secant bulk modulus fit [bar-1 degC-1 PSU-1] +real, parameter :: SP021 = 6.207323e-10 ! A coefficient in the secant bulk modulus fit [bar-1 degC-1 PSU-2] !>@} contains -!> This subroutine computes the in situ density of sea water (rho in -!! [kg m-3]) from salinity (S [PSU]), potential temperature -!! (T [degC]), and pressure [Pa], using the UNESCO (1981) equation of state. +!> This subroutine computes the in situ density of sea water (rho in [kg m-3]) +!! from salinity (S [PSU]), potential temperature (T [degC]), and pressure [Pa], +!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995). +!! If rho_ref is present, rho is an anomaly from rho_ref. subroutine calculate_density_scalar_UNESCO(T, S, pressure, rho, rho_ref) real, intent(in) :: T !< Potential temperature relative to the surface [degC]. real, intent(in) :: S !< Salinity [PSU]. @@ -64,8 +99,10 @@ subroutine calculate_density_scalar_UNESCO(T, S, pressure, rho, rho_ref) real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. ! Local variables - real, dimension(1) :: T0, S0, pressure0 - real, dimension(1) :: rho0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the potential temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the salinity [PSU] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: rho0 ! A 1-d array with a copy of the in situ density [kg m-3] T0(1) = T S0(1) = S @@ -76,9 +113,10 @@ subroutine calculate_density_scalar_UNESCO(T, S, pressure, rho, rho_ref) end subroutine calculate_density_scalar_UNESCO -!> This subroutine computes the in situ density of sea water (rho in -!! [kg m-3]) from salinity (S [PSU]), potential temperature -!! (T [degC]), and pressure [Pa], using the UNESCO (1981) equation of state. +!> This subroutine computes the in situ density of sea water (rho in [kg m-3]) +!! from salinity (S [PSU]), potential temperature (T [degC]) and pressure [Pa], +!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995). +!! If rho_ref is present, rho is an anomaly from rho_ref. subroutine calculate_density_array_UNESCO(T, S, pressure, rho, start, npts, rho_ref) real, dimension(:), intent(in) :: T !< potential temperature relative to the surface [degC]. real, dimension(:), intent(in) :: S !< salinity [PSU]. @@ -89,8 +127,12 @@ subroutine calculate_density_array_UNESCO(T, S, pressure, rho, start, npts, rho_ real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. ! Local variables - real :: t_local, t2, t3, t4, t5 ! Temperature to the 1st - 5th power [degC^n]. - real :: s_local, s32, s2 ! Salinity to the 1st, 3/2, & 2nd power [PSU^n]. + real :: t_local ! A copy of the temperature at a point [degC] + real :: t2, t3 ! Temperature squared [degC2] and cubed [degC3] + real :: t4, t5 ! Temperature to the 4th power [degC4] and 5th power [degC5] + real :: s_local ! A copy of the salinity at a point [PSU] + real :: s32 ! The square root of salinity cubed [PSU3/2] + real :: s2 ! Salinity squared [PSU2]. real :: p1, p2 ! Pressure (in bars) to the 1st and 2nd power [bar] and [bar2]. real :: rho0 ! Density at 1 bar pressure [kg m-3]. real :: sig0 ! The anomaly of rho0 from R00 [kg m-3]. @@ -103,9 +145,9 @@ subroutine calculate_density_array_UNESCO(T, S, pressure, rho, start, npts, rho_ cycle endif - p1 = pressure(j)*1.0e-5; p2 = p1*p1 - t_local = T(j); t2 = t_local*t_local; t3 = t_local*t2; t4 = t2*t2; t5 = t3*t2 - s_local = S(j); s2 = s_local*s_local; s32 = s_local*sqrt(s_local) + p1 = pressure(j)*1.0e-5 ; p2 = p1*p1 + t_local = T(j) ; t2 = t_local*t_local ; t3 = t_local*t2 ; t4 = t2*t2 ; t5 = t3*t2 + s_local = S(j) ; s2 = s_local*s_local ; s32 = s_local*sqrt(s_local) ! Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ). @@ -130,9 +172,9 @@ subroutine calculate_density_array_UNESCO(T, S, pressure, rho, start, npts, rho_ enddo end subroutine calculate_density_array_UNESCO -!> This subroutine computes the in situ specific volume of sea water (specvol in -!! [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC]) -!! and pressure [Pa], using the UNESCO (1981) equation of state. +!> This subroutine computes the in situ specific volume of sea water (specvol in [m3 kg-1]) +!! from salinity (S [PSU]), potential temperature (T [degC]) and pressure [Pa], +!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995). !! If spv_ref is present, specvol is an anomaly from spv_ref. subroutine calculate_spec_vol_scalar_UNESCO(T, S, pressure, specvol, spv_ref) real, intent(in) :: T !< potential temperature relative to the surface @@ -143,7 +185,10 @@ subroutine calculate_spec_vol_scalar_UNESCO(T, S, pressure, specvol, spv_ref) real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. ! Local variables - real, dimension(1) :: T0, S0, pressure0, spv0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the potential temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the salinity [PSU] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: spv0 ! A 1-d array with a copy of the specific volume [m3 kg-1] T0(1) = T ; S0(1) = S ; pressure0(1) = pressure @@ -151,9 +196,9 @@ subroutine calculate_spec_vol_scalar_UNESCO(T, S, pressure, specvol, spv_ref) specvol = spv0(1) end subroutine calculate_spec_vol_scalar_UNESCO -!> This subroutine computes the in situ specific volume of sea water (specvol in -!! [m3 kg-1]) from salinity (S [PSU]), potential temperature (T [degC]) -!! and pressure [Pa], using the UNESCO (1981) equation of state. +!> This subroutine computes the in situ specific volume of sea water (specvol in [m3 kg-1]) +!! from salinity (S [PSU]), potential temperature (T [degC]) and pressure [Pa], +!! using the UNESCO (1981) equation of state, as refit by Jackett and McDougall (1995). !! If spv_ref is present, specvol is an anomaly from spv_ref. subroutine calculate_spec_vol_array_UNESCO(T, S, pressure, specvol, start, npts, spv_ref) real, dimension(:), intent(in) :: T !< potential temperature relative to the surface @@ -166,8 +211,12 @@ subroutine calculate_spec_vol_array_UNESCO(T, S, pressure, specvol, start, npts, real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. ! Local variables - real :: t_local, t2, t3, t4, t5 ! Temperature to the 1st - 5th power [degC^n]. - real :: s_local, s32, s2 ! Salinity to the 1st, 3/2, & 2nd power [PSU^n]. + real :: t_local ! A copy of the temperature at a point [degC] + real :: t2, t3 ! Temperature squared [degC2] and cubed [degC3] + real :: t4, t5 ! Temperature to the 4th power [degC4] and 5th power [degC5] + real :: s_local ! A copy of the salinity at a point [PSU] + real :: s32 ! The square root of salinity cubed [PSU3/2] + real :: s2 ! Salinity squared [PSU2]. real :: p1, p2 ! Pressure (in bars) to the 1st and 2nd power [bar] and [bar2]. real :: rho0 ! Density at 1 bar pressure [kg m-3]. real :: ks ! The secant bulk modulus [bar]. @@ -180,9 +229,9 @@ subroutine calculate_spec_vol_array_UNESCO(T, S, pressure, specvol, start, npts, cycle endif - p1 = pressure(j)*1.0e-5; p2 = p1*p1 - t_local = T(j); t2 = t_local*t_local; t3 = t_local*t2; t4 = t2*t2; t5 = t3*t2 - s_local = S(j); s2 = s_local*s_local; s32 = s_local*sqrt(s_local) + p1 = pressure(j)*1.0e-5 ; p2 = p1*p1 + t_local = T(j) ; t2 = t_local*t_local ; t3 = t_local*t2 ; t4 = t2*t2 ; t5 = t3*t2 + s_local = S(j) ; s2 = s_local*s_local ; s32 = s_local*sqrt(s_local) ! Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ). @@ -222,8 +271,13 @@ subroutine calculate_density_derivs_UNESCO(T, S, pressure, drho_dT, drho_dS, sta integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: t_local, t2, t3, t4, t5 ! Temperature to the 1st - 5th power [degC^n]. - real :: s12, s_local, s32, s2 ! Salinity to the 1/2 - 2nd powers [PSU^n]. + real :: t_local ! A copy of the temperature at a point [degC] + real :: t2, t3 ! Temperature squared [degC2] and cubed [degC3] + real :: t4, t5 ! Temperature to the 4th power [degC4] and 5th power [degC5] + real :: s12 ! The square root of salinity [PSU1/2] + real :: s_local ! A copy of the salinity at a point [PSU] + real :: s32 ! The square root of salinity cubed [PSU3/2] + real :: s2 ! Salinity squared [PSU2]. real :: p1, p2 ! Pressure to the 1st & 2nd power [bar] and [bar2]. real :: rho0 ! Density at 1 bar pressure [kg m-3]. real :: ks ! The secant bulk modulus [bar]. @@ -240,9 +294,9 @@ subroutine calculate_density_derivs_UNESCO(T, S, pressure, drho_dT, drho_dS, sta cycle endif - p1 = pressure(j)*1.0e-5; p2 = p1*p1 - t_local = T(j); t2 = t_local*t_local; t3 = t_local*t2; t4 = t2*t2; t5 = t3*t2 - s_local = S(j); s2 = s_local*s_local; s12 = sqrt(s_local); s32 = s_local*s12 + p1 = pressure(j)*1.0e-5 ; p2 = p1*p1 + t_local = T(j) ; t2 = t_local*t_local ; t3 = t_local*t2 ; t4 = t2*t2 ; t5 = t3*t2 + s_local = S(j) ; s2 = s_local*s_local ; s12 = sqrt(s_local) ; s32 = s_local*s12 ! compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ) @@ -293,14 +347,20 @@ subroutine calculate_compress_UNESCO(T, S, pressure, rho, drho_dp, start, npts) integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: t_local, t2, t3, t4, t5 ! Temperature to the 1st - 5th power [degC^n]. - real :: s_local, s32, s2 ! Salinity to the 1st, 3/2, & 2nd power [PSU^n]. - real :: p1, p2 ! Pressure to the 1st & 2nd power [bar] and [bar2]. - real :: rho0 ! Density at 1 bar pressure [kg m-3]. - real :: ks ! The secant bulk modulus [bar]. - real :: ks_0, ks_1, ks_2 - real :: dks_dp ! The derivative of the secant bulk modulus - ! with pressure, nondimensional. + real :: t_local ! A copy of the temperature at a point [degC] + real :: t2, t3 ! Temperature squared [degC2] and cubed [degC3] + real :: t4, t5 ! Temperature to the 4th power [degC4] and 5th power [degC5] + real :: s_local ! A copy of the salinity at a point [PSU] + real :: s32 ! The square root of salinity cubed [PSU3/2] + real :: s2 ! Salinity squared [PSU2]. + real :: p1, p2 ! Pressure to the 1st & 2nd power [bar] and [bar2]. + real :: rho0 ! Density at 1 bar pressure [kg m-3]. + real :: ks ! The secant bulk modulus [bar]. + real :: ks_0 ! The secant bulk modulus at zero pressure [bar]. + real :: ks_1 ! The derivative of the secant bulk modulus with pressure at zero pressure [nondim]. + real :: ks_2 ! The second derivative of the secant bulk modulus with pressure at zero pressure [nondim]. + real :: dks_dp ! The derivative of the secant bulk modulus + ! with pressure [nondim] integer :: j do j=start,start+npts-1 @@ -309,9 +369,9 @@ subroutine calculate_compress_UNESCO(T, S, pressure, rho, drho_dp, start, npts) cycle endif - p1 = pressure(j)*1.0e-5; p2 = p1*p1 - t_local = T(j); t2 = t_local*t_local; t3 = t_local*t2; t4 = t2*t2; t5 = t3*t2 - s_local = S(j); s2 = s_local*s_local; s32 = s_local*sqrt(s_local) + p1 = pressure(j)*1.0e-5 ; p2 = p1*p1 + t_local = T(j) ; t2 = t_local*t_local ; t3 = t_local*t2 ; t4 = t2*t2 ; t5 = t3*t2 + s_local = S(j) ; s2 = s_local*s_local ; s32 = s_local*sqrt(s_local) ! Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ). diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index c2e50287b2..77e0d17ff3 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -27,15 +27,17 @@ module MOM_EOS_Wright !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to -!! a reference density, from salinity (in psu), potential temperature (in deg C), and pressure [Pa], -!! using the expressions from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. +!! a reference density, from salinity in practical salinity units ([PSU]), potential +!! temperature (in degrees Celsius [degC]), and pressure [Pa], using the expressions from +!! Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. interface calculate_density_wright module procedure calculate_density_scalar_wright, calculate_density_array_wright end interface calculate_density_wright !> Compute the in situ specific volume of sea water (in [m3 kg-1]), or an anomaly with respect -!! to a reference specific volume, from salinity (in psu), potential temperature (in deg C), and -!! pressure [Pa], using the expressions from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. +!! to a reference specific volume, from salinity in practical salinity units ([PSU]), potential +!! temperature (in degrees Celsius [degC]), and pressure [Pa], using the expressions from +!! Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. interface calculate_spec_vol_wright module procedure calculate_spec_vol_scalar_wright, calculate_spec_vol_array_wright end interface calculate_spec_vol_wright @@ -64,11 +66,25 @@ module MOM_EOS_Wright ! Following are the values for the reduced range formula. -real, parameter :: a0 = 7.057924e-4, a1 = 3.480336e-7, a2 = -1.112733e-7 ! a0/a1 ~= 2028 ; a0/a2 ~= -6343 -real, parameter :: b0 = 5.790749e8, b1 = 3.516535e6, b2 = -4.002714e4 ! b0/b1 ~= 165 ; b0/b4 ~= 974 -real, parameter :: b3 = 2.084372e2, b4 = 5.944068e5, b5 = -9.643486e3 -real, parameter :: c0 = 1.704853e5, c1 = 7.904722e2, c2 = -7.984422 ! c0/c1 ~= 216 ; c0/c4 ~= -740 -real, parameter :: c3 = 5.140652e-2, c4 = -2.302158e2, c5 = -3.079464 + ! Note that a0/a1 ~= 2028 [degC] ; a0/a2 ~= -6343 [PSU] + ! b0/b1 ~= 165 [degC] ; b0/b4 ~= 974 [PSU] + ! c0/c1 ~= 216 [degC] ; c0/c4 ~= -740 [PSU] + ! and also that (as always) [Pa] = [kg m-1 s-2] +real, parameter :: a0 = 7.057924e-4 ! A parameter in the Wright alpha_0 fit [m3 kg-1] +real, parameter :: a1 = 3.480336e-7 ! A parameter in the Wright alpha_0 fit [m3 kg-1 degC-1] +real, parameter :: a2 = -1.112733e-7 ! A parameter in the Wright alpha_0 fit [m3 kg-1 PSU-1] +real, parameter :: b0 = 5.790749e8 ! A parameter in the Wright p_0 fit [Pa] +real, parameter :: b1 = 3.516535e6 ! A parameter in the Wright p_0 fit [Pa degC-1] +real, parameter :: b2 = -4.002714e4 ! A parameter in the Wright p_0 fit [Pa degC-2] +real, parameter :: b3 = 2.084372e2 ! A parameter in the Wright p_0 fit [Pa degC-3] +real, parameter :: b4 = 5.944068e5 ! A parameter in the Wright p_0 fit [Pa PSU-1] +real, parameter :: b5 = -9.643486e3 ! A parameter in the Wright p_0 fit [Pa degC-1 PSU-1] +real, parameter :: c0 = 1.704853e5 ! A parameter in the Wright lambda fit [m2 s-2] +real, parameter :: c1 = 7.904722e2 ! A parameter in the Wright lambda fit [m2 s-2 degC-1] +real, parameter :: c2 = -7.984422 ! A parameter in the Wright lambda fit [m2 s-2 degC-2] +real, parameter :: c3 = 5.140652e-2 ! A parameter in the Wright lambda fit [m2 s-2 degC-3] +real, parameter :: c4 = -2.302158e2 ! A parameter in the Wright lambda fit [m2 s-2 PSU-1] +real, parameter :: c5 = -3.079464 ! A parameter in the Wright lambda fit [m2 s-2 degC-1 PSU-1] !>@} contains @@ -86,13 +102,16 @@ subroutine calculate_density_scalar_wright(T, S, pressure, rho, rho_ref) ! *====================================================================* ! * This subroutine computes the in situ density of sea water (rho in * -! * [kg m-3]) from salinity (S [PSU]), potential temperature * -! * (T [degC]), and pressure [Pa]. It uses the expression from * +! * [kg m-3]) from salinity (S [PSU]), potential temperature * +! * (T [degC]), and pressure [Pa]. It uses the expression from * ! * Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. * ! * Coded by R. Hallberg, 7/00 * ! *====================================================================* - real, dimension(1) :: T0, S0, pressure0, rho0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the potential temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the salinity [PSU] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: rho0 ! A 1-d array with a copy of the density [kg m-3] T0(1) = T S0(1) = S @@ -118,8 +137,13 @@ subroutine calculate_density_array_wright(T, S, pressure, rho, start, npts, rho_ ! Original coded by R. Hallberg, 7/00, anomaly coded in 3/18. ! Local variables - real :: al0, p0, lambda - real :: al_TS, p_TSp, lam_TS, pa_000 + real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1] + real :: p0 ! The pressure offset in the Wright EOS [Pa] + real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2] + real :: al_TS ! The contributions of temperature and salinity to al0 [m3 kg-1] + real :: p_TSp ! A combination of the pressure and the temperature and salinity contributions to p0 [Pa] + real :: lam_TS ! The contributions of temperature and salinity to lambda [m2 s-2] + real :: pa_000 ! A corrected offset to the pressure, including contributions from rho_ref [Pa] integer :: j if (present(rho_ref)) pa_000 = (b0*(1.0 - a0*rho_ref) - rho_ref*c0) @@ -155,7 +179,10 @@ subroutine calculate_spec_vol_scalar_wright(T, S, pressure, specvol, spv_ref) real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. ! Local variables - real, dimension(1) :: T0, S0, pressure0, spv0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the potential temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the salinity [PSU] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: spv0 ! A 1-d array with a copy of the specific volume [m3 kg-1] T0(1) = T ; S0(1) = S ; pressure0(1) = pressure @@ -170,7 +197,7 @@ end subroutine calculate_spec_vol_scalar_wright !! If spv_ref is present, specvol is an anomaly from spv_ref. subroutine calculate_spec_vol_array_wright(T, S, pressure, specvol, start, npts, spv_ref) real, dimension(:), intent(in) :: T !< potential temperature relative to the - !! surface [degC]. + !! surface [degC]. real, dimension(:), intent(in) :: S !< salinity [PSU]. real, dimension(:), intent(in) :: pressure !< pressure [Pa]. real, dimension(:), intent(inout) :: specvol !< in situ specific volume [m3 kg-1]. @@ -179,7 +206,9 @@ subroutine calculate_spec_vol_array_wright(T, S, pressure, specvol, start, npts, real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. ! Local variables - real :: al0, p0, lambda + real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1] + real :: p0 ! The pressure offset in the Wright EOS [Pa] + real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2] integer :: j do j=start,start+npts-1 @@ -209,7 +238,10 @@ subroutine calculate_density_derivs_array_wright(T, S, pressure, drho_dT, drho_d integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: al0, p0, lambda, I_denom2 + real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1] + real :: p0 ! The pressure offset in the Wright EOS [Pa] + real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2] + real :: I_denom2 ! The inverse of the square of the denominator of density in the Wright EOS [s4 m-4] integer :: j do j=start,start+npts-1 @@ -241,8 +273,11 @@ subroutine calculate_density_derivs_scalar_wright(T, S, pressure, drho_dT, drho_ !! in [kg m-3 PSU-1]. ! Local variables needed to promote the input/output scalars to 1-element arrays - real, dimension(1) :: T0, S0, P0 - real, dimension(1) :: drdt0, drds0 + real, dimension(1) :: T0 ! A 1-d array with a copy of the temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the salinity [PSU] + real, dimension(1) :: p0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: drdt0 ! The derivative of density with temperature [kg m-3 degC-1] + real, dimension(1) :: drds0 ! The derivative of density with salinity [kg m-3 PSU-1] T0(1) = T S0(1) = S @@ -261,19 +296,28 @@ subroutine calculate_density_second_derivs_array_wright(T, S, P, drho_ds_ds, drh real, dimension(:), intent(in ) :: P !< Pressure [Pa] real, dimension(:), intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect !! to S [kg m-3 PSU-2] - real, dimension(:), intent(inout) :: drho_ds_dt !< Partial derivative of beta with respcct + real, dimension(:), intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect !! to T [kg m-3 PSU-1 degC-1] real, dimension(:), intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect !! to T [kg m-3 degC-2] real, dimension(:), intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect - !! to pressure [kg m-3 PSU-1 Pa-1] + !! to pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1] real, dimension(:), intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect - !! to pressure [kg m-3 degC-1 Pa-1] + !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] integer, intent(in ) :: start !< Starting index in T,S,P integer, intent(in ) :: npts !< Number of points to loop over ! Local variables - real :: z0, z1, z2, z3, z4, z5, z6 ,z7, z8, z9, z10, z11, z2_2, z2_3 + real :: z0, z1 ! Local work variables [Pa] + real :: z2, z4 ! Local work variables [m2 s-2] + real :: z3, z5 ! Local work variables [Pa degC-1] + real :: z6, z8 ! Local work variables [m2 s-2 degC-1] + real :: z7 ! A local work variable [m2 s-2 PSU-1] + real :: z9 ! A local work variable [m3 kg-1] + real :: z10 ! A local work variable [Pa PSU-1] + real :: z11 ! A local work variable [Pa m2 s-2 PSU-1] = [kg m s-4 PSU-1] + real :: z2_2 ! A local work variable [m4 s-4] + real :: z2_3 ! A local work variable [m6 s-6] integer :: j ! Based on the above expression with common terms factored, there probably exists a more numerically stable ! and/or efficient expression @@ -313,17 +357,26 @@ subroutine calculate_density_second_derivs_scalar_wright(T, S, P, drho_ds_ds, dr real, intent(in ) :: P !< pressure [Pa] real, intent( out) :: drho_ds_ds !< Partial derivative of beta with respect !! to S [kg m-3 PSU-2] - real, intent( out) :: drho_ds_dt !< Partial derivative of beta with respcct + real, intent( out) :: drho_ds_dt !< Partial derivative of beta with respect !! to T [kg m-3 PSU-1 degC-1] real, intent( out) :: drho_dt_dt !< Partial derivative of alpha with respect !! to T [kg m-3 degC-2] real, intent( out) :: drho_ds_dp !< Partial derivative of beta with respect - !! to pressure [kg m-3 PSU-1 Pa-1] + !! to pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1] real, intent( out) :: drho_dt_dp !< Partial derivative of alpha with respect - !! to pressure [kg m-3 degC-1 Pa-1] + !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] ! Local variables - real, dimension(1) :: T0, S0, P0 - real, dimension(1) :: drdsds, drdsdt, drdtdt, drdsdp, drdtdp + real, dimension(1) :: T0 ! A 1-d array with a copy of the temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the salinity [PSU] + real, dimension(1) :: p0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: drdsds ! The second derivative of density with salinity [kg m-3 PSU-2] + real, dimension(1) :: drdsdt ! The second derivative of density with salinity and + ! temperature [kg m-3 PSU-1 degC-1] + real, dimension(1) :: drdtdt ! The second derivative of density with temperature [kg m-3 degC-2] + real, dimension(1) :: drdsdp ! The second derivative of density with salinity and + ! pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1] + real, dimension(1) :: drdtdp ! The second derivative of density with temperature and + ! pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] T0(1) = T S0(1) = S @@ -346,12 +399,14 @@ subroutine calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start real, intent(inout), dimension(:) :: dSV_dT !< The partial derivative of specific volume with !! potential temperature [m3 kg-1 degC-1]. real, intent(inout), dimension(:) :: dSV_dS !< The partial derivative of specific volume with - !! salinity [m3 kg-1 / Pa]. + !! salinity [m3 kg-1 PSU-1]. integer, intent(in) :: start !< The starting point in the arrays. integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real :: p0, lambda, I_denom + real :: p0 ! The pressure offset in the Wright EOS [Pa] + real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2] + real :: I_denom ! The inverse of the denominator of specific volume in the Wright EOS [Pa-1] integer :: j do j=start,start+npts-1 @@ -370,11 +425,10 @@ subroutine calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start end subroutine calculate_specvol_derivs_wright -!> This subroutine computes the in situ density of sea water (rho in -!! [kg m-3]) and the compressibility (drho/dp = C_sound^-2) -!! (drho_dp [s2 m-2]) from salinity (sal in psu), potential -!! temperature (T [degC]), and pressure [Pa]. It uses the expressions -!! from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. +!> This subroutine computes the in situ density of sea water (rho in [kg m-3]) +!! and the compressibility (drho/dp = C_sound^-2) (drho_dp [s2 m-2]) from +!! salinity (sal [PSU]), potential temperature (T [degC]), and pressure [Pa]. +!! It uses the expressions from Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. !! Coded by R. Hallberg, 1/01 subroutine calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts) real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface [degC]. @@ -389,7 +443,10 @@ subroutine calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts) ! Coded by R. Hallberg, 1/01 ! Local variables - real :: al0, p0, lambda, I_denom + real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1] + real :: p0 ! The pressure offset in the Wright EOS [Pa] + real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2] + real :: I_denom ! The inverse of the denominator of density in the Wright EOS [s2 m-2] integer :: j do j=start,start+npts-1 @@ -421,7 +478,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & intent(in) :: z_b !< Height at the top of the layer [Z ~> m]. real, intent(in) :: rho_ref !< A mean density [R ~> kg m-3], that is subtracted !! out to reduce the magnitude of each of the integrals. - !! (The pressure is calucated as p~=-z*rho_0*G_e.) + !! (The pressure is calculated as p~=-z*rho_0*G_e.) real, intent(in) :: rho_0 !< Density [R ~> kg m-3], that is used !! to calculate the pressure (as p~=-z*rho_0*G_e) !! used in the equation of state. @@ -454,7 +511,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real, optional, intent(in) :: temp_scale !< A multiplicative factor by which to scale !! temperature into degC [degC C-1 ~> 1] real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure - !! into ppt [ppt S-1 ~> 1]. + !! into PSU [PSU S-1 ~> 1]. real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! Local variables @@ -488,20 +545,20 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & ! pres_scale [R L2 T-2 Pa-1 ~> 1]. real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1] - real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 ppt-1] + real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1] real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1] real :: b2s ! Partly rescaled version of b2 [Pa C-2 ~> Pa degC-2] real :: b3s ! Partly rescaled version of b3 [Pa C-3 ~> Pa degC-3] - real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa ppt-1] - real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 ppt-1] + real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa PSU-1] + real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 PSU-1] real :: c1s ! Partly rescaled version of c1 [m2 s-2 C-1 ~> m2 s-2 degC-1] real :: c2s ! Partly rescaled version of c2 [m2 s-2 C-2 ~> m2 s-2 degC-2] real :: c3s ! Partly rescaled version of c3 [m2 s-2 C-3 ~> m2 s-2 degC-3] - real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 ppt-1] - real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 ppt-1] + real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] + real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. - real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants. - real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] + real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m ! These array bounds work for the indexing convention of the input arrays, but @@ -716,7 +773,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & real, optional, intent(in) :: temp_scale !< A multiplicative factor by which to scale !! temperature into degC [degC C-1 ~> 1] real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure - !! into ppt [ppt S-1 ~> 1]. + !! into PSU [PSU S-1 ~> 1]. ! Local variables real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [R-1 ~> m3 kg-1] @@ -743,27 +800,27 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & real :: intp(5) ! The integrals of specific volume with pressure at the ! 5 sub-column locations [L2 T-2 ~> m2 s-2]. real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1] - real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 ppt-1] + real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1] real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1] real :: b2s ! Partly rescaled version of b2 [Pa C-2 ~> Pa degC-2] real :: b3s ! Partly rescaled version of b3 [Pa C-3 ~> Pa degC-3] - real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa ppt-1] - real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 ppt-1] + real :: b4s ! Partly rescaled version of b4 [Pa S-1 ~> Pa PSU-1] + real :: b5s ! Partly rescaled version of b5 [Pa C-1 S-1 ~> Pa degC-1 PSU-1] real :: c1s ! Partly rescaled version of c1 [m2 s-2 C-1 ~> m2 s-2 degC-1] real :: c2s ! Partly rescaled version of c2 [m2 s-2 C-2 ~> m2 s-2 degC-2] real :: c3s ! Partly rescaled version of c3 [m2 s-2 C-3 ~> m2 s-2 degC-3] - real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 ppt-1] - real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 ppt-1] + real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] + real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. - real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants. - real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] + real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, 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 + 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 al0_scale = 1.0 ; if (present(SV_scale)) al0_scale = SV_scale @@ -842,7 +899,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & 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. + ! is linear, but for T and S it may be thickness weighted. al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i+1,j) p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i+1,j) lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j) @@ -883,7 +940,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & 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. + ! is linear, but for T and S it may be thickness weighted. al0 = wt_L*al0_2d(i,j) + wt_R*al0_2d(i,j+1) p0 = wt_L*p0_2d(i,j) + wt_R*p0_2d(i,j+1) lambda = wt_L*lambda_2d(i,j) + wt_R*lambda_2d(i,j+1) diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index 2b4f99adf0..07464861cb 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -21,16 +21,16 @@ module MOM_EOS_linear ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units ! vary with the Boussinesq approximation, the Boussinesq variant is given first. -!> Compute the density of sea water (in kg/m^3), or its anomaly from a reference density, -!! using a simple linear equation of state from salinity (in psu), potential temperature (in deg C) -!! and pressure [Pa]. +!> Compute the density of sea water (in [kg m-3]), or its anomaly from a reference density, +!! using a simple linear equation of state from salinity in practical salinity units ([PSU]), +!! potential temperature in degrees Celsius ([degC]) and pressure [Pa]. interface calculate_density_linear module procedure calculate_density_scalar_linear, calculate_density_array_linear end interface calculate_density_linear -!> Compute the specific volume of sea water (in m^3/kg), or its anomaly from a reference value, -!! using a simple linear equation of state from salinity (in psu), potential temperature (in deg C) -!! and pressure [Pa]. +!> Compute the specific volume of sea water (in [m3 kg-1]), or its anomaly from a reference value, +!! using a simple linear equation of state from salinity in practical salinity units ([PSU]), +!! potential temperature in degrees Celsius ([degC]) and pressure [Pa]. interface calculate_spec_vol_linear module procedure calculate_spec_vol_scalar_linear, calculate_spec_vol_array_linear end interface calculate_spec_vol_linear @@ -75,7 +75,7 @@ subroutine calculate_density_scalar_linear(T, S, pressure, rho, & end subroutine calculate_density_scalar_linear !> This subroutine computes the density of sea water with a trivial -!! linear equation of state (in kg/m^3) from salinity (sal in psu), +!! linear equation of state (in [kg m-3]) from salinity (sal [PSU]), !! potential temperature (T [degC]), and pressure [Pa]. subroutine calculate_density_array_linear(T, S, pressure, rho, start, npts, & Rho_T0_S0, dRho_dT, dRho_dS, rho_ref) @@ -561,8 +561,8 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & 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 + 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 @@ -612,7 +612,7 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & 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. + ! is linear, but for T and S it may be thickness weighted. dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i+1,j) - p_t(i+1,j)) dRho_TS = dRho_dT*(wtT_L*T(i,j) + wtT_R*T(i+1,j)) + & @@ -657,7 +657,7 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & 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. + ! is linear, but for T and S it may be thickness weighted. dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i,j+1) - p_t(i,j+1)) dRho_TS = dRho_dT*(wtT_L*T(i,j) + wtT_R*T(i,j+1)) + & diff --git a/src/equation_of_state/MOM_TFreeze.F90 b/src/equation_of_state/MOM_TFreeze.F90 index f0b22c8f4e..16a64c89ed 100644 --- a/src/equation_of_state/MOM_TFreeze.F90 +++ b/src/equation_of_state/MOM_TFreeze.F90 @@ -28,7 +28,7 @@ module MOM_TFreeze module procedure calculate_TFreeze_Millero_scalar, calculate_TFreeze_Millero_array end interface calculate_TFreeze_Millero -!> Compute the freezing point conservative temperature [degC] from absolute salinity [g/kg] +!> Compute the freezing point conservative temperature [degC] from absolute salinity [g kg-1] !! and pressure [Pa] using the TEOS10 package. interface calculate_TFreeze_teos10 module procedure calculate_TFreeze_teos10_scalar, calculate_TFreeze_teos10_array @@ -84,13 +84,15 @@ end subroutine calculate_TFreeze_linear_array !! expression for potential temperature (not in situ temperature), using a !! value that is correct at the freezing point at 35 PSU and 5e6 Pa (500 dbar). subroutine calculate_TFreeze_Millero_scalar(S, pres, T_Fr) - real, intent(in) :: S !< Salinity in PSU. - real, intent(in) :: pres !< Pressure [Pa]. - real, intent(out) :: T_Fr !< Freezing point potential temperature [degC]. + real, intent(in) :: S !< Salinity [PSU] + real, intent(in) :: pres !< Pressure [Pa] + real, intent(out) :: T_Fr !< Freezing point potential temperature [degC] ! Local variables - real, parameter :: cS1 = -0.0575, cS3_2 = 1.710523e-3, cS2 = -2.154996e-4 - real, parameter :: dTFr_dp = -7.75e-8 + real, parameter :: cS1 = -0.0575 ! A term in the freezing point fit [degC PSU-1] + real, parameter :: cS3_2 = 1.710523e-3 ! A term in the freezing point fit [degC PSU-3/2] + real, parameter :: cS2 = -2.154996e-4 ! A term in the freezing point fit [degC PSU-2] + real, parameter :: dTFr_dp = -7.75e-8 ! Derivative of freezing point with pressure [degC Pa-1] T_Fr = S*(cS1 + (cS3_2 * sqrt(max(S,0.0)) + cS2 * S)) + dTFr_dp*pres @@ -110,8 +112,10 @@ subroutine calculate_TFreeze_Millero_array(S, pres, T_Fr, start, npts) integer, intent(in) :: npts !< The number of values to calculate. ! Local variables - real, parameter :: cS1 = -0.0575, cS3_2 = 1.710523e-3, cS2 = -2.154996e-4 - real, parameter :: dTFr_dp = -7.75e-8 + real, parameter :: cS1 = -0.0575 ! A term in the freezing point fit [degC PSU-1] + real, parameter :: cS3_2 = 1.710523e-3 ! A term in the freezing point fit [degC PSU-3/2] + real, parameter :: cS2 = -2.154996e-4 ! A term in the freezing point fit [degC PSU-2] + real, parameter :: dTFr_dp = -7.75e-8 ! Derivative of freezing point with pressure [degC Pa-1] integer :: j do j=start,start+npts-1 @@ -121,17 +125,18 @@ subroutine calculate_TFreeze_Millero_array(S, pres, T_Fr, start, npts) end subroutine calculate_TFreeze_Millero_array -!> This subroutine computes the freezing point conservative temperature -!! [degC] from absolute salinity [g/kg], and pressure [Pa] using the +!> This subroutine computes the freezing point conservative temperature [degC] +!! from absolute salinity [g kg-1], and pressure [Pa] using the !! TEOS10 package. subroutine calculate_TFreeze_teos10_scalar(S, pres, T_Fr) - real, intent(in) :: S !< Absolute salinity [g/kg]. + real, intent(in) :: S !< Absolute salinity [g kg-1]. real, intent(in) :: pres !< Pressure [Pa]. real, intent(out) :: T_Fr !< Freezing point conservative temperature [degC]. ! Local variables - real, dimension(1) :: S0, pres0 - real, dimension(1) :: tfr0 + real, dimension(1) :: S0 ! Salinity at a point [g kg-1] + real, dimension(1) :: pres0 ! Pressure at a point [Pa] + real, dimension(1) :: tfr0 ! The freezing temperature [degC] S0(1) = S pres0(1) = pres @@ -141,22 +146,23 @@ subroutine calculate_TFreeze_teos10_scalar(S, pres, T_Fr) end subroutine calculate_TFreeze_teos10_scalar -!> This subroutine computes the freezing point conservative temperature -!! [degC] from absolute salinity [g/kg], and pressure [Pa] using the +!> This subroutine computes the freezing point conservative temperature [degC] +!! from absolute salinity [g kg-1], and pressure [Pa] using the !! TEOS10 package. subroutine calculate_TFreeze_teos10_array(S, pres, T_Fr, start, npts) - real, dimension(:), intent(in) :: S !< absolute salinity [g/kg]. + real, dimension(:), intent(in) :: S !< absolute salinity [g kg-1]. real, dimension(:), intent(in) :: pres !< pressure [Pa]. real, dimension(:), intent(out) :: T_Fr !< Freezing point conservative temperature [degC]. integer, intent(in) :: start !< the starting point in the arrays. integer, intent(in) :: npts !< the number of values to calculate. ! Local variables - real, parameter :: Pa2db = 1.e-4 ! The conversion factor from Pa to dbar. - real :: zs,zp + real, parameter :: Pa2db = 1.e-4 ! The conversion factor from Pa to dbar [dbar Pa-1] + real :: zs ! Salinity at a point [g kg-1] + real :: zp ! Pressures in [dbar] integer :: j ! Assume sea-water contains no dissolved air. - real, parameter :: saturation_fraction = 0.0 + real, parameter :: saturation_fraction = 0.0 ! Air saturation fraction in seawater [nondim] do j=start,start+npts-1 !Conversions