diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 9c4f355692..7b970f5686 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -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 @@ -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]. @@ -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 @@ -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. @@ -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) @@ -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. ) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 8d30ba28e3..03d5b6131c 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -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, & @@ -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] @@ -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] @@ -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 @@ -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 @@ -260,7 +268,8 @@ 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)) + (wt_R*z_t(i+1,j))) - & + ((wt_L*z0pres(i,j)) + (wt_R*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) @@ -326,7 +335,8 @@ 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)) + (wt_R*z_t(i,j+1))) - & + ((wt_L*z0pres(i,j)) + (wt_R*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) @@ -414,7 +424,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 @@ -464,7 +475,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] @@ -480,7 +491,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. @@ -517,7 +534,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) @@ -610,7 +627,8 @@ 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)) + (w_right*e(i+1,j,K))) - & + ((w_left*z0pres(i,j)) + (w_right*z0pres(i+1,j))) ) ! Pressure do n=2,5 @@ -706,7 +724,8 @@ 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)) + (w_right*e(i,j+1,K))) - & + ((w_left*z0pres(i,j)) + (w_right*z0pres(i,j+1))) ) ! Pressure do n=2,5 @@ -812,7 +831,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 @@ -864,7 +884,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] @@ -879,7 +899,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 @@ -924,7 +950,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) ) @@ -1011,7 +1037,8 @@ 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)) + (w_right*e(i+1,j,K))) - & + ((w_left*z0pres(i,j)) + (w_right*z0pres(i+1,j))) ) do n=2,5 p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_x(m,i) enddo @@ -1116,7 +1143,8 @@ 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)) + (w_right*e(i,j+1,K))) - & + ((w_left*z0pres(i,j)) + (w_right*z0pres(i,j+1))) ) do n=2,5 p15(pos+n) = p15(pos+n-1) + GxRho*0.25*dz_y(m,i) enddo diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index a66b53656c..74f540f64f 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -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 diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index f1893ae80f..6fbff16f11 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -435,12 +435,14 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! temperature into degC [degC C-1 ~> 1] real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure !! into PSU [PSU S-1 ~> 1]. - 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, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [m3 kg-1] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: p0_2d ! A term in the Wright EOS [Pa] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: lambda_2d ! A term in the Wright EOS [m2 s-2] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: al0 ! A term in the Wright EOS [m3 kg-1] real :: p0 ! A term in the Wright EOS [Pa] real :: lambda ! A term in the Wright EOS [m2 s-2] @@ -466,7 +468,6 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]. real :: Pa_to_RL2_T2 ! A conversion factor of pressures from Pa to the output units indicated by ! pres_scale [R L2 T-2 Pa-1 ~> 1]. - real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1] real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1] real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1] @@ -504,7 +505,13 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & else rho_ref_mks = rho_ref ; I_Rho = 1.0 / rho_0 endif - 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 a1s = a1 ; a2s = a2 b1s = b1 ; b2s = b2 ; b3s = b3 ; b4s = b4 ; b5s = b5 @@ -541,7 +548,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dz = z_t(i,j) - z_b(i,j) - p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres) + p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j)) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -585,7 +592,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j)) dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j))) - p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - z0pres) + p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - & + ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -626,7 +634,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1)) dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1))) - p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - z0pres) + p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - & + ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) diff --git a/src/equation_of_state/MOM_EOS_Wright_full.F90 b/src/equation_of_state/MOM_EOS_Wright_full.F90 index 576f390c70..1f03c1f4aa 100644 --- a/src/equation_of_state/MOM_EOS_Wright_full.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_full.F90 @@ -441,12 +441,14 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! temperature into degC [degC C-1 ~> 1] real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure !! into PSU [PSU S-1 ~> 1]. - 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, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [m3 kg-1] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: p0_2d ! A term in the Wright EOS [Pa] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: lambda_2d ! A term in the Wright EOS [m2 s-2] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: al0 ! A term in the Wright EOS [m3 kg-1] real :: p0 ! A term in the Wright EOS [Pa] real :: lambda ! A term in the Wright EOS [m2 s-2] @@ -472,7 +474,6 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]. real :: Pa_to_RL2_T2 ! A conversion factor of pressures from Pa to the output units indicated by ! pres_scale [R L2 T-2 Pa-1 ~> 1]. - real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1] real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1] real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1] @@ -510,7 +511,13 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & else rho_ref_mks = rho_ref ; I_Rho = 1.0 / rho_0 endif - 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 a1s = a1 ; a2s = a2 b1s = b1 ; b2s = b2 ; b3s = b3 ; b4s = b4 ; b5s = b5 @@ -547,7 +554,7 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dz = z_t(i,j) - z_b(i,j) - p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres) + p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j)) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0) @@ -590,7 +597,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j)) dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j))) - p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - z0pres) + p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - & + ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0) @@ -631,7 +639,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1)) dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1))) - p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - z0pres) + p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - & + ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0) diff --git a/src/equation_of_state/MOM_EOS_Wright_red.F90 b/src/equation_of_state/MOM_EOS_Wright_red.F90 index 58a9e2bce1..d46707e911 100644 --- a/src/equation_of_state/MOM_EOS_Wright_red.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_red.F90 @@ -443,12 +443,14 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! temperature into degC [degC C-1 ~> 1] real, optional, intent(in) :: saln_scale !< A multiplicative factor to convert pressure !! into PSU [PSU S-1 ~> 1]. - 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, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d ! A term in the Wright EOS [m3 kg-1] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: p0_2d ! A term in the Wright EOS [Pa] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: lambda_2d ! A term in the Wright EOS [m2 s-2] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: al0 ! A term in the Wright EOS [m3 kg-1] real :: p0 ! A term in the Wright EOS [Pa] real :: lambda ! A term in the Wright EOS [m2 s-2] @@ -474,7 +476,6 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa]. real :: Pa_to_RL2_T2 ! A conversion factor of pressures from Pa to the output units indicated by ! pres_scale [R L2 T-2 Pa-1 ~> 1]. - real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: a1s ! Partly rescaled version of a1 [m3 kg-1 C-1 ~> m3 kg-1 degC-1] real :: a2s ! Partly rescaled version of a2 [m3 kg-1 S-1 ~> m3 kg-1 PSU-1] real :: b1s ! Partly rescaled version of b1 [Pa C-1 ~> Pa degC-1] @@ -512,7 +513,13 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & else rho_ref_mks = rho_ref ; I_Rho = 1.0 / rho_0 endif - 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 a1s = a1 ; a2s = a2 b1s = b1 ; b2s = b2 ; b3s = b3 ; b4s = b4 ; b5s = b5 @@ -549,7 +556,7 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dz = z_t(i,j) - z_b(i,j) - p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres) + p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j)) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0) @@ -592,7 +599,8 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j)) dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j))) - p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - z0pres) + p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - & + ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0) @@ -633,7 +641,9 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1)) dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1))) - p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - z0pres) + p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - & + ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1)))) + I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)