Skip to content

Add flexibility for TFreeze form w/ different model definition of temp/salt #893

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marshallward merged 6 commits into from
May 8, 2025
Merged

Add flexibility for TFreeze form w/ different model definition of temp/salt #893

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
0c777ae
Change fatal error to warning with Roquet EOS and linear TFreeze.
Apr 30, 2025
996258b
Add flexibility for TFreeze and different definition of temp/salt fro…
May 1, 2025
55c6486
Remove spaces on blank line
May 5, 2025
9c28ecc
Initialize new EOS type logicals
May 5, 2025
4601982
Merge branch 'dev/gfdl' into TFreeze_error
marshallward May 8, 2025
ab4deb3
Merge branch 'dev/gfdl' into TFreeze_error
marshallward May 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3024,7 +3024,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
! initialization routine for tv.
if (use_EOS) then
allocate(CS%tv%eqn_of_state)
call EOS_init(param_file, CS%tv%eqn_of_state, US)
call EOS_init(param_file, CS%tv%eqn_of_state, US, use_conT_absS)
endif
if (use_temperature) then
allocate(CS%tv%TempxPmE(isd:ied,jsd:jed), source=0.0)
Expand Down
146 changes: 120 additions & 26 deletions src/equation_of_state/MOM_EOS.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module MOM_EOS
use MOM_EOS_Roquet_rho, only : Roquet_rho_EOS
use MOM_EOS_Roquet_SpV, only : Roquet_SpV_EOS
use MOM_EOS_TEOS10, only : TEOS10_EOS
use MOM_EOS_TEOS10, only : gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp
use MOM_EOS_TEOS10, only : gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp, gsw_ct_from_pt
use MOM_temperature_convert, only : poTemp_to_consTemp, consTemp_to_poTemp
use MOM_TFreeze, only : calculate_TFreeze_linear, calculate_TFreeze_Millero
use MOM_TFreeze, only : calculate_TFreeze_teos10, calculate_TFreeze_TEOS_poly
Expand Down Expand Up @@ -120,6 +120,13 @@ module MOM_EOS
real :: TFr_S0_P0 !< The freezing potential temperature at S=0, P=0 [degC]
real :: dTFr_dS !< The derivative of freezing point with salinity [degC ppt-1]
real :: dTFr_dp !< The derivative of freezing point with pressure [degC Pa-1]
! The following are logicals pertaining to definitions of the thermodynamic state variables
logical :: use_conT_absS =.false. !< True if the model internal temperature is the conservative temperature and
!! the salinity is absolute salinity. These could be separated into two flags,
!! but right now it is controlled by one input parameter and there is no known
!! need to have one True and one False.
logical :: TFreeze_S_is_pracS =.true. !< True if the freezing point expression is formulated from practical salinity
logical :: TFreeze_T_is_potT = .true. !< True if the freezing point expression yields a potential temperature

logical :: use_Wright_2nd_deriv_bug = .false. !< If true, use a separate subroutine that
!! retains a buggy version of the calculations of the second
Expand Down Expand Up @@ -529,28 +536,50 @@ subroutine calculate_TFreeze_scalar(S, pressure, T_fr, EOS, pres_scale, scale_fr
! Local variables
real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1]
real :: S_scale ! A factor to convert salinity to units of ppt [ppt S-1 ~> 1]
real :: iS_scale! A factor to convert salinity to units of S [S ppt-1 ~> 1]
real :: absS ! A salinity converted to absolute salinity, only used in specific scenarios [ppt]
real :: TFreeze_S ! The salinity for the freezing equation in model units [S ~> PSU or ppt]

p_scale = 1.0 ; S_scale = 1.0
p_scale = 1.0 ; S_scale = 1.0 ; iS_scale = 1.0
if (present(pres_scale)) p_scale = pres_scale
if (present(scale_from_EOS)) then ; if (scale_from_EOS) then
p_scale = EOS%RL2_T2_to_Pa
S_scale = EOS%S_to_ppt
iS_scale = EOS%ppt_to_S
endif ; endif

if (EOS%use_conT_absS) then
! Otherwise absS is unneeded and therefore unset
absS = S*S_scale
if (EOS%TFreeze_S_is_pracS) then
TFreeze_S = gsw_sp_from_sr(absS)*iS_scale
else
TFreeze_S = S
endif
else
TFreeze_S = S
endif

select case (EOS%form_of_TFreeze)
case (TFREEZE_LINEAR)
call calculate_TFreeze_linear(S_scale*S, p_scale*pressure, T_fr, EOS%TFr_S0_P0, &
call calculate_TFreeze_linear(S_scale*TFreeze_S, p_scale*pressure, T_fr, EOS%TFr_S0_P0, &
EOS%dTFr_dS, EOS%dTFr_dp)
case (TFREEZE_MILLERO)
call calculate_TFreeze_Millero(S_scale*S, p_scale*pressure, T_fr)
call calculate_TFreeze_Millero(S_scale*TFreeze_S, p_scale*pressure, T_fr)
case (TFREEZE_TEOSPOLY)
call calculate_TFreeze_TEOS_poly(S_scale*S, p_scale*pressure, T_fr)
call calculate_TFreeze_TEOS_poly(S_scale*TFreeze_S, p_scale*pressure, T_fr)
case (TFREEZE_TEOS10)
call calculate_TFreeze_teos10(S_scale*S, p_scale*pressure, T_fr)
call calculate_TFreeze_teos10(S_scale*TFreeze_S, p_scale*pressure, T_fr)
case default
call MOM_error(FATAL, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
end select

if (EOS%use_conT_absS .and. EOS%TFreeze_T_is_potT) then
! absS is set only if EOS%use_conT_absS is True
! absS and T_fr have physical units here and don't need converted
T_fr = gsw_ct_from_pt(absS,T_fr)
endif

if (present(scale_from_EOS)) then ; if (scale_from_EOS) then
T_fr = EOS%degC_to_C * T_fr
endif ; endif
Expand All @@ -561,8 +590,8 @@ end subroutine calculate_TFreeze_scalar
subroutine calculate_TFreeze_array(S, pressure, T_fr, start, npts, EOS, pres_scale)
real, dimension(:), intent(in) :: S !< Salinity [ppt]
real, dimension(:), intent(in) :: pressure !< Pressure, in [Pa] or [R L2 T-2 ~> Pa] depending on pres_scale
real, dimension(:), intent(inout) :: T_fr !< Freezing point potential temperature referenced
!! to the surface [degC]
real, dimension(:), intent(inout) :: T_fr !< Freezing point, either potential temperature referenced to the
!! surface or conservative temperature depending on settings [degC]
integer, intent(in) :: start !< Starting index within the array
integer, intent(in) :: npts !< The number of values to calculate
type(EOS_type), intent(in) :: EOS !< Equation of state structure
Expand All @@ -572,57 +601,80 @@ subroutine calculate_TFreeze_array(S, pressure, T_fr, start, npts, EOS, pres_sca
! Local variables
real, dimension(size(pressure)) :: pres ! Pressure converted to [Pa]
real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1]
real, dimension(size(S)) :: absS ! A salinity converted to absolute salinity, only used in specific scenarios [ppt]
real, dimension(size(S)) :: TFreeze_S ! The salinity for the freezing equation in model units [S ~> PSU or ppt]
integer :: j

p_scale = 1.0 ; if (present(pres_scale)) p_scale = pres_scale

if (EOS%use_conT_absS) then
! Otherwise absS is unneeded and therefore unset
absS(:) = S(:)
if (EOS%TFreeze_S_is_pracS) then
TFreeze_S(:) = gsw_sp_from_sr(absS(:))
else
TFreeze_S(:) = S(:)
endif
else
TFreeze_S(:) = S(:)
endif

if (p_scale == 1.0) then
select case (EOS%form_of_TFreeze)
case (TFREEZE_LINEAR)
call calculate_TFreeze_linear(S, pressure, T_fr, start, npts, &
call calculate_TFreeze_linear(TFreeze_S, pressure, T_fr, start, npts, &
EOS%TFr_S0_P0, EOS%dTFr_dS, EOS%dTFr_dp)
case (TFREEZE_MILLERO)
call calculate_TFreeze_Millero(S, pressure, T_fr, start, npts)
call calculate_TFreeze_Millero(TFreeze_S, pressure, T_fr, start, npts)
case (TFREEZE_TEOSPOLY)
call calculate_TFreeze_TEOS_poly(S, pressure, T_fr, start, npts)
call calculate_TFreeze_TEOS_poly(TFreeze_S, pressure, T_fr, start, npts)
case (TFREEZE_TEOS10)
call calculate_TFreeze_teos10(S, pressure, T_fr, start, npts)
call calculate_TFreeze_teos10(TFreeze_S, pressure, T_fr, start, npts)
case default
call MOM_error(FATAL, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
end select
else
do j=start,start+npts-1 ; pres(j) = p_scale * pressure(j) ; enddo
select case (EOS%form_of_TFreeze)
case (TFREEZE_LINEAR)
call calculate_TFreeze_linear(S, pres, T_fr, start, npts, &
call calculate_TFreeze_linear(TFreeze_S, pres, T_fr, start, npts, &
EOS%TFr_S0_P0, EOS%dTFr_dS, EOS%dTFr_dp)
case (TFREEZE_MILLERO)
call calculate_TFreeze_Millero(S, pres, T_fr, start, npts)
call calculate_TFreeze_Millero(TFreeze_S, pres, T_fr, start, npts)
case (TFREEZE_TEOS10)
call calculate_TFreeze_teos10(S, pres, T_fr, start, npts)
call calculate_TFreeze_teos10(TFreeze_S, pres, T_fr, start, npts)
case (TFREEZE_TEOSPOLY)
call calculate_TFreeze_TEOS_poly(S, pres, T_fr, start, npts)
call calculate_TFreeze_TEOS_poly(TFreeze_S, pres, T_fr, start, npts)
case default
call MOM_error(FATAL, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
end select
endif

if (EOS%use_conT_absS .and. EOS%TFreeze_T_is_potT) then
! absS is set only if EOS%use_conT_absS is True!
T_fr(:) = gsw_ct_from_pt(absS(:),T_fr(:))
endif


end subroutine calculate_TFreeze_array

!> Calls the appropriate subroutine to calculate the freezing point for a 1-D array, taking
!! dimensionally rescaled arguments with factors stored in EOS.
subroutine calculate_TFreeze_1d(S, pressure, T_fr, EOS, dom)
real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt]
real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa]
real, dimension(:), intent(inout) :: T_fr !< Freezing point potential temperature referenced
!! to the surface [C ~> degC]
real, dimension(:), intent(inout) :: T_fr !< Freezing point, either potential temperature referenced to the
!! surface or conservative temperature depending on settings
!! [C ~> degC]
type(EOS_type), intent(in) :: EOS !< Equation of state structure
integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking
!! into account that arrays start at 1.

! Local variables
real, dimension(size(T_fr)) :: pres ! Pressure converted to [Pa]
real, dimension(size(T_fr)) :: Sa ! Salinity converted to [ppt]
real, dimension(size(T_fr)) :: absS ! Salinity converted to absoluate salinity [ppt]
real, dimension(size(T_fr)) :: TFreeze_S ! The salinity for the freezing equation in model units [S ~> PSU or ppt]
integer :: i, is, ie, npts

if (present(dom)) then
Expand All @@ -631,24 +683,36 @@ subroutine calculate_TFreeze_1d(S, pressure, T_fr, EOS, dom)
is = 1 ; ie = size(T_Fr) ; npts = 1 + ie - is
endif

if (EOS%use_conT_absS) then
! Otherwise absS is unneeded and therefore unset
absS(:) = S(:)*EOS%S_to_ppt
if (EOS%TFreeze_S_is_pracS) then
TFreeze_S(:) = gsw_sp_from_sr(absS(:))*EOS%ppt_to_S
else
TFreeze_S(:) = S(:)
endif
else
TFreeze_S(:) = S(:)
endif

if ((EOS%RL2_T2_to_Pa == 1.0) .and. (EOS%S_to_ppt == 1.0)) then
select case (EOS%form_of_TFreeze)
case (TFREEZE_LINEAR)
call calculate_TFreeze_linear(S, pressure, T_fr, is, npts, &
call calculate_TFreeze_linear(TFreeze_S, pressure, T_fr, is, npts, &
EOS%TFr_S0_P0, EOS%dTFr_dS, EOS%dTFr_dp)
case (TFREEZE_MILLERO)
call calculate_TFreeze_Millero(S, pressure, T_fr, is, npts)
call calculate_TFreeze_Millero(TFreeze_S, pressure, T_fr, is, npts)
case (TFREEZE_TEOSPOLY)
call calculate_TFreeze_TEOS_poly(S, pressure, T_fr, is, npts)
call calculate_TFreeze_TEOS_poly(TFreeze_S, pressure, T_fr, is, npts)
case (TFREEZE_TEOS10)
call calculate_TFreeze_teos10(S, pressure, T_fr, is, npts)
call calculate_TFreeze_teos10(TFreeze_S, pressure, T_fr, is, npts)
case default
call MOM_error(FATAL, "calculate_TFreeze_scalar: form_of_TFreeze is not valid.")
end select
else
do i=is,ie
pres(i) = EOS%RL2_T2_to_Pa * pressure(i)
Sa(i) = EOS%S_to_ppt * S(i)
Sa(i) = EOS%S_to_ppt * TFreeze_S(i)
enddo
select case (EOS%form_of_TFreeze)
case (TFREEZE_LINEAR)
Expand All @@ -665,6 +729,13 @@ subroutine calculate_TFreeze_1d(S, pressure, T_fr, EOS, dom)
end select
endif

if (EOS%use_conT_absS .and. EOS%TFreeze_T_is_potT) then
! absS is set only if EOS%use_conT_absS is True!
! absS is in ppt and T_fr is in degC at this point.
T_fr(:) = gsw_ct_from_pt(absS(:),T_fr(:))
endif


if (EOS%degC_to_C /= 1.0) then
do i=is,ie ; T_fr(i) = EOS%degC_to_C * T_fr(i) ; enddo
endif
Expand Down Expand Up @@ -1462,17 +1533,19 @@ end function get_EOS_name

!> Initializes EOS_type by allocating and reading parameters. The scaling factors in
!! US are stored in EOS for later use.
subroutine EOS_init(param_file, EOS, US)
subroutine EOS_init(param_file, EOS, US, use_conT_absS)
type(param_file_type), intent(in) :: param_file !< Parameter file structure
type(EOS_type), intent(inout) :: EOS !< Equation of state structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
logical, intent(in), optional :: use_conT_absS !< True if the model is formulated for
!! conservative temp and absolute salinity
optional :: US
! Local variables
# include "version_variable.h"
character(len=40) :: mdl = "MOM_EOS" ! This module's name.
character(len=12) :: TFREEZE_DEFAULT ! The default freezing point expression
character(len=40) :: tmpstr
logical :: EOS_quad_default
logical :: EOS_quad_default, EOS_TS_default
real :: Rho_Tref_Sref ! Density at Tref degC and Sref ppt [kg m-3]
real :: Tref ! Reference temperature [degC]
real :: Sref ! Reference salinity [psu]
Expand Down Expand Up @@ -1551,6 +1624,12 @@ subroutine EOS_init(param_file, EOS, US)
call EOS_manual_init(EOS, form_of_EOS=EOS_WRIGHT, use_Wright_2nd_deriv_bug=EOS%use_Wright_2nd_deriv_bug)
endif

if (present(use_conT_absS)) then
EOS%use_conT_absS = use_conT_absS
else
EOS%use_conT_absS = .false. ! Assuming it is not needed, it is set to false
endif

EOS_quad_default = .not.((EOS%form_of_EOS == EOS_LINEAR) .or. &
(EOS%form_of_EOS == EOS_WRIGHT) .or. &
(EOS%form_of_EOS == EOS_WRIGHT_REDUCED) .or. &
Expand Down Expand Up @@ -1599,10 +1678,25 @@ subroutine EOS_init(param_file, EOS, US)
units="degC Pa-1", default=0.0)
endif

if ((EOS%form_of_TFreeze==TFREEZE_TEOSPOLY) .or. (EOS%form_of_TFreeze==TFREEZE_TEOS10)) then
! Which default is appropriate for Millero?
EOS_TS_default = .false.
else
EOS_TS_default = .true.
endif
call get_param(param_file, mdl, "TFREEZE_S_IS_PRACS", EOS%TFreeze_S_is_pracS, &
"When True, the model will check if the model internal salinity is "//&
"practical salinity. If the model uses absolute salinity, a "//&
"conversion will be applied.", default=EOS_TS_default)
call get_param(param_file, mdl, "TFREEZE_T_IS_POTT", EOS%TFreeze_T_is_potT, &
"When True, the model will check if the model internal temperature is "//&
"potential temperature. If the model uses conservative temperature, a "//&
"conversion will be applied.", default=EOS_TS_default)

if ((EOS%form_of_EOS == EOS_TEOS10 .or. EOS%form_of_EOS == EOS_ROQUET_RHO .or. &
EOS%form_of_EOS == EOS_ROQUET_SPV) .and. &
.not.((EOS%form_of_TFreeze == TFREEZE_TEOS10) .or. (EOS%form_of_TFreeze == TFREEZE_TEOSPOLY)) ) then
call MOM_error(FATAL, "interpret_eos_selection: EOS_TEOS10 or EOS_ROQUET_RHO or EOS_ROQUET_SPV "//&
call MOM_error(WARNING, "interpret_eos_selection: EOS_TEOS10 or EOS_ROQUET_RHO or EOS_ROQUET_SPV "//&
"should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 or TFREEZE_TEOSPOLY.")
endif

Expand Down
4 changes: 2 additions & 2 deletions src/equation_of_state/MOM_EOS_TEOS10.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,15 @@ module MOM_EOS_TEOS10

! This file is part of MOM6. See LICENSE.md for the license.

use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp
use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp, gsw_ct_from_pt
use gsw_mod_toolbox, only : gsw_rho, gsw_specvol
use gsw_mod_toolbox, only : gsw_rho_first_derivatives, gsw_specvol_first_derivatives
use gsw_mod_toolbox, only : gsw_rho_second_derivatives
use MOM_EOS_base_type, only : EOS_base

implicit none ; private

public gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp
public gsw_sp_from_sr, gsw_pt_from_ct, gsw_sr_from_sp, gsw_ct_from_pt
public TEOS10_EOS

real, parameter :: Pa2db = 1.e-4 !< The conversion factor from Pa to dbar [dbar Pa-1]
Expand Down