MOM_EOS_UNESCO.F90

!> The equation of state using the Jackett and McDougall fits to the UNESCO EOS
module mom_eos_unesco
 
! This file is part of MOM6. See LICENSE.md for the license.
 
!***********************************************************************
!*  The subroutines in this file implement the equation of state for   *
!*  sea water using the fit to the UNESCO equation of state given by   *
!*  the expressions from Jackett and McDougall, 1995, J. Atmos.        *
!*  Ocean. Tech., 12, 381-389.  Coded by J. Stephens, 9/99.            *
!***********************************************************************
 
implicit none ; private
 
public calculate_compress_unesco, calculate_density_unesco, calculate_spec_vol_unesco
public calculate_density_derivs_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.
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.
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
! 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
!!@}
 
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.
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].
  real,           intent(in)  :: pressure !< pressure [Pa].
  real,           intent(out) :: rho      !< In situ density [kg m-3].
  real, optional, intent(in)  :: rho_ref  !< A reference density [kg m-3].
 
  ! Local variables
  real, dimension(1) :: t0, s0, pressure0
  real, dimension(1) :: rho0
 
  t0(1) = t
  s0(1) = s
  pressure0(1) = pressure
 
  call calculate_density_array_unesco(t0, s0, pressure0, rho0, 1, 1, rho_ref)
  rho = rho0(1)
 
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.
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].
  real, dimension(:), intent(in)  :: pressure !< pressure [Pa].
  real, dimension(:), intent(out) :: rho      !< in situ density [kg m-3].
  integer,            intent(in)  :: start    !< the starting point in the arrays.
  integer,            intent(in)  :: npts     !< the number of values to calculate.
  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 :: 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].
  real :: ks          ! The secant bulk modulus [bar].
  integer :: j
 
  do j=start,start+npts-1
    if (s(j) < -1.0e-10) then !Can we assume safely that this is a missing value?
      rho(j) = 1000.0
      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)
 
!  Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ).
 
    sig0 = r10*t_local + r20*t2 + r30*t3 + r40*t4 + r50*t5 + &
           s_local*(r01 + r11*t_local + r21*t2 + r31*t3 + r41*t4) + &
           s32*(r032 + r132*t_local + r232*t2) + r02*s2
    rho0 = r00 + sig0
 
!  Compute rho(s,theta,p), first calculating the secant bulk modulus.
 
    ks = s00 + s10*t_local + s20*t2 + s30*t3 + s40*t4 + s_local*(s01 + s11*t_local + s21*t2 + s31*t3) + &
         s32*(s032 + s132*t_local + s232*t2) + &
         p1*(sp00 + sp10*t_local + sp20*t2 + sp30*t3 + &
             s_local*(sp01 + sp11*t_local + sp21*t2) + sp032*s32) + &
         p2*(sp000 + sp010*t_local + sp020*t2 + s_local*(sp001 + sp011*t_local + sp021*t2))
 
    if (present(rho_ref)) then
      rho(j) = ((r00 - rho_ref)*ks + (sig0*ks + p1*rho_ref)) / (ks - p1)
    else
      rho(j) = rho0*ks / (ks - p1)
    endif
  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.
!! 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
                                          !! [degC].
  real,           intent(in)  :: S        !< salinity [PSU].
  real,           intent(in)  :: pressure !< pressure [Pa].
  real,           intent(out) :: specvol  !< in situ specific volume [m3 kg-1].
  real, optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
 
  ! Local variables
  real, dimension(1) :: T0, S0, pressure0, spv0
 
  t0(1) = t ; s0(1) = s ; pressure0(1) = pressure
 
  call calculate_spec_vol_array_unesco(t0, s0, pressure0, spv0, 1, 1, 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.
!! 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
                                              !! [degC].
  real, dimension(:), intent(in)  :: S        !< salinity [PSU].
  real, dimension(:), intent(in)  :: pressure !< pressure [Pa].
  real, dimension(:), intent(out) :: specvol  !< in situ specific volume [m3 kg-1].
  integer,            intent(in)  :: start    !< the starting point in the arrays.
  integer,            intent(in)  :: npts     !< the number of values to calculate.
  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 :: 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].
  integer :: j
 
  do j=start,start+npts-1
    if (s(j) < -1.0e-10) then !Can we assume safely that this is a missing value?
      specvol(j) = 0.001
      if (present(spv_ref)) specvol(j) = 0.001 - spv_ref
      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)
 
!  Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ).
 
    rho0 = r00 + r10*t_local + r20*t2 + r30*t3 + r40*t4 + r50*t5 + &
           s_local*(r01 + r11*t_local + r21*t2 + r31*t3 + r41*t4) + &
           s32*(r032 + r132*t_local + r232*t2) + r02*s2
 
!  Compute rho(s,theta,p), first calculating the secant bulk modulus.
 
    ks = s00 + s10*t_local + s20*t2 + s30*t3 + s40*t4 + s_local*(s01 + s11*t_local + s21*t2 + s31*t3) + &
         s32*(s032 + s132*t_local + s232*t2) + &
         p1*(sp00 + sp10*t_local + sp20*t2 + sp30*t3 + &
             s_local*(sp01 + sp11*t_local + sp21*t2) + sp032*s32) + &
         p2*(sp000 + sp010*t_local + sp020*t2 + s_local*(sp001 + sp011*t_local + sp021*t2))
 
    if (present(spv_ref)) then
      specvol(j) = (ks*(1.0 - (rho0*spv_ref)) - p1) / (rho0*ks)
    else
      specvol(j) = (ks - p1) / (rho0*ks)
    endif
  enddo
end subroutine calculate_spec_vol_array_unesco
 
 
!> This subroutine calculates the partial derivatives of density
!! with potential temperature and salinity.
subroutine calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts)
  real,    intent(in),  dimension(:) :: t        !< Potential temperature relative to the surface
                                                 !! [degC].
  real,    intent(in),  dimension(:) :: s        !< Salinity [PSU].
  real,    intent(in),  dimension(:) :: pressure !< Pressure [Pa].
  real,    intent(out), dimension(:) :: drho_dt  !< The partial derivative of density with potential
                                                 !! temperature [kg m-3 degC-1].
  real,    intent(out), dimension(:) :: drho_ds  !< The partial derivative of density with salinity,
                                                 !! in [kg m-3 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 :: 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 :: 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 :: drho0_dt        ! Derivative of rho0 with T [kg m-3 degC-1].
  real :: drho0_ds        ! Derivative of rho0 with S [kg m-3 PSU-1].
  real :: dks_dt          ! Derivative of ks with T [bar degC-1].
  real :: dks_ds          ! Derivative of ks with S [bar psu-1].
  real :: denom           ! 1.0 / (ks - p1) [bar-1].
  integer :: j
 
  do j=start,start+npts-1
    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
      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
 
!       compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) )
 
    rho0 = r00 + r10*t_local + r20*t2 + r30*t3 + r40*t4 + r50*t5 + &
           s_local*(r01 + r11*t_local + r21*t2 + r31*t3 + r41*t4) + &
           s32*(r032 + r132*t_local + r232*t2) + r02*s2
    drho0_dt = r10 + 2.0*r20*t_local + 3.0*r30*t2 + 4.0*r40*t3 + 5.0*r50*t4 + &
               s_local*(r11 + 2.0*r21*t_local + 3.0*r31*t2 + 4.0*r41*t3) + &
               s32*(r132 + 2.0*r232*t_local)
    drho0_ds = (r01 + r11*t_local + r21*t2 + r31*t3 + r41*t4) + &
               1.5*s12*(r032 + r132*t_local + r232*t2) + 2.0*r02*s_local
 
!       compute rho(s,theta,p)
 
    ks = s00 + s10*t_local + s20*t2 + s30*t3 + s40*t4 + s_local*(s01 + s11*t_local + s21*t2 + s31*t3) + &
         s32*(s032 + s132*t_local + s232*t2) + &
         p1*(sp00 + sp10*t_local + sp20*t2 + sp30*t3 + &
             s_local*(sp01 + sp11*t_local + sp21*t2) + sp032*s32) + &
         p2*(sp000 + sp010*t_local + sp020*t2 + s_local*(sp001 + sp011*t_local + sp021*t2))
    dks_dt = s10 + 2.0*s20*t_local + 3.0*s30*t2 + 4.0*s40*t3 + &
             s_local*(s11 + 2.0*s21*t_local + 3.0*s31*t2) + s32*(s132 + 2.0*s232*t_local) + &
             p1*(sp10 + 2.0*sp20*t_local + 3.0*sp30*t2 + s_local*(sp11 + 2.0*sp21*t_local)) + &
             p2*(sp010 + 2.0*sp020*t_local + s_local*(sp011 + 2.0*sp021*t_local))
    dks_ds = (s01 + s11*t_local + s21*t2 + s31*t3) + 1.5*s12*(s032 + s132*t_local + s232*t2) + &
             p1*(sp01 + sp11*t_local + sp21*t2 + 1.5*sp032*s12) + &
             p2*(sp001 + sp011*t_local + sp021*t2)
 
    denom = 1.0 / (ks - p1)
    drho_dt(j) = denom*(ks*drho0_dt - rho0*p1*denom*dks_dt)
    drho_ds(j) = denom*(ks*drho0_ds - rho0*p1*denom*dks_ds)
  enddo
 
end subroutine calculate_density_derivs_unesco
 
!> This subroutine computes the in situ density of sea water (rho)
!! and the compressibility (drho/dp == C_sound^-2) at the given
!! salinity, potential temperature, and pressure.
subroutine calculate_compress_unesco(T, S, pressure, rho, drho_dp, start, npts)
  real,    intent(in),  dimension(:) :: t        !< Potential temperature relative to the surface
                                                 !! [degC].
  real,    intent(in),  dimension(:) :: s        !< Salinity [PSU].
  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
                                                 !! (also the inverse of the square of sound speed)
                                                 !! [s2 m-2].
  integer, intent(in)                :: start    !< The starting point in the arrays.
  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.
  integer :: j
 
  do j=start,start+npts-1
    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
      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)
 
!  Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ).
 
    rho0 = r00 + r10*t_local + r20*t2 + r30*t3 + r40*t4 + r50*t5 + &
           s_local*(r01 + r11*t_local + r21*t2 + r31*t3 + r41*t4) + &
           s32*(r032 + r132*t_local + r232*t2) + r02*s2
 
!  Compute rho(s,theta,p), first calculating the secant bulk modulus.
    ks_0 = s00 + s10*t_local + s20*t2 + s30*t3 + s40*t4 + &
           s_local*(s01 + s11*t_local + s21*t2 + s31*t3) + s32*(s032 + s132*t_local + s232*t2)
    ks_1 = sp00 + sp10*t_local + sp20*t2 + sp30*t3 + &
           s_local*(sp01 + sp11*t_local + sp21*t2) + sp032*s32
    ks_2 = sp000 + sp010*t_local + sp020*t2 + s_local*(sp001 + sp011*t_local + sp021*t2)
 
    ks = ks_0 + p1*ks_1 + p2*ks_2
    dks_dp = ks_1 + 2.0*p1*ks_2
 
    rho(j) = rho0*ks / (ks - p1)
! The factor of 1.0e-5 is because pressure here is in bars, not Pa.
    drho_dp(j) = 1.0e-5 * (rho(j) / (ks - p1)) * (1.0 - dks_dp*p1/ks)
  enddo
end subroutine calculate_compress_unesco
 
 
end module mom_eos_unesco