From 4d8d7afc7301188a97b542f0481e2cc683c3d4ac Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 05:01:06 -0400 Subject: [PATCH] +Add optional Z_0p argument to int_density_dz Added an optional Z_0p argument to the various Boussinesq density integral routines as a new intercept for the linear expression between the (arbitrary) interface heights and pressure used for the equation of state routines. If omitted, this is equivalent to setting this intercept to 0, and all answers are bitwise identical. There are new optional arguments to several publicly visible interfaces (that are not yet being exercised with this commit, but have been tested and will be used in an upcoming commit). --- src/core/MOM_density_integrals.F90 | 43 +++++++++++++++--------- src/equation_of_state/MOM_EOS.F90 | 8 +++-- src/equation_of_state/MOM_EOS_Wright.F90 | 20 ++++++----- 3 files changed, 43 insertions(+), 28 deletions(-) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 302ba0a714..04e151d5a7 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -38,7 +38,7 @@ module MOM_density_integrals !! required for calculating the finite-volume form pressure accelerations in a !! Boussinesq model. subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the arrays real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -78,13 +78,14 @@ 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] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, 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, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) else call analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) endif end subroutine int_density_dz @@ -93,8 +94,8 @@ end subroutine int_density_dz !> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which !! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, use_inaccurate_form) + EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dz_neglect, useMassWghtInterp, use_inaccurate_form, Z_0p) type(hor_index_type), intent(in) :: HI !< Horizontal index type for input variables. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature of the layer [degC] @@ -136,6 +137,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! interpolate T/S for top and bottom 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] + ! Local variables real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt] real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa] @@ -148,6 +151,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] real :: dz ! The layer thickness [Z ~> m] + real :: z0pres ! 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] @@ -173,6 +177,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p use_rho_ref = .true. if (present(use_inaccurate_form)) then if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form @@ -191,7 +196,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(n) = T(i,j) ; S5(n) = S(i,j) - p5(n) = -GxRho*(z_t(i,j) - 0.25*real(n-1)*dz) + p5(n) = -GxRho*((z_t(i,j) - z0pres) - 0.25*real(n-1)*dz) enddo if (use_rho_ref) then if (rho_scale /= 1.0) then @@ -245,7 +250,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) S5(1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) - p5(1) = -GxRho*(wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) + p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) - z0pres) do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) + GxRho*0.25*dz enddo @@ -300,7 +305,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) S5(1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) - p5(1) = -GxRho*(wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) + p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) - z0pres) do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) p5(n) = p5(n-1) + GxRho*0.25*dz @@ -336,7 +341,7 @@ end subroutine int_density_dz_generic_pcm subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, dpa, & intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, & - use_inaccurate_form) + use_inaccurate_form, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -379,6 +384,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & !! interpolate T/S for top and bottom 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] ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -427,6 +433,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 [degC] real :: Stl, Sbl, Str, Sbr ! Salinities at the velocity cell corners [ppt] + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thicknes 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] @@ -443,6 +450,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. if (present(useMassWghtInterp)) then if (useMassWghtInterp) massWeightToggle = 1. @@ -473,7 +481,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) - 0.25*real(n-1)*dz(i)) + p5(i*5+n) = -GxRho*((e(i,j,K) - z0pres) - 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) @@ -581,7 +589,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)) + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) ! Pressure do n=2,5 @@ -692,7 +700,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)) + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) ! Pressure do n=2,5 @@ -775,7 +783,7 @@ end subroutine int_density_dz_generic_plm !! are parabolic profiles subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, & - dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp) + dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -816,6 +824,7 @@ 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] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, 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 @@ -854,6 +863,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 [degC] real :: T_top, T_mn, T_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T real :: S_top, S_mn, S_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thicknes 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] @@ -868,6 +878,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. if (present(useMassWghtInterp)) then if (useMassWghtInterp) massWeightToggle = 1. @@ -900,7 +911,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(n) = -GxRho*(e(i,j,K) - 0.25*real(n-1)*dz) + p5(n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz) ! Salinity and temperature points are reconstructed with PPM S5(n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) ) T5(n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) ) @@ -978,7 +989,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! Pressure dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i+1,j,K) - e(i+1,j,K+1)) - p5(1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i+1,j,K)) + p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) do n=2,5 p5(n) = p5(n-1) + GxRho*0.25*dz enddo @@ -1066,7 +1077,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! Pressure dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i,j+1,K) - e(i,j+1,K+1)) - p5(1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i,j+1,K)) + p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) do n=2,5 p5(n) = p5(n-1) + GxRho*0.25*dz enddo diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 6e74c3ffa3..23f22d8a24 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1253,7 +1253,7 @@ end subroutine analytic_int_specific_vol_dp !! pressure anomalies across layers, which are required for calculating the !! finite-volume form pressure accelerations in a Boussinesq model. subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -1292,6 +1292,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] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, 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 ! desired units [R m3 kg-1 ~> 1] @@ -1322,11 +1324,11 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0)) then call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + dz_neglect, useMassWghtInterp, rho_scale, pres_scale, Z_0p=Z_0p) else call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp) + dz_neglect, useMassWghtInterp, Z_0p=Z_0p) endif case default call MOM_error(FATAL, "No analytic integration option is available with this EOS!") diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index e5cc9555b7..730687fbf6 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -408,7 +408,7 @@ end subroutine calculate_compress_wright !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -451,6 +451,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure !! into Pa [Pa T2 R-1 L-2 ~> 1]. + real, 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, p0_2d, lambda_2d @@ -461,7 +462,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: g_Earth ! The gravitational acceleration [m2 Z-1 s-2 ~> m s-2] real :: I_Rho ! The inverse of the Boussinesq density [m3 kg-1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] - real :: p_ave, I_al0, I_Lzz + real :: p_ave ! The layer averaged pressure [Pa] + real :: I_al0, I_Lzz real :: dz ! The layer thickness [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]. @@ -470,10 +472,11 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]. real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]. real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]. - real :: intz(5) ! The integrals of density with height at the - ! 5 sub-column locations [R L2 T-2 ~> Pa]. + real :: intz(5) ! The gravitational acceleration times the integrals of density + ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [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] or [1]. + real :: z0pres ! The height at which the pressure is zero [Z ~> m] logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants. real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants. @@ -499,6 +502,7 @@ 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 do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then @@ -517,7 +521,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 = -0.5*GxRho*(z_t(i,j)+z_b(i,j)) + p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -561,8 +565,7 @@ 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 = -0.5*GxRho*(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) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -603,8 +606,7 @@ 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 = -0.5*GxRho*(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) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)