Charnock parameter not set in ST6 NOAA-EMC#1307

Functional type is after Janssen (1991): CHARN = AHAT / SQRT(1-TAU_WAV/TAU) Changes ------- [x] Compute Charnock parameter in subroutine LFACTOR. [x] Add minimum Charock coefficient SIN6AHAT to namelist SIN6 (Eqv. to WAM-type ALPHA0). [x] Regression test `ww3_ts1` succeeded (output changed for ST6 gridded output). ToDo ---- [ ] Try alternatives w/ dynamically calculated AHAT.
stefanzieger · Oct 11, 2024 · 3780f17 · 3780f17
1 parent abc77b9
commit 3780f17
Show file tree

Hide file tree

Showing 9 changed files with 46 additions and 27 deletions.
diff --git a/model/src/gx_outp.F90 b/model/src/gx_outp.F90
@@ -1201,7 +1201,7 @@ SUBROUTINE GXEXPO
 #endif
 #ifdef W3_ST6
           CALL W3SIN6 (A, CG, WN2, UABS, USTAR, UDIRR, CD, &
-               DAIR, TAUWX, TAUWY, TAUWNX, TAUWNY, XWI, DIA )
+               DAIR, TAUWX, TAUWY, TAUWNX, TAUWNY, CHARN, XWI, DIA )
 #endif
         END IF
         IF ( FLSRCE(3) ) THEN

diff --git a/model/src/w3flx4md.F90 b/model/src/w3flx4md.F90
@@ -223,7 +223,7 @@ SUBROUTINE W3FLX4 ( ZWND, U10, U10D, UST, USTD, Z0, CD )
     ! Formats
     !
 1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3FLX4 : '/               &
-         '     HIGHT OF WIND SHOULD BE 10m IN THIS APPRACH '/   &
+         '     HIGHT OF WIND SHOULD BE 10m IN THIS APPROACH '/   &
          '     ZWND =',F8.2,'m'/)
     !/
     !/ End of W3FLX4 ----------------------------------------------------- /

diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90
@@ -915,7 +915,7 @@ MODULE W3GDATMD
     !
 #ifdef W3_ST6
     REAL                  :: SIN6A0, SDS6A1, SDS6A2, SWL6B1, &
-         SIN6WS, SIN6FC
+         SIN6WS, SIN6FC, SIN6AHAT
     INTEGER               :: SDS6P1, SDS6P2
     LOGICAL               :: SDS6ET, SWL6S6, SWL6CSTB1
 #endif
@@ -1333,7 +1333,7 @@ MODULE W3GDATMD
 #endif
 #ifdef W3_ST6
   REAL, POINTER           :: SIN6A0, SDS6A1, SDS6A2, SWL6B1, &
-       SIN6WS, SIN6FC
+       SIN6WS, SIN6FC, SIN6AHAT
   INTEGER, POINTER        :: SDS6P1, SDS6P2
   LOGICAL, POINTER        :: SDS6ET, SWL6S6, SWL6CSTB1
 #endif
@@ -2708,6 +2708,7 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     SIN6A0 => MPARS(IMOD)%SRCPS%SIN6A0
     SIN6WS => MPARS(IMOD)%SRCPS%SIN6WS
     SIN6FC => MPARS(IMOD)%SRCPS%SIN6FC
+    SIN6AHAT => MPARS(IMOD)%SRCPS%SIN6AHAT
     SDS6ET => MPARS(IMOD)%SRCPS%SDS6ET
     SDS6A1 => MPARS(IMOD)%SRCPS%SDS6A1
     SDS6P1 => MPARS(IMOD)%SRCPS%SDS6P1

diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90
@@ -866,7 +866,7 @@ MODULE W3GRIDMD
 #endif
   !
 #ifdef W3_ST6
-  REAL                    :: SINA0, SINWS, SINFC,                 &
+  REAL                    :: SINA0, SINWS, SINFC, SINAHAT,        &
        SDSA1, SDSA2, SWLB1
   INTEGER                 :: SDSP1, SDSP2
   LOGICAL                 :: SDSET, CSTB1
@@ -1053,7 +1053,7 @@ MODULE W3GRIDMD
 #endif
 
 #ifdef W3_ST6
-  NAMELIST /SIN6/ SINA0, SINWS, SINFC
+  NAMELIST /SIN6/ SINA0, SINWS, SINFC, SINAHAT
   NAMELIST /SDS6/ SDSET, SDSA1, SDSA2, SDSP1, SDSP2
   NAMELIST /SWL6/ SWLB1, CSTB1
 #endif
@@ -1744,6 +1744,7 @@ SUBROUTINE W3GRID()
     SINA0  = 0.09
     SINWS  = 32.0
     SINFC  = 6.0
+    SINAHAT  = 0.0095
 #endif
     !
 #ifdef W3_ST1
@@ -1839,9 +1840,10 @@ SUBROUTINE W3GRID()
     SIN6A0 = SINA0
     SIN6WS = SINWS
     SIN6FC = SINFC
+    SIN6AHAT = SINAHAT
     J = 1
     IF ( SIN6A0.LE.0. ) J = 2
-    WRITE (NDSO,921) YESXNO(J), SIN6A0, SIN6WS, SIN6FC
+    WRITE (NDSO,921) YESXNO(J), SIN6A0, SIN6WS, SIN6FC, SIN6AHAT
 #endif
     !
     ! 6.e Define Snl.
@@ -3236,7 +3238,7 @@ SUBROUTINE W3GRID()
            CAPCHA, CHAMIN, CHA0, UCAP, SIGMAUCAP
 #endif
 #ifdef W3_ST6
-      WRITE (NDSO,2920) SINA0, SINWS, SINFC
+      WRITE (NDSO,2920) SINA0, SINWS, SINFC, SINAHAT
 #endif
 #ifdef W3_NL1
       WRITE (NDSO,2922) LAMBDA, NLPROP, KDCONV, KDMIN,       &
@@ -6295,8 +6297,10 @@ SUBROUTINE W3GRID()
 921 FORMAT ( '  negative wind input active       :  ',A/    &
          '  attenuation factor               :  ',F6.2/ &
          '  wind speed scaling factor        :  ',F6.2/ &
-         '  frequency cut-off factor         :  ',F6.2/)
-2920 FORMAT ( '  &SIN6 SINA0 =', F6.3, ', SINWS =', F6.2, ', SINFC =', F6.2, ' /')
+         '  frequency cut-off factor         :  ',F6.2/ &
+         '  minimum Charnock coefficient     :  ',F8.5/ )
+2920 FORMAT ( '  &SIN6 SINA0 =', F6.3, ', SINWS =', F6.2, &
+         ', SINFC =', F6.2, ', SINAHAT =', F8.5, ' /')
 #endif
     !
 #ifdef W3_NL0

diff --git a/model/src/w3iogrmd.F90 b/model/src/w3iogrmd.F90
@@ -1573,20 +1573,20 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
     IF ( WRITE ) THEN
       WRITE (NDSM) SIN6A0, SDS6ET, SDS6A1, SDS6A2,           &
            SDS6P1, SDS6P2, SWL6S6, SWL6B1, SWL6CSTB1,        &
-           SIN6WS, SIN6FC
+           SIN6WS, SIN6FC, SIN6AHAT
 #ifdef W3_ASCII
       WRITE (NDSA,*) 'SIN6A0, SDS6ET, SDS6A1, SDS6A2,        &
            SDS6P1, SDS6P2, SWL6S6, SWL6B1, SWL6CSTB1,        &
-           SIN6WS, SIN6FC:',                                 &
+           SIN6WS, SIN6FC, AIN6AHAT:',                       &
                    SIN6A0, SDS6ET, SDS6A1, SDS6A2,           &
            SDS6P1, SDS6P2, SWL6S6, SWL6B1, SWL6CSTB1,        &
-           SIN6WS, SIN6FC
+           SIN6WS, SIN6FC, SIN6AHAT
 #endif
     ELSE
       READ (NDSM,END=801,ERR=802,IOSTAT=IERR)                &
            SIN6A0, SDS6ET, SDS6A1, SDS6A2,                   &
            SDS6P1, SDS6P2, SWL6S6, SWL6B1, SWL6CSTB1,        &
-           SIN6WS, SIN6FC
+           SIN6WS, SIN6FC, SIN6AHAT
     END IF
 #endif
     !

diff --git a/model/src/w3src6md.F90 b/model/src/w3src6md.F90
@@ -22,7 +22,7 @@
 !>
 !> @author S. Zieger
 !> @author Q. Liu
-!> @date   26-Jun-2018
+!> @date   11-Oct-2024
 !>
 MODULE W3SRC6MD
   !/
@@ -31,14 +31,15 @@ MODULE W3SRC6MD
   !/                  |           S. Zieger               |
   !/                  |           Q. Liu                  |
   !/                  |                        FORTRAN 90 |
-  !/                  | Last update :         26-Jun-2018 |
+  !/                  | Last update :         11-Oct-2024 |
   !/                  +-----------------------------------+
   !/
   !/    29-May-2009 : Origination (w3srcxmd.ftn)          ( version 3.14 )
   !/    10-Feb-2011 : Implementation of source terms      ( version 4.04 )
   !/                                                         (S. Zieger)
   !/    26-Jun-2017 : Recalibration of ST6                ( verison 6.06 )
   !/                                                         (Q. Liu   )
+  !/    11-Oct-2024 : Charnock parameter output added     ( verison 7.14 )
   !/
   !/    Copyright 2009 National Weather Service (NWS),
   !/       National Oceanic and Atmospheric Administration.  All rights
@@ -280,6 +281,7 @@ END SUBROUTINE W3SPR6
   !> @param[out] TAUWY    Component of the wave-supported stress.
   !> @param[out] TAUWNX   Component of the negative part of the stress.
   !> @param[out] TAUWNY   Component of the negative part of the stress.
+  !> @param[out] CHARN    Charnck parameter (sea-state dependent).
   !> @param[out] S        Source term.
   !> @param[out] D        Diagonal term of derivative.
   !>
@@ -288,7 +290,7 @@ END SUBROUTINE W3SPR6
   !> @date   13-Aug-2021
   !>
   SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
-       TAUWX, TAUWY, TAUNWX, TAUNWY, S, D )
+       TAUWX, TAUWY, TAUNWX, TAUNWY, CHARN, S, D )
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III      NOAA/NCEP/NOPP |
@@ -335,6 +337,7 @@ SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
     !      D¹       R.A. O  Diagonal term of derivative
     !      TAUWX-Y  Real O  Component of the wave-supported stress
     !      TAUNWX-Y Real O  Component of the negative part of the stress
+    !      CHARN    Real O  Charnock parameter
     !      ¹ Stored as 1-D array with dimension NTH*NK (column by column).
     !     ----------------------------------------------------------------
     !
@@ -389,7 +392,7 @@ SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
     !/ Parameter list
     REAL, INTENT(IN)       :: A (NSPEC), CG(NK), WN2(NSPEC)
     REAL, INTENT(IN)       :: UABS, USTAR, USDIR, CD, DAIR
-    REAL, INTENT(OUT)      :: TAUWX, TAUWY, TAUNWX, TAUNWY
+    REAL, INTENT(OUT)      :: TAUWX, TAUWY, TAUNWX, TAUNWY, CHARN
     REAL, INTENT(OUT)      :: S(NSPEC), D(NSPEC)
     !/
     !/ ------------------------------------------------------------------- /
@@ -487,7 +490,7 @@ SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
     CINV    = CINV2(IKN)
     SDENSIG = RESHAPE(S*SIG2/CG2,(/ NTH,NK /))
     CALL LFACTOR(SDENSIG, CINV, UABS, USTAR, USDIR, SIG, DSII, &
-         LFACT, TAUWX, TAUWY                          )
+         LFACT, TAUWX, TAUWY, CHARN                   )
     !
     !/ 6) --- apply reduction (LFACT) to the entire spectrum ------------- /
     IF (SUM(LFACT) .LT. NK) THEN
@@ -761,13 +764,14 @@ END SUBROUTINE W3SDS6
   !> @param[out] LFACT  Factor array.
   !> @param[out] TAUWX  Component of the wave-supported stress.
   !> @param[out] TAUWY  Component of the wave-supported stress.
+  !> @param[out] CHARN  Charnock parameter.
   !>
   !> @author S. Zieger
   !> @author Q. Liu
-  !> @date   26-Jun-2018
+  !> @date   11-Oct-2024
   !>
   SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
-       LFACT, TAUWX, TAUWY                    )
+       LFACT, TAUWX, TAUWY, CHARN             )
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III           NOAA/NCEP |
@@ -812,6 +816,10 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
     !        using reduction factor:
     !                          LFACT(F) = MIN(1,exp((1-U/C(F))*RTAU))
     !        Then alter RTAU and repeat 3) until our constraint is matched.
+    !     4) Charnock parameter after equation (3.47) (Komen el al, 1994):
+    !                                         SIN6AHAT
+    !                          CHARN = -----------------------
+    !                                  SQRT( 1.0 - TAU_W/TAU )
     !
     !  3. Parameters :
     !
@@ -826,6 +834,7 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
     !      DSII    R.A. I  Frequency bandwiths [in rad.]
     !      LFACTOR R.A. O  Factor array                       LFACT(sigma)
     !      TAUWX-Y Real O  Component of the wave-supported stress
+    !      CHARN   Real O  Charnock parameter
     !     ----------------------------------------------------------------
     !
     !  4. Subroutines used :
@@ -848,7 +857,7 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
     !/
     USE CONSTANTS, ONLY: DAIR, GRAV, TPI
     USE W3GDATMD,  ONLY: NK, NTH, NSPEC, DTH, XFR, ECOS, ESIN
-    USE W3GDATMD,  ONLY: SIN6WS
+    USE W3GDATMD,  ONLY: SIN6WS, SIN6AHAT
     USE W3ODATMD,  ONLY: NDST, NDSE, IAPROC, NAPERR
     USE W3TIMEMD,  ONLY: STME21
     USE W3WDATMD,  ONLY: TIME
@@ -866,6 +875,7 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
     REAL, INTENT(IN)  :: DSII(NK)       ! frequency bandwidths
     REAL, INTENT(OUT) :: LFACT(NK)      ! correction factor
     REAL, INTENT(OUT) :: TAUWX, TAUWY   ! normal stress components
+    REAL, INTENT(OUT) :: CHARN          ! Charnock parameter
     !
     !/    --- local parameters (in order of appearance) ------------------ /
 #ifdef W3_S
@@ -1020,6 +1030,9 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
     !
     LFACT(1:NK) = LF10Hz(1:NK)
     !
+    !/ 4) --- Sea-state depended Charnoc parameter ----------------------- /
+    CHARN = SIN6AHAT / SQRT(1.0 - MIN(TAU_WAV/TAU, 0.999))
+    !
 #ifdef W3_T6
     WRITE (NDST,273) 'Sin ', IDTIME(1:19), SDENS10Hz*TPI
     WRITE (NDST,273) 'SinR', IDTIME(1:19), SDENS10Hz*LF10Hz*TPI

diff --git a/model/src/w3srcemd.F90 b/model/src/w3srcemd.F90
@@ -14,7 +14,7 @@
 !>
 !> @author H. L. Tolman
 !> @author F. Ardhuin
-!> @date   22-Mar-2021
+!> @date   11-Oct-2024
 !>
 !> @copyright Copyright 2009-2022 National Weather Service (NWS),
 !>       National Oceanic and Atmospheric Administration.  All rights
@@ -28,7 +28,7 @@ MODULE W3SRCEMD
   !/                  |           H. L. Tolman            |
   !/                  |            F. Ardhuin             |
   !/                  |                        FORTRAN 90 |
-  !/                  | Last update :         22-Mar-2021 |
+  !/                  | Last update :         11-Oct-2024 |
   !/                  +-----------------------------------+
   !/
   !/    For updates see subroutine.
@@ -265,6 +265,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD,          &
     !/    22-Mar-2021 : Add extra fields used in coupling   ( version 7.13 )
     !/    07-Jun-2021 : S_{nl5} GKE NL5 (Q. Liu)            ( version 7.13 )
     !/    19-Jul-2021 : Momentum and air density support    ( version 7.14 )
+    !/    11-Oct-2024 : Provide CHARN to ST6 (S. Zieger)
     !/
     !/    Copyright 2009-2013 National Weather Service (NWS),
     !/       National Oceanic and Atmospheric Administration.  All rights
@@ -1214,7 +1215,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD,          &
 
 #ifdef W3_ST6
       CALL W3SIN6 ( SPEC, CG1, WN2, U10ABS, USTAR, USTDIR, CD, DAIR, &
-           TAUWX, TAUWY, TAUWAX, TAUWAY, VSIN, VDIN )
+           TAUWX, TAUWY, TAUWAX, TAUWAY, CHARN, VSIN, VDIN )
 #endif
       !
       ! 2.b Nonlinear interactions.

diff --git a/model/src/ww3_ounp.F90 b/model/src/ww3_ounp.F90
@@ -2437,7 +2437,7 @@ SUBROUTINE W3EXNC(I,NCID,NREQ,INDREQ,ORDER)
 #endif
 #ifdef W3_ST6
               CALL W3SIN6 (A, CG, WN2, UABS, USTAR, UDIRR, CD, DAIR, &
-                   TAUWX, TAUWY, TAUWNX, TAUWNY, XIN, DIA )
+                   TAUWX, TAUWY, TAUWNX, TAUWNY, CHARN, XIN, DIA )
 #endif
             END IF
             IF ( FLSRCE(3) ) THEN

diff --git a/model/src/ww3_outp.F90 b/model/src/ww3_outp.F90
@@ -1993,7 +1993,7 @@ SUBROUTINE W3EXPO
 #endif
 #ifdef W3_ST6
             CALL W3SIN6 (A, CG, WN2, UABS, USTAR, UDIRR, CD, DAIR, &
-                 TAUWX, TAUWY, TAUWNX, TAUWNY, XIN, DIA )
+                 TAUWX, TAUWY, TAUWNX, TAUWNY, CHARN, XIN, DIA )
 #endif
           END IF
           IF ( FLSRCE(3) ) THEN