Skip to content
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

(*)Correct two rescaling bugs in MOM_EOS #272

Merged
merged 3 commits into from
Dec 19, 2022
Merged

(*)Correct two rescaling bugs in MOM_EOS #272

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2efe1e5
(*)Correct two rescaling bugs in MOM_EOS
Hallberg-NOAA Dec 7, 2022
07b21f4
Merge branch 'dev/gfdl' into second_derivs_scalar_fix
marshallward Dec 18, 2022
632f07e
Merge branch 'dev/gfdl' into second_derivs_scalar_fix
marshallward Dec 19, 2022
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
6 changes: 5 additions & 1 deletion src/core/MOM_checksum_packages.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric)
intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1] or [m s-1]..
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type, which is
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type, which is
!! used to rescale u and v if present.
integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0).
logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully
Expand Down Expand Up @@ -130,6 +130,10 @@ subroutine MOM_thermo_chksum(mesg, tv, G, US, haloshift)
scale=US%Q_to_J_kg*US%R_to_kg_m3*US%Z_to_m)
if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs, &
scale=US%S_to_ppt*US%RZ_to_kg_m2)
if (associated(tv%varT)) call hchksum(tv%varT, mesg//" varT", G%HI, haloshift=hs, scale=US%C_to_degC**2)
if (associated(tv%varS)) call hchksum(tv%varS, mesg//" varS", G%HI, haloshift=hs, scale=US%S_to_ppt**2)
if (associated(tv%covarTS)) call hchksum(tv%covarTS, mesg//" covarTS", G%HI, haloshift=hs, &
scale=US%S_to_ppt*US%C_to_degC)

end subroutine MOM_thermo_chksum

Expand Down
65 changes: 43 additions & 22 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 @@ -126,16 +126,18 @@ module MOM_EOS
real :: dTFr_dp !< The derivative of freezing point with pressure [degC Pa-1]

! Unit conversion factors (normally used for dimensional testing but could also allow for
! change of units of arguments to functions)
real :: m_to_Z = 1. !< A constant that translates distances in meters to the units of depth.
real :: kg_m3_to_R = 1. !< A constant that translates kilograms per meter cubed to the units of density.
real :: R_to_kg_m3 = 1. !< A constant that translates the units of density to kilograms per meter cubed.
real :: RL2_T2_to_Pa = 1.!< Convert pressures from R L2 T-2 to Pa.
real :: L_T_to_m_s = 1. !< Convert lateral velocities from L T-1 to m s-1.
real :: degC_to_C = 1. !< A constant that translates degrees Celsius to the units of temperature.
real :: C_to_degC = 1. !< A constant that translates the units of temperature to degrees Celsius.
real :: ppt_to_S = 1. !< A constant that translates parts per thousand to the units of salinity.
real :: S_to_ppt = 1. !< A constant that translates the units of salinity to parts per thousand.
! change of units of arguments to functions
real :: m_to_Z = 1. !< A constant that translates distances in meters to the units of depth [Z m-1 ~> 1]
real :: kg_m3_to_R = 1. !< A constant that translates kilograms per meter cubed to the
!! units of density [R m3 kg-1 ~> 1]
real :: R_to_kg_m3 = 1. !< A constant that translates the units of density to
!! kilograms per meter cubed [kg m-3 R-1 ~> 1]
real :: RL2_T2_to_Pa = 1.!< Convert pressures from R L2 T-2 to Pa [Pa T2 R-1 L-2 ~> 1]
real :: L_T_to_m_s = 1. !< Convert lateral velocities from L T-1 to m s-1 [m T s-1 L-1 ~> 1]
real :: degC_to_C = 1. !< A constant that translates degrees Celsius to the units of temperature [C degC-1 ~> 1]
real :: C_to_degC = 1. !< A constant that translates the units of temperature to degrees Celsius [degC C-1 ~> 1]
real :: ppt_to_S = 1. !< A constant that translates parts per thousand to the units of salinity [S ppt-1 ~> 1]
real :: S_to_ppt = 1. !< A constant that translates the units of salinity to parts per thousand [ppt S-1 ~> 1]

! logical :: test_EOS = .true. ! If true, test the equation of state
end type EOS_type
Expand Down Expand Up @@ -219,7 +221,11 @@ subroutine calculate_stanley_density_scalar(T, S, pressure, Tvar, TScov, Svar, r
real, optional, intent(in) :: scale !< A multiplicative factor by which to scale output density in
!! combination with scaling stored in EOS [various]
! Local variables
real :: d2RdTT, d2RdST, d2RdSS, d2RdSp, d2RdTp ! Second derivatives of density wrt T,S,p
real :: d2RdTT ! Second derivative of density with temperature [kg m-3 degC-2]
real :: d2RdST ! Second derivative of density with temperature and salinity [kg m-3 degC-1 ppt-1]
real :: d2RdSS ! Second derivative of density with salinity [kg m-3 ppt-2]
real :: d2RdSp ! Second derivative of density with salinity and pressure [kg m-3 ppt-1 Pa-1]
real :: d2RdTp ! Second derivative of density with temperature and pressure [kg m-3 degC-1 Pa-1]
real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1]
real :: T_scale ! A factor to convert temperature to units of degC [degC C-1 ~> 1]
real :: S_scale ! A factor to convert salinity to units of ppt [ppt S-1 ~> 1]
Expand Down Expand Up @@ -309,7 +315,12 @@ subroutine calculate_stanley_density_array(T, S, pressure, Tvar, TScov, Svar, rh
real, optional, intent(in) :: scale !< A multiplicative factor by which to scale the output
!! density, perhaps to other units than kg m-3 [various]
! Local variables
real, dimension(size(T)) :: d2RdTT, d2RdST, d2RdSS, d2RdSp, d2RdTp ! Second derivatives of density wrt T,S,p
real, dimension(size(T)) :: &
d2RdTT, & ! Second derivative of density with temperature [kg m-3 degC-2]
d2RdST, & ! Second derivative of density with temperature and salinity [kg m-3 degC-1 ppt-1]
d2RdSS, & ! Second derivative of density with salinity [kg m-3 ppt-2]
d2RdSp, & ! Second derivative of density with salinity and pressure [kg m-3 ppt-1 Pa-1]
d2RdTp ! Second derivative of density with temperature and pressure [kg m-3 degC-1 Pa-1]
integer :: j

select case (EOS%form_of_EOS)
Expand Down Expand Up @@ -423,7 +434,12 @@ subroutine calculate_stanley_density_1d(T, S, pressure, Tvar, TScov, Svar, rho,
real, dimension(size(rho)) :: pres ! Pressure converted to [Pa]
real, dimension(size(rho)) :: Ta ! Temperature converted to [degC]
real, dimension(size(rho)) :: Sa ! Salinity converted to [ppt]
real, dimension(size(T)) :: d2RdTT, d2RdST, d2RdSS, d2RdSp, d2RdTp ! Second derivatives of density wrt T,S,p
real, dimension(size(T)) :: &
d2RdTT, & ! Second derivative of density with temperature [kg m-3 degC-2]
d2RdST, & ! Second derivative of density with temperature and salinity [kg m-3 degC-1 ppt-1]
d2RdSS, & ! Second derivative of density with salinity [kg m-3 ppt-2]
d2RdSp, & ! Second derivative of density with salinity and pressure [kg m-3 ppt-1 Pa-1]
d2RdTp ! Second derivative of density with temperature and pressure [kg m-3 degC-1 Pa-1]
integer :: i, is, ie, npts

if (present(dom)) then
Expand Down Expand Up @@ -670,7 +686,7 @@ 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.
real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1]
integer :: j

p_scale = 1.0 ; if (present(pres_scale)) p_scale = pres_scale
Expand Down Expand Up @@ -1028,7 +1044,6 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr
!! in combination with scaling stored in EOS [various]
! Local variables
real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1]
real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1]
real :: pres ! Pressure converted to [Pa]
real :: Ta ! Temperature converted to [degC]
real :: Sa ! Salinity converted to [ppt]
Expand Down Expand Up @@ -1061,9 +1076,9 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr
drho_dT_dP = rho_scale * drho_dT_dP
endif

if (p_scale /= 1.0) then
drho_dS_dP = p_scale * drho_dS_dP
drho_dT_dP = p_scale * drho_dT_dP
if (EOS%RL2_T2_to_Pa /= 1.0) then
drho_dS_dP = EOS%RL2_T2_to_Pa * drho_dS_dP
drho_dT_dP = EOS%RL2_T2_to_Pa * drho_dT_dP
endif

if (EOS%C_to_degC /= 1.0) then
Expand Down Expand Up @@ -1173,7 +1188,7 @@ subroutine calc_spec_vol_derivs_1d(T, S, pressure, dSV_dT, dSV_dS, EOS, dom, sca
if (present(scale)) spv_scale = spv_scale * scale
dSVdT_scale = spv_scale * EOS%C_to_degC
dSVdS_scale = spv_scale * EOS%S_to_ppt
if (spv_scale /= 1.0) then ; do i=is,ie
if ((dSVdT_scale /= 1.0) .or. (dSVdS_scale /= 1.0)) then ; do i=is,ie
dSV_dT(i) = dSVdT_scale * dSV_dT(i)
dSV_dS(i) = dSVdS_scale * dSV_dS(i)
enddo ; endif
Expand Down Expand Up @@ -1252,7 +1267,12 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS)

! Local variables
! These arrays use the same units as their counterparts in calcluate_compress_1d.
real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa
real, dimension(1) :: pa ! Pressure in a size-1 1d array [R L2 T-2 ~> Pa]
real, dimension(1) :: Ta ! Temperature in a size-1 1d array [C ~> degC]
real, dimension(1) :: Sa ! Salinity in a size-1 1d array [S ~> ppt]
real, dimension(1) :: rhoa ! In situ density in a size-1 1d array [R ~> kg m-3]
real, dimension(1) :: drho_dpa ! The partial derivative of density with pressure (also the
! inverse of the square of sound speed) in a 1d array [T2 L-2 ~> s2 m-2]

Ta(1) = T ; Sa(1) = S ; pa(1) = pressure

Expand Down Expand Up @@ -1629,11 +1649,12 @@ subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS)
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
intent(inout) :: S !< Salinity [S ~> ppt]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
intent(in) :: mask_z !< 3d mask regulating which points to convert.
intent(in) :: mask_z !< 3d mask regulating which points to convert [nondim]
type(EOS_type), intent(in) :: EOS !< Equation of state structure

real :: gsw_sr_from_sp ! Reference salinity after conversion from practical salinity [ppt]
real :: gsw_ct_from_pt ! Conservative temperature after conversion from potential temperature [degC]
integer :: i, j, k
real :: gsw_sr_from_sp, gsw_ct_from_pt

if ((EOS%form_of_EOS /= EOS_TEOS10) .and. (EOS%form_of_EOS /= EOS_NEMO)) return

Expand Down