From 1910a2012235dcdbdbb12eadba101e03e6709968 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 4 Apr 2024 16:12:37 -0400 Subject: [PATCH] Fix for using tidal amplitude in determining the BBL thickness 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. --- .../vertical/MOM_set_viscosity.F90 | 75 ++++++++++++------- 1 file changed, 47 insertions(+), 28 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index e601fdf2f7..e16a57d970 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -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. @@ -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]. @@ -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 @@ -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 @@ -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 @@ -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 @@ -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. @@ -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) @@ -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 @@ -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]. @@ -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 @@ -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 @@ -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)) @@ -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 @@ -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)) @@ -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 "//&