Skip to content

Commit

Permalink
+Add optional Z_0p argument to int_density_dz
Browse files Browse the repository at this point in the history 
  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).
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Aug 15, 2021
1 parent db43b35 commit 4d8d7af
Show file tree
Hide file tree
Showing 3 changed files with 43 additions and 28 deletions.
43 changes: 27 additions & 16 deletions src/core/MOM_density_integrals.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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]
Expand All @@ -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]
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand All @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand All @@ -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.
Expand Down Expand Up @@ -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) )
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
8 changes: 5 additions & 3 deletions src/equation_of_state/MOM_EOS.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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!")
Expand Down
20 changes: 11 additions & 9 deletions src/equation_of_state/MOM_EOS_Wright.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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].
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 4d8d7af

Please sign in to comment.