Functional form of Charnock calculation in ST4. (#1225)

Authored by: @ukmo-nievesvaliente
NOAA-EMC · May 10, 2024 · 14e6bed · 14e6bed
1 parent c13d0b1
commit 14e6bed
Show file tree

Hide file tree

Showing 14 changed files with 312 additions and 15 deletions.
diff --git a/model/inp/ww3_grid.inp b/model/inp/ww3_grid.inp
@@ -73,6 +73,12 @@ $                                    effective wind speed (!/STAB2).
 $   WAM4 and variants  : Namelist SIN3
 $                           ZWND    : Height of wind (m).
 $                           ALPHA0  : minimum value of Charnock coefficient
+$                           CAPCHA  : enable functional form for Cd (1=enabled),
+$                                       only valid when SINTABLE=0.
+$                             UCAP  : Wind speed threshold for cap (CAPCHA=1)
+$                        SIGMAUCAP  : Width of functional form (CAPCHA=1)
+$                             CHA0  : Initial Charnock coefficient (CAPCHA=1)
+$                           CHAMIN  : Minimum Charnock value (CAPCHA=1)
 $                           Z0MAX   : maximum value of air-side roughness z0
 $                           BETAMAX : maximum value of wind-wave coupling
 $                           SINTHP  : power of cosine in wind input

diff --git a/model/nml/namelists.nml b/model/nml/namelists.nml
@@ -62,6 +62,12 @@ $
 $   Janssen / Ardhuin   : Namelist SIN4
 $                             ZWND : Height of wind (m).
 $                           ALPHA0 : minimum value of Charnock coefficient
+$                           CAPCHA : enable functional form for Cd (1=enabled),
+$                                      only valid when SINTABLE=0.
+$                             UCAP : Wind speed threshold for cap (CAPCHA=1)
+$                        SIGMAUCAP : Width of functional form (CAPCHA=1)
+$                             CHA0 : Initial Charnock coefficient (CAPCHA=1)
+$                           CHAMIN : Minimum Charnock value (CAPCHA=1)
 $                            Z0MAX : maximum value of air-side roughness z0
 $                          BETAMAX : maximum value of wind-wave coupling
 $                           SINTHP : power of cosine in wind input

diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90
@@ -907,7 +907,8 @@ MODULE W3GDATMD
          SSDSP, WWNMEANP, SSTXFTF, SSTXFTWN,  &
          FFXPM, FFXFM, FFXFA,   &
          SSDSBRF1, SSDSBRF2, SSDSBINT,SSDSBCK,&
-         SSDSHCK, SSDSABK, SSDSPBK, SSINBR
+         SSDSHCK, SSDSABK, SSDSPBK, SSINBR,   &
+         CAPCHNK(1:10)
     REAL                  :: ZZWND
     REAL                  :: SSDSCOS, SSDSDTH, SSDSBT, SSDSBM(0:4)
 #endif
@@ -1327,7 +1328,8 @@ MODULE W3GDATMD
        SSDSPBK, SSINBR,SSINTHP,TTAUWSHELTER,&
        SINTAILPAR(:), SSWELLF(:), SSDSC(:), SSDSBR,        &
        SSDSP, WWNMEANP, SSTXFTF, SSTXFTWN,  &
-       SSDSBT, SSDSCOS, SSDSDTH, SSDSBM(:)
+       SSDSBT, SSDSCOS, SSDSDTH, SSDSBM(:), &
+       CAPCHNK(:) 
 #endif
 #ifdef W3_ST6
   REAL, POINTER           :: SIN6A0, SDS6A1, SDS6A2, SWL6B1, &
@@ -2668,6 +2670,7 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     ZZALP    => MPARS(IMOD)%SRCPS%ZZALP
     TTAUWSHELTER  => MPARS(IMOD)%SRCPS%TTAUWSHELTER
     SINTAILPAR  => MPARS(IMOD)%SRCPS%SINTAILPAR
+    CAPCHNK  => MPARS(IMOD)%SRCPS%CAPCHNK
     SSWELLFPAR  => MPARS(IMOD)%SRCPS%SSWELLFPAR
     SSWELLF  => MPARS(IMOD)%SRCPS%SSWELLF
     SSDSC    => MPARS(IMOD)%SRCPS%SSDSC

diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90
@@ -115,6 +115,7 @@ MODULE W3GRIDMD
   !/    19-Jul-2021 : Momentum and air density support    ( version 7.14 )
   !/    28-Feb-2023 : GQM as an alternative for NL1       ( version 7.15 )
   !/    11-Jan-2024 : New namelist parameters for IC4     ( version 7.15 )
+  !/    03-May-2024 : New CAPCHNK parameters for SIN4     ( version 7.15 )
   !/
   !/    Copyright 2009-2013 National Weather Service (NWS),
   !/       National Oceanic and Atmospheric Administration.  All rights
@@ -860,7 +861,8 @@ MODULE W3GRIDMD
        SDSBM0, SDSBM1, SDSBM2, SDSBM3,      &
        SDSBM4, SDSFACMTF, SDSCUMP,  SDSNUW, &
        SDSL, SDSMWD, SDSMWPOW, SPMSS, SDSNMTF, SINTAIL1, SINTAIL2, &
-       CUMSIGP, VISCSTRESS
+       CUMSIGP, VISCSTRESS,                 &
+       CAPCHA, CHAMIN, CHA0, UCAP, SIGMAUCAP
 #endif
   !
 #ifdef W3_ST6
@@ -1003,7 +1005,8 @@ MODULE W3GRIDMD
   NAMELIST /SIN4/ ZWND, ALPHA0, Z0MAX, BETAMAX, SINTHP, ZALP, &
        TAUWSHELTER, SWELLFPAR, SWELLF,                 &
        SWELLF2, SWELLF3, SWELLF4, SWELLF5, SWELLF6,    &
-       SWELLF7, Z0RAT, SINBR, SINTABLE, SINTAIL1, SINTAIL2, TAUWBUG, VISCSTRESS
+       SWELLF7, Z0RAT, SINBR, SINTABLE, SINTAIL1, SINTAIL2, TAUWBUG, VISCSTRESS, &
+       CAPCHA, CHAMIN, CHA0, UCAP, SIGMAUCAP
 #endif
 #ifdef W3_NL1
   NAMELIST /SNL1/ LAMBDA, NLPROP, KDCONV, KDMIN,                  &
@@ -1730,6 +1733,11 @@ SUBROUTINE W3GRID()
     TAUWBUG  = 1  !  TAUWBUG is 1 is the bug is kept:
     !  initializes TAUWX/Y to zero in W3SRCE
     VISCSTRESS =0
+    CAPCHA   = 0.     ! =1 indicates capping of drag is active
+    CHAMIN   = 0.0001 ! 
+    CHA0     = ALPHA0 ! initial value for charnock
+    UCAP     = 30.    ! U10 threshold from which drag capping is applied
+    SIGMAUCAP = 10.   ! Width for reduction of drag beyond UCAP
 #endif
     !
 #ifdef W3_ST6
@@ -1818,6 +1826,11 @@ SUBROUTINE W3GRID()
     SINTAILPAR(3) = SINTAIL2
     SINTAILPAR(4) = FLOAT(TAUWBUG)
     SINTAILPAR(5) = VISCSTRESS
+    CAPCHNK(1) = CAPCHA
+    CAPCHNK(2) = CHAMIN
+    CAPCHNK(3) = CHA0
+    CAPCHNK(4) = UCAP
+    CAPCHNK(5) = SIGMAUCAP
 #endif
     !
 #ifdef W3_ST6
@@ -3219,7 +3232,8 @@ SUBROUTINE W3GRID()
 #ifdef W3_ST4
       WRITE (NDSO,2920) ZWND, ALPHA0, Z0MAX, BETAMAX, SINTHP, ZALP,   &
            TAUWSHELTER, SWELLFPAR, SWELLF, SWELLF2, SWELLF3, SWELLF4, &
-           SWELLF5, SWELLF6, SWELLF7, Z0RAT, SINBR, SINTABLE, TAUWBUG, VISCSTRESS, SINTAIL1, SINTAIL2
+           SWELLF5, SWELLF6, SWELLF7, Z0RAT, SINBR, SINTABLE, TAUWBUG, VISCSTRESS, SINTAIL1, SINTAIL2, &
+           CAPCHA, CHAMIN, CHA0, UCAP, SIGMAUCAP
 #endif
 #ifdef W3_ST6
       WRITE (NDSO,2920) SINA0, SINWS, SINFC
@@ -6270,7 +6284,9 @@ SUBROUTINE W3GRID()
          '        SWELLF5 =',F8.5,', SWELLF6 =',F8.5,            &
          ', SWELLF7 =',F12.2,', Z0RAT =',F8.5,', SINBR =',F8.5,','/              &
          '        SINTABLE =',I2,', TAUWBUG =',I2,               &
-         ', VISCSTRESS =',F8.5,', SINTAIL1 =',F8.5,', SINTAIL2 =',F8.5,'  /')
+         ', VISCSTRESS =',F8.5,', SINTAIL1 =',F8.5,', SINTAIL2 =',F8.5,',' / &
+         ', CAPCHA =',F8.5,', CHAMIN =',F8.5,', CHA0 =',F8.5,', UCAP =',F5.1,', SIGMAUCAP =', &
+         F5.1,'  /')
 #endif
     !
 #ifdef W3_ST6

diff --git a/model/src/w3iogrmd.F90 b/model/src/w3iogrmd.F90
@@ -1508,7 +1508,7 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
            SSDSPBK, SSDSBINT, FFXPM, FFXFM, FFXFA,           &
            SSDSHCK,                                          &
            IKTAB, DCKI, QBI, SATINDICES, SATWEIGHTS,         &
-           DIKCUMUL, CUMULW, SINTAILPAR
+           DIKCUMUL, CUMULW, SINTAILPAR, CAPCHNK
 #ifdef W3_ASCII
       WRITE (NDSA,*)                                         &
            'ZZWND, AALPHA, ZZ0MAX, BBETA, SSINTHP, ZZALP,    &
@@ -1521,7 +1521,7 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
            SSDSPBK, SSDSBINT, FFXPM, FFXFM, FFXFA,           &
            SSDSHCK,                                          &
            IKTAB, DCKI, QBI, SATINDICES, SATWEIGHTS,         &
-           DIKCUMUL, CUMULW, SINTAILPAR:',                   &
+           DIKCUMUL, CUMULW, SINTAILPAR, CAPCHNK:',          &
            ZZWND, AALPHA, ZZ0MAX, BBETA, SSINTHP, ZZALP,     &
            TTAUWSHELTER, SSWELLFPAR, SSWELLF, SSINBR,        &
            ZZ0RAT, SSDSC,                                    &
@@ -1532,7 +1532,7 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
            SSDSPBK, SSDSBINT, FFXPM, FFXFM, FFXFA,           &
            SSDSHCK,                                          &
            IKTAB, DCKI, QBI, SATINDICES, SATWEIGHTS,         &
-           DIKCUMUL, CUMULW, SINTAILPAR
+           DIKCUMUL, CUMULW, SINTAILPAR, CAPCHNK
 #endif
       IF (SINTAILPAR(1).GT.0.5) THEN
         WRITE (NDSM) DELUST, DELTAIL, DELTAUW, DELU, DELALP, &
@@ -1558,7 +1558,7 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT &
            SSDSPBK, SSDSBINT, FFXPM, FFXFM, FFXFA,           &
            SSDSHCK,                                          &
            IKTAB, DCKI, QBI, SATINDICES, SATWEIGHTS,         &
-           DIKCUMUL, CUMULW, SINTAILPAR
+           DIKCUMUL, CUMULW, SINTAILPAR, CAPCHNK
       IF (SINTAILPAR(1).GT.0.5) THEN
         CALL INSIN4(.FALSE.)
         READ (NDSM,END=801,ERR=802,IOSTAT=IERR)              &

diff --git a/model/src/w3src4md.F90 b/model/src/w3src4md.F90
@@ -1849,6 +1849,8 @@ SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN)
     !/    14-Aug-2006 : Modified following Bidlot           ( version 2.22-SHOM )
     !/    18-Aug-2006 : Ported to version 3.09
     !/    03-Apr-2010 : Adding output of Charnock parameter ( version 3.14-IFREMER )
+    !/    03-May-2024 : Optional functional form of         ( version 7.15 )
+    !/                  Charnock coefficient and surface drag (UK Met Office).
     !
     !  1. Purpose :
     !
@@ -1895,7 +1897,7 @@ SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN)
     ! 10. Source code :
     !-----------------------------------------------------------------------------!
     USE CONSTANTS, ONLY: GRAV, KAPPA, NU_AIR
-    USE W3GDATMD,  ONLY: ZZWND, AALPHA, ZZ0MAX, SINTAILPAR
+    USE W3GDATMD,  ONLY: ZZWND, AALPHA, ZZ0MAX, SINTAILPAR, CAPCHNK
 #ifdef W3_T
     USE W3ODATMD, ONLY: NDST
 #endif
@@ -1908,6 +1910,7 @@ SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN)
     INTEGER          :: IND,J
     REAL             :: TAUW_LOCAL
     REAL             :: TAUOLD,CDRAG,WCD,USTOLD,X,UST,ZZ0,ZNU,ZZ00,F,DELF
+    REAL             :: CHATH, XMIN ! used for reduction of high winds
     INTEGER, PARAMETER      :: NITER=10
     REAL   , PARAMETER      :: XM=0.50, EPS1=0.00001
     INTEGER                 :: ITER
@@ -1918,6 +1921,7 @@ SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN)
     !      *EPS1*      REAL      SMALL NUMBER TO MAKE SURE THAT A SOLUTION
     !                            IS OBTAINED IN ITERATION WITH TAU>TAUW.
 
+    CHATH = AALPHA
     !
     IF (SINTAILPAR(1).GT.0.5) THEN
       TAUW_LOCAL=MAX(MIN(TAUW,TAUWMAX),0.)
@@ -1932,6 +1936,17 @@ SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN)
       USTAR=(TAUT(IND,J)*DELI2+TAUT(IND+1,J  )*DELI1)*DELJ2 &
            + (TAUT(IND,J+1)*DELI2+TAUT(IND+1,J+1)*DELI1)*DELJ1
     ELSE
+      IF (CAPCHNK(1).EQ.1.) THEN
+        ! Computation of sea surface roughness and charnock coefficient based
+        ! on Donelan (2018). Determines minimum charnock; reduction for winds
+        ! above a particular threshold
+        CHATH  = CAPCHNK(2) + 0.5 * (CAPCHNK(3) - CAPCHNK(2)) * (1 & 
+                 - TANH((WINDSPEED-CAPCHNK(4))/CAPCHNK(5)))
+        XMIN   = 0.15 * (CAPCHNK(3)-CHATH)
+      ELSE
+        XMIN  = 0.
+      END IF
+
       ! This max is for comparison ... to be removed later
       !        TAUW_LOCAL=MAX(MIN(TAUW,TAUWMAX),0.)
       TAUW_LOCAL=TAUW
@@ -1941,9 +1956,9 @@ SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN)
       TAUOLD  = MAX(USTOLD**2, TAUW_LOCAL+EPS1)
       ! Newton method to solve for ustar in U=ustar*log(Z/Z0)
       DO ITER=1,NITER
-        X   = TAUW_LOCAL/TAUOLD
+        X   = MAX(TAUW_LOCAL/TAUOLD, XMIN)
         UST = SQRT(TAUOLD)
-        ZZ00=AALPHA*TAUOLD/GRAV
+        ZZ00 = CHATH*TAUOLD/GRAV
         IF (ZZ0MAX.NE.0) ZZ00=MIN(ZZ00,ZZ0MAX)
         ! Corrects roughness ZZ00 for quasi-linear effect
         ZZ0 = ZZ00/(1.-X)**XM
@@ -1969,10 +1984,16 @@ SUBROUTINE CALC_USTAR(WINDSPEED,TAUW,USTAR,Z0,CHARN)
         SQRTCDM1  = MIN(WINDSPEED/USTAR,100.0)
         Z0  = ZZWND*EXP(-KAPPA*SQRTCDM1)
       ELSE
-        Z0 = AALPHA*0.001*0.001/GRAV
+        Z0 = CHATH*0.001*0.001/GRAV
       END IF
-      CHARN = AALPHA
+      CHARN = CHATH
     END IF
+    IF(CAPCHNK(1) .EQ. 1) THEN
+      ! Problem with large values of CHARN for low winds
+      CHARN = MIN( 0.09 , CHARN )
+      IF(CHARN.LT.CHATH) CHARN = CHATH
+    ENDIF
+
     !  WRITE(6,*) 'CALC_USTAR:',WINDSPEED,TAUW,AALPHA,CHARN,Z0,USTAR
     !
     RETURN

diff --git a/regtests/bin/matrix.base b/regtests/bin/matrix.base
@@ -918,6 +918,7 @@
     echo "$rtst -g ST4_T702 -w work_T702    -i input_10ms -N   $ww3 ww3_ts1" >> matrix.body
     echo "$rtst -g ST4_T707 -w work_T707GQM -i input_10ms -N  $ww3 ww3_ts1" >> matrix.body
     echo "$rtst -g ST4_T713 -w work_T713GQM -i input_10ms -N  $ww3 ww3_ts1" >> matrix.body
+    echo "$rtst -s ST4      -w work_ST4_CAP -i input_capcha -N  $ww3 ww3_ts1" >> matrix.body
   fi
 
 # fetch limited growth, no switch sharing here

diff --git a/regtests/ww3_ts1/input_capcha/namelists_CAPCHA.nml b/regtests/ww3_ts1/input_capcha/namelists_CAPCHA.nml
@@ -0,0 +1,12 @@
+&SIN4
+  SINTABLE=0,
+  BETAMAX=1.35,
+  CAPCHA=1,
+  CHA0=0.0095,
+  SIGMAUCAP=5.0,
+  SINTAIL1=0.8,
+  TAUWBUG = 1,
+  UCAP=28.0
+/
+
+END OF NAMELISTS
diff --git a/regtests/ww3_ts1/input_capcha/switch_ST4 b/regtests/ww3_ts1/input_capcha/switch_ST4
@@ -0,0 +1 @@
+NOGRB SHRD PR0 FLX0 LN1 ST4 NL1 BT1 DB1 TR0 BS0 IC0 IS0 REF0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7 O10 O11
diff --git a/regtests/ww3_ts1/input_capcha/ww3_grid.nml b/regtests/ww3_ts1/input_capcha/ww3_grid.nml
@@ -0,0 +1,70 @@
+! -------------------------------------------------------------------- !
+! WAVEWATCH III - ww3_grid.nml - Grid pre-processing                   !
+! -------------------------------------------------------------------- !
+
+! -------------------------------------------------------------------- !
+! Define the spectrum parameterization via SPECTRUM_NML namelist
+! -------------------------------------------------------------------- !
+&SPECTRUM_NML
+  SPECTRUM%XFR       =  1.10
+  SPECTRUM%FREQ1     =  0.0485
+  SPECTRUM%NK        =  36
+  SPECTRUM%NTH       =  24
+/
+
+! -------------------------------------------------------------------- !
+! Define the run parameterization via RUN_NML namelist
+! -------------------------------------------------------------------- !
+&RUN_NML
+  RUN%FLSOU        =  T
+/
+
+! -------------------------------------------------------------------- !
+! Define the timesteps parameterization via TIMESTEPS_NML namelist
+! -------------------------------------------------------------------- !
+&TIMESTEPS_NML
+  TIMESTEPS%DTMAX        =  900.
+  TIMESTEPS%DTXY         =  900.
+  TIMESTEPS%DTKTH        =  900.
+  TIMESTEPS%DTMIN        =  15.
+/
+
+! -------------------------------------------------------------------- !
+! Define the grid to preprocess via GRID_NML namelist
+! -------------------------------------------------------------------- !
+&GRID_NML
+  GRID%NAME         =  'HOMOGENEOUS SOURCE TERM TEST'
+  GRID%NML          =  '../input_capcha/namelists_CAPCHA.nml'
+  GRID%TYPE         =  'RECT'
+  GRID%COORD        =  'SPHE'
+  GRID%CLOS         =  'NONE'
+  GRID%ZLIM         =  -5.
+  GRID%DMIN         =  5.75
+/
+
+! -------------------------------------------------------------------- !
+! Define the rectilinear grid type via RECT_NML namelist
+! -------------------------------------------------------------------- !
+&RECT_NML
+  RECT%NX           =  3
+  RECT%NY           =  3
+  RECT%SX           =  1.
+  RECT%SY           =  1.
+  RECT%SF           =  1.E-2
+  RECT%X0           =  -1.
+  RECT%Y0           =  -1.
+  RECT%SF0          =  1.E-2
+/
+
+! -------------------------------------------------------------------- !
+! Define the depth to preprocess via DEPTH_NML namelist
+! -------------------------------------------------------------------- !
+&DEPTH_NML
+  DEPTH%SF        =  -2500.
+  DEPTH%FILENAME  =  '../input/HOMOGENEOUS.depth'
+  DEPTH%IDLA      =  3
+/
+
+! -------------------------------------------------------------------- !
+! WAVEWATCH III - end of namelist                                      !
+! -------------------------------------------------------------------- !
diff --git a/regtests/ww3_ts1/input_capcha/ww3_ounf.nml b/regtests/ww3_ts1/input_capcha/ww3_ounf.nml
@@ -0,0 +1,29 @@
+! -------------------------------------------------------------------- !
+! WAVEWATCH III ww3_ounf.nml - Grid output post-processing             !
+! -------------------------------------------------------------------- !
+
+! -------------------------------------------------------------------- !
+! Define the output fields to postprocess via FIELD_NML namelist
+! -------------------------------------------------------------------- !
+&FIELD_NML
+  FIELD%TIMESTART        =  '19680101 120000'
+  FIELD%TIMESTRIDE       =  '10'
+  FIELD%TIMECOUNT        =  '8000'
+  FIELD%LIST             =  'DPT WND ICE HS MSS MSD FAW WCC WCF WCH WCM FOC TAW CHA FWS WBT'
+  FIELD%PARTITION        =  '0 1 2'
+  FIELD%TYPE             =  4
+/
+
+! -------------------------------------------------------------------- !
+! Define the content of the output file via FILE_NML namelist
+! -------------------------------------------------------------------- !
+&FILE_NML
+  FILE%IX0           = 2
+  FILE%IXN           = 2
+  FILE%IY0           = 2
+  FILE%IYN           = 2
+/
+
+! -------------------------------------------------------------------- !
+! WAVEWATCH III - end of namelist                                      !
+! -------------------------------------------------------------------- !