From 63c451005e241f28a752a1adf6ab97b5a53d0ddd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 4 Apr 2020 09:22:58 -0400 Subject: [PATCH] +Rescaled Boussinesq pressure force calculations 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. --- src/core/MOM_PressureForce_Montgomery.F90 | 26 +++---- src/core/MOM_PressureForce_analytic_FV.F90 | 81 +++++++++----------- src/core/MOM_PressureForce_blocked_AFV.F90 | 88 ++++++++++------------ 3 files changed, 88 insertions(+), 107 deletions(-) diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 43de125701..e0177e35b9 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -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. @@ -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] @@ -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) @@ -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) @@ -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. @@ -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. @@ -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 @@ -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 @@ -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)) * & @@ -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., & diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index 75a2dfad7f..d0a6932810 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -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. @@ -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. @@ -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. @@ -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 @@ -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. @@ -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 @@ -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. @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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, & diff --git a/src/core/MOM_PressureForce_blocked_AFV.F90 b/src/core/MOM_PressureForce_blocked_AFV.F90 index faa7912f1e..60e1330aa6 100644 --- a/src/core/MOM_PressureForce_blocked_AFV.F90 +++ b/src/core/MOM_PressureForce_blocked_AFV.F90 @@ -36,7 +36,7 @@ module MOM_PressureForce_blk_AFV type, public :: PressureForce_blk_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. @@ -77,7 +77,7 @@ subroutine PressureForce_blk_AFV(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, !! 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. @@ -112,7 +112,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, !! 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. @@ -190,7 +190,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, use_EOS = associated(tv%eqn_of_state) 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 @@ -437,7 +437,7 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, !! 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. @@ -452,22 +452,21 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, Rho_cv_BL ! The coordinate potential density in the deepest variable ! density near-surface layer [R ~> kg m-3]. real, dimension(SZDI_(G%Block(1)),SZDJ_(G%Block(1))) :: & ! on block indices - dz_bk, & ! The change in geopotential thickness through a layer times some dimensional - ! rescaling factors [kg m-1 R-1 s-2 ~> m2 s-2]. + dz_bk, & ! The change in geopotential thickness through a layer [L2 T-2 ~> m2 s-2]. pa_bk, & ! 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_bk, & ! The change in pressure anomaly between the top and bottom - ! of a layer [Pa]. + ! of a layer [R L2 T-2 ~> Pa]. intz_dpa_bk ! 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(SZDIB_(G%Block(1)),SZDJ_(G%Block(1))) :: & ! on block indices intx_pa_bk, & ! The zonal integral of the pressure anomaly along the interface - ! atop a layer, divided by the grid spacing [Pa]. - intx_dpa_bk ! The change in intx_pa through a layer [Pa]. + ! atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa]. + intx_dpa_bk ! The change in intx_pa through a layer [R L2 T-2 ~> Pa]. real, dimension(SZDI_(G%Block(1)),SZDJB_(G%Block(1))) :: & ! on block indices inty_pa_bk, & ! The meridional integral of the pressure anomaly along the - ! interface atop a layer, divided by the grid spacing [Pa]. - inty_dpa_bk ! 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_bk ! 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 @@ -483,12 +482,9 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, 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 or kg m-2]. - real :: I_Rho0 ! 1/Rho0 times unit scaling factors [L2 m kg-1 s2 T-2 ~> m3 kg-1]. - 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 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]. real :: G_Rho0 ! G_Earth / Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. - real :: Rho_ref ! The reference density [R-1 ~> 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. @@ -518,12 +514,9 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, 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 @@ -629,9 +622,9 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, endif endif -!$OMP parallel do default(none) shared(use_p_atm,Rho_ref,Rho_ref_mks,G,GV,e,p_atm,nz,use_EOS,& -!$OMP use_ALE,T_t,T_b,S_t,S_b,CS,tv,tv_tmp,g_Earth_z_geo, & -!$OMP g_Earth_mks_z,h,PFu,I_Rho0,h_neglect,dz_neglect,PFv,dM)& +!$OMP parallel do default(none) shared(use_p_atm,rho_ref,G,GV,e,p_atm,nz,use_EOS,& +!$OMP use_ALE,T_t,T_b,S_t,S_b,CS,tv,tv_tmp, & +!$OMP h,PFu,I_Rho0,h_neglect,dz_neglect,PFv,dM)& !$OMP private(is_bk,ie_bk,js_bk,je_bk,Isq_bk,Ieq_bk,Jsq_bk, & !$OMP Jeq_bk,ioff_bk,joff_bk,pa_bk, & !$OMP intx_pa_bk,inty_pa_bk,dpa_bk,intz_dpa_bk, & @@ -650,12 +643,12 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, if (use_p_atm) then do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 i = ib+ioff_bk ; j = jb+joff_bk - pa_bk(ib,jb) = (Rho_ref*g_Earth_z_geo)*e(i,j,1) + p_atm(i,j) + pa_bk(ib,jb) = (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 do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 i = ib+ioff_bk ; j = jb+joff_bk - pa_bk(ib,jb) = (Rho_ref*g_Earth_z_geo)*e(i,j,1) + pa_bk(ib,jb) = (Rho_ref*GV%g_Earth)*e(i,j,1) enddo ; enddo endif do jb=js_bk,je_bk ; do Ib=Isq_bk,Ieq_bk @@ -677,24 +670,24 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ! 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%Block(n), & - tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & - 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%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & + 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%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, & - intx_dpa_bk, inty_dpa_bk) + tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, & + GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, & + intx_dpa_bk, inty_dpa_bk, & + 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%Block(n), tv%eqn_of_state, & + rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, & dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & - G%bathyT, dz_neglect, CS%useMassWghtInterp) + 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 do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 intz_dpa_bk(ib,jb) = intz_dpa_bk(ib,jb)*GV%Z_to_H @@ -702,7 +695,7 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, else do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 i = ib+ioff_bk ; j = jb+joff_bk - dz_bk(ib,jb) = g_Earth_z_geo*GV%H_to_Z*h(i,j,k) + dz_bk(ib,jb) = GV%g_Earth*GV%H_to_Z*h(i,j,k) dpa_bk(ib,jb) = (GV%Rlay(k) - Rho_ref)*dz_bk(ib,jb) intz_dpa_bk(ib,jb) = 0.5*(GV%Rlay(k) - Rho_ref) * dz_bk(ib,jb)*h(i,j,k) enddo ; enddo @@ -759,15 +752,14 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, 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 @@ -811,7 +803,7 @@ subroutine PressureForce_blk_AFV_init(Time, G, GV, US, param_file, diag, CS, tid "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, &