From 54aaffc65bb5f11c2ba5cd29d21152283a73fba4 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 29 Jul 2024 04:33:44 -0400 Subject: [PATCH 1/2] +Optionally use SSH in calculate density for PGF Added the option of including the atmospheric or ice pressure and sea surface height displacements from the global reference height in the pressures used in the density calculations for Boussinesq pressure gradient calculations. Note that the full pressures were already being used everywhere apart from the calculation of the equation of state. This capability is controlled by the new runtime parameter SSH_IN_EOS_PRESSURE_FOR_PGF. This commit changes the Z_0p argument to int_density_dz and 8 other routines (int_density_dz_generic_pcm, int_density_dz_generic_plm, int_density_dz_generic_ppm, analytic_int_density_dz, int_density_dz_wright, int_density_dz_wright_full, int_density_dz_wright_red and int_density_dz_linear) from scalars into 2-d arrays, as were the internal z0pres arrays in most of these routines. By default, all answers are bitwise identical, but there is a new runtime parameter in the MOM_parameter_doc files for Boussinesq cases. --- src/core/MOM_PressureForce_FV.F90 | 32 ++++++++- src/core/MOM_density_integrals.F90 | 66 +++++++++++++------ src/equation_of_state/MOM_EOS.F90 | 3 +- src/equation_of_state/MOM_EOS_Wright.F90 | 21 ++++-- src/equation_of_state/MOM_EOS_Wright_full.F90 | 21 ++++-- src/equation_of_state/MOM_EOS_Wright_red.F90 | 22 +++++-- 6 files changed, 124 insertions(+), 41 deletions(-) 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) From 2e3488fe764e5a8735c385d32b2548a6e59bc303 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 7 Sep 2024 13:53:55 -0400 Subject: [PATCH 2/2] (*)Refactor p_ave calculation Rearranged the calculation of the layer average pressure used in the density integrals to have one less multiply and be clearer to read. This could change answers at roundoff if REFERENCE_HEIGHT is nonzero (it is 0 by default and in all known cases) or if SSH_IN_EOS_PRESSURE_FOR_PGF is set to true. Because the former is mostly used for approximate self-consistency testing and the latter option has just been added, it is unlikely that answers will change for any production runs. --- src/core/MOM_density_integrals.F90 | 18 ++++++------------ src/equation_of_state/MOM_EOS_Wright.F90 | 8 ++++---- src/equation_of_state/MOM_EOS_Wright_full.F90 | 8 ++++---- src/equation_of_state/MOM_EOS_Wright_red.F90 | 9 ++++----- 4 files changed, 18 insertions(+), 25 deletions(-) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 03d5b6131c..8b820594aa 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -268,8 +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))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j))) ) + 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) @@ -335,8 +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))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1))) ) + 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) @@ -627,8 +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))) - & - ((w_left*z0pres(i,j)) + (w_right*z0pres(i+1,j))) ) + 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 @@ -724,8 +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))) - & - ((w_left*z0pres(i,j)) + (w_right*z0pres(i,j+1))) ) + 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 @@ -1037,8 +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))) - & - ((w_left*z0pres(i,j)) + (w_right*z0pres(i+1,j))) ) + 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 @@ -1143,8 +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))) - & - ((w_left*z0pres(i,j)) + (w_right*z0pres(i,j+1))) ) + 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 diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 6fbff16f11..11fa57644d 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -592,8 +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)))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j)))) + p_ave = -GxRho*((wt_L * (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + & + (wt_R * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -634,8 +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)))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1)))) + p_ave = -GxRho*((wt_L*(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + & + (wt_R*(0.5*(z_t(i,j+1)+z_b(i,j+1))-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 1f03c1f4aa..6dba8444dd 100644 --- a/src/equation_of_state/MOM_EOS_Wright_full.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_full.F90 @@ -597,8 +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)))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j)))) + p_ave = -GxRho*((wt_L * (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + & + (wt_R * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0) @@ -639,8 +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)))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1)))) + p_ave = -GxRho*((wt_L*(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + & + (wt_R*(0.5*(z_t(i,j+1)+z_b(i,j+1))-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 d46707e911..87d6a16dba 100644 --- a/src/equation_of_state/MOM_EOS_Wright_red.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_red.F90 @@ -599,8 +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)))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i+1,j)))) + p_ave = -GxRho*((wt_L * (0.5*(z_t(i,j)+z_b(i,j)) - z0pres(i,j))) + & + (wt_R * (0.5*(z_t(i+1,j)+z_b(i+1,j)) - z0pres(i+1,j)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0) @@ -641,9 +641,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,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)))) - & - ((wt_L*z0pres(i,j)) + (wt_R*z0pres(i,j+1)))) - + p_ave = -GxRho*((wt_L*(0.5*(z_t(i,j)+z_b(i,j))-z0pres(i,j))) + & + (wt_R*(0.5*(z_t(i,j+1)+z_b(i,j+1))-z0pres(i,j+1)))) I_al0 = 1.0 / al0 I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)