Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into vertFPmix_units_in_comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Apr 23, 2024
2 parents 8ef4d9c + 0efe83d commit 5d52b29
Show file tree
Hide file tree
Showing 6 changed files with 147 additions and 109 deletions.
13 changes: 10 additions & 3 deletions .testing/tools/parse_perf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,16 @@ def parse_perf_report(perf_data_path):

# get per-symbol count
else:
tokens = line.split()
symbol = tokens[2]
period = int(tokens[3])
try:
tokens = line.split()
symbol = tokens[2]
period = int(tokens[3])
except ValueError:
print("parse_perf.py: Error extracting symbol count",
file=sys.stderr)
print("line:", repr(line), file=sys.stderr)
print("tokens:", tokens, file=sys.stderr)
raise

profile[event_name]['symbol'][symbol] = period

Expand Down
7 changes: 7 additions & 0 deletions src/parameterizations/vertical/MOM_bkgnd_mixing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,8 @@ module MOM_bkgnd_mixing
real :: N0_2Omega !< ratio of the typical Buoyancy frequency to
!! twice the Earth's rotation period, used with the
!! Henyey scaling from the mixing [nondim]
real :: Henyey_max_lat !< A latitude poleward of which the Henyey profile
!! is returned to the minimum diffusivity [degN]
real :: prandtl_bkgnd !< Turbulent Prandtl number used to convert
!! vertical background diffusivity into viscosity [nondim]
real :: Kd_tanh_lat_scale !< A nondimensional scaling for the range of
Expand Down Expand Up @@ -282,6 +284,10 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS, physical_OBL
call get_param(param_file, mdl, "OMEGA", CS%omega, &
"The rotation rate of the earth.", &
units="s-1", default=7.2921e-5, scale=US%T_to_s)
call get_param(param_file, mdl, "HENYEY_MAX_LAT", CS%Henyey_max_lat, &
"A latitude poleward of which the Henyey profile "//&
"is returned to the minimum diffusivity", &
units="degN", default=95.0)
endif

call get_param(param_file, mdl, "KD_TANH_LAT_FN", CS%Kd_tanh_lat_fn, &
Expand Down Expand Up @@ -447,6 +453,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G,
I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg.
do i=is,ie
abs_sinlat = abs(sin(G%geoLatT(i,j)*deg_to_rad))
if (abs(G%geoLatT(i,j))>CS%Henyey_max_lat) abs_sinlat = min_sinlat
Kd_sfc(i) = max(CS%Kd_min, CS%Kd * &
((abs_sinlat * invcosh(CS%N0_2Omega / max(min_sinlat, abs_sinlat))) * I_x30) )
enddo
Expand Down
75 changes: 47 additions & 28 deletions src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ module MOM_set_visc
real :: drag_bg_vel !< An assumed unresolved background velocity for
!! calculating the bottom drag [L T-1 ~> m s-1].
!! Runtime parameter `DRAG_BG_VEL`.
!! Should not be used if BBL_USE_TIDAL_BG is True.
real :: BBL_thick_min !< The minimum bottom boundary layer thickness [Z ~> m].
!! This might be Kv / (cdrag * drag_bg_vel) to give
!! Kv as the minimum near-bottom viscosity.
Expand Down Expand Up @@ -229,11 +230,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
real :: BBL_thick_max ! A huge upper bound on the boundary layer thickness [Z ~> m].
real :: kv_bbl ! The bottom boundary layer viscosity [H Z T-1 ~> m2 s-1 or Pa s]
real :: C2f ! C2f = 2*f at velocity points [T-1 ~> s-1].

real :: U_bg_sq ! The square of an assumed background
! velocity, for calculating the mean
! magnitude near the bottom for use in the
! quadratic bottom drag [L2 T-2 ~> m2 s-2].
real :: u2_bg(SZIB_(G)) ! The square of an assumed background velocity, for calculating the mean
! magnitude near the bottom for use in the quadratic bottom drag [L2 T-2 ~> m2 s-2].
real :: hwtot ! Sum of the thicknesses used to calculate
! the near-bottom velocity magnitude [H ~> m or kg m-2].
real :: I_hwtot ! The Adcroft reciprocal of hwtot [H-1 ~> m-1 or m2 kg-1].
Expand Down Expand Up @@ -338,7 +336,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
use_BBL_EOS = associated(tv%eqn_of_state) .and. CS%BBL_use_EOS
OBC => CS%OBC

U_bg_sq = CS%drag_bg_vel * CS%drag_bg_vel
cdrag_sqrt = sqrt(CS%cdrag)
cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H
Expand Down Expand Up @@ -422,7 +419,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)

!$OMP parallel do default(private) shared(u,v,h,dz,tv,visc,G,GV,US,CS,Rml,nz,nkmb,nkml,K2, &
!$OMP Isq,Ieq,Jsq,Jeq,h_neglect,dz_neglect,Rho0x400_G, &
!$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, &
!$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, &
!$OMP cdrag_L_to_H,cdrag_RL_to_H,use_BBL_EOS,BBL_thick_max, &
!$OMP OBC,D_u,D_v,mask_u,mask_v,pbv)
do j=Jsq,Jeq ; do m=1,2
Expand Down Expand Up @@ -591,6 +588,19 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
dztot_vel = 0.0 ; dzwtot = 0.0
Thtot = 0.0 ; Shtot = 0.0 ; SpV_htot = 0.0

! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
if (CS%BBL_use_tidal_bg) then
if (m==1) then
u2_bg(I) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) )
else
u2_bg(i) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) )
endif
else
u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel
endif
do k=nz,1,-1

if (htot_vel>=CS%Hbbl) exit ! terminate the k loop
Expand All @@ -606,18 +616,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)

if ((.not.CS%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then
v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC)
if (CS%BBL_use_tidal_bg) then
U_bg_sq = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) )
endif
hutot = hutot + hweight * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + U_bg_sq)
hutot = hutot + hweight * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + u2_bg(I))
else
u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC)
if (CS%BBL_use_tidal_bg) then
U_bg_sq = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) )
endif
hutot = hutot + hweight * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + U_bg_sq)
hutot = hutot + hweight * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + u2_bg(i))
endif ; endif

if (use_BBL_EOS .and. (hweight >= 0.0)) then
Expand Down Expand Up @@ -798,6 +800,19 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
if (m==1) then ; C2f = G%CoriolisBu(I,J-1) + G%CoriolisBu(I,J)
else ; C2f = G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J) ; endif

! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
if (CS%BBL_use_tidal_bg) then
if (m==1) then
u2_bg(I) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
G%mask2dT(i+1,j)*(CS%tideamp(i+1,j)*CS%tideamp(i+1,j)) )
else
u2_bg(i) = 0.5*( G%mask2dT(i,j)*(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
G%mask2dT(i,j+1)*(CS%tideamp(i,j+1)*CS%tideamp(i,j+1)) )
endif
else
u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel
endif

! The thickness of a rotation limited BBL ignoring stratification is
! h_f ~ Cn u* / f (limit of KW99 eq. 2.20 for N->0).
! The buoyancy limit of BBL thickness (h_N) is already in the variable htot from above.
Expand All @@ -809,7 +824,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
! xp = 1/2 + sqrt( 1/4 + (2 f h_N/u*)^2 )
! To avoid dividing by zero if u*=0 then
! xp u* = 1/2 u* + sqrt( 1/4 u*^2 + (2 f h_N)^2 )
if (CS%cdrag * U_bg_sq <= 0.0) then
if (CS%cdrag * u2_bg(i) <= 0.0) then
! This avoids NaNs and overflows, and could be used in all cases,
! but is not bitwise identical to the current code.
ustH = ustar(i) ; root = sqrt(0.25*ustH**2 + (htot*C2f)**2)
Expand Down Expand Up @@ -957,12 +972,12 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
if (m==1) then
if (Rayleigh > 0.0) then
v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC)
visc%Ray_u(I,j,k) = Rayleigh * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + U_bg_sq)
visc%Ray_u(I,j,k) = Rayleigh * sqrt(u(I,j,k)*u(I,j,k) + v_at_u*v_at_u + u2_bg(I))
else ; visc%Ray_u(I,j,k) = 0.0 ; endif
else
if (Rayleigh > 0.0) then
u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC)
visc%Ray_v(i,J,k) = Rayleigh * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + U_bg_sq)
visc%Ray_v(i,J,k) = Rayleigh * sqrt(v(i,J,k)*v(i,J,k) + u_at_v*u_at_v + u2_bg(i))
else ; visc%Ray_v(i,J,k) = 0.0 ; endif
endif

Expand Down Expand Up @@ -1992,9 +2007,9 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
real :: frac_used ! The fraction of the present layer that contributes to Dh and Ddz [nondim]
real :: Dh ! The increment in layer thickness from the present layer [H ~> m or kg m-2].
real :: Ddz ! The increment in height change from the present layer [Z ~> m].
real :: U_bg_sq ! The square of an assumed background velocity, for
! calculating the mean magnitude near the top for use in
! the quadratic surface drag [L2 T-2 ~> m2 s-2].
real :: u2_bg(SZIB_(G)) ! The square of an assumed background velocity, for
! calculating the mean magnitude near the top for use in
! the quadratic surface drag [L2 T-2 ~> m2 s-2].
real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than
! h_tiny can not be the deepest in the viscous mixed layer.
real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1].
Expand Down Expand Up @@ -2025,7 +2040,6 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
associated(forces%frac_shelf_v)) ) return

Rho0x400_G = 400.0*(GV%H_to_RZ / (US%L_to_Z**2 * GV%g_Earth))
U_bg_sq = CS%drag_bg_vel * CS%drag_bg_vel
cdrag_sqrt = sqrt(CS%cdrag)
cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H
Expand Down Expand Up @@ -2099,7 +2113,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)

!$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
!$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
!$OMP js,je,OBC,Isq,Ieq,nz,nkml,U_star_2d,U_bg_sq,mask_v, &
!$OMP js,je,OBC,Isq,Ieq,nz,nkml,U_star_2d,mask_v, &
!$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
do j=js,je ! u-point loop
if (CS%dynamic_viscous_ML) then
Expand Down Expand Up @@ -2251,7 +2265,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)

if (.not.CS%linear_drag) then
v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC)
hutot = hutot + hweight * sqrt(u(I,j,k)**2 + v_at_u**2 + U_bg_sq)
hutot = hutot + hweight * sqrt(u(I,j,k)**2 + v_at_u**2 + u2_bg(I))
endif
if (use_EOS) then
Thtot(I) = Thtot(I) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
Expand Down Expand Up @@ -2369,7 +2383,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)

!$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
!$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
!$OMP is,ie,OBC,Jsq,Jeq,nz,nkml,U_bg_sq,U_star_2d,mask_u, &
!$OMP is,ie,OBC,Jsq,Jeq,nz,nkml,U_star_2d,mask_u, &
!$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
do J=Jsq,Jeq ! v-point loop
if (CS%dynamic_viscous_ML) then
Expand Down Expand Up @@ -2523,7 +2537,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)

if (.not.CS%linear_drag) then
u_at_v = set_u_at_v(u, h, G, GV, i, J, k, mask_u, OBC)
hutot = hutot + hweight * sqrt(v(i,J,k)**2 + u_at_v**2 + U_bg_sq)
hutot = hutot + hweight * sqrt(v(i,J,k)**2 + u_at_v**2 + u2_bg(i))
endif
if (use_EOS) then
Thtot(i) = Thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
Expand Down Expand Up @@ -2967,6 +2981,11 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS
call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, &
"The name of the tidal amplitude variable in the input file.", &
default="tideamp")
! This value is here only to detect whether it is inadvertently used. CS%drag_bg_vel should
! not be used if CS%BBL_use_tidal_bg is True. For this reason, we do not apply dimensions,
! nor dimensional testing in this mode. If we ever detect a dimensional sensitivity to
! this parameter, in this mode, then it means it is being used inappropriately.
CS%drag_bg_vel = 1.e30
else
call get_param(param_file, mdl, "DRAG_BG_VEL", CS%drag_bg_vel, &
"DRAG_BG_VEL is either the assumed bottom velocity (with "//&
Expand Down
Loading

0 comments on commit 5d52b29

Please sign in to comment.