diff --git a/.github/workflows/intel.yml b/.github/workflows/intel.yml
index 700553cea..1159e8599 100644
--- a/.github/workflows/intel.yml
+++ b/.github/workflows/intel.yml
@@ -1,3 +1,8 @@
+# This is a GitHub actions workflow for WW3.
+#
+# This workflow builds with the Intel compilers.
+#
+# Matt Masarik, Alex Richert, Ed Hartnett
 name: Intel Linux Build
 on: [push, pull_request, workflow_dispatch]
 
@@ -8,7 +13,7 @@ concurrency:
 
 # Set I_MPI_CC/F90 so Intel MPI wrapper uses icc/ifort instead of gcc/gfortran
 env:
-  cache_key: intel10
+  cache_key: intel12
   CC: icc
   FC: ifort
   CXX: icpc
@@ -21,7 +26,7 @@ env:
 
 jobs:
   setup:
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-latest
 
     steps:
 
@@ -51,12 +56,13 @@ jobs:
           sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB
           echo "deb https://apt.repos.intel.com/oneapi all main" | sudo tee /etc/apt/sources.list.d/oneAPI.list
           sudo apt-get update
-          sudo apt-get install intel-oneapi-dev-utilities intel-oneapi-mpi-devel intel-oneapi-openmp intel-oneapi-compiler-fortran intel-oneapi-compiler-dpcpp-cpp-and-cpp-classic
+          sudo apt-get install intel-oneapi-dev-utilities intel-oneapi-mpi-devel intel-oneapi-compiler-fortran-2023.2.1 intel-oneapi-compiler-dpcpp-cpp-and-cpp-classic-2023.2.1 intel-oneapi-openmp
 
       # Build WW3 spack environment
       - name: install-dependencies-with-spack
         if: steps.cache-env.outputs.cache-hit != 'true'
         run: |
+          sudo mv /usr/local /usr/local_mv
           # Install NetCDF, ESMF, g2, etc using Spack
           . /opt/intel/oneapi/setvars.sh
           git clone -c feature.manyFiles=true https://github.com/JCSDA/spack.git
@@ -67,7 +73,7 @@ jobs:
           spack compiler find
           sudo apt install cmake
           spack external find
-          spack add intel-oneapi-mpi
+          spack config add "packages:mpi:require:'intel-oneapi-mpi'"
           spack config add "packages:all:require:['%intel']"
           spack concretize
           spack install --dirty -v --fail-fast
@@ -92,7 +98,7 @@ jobs:
     strategy:
       matrix:
         switch: [Ifremer1, NCEP_st2, NCEP_st4, ite_pdlib, NCEP_st4sbs, NCEP_glwu, OASACM, UKMO, MULTI_ESMF]
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-latest
 
     steps:
       - name: checkout-ww3
@@ -113,6 +119,8 @@ jobs:
 
       - name: build-ww3
         run: |
+          sudo mv /usr/local /usr/local_mv
+          sudo apt install cmake
           . /opt/intel/oneapi/setvars.sh
           source spack/share/spack/setup-env.sh
           spack env activate ww3-intel
diff --git a/model/inp/ww3_grid.inp b/model/inp/ww3_grid.inp
index b802c6717..55b75a598 100644
--- a/model/inp/ww3_grid.inp
+++ b/model/inp/ww3_grid.inp
@@ -305,25 +305,25 @@ $   limitation and the GSE alleviation.
 $
 $ Unstructured grids ------------------------------------------------ $
 $   UNST parameters       : Namelist UNST
-$                           UGBCCFL    : Turns on/off (TRUE/FALSE) the computation of the CFL number on the physical domain boundary. 
+$                           UGBCCFL    : Turns on/off (TRUE/FALSE) the computation of the CFL number on the physical domain boundary.
 $                                        If FALSE the explicit scheme can be much faster though stability is not guaranteed
 $                                        (default TRUE)
 $                           UGOBCAUTO  : TRUE: OBC points are taken from  type 15 elements (default)
 $                                        FALSE: OBC points must be listed in ww3_grid.inp
 $                           UGOBCDEPTH : Threshold ( < 0) depth (default -10) for OBC points if UGOBCAUTO is TRUE
-$                           UGOBCFILE  : File name of file for reading boudary (default 'unset') 
+$                           UGOBCFILE  : File name of file for reading boudary (default 'unset')
 $
-$                           The following are TRUE/FALSE variables and only one can be TRUE. 
-$                           By default, EXPFSN is TRUE and must be set to false to activate another option  
+$                           The following are TRUE/FALSE variables and only one can be TRUE.
+$                           By default, EXPFSN is TRUE and must be set to false to activate another option
 $                           EXPFSN    : Activation of N scheme (default option)
 $                           EXPFSPSI  : Activation of PSI scheme
 $                           EXPFSFCT  : Activation of FCT scheme
 $                           IMPFSN    : Activation of N implicit scheme
 $                           EXPTOTAL  : Activation of Block explicit N scheme solver
 $                           IMPTOTAL  : Activation of fully implicit scheme, non splitting
-$                             
+$
 $                           The following TRUE/FALSE variables are only for IMPTOTAL=TRUE
-$                           IMPREFRACTION : Turn on implicit freq. shift (default FALSE) 
+$                           IMPREFRACTION : Turn on implicit freq. shift (default FALSE)
 $                           IMPFREQSHIFT  : Turn on implicit freq. shift terms (default FALSE)
 $                           IMPSOURCE     : Turn on implicit source terms (default FALSE)
 $
@@ -332,17 +332,17 @@ $                           JGS_TERMINATE_DIFFERENCE  : Terminate based on the t
 $                           JGS_TERMINATE_NORM  : Terminate based on the norm of the solution (TRUE/FALSE, default FALSE)
 $                           JGS_USE_JACOBI  : Use Jacobi solver family (TRUE/FALSE, default TRUE)
 $                           JGS_BLOCK_GAUSS_SEIDEL  : Use Block Gauss Seidel method for imptotal instead of the conservative jacobi iterator. (TRUE/FALSE, default TRUE)
-$                           JGS_MAXITER  : Max. Number of solver iterations for JGS_TERMINATE_MAXITER (integer, default 100) 
-$                           JGS_PMIN  : % of grid points that do not need to converge during solver iteration (real, default 1) 
+$                           JGS_MAXITER  : Max. Number of solver iterations for JGS_TERMINATE_MAXITER (integer, default 100)
+$                           JGS_PMIN  : % of grid points that do not need to converge during solver iteration (real, default 1)
 $                           JGS_DIFF_THR  : Implicit solver threshold for JGS_TERMINATE_DIFFERENCE (real, default 1.0e-10)
 $                           JGS_NORM_THR  : Norm of the solution for JGS_TERMINATE_NORM (real, default 1.0e-20)
 $                           JGS_LIMITER  : TRUE: Use total (quasi-steady: limits whole equation) instead of local limiter (un-steady: limits only source terms)
-$                                          FALSE: default 
+$                                          FALSE: default
 $                           JGS_LIMITER_FUNC  : 1 - old limiter (default)
 $                                               2 - alternatnive limiter
-$                           SETUP_APPLY_WLV  : Compute wave setup (TRUE/FALSE, default TRUE) 
+$                           SETUP_APPLY_WLV  : Compute wave setup (TRUE/FALSE, default TRUE)
 $                           SOLVERTHR_SETUP  : Solver threshold for setup computations (default 1E-6)
-$                           CRIT_DEP_SETUP  : Critical depth for setup computations (default 0.1) 
+$                           CRIT_DEP_SETUP  : Critical depth for setup computations (default 0.1)
 
 $
 $ SMC grid propagation    : Namelist PSMC and default values
@@ -434,6 +434,33 @@ $                                 provides 3 partitions of both x & y component,
 $                                 with a reasonable accuracy for using in
 $                                 a climate model.
 $
+$ Langmuir turbulence parameterization ------------------------------- $
+$   The LMPN namelist group configures the Langmuir Mixing
+$   Parameterization used within the context of climate modeling and
+$   when WW3 is coupled to an active ocean model (.e.g., POP2, MOM6).
+$   The LMPENABLED flag turns on the parameterization computations.
+$   The SDTAIL flag enables the spectral tail contribution. HSLMODE
+$   parameter determines the surface layer depth (HSL).
+$
+$   Input parameters  : Namelist LMPN
+$                       LMPENABLED : Logical controlling Langmuir
+$                                    related computations.
+$                                     T : Activates calculations
+$                                     F : Calculations off (default)
+$                       SDTAIL     : Logical controlling addition of
+$                                    tail contribution.
+$                                     T : Tail added to Stokes
+$                                         drift (LMPENABLED=T also).
+$                                     F : No tail added (default)
+$                       HSLMODE    : Integer indicating surface layer
+$                                    depth (HSL) type.
+$                                     0 : Uniform HSL=10m (default)
+$                                     1 : HSL from coupler
+$
+$   For more info on the Langmuir parameterization, see:
+$       Li, Qing, et al. "Langmuir mixing effects on global climate:
+$       WAVEWATCH III in CESM." Ocean Modelling 103 (2016): 145-160.
+$
 $ Miscellaneous ------------------------------------------------------ $
 $   Misc. parameters    : Namelist MISC
 $                           CICE0  : Ice concentration cut-off.
diff --git a/model/inp/ww3_ounf.inp b/model/inp/ww3_ounf.inp
index 7bde30754..8505a4233 100644
--- a/model/inp/ww3_ounf.inp
+++ b/model/inp/ww3_ounf.inp
@@ -15,7 +15,7 @@ $
 $ DPT CUR WND AST WLV ICE IBG TAU RHO D50 IC1 IC5 HS LM T02 T0M1 T01 FP
 $ DIR SPR DP HIG EF TH1M STH1M TH2M STH2M WN PHS PTP PLP PDIR PSPR PWS PDP
 $ PQP PPE PGW PSW PTM10 PT01 PT02 PEP TWS PNR UST CHA CGE FAW TAW TWA WCC
-$ WCF WCH WCM SXY TWO BHD FOC TUS USS P2S USF P2L TWI FIC ABR UBR BED
+$ WCF WCH WCM SXY TWO BHD FOC TUS USS USSH P2S USF P2L TWI FIC ABR UBR BED
 $ FBB TBB MSS MSC DTD FC CFX CFD CFK U1 U2 WNM TOC
 $
  N
diff --git a/model/inp/ww3_shel.inp b/model/inp/ww3_shel.inp
index 8eca9b58a..1fc2220c1 100644
--- a/model/inp/ww3_shel.inp
+++ b/model/inp/ww3_shel.inp
@@ -190,7 +190,7 @@ $  F  F  6     2   TAUO[X,Y]  TWO   Wave to ocean momentum flux
 $  F  F  6     3   BHD        BHD   Bernoulli head (J term)
 $  F  F  6     4   PHIOC      FOC   Wave to ocean energy flux
 $  F  F  6     5   TUS[X,Y]   TUS   Stokes transport
-$  F  F  6     6   USS[X,Y]   USS   Surface Stokes drift
+$  F  F  6     6   USS[X,Y]   USS   Surface layer Stokes drift
 $  F  F  6     7   [PR,TP]MS  P2S   Second-order sum pressure
 $  F  F  6     8   US3D       USF   Spectrum of surface Stokes drift
 $  F  F  6     9   P2SMS      P2L   Micro seism  source term
@@ -198,6 +198,7 @@ $  F  F  6    10   TAUICE     TWI   Wave to sea ice stress
 $  F  F  6    11   PHICE      FIC   Wave to sea ice energy flux
 $  F  F  6    12   USSP       USP   Partitioned surface Stokes drift
 $  F  F  6    13   TAUOC[X,Y] TOC   Total momentum to the ocean
+$  F  F  6    14   USSH[X,Y]  USSH  Surface layer averaged Stokes drift
 $   -------------------------------------------------
 $        7                          Wave-bottom layer
 $   -------------------------------------------------
@@ -257,8 +258,8 @@ $ UST CHA CGE FAW TAW TWA WCC WCF WCH WCM FWS
   T   T   T   T   T   T   T   T   T   T   T
 $ (6) Wave-Ocean layer
   T
-$ SXY TWO BHD FOC TUS USS P2S USF P2L TWI FIC USP
-  T   T   T   T   T   T   T   F   F   F   F   T
+$ SXY TWO BHD FOC TUS USS P2S USF P2L TWI FIC USP TOC USSH
+  T   T   T   T   T   T   T   F   F   F   F   T   T   T
 $ (7) Wave-bottom layer
   T
 $ ABR UBR BED FBB TBB
diff --git a/model/src/w3adatmd.F90 b/model/src/w3adatmd.F90
index a4aea27d0..a48af3199 100644
--- a/model/src/w3adatmd.F90
+++ b/model/src/w3adatmd.F90
@@ -174,6 +174,8 @@ MODULE W3ADATMD
   !      US3D      R.A.  Public   3D Stokes drift.
   !      USSP      R.A.  Public   Partitioned Surface Stokes drift
   !
+  !      USSHX/Y   R.A.  Public   Surface layer averaged Stokes drift.
+  !
   !      ABA       R.A.  Public   Near-bottom rms wave ex. amplitude.
   !      ABD       R.A.  Public   Corresponding direction.
   !      UBA       R.A.  Public   Near-bottom rms wave velocity.
@@ -461,9 +463,7 @@ MODULE W3ADATMD
          XPRMS(:), XTPMS(:), XPHICE(:),       &
          XTAUICE(:,:)
     REAL, POINTER         :: XP2SMS(:,:), XUS3D(:,:), XUSSP(:,:)
-#ifdef W3_CESMCOUPLED
-    REAL, POINTER         ::  XLANGMT(:)
-#endif
+    REAL, POINTER         :: XUSSHX(:), XUSSHY(:)
     !
     ! Output fields group 7)
     !
@@ -492,12 +492,8 @@ MODULE W3ADATMD
     !
     REAL, POINTER         ::  USERO(:,:)
     REAL, POINTER         :: XUSERO(:,:)
-#ifdef W3_CESMCOUPLED
-    ! Output fileds for Langmuir mixing in group
-    REAL, POINTER         :: LANGMT(:), LAPROJ(:), LASL(:),       &
-         LASLPJ(:), LAMULT(:), ALPHAL(:),     &
-         ALPHALS(:), USSXH(:), USSYH(:)
-#endif
+    ! Output fileds for Langmuir mixing parameterization
+    REAL, POINTER         :: USSHX(:), USSHY(:)
     !
     ! Spatial derivatives
     !
@@ -578,11 +574,8 @@ MODULE W3ADATMD
   !/
   !/ Data aliases for structure WADAT(S)
   !/
-#ifdef W3_CESMCOUPLED
-  REAL, POINTER           :: LANGMT(:), LAPROJ(:), ALPHAL(:),      &
-       ALPHALS(:), LAMULT(:), LASL(:),       &
-       LASLPJ(:), USSXH(:), USSYH(:)
-#endif
+  REAL, POINTER :: USSHX(:), USSHY(:)
+  !
   REAL, POINTER           :: CG(:,:), WN(:,:)
   REAL, POINTER           :: IC3WN_R(:,:), IC3WN_I(:,:), IC3CG(:,:)
   !
@@ -1064,19 +1057,6 @@ SUBROUTINE W3DIMA  ( IMOD, NDSE, NDST, D_ONLY )
          STAT=ISTAT )
     CHECK_ALLOC_STATUS ( ISTAT )
 
-#ifdef W3_CESMCOUPLED
-    ALLOCATE ( WADATS(IMOD)%USSXH(NSEALM)   , &
-         WADATS(IMOD)%USSYH(NSEALM)   , &
-         WADATS(IMOD)%LANGMT(NSEALM)  , &
-         WADATS(IMOD)%LAPROJ(NSEALM)  , &
-         WADATS(IMOD)%LASL(NSEALM)    , &
-         WADATS(IMOD)%LASLPJ(NSEALM)  , &
-         WADATS(IMOD)%ALPHAL(NSEALM)  , &
-         WADATS(IMOD)%ALPHALS(NSEALM) , &
-         WADATS(IMOD)%LAMULT(NSEALM)  , &
-         STAT=ISTAT )
-    CHECK_ALLOC_STATUS ( ISTAT )
-#endif
     !
     WADATS(IMOD)%HS     = UNDEF
     WADATS(IMOD)%WLM    = UNDEF
@@ -1226,6 +1206,8 @@ SUBROUTINE W3DIMA  ( IMOD, NDSE, NDST, D_ONLY )
          WADATS(IMOD)%TPMS  (NSEALM) ,                        &
          WADATS(IMOD)%PHICE (NSEALM) ,                        &
          WADATS(IMOD)%TAUICE(NSEALM,2),                       &
+         WADATS(IMOD)%USSHX(NSEALM),                          &
+         WADATS(IMOD)%USSHY(NSEALM),                          &
          STAT=ISTAT )
     CHECK_ALLOC_STATUS ( ISTAT )
     !
@@ -1262,9 +1244,8 @@ SUBROUTINE W3DIMA  ( IMOD, NDSE, NDST, D_ONLY )
     WADATS(IMOD)%TPMS   = UNDEF
     WADATS(IMOD)%PHICE  = UNDEF
     WADATS(IMOD)%TAUICE = UNDEF
-#ifdef W3_CESMCOUPLED
-    WADATS(IMOD)%LANGMT = UNDEF
-#endif
+    WADATS(IMOD)%USSHX  = UNDEF
+    WADATS(IMOD)%USSHY  = UNDEF
     IF (  P2MSF(1).GT.0 ) WADATS(IMOD)%P2SMS  = UNDEF
     IF (  US3DF(1).GT.0 ) WADATS(IMOD)%US3D   = UNDEF
     IF (  USSPF(1).GT.0 ) WADATS(IMOD)%USSP   = UNDEF
@@ -2180,15 +2161,17 @@ SUBROUTINE W3XDMA  ( IMOD, NDSE, NDST, OUTFLAGS )
       ALLOCATE ( WADATS(IMOD)%XTAUOCY(1), STAT=ISTAT )
       CHECK_ALLOC_STATUS ( ISTAT )
     END IF
-#ifdef W3_CESMCOUPLED
     IF ( OUTFLAGS( 6, 14) ) THEN
-      ALLOCATE ( WADATS(IMOD)%XLANGMT(NXXX), STAT=ISTAT )
+      ALLOCATE ( WADATS(IMOD)%XUSSHX(NXXX), STAT=ISTAT )
+      CHECK_ALLOC_STATUS ( ISTAT )
+      ALLOCATE ( WADATS(IMOD)%XUSSHY(NXXX), STAT=ISTAT )
       CHECK_ALLOC_STATUS ( ISTAT )
     ELSE
-      ALLOCATE ( WADATS(IMOD)%XLANGMT(1), STAT=ISTAT )
+      ALLOCATE ( WADATS(IMOD)%XUSSHX(1), STAT=ISTAT )
+      CHECK_ALLOC_STATUS ( ISTAT )
+      ALLOCATE ( WADATS(IMOD)%XUSSHY(1), STAT=ISTAT )
       CHECK_ALLOC_STATUS ( ISTAT )
     END IF
-#endif
     !
     WADATS(IMOD)%XSXX    = UNDEF
     WADATS(IMOD)%XSYY    = UNDEF
@@ -2210,9 +2193,8 @@ SUBROUTINE W3XDMA  ( IMOD, NDSE, NDST, OUTFLAGS )
     WADATS(IMOD)%XUSSP   = UNDEF
     WADATS(IMOD)%XTAUOCX = UNDEF
     WADATS(IMOD)%XTAUOCY = UNDEF
-#ifdef W3_CESMCOUPLED
-    WADATS(IMOD)%XLANGMT = UNDEF
-#endif
+    WADATS(IMOD)%XUSSHX   = UNDEF
+    WADATS(IMOD)%XUSSHY   = UNDEF
     !
     IF ( OUTFLAGS( 7, 1) ) THEN
       ALLOCATE ( WADATS(IMOD)%XABA(NXXX), STAT=ISTAT )
@@ -2946,18 +2928,8 @@ SUBROUTINE W3SETA ( IMOD, NDSE, NDST )
       USERO  => WADATS(IMOD)%USERO
       !
       WN     => WADATS(IMOD)%WN
-#ifdef W3_CESMCOUPLED
-      ! USSX and USSY are already set
-      LANGMT => WADATS(IMOD)%LANGMT
-      LAPROJ => WADATS(IMOD)%LAPROJ
-      LASL   => WADATS(IMOD)%LASL
-      LASLPJ => WADATS(IMOD)%LASLPJ
-      ALPHAL => WADATS(IMOD)%ALPHAL
-      ALPHALS=> WADATS(IMOD)%ALPHALS
-      USSXH  => WADATS(IMOD)%USSXH
-      USSYH  => WADATS(IMOD)%USSYH
-      LAMULT => WADATS(IMOD)%LAMULT
-#endif
+      USSHX  => WADATS(IMOD)%USSHX
+      USSHY  => WADATS(IMOD)%USSHY
 #ifdef W3_IC3
       IC3WN_R=> WADATS(IMOD)%IC3WN_R
       IC3WN_I=> WADATS(IMOD)%IC3WN_I
@@ -3279,9 +3251,6 @@ SUBROUTINE W3XETA ( IMOD, NDSE, NDST )
       BEDFORMS=> WADATS(IMOD)%XBEDFORMS
       PHIBBL => WADATS(IMOD)%XPHIBBL
       TAUBBL => WADATS(IMOD)%XTAUBBL
-#ifdef W3_CESMCOUPLED
-      LANGMT => WADATS(IMOD)%XLANGMT
-#endif
       !
       MSSX   => WADATS(IMOD)%XMSSX
       MSSY   => WADATS(IMOD)%XMSSY
@@ -3298,6 +3267,9 @@ SUBROUTINE W3XETA ( IMOD, NDSE, NDST )
       !
       USERO  => WADATS(IMOD)%XUSERO
       !
+      USSHX   => WADATS(IMOD)%XUSSHX
+      USSHY   => WADATS(IMOD)%XUSSHY
+      !
     END IF
     !
     RETURN
diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90
index 6cd6e91d8..add5db979 100644
--- a/model/src/w3gdatmd.F90
+++ b/model/src/w3gdatmd.F90
@@ -728,6 +728,9 @@ MODULE W3GDATMD
 #ifdef W3_IS2
     REAL,    POINTER      :: IS2PARS(:)
 #endif
+    LOGICAL               :: LMPENABLED ! flag to enable Li et al. Langmuir parameterization
+    LOGICAL               :: SDTAIL ! flag to enable high-freq tail in Li et al. Stokes Drift computations
+    INTEGER               :: HSLMODE ! 0 for test (HSL=10m everywhere, 1 for coupler-based HSL)
     !
     ! unstructured data
     !
@@ -1084,6 +1087,10 @@ MODULE W3GDATMD
   INTEGER, POINTER        :: NBEDGE
   INTEGER, POINTER        :: EDGES(:,:), NEIGH(:,:)
   !
+  LOGICAL, POINTER        :: LMPENABLED
+  LOGICAL, POINTER        :: SDTAIL
+  INTEGER, POINTER        :: HSLMODE
+  !
   ! Variables for unstructured grids
   !
   INTEGER, POINTER        :: NTRI,COUNTRI,COUNTOT,NNZ
@@ -2274,6 +2281,10 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     USSPF  => GRIDS(IMOD)%USSPF
     USSP_WN => GRIDS(IMOD)%USSP_WN
     FFACBERG => GRIDS(IMOD)%FFACBERG
+    !
+    LMPENABLED => GRIDS(IMOD)%LMPENABLED
+    SDTAIL => GRIDS(IMOD)%SDTAIL
+    HSLMODE => GRIDS(IMOD)%HSLMODE
 #ifdef W3_REF1
     REFLC  => GRIDS(IMOD)%REFLC
     REFLD  => GRIDS(IMOD)%REFLD
diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90
index 281ed7a3f..052b1f2c9 100644
--- a/model/src/w3gridmd.F90
+++ b/model/src/w3gridmd.F90
@@ -1105,6 +1105,7 @@ MODULE W3GRIDMD
        STH1MF, I1STH1M, I2STH1M,                &
        TH2MF, I1TH2M, I2TH2M,                   &
        STH2MF, I1STH2M, I2STH2M
+  NAMELIST /LMPN/ LMPENABLED, SDTAIL, HSLMODE
 #ifdef W3_IS1
   NAMELIST /SIS1/ ISC1, ISC2
 #endif
@@ -2731,6 +2732,11 @@ SUBROUTINE W3GRID()
     I2STH2M=NK
     !
     FACBERG=1.
+    !
+    LMPENABLED = .false.
+    SDTAIL = .false.
+    HSLMODE = 0   ! 0 for test (HSL=10m everywhere, 1 for coupler-based HSL)
+    !
 #ifdef W3_IS0
     WRITE (NDSO,944)
 #endif
@@ -2882,6 +2888,10 @@ SUBROUTINE W3GRID()
          IC5MAXKI, IC5MINHW, IC5MAXITER, IC5RKICK, &
          IC5KFILTER, IC5MSTR(NINT(IC5VEMOD))
 #endif
+    !
+    CALL READNL ( NDSS, 'LMPN', STATUS )
+    WRITE (NDSO,4960) STATUS
+    WRITE (NDSO,4961) LMPENABLED, SDTAIL, HSLMODE
     !
     CALL READNL ( NDSS, 'OUTS', STATUS )
     WRITE (NDSO,4970) STATUS
@@ -6673,7 +6683,11 @@ SUBROUTINE W3GRID()
          /'         (0.0==> no reduction and 1.0==> no wind', &
          /'         input with 100% ice cover)')
     !
-    !
+!
+4960 FORMAT (/'  Langmuir Mixing Parameterization ',A/                   &
+              ' --------------------------------------------------')
+4961 FORMAT ('  &LMPN LMPENABLED = ',L, 'SDTAIL = ', L, ' HSLMODE = ', I2 '/' )
+!
 4970 FORMAT (/'  Spectral output on full grid ',A/                   &
          ' --------------------------------------------------')
 4971 FORMAT ( '       Second order pressure at K=0:',3I4)
@@ -7344,6 +7358,8 @@ SUBROUTINE READNL ( NDS, NAME, STATUS )
                 READ (NDS,NML=UNST,END=801,ERR=802,IOSTAT=J)
               CASE('OUTS')
                 READ (NDS,NML=OUTS,END=801,ERR=802,IOSTAT=J)
+              CASE('LMPN')
+                READ (NDS,NML=LMPN,END=801,ERR=802,IOSTAT=J)
               CASE('MISC')
                 READ (NDS,NML=MISC,END=801,ERR=802,IOSTAT=J)
               CASE DEFAULT
diff --git a/model/src/w3idatmd.F90 b/model/src/w3idatmd.F90
index 83d9af3ee..52035dafd 100644
--- a/model/src/w3idatmd.F90
+++ b/model/src/w3idatmd.F90
@@ -94,9 +94,8 @@ MODULE W3IDATMD
   !      FLCUR     Log.  Public   Flag for current input.
   !      FLWIND    Log.  Public   Flag for wind input.
   !      FLICE     Log.  Public   Flag for ice input.
-#ifdef W3_CESMCOUPLED
-  !      HML       R.A.  Public   Mixed layer depth
-#endif
+  !      HSL       R.A.  Public   Depth of a surface layer over which Stokes
+  !                               drift is averaged
   !      FLTAUA    Log.  Public   Flag for atmospheric momentum input
   !      FLRHOA    Log.  Public   Flag for air density input
   !      INFLAGS1  L.A.  Public   Array consolidating the above six
@@ -219,9 +218,7 @@ MODULE W3IDATMD
     REAL, POINTER         :: CYTIDE(:,:,:,:)
     REAL, POINTER         :: WLTIDE(:,:,:,:)
 #endif
-#ifdef W3_CESMCOUPLED
-    REAL, POINTER         :: HML(:,:)
-#endif
+    REAL, POINTER         :: HSL(:,:)
     LOGICAL               :: IINIT
 #ifdef W3_WRST
     LOGICAL               :: WRSTIINIT=.FALSE.
@@ -272,9 +269,7 @@ MODULE W3IDATMD
   LOGICAL, POINTER        ::  FLLEVTIDE, FLCURTIDE,  &
        FLLEVRESI, FLCURRESI
 #endif
-#ifdef W3_CESMCOUPLED
-  REAL   , POINTER        :: HML(:,:)
-#endif
+  REAL , POINTER :: HSL(:,:)
   !/
 CONTAINS
   !/ ------------------------------------------------------------------- /
@@ -743,10 +738,8 @@ SUBROUTINE W3DIMI  ( IMOD, NDSE, NDST, FLAGSTIDEIN )
       CHECK_ALLOC_STATUS ( ISTAT )
     END IF
     !
-#ifdef W3_CESMCOUPLED
-    ALLOCATE ( INPUTS(IMOD)%HML(NX,NY), STAT=ISTAT )
+    ALLOCATE ( INPUTS(IMOD)%HSL(NX,NY), STAT=ISTAT )
     CHECK_ALLOC_STATUS ( ISTAT )
-#endif
     !
     INPUTS(IMOD)%IINIT  = .TRUE.
     !
@@ -1061,9 +1054,7 @@ SUBROUTINE W3SETI ( IMOD, NDSE, NDST )
         ICEI   => INPUTS(IMOD)%ICEI
         BERGI  => INPUTS(IMOD)%BERGI
       END IF
-#ifdef W3_CESMCOUPLED
-      HML    => INPUTS(IMOD)%HML
-#endif
+      HSL => INPUTS(IMOD)%HSL
       !
       IF ( FLTAUA  ) THEN
         UX0    => INPUTS(IMOD)%UX0
diff --git a/model/src/w3initmd.F90 b/model/src/w3initmd.F90
index ce7e2ddf8..9f7c62de8 100644
--- a/model/src/w3initmd.F90
+++ b/model/src/w3initmd.F90
@@ -2152,10 +2152,7 @@ SUBROUTINE W3MPIO ( IMOD )
          TAUOCX, TAUOCY, WNMEAN
 #endif
 
-#ifdef W3_CESMCOUPLED
-    USE W3ADATMD, ONLY: LANGMT, LAPROJ, ALPHAL, LASL, LASLPJ,  &
-         ALPHALS, LAMULT
-#endif
+    USE W3ADATMD, ONLY: USSHX, USSHY
 
 #ifdef W3_MPI
     USE W3GDATMD, ONLY: NK
@@ -2248,7 +2245,7 @@ SUBROUTINE W3MPIO ( IMOD )
            0                 +    0  +    0  +  &  ! group 3 (extra contributions below)
            2+(NOGE(4)-2)*(NOSWLL+1) +    0  +    0  +  &  ! group 4
            11                +    3  +    1  +  &  ! group 5
-           12                +    7  +    1  +  &  ! group 6 (extra contributions below)
+           10                +    7  +    1  +  &  ! group 6 (extra contributions below)
            5                 +    4  +    1  +  &  ! group 7
            5                 +    2  +    0  +  &  ! group 8
            5                 +    0  +    0  +  &  ! group 9
@@ -2262,6 +2259,7 @@ SUBROUTINE W3MPIO ( IMOD )
       IF ( FLGRDALL( 6,9)) NRQMAX = NRQMAX + P2MSF(3) - P2MSF(2) + 1
       IF ( FLGRDALL( 6, 8) ) NRQMAX = NRQMAX + 2*NK
       IF ( FLGRDALL( 6,12) ) NRQMAX = NRQMAX + 2*NK
+      IF ( FLGRDALL( 6,14) ) NRQMAX = NRQMAX + 2
       !
       IF ( NRQMAX .GT. 0 ) THEN
         ALLOCATE ( OUTPTS(IMOD)%OUT1%IRQGO(NRQMAX) )
@@ -3218,17 +3216,22 @@ SUBROUTINE W3MPIO ( IMOD )
 #ifdef W3_MPI
         END IF
         !
-#ifdef W3_CESMCOUPLED
         IF ( FLGRDALL( 6, 14) ) THEN
           IH     = IH + 1
           IT     = IT + 1
-          CALL MPI_SEND_INIT (LANGMT(1),NSEALM , MPI_REAL, IROOT,   &
+          CALL MPI_SEND_INIT (USSHX (1),NSEALM , MPI_REAL, IROOT,   &
+               IT, MPI_COMM_WAVE, IRQGO(IH), IERR)
+#ifdef W3_MPIT
+          WRITE (NDST,9011) IH, ' 6/14', IROOT, IT, IRQGO(IH), IERR
+#endif
+          IH     = IH + 1
+          IT     = IT + 1
+          CALL MPI_SEND_INIT (USSHY (1),NSEALM , MPI_REAL, IROOT,   &
                IT, MPI_COMM_WAVE, IRQGO(IH), IERR)
 #ifdef W3_MPIT
           WRITE (NDST,9011) IH, ' 6/14', IROOT, IT, IRQGO(IH), IERR
 #endif
         END IF
-#endif !W3_CESMCOUPLED
         IF ( FLGRDALL( 7, 1) ) THEN
           IH     = IH + 1
           IT     = IT + 1
@@ -4462,17 +4465,23 @@ SUBROUTINE W3MPIO ( IMOD )
 #ifdef W3_MPI
           END IF
           !
-#ifdef W3_CESMCOUPLED
           IF ( FLGRDALL( 6, 14) ) THEN
             IH     = IH + 1
             IT     = IT + 1
-            CALL MPI_RECV_INIT (LANGMT(I0),1,WW3_FIELD_VEC, IFROM, IT,  &
-                 MPI_COMM_WAVE, IRQGO2(IH), IERR)
+            CALL MPI_RECV_INIT (USSHX (I0),1,WW3_FIELD_VEC, IFROM, IT,  &
+                 MPI_COMM_WAVE, IRQGO2(IH), IERR )
+#ifdef W3_MPIT
+            WRITE (NDST,9011) IH, ' 6/14', IFROM, IT, IRQGO2(IH), IERR
+#endif
+            IH     = IH + 1
+            IT     = IT + 1
+            CALL MPI_RECV_INIT (USSHY (I0),1,WW3_FIELD_VEC, IFROM, IT,  &
+                 MPI_COMM_WAVE, IRQGO2(IH), IERR )
 #ifdef W3_MPIT
             WRITE (NDST,9011) IH, ' 6/14', IFROM, IT, IRQGO2(IH), IERR
 #endif
           END IF
-#endif ! W3_CESMCOUPLED
+          !
           IF ( FLGRDALL( 7, 1) ) THEN
             IH     = IH + 1
             IT     = IT + 1
diff --git a/model/src/w3iogomd.F90 b/model/src/w3iogomd.F90
index 43dd1bf64..61495a6fe 100644
--- a/model/src/w3iogomd.F90
+++ b/model/src/w3iogomd.F90
@@ -1084,6 +1084,9 @@ SUBROUTINE W3FLDTOIJ(FLD, I, J, IAPROC, NAPOUT, NDSEN)
     CASE('TOC')
       I = 6
       J = 13
+    CASE('USSH')
+      I = 6
+      J = 14
       !
       ! Group 7
       !
@@ -1296,23 +1299,11 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
          ICPRT, DTPRT, WSCUT, NOSWLL, FLOGRD, FLOGR2,&
          NOGRP, NGRPP
     USE W3ADATMD, ONLY: NSEALM
-#ifdef W3_CESMCOUPLED
-    ! USSX, USSY   : surface Stokes drift (SD)
-    ! USSXH, USSYH : surface layer (SL) averaged SD
-    ! LANGMT       : La_t
-    ! LAPROJ       : La_{Proj}
-    ! LASL         : La_{SL}
-    ! LASLPJ       : La_{SL,Proj}
-    ! ALPHAL       : angle between wind and Langmuir cells (SL averaged)
-    ! ALPHALS      : angle between wind and Langmuir cells (surface)
-    ! UD           : wind direction
-    ! LAMULT       : enhancement factor
-    ! HML          : mixing layer depth (from coupler)
-    USE W3ADATMD, ONLY: LAMULT, USSXH, USSYH, LANGMT, LAPROJ, &
-         ALPHAL, ALPHALS, LASL, UD, LASLPJ
-    USE W3IDATMD, ONLY: HML
-    USE W3WDATMD, ONLY: ASF
-#endif
+      ! USSHX, USSHY : surface layer (SL) averaged SD
+      ! HSL          : surface layer depth (1/5 of the mixed layer depth
+      !                from the coupler)
+    USE W3ADATMD, ONLY: USSHX, USSHY
+    USE W3IDATMD, ONLY: HSL
 #ifdef W3_S
     USE W3SERVMD, ONLY: STRACE
 #endif
@@ -1370,17 +1361,12 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
     REAL                       USSCO, FT1
     REAL, SAVE              :: HSMIN = 0.01
     LOGICAL                 :: FLOLOC(NOGRP,NGRPP)
-#ifdef W3_CESMCOUPLED
     ! SWW: angle between wind and waves
-    ! HSL: surface layer depth (=0.2*HML)
-    REAL                    :: SWW !angle between wind and waves
-    REAL                    :: HSL !surface layer depth (=0.2*HML)
-    ! tmp variables for surface and SL averaged SD
-    REAL                    :: ETUSSX(NSEAL),        &
-         ETUSSY(NSEAL),        &
-         ETUSSXH(NSEAL),       &
-         ETUSSYH(NSEAL)
-#endif
+    ! LHSL: local surface layer depth
+    REAL                    :: SWW
+    REAL                    :: LHSL
+    ! tmp variable for surface layer averaged Stokes drift
+    REAL                    :: USSCOH
     !/
     !/ ------------------------------------------------------------------- /
     !/
@@ -1487,25 +1473,10 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
     QP    = UNDEF
     WBT    = UNDEF
     !
-#ifdef W3_CESMCOUPLED
-    ETUSSX  = 0.
-    ETUSSY  = 0.
     ETUSCX  = 0.
     ETUSCY  = 0.
-    ETUSSXH  = 0.
-    ETUSSYH  = 0
-    LANGMT = UNDEF
-    LAPROJ = UNDEF
-    LASL   = UNDEF
-    LASLPJ = UNDEF
-    ALPHAL = UNDEF
-    ALPHALS = UNDEF
-    USSX   = 0.
-    USSY   = 0.
-    USSXH  = 0.
-    USSYH  = 0.
-    LAMULT  = 1.
-#endif
+    USSHX  = 0.
+    USSHY  = 0.
     !
     ! 2.  Integral over discrete part of spectrum ------------------------ *
     !
@@ -1631,12 +1602,16 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
           TPMS(JSEA) = TPI/SIG(IK)
         END IF
 
-#ifdef W3_CESMCOUPLED
-        ! Get surface layer depth
-        IX    = MAPSF(ISEA,1)
-        IY    = MAPSF(ISEA,2)
-        HSL   = HML(IX,IY)/5.     ! depth over which SD is averaged
-#endif
+        IF (LMPENABLED) then
+            IF (HSLMODE.EQ.0) then
+              LHSL = 10.0 ! a constant value for testing purposes
+            ELSE
+              ! Get surface layer depth from coupler
+              IX    = MAPSF(ISEA,1)
+              IY    = MAPSF(ISEA,2)
+              LHSL  = HSL(IX,IY)      ! depth over which SD is averaged
+            END IF
+        END IF
 
         !
         ! Directional moments in the last freq. band
@@ -1677,39 +1652,14 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
           USSCO=FKD*SIG(IK)*WN(IK,ISEA)*COSH(2.*KD)
           BHD(JSEA) = BHD(JSEA) +                             &
                GRAV*WN(IK,ISEA) * EBD(IK,JSEA) / (SINH(2.*KD))
-#ifdef W3_CESMCOUPLED
-          ! Surface Stokes Drift
-          ETUSSX(JSEA)  = ETUSSX(JSEA) + ABX(JSEA)*FACTOR*SIG(IK) &
-               *WN(IK,ISEA)*COSH(2*WN(IK,ISEA)*DW(ISEA))          &
-               /(SINH(WN(IK,ISEA)*DW(ISEA)))**2
-          ETUSSY(JSEA)  = ETUSSY(JSEA) + ABY(JSEA)*FACTOR*SIG(IK) &
-               *WN(IK,ISEA)*COSH(2*WN(IK,ISEA)*DW(ISEA))          &
-               /(SINH(WN(IK,ISEA)*DW(ISEA)))**2
-          ! Depth averaged Stokes Drift
-          ETUSSXH(JSEA)  = ETUSSXH(JSEA) + ABX(JSEA)*FACTOR*SIG(IK) &
-               *(1.-EXP(-2.*WN(IK,ISEA)*HSL))/2./HSL                &
-               *COSH(2*WN(IK,ISEA)*DW(ISEA))                        &
-               /(SINH(WN(IK,ISEA)*DW(ISEA)))**2
-          ETUSSYH(JSEA)  = ETUSSYH(JSEA) + ABY(JSEA)*FACTOR*SIG(IK) &
-               *(1.-EXP(-2.*WN(IK,ISEA)*HSL))/2./HSL                &
-               *COSH(2*WN(IK,ISEA)*DW(ISEA))                        &
-               /(SINH(WN(IK,ISEA)*DW(ISEA)))**2
-#endif
+          IF (LMPENABLED) THEN
+            USSCOH=0.5*FKD*SIG(IK)*(1.-EXP(-2.*WN(IK,ISEA)*LHSL))/LHSL*COSH(2.*KD)
+          ENDIF
         ELSE
           USSCO=FACTOR*SIG(IK)*2.*WN(IK,ISEA)
-#ifdef W3_CESMCOUPLED
-          ! deep water limit
-          ! Surface Stokes Drift
-          ETUSSX(JSEA)  = ETUSSX(JSEA) + ABX(JSEA)*FACTOR*SIG(IK) &
-               *2.*WN(IK,ISEA)
-          ETUSSY(JSEA)  = ETUSSY(JSEA) + ABY(JSEA)*FACTOR*SIG(IK) &
-               *2.*WN(IK,ISEA)
-          ! Depth averaged Stokes Drift
-          ETUSSXH(JSEA)  = ETUSSXH(JSEA) + ABX(JSEA)*FACTOR*SIG(IK) &
-               *(1.-EXP(-2.*WN(IK,ISEA)*HSL))/HSL
-          ETUSSYH(JSEA)  = ETUSSYH(JSEA) + ABY(JSEA)*FACTOR*SIG(IK) &
-               *(1.-EXP(-2.*WN(IK,ISEA)*HSL))/HSL
-#endif
+          IF (LMPENABLED) THEN
+            USSCOH=FACTOR*SIG(IK)*(1.-EXP(-2.*WN(IK,ISEA)*LHSL))/LHSL
+          ENDIF
         END IF
         !
         ABXX(JSEA)   = MAX ( 0. , ABXX(JSEA) ) * FACTOR
@@ -1725,6 +1675,10 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
         !
         USSX(JSEA)  = USSX(JSEA) + ABX(JSEA)*USSCO
         USSY(JSEA)  = USSY(JSEA) + ABY(JSEA)*USSCO
+        IF (LMPENABLED) THEN
+          USSHX(JSEA) = USSHX(JSEA) + ABX(JSEA)*USSCOH
+          USSHY(JSEA) = USSHY(JSEA) + ABY(JSEA)*USSCOH
+        ENDIF
         !
         ! Fills the 3D Stokes drift spectrum array
         !  ! The US3D Stokes drift specrum array is now calculated in a
@@ -1998,11 +1952,17 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
     !
     DO JSEA=1, NSEAL
       CALL INIT_GET_ISEA(ISEA, JSEA)
-#ifdef W3_CESMCOUPLED
-      IX = MAPSF(ISEA,1)
-      IY = MAPSF(ISEA,2)
-      HS = HML(IX,IY)/5.     ! depth over which SD is averaged
-#endif
+
+      IF (LMPENABLED) then
+        IF (HSLMODE.EQ.0) then
+          LHSL = 10.0 ! a constant value for testing purposes
+        ELSE
+          ! Get surface layer depth from coupler
+          IX    = MAPSF(ISEA,1)
+          IY    = MAPSF(ISEA,2)
+          LHSL  = HSL(IX,IY)      ! depth over which SD is averaged
+        END IF
+      END IF
       !
       ! 3.a Directional mss parameters
       !     NB: the slope PDF is proportional to ell1=ETYY*EC2-2*ETXY*ECS+ETXX*ES2 = C*EC2-2*B*ECS+A*ES2
@@ -2032,16 +1992,23 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
       SXX(JSEA) = SXX(JSEA) + FTE * ABXX(JSEA) / CG(NK,ISEA)
       SYY(JSEA) = SYY(JSEA) + FTE * ABYY(JSEA) / CG(NK,ISEA)
       SXY(JSEA) = SXY(JSEA) + FTE * ABXY(JSEA) / CG(NK,ISEA)
-#ifdef W3_CESMCOUPLED
-      ! tail for SD
-      ETUSSX(JSEA)  = ETUSSX(JSEA) + 2*GRAV*ETUSCX(JSEA)/SIG(NK)
-      ETUSSY(JSEA)  = ETUSSY(JSEA) + 2*GRAV*ETUSCY(JSEA)/SIG(NK)
-#endif
       !
       ! Tail for surface stokes drift is commented out: very sensitive to tail power
       !
       !       USSX(JSEA)  = USSX(JSEA) + 2*GRAV*ETUSCX(JSEA)/SIG(NK)
       !       USSY(JSEA)  = USSY(JSEA) + 2*GRAV*ETUSCY(JSEA)/SIG(NK)
+
+      ! Add tail contribution for surface and layer averaged Stokes drift
+      IF (LMPENABLED.and.SDTAIL) then
+        USSX(JSEA)  = USSX(JSEA) + 2*GRAV*ETUSCX(JSEA)/SIG(NK)
+        USSY(JSEA)  = USSY(JSEA) + 2*GRAV*ETUSCY(JSEA)/SIG(NK)
+        USSHX(JSEA) = USSHX(JSEA) + 2*GRAV*ETUSCX(JSEA)/SIG(NK)     &
+          *(1.-(1.-4.*LHSL*WN(NK,ISEA))*EXP(-2.*WN(NK,ISEA)*LHSL))    &
+          /6./WN(NK,ISEA)/LHSL
+        USSHY(JSEA)  = USSHY(JSEA) + 2*GRAV*ETUSCY(JSEA)/SIG(NK)    &
+          *(1.-(1.-4.*LHSL*WN(NK,ISEA))*EXP(-2.*WN(NK,ISEA)*LHSL))    &
+          /6./WN(NK,ISEA)/LHSL
+      END IF
       UBS(JSEA) = UBS(JSEA) + FTWL * EBAND/GRAV
     END DO
     !
@@ -2110,87 +2077,6 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
           T02(JSEA) = TPI / SIG(NK)
           T01(JSEA)= T02(JSEA)
         ENDIF
-#ifdef W3_CESMCOUPLED
-        !TODO is this affected by the NXXX vs. NSEALM?
-        ! Should LAMULT, etc. be NSEAML length?
-        ! Output Stokes drift and Langmuir numbers
-        ! USERO(JSEA,1) = HS(JSEA) / MAX ( 0.001 , DW(JSEA) )
-        ! USERO(JSEA,2) = ASF(ISEA)
-        IF (ETUSSX(JSEA) .NE. 0. .OR. ETUSSY(JSEA) .NE. 0.) THEN
-
-          USSX(JSEA) = ETUSSX(JSEA)
-          USSY(JSEA) = ETUSSY(JSEA)
-          USSXH(JSEA) = ETUSSXH(JSEA)
-          USSYH(JSEA) = ETUSSYH(JSEA)
-
-          ! this check is to divide by zeror error with gx17
-          ! is there a better way to do this check?
-          IF( SQRT(USSX(JSEA)**2 + USSY(JSEA)**2) .GT. 0) THEN
-            IF( SQRT(USSXH(JSEA)**2+USSYH(JSEA)**2) .GT. 0) THEN
-
-              LANGMT(JSEA) = SQRT ( UST(ISEA) * ASF(ISEA)        &
-                   * SQRT ( DAIR / DWAT )                   &
-                   / SQRT ( USSX(JSEA)**2 + USSY(JSEA)**2 ) )
-              ! Calculating Langmuir Number for misaligned wind and waves
-              ! see Van Roekel et al., 2012
-              ! take z1 = 4 * HS
-              ! SWW: angle between Stokes drift and wind
-
-              ! no Stokes depth
-              SWW = ATAN2(USSY(JSEA),USSX(JSEA)) - UD(ISEA)
-              ! ALPHALS: angle between wind and LC direction, Surface
-              ! Stokes drift
-              ! LR check for divide by zero
-              if ((LANGMT(JSEA)**2  &
-                   /0.4*LOG(MAX(ABS(HML(IX,IY)/4./HS(JSEA)),1.0))+COS(SWW)).eq.0.) then
-                print *, 'LR warning A denom 0.'
-                ! This appears to be a decimal precision error
-                ! The first term equals minus the second term to 6 decimal places
-                ! The denominator should be a very small number (e-7)
-                ! ATAN(sin(sww)/small number) tends to pi/2
-                ! So I hardcoded this here.
-                ALPHALS(JSEA) = -1.5707956594501575
-              else
-
-                ALPHALS(JSEA) = ATAN(SIN(SWW) / (LANGMT(JSEA)**2  &
-                     /0.4*LOG(MAX(ABS(HML(IX,IY)/4./HS(JSEA)),1.0))+COS(SWW)))
-              end if
-
-
-              ALPHALS(JSEA) = ATAN( SIN(SWW) / ( LANGMT(JSEA)**2  &
-                   /0.4*LOG(MAX(ABS(HML(IX,IY)/4./HS(JSEA)),1.0))+COS(SWW)))
-              LAPROJ(JSEA) = LANGMT(JSEA) &
-                   * SQRT(ABS(COS(ALPHALS(JSEA))) &
-                   / ABS(COS(SWW-ALPHALS(JSEA))))
-              ! Stokes depth
-              SWW = ATAN2(USSYH(JSEA),USSXH(JSEA)) - UD(ISEA)
-              ! ALPHAL: angle between wind and LC direction
-
-              ! LR check for divide by zero (same as above)
-              if ((LANGMT(JSEA)**2  &
-                   /0.4*LOG(MAX(ABS(HML(IX,IY)/4./HS(JSEA)),1.0))+COS(SWW)).eq.0.) then
-                print *, 'LR warning B denom 0.'
-                ALPHAL(JSEA) = -1.5707956594501575
-              else
-
-                ALPHAL(JSEA) = ATAN(SIN(SWW) / (LANGMT(JSEA)**2  &
-                     /0.4*LOG(MAX(ABS(HML(IX,IY)/4./HS(JSEA)),1.0))+COS(SWW)))
-              end if
-              LASL(JSEA) = SQRT(UST(ISEA)*ASF(ISEA)         &
-                   * SQRT(DAIR/DWAT)                       &
-                   / SQRT(USSXH(JSEA)**2+USSYH(JSEA)**2))
-              LASLPJ(JSEA) = LASL(JSEA) * SQRT(ABS(COS(ALPHAL(JSEA))) &
-                   / ABS(COS(SWW-ALPHAL(JSEA))))
-              ! LAMULT
-              LAMULT(JSEA) = MIN(5.0, ABS(COS(ALPHAL(JSEA))) * &
-                   SQRT(1.0+(1.5*LASLPJ(JSEA))**(-2)+(5.4*real(LASLPJ(JSEA),kind=8))**(-4)))
-              ! user defined output
-              USERO(JSEA,1) = HML(IX,IY)
-              !USERO(JSEA,2) = COS(ALPHAL(JSEA)
-            END IF
-          END IF
-        END IF
-#endif
         !
         !  Add here USERO(JSEA,1) ...
         !
@@ -2669,6 +2555,7 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD )
          TH1M, STH1M, TH2M, STH2M, HSIG, PHICE, TAUICE,&
          STMAXE, STMAXD, HMAXE, HCMAXE, HMAXD, HCMAXD,&
          USSP, TAUOCX, TAUOCY
+    USE W3ADATMD, ONLY: USSHX, USSHY
     !/
     USE W3ODATMD, ONLY: NOGRP, NGRPP, IDOUT, UNDEF, NDST, NDSE,     &
          FLOGRD, IPASS => IPASS1, WRITE => WRITE1,   &
@@ -3043,6 +2930,10 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD )
             TAUOCX(ISEA) = UNDEF
             TAUOCY(ISEA) = UNDEF
           END IF
+          IF ( FLOGRD( 6, 14) ) THEN
+            USSHX (ISEA) = UNDEF
+            USSHY (ISEA) = UNDEF
+          END IF
           !
           IF ( FLOGRD( 7, 1) ) THEN
             ABA   (ISEA) = UNDEF
@@ -3366,6 +3257,9 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD )
             ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 13 ) THEN
               WRITE ( NDSOG ) TAUOCX(1:NSEA)
               WRITE ( NDSOG ) TAUOCY(1:NSEA)
+            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 14 ) THEN
+              WRITE ( NDSOG ) USSHX(1:NSEA)
+              WRITE ( NDSOG ) USSHY(1:NSEA)
               !
               !     Section 7)
               !
@@ -3707,6 +3601,11 @@ SUBROUTINE W3IOGO ( INXOUT, NDSOG, IOTST, IMOD )
                    TAUOCX(1:NSEA)
               READ (NDSOG,END=801,ERR=802,IOSTAT=IERR)         &
                    TAUOCY(1:NSEA)
+            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 14 ) THEN
+              READ (NDSOG,END=801,ERR=802,IOSTAT=IERR)         &
+                   USSHX(1:NSEA)
+              READ (NDSOG,END=801,ERR=802,IOSTAT=IERR)         &
+                   USSHY(1:NSEA)
 
               !
               !     Section 7)
diff --git a/model/src/w3iogoncdmd.F90 b/model/src/w3iogoncdmd.F90
index 7520a0aa2..813aa28d2 100644
--- a/model/src/w3iogoncdmd.F90
+++ b/model/src/w3iogoncdmd.F90
@@ -60,9 +60,7 @@ subroutine w3iogoncd ()
     use w3adatmd   , only : cflxymax, cflthmax, cflkmax, p2sms, us3d
     use w3adatmd   , only : th1m, sth1m, th2m, sth2m, hsig, phice, tauice
     use w3adatmd   , only : stmaxe, stmaxd, hmaxe, hcmaxe, hmaxd, hcmaxd, ussp, tauocx, tauocy
-#ifdef W3_CESMCOUPLED
-    use w3adatmd   , only : langmt
-#endif
+    use w3adatmd , only : usshx, usshy
     use wav_grdout , only : varatts, outvars
     use w3timemd   , only : set_user_timestring
     use w3odatmd   , only : time_origin, calendar_name, elapsed_secs
@@ -355,9 +353,8 @@ subroutine w3iogoncd ()
         if (vname .eq.   'PHICE') call write_var2d(vname, phice    (1:nsea) )
         if (vname .eq.  'TAUOCX') call write_var2d(vname, tauocx   (1:nsea) )
         if (vname .eq.  'TAUOCY') call write_var2d(vname, tauocy   (1:nsea) )
-#ifdef W3_CESMCOUPLED
-        if (vname .eq.  'LANGMT') call write_var2d(vname, langmt   (1:nsea) )
-#endif
+        if (vname .eq.   'USSHX') call write_var2d(vname, usshx    (1:nsea) )
+        if (vname .eq.   'USSHY') call write_var2d(vname, usshy    (1:nsea) )
         ! Group 7
         if (vname .eq.    'ABAX') call write_var2d(vname, aba      (1:nsea), dir=cos(abd(1:nsea)) )
         if (vname .eq.    'ABAY') call write_var2d(vname, aba      (1:nsea), dir=sin(abd(1:nsea)) )
diff --git a/model/src/w3iogrmd.F90 b/model/src/w3iogrmd.F90
index 3aa2688ab..417428591 100644
--- a/model/src/w3iogrmd.F90
+++ b/model/src/w3iogrmd.F90
@@ -993,6 +993,16 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT )
          FTWL, FACTI1, FACTI2, FACHFA, FACHFE
 #endif
     !
+    ! Langmuir mixing parameterization --------------
+    IF ( WRITE ) THEN
+      WRITE (NDSM)                                                &
+        LMPENABLED, SDTAIL, HSLMODE
+      ELSE
+        READ (NDSM,END=801,ERR=802,IOSTAT=IERR)                     &
+             LMPENABLED, SDTAIL, HSLMODE
+    END IF
+
+    ! --------------
     !
     ! Output flags for 3D parameters ------------------------------------- *
     !                                                 Module W3GDATMD
diff --git a/model/src/w3iorsmd.F90 b/model/src/w3iorsmd.F90
index 40d65f413..24d9a280c 100644
--- a/model/src/w3iorsmd.F90
+++ b/model/src/w3iorsmd.F90
@@ -328,6 +328,7 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
     USE W3PARALL, ONLY: PRINT_MY_TIME
 #endif
     USE w3odatmd, ONLY : RUNTYPE, INITFILE
+    USE w3adatmd, ONLY : USSHX, USSHY
 #ifdef W3_PDLIB
     USE PDLIB_FIELD_VEC
 #endif
@@ -1054,6 +1055,10 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
               WRITE(NDSR,ERR=803,IOSTAT=IERR) TAUOCX(1:NSEA)
               WRITE(NDSR,ERR=803,IOSTAT=IERR) TAUOCY(1:NSEA)
             ENDIF
+            IF ( FLOGRR(6,14) ) THEN
+              WRITE(NDSR,ERR=803,IOSTAT=IERR) USSHX(1:NSEA)
+              WRITE(NDSR,ERR=803,IOSTAT=IERR) USSHY(1:NSEA)
+            ENDIF
             IF ( FLOGRR(7,2) ) THEN
               WRITE(NDSR,ERR=803,IOSTAT=IERR) UBA(1:NSEA)
               WRITE(NDSR,ERR=803,IOSTAT=IERR) UBD(1:NSEA)
@@ -1321,6 +1326,17 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
               ENDIF
             ENDDO
           ENDIF
+          IF ( FLOGOA(6,14) ) THEN
+            READ (NDSR,ERR=802,IOSTAT=IERR) TMP(1:NSEA)
+            READ (NDSR,ERR=802,IOSTAT=IERR) TMP2(1:NSEA)
+            DO I=1, NSEALM
+              J = IAPROC + (I-1)*NAPROC
+              IF (J .LE. NSEA) THEN
+                USSHX(I) = TMP(J)
+                USSHY(I) = TMP2(J)
+              ENDIF
+            ENDDO
+          ENDIF
           IF ( FLOGOA(7,2) ) THEN
             READ (NDSR,ERR=802,IOSTAT=IERR) TMP(1:NSEA)
             READ (NDSR,ERR=802,IOSTAT=IERR) TMP2(1:NSEA)
@@ -1404,6 +1420,8 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
           UBD     = 0.
           PHIBBL  = 0.
           TAUBBL  = 0.
+          USSHX   = 0.
+          USSHY   = 0.
         ENDIF
 #ifdef W3_T
         WRITE (NDST,9008)
diff --git a/model/src/w3odatmd.F90 b/model/src/w3odatmd.F90
index 28e97c799..5c1c7d239 100644
--- a/model/src/w3odatmd.F90
+++ b/model/src/w3odatmd.F90
@@ -876,7 +876,7 @@ SUBROUTINE W3NOUT ( NDSERR, NDSTST )
     !
     ! 6) Wave-ocean layer
     !
-    NOGE(6) = 13
+    NOGE(6) = 14
     !
     IDOUT( 6, 1)  = 'Radiation stresses  '
     IDOUT( 6, 2)  = 'Wave-ocean mom. flux'
@@ -891,9 +891,7 @@ SUBROUTINE W3NOUT ( NDSERR, NDSTST )
     IDOUT( 6,11)  = 'Wave-ice energy flux'
     IDOUT( 6,12)  = 'Split Surface Stokes'
     IDOUT( 6,13)  = 'Tot wav-ocn mom flux'
-#ifdef W3_CESMCOUPLED
-    IDOUT( 6,14)  = 'Turbulent Langmuir number'
-#endif
+    IDOUT( 6,14)  = 'Stokes drift sfc ave'
     !
     ! 7) Wave-bottom layer
     !
diff --git a/model/src/w3ounfmetamd.F90 b/model/src/w3ounfmetamd.F90
index 44704c375..03e8e1b62 100644
--- a/model/src/w3ounfmetamd.F90
+++ b/model/src/w3ounfmetamd.F90
@@ -3744,6 +3744,30 @@ SUBROUTINE DEFAULT_META()
     META(2)%VARNS=''
     META(2)%VARNG=''
     META(2)%VARNC='toc=sqrt(utoc**2+vtoc**2)'
+
+    ! IFI=6, IFJ=14, USSH
+    META => GROUP(6)%FIELD(14)%META
+    META(1)%FSC    = 0.0005
+    META(1)%UNITS  = 'm s-1'
+    META(1)%ENAME  = '.ussh'
+
+    ! First component
+    META(1)%VARNM='uussh'
+    META(1)%VARNL='eastward surface averaged stokes drift'
+    META(1)%VARNS='sea_surface_wave_averaged_stokes_drift_eastward_velocity'
+    META(1)%VARNG='eastward_surface_averaged_stokes_drift'
+    META(1)%VARND=DIRCOM
+    META(1)%VMIN = -4.95
+    META(1)%VMAX =  4.95
+
+    ! Second component
+    META(2) = META(1)
+    META(2)%VARNM='vussh'
+    META(2)%VARNL='northward surface averaged stokes drift'
+    META(2)%VARNS='sea_surface_wave_averaged_stokes_drift_northward_velocity'
+    META(2)%VARNG='northward_surface_averaged_stokes_drift'
+
+
     !
     !---------- GROUP 7 ----------------
     !
@@ -3956,6 +3980,7 @@ SUBROUTINE DEFAULT_META()
     META(1)%VARNC='Goda wave peakedness parameter'
     META(1)%VMIN = 0
     META(1)%VMAX = 32
+
     !
     !---------- GROUP 9 ----------------
     !
diff --git a/model/src/wav_comp_nuopc.F90 b/model/src/wav_comp_nuopc.F90
index 6f3eeef5a..c7bb14ef2 100644
--- a/model/src/wav_comp_nuopc.F90
+++ b/model/src/wav_comp_nuopc.F90
@@ -369,7 +369,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
          isSet=isSet, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
     if (isPresent .and. isSet) then
-      wav_coupling_to_cice=(trim(cvalue)=="true")
+       read(cvalue,*) wav_coupling_to_cice
     end if
     write(logmsg,'(A,l)') trim(subname)//': Wave wav_coupling_to_cice setting is ',wav_coupling_to_cice
     call ESMF_LogWrite(trim(logmsg), ESMF_LOGMSG_INFO)
@@ -567,6 +567,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
     if ( root_task ) then
       write(stdout,'(a)')'      *** WAVEWATCH III Program shell ***      '
       write(stdout,'(a)')'==============================================='
+      write(stdout,'(a,l)')' Wave wav_coupling_to_cice setting is ',wav_coupling_to_cice
     end if
 
     !--------------------------------------------------------------------
@@ -917,6 +918,7 @@ subroutine DataInitialize(gcomp, rc)
     type(ESMF_State)  :: exportState
     real(r8), pointer :: z0rlen(:)
     real(r8), pointer :: sw_lamult(:)
+    real(r8), pointer :: sw_lasl(:)
     real(r8), pointer :: sw_ustokes(:)
     real(r8), pointer :: sw_vstokes(:)
     real(r8), pointer :: wave_elevation_spectrum(:,:)
@@ -938,6 +940,14 @@ subroutine DataInitialize(gcomp, rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
       sw_lamult (:) = 1.
     endif
+    if (state_fldchk(exportState, 'Sw_lasl')) then
+      call state_getfldptr(exportState, 'Sw_lasl', sw_lasl, rc=rc)
+      if (ChkErr(rc,__LINE__,u_FILE_u)) return
+      ! note: the default value of this surface layer averaged Langmuir number
+      ! should be a large number to be consistent with lamult=1., ustokes=0.,
+      ! and vstokes=0.
+      sw_lasl (:) = 1.e6
+    endif
     if (state_fldchk(exportState, 'Sw_ustokes')) then
       call state_getfldptr(exportState, 'Sw_ustokes', sw_ustokes, rc=rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -954,7 +964,7 @@ subroutine DataInitialize(gcomp, rc)
       call CalcRoughl(z0rlen)
     endif
     if (wav_coupling_to_cice) then
-      call state_getfldptr(exportState, 'wave_elevation_spectrum', wave_elevation_spectrum, rc=rc)
+      call state_getfldptr(exportState, 'Sw_elevation_spectrum', wave_elevation_spectrum, rc=rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
       wave_elevation_spectrum(:,:) = 0.
     endif
diff --git a/model/src/wav_grdout.F90 b/model/src/wav_grdout.F90
index 7e592e618..4583070d7 100644
--- a/model/src/wav_grdout.F90
+++ b/model/src/wav_grdout.F90
@@ -4,7 +4,7 @@ module wav_grdout
 
   implicit none
 
-  integer, parameter   :: maxvars = 24         ! maximum number of variables/group
+  integer, parameter   :: maxvars = 25         ! maximum number of variables/group
 
   private ! except
 
@@ -222,7 +222,7 @@ subroutine initialize_gridout
          ]
 
     !  6   Wave-ocean layer
-    gridoutdefs(6,1:24) = [ &
+    gridoutdefs(6,1:25) = [ &
          varatts( "SXY  ", "SXX       ", "Radiation stresses xx                           ", "N m-1     ", "  ", .false.) , &
          varatts( "SXY  ", "SYY       ", "Radiation stresses yy                           ", "N m-1     ", "  ", .false.) , &
          varatts( "SXY  ", "SXY       ", "Radiation stresses xy                           ", "N m-1     ", "  ", .false.) , &
@@ -246,7 +246,8 @@ subroutine initialize_gridout
          varatts( "USP  ", "USSPY     ", "Partitioned surface Stokes drift y              ", "m s-1     ", "p ", .false.) , &
          varatts( "TWC  ", "TAUOCX    ", "Total wave to ocean stress x                    ", "Pa        ", "  ", .false.) , &
          varatts( "TWC  ", "TAUOCY    ", "Total wave to ocean stress y                    ", "Pa        ", "  ", .false.) , &
-         varatts( "LAN  ", "LANGMT    ", "Turbulent Langmuir number (La_t)                ", "nd        ", "  ", .false.)   &
+         varatts( "USSH ", "USSHX     ", "Surface layer averaged Stokes drift x           ", "m s-1     ", "  ", .false.) , &
+         varatts( "USSH ", "USSHY     ", "Surface layer averaged Stokes drift y           ", "m s-1     ", "  ", .false.)   &
          ]
 
     !  7   Wave-bottom layer
diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90
index 025edb630..e6caa1a2d 100644
--- a/model/src/wav_import_export.F90
+++ b/model/src/wav_import_export.F90
@@ -20,7 +20,7 @@ module wav_import_export
   use wav_shr_mod  , only : chkerr
   use wav_shr_mod  , only : state_diagnose, state_reset, state_getfldptr, state_fldchk
   use wav_shr_mod  , only : wav_coupling_to_cice, nwav_elev_spectrum, merge_import, dbug_flag, multigrid, unstr_mesh
-  use constants    , only : grav, tpi, dwat
+  use constants    , only : grav, tpi, dwat, dair
   use w3parall     , only : init_get_isea
 
   implicit none
@@ -140,9 +140,9 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc)
     end if
     if (cesmcoupled) then
       call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_lamult' )
+      call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_lasl' )
       call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_ustokes')
       call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_vstokes')
-      !call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_hstokes')
     else
       call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_z0')
     end if
@@ -154,7 +154,7 @@ subroutine advertise_fields(importState, ExportState, flds_scalar_name, rc)
     ! will be implemented soon based on receiving USSP and USSPF from the coupler instead of the mod_def file. This will
     ! also ensure compatibility with the ocean component since ocean will also receive these from the coupler.
     if (wav_coupling_to_cice) then
-      call fldlist_add(fldsFrWav_num, fldsFrWav, 'wave_elevation_spectrum', &
+      call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_elevation_spectrum', &
            ungridded_lbound=1, ungridded_ubound=nwav_elev_spectrum)
     end if
 
@@ -267,7 +267,7 @@ subroutine import_fields( gcomp, time0, timen, rc )
     use w3odatmd    , only: w3seto
     use w3wdatmd    , only: time, w3setw
 #ifdef W3_CESMCOUPLED
-    use w3idatmd    , only: HML
+    use w3idatmd    , only: HSL
 #else
     use wmupdtmd    , only: wmupd2
     use wmmdatmd    , only: wmsetm
@@ -480,8 +480,8 @@ subroutine import_fields( gcomp, time0, timen, rc )
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
     ! ocn mixing layer depth
-    global_data = max(global_data, 5.)
-    call FillGlobalInput(global_data, HML)
+    global_data = max(global_data, 5.)*0.2
+    call FillGlobalInput(global_data, HSL)
 #endif
     ! ---------------
     ! INFLAGS1(5) - atm momentum fields
@@ -588,7 +588,9 @@ subroutine export_fields (gcomp, rc)
     use w3gdatmd      , only : nseal, mapsf, MAPSTA, USSPF, NK, w3setg
     use w3iogomd      , only : CALC_U3STOKES
 #ifdef W3_CESMCOUPLED
-    use w3adatmd      , only : LAMULT
+    use w3wdatmd      , only : ASF, UST
+    use w3adatmd      , only : USSHX, USSHY, UD, HS
+    use w3idatmd      , only : HSL
 #else
     use wmmdatmd      , only : mdse, mdst, wmsetm
 #endif
@@ -600,6 +602,7 @@ subroutine export_fields (gcomp, rc)
     ! Local variables
 #ifdef W3_CESMCOUPLED
     real(R8)          :: fillvalue = 1.0e30_R8                 ! special missing value
+    real              :: sww, langmt, lasl, laslpj, alphal
 #else
     real(R8)          :: fillvalue = zero                      ! special missing value
 #endif
@@ -616,6 +619,7 @@ subroutine export_fields (gcomp, rc)
     real(r8), pointer :: syyn(:)
 
     real(r8), pointer :: sw_lamult(:)
+    real(r8), pointer :: sw_lasl(:)
     real(r8), pointer :: sw_ustokes(:)
     real(r8), pointer :: sw_vstokes(:)
 
@@ -653,13 +657,45 @@ subroutine export_fields (gcomp, rc)
         call init_get_isea(isea, jsea)
         ix  = mapsf(isea,1)
         iy  = mapsf(isea,2)
-        if (mapsta(iy,ix) == 1) then
-          sw_lamult(jsea) = LAMULT(jsea)
+        if (mapsta(iy,ix) == 1 .and. HS(jsea) > zero .and. &
+            sqrt(USSX(jsea)**2+USSY(jsea)**2)>zero .and. sqrt(USSHX(jsea)**2+USSHY(jsea)**2)>zero ) then
+           sww = atan2(USSHY(jsea),USSHX(jsea)) - UD(isea)
+           alphal = atan( sin(sww) / (                                       &
+                          2.5 * UST(isea)*ASF(isea)*sqrt(dair/dwat)          &
+                        / max(1.e-14_r8, sqrt(USSX(jsea)**2+USSY(jsea)**2))     &
+                        * log(max(1.0, abs(1.25*HSL(ix,iy)/HS(jsea))))       &
+                        + cos(sww)   )                                       &
+                        )
+           lasl = sqrt(ust(isea) * asf(isea) * sqrt(dair/dwat) &
+                                 / sqrt(usshx(jsea)**2 + usshy(jsea)**2 ))
+           laslpj = lasl * sqrt(abs(cos(alphal)) &
+               / abs(cos(sww-alphal)))
+           sw_lamult(jsea) = min(5.0, abs(cos(alphal)) * &
+                              sqrt(1.0+(1.5*laslpj)**(-2)+(5.4_r8*laslpj)**(-4)))
         else
           sw_lamult(jsea)  = 1.
         endif
       enddo
     end if
+    if (state_fldchk(exportState, 'Sw_lasl')) then
+      call state_getfldptr(exportState, 'Sw_lasl', sw_lasl, rc=rc)
+      if (ChkErr(rc,__LINE__,u_FILE_u)) return
+      sw_lasl(:) = fillvalue
+      do jsea=1, nseal
+         isea = iaproc + (jsea-1)*naproc
+         ix  = mapsf(isea,1)
+         iy  = mapsf(isea,2)
+         if (mapsta(iy,ix) == 1) then
+            ! note: an arbitrary minimum value of 0.2 is set to avoid zero
+            !       Langmuir number which may result from zero surface friction
+            !       velocity but may cause unphysically strong Langmuir mixing
+            sw_lasl(jsea) = max(0.2, sqrt(UST(isea)*ASF(isea)*sqrt(dair/dwat) &
+                          / max(1.e-14, sqrt(USSHX(jsea)**2+USSHY(jsea)**2))))
+         else
+            sw_lasl(jsea)  = 1.e6
+         endif
+      enddo
+    end if
 #endif
     ! surface stokes drift
     if (state_fldchk(exportState, 'Sw_ustokes')) then
@@ -729,7 +765,7 @@ subroutine export_fields (gcomp, rc)
       call CalcRadstr2D( va, sxxn, sxyn, syyn)
     end if
     if (wav_coupling_to_cice) then
-      call state_getfldptr(exportState, 'wave_elevation_spectrum', wave_elevation_spectrum, rc=rc)
+      call state_getfldptr(exportState, 'Sw_elevation_spectrum', wave_elevation_spectrum, rc=rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
       ! Initialize wave elevation spectrum
       wave_elevation_spectrum(:,:) = fillvalue
diff --git a/model/src/ww3_ounf.F90 b/model/src/ww3_ounf.F90
index 0a2cadfca..3c8ebb4ec 100644
--- a/model/src/ww3_ounf.F90
+++ b/model/src/ww3_ounf.F90
@@ -186,7 +186,7 @@ PROGRAM W3OUNF
        TAUOX, TAUOY, TAUWIX,                        &
        TAUWIY, PHIAW, PHIOC, TUSX, TUSY, PRMS, TPMS,&
        USSX, USSY, MSSX, MSSY, MSSD, MSCX, MSCY,    &
-       MSCD, CHARN, TWS, TAUA, TAUADIR,             &
+       MSCD, CHARN, TWS, TAUA, TAUADIR,USSHX,USSHY, &
        TAUWNX, TAUWNY, BHD, T02, HSIG, CGE,         &
        T01, BEDFORMS, WHITECAP, TAUBBL, PHIBBL,     &
        CFLTHMAX, CFLXYMAX, CFLKMAX, TAUICE, PHICE,  &
@@ -1854,6 +1854,23 @@ SUBROUTINE W3EXNC ( NX, NY, IX1, IXN, IY1, IYN, NSEA,             &
               CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUOCY(1:NSEA)     &
                    , MAPSF, XY )
             ENDIF ! SMCGRD
+            !
+            NFIELD=2
+
+            !
+            ! Surface averaged stokes drift
+          ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 14 ) THEN
+            DO ISEA=1, NSEA
+              USSHX(ISEA)=MAX(-0.9998,MIN(0.9998,USSHX(ISEA)))
+              USSHY(ISEA)=MAX(-0.9998,MIN(0.9998,USSHY(ISEA)))
+            END DO
+#ifdef W3_RTD
+            ! Rotate x,y vector back to standard pole
+            IF ( FLAGUNR ) CALL W3XYRTN(NSEA, USSHX(1:NSEA), USSHY(1:NSEA), AnglD)
+#endif
+            CALL S2GRID(USSHX(1:NSEA), XX)
+            CALL S2GRID(USSHY(1:NSEA), XY)
+            !
             NFIELD=2
             !
             ! RMS of bottom displacement amplitude
diff --git a/model/src/ww3_outf.F90 b/model/src/ww3_outf.F90
index e4c1affed..12f0a0f06 100644
--- a/model/src/ww3_outf.F90
+++ b/model/src/ww3_outf.F90
@@ -167,7 +167,7 @@ PROGRAM W3OUTF
        CFLTHMAX, CFLKMAX, BEDFORMS, WHITECAP, T02, &
        CGE, T01, HSIG, STMAXE, STMAXD, HMAXE,      &
        HCMAXE, HMAXD, HCMAXD, MSSD, MSCD, WBT,     &
-       WNMEAN, TAUA, TAUADIR
+       WNMEAN, TAUA, TAUADIR, USSHX, USSHY
   USE W3ODATMD, ONLY: NDSO, NDSE, NDST, NOGRP, NGRPP, IDOUT,      &
        UNDEF, FLOGRD, FNMPRE, NOSWLL, NOGE
   !
@@ -1940,6 +1940,57 @@ SUBROUTINE W3EXGO ( NX, NY, NSEA )
               CALL W3S2XY (NSEA,NSEA,NX+1,NY, TAUOCY(1:NSEA)        &
                    , MAPSF, X2 )
             ENDIF
+
+            !
+          ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 14 ) THEN
+            IF ( VECTOR ) THEN
+              FLTWO  = .TRUE.
+            ELSE
+              FLDIR  = .TRUE.
+            END IF
+            FSC    = 0.001
+            UNITS  = 'm s-1'
+            ENAME  = '.ussh'
+            DO ISEA=1, NSEA
+              IF (USSHX(ISEA) .NE. UNDEF) THEN
+                USSHX(ISEA)=MAX(-0.9998,MIN(0.9998,USSHX(ISEA)))
+                USSHY(ISEA)=MAX(-0.9998,MIN(0.9998,USSHY(ISEA)))
+              END IF
+            END DO
+#ifdef W3_RTD
+            ! Rotate x,y vector back to standard pole
+            IF ( FLAGUNR ) CALL W3XYRTN(NSEA, USSHX, USSHY, AnglD)
+#endif
+            IF ( ITYPE .EQ. 4 ) THEN
+              XS1    = USSHX(1:NSEA)
+              XS2    = USSHY(1:NSEA)
+            ELSE
+              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, USSHX(1:NSEA)     &
+                   , MAPSF, XX )
+              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, USSHY(1:NSEA)     &
+                   , MAPSF, XY )
+            ENDIF
+            DO ISEA=1, NSEA
+              CABS   = SQRT(USSHX(ISEA)**2+USSHY(ISEA)**2)
+              IF ( USSHX(ISEA)  .NE. UNDEF ) THEN
+                USSHY(ISEA) = MOD ( 630. -                      &
+                     RADE*ATAN2(USSHY(ISEA),USSHX(ISEA)) , 360. )
+              ELSE
+                USSHY(ISEA) = UNDEF
+                CABS        = UNDEF
+              END IF
+              USSHX(ISEA) = CABS
+            END DO
+            IF ( ITYPE .EQ. 4 ) THEN
+              XS3    = USSHX(1:NSEA)
+              XS4    = USSHY(1:NSEA)
+            ELSE
+              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, USSHX(1:NSEA),     &
+                   MAPSF, X1 )
+              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, USSHY(1:NSEA),     &
+                   MAPSF, X2 )
+            ENDIF
+
             !
           ELSE IF ( IFI .EQ. 7 .AND. IFJ .EQ. 1 ) THEN
             IF ( VECTOR ) THEN
diff --git a/regtests/bin/matrix.base b/regtests/bin/matrix.base
index a273372a9..28d78634a 100755
--- a/regtests/bin/matrix.base
+++ b/regtests/bin/matrix.base
@@ -31,7 +31,7 @@
              prop1D prop2D time fetch hur1mg \
              multi01 multi02 multi03 multi04 multi05 \
              hybd shwtr unstr pdlib smcgr rtd mudice infgrv \
-             uost assim calendar confignc oasis \
+             lmp uost assim calendar confignc oasis \
              esmf ufs ufscoarse rstrt_b4b nth_b4b npl_b4b grib \
              multi06 multi07 multi08 multi09
   do
@@ -74,6 +74,7 @@
   echo "  echo '      Growth curves (fetch)              : $fetch'"            >> matrix.head
   echo "  echo '      Shallow water tests                : $shwtr'"            >> matrix.head
   echo "  echo '      Unstructured grid tests            : $unstr'"            >> matrix.head
+  echo "  echo '      Langmuir mixing parameterization   : $lmp'"              >> matrix.head
   echo "  echo '      SMC grid tests                     : $smcgr'"            >> matrix.head
   echo "  echo '      Rotated pole tests                 : $rtd'"              >> matrix.head
   echo "  echo '      Hurricane with one moving grid     : $hur1mg'"           >> matrix.head
@@ -171,7 +172,7 @@
       echo "$rtst -s PR1 -w work_PR1_f -m grdset_f $ww3 mww3_test_04" >> matrix.body
       echo "$rtst -s PR1 -w work_PR1_g -m grdset_g $ww3 mww3_test_04" >> matrix.body
     fi
-  
+
   fi
 
   if [ "$dist" = 'y' ]
@@ -224,7 +225,7 @@
       echo "$rtst -s PR1_MPI -w work_PR1_MPI_d  -m grdset_d  -f -p $mpi -n $np $ww3 mww3_test_04" >> matrix.body
       echo "$rtst -s PR1_MPI -w work_PR1_MPI_e  -m grdset_e  -f -p $mpi -n $np $ww3 mww3_test_04" >> matrix.body
     fi
-  
+
   fi
 
   if [ "$prop2D" = 'y' ] && [ "$dist" = 'y' ]
@@ -812,6 +813,13 @@
       echo "$rtst -s MPI                     -w work_MPI_5km         -g 5km      -f -p $mpi -n $np $ww3 ww3_tp2.15" >> matrix.body
       echo "$rtst -s MPI        -i input_rho -w work_PR3_UQ_RHO_MPI              -f -p $mpi -n $np $ww3 ww3_tp2.15" >> matrix.body
       echo "$rtst -s TIDE_MPI                -w work_TIDE_MPI                    -f -p $mpi -n $np $ww3 ww3_tp2.18" >> matrix.body
+
+      if [ "$lmp" == "y" ];      # Langmuir mixing parameterization
+      then
+        echo ' ' >> matrix.body
+        echo "$rtst                          -w work_LMPN                        -f -p $mpi -n $np $ww3 ww3_tp2.22" >> matrix.body
+        echo "$rtst             -g sdtail    -w work_LMPN_SDTAIL                 -f -p $mpi -n $np $ww3 ww3_tp2.22" >> matrix.body
+      fi
     fi
 
     if [ "$multi01" = 'y' ]
@@ -985,9 +993,9 @@
   if [ "$dist" = 'y' ]
   then
 
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -998,7 +1006,7 @@
       echo "$rtst -s ST1_PR1_MPI     -w work_ST1_PR1_MPI     -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1009,7 +1017,7 @@
       echo "$rtst -s ST1_PR2_UQ_MPI  -w work_ST1_PR2_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1020,7 +1028,7 @@
       echo "$rtst -s ST1_PR2_UNO_MPI -w work_ST1_PR2_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1031,7 +1039,7 @@
       echo "$rtst -s ST1_PR3_UQ_MPI  -w work_ST1_PR3_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1041,9 +1049,9 @@
     then
       echo "$rtst -s ST1_PR3_UNO_MPI -w work_ST1_PR3_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1054,7 +1062,7 @@
       echo "$rtst -s ST2_PR1_MPI     -w work_ST2_PR1_MPI     -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1065,7 +1073,7 @@
       echo "$rtst -s ST2_PR2_UQ_MPI  -w work_ST2_PR2_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1076,7 +1084,7 @@
       echo "$rtst -s ST2_PR2_UNO_MPI -w work_ST2_PR2_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1087,7 +1095,7 @@
       echo "$rtst -s ST2_PR3_UQ_MPI  -w work_ST2_PR3_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1097,9 +1105,9 @@
     then
       echo "$rtst -s ST2_PR3_UNO_MPI -w work_ST2_PR3_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1110,7 +1118,7 @@
       echo "$rtst -s ST3_PR1_MPI     -w work_ST3_PR1_MPI     -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1121,7 +1129,7 @@
       echo "$rtst -s ST3_PR2_UQ_MPI  -w work_ST3_PR2_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1132,7 +1140,7 @@
       echo "$rtst -s ST3_PR2_UNO_MPI -w work_ST3_PR2_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1143,7 +1151,7 @@
       echo "$rtst -s ST3_PR3_UQ_MPI  -w work_ST3_PR3_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1153,9 +1161,9 @@
     then
       echo "$rtst -s ST3_PR3_UNO_MPI -w work_ST3_PR3_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1168,7 +1176,7 @@
       echo "$rtst -s ST4_PR1_MPI     -w work_ST4_PR1_MPI     -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1179,7 +1187,7 @@
       echo "$rtst -s ST4_PR2_UQ_MPI  -w work_ST4_PR2_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1190,7 +1198,7 @@
       echo "$rtst -s ST4_PR2_UNO_MPI -w work_ST4_PR2_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1201,7 +1209,7 @@
       echo "$rtst -s ST4_PR3_UQ_MPI  -w work_ST4_PR3_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1211,9 +1219,9 @@
     then
       echo "$rtst -s ST4_PR3_UNO_MPI -w work_ST4_PR3_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1224,7 +1232,7 @@
       echo "$rtst -s ST6_PR1_MPI     -w work_ST6_PR1_MPI     -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1235,7 +1243,7 @@
       echo "$rtst -s ST6_PR2_UQ_MPI  -w work_ST6_PR2_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1246,7 +1254,7 @@
       echo "$rtst -s ST6_PR2_UNO_MPI -w work_ST6_PR2_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1257,7 +1265,7 @@
       echo "$rtst -s ST6_PR3_UQ_MPI  -w work_ST6_PR3_UQ_MPI  -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1268,7 +1276,7 @@
       echo "$rtst -s ST6_PR3_UNO_MPI -w work_ST6_PR3_UNO_MPI -m grdset_a -f -p $mpi -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
   fi
 
@@ -1278,9 +1286,9 @@
   if [ "$hybd" = 'y' ]
   then
 
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1291,7 +1299,7 @@
       echo "$rtst -s ST1_PR1_MPI_OMPH     -w work_ST1_PR1_MPI_OMPH     -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1302,7 +1310,7 @@
       echo "$rtst -s ST1_PR2_UQ_MPI_OMPH  -w work_ST1_PR2_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1313,7 +1321,7 @@
       echo "$rtst -s ST1_PR2_UNO_MPI_OMPH -w work_ST1_PR2_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1324,7 +1332,7 @@
       echo "$rtst -s ST1_PR3_UQ_MPI_OMPH  -w work_ST1_PR3_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1334,9 +1342,9 @@
     then
       echo "$rtst -s ST1_PR3_UNO_MPI_OMPH -w work_ST1_PR3_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1347,7 +1355,7 @@
       echo "$rtst -s ST2_PR1_MPI_OMPH     -w work_ST2_PR1_MPI_OMPH     -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1358,7 +1366,7 @@
       echo "$rtst -s ST2_PR2_UQ_MPI_OMPH  -w work_ST2_PR2_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1369,7 +1377,7 @@
       echo "$rtst -s ST2_PR2_UNO_MPI_OMPH -w work_ST2_PR2_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1380,7 +1388,7 @@
       echo "$rtst -s ST2_PR3_UQ_MPI_OMPH  -w work_ST2_PR3_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1390,9 +1398,9 @@
     then
       echo "$rtst -s ST2_PR3_UNO_MPI_OMPH -w work_ST2_PR3_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1403,7 +1411,7 @@
       echo "$rtst -s ST3_PR1_MPI_OMPH     -w work_ST3_PR1_MPI_OMPH     -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1414,7 +1422,7 @@
       echo "$rtst -s ST3_PR2_UQ_MPI_OMPH  -w work_ST3_PR2_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1425,7 +1433,7 @@
       echo "$rtst -s ST3_PR2_UNO_MPI_OMPH -w work_ST3_PR2_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1436,7 +1444,7 @@
       echo "$rtst -s ST3_PR3_UQ_MPI_OMPH  -w work_ST3_PR3_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1446,9 +1454,9 @@
     then
       echo "$rtst -s ST3_PR3_UNO_MPI_OMPH -w work_ST3_PR3_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1459,7 +1467,7 @@
       echo "$rtst -s ST4_PR1_MPI_OMPH     -w work_ST4_PR1_MPI_OMPH     -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1470,7 +1478,7 @@
       echo "$rtst -s ST4_PR2_UQ_MPI_OMPH  -w work_ST4_PR2_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1481,7 +1489,7 @@
       echo "$rtst -s ST4_PR2_UNO_MPI_OMPH -w work_ST4_PR2_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1492,7 +1500,7 @@
       echo "$rtst -s ST4_PR3_UQ_MPI_OMPH  -w work_ST4_PR3_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1502,9 +1510,9 @@
     then
       echo "$rtst -s ST4_PR3_UNO_MPI_OMPH -w work_ST4_PR3_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1515,7 +1523,7 @@
       echo "$rtst -s ST6_PR1_MPI_OMPH     -w work_ST6_PR1_MPI_OMPH     -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1526,7 +1534,7 @@
       echo "$rtst -s ST6_PR2_UQ_MPI_OMPH  -w work_ST6_PR2_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1537,7 +1545,7 @@
       echo "$rtst -s ST6_PR2_UNO_MPI_OMPH -w work_ST6_PR2_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1548,7 +1556,7 @@
       echo "$rtst -s ST6_PR3_UQ_MPI_OMPH  -w work_ST6_PR3_UQ_MPI_OMPH  -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1559,7 +1567,7 @@
       echo "$rtst -s ST6_PR3_UNO_MPI_OMPH -w work_ST6_PR3_UNO_MPI_OMPH -m grdset_a -f -p $mpi -n $nr -t $nth $ww3 mww3_test_05" >> matrix.body
     fi
 
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
   fi
 # Moving grid cases for OpenMP, ww3_ts3 and mww3_test_05 share switches
@@ -1568,9 +1576,9 @@
   if [ "$omp" = 'y' ]
   then
 
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1581,7 +1589,7 @@
       echo "$rtst -s ST1_PR1_OMP     -w work_ST1_PR1_OMP     -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1592,7 +1600,7 @@
       echo "$rtst -s ST1_PR2_UQ_OMP  -w work_ST1_PR2_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1603,7 +1611,7 @@
       echo "$rtst -s ST1_PR2_UNO_OMP -w work_ST1_PR2_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1614,7 +1622,7 @@
       echo "$rtst -s ST1_PR3_UQ_OMP  -w work_ST1_PR3_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1624,9 +1632,9 @@
     then
       echo "$rtst -s ST1_PR3_UNO_OMP -w work_ST1_PR3_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1637,7 +1645,7 @@
       echo "$rtst -s ST2_PR1_OMP     -w work_ST2_PR1_OMP     -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1648,7 +1656,7 @@
       echo "$rtst -s ST2_PR2_UQ_OMP  -w work_ST2_PR2_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1659,7 +1667,7 @@
       echo "$rtst -s ST2_PR2_UNO_OMP -w work_ST2_PR2_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1670,7 +1678,7 @@
       echo "$rtst -s ST2_PR3_UQ_OMP  -w work_ST2_PR3_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1680,9 +1688,9 @@
     then
       echo "$rtst -s ST2_PR3_UNO_OMP -w work_ST2_PR3_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1693,7 +1701,7 @@
       echo "$rtst -s ST3_PR1_OMP     -w work_ST3_PR1_OMP     -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1704,7 +1712,7 @@
       echo "$rtst -s ST3_PR2_UQ_OMP  -w work_ST3_PR2_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1715,7 +1723,7 @@
       echo "$rtst -s ST3_PR2_UNO_OMP -w work_ST3_PR2_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1726,7 +1734,7 @@
       echo "$rtst -s ST3_PR3_UQ_OMP  -w work_ST3_PR3_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1736,9 +1744,9 @@
     then
       echo "$rtst -s ST3_PR3_UNO_OMP -w work_ST3_PR3_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1749,7 +1757,7 @@
       echo "$rtst -s ST4_PR1_OMP     -w work_ST4_PR1_OMP     -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1760,7 +1768,7 @@
       echo "$rtst -s ST4_PR2_UQ_OMP  -w work_ST4_PR2_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1771,7 +1779,7 @@
       echo "$rtst -s ST4_PR2_UNO_OMP -w work_ST4_PR2_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1782,7 +1790,7 @@
       echo "$rtst -s ST4_PR3_UQ_OMP  -w work_ST4_PR3_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1792,9 +1800,9 @@
     then
       echo "$rtst -s ST4_PR3_UNO_OMP -w work_ST4_PR3_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1805,7 +1813,7 @@
       echo "$rtst -s ST6_PR1_OMP     -w work_ST6_PR1_OMP     -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1816,7 +1824,7 @@
       echo "$rtst -s ST6_PR2_UQ_OMP  -w work_ST6_PR2_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1827,7 +1835,7 @@
       echo "$rtst -s ST6_PR2_UNO_OMP -w work_ST6_PR2_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1838,7 +1846,7 @@
       echo "$rtst -s ST6_PR3_UQ_OMP  -w work_ST6_PR3_UQ_OMP  -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
 
-    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then 
+    if [ "$hur1mg" = 'y' ] || [ "$multi05" = 'y' ] ; then
       echo ' ' >> matrix.body ; fi
     if [ "$hur1mg" = 'y' ]
     then
@@ -1848,13 +1856,13 @@
     then
       echo "$rtst -s ST6_PR3_UNO_OMP -w work_ST6_PR3_UNO_OMP -m grdset_a -f -O -n $np $ww3 mww3_test_05" >> matrix.body
     fi
-# - - - - - - - - - - - - 
+# - - - - - - - - - - - -
 
   fi
 
 # Second harmonic test case
 
-  if [ "$infgrv" = 'y' ] && [ "$shrd" = 'y' ] 
+  if [ "$infgrv" = 'y' ] && [ "$shrd" = 'y' ]
   then
     echo ' ' >> matrix.body
     echo "$rtst $ww3 -w work_IG1 ww3_tig1.1" >> matrix.body
@@ -1974,7 +1982,7 @@
     echo "$rtst -s IC1IS2        -w work_IC1IS2_1000      -g 1000m                  $ww3 ww3_tic1.4" >> matrix.body
     echo "$rtst -s IC2IS2        -w work_IC2IS2_IC2b      -g IC2b_1000m             $ww3 ww3_tic1.4" >> matrix.body
     echo "$rtst -s IC2IS2        -w work_IC2IS2_IC2d      -g IC2d_1000m             $ww3 ww3_tic1.4" >> matrix.body
-    echo "$rtst -s IC2IS2        -w work_IC2IS2scat       -g scat                   $ww3 ww3_tic2.3" >> matrix.body    
+    echo "$rtst -s IC2IS2        -w work_IC2IS2scat       -g scat                   $ww3 ww3_tic2.3" >> matrix.body
     echo "$rtst -s IC2IS2        -w work_IC2IS2creep      -g creepOnly              $ww3 ww3_tic2.3" >> matrix.body
     echo "$rtst -s IC2IS2        -w work_IC2IS2dissip     -g dissipOnly             $ww3 ww3_tic2.3" >> matrix.body
   fi
@@ -2006,7 +2014,7 @@
    echo "$rtst -s MPI -s NO_PDLIB -w work_ma    -m grdset_a -f -p $mpi -n $np $ww3 ww3_tp2.17"     >> matrix.body
    echo "$rtst -s MPI -s PDLIB    -w work_b     -g b        -f -p $mpi -n $np $ww3 ww3_tp2.17"     >> matrix.body
    echo "$rtst -s MPI -s PDLIB    -w work_c     -g c        -f -p $mpi -n $np $ww3 ww3_tp2.17"     >> matrix.body
-   echo "$rtst -s MPI -s PDLIB    -w work_pdlib -g pdlib    -f -p $mpi -n $np $ww3 ww3_tp2.6"      >> matrix.body 
+   echo "$rtst -s MPI -s PDLIB    -w work_pdlib -g pdlib    -f -p $mpi -n $np $ww3 ww3_tp2.6"      >> matrix.body
    echo "$rtst -s MPI -s PDLIB    -w work_mb    -m grdset_b -f -p $mpi -n $np $ww3 ww3_tp2.17"     >> matrix.body
    echo "$rtst -s MPI -s PDLIB    -w work_mc    -m grdset_c -f -p $mpi -n $np $ww3 ww3_tp2.17"     >> matrix.body
      if [ "$rstrt_b4b" = 'y' ]
@@ -2023,7 +2031,7 @@
  fi
 
  # unstructured grid with PDLIB for domain decomposition and implicit solver, Neumann Bnd, depth breaking & triad interaction ST
- if [ "$pdlib" = 'y' ] && [ "$dist" = 'y' ] 
+ if [ "$pdlib" = 'y' ] && [ "$dist" = 'y' ]
  then
    echo ' ' >> matrix.body
    echo "$rtst -s MPI -s PDLIB -w work_1A_a -f -g a -p $mpi -n $np -o netcdf -i input_Case1A $ww3 ww3_tp2.19" >> matrix.body
@@ -2032,23 +2040,23 @@
  fi
 
 
- #Unresolved Obstacles Source Term (UOST) 
- if [ "$uost" = 'y' ] 
+ #Unresolved Obstacles Source Term (UOST)
+ if [ "$uost" = 'y' ]
  then
    echo ' ' >> matrix.body
    if [ "$dist" = 'y' ]
-   then 
+   then
      echo "$rtst -s MPI -w work_rg_shel_MPI  -i input_rg_shel            -f -p $mpi -n $np $ww3 ww3_ts4" >> matrix.body
      echo "$rtst -s MPI -w work_rg_multi_MPI -i input_rg_multi -m grdset -f -p $mpi -n $np $ww3 ww3_ts4" >> matrix.body
      echo "$rtst -s MPI -w work_ug_MPI       -i input_ug                 -f -p $mpi -n $np $ww3 ww3_ts4" >> matrix.body
-   else 
+   else
      echo "$rtst        -w work_rg_shel      -i input_rg_shel                              $ww3 ww3_ts4" >> matrix.body
      echo "$rtst        -w work_rg_multi     -i input_rg_multi -m grdset                   $ww3 ww3_ts4" >> matrix.body
      echo "$rtst        -w work_ug           -i input_ug                                   $ww3 ww3_ts4" >> matrix.body
    fi
  fi
 
- # Global unstructured grid and Unresolved Obstacles Source Term (UOST) 
+ # Global unstructured grid and Unresolved Obstacles Source Term (UOST)
  if [ "$uost" = 'y' ] && [ "$dist" = 'y' ]
  then
    echo ' ' >> matrix.body
@@ -2062,7 +2070,10 @@
    fi
  fi
 
- #Test of updating the restart spectra 
+
+
+
+ #Test of updating the restart spectra
  if [ "$assim" = 'y' ]
  then
    echo ' ' >> matrix.body
@@ -2082,7 +2093,7 @@
    echo "$rtst -s ST4 -w work_UPD6_U_cap -i input_UPD6_U_cap $ww3 ww3_ta1" >> matrix.body
   fi
 
- #Test of atmosphere, ocean, and ice coupling using OASIS 
+ #Test of atmosphere, ocean, and ice coupling using OASIS
  if [ "$oasis" = 'y' ] && [ "$dist" = 'y' ]
  then
    echo ' ' >> matrix.body
@@ -2204,15 +2215,15 @@
     fi
   fi
 
- #Test gint for restarts 
+ #Test gint for restarts
   if [ "$ufs" = 'y' ] && [ "$grib" = 'y' ]
   then
-    echo ' ' >> matrix.body    
+    echo ' ' >> matrix.body
     echo "$rtst -s MPI_OMPH -w work_c -m grdset_c -f -p $mpi -n $npl  -t $nth1 $ww3 ww3_ufs1.2"            >> matrix.body
     echo "mkdir -p ww3_ufs1.2/work_l"                                                                      >> matrix.body
     echo "cp ww3_ufs1.2/work_c/restart.hafsl ww3_ufs1.2/work_l/restart.hafsl"                              >> matrix.body
     echo "$rtst -s MPI_OMPH -w work_l -m grdset_l -f -p $mpi -n $npl  -t $nth1 $ww3 ww3_ufs1.2"            >> matrix.body
-  fi 
+  fi
 
 
  #GEFSv12 setup with ww3_multi and grib2 output
diff --git a/regtests/bin/matrix_cmake_datarmor b/regtests/bin/matrix_cmake_datarmor
index f8ffaaa63..a197d3a03 100755
--- a/regtests/bin/matrix_cmake_datarmor
+++ b/regtests/bin/matrix_cmake_datarmor
@@ -27,7 +27,7 @@ usage ()
 {
 cat 2>&1 << EOF
 
-Usage: $myname model_dir 
+Usage: $myname model_dir
 Required:
   model_dir : path to model dir of WW3 source
 EOF
@@ -51,7 +51,7 @@ main_dir="`cd $main_dir 1>/dev/null 2>&1 && pwd`"
 
 
 # 1. Set up
-   
+
   export np='28'      #number of mpi tasks
   export npl='28'    #number of mpi tasks for ufs applications and large setups
   export npl1='20'   #number of mpi tasks for ufs/large setups (b4b check)
@@ -151,9 +151,10 @@ main_dir="`cd $main_dir 1>/dev/null 2>&1 && pwd`"
   export      shwtr='y' # shallow water tests
   export      unstr='y' # unstructured grid tests
   export      pdlib='y' # unstr with pdlib for domain decomposition and implicit solver
+  export        lmp='y' # Langmuir mixing parameterization
   export      smcgr='y' # SMC grid test
   export        rtd='y' # Rotated pole test
-  export     mudice='y' # Mud/Ice and wave interaction tests 
+  export     mudice='y' # Mud/Ice and wave interaction tests
   export     infgrv='y' # Second harmonic generation tests
   export       uost='y' # ww3_ts4 Unresolved Obstacles Source Term (UOST)
   export      assim='y' # Restart spectra update
@@ -178,17 +179,17 @@ main_dir="`cd $main_dir 1>/dev/null 2>&1 && pwd`"
   export    nth_b4b='y' # Thread Reproducibility
   export       esmf='n' # ESMF coupling
 # export   filter='PR3 ST2 UQ'
-                      # The filter does a set of consecutive greps on the 
+                      # The filter does a set of consecutive greps on the
                       # command lines generated by filter.base with the above
                       # selected options.
 
 # --------------------------------------------------------------------------- #
 # 2.  Execute matrix.base ...                                                 #
 # --------------------------------------------------------------------------- #
-             
+
   $main_dir/../regtests/bin/matrix.base
 
-  $main_dir/../regtests/bin/matrix_divider_cmake.sh 
+  $main_dir/../regtests/bin/matrix_divider_cmake.sh
 
 # --------------------------------------------------------------------------- #
 # End to the matrix                                                           #
diff --git a/regtests/bin/matrix_cmake_milhydro b/regtests/bin/matrix_cmake_milhydro
index 8213b9b33..6be47ab24 100755
--- a/regtests/bin/matrix_cmake_milhydro
+++ b/regtests/bin/matrix_cmake_milhydro
@@ -19,10 +19,10 @@
 #                                                                             #
 # --------------------------------------------------------------------------- #
 usage ()
-{ 
+{
   cat 2>&1 << EOF
 
-  Usage: $myname model_dir 
+  Usage: $myname model_dir
   Required:
     model_dir : path to model dir of WW3 source
 EOF
@@ -41,7 +41,7 @@ EOF
 # Convert main_dir to absolute path
   main_dir="`cd $main_dir 1>/dev/null 2>&1 && pwd`"
 
-# Set batchq queue, choose modules and other custom variables to fit system and 
+# Set batchq queue, choose modules and other custom variables to fit system and
 # to define headers etc (default to original version if empty)
 # 1. Set up
 # 1.a Computer/ user dependent set up
@@ -58,7 +58,7 @@ EOF
 echo "  export path_build_root=$(dirname $main_dir)/regtests/buildmatrix" >> matrix.head
   echo '  [[ -d ${path_build_root} ]] && rm -rf ${path_build_root}'         >> matrix.head
 
-  echo ' ' 
+  echo ' '
 
   if [ "$batchq" = 'slurm' ]
   then
@@ -97,9 +97,10 @@ echo "  export path_build_root=$(dirname $main_dir)/regtests/buildmatrix" >> mat
   export      shwtr='y' # shallow water tests
   export      unstr='y' # unstructured grid tests
   export      pdlib='y' # unstr with pdlib for domain decomposition and implicit solver
+  export        lmp='y' # Langmuir mixing parameterization
   export      smcgr='y' # SMC grid test
   export        rtd='y' # Rotated pole test
-  export     mudice='y' # Mud/Ice and wave interaction tests 
+  export     mudice='y' # Mud/Ice and wave interaction tests
   export     infgrv='y' # Second harmonic generation tests
   export       uost='y' # ww3_ts4 Unresolved Obstacles Source Term (UOST)
   export      assim='y' # Restart spectra update
@@ -124,14 +125,14 @@ echo "  export path_build_root=$(dirname $main_dir)/regtests/buildmatrix" >> mat
   export    npl_b4b='n' # MPI task Reproducibility
   export    nth_b4b='n' # Thread Reproducibility
   export       esmf='n' # ESMF coupling
-                      # The filter does a set of consecutinve greps on the 
+                      # The filter does a set of consecutinve greps on the
                       # command lines generated by filter.base with the above
                       # selected options.
 
 # --------------------------------------------------------------------------- #
 # 2.  Execute matrix.base ...                                                 #
 # --------------------------------------------------------------------------- #
-             
+
 
   $main_dir/../regtests/bin/matrix.base
 
@@ -146,7 +147,7 @@ echo "  export path_build_root=$(dirname $main_dir)/regtests/buildmatrix" >> mat
     do
       echo "sbatch < $file"  >> msuball.sh
     done
-  fi 
+  fi
 
 # --------------------------------------------------------------------------- #
 # End to the matrix                                                           #
diff --git a/regtests/bin/matrix_cmake_ncep b/regtests/bin/matrix_cmake_ncep
index 771b6f96e..8d1ca00e5 100755
--- a/regtests/bin/matrix_cmake_ncep
+++ b/regtests/bin/matrix_cmake_ncep
@@ -203,6 +203,7 @@ EOF
   export      shwtr='y' # shallow water tests
   export      unstr='y' # unstructured grid tests
   export      pdlib='y' # unstr with pdlib for domain decomposition and implicit solver
+  export        lmp='y' # Langmuir mixing parameterization
   export      smcgr='y' # SMC grid test
   export        rtd='y' # Rotated pole test
   export     mudice='y' # Mud/Ice and wave interaction tests
diff --git a/regtests/bin/matrix_cmake_ukmo_cray b/regtests/bin/matrix_cmake_ukmo_cray
index fc6cf64b9..a4df5da0a 100755
--- a/regtests/bin/matrix_cmake_ukmo_cray
+++ b/regtests/bin/matrix_cmake_ukmo_cray
@@ -29,7 +29,7 @@ usage ()
 {
   cat 2>&1 << EOF
 
-  Usage: $myname model_dir 
+  Usage: $myname model_dir
   Required:
     model_dir : path to model dir of WW3 source
 EOF
@@ -139,14 +139,15 @@ fi
   export      shwtr='y' # shallow water tests
   export      unstr='y' # unstructured grid tests
   export      pdlib='y' # unstr with pdlib for domain decomposition and implicit solver
+  export        lmp='y' # Langmuir mixing parameterization
   export      smcgr='y' # SMC grid test
   export        rtd='y' # Rotated pole test
-  export     mudice='y' # Mud/Ice and wave interaction tests 
+  export     mudice='y' # Mud/Ice and wave interaction tests
   export     infgrv='y' # Second harmonic generation tests
   export       uost='y' # ww3_ts4 Unresolved Obstacles Source Term (UOST)
   export      assim='y' # Restart spectra update
   export      oasis='y' # Atmosphere, ocean, and ice coupling using oasis
-  export   calendar='y' # Calendar type 
+  export   calendar='y' # Calendar type
   export   confignc='y' # Configurable netCDF meta data (ww3_ounf)
 
   export    multi01='y' # mww3_test_01 (wetting and drying)
@@ -166,14 +167,14 @@ fi
   export    nth_b4b='y' # Thread Reproducibility
   export       esmf='n' # ESMF coupling
 # export   filter='PR3 ST2 UQ'
-                      # The filter does a set of consecutinve greps on the 
+                      # The filter does a set of consecutinve greps on the
                       # command lines generated by filter.base with the above
                       # selected options.
 
 # --------------------------------------------------------------------------- #
 # 2.  Execute matrix.base ...                                                 #
 # --------------------------------------------------------------------------- #
-             
+
   $main_dir/../regtests/bin/matrix.base
 
   $main_dir/../regtests/bin/matrix_divider_cmake.sh
diff --git a/regtests/ww3_tp2.22/info b/regtests/ww3_tp2.22/info
new file mode 100644
index 000000000..c5e483a3c
--- /dev/null
+++ b/regtests/ww3_tp2.22/info
@@ -0,0 +1,58 @@
+#############################################################################
+#                                                                           #
+# ww3_tp2.22 Test script for WW-III demonstrating Langmuir Mixing           #
+#            Parmeterization for the case of two-dimensional propagation    #
+#            over half the globe. Based on regtest ww3_tp2.2 with           #
+#            switch_PR3_UQ_MPI used as a sample switch.                     #
+#                                                                           #
+#            Two test cases are run. Both have LMPENABLED, with one which   #
+#            includes the tail contribution, the other does not.            #
+#                                                                           #
+# Remarks :                                                                 #
+# - Test case input parameters:                                             #
+#   I. NO TAIL                                                              #
+#   ----------                                                              #
+#     + LMPENABLED = T:  Langmuir mixing parameterization enabled.          #
+#     + SDTAIL     = F:  tail contribution NOT enabled.                     #
+#     + HSLMODE    = 0:  uniform (10m) surface layer depth is assumed.      #
+#                                                                           #
+#   II. TAIL INCLUDED                                                       #
+#   -----------------                                                       #
+#     + LMPENABLED = T:  Langmuir mixing parameterization enabled.          #
+#     + SDTAIL     = T:  tail contribution enabled.                         #
+#     + HSLMODE    = 0:  uniform (10m) surface layer depth is assumed.      #
+#                                                                           #
+# - Test case output parameters:                                            #
+#     + USSX/USSY:   surface Stokes drift           (netCDF:   uuss/vuss)   #
+#     + USSHX/USSHY: surface averaged Stokes drift  (netCDF: uussh/vussh)   #
+#                                                                           #
+# Sample run_test commands :                                                #
+#   (Note: mpirun commands differ by local system)                          #
+# ./bin/run_cmake_test -n 4 -p mpirun -f           ../model ww3_tp2.22      #
+# ./bin/run_cmake_test -n 4 -p mpirun -f -g sdtail ../model ww3_tp2.22      #
+#                                                                           #
+# Notes on LMPN namelist and Langmuir Mixing Parameterization:              #
+#                                                                           #
+#    The LMPN namelist group configures the Langmuir Mixing                 #
+#    Parameterization used within the context of climate modeling and       #
+#    when WW3 is coupled with an active ocean model (.e.g., POP2, MOM6).    #
+#    The LMPENABLED flag turns on the parameterization computations.        #
+#    The SDTAIL flag enables the spectral tail contribution.                #
+#    HSLMODE parameter determines the surface layer depth (HSL).            #
+#    When HSLMODE == 0, a uniform 10m HSL is assumed for testing purposes.  #
+#    When HSLMODE == 1, the HSL field is to be received from the active     #
+#    ocean model via coupler.                                               #
+#                                                                           #
+#    For more info on the Langmuir parameterization, see:                   #
+#       Li, Qing, et al. "Langmuir mixing effects on global climate:        #
+#       WAVEWATCH III in CESM." Ocean Modelling 103 (2016): 145-160.        #
+#                                                                           #
+#                                              Hendrik Tolman, Jun 2002     #
+#                                   Last Mod : Li / Altuntas,  Nov 2023     #
+#                                                                           #
+#    Copyright 2009-2023 National Weather Service (NWS),                    #
+#       National Oceanic and Atmospheric Administration.  All rights        #
+#       reserved.  WAVEWATCH III is a trademark of the NWS.                 #
+#       No unauthorized use without permission.                             #
+#                                                                           #
+#############################################################################
diff --git a/regtests/ww3_tp2.22/input/switch b/regtests/ww3_tp2.22/input/switch
new file mode 100644
index 000000000..0f07d4935
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/switch
@@ -0,0 +1 @@
+NOGRB DIST MPI PR3 UQ FLX2 LN0 ST0 NL0 BT0 DB0 TR0 BS0 IC0 IS0 REF0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7 O10 O11
diff --git a/regtests/ww3_tp2.22/input/track_i.ww3 b/regtests/ww3_tp2.22/input/track_i.ww3
new file mode 100644
index 000000000..ae2154a44
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/track_i.ww3
@@ -0,0 +1,6 @@
+WAVEWATCH III TRACK LOCATIONS DATA 
+19680606 000000       0       0    S1A
+19680606 040000       1       0    S1B
+19680606 060000       1       0    S1C
+19680606 080000       2       0    S1D
+19680606 120000       0.5     0    S1E
diff --git a/regtests/ww3_tp2.22/input/ww3_grid.inp b/regtests/ww3_tp2.22/input/ww3_grid.inp
new file mode 100644
index 000000000..83213e0e1
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_grid.inp
@@ -0,0 +1,33 @@
+$ WAVEWATCH III Grid preprocessor input file
+$ ------------------------------------------
+  'LANGMUIR MIXING PARAMETERIZATION FOR 2-D PROP. AROUND GLOBE'
+$
+   1.1 0.040 9  24  0.
+$
+   F T T T F F
+  2200.  2200. 2200.  2200.
+$
+  &PRO2 DTIME = 345600. /
+  &PRO3 WDTHCG = 1.50, WDTHTH =1.50 /
+  &OUTS E3D = 1 /
+  &LMPN LMPENABLED = T, SDTAIL = F, HSLMODE = 0 /
+END OF NAMELISTS
+$
+  'RECT' T 'NONE'
+  193   93
+    1.   1.  1.
+   -6. -46.  1.
+$
+  -5. 5.75  10  -2500. 2 1 '(....)' 'UNIT' 'input'
+$
+   17949*1
+$
+  10  1 1 '(....)' 'PART' 'input'
+$
+   0  0  F
+   0  0  F
+   0  0
+$
+   0. -40.  0.5  0.866  100
+   0.  40.  0.5 -0.866 -100
+   0. 0. 0. 0.  0
diff --git a/regtests/ww3_tp2.22/input/ww3_grid_sdtail.inp b/regtests/ww3_tp2.22/input/ww3_grid_sdtail.inp
new file mode 100644
index 000000000..e021a094a
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_grid_sdtail.inp
@@ -0,0 +1,33 @@
+$ WAVEWATCH III Grid preprocessor input file
+$ ------------------------------------------
+  'LANGMUIR MIXING PARAMETERIZATION FOR 2-D PROP. AROUND GLOBE'
+$
+   1.1 0.040 9  24  0.
+$
+   F T T T F F
+  2200.  2200. 2200.  2200.
+$
+  &PRO2 DTIME = 345600. /
+  &PRO3 WDTHCG = 1.50, WDTHTH =1.50 /
+  &OUTS E3D = 1 /
+  &LMPN LMPENABLED = T, SDTAIL = T, HSLMODE = 0 /
+END OF NAMELISTS
+$
+  'RECT' T 'NONE'
+  193   93
+    1.   1.  1.
+   -6. -46.  1.
+$
+  -5. 5.75  10  -2500. 2 1 '(....)' 'UNIT' 'input'
+$
+   17949*1
+$
+  10  1 1 '(....)' 'PART' 'input'
+$
+   0  0  F
+   0  0  F
+   0  0
+$
+   0. -40.  0.5  0.866  100
+   0.  40.  0.5 -0.866 -100
+   0. 0. 0. 0.  0
diff --git a/regtests/ww3_tp2.22/input/ww3_ounf.inp b/regtests/ww3_tp2.22/input/ww3_ounf.inp
new file mode 100644
index 000000000..6d1998f36
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_ounf.inp
@@ -0,0 +1,43 @@
+$ -------------------------------------------------------------------- $
+$ WAVEWATCH III Grid output post-processing                            $
+$--------------------------------------------------------------------- $
+$ First output time (yyyymmdd hhmmss), increment of output (s),
+$ and number of output times.
+$
+  19680606 000000 3600. 12
+$
+$ Fields requested --------------------------------------------------- $
+$
+$ Output request flags identifying fields as in ww3_shel.inp. See that
+$ file for a full documentation of field output options. Namelist type
+$ selection is used here (for alternative F/T flags, see ww3_shel.inp).
+$
+ N
+ HS T01 DIR SPR USS USSH
+$
+$--------------------------------------------------------------------- $
+$ netCDF version [3,4]
+$        and variable type 4 [2 = SHORT, 3 = it depends , 4 = REAL]
+$ swell partitions [0 1 2 3 4 5]
+$ variables in same file [T] or not [F]
+$
+ 4 4
+ 0 1 2
+ T
+$
+$ -------------------------------------------------------------------- $
+$ File prefix
+$ number of characters in date [0(nodate),4(yearly),6(monthly),8(daily),10(hourly)]
+$ IX and IY ranges [regular:IX NX IY NY, unstructured:IP NP 1 1]
+$
+ ww3.
+ 6
+ 1 100000 1 100000
+$
+$ For each field and time a new file is generated with the file name
+$ ww3.date_xxx.nc , where date is a conventional time indicator with S3
+$ characters, and xxx is a field identifier.
+$
+$ -------------------------------------------------------------------- $
+$ End of input file                                                    $
+$ -------------------------------------------------------------------- $
diff --git a/regtests/ww3_tp2.22/input/ww3_ounp.inp b/regtests/ww3_tp2.22/input/ww3_ounp.inp
new file mode 100644
index 000000000..4358accbd
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_ounp.inp
@@ -0,0 +1,119 @@
+$ -------------------------------------------------------------------- $
+$ WAVEWATCH III NETCDF Point output post-processing                    $
+$--------------------------------------------------------------------- $
+$ First output time (yyyymmdd hhmmss), increment of output (s),
+$ and number of output times.
+$
+  19680606 000000 21600.  1000
+$
+$ Points requested --------------------------------------------------- $
+$
+$ Define points index for which output is to be generated.
+$ If no one defined, all points are selected
+$ One index number per line, negative number identifies end of list.
+ 1
+ 2
+ 3
+$ mandatory end of list
+ -1
+$
+$--------------------------------------------------------------------- $
+$ file prefix
+$ number of characters in date [4(yearly),6(monthly),8(daily),10(hourly)]
+$ netCDF version [3,4]
+$ points in same file [T] or not [F]
+$                    and max number of points to be processed in one pass
+$ output type ITYPE [0,1,2,3]
+$ flag for global attributes WW3 [0] or variable version [1-2-3-4]
+$ flag for dimensions order time,station [T] or station,time [F]
+$
+ ww3.
+ 6
+ 4
+ T 150
+ 1
+ 0
+ T
+$
+$ -------------------------------------------------------------------- $
+$ ITYPE = 0, inventory of file.
+$            No additional input, the above time range is ignored.
+$
+$ -------------------------------------------------------------------- $
+$ ITYPE = 1, netCDF Spectra.
+$          - Sub-type OTYPE :  1 : Print plots.
+$                              2 : Table of 1-D spectra
+$                              3 : Transfer file.
+$                              4 : Spectral partitioning.
+$          - Scaling factors for 1-D and 2-D spectra Negative factor
+$            disables, output, factor = 0. gives normalized spectrum.
+$          - Netcdf variable type [2=SHORT, 3=it depends, 4=REAL]
+$
+  3  1  0  4
+$
+$ The transfer file contains records with the following contents.
+$
+$ - File ID in quotes, number of frequencies, directions and points.
+$   grid name in quotes (for unformatted file C*21,3I,C*30).
+$ - Bin frequencies in Hz for all bins.
+$ - Bin directions in radians for all bins (Oceanographic conv.).
+$                                                         -+
+$ - Time in yyyymmdd hhmmss format                         | loop
+$                                             -+           |
+$ - Point name (C*40), lat, lon, d, U10 and    |  loop     | over
+$   direction, current speed and direction     |   over    |
+$ - E(f,theta)                                 |  points   | times
+$                                             -+          -+
+$
+$ -------------------------------------------------------------------- $
+$ ITYPE = 2, netCDF Tables of (mean) parameter
+$          - Sub-type OTYPE :  1 : Depth, current, wind
+$                              2 : Mean wave pars.
+$                              3 : Nondimensional pars. (U*)
+$                              4 : Nondimensional pars. (U10)
+$                              5 : 'Validation table'
+$                              6 : WMO standard output
+$  4
+$
+$ -------------------------------------------------------------------- $
+$ ITYPE = 3, netCDF Source terms
+$          - Sub-type OTYPE :  1 : Print plots.
+$                              2 : Table of 1-D S(f).
+$                              3 : Table of 1-D inverse time scales
+$                                  (1/T = S/F).
+$                              4 : Transfer file
+$          - Scaling factors for 1-D and 2-D source terms. Negative
+$            factor disables print plots, factor = 0. gives normalized
+$            print plots.
+$          - Flags for spectrum, input, interactions, dissipation,
+$            bottom and total source term.
+$          - scale ISCALE for OTYPE=2,3
+$                              0 : Dimensional.
+$                              1 : Nondimensional in terms of U10
+$                              2 : Nondimensional in terms of U*
+$                             3-5: like 0-2 with f normalized with fp.
+$
+$  4  0  0  T T T T T T  0
+$
+$ The transfer file contains records with the following contents.
+$
+$ - File ID in quotes, nubmer of frequencies, directions and points,
+$   flags for spectrum and source terms (C*21, 3I, 6L)
+$ - Bin frequencies in Hz for all bins.
+$ - Bin directions in radians for all bins (Oceanographic conv.).
+$                                                         -+
+$ - Time in yyyymmdd hhmmss format                         | loop
+$                                             -+           |
+$ - Point name (C*40), depth, wind speed and   |  loop     | over
+$   direction, current speed and direction     |   over    |
+$ - E(f,theta) if requested                    |  points   | times
+$ - Sin(f,theta) if requested                  |           |
+$ - Snl(f,theta) if requested                  |           |
+$ - Sds(f,theta) if requested                  |           |
+$ - Sbt(f,theta) if requested                  |           |
+$ - Stot(f,theta) if requested                 |           |
+$                                             -+          -+
+$
+$ -------------------------------------------------------------------- $
+$ End of input file                                                    $
+$ -------------------------------------------------------------------- $
diff --git a/regtests/ww3_tp2.22/input/ww3_outf.inp b/regtests/ww3_tp2.22/input/ww3_outf.inp
new file mode 100644
index 000000000..2ee125996
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_outf.inp
@@ -0,0 +1,10 @@
+$ WAVEWATCH III Grid output post-processing
+$ -----------------------------------------
+  19680606 000000 3600. 12
+$
+  N
+  HS T01 DIR SPR USS USSH
+$
+  1 0
+$
+  1 193 1 1 93 1 T F F
diff --git a/regtests/ww3_tp2.22/input/ww3_outp_spec.inp b/regtests/ww3_tp2.22/input/ww3_outp_spec.inp
new file mode 100644
index 000000000..d4c452b1f
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_outp_spec.inp
@@ -0,0 +1,11 @@
+$ WAVEWATCH III Point output post-processing
+$ ------------------------------------------
+  19680606 000000 43200. 25
+$
+  1
+  2
+  3
+ -1
+$
+  1
+  1  -1.  0.  33   F
diff --git a/regtests/ww3_tp2.22/input/ww3_outp_tab51.inp b/regtests/ww3_tp2.22/input/ww3_outp_tab51.inp
new file mode 100644
index 000000000..d40f4a55d
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_outp_tab51.inp
@@ -0,0 +1,10 @@
+$ WAVEWATCH III Point output post-processing
+$ ------------------------------------------
+  19680606 000000 43200. 25
+$
+  1
+ -1
+$
+  2
+$
+  2  51
diff --git a/regtests/ww3_tp2.22/input/ww3_outp_tab52.inp b/regtests/ww3_tp2.22/input/ww3_outp_tab52.inp
new file mode 100644
index 000000000..06431db5d
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_outp_tab52.inp
@@ -0,0 +1,10 @@
+$ WAVEWATCH III Point output post-processing
+$ ------------------------------------------
+  19680606 000000 43200. 25
+$
+  2
+ -1
+$
+  2
+$
+  2  52
diff --git a/regtests/ww3_tp2.22/input/ww3_outp_tab53.inp b/regtests/ww3_tp2.22/input/ww3_outp_tab53.inp
new file mode 100644
index 000000000..209a605b0
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_outp_tab53.inp
@@ -0,0 +1,10 @@
+$ WAVEWATCH III Point output post-processing
+$ ------------------------------------------
+  19680606 000000 43200. 25
+$
+  3
+ -1
+$
+  2
+$
+  2  53
diff --git a/regtests/ww3_tp2.22/input/ww3_shel.inp b/regtests/ww3_tp2.22/input/ww3_shel.inp
new file mode 100644
index 000000000..8a32e5584
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_shel.inp
@@ -0,0 +1,42 @@
+$ WAVEWATCH III shell input file
+$ ------------------------------
+   F T
+   T T
+   T T
+   F F
+   F F
+   F F
+   F
+   F
+   F
+$
+   19680606 000000
+   19680606 120000
+$
+   1
+$
+   19680606 000000  3600  19680608 000000
+$
+  N
+  HS T01 DIR SPR USS USSH
+$
+   19680606 000000  21600  19680608 000000
+     0.0   0.0   'LEFT'
+    90.0   0.0   'CENTER'
+   180.0   0.0   'RIGHT'
+     0.0   0.0   'STOPSTRING'
+   19680606 000000  14400  19680608 000000
+   T
+   19680606 000000      0  19680608 000000
+   19680606 000000      0  19680608 000000
+   19680606 000000      0  19680608 000000
+$
+   'CUR' 19680606 030000   2.0  45.
+   'WND' 19680606 000000  20.0  180.  2.
+   'WND' 19680606 040000  15.0  130.  1.
+   'WND' 19680606 080000  25.0   90.  3.
+   'STP'
+$
+$ -------------------------------------------------------------------- $
+$ End of input file                                                    $
+$ -------------------------------------------------------------------- $
diff --git a/regtests/ww3_tp2.22/input/ww3_strt.inp b/regtests/ww3_tp2.22/input/ww3_strt.inp
new file mode 100644
index 000000000..62e9be21a
--- /dev/null
+++ b/regtests/ww3_tp2.22/input/ww3_strt.inp
@@ -0,0 +1,5 @@
+$ WAVEWATCH III Initial conditions input file
+$ -------------------------------------------
+  1
+$  fp       sip    thm  ncos xm  six   ym  siy hmax
+   0.06     0.01   310. 8    0.0 -1.0  0.0 1.0 2.5