Skip to content

Commit

Permalink
Fix for using tidal amplitude in determining the BBL thickness
Browse files Browse the repository at this point in the history 
Setting `BBL_USE_TIDAL_BG` in OM5_b003 revealed non-conservation and
non-reproducibility across layouts. A scalar variable was correctly
being set based on the tidal amplitude but was being used after a break
between loops.

- Replaced the scalar variable U_bg_sq with a vector u2_bg(:) that
  is set to either a constant or the square of the tidal amplitude.
- The module parameter CS%drag_bg_vel was not set if using the tidal
  amplitude but WAS being used for conditionals. We now set it to
  1e30 because it should NOT be used or impact the solution. Setting
  it to zero, or anything meaningful would allow code to use it
  inadvertently. This "constant" does not need to satisfy any
  dimensional testing.
  • Loading branch information
@adcroft
adcroft committed Apr 4, 2024
1 parent f9372f3 commit d71ed45
Showing 1 changed file with 42 additions and 28 deletions.
70 changes: 42 additions & 28 deletions src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -229,11 +229,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 +335,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 +418,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 +587,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 +615,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 +799,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 +823,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 +971,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 +2006,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 +2039,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 +2112,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 +2264,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 +2382,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 +2536,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 +2980,7 @@ 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")
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

0 comments on commit d71ed45

Please sign in to comment.