+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

@@ -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,14 +1739,17 @@ 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]
   real :: uS0_l, uS0_r, vS0_l, vS0_r ! Surface Stokes velocity components (left/right of point) [L T-1 ~> m s-1]
   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