From @KarinaAsmar-NOAA and @AndrewHazelton:

Add new gust factor to HAFS products NOAA-EMC#1422
The CIMAS/AOML/HRD team has developed a new gust factor as part of HAFS
development. This PR implements the new calculations into the UPP source
code and adds the new GUST ID to the HAFS products list.

Author: BinLiu-NOAA

parm/hafs/postcntrl_hafs_nosat.xml

@@ -554,6 +554,11 @@
        <scale>3.0</scale>
        </param>
 
+       <param>
+       <shortname>HAFS_GUST_ON_SPEC_ALT_ABOVE_MEAN_SEA_LVL</shortname>
+       <scale>3.0</scale>
+       </param>
+
        <param>
        <shortname>MAX_MAXUVV_ON_ISOBARIC_SFC_100-1000hpa</shortname>
        <pname>MAXUVV</pname>

parm/hafs/postxconfig-NT-hafs_nosat.txt

@@ -1,5 +1,5 @@
 1
-84
+85
 HURPRS
 32769
 ncep_nco
@@ -3376,6 +3376,48 @@ surface
 ?
 ?
 ?
+1026
+HAFS_GUST_ON_SPEC_ALT_ABOVE_MEAN_SEA_LVL
+?
+1
+tmpl4_0
+GUST
+?
+?
+spec_hgt_lvl_above_grnd
+0
+?
+1
+10.
+?
+0
+?
+0
+?
+?
+?
+?
+0
+0.0
+0
+0.0
+?
+0
+0.0
+0
+0.0
+0
+0.0
+0
+0.0
+1
+3.0
+0
+0
+0
+?
+?
+?
 423
 MAX_MAXUVV_ON_ISOBARIC_SFC_100-1000hpa
 hourly maximum Upward Vertical Velocity between 100-1000hpa

parm/post_avblflds.xml

@@ -8660,5 +8660,14 @@
        <scale>3.0</scale>
     </param>
 
+    <param>
+       <post_avblfldidx>1026</post_avblfldidx>
+       <shortname>HAFS_GUST_ON_SPEC_ALT_ABOVE_MEAN_SEA_LVL</shortname>
+       <pname>GUST</pname>
+       <fixed_sfc1_type>spec_hgt_lvl_above_grnd</fixed_sfc1_type>
+       <level>10.</level>
+       <scale>3.0</scale>
+    </param>
+
   </post_avblflds>
 </postxml>

sorc/ncep_post.fd/CALGUSTCONV.f

@@ -0,0 +1,82 @@
+!> @file
+!> @brief Subroutine that computes wind gust from turbulence and convective components
+!> for HAFS model.
+!>
+!> This subroutine is based on the Python code provided by Andrew Hazelton (HRD),
+!> which computes a new wind gust factor based on turbulence and convective components.
+!>
+!> ### Program history log:
+!> Date | Programmer | Comments
+!> -----|------------|---------
+!> 2026-01-06 | Karina Asmar | Initial
+!>
+!> @author Karina Asmar NCEP/EMC @date 2026-01-06
+! ------------------------------------------------------------------------------------------
+!>
+!> @param[in] SPEED850 Wind speed at 850 mb.
+!> @param[in] SPEED950 Wind speed at 950 mb.
+!> @param[inout] GUSTCONV Gust factor.
+!>
+      SUBROUTINE CALGUSTCONV(SPEED850,SPEED950,GUSTCONV)
+!
+!
+      use vrbls2d , only: u10,v10, ustar
+      use ctlblk_mod, only: ista, iend, jsta, jend, ista_2l, iend_2u, jsta_2l, jend_2u, spval
+!
+      implicit none
+!
+      INCLUDE "mpif.h"
+!
+!     INCLUDE ETA GRID DIMENSIONS.  SET/DERIVE PARAMETERS.
+!
+!     DECLARE VARIABLES.
+!
+      REAL,intent(in)    :: SPEED850(ista_2l:iend_2u,jsta_2l:jend_2u)
+      REAL,intent(in)    :: SPEED950(ista_2l:iend_2u,jsta_2l:jend_2u)
+      REAL,intent(inout) :: GUSTCONV(ista_2l:iend_2u,jsta_2l:jend_2u)
+!
+      REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: WIND10, WSD, GUST1_NEW, GUST2_NEW, WSTT1_NEW, &
+                                                            WSTT2_NEW, GF1_NEW, GF2_NEW
+!
+!
+      integer I,J
+!
+!
+!*****************************************************************************
+!> CALGUSTCONV computes new gust wind factor for HAFS based on convective and
+!> turbulent components.
+!     START CALGUSTCONV HERE.
+!
+!
+!$omp parallel do private(i,j)
+      DO J=JSTA,JEND
+        DO I=ISTA,IEND
+          GUSTCONV(I,J) = SPVAL
+        ENDDO
+      ENDDO
+
+! CALCULATE GUST FACTORS
+!$omp parallel do private(i,j)
+     DO J=JSTA,JEND
+       DO I=ISTA,IEND
+       IF(U10(I,J)<SPVAL.AND.V10(I,J)<SPVAL.AND.SPEED850(I,J)<SPVAL.AND.SPEED950(I,J)<SPVAL) THEN
+       ! CALCULATE WINDS AND WIND SPEED DIFFERENCE
+       WIND10(I,J) = SQRT(U10(I,J)**2 + V10(I,J)**2)
+       WSD(I,J) = SPEED850(I,J) - SPEED950(I,J)
+       IF (WSD(I,J) < 0.0) WSD(I,J) = 0.0      ! MAKE SURE WIND SPEED IS >= ZERO
+
+       ! CALCULATE NEW GUST FACTOR
+       GUST1_NEW(I,J) = 3.0 * USTAR(I,J)  ! TURBULENT COMPONENT
+       WSTT1_NEW(I,J) = WIND10(I,J) + GUST1_NEW(I,J)
+       GF1_NEW(I,J) = WSTT1_NEW(I,J) / WIND10(I,J)
+       GUST2_NEW(I,J) = 0.3 * WSD(I,J)  ! CONVECTIVE MIXING COMPONENT
+
+       GUSTCONV(I,J) = (WSTT1_NEW(I,J) + GUST2_NEW(I,J))
+       ENDIF
+       ENDDO
+     ENDDO
+!
+!     END OF ROUTINE.
+!
+      RETURN
+      END

sorc/ncep_post.fd/CMakeLists.txt

@@ -7,6 +7,7 @@ list(APPEND LIB_SRC
     CALDRG.f
     CALDWP.f
     CALGUST.f
+    CALGUSTCONV.f
     CALHEL.f
     CALHEL2.f
     CALHEL3.f

sorc/ncep_post.fd/MDL2P.f

@@ -39,12 +39,13 @@
 !> 2023-08-24 | Y Mao           | Add gtg_on option for GTG interpolation
 !> 2023-09-12 | J Kenyon        | Prevent spurious supercooled rain and cloud water
 !> 2024-04-23 | E James         | Adding smoke emissions (ebb) from RRFS
-!> 2024-09-23 | K Asmar		| Add velocity potential and streamfunction from wind vectors
+!> 2024-09-23 | K Asmar		    | Add velocity potential and streamfunction from wind vectors
 !> 2024-12-12 | J Meng          | Adding UUtah 2024 SLR algorithm
 !> 2025-01-17 | J Kenyon        | Add graupel number concentration (QQNG)
 !> 2025-11-13 | L Pan           | enable aerosols to be output on isobaric surfaces
 !> 2025-11-17 | W Meng          | Correct variable allocation
 !> 2025-11-19 | W Meng          | Relocate dxm calculation
+!> 2026-01-06 | K Asmar			| Add convective wind gust calculation
 !>
 !> @author T Black W/NP2 @date 1999-09-23
 !--------------------------------------------------------------------------------------
@@ -114,6 +115,7 @@ SUBROUTINE MDL2P(iostatusD3D)
      &,                      SALTSL(:,:,:), SOOTSL(:,:,:), WASOSL(:,:,:)           &
      &,                      SUSOSL(:,:,:)
       REAL, allocatable :: GTGSL(:,:),CATSL(:,:),MWTSL(:,:)
+      REAL, allocatable :: WS850(:,:), WS950(:,:), GUSTCONV(:,:)
 !
       integer,intent(in) :: iostatusD3D
       INTEGER, dimension(ista_2l:iend_2u,jsta_2l:jend_2u)  :: NL1X, NL1XF
@@ -277,6 +279,10 @@ SUBROUTINE MDL2P(iostatusD3D)
         if (.not. allocated(CATSL)) allocate(CATSL(ista_2l:iend_2u,jsta_2l:jend_2u))
         if (.not. allocated(MWTSL)) allocate(MWTSL(ista_2l:iend_2u,jsta_2l:jend_2u))
       endif
+! Allocate HAFS wind speeds and gust factor
+      if (.not. allocated(ws850)) allocate(ws850(ista_2l:iend_2u,jsta_2l:jend_2u))
+      if (.not. allocated(ws950)) allocate(ws950(ista_2l:iend_2u,jsta_2l:jend_2u))
+      if (.not. allocated(gustconv)) allocate(gustconv(ista_2l:iend_2u,jsta_2l:jend_2u))
 !     
 !     SET TOTAL NUMBER OF POINTS ON OUTPUT GRID.
 !
@@ -301,7 +307,8 @@ SUBROUTINE MDL2P(iostatusD3D)
          (IGET(257) > 0) .OR. (IGET(258) > 0) .OR.      &
          (IGET(294) > 0) .OR. (IGET(268) > 0) .OR.      &
          (IGET(331) > 0) .OR. (IGET(326) > 0) .OR.      &
-	 (IGET(1021) > 0) .OR. (IGET(1022) > 0) .OR.	&
+	     (IGET(1021) > 0) .OR. (IGET(1022) > 0) .OR.	&
+		 (IGET(1026) > 0) .OR.							&
 ! add D3D fields
          (IGET(354) > 0) .OR. (IGET(355) > 0) .OR.      &
          (IGET(356) > 0) .OR. (IGET(357) > 0) .OR.      &
@@ -350,6 +357,10 @@ SUBROUTINE MDL2P(iostatusD3D)
         if(gridtype == 'B' .or. gridtype == 'E')                         &
           call exch(PINT(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LP1)) 
 
+! wind speeds at 850 and 950 mb for HAFS wind gust
+        WS850 = SPVAL
+        WS950 = SPVAL
+
         DO LP=1,LSM
 
 !         if(me == 0) print *,'in LP loop me=',me,'UH=',UH(1:10,JSTA,LP), &
@@ -1864,11 +1875,11 @@ SUBROUTINE MDL2P(iostatusD3D)
 !
         IF(IGET(018) > 0.OR.IGET(019) > 0)THEN
           log1=.false.
-          IF(IGET(018) > 0.) then
-             if(LVLS(LP,IGET(018)) > 0 ) log1=.true.
+          IF(IGET(018) > 0) then
+             if(LVLS(LP,IGET(018)) > 0) log1=.true.
           endif
           IF(IGET(019) > 0) then
-             if(LVLS(LP,IGET(019)) > 0 ) log1=.true.
+             if(LVLS(LP,IGET(019)) > 0) log1=.true.
           endif
           if ( log1 ) then
 !$omp  parallel do private(i,j)
@@ -1917,6 +1928,19 @@ SUBROUTINE MDL2P(iostatusD3D)
                 enddo
               enddo
             endif
+
+            ! *** WIND SPEED AT 850 AND 950 MB FOR HAFS WIND GUST
+            IF(IGET(1026)>0)THEN
+              DO J=JSTA,JEND
+              DO I=ISTA,IEND
+                IF(GRID1(I,J)<SPVAL .and.  GRID2(I,J)<SPVAL) then
+                  IF (ABS(SPL(LP)-85000.)<SMALL) WS850(I,J) = SQRT(GRID1(I,J)**2 + GRID2(I,J)**2)
+                  IF (ABS(SPL(LP)-95000.)<SMALL) WS950(I,J) = SQRT(GRID1(I,J)**2 + GRID2(I,J)**2)
+                ENDIF
+              ENDDO
+              ENDDO
+            ENDIF
+
           ENDIF
         ENDIF
 !     
@@ -4967,6 +4991,31 @@ SUBROUTINE MDL2P(iostatusD3D)
             enddo
          endif
       ENDIF
+
+!
+! *** HAFS WIND GUST
+!
+      IF (IGET(1026) > 0) THEN
+ 		CALL CALGUSTCONV(WS850,WS950,GUSTCONV)
+!$omp  parallel do private(i,j)
+         DO J=JSTA,JEND
+           DO I=ISTA,IEND
+             GRID1(I,J) = GUSTCONV(I,J)
+           ENDDO
+         ENDDO
+		if(grib=='grib2') then
+        cfld=cfld+1
+        fld_info(cfld)%ifld=IAVBLFLD(IGET(1026))
+!$omp parallel do private(i,j,ii,jj)
+        do j=1,jend-jsta+1
+          jj = jsta+j-1
+          do i=1,iend-ista+1
+            ii = ista+i-1
+            datapd(i,j,cfld) = GRID1(ii,jj)
+          enddo
+        enddo
+       endif
+      ENDIF
 !
 if(allocated(d3dsl))   deallocate(d3dsl)
 if(allocated(smokesl)) deallocate(smokesl)
@@ -4982,6 +5031,10 @@ SUBROUTINE MDL2P(iostatusD3D)
 if(allocated(GTGSL)) deallocate(GTGSL)
 if(allocated(CATSL)) deallocate(CATSL)
 if(allocated(MWTSL)) deallocate(MWTSL)
+! HAFS winds and gust factor
+if(allocated(WS850)) deallocate(WS850)
+if(allocated(WS950)) deallocate(WS950)
+if(allocated(GUSTCONV)) deallocate(GUSTCONV)
 !     END OF ROUTINE.
 !
       RETURN