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

Jump to bottom

+Optionally use SSH in calculate density for PGF #696

Merged
merged 3 commits into from
Sep 11, 2024
Merged

+Optionally use SSH in calculate density for PGF #696

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
54aaffc
+Optionally use SSH in calculate density for PGF
Hallberg-NOAA Jul 29, 2024
2e3488f
(*)Refactor p_ave calculation
Hallberg-NOAA Sep 7, 2024
7707230
Merge branch 'dev/gfdl' into add_SSH_IN_EOS_PRESSURE_FOR_PGF
marshallward Sep 10, 2024
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
32 changes: 29 additions & 3 deletions src/core/MOM_PressureForce_FV.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,9 @@ module MOM_PressureForce_FV
!! for the finite volume pressure gradient calculation.
!! By the default (1) is for a piecewise linear method

logical :: use_SSH_in_Z0p !< If true, adjust the height at which the pressure used in the
!! equation of state is 0 to account for the displacement of the sea
!! surface including adjustments for atmospheric or sea-ice pressure.
logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF
integer :: tides_answer_date !< Recover old answers with tides in Boussinesq mode
integer :: id_e_tide = -1 !< Diagnostic identifier
Expand Down Expand Up @@ -522,6 +525,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
! [Z ~> m].
e_tide_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
! specific to tides [Z ~> m].
Z_0p, & ! The height at which the pressure used in the equation of state is 0 [Z ~> m]
SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
dM ! The barotropic adjustment to the Montgomery potential to
! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
Expand Down Expand Up @@ -577,6 +581,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
! in roundoff and can be neglected [H ~> m].
real :: I_Rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1].
real :: G_Rho0 ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1]
real :: rho_ref ! The reference density [R ~> kg m-3].
real :: dz_neglect ! A minimal thickness [Z ~> m], like e.
real :: H_to_RL2_T2 ! A factor to convert from thickness units (H) to pressure
Expand Down Expand Up @@ -753,6 +758,21 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
enddo ; enddo
endif

if (CS%use_SSH_in_Z0p .and. use_p_atm) then
I_g_rho = 1.0 / (CS%rho0*GV%g_Earth)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Z_0p(i,j) = e(i,j,1) + p_atm(i,j) * I_g_rho
enddo ; enddo
elseif (CS%use_SSH_in_Z0p) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Z_0p(i,j) = e(i,j,1)
enddo ; enddo
else
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Z_0p(i,j) = G%Z_ref
enddo ; enddo
endif

do k=1,nz
! Calculate 4 integrals through the layer that are required in the
! subsequent calculation.
Expand All @@ -769,19 +789,19 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &
intx_dpa(:,:,k), inty_dpa(:,:,k), &
MassWghtInterp=CS%MassWghtInterp, &
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref)
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=Z_0p)
elseif ( CS%Recon_Scheme == 2 ) then
call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &
intx_dpa(:,:,k), inty_dpa(:,:,k), &
MassWghtInterp=CS%MassWghtInterp, Z_0p=G%Z_ref)
MassWghtInterp=CS%MassWghtInterp, Z_0p=Z_0p)
endif
else
call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref, CS%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa(:,:,k), &
intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), G%bathyT, dz_neglect, &
CS%MassWghtInterp, Z_0p=G%Z_ref)
CS%MassWghtInterp, Z_0p=Z_0p)
endif
if (GV%Z_to_H /= 1.0) then
!$OMP parallel do default(shared)
Expand Down Expand Up @@ -1042,6 +1062,12 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
endif
call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
"If true, calculate self-attraction and loading.", default=CS%tides)
call get_param(param_file, mdl, "SSH_IN_EOS_PRESSURE_FOR_PGF", CS%use_SSH_in_Z0p, &
"If true, include contributions from the sea surface height in the height-based "//&
"pressure used in the equation of state calculations for the Boussinesq pressure "//&
"gradient forces, including adjustments for atmospheric or sea-ice pressure.", &
default=.false., do_not_log=.not.GV%Boussinesq)

call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, &
"If True, use the ALE algorithm (regridding/remapping). "//&
"If False, use the layered isopycnal algorithm.", default=.false. )
Expand Down
60 changes: 41 additions & 19 deletions src/core/MOM_density_integrals.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,8 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa,
real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m]
integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
!! mass weighting to interpolate T/S in integrals
real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]

if (EOS_quadrature(EOS)) then
call int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, &
Expand Down Expand Up @@ -139,7 +140,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
!! mass weighting to interpolate T/S in integrals
logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of
!! density anomalies, as was used prior to March 2018.
real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]

! Local variables
real :: T5((5*HI%iscB+1):(5*(HI%iecB+2))) ! Temperatures along a line of subgrid locations [C ~> degC]
Expand All @@ -157,7 +159,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
real :: dz ! The layer thickness [Z ~> m]
real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m]
real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m]
real :: z0pres ! The height at which the pressure is zero [Z ~> m]
real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
real :: hWght ! A pressure-thickness below topography [Z ~> m]
real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m]
real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2]
Expand All @@ -184,7 +186,13 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &

GxRho = G_e * rho_0
I_Rho = 1.0 / rho_0
z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p
if (present(Z_0p)) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
z0pres(i,j) = Z_0p(i,j)
enddo ; enddo
else
z0pres(:,:) = 0.0
endif
use_rho_ref = .true.
if (present(use_inaccurate_form)) then
if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form
Expand All @@ -209,7 +217,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
dz = z_t(i,j) - z_b(i,j)
do n=1,5
T5(i*5+n) = T(i,j) ; S5(i*5+n) = S(i,j)
p5(i*5+n) = -GxRho*((z_t(i,j) - z0pres) - 0.25*real(n-1)*dz)
p5(i*5+n) = -GxRho*((z_t(i,j) - z0pres(i,j)) - 0.25*real(n-1)*dz)
enddo
enddo

Expand Down Expand Up @@ -260,7 +268,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
pos = i*15+(m-2)*5
T15(pos+1) = (wtT_L*T(i,j)) + (wtT_R*T(i+1,j))
S15(pos+1) = (wtT_L*S(i,j)) + (wtT_R*S(i+1,j))
p15(pos+1) = -GxRho*(((wt_L*z_t(i,j)) + (wt_R*z_t(i+1,j))) - z0pres)
p15(pos+1) = -GxRho * ((wt_L*(z_t(i,j)-z0pres(i,j))) + (wt_R*(z_t(i+1,j)-z0pres(i+1,j))))
do n=2,5
T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i)
Expand Down Expand Up @@ -326,7 +334,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
pos = i*15+(m-2)*5
T15(pos+1) = (wtT_L*T(i,j)) + (wtT_R*T(i,j+1))
S15(pos+1) = (wtT_L*S(i,j)) + (wtT_R*S(i,j+1))
p15(pos+1) = -GxRho*(((wt_L*z_t(i,j)) + (wt_R*z_t(i,j+1))) - z0pres)
p15(pos+1) = -GxRho * ((wt_L*(z_t(i,j)-z0pres(i,j))) + (wt_R*(z_t(i,j+1)-z0pres(i,j+1))))
do n=2,5
T15(pos+n) = T15(pos+1) ; S15(pos+n) = S15(pos+1)
p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_y(m,i)
Expand Down Expand Up @@ -414,7 +422,8 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
!! mass weighting to interpolate T/S in integrals
logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of
!! density anomalies, as was used prior to March 2018.
real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]

! This subroutine calculates (by numerical quadrature) integrals of
! pressure anomalies across layers, which are required for calculating the
Expand Down Expand Up @@ -464,7 +473,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
real :: massWeightToggle ! A non-dimensional toggle factor (0 or 1) [nondim]
real :: Ttl, Tbl, Ttr, Tbr ! Temperatures at the velocity cell corners [C ~> degC]
real :: Stl, Sbl, Str, Sbr ! Salinities at the velocity cell corners [S ~> ppt]
real :: z0pres ! The height at which the pressure is zero [Z ~> m]
real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
real :: hWght ! A topographically limited thickness weight [Z ~> m]
real :: hL, hR ! Thicknesses to the left and right [Z ~> m]
real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
Expand All @@ -480,7 +489,13 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &

GxRho = G_e * rho_0
I_Rho = 1.0 / rho_0
z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p
if (present(Z_0p)) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
z0pres(i,j) = Z_0p(i,j)
enddo ; enddo
else
z0pres(:,:) = 0.0
endif
massWeightToggle = 0.
if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif
use_rho_ref = .true.
Expand Down Expand Up @@ -517,7 +532,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
do i = Isq,Ieq+1
dz(i) = e(i,j,K) - e(i,j,K+1)
do n=1,5
p5(i*5+n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz(i))
p5(i*5+n) = -GxRho*((e(i,j,K) - z0pres(i,j)) - 0.25*real(n-1)*dz(i))
! Salinity and temperature points are linearly interpolated
S5(i*5+n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * S_b(i,j,k)
T5(i*5+n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * T_b(i,j,k)
Expand Down Expand Up @@ -610,7 +625,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
S15(pos+1) = (w_left*Stl) + (w_right*Str)
S15(pos+5) = (w_left*Sbl) + (w_right*Sbr)

p15(pos+1) = -GxRho*(((w_left*e(i,j,K)) + (w_right*e(i+1,j,K))) - z0pres)
p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i+1,j,K)-z0pres(i+1,j))))

! Pressure
do n=2,5
Expand Down Expand Up @@ -706,7 +721,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, &
S15(pos+1) = (w_left*Stl) + (w_right*Str)
S15(pos+5) = (w_left*Sbl) + (w_right*Sbr)

p15(pos+1) = -GxRho*(((w_left*e(i,j,K)) + (w_right*e(i,j+1,K))) - z0pres)
p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i,j+1,K)-z0pres(i,j+1))))

! Pressure
do n=2,5
Expand Down Expand Up @@ -812,7 +827,8 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
!! divided by the y grid spacing [R L2 T-2 ~> Pa]
integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
!! mass weighting to interpolate T/S in integrals
real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]

! This subroutine calculates (by numerical quadrature) integrals of
! pressure anomalies across layers, which are required for calculating the
Expand Down Expand Up @@ -864,7 +880,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
real :: t6 ! PPM curvature coefficient for T [C ~> degC]
real :: T_top, T_mn, T_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T [C ~> degC]
real :: S_top, S_mn, S_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S [S ~> ppt]
real :: z0pres ! The height at which the pressure is zero [Z ~> m]
real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m]
real :: hWght ! A topographically limited thickness weight [Z ~> m]
real :: hL, hR ! Thicknesses to the left and right [Z ~> m]
real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2]
Expand All @@ -879,7 +895,13 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &

GxRho = G_e * rho_0
I_Rho = 1.0 / rho_0
z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p
if (present(Z_0p)) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
z0pres(i,j) = Z_0p(i,j)
enddo ; enddo
else
z0pres(:,:) = 0.0
endif
massWeightToggle = 0.
if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif

Expand Down Expand Up @@ -924,7 +946,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
endif
dz = e(i,j,K) - e(i,j,K+1)
do n=1,5
p5(I*5+n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz)
p5(I*5+n) = -GxRho*((e(i,j,K) - z0pres(i,j)) - 0.25*real(n-1)*dz)
! Salinity and temperature points are reconstructed with PPM
S5(I*5+n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) )
T5(I*5+n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) )
Expand Down Expand Up @@ -1011,7 +1033,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
dz_x(m,i) = (w_left*(e(i,j,K) - e(i,j,K+1))) + (w_right*(e(i+1,j,K) - e(i+1,j,K+1)))

pos = i*15+(m-2)*5
p15(pos+1) = -GxRho*(((w_left*e(i,j,K)) + (w_right*e(i+1,j,K))) - z0pres)
p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i+1,j,K)-z0pres(i+1,j))))
do n=2,5
p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i)
enddo
Expand Down Expand Up @@ -1116,7 +1138,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
dz_y(m,i) = (w_left*(e(i,j,K) - e(i,j,K+1))) + (w_right*(e(i,j+1,K) - e(i,j+1,K+1)))

pos = i*15+(m-2)*5
p15(pos+1) = -GxRho*(((w_left*e(i,j,K)) + (w_right*e(i,j+1,K))) - z0pres)
p15(pos+1) = -GxRho * ((w_left*(e(i,j,K)-z0pres(i,j))) + (w_right*(e(i,j+1,K)-z0pres(i,j+1))))
do n=2,5
p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_y(m,i)
enddo
Expand Down
3 changes: 2 additions & 1 deletion src/equation_of_state/MOM_EOS.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1345,7 +1345,8 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS,
real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]
integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
!! mass weighting to interpolate T/S in integrals
real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m]

! Local variables
real :: rho_scale ! A multiplicative factor by which to scale density from kg m-3 to the
Expand Down
Loading
Loading