Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into f2023
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Sep 17, 2024
2 parents 6b632ac + 0363d2b commit 6bed957
Show file tree
Hide file tree
Showing 10 changed files with 1,165 additions and 251 deletions.
6 changes: 6 additions & 0 deletions .testing/tc1/MOM_input
View file
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
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
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

Large diffs are not rendered by default.

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

Large diffs are not rendered by default.

32 changes: 18 additions & 14 deletions src/equation_of_state/MOM_EOS.F90
View file
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
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

0 comments on commit 6bed957

Please sign in to comment.