MOM_EOS_UNESCO.F90

Go to the documentation of this file.

module mom_eos_unesco
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!***********************************************************************
!*  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
public calculate_density_derivs_unesco
public calculate_density_scalar_unesco, calculate_density_array_unesco

interface calculate_density_unesco
  module procedure calculate_density_scalar_unesco, calculate_density_array_unesco
end interface calculate_density_unesco

! 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
!! units of kg/m^3) from salinity (S in psu), potential temperature
!! (T in deg C), and pressure in Pa.  It uses the expression from
!! Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740.
!! Coded by R. Hallberg, 7/00
subroutine calculate_density_scalar_unesco(T, S, pressure, rho)
real,    intent(in)  :: T        !< Potential temperature relative to the surface in C.
real,    intent(in)  :: S        !< Salinity in PSU.
real,    intent(in)  :: pressure !< Pressure in Pa.
real,    intent(out) :: rho      !< In situ density in kg m-3.

! * Arguments: T - potential temperature relative to the surface in C. *
! *  (in)      S - salinity in PSU.                                    *
! *  (in)      pressure - pressure in Pa.                              *
! *  (out)     rho - in situ density in kg m-3.                        *
! *  (in)      start - the starting point in the arrays.               *
! *  (in)      npts - the number of values to calculate.               *

! *====================================================================*
! *  This subroutine computes the in situ density of sea water (rho in *
! *  units of kg/m^3) from salinity (S in psu), potential temperature  *
! *  (T in deg C), and pressure in 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
  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 = rho0(1)

end subroutine calculate_density_scalar_unesco

!> This subroutine computes the in situ density of sea water (rho in
!! units of kg/m^3) from salinity (S in psu), potential temperature
!! (T in deg C), and pressure in Pa.
subroutine calculate_density_array_unesco(T, S, pressure, rho, start, npts)
  real,    intent(in),  dimension(:) :: T        !< Potential temperature relative to the surface
                                                 !! in C.
  real,    intent(in),  dimension(:) :: S        !< Salinity in PSU.
  real,    intent(in),  dimension(:) :: pressure !< Pressure in Pa.
  real,    intent(out), dimension(:) :: rho      !< In situ density in kg m-3.
  integer, intent(in)                :: start    !< The starting point in the arrays.
  integer, intent(in)                :: npts     !< The number of values to calculate.

! *  This subroutine computes the in situ density of sea water (rho in *
! *  units of kg/m^3) from salinity (S in psu), potential temperature  *
! *  (T in deg C), and pressure in Pa.                                 *

! * Arguments: T - potential temperature relative to the surface in C. *
! *  (in)      S - salinity in PSU.                                    *
! *  (in)      pressure - pressure in Pa.                              *
! *  (out)     rho - in situ density in kg m-3.                        *
! *  (in)      start - the starting point in the arrays.               *
! *  (in)      npts - the number of values to calculate.               *

  real :: t_local, t2, t3, t4, t5; ! Temperature to the 1st - 5th power.
  real :: s_local, s32, s2;        ! Salinity to the 1st, 3/2, & 2nd power.
  real :: p1, p2;      ! Pressure (in bars) to the 1st and 2nd power.
  real :: rho0;        ! Density at 1 bar pressure, in kg m-3.
  real :: ks;          ! The secant bulk modulus in bar.
  integer :: j

  do j=start,start+npts-1
    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));

    rho(j) = rho0*ks / (ks - p1);
  enddo
end subroutine calculate_density_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
                                                 !! in C.
  real,    intent(in),  dimension(:) :: S        !< Salinity in PSU.
  real,    intent(in),  dimension(:) :: pressure !< Pressure in Pa.
  real,    intent(out), dimension(:) :: drho_dT  !< The partial derivative of density with potential
                                                 !! temperature, in kg m-3 K-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.

! *   This subroutine calculates the partial derivatives of density    *
! * with potential temperature and salinity.                           *
! *                                                                    *
! * Arguments: T - potential temperature relative to the surface in C. *
! *  (in)      S - salinity in PSU.                                    *
! *  (in)      pressure - pressure in Pa.                              *
! *  (out)     drho_dT - the partial derivative of density with        *
! *                      potential temperature, in kg m-3 K-1.         *
! *  (out)     drho_dS - the partial derivative of density with        *
! *                      salinity, in kg m-3 psu-1.                    *
! *  (in)      start - the starting point in the arrays.               *
! *  (in)      npts - the number of values to calculate.               *
  real :: t_local, t2, t3, t4, t5; ! Temperature to the 1st - 5th power.
  real :: s12, s_local, s32, s2;   ! Salinity to the 1/2 - 2nd powers.
  real :: p1, p2;         ! Pressure (in bars) to the 1st & 2nd power.
  real :: rho0;           ! Density at 1 bar pressure, in kg m-3.
  real :: ks;             ! The secant bulk modulus, in bar.
  real :: drho0_dT;       ! Derivative of rho0 with T, in kg m-3 K-1.
  real :: drho0_dS;       ! Derivative of rho0 with S, kg m-3 psu-1.
  real :: dks_dT;         ! Derivative of ks with T, in bar K-1.
  real :: dks_dS;         ! Derivative of ks with S, in bar psu-1.
  real :: denom;          ! 1.0 / (ks - p1) in bar-1.
  integer :: j

  do j=start,start+npts-1
    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
                                                 !! in C.
  real,    intent(in),  dimension(:) :: S        !< Salinity in PSU.
  real,    intent(in),  dimension(:) :: pressure !< Pressure in Pa.
  real,    intent(out), dimension(:) :: rho      !< In situ density in 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)
                                                 !! in s2 m-2.
  integer, intent(in)                :: start    !< The starting point in the arrays.
  integer, intent(in)                :: npts     !< The number of values to calculate.

! *  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.                    *
! *                                                                    *
! * Arguments: T - potential temperature relative to the surface in C. *
! *  (in)      S - salinity in PSU.                                    *
! *  (in)      pressure - pressure in 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.                       *
! *  (in)      start - the starting point in the arrays.               *
! *  (in)      npts - the number of values to calculate.               *

  real :: t_local, t2, t3, t4, t5; ! Temperature to the 1st - 5th power.
  real :: s_local, s32, s2;        ! Salinity to the 1st, 3/2, & 2nd power.
  real :: p1, p2;      ! Pressure (in bars) to the 1st and 2nd power.
  real :: rho0;        ! Density at 1 bar pressure, in kg m-3.
  real :: ks;          ! The secant bulk modulus in 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
    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