Merge pull request #13881 from mcgratta/master

FDS Source: Simplify MPI exchange of global variables

Author: mcgratta

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -14404,9 +14404,8 @@ \section{Diagnostic Output ({\tt .out})}
 Min divergence: -0.20E+00 at ( 66, 13,  1)
 Max VN number:   0.51E+00 at (  1, 25, 18)
 No. of Lagrangian Particles:            27
-Radiation Loss to Boundaries:        13.830 kW
 \end{lstlisting}
-The \ct{Time Step} indicates the total number of iterations. The date and time indicate the current wall clock time. The \ct{STEP SIZE} indicates the size of the numerical time step. The \ct{Total Time} indicates the total simulation time calculated up to that point. The \ct{Pressure Iterations} are the number of iterations of the pressure solver for the corrector (second) half of the time step. The pressure solver iterations are designed to minimize the error in the normal component of velocity at solid walls or the interface of two meshes. The \ct{Maximum Velocity Error} indicates this error and in which grid cell it occurs.  \ct{Max/Min divergence} is the max/min value of the function $\nabla \cdot \bu$ and is used as a diagnostic when the flow is incompressible (i.e., no heating); \ct{Max CFL number} is the maximum value of the CFL number, the primary time step constraint; \ct{Max VN number} is the maximum value of the Von Neumann number, the secondary time step constraint. The \ct{No. of Lagrangian Particles} refers to the number of particles in the current mesh. The \ct{Radiation Loss to Boundaries} is the amount of energy that is being radiated to the boundaries. As compartments heat up, the energy lost to the boundaries can grow to be an appreciable fraction of the \ct{Total Heat Release Rate}.
+The \ct{Time Step} indicates the total number of iterations. The date and time indicate the current wall clock time. The \ct{STEP SIZE} indicates the size of the numerical time step. The \ct{Total Time} indicates the total simulation time calculated up to that point. The \ct{Pressure Iterations} are the number of iterations of the pressure solver for the corrector (second) half of the time step. The pressure solver iterations are designed to minimize the error in the normal component of velocity at solid walls or the interface of two meshes. The \ct{Maximum Velocity Error} indicates this error and in which grid cell it occurs.  \ct{Max/Min divergence} is the max/min value of the function $\nabla \cdot \bu$ and is used as a diagnostic when the flow is incompressible (i.e., no heating); \ct{Max CFL number} is the maximum value of the CFL number, the primary time step constraint; \ct{Max VN number} is the maximum value of the Von Neumann number, the secondary time step constraint. The \ct{No. of Lagrangian Particles} refers to the number of particles in the current mesh.
 
 Following the completion of a successful run, a summary of the CPU usage per subroutine is listed in the file called \ct{CHID_cpu.csv} (Section~\ref{out:CPU}). This is useful in determining where most of the computational effort is being placed.
 

Source/data.f90

@@ -16,7 +16,7 @@ MODULE OUTPUT_DATA
 INTEGER :: HVAC_SMV_EQUIVALENCE(300:350)
 REAL(EB) :: T_LAST_DUMP_HRR,T_LAST_DUMP_MASS,T_LAST_DUMP_MOM
 REAL(EB),ALLOCATABLE, DIMENSION(:) :: ENTHALPY_SUM
-REAL(EB),ALLOCATABLE, DIMENSION(:,:) :: MASS_DT,Q_DOT,Q_DOT_SUM,M_DOT,M_DOT_SUM
+REAL(EB),ALLOCATABLE, DIMENSION(:) :: MASS_DT,Q_DOT,Q_DOT_SUM,M_DOT,M_DOT_SUM
 
 TYPE SMOKE3D_TYPE
    INTEGER :: QUANTITY_INDEX,Y_INDEX=0,Z_INDEX=0

Source/dump.f90

@@ -3446,7 +3446,7 @@ SUBROUTINE DUMP_RESTART(T,DT,NM)
              RTE_SOURCE_CORRECTION_FACTOR,RAD_Q_SUM,KFST4_SUM,ENTHALPY_SUM(NM)
 WRITE(LU_CORE(NM)) DT_BNDF,DT_CPU,DT_CTRL,DT_DEVC,DT_FLUSH,DT_GEOM,DT_HRR,DT_ISOF,DT_MASS,DT_PART,DT_PL3D,DT_PROF,DT_RADF,&
                    DT_SLCF,DT_SL3D,DT_SMOKE3D,DT_UVW
-WRITE(LU_CORE(NM)) Q_DOT_SUM(1:N_Q_DOT,NM),M_DOT_SUM(1:N_TRACKED_SPECIES,NM),MASS_DT(0:N_SPECIES+N_TRACKED_SPECIES,NM)
+WRITE(LU_CORE(NM)) Q_DOT_SUM(1:N_Q_DOT),M_DOT_SUM(1:N_TRACKED_SPECIES),MASS_DT(0:N_SPECIES+N_TRACKED_SPECIES)
 DO N=1,N_DEVC
    DV => DEVICE(N)
    WRITE(LU_CORE(NM)) DV%T,DV%T_CHANGE,DV%TMP_L,DV%Y_C,DV%CURRENT_STATE,DV%PRIOR_STATE,&
@@ -3627,7 +3627,7 @@ SUBROUTINE READ_RESTART(T,DT,NM)
 
 READ(LU_RESTART(NM)) DT_BNDF,DT_CPU,DT_CTRL,DT_DEVC,DT_FLUSH,DT_GEOM,DT_HRR,DT_ISOF,DT_MASS,DT_PART,DT_PL3D,DT_PROF,DT_RADF,&
                      DT_SLCF,DT_SL3D,DT_SMOKE3D,DT_UVW
-READ(LU_RESTART(NM)) Q_DOT_SUM(1:N_Q_DOT,NM),M_DOT_SUM(1:N_TRACKED_SPECIES,NM),MASS_DT(0:N_SPECIES+N_TRACKED_SPECIES,NM)
+READ(LU_RESTART(NM)) Q_DOT_SUM(1:N_Q_DOT),M_DOT_SUM(1:N_TRACKED_SPECIES),MASS_DT(0:N_SPECIES+N_TRACKED_SPECIES)
 DO N=1,N_DEVC
    DV => DEVICE(N)
    READ(LU_RESTART(NM)) DV%T,DV%T_CHANGE,DV%TMP_L,DV%Y_C,DV%CURRENT_STATE,DV%PRIOR_STATE,&
@@ -3922,8 +3922,6 @@ SUBROUTINE WRITE_DIAGNOSTICS(T,DT)
    IF (CHECK_POISSON) WRITE(LU_OUTPUT,'(A,E9.2)') '       Poisson Error : ',M%POIS_ERR
    IF (SIM_MODE==DNS_MODE .OR. CHECK_VN) WRITE(LU_OUTPUT,230) M%VN,M%I_VN,M%J_VN,M%K_VN
    IF (M%NLP>0) WRITE(LU_OUTPUT,141) M%NLP
-   IF (ABS(Q_DOT(1,NM))>1._EB) WRITE(LU_OUTPUT,119) Q_DOT(1,NM)/1000._EB
-   IF (ABS(Q_DOT(3,NM))>1._EB) WRITE(LU_OUTPUT,120) Q_DOT(3,NM)/1000._EB
    IF (M%DT_RESTRICT_STORE>0 ) THEN
       WRITE(LU_OUTPUT,121) M%DT_RESTRICT_STORE
       M%DT_RESTRICT_STORE=0
@@ -3939,8 +3937,6 @@ SUBROUTINE WRITE_DIAGNOSTICS(T,DT)
            6X,' Min divergence: ',E9.2,' at (',I0,',',I0,',',I0,')')
 133 FORMAT(6X,' Max div. error: ',E9.2,' at (',I0,',',I0,',',I0,')')
 230 FORMAT(6X,' Max VN number : ',E9.2,' at (',I0,',',I0,',',I0,')')
-119 FORMAT(6X,' Total Heat Release Rate:      ',F13.3,' kW')
-120 FORMAT(6X,' Radiation Loss to Boundaries: ',F13.3,' kW')
 141 FORMAT(6X,' No. of Lagrangian Particles:  ',I0)
 121 FORMAT(6X,' No. of CLIP DT restrictions:  ',I0)
 
@@ -9950,15 +9946,15 @@ END SUBROUTINE DUMP_HVAC
 !> \param DT Current time step size (s)
 !> \param NM Mesh number
 !> \details
-!> Q_DOT(1,NM) = \f$ \int \dot{q}''' \, dV \f$
-!> Q_DOT(2,NM) = \f$ \int \dot{q}_{\rm ox}'' \, dS \f$
-!> Q_DOT(3,NM) = \f$ \int \nabla \cdot \mathbf{q}_{\rm r}'' \, dV \f$
-!> Q_DOT(4,NM) = \f$ \int \mathbf{u} \rho h_{\rm s} \cdot \, d\mathbf{S} \f$
-!> Q_DOT(5,NM) = \f$ \int k \nabla T \cdot d\mathbf{S} \f$
-!> Q_DOT(6,NM) = \f$ \int \sum_\alpha h_{{\rm s},\alpha} \rho D_\alpha \nabla Z_\alpha \cdot d\mathbf{S} \f$
-!> Q_DOT(7,NM) = \f$ \int dp/dt \, dV \f$
-!> Q_DOT(8,NM) = \f$ \sum \dot{q}_{\rm p} \f$
-!> Q_DOT(9,NM) = \f$ \int d(\rho h_{\rm s})/dt \, dV \f$
+!> Q_DOT(1) = \f$ \int \dot{q}''' \, dV \f$
+!> Q_DOT(2) = \f$ \int \dot{q}_{\rm ox}'' \, dS \f$
+!> Q_DOT(3) = \f$ \int \nabla \cdot \mathbf{q}_{\rm r}'' \, dV \f$
+!> Q_DOT(4) = \f$ \int \mathbf{u} \rho h_{\rm s} \cdot \, d\mathbf{S} \f$
+!> Q_DOT(5) = \f$ \int k \nabla T \cdot d\mathbf{S} \f$
+!> Q_DOT(6) = \f$ \int \sum_\alpha h_{{\rm s},\alpha} \rho D_\alpha \nabla Z_\alpha \cdot d\mathbf{S} \f$
+!> Q_DOT(7) = \f$ \int dp/dt \, dV \f$
+!> Q_DOT(8) = \f$ \sum \dot{q}_{\rm p} \f$
+!> Q_DOT(9) = \f$ \int d(\rho h_{\rm s})/dt \, dV \f$
 
 SUBROUTINE UPDATE_HRR(DT,NM)
 
@@ -9989,22 +9985,22 @@ SUBROUTINE UPDATE_HRR(DT,NM)
          IF (TWO_D)       VC = VC/DY(J)
          IF (CYLINDRICAL) VC = VC*2._EB*PI
 
-         Q_DOT(1,NM) = Q_DOT(1,NM) + Q(I,J,K)*VC
-         Q_DOT(3,NM) = Q_DOT(3,NM) + QR(I,J,K)*VC
-         Q_DOT(7,NM) = Q_DOT(7,NM) + 0.5_EB*(D_PBAR_DT_S(PRESSURE_ZONE(I,J,K))+D_PBAR_DT(PRESSURE_ZONE(I,J,K)))*VC
+         Q_DOT(1) = Q_DOT(1) + Q(I,J,K)*VC
+         Q_DOT(3) = Q_DOT(3) + QR(I,J,K)*VC
+         Q_DOT(7) = Q_DOT(7) + 0.5_EB*(D_PBAR_DT_S(PRESSURE_ZONE(I,J,K))+D_PBAR_DT(PRESSURE_ZONE(I,J,K)))*VC
          ZZ_GET(1:N_TRACKED_SPECIES) = ZZ(I,J,K,1:N_TRACKED_SPECIES)
          CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMP(I,J,K))
          ENTHALPY_SUM(NM) = ENTHALPY_SUM(NM) + RHO(I,J,K)*H_S*VC
       ENDDO
    ENDDO
 ENDDO
 
-IF (CC_IBM) CALL ADD_Q_DOT_CUTCELLS(NM,Q_DOT(1,NM),Q_DOT(3,NM),Q_DOT(7,NM),ENTHALPY_SUM(NM))
+IF (CC_IBM) CALL ADD_Q_DOT_CUTCELLS(NM,Q_DOT(1),Q_DOT(3),Q_DOT(7),ENTHALPY_SUM(NM))
 
 IF (ICYC>0) THEN
-   Q_DOT(9,NM) = (ENTHALPY_SUM(NM)-ENTHALPY_SUM_OLD)/DT
+   Q_DOT(9) = Q_DOT(9) + (ENTHALPY_SUM(NM)-ENTHALPY_SUM_OLD)/DT
 ELSE
-   Q_DOT(9,NM) = 0._EB
+   Q_DOT(9) = 0._EB
 ENDIF
 
 ! Compute the surface integral of all Del Dot terms
@@ -10053,10 +10049,10 @@ SUBROUTINE UPDATE_HRR(DT,NM)
    AREA_F = B1%AREA
    IF (TWO_D)       AREA_F = AREA_F/DY(BC%JJG)
    IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI
-   Q_DOT(2,NM) = Q_DOT(2,NM) + B1%Q_DOT_O2_PP*AREA_F
-   Q_DOT(4,NM) = Q_DOT(4,NM) - U_N*B1%RHO_F*H_S*AREA_F
-   Q_DOT(5,NM) = Q_DOT(5,NM) - B1%Q_CON_F*AREA_F
-   Q_DOT(6,NM) = Q_DOT(6,NM) - H_S_J_ALPHA*AREA_F
+   Q_DOT(2) = Q_DOT(2) + B1%Q_DOT_O2_PP*AREA_F
+   Q_DOT(4) = Q_DOT(4) - U_N*B1%RHO_F*H_S*AREA_F
+   Q_DOT(5) = Q_DOT(5) - B1%Q_CON_F*AREA_F
+   Q_DOT(6) = Q_DOT(6) - H_S_J_ALPHA*AREA_F
 ENDDO WALL_LOOP
 
 CFACE_LOOP : DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
@@ -10085,9 +10081,9 @@ SUBROUTINE UPDATE_HRR(DT,NM)
    AREA_F = B1%AREA
    IF (TWO_D)       AREA_F = AREA_F/DY(BC%JJG)
    IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI
-   Q_DOT(4,NM) = Q_DOT(4,NM) - U_N*B1%RHO_F*H_S*AREA_F
-   Q_DOT(5,NM) = Q_DOT(5,NM) - B1%Q_CON_F*AREA_F
-   Q_DOT(6,NM) = Q_DOT(6,NM) - H_S_J_ALPHA*AREA_F
+   Q_DOT(4) = Q_DOT(4) - U_N*B1%RHO_F*H_S*AREA_F
+   Q_DOT(5) = Q_DOT(5) - B1%Q_CON_F*AREA_F
+   Q_DOT(6) = Q_DOT(6) - H_S_J_ALPHA*AREA_F
 
 ENDDO CFACE_LOOP
 
@@ -10101,7 +10097,7 @@ SUBROUTINE UPDATE_HRR(DT,NM)
       AREA_F = B1%AREA
       IF (TWO_D)       AREA_F = AREA_F/DY(BC%JJG)
       IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI
-      Q_DOT(2,NM) = Q_DOT(2,NM) + LP%PWT*B1%Q_DOT_O2_PP*AREA_F
+      Q_DOT(2) = Q_DOT(2) + LP%PWT*B1%Q_DOT_O2_PP*AREA_F
    ENDDO PARTICLE_LOOP
 ENDIF
 
@@ -10116,7 +10112,7 @@ SUBROUTINE UPDATE_HRR(DT,NM)
    IF (TWO_D)       AREA_F = AREA_F/DY(BC%JJG)
    IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI
    DO N=1,N_TRACKED_SPECIES
-      M_DOT(N,NM) = M_DOT(N,NM) + B1%M_DOT_G_PP_ADJUST(N)*AREA_F
+      M_DOT(N) = M_DOT(N) + B1%M_DOT_G_PP_ADJUST(N)*AREA_F
    ENDDO
 ENDDO WALL_LOOP2
 
@@ -10129,12 +10125,12 @@ SUBROUTINE UPDATE_HRR(DT,NM)
    IF (TWO_D)       AREA_F = AREA_F/DY(BC%JJG)
    IF (CYLINDRICAL) AREA_F = AREA_F*2._EB*PI
    DO N=1,N_TRACKED_SPECIES
-      M_DOT(N,NM) = M_DOT(N,NM) + B1%M_DOT_G_PP_ADJUST(N)*AREA_F*B1%AREA_ADJUST
+      M_DOT(N) = M_DOT(N) + B1%M_DOT_G_PP_ADJUST(N)*AREA_F*B1%AREA_ADJUST
    ENDDO
 ENDDO CFACE_LOOP_2
 
-Q_DOT_SUM(:,NM) = Q_DOT_SUM(:,NM) + DT*Q_DOT(:,NM)
-M_DOT_SUM(:,NM) = M_DOT_SUM(:,NM) + DT*M_DOT(:,NM)
+Q_DOT_SUM = Q_DOT_SUM + DT*Q_DOT
+M_DOT_SUM = M_DOT_SUM + DT*M_DOT
 
 END SUBROUTINE UPDATE_HRR
 
@@ -10147,18 +10143,13 @@ SUBROUTINE DUMP_HRR(T,DT)
 
 REAL(EB), INTENT(IN) :: T,DT
 REAL(FB) :: STIME
-INTEGER :: NM,I,N_ZONE_TMP
-REAL(EB) :: Q_DOT_TOTAL(N_Q_DOT),M_DOT_TOTAL(N_TRACKED_SPECIES)
+INTEGER :: I,N_ZONE_TMP
 REAL(EB), DIMENSION(:), ALLOCATABLE ::  P_ZONE_P
 
 STIME = REAL(T_BEGIN + (T-T_BEGIN)*TIME_SHRINK_FACTOR,FB)
-Q_DOT_TOTAL = 0._EB
-M_DOT_TOTAL = 0._EB
 
-DO NM=1,NMESHES
-   Q_DOT_TOTAL(:) = Q_DOT_TOTAL(:) + Q_DOT_SUM(:,NM)/MAX(DT,T-T_LAST_DUMP_HRR)
-   M_DOT_TOTAL(:) = M_DOT_TOTAL(:) + M_DOT_SUM(:,NM)/MAX(DT,T-T_LAST_DUMP_HRR)
-ENDDO
+Q_DOT_SUM = Q_DOT_SUM/MAX(DT,T-T_LAST_DUMP_HRR)
+M_DOT_SUM = M_DOT_SUM/MAX(DT,T-T_LAST_DUMP_HRR)
 
 N_ZONE_TMP = 0
 IF (N_ZONE>0) THEN
@@ -10171,11 +10162,11 @@ SUBROUTINE DUMP_HRR(T,DT)
 
 WRITE(TCFORM,'(A,I0,5A)') "(",N_Q_DOT+1+N_TRACKED_SPECIES+N_ZONE_TMP,"(",FMT_R,",','),",FMT_R,")"
 IF (N_ZONE_TMP>0) THEN
-   WRITE(LU_HRR,TCFORM) STIME,0.001_EB*Q_DOT_TOTAL(1:N_Q_DOT),0.001_EB*SUM(Q_DOT_TOTAL(1:N_Q_DOT-1)),&
-                        M_DOT_TOTAL(1:N_TRACKED_SPECIES),(P_ZONE_P(I),I=1,N_ZONE_TMP)
+   WRITE(LU_HRR,TCFORM) STIME,0.001_EB*Q_DOT_SUM(1:N_Q_DOT),0.001_EB*SUM(Q_DOT_SUM(1:N_Q_DOT-1)),&
+                        M_DOT_SUM(1:N_TRACKED_SPECIES),(P_ZONE_P(I),I=1,N_ZONE_TMP)
 ELSE
-   WRITE(LU_HRR,TCFORM) STIME,0.001_EB*Q_DOT_TOTAL(1:N_Q_DOT),0.001_EB*SUM(Q_DOT_TOTAL(1:N_Q_DOT-1)),&
-                        M_DOT_TOTAL(1:N_TRACKED_SPECIES)
+   WRITE(LU_HRR,TCFORM) STIME,0.001_EB*Q_DOT_SUM(1:N_Q_DOT),0.001_EB*SUM(Q_DOT_SUM(1:N_Q_DOT-1)),&
+                        M_DOT_SUM(1:N_TRACKED_SPECIES)
 ENDIF
 
 IF (N_ZONE>0) DEALLOCATE(P_ZONE_P)
@@ -10234,7 +10225,7 @@ SUBROUTINE UPDATE_MASS(DT,NM)
    ENDDO
 ENDDO
 
-MASS_DT(:,NM) = MASS_DT(:,NM) + DT*MASS_INTEGRAL(:)
+MASS_DT = MASS_DT + DT*MASS_INTEGRAL
 
 END SUBROUTINE UPDATE_MASS
 
@@ -10247,21 +10238,14 @@ SUBROUTINE DUMP_MASS(T,DT)
 
 REAL(EB), INTENT(IN) :: T,DT
 REAL(FB) :: STIME
-REAL(EB) :: MASS_TOTAL(0:N_SPECIES+N_TRACKED_SPECIES)
-INTEGER :: NM,N,N_TOTAL_SPECIES
+INTEGER :: N,N_TOTAL_SPECIES
 
 IF (.NOT.MASS_FILE) RETURN
-
 STIME = REAL(T_BEGIN + (T-T_BEGIN)*TIME_SHRINK_FACTOR,FB)
-MASS_TOTAL(:) = 0._EB
-
 N_TOTAL_SPECIES=N_SPECIES+N_TRACKED_SPECIES
-DO NM=1,NMESHES
-   MASS_TOTAL(0:N_TOTAL_SPECIES) = MASS_TOTAL(0:N_TOTAL_SPECIES) + MASS_DT(0:N_TOTAL_SPECIES,NM)/MAX(DT,T-T_LAST_DUMP_MASS)
-ENDDO
 
 WRITE(TCFORM,'(A,I0,5A)') "(",N_TOTAL_SPECIES+1,"(",FMT_R,",','),",FMT_R,")"
-WRITE(LU_MASS,TCFORM) STIME,(MASS_TOTAL(N),N=0,N_TOTAL_SPECIES)
+WRITE(LU_MASS,TCFORM) STIME,(MASS_DT(N)/MAX(DT,T-T_LAST_DUMP_MASS),N=0,N_TOTAL_SPECIES)
 
 END SUBROUTINE DUMP_MASS
 

Source/init.f90

@@ -2832,22 +2832,11 @@ SUBROUTINE INITIALIZE_GLOBAL_VARIABLES
 ALLOCATE(ENTHALPY_SUM(NMESHES),STAT=IZERO)
 CALL ChkMemErr('INIT','ENTHALPY_SUM',IZERO)
 ENTHALPY_SUM = 0._EB
-ALLOCATE(Q_DOT(N_Q_DOT,NMESHES),STAT=IZERO)
-CALL ChkMemErr('INIT','Q_DOT',IZERO)
-Q_DOT = 0._EB
-ALLOCATE(Q_DOT_SUM(N_Q_DOT,NMESHES),STAT=IZERO)
-CALL ChkMemErr('INIT','Q_DOT_SUM',IZERO)
-Q_DOT_SUM = 0._EB
-ALLOCATE(M_DOT(N_TRACKED_SPECIES,NMESHES),STAT=IZERO)
-CALL ChkMemErr('INIT','M_DOT',IZERO)
-M_DOT = 0._EB
-ALLOCATE(M_DOT_SUM(N_TRACKED_SPECIES,NMESHES),STAT=IZERO)
-CALL ChkMemErr('INIT','M_DOT_SUM',IZERO)
-M_DOT_SUM=0._EB
-
-ALLOCATE(MASS_DT(0:N_SPECIES+N_TRACKED_SPECIES,NMESHES),STAT=IZERO)
-CALL ChkMemErr('INIT','MASS_DT',IZERO)
-MASS_DT=0._EB
+ALLOCATE(Q_DOT(N_Q_DOT),STAT=IZERO)     ; CALL ChkMemErr('INIT','Q_DOT',IZERO)     ; Q_DOT = 0._EB
+ALLOCATE(Q_DOT_SUM(N_Q_DOT),STAT=IZERO) ; CALL ChkMemErr('INIT','Q_DOT_SUM',IZERO) ; Q_DOT_SUM = 0._EB
+ALLOCATE(M_DOT(N_TRACKED_SPECIES),STAT=IZERO)     ; CALL ChkMemErr('INIT','M_DOT',IZERO)     ; M_DOT = 0._EB
+ALLOCATE(M_DOT_SUM(N_TRACKED_SPECIES),STAT=IZERO) ; CALL ChkMemErr('INIT','M_DOT_SUM',IZERO) ; M_DOT_SUM=0._EB
+ALLOCATE(MASS_DT(0:N_SPECIES+N_TRACKED_SPECIES),STAT=IZERO) ; CALL ChkMemErr('INIT','MASS_DT',IZERO) ; MASS_DT=0._EB
 
 ALLOCATE(PRESSURE_ERROR_MAX(NMESHES),STAT=IZERO)
 CALL ChkMemErr('INIT','PRESSURE_ERROR_MAX',IZERO)