diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 0932758432..bd5965907c 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -28,9 +28,13 @@ module MOM_EOS use MOM_EOS_UNESCO, only : calculate_density_second_derivs_UNESCO use MOM_EOS_UNESCO, only : calculate_compress_unesco use MOM_EOS_NEMO, only : calculate_density_nemo -use MOM_EOS_NEMO, only : calculate_density_derivs_nemo, calculate_density_nemo +use MOM_EOS_NEMO, only : calculate_density_derivs_nemo use MOM_EOS_NEMO, only : calculate_density_second_derivs_NEMO use MOM_EOS_NEMO, only : calculate_compress_nemo +use MOM_EOS_Roquet_SpV, only : calculate_density_Roquet_SpV, calculate_spec_vol_Roquet_SpV +use MOM_EOS_Roquet_SpV, only : calculate_density_derivs_Roquet_SpV, calculate_specvol_derivs_Roquet_SpV +use MOM_EOS_Roquet_SpV, only : calculate_compress_Roquet_SpV +use MOM_EOS_Roquet_SpV, only : calculate_density_second_derivs_Roquet_SpV use MOM_EOS_TEOS10, only : calculate_density_teos10, calculate_spec_vol_teos10 use MOM_EOS_TEOS10, only : calculate_density_derivs_teos10 use MOM_EOS_TEOS10, only : calculate_specvol_derivs_teos10 @@ -169,14 +173,18 @@ module MOM_EOS integer, parameter, public :: EOS_WRIGHT_RED = 5 !< A named integer specifying an equation of state integer, parameter, public :: EOS_TEOS10 = 6 !< A named integer specifying an equation of state integer, parameter, public :: EOS_NEMO = 7 !< A named integer specifying an equation of state +integer, parameter, public :: EOS_ROQUET_SPV = 8 !< A named integer specifying an equation of state -character*(10), parameter :: EOS_LINEAR_STRING = "LINEAR" !< A string for specifying the equation of state -character*(10), parameter :: EOS_UNESCO_STRING = "UNESCO" !< A string for specifying the equation of state -character*(10), parameter :: EOS_WRIGHT_STRING = "WRIGHT" !< A string for specifying the equation of state +character*(12), parameter :: EOS_LINEAR_STRING = "LINEAR" !< A string for specifying the equation of state +character*(12), parameter :: EOS_UNESCO_STRING = "UNESCO" !< A string for specifying the equation of state +character*(12), parameter :: EOS_JACKETT_STRING = "JACKETT_MCD" !< A string for specifying the equation of state +character*(12), parameter :: EOS_WRIGHT_STRING = "WRIGHT" !< A string for specifying the equation of state character*(12), parameter :: EOS_WRIGHT_RED_STRING = "WRIGHT_RED" !< A string for specifying the equation of state character*(12), parameter :: EOS_WRIGHT_FULL_STRING = "WRIGHT_FULL" !< A string for specifying the equation of state -character*(10), parameter :: EOS_TEOS10_STRING = "TEOS10" !< A string for specifying the equation of state -character*(10), parameter :: EOS_NEMO_STRING = "NEMO" !< A string for specifying the equation of state +character*(12), parameter :: EOS_TEOS10_STRING = "TEOS10" !< A string for specifying the equation of state +character*(12), parameter :: EOS_NEMO_STRING = "NEMO" !< A string for specifying the equation of state +character*(12), parameter :: EOS_ROQUET_RHO_STRING = "ROQUET_RHO" !< A string for specifying the equation of state +character*(12), parameter :: EOS_ROQUET_SPV_STRING = "ROQUET_SPV" !< A string for specifying the equation of state character*(12), parameter :: EOS_DEFAULT = EOS_WRIGHT_STRING !< The default equation of state integer, parameter :: TFREEZE_LINEAR = 1 !< A named integer specifying a freezing point expression @@ -281,6 +289,9 @@ subroutine calculate_stanley_density_scalar(T, S, pressure, Tvar, TScov, Svar, r case (EOS_NEMO) call calculate_density_second_derivs_NEMO(T_scale*T, S_scale*S, p_scale*pressure, & d2RdSS, d2RdST, d2RdTT, d2RdSp, d2RdTP) + case (EOS_ROQUET_SPV) + call calculate_density_second_derivs_Roquet_SpV(T_scale*T, S_scale*S, p_scale*pressure, & + d2RdSS, d2RdST, d2RdTT, d2RdSp, d2RdTP) case (EOS_TEOS10) call calculate_density_second_derivs_teos10(T_scale*T, S_scale*S, p_scale*pressure, & d2RdSS, d2RdST, d2RdTT, d2RdSp, d2RdTP) @@ -327,7 +338,9 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re case (EOS_TEOS10) call calculate_density_teos10(T, S, pressure, rho, start, npts, rho_ref) case (EOS_NEMO) - call calculate_density_nemo(T, S, pressure, rho, start, npts, rho_ref) + call calculate_density_nemo(T, S, pressure, rho, start, npts, rho_ref) + case (EOS_ROQUET_SPV) + call calculate_density_Roquet_SpV(T, S, pressure, rho, start, npts, rho_ref) case default call MOM_error(FATAL, "calculate_density_array: EOS%form_of_EOS is not valid.") end select @@ -397,6 +410,10 @@ subroutine calculate_stanley_density_array(T, S, pressure, Tvar, TScov, Svar, rh call calculate_density_NEMO(T, S, pressure, rho, start, npts, rho_ref) call calculate_density_second_derivs_NEMO(T, S, pressure, d2RdSS, d2RdST, & d2RdTT, d2RdSp, d2RdTP, start, npts) + case (EOS_ROQUET_SPV) + call calculate_density_Roquet_SpV(T, S, pressure, rho, start, npts, rho_ref) + call calculate_density_second_derivs_Roquet_SpV(T, S, pressure, d2RdSS, d2RdST, & + d2RdTT, d2RdSp, d2RdTP, start, npts) case (EOS_TEOS10) call calculate_density_teos10(T, S, pressure, rho, start, npts, rho_ref) call calculate_density_second_derivs_teos10(T, S, pressure, d2RdSS, d2RdST, & @@ -557,6 +574,10 @@ subroutine calculate_stanley_density_1d(T, S, pressure, Tvar, TScov, Svar, rho, call calculate_density_NEMO(Ta, Sa, pres, rho, is, npts, rho_reference) call calculate_density_second_derivs_NEMO(Ta, Sa, pres, d2RdSS, d2RdST, & d2RdTT, d2RdSp, d2RdTP, is, npts) + case (EOS_ROQUET_SPV) + call calculate_density_Roquet_SpV(Ta, Sa, pres, rho, is, npts, rho_reference) + call calculate_density_second_derivs_Roquet_SpV(Ta, Sa, pres, d2RdSS, d2RdST, & + d2RdTT, d2RdSp, d2RdTP, is, npts) case (EOS_TEOS10) call calculate_density_teos10(Ta, Sa, pres, rho, is, npts, rho_reference) call calculate_density_second_derivs_teos10(Ta, Sa, pres, d2RdSS, d2RdST, & @@ -618,6 +639,8 @@ subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, s else specvol(:) = 1.0 / rho(:) endif + case (EOS_ROQUET_SpV) + call calculate_spec_vol_Roquet_SpV(T, S, pressure, specvol, start, npts, spv_ref) case default call MOM_error(FATAL, "calculate_spec_vol_array: EOS%form_of_EOS is not valid.") end select @@ -904,6 +927,8 @@ subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, star 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 (EOS_ROQUET_SPV) + call calculate_density_derivs_Roquet_SpV(T, S, pressure, drho_dT, drho_dS, start, npts) case default call MOM_error(FATAL, "calculate_density_derivs_array: EOS%form_of_EOS is not valid.") end select @@ -1085,6 +1110,9 @@ subroutine calculate_density_second_derivs_1d(T, S, pressure, drho_dS_dS, drho_d case (EOS_NEMO) call calculate_density_second_derivs_NEMO(T, S, pressure, drho_dS_dS, drho_dS_dT, & drho_dT_dT, drho_dS_dP, drho_dT_dP, is, npts) + case (EOS_ROQUET_SPV) + call calculate_density_second_derivs_Roquet_SpV(T, S, pressure, drho_dS_dS, drho_dS_dT, & + drho_dT_dT, drho_dS_dP, drho_dT_dP, is, npts) case (EOS_TEOS10) call calculate_density_second_derivs_teos10(T, S, pressure, drho_dS_dS, drho_dS_dT, & drho_dT_dT, drho_dS_dP, drho_dT_dP, is, npts) @@ -1121,6 +1149,9 @@ subroutine calculate_density_second_derivs_1d(T, S, pressure, drho_dS_dS, drho_d case (EOS_NEMO) call calculate_density_second_derivs_NEMO(Ta, Sa, pres, drho_dS_dS, drho_dS_dT, & drho_dT_dT, drho_dS_dP, drho_dT_dP, is, npts) + case (EOS_ROQUET_SpV) + call calculate_density_second_derivs_Roquet_SpV(Ta, Sa, pres, drho_dS_dS, drho_dS_dT, & + drho_dT_dT, drho_dS_dP, drho_dT_dP, is, npts) case (EOS_TEOS10) call calculate_density_second_derivs_teos10(Ta, Sa, pres, drho_dS_dS, drho_dS_dT, & drho_dT_dT, drho_dS_dP, drho_dT_dP, is, npts) @@ -1211,6 +1242,9 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr case (EOS_NEMO) call calculate_density_second_derivs_NEMO(Ta, Sa, pres, drho_dS_dS, drho_dS_dT, & drho_dT_dT, drho_dS_dP, drho_dT_dP) + case (EOS_ROQUET_SPV) + call calculate_density_second_derivs_Roquet_SpV(Ta, Sa, pres, drho_dS_dS, drho_dS_dT, & + drho_dT_dT, drho_dS_dP, drho_dT_dP) case (EOS_TEOS10) call calculate_density_second_derivs_teos10(Ta, Sa, pres, drho_dS_dS, drho_dS_dT, & drho_dT_dT, drho_dS_dP, drho_dT_dP) @@ -1292,6 +1326,8 @@ subroutine calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, start dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) enddo + case (EOS_ROQUET_SPV) + call calculate_specvol_derivs_Roquet_SpV(T, S, pressure, dSV_dT, dSV_dS, start, npts) case default call MOM_error(FATAL, "calculate_spec_vol_derivs_array: EOS%form_of_EOS is not valid.") end select @@ -1400,6 +1436,8 @@ subroutine calculate_compress_1d(T, S, pressure, rho, drho_dp, EOS, dom) call calculate_compress_teos10(Ta, Sa, pres, rho, drho_dp, is, npts) case (EOS_NEMO) call calculate_compress_nemo(Ta, Sa, pres, rho, drho_dp, is, npts) + case (EOS_ROQUET_SpV) + call calculate_compress_Roquet_SpV(Ta, Sa, pres, rho, drho_dp, is, npts) case default call MOM_error(FATAL, "calculate_compress: EOS%form_of_EOS is not valid.") end select @@ -1511,7 +1549,6 @@ subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & real :: dRdS_scale ! A factor to convert drho_dS to the desired units [R ppt m3 S-1 kg-1 ~> 1] - ! We should never reach this point with quadrature. EOS_quadrature indicates that numerical ! integration be used instead of analytic. This is a safety check. if (EOS%EOS_quadrature) call MOM_error(FATAL, "EOS_quadrature is set!") @@ -1686,14 +1723,16 @@ subroutine EOS_init(param_file, EOS, US) call get_param(param_file, mdl, "EQN_OF_STATE", tmpstr, & "EQN_OF_STATE determines which ocean equation of state should be used. "//& - 'Currently, the valid choices are "LINEAR", "UNESCO", '//& - '"WRIGHT", "WRIGHT_RED", "WRIGHT_FULL", "NEMO" and "TEOS10". '//& - "This is only used if USE_EOS is true.", default=EOS_DEFAULT) + 'Currently, the valid choices are "LINEAR", "UNESCO", "JACKETT_MCD", '//& + '"WRIGHT", "WRIGHT_RED", "WRIGHT_FULL", "NEMO", "ROQUET_RHO", "ROQUET_SPV" '//& + 'and "TEOS10". This is only used if USE_EOS is true.', default=EOS_DEFAULT) select case (uppercase(tmpstr)) case (EOS_LINEAR_STRING) EOS%form_of_EOS = EOS_LINEAR case (EOS_UNESCO_STRING) EOS%form_of_EOS = EOS_UNESCO + case (EOS_JACKETT_STRING) + EOS%form_of_EOS = EOS_UNESCO case (EOS_WRIGHT_STRING) EOS%form_of_EOS = EOS_WRIGHT case (EOS_WRIGHT_RED_STRING) @@ -1704,6 +1743,10 @@ subroutine EOS_init(param_file, EOS, US) EOS%form_of_EOS = EOS_TEOS10 case (EOS_NEMO_STRING) EOS%form_of_EOS = EOS_NEMO + case (EOS_ROQUET_RHO_STRING) + EOS%form_of_EOS = EOS_NEMO + case (EOS_ROQUET_SPV_STRING) + EOS%form_of_EOS = EOS_ROQUET_SPV case default call MOM_error(FATAL, "interpret_eos_selection: EQN_OF_STATE "//& trim(tmpstr) // " in input file is invalid.") @@ -1741,7 +1784,8 @@ subroutine EOS_init(param_file, EOS, US) "code for the integrals of density.", default=EOS_quad_default) TFREEZE_DEFAULT = TFREEZE_LINEAR_STRING - if ((EOS%form_of_EOS == EOS_TEOS10 .or. EOS%form_of_EOS == EOS_NEMO)) & + if ((EOS%form_of_EOS == EOS_TEOS10 .or. EOS%form_of_EOS == EOS_NEMO .or. & + EOS%form_of_EOS == EOS_ROQUET_SPV)) & TFREEZE_DEFAULT = TFREEZE_TEOS10_STRING call get_param(param_file, mdl, "TFREEZE_FORM", tmpstr, & "TFREEZE_FORM determines which expression should be "//& @@ -1777,9 +1821,9 @@ subroutine EOS_init(param_file, EOS, US) units="deg C Pa-1", default=0.0) endif - if ((EOS%form_of_EOS == EOS_TEOS10 .or. EOS%form_of_EOS == EOS_NEMO) .and. & + if ((EOS%form_of_EOS == EOS_TEOS10 .or. EOS%form_of_EOS == EOS_NEMO .or. EOS%form_of_EOS == EOS_ROQUET_SPV) .and. & (EOS%form_of_TFreeze /= TFREEZE_TEOS10)) then - call MOM_error(FATAL, "interpret_eos_selection: EOS_TEOS10 or EOS_NEMO "//& + call MOM_error(FATAL, "interpret_eos_selection: EOS_TEOS10 or EOS_NEMO or EOS_ROQUET_SPV "//& "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 .") endif @@ -1870,7 +1914,8 @@ subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS) real :: gsw_ct_from_pt ! Conservative temperature after conversion from potential temperature [degC] integer :: i, j, k - if ((EOS%form_of_EOS /= EOS_TEOS10) .and. (EOS%form_of_EOS /= EOS_NEMO)) return + if ((EOS%form_of_EOS /= EOS_TEOS10) .and. (EOS%form_of_EOS /= EOS_NEMO) .and. & + (EOS%form_of_EOS /= EOS_ROQUET_SPV)) return do k=1,kd ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec if (mask_z(i,j,k) >= 1.0) then @@ -1886,7 +1931,7 @@ end subroutine convert_temp_salt_for_TEOS10 !> Converts an array of conservative temperatures to potential temperatures. The input arguments -!! use the dimesionally rescaling as specified within the EOS type. The output potential +!! use the dimensionally rescaling as specified within the EOS type. The output potential !! temperature uses this same scaling, but this can be replaced by the factor given by scale. subroutine cons_temp_to_pot_temp(T, S, poTemp, EOS, dom, scale) real, dimension(:), intent(in) :: T !< Conservative temperature [C ~> degC] @@ -1933,7 +1978,7 @@ end subroutine cons_temp_to_pot_temp !> Converts an array of absolute salinity to practical salinity. The input arguments -!! use the dimesionally rescaling as specified within the EOS type. The output potential +!! use the dimensionally rescaling as specified within the EOS type. The output potential !! temperature uses this same scaling, but this can be replaced by the factor given by scale. subroutine abs_saln_to_prac_saln(S, prSaln, EOS, dom, scale) real, dimension(:), intent(in) :: S !< Absolute salinity [S ~> ppt] @@ -2029,44 +2074,65 @@ logical function EOS_unit_tests(verbose) call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_UNESCO) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "UNESCO", & - rho_check=1027.5434579611974*EOS_tmp%kg_m3_to_R) + rho_check=1027.54345796120*EOS_tmp%kg_m3_to_R) if (verbose .and. fail) call MOM_error(WARNING, "UNESCO EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT_FULL) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT_FULL", & - rho_check=1027.5517744761617*EOS_tmp%kg_m3_to_R) + rho_check=1027.55177447616*EOS_tmp%kg_m3_to_R) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT_FULL EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT_RED) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT_RED", & - rho_check=1027.5430359634624*EOS_tmp%kg_m3_to_R) + rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT_RED EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT", & - rho_check=1027.5430359634624*EOS_tmp%kg_m3_to_R) + rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_NEMO) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "NEMO", & - rho_check=1027.4238566366823*EOS_tmp%kg_m3_to_R) + rho_check=1027.42385663668*EOS_tmp%kg_m3_to_R) if (verbose .and. fail) call MOM_error(WARNING, "NEMO EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail + call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_ROQUET_SPV) + fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "ROQUET_SPV", & + rho_check=1027.42387475199*EOS_tmp%kg_m3_to_R) + if (verbose .and. fail) call MOM_error(WARNING, "ROQUET_SPV EOS has failed some self-consistency tests.") + EOS_unit_tests = EOS_unit_tests .or. fail + ! The TEOS10 equation of state is not passing the self consistency tests for dho_dS_dp due ! to a bug (a missing division by the square root of salinity) on line 109 of ! pkg/GSW-Fortan/toolbox/gsw_specvol_second_derivatives.f90. This bug has been highlighted in an ! issue posted to the TEOS-10/GSW-Fortran page at github.com/TEOS-10/GSW-Fortran/issues/26. call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_TEOS10) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "TEOS10", skip_2nd=.true., & - rho_check=1027.4235596149185*EOS_tmp%kg_m3_to_R) + rho_check=1027.42355961492*EOS_tmp%kg_m3_to_R) if (verbose .and. fail) call MOM_error(WARNING, "TEOS10 EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail + call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_NEMO) + fail = test_EOS_consistency(10.0, 30.0, 1.0e7, EOS_tmp, verbose, "NEMO", & + rho_check=1027.45140117152*EOS_tmp%kg_m3_to_R) + ! The corresponding check value published by Roquet et al. (2015) is 1027.45140 [kg m-3]. + if (verbose .and. fail) call MOM_error(WARNING, "NEMO EOS has failed some self-consistency tests.") + EOS_unit_tests = EOS_unit_tests .or. fail + + call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_ROQUET_SPV) + fail = test_EOS_consistency(10.0, 30.0, 1.0e7, EOS_tmp, verbose, "ROQUET_SPV", & + spv_check=9.73282046614623e-04*EOS_tmp%R_to_kg_m3) + ! The corresponding check value here published by Roquet et al. (2015) is 9.732819628e-04 [m3 kg-1], + ! but the order of arithmetic there was not completely specified with parentheses. + if (verbose .and. fail) call MOM_error(WARNING, "ROQUET_SPV EOS has failed some self-consistency tests.") + EOS_unit_tests = EOS_unit_tests .or. fail + call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_LINEAR, Rho_T0_S0=1000.0, drho_dT=-0.2, dRho_dS=0.8) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "LINEAR", & rho_check=1023.0*EOS_tmp%kg_m3_to_R) @@ -2094,7 +2160,7 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & ! Local variables real, dimension(-3:3,-3:3,-3:3) :: T ! Temperatures at the test value and perturbed points [C ~> degC] - real, dimension(-3:3,-3:3,-3:3) :: S ! Salinites at the test value and perturbed points [S ~> ppt] + real, dimension(-3:3,-3:3,-3:3) :: S ! Salinities at the test value and perturbed points [S ~> ppt] real, dimension(-3:3,-3:3,-3:3) :: P ! Pressures at the test value and perturbed points [R L2 T-2 ~> Pa] real, dimension(-3:3,-3:3,-3:3,2) :: rho ! Densities relative to rho_ref at the test value and ! perturbed points [R ~> kg m-3] @@ -2176,7 +2242,7 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & do n=1,2 ! Calculate density values with a wide enough stencil to estimate first and second derivatives - ! with up to 6th order accuracy. Doing this twice with different sizes of pertubations allows + ! with up to 6th order accuracy. Doing this twice with different sizes of perturbations allows ! the evaluation of whether the finite differences are converging to the calculated values at a ! rate that is consistent with the order of accuracy of the finite difference forms, and hence ! the consistency of the calculated values. @@ -2300,9 +2366,9 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & contains - !> Return a finite difference estimate of the first derivative of a field in arbitary units [A B-1] + !> Return a finite difference estimate of the first derivative of a field in arbitrary units [A B-1] real function first_deriv(R, dx, order) - real, intent(in) :: R(-3:3) !< The field whose derivative is being taken, in abitrary units [A] + real, intent(in) :: R(-3:3) !< The field whose derivative is being taken, in arbitrary units [A] real, intent(in) :: dx !< The spacing in parameter space, in different arbitrary units [B] integer, intent(in) :: order !< The order of accuracy of the centered finite difference estimates (2, 4 or 6) @@ -2315,9 +2381,9 @@ real function first_deriv(R, dx, order) endif end function first_deriv - !> Return a finite difference estimate of the second derivative of a field in arbitary units [A B-2] + !> Return a finite difference estimate of the second derivative of a field in arbitrary units [A B-2] real function second_deriv(R, dx, order) - real, intent(in) :: R(-3:3) !< The field whose derivative is being taken, in abitrary units [A] + real, intent(in) :: R(-3:3) !< The field whose derivative is being taken, in arbitrary units [A] real, intent(in) :: dx !< The spacing in parameter space, in different arbitrary units [B] integer, intent(in) :: order !< The order of accuracy of the centered finite difference estimates (2, 4 or 6) @@ -2331,9 +2397,9 @@ real function second_deriv(R, dx, order) end function second_deriv !> Return a finite difference estimate of the second derivative with respect to two different - !! parameters of a field in arbitary units [A B-2] + !! parameters of a field in arbitrary units [A B-1 C-1] real function derivs_2d(R, dxdy, order) - real, intent(in) :: R(-3:3,-3:3) !< The field whose derivative is being taken in abitrary units [A] + real, intent(in) :: R(-3:3,-3:3) !< The field whose derivative is being taken in arbitrary units [A] real, intent(in) :: dxdy !< The spacing in two directions in parameter space in different arbitrary units [B C] integer, intent(in) :: order !< The order of accuracy of the centered finite difference estimates (2, 4 or 6) @@ -2386,6 +2452,6 @@ end module MOM_EOS !> \namespace mom_eos !! -!! The MOM_EOS module is a wrapper for various equations of state (e.g. Linear, -!! Wright, UNESCO, TEOS10 or NEMO) and provides a uniform interface to the rest of the model -!! independent of which equation of state is being used. +!! The MOM_EOS module is a wrapper for various equations of state (i.e. Linear, Wright, +!! Wright_full, Wright_red, UNESCO, TEOS10, Roquet_SpV or NEMO) and provides a uniform +!! interface to the rest of the model independent of which equation of state is being used. diff --git a/src/equation_of_state/MOM_EOS_Roquet_SpV.F90 b/src/equation_of_state/MOM_EOS_Roquet_SpV.F90 new file mode 100644 index 0000000000..5a276065dd --- /dev/null +++ b/src/equation_of_state/MOM_EOS_Roquet_SpV.F90 @@ -0,0 +1,790 @@ +!> The equation of state for specific volume (SpV) using the expressions of Roquet et al. 2015 +module MOM_EOS_Roquet_Spv + +! This file is part of MOM6. See LICENSE.md for the license. + +!use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt + +implicit none ; private + +public calculate_compress_Roquet_SpV, calculate_density_Roquet_SpV, calculate_spec_vol_Roquet_SpV +public calculate_density_derivs_Roquet_SpV, calculate_specvol_derivs_Roquet_SpV +public calculate_density_scalar_Roquet_SpV, calculate_density_array_Roquet_SpV +public calculate_density_second_derivs_Roquet_SpV + +!> Compute the in situ density of sea water [kg m-3], or its anomaly with respect to +!! a reference density, from absolute salinity [g kg-1], conservative temperature [degC], +!! and pressure [Pa], using the specific volume polynomial fit from Roquet et al. (2015) +interface calculate_density_Roquet_SpV + module procedure calculate_density_scalar_Roquet_SpV, calculate_density_array_Roquet_SpV +end interface calculate_density_Roquet_SpV + +!> Compute the in situ specific volume of sea water (in [m3 kg-1]), or an anomaly with respect +!! to a reference specific volume, from absolute salinity ([g kg-1]), conservative +!! temperature (in degrees Celsius [degC]), and pressure [Pa], using the specific volume +!! polynomial fit from Roquet et al. (2015) +interface calculate_spec_vol_Roquet_SpV + module procedure calculate_spec_vol_scalar_Roquet_SpV, calculate_spec_vol_array_Roquet_SpV +end interface calculate_spec_vol_Roquet_SpV + +!> For a given thermodynamic state, return the derivatives of density with conservative temperature +!! and absolute salinity, using the specific volume polynomial fit from Roquet et al. (2015) +interface calculate_density_derivs_Roquet_SpV + module procedure calculate_density_derivs_scalar_Roquet_SpV, calculate_density_derivs_array_Roquet_SpV +end interface calculate_density_derivs_Roquet_SpV + +!> Compute the second derivatives of density with various combinations of temperature, salinity +!! and pressure using the specific volume polynomial fit from Roquet et al. (2015) +interface calculate_density_second_derivs_Roquet_SpV + module procedure calculate_density_second_derivs_scalar_Roquet_SpV + module procedure calculate_density_second_derivs_array_Roquet_SpV +end interface calculate_density_second_derivs_Roquet_SpV + +real, parameter :: Pa2db = 1.e-4 !< Conversion factor between Pa and dbar [dbar Pa-1] +!>@{ Parameters in the Roquet specific volume polynomial equation of state +real, parameter :: rdeltaS = 24. ! An offset to salinity before taking its square root [g kg-1] +real, parameter :: r1_S0 = 0.875/35.16504 ! The inverse of a plausible range of oceanic salinities [kg g-1] +real, parameter :: r1_T0 = 1./40. ! The inverse of a plausible range of oceanic temperatures [degC-1] +real, parameter :: r1_P0 = 1.e-4 ! The inverse of a plausible range of oceanic pressures [dbar-1] +real, parameter :: V00 = -4.4015007269e-05 ! Contribution to SpV00p proportional to zp [m3 kg-1] +real, parameter :: V01 = 6.9232335784e-06 ! Contribution to SpV00p proportional to zp**2 [m3 kg-1] +real, parameter :: V02 = -7.5004675975e-07 ! Contribution to SpV00p proportional to zp**3 [m3 kg-1] +real, parameter :: V03 = 1.7009109288e-08 ! Contribution to SpV00p proportional to zp**4 [m3 kg-1] +real, parameter :: V04 = -1.6884162004e-08 ! Contribution to SpV00p proportional to zp**5 [m3 kg-1] +real, parameter :: V05 = 1.9613503930e-09 ! Contribution to SpV00p proportional to zp**6 [m3 kg-1] + +! The following terms are contributions to specific volume as a function of the normalized square root of salinity +! with an offset (zs), temperature (zt) and pressure (zp), with a contribution SPVabc * zs**a * zt**b * zp**c +real, parameter :: SPV000 = 1.0772899069e-03 ! A constant specific volume (SpV) contribution [m3 kg-1] +real, parameter :: SPV100 = -3.1263658781e-04 ! Coefficient of SpV proportional to zs [m3 kg-1] +real, parameter :: SPV200 = 6.7615860683e-04 ! Coefficient of SpV proportional to zs**2 [m3 kg-1] +real, parameter :: SPV300 = -8.6127884515e-04 ! Coefficient of SpV proportional to zs**3 [m3 kg-1] +real, parameter :: SPV400 = 5.9010812596e-04 ! Coefficient of SpV proportional to zs**4 [m3 kg-1] +real, parameter :: SPV500 = -2.1503943538e-04 ! Coefficient of SpV proportional to zs**5 [m3 kg-1] +real, parameter :: SPV600 = 3.2678954455e-05 ! Coefficient of SpV proportional to zs**6 [m3 kg-1] +real, parameter :: SPV010 = -1.4949652640e-05 ! Coefficient of SpV proportional to zt [m3 kg-1] +real, parameter :: SPV110 = 3.1866349188e-05 ! Coefficient of SpV proportional to zs * zt [m3 kg-1] +real, parameter :: SPV210 = -3.8070687610e-05 ! Coefficient of SpV proportional to zs**2 * zt [m3 kg-1] +real, parameter :: SPV310 = 2.9818473563e-05 ! Coefficient of SpV proportional to zs**3 * zt [m3 kg-1] +real, parameter :: SPV410 = -1.0011321965e-05 ! Coefficient of SpV proportional to zs**4 * zt [m3 kg-1] +real, parameter :: SPV510 = 1.0751931163e-06 ! Coefficient of SpV proportional to zs**5 * zt [m3 kg-1] +real, parameter :: SPV020 = 2.7546851539e-05 ! Coefficient of SpV proportional to zt**2 [m3 kg-1] +real, parameter :: SPV120 = -3.6597334199e-05 ! Coefficient of SpV proportional to zs * zt**2 [m3 kg-1] +real, parameter :: SPV220 = 3.4489154625e-05 ! Coefficient of SpV proportional to zs**2 * zt**2 [m3 kg-1] +real, parameter :: SPV320 = -1.7663254122e-05 ! Coefficient of SpV proportional to zs**3 * zt**2 [m3 kg-1] +real, parameter :: SPV420 = 3.5965131935e-06 ! Coefficient of SpV proportional to zs**4 * zt**2 [m3 kg-1] +real, parameter :: SPV030 = -1.6506828994e-05 ! Coefficient of SpV proportional to zt**3 [m3 kg-1] +real, parameter :: SPV130 = 2.4412359055e-05 ! Coefficient of SpV proportional to zs * zt**3 [m3 kg-1] +real, parameter :: SPV230 = -1.4606740723e-05 ! Coefficient of SpV proportional to zs**2 * zt**3 [m3 kg-1] +real, parameter :: SPV330 = 2.3293406656e-06 ! Coefficient of SpV proportional to zs**3 * zt**3 [m3 kg-1] +real, parameter :: SPV040 = 6.7896174634e-06 ! Coefficient of SpV proportional to zt**4 [m3 kg-1] +real, parameter :: SPV140 = -8.7951832993e-06 ! Coefficient of SpV proportional to zs * zt**4 [m3 kg-1] +real, parameter :: SPV240 = 4.4249040774e-06 ! Coefficient of SpV proportional to zs**2 * zt**4 [m3 kg-1] +real, parameter :: SPV050 = -7.2535743349e-07 ! Coefficient of SpV proportional to zt**5 [m3 kg-1] +real, parameter :: SPV150 = -3.4680559205e-07 ! Coefficient of SpV proportional to zs * zt**5 [m3 kg-1] +real, parameter :: SPV060 = 1.9041365570e-07 ! Coefficient of SpV proportional to zt**6 [m3 kg-1] +real, parameter :: SPV001 = -1.6889436589e-05 ! Coefficient of SpV proportional to zp [m3 kg-1] +real, parameter :: SPV101 = 2.1106556158e-05 ! Coefficient of SpV proportional to zs * zp [m3 kg-1] +real, parameter :: SPV201 = -2.1322804368e-05 ! Coefficient of SpV proportional to zs**2 * zp [m3 kg-1] +real, parameter :: SPV301 = 1.7347655458e-05 ! Coefficient of SpV proportional to zs**3 * zp [m3 kg-1] +real, parameter :: SPV401 = -4.3209400767e-06 ! Coefficient of SpV proportional to zs**4 * zp [m3 kg-1] +real, parameter :: SPV011 = 1.5355844621e-05 ! Coefficient of SpV proportional to zt * zp [m3 kg-1] +real, parameter :: SPV111 = 2.0914122241e-06 ! Coefficient of SpV proportional to zs * zt * zp [m3 kg-1] +real, parameter :: SPV211 = -5.7751479725e-06 ! Coefficient of SpV proportional to zs**2 * zt * zp [m3 kg-1] +real, parameter :: SPV311 = 1.0767234341e-06 ! Coefficient of SpV proportional to zs**3 * zt * zp [m3 kg-1] +real, parameter :: SPV021 = -9.6659393016e-06 ! Coefficient of SpV proportional to zt**2 * zp [m3 kg-1] +real, parameter :: SPV121 = -7.0686982208e-07 ! Coefficient of SpV proportional to zs * zt**2 * zp [m3 kg-1] +real, parameter :: SPV221 = 1.4488066593e-06 ! Coefficient of SpV proportional to zs**2 * zt**2 * zp [m3 kg-1] +real, parameter :: SPV031 = 3.1134283336e-06 ! Coefficient of SpV proportional to zt**3 * zp [m3 kg-1] +real, parameter :: SPV131 = 7.9562529879e-08 ! Coefficient of SpV proportional to zs * zt**3 * zp [m3 kg-1] +real, parameter :: SPV041 = -5.6590253863e-07 ! Coefficient of SpV proportional to zt * zp [m3 kg-1] +real, parameter :: SPV002 = 1.0500241168e-06 ! Coefficient of SpV proportional to zp**2 [m3 kg-1] +real, parameter :: SPV102 = 1.9600661704e-06 ! Coefficient of SpV proportional to zs * zp**2 [m3 kg-1] +real, parameter :: SPV202 = -2.1666693382e-06 ! Coefficient of SpV proportional to zs**2 * zp**2 [m3 kg-1] +real, parameter :: SPV012 = -3.8541359685e-06 ! Coefficient of SpV proportional to zt * zp**2 [m3 kg-1] +real, parameter :: SPV112 = 1.0157632247e-06 ! Coefficient of SpV proportional to zs * zt * zp**2 [m3 kg-1] +real, parameter :: SPV022 = 1.7178343158e-06 ! Coefficient of SpV proportional to zt**2 * zp**2 [m3 kg-1] +real, parameter :: SPV003 = -4.1503454190e-07 ! Coefficient of SpV proportional to zp**3 [m3 kg-1] +real, parameter :: SPV103 = 3.5627020989e-07 ! Coefficient of SpV proportional to zs * zp**3 [m3 kg-1] +real, parameter :: SPV013 = -1.1293871415e-07 ! Coefficient of SpV proportional to zt * zp**3 [m3 kg-1] + +real, parameter :: ALP000 = SPV010*r1_T0 ! Constant in the dSpV_dT fit [m3 kg-1 degC-1] +real, parameter :: ALP100 = SPV110*r1_T0 ! Coefficient of the dSpV_dT fit zs term [m3 kg-1 degC-1] +real, parameter :: ALP200 = SPV210*r1_T0 ! Coefficient of the dSpV_dT fit zs**2 term [m3 kg-1 degC-1] +real, parameter :: ALP300 = SPV310*r1_T0 ! Coefficient of the dSpV_dT fit zs**3 term [m3 kg-1 degC-1] +real, parameter :: ALP400 = SPV410*r1_T0 ! Coefficient of the dSpV_dT fit zs**4 term [m3 kg-1 degC-1] +real, parameter :: ALP500 = SPV510*r1_T0 ! Coefficient of the dSpV_dT fit zs**5 term [m3 kg-1 degC-1] +real, parameter :: ALP010 = 2.*SPV020*r1_T0 ! Coefficient of the dSpV_dT fit zt term [m3 kg-1 degC-1] +real, parameter :: ALP110 = 2.*SPV120*r1_T0 ! Coefficient of the dSpV_dT fit zs * zt term [m3 kg-1 degC-1] +real, parameter :: ALP210 = 2.*SPV220*r1_T0 ! Coefficient of the dSpV_dT fit zs**2 * zt term [m3 kg-1 degC-1] +real, parameter :: ALP310 = 2.*SPV320*r1_T0 ! Coefficient of the dSpV_dT fit zs**3 * zt term [m3 kg-1 degC-1] +real, parameter :: ALP410 = 2.*SPV420*r1_T0 ! Coefficient of the dSpV_dT fit zs**4 * zt term [m3 kg-1 degC-1] +real, parameter :: ALP020 = 3.*SPV030*r1_T0 ! Coefficient of the dSpV_dT fit zt**2 term [m3 kg-1 degC-1] +real, parameter :: ALP120 = 3.*SPV130*r1_T0 ! Coefficient of the dSpV_dT fit zs * zt**2 term [m3 kg-1 degC-1] +real, parameter :: ALP220 = 3.*SPV230*r1_T0 ! Coefficient of the dSpV_dT fit zs**2 * zt**2 term [m3 kg-1 degC-1] +real, parameter :: ALP320 = 3.*SPV330*r1_T0 ! Coefficient of the dSpV_dT fit zs**3 * zt**2 term [m3 kg-1 degC-1] +real, parameter :: ALP030 = 4.*SPV040*r1_T0 ! Coefficient of the dSpV_dT fit zt**3 term [m3 kg-1 degC-1] +real, parameter :: ALP130 = 4.*SPV140*r1_T0 ! Coefficient of the dSpV_dT fit zs * zt**3 term [m3 kg-1 degC-1] +real, parameter :: ALP230 = 4.*SPV240*r1_T0 ! Coefficient of the dSpV_dT fit zs**2 * zt**3 term [m3 kg-1 degC-1] +real, parameter :: ALP040 = 5.*SPV050*r1_T0 ! Coefficient of the dSpV_dT fit zt**4 term [m3 kg-1 degC-1] +real, parameter :: ALP140 = 5.*SPV150*r1_T0 ! Coefficient of the dSpV_dT fit zs* * zt**4 term [m3 kg-1 degC-1] +real, parameter :: ALP050 = 6.*SPV060*r1_T0 ! Coefficient of the dSpV_dT fit zt**5 term [m3 kg-1 degC-1] +real, parameter :: ALP001 = SPV011*r1_T0 ! Coefficient of the dSpV_dT fit zp term [m3 kg-1 degC-1] +real, parameter :: ALP101 = SPV111*r1_T0 ! Coefficient of the dSpV_dT fit zs * zp term [m3 kg-1 degC-1] +real, parameter :: ALP201 = SPV211*r1_T0 ! Coefficient of the dSpV_dT fit zs**2 * zp term [m3 kg-1 degC-1] +real, parameter :: ALP301 = SPV311*r1_T0 ! Coefficient of the dSpV_dT fit zs**3 * zp term [m3 kg-1 degC-1] +real, parameter :: ALP011 = 2.*SPV021*r1_T0 ! Coefficient of the dSpV_dT fit zt * zp term [m3 kg-1 degC-1] +real, parameter :: ALP111 = 2.*SPV121*r1_T0 ! Coefficient of the dSpV_dT fit zs * zt * zp term [m3 kg-1 degC-1] +real, parameter :: ALP211 = 2.*SPV221*r1_T0 ! Coefficient of the dSpV_dT fit zs**2 * zt * zp term [m3 kg-1 degC-1] +real, parameter :: ALP021 = 3.*SPV031*r1_T0 ! Coefficient of the dSpV_dT fit zt**2 * zp term [m3 kg-1 degC-1] +real, parameter :: ALP121 = 3.*SPV131*r1_T0 ! Coefficient of the dSpV_dT fit zs * zt**2 * zp term [m3 kg-1 degC-1] +real, parameter :: ALP031 = 4.*SPV041*r1_T0 ! Coefficient of the dSpV_dT fit zt**3 * zp term [m3 kg-1 degC-1] +real, parameter :: ALP002 = SPV012*r1_T0 ! Coefficient of the dSpV_dT fit zp**2 term [m3 kg-1 degC-1] +real, parameter :: ALP102 = SPV112*r1_T0 ! Coefficient of the dSpV_dT fit zs * zp**2 term [m3 kg-1 degC-1] +real, parameter :: ALP012 = 2.*SPV022*r1_T0 ! Coefficient of the dSpV_dT fit zt * zp**2 term [m3 kg-1 degC-1] +real, parameter :: ALP003 = SPV013*r1_T0 ! Coefficient of the dSpV_dT fit zp**3 term [m3 kg-1 degC-1] + +real, parameter :: BET000 = 0.5*SPV100*r1_S0 ! Constant in the dSpV_dS fit [m3 kg-1 ppt-1] +real, parameter :: BET100 = SPV200*r1_S0 ! Coefficient of the dSpV_dS fit zs term [m3 kg-1 ppt-1] +real, parameter :: BET200 = 1.5*SPV300*r1_S0 ! Coefficient of the dSpV_dS fit zs**2 term [m3 kg-1 ppt-1] +real, parameter :: BET300 = 2.0*SPV400*r1_S0 ! Coefficient of the dSpV_dS fit zs**3 term [m3 kg-1 ppt-1] +real, parameter :: BET400 = 2.5*SPV500*r1_S0 ! Coefficient of the dSpV_dS fit zs**4 term [m3 kg-1 ppt-1] +real, parameter :: BET500 = 3.0*SPV600*r1_S0 ! Coefficient of the dSpV_dS fit zs**5 term [m3 kg-1 ppt-1] +real, parameter :: BET010 = 0.5*SPV110*r1_S0 ! Coefficient of the dSpV_dS fit zt term [m3 kg-1 ppt-1] +real, parameter :: BET110 = SPV210*r1_S0 ! Coefficient of the dSpV_dS fit zs * zt term [m3 kg-1 ppt-1] +real, parameter :: BET210 = 1.5*SPV310*r1_S0 ! Coefficient of the dSpV_dS fit zs**2 * zt term [m3 kg-1 ppt-1] +real, parameter :: BET310 = 2.0*SPV410*r1_S0 ! Coefficient of the dSpV_dS fit zs**3 * zt term [m3 kg-1 ppt-1] +real, parameter :: BET410 = 2.5*SPV510*r1_S0 ! Coefficient of the dSpV_dS fit zs**4 * zt term [m3 kg-1 ppt-1] +real, parameter :: BET020 = 0.5*SPV120*r1_S0 ! Coefficient of the dSpV_dS fit zt**2 term [m3 kg-1 ppt-1] +real, parameter :: BET120 = SPV220*r1_S0 ! Coefficient of the dSpV_dS fit zs * zt**2 term [m3 kg-1 ppt-1] +real, parameter :: BET220 = 1.5*SPV320*r1_S0 ! Coefficient of the dSpV_dS fit zs**2 * zt**2 term [m3 kg-1 ppt-1] +real, parameter :: BET320 = 2.0*SPV420*r1_S0 ! Coefficient of the dSpV_dS fit zs**3 * zt**2 term [m3 kg-1 ppt-1] +real, parameter :: BET030 = 0.5*SPV130*r1_S0 ! Coefficient of the dSpV_dS fit zt**3 term [m3 kg-1 ppt-1] +real, parameter :: BET130 = SPV230*r1_S0 ! Coefficient of the dSpV_dS fit zs * zt**3 term [m3 kg-1 ppt-1] +real, parameter :: BET230 = 1.5*SPV330*r1_S0 ! Coefficient of the dSpV_dS fit zs**2 * zt**3 term [m3 kg-1 ppt-1] +real, parameter :: BET040 = 0.5*SPV140*r1_S0 ! Coefficient of the dSpV_dS fit zt**4 term [m3 kg-1 ppt-1] +real, parameter :: BET140 = SPV240*r1_S0 ! Coefficient of the dSpV_dS fit zs * zt**4 term [m3 kg-1 ppt-1] +real, parameter :: BET050 = 0.5*SPV150*r1_S0 ! Coefficient of the dSpV_dS fit zt**5 term [m3 kg-1 ppt-1] +real, parameter :: BET001 = 0.5*SPV101*r1_S0 ! Coefficient of the dSpV_dS fit zp term [m3 kg-1 ppt-1] +real, parameter :: BET101 = SPV201*r1_S0 ! Coefficient of the dSpV_dS fit zs * zp term [m3 kg-1 ppt-1] +real, parameter :: BET201 = 1.5*SPV301*r1_S0 ! Coefficient of the dSpV_dS fit zs**2 * zp term [m3 kg-1 ppt-1] +real, parameter :: BET301 = 2.0*SPV401*r1_S0 ! Coefficient of the dSpV_dS fit zs**3 * zp term [m3 kg-1 ppt-1] +real, parameter :: BET011 = 0.5*SPV111*r1_S0 ! Coefficient of the dSpV_dS fit zt * zp term [m3 kg-1 ppt-1] +real, parameter :: BET111 = SPV211*r1_S0 ! Coefficient of the dSpV_dS fit zs * zt * zp term [m3 kg-1 ppt-1] +real, parameter :: BET211 = 1.5*SPV311*r1_S0 ! Coefficient of the dSpV_dS fit zs**2 * zt * zp term [m3 kg-1 ppt-1] +real, parameter :: BET021 = 0.5*SPV121*r1_S0 ! Coefficient of the dSpV_dS fit zt**2 * zp term [m3 kg-1 ppt-1] +real, parameter :: BET121 = SPV221*r1_S0 ! Coefficient of the dSpV_dS fit zs * zt**2 * zp term [m3 kg-1 ppt-1] +real, parameter :: BET031 = 0.5*SPV131*r1_S0 ! Coefficient of the dSpV_dS fit zt**3 * zp term [m3 kg-1 ppt-1] +real, parameter :: BET002 = 0.5*SPV102*r1_S0 ! Coefficient of the dSpV_dS fit zp**2 term [m3 kg-1 ppt-1] +real, parameter :: BET102 = SPV202*r1_S0 ! Coefficient of the dSpV_dS fit zs * zp**2 term [m3 kg-1 ppt-1] +real, parameter :: BET012 = 0.5*SPV112*r1_S0 ! Coefficient of the dSpV_dS fit zt * zp**2 term [m3 kg-1 ppt-1] +real, parameter :: BET003 = 0.5*SPV103*r1_S0 ! Coefficient of the dSpV_dS fit zp**3 term [m3 kg-1 ppt-1] +!>@} + +contains + +!> Computes the Roquet et al. in situ specific volume of sea water for scalar inputs and outputs. +!! +!! Returns the in situ specific volume of sea water (specvol in [m3 kg-1]) from absolute salinity (S [g kg-1]), +!! conservative temperature (T [degC]) and pressure [Pa]. It uses the specific volume polynomial +!! fit from Roquet et al. (2015). +!! If spv_ref is present, specvol is an anomaly from spv_ref. +subroutine calculate_spec_vol_scalar_Roquet_SpV(T, S, pressure, specvol, spv_ref) + real, intent(in) :: T !< Conservative temperature [degC] + real, intent(in) :: S !< Absolute salinity [g kg-1] + real, intent(in) :: pressure !< Pressure [Pa] + real, intent(out) :: specvol !< In situ specific volume [m3 kg-1] + real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1] + + ! Local variables + real, dimension(1) :: T0 ! A 1-d array with a copy of the conservative temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the absolutes salinity [g kg-1] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: spv0 ! A 1-d array with a copy of the specific volume [m3 kg-1] + + T0(1) = T ; S0(1) = S ; pressure0(1) = pressure + + call calculate_spec_vol_array_Roquet_SpV(T0, S0, pressure0, spv0, 1, 1, spv_ref) + specvol = spv0(1) + +end subroutine calculate_spec_vol_scalar_Roquet_SpV + +!> Computes the Roquet et al. in situ specific volume of sea water for 1-d array inputs and outputs. +!! +!! Returns the in situ specific volume of sea water (specvol in [m3 kg-1]) from absolute salinity (S [g kg-1]), +!! conservative temperature (T [degC]) and pressure [Pa]. It uses the specific volume polynomial +!! fit from Roquet et al. (2015). +!! If spv_ref is present, specvol is an anomaly from spv_ref. +subroutine calculate_spec_vol_array_Roquet_SpV(T, S, pressure, specvol, start, npts, spv_ref) + real, dimension(:), intent(in) :: T !< Conservative temperature [degC] + real, dimension(:), intent(in) :: S !< Absolute salinity [g kg-1] + real, dimension(:), intent(in) :: pressure !< pressure [Pa] + real, dimension(:), intent(inout) :: specvol !< in situ specific volume [m3 kg-1] + integer, intent(in) :: start !< The starting index for calculations + integer, intent(in) :: npts !< the number of values to calculate + real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1] + + ! Local variables + real :: zp ! Pressure, first in [dbar], then normalized by an assumed pressure range [nondim] + real :: zt ! Conservative temperature, first in [degC], then normalized by an assumed temperature range [nondim] + real :: zs ! Absolute salinity, first in [g kg-1], then the square root of salinity with an offset normalized + ! by an assumed salinity range [nondim] + real :: SV_00p ! A pressure-dependent but temperature and salinity independent contribution to + ! specific volume at the reference temperature and salinity [m3 kg-1] + real :: SV_TS ! Specific volume without a pressure-dependent contribution [m3 kg-1] + real :: SV_TS0 ! A contribution to specific volume from temperature and salinity anomalies at + ! the surface pressure [m3 kg-1] + real :: SV_TS1 ! A temperature and salinity dependent specific volume contribution that is + ! proportional to pressure [m3 kg-1] + real :: SV_TS2 ! A temperature and salinity dependent specific volume contribution that is + ! proportional to pressure**2 [m3 kg-1] + real :: SV_TS3 ! A temperature and salinity dependent specific volume contribution that is + ! proportional to pressure**3 [m3 kg-1] + real :: SV_0S0 ! Salinity dependent specific volume at the surface pressure and zero temperature [m3 kg-1] + integer :: j + + ! The following algorithm was published by Roquet et al. (2015), intended for use in non-Boussinesq ocean models. + do j=start,start+npts-1 + ! Conversions to the units used here. + zt = T(j) * r1_T0 ! Conservative temperature normalized by a plausible oceanic range [nondim] + zs = SQRT( ABS( S(j) + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim] + zp = pressure(j) * (Pa2db*r1_P0) ! Convert pressure from Pascals to kilobars to normalize it [nondim] + + ! The next two lines should be used if it is necessary to convert potential temperature and + ! practical salinity to conservative temperature and absolute salinity. + ! zt = r1_T0 * gsw_ct_from_pt(S(j),T(j)) ! Convert potential temp to conservative temp [degC] + ! zs = SQRT( ABS( gsw_sr_from_sp(S(j)) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity. + + SV_TS3 = SPV003 + (zs*SPV103 + zt*SPV013) + SV_TS2 = SPV002 + (zs*(SPV102 + zs*SPV202) & + + zt*(SPV012 + (zs*SPV112 + zt*SPV022)) ) + SV_TS1 = SPV001 + (zs*(SPV101 + zs*(SPV201 + zs*(SPV301 + zs*SPV401))) & + + zt*(SPV011 + (zs*(SPV111 + zs*(SPV211 + zs*SPV311)) & + + zt*(SPV021 + (zs*(SPV121 + zs*SPV221) & + + zt*(SPV031 + (zs*SPV131 + zt*SPV041)) )) )) ) + SV_TS0 = zt*(SPV010 & + + (zs*(SPV110 + zs*(SPV210 + zs*(SPV310 + zs*(SPV410 + zs*SPV510)))) & + + zt*(SPV020 + (zs*(SPV120 + zs*(SPV220 + zs*(SPV320 + zs*SPV420))) & + + zt*(SPV030 + (zs*(SPV130 + zs*(SPV230 + zs*SPV330)) & + + zt*(SPV040 + (zs*(SPV140 + zs*SPV240) & + + zt*(SPV050 + (zs*SPV150 + zt*SPV060)) )) )) )) ) ) + + SV_0S0 = SPV000 + zs*(SPV100 + zs*(SPV200 + zs*(SPV300 + zs*(SPV400 + zs*(SPV500 + zs*SPV600))))) + + SV_00p = zp*(V00 + zp*(V01 + zp*(V02 + zp*(V03 + zp*(V04 + zp*V05))))) + + if (present(spv_ref)) SV_0S0 = SV_0S0 - spv_ref + + SV_TS = (SV_TS0 + SV_0S0) + zp*(SV_TS1 + zp*(SV_TS2 + zp*SV_TS3)) + specvol(j) = SV_TS + SV_00p ! In situ specific volume [m3 kg-1] + enddo + +end subroutine calculate_spec_vol_array_Roquet_SpV + + +!> Compute the in situ density of sea water at a point (rho in [kg m-3]) from absolute +!! salinity (S [g kg-1]), conservative temperature (T [degC]) and pressure [Pa], using the +!! specific volume polynomial fit from Roquet et al. (2015). +subroutine calculate_density_scalar_Roquet_SpV(T, S, pressure, rho, rho_ref) + real, intent(in) :: T !< Conservative temperature [degC] + real, intent(in) :: S !< Absolute salinity [g kg-1] + real, intent(in) :: pressure !< Pressure [Pa] + real, intent(out) :: rho !< In situ density [kg m-3] + real, optional, intent(in) :: rho_ref !< A reference density [kg m-3] + + real, dimension(1) :: T0 ! A 1-d array with a copy of the conservative temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the absolute salinity [g kg-1] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: spv ! A 1-d array with the specific volume [m3 kg-1] + + T0(1) = T + S0(1) = S + pressure0(1) = pressure + + if (present(rho_ref)) then + call calculate_spec_vol_array_Roquet_SpV(T0, S0, pressure0, spv, 1, 1, spv_ref=1.0/rho_ref) + rho = -rho_ref**2*spv(1) / (rho_ref*spv(1) + 1.0) ! In situ density [kg m-3] + else + call calculate_spec_vol_array_Roquet_SpV(T0, S0, pressure0, spv, 1, 1) + rho = 1.0 / spv(1) + endif + +end subroutine calculate_density_scalar_Roquet_SpV + +!> Compute an array of in situ densities of sea water (rho in [kg m-3]) from absolute +!! salinity (S [g kg-1]), conservative temperature (T [degC]) and pressure [Pa], +!! using the specific volume polynomial fit from Roquet et al. (2015). +subroutine calculate_density_array_Roquet_SpV(T, S, pressure, rho, start, npts, rho_ref) + real, dimension(:), intent(in) :: T !< Conservative temperature [degC] + real, dimension(:), intent(in) :: S !< Absolute salinity [g kg-1] + real, dimension(:), intent(in) :: pressure !< Pressure [Pa] + real, dimension(:), intent(out) :: rho !< In situ density [kg m-3] + integer, intent(in) :: start !< The starting index for calculations + integer, intent(in) :: npts !< The number of values to calculate + real, optional, intent(in) :: rho_ref !< A reference density [kg m-3] + + ! Local variables + real, dimension(size(T)) :: spv ! The specific volume [m3 kg-1] + integer :: j + + if (present(rho_ref)) then + call calculate_spec_vol_array_Roquet_SpV(T, S, pressure, spv, start, npts, spv_ref=1.0/rho_ref) + do j=start,start+npts-1 + rho(j) = -rho_ref**2*spv(j) / (rho_ref*spv(j) + 1.0) ! In situ density [kg m-3] + enddo + else + call calculate_spec_vol_array_Roquet_SpV(T, S, pressure, spv, start, npts) + do j=start,start+npts-1 + rho(j) = 1.0 / spv(j) ! In situ density [kg m-3] + enddo + endif + +end subroutine calculate_density_array_Roquet_SpV + +!> Return the partial derivatives of specific volume with temperature and salinity for 1-d array +!! inputs and outputs, using the specific volume polynomial fit from Roquet et al. (2015). +subroutine calculate_specvol_derivs_Roquet_SpV(T, S, pressure, dSV_dT, dSV_dS, start, npts) + real, intent(in), dimension(:) :: T !< Conservative temperature [degC] + real, intent(in), dimension(:) :: S !< Absolute salinity [g kg-1] + real, intent(in), dimension(:) :: pressure !< Pressure [Pa] + real, intent(inout), dimension(:) :: dSV_dT !< The partial derivative of specific volume with + !! conservative temperature [m3 kg-1 degC-1] + real, intent(inout), dimension(:) :: dSV_dS !< The partial derivative of specific volume with + !! absolute salinity [m3 kg-1 ppt-1] + integer, intent(in) :: start !< The starting index for calculations + integer, intent(in) :: npts !< The number of values to calculate + + real :: zp ! Pressure, first in [dbar], then normalized by an assumed pressure range [nondim] + real :: zt ! Conservative temperature, first in [degC], then normalized by an assumed temperature range [nondim] + real :: zs ! Absolute salinity, first in [g kg-1], then the square root of salinity with an offset normalized + ! by an assumed salinity range [nondim] + real :: dSVdzt0 ! A contribution to the partial derivative of specific volume with temperature [m3 kg-1 degC-1] + ! from temperature anomalies at the surface pressure + real :: dSVdzt1 ! A contribution to the partial derivative of specific volume with temperature [m3 kg-1 degC-1] + ! that is proportional to pressure + real :: dSVdzt2 ! A contribution to the partial derivative of specific volume with temperature [m3 kg-1 degC-1] + ! that is proportional to pressure^2 + real :: dSVdzt3 ! A contribution to the partial derivative of specific volume with temperature [m3 kg-1 degC-1] + ! that is proportional to pressure^3 + real :: dSVdzs0 ! A contribution to the partial derivative of specific volume with + ! salinity [m3 kg-1 ppt-1] from temperature anomalies at the surface pressure + real :: dSVdzs1 ! A contribution to the partial derivative of specific volume with + ! salinity [m3 kg-1 ppt-1] proportional to pressure + real :: dSVdzs2 ! A contribution to the partial derivative of specific volume with + ! salinity [m3 kg-1 ppt-1] proportional to pressure^2 + real :: dSVdzs3 ! A contribution to the partial derivative of specific volume with + ! salinity [m3 kg-1 ppt-1] proportional to pressure^3 + integer :: j + + do j=start,start+npts-1 + ! Conversions to the units used here. + zt = T(j) * r1_T0 ! Conservative temperature normalized by a plausible oceanic range [nondim] + zs = SQRT( ABS( S(j) + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim] + zp = pressure(j) * (Pa2db*r1_P0) ! Convert pressure from Pascals to kilobars to normalize it [nondim] + + ! The next two lines should be used if it is necessary to convert potential temperature and + ! practical salinity to conservative temperature and absolute salinity. + ! zt = r1_T0 * gsw_ct_from_pt(S(j),T(j)) ! Convert potential temp to conservative temp [degC] + ! zs = SQRT( ABS( gsw_sr_from_sp(S(j)) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity. + + ! Find the partial derivative of specific volume with temperature + dSVdzt3 = ALP003 + dSVdzt2 = ALP002 + (zs*ALP102 + zt*ALP012) + dSVdzt1 = ALP001 + (zs*(ALP101 + zs*(ALP201 + zs*ALP301)) & + + zt*(ALP011 + (zs*(ALP111 + zs*ALP211) & + + zt*(ALP021 + (zs*ALP121 + zt*ALP031)) )) ) + dSVdzt0 = ALP000 + (zs*(ALP100 + zs*(ALP200 + zs*(ALP300 + zs*(ALP400 + zs*ALP500)))) & + + zt*(ALP010 + (zs*(ALP110 + zs*(ALP210 + zs*(ALP310 + zs*ALP410))) & + + zt*(ALP020 + (zs*(ALP120 + zs*(ALP220 + zs*ALP320)) & + + zt*(ALP030 + (zt*(ALP040 + (zs*ALP140 + zt*ALP050)) & + + zs*(ALP130 + zs*ALP230) )) )) )) ) + + dSV_dT(j) = dSVdzt0 + zp*(dSVdzt1 + zp*(dSVdzt2 + zp*dSVdzt3)) + + ! Find the partial derivative of specific volume with salinity + dSVdzs3 = BET003 + dSVdzs2 = BET002 + (zs*BET102 + zt*BET012) + dSVdzs1 = BET001 + (zs*(BET101 + zs*(BET201 + zs*BET301)) & + + zt*(BET011 + (zs*(BET111 + zs*BET211) & + + zt*(BET021 + (zs*BET121 + zt*BET031)) )) ) + dSVdzs0 = BET000 + (zs*(BET100 + zs*(BET200 + zs*(BET300 + zs*(BET400 + zs*BET500)))) & + + zt*(BET010 + (zs*(BET110 + zs*(BET210 + zs*(BET310 + zs*BET410))) & + + zt*(BET020 + (zs*(BET120 + zs*(BET220 + zs*BET320)) & + + zt*(BET030 + (zt*(BET040 + (zs*BET140 + zt*BET050)) & + + zs*(BET130 + zs*BET230) )) )) )) ) + + ! The division by zs here is because zs = sqrt(S + S0), so dSV_dS = dzs_dS * dSV_dzs = (0.5 / zs) * dSV_dzs + dSV_dS(j) = (dSVdzs0 + zp*(dSVdzs1 + zp*(dSVdzs2 + zp * dSVdzs3))) / zs + enddo + +end subroutine calculate_specvol_derivs_Roquet_SpV + + +!> Compute an array of derivatives of densities of sea water with temperature (drho_dT in [kg m-3 degC-1]) +!! and salinity (drho_dS in [kg m-3 ppt-1]) from absolute salinity (S [g kg-1]), conservative temperature +!! (T [degC]) and pressure [Pa], using the specific volume polynomial fit from Roquet et al. (2015). +subroutine calculate_density_derivs_array_Roquet_SpV(T, S, pressure, drho_dT, drho_dS, start, npts) + real, intent(in), dimension(:) :: T !< Conservative temperature [degC] + real, intent(in), dimension(:) :: S !< Absolute salinity [g kg-1] + real, intent(in), dimension(:) :: pressure !< pressure [Pa] + real, intent(out), dimension(:) :: drho_dT !< The partial derivative of density with + !! conservative temperature [kg m-3 degC-1] + real, intent(out), dimension(:) :: drho_dS !< The partial derivative of density with + !! absolute salinity [kg m-3 ppt-1] + integer, intent(in) :: start !< The starting index for calculations + integer, intent(in) :: npts !< The number of values to calculate + + ! Local variables + real, dimension(size(T)) :: specvol ! The specific volume [m3 kg-1] + real, dimension(size(T)) :: dSV_dT ! The partial derivative of specific volume with + ! conservative temperature [m3 kg-1 degC-1] + real, dimension(size(T)) :: dSV_dS ! The partial derivative of specific volume with + ! absolute salinity [m3 kg-1 ppt-1] + real :: rho ! The in situ density [kg m-3] + integer :: j + + call calculate_spec_vol_array_Roquet_SpV(T, S, pressure, specvol, start, npts) + call calculate_specvol_derivs_Roquet_SpV(T, S, pressure, dSV_dT, dSV_dS, start, npts) + + do j=start,start+npts-1 + rho = 1.0 / specvol(j) + drho_dT(j) = -dSv_dT(j) * rho**2 + drho_dS(j) = -dSv_dS(j) * rho**2 + enddo + +end subroutine calculate_density_derivs_array_Roquet_SpV + +!> Wrapper to calculate_density_derivs_array_Roquet_SpV for scalar inputs +subroutine calculate_density_derivs_scalar_Roquet_SpV(T, S, pressure, drho_dt, drho_ds) + real, intent(in) :: T !< Conservative temperature [degC] + real, intent(in) :: S !< Absolute salinity [g kg-1] + real, intent(in) :: pressure !< Pressure [Pa] + real, intent(out) :: drho_dT !< The partial derivative of density with + !! conservative temperature [kg m-3 degC-1] + real, intent(out) :: drho_dS !< The partial derivative of density with + !! absolute salinity [kg m-3 ppt-1] + ! Local variables + real, dimension(1) :: T0 ! A 1-d array with a copy of the conservative temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the absolute salinity [g kg-1] + real, dimension(1) :: pressure0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: drdt0 ! A 1-d array with a copy of the derivative of density + ! with conservative temperature [kg m-3 degC-1] + real, dimension(1) :: drds0 ! A 1-d array with a copy of the derivative of density + ! with absolute salinity [kg m-3 ppt-1] + + T0(1) = T + S0(1) = S + pressure0(1) = pressure + + call calculate_density_derivs_array_Roquet_SpV(T0, S0, pressure0, drdt0, drds0, 1, 1) + drho_dt = drdt0(1) + drho_ds = drds0(1) +end subroutine calculate_density_derivs_scalar_Roquet_SpV + +!> Compute the in situ density of sea water (rho in [kg m-3]) and the compressibility +!! (drho/dp = C_sound^-2, stored as drho_dp [s2 m-2]) from absolute salinity (sal [g kg-1]), +!! conservative temperature (T [degC]), and pressure [Pa], using the specific volume +!! polynomial fit from Roquet et al. (2015). +subroutine calculate_compress_Roquet_SpV(T, S, pressure, rho, drho_dp, start, npts) + real, intent(in), dimension(:) :: T !< Conservative temperature [degC] + real, intent(in), dimension(:) :: S !< Absolute salinity [g kg-1] + real, intent(in), dimension(:) :: pressure !< pressure [Pa] + real, intent(out), dimension(:) :: rho !< In situ density [kg m-3] + real, intent(out), dimension(:) :: drho_dp !< The partial derivative of density with pressure + !! (also the inverse of the square of sound speed) + !! [s2 m-2] + integer, intent(in) :: start !< The starting index for calculations + integer, intent(in) :: npts !< The number of values to calculate. + + ! Local variables + real :: zp ! Pressure normalized by an assumed pressure range [nondim] + real :: zt ! Conservative temperature normalized by an assumed temperature range [nondim] + real :: zs ! The square root of absolute salinity with an offset normalized + ! by an assumed salinity range [nondim] + real :: dSV_00p_dp ! Derivative of the pressure-dependent reference specific volume profile with + ! normalized pressure [m3 kg-1] + real :: dSV_TS_dp ! Derivative of the specific volume anomaly from the reference profile with + ! normalized pressure [m3 kg-1] + real :: SV_00p ! A pressure-dependent but temperature and salinity independent contribution to + ! specific volume at the reference temperature and salinity [m3 kg-1] + real :: SV_TS ! specific volume without a pressure-dependent contribution [m3 kg-1] + real :: SV_TS0 ! A contribution to specific volume from temperature and salinity anomalies at + ! the surface pressure [m3 kg-1] + real :: SV_TS1 ! A temperature and salinity dependent specific volume contribution that is + ! proportional to pressure [m3 kg-1] + real :: SV_TS2 ! A temperature and salinity dependent specific volume contribution that is + ! proportional to pressure**2 [m3 kg-1] + real :: SV_TS3 ! A temperature and salinity dependent specific volume contribution that is + ! proportional to pressure**3 [m3 kg-1] + real :: SV_0S0 ! Salinity dependent specific volume at the surface pressure and zero temperature [m3 kg-1] + real :: dSpecVol_dp ! The partial derivative of specific volume with pressure [m3 kg-1 Pa-1] + integer :: j + + ! The following algorithm was published by Roquet et al. (2015), intended for use + ! with NEMO, but it is not necessarily the algorithm used in NEMO ocean model. + do j=start,start+npts-1 + ! Conversions to the units used here. + zt = T(j) * r1_T0 ! Conservative temperature normalized by a plausible oceanic range [nondim] + zs = SQRT( ABS( S(j) + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim] + zp = pressure(j) * (Pa2db*r1_P0) ! Convert pressure from Pascals to kilobars to normalize it [nondim] + + ! The next two lines should be used if it is necessary to convert potential temperature and + ! practical salinity to conservative temperature and absolute salinity. + ! zt = r1_T0 * gsw_ct_from_pt(S(j),T(j)) ! Convert potential temp to conservative temp [degC] + ! zs = SQRT( ABS( gsw_sr_from_sp(S(j)) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity. + + SV_TS3 = SPV003 + (zs*SPV103 + zt*SPV013) + SV_TS2 = SPV002 + (zs*(SPV102 + zs*SPV202) & + + zt*(SPV012 + (zs*SPV112 + zt*SPV022)) ) + SV_TS1 = SPV001 + (zs*(SPV101 + zs*(SPV201 + zs*(SPV301 + zs*SPV401))) & + + zt*(SPV011 + (zs*(SPV111 + zs*(SPV211 + zs*SPV311)) & + + zt*(SPV021 + (zs*(SPV121 + zs*SPV221) & + + zt*(SPV031 + (zs*SPV131 + zt*SPV041)) )) )) ) + + SV_TS0 = zt*(SPV010 & + + (zs*(SPV110 + zs*(SPV210 + zs*(SPV310 + zs*(SPV410 + zs*SPV510)))) & + + zt*(SPV020 + (zs*(SPV120 + zs*(SPV220 + zs*(SPV320 + zs*SPV420))) & + + zt*(SPV030 + (zs*(SPV130 + zs*(SPV230 + zs*SPV330)) & + + zt*(SPV040 + (zs*(SPV140 + zs*SPV240) & + + zt*(SPV050 + (zs*SPV150 + zt*SPV060)) )) )) )) ) ) + + SV_0S0 = SPV000 + zs*(SPV100 + zs*(SPV200 + zs*(SPV300 + zs*(SPV400 + zs*(SPV500 + zs*SPV600))))) + + SV_00p = zp*(V00 + zp*(V01 + zp*(V02 + zp*(V03 + zp*(V04 + zp*V05))))) + + SV_TS = (SV_TS0 + SV_0S0) + zp*(SV_TS1 + zp*(SV_TS2 + zp*SV_TS3)) + ! specvol = SV_TS + SV_00p ! In situ specific volume [m3 kg-1] + rho(j) = 1.0 / (SV_TS + SV_00p) ! In situ density [kg m-3] + + dSV_00p_dp = V00 + zp*(2.*V01 + zp*(3.*V02 + zp*(4.*V03 + zp*(5.*V04 + zp*(6.*V05))))) + dSV_TS_dp = SV_TS1 + zp*(2.*SV_TS2 + zp*(3.*SV_TS3)) + dSpecVol_dp = (dSV_TS_dp + dSV_00p_dp) * (Pa2db*r1_P0) ! [m3 kg-1 Pa-1] + drho_dp(j) = -dSpecVol_dp * rho(j)**2 ! Compressibility [s2 m-2] + + enddo +end subroutine calculate_compress_Roquet_SpV + + +!> Second derivatives of specific volume with respect to temperature, salinity, and pressure for a +!! 1-d array inputs and outputs using the specific volume polynomial fit from Roquet et al. (2015). +subroutine calc_spec_vol_second_derivs_array_Roquet_SpV(T, S, P, dSV_ds_ds, dSV_ds_dt, dSV_dt_dt, & + dSV_ds_dp, dSV_dt_dp, start, npts) + real, dimension(:), intent(in ) :: T !< Conservative temperature [degC] + real, dimension(:), intent(in ) :: S !< Absolute salinity [g kg-1] + real, dimension(:), intent(in ) :: P !< Pressure [Pa] + real, dimension(:), intent(inout) :: dSV_ds_ds !< Second derivative of specific volume with respect + !! to salinity [m3 kg-1 ppt-2] + real, dimension(:), intent(inout) :: dSV_ds_dt !< Second derivative of specific volume with respect + !! to salinity and temperature [m3 kg-1 ppt-1 degC-1] + real, dimension(:), intent(inout) :: dSV_dt_dt !< Second derivative of specific volume with respect + !! to temperature [m3 kg-1 degC-2] + real, dimension(:), intent(inout) :: dSV_ds_dp !< Second derivative of specific volume with respect to pressure + !! and salinity [m3 kg-1 ppt-1 Pa-1] + real, dimension(:), intent(inout) :: dSV_dt_dp !< Second derivative of specific volume with respect to pressure + !! and temperature [m3 kg-1 degC-1 Pa-1] + integer, intent(in ) :: start !< Starting index in T,S,P + integer, intent(in ) :: npts !< Number of points to loop over + + ! Local variables + real :: zp ! Pressure normalized by an assumed pressure range [nondim] + real :: zt ! Conservative temperature normalized by an assumed temperature range [nondim] + real :: zs ! The square root of absolute salinity with an offset normalized + ! by an assumed salinity range [nondim] + real :: I_s ! The inverse of zs [nondim] + real :: d2SV_p0 ! A contribution to one of the second derivatives that is independent of pressure [various] + real :: d2SV_p1 ! A contribution to one of the second derivatives that is proportional to pressure [various] + real :: d2SV_p2 ! A contribution to one of the second derivatives that is proportional to pressure^2 [various] + real :: d2SV_p3 ! A contribution to one of the second derivatives that is proportional to pressure^3 [various] + integer :: j + + do j = start,start+npts-1 + ! Conversions to the units used here. + zt = T(j) * r1_T0 ! Conservative temperature normalized by a plausible oceanic range [nondim] + zs = SQRT( ABS( S(j) + rdeltaS ) * r1_S0 ) ! square root of normalized salinity plus an offset [nondim] + zp = P(j) * (Pa2db*r1_P0) ! Convert pressure from Pascals to kilobars to normalize it [nondim] + + ! The next two lines should be used if it is necessary to convert potential temperature and + ! practical salinity to conservative temperature and absolute salinity. + ! zt = r1_T0 * gsw_ct_from_pt(S(j),T(j)) ! Convert potential temp to conservative temp [degC] + ! zs = SQRT( ABS( gsw_sr_from_sp(S(j)) + rdeltaS ) * r1_S0 ) ! Convert S from practical to absolute salinity. + + I_s = 1.0 / zs + + ! Find dSV_ds_ds + d2SV_p3 = -SPV103*I_s**2 + d2SV_p2 = -(SPV102 + zt*SPV112)*I_s**2 + d2SV_p1 = (3.*SPV301 + (zt*(3.*SPV311) + zs*(8.*SPV401))) & + - ( SPV101 + zt*(SPV111 + zt*(SPV121 + zt*SPV131)) )*I_s**2 + d2SV_p0 = (3.*SPV300 + (zs*(8.*SPV400 + zs*(15.*SPV500 + zs*(24.*SPV600))) & + + zt*(3.*SPV310 + (zs*(8.*SPV410 + zs*(15.*SPV510)) & + + zt*(3.*SPV320 + (zs*(8.*SPV420) + zt*(3.*SPV330))) )) )) & + - (SPV100 + zt*(SPV110 + zt*(SPV120 + zt*(SPV130 + zt*(SPV140 + zt*SPV150)))) )*I_s**2 + dSV_dS_dS(j) = (0.5*r1_S0)**2 * ((d2SV_p0 + zp*(d2SV_p1 + zp*(d2SV_p2 + zp*d2SV_p3))) * I_s) + + ! Find dSV_ds_dt + d2SV_p2 = SPV112 + d2SV_p1 = SPV111 + (zs*(2.*SPV211 + zs*(3.*SPV311)) & + + zt*(2.*SPV121 + (zs*(4.*SPV221) + zt*(3.*SPV131))) ) + d2SV_p0 = SPV110 + (zs*(2.*SPV210 + zs*(3.*SPV310 + zs*(4.*SPV410 + zs*(5.*SPV510)))) & + + zt*(2.*SPV120 + (zs*(4.*SPV220 + zs*(6.*SPV320 + zs*(8.*SPV420))) & + + zt*(3.*SPV130 + (zs*(6.*SPV230 + zs*(9.*SPV330)) & + + zt*(4.*SPV140 + (zs*(8.*SPV240) & + + zt*(5.*SPV150))) )) )) ) + dSV_ds_dt(j) = (0.5*r1_S0*r1_T0) * ((d2SV_p0 + zp*(d2SV_p1 + zp*d2SV_p2)) * I_s) + + ! Find dSV_dt_dt + d2SV_p2 = 2.*SPV022 + d2SV_p1 = 2.*SPV021 + (zs*(2.*SPV121 + zs*(2.*SPV221)) & + + zt*(6.*SPV031 + (zs*(6.*SPV131) + zt*(12.*SPV041))) ) + d2SV_p0 = 2.*SPV020 + (zs*(2.*SPV120 + zs*( 2.*SPV220 + zs*( 2.*SPV320 + zs * (2.*SPV420)))) & + + zt*(6.*SPV030 + (zs*( 6.*SPV130 + zs*( 6.*SPV230 + zs * (6.*SPV330))) & + + zt*(12.*SPV040 + (zs*(12.*SPV140 + zs *(12.*SPV240)) & + + zt*(20.*SPV050 + (zs*(20.*SPV150) & + + zt*(30.*SPV060) )) )) )) ) + dSV_dt_dt(j) = (d2SV_p0 + zp*(d2SV_p1 + zp*d2SV_p2)) * r1_T0**2 + + ! Find dSV_ds_dp + d2SV_p2 = 3.*SPV103 + d2SV_p1 = 2.*SPV102 + (zs*(4.*SPV202) + zt*(2.*SPV112)) + d2SV_p0 = SPV101 + (zs*(2.*SPV201 + zs*(3.*SPV301 + zs*(4.*SPV401))) & + + zt*(SPV111 + (zs*(2.*SPV211 + zs*(3.*SPV311)) & + + zt*( SPV121 + (zs*(2.*SPV221) + zt*SPV131)) )) ) + dSV_ds_dp(j) = ((d2SV_p0 + zp*(d2SV_p1 + zp*d2SV_p2)) * I_s) * (0.5*r1_S0 * Pa2db*r1_P0) + + ! Find dSV_dt_dp + d2SV_p2 = 3.*SPV013 + d2SV_p1 = 2.*SPV012 + (zs*(2.*SPV112) + zt*(4.*SPV022)) + d2SV_p0 = SPV011 + (zs*(SPV111 + zs*( SPV211 + zs* SPV311)) & + + zt*(2.*SPV021 + (zs*(2.*SPV121 + zs*(2.*SPV221)) & + + zt*(3.*SPV031 + (zs*(3.*SPV131) + zt*(4.*SPV041))) )) ) + dSV_dt_dp(j) = (d2SV_p0 + zp*(d2SV_p1 + zp*d2SV_p2)) * (Pa2db*r1_P0* r1_T0) + enddo + +end subroutine calc_spec_vol_second_derivs_array_Roquet_SpV + + +!> Second derivatives of density with respect to temperature, salinity, and pressure for a +!! 1-d array inputs and outputs using the specific volume polynomial fit from Roquet et al. (2015). +subroutine calculate_density_second_derivs_array_Roquet_SpV(T, S, P, drho_ds_ds, drho_ds_dt, drho_dt_dt, & + drho_ds_dp, drho_dt_dp, start, npts) + real, dimension(:), intent(in ) :: T !< Conservative temperature [degC] + real, dimension(:), intent(in ) :: S !< Absolute salinity [g kg-1] + real, dimension(:), intent(in ) :: P !< Pressure [Pa] + real, dimension(:), intent(inout) :: drho_ds_ds !< Second derivative of density with respect + !! to salinity [kg m-3 ppt-2] + real, dimension(:), intent(inout) :: drho_ds_dt !< Second derivative of density with respect + !! to salinity and temperature [kg m-3 ppt-1 degC-1] + real, dimension(:), intent(inout) :: drho_dt_dt !< Second derivative of density with respect + !! to temperature [kg m-3 degC-2] + real, dimension(:), intent(inout) :: drho_ds_dp !< Second derivative of density with respect to pressure + !! and salinity [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1] + real, dimension(:), intent(inout) :: drho_dt_dp !< Second derivative of density with respect to pressure + !! and temperature [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] + integer, intent(in ) :: start !< Starting index in T,S,P + integer, intent(in ) :: npts !< Number of points to loop over + + ! Local variables + real, dimension(size(T)) :: rho ! The in situ density [kg m-3] + real, dimension(size(T)) :: drho_dp ! The partial derivative of density with pressure + ! (also the inverse of the square of sound speed) [s2 m-2] + real, dimension(size(T)) :: dSV_dT ! The partial derivative of specific volume with + ! conservative temperature [m3 kg-1 degC-1] + real, dimension(size(T)) :: dSV_dS ! The partial derivative of specific volume with + ! absolute salinity [m3 kg-1 ppt-1] + real, dimension(size(T)) :: dSV_ds_ds ! Second derivative of specific volume with respect + ! to salinity [m3 kg-1 ppt-2] + real, dimension(size(T)) :: dSV_ds_dt ! Second derivative of specific volume with respect + ! to salinity and temperature [m3 kg-1 ppt-1 degC-1] + real, dimension(size(T)) :: dSV_dt_dt ! Second derivative of specific volume with respect + ! to temperature [m3 kg-1 degC-2] + real, dimension(size(T)) :: dSV_ds_dp ! Second derivative of specific volume with respect to pressure + ! and salinity [m3 kg-1 ppt-1 Pa-1] + real, dimension(size(T)) :: dSV_dt_dp ! Second derivative of specific volume with respect to pressure + ! and temperature [m3 kg-1 degC-1 Pa-1] + integer :: j + + call calc_spec_vol_second_derivs_array_Roquet_SpV(T, S, P, dSV_ds_ds, dSV_ds_dt, dSV_dt_dt, & + dSV_ds_dp, dSV_dt_dp, start, npts) + call calculate_specvol_derivs_Roquet_SpV(T, S, P, dSV_dT, dSV_dS, start, npts) + call calculate_compress_Roquet_SpV(T, S, P, rho, drho_dp, start, npts) + + do j = start,start+npts-1 + ! Find drho_ds_ds + drho_dS_dS(j) = rho(j)**2 * (2.0*rho(j)*dSV_dS(j)**2 - dSV_dS_dS(j)) + + ! Find drho_ds_dt + drho_ds_dt(j) = rho(j)**2 * (2.0*rho(j)*(dSV_dT(j)*dSV_dS(j)) - dSV_dS_dT(j)) + + ! Find drho_dt_dt + drho_dT_dT(j) = rho(j)**2 * (2.0*rho(j)*dSV_dT(j)**2 - dSV_dT_dT(j)) + + ! Find drho_ds_dp + drho_ds_dp(j) = -rho(j) * (2.0*dSV_dS(j) * drho_dp(j) + rho(j) * dSV_dS_dp(j)) + + ! Find drho_dt_dp + drho_dt_dp(j) = -rho(j) * (2.0*dSV_dT(j) * drho_dp(j) + rho(j) * dSV_dT_dp(j)) + enddo + +end subroutine calculate_density_second_derivs_array_Roquet_SpV + +!> Second derivatives of density with respect to temperature, salinity, and pressure for scalar inputs. +!! +!! The scalar version of calculate_density_second_derivs promotes scalar inputs to 1-element array +!! and then demotes the output back to a scalar +subroutine calculate_density_second_derivs_scalar_Roquet_SpV(T, S, P, drho_ds_ds, drho_ds_dt, drho_dt_dt, & + drho_ds_dp, drho_dt_dp) + real, intent(in ) :: T !< Conservative temperature [degC] + real, intent(in ) :: S !< Absolute salinity [g kg-1] + real, intent(in ) :: P !< pressure [Pa] + real, intent( out) :: drho_ds_ds !< Second derivative of density with respect + !! to salinity [kg m-3 ppt-2] + real, intent( out) :: drho_ds_dt !< Second derivative of density with respect + !! to salinity and temperature [kg m-3 ppt-1 degC-1] + real, intent( out) :: drho_dt_dt !< Second derivative of density with respect + !! to temperature [kg m-3 degC-2] + real, intent( out) :: drho_ds_dp !< Second derivative of density with respect to pressure + !! and salinity [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1] + real, intent( out) :: drho_dt_dp !< Second derivative of density with respect to pressure + !! and temperature [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] + ! Local variables + real, dimension(1) :: T0 ! A 1-d array with a copy of the temperature [degC] + real, dimension(1) :: S0 ! A 1-d array with a copy of the salinity [g kg-1] + real, dimension(1) :: p0 ! A 1-d array with a copy of the pressure [Pa] + real, dimension(1) :: drdsds ! The second derivative of density with salinity [kg m-3 ppt-2] + real, dimension(1) :: drdsdt ! The second derivative of density with salinity and + ! temperature [kg m-3 ppt-1 degC-1] + real, dimension(1) :: drdtdt ! The second derivative of density with temperature [kg m-3 degC-2] + real, dimension(1) :: drdsdp ! The second derivative of density with salinity and + ! pressure [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1] + real, dimension(1) :: drdtdp ! The second derivative of density with temperature and + ! pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] + + T0(1) = T + S0(1) = S + P0(1) = P + call calculate_density_second_derivs_array_Roquet_SpV(T0, S0, P0, drdsds, drdsdt, drdtdt, drdsdp, drdtdp, 1, 1) + drho_ds_ds = drdsds(1) + drho_ds_dt = drdsdt(1) + drho_dt_dt = drdtdt(1) + drho_ds_dp = drdsdp(1) + drho_dt_dp = drdtdp(1) + +end subroutine calculate_density_second_derivs_scalar_Roquet_SpV + +!> \namespace mom_eos_Roquet_SpV +!! +!! \section section_EOS_Roquet_SpV NEMO equation of state +!! +!! Fabien Roquet and colleagues developed this equation of state using a simple polynomial fit +!! to the TEOS-10 equation of state expressions for specific, for efficiency when used with a +!! non-Boussinesq ocean model. This particular equation of state is a balance between an +!! accuracy that matches the TEOS-10 density to better than observational uncertainty with a +!! polynomial form that can be evaluated quickly despite having 55 terms. +!! +!! \subsection section_EOS_Roquet_Spv_references References +!! +!! 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. + +end module MOM_EOS_Roquet_Spv