MOM_EOS_NEMO.F90

Go to the documentation of this file.

module mom_eos_nemo
!***********************************************************************
!*                   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 formulae provided by NEMO developer Roquet     *
!*  in a private communication , Roquet et al, Ocean Modelling (2015)  *
!*  Roquet, F., Madec, G., McDougall, T. J., and Barker, P. M., 2015.  *
!*  Accurate polynomial expressions for the density and specific volume*
!*  of seawater using the TEOS-10 standard. Ocean Modelling, 90:29-43. *
!*  These algorithms are NOT from NEMO package!!                       *
!***********************************************************************

!use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt
use gsw_mod_toolbox, only : gsw_rho_first_derivatives

implicit none ; private

public calculate_compress_nemo, calculate_density_nemo
public calculate_density_derivs_nemo
public calculate_density_scalar_nemo, calculate_density_array_nemo

interface calculate_density_nemo
  module procedure calculate_density_scalar_nemo, calculate_density_array_nemo
end interface calculate_density_nemo

   real, parameter :: pa2db  = 1.e-4
   real, parameter :: rdeltas = 32.
   real, parameter :: r1_s0  = 0.875/35.16504
   real, parameter :: r1_t0  = 1./40.
   real, parameter :: r1_p0  = 1.e-4
   real, parameter :: r00 = 4.6494977072e+01
   real, parameter :: r01 = -5.2099962525
   real, parameter :: r02 = 2.2601900708e-01
   real, parameter :: r03 = 6.4326772569e-02
   real, parameter :: r04 = 1.5616995503e-02
   real, parameter :: r05 = -1.7243708991e-03
   real, parameter :: eos000 = 8.0189615746e+02
   real, parameter :: eos100 = 8.6672408165e+02
   real, parameter :: eos200 = -1.7864682637e+03
   real, parameter :: eos300 = 2.0375295546e+03
   real, parameter :: eos400 = -1.2849161071e+03
   real, parameter :: eos500 = 4.3227585684e+02
   real, parameter :: eos600 = -6.0579916612e+01
   real, parameter :: eos010 = 2.6010145068e+01
   real, parameter :: eos110 = -6.5281885265e+01
   real, parameter :: eos210 = 8.1770425108e+01
   real, parameter :: eos310 = -5.6888046321e+01
   real, parameter :: eos410 = 1.7681814114e+01
   real, parameter :: eos510 = -1.9193502195
   real, parameter :: eos020 = -3.7074170417e+01
   real, parameter :: eos120 = 6.1548258127e+01
   real, parameter :: eos220 = -6.0362551501e+01
   real, parameter :: eos320 = 2.9130021253e+01
   real, parameter :: eos420 = -5.4723692739
   real, parameter :: eos030 = 2.1661789529e+01
   real, parameter :: eos130 = -3.3449108469e+01
   real, parameter :: eos230 = 1.9717078466e+01
   real, parameter :: eos330 = -3.1742946532
   real, parameter :: eos040 = -8.3627885467
   real, parameter :: eos140 = 1.1311538584e+01
   real, parameter :: eos240 = -5.3563304045
   real, parameter :: eos050 = 5.4048723791e-01
   real, parameter :: eos150 = 4.8169980163e-01
   real, parameter :: eos060 = -1.9083568888e-01
   real, parameter :: eos001 = 1.9681925209e+01
   real, parameter :: eos101 = -4.2549998214e+01
   real, parameter :: eos201 = 5.0774768218e+01
   real, parameter :: eos301 = -3.0938076334e+01
   real, parameter :: eos401 = 6.6051753097
   real, parameter :: eos011 = -1.3336301113e+01
   real, parameter :: eos111 = -4.4870114575
   real, parameter :: eos211 = 5.0042598061
   real, parameter :: eos311 = -6.5399043664e-01
   real, parameter :: eos021 = 6.7080479603
   real, parameter :: eos121 = 3.5063081279
   real, parameter :: eos221 = -1.8795372996
   real, parameter :: eos031 = -2.4649669534
   real, parameter :: eos131 = -5.5077101279e-01
   real, parameter :: eos041 = 5.5927935970e-01
   real, parameter :: eos002 = 2.0660924175
   real, parameter :: eos102 = -4.9527603989
   real, parameter :: eos202 = 2.5019633244
   real, parameter :: eos012 = 2.0564311499
   real, parameter :: eos112 = -2.1311365518e-01
   real, parameter :: eos022 = -1.2419983026
   real, parameter :: eos003 = -2.3342758797e-02
   real, parameter :: eos103 = -1.8507636718e-02
   real, parameter :: eos013 = 3.7969820455e-01
   real, parameter :: alp000 = -6.5025362670e-01
   real, parameter :: alp100 = 1.6320471316
   real, parameter :: alp200 = -2.0442606277
   real, parameter :: alp300 = 1.4222011580
   real, parameter :: alp400 = -4.4204535284e-01
   real, parameter :: alp500 = 4.7983755487e-02
   real, parameter :: alp010 = 1.8537085209
   real, parameter :: alp110 = -3.0774129064
   real, parameter :: alp210 = 3.0181275751
   real, parameter :: alp310 = -1.4565010626
   real, parameter :: alp410 = 2.7361846370e-01
   real, parameter :: alp020 = -1.6246342147
   real, parameter :: alp120 = 2.5086831352
   real, parameter :: alp220 = -1.4787808849
   real, parameter :: alp320 = 2.3807209899e-01
   real, parameter :: alp030 = 8.3627885467e-01
   real, parameter :: alp130 = -1.1311538584
   real, parameter :: alp230 = 5.3563304045e-01
   real, parameter :: alp040 = -6.7560904739e-02
   real, parameter :: alp140 = -6.0212475204e-02
   real, parameter :: alp050 = 2.8625353333e-02
   real, parameter :: alp001 = 3.3340752782e-01
   real, parameter :: alp101 = 1.1217528644e-01
   real, parameter :: alp201 = -1.2510649515e-01
   real, parameter :: alp301 = 1.6349760916e-02
   real, parameter :: alp011 = -3.3540239802e-01
   real, parameter :: alp111 = -1.7531540640e-01
   real, parameter :: alp211 = 9.3976864981e-02
   real, parameter :: alp021 = 1.8487252150e-01
   real, parameter :: alp121 = 4.1307825959e-02
   real, parameter :: alp031 = -5.5927935970e-02
   real, parameter :: alp002 = -5.1410778748e-02
   real, parameter :: alp102 = 5.3278413794e-03
   real, parameter :: alp012 = 6.2099915132e-02
   real, parameter :: alp003 = -9.4924551138e-03
   real, parameter :: bet000 = 1.0783203594e+01
   real, parameter :: bet100 = -4.4452095908e+01
   real, parameter :: bet200 = 7.6048755820e+01
   real, parameter :: bet300 = -6.3944280668e+01
   real, parameter :: bet400 = 2.6890441098e+01
   real, parameter :: bet500 = -4.5221697773
   real, parameter :: bet010 = -8.1219372432e-01
   real, parameter :: bet110 = 2.0346663041
   real, parameter :: bet210 = -2.1232895170
   real, parameter :: bet310 = 8.7994140485e-01
   real, parameter :: bet410 = -1.1939638360e-01
   real, parameter :: bet020 = 7.6574242289e-01
   real, parameter :: bet120 = -1.5019813020
   real, parameter :: bet220 = 1.0872489522
   real, parameter :: bet320 = -2.7233429080e-01
   real, parameter :: bet030 = -4.1615152308e-01
   real, parameter :: bet130 = 4.9061350869e-01
   real, parameter :: bet230 = -1.1847737788e-01
   real, parameter :: bet040 = 1.4073062708e-01
   real, parameter :: bet140 = -1.3327978879e-01
   real, parameter :: bet050 = 5.9929880134e-03
   real, parameter :: bet001 = -5.2937873009e-01
   real, parameter :: bet101 = 1.2634116779
   real, parameter :: bet201 = -1.1547328025
   real, parameter :: bet301 = 3.2870876279e-01
   real, parameter :: bet011 = -5.5824407214e-02
   real, parameter :: bet111 = 1.2451933313e-01
   real, parameter :: bet211 = -2.4409539932e-02
   real, parameter :: bet021 = 4.3623149752e-02
   real, parameter :: bet121 = -4.6767901790e-02
   real, parameter :: bet031 = -6.8523260060e-03
   real, parameter :: bet002 = -6.1618945251e-02
   real, parameter :: bet102 = 6.2255521644e-02
   real, parameter :: bet012 = -2.6514181169e-03
   real, parameter :: bet003 = -2.3025968587e-04



contains

subroutine calculate_density_scalar_nemo(T, S, pressure, rho)
real,    intent(in)  :: T, S, pressure
real,    intent(out) :: rho
! * Arguments: T - conservative temperature in C.                      *
! *  (in)      S - absoulte salinity in g/Kg.                          *
! *  (in)      pressure - pressure in Pa.                              *
! *  (out)     rho - in situ density in kg m-3.                        *

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

  real :: al0, p0, lambda
  integer :: j
  real, dimension(1) :: T0, S0, pressure0
  real, dimension(1) :: rho0

  t0(1) = t
  s0(1) = s
  pressure0(1) = pressure

  call calculate_density_array_nemo(t0, s0, pressure0, rho0, 1, 1)
  rho = rho0(1)
end subroutine calculate_density_scalar_nemo
!> This subroutine computes the in situ density of sea water (rho in
!! units of kg/m^3) from absolute salinity (S in g/Kg),
!! conservative temperature (T in deg C), and pressure in Pa.
subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts)
  real,    intent(in),  dimension(:) :: T     !< Conservative temperature in C.
  real,    intent(in),  dimension(:) :: S     !< Absoulte salinity in g/Kg.
  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.

! * Arguments: T - conservative temperature in C.                      *
! *  (in)      S - absoulte salinity in g/Kg.                          *
! *  (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 absolute salinity (S in g/Kg),              *
! *  conservative temperature (T in deg C), and pressure in Pa.        *
! *====================================================================*
  real :: zp,zt , zh , zs , zr0, zn , zn0, zn1, zn2, zn3
  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
    zp = pressure(j)* pa2db         !Convert pressure from Pascal to decibar

    !The following algorithm was provided by Roquet in a private communication.
    !It is not necessarily the algorithm used in NEMO ocean!
    zp  = zp * r1_p0 !pressure
    zt  = zt * r1_t0 !temperature
    zs  = sqrt( abs( zs + rdeltas ) * r1_s0 )   ! square root salinity
    !
    zn3 = eos013*zt   &
       &   + eos103*zs+eos003
       !
    zn2 = (eos022*zt   &
       &   + eos112*zs+eos012)*zt   &
       &   + (eos202*zs+eos102)*zs+eos002
       !
    zn1 = (((eos041*zt   &
       &   + eos131*zs+eos031)*zt   &
       &   + (eos221*zs+eos121)*zs+eos021)*zt   &
       &   + ((eos311*zs+eos211)*zs+eos111)*zs+eos011)*zt   &
       &   + (((eos401*zs+eos301)*zs+eos201)*zs+eos101)*zs+eos001
       !
    zn0 = (((((eos060*zt   &
       &   + eos150*zs+eos050)*zt   &
       &   + (eos240*zs+eos140)*zs+eos040)*zt   &
       &   + ((eos330*zs+eos230)*zs+eos130)*zs+eos030)*zt   &
       &   + (((eos420*zs+eos320)*zs+eos220)*zs+eos120)*zs+eos020)*zt   &
       &   + ((((eos510*zs+eos410)*zs+eos310)*zs+eos210)*zs+eos110)*zs+eos010)*zt   &
       &   + (((((eos600*zs+eos500)*zs+eos400)*zs+eos300)*zs+eos200)*zs+eos100)*zs+eos000
       !
    zn  = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0
    !
    zr0 = (((((r05 * zp+r04) * zp+r03 ) * zp+r02 ) * zp+r01) * zp+r00) * zp
    !
    rho(j) =  ( zn + zr0 ) ! density

 enddo
end subroutine calculate_density_array_nemo

subroutine calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts)
  real,    intent(in),  dimension(:) :: T        !< Conservative temperature in C.
  real,    intent(in),  dimension(:) :: S        !< Absolute salinity in g/kg.
  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.
! * Arguments: T - conservative temperature in C.                      *
! *  (in)      S - absolute salinity in g/kg.                          *
! *  (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 :: zp,zt , zh , zs , zr0, zn , zn0, zn1, zn2, zn3
  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
    zp = pressure(j)* pa2db         !Convert pressure from Pascal to decibar

    !The following algorithm was provided by Roquet in a private communication.
    !It is not necessarily the algorithm used in NEMO ocean!
    zp  = zp * r1_p0  ! pressure (first converted to decibar)
    zt  = zt * r1_t0                ! temperature
    zs  = sqrt( abs( zs + rdeltas ) * r1_s0 )   ! square root salinity
    !
    ! alpha
    zn3 = alp003
    !
    zn2 = alp012*zt + alp102*zs+alp002
    !
    zn1 = ((alp031*zt   &
       &   + alp121*zs+alp021)*zt   &
       &   + (alp211*zs+alp111)*zs+alp011)*zt   &
       &   + ((alp301*zs+alp201)*zs+alp101)*zs+alp001
       !
    zn0 = ((((alp050*zt   &
       &   + alp140*zs+alp040)*zt   &
       &   + (alp230*zs+alp130)*zs+alp030)*zt   &
       &   + ((alp320*zs+alp220)*zs+alp120)*zs+alp020)*zt   &
       &   + (((alp410*zs+alp310)*zs+alp210)*zs+alp110)*zs+alp010)*zt   &
       &   + ((((alp500*zs+alp400)*zs+alp300)*zs+alp200)*zs+alp100)*zs+alp000
       !
    zn  = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0
    !
    drho_dt(j) = -zn
    !
    ! beta
    !
    zn3 = bet003
    !
    zn2 = bet012*zt + bet102*zs+bet002
    !
    zn1 = ((bet031*zt   &
       &   + bet121*zs+bet021)*zt   &
       &   + (bet211*zs+bet111)*zs+bet011)*zt   &
       &   + ((bet301*zs+bet201)*zs+bet101)*zs+bet001
       !
    zn0 = ((((bet050*zt   &
       &   + bet140*zs+bet040)*zt   &
       &   + (bet230*zs+bet130)*zs+bet030)*zt   &
       &   + ((bet320*zs+bet220)*zs+bet120)*zs+bet020)*zt   &
       &   + (((bet410*zs+bet310)*zs+bet210)*zs+bet110)*zs+bet010)*zt   &
       &   + ((((bet500*zs+bet400)*zs+bet300)*zs+bet200)*zs+bet100)*zs+bet000
       !
    zn  = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0
    !
    drho_ds(j) = zn / zs
  enddo

end subroutine calculate_density_derivs_nemo

subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts)
  real,    intent(in),  dimension(:) :: T        !< Conservative temperature in C.
  real,    intent(in),  dimension(:) :: S        !< Absolute salinity in g/kg.
  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.
! * Arguments: T - conservative temperature in C.                      *
! *  (in)      S - absolute salinity in g/kg.                          *
! *  (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 ::  zs,zt,zp
  integer :: j

  call calculate_density_array_nemo(t, s, pressure, rho, start, npts)
  !
  !NOTE: The following calculates the TEOS10 approximation to compressibility
  !      since the corresponding NEMO approximation is not available yet.
  !
  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
    zp = pressure(j)* pa2db         !Convert pressure from Pascal to decibar
    call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j))
 enddo
end subroutine calculate_compress_nemo

end module mom_eos_nemo