Clean up and refine Langmuir turbulence parameterization computations. (

NOAA-EMC#1034)

Authors: @qingli411, @alperaltuntas

Co-authored-by: qingli411 <lqingpku@gmail.com>
Co-authored-by: Matthew Masarik <Matthew.Masarik@noaa.gov>

.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

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.

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

model/inp/ww3_shel.inp

@@ -190,14 +190,15 @@ $  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
 $  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

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

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