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

Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting in PGF #810

Merged
merged 3 commits into from
Feb 13, 2025
Merged

Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting in PGF #810

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9944ebc
Add MASS_WEIGHT_IN_PGF_VANISHED_ONLY to modify mass weighting
claireyung Jan 21, 2025
29dd0b2
Minor style updates to previous commit of Add MASS_WEIGHT_IN_PGF_VANI…
claireyung Jan 21, 2025
6630441
Address comments from Bob by adding scaling to dimensional numbers, r…
claireyung Feb 6, 2025
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
32 changes: 25 additions & 7 deletions src/core/MOM_PressureForce_FV.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ module MOM_PressureForce_FV
logical :: reset_intxpa_integral !< If true and the surface displacement between adjacent cells
!! exceeds the vertical grid spacing, reset intxpa at the interface below
!! a trusted interior cell. (This often applies in ice shelf cavities.)
logical :: MassWghtInterpVanOnly !< If true, don't do mass weighting of T/S interpolation unless vanished
real :: h_nonvanished !< A minimal layer thickness that indicates that a layer is thick enough
!! to usefully reestimate the pressure integral across the interface
!! below it [H ~> m or kg m-2]
Expand Down Expand Up @@ -224,6 +225,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_

real :: dp_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [R L2 T-2 ~> Pa].
real :: p_nonvanished ! nonvanshed pressure [R L2 T-2 ~> Pa]
real :: I_gEarth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1]
real :: alpha_anom ! The in-situ specific volume, averaged over a
! layer, less alpha_ref [R-1 ~> m3 kg-1].
Expand Down Expand Up @@ -278,6 +280,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
dp_neglect = GV%g_Earth*GV%H_to_RZ * GV%H_subroundoff
alpha_ref = 1.0 / CS%Rho0
I_gEarth = 1.0 / GV%g_Earth
p_nonvanished = GV%g_Earth*GV%H_to_RZ*CS%h_nonvanished

if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then
MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0
Expand Down Expand Up @@ -354,7 +357,8 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
call int_spec_vol_dp_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), &
p(:,:,K), p(:,:,K+1), alpha_ref, dp_neglect, p(:,:,nz+1), G%HI, &
tv%eqn_of_state, US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), &
P_surf=p(:,:,1), MassWghtInterp=CS%MassWghtInterp)
P_surf=p(:,:,1), MassWghtInterp=CS%MassWghtInterp, &
MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, p_nv=p_nonvanished)
elseif ( CS%Recon_Scheme == 2 ) then
call MOM_error(FATAL, "PressureForce_FV_nonBouss: "//&
"int_spec_vol_dp_generic_ppm does not exist yet.")
Expand All @@ -368,11 +372,13 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, &
US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
inty_dza(:,:,k), bathyP=p(:,:,nz+1), P_surf=p(:,:,1), dP_tiny=dp_neglect, &
MassWghtInterp=CS%MassWghtInterp)
MassWghtInterp=CS%MassWghtInterp, &
MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, p_nv=p_nonvanished)
endif
if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) &
call diagnose_mass_weight_p(p(:,:,K), p(:,:,K+1), p(:,:,nz+1), p(:,:,1), dp_neglect, CS%MassWghtInterp, &
G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k))
G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k), &
MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, p_nv=p_nonvanished)
else
alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand Down Expand Up @@ -979,6 +985,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
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 :: dz_nonvanished ! A small thickness considered to be vanished for mass weighting [Z ~> m]
real :: H_to_RL2_T2 ! A factor to convert from thickness units (H) to pressure
! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].
real :: T5(5) ! Temperatures and salinities at five quadrature points [C ~> degC]
Expand Down Expand Up @@ -1019,6 +1026,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm

h_neglect = GV%H_subroundoff
dz_neglect = GV%dZ_subroundoff
dz_nonvanished = GV%H_to_Z*CS%h_nonvanished
I_Rho0 = 1.0 / GV%Rho0
G_Rho0 = GV%g_Earth / GV%Rho0
GxRho = GV%g_Earth * GV%Rho0
Expand Down Expand Up @@ -1173,19 +1181,22 @@ 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=Z_0p)
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=Z_0p, &
MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=dz_nonvanished)
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=Z_0p)
MassWghtInterp=CS%MassWghtInterp, Z_0p=Z_0p, &
MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=dz_nonvanished)
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, e(:,:,1), dz_neglect, &
CS%MassWghtInterp, Z_0p=Z_0p)
CS%MassWghtInterp, Z_0p=Z_0p, &
MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=dz_nonvanished)
endif
if (GV%Z_to_H /= 1.0) then
!$OMP parallel do default(shared)
Expand All @@ -1195,7 +1206,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
endif
if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) &
call diagnose_mass_weight_Z(e(:,:,K), e(:,:,K+1), G%bathyT, e(:,:,1), dz_neglect, CS%MassWghtInterp, &
G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k))
G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k), &
MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=CS%h_nonvanished)
else
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand Down Expand Up @@ -1806,6 +1818,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
logical :: useMassWghtInterp ! If true, use near-bottom mass weighting for T and S
logical :: MassWghtInterpTop ! If true, use near-surface mass weighting for T and S under ice shelves
logical :: MassWghtInterp_NonBous_bug ! If true, use a buggy mass weighting when non-Boussinesq
logical :: MassWghtInterpVanOnly ! If true, turn of mass weighting unless one side is vanished
! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=40) :: mdl ! This module's name.
Expand Down Expand Up @@ -1888,6 +1901,11 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
"when interpolating T/S for integrals near the bathymetry in FV pressure "//&
"gradient calculations.", &
default=.false., do_not_log=(GV%Boussinesq .or. (.not.useMassWghtInterp)))
call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_VANISHED_ONLY", CS%MassWghtInterpVanOnly, &
"If true, use mass weighting when interpolating T/S for integrals "//&
"only if one side is vanished according to RESET_INTXPA_H_NONVANISHED. ", &
default=.false.)

CS%MassWghtInterp = 0
if (useMassWghtInterp) &
CS%MassWghtInterp = ibset(CS%MassWghtInterp, 0) ! Same as CS%MassWghtInterp + 1
Expand Down
Loading