Merge 07e0602 into 1444864

src/core/MOM_unit_tests.F90

@@ -11,6 +11,7 @@ module MOM_unit_tests
 use MOM_random,                     only : random_unit_tests
 use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests
 use MOM_CFC_cap,                    only : CFC_cap_unit_tests
+use MOM_EOS,                        only : EOS_unit_tests
 implicit none ; private
 
 public unit_tests
@@ -30,6 +31,8 @@ subroutine unit_tests(verbosity)
   if (is_root_pe()) then ! The following need only be tested on 1 PE
     if (string_functions_unit_tests(verbose)) call MOM_error(FATAL, &
        "MOM_unit_tests: string_functions_unit_tests FAILED")
+    if (EOS_unit_tests(verbose)) call MOM_error(FATAL, &
+       "MOM_unit_tests: EOS_unit_tests FAILED")
     if (remapping_unit_tests(verbose)) call MOM_error(FATAL, &
        "MOM_unit_tests: remapping_unit_tests FAILED")
     if (neutral_diffusion_unit_tests(verbose)) call MOM_error(FATAL, &

src/equation_of_state/MOM_EOS_TEOS10.F90

@@ -17,9 +17,8 @@ module MOM_EOS_TEOS10
 implicit none ; private
 
 public calculate_compress_teos10, calculate_density_teos10, calculate_spec_vol_teos10
-public calculate_density_derivs_teos10
-public calculate_specvol_derivs_teos10
-public calculate_density_second_derivs_teos10
+public calculate_density_derivs_teos10, calculate_specvol_derivs_teos10
+public calculate_density_second_derivs_teos10, EoS_fit_range_teos10
 public gsw_sp_from_sr, gsw_pt_from_ct
 
 !> Compute the in situ density of sea water ([kg m-3]), or its anomaly with respect to
@@ -369,4 +368,25 @@ subroutine calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts)
   enddo
 end subroutine calculate_compress_teos10
 
+
+!> Return the range of temperatures, salinities and pressures for which the TEOS-10
+!! equation of state has been fitted to observations.  Care should be taken when
+!! applying this equation of state outside of its fit range.
+subroutine EoS_fit_range_teos10(T_min, T_max, S_min, S_max, p_min, p_max)
+  real, optional, intent(out) :: T_min !< The minimum conservative temperature over which this EoS is fitted [degC]
+  real, optional, intent(out) :: T_max !< The maximum conservative temperature over which this EoS is fitted [degC]
+  real, optional, intent(out) :: S_min !< The minimum absolute salinity over which this EoS is fitted [g kg-1]
+  real, optional, intent(out) :: S_max !< The maximum absolute salinity over which this EoS is fitted [g kg-1]
+  real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
+  real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
+
+  if (present(T_min)) T_min = -6.0
+  if (present(T_max)) T_max = 40.0
+  if (present(S_min)) S_min =  0.0
+  if (present(S_max)) S_max = 42.0
+  if (present(p_min)) p_min = 0.0
+  if (present(p_max)) p_max = 1.0e8
+
+end subroutine EoS_fit_range_teos10
+
 end module MOM_EOS_TEOS10

src/equation_of_state/MOM_EOS_linear.F90

@@ -7,13 +7,11 @@ module MOM_EOS_linear
 
 implicit none ; private
 
-#include <MOM_memory.h>
-
 public calculate_compress_linear, calculate_density_linear, calculate_spec_vol_linear
 public calculate_density_derivs_linear, calculate_density_derivs_scalar_linear
 public calculate_specvol_derivs_linear
 public calculate_density_scalar_linear, calculate_density_array_linear
-public calculate_density_second_derivs_linear
+public calculate_density_second_derivs_linear, EoS_fit_range_linear
 public int_density_dz_linear, int_spec_vol_dp_linear
 
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
@@ -119,7 +117,7 @@ subroutine calculate_spec_vol_scalar_linear(T, S, pressure, specvol, &
   real, optional, intent(in)  :: spv_ref  !< A reference specific volume [m3 kg-1].
 
   if (present(spv_ref)) then
-    specvol = ((1.0 - Rho_T0_S0*spv_ref) + spv_ref*(dRho_dT*T + dRho_dS*S)) / &
+    specvol = ((1.0 - Rho_T0_S0*spv_ref) - spv_ref*(dRho_dT*T + dRho_dS*S)) / &
              ( Rho_T0_S0 + (dRho_dT*T + dRho_dS*S))
   else
     specvol = 1.0 / ( Rho_T0_S0 + (dRho_dT*T + dRho_dS*S))
@@ -148,7 +146,7 @@ subroutine calculate_spec_vol_array_linear(T, S, pressure, specvol, start, npts,
   integer :: j
 
   if (present(spv_ref)) then ; do j=start,start+npts-1
-    specvol(j) = ((1.0 - Rho_T0_S0*spv_ref) + spv_ref*(dRho_dT*T(j) + dRho_dS*S(j))) / &
+    specvol(j) = ((1.0 - Rho_T0_S0*spv_ref) - spv_ref*(dRho_dT*T(j) + dRho_dS*S(j))) / &
                  ( Rho_T0_S0 + (dRho_dT*T(j) + dRho_dS*S(j)))
   enddo ; else ; do j=start,start+npts-1
     specvol(j) = 1.0 / ( Rho_T0_S0 + (dRho_dT*T(j) + dRho_dS*S(j)))
@@ -320,6 +318,26 @@ subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts,&
   enddo
 end subroutine calculate_compress_linear
 
+!> Return the range of temperatures, salinities and pressures for which the reduced-range equation
+!! of state from Wright (1997) has been fitted to observations.  Care should be taken when applying
+!! this equation of state outside of its fit range.
+subroutine EoS_fit_range_linear(T_min, T_max, S_min, S_max, p_min, p_max)
+  real, optional, intent(out) :: T_min !< The minimum potential temperature over which this EoS is fitted [degC]
+  real, optional, intent(out) :: T_max !< The maximum potential temperature over which this EoS is fitted [degC]
+  real, optional, intent(out) :: S_min !< The minimum salinity over which this EoS is fitted [ppt]
+  real, optional, intent(out) :: S_max !< The maximum salinity over which this EoS is fitted [ppt]
+  real, optional, intent(out) :: p_min !< The minimum pressure over which this EoS is fitted [Pa]
+  real, optional, intent(out) :: p_max !< The maximum pressure over which this EoS is fitted [Pa]
+
+  if (present(T_min)) T_min = -273.0
+  if (present(T_max)) T_max = 100.0
+  if (present(S_min)) S_min = 0.0
+  if (present(S_max)) S_max = 1000.0
+  if (present(p_min)) p_min = 0.0
+  if (present(p_max)) p_max = 1.0e9
+
+end subroutine EoS_fit_range_linear
+
 !>   This subroutine calculates analytical and nearly-analytical integrals of
 !! pressure anomalies across layers, which are required for calculating the
 !! finite-volume form pressure accelerations in a Boussinesq model.

src/equation_of_state/MOM_TFreeze.F90

@@ -5,13 +5,14 @@ module MOM_TFreeze
 
 !********+*********+*********+*********+*********+*********+*********+**
 !*  The subroutines in this file determine the potential temperature   *
-!* at which sea-water freezes.                                         *
+!* or conservative temperature at which sea-water freezes.             *
 !********+*********+*********+*********+*********+*********+*********+**
 use gsw_mod_toolbox, only : gsw_ct_freezing_exact
 
 implicit none ; private
 
 public calculate_TFreeze_linear, calculate_TFreeze_Millero, calculate_TFreeze_teos10
+public calculate_TFreeze_TEOS_poly
 
 !> Compute the freezing point potential temperature [degC] from salinity [ppt] and
 !! pressure [Pa] using a simple linear expression, with coefficients passed in as arguments.
@@ -34,11 +35,17 @@ module MOM_TFreeze
   module procedure calculate_TFreeze_teos10_scalar, calculate_TFreeze_teos10_array
 end interface calculate_TFreeze_teos10
 
+!> Compute the freezing point conservative temperature [degC] from absolute salinity [g kg-1] and
+!! pressure [Pa] using a rescaled and refactored version of the expressions from the TEOS10 package.
+interface calculate_TFreeze_TEOS_poly
+  module procedure calculate_TFreeze_TEOS_poly_scalar, calculate_TFreeze_TEOS_poly_array
+end interface calculate_TFreeze_TEOS_poly
+
 contains
 
-!>  This subroutine computes the freezing point potential temperature
-!!  [degC] from salinity [ppt], and pressure [Pa] using a simple
-!!  linear expression, with coefficients passed in as arguments.
+!>  This subroutine computes the freezing point potential temperature [degC] from
+!!  salinity [ppt], and pressure [Pa] using a simple linear expression,
+!!  with coefficients passed in as arguments.
 subroutine calculate_TFreeze_linear_scalar(S, pres, T_Fr, TFr_S0_P0, &
                                            dTFr_dS, dTFr_dp)
   real,  intent(in)  :: S         !< salinity [ppt].
@@ -66,7 +73,7 @@ subroutine calculate_TFreeze_linear_array(S, pres, T_Fr, start, npts, &
   integer,             intent(in)  :: npts      !< the number of values to calculate.
   real,                intent(in)  :: TFr_S0_P0 !< The freezing point at S=0, p=0, [degC].
   real,                intent(in)  :: dTFr_dS   !< The derivative of freezing point with salinity,
-                                                !! [degC PSU-1].
+                                                !! [degC ppt-1].
   real,                intent(in)  :: dTFr_dp   !< The derivative of freezing point with pressure,
                                                 !! [degC Pa-1].
   integer :: j
@@ -94,13 +101,13 @@ subroutine calculate_TFreeze_Millero_scalar(S, pres, T_Fr)
   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
+  T_Fr = S*(cS1 + (cS3_2 * sqrt(max(S, 0.0)) + cS2 * S)) + dTFr_dp*pres
 
 end subroutine calculate_TFreeze_Millero_scalar
 
 !> This subroutine computes the freezing point potential temperature
 !! [degC] from salinity [ppt], and pressure [Pa] using the expression
-!! from Millero (1978) (and in appendix A of Gill 1982), but with the of the
+!! from Millero (1978) (and in appendix A of Gill 1982), but with the
 !! pressure dependence changed from 7.53e-8 to 7.75e-8 to make this an
 !! 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).
@@ -119,12 +126,82 @@ subroutine calculate_TFreeze_Millero_array(S, pres, T_Fr, start, npts)
   integer :: j
 
   do j=start,start+npts-1
-    T_Fr(j) = S(j)*(cS1 + (cS3_2 * sqrt(max(S(j),0.0)) + cS2 * S(j))) + &
+    T_Fr(j) = S(j)*(cS1 + (cS3_2 * sqrt(max(S(j), 0.0)) + cS2 * S(j))) + &
               dTFr_dp*pres(j)
   enddo
 
 end subroutine calculate_TFreeze_Millero_array
 
+!> This subroutine computes the freezing point conservative temperature [degC]
+!! from absolute salinity [g kg-1], and pressure [Pa] using a rescaled and
+!! refactored version of the polynomial expressions from the TEOS10 package.
+subroutine calculate_TFreeze_TEOS_poly_scalar(S, pres, T_Fr)
+  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    ! 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
+
+  call calculate_TFreeze_TEOS_poly_array(S0, pres0, tfr0, 1, 1)
+  T_Fr = tfr0(1)
+
+end subroutine calculate_TFreeze_TEOS_poly_scalar
+
+!> This subroutine computes the freezing point conservative temperature [degC]
+!! from absolute salinity [g kg-1], and pressure [Pa] using a rescaled and
+!! refactored version of the polynomial expressions from the TEOS10 package.
+subroutine calculate_TFreeze_TEOS_poly_array(S, pres, T_Fr, start, npts)
+  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 :: Sa    ! Absolute salinity [g kg-1] = [ppt]
+  real :: rS    ! Square root of salinity [ppt1/2]
+  ! The coefficients here use the notation TFab for contributions proportional to S**a/2 * P**b.
+  real, parameter :: TF00 =  0.017947064327968736  ! Freezing point coefficient [degC]
+  real, parameter :: TF20 = -6.076099099929818e-2  ! Freezing point coefficient [degC ppt-1]
+  real, parameter :: TF30 =  4.883198653547851e-3  ! Freezing point coefficient [degC ppt-3/2]
+  real, parameter :: TF40 = -1.188081601230542e-3  ! Freezing point coefficient [degC ppt-2]
+  real, parameter :: TF50 =  1.334658511480257e-4  ! Freezing point coefficient [degC ppt-5/2]
+  real, parameter :: TF60 = -8.722761043208607e-6  ! Freezing point coefficient [degC ppt-3]
+  real, parameter :: TF70 =  2.082038908808201e-7  ! Freezing point coefficient [degC ppt-7/2]
+  real, parameter :: TF01 = -7.389420998107497e-8  ! Freezing point coefficient [degC Pa-1]
+  real, parameter :: TF21 = -9.891538123307282e-11 ! Freezing point coefficient [degC ppt-1 Pa-1]
+  real, parameter :: TF31 = -8.987150128406496e-13 ! Freezing point coefficient [degC ppt-3/2 Pa-1]
+  real, parameter :: TF41 =  1.054318231187074e-12 ! Freezing point coefficient [degC ppt-2 Pa-1]
+  real, parameter :: TF51 =  3.850133554097069e-14 ! Freezing point coefficient [degC ppt-5/2 Pa-1]
+  real, parameter :: TF61 = -2.079022768390933e-14 ! Freezing point coefficient [degC ppt-3 Pa-1]
+  real, parameter :: TF71 =  1.242891021876471e-15 ! Freezing point coefficient [degC ppt-7/2 Pa-1]
+  real, parameter :: TF02 = -2.110913185058476e-16 ! Freezing point coefficient [degC Pa-2]
+  real, parameter :: TF22 =  3.831132432071728e-19 ! Freezing point coefficient [degC ppt-1 Pa-2]
+  real, parameter :: TF32 =  1.065556599652796e-19 ! Freezing point coefficient [degC ppt-3/2 Pa-2]
+  real, parameter :: TF42 = -2.078616693017569e-20 ! Freezing point coefficient [degC ppt-2 Pa-2]
+  real, parameter :: TF52 =  1.596435439942262e-21 ! Freezing point coefficient [degC ppt-5/2 Pa-2]
+  real, parameter :: TF03 =  2.295491578006229e-25 ! Freezing point coefficient [degC Pa-3]
+  real, parameter :: TF23 = -7.997496801694032e-27 ! Freezing point coefficient [degC ppt-1 Pa-3]
+  real, parameter :: TF33 =  8.756340772729538e-28 ! Freezing point coefficient [degC ppt-3/2 Pa-3]
+  real, parameter :: TF43 =  1.338002171109174e-29 ! Freezing point coefficient [degC ppt-2 Pa-3]
+  integer :: j
+
+  do j=start,start+npts-1
+    rS = sqrt(max(S(j), 0.0))
+    T_Fr(j) =       (TF00 + S(j)*(TF20 + rS*(TF30 + rS*(TF40 + rS*(TF50 + rS*(TF60 + rS*TF70)))))) &
+        + pres(j)*( (TF01 + S(j)*(TF21 + rS*(TF31 + rS*(TF41 + rS*(TF51 + rS*(TF61 + rS*TF71)))))) &
+         + pres(j)*((TF02 + S(j)*(TF22 + rS*(TF32 + rS*(TF42 + rS* TF52)))) &
+          + pres(j)*(TF03 + S(j)*(TF23 + rS*(TF33 + rS* TF43))) ) )
+  enddo
+
+end subroutine calculate_TFreeze_TEOS_poly_array
+
 !> This subroutine computes the freezing point conservative temperature [degC]
 !! from absolute salinity [g kg-1], and pressure [Pa] using the
 !! TEOS10 package.
@@ -158,19 +235,17 @@ subroutine calculate_TFreeze_teos10_array(S, pres, T_Fr, start, npts)
 
   ! Local variables
   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 ! Air saturation fraction in seawater [nondim]
 
   do j=start,start+npts-1
     !Conversions
-    zs = S(j)
     zp = pres(j)* Pa2db         !Convert pressure from Pascal to decibar
 
     if (S(j) < -1.0e-10) cycle !Can we assume safely that this is a missing value?
-    T_Fr(j) = gsw_ct_freezing_exact(zs,zp,saturation_fraction)
+    T_Fr(j) = gsw_ct_freezing_exact(S(j), zp, saturation_fraction)
   enddo
 
 end subroutine calculate_TFreeze_teos10_array