Skip to content

Commit

Permalink
+Rescaled Boussinesq pressure force calculations
Browse files Browse the repository at this point in the history 
  Rescaled the pressures in Boussinesq pressure force calculations, including
changing the units of the densities passed to set_pbce_Bouss and using the new
rho_scale and pres_scale arguments to the equation of state routines.  All
answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 4, 2020
1 parent 5a56df7 commit 63c4510
Show file tree
Hide file tree
Showing 3 changed files with 88 additions and 107 deletions.
26 changes: 13 additions & 13 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ module MOM_PressureForce_Mont
type, public :: PressureForce_Mont_CS ; private
logical :: tides !< If true, apply tidal momentum forcing.
real :: Rho0 !< The density used in the Boussinesq
!! approximation [kg m-3].
!! approximation [R ~> kg m-3].
real :: GFS_scale !< Ratio between gravity applied to top interface and the
!! gravitational acceleration of the planet [nondim].
!! Usually this ratio is 1.
Expand Down Expand Up @@ -401,7 +401,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
! attraction and loading, in depth units [Z ~> m].
real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
! density [Pa] (usually 2e7 Pa = 2000 dbar).
real :: I_Rho0 ! 1/Rho0 [m3 kg-1].
real :: I_Rho0 ! 1/Rho0 [R-1 ~> m3 kg-1].
real :: G_Rho0 ! G_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
real :: PFu_bc, PFv_bc ! The pressure gradient force due to along-layer
! compensated density gradients [L T-2 ~> m s-2]
Expand Down Expand Up @@ -520,7 +520,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
M(i,j,1) = CS%GFS_scale * (rho_star(i,j,1) * e(i,j,1))
if (use_p_atm) M(i,j,1) = M(i,j,1) + US%m_s_to_L_T**2*p_atm(i,j) * I_Rho0
if (use_p_atm) M(i,j,1) = M(i,j,1) + US%kg_m3_to_R*US%m_s_to_L_T**2*p_atm(i,j) * I_Rho0
enddo
do k=2,nz ; do i=Isq,Ieq+1
M(i,j,k) = M(i,j,k-1) + (rho_star(i,j,k) - rho_star(i,j,k-1)) * e(i,j,K)
Expand All @@ -531,7 +531,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
do j=Jsq,Jeq+1
do i=Isq,Ieq+1
M(i,j,1) = GV%g_prime(1) * e(i,j,1)
if (use_p_atm) M(i,j,1) = M(i,j,1) + US%m_s_to_L_T**2*p_atm(i,j) * I_Rho0
if (use_p_atm) M(i,j,1) = M(i,j,1) + US%kg_m3_to_R*US%m_s_to_L_T**2*p_atm(i,j) * I_Rho0
enddo
do k=2,nz ; do i=Isq,Ieq+1
M(i,j,k) = M(i,j,k-1) + GV%g_prime(K) * e(i,j,K)
Expand Down Expand Up @@ -609,7 +609,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface height [Z ~> m].
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, intent(in) :: Rho0 !< The "Boussinesq" ocean density [kg m-3].
real, intent(in) :: Rho0 !< The "Boussinesq" ocean density [R ~> kg m-3].
real, intent(in) :: GFS_scale !< Ratio between gravity applied to top
!! interface and the gravitational acceleration of
!! the planet [nondim]. Usually this ratio is 1.
Expand All @@ -623,13 +623,13 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)

! Local variables
real :: Ihtot(SZI_(G)) ! The inverse of the sum of the layer thicknesses [H-1 ~> m-1 or m2 kg-1].
real :: press(SZI_(G)) ! Interface pressure [Pa].
real :: press(SZI_(G)) ! Interface pressure [R L2 T-2 ~> Pa].
real :: T_int(SZI_(G)) ! Interface temperature [degC].
real :: S_int(SZI_(G)) ! Interface salinity [ppt].
real :: dR_dT(SZI_(G)) ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
real :: dR_dS(SZI_(G)) ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
real :: rho_in_situ(SZI_(G)) ! In-situ density at the top of a layer [R ~> kg m-3].
real :: G_Rho0 ! A scaled version of g_Earth / Rho0 [L2 m3 Z-1 T-2 kg-1 ~> m4 s-2 kg-1]
real :: G_Rho0 ! A scaled version of g_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
real :: Rho0xG ! g_Earth * Rho0 [kg s-2 m-1 Z-1 ~> kg s-2 m-2]
logical :: use_EOS ! If true, density is calculated from T & S using
! an equation of state.
Expand All @@ -639,7 +639,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)

Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke

Rho0xG = Rho0*US%L_T_to_m_s**2 * GV%g_Earth
Rho0xG = Rho0 * GV%g_Earth
G_Rho0 = GV%g_Earth / GV%Rho0
use_EOS = associated(tv%eqn_of_state)
z_neglect = GV%H_subroundoff*GV%H_to_Z
Expand All @@ -664,8 +664,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
press(i) = -Rho0xG*e(i,j,1)
enddo
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, Isq, Ieq-Isq+2, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%R_to_kg_m3*US%L_T_to_m_s**2)
do i=Isq,Ieq+1
pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z
enddo
Expand All @@ -675,8 +675,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
enddo
call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, &
Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, Isq, Ieq-Isq+2, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%R_to_kg_m3*US%L_T_to_m_s**2)
do i=Isq,Ieq+1
pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * &
((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * &
Expand Down Expand Up @@ -851,7 +851,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_
"calculate accelerations and the mass for conservation "//&
"properties, or with BOUSSINSEQ false to convert some "//&
"parameters from vertical units of m to kg m-2.", &
units="kg m-3", default=1035.0)
units="kg m-3", default=1035.0, scale=US%R_to_kg_m3)
call get_param(param_file, mdl, "TIDES", CS%tides, &
"If true, apply tidal momentum forcing.", default=.false.)
call get_param(param_file, mdl, "USE_EOS", use_EOS, default=.true., &
Expand Down
81 changes: 35 additions & 46 deletions src/core/MOM_PressureForce_analytic_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ module MOM_PressureForce_AFV
type, public :: PressureForce_AFV_CS ; private
logical :: tides !< If true, apply tidal momentum forcing.
real :: Rho0 !< The density used in the Boussinesq
!! approximation [kg m-3].
!! approximation [R ~> kg m-3].
real :: GFS_scale !< A scaling of the surface pressure gradients to
!! allow the use of a reduced gravity model [nondim].
type(time_type), pointer :: Time !< A pointer to the ocean model's clock.
Expand Down Expand Up @@ -77,7 +77,7 @@ subroutine PressureForce_AFV(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbc
!! or atmosphere-ocean interface [Pa].
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
!! anomaly in each layer due to eta anomalies
!! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
!! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
!! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
!! contributions or compressibility compensation.
Expand Down Expand Up @@ -113,7 +113,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
!! or atmosphere-ocean interface [Pa].
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
!! anomaly in each layer due to eta anomalies
!! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
!! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
!! calculate PFu and PFv [H ~> m or kg m-2], with any tidal
!! contributions or compressibility compensation.
Expand Down Expand Up @@ -194,7 +194,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS

dp_neglect = GV%H_to_Pa * GV%H_subroundoff
alpha_ref = 1.0/CS%Rho0
alpha_ref = 1.0/(US%R_to_kg_m3*CS%Rho0)
g_Earth_z = US%L_T_to_m_s**2 * GV%g_Earth
I_gEarth = 1.0 / g_Earth_z

Expand Down Expand Up @@ -456,7 +456,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
!! or atmosphere-ocean interface [Pa].
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: pbce !< The baroclinic pressure
!! anomaly in each layer due to eta anomalies
!! [m2 s-2 H-1 ~> m s-2 or m4 s-2 kg-1].
!! [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The bottom mass used to
!! calculate PFu and PFv [H ~> m or kg m-2], with any
!! tidal contributions or compressibility compensation.
Expand All @@ -471,22 +471,21 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
Rho_cv_BL ! The coordinate potential density in the deepest variable
! density near-surface layer [R ~> kg m-3].
real, dimension(SZI_(G),SZJ_(G)) :: &
dz_geo, & ! The change in geopotential thickness through a layer times some dimensional
! rescaling factors [kg m-1 R-1 s-2 ~> m2 s-2].
dz_geo, & ! The change in geopotential thickness through a layer [L2 T-2 ~> m2 s-2].
pa, & ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the
! the interface atop a layer [Pa].
! the interface atop a layer [R L2 T-2 ~> Pa].
dpa, & ! The change in pressure anomaly between the top and bottom
! of a layer [Pa].
! of a layer [R L2 T-2 ~> Pa].
intz_dpa ! The vertical integral in depth of the pressure anomaly less the
! pressure anomaly at the top of the layer [H Pa ~> m Pa or kg m-2 Pa].
! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa or kg m-2 Pa].
real, dimension(SZIB_(G),SZJ_(G)) :: &
intx_pa, & ! The zonal integral of the pressure anomaly along the interface
! atop a layer, divided by the grid spacing [Pa].
intx_dpa ! The change in intx_pa through a layer [Pa].
! atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].
intx_dpa ! The change in intx_pa through a layer [R L2 T-2 ~> Pa].
real, dimension(SZI_(G),SZJB_(G)) :: &
inty_pa, & ! The meridional integral of the pressure anomaly along the
! interface atop a layer, divided by the grid spacing [Pa].
inty_dpa ! The change in inty_pa through a layer [Pa].
! interface atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].
inty_dpa ! The change in inty_pa through a layer [R L2 T-2 ~> Pa].

real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
T_tmp, & ! Temporary array of temperatures where layers that are lighter
Expand All @@ -502,12 +501,9 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
real :: p0(SZI_(G)) ! An array of zeros to use for pressure [Pa].
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m].
real :: g_Earth_mks_z ! A scaled version of g_Earth [m2 Z-1 s-2 ~> m s-2].
real :: g_Earth_z_geo ! Another scaled version of g_Earth [R m5 kg-1 Z-1 s-2 ~> m s-2].
real :: I_Rho0 ! 1/Rho0 times unit scaling factors [L2 m kg-1 s2 T-2 ~> m3 kg-1].
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 :: Rho_ref ! The reference density [R ~> kg m-3].
real :: Rho_ref_mks ! The reference density in mks units [kg m-3].
real :: rho_ref ! The reference density [R ~> kg m-3].
real :: dz_neglect ! A minimal thickness [Z ~> m], like e.
logical :: use_p_atm ! If true, use the atmospheric pressure.
logical :: use_ALE ! If true, use an ALE pressure reconstruction.
Expand All @@ -534,12 +530,9 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at

h_neglect = GV%H_subroundoff
dz_neglect = GV%H_subroundoff * GV%H_to_Z
I_Rho0 = US%m_s_to_L_T**2 / (US%R_to_kg_m3*GV%Rho0)
g_Earth_mks_z = US%L_T_to_m_s**2 * GV%g_Earth
g_Earth_z_geo = US%R_to_kg_m3*US%L_T_to_m_s**2 * GV%g_Earth
I_Rho0 = 1.0 / GV%Rho0
G_Rho0 = GV%g_Earth / GV%Rho0
rho_ref_mks = CS%Rho0
rho_ref = rho_ref_mks*US%kg_m3_to_R
rho_ref = CS%Rho0

if (CS%tides) then
! Determine the surface height anomaly for calculating self attraction
Expand Down Expand Up @@ -651,12 +644,12 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
if (use_p_atm) then
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pa(i,j) = (rho_ref*g_Earth_z_geo)*e(i,j,1) + p_atm(i,j)
pa(i,j) = (rho_ref*GV%g_Earth)*e(i,j,1) + US%kg_m3_to_R*US%m_s_to_L_T**2 * p_atm(i,j)
enddo ; enddo
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
pa(i,j) = (rho_ref*g_Earth_z_geo)*e(i,j,1)
pa(i,j) = (rho_ref*GV%g_Earth)*e(i,j,1)
enddo ; enddo
endif
!$OMP parallel do default(shared)
Expand All @@ -680,24 +673,21 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
! where the layers are located.
if ( use_ALE ) then
if ( CS%Recon_Scheme == 1 ) then
call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), &
S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref_mks, CS%Rho0, g_Earth_mks_z, &
dz_neglect, G%bathyT, G%HI, G%HI, &
tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, &
useMassWghtInterp = CS%useMassWghtInterp)
call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k),&
e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, G%HI, G%HI, &
tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp=CS%useMassWghtInterp, &
rho_scale=US%kg_m3_to_R, pres_scale=US%R_to_kg_m3*US%L_T_to_m_s**2)
elseif ( CS%Recon_Scheme == 2 ) then
call int_density_dz_generic_ppm( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref_mks, CS%Rho0, g_Earth_mks_z, &
G%HI, G%HI, tv%eqn_of_state, dpa, intz_dpa, &
intx_dpa, inty_dpa)
tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, &
GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, &
rho_scale=US%kg_m3_to_R, pres_scale=US%R_to_kg_m3*US%L_T_to_m_s**2)
endif
else
call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref_mks, CS%Rho0, g_Earth_mks_z, G%HI, G%HI, tv%eqn_of_state, &
dpa, intz_dpa, intx_dpa, inty_dpa, &
G%bathyT, dz_neglect, CS%useMassWghtInterp)
rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, &
dpa, intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, &
rho_scale=US%kg_m3_to_R, pres_scale=US%R_to_kg_m3*US%L_T_to_m_s**2)
endif
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand All @@ -706,7 +696,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
dz_geo(i,j) = g_Earth_z_geo * GV%H_to_Z*h(i,j,k)
dz_geo(i,j) = GV%g_Earth * GV%H_to_Z*h(i,j,k)
dpa(i,j) = (GV%Rlay(k) - rho_ref) * dz_geo(i,j)
intz_dpa(i,j) = 0.5*(GV%Rlay(k) - rho_ref) * dz_geo(i,j)*h(i,j,k)
enddo ; enddo
Expand Down Expand Up @@ -767,15 +757,14 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at

if (present(eta)) then
if (CS%tides) then
! eta is the sea surface height relative to a time-invariant geoid, for
! comparison with what is used for eta in btstep. See how e was calculated
! about 200 lines above.
!$OMP parallel do default(shared)
! eta is the sea surface height relative to a time-invariant geoid, for comparison with
! what is used for eta in btstep. See how e was calculated about 200 lines above.
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta(i,j) = e(i,j,1)*GV%Z_to_H + e_tidal(i,j)*GV%Z_to_H
enddo ; enddo
else
!$OMP parallel do default(shared)
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta(i,j) = e(i,j,1)*GV%Z_to_H
enddo ; enddo
Expand Down Expand Up @@ -819,7 +808,7 @@ subroutine PressureForce_AFV_init(Time, G, GV, US, param_file, diag, CS, tides_C
"calculate accelerations and the mass for conservation "//&
"properties, or with BOUSSINSEQ false to convert some "//&
"parameters from vertical units of m to kg m-2.", &
units="kg m-3", default=1035.0)
units="kg m-3", default=1035.0, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "TIDES", CS%tides, &
"If true, apply tidal momentum forcing.", default=.false.)
call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, &
Expand Down
Loading

0 comments on commit 63c4510

Please sign in to comment.