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/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 345d0dd5eb..8e199cf2e6 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 @@ -240,7 +241,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, & @@ -495,7 +496,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_absS) call callTree_leave("update_ocean_model()") end subroutine update_ocean_model @@ -741,10 +742,11 @@ 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_absS, patm, press_to_z) 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_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 @@ -769,9 +771,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_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 + 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 @@ -825,7 +841,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_absS) end subroutine ocean_model_init_sfc ! 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/core/MOM.F90 b/src/core/MOM.F90 index 3a6203dccb..d5c20f3f96 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_EOS, only : calculate_density use MOM_error_checking, only : check_redundant use MOM_grid, only : ocean_grid_type, set_first_direction @@ -125,7 +126,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 - ! Offline modules use MOM_offline_transport, only : offline_transport_CS use MOM_offline_transport, only : transport_by_files, next_modulo_time @@ -192,6 +192,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_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. logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics. logical :: mixedlayer_restrat !< If true, use submesoscale mixed layer restratifying scheme. @@ -285,6 +289,8 @@ module MOM integer :: id_h = -1 integer :: id_T = -1 integer :: id_S = -1 + integer :: id_Tcon = -1 + integer :: id_Sabs = -1 ! 2-d surface and bottom fields integer :: id_zos = -1 @@ -302,6 +308,8 @@ module MOM integer :: id_ssh_inst = -1 integer :: id_tob = -1 integer :: id_sob = -1 + integer :: id_sstcon = -1 + integer :: id_sssabs = -1 ! heat and salt flux fields integer :: id_fraz = -1 @@ -487,6 +495,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 !TEOS10 Diagnostics real, dimension(SZIB_(CS%G), SZJ_(CS%G)) :: umo2d ! Diagnostics real, dimension(SZI_(CS%G), SZJB_(CS%G)) :: vmo2d ! Diagnostics real, dimension(SZIB_(CS%G), SZJ_(CS%G), SZK_(CS%G)) :: umo ! Diagnostics @@ -1212,11 +1221,29 @@ 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) + 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) + 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) + !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)) + 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) + 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) @@ -1406,24 +1433,38 @@ 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) + 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) + !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)) + 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) & @@ -2067,6 +2108,12 @@ 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_ABSSAL", CS%use_conT_absS, & + "If true, , the prognostics T&S are the conservative temperature \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.) 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"//& @@ -2916,6 +2963,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_Sabs = register_diag_field('ocean_model', 'abssalt', diag%axesTL, Time, & + long_name='Absolute Salinity', units='g/Kg') + CS%id_sstcon = register_diag_field('ocean_model', 'conSST', diag%axesT1, Time, & + 'Sea Surface Conservative Temperature', 'Celsius', CS%missing) + CS%id_sssabs = register_diag_field('ocean_model', 'absSSS', diag%axesT1, Time, & + 'Sea Surface Absolute Salinity', 'g/Kg', CS%missing) + endif if (CS%use_temperature .and. CS%use_frazil) then diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 980074db69..54478a4dea 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -3,9 +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_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 @@ -26,6 +39,8 @@ 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 +public gsw_sp_from_sr, gsw_pt_from_ct !> Calculates density of sea water from T, S and P interface calculate_density @@ -61,10 +76,14 @@ module MOM_EOS integer, parameter :: EOS_LINEAR = 1 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 @@ -94,6 +113,10 @@ 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 (EOS_NEMO) + call calculate_density_scalar_nemo(T, S, pressure, rho) case default call MOM_error(FATAL, & "calculate_density_scalar: EOS is not valid.") @@ -122,10 +145,15 @@ 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 (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. @@ -186,7 +214,7 @@ 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 - + !! if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") @@ -198,6 +226,10 @@ 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 (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.") @@ -235,6 +267,15 @@ 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 (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.") @@ -265,6 +306,10 @@ 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 (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.") @@ -295,6 +340,10 @@ 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 (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.") @@ -451,7 +500,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", "NEMO" 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 +509,10 @@ 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 (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.") @@ -2091,6 +2144,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) .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 + 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_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 new file mode 100644 index 0000000000..ff7fbd8177 --- /dev/null +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -0,0 +1,445 @@ +module MOM_EOS_NEMO +!*********************************************************************** +!* GNU General Public License * +!* This file is a part of MOM. * +!* * +!* MOM is free software; you can redistribute it and/or modify it and * +!* are expected to follow the terms of the GNU General Public License * +!* as published by the Free Software Foundation; either version 2 of * +!* the License, or (at your option) any later version. * +!* * +!* MOM is distributed in the hope that it will be useful, but WITHOUT * +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * +!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * +!* License for more details. * +!* * +!* For the full text of the GNU General Public License, * +!* write to: Free Software Foundation, Inc., * +!* 675 Mass Ave, Cambridge, MA 02139, USA. * +!* or see: http://www.gnu.org/licenses/gpl.html * +!*********************************************************************** + +!*********************************************************************** +!* The subroutines in this file implement the equation of state for * +!* sea water using the formulae provided by NEMO developer Roquet * +!* in a private communication , Roquet et al, Ocean Modelling (2015) * +!* Roquet, F., Madec, G., McDougall, T. J., and Barker, P. M., 2015. * +!* Accurate polynomial expressions for the density and specific volume* +!* of seawater using the TEOS-10 standard. Ocean Modelling, 90:29-43. * +!* These algorithms are NOT from NEMO package!! * +!*********************************************************************** + +!use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt +use gsw_mod_toolbox, only : gsw_rho_first_derivatives + + +implicit none ; private + +public calculate_compress_nemo, calculate_density_nemo +public calculate_density_derivs_nemo, 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 + !Conversions + zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + + !The following algorithm was provided by Roquet in a private communication. + !It is not necessarily the algorithm used in NEMO ocean! + zp = zp * r1_P0 !pressure + zt = zt * r1_T0 !temperature + zs = SQRT( ABS( zs + rdeltaS ) * r1_S0 ) ! square root salinity + ! + zn3 = EOS013*zt & + & + EOS103*zs+EOS003 + ! + zn2 = (EOS022*zt & + & + EOS112*zs+EOS012)*zt & + & + (EOS202*zs+EOS102)*zs+EOS002 + ! + zn1 = (((EOS041*zt & + & + EOS131*zs+EOS031)*zt & + & + (EOS221*zs+EOS121)*zs+EOS021)*zt & + & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & + & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 + ! + zn0 = (((((EOS060*zt & + & + EOS150*zs+EOS050)*zt & + & + (EOS240*zs+EOS140)*zs+EOS040)*zt & + & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & + & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & + & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & + & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 + ! + zn = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0 + ! + zr0 = (((((R05 * zp+R04) * zp+R03 ) * zp+R02 ) * zp+R01) * zp+R00) * zp + ! + rho(j) = ( zn + zr0 ) ! density + + enddo +end subroutine calculate_density_array_nemo + +subroutine calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) + real, intent(in), dimension(:) :: T, 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 + !Conversions + zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + + !The following algorithm was provided by Roquet in a private communication. + !It is not necessarily the algorithm used in NEMO ocean! + zp = zp * r1_P0 ! pressure (first converted to decibar) + zt = zt * r1_T0 ! temperature + zs = SQRT( ABS( zs + rdeltaS ) * r1_S0 ) ! square root salinity + ! + ! alpha + zn3 = ALP003 + ! + zn2 = ALP012*zt + ALP102*zs+ALP002 + ! + zn1 = ((ALP031*zt & + & + ALP121*zs+ALP021)*zt & + & + (ALP211*zs+ALP111)*zs+ALP011)*zt & + & + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 + ! + zn0 = ((((ALP050*zt & + & + ALP140*zs+ALP040)*zt & + & + (ALP230*zs+ALP130)*zs+ALP030)*zt & + & + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt & + & + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt & + & + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 + ! + zn = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0 + ! + drho_dT(j) = -zn + ! + ! beta + ! + zn3 = BET003 + ! + zn2 = BET012*zt + BET102*zs+BET002 + ! + zn1 = ((BET031*zt & + & + BET121*zs+BET021)*zt & + & + (BET211*zs+BET111)*zs+BET011)*zt & + & + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 + ! + zn0 = ((((BET050*zt & + & + BET140*zs+BET040)*zt & + & + (BET230*zs+BET130)*zs+BET030)*zt & + & + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt & + & + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt & + & + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 + ! + zn = ( ( zn3 * zp + zn2 ) * zp + zn1 ) * zp + zn0 + ! + drho_dS(j) = zn / zs + enddo + +end subroutine calculate_density_derivs_nemo + +subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) + real, intent(in), dimension(:) :: T, 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 :: zs,zt,zp + integer :: j + + call calculate_density_array_nemo(T, S, pressure, rho, start, npts) + ! + !NOTE: The following calculates the TEOS10 approximation to compressibility + ! since the corresponding NEMO approximation is not available yet. + ! + do j=start,start+npts-1 + !Conversions + zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp(j)) + enddo +end subroutine calculate_compress_nemo + +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 + + 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 = 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 + + !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) + 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 + ! + 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 new file mode 100644 index 0000000000..38bb0a0df8 --- /dev/null +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -0,0 +1,229 @@ + +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 TEOS10 functions * +!*********************************************************************** + +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 + +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 +public gsw_sp_from_sr, gsw_pt_from_ct + +interface calculate_density_teos10 + module procedure calculate_density_scalar_teos10, calculate_density_array_teos10 +end interface calculate_density_teos10 + + real, parameter :: Pa2db = 1.e-4 + +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 * +! * TEOS10 website. * +! *====================================================================* + + 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 :: zs,zt,zp + integer :: j + + do j=start,start+npts-1 + !Conversions + zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + + 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) + 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 :: zs,zt,zp + integer :: j + + do j=start,start+npts-1 + !Conversions + zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + 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)) + 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 :: zs, zt, zp + integer :: j + + do j=start,start+npts-1 + !Conversions + zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + 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)) + 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 :: zs,zt,zp + integer :: j + + do j=start,start+npts-1 + !Conversions + zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + 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) + call gsw_rho_first_derivatives(zs,zt,zp, 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 :: zp1, zp2, zt, zs + 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 = 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) + 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 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 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 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 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 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 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 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 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 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 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 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 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 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 + +!-------------------------------------------------------------------------- + + + diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 5b133d49fd..c6a7728110 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -38,7 +38,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_OBC_data @@ -1771,6 +1771,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 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 call calculate_density(temp_z(:,j,k),salt_z(:,j,k), press, rho_z(:,j,k), is, ie, eos)