Merge remote-tracking branch 'firemodels/master' into PMMA

Build/makefile

@@ -59,6 +59,8 @@ FFLAGSMKL_CUSTOM =
 LFLAGSMKL_CUSTOM =
 FFLAGS_SUNDIALS =
 LFLAGS_SUNDIALS =
+FFLAGS_HYPRE =
+LFLAGS_HYPRE =
 ifdef MKLROOT # This assumes the MKL library is available.
 ifeq ($(shell echo "check_quotes"),"check_quotes")
 # windows
@@ -109,6 +111,11 @@ ifdef SUNDIALS_HOME # This assumes the SUNDIALS library is available.
    LFLAGS_SUNDIALS_WIN = ${SUNDIALS_HOME}/lib/sundials_fcvode_mod.lib ${SUNDIALS_HOME}/lib/sundials_fnvecserial_mod.lib ${SUNDIALS_HOME}/lib/sundials_cvode.lib  /link /NODEFAULTLIB:MSVCRTD /NODEFAULTLIB:libcmtd.lib
 endif
 
+ifdef HYPRE_HOME # This assumes the HYPRE library is available.
+   FFLAGS_HYPRE = -DWITH_HYPRE -I${HYPRE_HOME}/include
+   LFLAGS_HYPRE = -L${HYPRE_HOME}/lib -lHYPRE -lm
+endif
+
 obj_mpi = prec.o cons.o chem.o prop.o devc.o type.o data.o mesh.o func.o gsmv.o smvv.o rcal.o turb.o soot.o \
           pois.o geom.o ccib.o radi.o part.o vege.o ctrl.o hvac.o mass.o \
           wall.o fire.o velo.o pres.o init.o dump.o read.o divg.o main.o
@@ -337,6 +344,14 @@ ompi_gnu_linux_db : obj = fds_ompi_gnu_linux_db
 ompi_gnu_linux_db : setup $(obj_mpi)
 	$(FCOMPL) $(FFLAGS) $(FOPENMPFLAGS) -o $(obj) $(obj_mpi) $(LFLAGSMKL)
 
+ompi_gnu_linux_dv : FFLAGS = -m64 -O1 -fbacktrace -std=f2018 -frecursive -ffpe-summary=none -fall-intrinsics $(GITINFOGNU) $(FFLAGSMKL_GNU_OPENMPI) $(GFORTRAN_OPTIONS)
+ompi_gnu_linux_dv : LFLAGSMKL = $(LFLAGSMKL_GNU_OPENMPI)
+ompi_gnu_linux_dv : FCOMPL = mpifort
+ompi_gnu_linux_dv : FOPENMPFLAGS = -fopenmp
+ompi_gnu_linux_dv : obj = fds_ompi_gnu_linux_dv
+ompi_gnu_linux_dv : setup $(obj_mpi)
+	$(FCOMPL) $(FFLAGS) $(FOPENMPFLAGS) -o $(obj) $(obj_mpi) $(LFLAGSMKL)
+
 ompi_gnu_osx : FFLAGS = -m64 -O2 -std=f2018 -frecursive -ffpe-summary=none -fall-intrinsics $(GITINFOGNU) $(FFLAGSMKL_GNU_CUSTOM) $(GFORTRAN_OPTIONS) $(FFLAGS_SUNDIALS)
 ompi_gnu_osx : LFLAGSMKL = $(LFLAGSMKL_GNU_CUSTOM) $(CLT_VERSION) $(LFLAGS_SUNDIALS_OSX)
 ompi_gnu_osx : FCOMPL = mpifort

Build/ompi_gnu_linux_dv/make_fds.sh

@@ -0,0 +1,6 @@
+#!/bin/bash
+dir=`pwd`
+target=${dir##*/}
+
+echo Building $target
+make -j4 VPATH="../../Source" -f ../makefile $target

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -3066,11 +3066,12 @@ \subsection{Liquid Fuels}
 \subsubsection{Evaporation of a Pure Liquid}
 \label{methanol_evaporation}
 
-An example of liquid evaporation is given by the sample case found in the {\ct Pyrolysis} folder called {\ct methanol\_evaporation.fds}. A 1~m by 1~m pan filled with methanol at $T_\infty=20$~$^\circ$C is exposed to a uniform heat flux, $\dot{q}''=20$~\unit{kW/m^2}. The boiling temperature of methanol is $T_{\rm b}=64.65$~$^\circ$C, its specific heat, $c=2.48$~kJ/(kg$\cdot$K), and heat of vaporization, $h_{\rm v}=1099$~kJ/kg. The evaporation rate of a burning liquid in steady state is approximately
+An example of liquid evaporation is given by the sample case found in the {\ct Pyrolysis} folder called {\ct methanol\_evaporation.fds}. A 1~m by 1~m pan filled with methanol at $T_\infty=20$~$^\circ$C is exposed to a uniform heat flux, $\dot{q}''=20$~\unit{kW/m^2}. The boiling temperature of methanol is $T_{\rm b}=64.65$~$^\circ$C, its specific heat, $c=2.48$~kJ/(kg$\cdot$K), and heat of vaporization, $h_{\rm v}=1099$~kJ/kg. At steady state, the heat balance at the pool surface is
 \be
-   \dot{m}'' \approx \frac{\dot{q}''}{h_{\rm g}}  \quad ; \quad  h_{\rm g} = c (T_{\rm b}-T_\infty) + h_{\rm v}
+   \dot{q}''_{total} - \dot{q}''_{c}= \dot{m}'' h_{\rm v}(T_s)
 \ee
-In this example, the methanol evaporates in an oxygen-depleted atmosphere and no burning occurs. The left hand plot in Fig.~\ref{methanol_evaporation_plot} displays the computed evaporation rate, $\dot{m}''$, versus the ideal, $\dot{q}''/h_{\rm g}$. The former approaches the latter as all of the absorbed energy is used to evaporate the liquid. The right hand plot shows the computed liquid surface temperature versus the liquid boiling temperature.
+
+where $ \dot{q}''_{c}$ is the heat being conducted away from the surface. If $ \dot{q}''_{c}$ is made zero, which can be done by using {\ct BACKING='INSULATED'} and a high thermal conductivity, then the pool surface temperature,$T_s$, will approach the boiling temperature,$T_b$. In this example, the methanol evaporates in an oxygen-depleted atmosphere and no burning occurs. The left hand plot in Fig.~\ref{methanol_evaporation_plot} displays the computed evaporation rate, $\dot{m}''$, versus the ideal, $\dot{q}''_{}total}/h_{\rm v}(T_b)$. The former approaches the latter as all of the absorbed energy is used to evaporate the liquid. The right hand plot shows the computed liquid surface temperature versus the liquid boiling temperature.
 \begin{figure}[!ht]
 \includegraphics[width=3.2in]{SCRIPT_FIGURES/methanol_evaporation_mdot}
 \includegraphics[width=3.2in]{SCRIPT_FIGURES/methanol_evaporation_temp}

Source/cons.f90

@@ -211,6 +211,7 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: CHECK_VN=.TRUE.                  !< Check the Von Neumann number
 LOGICAL :: CHECK_FO=.FALSE.                 !< Check the solid phase Fourier number
 LOGICAL :: SOLID_PARTICLES=.FALSE.          !< Indicates the existence of solid particles
+LOGICAL :: ORIENTED_PARTICLES=.FALSE.       !< Indicates the existence of particles with a specified orientation
 LOGICAL :: HVAC=.FALSE.                     !< Perform an HVAC calculation
 LOGICAL :: BAROCLINIC=.TRUE.                !< Include the baroclinic terms in the momentum equation
 LOGICAL :: GRAVITATIONAL_DEPOSITION=.TRUE.  !< Allow aerosol gravitational deposition
@@ -248,6 +249,7 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: DUCT_HT_INSERTED=.FALSE.
 LOGICAL :: HVAC_QFAN=.FALSE.
 LOGICAL :: USE_ATMOSPHERIC_INTERPOLATION=.FALSE.
+LOGICAL :: NEAR_WALL_PARTICLE_INTERPOLATION=.FALSE.
 LOGICAL :: POSITIVE_ERROR_TEST=.FALSE.
 LOGICAL :: OBST_SHAPE_AREA_ADJUST=.FALSE.
 LOGICAL :: STORE_SPECIES_FLUX=.FALSE.
@@ -725,8 +727,8 @@ MODULE GLOBAL_CONSTANTS
 
 REAL(EB), POINTER, DIMENSION(:,:) :: ORIENTATION_VECTOR       !< Global array of orientation vectors
 INTEGER, ALLOCATABLE, DIMENSION(:) :: NEAREST_RADIATION_ANGLE !< Index of the rad angle most opposite the given ORIENTATION_VECTOR
-REAL(EB), POINTER, DIMENSION(:) :: ORIENTATION_VIEW_ANGLE     !< View angle of the given ORIENTATION_VECTOR
-REAL(EB), ALLOCATABLE, DIMENSION(:) :: VIEW_ANGLE_AREA        !< View angle area ORIENTATION_VECTOR
+REAL(EB), POINTER, DIMENSION(:) :: COS_HALF_VIEW_ANGLE     !< View angle of the given ORIENTATION_VECTOR
+REAL(EB), ALLOCATABLE, DIMENSION(:) :: VIEW_ANGLE_FACTOR        !< View angle area ORIENTATION_VECTOR
 INTEGER :: N_ORIENTATION_VECTOR                               !< Number of ORIENTATION_VECTORs
 
 INTEGER :: TGA_MESH_INDEX=HUGE(INTEGER_ONE)  !< Mesh for the special TGA calculation

Source/dump.f90

@@ -3011,7 +3011,10 @@ SUBROUTINE INITIALIZE_DIAGNOSTIC_FILE(DT)
          WRITE(LU_OUTPUT,'(A,A,A,F8.2)')'        ',SPECIES_MIXTURE(NS)%ID,': ',ML%NU_GAS(NS,1)
       ENDDO
       WRITE(LU_OUTPUT,'(A,F8.2)') '        Boiling temperature (C): ',ML%TMP_BOIL-TMPM
-      WRITE(LU_OUTPUT,'(A,ES10.3)')'        H_R (kJ/kg)            : ',ML%H_R(1,NINT(TMPA))/1000._EB
+      ITMP = NINT(TMPA)
+      WRITE(LU_OUTPUT,'(A,I4,A,ES10.3)') '        H_R (kJ/kg) TMPA,    ',ITMP,' K: ',ML%H_R(1,ITMP)/1000._EB
+      ITMP = NINT(ML%TMP_REF(1))
+      WRITE(LU_OUTPUT,'(A,I4,A,ES10.3)') '        H_R (kJ/kg) TMP_REF, ',ITMP,' K: ',ML%H_R(1,ITMP)/1000._EB
    ENDIF
 
 ENDDO MATL_LOOP

Source/func.f90

@@ -1477,7 +1477,6 @@ SUBROUTINE PACK_PARTICLE(NM,OS,LP,LPC_INDEX,RC,IC,LC,UNPACK_IT,COUNT_ONLY)
 
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%TAG,UNPACK_IT)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%CLASS_INDEX,UNPACK_IT)
-IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%INITIALIZATION_INDEX,UNPACK_IT)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%ORIENTATION_INDEX,UNPACK_IT)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%WALL_INDEX,UNPACK_IT)
 IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC),LP%DUCT_INDEX,UNPACK_IT)
@@ -3497,100 +3496,6 @@ SUBROUTINE UPDATE_HISTOGRAM(NBINS,LIMITS,COUNTS,VAL,WEIGHT)
 COUNTS(IND)=COUNTS(IND)+WEIGHT
 END SUBROUTINE UPDATE_HISTOGRAM
 
-
-!> \brief Linearly interpolate the a mesh quantity onto a point
-!> \param X The interpolated value of the 3D array
-!> \param A The 3D array of values
-!> \param I The lower x index of the array
-!> \param J The lower y index of the array
-!> \param K The lower z index of the array
-!> \param P Fraction of the distance from the lower to upper x coordinate
-!> \param R Fraction of the distance from the lower to upper y coordinate
-!> \param S Fraction of the distance from the lower to upper z coordinate
-
-SUBROUTINE MESH_TO_PARTICLE(X,A,I,J,K,P,R,S)
-
-REAL(EB), INTENT(IN), DIMENSION(0:,0:,0:) :: A
-INTEGER, INTENT(IN) :: I,J,K
-REAL(EB), INTENT(IN) :: P,R,S
-REAL(EB), INTENT(OUT) :: X
-REAL(EB) :: PP,RR,SS
-
-PP = 1._EB-P
-RR = 1._EB-R
-SS = 1._EB-S
-X  = ((PP*A(I,J,K)  +P*A(I+1,J,K)  )*RR+(PP*A(I,J+1,K)  +P*A(I+1,J+1,K)  )*R)*SS + &
-     ((PP*A(I,J,K+1)+P*A(I+1,J,K+1))*RR+(PP*A(I,J+1,K+1)+P*A(I+1,J+1,K+1))*R)*S
-
-END SUBROUTINE MESH_TO_PARTICLE
-
-
-!> \brief Linearly interpolate the value at a point onto the mesh
-!> \param X The interpolated value of the 3D array
-!> \param A The 3D array of values
-!> \param I The lower x index of the array
-!> \param J The lower y index of the array
-!> \param K The lower z index of the array
-!> \param P Fraction of the distance from the lower to upper x coordinate
-!> \param R Fraction of the distance from the lower to upper y coordinate
-!> \param S Fraction of the distance from the lower to upper z coordinate
-
-SUBROUTINE PARTICLE_TO_MESH(X,A,I,J,K,P,R,S)
-
-REAL(EB), INTENT(INOUT), DIMENSION(0:,0:,0:) :: A
-INTEGER, INTENT(IN) :: I,J,K
-REAL(EB), INTENT(IN) :: P,R,S
-REAL(EB), INTENT(IN) :: X
-REAL(EB) :: PP,RR,SS
-
-PP = 1._EB-P
-RR = 1._EB-R
-SS = 1._EB-S
-A(I  ,J  ,K  ) = A(I  ,J  ,K  ) - X*PP*RR*SS
-A(I+1,J  ,K  ) = A(I+1,J  ,K  ) - X*P *RR*SS
-A(I  ,J+1,K  ) = A(I  ,J+1,K  ) - X*PP*R *SS
-A(I+1,J+1,K  ) = A(I+1,J+1,K  ) - X*P *R *SS
-A(I  ,J  ,K+1) = A(I  ,J  ,K+1) - X*PP*RR*S
-A(I+1,J  ,K+1) = A(I+1,J  ,K+1) - X*P *RR*S
-A(I  ,J+1,K+1) = A(I  ,J+1,K+1) - X*PP*R *S
-A(I+1,J+1,K+1) = A(I+1,J+1,K+1) - X*P *R *S
-
-END SUBROUTINE PARTICLE_TO_MESH
-
-
-!> \brief Trilinear interpolation https://paulbourke.net/miscellaneous/interpolation/
-!> \param V The interpolated value of the 3D array
-!> \param A The 3D array of box corner values
-!> \param X Fractional distance from the lower to upper x coordinate
-!> \param Y Fractional distance from the lower to upper y coordinate
-!> \param Z Fractional distance from the lower to upper z coordinate
-
-SUBROUTINE TRILIN_INTERP(V,A,X,Y,Z)
-
-REAL(EB), INTENT(IN), DIMENSION(0:1,0:1,0:1) :: A
-REAL(EB), INTENT(IN) :: X,Y,Z
-REAL(EB), INTENT(OUT) :: V
-REAL(EB), DIMENSION(0:1,0:1,0:1) :: WGT
-REAL(EB) :: XX,YY,ZZ
-
-XX = 1._EB-X
-YY = 1._EB-Y
-ZZ = 1._EB-Z
-
-WGT(0,0,0) = XX * YY * ZZ
-WGT(1,0,0) = X  * YY * ZZ
-WGT(0,1,0) = XX * Y  * ZZ
-WGT(0,0,1) = XX * YY * Z
-WGT(1,0,1) = X  * YY * Z
-WGT(0,1,1) = XX * Y  * Z
-WGT(1,1,0) = X  * Y  * ZZ
-WGT(1,1,1) = X  * Y  * Z
-
-V = SUM(A*WGT)
-
-END SUBROUTINE TRILIN_INTERP
-
-
 !> \brief Calculate the value of polynomial function.
 !> \param N Number of coefficients in the polynomial
 !> \param TEMP The independent variable

Source/main.f90

@@ -52,10 +52,9 @@ PROGRAM FDS
 INTEGER  :: LO10,NM,IZERO,ANG_INC_COUNTER
 REAL(EB) :: T,DT,TNOW
 REAL :: CPUTIME
-REAL(EB), ALLOCATABLE, DIMENSION(:) ::  TC_GLB,TC_LOC,DT_NEW,TI_LOC,TI_GLB, &
-                                        DSUM_ALL,PSUM_ALL,USUM_ALL,DSUM_ALL_LOCAL,PSUM_ALL_LOCAL,USUM_ALL_LOCAL
+REAL(EB), ALLOCATABLE, DIMENSION(:) ::  TC_GLB,TC_LOC,DT_NEW,TI_LOC,TI_GLB,DSUM_ALL,PSUM_ALL,USUM_ALL
 REAL(EB), ALLOCATABLE, DIMENSION(:,:) ::  TC2_GLB,TC2_LOC
-LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES_GLOBAL,CONNECTED_ZONES_LOCAL
+LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES_ALL
 LOGICAL, ALLOCATABLE, DIMENSION(:) ::  STATE_GLB,STATE_LOC
 INTEGER :: ITER
 TYPE (MESH_TYPE), POINTER :: M,M4
@@ -519,7 +518,7 @@ PROGRAM FDS
 
 ! If there are zones and HVAC pass PSUM
 
-IF (HVAC_SOLVE .AND. N_ZONE>0) CALL EXCHANGE_DIVERGENCE_INFO
+IF (HVAC_SOLVE .AND. N_ZONE>0 .AND. .NOT.SOLID_PHASE_ONLY) CALL EXCHANGE_DIVERGENCE_INFO
 
 ! Make an initial dump of global output quantities
 
@@ -687,7 +686,7 @@ PROGRAM FDS
 
       ! If there are pressure ZONEs, exchange integrated quantities mesh to mesh for use in the divergence calculation
 
-      IF (N_ZONE>0) CALL EXCHANGE_DIVERGENCE_INFO
+      IF (N_ZONE>0 .AND. .NOT.SOLID_PHASE_ONLY) CALL EXCHANGE_DIVERGENCE_INFO
 
       ! Update global pressure matrices after zone connections
 
@@ -886,7 +885,7 @@ PROGRAM FDS
 
    ! Exchange global pressure zone information
 
-   IF (N_ZONE>0) CALL EXCHANGE_DIVERGENCE_INFO
+   IF (N_ZONE>0 .AND. .NOT.SOLID_PHASE_ONLY) CALL EXCHANGE_DIVERGENCE_INFO
 
    ! Update global pressure matrices after zone connections
 
@@ -1390,11 +1389,7 @@ SUBROUTINE MPI_INITIALIZATION_CHORES(TASK_NUMBER)
          ALLOCATE(DSUM_ALL(N_ZONE),STAT=IZERO)
          ALLOCATE(PSUM_ALL(N_ZONE),STAT=IZERO)
          ALLOCATE(USUM_ALL(N_ZONE),STAT=IZERO)
-         ALLOCATE(CONNECTED_ZONES_GLOBAL(0:N_ZONE,0:N_ZONE),STAT=IZERO)
-         ALLOCATE(DSUM_ALL_LOCAL(N_ZONE),STAT=IZERO)
-         ALLOCATE(PSUM_ALL_LOCAL(N_ZONE),STAT=IZERO)
-         ALLOCATE(USUM_ALL_LOCAL(N_ZONE),STAT=IZERO)
-         ALLOCATE(CONNECTED_ZONES_LOCAL(0:N_ZONE,0:N_ZONE),STAT=IZERO)
+         ALLOCATE(CONNECTED_ZONES_ALL(0:N_ZONE,0:N_ZONE),STAT=IZERO)
       ENDIF
 
       ALLOCATE(CONNECTED_ZONES(0:N_ZONE,0:N_ZONE,NMESHES),STAT=IZERO)
@@ -1818,40 +1813,35 @@ SUBROUTINE EXCHANGE_DIVERGENCE_INFO
 
 TNOW = CURRENT_TIME()
 
-CONNECTED_ZONES_LOCAL = .FALSE.
+CONNECTED_ZONES_ALL = .FALSE.
 
 DO IPZ=1,N_ZONE
-   DSUM_ALL_LOCAL(IPZ) = 0._EB
-   PSUM_ALL_LOCAL(IPZ) = 0._EB
-   USUM_ALL_LOCAL(IPZ) = 0._EB
+   DSUM_ALL(IPZ) = 0._EB
+   PSUM_ALL(IPZ) = 0._EB
+   USUM_ALL(IPZ) = 0._EB
    DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
-      DSUM_ALL_LOCAL(IPZ) = DSUM_ALL_LOCAL(IPZ) + DSUM(IPZ,NM)
-      PSUM_ALL_LOCAL(IPZ) = PSUM_ALL_LOCAL(IPZ) + PSUM(IPZ,NM)
-      USUM_ALL_LOCAL(IPZ) = USUM_ALL_LOCAL(IPZ) + USUM(IPZ,NM)
+      DSUM_ALL(IPZ) = DSUM_ALL(IPZ) + DSUM(IPZ,NM)
+      PSUM_ALL(IPZ) = PSUM_ALL(IPZ) + PSUM(IPZ,NM)
+      USUM_ALL(IPZ) = USUM_ALL(IPZ) + USUM(IPZ,NM)
       DO IOPZ=0,N_ZONE
-         IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_LOCAL(IPZ,IOPZ) = .TRUE.
+         IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_ALL(IPZ,IOPZ) = .TRUE.
       ENDDO
    ENDDO
 ENDDO
 
 IF (N_MPI_PROCESSES>1) THEN
-   CALL MPI_ALLREDUCE(DSUM_ALL_LOCAL(1),DSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
-   CALL MPI_ALLREDUCE(PSUM_ALL_LOCAL(1),PSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
-   CALL MPI_ALLREDUCE(USUM_ALL_LOCAL(1),USUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
-   CALL MPI_ALLREDUCE(CONNECTED_ZONES_LOCAL(0,0),CONNECTED_ZONES_GLOBAL(0,0),(N_ZONE+1)**2,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR)
-ELSE
-   DSUM_ALL = DSUM_ALL_LOCAL
-   PSUM_ALL = PSUM_ALL_LOCAL
-   USUM_ALL = USUM_ALL_LOCAL
-   CONNECTED_ZONES_GLOBAL = CONNECTED_ZONES_LOCAL
+   CALL MPI_ALLREDUCE(MPI_IN_PLACE,DSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
+   CALL MPI_ALLREDUCE(MPI_IN_PLACE,PSUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
+   CALL MPI_ALLREDUCE(MPI_IN_PLACE,USUM_ALL(1),N_ZONE,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
+   CALL MPI_ALLREDUCE(MPI_IN_PLACE,CONNECTED_ZONES_ALL(0,0),(N_ZONE+1)**2,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR)
 ENDIF
 
 DO IPZ=1,N_ZONE
    DO NM=1,NMESHES
       DSUM(IPZ,NM) = DSUM_ALL(IPZ)
       PSUM(IPZ,NM) = PSUM_ALL(IPZ)
       USUM(IPZ,NM) = USUM_ALL(IPZ)
-      CONNECTED_ZONES(IPZ,:,NM) = CONNECTED_ZONES_GLOBAL(IPZ,:)
+      CONNECTED_ZONES(IPZ,:,NM) = CONNECTED_ZONES_ALL(IPZ,:)
    ENDDO
 ENDDO