From 042a3c1a86d975562812f9e3e194f02b14cebe63 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 30 Jun 2016 19:03:38 -0400 Subject: [PATCH 1/9] Initial support for TEOS10 equation of state - Added TEOS10 GSW-Fortran as a pakage submodule https://github.com/TEOS-10/GSW-Fortran.git - Added interfaces to use TEOS10 equation of state via #override EQN_OF_STATE = "TEOS10" - Added links to a subset of GSW-Fortran modules to be compiled These are under src/equation_of_state/TEOS10 - Added a conversion from ptemp and salt to ctemp and ref. salinity after initialization in src/initialization/MOM_state_initialization.F90 - ocean_solo/ALE/z could run for 1 month using TEOS10. - Still needs to be done - Convert ctemp to ptemp at the surface before use in flux calculations - Convert from ctemp,refsalt to ptemp,salt for diagnostics - Find out why conversion of salt to absolute salinity produces junk - Find out how much things change when converting to absolute salinity instead of to reference salinity - Find how to pass to and use land_mask in TEOS10 functions instead of relying on salt<-1.0e33 --- .gitmodules | 3 + pkg/GSW-Fortran | 1 + src/equation_of_state/MOM_EOS.F90 | 59 ++++- src/equation_of_state/MOM_EOS_TEOS10.F90 | 233 ++++++++++++++++++ .../TEOS10/gsw_ct_from_pt.f90 | 1 + .../TEOS10/gsw_mod_kinds.f90 | 1 + .../TEOS10/gsw_mod_specvol_coefficients.f90 | 1 + .../TEOS10/gsw_mod_teos10_constants.f90 | 1 + .../TEOS10/gsw_mod_toolbox.f90 | 1 + src/equation_of_state/TEOS10/gsw_rho.f90 | 1 + .../TEOS10/gsw_rho_first_derivatives.f90 | 1 + src/equation_of_state/TEOS10/gsw_specvol.f90 | 1 + .../TEOS10/gsw_specvol_first_derivatives.f90 | 1 + .../TEOS10/gsw_sr_from_sp.f90 | 1 + .../MOM_state_initialization.F90 | 5 +- 15 files changed, 309 insertions(+), 2 deletions(-) create mode 160000 pkg/GSW-Fortran create mode 100644 src/equation_of_state/MOM_EOS_TEOS10.F90 create mode 120000 src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_mod_kinds.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_rho.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_specvol.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 diff --git a/.gitmodules b/.gitmodules index ae41d00a4d..637f1188ed 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,6 @@ [submodule "pkg/CVMix-src"] path = pkg/CVMix-src url = https://github.com/CVMix/CVMix-src.git +[submodule "pkg/GSW-Fortran"] + path = pkg/GSW-Fortran + url = https://github.com/TEOS-10/GSW-Fortran.git diff --git a/pkg/GSW-Fortran b/pkg/GSW-Fortran new file mode 160000 index 0000000000..29e64d6527 --- /dev/null +++ b/pkg/GSW-Fortran @@ -0,0 +1 @@ +Subproject commit 29e64d652786e1d076a05128c920f394202bfe10 diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index e05d868c88..45303a0eef 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -6,6 +6,7 @@ module MOM_EOS use MOM_EOS_linear use MOM_EOS_Wright use MOM_EOS_UNESCO +use MOM_EOS_TEOS10 use MOM_TFreeze, only : calculate_TFreeze_linear, calculate_TFreeze_Millero use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_file_parser, only : get_param, log_version, param_file_type @@ -26,6 +27,7 @@ module MOM_EOS public int_density_dz_generic_plm_analytic public find_depth_of_pressure_in_cell public calculate_TFreeze +public convert_temp_salt_for_TEOS10 !> Calculates density of sea water from T, S and P interface calculate_density @@ -61,10 +63,12 @@ module MOM_EOS integer, parameter :: EOS_LINEAR = 1 integer, parameter :: EOS_UNESCO = 2 integer, parameter :: EOS_WRIGHT = 3 +integer, parameter :: EOS_TEOS10 = 4 character*(10), parameter :: EOS_LINEAR_STRING = "LINEAR" character*(10), parameter :: EOS_UNESCO_STRING = "UNESCO" character*(10), parameter :: EOS_WRIGHT_STRING = "WRIGHT" +character*(10), parameter :: EOS_TEOS10_STRING = "TEOS10" character*(10), parameter :: EOS_DEFAULT = EOS_WRIGHT_STRING integer, parameter :: TFREEZE_LINEAR = 1 @@ -94,6 +98,8 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS) call calculate_density_scalar_unesco(T, S, pressure, rho) case (EOS_WRIGHT) call calculate_density_scalar_wright(T, S, pressure, rho) + case (EOS_TEOS10) + call calculate_density_scalar_teos10(T, S, pressure, rho) case default call MOM_error(FATAL, & "calculate_density_scalar: EOS is not valid.") @@ -122,6 +128,8 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS) call calculate_density_array_unesco(T, S, pressure, rho, start, npts) case (EOS_WRIGHT) call calculate_density_array_wright(T, S, pressure, rho, start, npts) + case (EOS_TEOS10) + call calculate_density_array_teos10(T, S, pressure, rho, start, npts) case default call MOM_error(FATAL, & "calculate_density_array: EOS%form_of_EOS is not valid.") @@ -198,6 +206,8 @@ subroutine calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npt call calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts) case (EOS_WRIGHT) call calculate_density_derivs_wright(T, S, pressure, drho_dT, drho_dS, start, npts) + case (EOS_TEOS10) + call calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, start, npts) case default call MOM_error(FATAL, & "calculate_density_derivs: EOS%form_of_EOS is not valid.") @@ -235,6 +245,8 @@ subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, enddo case (EOS_WRIGHT) call calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts) + case (EOS_TEOS10) + call calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts) case default call MOM_error(FATAL, & "calculate_density_derivs: EOS%form_of_EOS is not valid.") @@ -265,6 +277,8 @@ subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS) call calculate_compress_unesco(T, S, pressure, rho, drho_dp, start, npts) case (EOS_WRIGHT) call calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts) + case (EOS_TEOS10) + call calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts) case default call MOM_error(FATAL, & "calculate_compress: EOS%form_of_EOS is not valid.") @@ -295,6 +309,8 @@ subroutine calculate_2_densities( T, S, pressure1, pressure2, rho1, rho2, start, call calculate_2_densities_unesco(T, S, pressure1, pressure2, rho1, rho2, start, npts) case (EOS_WRIGHT) call calculate_2_densities_wright(T, S, pressure1, pressure2, rho1, rho2, start, npts) + case (EOS_TEOS10) + call calculate_2_densities_teos10(T, S, pressure1, pressure2, rho1, rho2, start, npts) case default call MOM_error(FATAL, & "calculate_2_densities: EOS%form_of_EOS is not valid.") @@ -451,7 +467,7 @@ subroutine EOS_init(param_file, EOS) call get_param(param_file, mod, "EQN_OF_STATE", tmpstr, & "EQN_OF_STATE determines which ocean equation of state \n"//& "should be used. Currently, the valid choices are \n"//& - '"LINEAR", "UNESCO", and "WRIGHT". \n'//& + '"LINEAR", "UNESCO", "WRIGHT" and "TEOS10". \n'//& "This is only used if USE_EOS is true.", default=EOS_DEFAULT) select case (uppercase(tmpstr)) case (EOS_LINEAR_STRING) @@ -460,6 +476,8 @@ subroutine EOS_init(param_file, EOS) EOS%form_of_EOS = EOS_UNESCO case (EOS_WRIGHT_STRING) EOS%form_of_EOS = EOS_WRIGHT + case (EOS_TEOS10_STRING) + EOS%form_of_EOS = EOS_TEOS10 case default call MOM_error(FATAL, "interpret_eos_selection: EQN_OF_STATE "//& trim(tmpstr) // "in input file is invalid.") @@ -2068,6 +2086,45 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & end subroutine int_spec_vol_dp_generic +!> Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 +subroutine convert_temp_salt_for_TEOS10(T, S, press, G, kd, mask_z, EOS) + use MOM_grid, only : ocean_grid_type + !> The horizontal index structure + type(ocean_grid_type), intent(in) :: G + + !> Potential temperature referenced to the surface (degC) + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: T + !> Salinity (PSU) + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: S + !> Pressure at the top of the layer in Pa. + real, dimension(:), intent(in) :: press + !> Equation of state structure + type(EOS_type), pointer :: EOS + !> 3d mask + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: mask_z + integer, intent(in) :: kd + ! + integer :: i,j,k + real :: gsw_sr_from_sp, gsw_ct_from_pt, gsw_sa_from_sp + real :: p + + if (.not.associated(EOS)) call MOM_error(FATAL, & + "convert_temp_salt_to_TEOS10 called with an unassociated EOS_type EOS.") + + if (EOS%form_of_EOS .ne. EOS_TEOS10) return + + do k=1,kd ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (mask_z(i,j,k) .ge. 1.0) then + S(i,j,k) = gsw_sr_from_sp(S(i,j,k)) +! p=press(k)/10000. !convert pascal to dbar +! S(i,j,k) = gsw_sa_from_sp(S(i,j,k),p,G%geoLonT(i,j),G%geoLatT(i,j)) + T(i,j,k) = gsw_ct_from_pt(S(i,j,k),T(i,j,k)) + endif + enddo; enddo; enddo +end subroutine convert_temp_salt_for_TEOS10 + + + end module MOM_EOS !> \namespace mom_eos diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 new file mode 100644 index 0000000000..1b4f0487aa --- /dev/null +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -0,0 +1,233 @@ +module MOM_EOS_TEOS10 +!*********************************************************************** +!* 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 given by Wright, 1997, J. Atmos. * +!* Ocean. Tech., 14, 735-740. Coded by R. Hallberg, 7/00. * +!*********************************************************************** + +use MOM_hor_index, only : hor_index_type +use gsw_mod_toolbox, only : gsw_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives + +implicit none ; private + +#include + +public calculate_compress_teos10, calculate_density_teos10 +public calculate_density_derivs_teos10, calculate_2_densities_teos10 +public calculate_specvol_derivs_teos10 +public calculate_density_scalar_teos10, calculate_density_array_teos10 + +interface calculate_density_teos10 + module procedure calculate_density_scalar_teos10, calculate_density_array_teos10 +end interface calculate_density_teos10 + +!real :: a0, a1, a2, b0, b1, b2, b3, b4, b5, c0, c1, c2, c3, c4, c5 +! One of the two following blocks of values should be commented out. +! Following are the values for the full range formula. +! +!real, parameter :: a0 = 7.133718e-4, a1 = 2.724670e-7, a2 = -1.646582e-7 +!real, parameter :: b0 = 5.613770e8, b1 = 3.600337e6, b2 = -3.727194e4 +!real, parameter :: b3 = 1.660557e2, b4 = 6.844158e5, b5 = -8.389457e3 +!real, parameter :: c0 = 1.609893e5, c1 = 8.427815e2, c2 = -6.931554 +!real, parameter :: c3 = 3.869318e-2, c4 = -1.664201e2, c5 = -2.765195 + + +! Following are the values for the reduced range formula. +real, parameter :: a0 = 7.057924e-4, a1 = 3.480336e-7, a2 = -1.112733e-7 +real, parameter :: b0 = 5.790749e8, b1 = 3.516535e6, b2 = -4.002714e4 +real, parameter :: b3 = 2.084372e2, b4 = 5.944068e5, b5 = -9.643486e3 +real, parameter :: c0 = 1.704853e5, c1 = 7.904722e2, c2 = -7.984422 +real, parameter :: c3 = 5.140652e-2, c4 = -2.302158e2, c5 = -3.079464 + +contains + +subroutine calculate_density_scalar_teos10(T, S, pressure, rho) +real, intent(in) :: T, S, pressure +real, intent(out) :: rho +! * 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 :: 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_teos10(T0, S0, pressure0, rho0, 1, 1) + rho = rho0(1) + +end subroutine calculate_density_scalar_teos10 + +subroutine calculate_density_array_teos10(T, S, pressure, rho, start, npts) + real, intent(in), dimension(:) :: T, S, pressure + real, intent(out), dimension(:) :: rho + integer, intent(in) :: start, npts +! * 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. * +! * (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. * +! * It uses the functions from TEOS10 website * +! *====================================================================* + real :: gsw_rho !external TEOS10 function + real :: p_dbar + integer :: j + + do j=start,start+npts-1 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar + rho(j) = gsw_rho(S(j),T(j),p_dbar) + enddo +end subroutine calculate_density_array_teos10 + +subroutine calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, start, npts) + real, intent(in), dimension(:) :: T, S, pressure + real, intent(out), dimension(:) :: drho_dT, drho_dS + integer, intent(in) :: start, npts +! * 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 :: p_dbar + integer :: j + + do j=start,start+npts-1 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar + call gsw_rho_first_derivatives(S(j),T(j), p_dbar, drho_dsa=drho_dS(j), drho_dct=drho_dT(j)) + enddo + +end subroutine calculate_density_derivs_teos10 + +subroutine calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts) + real, intent(in), dimension(:) :: T, S, pressure + real, intent(out), dimension(:) :: dSV_dT, dSV_dS + integer, intent(in) :: start, npts +! * Arguments: T - conservative temperature in C. * +! * (in) S - absolute salinity in g/kg. * +! * (in) pressure - pressure in Pa. * +! * (out) dSV_dT - the partial derivative of specific volume with * +! * potential temperature, in m3 kg-1 K-1. * +! * (out) dSV_dS - the partial derivative of specific volume with * +! * salinity, in m3 kg-1 / (g/kg). * +! * (in) start - the starting point in the arrays. * +! * (in) npts - the number of values to calculate. * + real :: p_dbar + integer :: j + + do j=start,start+npts-1 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar + call gsw_specvol_first_derivatives(S(j),T(j), p_dbar, v_sa=dSV_dS(j), v_ct=dSV_dT(j)) + enddo + +end subroutine calculate_specvol_derivs_teos10 + +subroutine calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts) + real, intent(in), dimension(:) :: T, S, pressure + real, intent(out), dimension(:) :: rho, drho_dp + integer, intent(in) :: start, npts +! * 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. * +! *====================================================================* +! * This subroutine computes the in situ density of sea water (rho in * +! * units of kg/m^3) and the compressibility (drho/dp = C_sound^-2) * +! * (drho_dp in units of s2 m-2) from salinity (sal in psu), potential* +! * temperature (T in deg C), and pressure in Pa. It uses the * +! * subroutines from TEOS10 website * +! *====================================================================* + real :: gsw_rho !external TEOS10 function + real :: p_dbar + integer :: j + + do j=start,start+npts-1 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar + rho(j) = gsw_rho(S(j),T(j),p_dbar) + call gsw_rho_first_derivatives(S(j),T(j), p_dbar, drho_dp=drho_dp(j)) + enddo +end subroutine calculate_compress_teos10 + +subroutine calculate_2_densities_teos10( T, S, pressure1, pressure2, rho1, rho2, start, npts) + real, intent(in), dimension(:) :: T, S + real, intent(in) :: pressure1, pressure2 + real, intent(out), dimension(:) :: rho1, rho2 + integer, intent(in) :: start, npts +! * Arguments: T - conservative temperature in C. * +! * (in) S - absolute salinity in g/kg. * +! * (in) pressure1 - the first pressure in Pa. * +! * (in) pressure2 - the second pressure in Pa. * +! * (out) rho1 - density at pressure1 in kg m-3. * +! * (out) rho2 - density at pressure2 in kg m-3. * +! * (in) start - the starting point in the arrays. * +! * (in) npts - the number of values to calculate. * + + real :: gsw_rho !external TEOS10 function + real :: p1_dbar,p2_dbar + integer :: j + + p1_dbar = pressure1*1.0e-4 !convert pressure to dbar + p2_dbar = pressure2*1.0e-4 !convert pressure to dbar + + do j=start,start+npts-1 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + rho1(j) = gsw_rho(S(j),T(j),p1_dbar) + rho2(j) = gsw_rho(S(j),T(j),p2_dbar) + enddo +end subroutine calculate_2_densities_teos10 + + +end module MOM_EOS_TEOS10 diff --git a/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 b/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 new file mode 120000 index 0000000000..d67d2df3e2 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_ct_from_pt.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 b/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 new file mode 120000 index 0000000000..fa0926e540 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/modules/gsw_mod_kinds.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 b/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 new file mode 120000 index 0000000000..934f689c20 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/modules/gsw_mod_specvol_coefficients.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 b/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 new file mode 120000 index 0000000000..17dec5add5 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/modules/gsw_mod_teos10_constants.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 b/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 new file mode 120000 index 0000000000..f2f4761ec4 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/modules/gsw_mod_toolbox.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_rho.f90 b/src/equation_of_state/TEOS10/gsw_rho.f90 new file mode 120000 index 0000000000..22eea6219a --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_rho.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_rho.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 new file mode 120000 index 0000000000..3a8ba38824 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_rho_first_derivatives.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_specvol.f90 b/src/equation_of_state/TEOS10/gsw_specvol.f90 new file mode 120000 index 0000000000..7a41a5cea0 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_specvol.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_specvol.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 new file mode 120000 index 0000000000..ee6ee1f906 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_specvol_first_derivatives.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 b/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 new file mode 120000 index 0000000000..eda229ff66 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_sr_from_sp.f90 \ No newline at end of file diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 08825eaff5..218e1d1303 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -36,7 +36,7 @@ module MOM_state_initialization use MOM_verticalGrid, only : setVerticalGridAxes, verticalGrid_type use MOM_ALE, only : pressure_gradient_plm use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type -use MOM_EOS, only : int_specific_vol_dp +use MOM_EOS, only : int_specific_vol_dp, convert_temp_salt_for_TEOS10 use user_initialization, only : user_initialize_thickness, user_initialize_velocity use user_initialization, only : user_init_temperature_salinity use user_initialization, only : user_set_Open_Bdry_Conds, user_initialize_sponges @@ -2249,6 +2249,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) press(:)=tv%p_ref + !Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 + call convert_temp_salt_for_TEOS10(temp_z,salt_z, press, G, kd, mask_z, eos) + do k=1, kd do j=js,je call calculate_density(temp_z(:,j,k),salt_z(:,j,k), press, rho_z(:,j,k), is, ie, eos) From 3dd705d68d7f0da130fe18bc5686fef9a0ca271a Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Tue, 9 Aug 2016 10:26:46 -0400 Subject: [PATCH 2/9] Adding NEMO as an EOS option - The approximate polynomials for density, alpha and beta were provided by Fabien Roquet in the form of a code snippet. He gave the reference 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. - The corresponding approximation for compressibility is lacking and it is calculated via a TEOS10 function call. - The corresponding approximation for specvol deriverties is lacking and these are calculated explicitly using the density derivatives. --- src/equation_of_state/MOM_EOS.F90 | 43 ++- src/equation_of_state/MOM_EOS_NEMO.F90 | 424 +++++++++++++++++++++++ src/equation_of_state/MOM_EOS_TEOS10.F90 | 26 +- 3 files changed, 467 insertions(+), 26 deletions(-) create mode 100644 src/equation_of_state/MOM_EOS_NEMO.F90 diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 45303a0eef..1fd055aa4c 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -7,6 +7,7 @@ module MOM_EOS use MOM_EOS_Wright use MOM_EOS_UNESCO use MOM_EOS_TEOS10 +use MOM_EOS_NEMO use MOM_TFreeze, only : calculate_TFreeze_linear, calculate_TFreeze_Millero use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_file_parser, only : get_param, log_version, param_file_type @@ -64,11 +65,13 @@ module MOM_EOS integer, parameter :: EOS_UNESCO = 2 integer, parameter :: EOS_WRIGHT = 3 integer, parameter :: EOS_TEOS10 = 4 +integer, parameter :: EOS_NEMO = 5 character*(10), parameter :: EOS_LINEAR_STRING = "LINEAR" character*(10), parameter :: EOS_UNESCO_STRING = "UNESCO" character*(10), parameter :: EOS_WRIGHT_STRING = "WRIGHT" character*(10), parameter :: EOS_TEOS10_STRING = "TEOS10" +character*(10), parameter :: EOS_NEMO_STRING = "NEMO" character*(10), parameter :: EOS_DEFAULT = EOS_WRIGHT_STRING integer, parameter :: TFREEZE_LINEAR = 1 @@ -100,6 +103,8 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS) call calculate_density_scalar_wright(T, S, pressure, rho) case (EOS_TEOS10) call calculate_density_scalar_teos10(T, S, pressure, rho) + case (EOS_NEMO) + call calculate_density_scalar_nemo(T, S, pressure, rho) case default call MOM_error(FATAL, & "calculate_density_scalar: EOS is not valid.") @@ -130,10 +135,13 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS) call calculate_density_array_wright(T, S, pressure, rho, start, npts) case (EOS_TEOS10) call calculate_density_array_teos10(T, S, pressure, rho, start, npts) + case (EOS_nemo) + call calculate_density_array_nemo (T, S, pressure, rho, start, npts) case default call MOM_error(FATAL, & "calculate_density_array: EOS%form_of_EOS is not valid.") end select + end subroutine calculate_density_array !> Calls the appropriate subroutine to calculate the freezing point for scalar inputs. @@ -194,7 +202,23 @@ subroutine calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npt integer, intent(in) :: start !< Starting index within the array integer, intent(in) :: npts !< The number of values to calculate type(EOS_type), pointer :: EOS !< Equation of state structure - + !! + !!Testing section + !!NEMO Check value: rho = 1027.45140 kg/m^3 for p=1000 dbar, ct=10 Celcius, sa=30 g/kg + !real :: T0(1),S0(1),pressure0(1) + !T0(1) = 10. + !S0(1) = 30. + !pressure0(1) = 1000. * 1e4 !pa + !call calculate_density_derivs_nemo(T0, S0, pressure0, drho_dT, drho_dS, 1,1) + !print*, 'NEMO drho: ', drho_dT(1), drho_dS(1) + !call calculate_density_derivs_teos10(T0, S0, pressure0, drho_dT, drho_dS, 1,1) + !print*, 'TEOS10 drho: ', drho_dT(1), drho_dS(1) + !call calculate_density_derivs_wright(T0, S0, pressure0, drho_dT, drho_dS, 1,1) + !print*, 'WRIGHT drho: ', drho_dT(1), drho_dS(1) + !call calculate_density_derivs_unesco(T0, S0, pressure0, drho_dT, drho_dS, 1,1) + !print*, 'UNESCO drho: ', drho_dT(1), drho_dS(1) + !stop + !! if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") @@ -208,6 +232,8 @@ subroutine calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npt call calculate_density_derivs_wright(T, S, pressure, drho_dT, drho_dS, start, npts) case (EOS_TEOS10) call calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, start, npts) + case (EOS_NEMO) + call calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) case default call MOM_error(FATAL, & "calculate_density_derivs: EOS%form_of_EOS is not valid.") @@ -247,6 +273,13 @@ subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, call calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts) case (EOS_TEOS10) call calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pressure, rho, start, npts) + call calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) + do j=start,start+npts-1 + dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) + dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) + enddo case default call MOM_error(FATAL, & "calculate_density_derivs: EOS%form_of_EOS is not valid.") @@ -279,6 +312,8 @@ subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS) call calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts) case (EOS_TEOS10) call calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts) + case (EOS_NEMO) + call calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) case default call MOM_error(FATAL, & "calculate_compress: EOS%form_of_EOS is not valid.") @@ -311,6 +346,8 @@ subroutine calculate_2_densities( T, S, pressure1, pressure2, rho1, rho2, start, call calculate_2_densities_wright(T, S, pressure1, pressure2, rho1, rho2, start, npts) case (EOS_TEOS10) call calculate_2_densities_teos10(T, S, pressure1, pressure2, rho1, rho2, start, npts) + case (EOS_NEMO) + call calculate_2_densities_nemo(T, S, pressure1, pressure2, rho1, rho2, start, npts) case default call MOM_error(FATAL, & "calculate_2_densities: EOS%form_of_EOS is not valid.") @@ -467,7 +504,7 @@ subroutine EOS_init(param_file, EOS) call get_param(param_file, mod, "EQN_OF_STATE", tmpstr, & "EQN_OF_STATE determines which ocean equation of state \n"//& "should be used. Currently, the valid choices are \n"//& - '"LINEAR", "UNESCO", "WRIGHT" and "TEOS10". \n'//& + '"LINEAR", "UNESCO", "WRIGHT", "NEMO" and "TEOS10". \n'//& "This is only used if USE_EOS is true.", default=EOS_DEFAULT) select case (uppercase(tmpstr)) case (EOS_LINEAR_STRING) @@ -478,6 +515,8 @@ subroutine EOS_init(param_file, EOS) EOS%form_of_EOS = EOS_WRIGHT case (EOS_TEOS10_STRING) EOS%form_of_EOS = EOS_TEOS10 + case (EOS_NEMO_STRING) + EOS%form_of_EOS = EOS_NEMO case default call MOM_error(FATAL, "interpret_eos_selection: EQN_OF_STATE "//& trim(tmpstr) // "in input file is invalid.") diff --git a/src/equation_of_state/MOM_EOS_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 new file mode 100644 index 0000000000..ba994c4369 --- /dev/null +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -0,0 +1,424 @@ +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_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives + + +implicit none ; private + +public calculate_compress_nemo, calculate_density_nemo +public calculate_density_derivs_nemo, calculate_2_densities_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 + +subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts) + real, intent(in), dimension(:) :: T, S, pressure + real, intent(out), dimension(:) :: rho + integer, intent(in) :: start, npts +! * 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 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + !The following algorithm was provided by Roquet in a private communication. + !It is not necessarily the algorithm used in NEMO ocean! + zp = pressure(j)* r1_P0 *Pa2db !pressure (first converted to decibar) + zt = T(j) * r1_T0 ! temperature + zs = SQRT( ABS( S(j) + 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, S, pressure + real, intent(out), dimension(:) :: drho_dT, drho_dS + integer, intent(in) :: start, npts +! * 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 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + !The following algorithm was provided by Roquet in a private communication. + !It is not necessarily the algorithm used in NEMO ocean! + zp = pressure(j)*Pa2db * r1_P0 ! pressure (first converted to decibar) + zt = T(j) * r1_T0 ! temperature + zs = SQRT( ABS( S(j) + 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, S, pressure + real, intent(out), dimension(:) :: rho, drho_dp + integer, intent(in) :: start, npts +! * 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 :: p_dbar + integer :: j + + do j=start,start+npts-1 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar + rho(j) = gsw_rho(S(j),T(j),p_dbar) + call gsw_rho_first_derivatives(S(j),T(j), p_dbar, drho_dp=drho_dp(j)) + enddo +end subroutine calculate_compress_nemo + +subroutine calculate_2_densities_nemo( T, S, pressure1, pressure2, rho1, rho2, start, npts) + real, intent(in), dimension(:) :: T, S + real, intent(in) :: pressure1, pressure2 + real, intent(out), dimension(:) :: rho1, rho2 + integer, intent(in) :: start, npts +! * Arguments: T - conservative temperature in C. * +! * (in) S - absolute salinity in g/kg. * +! * (in) pressure1 - the first pressure in Pa. * +! * (in) pressure2 - the second pressure in Pa. * +! * (out) rho1 - density at pressure1 in kg m-3. * +! * (out) rho2 - density at pressure2 in kg m-3. * +! * (in) start - the starting point in the arrays. * +! * (in) npts - the number of values to calculate. * + + real :: zp1, zp2, zt , zh , zs , zr0, zn , zn0, zn1, zn2, zn3 + integer :: j + + do j=start,start+npts-1 + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + !The following algorithm was provided by Roquet in a private communication. + !It is not necessarily the algorithm used in NEMO ocean! + zp1 = pressure1*Pa2db * r1_P0 ! pressure (first converted to decibar) + zp2 = pressure2*Pa2db * r1_P0 ! pressure (first converted to decibar) + zt = T(j) * r1_T0 ! temperature + zs = SQRT( ABS( S(j) + 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 * zp1 + zn2 ) * zp1 + zn1 ) * zp1 + zn0 + ! + zr0 = (((((R05 * zp1+R04) * zp1+R03 ) * zp1+R02 ) * zp1+R01) * zp1+R00) * zp1 + ! + rho1(j) = ( zn + zr0 ) ! density + ! + zn = ( ( zn3 * zp2 + zn2 ) * zp2 + zn1 ) * zp2 + zn0 + ! + zr0 = (((((R05 * zp2+R04) * zp2+R03 ) * zp2+R02 ) * zp2+R01) * zp2+R00) * zp2 + ! + rho2(j) = ( zn + zr0 ) ! density + enddo +end subroutine calculate_2_densities_nemo + + +end module MOM_EOS_NEMO diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 index 1b4f0487aa..d08b9ddf42 100644 --- a/src/equation_of_state/MOM_EOS_TEOS10.F90 +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -21,17 +21,13 @@ module MOM_EOS_TEOS10 !*********************************************************************** !* The subroutines in this file implement the equation of state for * -!* sea water using the formulae given by Wright, 1997, J. Atmos. * -!* Ocean. Tech., 14, 735-740. Coded by R. Hallberg, 7/00. * +!* sea water using the TEOS10 functions * !*********************************************************************** -use MOM_hor_index, only : hor_index_type use gsw_mod_toolbox, only : gsw_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives implicit none ; private -#include - public calculate_compress_teos10, calculate_density_teos10 public calculate_density_derivs_teos10, calculate_2_densities_teos10 public calculate_specvol_derivs_teos10 @@ -41,23 +37,6 @@ module MOM_EOS_TEOS10 module procedure calculate_density_scalar_teos10, calculate_density_array_teos10 end interface calculate_density_teos10 -!real :: a0, a1, a2, b0, b1, b2, b3, b4, b5, c0, c1, c2, c3, c4, c5 -! One of the two following blocks of values should be commented out. -! Following are the values for the full range formula. -! -!real, parameter :: a0 = 7.133718e-4, a1 = 2.724670e-7, a2 = -1.646582e-7 -!real, parameter :: b0 = 5.613770e8, b1 = 3.600337e6, b2 = -3.727194e4 -!real, parameter :: b3 = 1.660557e2, b4 = 6.844158e5, b5 = -8.389457e3 -!real, parameter :: c0 = 1.609893e5, c1 = 8.427815e2, c2 = -6.931554 -!real, parameter :: c3 = 3.869318e-2, c4 = -1.664201e2, c5 = -2.765195 - - -! Following are the values for the reduced range formula. -real, parameter :: a0 = 7.057924e-4, a1 = 3.480336e-7, a2 = -1.112733e-7 -real, parameter :: b0 = 5.790749e8, b1 = 3.516535e6, b2 = -4.002714e4 -real, parameter :: b3 = 2.084372e2, b4 = 5.944068e5, b5 = -9.643486e3 -real, parameter :: c0 = 1.704853e5, c1 = 7.904722e2, c2 = -7.984422 -real, parameter :: c3 = 5.140652e-2, c4 = -2.302158e2, c5 = -3.079464 contains @@ -75,8 +54,7 @@ subroutine calculate_density_scalar_teos10(T, S, pressure, rho) ! * 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 * +! * TEOS10 website. * ! *====================================================================* real :: al0, p0, lambda From b4512c0b2126a9f966a87611b85c08fe85a0cfb8 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 18 Aug 2016 13:44:49 -0400 Subject: [PATCH 3/9] Assuming model T&S are ptemp,psalt - This commit assumes that the current MOM6 prognostics T&S variables are the potential temp and practical salinity. It converts them to conservative temp and reference (absolute) salinity at each call to NEMO or TEOS10 subroutines. --- src/equation_of_state/MOM_EOS.F90 | 2 +- src/equation_of_state/MOM_EOS_NEMO.F90 | 56 +++++++++++++----- src/equation_of_state/MOM_EOS_TEOS10.F90 | 57 ++++++++++++------- .../MOM_state_initialization.F90 | 2 +- 4 files changed, 79 insertions(+), 38 deletions(-) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 46254118de..feb95dd488 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -2173,7 +2173,7 @@ subroutine convert_temp_salt_for_TEOS10(T, S, press, G, kd, mask_z, EOS) if (.not.associated(EOS)) call MOM_error(FATAL, & "convert_temp_salt_to_TEOS10 called with an unassociated EOS_type EOS.") - if (EOS%form_of_EOS .ne. EOS_TEOS10) return + if ((EOS%form_of_EOS .ne. EOS_TEOS10) .and. (EOS%form_of_EOS .ne. EOS_NEMO)) return do k=1,kd ; do j=G%jsc,G%jec ; do i=G%isc,G%iec if (mask_z(i,j,k) .ge. 1.0) then diff --git a/src/equation_of_state/MOM_EOS_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 index ba994c4369..abeb5b86a7 100644 --- a/src/equation_of_state/MOM_EOS_NEMO.F90 +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -29,7 +29,8 @@ module MOM_EOS_NEMO !* These algorithms are NOT from NEMO package!! * !*********************************************************************** -use gsw_mod_toolbox, only : gsw_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives +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 @@ -227,12 +228,17 @@ subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts) integer :: j do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? !The following algorithm was provided by Roquet in a private communication. !It is not necessarily the algorithm used in NEMO ocean! - zp = pressure(j)* r1_P0 *Pa2db !pressure (first converted to decibar) - zt = T(j) * r1_T0 ! temperature - zs = SQRT( ABS( S(j) + rdeltaS ) * r1_S0 ) ! square root salinity + 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 @@ -281,12 +287,18 @@ subroutine calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start integer :: j do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? + !The following algorithm was provided by Roquet in a private communication. !It is not necessarily the algorithm used in NEMO ocean! - zp = pressure(j)*Pa2db * r1_P0 ! pressure (first converted to decibar) - zt = T(j) * r1_T0 ! temperature - zs = SQRT( ABS( S(j) + rdeltaS ) * r1_S0 ) ! square root salinity + 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 @@ -348,14 +360,21 @@ subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) ! * (in) start - the starting point in the arrays. * ! * (in) npts - the number of values to calculate. * ! *====================================================================* - real :: p_dbar + 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 = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? - p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar - rho(j) = gsw_rho(S(j),T(j),p_dbar) - call gsw_rho_first_derivatives(S(j),T(j), p_dbar, drho_dp=drho_dp(j)) + call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j)) enddo end subroutine calculate_compress_nemo @@ -376,14 +395,21 @@ subroutine calculate_2_densities_nemo( T, S, pressure1, pressure2, rho1, rho2, s real :: zp1, zp2, zt , zh , zs , zr0, zn , zn0, zn1, zn2, zn3 integer :: j + zp1 = pressure1 * Pa2db !Convert pressure from Pascal to decibar + zp2 = pressure2 * Pa2db !Convert pressure from Pascal to decibar + do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? !The following algorithm was provided by Roquet in a private communication. !It is not necessarily the algorithm used in NEMO ocean! - zp1 = pressure1*Pa2db * r1_P0 ! pressure (first converted to decibar) - zp2 = pressure2*Pa2db * r1_P0 ! pressure (first converted to decibar) - zt = T(j) * r1_T0 ! temperature - zs = SQRT( ABS( S(j) + rdeltaS ) * r1_S0 ) ! square root salinity + zp1 = pressure1 * r1_P0 ! pressure (first converted to decibar) + zp2 = pressure2 * r1_P0 ! pressure (first converted to decibar) + zt = zt * r1_T0 ! temperature + zs = SQRT( ABS( zs + rdeltaS ) * r1_S0 ) ! square root salinity ! zn3 = EOS013*zt & & + EOS103*zs+EOS003 diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 index d08b9ddf42..5e67780790 100644 --- a/src/equation_of_state/MOM_EOS_TEOS10.F90 +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -24,6 +24,7 @@ module MOM_EOS_TEOS10 !* sea water using the TEOS10 functions * !*********************************************************************** +use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt use gsw_mod_toolbox, only : gsw_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives implicit none ; private @@ -37,6 +38,7 @@ module MOM_EOS_TEOS10 module procedure calculate_density_scalar_teos10, calculate_density_array_teos10 end interface calculate_density_teos10 + real, parameter :: Pa2db = 1.e-4 contains @@ -88,14 +90,17 @@ subroutine calculate_density_array_teos10(T, S, pressure, rho, start, npts) ! * conservative temperature (T in deg C), and pressure in Pa. * ! * It uses the functions from TEOS10 website * ! *====================================================================* - real :: gsw_rho !external TEOS10 function - real :: p_dbar + real :: zs,zt,zp integer :: j do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? - p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar - rho(j) = gsw_rho(S(j),T(j),p_dbar) + rho(j) = gsw_rho(zs,zt,zp) enddo end subroutine calculate_density_array_teos10 @@ -112,13 +117,16 @@ subroutine calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, sta ! * salinity, in kg m-3 psu-1. * ! * (in) start - the starting point in the arrays. * ! * (in) npts - the number of values to calculate. * - real :: p_dbar + real :: zs,zt,zp integer :: j do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? - p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar - call gsw_rho_first_derivatives(S(j),T(j), p_dbar, drho_dsa=drho_dS(j), drho_dct=drho_dT(j)) + call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS(j), drho_dct=drho_dT(j)) enddo end subroutine calculate_density_derivs_teos10 @@ -136,13 +144,16 @@ subroutine calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start ! * salinity, in m3 kg-1 / (g/kg). * ! * (in) start - the starting point in the arrays. * ! * (in) npts - the number of values to calculate. * - real :: p_dbar + real :: zs, zt, zp integer :: j do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? - p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar - call gsw_specvol_first_derivatives(S(j),T(j), p_dbar, v_sa=dSV_dS(j), v_ct=dSV_dT(j)) + call gsw_specvol_first_derivatives(zs,zt,zp, v_sa=dSV_dS(j), v_ct=dSV_dT(j)) enddo end subroutine calculate_specvol_derivs_teos10 @@ -167,15 +178,17 @@ subroutine calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts) ! * temperature (T in deg C), and pressure in Pa. It uses the * ! * subroutines from TEOS10 website * ! *====================================================================* - real :: gsw_rho !external TEOS10 function - real :: p_dbar + real :: zs,zt,zp integer :: j do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? - p_dbar = pressure(j)*1.0e-4 !convert pressure to dbar - rho(j) = gsw_rho(S(j),T(j),p_dbar) - call gsw_rho_first_derivatives(S(j),T(j), p_dbar, drho_dp=drho_dp(j)) + rho(j) = gsw_rho(zs,zt,zp) + call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j)) enddo end subroutine calculate_compress_teos10 @@ -193,17 +206,19 @@ subroutine calculate_2_densities_teos10( T, S, pressure1, pressure2, rho1, rho2, ! * (in) start - the starting point in the arrays. * ! * (in) npts - the number of values to calculate. * - real :: gsw_rho !external TEOS10 function - real :: p1_dbar,p2_dbar + real :: zp1, zp2, zt, zs integer :: j - p1_dbar = pressure1*1.0e-4 !convert pressure to dbar - p2_dbar = pressure2*1.0e-4 !convert pressure to dbar + zp1 = pressure1 * Pa2db !Convert pressure from Pascal to decibar + zp2 = pressure2 * Pa2db !Convert pressure from Pascal to decibar do j=start,start+npts-1 + !Conversions + zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? - rho1(j) = gsw_rho(S(j),T(j),p1_dbar) - rho2(j) = gsw_rho(S(j),T(j),p2_dbar) + rho1(j) = gsw_rho(zs,zt,zp1) + rho2(j) = gsw_rho(zs,zt,zp2) enddo end subroutine calculate_2_densities_teos10 diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index a43c4ab085..0975c2c0ee 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2284,7 +2284,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) press(:)=tv%p_ref !Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 - call convert_temp_salt_for_TEOS10(temp_z,salt_z, press, G, kd, mask_z, eos) +! call convert_temp_salt_for_TEOS10(temp_z,salt_z, press, G, kd, mask_z, eos) do k=1, kd do j=js,je From 12f82e7acf7e8de0f5a80de74134f84c46776542 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Tue, 23 Aug 2016 13:40:24 -0400 Subject: [PATCH 4/9] Assume prognostic T&S are conservative&reference - This commit corresponds to the experiment OMp5_EOS_NEMO_modelStepsCtRs3 - In this commit we assume that the internal model variables for T&S tracers are the conservative temperature and reference salinity. - This assumption is justified by the lack of source terms in advection equations for T&S - The T&S variables are being converted at 3 different places: 1. At initialization of T&S from WOA they are converted from potential&practical to conservative&reference before the model starts. No need to convert when reading restart files but we have to ensure that restart T&S are conservative&reference. 2. Before exposing T&S to coupler/exchange we have to convert from conservative&reference to potential&practical since the coupler assumes that. 3. The archived diagnostics for T&S should be strictly potential&practical, so we need to convert from model conservative&reference to potential&practical before dianostics sends. - Still to do: Ensure SSS and SST diagnostics are converted to potential&practical ?! --- config_src/coupled_driver/ocean_model_MOM.F90 | 28 +++++++--- src/core/MOM.F90 | 53 +++++++++++++++++-- src/equation_of_state/MOM_EOS_NEMO.F90 | 18 +++---- src/equation_of_state/MOM_EOS_TEOS10.F90 | 23 ++++---- .../TEOS10/gsw_gibbs_pt0_pt0.f90 | 1 + .../TEOS10/gsw_pt_from_ct.f90 | 1 + .../TEOS10/gsw_sp_from_sr.f90 | 1 + .../MOM_state_initialization.F90 | 4 +- 8 files changed, 97 insertions(+), 32 deletions(-) create mode 120000 src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 create mode 120000 src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index ea4942cbe8..f4aa4e1c20 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -229,7 +229,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in) OS%grid => OS%MOM_CSp%G ; OS%GV => OS%MOM_CSp%GV OS%C_p = OS%MOM_CSp%tv%C_p OS%fluxes%C_p = OS%MOM_CSp%tv%C_p - + ! Read all relevant parameters and write them to the model log. call log_version(param_file, mod, version, "") call get_param(param_file, mod, "RESTART_CONTROL", OS%Restart_control, & @@ -454,7 +454,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! Translate state into Ocean. ! call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & ! Ice_ocean_boundary%p, OS%press_to_z) - call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid) + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, OS%MOM_CSp%use_conT_refS) call callTree_leave("update_ocean_model()") end subroutine update_ocean_model @@ -606,10 +606,12 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, maskmap) end subroutine initialize_ocean_public_type -subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, patm, press_to_z) +subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_refS, patm, press_to_z) + use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct type(surface), intent(inout) :: state type(ocean_public_type), target, intent(inout) :: Ocean_sfc type(ocean_grid_type), intent(inout) :: G + logical, intent(in) :: use_conT_refS real, optional, intent(in) :: patm(:,:) real, optional, intent(in) :: press_to_z ! This subroutine translates the coupler's ocean_data_type into MOM's @@ -634,9 +636,23 @@ subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, patm, press_to_z) endif i0 = is - isc_bnd ; j0 = js - jsc_bnd + !If directed convert the surface T&S + !from conservative T to potential T and + !from absolute (reference) salinity to practical salinity + ! + if(use_conT_refS) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%s_surf(i,j) = gsw_sp_from_sr(state%SSS(i+i0,j+j0)) + Ocean_sfc%t_surf(i,j) = gsw_pt_from_ct(state%SSS(i+i0,j+j0),state%SST(i+i0,j+j0)) + CELSIUS_KELVIN_OFFSET + enddo ; enddo + else + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%t_surf(i,j) = state%SST(i+i0,j+j0) + CELSIUS_KELVIN_OFFSET + Ocean_sfc%s_surf(i,j) = state%SSS(i+i0,j+j0) + enddo ; enddo + endif + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%t_surf(i,j) = state%SST(i+i0,j+j0) + CELSIUS_KELVIN_OFFSET - Ocean_sfc%s_surf(i,j) = state%SSS(i+i0,j+j0) Ocean_sfc%sea_lev(i,j) = state%sea_lev(i+i0,j+j0) if (present(patm)) & Ocean_sfc%sea_lev(i,j) = Ocean_sfc%sea_lev(i,j) + patm(i,j) * press_to_z @@ -690,7 +706,7 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc) OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%ave_ssh,& OS%grid, OS%GV, OS%MOM_CSp) - call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid) + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, OS%MOM_CSp%use_conT_refS) end subroutine ocean_model_init_sfc ! diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f888599b6e..9abd893e0a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -120,6 +120,7 @@ module MOM use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd use MOM_verticalGrid, only : get_thickness_units, get_flux_units, get_tr_flux_units use MOM_wave_speed, only : wave_speed_init, wave_speed_CS +use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct implicit none ; private @@ -177,6 +178,10 @@ module MOM !! with nkml sublayers and nkbl buffer layer. logical :: diabatic_first !< If true, apply diabatic and thermodynamic !! processes before time stepping the dynamics. + logical :: use_conT_refS !< If true, , the prognostics T&S are the conservative temperature + !! and reference salinity. Care should be taken to convert them + !! to potential temperature and practical salinity before + !! exchanging them with the coupler and/or reporting T&S diagnostics. logical :: thickness_diffuse !< If true, diffuse interface height w/ a diffusivity KHTH. logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics. logical :: mixedlayer_restrat !< If true, use submesoscale mixed layer restratifying scheme. @@ -265,6 +270,8 @@ module MOM integer :: id_h = -1 integer :: id_T = -1 integer :: id_S = -1 + integer :: id_Tcon = -1 + integer :: id_Sref = -1 ! 2-d surface and bottom fields integer :: id_zos = -1 @@ -282,6 +289,8 @@ module MOM integer :: id_ssh_inst = -1 integer :: id_tob = -1 integer :: id_sob = -1 + integer :: id_consst = -1 + integer :: id_refsss = -1 ! heat and salt flux fields integer :: id_fraz = -1 @@ -462,6 +471,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) h ! h : layer thickness (meter (Bouss) or kg/m2 (non-Bouss)) real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_predia, eta_preale + real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: potTemp, pracSal real :: tot_wt_ssh, Itot_wt_ssh, I_time_int real :: zos_area_mean, volo, ssh_ga @@ -1118,11 +1128,31 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! post some diagnostics - if (CS%id_T > 0) call post_data(CS%id_T, CS%tv%T, CS%diag) - if (CS%id_S > 0) call post_data(CS%id_S, CS%tv%S, CS%diag) + ! + !In case the internal representation of the model is conservative temperature and reference salinity + !convert, for diagnostics + !from conservative temp to potential temp and + !from reference salinity to practical salinity + ! + if(CS%use_conT_refS) then + if (CS%id_Tcon > 0) call post_data(CS%id_Tcon, CS%tv%T, CS%diag) + if (CS%id_Sref > 0) call post_data(CS%id_Sref, CS%tv%S, CS%diag) + !Conversions + do k=1,nz ; do j=js,je ; do i=is,ie + pracSal(i,j,k) = gsw_sp_from_sr(CS%tv%S(i,j,k)) + potTemp(i,j,k) = gsw_pt_from_ct(CS%tv%S(i,j,k),CS%tv%T(i,j,k)) + enddo; enddo ; enddo + if (CS%id_T > 0) call post_data(CS%id_T, potTemp, CS%diag) + if (CS%id_S > 0) call post_data(CS%id_S, pracSal, CS%diag) + if (CS%id_tob > 0) call post_data(CS%id_tob, potTemp(:,:,G%ke), CS%diag, mask=G%mask2dT) + if (CS%id_sob > 0) call post_data(CS%id_sob, pracSal(:,:,G%ke), CS%diag, mask=G%mask2dT) + else + if (CS%id_T > 0) call post_data(CS%id_T, CS%tv%T, CS%diag) + if (CS%id_S > 0) call post_data(CS%id_S, CS%tv%S, CS%diag) - if (CS%id_tob > 0) call post_data(CS%id_tob, CS%tv%T(:,:,G%ke), CS%diag, mask=G%mask2dT) - if (CS%id_sob > 0) call post_data(CS%id_sob, CS%tv%S(:,:,G%ke), CS%diag, mask=G%mask2dT) + if (CS%id_tob > 0) call post_data(CS%id_tob, CS%tv%T(:,:,G%ke), CS%diag, mask=G%mask2dT) + if (CS%id_sob > 0) call post_data(CS%id_sob, CS%tv%S(:,:,G%ke), CS%diag, mask=G%mask2dT) + endif if (CS%id_Tadx > 0) call post_data(CS%id_Tadx, CS%T_adx, CS%diag) if (CS%id_Tady > 0) call post_data(CS%id_Tady, CS%T_ady, CS%diag) @@ -1486,6 +1516,12 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) "If true, apply diabatic and thermodynamic processes, \n"//& "including buoyancy forcing and mass gain or loss, \n"//& "before stepping the dynamics forward.", default=.false.) + call get_param(param_file, "MOM", "USE_CONTEMP_REFSAL", CS%use_conT_refS, & + "If true, , the prognostics T&S are the conservative temperature \n"//& + "and reference salinity. Care should be taken to convert them \n"//& + "to potential temperature and practical salinity before \n"//& + "exchanging them with the coupler and/or reporting T&S diagnostics. \n"& + , default=.false.) call get_param(param_file, "MOM", "ADIABATIC", CS%adiabatic, & "There are no diapycnal mass fluxes if ADIABATIC is \n"//& "true. This assumes that KD = KDML = 0.0 and that \n"//& @@ -2224,6 +2260,15 @@ subroutine register_diags(Time, G, GV, CS, ADp) cmor_long_name='Square of Sea Surface Salinity ', cmor_units='ppt^2', & cmor_standard_name='square_of_sea_surface_salinity') if (CS%id_sss_sq > 0) call safe_alloc_ptr(CS%SSS_sq,isd,ied,jsd,jed) + CS%id_Tcon = register_diag_field('ocean_model', 'contemp', diag%axesTL, Time, & + 'Conservative Temperature', 'Celsius') + CS%id_Sref = register_diag_field('ocean_model', 'refsalt', diag%axesTL, Time, & + long_name='Reference Salinity', units='g/Kg') + CS%id_consst = register_diag_field('ocean_model', 'conSST', diag%axesT1, Time, & + 'Sea Surface Conservative Temperature', 'Celsius', CS%missing) + CS%id_refsss = register_diag_field('ocean_model', 'refSSS', diag%axesT1, Time, & + 'Sea Surface Reference Salinity', 'g/Kg', CS%missing) + endif if (CS%use_temperature .and. CS%use_frazil) then diff --git a/src/equation_of_state/MOM_EOS_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 index abeb5b86a7..6edbe3d352 100644 --- a/src/equation_of_state/MOM_EOS_NEMO.F90 +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -29,7 +29,7 @@ module MOM_EOS_NEMO !* 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_sr_from_sp, gsw_ct_from_pt use gsw_mod_toolbox, only : gsw_rho_first_derivatives @@ -229,8 +229,8 @@ subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts) do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? @@ -288,8 +288,8 @@ subroutine calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? @@ -370,8 +370,8 @@ subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) ! do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j)) @@ -400,8 +400,8 @@ subroutine calculate_2_densities_nemo( T, S, pressure1, pressure2, rho1, rho2, s do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? !The following algorithm was provided by Roquet in a private communication. diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 index 5e67780790..c45e41e917 100644 --- a/src/equation_of_state/MOM_EOS_TEOS10.F90 +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -1,3 +1,4 @@ + module MOM_EOS_TEOS10 !*********************************************************************** !* GNU General Public License * @@ -24,7 +25,7 @@ module MOM_EOS_TEOS10 !* sea water using the TEOS10 functions * !*********************************************************************** -use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt +!use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt use gsw_mod_toolbox, only : gsw_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives implicit none ; private @@ -95,8 +96,8 @@ subroutine calculate_density_array_teos10(T, S, pressure, rho, start, npts) do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? @@ -122,8 +123,8 @@ subroutine calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, sta do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS(j), drho_dct=drho_dT(j)) @@ -149,8 +150,8 @@ subroutine calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? call gsw_specvol_first_derivatives(zs,zt,zp, v_sa=dSV_dS(j), v_ct=dSV_dT(j)) @@ -183,8 +184,8 @@ subroutine calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts) do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? rho(j) = gsw_rho(zs,zt,zp) @@ -214,8 +215,8 @@ subroutine calculate_2_densities_teos10( T, S, pressure1, pressure2, rho1, rho2, do j=start,start+npts-1 !Conversions - zs = gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity - zt = gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + 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 if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? rho1(j) = gsw_rho(zs,zt,zp1) rho2(j) = gsw_rho(zs,zt,zp2) diff --git a/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 b/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 new file mode 120000 index 0000000000..e345379f5d --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_gibbs_pt0_pt0.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 b/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 new file mode 120000 index 0000000000..cd794a1316 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_pt_from_ct.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 b/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 new file mode 120000 index 0000000000..d8cd41f4bf --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 @@ -0,0 +1 @@ +../../../pkg/GSW-Fortran/toolbox/gsw_sp_from_sr.f90 \ No newline at end of file diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0975c2c0ee..1b2fa4ed1c 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2283,8 +2283,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) press(:)=tv%p_ref - !Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 -! call convert_temp_salt_for_TEOS10(temp_z,salt_z, press, G, kd, mask_z, eos) + !Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO + call convert_temp_salt_for_TEOS10(temp_z,salt_z, press, G, kd, mask_z, eos) do k=1, kd do j=js,je From 78a23531490f16f860e710ee4f814413f42144a2 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Tue, 18 Oct 2016 14:36:08 -0400 Subject: [PATCH 5/9] change symbolic link files to copies - fremake list_paths tool refuses to discover symbolic link files. Make copies of TEOS-10 files instead of using links --- .../TEOS10/gsw_ct_from_pt.f90 | 53 +- .../TEOS10/gsw_gibbs_pt0_pt0.f90 | 48 +- .../TEOS10/gsw_mod_kinds.f90 | 17 +- .../TEOS10/gsw_mod_specvol_coefficients.f90 | 314 +++- .../TEOS10/gsw_mod_teos10_constants.f90 | 72 +- .../TEOS10/gsw_mod_toolbox.f90 | 1494 ++++++++++++++++- .../TEOS10/gsw_pt_from_ct.f90 | 73 +- src/equation_of_state/TEOS10/gsw_rho.f90 | 37 +- .../TEOS10/gsw_rho_first_derivatives.f90 | 111 +- .../TEOS10/gsw_sp_from_sr.f90 | 31 +- src/equation_of_state/TEOS10/gsw_specvol.f90 | 53 +- .../TEOS10/gsw_specvol_first_derivatives.f90 | 105 +- .../TEOS10/gsw_sr_from_sp.f90 | 31 +- 13 files changed, 2426 insertions(+), 13 deletions(-) mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_mod_kinds.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_rho.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_specvol.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 mode change 120000 => 100644 src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 diff --git a/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 b/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 deleted file mode 120000 index d67d2df3e2..0000000000 --- a/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_ct_from_pt.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 b/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 new file mode 100644 index 0000000000..c4a624ed37 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_ct_from_pt.f90 @@ -0,0 +1,52 @@ +!========================================================================== +elemental function gsw_ct_from_pt (sa, pt) +!========================================================================== +! +! Calculates Conservative Temperature from potential temperature of seawater +! +! sa : Absolute Salinity [g/kg] +! pt : potential temperature with [deg C] +! reference pressure of 0 dbar +! +! gsw_ct_from_pt : Conservative Temperature [deg C] +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : gsw_cp0, gsw_sfac + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, pt + +real (r8) :: gsw_ct_from_pt + +real (r8) :: pot_enthalpy, x2, x, y + +x2 = gsw_sfac*sa +x = sqrt(x2) +y = pt*0.025_r8 ! normalize for F03 and F08 + +pot_enthalpy = 61.01362420681071_r8 + y*(168776.46138048015_r8 + & + y*(-2735.2785605119625_r8 + y*(2574.2164453821433_r8 + & + y*(-1536.6644434977543_r8 + y*(545.7340497931629_r8 + & + (-50.91091728474331_r8 - 18.30489878927802_r8*y)*y))))) + & + x2*(268.5520265845071_r8 + y*(-12019.028203559312_r8 + & + y*(3734.858026725145_r8 + y*(-2046.7671145057618_r8 + & + y*(465.28655623826234_r8 + (-0.6370820302376359_r8 - & + 10.650848542359153_r8*y)*y)))) + & + x*(937.2099110620707_r8 + y*(588.1802812170108_r8 + & + y*(248.39476522971285_r8 + (-3.871557904936333_r8 - & + 2.6268019854268356_r8*y)*y)) + & + x*(-1687.914374187449_r8 + x*(246.9598888781377_r8 + & + x*(123.59576582457964_r8 - 48.5891069025409_r8*x)) + & + y*(936.3206544460336_r8 + & + y*(-942.7827304544439_r8 + y*(369.4389437509002_r8 + & + (-33.83664947895248_r8 - 9.987880382780322_r8*y)*y)))))) + +gsw_ct_from_pt = pot_enthalpy/gsw_cp0 + +return +end function + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 b/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 deleted file mode 120000 index e345379f5d..0000000000 --- a/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_gibbs_pt0_pt0.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 b/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 new file mode 100644 index 0000000000..6e8bcfc779 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_gibbs_pt0_pt0.f90 @@ -0,0 +1,47 @@ +!========================================================================== +elemental function gsw_gibbs_pt0_pt0 (sa, pt0) +!========================================================================== +! +! gibbs_tt at (sa,pt,0) +! +! sa : Absolute Salinity [g/kg] +! pt0 : potential temperature [deg C] +! +! gsw_gibbs_pt0_pt0 : gibbs_tt at (sa,pt,0) +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : gsw_sfac + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, pt0 + +real (r8) :: gsw_gibbs_pt0_pt0 + +real (r8) :: x2, x, y, g03, g08 + +x2 = gsw_sfac*sa +x = sqrt(x2) +y = pt0*0.025_r8 + +g03 = -24715.571866078_r8 + & + y*(4420.4472249096725_r8 + & + y*(-1778.231237203896_r8 + & + y*(1160.5182516851419_r8 + & + y*(-569.531539542516_r8 + y*128.13429152494615_r8)))) + +g08 = x2*(1760.062705994408_r8 + x*(-86.1329351956084_r8 + & + x*(-137.1145018408982_r8 + y*(296.20061691375236_r8 + & + y*(-205.67709290374563_r8 + 49.9394019139016_r8*y))) + & + y*(-60.136422517125_r8 + y*10.50720794170734_r8)) + & + y*(-1351.605895580406_r8 + y*(1097.1125373015109_r8 + & + y*(-433.20648175062206_r8 + 63.905091254154904_r8*y)))) + +gsw_gibbs_pt0_pt0 = (g03 + g08)*0.000625_r8 + +return +end function + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 b/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 deleted file mode 120000 index fa0926e540..0000000000 --- a/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/modules/gsw_mod_kinds.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 b/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 new file mode 100644 index 0000000000..7a2a80891f --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_kinds.f90 @@ -0,0 +1,16 @@ +!========================================================================== +module gsw_mod_kinds +!========================================================================== + +implicit none + +integer, parameter :: r4 = selected_real_kind(6,30) + +integer, parameter :: r8 = selected_real_kind(14,30) + +end module + +!-------------------------------------------------------------------------- + + + diff --git a/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 b/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 deleted file mode 120000 index 934f689c20..0000000000 --- a/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/modules/gsw_mod_specvol_coefficients.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 b/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 new file mode 100644 index 0000000000..7bc89c7b5e --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_specvol_coefficients.f90 @@ -0,0 +1,313 @@ +!========================================================================== +module gsw_mod_specvol_coefficients +!========================================================================== + +use gsw_mod_kinds + +implicit none + +real (r8), parameter :: a000 = -1.56497346750e-5_r8 +real (r8), parameter :: a001 = 1.85057654290e-5_r8 +real (r8), parameter :: a002 = -1.17363867310e-6_r8 +real (r8), parameter :: a003 = -3.65270065530e-7_r8 +real (r8), parameter :: a004 = 3.14540999020e-7_r8 +real (r8), parameter :: a010 = 5.55242129680e-5_r8 +real (r8), parameter :: a011 = -2.34332137060e-5_r8 +real (r8), parameter :: a012 = 4.26100574800e-6_r8 +real (r8), parameter :: a013 = 5.73918103180e-7_r8 +real (r8), parameter :: a020 = -4.95634777770e-5_r8 +real (r8), parameter :: a021 = 2.37838968519e-5_r8 +real (r8), parameter :: a022 = -1.38397620111e-6_r8 +real (r8), parameter :: a030 = 2.76445290808e-5_r8 +real (r8), parameter :: a031 = -1.36408749928e-5_r8 +real (r8), parameter :: a032 = -2.53411666056e-7_r8 +real (r8), parameter :: a040 = -4.02698077700e-6_r8 +real (r8), parameter :: a041 = 2.53683834070e-6_r8 +real (r8), parameter :: a050 = 1.23258565608e-6_r8 +real (r8), parameter :: a100 = 3.50095997640e-5_r8 +real (r8), parameter :: a101 = -9.56770881560e-6_r8 +real (r8), parameter :: a102 = -5.56991545570e-6_r8 +real (r8), parameter :: a103 = -2.72956962370e-7_r8 +real (r8), parameter :: a110 = -7.48716846880e-5_r8 +real (r8), parameter :: a111 = -4.73566167220e-7_r8 +real (r8), parameter :: a112 = 7.82747741600e-7_r8 +real (r8), parameter :: a120 = 7.24244384490e-5_r8 +real (r8), parameter :: a121 = -1.03676320965e-5_r8 +real (r8), parameter :: a122 = 2.32856664276e-8_r8 +real (r8), parameter :: a130 = -3.50383492616e-5_r8 +real (r8), parameter :: a131 = 5.18268711320e-6_r8 +real (r8), parameter :: a140 = -1.65263794500e-6_r8 +real (r8), parameter :: a200 = -4.35926785610e-5_r8 +real (r8), parameter :: a201 = 1.11008347650e-5_r8 +real (r8), parameter :: a202 = 5.46207488340e-6_r8 +real (r8), parameter :: a210 = 7.18156455200e-5_r8 +real (r8), parameter :: a211 = 5.85666925900e-6_r8 +real (r8), parameter :: a212 = -1.31462208134e-6_r8 +real (r8), parameter :: a220 = -4.30608991440e-5_r8 +real (r8), parameter :: a221 = 9.49659182340e-7_r8 +real (r8), parameter :: a230 = 1.74814722392e-5_r8 +real (r8), parameter :: a300 = 3.45324618280e-5_r8 +real (r8), parameter :: a301 = -9.84471178440e-6_r8 +real (r8), parameter :: a302 = -1.35441856270e-6_r8 +real (r8), parameter :: a310 = -3.73971683740e-5_r8 +real (r8), parameter :: a311 = -9.76522784000e-7_r8 +real (r8), parameter :: a320 = 6.85899736680e-6_r8 +real (r8), parameter :: a400 = -1.19594097880e-5_r8 +real (r8), parameter :: a401 = 2.59092252600e-6_r8 +real (r8), parameter :: a410 = 7.71906784880e-6_r8 +real (r8), parameter :: a500 = 1.38645945810e-6_r8 + +real (r8), parameter :: b000 = -3.10389819760e-4_r8 +real (r8), parameter :: b003 = 3.63101885150e-7_r8 +real (r8), parameter :: b004 = -1.11471254230e-7_r8 +real (r8), parameter :: b010 = 3.50095997640e-5_r8 +real (r8), parameter :: b013 = -2.72956962370e-7_r8 +real (r8), parameter :: b020 = -3.74358423440e-5_r8 +real (r8), parameter :: b030 = 2.41414794830e-5_r8 +real (r8), parameter :: b040 = -8.75958731540e-6_r8 +real (r8), parameter :: b050 = -3.30527589000e-7_r8 +real (r8), parameter :: b100 = 1.33856134076e-3_r8 +real (r8), parameter :: b103 = 3.34926075600e-8_r8 +real (r8), parameter :: b110 = -8.71853571220e-5_r8 +real (r8), parameter :: b120 = 7.18156455200e-5_r8 +real (r8), parameter :: b130 = -2.87072660960e-5_r8 +real (r8), parameter :: b140 = 8.74073611960e-6_r8 +real (r8), parameter :: b200 = -2.55143801811e-3_r8 +real (r8), parameter :: b210 = 1.03597385484e-4_r8 +real (r8), parameter :: b220 = -5.60957525610e-5_r8 +real (r8), parameter :: b230 = 6.85899736680e-6_r8 +real (r8), parameter :: b300 = 2.32344279772e-3_r8 +real (r8), parameter :: b310 = -4.78376391520e-5_r8 +real (r8), parameter :: b320 = 1.54381356976e-5_r8 +real (r8), parameter :: b400 = -1.05461852535e-3_r8 +real (r8), parameter :: b410 = 6.93229729050e-6_r8 +real (r8), parameter :: b500 = 1.91594743830e-4_r8 +real (r8), parameter :: b001 = 2.42624687470e-5_r8 +real (r8), parameter :: b011 = -9.56770881560e-6_r8 +real (r8), parameter :: b021 = -2.36783083610e-7_r8 +real (r8), parameter :: b031 = -3.45587736550e-6_r8 +real (r8), parameter :: b041 = 1.29567177830e-6_r8 +real (r8), parameter :: b101 = -6.95849219480e-5_r8 +real (r8), parameter :: b111 = 2.22016695300e-5_r8 +real (r8), parameter :: b121 = 5.85666925900e-6_r8 +real (r8), parameter :: b131 = 6.33106121560e-7_r8 +real (r8), parameter :: b201 = 1.12412331915e-4_r8 +real (r8), parameter :: b211 = -2.95341353532e-5_r8 +real (r8), parameter :: b221 = -1.46478417600e-6_r8 +real (r8), parameter :: b301 = -6.92888744480e-5_r8 +real (r8), parameter :: b311 = 1.03636901040e-5_r8 +real (r8), parameter :: b401 = 1.54637136265e-5_r8 +real (r8), parameter :: b002 = -5.84844329840e-7_r8 +real (r8), parameter :: b012 = -5.56991545570e-6_r8 +real (r8), parameter :: b022 = 3.91373870800e-7_r8 +real (r8), parameter :: b032 = 7.76188880920e-9_r8 +real (r8), parameter :: b102 = -9.62445031940e-6_r8 +real (r8), parameter :: b112 = 1.09241497668e-5_r8 +real (r8), parameter :: b122 = -1.31462208134e-6_r8 +real (r8), parameter :: b202 = 1.47789320994e-5_r8 +real (r8), parameter :: b212 = -4.06325568810e-6_r8 +real (r8), parameter :: b302 = -7.12478989080e-6_r8 + +real (r8), parameter :: c000 = -6.07991438090e-5_r8 +real (r8), parameter :: c001 = 1.99712338438e-5_r8 +real (r8), parameter :: c002 = -3.39280843110e-6_r8 +real (r8), parameter :: c003 = 4.21246123200e-7_r8 +real (r8), parameter :: c004 = -6.32363064300e-8_r8 +real (r8), parameter :: c005 = 1.17681023580e-8_r8 +real (r8), parameter :: c010 = 1.85057654290e-5_r8 +real (r8), parameter :: c011 = -2.34727734620e-6_r8 +real (r8), parameter :: c012 = -1.09581019659e-6_r8 +real (r8), parameter :: c013 = 1.25816399608e-6_r8 +real (r8), parameter :: c020 = -1.17166068530e-5_r8 +real (r8), parameter :: c021 = 4.26100574800e-6_r8 +real (r8), parameter :: c022 = 8.60877154770e-7_r8 +real (r8), parameter :: c030 = 7.92796561730e-6_r8 +real (r8), parameter :: c031 = -9.22650800740e-7_r8 +real (r8), parameter :: c040 = -3.41021874820e-6_r8 +real (r8), parameter :: c041 = -1.26705833028e-7_r8 +real (r8), parameter :: c050 = 5.07367668140e-7_r8 +real (r8), parameter :: c100 = 2.42624687470e-5_r8 +real (r8), parameter :: c101 = -1.16968865968e-6_r8 +real (r8), parameter :: c102 = 1.08930565545e-6_r8 +real (r8), parameter :: c103 = -4.45885016920e-7_r8 +real (r8), parameter :: c110 = -9.56770881560e-6_r8 +real (r8), parameter :: c111 = -1.11398309114e-5_r8 +real (r8), parameter :: c112 = -8.18870887110e-7_r8 +real (r8), parameter :: c120 = -2.36783083610e-7_r8 +real (r8), parameter :: c121 = 7.82747741600e-7_r8 +real (r8), parameter :: c130 = -3.45587736550e-6_r8 +real (r8), parameter :: c131 = 1.55237776184e-8_r8 +real (r8), parameter :: c140 = 1.29567177830e-6_r8 +real (r8), parameter :: c200 = -3.47924609740e-5_r8 +real (r8), parameter :: c201 = -9.62445031940e-6_r8 +real (r8), parameter :: c202 = 5.02389113400e-8_r8 +real (r8), parameter :: c210 = 1.11008347650e-5_r8 +real (r8), parameter :: c211 = 1.09241497668e-5_r8 +real (r8), parameter :: c220 = 2.92833462950e-6_r8 +real (r8), parameter :: c221 = -1.31462208134e-6_r8 +real (r8), parameter :: c230 = 3.16553060780e-7_r8 +real (r8), parameter :: c300 = 3.74707773050e-5_r8 +real (r8), parameter :: c301 = 9.85262139960e-6_r8 +real (r8), parameter :: c310 = -9.84471178440e-6_r8 +real (r8), parameter :: c311 = -2.70883712540e-6_r8 +real (r8), parameter :: c320 = -4.88261392000e-7_r8 +real (r8), parameter :: c400 = -1.73222186120e-5_r8 +real (r8), parameter :: c401 = -3.56239494540e-6_r8 +real (r8), parameter :: c410 = 2.59092252600e-6_r8 +real (r8), parameter :: c500 = 3.09274272530e-6_r8 + +real (r8), parameter :: h001 = 1.07699958620e-3_r8 +real (r8), parameter :: h002 = -3.03995719050e-5_r8 +real (r8), parameter :: h003 = 3.32853897400e-6_r8 +real (r8), parameter :: h004 = -2.82734035930e-7_r8 +real (r8), parameter :: h005 = 2.10623061600e-8_r8 +real (r8), parameter :: h006 = -2.10787688100e-9_r8 +real (r8), parameter :: h007 = 2.80192913290e-10_r8 +real (r8), parameter :: h011 = -1.56497346750e-5_r8 +real (r8), parameter :: h012 = 9.25288271450e-6_r8 +real (r8), parameter :: h013 = -3.91212891030e-7_r8 +real (r8), parameter :: h014 = -9.13175163830e-8_r8 +real (r8), parameter :: h015 = 6.29081998040e-8_r8 +real (r8), parameter :: h021 = 2.77621064840e-5_r8 +real (r8), parameter :: h022 = -5.85830342650e-6_r8 +real (r8), parameter :: h023 = 7.10167624670e-7_r8 +real (r8), parameter :: h024 = 7.17397628980e-8_r8 +real (r8), parameter :: h031 = -1.65211592590e-5_r8 +real (r8), parameter :: h032 = 3.96398280870e-6_r8 +real (r8), parameter :: h033 = -1.53775133460e-7_r8 +real (r8), parameter :: h042 = -1.70510937410e-6_r8 +real (r8), parameter :: h043 = -2.11176388380e-8_r8 +real (r8), parameter :: h041 = 6.91113227020e-6_r8 +real (r8), parameter :: h051 = -8.05396155400e-7_r8 +real (r8), parameter :: h052 = 2.53683834070e-7_r8 +real (r8), parameter :: h061 = 2.05430942680e-7_r8 +real (r8), parameter :: h101 = -3.10389819760e-4_r8 +real (r8), parameter :: h102 = 1.21312343735e-5_r8 +real (r8), parameter :: h103 = -1.94948109950e-7_r8 +real (r8), parameter :: h104 = 9.07754712880e-8_r8 +real (r8), parameter :: h105 = -2.22942508460e-8_r8 +real (r8), parameter :: h111 = 3.50095997640e-5_r8 +real (r8), parameter :: h112 = -4.78385440780e-6_r8 +real (r8), parameter :: h113 = -1.85663848520e-6_r8 +real (r8), parameter :: h114 = -6.82392405930e-8_r8 +real (r8), parameter :: h121 = -3.74358423440e-5_r8 +real (r8), parameter :: h122 = -1.18391541805e-7_r8 +real (r8), parameter :: h123 = 1.30457956930e-7_r8 +real (r8), parameter :: h131 = 2.41414794830e-5_r8 +real (r8), parameter :: h132 = -1.72793868275e-6_r8 +real (r8), parameter :: h133 = 2.58729626970e-9_r8 +real (r8), parameter :: h141 = -8.75958731540e-6_r8 +real (r8), parameter :: h142 = 6.47835889150e-7_r8 +real (r8), parameter :: h151 = -3.30527589000e-7_r8 +real (r8), parameter :: h201 = 6.69280670380e-4_r8 +real (r8), parameter :: h202 = -1.73962304870e-5_r8 +real (r8), parameter :: h203 = -1.60407505320e-6_r8 +real (r8), parameter :: h204 = 4.18657594500e-9_r8 +real (r8), parameter :: h211 = -4.35926785610e-5_r8 +real (r8), parameter :: h212 = 5.55041738250e-6_r8 +real (r8), parameter :: h213 = 1.82069162780e-6_r8 +real (r8), parameter :: h221 = 3.59078227600e-5_r8 +real (r8), parameter :: h222 = 1.46416731475e-6_r8 +real (r8), parameter :: h223 = -2.19103680220e-7_r8 +real (r8), parameter :: h231 = -1.43536330480e-5_r8 +real (r8), parameter :: h232 = 1.58276530390e-7_r8 +real (r8), parameter :: h241 = 4.37036805980e-6_r8 +real (r8), parameter :: h301 = -8.50479339370e-4_r8 +real (r8), parameter :: h302 = 1.87353886525e-5_r8 +real (r8), parameter :: h303 = 1.64210356660e-6_r8 +real (r8), parameter :: h311 = 3.45324618280e-5_r8 +real (r8), parameter :: h312 = -4.92235589220e-6_r8 +real (r8), parameter :: h313 = -4.51472854230e-7_r8 +real (r8), parameter :: h321 = -1.86985841870e-5_r8 +real (r8), parameter :: h322 = -2.44130696000e-7_r8 +real (r8), parameter :: h331 = 2.28633245560e-6_r8 +real (r8), parameter :: h401 = 5.80860699430e-4_r8 +real (r8), parameter :: h402 = -8.66110930600e-6_r8 +real (r8), parameter :: h403 = -5.93732490900e-7_r8 +real (r8), parameter :: h411 = -1.19594097880e-5_r8 +real (r8), parameter :: h421 = 3.85953392440e-6_r8 +real (r8), parameter :: h412 = 1.29546126300e-6_r8 +real (r8), parameter :: h501 = -2.10923705070e-4_r8 +real (r8), parameter :: h502 = 1.54637136265e-6_r8 +real (r8), parameter :: h511 = 1.38645945810e-6_r8 +real (r8), parameter :: h601 = 3.19324573050e-5_r8 + +real (r8), parameter :: v000 = 1.0769995862e-3_r8 +real (r8), parameter :: v001 = -6.0799143809e-5_r8 +real (r8), parameter :: v002 = 9.9856169219e-6_r8 +real (r8), parameter :: v003 = -1.1309361437e-6_r8 +real (r8), parameter :: v004 = 1.0531153080e-7_r8 +real (r8), parameter :: v005 = -1.2647261286e-8_r8 +real (r8), parameter :: v006 = 1.9613503930e-9_r8 +real (r8), parameter :: v010 = -3.1038981976e-4_r8 +real (r8), parameter :: v011 = 2.4262468747e-5_r8 +real (r8), parameter :: v012 = -5.8484432984e-7_r8 +real (r8), parameter :: v013 = 3.6310188515e-7_r8 +real (r8), parameter :: v014 = -1.1147125423e-7_r8 +real (r8), parameter :: v020 = 6.6928067038e-4_r8 +real (r8), parameter :: v021 = -3.4792460974e-5_r8 +real (r8), parameter :: v022 = -4.8122251597e-6_r8 +real (r8), parameter :: v023 = 1.6746303780e-8_r8 +real (r8), parameter :: v030 = -8.5047933937e-4_r8 +real (r8), parameter :: v031 = 3.7470777305e-5_r8 +real (r8), parameter :: v032 = 4.9263106998e-6_r8 +real (r8), parameter :: v040 = 5.8086069943e-4_r8 +real (r8), parameter :: v041 = -1.7322218612e-5_r8 +real (r8), parameter :: v042 = -1.7811974727e-6_r8 +real (r8), parameter :: v050 = -2.1092370507e-4_r8 +real (r8), parameter :: v051 = 3.0927427253e-6_r8 +real (r8), parameter :: v060 = 3.1932457305e-5_r8 +real (r8), parameter :: v100 = -1.5649734675e-5_r8 +real (r8), parameter :: v101 = 1.8505765429e-5_r8 +real (r8), parameter :: v102 = -1.1736386731e-6_r8 +real (r8), parameter :: v103 = -3.6527006553e-7_r8 +real (r8), parameter :: v104 = 3.1454099902e-7_r8 +real (r8), parameter :: v110 = 3.5009599764e-5_r8 +real (r8), parameter :: v111 = -9.5677088156e-6_r8 +real (r8), parameter :: v112 = -5.5699154557e-6_r8 +real (r8), parameter :: v113 = -2.7295696237e-7_r8 +real (r8), parameter :: v120 = -4.3592678561e-5_r8 +real (r8), parameter :: v121 = 1.1100834765e-5_r8 +real (r8), parameter :: v122 = 5.4620748834e-6_r8 +real (r8), parameter :: v130 = 3.4532461828e-5_r8 +real (r8), parameter :: v131 = -9.8447117844e-6_r8 +real (r8), parameter :: v132 = -1.3544185627e-6_r8 +real (r8), parameter :: v140 = -1.1959409788e-5_r8 +real (r8), parameter :: v141 = 2.5909225260e-6_r8 +real (r8), parameter :: v150 = 1.3864594581e-6_r8 +real (r8), parameter :: v200 = 2.7762106484e-5_r8 +real (r8), parameter :: v201 = -1.1716606853e-5_r8 +real (r8), parameter :: v202 = 2.1305028740e-6_r8 +real (r8), parameter :: v203 = 2.8695905159e-7_r8 +real (r8), parameter :: v210 = -3.7435842344e-5_r8 +real (r8), parameter :: v211 = -2.3678308361e-7_r8 +real (r8), parameter :: v212 = 3.9137387080e-7_r8 +real (r8), parameter :: v220 = 3.5907822760e-5_r8 +real (r8), parameter :: v221 = 2.9283346295e-6_r8 +real (r8), parameter :: v222 = -6.5731104067e-7_r8 +real (r8), parameter :: v230 = -1.8698584187e-5_r8 +real (r8), parameter :: v231 = -4.8826139200e-7_r8 +real (r8), parameter :: v240 = 3.8595339244e-6_r8 +real (r8), parameter :: v300 = -1.6521159259e-5_r8 +real (r8), parameter :: v301 = 7.9279656173e-6_r8 +real (r8), parameter :: v302 = -4.6132540037e-7_r8 +real (r8), parameter :: v310 = 2.4141479483e-5_r8 +real (r8), parameter :: v311 = -3.4558773655e-6_r8 +real (r8), parameter :: v312 = 7.7618888092e-9_r8 +real (r8), parameter :: v320 = -1.4353633048e-5_r8 +real (r8), parameter :: v321 = 3.1655306078e-7_r8 +real (r8), parameter :: v330 = 2.2863324556e-6_r8 +real (r8), parameter :: v400 = 6.9111322702e-6_r8 +real (r8), parameter :: v401 = -3.4102187482e-6_r8 +real (r8), parameter :: v402 = -6.3352916514e-8_r8 +real (r8), parameter :: v410 = -8.7595873154e-6_r8 +real (r8), parameter :: v411 = 1.2956717783e-6_r8 +real (r8), parameter :: v420 = 4.3703680598e-6_r8 +real (r8), parameter :: v500 = -8.0539615540e-7_r8 +real (r8), parameter :: v501 = 5.0736766814e-7_r8 +real (r8), parameter :: v510 = -3.3052758900e-7_r8 +real (r8), parameter :: v600 = 2.0543094268e-7_r8 + +end module + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 b/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 deleted file mode 120000 index 17dec5add5..0000000000 --- a/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/modules/gsw_mod_teos10_constants.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 b/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 new file mode 100644 index 0000000000..e3c6afbce0 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_teos10_constants.f90 @@ -0,0 +1,71 @@ +!========================================================================== +module gsw_mod_teos10_constants +!========================================================================== + +use gsw_mod_kinds + +implicit none + +real (r8), parameter :: db2pa = 1.0e4_r8 +real (r8), parameter :: rec_db2pa = 1.0e-4_r8 + +real (r8), parameter :: pa2db = 1.0e-4_r8 +real (r8), parameter :: rec_pa2db = 1.0e4_r8 + +real (r8), parameter :: pi = 3.141592653589793_r8 +real (r8), parameter :: deg2rad = pi/180.0_r8 +real (r8), parameter :: rad2deg = 180.0_r8/pi + +real (r8), parameter :: gamma = 2.26e-7_r8 + +! cp0 = The "specific heat" for use [ J/(kg K) ] +! with Conservative Temperature + +real (r8), parameter :: gsw_cp0 = 3991.86795711963_r8 + +! T0 = the Celcius zero point. [ K ] + +real (r8), parameter :: gsw_t0 = 273.15_r8 + +! P0 = Absolute Pressure of one standard atmosphere. [ Pa ] + +real (r8), parameter :: gsw_p0 = 101325.0_r8 + +! SSO = Standard Ocean Reference Salinity. [ g/kg ] + +real (r8), parameter :: gsw_sso = 35.16504_r8 +real (r8), parameter :: gsw_sqrtsso = 5.930011804372737_r8 + +! uPS = unit conversion factor for salinities [ g/kg ] + +real (r8), parameter :: gsw_ups = gsw_sso/35.0_r8 + +! sfac = 1/(40*gsw_ups) + +real (r8), parameter :: gsw_sfac = 0.0248826675584615_r8 + +! deltaS = 24, offset = deltaS*gsw_sfac + +real (r8), parameter :: offset = 5.971840214030754e-1_r8 + +! C3515 = Conductivity at (SP=35, t_68=15, p=0) [ mS/cm ] + +real (r8), parameter :: gsw_c3515 = 42.9140_r8 + +! SonCl = SP to Chlorinity ratio [ (g/kg)^-1 ] + +real (r8), parameter :: gsw_soncl = 1.80655_r8 + +! valence_factor = valence factor of sea salt of Reference Composition +! [ unitless ] + +real (r8), parameter :: gsw_valence_factor = 1.2452898_r8 + +! atomic_weight = mole-weighted atomic weight of sea salt of Reference +! Composition [ g/mol ] + +real (r8), parameter :: gsw_atomic_weight = 31.4038218_r8 + +end module + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 b/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 deleted file mode 120000 index f2f4761ec4..0000000000 --- a/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/modules/gsw_mod_toolbox.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 b/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 new file mode 100644 index 0000000000..a8012e1274 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_mod_toolbox.f90 @@ -0,0 +1,1493 @@ +module gsw_mod_toolbox + +use gsw_mod_kinds + +implicit none + +public :: gsw_add_barrier +public :: gsw_add_mean +public :: gsw_adiabatic_lapse_rate_from_ct +public :: gsw_adiabatic_lapse_rate_ice +public :: gsw_alpha +public :: gsw_alpha_on_beta +public :: gsw_alpha_wrt_t_exact +public :: gsw_alpha_wrt_t_ice +public :: gsw_beta_const_t_exact +public :: gsw_beta +public :: gsw_cabbeling +public :: gsw_c_from_sp +public :: gsw_chem_potential_water_ice +public :: gsw_chem_potential_water_t_exact +public :: gsw_cp_ice +public :: gsw_ct_first_derivatives +public :: gsw_ct_first_derivatives_wrt_t_exact +public :: gsw_ct_freezing_exact +public :: gsw_ct_freezing +public :: gsw_ct_freezing_first_derivatives +public :: gsw_ct_freezing_first_derivatives_poly +public :: gsw_ct_freezing_poly +public :: gsw_ct_from_enthalpy_exact +public :: gsw_ct_from_enthalpy +public :: gsw_ct_from_entropy +public :: gsw_ct_from_pt +public :: gsw_ct_from_rho +public :: gsw_ct_from_t +public :: gsw_ct_maxdensity +public :: gsw_ct_second_derivatives +public :: gsw_deltasa_atlas +public :: gsw_deltasa_from_sp +public :: gsw_dilution_coefficient_t_exact +public :: gsw_dynamic_enthalpy +public :: gsw_enthalpy_ct_exact +public :: gsw_enthalpy_diff +public :: gsw_enthalpy +public :: gsw_enthalpy_first_derivatives_ct_exact +public :: gsw_enthalpy_first_derivatives +public :: gsw_enthalpy_ice +public :: gsw_enthalpy_second_derivatives_ct_exact +public :: gsw_enthalpy_second_derivatives +public :: gsw_enthalpy_sso_0 +public :: gsw_enthalpy_t_exact +public :: gsw_entropy_first_derivatives +public :: gsw_entropy_from_pt +public :: gsw_entropy_from_t +public :: gsw_entropy_ice +public :: gsw_entropy_part +public :: gsw_entropy_part_zerop +public :: gsw_entropy_second_derivatives +public :: gsw_fdelta +public :: gsw_frazil_properties +public :: gsw_frazil_properties_potential +public :: gsw_frazil_properties_potential_poly +public :: gsw_frazil_ratios_adiabatic +public :: gsw_frazil_ratios_adiabatic_poly +public :: gsw_geo_strf_dyn_height +public :: gsw_geo_strf_dyn_height_pc +public :: gsw_gibbs +public :: gsw_gibbs_ice +public :: gsw_gibbs_ice_part_t +public :: gsw_gibbs_ice_pt0 +public :: gsw_gibbs_ice_pt0_pt0 +public :: gsw_gibbs_pt0_pt0 +public :: gsw_grav +public :: gsw_helmholtz_energy_ice +public :: gsw_hill_ratio_at_sp2 +public :: gsw_ice_fraction_to_freeze_seawater +public :: gsw_internal_energy +public :: gsw_internal_energy_ice +public :: gsw_ipv_vs_fnsquared_ratio +public :: gsw_kappa_const_t_ice +public :: gsw_kappa +public :: gsw_kappa_ice +public :: gsw_kappa_t_exact +public :: gsw_latentheat_evap_ct +public :: gsw_latentheat_evap_t +public :: gsw_latentheat_melting +public :: gsw_linear_interp_sa_ct +public :: gsw_melting_ice_equilibrium_sa_ct_ratio +public :: gsw_melting_ice_equilibrium_sa_ct_ratio_poly +public :: gsw_melting_ice_into_seawater +public :: gsw_melting_ice_sa_ct_ratio +public :: gsw_melting_ice_sa_ct_ratio_poly +public :: gsw_melting_seaice_equilibrium_sa_ct_ratio +public :: gsw_melting_seaice_equilibrium_sa_ct_ratio_poly +public :: gsw_melting_seaice_into_seawater +public :: gsw_melting_seaice_sa_ct_ratio +public :: gsw_melting_seaice_sa_ct_ratio_poly +public :: gsw_nsquared +public :: gsw_pot_enthalpy_from_pt_ice +public :: gsw_pot_enthalpy_from_pt_ice_poly +public :: gsw_pot_enthalpy_ice_freezing +public :: gsw_pot_enthalpy_ice_freezing_first_derivatives +public :: gsw_pot_enthalpy_ice_freezing_first_derivatives_poly +public :: gsw_pot_enthalpy_ice_freezing_poly +public :: gsw_pot_rho_t_exact +public :: gsw_pressure_coefficient_ice +public :: gsw_pressure_freezing_ct +public :: gsw_pt0_cold_ice_poly +public :: gsw_pt0_from_t +public :: gsw_pt0_from_t_ice +public :: gsw_pt_first_derivatives +public :: gsw_pt_from_ct +public :: gsw_pt_from_entropy +public :: gsw_pt_from_pot_enthalpy_ice +public :: gsw_pt_from_pot_enthalpy_ice_poly_dh +public :: gsw_pt_from_pot_enthalpy_ice_poly +public :: gsw_pt_from_t +public :: gsw_pt_from_t_ice +public :: gsw_pt_second_derivatives +public :: gsw_rho_alpha_beta +public :: gsw_rho +public :: gsw_rho_first_derivatives +public :: gsw_rho_first_derivatives_wrt_enthalpy +public :: gsw_rho_ice +public :: gsw_rho_second_derivatives +public :: gsw_rho_second_derivatives_wrt_enthalpy +public :: gsw_rho_t_exact +public :: gsw_rr68_interp_sa_ct +public :: gsw_saar +public :: gsw_sa_freezing_estimate +public :: gsw_sa_freezing_from_ct +public :: gsw_sa_freezing_from_ct_poly +public :: gsw_sa_freezing_from_t +public :: gsw_sa_freezing_from_t_poly +public :: gsw_sa_from_rho +public :: gsw_sa_from_sp_baltic +public :: gsw_sa_from_sp +public :: gsw_sa_from_sstar +public :: gsw_sa_p_inrange +public :: gsw_seaice_fraction_to_freeze_seawater +public :: gsw_sigma0 +public :: gsw_sigma1 +public :: gsw_sigma2 +public :: gsw_sigma3 +public :: gsw_sigma4 +public :: gsw_sound_speed +public :: gsw_sound_speed_ice +public :: gsw_sound_speed_t_exact +public :: gsw_specvol_alpha_beta +public :: gsw_specvol_anom_standard +public :: gsw_specvol +public :: gsw_specvol_first_derivatives +public :: gsw_specvol_first_derivatives_wrt_enthalpy +public :: gsw_specvol_ice +public :: gsw_specvol_second_derivatives +public :: gsw_specvol_second_derivatives_wrt_enthalpy +public :: gsw_specvol_sso_0 +public :: gsw_specvol_t_exact +public :: gsw_sp_from_c +public :: gsw_sp_from_sa_baltic +public :: gsw_sp_from_sa +public :: gsw_sp_from_sk +public :: gsw_sp_from_sr +public :: gsw_sp_from_sstar +public :: gsw_spiciness0 +public :: gsw_spiciness1 +public :: gsw_spiciness2 +public :: gsw_sr_from_sp +public :: gsw_sstar_from_sa +public :: gsw_sstar_from_sp +public :: gsw_t_deriv_chem_potential_water_t_exact +public :: gsw_t_freezing_exact +public :: gsw_t_freezing +public :: gsw_t_freezing_first_derivatives +public :: gsw_t_freezing_first_derivatives_poly +public :: gsw_t_freezing_poly +public :: gsw_t_from_ct +public :: gsw_t_from_pt0_ice +public :: gsw_thermobaric +public :: gsw_turner_rsubrho +public :: gsw_util_indx +public :: gsw_util_interp1q_int +public :: gsw_util_sort_real +public :: gsw_util_xinterp1 +public :: gsw_z_from_p + +interface + + pure subroutine gsw_add_barrier (input_data, long, lat, long_grid, & + lat_grid, dlong_grid, dlat_grid, output_data) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: long, lat, long_grid, lat_grid, dlong_grid + real (r8), intent(in) :: dlat_grid + real (r8), intent(in), dimension(4) :: input_data + real (r8), intent(out), dimension(4) :: output_data + end subroutine gsw_add_barrier + + pure subroutine gsw_add_mean (data_in, data_out) + use gsw_mod_kinds + implicit none + real (r8), intent(in), dimension(4) :: data_in + real (r8), intent(out), dimension(4) :: data_out + end subroutine gsw_add_mean + + elemental function gsw_adiabatic_lapse_rate_from_ct (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_adiabatic_lapse_rate_from_ct + end function gsw_adiabatic_lapse_rate_from_ct + + elemental function gsw_adiabatic_lapse_rate_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_adiabatic_lapse_rate_ice + end function gsw_adiabatic_lapse_rate_ice + + elemental function gsw_alpha (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_alpha + end function gsw_alpha + + elemental function gsw_alpha_on_beta (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_alpha_on_beta + end function gsw_alpha_on_beta + + elemental function gsw_alpha_wrt_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_alpha_wrt_t_exact + end function gsw_alpha_wrt_t_exact + + elemental function gsw_alpha_wrt_t_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_alpha_wrt_t_ice + end function gsw_alpha_wrt_t_ice + + elemental function gsw_beta_const_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_beta_const_t_exact + end function gsw_beta_const_t_exact + + elemental function gsw_beta (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_beta + end function gsw_beta + + elemental function gsw_cabbeling (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_cabbeling + end function gsw_cabbeling + + elemental function gsw_c_from_sp (sp, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sp, t, p + real (r8) :: gsw_c_from_sp + end function gsw_c_from_sp + + elemental function gsw_chem_potential_water_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_chem_potential_water_ice + end function gsw_chem_potential_water_ice + + elemental function gsw_chem_potential_water_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_chem_potential_water_t_exact + end function gsw_chem_potential_water_t_exact + + elemental function gsw_cp_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_cp_ice + end function gsw_cp_ice + + elemental subroutine gsw_ct_first_derivatives (sa, pt, ct_sa, ct_pt) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, pt + real (r8), intent(out), optional :: ct_sa, ct_pt + end subroutine gsw_ct_first_derivatives + + elemental subroutine gsw_ct_first_derivatives_wrt_t_exact (sa, t, p, & + ct_sa_wrt_t, ct_t_wrt_t, ct_p_wrt_t) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8), intent(out), optional :: ct_p_wrt_t, ct_sa_wrt_t, ct_t_wrt_t + end subroutine gsw_ct_first_derivatives_wrt_t_exact + + elemental function gsw_ct_freezing_exact (sa, p, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + real (r8) :: gsw_ct_freezing_exact + end function gsw_ct_freezing_exact + + elemental function gsw_ct_freezing (sa, p, saturation_fraction, poly) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + logical, intent(in), optional :: poly + real (r8) :: gsw_ct_freezing + end function gsw_ct_freezing + + elemental subroutine gsw_ct_freezing_first_derivatives (sa, p, & + saturation_fraction, ctfreezing_sa, ctfreezing_p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + real (r8), intent(out), optional :: ctfreezing_sa, ctfreezing_p + end subroutine gsw_ct_freezing_first_derivatives + + elemental subroutine gsw_ct_freezing_first_derivatives_poly (sa, p, & + saturation_fraction, ctfreezing_sa, ctfreezing_p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + real (r8), intent(out), optional :: ctfreezing_sa, ctfreezing_p + end subroutine gsw_ct_freezing_first_derivatives_poly + + elemental function gsw_ct_freezing_poly (sa, p, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + real (r8) :: gsw_ct_freezing_poly + end function gsw_ct_freezing_poly + + elemental function gsw_ct_from_enthalpy_exact (sa, h, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, h, p + real (r8) :: gsw_ct_from_enthalpy_exact + end function gsw_ct_from_enthalpy_exact + + elemental function gsw_ct_from_enthalpy (sa, h, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, h, p + real (r8) :: gsw_ct_from_enthalpy + end function gsw_ct_from_enthalpy + + elemental function gsw_ct_from_entropy (sa, entropy) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, entropy + real (r8) :: gsw_ct_from_entropy + end function gsw_ct_from_entropy + + elemental function gsw_ct_from_pt (sa, pt) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, pt + real (r8) :: gsw_ct_from_pt + end function gsw_ct_from_pt + + elemental subroutine gsw_ct_from_rho (rho, sa, p, ct, ct_multiple) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: rho, sa, p + real (r8), intent(out) :: ct + real (r8), intent(out), optional :: ct_multiple + end subroutine gsw_ct_from_rho + + elemental function gsw_ct_from_t (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_ct_from_t + end function gsw_ct_from_t + + elemental function gsw_ct_maxdensity (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_ct_maxdensity + end function gsw_ct_maxdensity + + elemental subroutine gsw_ct_second_derivatives (sa, pt, ct_sa_sa, ct_sa_pt, & + ct_pt_pt) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, pt + real (r8), intent(out), optional :: ct_sa_sa, ct_sa_pt, ct_pt_pt + end subroutine gsw_ct_second_derivatives + + elemental function gsw_deltasa_atlas (p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: p, long, lat + real (r8) :: gsw_deltasa_atlas + end function gsw_deltasa_atlas + + elemental function gsw_deltasa_from_sp (sp, p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sp, p, long, lat + real (r8) :: gsw_deltasa_from_sp + end function gsw_deltasa_from_sp + + elemental function gsw_dilution_coefficient_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_dilution_coefficient_t_exact + end function gsw_dilution_coefficient_t_exact + + elemental function gsw_dynamic_enthalpy (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_dynamic_enthalpy + end function gsw_dynamic_enthalpy + + elemental function gsw_enthalpy_ct_exact (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_enthalpy_ct_exact + end function gsw_enthalpy_ct_exact + + elemental function gsw_enthalpy_diff (sa, ct, p_shallow, p_deep) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p_shallow, p_deep + real (r8) :: gsw_enthalpy_diff + end function gsw_enthalpy_diff + + elemental function gsw_enthalpy (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_enthalpy + end function gsw_enthalpy + + elemental subroutine gsw_enthalpy_first_derivatives_ct_exact (sa, ct, p, & + h_sa, h_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: h_sa, h_ct + end subroutine gsw_enthalpy_first_derivatives_ct_exact + + elemental subroutine gsw_enthalpy_first_derivatives (sa, ct, p, h_sa, h_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: h_sa, h_ct + end subroutine gsw_enthalpy_first_derivatives + + elemental function gsw_enthalpy_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_enthalpy_ice + end function gsw_enthalpy_ice + + elemental subroutine gsw_enthalpy_second_derivatives_ct_exact (sa, ct, p, & + h_sa_sa, h_sa_ct, h_ct_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: h_sa_sa, h_sa_ct, h_ct_ct + end subroutine gsw_enthalpy_second_derivatives_ct_exact + + elemental subroutine gsw_enthalpy_second_derivatives (sa, ct, p, h_sa_sa, & + h_sa_ct, h_ct_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: h_sa_sa, h_sa_ct, h_ct_ct + end subroutine gsw_enthalpy_second_derivatives + + elemental function gsw_enthalpy_sso_0 (p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: p + real (r8) :: gsw_enthalpy_sso_0 + end function gsw_enthalpy_sso_0 + + elemental function gsw_enthalpy_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_enthalpy_t_exact + end function gsw_enthalpy_t_exact + + elemental subroutine gsw_entropy_first_derivatives (sa, ct, eta_sa, eta_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8), intent(out), optional :: eta_sa, eta_ct + end subroutine gsw_entropy_first_derivatives + + elemental function gsw_entropy_from_pt (sa, pt) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, pt + real (r8) :: gsw_entropy_from_pt + end function gsw_entropy_from_pt + + elemental function gsw_entropy_from_t (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_entropy_from_t + end function gsw_entropy_from_t + + elemental function gsw_entropy_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_entropy_ice + end function gsw_entropy_ice + + elemental function gsw_entropy_part (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_entropy_part + end function gsw_entropy_part + + elemental function gsw_entropy_part_zerop (sa, pt0) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, pt0 + real (r8) :: gsw_entropy_part_zerop + end function gsw_entropy_part_zerop + + elemental subroutine gsw_entropy_second_derivatives (sa, ct, eta_sa_sa, & + eta_sa_ct, eta_ct_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8), intent(out), optional :: eta_sa_sa, eta_sa_ct, eta_ct_ct + end subroutine gsw_entropy_second_derivatives + + elemental function gsw_fdelta (p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: p, long, lat + real (r8) :: gsw_fdelta + end function gsw_fdelta + + elemental subroutine gsw_frazil_properties (sa_bulk, h_bulk, p, & + sa_final, ct_final, w_ih_final) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa_bulk, h_bulk, p + real (r8), intent(out) :: sa_final, ct_final, w_ih_final + end subroutine gsw_frazil_properties + + elemental subroutine gsw_frazil_properties_potential (sa_bulk, h_pot_bulk,& + p, sa_final, ct_final, w_ih_final) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa_bulk, h_pot_bulk, p + real (r8), intent(out) :: sa_final, ct_final, w_ih_final + end subroutine gsw_frazil_properties_potential + + elemental subroutine gsw_frazil_properties_potential_poly (sa_bulk, & + h_pot_bulk, p, sa_final, ct_final, w_ih_final) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa_bulk, h_pot_bulk, p + real (r8), intent(out) :: sa_final, ct_final, w_ih_final + end subroutine gsw_frazil_properties_potential_poly + + elemental subroutine gsw_frazil_ratios_adiabatic (sa, p, w_ih, & + dsa_dct_frazil, dsa_dp_frazil, dct_dp_frazil) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, w_ih + real (r8), intent(out) :: dsa_dct_frazil, dsa_dp_frazil, dct_dp_frazil + end subroutine gsw_frazil_ratios_adiabatic + + elemental subroutine gsw_frazil_ratios_adiabatic_poly (sa, p, w_ih, & + dsa_dct_frazil, dsa_dp_frazil, dct_dp_frazil) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, w_ih + real (r8), intent(out) :: dsa_dct_frazil, dsa_dp_frazil, dct_dp_frazil + end subroutine gsw_frazil_ratios_adiabatic_poly + + pure function gsw_geo_strf_dyn_height (sa, ct, p, p_ref) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa(:), ct(:), p(:), p_ref + real (r8) :: gsw_geo_strf_dyn_height(size(sa)) + end function gsw_geo_strf_dyn_height + + pure subroutine gsw_geo_strf_dyn_height_pc (sa, ct, delta_p, & + geo_strf_dyn_height_pc, p_mid) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa(:), ct(:), delta_p(:) + real (r8), intent(out) :: geo_strf_dyn_height_pc(:), p_mid(:) + end subroutine gsw_geo_strf_dyn_height_pc + + elemental function gsw_gibbs (ns, nt, np, sa, t, p) + use gsw_mod_kinds + implicit none + integer, intent(in) :: ns, nt, np + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_gibbs + end function gsw_gibbs + + elemental function gsw_gibbs_ice (nt, np, t, p) + use gsw_mod_kinds + implicit none + integer, intent(in) :: nt, np + real (r8), intent(in) :: t, p + real (r8) :: gsw_gibbs_ice + end function gsw_gibbs_ice + + elemental function gsw_gibbs_ice_part_t (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_gibbs_ice_part_t + end function gsw_gibbs_ice_part_t + + elemental function gsw_gibbs_ice_pt0 (pt0) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pt0 + real (r8) :: gsw_gibbs_ice_pt0 + end function gsw_gibbs_ice_pt0 + + elemental function gsw_gibbs_ice_pt0_pt0 (pt0) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pt0 + real (r8) :: gsw_gibbs_ice_pt0_pt0 + end function gsw_gibbs_ice_pt0_pt0 + + elemental function gsw_gibbs_pt0_pt0 (sa, pt0) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, pt0 + real (r8) :: gsw_gibbs_pt0_pt0 + end function gsw_gibbs_pt0_pt0 + + elemental function gsw_grav (lat, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: lat, p + real (r8) :: gsw_grav + end function gsw_grav + + elemental function gsw_helmholtz_energy_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_helmholtz_energy_ice + end function gsw_helmholtz_energy_ice + + elemental function gsw_hill_ratio_at_sp2 (t) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t + real (r8) :: gsw_hill_ratio_at_sp2 + end function gsw_hill_ratio_at_sp2 + + elemental subroutine gsw_ice_fraction_to_freeze_seawater (sa, ct, p, & + t_ih, sa_freeze, ct_freeze, w_ih) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, t_ih + real (r8), intent(out) :: sa_freeze, ct_freeze, w_ih + end subroutine gsw_ice_fraction_to_freeze_seawater + + elemental function gsw_internal_energy (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_internal_energy + end function gsw_internal_energy + + elemental function gsw_internal_energy_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_internal_energy_ice + end function gsw_internal_energy_ice + + pure subroutine gsw_ipv_vs_fnsquared_ratio (sa, ct, p, p_ref, & + ipv_vs_fnsquared_ratio, p_mid) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa(:), ct(:), p(:), p_ref + real (r8), intent(out) :: ipv_vs_fnsquared_ratio(:), p_mid(:) + end subroutine gsw_ipv_vs_fnsquared_ratio + + elemental function gsw_kappa_const_t_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_kappa_const_t_ice + end function gsw_kappa_const_t_ice + + elemental function gsw_kappa (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_kappa + end function gsw_kappa + + elemental function gsw_kappa_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_kappa_ice + end function gsw_kappa_ice + + elemental function gsw_kappa_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_kappa_t_exact + end function gsw_kappa_t_exact + + elemental function gsw_latentheat_evap_ct (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_latentheat_evap_ct + end function gsw_latentheat_evap_ct + + elemental function gsw_latentheat_evap_t (sa, t) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t + real (r8) :: gsw_latentheat_evap_t + end function gsw_latentheat_evap_t + + elemental function gsw_latentheat_melting (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_latentheat_melting + end function gsw_latentheat_melting + + pure subroutine gsw_linear_interp_sa_ct (sa, ct, p, p_i, sa_i, ct_i) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa(:), ct(:), p(:), p_i(:) + real (r8), intent(out) :: sa_i(:), ct_i(:) + end subroutine gsw_linear_interp_sa_ct + + elemental function gsw_melting_ice_equilibrium_sa_ct_ratio (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_melting_ice_equilibrium_sa_ct_ratio + end function gsw_melting_ice_equilibrium_sa_ct_ratio + + elemental function gsw_melting_ice_equilibrium_sa_ct_ratio_poly (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_melting_ice_equilibrium_sa_ct_ratio_poly + end function gsw_melting_ice_equilibrium_sa_ct_ratio_poly + + elemental subroutine gsw_melting_ice_into_seawater (sa, ct, p, w_ih, t_ih,& + sa_final, ct_final, w_ih_final) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, w_ih, t_ih + real (r8), intent(out) :: sa_final, ct_final, w_ih_final + end subroutine gsw_melting_ice_into_seawater + + elemental function gsw_melting_ice_sa_ct_ratio (sa, ct, p, t_ih) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, t_ih + real (r8) :: gsw_melting_ice_sa_ct_ratio + end function gsw_melting_ice_sa_ct_ratio + + elemental function gsw_melting_ice_sa_ct_ratio_poly (sa, ct, p, t_ih) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, t_ih + real (r8) :: gsw_melting_ice_sa_ct_ratio_poly + end function gsw_melting_ice_sa_ct_ratio_poly + + elemental function gsw_melting_seaice_equilibrium_sa_ct_ratio (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_melting_seaice_equilibrium_sa_ct_ratio + end function gsw_melting_seaice_equilibrium_sa_ct_ratio + + elemental function gsw_melting_seaice_equilibrium_sa_ct_ratio_poly (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_melting_seaice_equilibrium_sa_ct_ratio_poly + end function gsw_melting_seaice_equilibrium_sa_ct_ratio_poly + + elemental subroutine gsw_melting_seaice_into_seawater (sa, ct, p, & + w_seaice, sa_seaice, t_seaice, sa_final, ct_final) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, w_seaice, sa_seaice, t_seaice + real (r8), intent(out) :: sa_final, ct_final + end subroutine gsw_melting_seaice_into_seawater + + elemental function gsw_melting_seaice_sa_ct_ratio (sa, ct, p, sa_seaice, & + t_seaice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, sa_seaice, t_seaice + real (r8) :: gsw_melting_seaice_sa_ct_ratio + end function gsw_melting_seaice_sa_ct_ratio + + elemental function gsw_melting_seaice_sa_ct_ratio_poly (sa, ct, p, & + sa_seaice, t_seaice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, sa_seaice, t_seaice + real (r8) :: gsw_melting_seaice_sa_ct_ratio_poly + end function gsw_melting_seaice_sa_ct_ratio_poly + + pure subroutine gsw_nsquared (sa, ct, p, lat, n2, p_mid) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa(:), ct(:), p(:), lat(:) + real (r8), intent(out) :: n2(:), p_mid(:) + end subroutine gsw_nsquared + + elemental function gsw_pot_enthalpy_from_pt_ice (pt0_ice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pt0_ice + real (r8) :: gsw_pot_enthalpy_from_pt_ice + end function gsw_pot_enthalpy_from_pt_ice + + elemental function gsw_pot_enthalpy_from_pt_ice_poly (pt0_ice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pt0_ice + real (r8) :: gsw_pot_enthalpy_from_pt_ice_poly + end function gsw_pot_enthalpy_from_pt_ice_poly + + elemental function gsw_pot_enthalpy_ice_freezing (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_pot_enthalpy_ice_freezing + end function gsw_pot_enthalpy_ice_freezing + + elemental subroutine gsw_pot_enthalpy_ice_freezing_first_derivatives (sa, & + p, pot_enthalpy_ice_freezing_sa, pot_enthalpy_ice_freezing_p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8), intent(out), optional :: pot_enthalpy_ice_freezing_sa + real (r8), intent(out), optional :: pot_enthalpy_ice_freezing_p + end subroutine gsw_pot_enthalpy_ice_freezing_first_derivatives + + elemental subroutine gsw_pot_enthalpy_ice_freezing_first_derivatives_poly(& + sa, p, pot_enthalpy_ice_freezing_sa, pot_enthalpy_ice_freezing_p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8), intent(out), optional :: pot_enthalpy_ice_freezing_sa + real (r8), intent(out), optional :: pot_enthalpy_ice_freezing_p + end subroutine gsw_pot_enthalpy_ice_freezing_first_derivatives_poly + + elemental function gsw_pot_enthalpy_ice_freezing_poly (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8) :: gsw_pot_enthalpy_ice_freezing_poly + end function gsw_pot_enthalpy_ice_freezing_poly + + elemental function gsw_pot_rho_t_exact (sa, t, p, p_ref) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p, p_ref + real (r8) :: gsw_pot_rho_t_exact + end function gsw_pot_rho_t_exact + + elemental function gsw_pressure_coefficient_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_pressure_coefficient_ice + end function gsw_pressure_coefficient_ice + + elemental function gsw_pressure_freezing_ct (sa, ct, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, saturation_fraction + real (r8) :: gsw_pressure_freezing_ct + end function gsw_pressure_freezing_ct + + elemental function gsw_pt0_cold_ice_poly (pot_enthalpy_ice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pot_enthalpy_ice + real (r8) :: gsw_pt0_cold_ice_poly + end function gsw_pt0_cold_ice_poly + + elemental function gsw_pt0_from_t (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_pt0_from_t + end function gsw_pt0_from_t + + elemental function gsw_pt0_from_t_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_pt0_from_t_ice + end function gsw_pt0_from_t_ice + + elemental subroutine gsw_pt_first_derivatives (sa, ct, pt_sa, pt_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8), intent(out), optional :: pt_sa, pt_ct + end subroutine gsw_pt_first_derivatives + + elemental function gsw_pt_from_ct (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_pt_from_ct + end function gsw_pt_from_ct + + elemental function gsw_pt_from_entropy (sa, entropy) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, entropy + real (r8) :: gsw_pt_from_entropy + end function gsw_pt_from_entropy + + elemental function gsw_pt_from_pot_enthalpy_ice (pot_enthalpy_ice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pot_enthalpy_ice + real (r8) :: gsw_pt_from_pot_enthalpy_ice + end function gsw_pt_from_pot_enthalpy_ice + + elemental function gsw_pt_from_pot_enthalpy_ice_poly_dh (pot_enthalpy_ice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pot_enthalpy_ice + real (r8) :: gsw_pt_from_pot_enthalpy_ice_poly_dh + end function gsw_pt_from_pot_enthalpy_ice_poly_dh + + elemental function gsw_pt_from_pot_enthalpy_ice_poly (pot_enthalpy_ice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pot_enthalpy_ice + real (r8) :: gsw_pt_from_pot_enthalpy_ice_poly + end function gsw_pt_from_pot_enthalpy_ice_poly + + elemental function gsw_pt_from_t (sa, t, p, p_ref) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p, p_ref + real (r8) :: gsw_pt_from_t + end function gsw_pt_from_t + + elemental function gsw_pt_from_t_ice (t, p, p_ref) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p, p_ref + real (r8) :: gsw_pt_from_t_ice + end function gsw_pt_from_t_ice + + elemental subroutine gsw_pt_second_derivatives (sa, ct, pt_sa_sa, & + pt_sa_ct, pt_ct_ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8), intent(out), optional :: pt_sa_sa, pt_sa_ct, pt_ct_ct + end subroutine gsw_pt_second_derivatives + + elemental subroutine gsw_rho_alpha_beta (sa, ct, p, rho, alpha, beta) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: rho, alpha, beta + end subroutine gsw_rho_alpha_beta + + elemental function gsw_rho (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_rho + end function gsw_rho + + elemental subroutine gsw_rho_first_derivatives (sa, ct, p, drho_dsa, & + drho_dct, drho_dp) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: drho_dsa, drho_dct, drho_dp + end subroutine gsw_rho_first_derivatives + + elemental subroutine gsw_rho_first_derivatives_wrt_enthalpy (sa, ct, p, & + rho_sa, rho_h) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: rho_sa, rho_h + end subroutine gsw_rho_first_derivatives_wrt_enthalpy + + elemental function gsw_rho_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_rho_ice + end function gsw_rho_ice + + elemental subroutine gsw_rho_second_derivatives (sa, ct, p, rho_sa_sa, & + rho_sa_ct, rho_ct_ct, rho_sa_p, rho_ct_p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: rho_sa_sa, rho_sa_ct, rho_ct_ct + real (r8), intent(out), optional :: rho_sa_p, rho_ct_p + end subroutine gsw_rho_second_derivatives + + elemental subroutine gsw_rho_second_derivatives_wrt_enthalpy (sa, ct, p, & + rho_sa_sa, rho_sa_h, rho_h_h) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: rho_sa_sa, rho_sa_h, rho_h_h + end subroutine gsw_rho_second_derivatives_wrt_enthalpy + + elemental function gsw_rho_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_rho_t_exact + end function gsw_rho_t_exact + + pure subroutine gsw_rr68_interp_sa_ct (sa, ct, p, p_i, sa_i, ct_i) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa(:), ct(:), p(:), p_i(:) + real (r8), intent(out) :: sa_i(:), ct_i(:) + end subroutine gsw_rr68_interp_sa_ct + + elemental function gsw_saar (p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: p, long, lat + real (r8) :: gsw_saar + end function gsw_saar + + elemental function gsw_sa_freezing_estimate (p, saturation_fraction, ct, t) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: p, saturation_fraction + real (r8), intent(in), optional :: ct, t + real (r8) :: gsw_sa_freezing_estimate + end function gsw_sa_freezing_estimate + + elemental function gsw_sa_freezing_from_ct (ct, p, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: ct, p, saturation_fraction + real (r8) :: gsw_sa_freezing_from_ct + end function gsw_sa_freezing_from_ct + + elemental function gsw_sa_freezing_from_ct_poly (ct, p, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: ct, p, saturation_fraction + real (r8) :: gsw_sa_freezing_from_ct_poly + end function gsw_sa_freezing_from_ct_poly + + elemental function gsw_sa_freezing_from_t (t, p, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p, saturation_fraction + real (r8) :: gsw_sa_freezing_from_t + end function gsw_sa_freezing_from_t + + elemental function gsw_sa_freezing_from_t_poly (t, p, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p, saturation_fraction + real (r8) :: gsw_sa_freezing_from_t_poly + end function gsw_sa_freezing_from_t_poly + + elemental function gsw_sa_from_rho (rho, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: rho, ct, p + real (r8) :: gsw_sa_from_rho + end function gsw_sa_from_rho + + elemental function gsw_sa_from_sp_baltic (sp, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sp, long, lat + real (r8) :: gsw_sa_from_sp_baltic + end function gsw_sa_from_sp_baltic + + elemental function gsw_sa_from_sp (sp, p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sp, p, long, lat + real (r8) :: gsw_sa_from_sp + end function gsw_sa_from_sp + + elemental function gsw_sa_from_sstar (sstar, p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sstar, p, long, lat + real (r8) :: gsw_sa_from_sstar + end function gsw_sa_from_sstar + + elemental function gsw_sa_p_inrange (sa, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + logical :: gsw_sa_p_inrange + end function gsw_sa_p_inrange + + elemental subroutine gsw_seaice_fraction_to_freeze_seawater (sa, ct, p, & + sa_seaice, t_seaice, sa_freeze, ct_freeze, w_seaice) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p, sa_seaice, t_seaice + real (r8), intent(out) :: sa_freeze, ct_freeze, w_seaice + end subroutine gsw_seaice_fraction_to_freeze_seawater + + elemental function gsw_sigma0 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_sigma0 + end function gsw_sigma0 + + elemental function gsw_sigma1 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_sigma1 + end function gsw_sigma1 + + elemental function gsw_sigma2 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_sigma2 + end function gsw_sigma2 + + elemental function gsw_sigma3 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_sigma3 + end function gsw_sigma3 + + elemental function gsw_sigma4 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_sigma4 + end function gsw_sigma4 + + elemental function gsw_sound_speed (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_sound_speed + end function gsw_sound_speed + + elemental function gsw_sound_speed_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_sound_speed_ice + end function gsw_sound_speed_ice + + elemental function gsw_sound_speed_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_sound_speed_t_exact + end function gsw_sound_speed_t_exact + + elemental subroutine gsw_specvol_alpha_beta (sa, ct, p, specvol, alpha, & + beta) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8), intent(out), optional :: specvol, alpha, beta + end subroutine gsw_specvol_alpha_beta + + elemental function gsw_specvol_anom_standard (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_specvol_anom_standard + end function gsw_specvol_anom_standard + + elemental function gsw_specvol (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_specvol + end function gsw_specvol + + elemental subroutine gsw_specvol_first_derivatives (sa, ct, p, v_sa, v_ct, & + v_p, iflag) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + integer, intent(in), optional :: iflag + real (r8), intent(out), optional :: v_sa, v_ct, v_p + end subroutine gsw_specvol_first_derivatives + + elemental subroutine gsw_specvol_first_derivatives_wrt_enthalpy (sa, ct, & + p, v_sa, v_h, iflag) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + integer, intent(in), optional :: iflag + real (r8), intent(out), optional :: v_sa, v_h + end subroutine gsw_specvol_first_derivatives_wrt_enthalpy + + elemental function gsw_specvol_ice (t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: t, p + real (r8) :: gsw_specvol_ice + end function gsw_specvol_ice + + elemental subroutine gsw_specvol_second_derivatives (sa, ct, p, v_sa_sa, & + v_sa_ct, v_ct_ct, v_sa_p, v_ct_p, iflag) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + integer, intent(in), optional :: iflag + real (r8), intent(out), optional :: v_sa_sa, v_sa_ct, v_ct_ct, v_sa_p, v_ct_p + end subroutine gsw_specvol_second_derivatives + + elemental subroutine gsw_specvol_second_derivatives_wrt_enthalpy (sa, ct, & + p, v_sa_sa, v_sa_h, v_h_h, iflag) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + integer, intent(in), optional :: iflag + real (r8), intent(out), optional :: v_sa_sa, v_sa_h, v_h_h + end subroutine gsw_specvol_second_derivatives_wrt_enthalpy + + elemental function gsw_specvol_sso_0 (p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: p + real (r8) :: gsw_specvol_sso_0 + end function gsw_specvol_sso_0 + + elemental function gsw_specvol_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_specvol_t_exact + end function gsw_specvol_t_exact + + elemental function gsw_sp_from_c (c, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: c, t, p + real (r8) :: gsw_sp_from_c + end function gsw_sp_from_c + + elemental function gsw_sp_from_sa_baltic (sa, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, long, lat + real (r8) :: gsw_sp_from_sa_baltic + end function gsw_sp_from_sa_baltic + + elemental function gsw_sp_from_sa (sa, p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, long, lat + real (r8) :: gsw_sp_from_sa + end function gsw_sp_from_sa + + elemental function gsw_sp_from_sk (sk) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sk + real (r8) :: gsw_sp_from_sk + end function gsw_sp_from_sk + + elemental function gsw_sp_from_sr (sr) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sr + real (r8) :: gsw_sp_from_sr + end function gsw_sp_from_sr + + elemental function gsw_sp_from_sstar (sstar, p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sstar, p, long, lat + real (r8) :: gsw_sp_from_sstar + end function gsw_sp_from_sstar + + elemental function gsw_spiciness0 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_spiciness0 + end function gsw_spiciness0 + + elemental function gsw_spiciness1 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_spiciness1 + end function gsw_spiciness1 + + elemental function gsw_spiciness2 (sa, ct) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct + real (r8) :: gsw_spiciness2 + end function gsw_spiciness2 + + elemental function gsw_sr_from_sp (sp) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sp + real (r8) :: gsw_sr_from_sp + end function gsw_sr_from_sp + + elemental function gsw_sstar_from_sa (sa, p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, long, lat + real (r8) :: gsw_sstar_from_sa + end function gsw_sstar_from_sa + + elemental function gsw_sstar_from_sp (sp, p, long, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sp, p, long, lat + real (r8) :: gsw_sstar_from_sp + end function gsw_sstar_from_sp + + elemental function gsw_t_deriv_chem_potential_water_t_exact (sa, t, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, t, p + real (r8) :: gsw_t_deriv_chem_potential_water_t_exact + end function gsw_t_deriv_chem_potential_water_t_exact + + elemental function gsw_t_freezing_exact (sa, p, saturation_fraction) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + real (r8) :: gsw_t_freezing_exact + end function gsw_t_freezing_exact + + elemental function gsw_t_freezing (sa, p, saturation_fraction, poly) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + logical, intent(in), optional :: poly + real (r8) :: gsw_t_freezing + end function gsw_t_freezing + + elemental subroutine gsw_t_freezing_first_derivatives (sa, p, & + saturation_fraction, tfreezing_sa, tfreezing_p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + real (r8), intent(out), optional :: tfreezing_sa, tfreezing_p + end subroutine gsw_t_freezing_first_derivatives + + elemental subroutine gsw_t_freezing_first_derivatives_poly (sa, p, & + saturation_fraction, tfreezing_sa, tfreezing_p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p, saturation_fraction + real (r8), intent(out), optional :: tfreezing_sa, tfreezing_p + end subroutine gsw_t_freezing_first_derivatives_poly + + elemental function gsw_t_freezing_poly (sa, p, saturation_fraction, polynomial) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, p + real (r8), intent(in), optional :: saturation_fraction + logical, intent(in), optional :: polynomial + real (r8) :: gsw_t_freezing_poly + end function gsw_t_freezing_poly + + elemental function gsw_t_from_ct (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_t_from_ct + end function gsw_t_from_ct + + elemental function gsw_t_from_pt0_ice (pt0_ice, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: pt0_ice, p + real (r8) :: gsw_t_from_pt0_ice + end function gsw_t_from_pt0_ice + + elemental function gsw_thermobaric (sa, ct, p) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa, ct, p + real (r8) :: gsw_thermobaric + end function gsw_thermobaric + + pure subroutine gsw_turner_rsubrho (sa, ct, p, tu, rsubrho, p_mid) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: sa(:), ct(:), p(:) + real (r8), intent(out) :: tu(:), rsubrho(:), p_mid(:) + end subroutine gsw_turner_rsubrho + + pure subroutine gsw_util_indx (x, n, z, k) + use gsw_mod_kinds + integer, intent(in) :: n + integer, intent(out) :: k + real (r8), intent(in), dimension(n) :: x + real (r8), intent(in) :: z + end subroutine gsw_util_indx + + pure function gsw_util_interp1q_int (x, iy, x_i) result(y_i) + use gsw_mod_kinds + implicit none + integer, intent(in) :: iy(:) + real (r8), intent(in) :: x(:), x_i(:) + real (r8) :: y_i(size(x_i)) + end function gsw_util_interp1q_int + + pure function gsw_util_sort_real (rarray) result(iarray) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: rarray(:) ! Values to be sorted + integer :: iarray(size(rarray)) ! Sorted ids + end function gsw_util_sort_real + + pure function gsw_util_xinterp1 (x, y, n, x0) + use gsw_mod_kinds + implicit none + integer, intent(in) :: n + real (r8), intent(in) :: x0 + real (r8), dimension(n), intent(in) :: x, y + real (r8) :: gsw_util_xinterp1 + end function gsw_util_xinterp1 + + elemental function gsw_z_from_p (p, lat) + use gsw_mod_kinds + implicit none + real (r8), intent(in) :: p, lat + real (r8) :: gsw_z_from_p + end function gsw_z_from_p + +end interface + +end module gsw_mod_toolbox diff --git a/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 b/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 deleted file mode 120000 index cd794a1316..0000000000 --- a/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_pt_from_ct.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 b/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 new file mode 100644 index 0000000000..b856b923c8 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_pt_from_ct.f90 @@ -0,0 +1,72 @@ +!========================================================================== +elemental function gsw_pt_from_ct (sa, ct) +!========================================================================== +! +! potential temperature of seawater from conservative temperature +! +! sa : Absolute Salinity [g/kg] +! ct : Conservative Temperature [deg C] +! p : sea pressure [dbar] +! +! gsw_pt_from_ct : potential temperature with [deg C] +! reference pressure of 0 dbar +!-------------------------------------------------------------------------- + +use gsw_mod_toolbox, only : gsw_ct_from_pt, gsw_gibbs_pt0_pt0 + +use gsw_mod_teos10_constants, only : gsw_cp0, gsw_ups, gsw_t0 + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, ct + +real (r8) :: gsw_pt_from_ct + +real (r8) :: a5ct, b3ct, ct_factor, pt_num, pt_recden, ct_diff +real (r8) :: ct0, pt, pt_old, ptm, dct, dpt_dct, s1 + +real (r8), parameter :: a0 = -1.446013646344788e-2_r8 +real (r8), parameter :: a1 = -3.305308995852924e-3_r8 +real (r8), parameter :: a2 = 1.062415929128982e-4_r8 +real (r8), parameter :: a3 = 9.477566673794488e-1_r8 +real (r8), parameter :: a4 = 2.166591947736613e-3_r8 +real (r8), parameter :: a5 = 3.828842955039902e-3_r8 + +real (r8), parameter :: b0 = 1.0_r8 +real (r8), parameter :: b1 = 6.506097115635800e-4_r8 +real (r8), parameter :: b2 = 3.830289486850898e-3_r8 +real (r8), parameter :: b3 = 1.247811760368034e-6_r8 + +s1 = sa/gsw_ups + +a5ct = a5*ct +b3ct = b3*ct + +ct_factor = (a3 + a4*s1 + a5ct) +pt_num = a0 + s1*(a1 + a2*s1) + ct*ct_factor +pt_recden = 1.0_r8/(b0 + b1*s1 + ct*(b2 + b3ct)) +pt = pt_num*pt_recden + +dpt_dct = (ct_factor + a5ct - (b2 + b3ct + b3ct)*pt)*pt_recden + +! Start the 1.5 iterations through the modified Newton-Rapshon iterative, +! method, which is also known as the Newton-McDougall method. + +ct_diff = gsw_ct_from_pt(sa,pt) - ct +pt_old = pt +pt = pt_old - ct_diff*dpt_dct +ptm = 0.5_r8*(pt + pt_old) + +dpt_dct = -gsw_cp0/((ptm + gsw_t0)*gsw_gibbs_pt0_pt0(sa,ptm)) + +pt = pt_old - ct_diff*dpt_dct +ct_diff = gsw_ct_from_pt(sa,pt) - ct +pt_old = pt +gsw_pt_from_ct = pt_old - ct_diff*dpt_dct + +return +end function + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_rho.f90 b/src/equation_of_state/TEOS10/gsw_rho.f90 deleted file mode 120000 index 22eea6219a..0000000000 --- a/src/equation_of_state/TEOS10/gsw_rho.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_rho.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_rho.f90 b/src/equation_of_state/TEOS10/gsw_rho.f90 new file mode 100644 index 0000000000..3daa65746e --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_rho.f90 @@ -0,0 +1,36 @@ +!========================================================================== +elemental function gsw_rho (sa, ct, p) +!========================================================================== +! +! Calculates in-situ density from Absolute Salinity and Conservative +! Temperature, using the computationally-efficient expression for +! specific volume in terms of SA, CT and p (Roquet et al., 2014). +! +! Note that potential density with respect to reference pressure, pr, is +! obtained by calling this function with the pressure argument being pr +! (i.e. "gsw_rho(SA,CT,pr)"). +! +! SA = Absolute Salinity [ g/kg ] +! CT = Conservative Temperature (ITS-90) [ deg C ] +! p = sea pressure [ dbar ] +! ( i.e. absolute pressure - 10.1325 dbar ) +! +! rho = in-situ density [ kg/m ] +!-------------------------------------------------------------------------- + +use gsw_mod_toolbox, only : gsw_specvol + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, ct, p + +real (r8) :: gsw_rho + +gsw_rho = 1.0_r8/gsw_specvol(sa,ct,p) + +return +end function + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 deleted file mode 120000 index 3a8ba38824..0000000000 --- a/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_rho_first_derivatives.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 new file mode 100644 index 0000000000..b4ee696a1d --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_rho_first_derivatives.f90 @@ -0,0 +1,110 @@ +!========================================================================== +elemental subroutine gsw_rho_first_derivatives (sa, ct, p, drho_dsa, & + drho_dct, drho_dp) +!========================================================================== +! +! Calculates the three (3) partial derivatives of in-situ density with +! respect to Absolute Salinity, Conservative Temperature and pressure. +! Note that the pressure derivative is done with respect to pressure in +! Pa, not dbar. This function uses the computationally-efficient expression +! for specific volume in terms of SA, CT and p (Roquet et al., 2014). +! +! SA = Absolute Salinity [ g/kg ] +! CT = Conservative Temperature (ITS-90) [ deg C ] +! p = sea pressure [ dbar ] +! ( i.e. absolute pressure - 10.1325 dbar ) +! +! drho_dSA = partial derivatives of density [ kg^2/(g m^3) ] +! with respect to Absolute Salinity +! drho_dCT = partial derivatives of density [ kg/(K m^3) ] +! with respect to Conservative Temperature +! drho_dP = partial derivatives of density [ kg/(Pa m^3) ] +! with respect to pressure in Pa +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : pa2db, gsw_sfac, offset + +use gsw_mod_specvol_coefficients + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, ct, p +real (r8), intent(out), optional :: drho_dsa, drho_dct, drho_dp + +real (r8) :: rho2, v_ct, v_p, v_sa, xs, ys, z, v + +xs = sqrt(gsw_sfac*sa + offset) +ys = ct*0.025_r8 +z = p*1e-4_r8 + +v = v000 + xs*(v010 + xs*(v020 + xs*(v030 + xs*(v040 + xs*(v050 & + + v060*xs))))) + ys*(v100 + xs*(v110 + xs*(v120 + xs*(v130 + xs*(v140 & + + v150*xs)))) + ys*(v200 + xs*(v210 + xs*(v220 + xs*(v230 + v240*xs))) & + + ys*(v300 + xs*(v310 + xs*(v320 + v330*xs)) + ys*(v400 + xs*(v410 & + + v420*xs) + ys*(v500 + v510*xs + v600*ys))))) + z*(v001 + xs*(v011 & + + xs*(v021 + xs*(v031 + xs*(v041 + v051*xs)))) + ys*(v101 + xs*(v111 & + + xs*(v121 + xs*(v131 + v141*xs))) + ys*(v201 + xs*(v211 + xs*(v221 & + + v231*xs)) + ys*(v301 + xs*(v311 + v321*xs) + ys*(v401 + v411*xs & + + v501*ys)))) + z*(v002 + xs*(v012 + xs*(v022 + xs*(v032 + v042*xs))) & + + ys*(v102 + xs*(v112 + xs*(v122 + v132*xs)) + ys*(v202 + xs*(v212 & + + v222*xs) + ys*(v302 + v312*xs + v402*ys))) + z*(v003 + xs*(v013 & + + v023*xs) + ys*(v103 + v113*xs + v203*ys) + z*(v004 + v014*xs + v104*ys & + + z*(v005 + v006*z))))) + +rho2 = (1.0_r8/v)**2 + +if (present(drho_dsa)) then + + v_sa = b000 + xs*(b100 + xs*(b200 + xs*(b300 + xs*(b400 + b500*xs)))) & + + ys*(b010 + xs*(b110 + xs*(b210 + xs*(b310 + b410*xs))) & + + ys*(b020 + xs*(b120 + xs*(b220 + b320*xs)) + ys*(b030 & + + xs*(b130 + b230*xs) + ys*(b040 + b140*xs + b050*ys)))) & + + z*(b001 + xs*(b101 + xs*(b201 + xs*(b301 + b401*xs))) & + + ys*(b011 + xs*(b111 + xs*(b211 + b311*xs)) + ys*(b021 & + + xs*(b121 + b221*xs) + ys*(b031 + b131*xs + b041*ys))) & + + z*(b002 + xs*(b102 + xs*(b202 + b302*xs))+ ys*(b012 & + + xs*(b112 + b212*xs) + ys*(b022 + b122*xs + b032*ys)) & + + z*(b003 + b103*xs + b013*ys + b004*z))) + + drho_dsa = -rho2*0.5_r8*gsw_sfac*v_sa/xs + +end if + +if (present(drho_dct)) then + + v_ct = a000 + xs*(a100 + xs*(a200 + xs*(a300 + xs*(a400 + a500*xs)))) & + + ys*(a010 + xs*(a110 + xs*(a210 + xs*(a310 + a410*xs))) & + + ys*(a020 + xs*(a120 + xs*(a220 + a320*xs)) + ys*(a030 & + + xs*(a130 + a230*xs) + ys*(a040 + a140*xs + a050*ys )))) & + + z*(a001 + xs*(a101 + xs*(a201 + xs*(a301 + a401*xs))) & + + ys*(a011 + xs*(a111 + xs*(a211 + a311*xs)) + ys*(a021 & + + xs*(a121 + a221*xs) + ys*(a031 + a131*xs + a041*ys))) & + + z*(a002 + xs*(a102 + xs*(a202 + a302*xs)) + ys*(a012 & + + xs*(a112 + a212*xs) + ys*(a022 + a122*xs + a032*ys)) & + + z*(a003 + a103*xs + a013*ys + a004*z))) + + drho_dct = -rho2*0.025_r8*v_ct + +end if + +if (present(drho_dp)) then + + v_p = c000 + xs*(c100 + xs*(c200 + xs*(c300 + xs*(c400 + c500*xs)))) & + + ys*(c010 + xs*(c110 + xs*(c210 + xs*(c310 + c410*xs))) + ys*(c020 & + + xs*(c120 + xs*(c220 + c320*xs)) + ys*(c030 + xs*(c130 + c230*xs) & + + ys*(c040 + c140*xs + c050*ys)))) + z*(c001 + xs*(c101 + xs*(c201 & + + xs*(c301 + c401*xs))) + ys*(c011 + xs*(c111 + xs*(c211 + c311*xs)) & + + ys*(c021 + xs*(c121 + c221*xs) + ys*(c031 + c131*xs + c041*ys))) & + + z*(c002 + xs*(c102 + c202*xs) + ys*(c012 + c112*xs + c022*ys) & + + z*(c003 + c103*xs + c013*ys + z*(c004 + c005*z)))) + + drho_dp = -rho2*1e-4_r8*pa2db*v_p + +end if + +return +end subroutine + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 b/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 deleted file mode 120000 index d8cd41f4bf..0000000000 --- a/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_sp_from_sr.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 b/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 new file mode 100644 index 0000000000..c01377546c --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_sp_from_sr.f90 @@ -0,0 +1,30 @@ +!========================================================================== +elemental function gsw_sp_from_sr (sr) +!========================================================================== +! +! Calculates Practical Salinity, sp, from Reference Salinity, sr. +! +! sr : Reference Salinity [g/kg] +! +! gsw_sp_from_sr : Practical Salinity [unitless] +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : gsw_ups + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sr + +real (r8) :: gsw_sp_from_sr + +gsw_sp_from_sr = sr/gsw_ups + +return +end function + +!-------------------------------------------------------------------------- + + + diff --git a/src/equation_of_state/TEOS10/gsw_specvol.f90 b/src/equation_of_state/TEOS10/gsw_specvol.f90 deleted file mode 120000 index 7a41a5cea0..0000000000 --- a/src/equation_of_state/TEOS10/gsw_specvol.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_specvol.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_specvol.f90 b/src/equation_of_state/TEOS10/gsw_specvol.f90 new file mode 100644 index 0000000000..00cfaab125 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_specvol.f90 @@ -0,0 +1,52 @@ +!========================================================================== +elemental function gsw_specvol (sa, ct, p) +!========================================================================== +! +! Calculates specific volume from Absolute Salinity, Conservative +! Temperature and pressure, using the computationally-efficient +! polynomial expression for specific volume (Roquet et al., 2014). +! +! SA = Absolute Salinity [ g/kg ] +! CT = Conservative Temperature (ITS-90) [ deg C ] +! p = sea pressure [ dbar ] +! ( i.e. absolute pressure - 10.1325 dbar ) +! +! specvol = specific volume [ m^3/kg ] +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : gsw_sfac, offset + +use gsw_mod_specvol_coefficients + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, ct, p + +real (r8) :: gsw_specvol + +real (r8) :: xs, ys, z + +xs = sqrt(gsw_sfac*sa + offset) +ys = ct*0.025_r8 +z = p*1e-4_r8 + +gsw_specvol = v000 + xs*(v010 + xs*(v020 + xs*(v030 + xs*(v040 + xs*(v050 & + + v060*xs))))) + ys*(v100 + xs*(v110 + xs*(v120 + xs*(v130 + xs*(v140 & + + v150*xs)))) + ys*(v200 + xs*(v210 + xs*(v220 + xs*(v230 + v240*xs))) & + + ys*(v300 + xs*(v310 + xs*(v320 + v330*xs)) + ys*(v400 + xs*(v410 & + + v420*xs) + ys*(v500 + v510*xs + v600*ys))))) + z*(v001 + xs*(v011 & + + xs*(v021 + xs*(v031 + xs*(v041 + v051*xs)))) + ys*(v101 + xs*(v111 & + + xs*(v121 + xs*(v131 + v141*xs))) + ys*(v201 + xs*(v211 + xs*(v221 & + + v231*xs)) + ys*(v301 + xs*(v311 + v321*xs) + ys*(v401 + v411*xs & + + v501*ys)))) + z*(v002 + xs*(v012 + xs*(v022 + xs*(v032 + v042*xs))) & + + ys*(v102 + xs*(v112 + xs*(v122 + v132*xs)) + ys*(v202 + xs*(v212 & + + v222*xs) + ys*(v302 + v312*xs + v402*ys))) + z*(v003 + xs*(v013 & + + v023*xs) + ys*(v103 + v113*xs + v203*ys) + z*(v004 + v014*xs + v104*ys & + + z*(v005 + v006*z))))) + +return +end function + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 deleted file mode 120000 index ee6ee1f906..0000000000 --- a/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_specvol_first_derivatives.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 new file mode 100644 index 0000000000..2f2a006b17 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_specvol_first_derivatives.f90 @@ -0,0 +1,104 @@ +!========================================================================== +elemental subroutine gsw_specvol_first_derivatives (sa, ct, p, v_sa, v_ct, & + v_p, iflag) +! ========================================================================= +! +! Calculates three first-order derivatives of specific volume (v). +! Note that this function uses the computationally-efficient +! expression for specific volume (Roquet et al., 2014). +! +! SA = Absolute Salinity [ g/kg ] +! CT = Conservative Temperature (ITS-90) [ deg C ] +! p = sea pressure [ dbar ] +! ( i.e. absolute pressure - 10.1325 dbar ) +! +! v_SA = The first derivative of specific volume with respect to +! Absolute Salinity at constant CT & p. [ J/(kg (g/kg)^2) ] +! v_CT = The first derivative of specific volume with respect to +! CT at constant SA and p. [ J/(kg K(g/kg)) ] +! v_P = The first derivative of specific volume with respect to +! P at constant SA and CT. [ J/(kg K^2) ] +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : gsw_sfac, offset + +use gsw_mod_specvol_coefficients + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, ct, p +integer, intent(in), optional :: iflag +real (r8), intent(out), optional :: v_sa, v_ct, v_p + +integer :: i +logical :: flags(3) +real (r8) :: v_ct_part, v_p_part, v_sa_part, xs, ys, z + +xs = sqrt(gsw_sfac*sa + offset) +ys = ct*0.025_r8 +z = p*1e-4_r8 + +if (present(iflag)) then + do i = 1, 3 + flags(i) = btest(iflag,i) + end do +else + flags = .true. +end if + +if (present(v_sa) .and. flags(1)) then + + v_sa_part = b000 + xs*(b100 + xs*(b200 + xs*(b300 + xs*(b400 + b500*xs)))) & + + ys*(b010 + xs*(b110 + xs*(b210 + xs*(b310 + b410*xs))) & + + ys*(b020 + xs*(b120 + xs*(b220 + b320*xs)) + ys*(b030 & + + xs*(b130 + b230*xs) + ys*(b040 + b140*xs + b050*ys)))) & + + z*(b001 + xs*(b101 + xs*(b201 + xs*(b301 + b401*xs))) & + + ys*(b011 + xs*(b111 + xs*(b211 + b311*xs)) + ys*(b021 & + + xs*(b121 + b221*xs) + ys*(b031 + b131*xs + b041*ys))) & + + z*(b002 + xs*(b102 + xs*(b202 + b302*xs))+ ys*(b012 & + + xs*(b112 + b212*xs) + ys*(b022 + b122*xs + b032*ys)) & + + z*(b003 + b103*xs + b013*ys + b004*z))) + + v_sa = 0.5_r8*gsw_sfac*v_sa_part/xs + +end if + + +if (present(v_ct) .and. flags(2)) then + + v_ct_part = a000 + xs*(a100 + xs*(a200 + xs*(a300 + xs*(a400 + a500*xs)))) & + + ys*(a010 + xs*(a110 + xs*(a210 + xs*(a310 + a410*xs))) & + + ys*(a020 + xs*(a120 + xs*(a220 + a320*xs)) + ys*(a030 & + + xs*(a130 + a230*xs) + ys*(a040 + a140*xs + a050*ys )))) & + + z*(a001 + xs*(a101 + xs*(a201 + xs*(a301 + a401*xs))) & + + ys*(a011 + xs*(a111 + xs*(a211 + a311*xs)) + ys*(a021 & + + xs*(a121 + a221*xs) + ys*(a031 + a131*xs + a041*ys))) & + + z*(a002 + xs*(a102 + xs*(a202 + a302*xs)) + ys*(a012 & + + xs*(a112 + a212*xs) + ys*(a022 + a122*xs + a032*ys)) & + + z*(a003 + a103*xs + a013*ys + a004*z))) + + v_ct = 0.025_r8*v_ct_part + +end if + +if (present(v_p) .and. flags(3)) then + + v_p_part = c000 + xs*(c100 + xs*(c200 + xs*(c300 + xs*(c400 + c500*xs)))) & + + ys*(c010 + xs*(c110 + xs*(c210 + xs*(c310 + c410*xs))) + ys*(c020 & + + xs*(c120 + xs*(c220 + c320*xs)) + ys*(c030 + xs*(c130 + c230*xs) & + + ys*(c040 + c140*xs + c050*ys)))) + z*(c001 + xs*(c101 + xs*(c201 & + + xs*(c301 + c401*xs))) + ys*(c011 + xs*(c111 + xs*(c211 + c311*xs)) & + + ys*(c021 + xs*(c121 + c221*xs) + ys*(c031 + c131*xs + c041*ys))) & + + z*( c002 + xs*(c102 + c202*xs) + ys*(c012 + c112*xs + c022*ys) & + + z*(c003 + c103*xs + c013*ys + z*(c004 + c005*z)))) + + v_p = 1e-8_r8*v_p_part + +end if + +return +end subroutine + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 b/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 deleted file mode 120000 index eda229ff66..0000000000 --- a/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 +++ /dev/null @@ -1 +0,0 @@ -../../../pkg/GSW-Fortran/toolbox/gsw_sr_from_sp.f90 \ No newline at end of file diff --git a/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 b/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 new file mode 100644 index 0000000000..cbcc4fea0b --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_sr_from_sp.f90 @@ -0,0 +1,30 @@ +!========================================================================== +elemental function gsw_sr_from_sp (sp) +!========================================================================== +! +! Calculates Reference Salinity, SR, from Practical Salinity, SP. +! +! sp : Practical Salinity [unitless] +! +! gsw_sr_from_sp : Reference Salinity [g/kg] +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : gsw_ups + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sp + +real (r8) :: gsw_sr_from_sp + +gsw_sr_from_sp = sp*gsw_ups + +return +end function + +!-------------------------------------------------------------------------- + + + From eff7bf2a62ace2c4924089952bf275454d0f1d0b Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 16 Dec 2016 15:28:14 -0500 Subject: [PATCH 6/9] Changed all stringd "reference salinity" to "absolute salinity" - To avoid further confusions I changed all references to "reference salinity" to "absolute salinity". --- config_src/coupled_driver/ocean_model_MOM.F90 | 10 +++---- src/core/MOM.F90 | 28 +++++++++---------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index c73d7e52ed..00146a11b2 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -488,7 +488,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! Translate state into Ocean. ! call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & ! Ice_ocean_boundary%p, OS%press_to_z) - call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, OS%MOM_CSp%use_conT_refS) + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, OS%MOM_CSp%use_conT_absS) call callTree_leave("update_ocean_model()") end subroutine update_ocean_model @@ -731,12 +731,12 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, maskmap) end subroutine initialize_ocean_public_type -subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_refS, patm, press_to_z) +subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_absS, patm, press_to_z) use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct type(surface), intent(inout) :: state type(ocean_public_type), target, intent(inout) :: Ocean_sfc type(ocean_grid_type), intent(inout) :: G - logical, intent(in) :: use_conT_refS + logical, intent(in) :: use_conT_absS real, optional, intent(in) :: patm(:,:) real, optional, intent(in) :: press_to_z ! This subroutine translates the coupler's ocean_data_type into MOM's @@ -765,7 +765,7 @@ subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_refS, patm, !from conservative T to potential T and !from absolute (reference) salinity to practical salinity ! - if(use_conT_refS) then + if(use_conT_absS) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd Ocean_sfc%s_surf(i,j) = gsw_sp_from_sr(state%SSS(i+i0,j+j0)) Ocean_sfc%t_surf(i,j) = gsw_pt_from_ct(state%SSS(i+i0,j+j0),state%SST(i+i0,j+j0)) + CELSIUS_KELVIN_OFFSET @@ -831,7 +831,7 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc) OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%ave_ssh,& OS%grid, OS%GV, OS%MOM_CSp) - call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, OS%MOM_CSp%use_conT_refS) + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, OS%MOM_CSp%use_conT_absS) end subroutine ocean_model_init_sfc ! diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index b6ccc1abfb..bb6861d02e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -190,8 +190,8 @@ module MOM !! with nkml sublayers and nkbl buffer layer. logical :: diabatic_first !< If true, apply diabatic and thermodynamic !! processes before time stepping the dynamics. - logical :: use_conT_refS !< If true, , the prognostics T&S are the conservative temperature - !! and reference salinity. Care should be taken to convert them + logical :: use_conT_absS !< If true, , the prognostics T&S are the conservative temperature + !! and absolute salinity. Care should be taken to convert them !! to potential temperature and practical salinity before !! exchanging them with the coupler and/or reporting T&S diagnostics. logical :: thickness_diffuse !< If true, diffuse interface height w/ a diffusivity KHTH. @@ -288,7 +288,7 @@ module MOM integer :: id_T = -1 integer :: id_S = -1 integer :: id_Tcon = -1 - integer :: id_Sref = -1 + integer :: id_Sabs = -1 ! 2-d surface and bottom fields integer :: id_zos = -1 @@ -307,7 +307,7 @@ module MOM integer :: id_tob = -1 integer :: id_sob = -1 integer :: id_consst = -1 - integer :: id_refsss = -1 + integer :: id_abssss = -1 ! heat and salt flux fields integer :: id_fraz = -1 @@ -1191,14 +1191,14 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! post some diagnostics ! - !In case the internal representation of the model is conservative temperature and reference salinity + !In case the internal representation of the model is conservative temperature and absolute salinity !convert, for diagnostics !from conservative temp to potential temp and - !from reference salinity to practical salinity + !from absolute salinity to practical salinity ! - if(CS%use_conT_refS) then + if(CS%use_conT_absS) then if (CS%id_Tcon > 0) call post_data(CS%id_Tcon, CS%tv%T, CS%diag) - if (CS%id_Sref > 0) call post_data(CS%id_Sref, CS%tv%S, CS%diag) + if (CS%id_Sabs > 0) call post_data(CS%id_Sabs, CS%tv%S, CS%diag) !Conversions do k=1,nz ; do j=js,je ; do i=is,ie pracSal(i,j,k) = gsw_sp_from_sr(CS%tv%S(i,j,k)) @@ -2087,9 +2087,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo "If true, apply diabatic and thermodynamic processes, \n"//& "including buoyancy forcing and mass gain or loss, \n"//& "before stepping the dynamics forward.", default=.false.) - call get_param(param_file, "MOM", "USE_CONTEMP_REFSAL", CS%use_conT_refS, & + call get_param(param_file, "MOM", "USE_CONTEMP_ABSSAL", CS%use_conT_absS, & "If true, , the prognostics T&S are the conservative temperature \n"//& - "and reference salinity. Care should be taken to convert them \n"//& + "and absolute salinity. Care should be taken to convert them \n"//& "to potential temperature and practical salinity before \n"//& "exchanging them with the coupler and/or reporting T&S diagnostics. \n"& , default=.false.) @@ -2944,12 +2944,12 @@ subroutine register_diags(Time, G, GV, CS, ADp) if (CS%id_sss_sq > 0) call safe_alloc_ptr(CS%SSS_sq,isd,ied,jsd,jed) CS%id_Tcon = register_diag_field('ocean_model', 'contemp', diag%axesTL, Time, & 'Conservative Temperature', 'Celsius') - CS%id_Sref = register_diag_field('ocean_model', 'refsalt', diag%axesTL, Time, & - long_name='Reference Salinity', units='g/Kg') + CS%id_Sabs = register_diag_field('ocean_model', 'abssalt', diag%axesTL, Time, & + long_name='Absolute Salinity', units='g/Kg') CS%id_consst = register_diag_field('ocean_model', 'conSST', diag%axesT1, Time, & 'Sea Surface Conservative Temperature', 'Celsius', CS%missing) - CS%id_refsss = register_diag_field('ocean_model', 'refSSS', diag%axesT1, Time, & - 'Sea Surface Reference Salinity', 'g/Kg', CS%missing) + CS%id_abssss = register_diag_field('ocean_model', 'absSSS', diag%axesT1, Time, & + 'Sea Surface Absolute Salinity', 'g/Kg', CS%missing) endif From 90958311ff4aa9d84424a6f22162e030d3612029 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 16 Dec 2016 17:55:14 -0500 Subject: [PATCH 7/9] Fixed surface T&S diagnostics and use statements - Surface T&S diagnostics need to be converted to potTemp&pracSal - Added T&S diagnostics for conTemp and absSal - published TEOS10 functions used in MOM.F90 from MOM_EOS.F90 - Added ,only to use statements in MOM_EOS.F90 --- src/core/MOM.F90 | 35 ++++++++++++++-------- src/equation_of_state/MOM_EOS.F90 | 38 +++++++++++------------- src/equation_of_state/MOM_EOS_TEOS10.F90 | 4 ++- 3 files changed, 43 insertions(+), 34 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index bb6861d02e..7e23b6198f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -88,6 +88,7 @@ module MOM use MOM_dynamics_legacy_split, only : adjustments_dyn_legacy_split, MOM_dyn_legacy_split_CS use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid use MOM_EOS, only : EOS_init +use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_error_checking, only : check_redundant use MOM_grid, only : ocean_grid_type, set_first_direction use MOM_grid, only : MOM_grid_init, MOM_grid_end @@ -124,7 +125,6 @@ module MOM use MOM_vert_friction, only : vertvisc_limit_vel, vertvisc_init use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd use MOM_verticalGrid, only : get_thickness_units, get_flux_units, get_tr_flux_units -use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct ! Offline modules use MOM_offline_transport, only : offline_transport_CS use MOM_offline_transport, only : transport_by_files, next_modulo_time @@ -306,8 +306,8 @@ module MOM integer :: id_ssh_inst = -1 integer :: id_tob = -1 integer :: id_sob = -1 - integer :: id_consst = -1 - integer :: id_abssss = -1 + integer :: id_sstcon = -1 + integer :: id_sssabs = -1 ! heat and salt flux fields integer :: id_fraz = -1 @@ -493,7 +493,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) h ! h : layer thickness (meter (Bouss) or kg/m2 (non-Bouss)) real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_predia, eta_preale - real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: potTemp, pracSal + real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: potTemp, pracSal !TEOS10 Diagnostics real, dimension(SZIB_(CS%G), SZJ_(CS%G)) :: umo2d ! Diagnostics real, dimension(SZI_(CS%G), SZJB_(CS%G)) :: vmo2d ! Diagnostics @@ -1430,24 +1430,35 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call enable_averaging(dt*n_max,Time_local, CS%diag) - if (CS%id_sst > 0) & - call post_data(CS%id_sst, state%SST, CS%diag, mask=G%mask2dT) + !We may need to convert surface T&S diagnostics from conservative&absolute to potential&practical + if(.NOT. CS%use_conT_absS) then + if (CS%id_sst > 0) call post_data(CS%id_sst, state%SST, CS%diag, mask=G%mask2dT) + if (CS%id_sss > 0) call post_data(CS%id_sss, state%SSS, CS%diag, mask=G%mask2dT) + else + if (CS%id_sstcon > 0) call post_data(CS%id_sstcon, state%SST, CS%diag, mask=G%mask2dT) + if (CS%id_sssabs > 0) call post_data(CS%id_sssabs, state%SSS, CS%diag, mask=G%mask2dT) + !Conversions + do j=js,je ; do i=is,ie + pracSal(i,j,1) = gsw_sp_from_sr(state%SSS(i,j)) + potTemp(i,j,1) = gsw_pt_from_ct(state%SSS(i,j),state%SST(i,j)) + enddo ; enddo + if (CS%id_sst > 0) call post_data(CS%id_sst, potTemp(:,:,1), CS%diag, mask=G%mask2dT) + if (CS%id_sss > 0) call post_data(CS%id_sss, pracSal(:,:,1), CS%diag, mask=G%mask2dT) + endif + if (CS%id_sst_sq > 0) then do j=js,je ; do i=is,ie CS%SST_sq(i,j) = state%SST(i,j)*state%SST(i,j) enddo ; enddo call post_data(CS%id_sst_sq, CS%SST_sq, CS%diag, mask=G%mask2dT) endif - - if (CS%id_sss > 0) & - call post_data(CS%id_sss, state%SSS, CS%diag, mask=G%mask2dT) if (CS%id_sss_sq > 0) then do j=js,je ; do i=is,ie CS%SSS_sq(i,j) = state%SSS(i,j)*state%SSS(i,j) enddo ; enddo call post_data(CS%id_sss_sq, CS%SSS_sq, CS%diag, mask=G%mask2dT) endif - + if (CS%id_ssu > 0) & call post_data(CS%id_ssu, state%u, CS%diag, mask=G%mask2dCu) if (CS%id_ssv > 0) & @@ -2946,9 +2957,9 @@ subroutine register_diags(Time, G, GV, CS, ADp) 'Conservative Temperature', 'Celsius') CS%id_Sabs = register_diag_field('ocean_model', 'abssalt', diag%axesTL, Time, & long_name='Absolute Salinity', units='g/Kg') - CS%id_consst = register_diag_field('ocean_model', 'conSST', diag%axesT1, Time, & + CS%id_sstcon = register_diag_field('ocean_model', 'conSST', diag%axesT1, Time, & 'Sea Surface Conservative Temperature', 'Celsius', CS%missing) - CS%id_abssss = register_diag_field('ocean_model', 'absSSS', diag%axesT1, Time, & + CS%id_sssabs = register_diag_field('ocean_model', 'absSSS', diag%axesT1, Time, & 'Sea Surface Absolute Salinity', 'g/Kg', CS%missing) endif diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index feb95dd488..54478a4dea 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -3,11 +3,22 @@ module MOM_EOS ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_EOS_linear -use MOM_EOS_Wright -use MOM_EOS_UNESCO -use MOM_EOS_TEOS10 -use MOM_EOS_NEMO +use MOM_EOS_linear, only : calculate_density_scalar_linear, calculate_density_array_linear +use MOM_EOS_linear, only : calculate_density_derivs_linear, calculate_specvol_derivs_linear, int_density_dz_linear +use MOM_EOS_linear, only : calculate_compress_linear, calculate_2_densities_linear, int_spec_vol_dp_linear +use MOM_EOS_Wright, only : calculate_density_scalar_wright, calculate_density_array_wright +use MOM_EOS_Wright, only : calculate_density_derivs_wright, calculate_specvol_derivs_wright, int_density_dz_wright +use MOM_EOS_Wright, only : calculate_compress_wright, calculate_2_densities_wright, int_spec_vol_dp_wright +use MOM_EOS_UNESCO, only : calculate_density_scalar_unesco, calculate_density_array_unesco +use MOM_EOS_UNESCO, only : calculate_density_derivs_unesco, calculate_density_unesco +use MOM_EOS_UNESCO, only : calculate_compress_unesco, calculate_2_densities_unesco +use MOM_EOS_NEMO, only : calculate_density_scalar_nemo, calculate_density_array_nemo +use MOM_EOS_NEMO, only : calculate_density_derivs_nemo, calculate_density_nemo +use MOM_EOS_NEMO, only : calculate_compress_nemo, calculate_2_densities_nemo +use MOM_EOS_TEOS10, only : calculate_density_scalar_teos10, calculate_density_array_teos10 +use MOM_EOS_TEOS10, only : calculate_density_derivs_teos10, calculate_specvol_derivs_teos10 +use MOM_EOS_TEOS10, only : calculate_compress_teos10, calculate_2_densities_teos10 +use MOM_EOS_TEOS10, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_TFreeze, only : calculate_TFreeze_linear, calculate_TFreeze_Millero use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_file_parser, only : get_param, log_version, param_file_type @@ -29,6 +40,7 @@ module MOM_EOS public find_depth_of_pressure_in_cell public calculate_TFreeze public convert_temp_salt_for_TEOS10 +public gsw_sp_from_sr, gsw_pt_from_ct !> Calculates density of sea water from T, S and P interface calculate_density @@ -203,22 +215,6 @@ subroutine calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npt integer, intent(in) :: npts !< The number of values to calculate type(EOS_type), pointer :: EOS !< Equation of state structure !! - !!Testing section - !!NEMO Check value: rho = 1027.45140 kg/m^3 for p=1000 dbar, ct=10 Celcius, sa=30 g/kg - !real :: T0(1),S0(1),pressure0(1) - !T0(1) = 10. - !S0(1) = 30. - !pressure0(1) = 1000. * 1e4 !pa - !call calculate_density_derivs_nemo(T0, S0, pressure0, drho_dT, drho_dS, 1,1) - !print*, 'NEMO drho: ', drho_dT(1), drho_dS(1) - !call calculate_density_derivs_teos10(T0, S0, pressure0, drho_dT, drho_dS, 1,1) - !print*, 'TEOS10 drho: ', drho_dT(1), drho_dS(1) - !call calculate_density_derivs_wright(T0, S0, pressure0, drho_dT, drho_dS, 1,1) - !print*, 'WRIGHT drho: ', drho_dT(1), drho_dS(1) - !call calculate_density_derivs_unesco(T0, S0, pressure0, drho_dT, drho_dS, 1,1) - !print*, 'UNESCO drho: ', drho_dT(1), drho_dS(1) - !stop - !! if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 index c45e41e917..38bb0a0df8 100644 --- a/src/equation_of_state/MOM_EOS_TEOS10.F90 +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -25,8 +25,9 @@ module MOM_EOS_TEOS10 !* sea water using the TEOS10 functions * !*********************************************************************** -!use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt +use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct use gsw_mod_toolbox, only : gsw_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives +!use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt implicit none ; private @@ -34,6 +35,7 @@ module MOM_EOS_TEOS10 public calculate_density_derivs_teos10, calculate_2_densities_teos10 public calculate_specvol_derivs_teos10 public calculate_density_scalar_teos10, calculate_density_array_teos10 +public gsw_sp_from_sr, gsw_pt_from_ct interface calculate_density_teos10 module procedure calculate_density_scalar_teos10, calculate_density_array_teos10 From 8274beb41ace521d032367e53bbfc44c14b6e2b3 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 16 Dec 2016 18:33:00 -0500 Subject: [PATCH 8/9] Corrections and comments - use TEOS-10 functions published from MOM_EOS rather than directly - Added comments and made diffs closer to the original --- config_src/coupled_driver/ocean_model_MOM.F90 | 2 +- src/core/MOM.F90 | 33 ++++++++++--------- 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 00146a11b2..09827882ef 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -71,6 +71,7 @@ module ocean_model_mod use fms_mod, only : stdout use mpp_mod, only : mpp_chksum use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE +use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct #include @@ -732,7 +733,6 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, maskmap) end subroutine initialize_ocean_public_type subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_absS, patm, press_to_z) - use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct type(surface), intent(inout) :: state type(ocean_public_type), target, intent(inout) :: Ocean_sfc type(ocean_grid_type), intent(inout) :: G diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7e23b6198f..cbd590c971 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1190,16 +1190,20 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! post some diagnostics - ! - !In case the internal representation of the model is conservative temperature and absolute salinity - !convert, for diagnostics - !from conservative temp to potential temp and - !from absolute salinity to practical salinity - ! - if(CS%use_conT_absS) then + if(.NOT. CS%use_conT_absS) then + !Internal T&S variables are assumed to be potential&practical + if (CS%id_T > 0) call post_data(CS%id_T, CS%tv%T, CS%diag) + if (CS%id_S > 0) call post_data(CS%id_S, CS%tv%S, CS%diag) + + if (CS%id_tob > 0) call post_data(CS%id_tob, CS%tv%T(:,:,G%ke), CS%diag, mask=G%mask2dT) + if (CS%id_sob > 0) call post_data(CS%id_sob, CS%tv%S(:,:,G%ke), CS%diag, mask=G%mask2dT) + else + !Internal T&S variables are assumed to be conservative&absolute if (CS%id_Tcon > 0) call post_data(CS%id_Tcon, CS%tv%T, CS%diag) if (CS%id_Sabs > 0) call post_data(CS%id_Sabs, CS%tv%S, CS%diag) - !Conversions + !Using TEOS-10 function calls convert T&S diagnostics + !from conservative temp to potential temp and + !from absolute salinity to practical salinity do k=1,nz ; do j=js,je ; do i=is,ie pracSal(i,j,k) = gsw_sp_from_sr(CS%tv%S(i,j,k)) potTemp(i,j,k) = gsw_pt_from_ct(CS%tv%S(i,j,k),CS%tv%T(i,j,k)) @@ -1208,12 +1212,6 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (CS%id_S > 0) call post_data(CS%id_S, pracSal, CS%diag) if (CS%id_tob > 0) call post_data(CS%id_tob, potTemp(:,:,G%ke), CS%diag, mask=G%mask2dT) if (CS%id_sob > 0) call post_data(CS%id_sob, pracSal(:,:,G%ke), CS%diag, mask=G%mask2dT) - else - if (CS%id_T > 0) call post_data(CS%id_T, CS%tv%T, CS%diag) - if (CS%id_S > 0) call post_data(CS%id_S, CS%tv%S, CS%diag) - - if (CS%id_tob > 0) call post_data(CS%id_tob, CS%tv%T(:,:,G%ke), CS%diag, mask=G%mask2dT) - if (CS%id_sob > 0) call post_data(CS%id_sob, CS%tv%S(:,:,G%ke), CS%diag, mask=G%mask2dT) endif if (CS%id_Tadx > 0) call post_data(CS%id_Tadx, CS%T_adx, CS%diag) @@ -1430,14 +1428,17 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) call enable_averaging(dt*n_max,Time_local, CS%diag) - !We may need to convert surface T&S diagnostics from conservative&absolute to potential&practical if(.NOT. CS%use_conT_absS) then + !Internal T&S variables are assumed to be potential&practical if (CS%id_sst > 0) call post_data(CS%id_sst, state%SST, CS%diag, mask=G%mask2dT) if (CS%id_sss > 0) call post_data(CS%id_sss, state%SSS, CS%diag, mask=G%mask2dT) else + !Internal T&S variables are assumed to be conservative&absolute if (CS%id_sstcon > 0) call post_data(CS%id_sstcon, state%SST, CS%diag, mask=G%mask2dT) if (CS%id_sssabs > 0) call post_data(CS%id_sssabs, state%SSS, CS%diag, mask=G%mask2dT) - !Conversions + !Using TEOS-10 function calls convert T&S diagnostics + !from conservative temp to potential temp and + !from absolute salinity to practical salinity do j=js,je ; do i=is,ie pracSal(i,j,1) = gsw_sp_from_sr(state%SSS(i,j)) potTemp(i,j,1) = gsw_pt_from_ct(state%SSS(i,j),state%SST(i,j)) From 1fee3ff9b23f5a6d3944c61a557512c5e1014190 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Tue, 3 Jan 2017 17:48:14 -0500 Subject: [PATCH 9/9] Removing unnecessary if statement - Input T&S are masked so it is unnecessary to check - Note that this changes answers so it means that input S does become slightly negative at some places! - But bypassing negative S via if statement is not safe. It causes restart issues when the model uses a processor mask_table. - Slighlty negative S is tolerated by NEMO formulation since S appears in sqrt as sqrt(S + epsilon) --- src/equation_of_state/MOM_EOS_NEMO.F90 | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/equation_of_state/MOM_EOS_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 index 6edbe3d352..ff7fbd8177 100644 --- a/src/equation_of_state/MOM_EOS_NEMO.F90 +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -233,7 +233,6 @@ subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts) 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 - if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? !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 @@ -292,8 +291,6 @@ subroutine calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start 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 - if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? - !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) @@ -373,7 +370,6 @@ subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) 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 - if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j)) enddo end subroutine calculate_compress_nemo @@ -403,7 +399,6 @@ subroutine calculate_2_densities_nemo( T, S, pressure1, pressure2, rho1, rho2, s 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 - if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value? !The following algorithm was provided by Roquet in a private communication. !It is not necessarily the algorithm used in NEMO ocean! zp1 = pressure1 * r1_P0 ! pressure (first converted to decibar)