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

Fix for using tidal amplitude when determining the BBL thickness #595

Merged
merged 2 commits into from
Apr 18, 2024
Merged

Fix for using tidal amplitude when determining the BBL thickness #595

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1910a20
Fix for using tidal amplitude in determining the BBL thickness
adcroft Apr 4, 2024
8734901
Merge branch 'dev/gfdl' into bbl-tideamp-fix
marshallward Apr 17, 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
75 changes: 47 additions & 28 deletions src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Open in desktop
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
Hallberg-NOAA marked this conversation as resolved.
Show resolved Hide resolved
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