From f182f3dba92a55ab0976d3bb414d03532f08e330 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 29 May 2024 18:08:21 -0400 Subject: [PATCH] +Add ROBUST_STOKES_PGF and LA_MISALIGNMENT_BUG Added the new runtime parameters ROBUST_STOKES_PGF and LA_MISALIGNMENT_BUG to allow for the selection of Stokes pressure gradient force calculations that work properly in the limit of vanishingly thin layers and to correct a sign error in the calculations of the misalignment between the waves and shears in the Langmuir number calculations when LA_MISALIGNMENT is true. By default, all answers are bitwise identical, but as these options are not yet widely used it might make sense to move aggressively to obsolete the previous code once more extensive testing has take place. There will be new entries in the MOM_parameter_doc files for some cases, but there are not yet any such cases in the MOM-examples regression suite. --- src/user/MOM_wave_interface.F90 | 191 ++++++++++++++++++++++++-------- 1 file changed, 143 insertions(+), 48 deletions(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 656ff5b569..bcd66843ac 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -59,6 +59,9 @@ module MOM_wave_interface logical, public :: Stokes_VF = .false. !< True if Stokes vortex force is used logical, public :: Passive_Stokes_VF = .false. !< Computes Stokes VF, but doesn't affect dynamics logical, public :: Stokes_PGF = .false. !< True if Stokes shear pressure Gradient force is used + logical, public :: robust_Stokes_PGF = .false. !< If true, use expressions to calculate the + !! Stokes-induced pressure gradient anomalies that are + !! more accurate in the limit of thin layers. logical, public :: Passive_Stokes_PGF = .false. !< Keeps Stokes_PGF on, but doesn't affect dynamics logical, public :: Stokes_DDT = .false. !< Developmental: !! True if Stokes d/dt is used @@ -164,6 +167,8 @@ module MOM_wave_interface real :: LA_FracHBL !< Fraction of OSBL for averaging Langmuir number [nondim] real :: LA_HBL_min !< Minimum boundary layer depth for averaging Langmuir number [Z ~> m] logical :: LA_Misalignment = .false. !< Flag to use misalignment in Langmuir number + logical :: LA_misalign_bug = .false. !< Flag to use code with a sign error when calculating the + !! misalignment between the shear and waves in the Langmuir number calculation. real :: g_Earth !< The gravitational acceleration, equivalent to GV%g_Earth but with !! different dimensional rescaling appropriate for deep-water gravity !! waves [Z T-2 ~> m s-2] @@ -377,22 +382,27 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) call get_param(param_file, mdl, "STOKES_VF", CS%Stokes_VF, & "Flag to use Stokes vortex force", & - Default=.false.) + default=.false.) call get_param(param_file, mdl, "PASSIVE_STOKES_VF", CS%Passive_Stokes_VF, & "Flag to make Stokes vortex force diagnostic only.", & - Default=.false.) + default=.false.) call get_param(param_file, mdl, "STOKES_PGF", CS%Stokes_PGF, & "Flag to use Stokes-induced pressure gradient anomaly", & - Default=.false.) + default=.false.) + call get_param(param_file, mdl, "ROBUST_STOKES_PGF", CS%robust_Stokes_PGF, & + "If true, use expressions to calculate the Stokes-induced pressure gradient "//& + "anomalies that are more accurate in the limit of thin layers.", & + default=.false., do_not_log=.not.CS%Stokes_PGF) + !### Change the default for ROBUST_STOKES_PGF to True. call get_param(param_file, mdl, "PASSIVE_STOKES_PGF", CS%Passive_Stokes_PGF, & "Flag to make Stokes-induced pressure gradient anomaly diagnostic only.", & - Default=.false.) + default=.false.) call get_param(param_file, mdl, "STOKES_DDT", CS%Stokes_DDT, & "Flag to use Stokes d/dt", & - Default=.false.) + default=.false.) call get_param(param_file, mdl, "PASSIVE_STOKES_DDT", CS%Passive_Stokes_DDT, & "Flag to make Stokes d/dt diagnostic only", & - Default=.false.) + default=.false.) ! Get Wave Method and write to integer WaveMethod call get_param(param_file,mdl,"WAVE_METHOD",TMPSTRING1, & @@ -526,6 +536,11 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) call get_param(param_file, mdl, "LA_MISALIGNMENT", CS%LA_Misalignment, & "Flag (logical) if using misalignment between shear and waves in LA", & default=.false.) + call get_param(param_file, mdl, "LA_MISALIGNMENT_BUG", CS%LA_misalign_bug, & + "If true, use a code with a sign error when calculating the misalignment between "//& + "the shear and waves when LA_MISALIGNMENT is true.", & + default=CS%LA_Misalignment, do_not_log=.not.CS%LA_Misalignment) + !### Change the default for LA_MISALIGNMENT_BUG to .false. call get_param(param_file, mdl, "MIN_LANGMUIR", CS%La_min, & "A minimum value for all Langmuir numbers that is not physical, "//& "but is likely only encountered when the wind is very small and "//& @@ -1030,6 +1045,18 @@ real function one_minus_exp_x(x) endif end function one_minus_exp_x +!> Return the value of (1 - exp(-x)), using an accurate expression for small values of x. +real function one_minus_exp(x) + real, intent(in) :: x !< The argument of the function ((1 - exp(-x))/x) [nondim] + real, parameter :: C1_6 = 1.0/6.0 ! A rational fraction [nondim] + if (abs(x) <= 2.0e-5) then + ! The Taylor series expression for exp(-x) gives a more accurate expression for 64-bit reals. + one_minus_exp = x * (1.0 - x * (0.5 - C1_6*x)) + else + one_minus_exp = 1.0 - exp(-x) + endif +end function one_minus_exp + !> A subroutine to fill the Stokes drift from a NetCDF file !! using the data_override procedures. subroutine Surface_Bands_by_data_override(Time, G, GV, US, CS) @@ -1199,12 +1226,18 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, dz, Waves, & Top = Bottom MidPoint = Bottom + 0.5*dz(k) Bottom = Bottom + dz(k) - !### Given the sign convention that Dpt_LASL is negative, the next line seems to have a bug. - ! To correct this bug, this line should be changed to: - ! if (MidPoint > abs(Dpt_LASL) .and. (k > 1) .and. ContinueLoop) then - if (MidPoint > Dpt_LASL .and. k > 1 .and. ContinueLoop) then - ShearDirection = atan2(V_H(1)-V_H(k), U_H(1)-U_H(k)) - ContinueLoop = .false. + + if (Waves%LA_Misalign_bug) then + ! Given the sign convention that Dpt_LASL is negative, the next line has a bug. + if (MidPoint > Dpt_LASL .and. k > 1 .and. ContinueLoop) then + ShearDirection = atan2(V_H(1)-V_H(k), U_H(1)-U_H(k)) + ContinueLoop = .false. + endif + else ! This version avoids the bug in the version above. + if (MidPoint > abs(Dpt_LASL) .and. (k > 1) .and. ContinueLoop) then + ShearDirection = atan2(V_H(1)-V_H(k), U_H(1)-U_H(k)) + ContinueLoop = .false. + endif endif enddo endif @@ -1706,7 +1739,9 @@ subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) real :: P_Stokes_l0, P_Stokes_r0 ! Stokes-induced pressure anomaly at interface ! (left/right of point) [L2 T-2 ~> m2 s-2] real :: dP_Stokes_l_dz, dP_Stokes_r_dz ! Contribution of layer to integrated Stokes pressure anomaly for summation - ! (left/right of point) [L3 T-2 ~> m3 s-2] + ! (left/right of point) [Z L2 T-2 ~> m3 s-2] + real :: dP_lay_Stokes_l, dP_lay_Stokes_r ! Contribution of layer to integrated Stokes pressure anomaly for summation + ! (left/right of point) [L2 T-2 ~> m2 s-2] real :: dP_Stokes_l, dP_Stokes_r ! Net increment of Stokes pressure anomaly across layer for summation ! (left/right of point) [L2 T-2 ~> m2 s-2] real :: uE_l, uE_r, vE_l, vE_r ! Eulerian velocity components (left/right of point) [L T-1 ~> m s-1] @@ -1714,6 +1749,7 @@ subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) real :: zi_l(SZK_(G)+1), zi_r(SZK_(G)+1) ! The height of the edges of the cells (left/right of point) [Z ~> m]. real :: idz_l(SZK_(G)), idz_r(SZK_(G)) ! The inverse thickness of the cells (left/right of point) [Z-1 ~> m-1] real :: h_l, h_r ! The thickness of the cell (left/right of point) [Z ~> m]. + real :: exp_top ! The decay of the surface stokes drift to the interface atop a layer [nondim] real :: dexp2kzL, dexp4kzL, dexp2kzR, dexp4kzR ! Analytical evaluation of multi-exponential decay ! contribution to Stokes pressure anomalies [nondim]. real :: TwoK, FourK ! Wavenumbers multiplied by a factor [Z-1 ~> m-1] @@ -1762,9 +1798,11 @@ subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) h_r = dz(i+1,j,k) zi_l(k+1) = zi_l(k) - h_l zi_r(k+1) = zi_r(k) - h_r - !### If the code were properly refactored, the following hard-coded constants would be unnecessary. - Idz_l(k) = 1./max(0.1*US%m_to_Z, h_l) - Idz_r(k) = 1./max(0.1*US%m_to_Z, h_r) + if (.not.CS%robust_Stokes_PGF) then + ! When the code is properly refactored, the following hard-coded constants are unnecessary. + Idz_l(k) = 1./max(0.1*US%m_to_Z, h_l) + Idz_r(k) = 1./max(0.1*US%m_to_Z, h_r) + endif enddo do k = 1,G%ke ! Computing (left/right) Eulerian velocities assuming the velocity passed to this routine is the @@ -1798,31 +1836,59 @@ subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) ! Wavenumber terms that are useful to simplify the pressure calculations TwoK = 2.*CS%WaveNum_Cen(l) FourK = 2.*TwoK - iTwoK = 1./TwoK - iFourK = 1./(FourK) - dexp2kzL = exp(TwoK*zi_l(k))-exp(TwoK*zi_l(k+1)) - dexp2kzR = exp(TwoK*zi_r(k))-exp(TwoK*zi_r(k+1)) - dexp4kzL = exp(FourK*zi_l(k))-exp(FourK*zi_l(k+1)) - dexp4kzR = exp(FourK*zi_r(k))-exp(FourK*zi_r(k+1)) + if (.not.CS%robust_Stokes_PGF) then + iTwoK = 1. / TwoK + iFourK = 1. / FourK + endif ! Compute Pressure at interface and integrated over layer on left/right bounding points. ! These are summed over wavenumber bands. if (G%mask2dT(i,j)>0.5) then - dP_Stokes_l_dz = dP_Stokes_l_dz + & - ((uE_l*uS0_l+vE_l*vS0_l)*iTwoK*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzL) - dP_Stokes_l = dP_Stokes_l + (uE_l*uS0_l+vE_l*vS0_l)*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzL + if (.not.CS%robust_Stokes_PGF) then + dexp2kzL = exp(TwoK*zi_l(k))-exp(TwoK*zi_l(k+1)) + dexp4kzL = exp(FourK*zi_l(k))-exp(FourK*zi_l(k+1)) + dP_Stokes_l_dz = dP_Stokes_l_dz + & + ((uE_l*uS0_l+vE_l*vS0_l)*iTwoK*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzL) + dP_Stokes_l = dP_Stokes_l + (uE_l*uS0_l+vE_l*vS0_l)*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzL + else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers. + exp_top = exp(TwoK*zi_l(k)) + dP_lay_Stokes_l = dP_lay_Stokes_l + & + ((((uE_l*uS0_l)+(vE_l*vS0_l)) * exp_top) * one_minus_exp_x(TwoK*dz(i,j,k)) + & + (0.5*((uS0_l**2)+(vS0_l**2)) * exp_top**2) * one_minus_exp_x(FourK*dz(i,j,k)) ) + dP_Stokes_l = dP_Stokes_l + & + ((((uE_l*uS0_l)+(vE_l*vS0_l)) * exp_top) * one_minus_exp(TwoK*dz(i,j,k)) + & + (0.5*((uS0_l**2)+(vS0_l**2)) * exp_top**2) * one_minus_exp(FourK*dz(i,j,k)) ) + endif endif if (G%mask2dT(i+1,j)>0.5) then - dP_Stokes_r_dz = dP_Stokes_r_dz + & - ((uE_r*uS0_r+vE_r*vS0_r)*iTwoK*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzR) - dP_Stokes_r = dP_Stokes_r + (uE_r*uS0_r+vE_r*vS0_r)*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzR + if (.not.CS%robust_Stokes_PGF) then + dexp2kzR = exp(TwoK*zi_r(k))-exp(TwoK*zi_r(k+1)) + dexp4kzR = exp(FourK*zi_r(k))-exp(FourK*zi_r(k+1)) + dP_Stokes_r_dz = dP_Stokes_r_dz + & + ((uE_r*uS0_r+vE_r*vS0_r)*iTwoK*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzR) + dP_Stokes_r = dP_Stokes_r + (uE_r*uS0_r+vE_r*vS0_r)*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzR + else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers. + exp_top = exp(TwoK*zi_r(k)) + dP_lay_Stokes_r = dP_lay_Stokes_r + & + ((((uE_r*uS0_r)+(vE_r*vS0_r)) * exp_top) * one_minus_exp_x(TwoK*dz(i+1,j,k)) + & + (0.5*((uS0_r**2)+(vS0_r**2)) * exp_top**2) * one_minus_exp_x(FourK*dz(i+1,j,k)) ) + dP_Stokes_r = dP_Stokes_r + & + ((((uE_r*uS0_r)+(vE_r*vS0_r)) * exp_top) * one_minus_exp(TwoK*dz(i+1,j,k)) + & + (0.5*((uS0_r**2)+(vS0_r**2)) * exp_top**2) * one_minus_exp(FourK*dz(i+1,j,k)) ) + endif endif enddo ! Summing PF over bands ! > Increment the Layer averaged pressure - P_Stokes_l = P_Stokes_l0 + dP_Stokes_l_dz*Idz_l(k) - P_Stokes_r = P_Stokes_r0 + dP_Stokes_r_dz*Idz_r(k) + if (.not.CS%robust_Stokes_PGF) then + P_Stokes_l = P_Stokes_l0 + dP_Stokes_l_dz*Idz_l(k) + P_Stokes_r = P_Stokes_r0 + dP_Stokes_r_dz*Idz_r(k) + else + P_Stokes_l = P_Stokes_l0 + dP_lay_Stokes_l + P_Stokes_r = P_Stokes_r0 + dP_lay_Stokes_r + endif + ! > Increment the Interface pressure P_Stokes_l0 = P_Stokes_l0 + dP_Stokes_l P_Stokes_r0 = P_Stokes_r0 + dP_Stokes_r @@ -1856,9 +1922,11 @@ subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) h_r = dz(i,j+1,k) zi_l(k+1) = zi_l(k) - h_l zi_r(k+1) = zi_r(k) - h_r - !### If the code were properly refactored, the following hard-coded constants would be unnecessary. - Idz_l(k) = 1. / max(0.1*US%m_to_Z, h_l) - Idz_r(k) = 1. / max(0.1*US%m_to_Z, h_r) + if (.not.CS%robust_Stokes_PGF) then + ! When the code is properly refactored, the following hard-coded constants are unnecessary. + Idz_l(k) = 1. / max(0.1*US%m_to_Z, h_l) + Idz_r(k) = 1. / max(0.1*US%m_to_Z, h_r) + endif enddo do k = 1,G%ke ! Computing (left/right) Eulerian velocities assuming the velocity passed to this routine is the @@ -1892,31 +1960,59 @@ subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) ! Wavenumber terms that are useful to simplify the pressure calculations TwoK = 2.*CS%WaveNum_Cen(l) FourK = 2.*TwoK - iTwoK = 1./TwoK - iFourK = 1./(FourK) - dexp2kzL = exp(TwoK*zi_l(k))-exp(TwoK*zi_l(k+1)) - dexp2kzR = exp(TwoK*zi_r(k))-exp(TwoK*zi_r(k+1)) - dexp4kzL = exp(FourK*zi_l(k))-exp(FourK*zi_l(k+1)) - dexp4kzR = exp(FourK*zi_r(k))-exp(FourK*zi_r(k+1)) + if (.not.CS%robust_Stokes_PGF) then + iTwoK = 1. / TwoK + iFourK = 1. / FourK + endif ! Compute Pressure at interface and integrated over layer on left/right bounding points. ! These are summed over wavenumber bands. if (G%mask2dT(i,j)>0.5) then - dP_Stokes_l_dz = dP_Stokes_l_dz + & - ((uE_l*uS0_l+vE_l*vS0_l)*iTwoK*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzL) - dP_Stokes_l = dP_Stokes_l + (uE_l*uS0_l+vE_l*vS0_l)*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzL + if (.not.CS%robust_Stokes_PGF) then + dexp2kzL = exp(TwoK*zi_l(k))-exp(TwoK*zi_l(k+1)) + dexp4kzL = exp(FourK*zi_l(k))-exp(FourK*zi_l(k+1)) + dP_Stokes_l_dz = dP_Stokes_l_dz + & + ((uE_l*uS0_l+vE_l*vS0_l)*iTwoK*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzL) + dP_Stokes_l = dP_Stokes_l + (uE_l*uS0_l+vE_l*vS0_l)*dexp2kzL + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzL + else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers. + exp_top = exp(TwoK*zi_l(k)) + dP_lay_Stokes_l = dP_lay_Stokes_l + & + ((((uE_l*uS0_l)+(vE_l*vS0_l)) * exp_top) * one_minus_exp_x(TwoK*dz(i,j,k)) + & + (0.5*((uS0_l**2)+(vS0_l**2)) * exp_top**2) * one_minus_exp_x(FourK*dz(i,j,k)) ) + dP_Stokes_l = dP_Stokes_l + & + ((((uE_l*uS0_l)+(vE_l*vS0_l)) * exp_top) * one_minus_exp(TwoK*dz(i,j,k)) + & + (0.5*((uS0_l**2)+(vS0_l**2)) * exp_top**2) * one_minus_exp(FourK*dz(i,j,k)) ) + endif endif if (G%mask2dT(i,j+1)>0.5) then - dP_Stokes_r_dz = dP_Stokes_r_dz + & - ((uE_r*uS0_r+vE_r*vS0_r)*iTwoK*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzR) - dP_Stokes_r = dP_Stokes_r + (uE_r*uS0_r+vE_r*vS0_r)*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzR + if (.not.CS%robust_Stokes_PGF) then + dexp2kzR = exp(TwoK*zi_r(k))-exp(TwoK*zi_r(k+1)) + dexp4kzR = exp(FourK*zi_r(k))-exp(FourK*zi_r(k+1)) + dP_Stokes_r_dz = dP_Stokes_r_dz + & + ((uE_r*uS0_r+vE_r*vS0_r)*iTwoK*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*iFourK*dexp4kzR) + dP_Stokes_r = dP_Stokes_r + (uE_r*uS0_r+vE_r*vS0_r)*dexp2kzR + 0.5*(uS0_l*uS0_l+vS0_l*vS0_l)*dexp4kzR + else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers. + exp_top = exp(TwoK*zi_r(k)) + dP_lay_Stokes_r = dP_lay_Stokes_r + & + ((((uE_r*uS0_r)+(vE_r*vS0_r)) * exp_top) * one_minus_exp_x(TwoK*dz(i,j+1,k)) + & + (0.5*((uS0_r**2)+(vS0_r**2)) * exp_top**2) * one_minus_exp_x(FourK*dz(i,j+1,k)) ) + dP_Stokes_r = dP_Stokes_r + & + ((((uE_r*uS0_r)+(vE_r*vS0_r)) * exp_top) * one_minus_exp(TwoK*dz(i,j+1,k)) + & + (0.5*((uS0_r**2)+(vS0_r**2)) * exp_top**2) * one_minus_exp(FourK*dz(i,j+1,k)) ) + endif endif enddo ! Summing PF over bands ! > Increment the Layer averaged pressure - P_Stokes_l = P_Stokes_l0 + dP_Stokes_l_dz*Idz_l(k) - P_Stokes_r = P_Stokes_r0 + dP_Stokes_r_dz*Idz_r(k) + if (.not.CS%robust_Stokes_PGF) then + P_Stokes_l = P_Stokes_l0 + dP_Stokes_l_dz*Idz_l(k) + P_Stokes_r = P_Stokes_r0 + dP_Stokes_r_dz*Idz_r(k) + else + P_Stokes_l = P_Stokes_l0 + dP_lay_Stokes_l + P_Stokes_r = P_Stokes_r0 + dP_lay_Stokes_r + endif + ! > Increment the Interface pressure P_Stokes_l0 = P_Stokes_l0 + dP_Stokes_l P_Stokes_r0 = P_Stokes_r0 + dP_Stokes_r @@ -1939,7 +2035,6 @@ subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) end subroutine Stokes_PGF - !> Computes wind speed from ustar_air based on COARE 3.5 Cd relationship !! Probably doesn't belong in this module, but it is used here to estimate !! wind speed for wind-wave relationships. Should be a fine way to estimate