Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

+(*)Pressure gradient corrections under ice shelves #713

Merged
merged 11 commits into from
Sep 17, 2024
Merged

+(*)Pressure gradient corrections under ice shelves #713

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
ecb365e
+Add MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP
claireyung Apr 30, 2024
4b9a8c0
+Add top mass_weight_in_PGF option to 13 integrals
Hallberg-NOAA Jul 24, 2024
94e4db0
+Add CORRECTION_INTXPA
claireyung Apr 30, 2024
e0870a8
+Add CORRECTION_INTXPA_5PT
claireyung Apr 30, 2024
37418fc
+Add RESET_INTXPA_INTEGRAL
claireyung Apr 30, 2024
3c9bdfe
Revisions of sub-ice pressure gradient fixes
Hallberg-NOAA Jul 19, 2024
65d9988
*Refactor CORRECTION_INTX_PA
Hallberg-NOAA Jul 29, 2024
c6ec948
*Non-Boussinesq code for RESET_INTXPA_INTEGRAL
Hallberg-NOAA Aug 5, 2024
3a7e4a1
+(*)Add 5-point quadrature in RESET_INTXPA_INTEGRAL
Hallberg-NOAA Aug 20, 2024
64628f7
+(*)Eliminate CORRECTION_INTXPA_5PT
Hallberg-NOAA Sep 16, 2024
c69581e
*Set MASS_WEIGHT_IN_PRESSURE_GRADIENT in .testing
Hallberg-NOAA Sep 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions .testing/tc1/MOM_input
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -278,6 +278,12 @@ BOUND_CORIOLIS = True ! [Boolean] default = False
! have no effect on the SADOURNY Coriolis scheme if it
! were possible to use centered difference thickness fluxes.

! === module MOM_PressureForce_FV ===
MASS_WEIGHT_IN_PRESSURE_GRADIENT = True ! [Boolean] default = False
! If true, use mass weighting when interpolating T/S for integrals
! near the bathymetry in FV pressure gradient calculations.


! === module MOM_hor_visc ===
AH_VEL_SCALE = 0.05 ! [m s-1] default = 0.0
! The velocity scale which is multiplied by the cube of
Expand Down
7 changes: 6 additions & 1 deletion .testing/tc2/MOM_input
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -302,11 +302,16 @@ PGF_STANLEY_T2_DET_COEFF = -1.0 ! [nondim] default = -1.0
! gradient in the deterministic part of the Stanley form of the Brankart
! correction. Negative values disable the scheme.

! === module MOM_PressureForce_FV ===
MASS_WEIGHT_IN_PRESSURE_GRADIENT = True ! [Boolean] default = False
! If true, use mass weighting when interpolating T/S for integrals
! near the bathymetry in FV pressure gradient calculations.

! === module MOM_hor_visc ===
LAPLACIAN = True
KH_VEL_SCALE = 0.05
SMAGORINSKY_KH = True ! [Boolean] default = False
SMAG_LAP_CONST = 0.06 ! [nondim] default = 0.0
SMAG_LAP_CONST = 0.06 ! [nondim] default = 0.0
AH_VEL_SCALE = 0.05 ! [m s-1] default = 0.0
! The velocity scale which is multiplied by the cube of
! the grid spacing to calculate the Laplacian viscosity.
Expand Down
3 changes: 3 additions & 0 deletions .testing/tc4/MOM_input
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -269,6 +269,9 @@ BOUND_CORIOLIS = True ! [Boolean] default = False
! === module MOM_PressureForce ===

! === module MOM_PressureForce_FV ===
MASS_WEIGHT_IN_PRESSURE_GRADIENT = True ! [Boolean] default = False
! If true, use mass weighting when interpolating T/S for integrals
! near the bathymetry in FV pressure gradient calculations.
RECONSTRUCT_FOR_PRESSURE = False ! [Boolean] default = True
! If True, use vertical reconstruction of T & S within the integrals of the FV
! pressure gradient calculation. If False, use the constant-by-layer algorithm.
Expand Down
859 changes: 805 additions & 54 deletions src/core/MOM_PressureForce_FV.F90
View file
Open in desktop

Large diffs are not rendered by default.

325 changes: 219 additions & 106 deletions src/core/MOM_density_integrals.F90
View file
Open in desktop

Large diffs are not rendered by default.

32 changes: 18 additions & 14 deletions src/equation_of_state/MOM_EOS.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1226,7 +1226,7 @@ end function EOS_domain
!! series for log(1-eps/1+eps) that assumes that |eps| < 0.34.
subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, &
dza, intp_dza, intx_dza, inty_dza, halo_size, &
bathyP, dP_tiny, MassWghtInterp)
bathyP, P_surf, dP_tiny, MassWghtInterp)
type(hor_index_type), intent(in) :: HI !< The horizontal index structure
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC]
Expand Down Expand Up @@ -1259,6 +1259,8 @@ subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, &
integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza.
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
real, optional, intent(in) :: dP_tiny !< A miniscule pressure change with
!! the same units as p_t [R L2 T-2 ~> Pa]
integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
Expand All @@ -1280,20 +1282,20 @@ subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, &
call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%kg_m3_to_R*EOS%Rho_T0_S0, &
dRdT_scale*EOS%dRho_dT, dRdS_scale*EOS%dRho_dS, dza, &
intp_dza, intx_dza, inty_dza, halo_size, &
bathyP, dP_tiny, MassWghtInterp)
bathyP, P_surf, dP_tiny, MassWghtInterp)
case (EOS_WRIGHT)
call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, &
inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, &
inty_dza, halo_size, bathyP, P_surf, dP_tiny, MassWghtInterp, &
SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, &
temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt)
case (EOS_WRIGHT_FULL)
call int_spec_vol_dp_wright_full(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, &
inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, &
inty_dza, halo_size, bathyP, P_surf, dP_tiny, MassWghtInterp, &
SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, &
temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt)
case (EOS_WRIGHT_REDUCED)
call int_spec_vol_dp_wright_red(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, &
inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, &
inty_dza, halo_size, bathyP, P_surf, dP_tiny, MassWghtInterp, &
SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, &
temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt)
case default
Expand All @@ -1306,7 +1308,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, MassWghtInterp, Z_0p)
intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, 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 [C ~> degC]
Expand Down Expand Up @@ -1342,6 +1344,8 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS,
!! layer divided by the y grid spacing [R L2 T-2 ~> Pa]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: SSH !< The sea surface height [Z ~> m]
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
Expand All @@ -1367,49 +1371,49 @@ 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. (dRdT_scale /= 1.0) .or. (dRdS_scale /= 1.0)) then
call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
rho_scale*EOS%Rho_T0_S0, dRdT_scale*EOS%dRho_dT, dRdS_scale*EOS%dRho_dS, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp)
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp)
else
call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp)
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp)
endif
case (EOS_WRIGHT)
rho_scale = EOS%kg_m3_to_R
pres_scale = EOS%RL2_T2_to_Pa
if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 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, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &
dz_neglect, MassWghtInterp, rho_scale, pres_scale, &
temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, 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, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &
dz_neglect, MassWghtInterp, Z_0p=Z_0p)
endif
case (EOS_WRIGHT_FULL)
rho_scale = EOS%kg_m3_to_R
pres_scale = EOS%RL2_T2_to_Pa
if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then
call int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &
dz_neglect, MassWghtInterp, rho_scale, pres_scale, &
temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p)
else
call int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &
dz_neglect, MassWghtInterp, Z_0p=Z_0p)
endif
case (EOS_WRIGHT_REDUCED)
rho_scale = EOS%kg_m3_to_R
pres_scale = EOS%RL2_T2_to_Pa
if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then
call int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &
dz_neglect, MassWghtInterp, rho_scale, pres_scale, &
temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p)
else
call int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, &
dz_neglect, MassWghtInterp, Z_0p=Z_0p)
endif
case default
Expand Down
46 changes: 27 additions & 19 deletions src/equation_of_state/MOM_EOS_Wright.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -387,7 +387,7 @@ end subroutine EoS_fit_range_buggy_Wright
!! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps)
!! that assumes that |eps| < 0.34.
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, &
dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, &
MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_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), &
Expand Down Expand Up @@ -424,6 +424,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
!! layer divided by the y grid spacing [R L2 T-2 ~> Pa].
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m].
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: SSH !< The sea surface height [Z ~> m]
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
Expand Down Expand Up @@ -481,6 +483,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1]
real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1]
logical :: do_massWeight ! Indicates whether to do mass weighting.
logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim]
real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim]
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m
Expand Down Expand Up @@ -531,14 +534,11 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
c4s = c4s * saln_scale ; c5s = c5s * saln_scale
endif ; endif

do_massWeight = .false.
if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
! if (do_massWeight) then
! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//&
! "bathyT must be present if MassWghtInterp is present and true.")
! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//&
! "dz_neglect must be present if MassWghtInterp is present and true.")
! endif
do_massWeight = .false. ; top_massWeight = .false.
if (present(MassWghtInterp)) then
do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
endif

do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
al0_2d(i,j) = (a0 + a1s*T(i,j)) + a2s*S(i,j)
Expand Down Expand Up @@ -571,6 +571,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hWght = 0.0
if (do_massWeight) &
hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j))
if (top_massWeight) &
hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j))
if (hWght > 0.) then
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
Expand Down Expand Up @@ -613,6 +615,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hWght = 0.0
if (do_massWeight) &
hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j))
if (top_massWeight) &
hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1))
if (hWght > 0.) then
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
Expand Down Expand Up @@ -656,7 +660,7 @@ end subroutine int_density_dz_wright
!! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps)
!! that assumes that |eps| < 0.34.
subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, &
intp_dza, intx_dza, inty_dza, halo_size, bathyP, P_surf, dP_neglect, &
MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale)
type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type.
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
Expand Down Expand Up @@ -693,6 +697,8 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
!! dza.
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa]
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with
!! the same units as p_t [R L2 T-2 ~> Pa]
integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use
Expand Down Expand Up @@ -743,6 +749,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1]
real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1]
logical :: do_massWeight ! Indicates whether to do mass weighting.
logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface
logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting
real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim]
real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim]
Expand Down Expand Up @@ -780,15 +787,12 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
c4s = c4s * saln_scale ; c5s = c5s * saln_scale
endif ; endif

do_massWeight = .false. ; massWeight_bug = .false.
if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set
! if (do_massWeight) then
! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//&
! "bathyP must be present if MassWghtInterp is present and true.")
! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//&
! "dP_neglect must be present if MassWghtInterp is present and true.")
! endif
do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false.
if (present(MassWghtInterp)) then
do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values
top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set
massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set
endif

! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0)
do j=jsh,jeh ; do i=ish,ieh
Expand Down Expand Up @@ -818,6 +822,8 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
elseif (do_massWeight) then
hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j))
endif
if (top_massWeight) &
hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j))
if (hWght > 0.) then
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
Expand Down Expand Up @@ -862,6 +868,8 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
elseif (do_massWeight) then
hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1))
endif
if (top_massWeight) &
hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j))
if (hWght > 0.) then
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
Expand Down
Loading
Loading