From 64760649a49df4323778ed7a401c69b21969a45c Mon Sep 17 00:00:00 2001
From: ericvmueller <ericvmueller@gmail.com>
Date: Thu, 30 Jul 2026 13:54:52 -0400
Subject: [PATCH] FDS Source: scale SEM eddy contributions using all mesh
 segments of VENT to preserve mean flow and field continuity

---
 Source/turb.f90 | 52 ++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 49 insertions(+), 3 deletions(-)

diff --git a/Source/turb.f90 b/Source/turb.f90
index 29f15e8add..4056527a55 100644
--- a/Source/turb.f90
+++ b/Source/turb.f90
@@ -1692,15 +1692,16 @@ SUBROUTINE SYNTHETIC_TURBULENCE(DT,T)
 
 USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
 USE PHYSICAL_FUNCTIONS, ONLY: GET_WIND_AT_HEIGHT
+USE MPI_F08
 
 REAL(EB), INTENT(IN) :: DT,T
 INTEGER :: NE,NV,NT,II,JJ,KK,NM
-INTEGER :: IN,IA,IC
+INTEGER :: IN,IA,IC,IERR
 TYPE(VENTS_TYPE), POINTER :: VT
 TYPE(SURFACE_TYPE), POINTER :: SF
 REAL(EB) :: RAMP_T,TSI,SIGMA_MAX,ADV_VEL(3),SGN,UW,VW,WW,U1,V1,W1
-REAL(EB) :: U_ADD,V_ADD,W_ADD
-REAL(EB), ALLOCATABLE :: VENT_VEL(:,:)
+REAL(EB) :: U_ADD,V_ADD,W_ADD,U_MEAN
+REAL(EB), ALLOCATABLE :: VENT_VEL(:,:),VENT_MEAN_SUM(:),VENT_MEAN_CNT(:)
 LOGICAL, ALLOCATABLE :: VENT_VEL_SET(:)
 
 ! Reference:
@@ -1860,6 +1861,51 @@ SUBROUTINE SYNTHETIC_TURBULENCE(DT,T)
    ENDDO VENT_APPLY_LOOP
 ENDDO MESH_APPLY_LOOP
 
+! Subtract the vent mean of the normal component so fluctuations do not
+! change volume flow. Must be global over TOTAL_INDEX (all mesh sections) 
+! to avod jumps at mesh interfaces
+ALLOCATE(VENT_MEAN_SUM(N_VENT_TOTAL),SOURCE=0._EB)
+ALLOCATE(VENT_MEAN_CNT(N_VENT_TOTAL),SOURCE=0._EB)
+DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
+   CALL POINT_TO_MESH(NM)
+   DO NV=1,N_VENT
+      VT => VENTS(NV)
+      IF (VT%N_EDDY==0) CYCLE
+      NT = VT%TOTAL_INDEX
+      SELECT CASE (ABS(VT%IOR))
+         CASE(1)
+            VENT_MEAN_SUM(NT) = VENT_MEAN_SUM(NT) + SUM(VT%U_EDDY)
+            VENT_MEAN_CNT(NT) = VENT_MEAN_CNT(NT) + REAL(SIZE(VT%U_EDDY),EB)
+         CASE(2)
+            VENT_MEAN_SUM(NT) = VENT_MEAN_SUM(NT) + SUM(VT%V_EDDY)
+            VENT_MEAN_CNT(NT) = VENT_MEAN_CNT(NT) + REAL(SIZE(VT%V_EDDY),EB)
+         CASE(3)
+            VENT_MEAN_SUM(NT) = VENT_MEAN_SUM(NT) + SUM(VT%W_EDDY)
+            VENT_MEAN_CNT(NT) = VENT_MEAN_CNT(NT) + REAL(SIZE(VT%W_EDDY),EB)
+      END SELECT
+   ENDDO
+ENDDO
+IF (N_MPI_PROCESSES>1) THEN
+   CALL MPI_ALLREDUCE(MPI_IN_PLACE,VENT_MEAN_SUM,N_VENT_TOTAL,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
+   CALL MPI_ALLREDUCE(MPI_IN_PLACE,VENT_MEAN_CNT,N_VENT_TOTAL,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
+ENDIF
+DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
+   CALL POINT_TO_MESH(NM)
+   DO NV=1,N_VENT
+      VT => VENTS(NV)
+      IF (VT%N_EDDY==0) CYCLE
+      NT = VT%TOTAL_INDEX
+      IF (VENT_MEAN_CNT(NT)<=TWO_EPSILON_EB) CYCLE
+      U_MEAN = VENT_MEAN_SUM(NT)/VENT_MEAN_CNT(NT)
+      SELECT CASE (ABS(VT%IOR))
+         CASE(1); VT%U_EDDY = VT%U_EDDY - U_MEAN
+         CASE(2); VT%V_EDDY = VT%V_EDDY - U_MEAN
+         CASE(3); VT%W_EDDY = VT%W_EDDY - U_MEAN
+      END SELECT
+   ENDDO
+ENDDO
+DEALLOCATE(VENT_MEAN_SUM,VENT_MEAN_CNT)
+
 END SUBROUTINE SYNTHETIC_TURBULENCE