main_mpi.f90 - fds-smv - Project Hosting on Google Code

‹r5359

r5421

Source path: svn/ trunk/ FDS/ trunk/ FDS_Source/ main_mpi.f90

PROGRAM FDS  

! Fire Dynamics Simulator, Main Program, Multiple CPU version.

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
USE TRAN
USE DUMP
USE READ_INPUT
USE INIT
USE DIVG
USE PRES
USE MASS
USE PART
USE VEGE
USE VELO
USE RAD
USE MEMORY_FUNCTIONS
USE COMP_FUNCTIONS, ONLY : SECOND, WALL_CLOCK_TIME, SHUTDOWN
USE MATH_FUNCTIONS, ONLY : GAUSSJ
USE DEVICE_VARIABLES
USE WALL_ROUTINES
USE FIRE
USE RADCONS, ONLY: ANGLE_INCREMENT,TIME_STEP_INCREMENT
USE CONTROL_FUNCTIONS
USE EVAC
USE TURBULENCE, ONLY: ANALYTICAL_SOLUTION,SANDIA_DAT,INIT_TURB_ARRAYS,SPECTRAL_OUTPUT
USE SOOT
USE OPENMP

IMPLICIT NONE

! Miscellaneous declarations

CHARACTER(255), PARAMETER :: mainmpiid='$Id$'
CHARACTER(255), PARAMETER :: mainmpirev='$Revision$'
CHARACTER(255), PARAMETER :: mainmpidate='$Date$'
LOGICAL  :: EX,EX_GLOBAL,DIAGNOSTICS,EXCHANGE_RADIATION,EXCHANGE_EVACUATION=.FALSE.,FIRST_PASS
INTEGER  :: LO10,NM,IZERO,REVISION_NUMBER,IOS
CHARACTER(255) :: REVISION_DATE
REAL(EB) :: T_MAX,T_MIN
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: AINV
REAL(EB), ALLOCATABLE, DIMENSION(:) ::  T,TC_GLB,TC_LOC,DT_SYNC,DT_NEXT_SYNC,TI_LOC,TI_GLB, &
                                        DSUM_ALL,PSUM_ALL,USUM_ALL,DSUM_ALL_LOCAL,PSUM_ALL_LOCAL,USUM_ALL_LOCAL
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: ASUM_ALL,ASUM_ALL_LOCAL,FDS_LEAK_AREA_LOCAL,FDS_LEAK_AREA_GLOBAL
LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: CONNECTED_ZONES_GLOBAL,CONNECTED_ZONES_LOCAL
INTEGER, ALLOCATABLE, DIMENSION(:) ::  MESH_STOP_STATUS
LOGICAL, ALLOCATABLE, DIMENSION(:) ::  ACTIVE_MESH,STATE_GLB,STATE_LOC
INTEGER :: NOM,IWW,IW,STOP_STATUS=0
INTEGER, PARAMETER :: N_DROP_ADOPT_MAX=10000
TYPE (MESH_TYPE), POINTER :: M,M4
TYPE (OMESH_TYPE), POINTER :: M2,M3
 
! MPI stuff

! INCLUDE '/usr/local/include/mpif.h'  ! Uncomment this line only if the compiler needs help finding mpif.h
INCLUDE 'mpif.h'
INTEGER :: N,I,IERR,STATUS(MPI_STATUS_SIZE)
INTEGER :: BUFFER_SIZE,TAG,PNAMELEN,DISP
INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: TAGS
INTEGER, ALLOCATABLE, DIMENSION(:) :: REQ,COUNTS,DISPLS,COUNTS2D,DISPLS2D,COUNTS_TIMERS,DISPLS_TIMERS, &
                                      COUNTS_MASS,DISPLS_MASS,DISPLS_PC,COUNTS_PC,DISPLS2D_PC,COUNTS2D_PC
INTEGER, ALLOCATABLE, DIMENSION(:) :: COUNTS_TREE,DISPLS_TREE
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ARRAY_OF_STATUSES
INTEGER :: N_REQ
CHARACTER(MPI_MAX_PROCESSOR_NAME) PNAME
REAL(EB) :: MPI_WALL_TIME,MPI_WALL_TIME_START,RHO_UPPER_LOCAL,RHO_LOWER_LOCAL
 
! Initialize MPI (First executable lines of code)
 
CALL MPI_INIT(IERR)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NUMPROCS, IERR)
CALL MPI_GET_PROCESSOR_NAME(PNAME, PNAMELEN, IERR)
MPI_WALL_TIME_START = MPI_WTIME()
 
WRITE(LU_ERR,'(A,I3,A,I3,A,A)') 'Process ',MYID,' of ',NUMPROCS-1,' is running on ',PNAME(1:PNAMELEN)
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Check for OpenMP

CALL OPENMP_CHECK

IF (USE_OPENMP) THEN
   WRITE(LU_ERR,'(A,A,I3,A)') PNAME(1:PNAMELEN),' has ',OPENMP_AVAILABLE_THREADS,' threads available'
ENDIF
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)

! Start wall clock timing

WALL_CLOCK_START = WALL_CLOCK_TIME()
 
! Assign a compilation date (All Nodes)

WRITE(VERSION_STRING,'(A)') '5.4.3'
PARALLEL = .TRUE.

IF (INDEX(mainmpirev,':',BACK=.TRUE.)>0) THEN
   WRITE(REVISION_DATE,'(A)',IOSTAT=IOS,ERR=5) mainmpirev(INDEX(mainmpirev,':')+1:LEN_TRIM(mainmpirev)-2)
   5 REVISION_NUMBER = 0
   IF (IOS==0) READ(REVISION_DATE,'(I5)') REVISION_NUMBER
   WRITE(REVISION_DATE,'(A)') mainmpidate
   CALL GET_REVISION_NUMBER(REVISION_NUMBER,REVISION_DATE)
   SVN_REVISION_NUMBER = REVISION_NUMBER
   WRITE(COMPILE_DATE,'(A)',IOSTAT=IOS,ERR=10) REVISION_DATE(INDEX(REVISION_DATE,'(')+1:INDEX(REVISION_DATE,')')-1)
   10 IF (IOS>0) COMPILE_DATE = 'null'
ENDIF

! Set some constants

CALL SET_OFTEN_USED

! Read input from CHID.data file (All Nodes)
 
CALL READ_DATA
 
! Make sure that the appropriate number of processes (NUMPROCS) have been assigned

IF (RESTRICT_MESH_ASSIGNMENT .AND. NMESHES/=NUMPROCS .AND. .NOT.ANY(EVACUATION_GRID)) THEN
   IF (MYID==0) WRITE(LU_ERR,'(A)') 'ERROR: The number of meshes is not equal to the number of MPI processes'
   STOP
ENDIF

! Allocate save array for PRESSURE_CORRECTION

IF (PRESSURE_CORRECTION) THEN
   ALLOCATE(AINV(NCGC,NCGC),STAT=IZERO)
   CALL ChkMemErr('MAIN','AINV',IZERO)
   AINV = 0._EB
ENDIF

! Determine if radiation exchanges should be done by default

IF (RADIATION) THEN
   EXCHANGE_RADIATION = .TRUE.
ELSE
   EXCHANGE_RADIATION = .FALSE.
ENDIF

! Set up send and receive buffer counts and displacements

ALLOCATE(COUNTS(0:NUMPROCS-1))
ALLOCATE(COUNTS_PC(0:NUMPROCS-1))
ALLOCATE(COUNTS2D(0:NUMPROCS-1))
ALLOCATE(COUNTS2D_PC(0:NUMPROCS-1))
ALLOCATE(COUNTS_TIMERS(0:NUMPROCS-1))
ALLOCATE(COUNTS_MASS(0:NUMPROCS-1))
ALLOCATE(DISPLS(0:NUMPROCS-1))
ALLOCATE(DISPLS2D(0:NUMPROCS-1))
ALLOCATE(DISPLS_MASS(0:NUMPROCS-1))
ALLOCATE(DISPLS_TIMERS(0:NUMPROCS-1))
ALLOCATE(DISPLS_PC(0:NUMPROCS-1))
ALLOCATE(DISPLS2D_PC(0:NUMPROCS-1))

ALLOCATE(COUNTS_TREE(0:NUMPROCS-1))
ALLOCATE(DISPLS_TREE(0:NUMPROCS-1))

COUNTS    = 0
COUNTS_PC = 0
DO N=0,NUMPROCS-1
   DO NM=1,NMESHES
      IF (PROCESS(NM)==N) COUNTS(N)    = COUNTS(N)    + 1
      IF (PROCESS(NM)==N) COUNTS_PC(N) = COUNTS_PC(N) + MESHES(NM)%IBAR2*MESHES(NM)%JBAR2*MESHES(NM)%KBAR2
   ENDDO
ENDDO
DISPLS(0)    = 0
DISPLS_PC(0) = 0
DO N=1,NUMPROCS-1
   DISPLS(N)    = COUNTS(N-1)    + DISPLS(N-1)
   DISPLS_PC(N) = COUNTS_PC(N-1) + DISPLS_PC(N-1)
ENDDO
COUNTS2D      = COUNTS*NMESHES
DISPLS2D      = DISPLS*NMESHES
COUNTS2D_PC   = COUNTS_PC*NCGC
DISPLS2D_PC   = DISPLS_PC*NCGC
COUNTS_TIMERS = COUNTS*N_TIMERS_DIM
DISPLS_TIMERS = DISPLS*N_TIMERS_DIM
COUNTS_MASS   = COUNTS*(MAX_SPECIES+1)
DISPLS_MASS   = DISPLS*(MAX_SPECIES+1)
 
! Open and write to Smokeview and status file (Master Node Only)
 
CALL ASSIGN_FILE_NAMES

DO N=0,NUMPROCS-1
   IF (MYID==N) CALL WRITE_SMOKEVIEW_FILE
   IF (N==NUMPROCS-1) EXIT
   CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
IF (MYID==0) THEN
   OPEN(LU_SMV,FILE=FN_SMV,FORM='FORMATTED', STATUS='OLD',POSITION='APPEND')
   CALL WRITE_STATUS_FILES
ENDIF

! Stop all the processes if this is just a set-up run
 
IF (SET_UP) CALL SHUTDOWN('Stop FDS, Set-up only')
 
! Set up Time arrays (All Nodes)
 
ALLOCATE(ACTIVE_MESH(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','ACTIVE_MESH',IZERO)
ACTIVE_MESH = .TRUE.
ALLOCATE(T(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','T',IZERO)
ALLOCATE(DT_SYNC(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','DT_SYNC',IZERO)
ALLOCATE(DT_NEXT_SYNC(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','DT_NEXT_SYNC',IZERO)
ALLOCATE(MESH_STOP_STATUS(NMESHES),STAT=IZERO)
CALL ChkMemErr('MAIN','MESH_STOP_STATUS',IZERO)


! Set up dummy arrays to hold various arrays that must be exchanged among meshes

ALLOCATE(TI_LOC(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','TI_LOC',IZERO) 
ALLOCATE(TI_GLB(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','TI_GLB',IZERO) 
ALLOCATE(STATE_GLB(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','STATE_GLB',IZERO) 
ALLOCATE(STATE_LOC(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','STATE_LOC',IZERO) 
ALLOCATE(TC_GLB(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','TC_GLB',IZERO)
ALLOCATE(TC_LOC(N_DEVC),STAT=IZERO)
CALL ChkMemErr('MAIN','TC_LOC',IZERO)

! Allocate a few arrays needed to exchange divergence and pressure info among meshes

IF (N_ZONE > 0) THEN
   ALLOCATE(DSUM_ALL(N_ZONE),STAT=IZERO)
   ALLOCATE(PSUM_ALL(N_ZONE),STAT=IZERO)
   ALLOCATE(USUM_ALL(N_ZONE),STAT=IZERO)
   ALLOCATE(ASUM_ALL(N_ZONE,0:N_ZONE),STAT=IZERO)
   ALLOCATE(FDS_LEAK_AREA_LOCAL(0:N_ZONE,0: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(ASUM_ALL_LOCAL(N_ZONE,0:N_ZONE),STAT=IZERO)
   ALLOCATE(FDS_LEAK_AREA_GLOBAL(0:N_ZONE,0:N_ZONE),STAT=IZERO)
   ALLOCATE(CONNECTED_ZONES_LOCAL(0:N_ZONE,0:N_ZONE),STAT=IZERO)
ENDIF

! Start the clock

T     = T_BEGIN
MESH_STOP_STATUS = NO_STOP
 
! Create unique tags for all mesh exchanges
 
ALLOCATE(REQ(NMESHES*NMESHES*20)) 
REQ = MPI_REQUEST_NULL
ALLOCATE(ARRAY_OF_STATUSES(MPI_STATUS_SIZE,NMESHES*NMESHES*20))
ALLOCATE(TAGS(NMESHES,NMESHES,0:6))
TAG = 0
DO NM=1,NMESHES
   DO NOM=NM,NMESHES
      TAG = TAG+1
      TAGS(NM,NOM,0) = TAG
      TAGS(NOM,NM,0) = TAG
   ENDDO
ENDDO
TAGS(:,:,1) = TAGS(:,:,0) + 1000
TAGS(:,:,2) = TAGS(:,:,0) + 2000
TAGS(:,:,3) = TAGS(:,:,0) + 3000
TAGS(:,:,4) = TAGS(:,:,0) + 4000
TAGS(:,:,5) = TAGS(:,:,0) + 5000
TAGS(:,:,6) = TAGS(:,:,0) + 6000

! Initialize global parameters
 
CALL INITIALIZE_GLOBAL_VARIABLES
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Initialize radiation 
 
IF (RADIATION) CALL INIT_RADIATION
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Allocate and initialize mesh-specific variables
 
DO NM=1,NMESHES
   IF (PROCESS(NM)==MYID) CALL INITIALIZE_MESH_VARIABLES(NM)
ENDDO
CALL MPI_ALLREDUCE(PROCESS_STOP_STATUS,STOP_STATUS,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,IERR)
IF (STOP_STATUS/=NO_STOP) CALL END_FDS
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Allocate and initialize mesh variable exchange arrays
 
DO NM=1,NMESHES
   IF (PROCESS(NM)==MYID) CALL INITIALIZE_MESH_EXCHANGE(NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
CALL MPI_ALLGATHERV(I_MIN(1,DISPLS(MYID)+1),COUNTS2D(MYID),MPI_INTEGER,I_MIN,COUNTS2D,DISPLS2D,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHERV(I_MAX(1,DISPLS(MYID)+1),COUNTS2D(MYID),MPI_INTEGER,I_MAX,COUNTS2D,DISPLS2D,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHERV(J_MIN(1,DISPLS(MYID)+1),COUNTS2D(MYID),MPI_INTEGER,J_MIN,COUNTS2D,DISPLS2D,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHERV(J_MAX(1,DISPLS(MYID)+1),COUNTS2D(MYID),MPI_INTEGER,J_MAX,COUNTS2D,DISPLS2D,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHERV(K_MIN(1,DISPLS(MYID)+1),COUNTS2D(MYID),MPI_INTEGER,K_MIN,COUNTS2D,DISPLS2D,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHERV(K_MAX(1,DISPLS(MYID)+1),COUNTS2D(MYID),MPI_INTEGER,K_MAX,COUNTS2D,DISPLS2D,MPI_INTEGER,MPI_COMM_WORLD,IERR)
CALL MPI_ALLGATHERV(NIC(1,DISPLS(MYID)+1),  COUNTS2D(MYID),MPI_INTEGER,NIC,  COUNTS2D,DISPLS2D,MPI_INTEGER,MPI_COMM_WORLD,IERR)
 
I_MIN = TRANSPOSE(I_MIN)
I_MAX = TRANSPOSE(I_MAX)
J_MIN = TRANSPOSE(J_MIN)
J_MAX = TRANSPOSE(J_MAX)
K_MIN = TRANSPOSE(K_MIN)
K_MAX = TRANSPOSE(K_MAX)
NIC   = TRANSPOSE(NIC)
 
DO NM=1,NMESHES
   IF (PROCESS(NM)==MYID) CALL DOUBLE_CHECK(NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Initialize Mesh Exchange Arrays (All Nodes)

N_REQ=0  ! Counter for MPI requests
CALL POST_RECEIVES(0)
CALL MESH_EXCHANGE(0)
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)

! Initialize turb arrays

IF (PERIODIC_TEST==2 .OR. DYNSMAG .OR. CHECK_KINETIC_ENERGY) THEN
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      IF (.NOT.EVACUATION_ONLY(NM)) CALL INIT_TURB_ARRAYS(NM)
   ENDDO
ENDIF

! Initialize the flow field with random noise to eliminate false symmetries

IF (NOISE .OR. PERIODIC_TEST>0) THEN
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      IF (NOISE) CALL INITIAL_NOISE(NM)
      IF (PERIODIC_TEST==1) CALL ANALYTICAL_SOLUTION(NM)
      IF (PERIODIC_TEST==2) THEN
         CALL SANDIA_DAT(NM)
         IF (NM==1) CALL SPECTRAL_OUTPUT(1)
      ENDIF
   ENDDO
   CALL POST_RECEIVES(6)
   CALL MESH_EXCHANGE(6)
   PREDICTOR = .FALSE.
   CORRECTOR = .TRUE.
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      CALL MATCH_VELOCITY(NM)
      CALL VELOCITY_BC(T_BEGIN,NM)
   ENDDO
ENDIF

! Potentially read data from a previous calculation 
 
DO NM=1,NMESHES
   IF (RESTART .AND. PROCESS(NM)==MYID) CALL READ_RESTART(T(NM),NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Initialize output files containing global data (Master Node Only)
 
IF (MYID==0) CALL INITIALIZE_GLOBAL_DUMPS
CALL INIT_EVAC_DUMPS
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
 
! Initialize output files that are mesh-specific
 
DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID) CYCLE
   CALL INITIALIZE_MESH_DUMPS(NM)
   CALL INITIALIZE_DROPLETS(NM)
   CALL INSERT_DROPLETS_AND_PARTICLES(T_BEGIN,NM)
   CALL INITIALIZE_RAISED_VEG(NM)
ENDDO
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)

! Initialize EVACuation routines

IF (ANY(EVACUATION_GRID)) THEN
   CALL INITIALIZE_EVAC
   IF (MYID==EVAC_PROCESS) CALL INIT_EVAC_GROUPS
ENDIF

! Write out character strings to .smv file
 
CALL WRITE_STRINGS
 
! Make an initial dump of ambient values

IF (.NOT.RESTART) THEN
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      CALL UPDATE_GLOBAL_OUTPUTS(T(NM),NM)      
      CALL DUMP_MESH_OUTPUTS(T(NM),NM)
   ENDDO
ENDIF

! Ensure the time is known to all meshes

CALL MPI_ALLGATHERV(T(DISPLS(MYID)+1),COUNTS(MYID),MPI_DOUBLE_PRECISION,T,COUNTS,DISPLS,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)

! Make an initial dump of global output quantities

IF (.NOT.RESTART) THEN
   CALL UPDATE_CONTROLS(T)
   CALL DUMP_GLOBAL_OUTPUTS(T(1))
ENDIF
 
! Check for changes in VENT or OBSTruction control and device status at t=T_BEGIN


IF (.NOT.RESTART) THEN
   DO NM=1,NMESHES
      IF (PROCESS(NM)==MYID) CALL OPEN_AND_CLOSE(T(NM),NM)  
   ENDDO
ENDIF

IF (ANY(EVACUATION_GRID)) THEN
   CALL MPI_ALLREDUCE(PROCESS_STOP_STATUS,STOP_STATUS,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,IERR)
   IF (STOP_STATUS/=NO_STOP) CALL END_FDS
   CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
ENDIF

! Start the clock for time stepping

WALL_CLOCK_START_ITERATIONS = WALL_CLOCK_TIME()

!***********************************************************************************************************************************
!                                                           MAIN TIMESTEPPING LOOP
!***********************************************************************************************************************************
 
MAIN_LOOP: DO  
 
   ICYC  = ICYC + 1   ! Time step iterations
   N_REQ = 0          ! Counter for MPI requests

   IF (MOD(ICYC,3)==0 .AND. TIMING) THEN
      MPI_WALL_TIME = MPI_WTIME() - MPI_WALL_TIME_START
      WRITE(LU_ERR,'(A,I3,A,I6,A,F12.4)')  ' Process ',MYID,' starts iteration',ICYC,' at ', MPI_WALL_TIME
   ENDIF

   ! Determine if radiation or evacuation info is to be exchanged this time step
 
   EXCHANGE_RADIATION = .FALSE.
   IF (RADIATION) THEN
      IF (MOD(ICYC,ANGLE_INCREMENT*TIME_STEP_INCREMENT)==0) EXCHANGE_RADIATION = .TRUE.
   ENDIF

   EXCHANGE_EVACUATION = .FALSE.

   ! Synchronize clocks
   
   CALL MPI_ALLGATHERV(T(DISPLS(MYID)+1),COUNTS(MYID),MPI_DOUBLE_PRECISION,T,COUNTS,DISPLS,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)

   ! Check for program stops
 
   CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)
   INQUIRE(FILE=FN_STOP,EXIST=EX)
   CALL MPI_ALLREDUCE(EX,EX_GLOBAL,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR)
   EX = EX_GLOBAL
   IF (EX) MESH_STOP_STATUS = USER_STOP
   IF (EX) WRITE(LU_ERR,'(A,A)') PNAME(1:PNAMELEN),' read the stop file'

   ! Figure out fastest and slowest meshes
 
   T_MAX = MAXVAL(T)
   T_MIN = MINVAL(T)
   DO NM=1,NMESHES
      IF (MESH_STOP_STATUS(NM)/=NO_STOP) PROCESS_STOP_STATUS = MESH_STOP_STATUS(NM)
   ENDDO
 
   ! Determine time step
 
   IF (SYNCHRONIZE) THEN
      DO NM=1,NMESHES
         IF (PROCESS(NM)==MYID) DT_SYNC(NM) = MESHES(NM)%DT_NEXT
      ENDDO

      CALL MPI_ALLGATHERV(DT_SYNC(DISPLS(MYID)+1),COUNTS(MYID),MPI_DOUBLE_PRECISION,DT_SYNC,COUNTS,DISPLS, &
                          MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)

      DO NM=1,NMESHES
         IF (PROCESS(NM)/=MYID) CYCLE
         IF (SYNC_TIME_STEP(NM)) THEN
            MESHES(NM)%DT_NEXT = MINVAL(DT_SYNC,MASK=SYNC_TIME_STEP)
            T(NM) = MINVAL(T,MASK=SYNC_TIME_STEP)
            ACTIVE_MESH(NM) = .TRUE.
         ELSE
            ACTIVE_MESH(NM) = .FALSE.
            IF (T(NM)+MESHES(NM)%DT_NEXT <= T_MAX)  ACTIVE_MESH(NM) = .TRUE.
            IF (PROCESS_STOP_STATUS/=NO_STOP) ACTIVE_MESH(NM) = .TRUE.
         ENDIF
      ENDDO
   ELSE
      ACTIVE_MESH = .FALSE.
      DO NM=1,NMESHES
         IF (PROCESS(NM)/=MYID) CYCLE
         IF (T(NM)+MESHES(NM)%DT_NEXT <= T_MAX) ACTIVE_MESH(NM) = .TRUE.
         IF (PROCESS_STOP_STATUS/=NO_STOP) ACTIVE_MESH(NM) = .TRUE.
      ENDDO
   ENDIF
 
   ! Determine when to dump out diagnostics to the .out file
 
   DIAGNOSTICS = .FALSE.
   LO10 = LOG10(REAL(MAX(1,ABS(ICYC)),EB))
   IF (MOD(ICYC,10**LO10)==0 .OR. MOD(ICYC,100)==0 .OR. T_MIN>=T_END .OR. PROCESS_STOP_STATUS/=NO_STOP) DIAGNOSTICS = .TRUE.

   ! Give every processor the full ACTIVE_MESH array
 
   CALL MPI_ALLGATHERV(ACTIVE_MESH(DISPLS(MYID)+1),COUNTS(MYID),MPI_LOGICAL,ACTIVE_MESH,COUNTS,DISPLS,MPI_LOGICAL, &
                       MPI_COMM_WORLD,IERR)

   ! If no meshes are due to be updated, update them all
 
   IF (ALL(.NOT.ACTIVE_MESH)) ACTIVE_MESH = .TRUE.

   ! If evacuation, set up special time iteration parameters

   CALL EVAC_MAIN_LOOP

   !============================================================================================================================
   !                                          Start of Predictor part of time step
   !============================================================================================================================
 
   PREDICTOR = .TRUE.
   CORRECTOR = .FALSE.

   ! Diagnostic calls

   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      IF (DEBUG) WRITE(LU_ERR,'(A,I8,A,I3,A,L2)') 'Cycle ',ICYC,', Mesh ',NM, ' status=',ACTIVE_MESH(NM)

      IF (MOD(ICYC,3)==0 .AND. TIMING .AND. ACTIVE_MESH(NM)) THEN
         MPI_WALL_TIME = MPI_WTIME() - MPI_WALL_TIME_START
         WRITE(LU_ERR,'(A,I3,A,F12.4)')  ' Mesh ',NM,' is active at ', MPI_WALL_TIME
      ENDIF
   ENDDO
 
   ! Begin the finite differencing of the PREDICTOR step
 
   COMPUTE_FINITE_DIFFERENCES_1: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE COMPUTE_FINITE_DIFFERENCES_1
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_FINITE_DIFFERENCES_1
      MESHES(NM)%DT = MESHES(NM)%DT_NEXT
      NTCYC(NM)     = NTCYC(NM) + 1
 
      CALL INSERT_DROPLETS_AND_PARTICLES(T(NM),NM)
      CALL COMPUTE_VELOCITY_FLUX(T(NM),NM,1)
      IF (FLUX_LIMITER<0 .AND. (.NOT.ISOTHERMAL .OR. N_SPECIES>0)) CALL MASS_FINITE_DIFFERENCES(NM)
   ENDDO COMPUTE_FINITE_DIFFERENCES_1

   ! Estimate quantities at next time step, and decrease/increase time step if necessary based on CFL condition

   FIRST_PASS = .TRUE.
   
   ! Initialize density limits

   RHO_LOWER_GLOBAL = RHOA
   RHO_UPPER_GLOBAL = RHOA
 
   CHANGE_TIME_STEP_LOOP: DO

      IF (FIRST_PASS .OR. SYNCHRONIZE) CALL POST_RECEIVES(1)  

      ! Predict density and mass fractions at next time step, and then start the divergence calculation
 
      COMPUTE_DENSITY_LOOP: DO NM=1,NMESHES
         IF (PROCESS(NM)/=MYID)    CYCLE COMPUTE_DENSITY_LOOP
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_DENSITY_LOOP
         IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) THEN
            IF (FLUX_LIMITER>=0) THEN
               CALL SCALARF(NM) ! leave this here
               CALL DENSITY_TVD(NM)
            ELSE
               CALL DENSITY(NM)
            ENDIF
         ENDIF
      ENDDO COMPUTE_DENSITY_LOOP
      
      IF (BAROCLINIC) THEN
         RHO_UPPER_LOCAL = RHO_UPPER_GLOBAL
         RHO_LOWER_LOCAL = RHO_LOWER_GLOBAL
         CALL MPI_ALLREDUCE(RHO_UPPER_LOCAL,RHO_UPPER_GLOBAL,1,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,IERR)
         CALL MPI_ALLREDUCE(RHO_LOWER_LOCAL,RHO_LOWER_GLOBAL,1,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,IERR)
      ENDIF

      ! Exchange density and species mass fractions in interpolated boundaries

      IF (FIRST_PASS .OR. SYNCHRONIZE) THEN
         CALL MESH_EXCHANGE(1)
         IF (PRESSURE_CORRECTION) CALL POST_RECEIVES(2) 
      ENDIF

      ! Do mass and energy boundary conditions, and begin divergence calculation

      COMPUTE_DIVERGENCE_LOOP: DO NM=1,NMESHES
         IF (PROCESS(NM)/=MYID)    CYCLE COMPUTE_DIVERGENCE_LOOP
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_DIVERGENCE_LOOP                
         IF (FIRST_PASS) CALL COMPUTE_VELOCITY_FLUX(T(NM),NM,2)
         IF (FIRST_PASS) CALL UPDATE_PARTICLES(T(NM),NM)
         IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) CALL WALL_BC(T(NM),NM)
         CALL DIVERGENCE_PART_1(T(NM),NM)
      ENDDO COMPUTE_DIVERGENCE_LOOP

      ! 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

      ! Finish the divergence calculation and solve for the pressure at the current time step


      COMPUTE_PRESSURE_LOOP: DO NM=1,NMESHES
         IF (PROCESS(NM)/=MYID)    CYCLE COMPUTE_PRESSURE_LOOP
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_PRESSURE_LOOP
         CALL DIVERGENCE_PART_2(NM)
         CALL PRESSURE_SOLVER(T(NM),NM)
         CALL EVAC_PRESSURE_LOOP(NM)
      ENDDO COMPUTE_PRESSURE_LOOP

      ! Optional pressure correction

      IF (PRESSURE_CORRECTION .AND. (FIRST_PASS .OR. SYNCHRONIZE)) THEN
         CALL MESH_EXCHANGE(2)
         CALL CORRECT_PRESSURE
      ENDIF

      ! Predict the velocity components at the next time step

      PREDICT_VELOCITY_LOOP: DO NM=1,NMESHES
         IF (PROCESS(NM)/=MYID)    CYCLE PREDICT_VELOCITY_LOOP
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE PREDICT_VELOCITY_LOOP
         CALL VELOCITY_PREDICTOR(NM,MESH_STOP_STATUS(NM))
      ENDDO PREDICT_VELOCITY_LOOP

      ! Exchange information about the time step status, and if need be, repeat the CHANGE_TIME_STEP_LOOP

      SYNCHRONIZE_ONLY: IF (SYNCHRONIZE) THEN
         DISP = DISPLS(MYID)+1
         CALL MPI_ALLGATHERV(CHANGE_TIME_STEP(DISP),COUNTS(MYID),MPI_LOGICAL,CHANGE_TIME_STEP,COUNTS,DISPLS, &
                             MPI_LOGICAL,MPI_COMM_WORLD,IERR)
         CALL MPI_ALLGATHERV(MESH_STOP_STATUS(DISP),COUNTS(MYID),MPI_INTEGER,MESH_STOP_STATUS,COUNTS,DISPLS, &
                            MPI_INTEGER,MPI_COMM_WORLD,IERR)
         IF (ANY(MESH_STOP_STATUS/=NO_STOP)) THEN
            IF (ANY(MESH_STOP_STATUS==INSTABILITY_STOP)) THEN
               PROCESS_STOP_STATUS = INSTABILITY_STOP
               DIAGNOSTICS = .TRUE.
            ENDIF
            EXIT CHANGE_TIME_STEP_LOOP
         ENDIF
         IF (ANY(CHANGE_TIME_STEP)) THEN
            CHANGE_TIME_STEP = .TRUE.
            DO NM=1,NMESHES
               IF (PROCESS(NM)/=MYID) CYCLE
               DT_SYNC(NM)     = MESHES(NM)%DT
               DT_NEXT_SYNC(NM) = MESHES(NM)%DT_NEXT
            ENDDO
            CALL MPI_ALLGATHERV(DT_SYNC(DISP),COUNTS(MYID),MPI_DOUBLE_PRECISION,DT_SYNC,COUNTS,DISPLS, &
                                MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
            CALL MPI_ALLGATHERV(DT_NEXT_SYNC(DISP),COUNTS(MYID),MPI_DOUBLE_PRECISION,DT_NEXT_SYNC,COUNTS,DISPLS, &
                                MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
            DO NM=1,NMESHES
               IF (PROCESS(NM)/=MYID) CYCLE
               IF (EVACUATION_ONLY(NM)) CHANGE_TIME_STEP(NM) = .FALSE.
               MESHES(NM)%DT_NEXT = MINVAL(DT_NEXT_SYNC,MASK=SYNC_TIME_STEP)
               MESHES(NM)%DT     = MINVAL(DT_SYNC,MASK=SYNC_TIME_STEP)
            ENDDO
         ENDIF
      ENDIF SYNCHRONIZE_ONLY

      DO NM=1,NMESHES
         IF (PROCESS(NM)==MYID .AND. MESH_STOP_STATUS(NM)/=NO_STOP) THEN
            PROCESS_STOP_STATUS = MESH_STOP_STATUS(NM)
            IF (PROCESS_STOP_STATUS==INSTABILITY_STOP) DIAGNOSTICS = .TRUE.
            EXIT CHANGE_TIME_STEP_LOOP
         ENDIF
      ENDDO
 
      IF (.NOT.ANY(CHANGE_TIME_STEP)) EXIT CHANGE_TIME_STEP_LOOP
 
      FIRST_PASS = .FALSE.

   ENDDO CHANGE_TIME_STEP_LOOP
 
   CHANGE_TIME_STEP = .FALSE.
 
   ! Exchange velocity and pressures at interpolated boundaries

   CALL POST_RECEIVES(3)
   CALL MESH_EXCHANGE(3)

   ! Force normal components of velocity to match at interpolated boundaries

   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID .OR. .NOT.ACTIVE_MESH(NM)) CYCLE
      CALL MATCH_VELOCITY(NM)
   ENDDO

   ! Apply tangential velocity boundary conditions

   VELOCITY_BC_LOOP: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE VELOCITY_BC_LOOP
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE VELOCITY_BC_LOOP
      CALL VELOCITY_BC(T(NM),NM)
   ENDDO VELOCITY_BC_LOOP

   ! Advance the time to start the CORRECTOR step
 
   UPDATE_TIME: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE UPDATE_TIME
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE UPDATE_TIME
      T(NM) = T(NM) + MESHES(NM)%DT  
   ENDDO UPDATE_TIME

   ! Ensure the time is known to all meshes

   CALL MPI_ALLGATHERV(T(DISPLS(MYID)+1),COUNTS(MYID),MPI_DOUBLE_PRECISION,T,COUNTS,DISPLS,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
   T_MIN = MINVAL(T)
   T_MAX = MAXVAL(T)
 
   !============================================================================================================================
   !                                          Start of Corrector part of time step
   !============================================================================================================================
 
   CORRECTOR = .TRUE.
   PREDICTOR = .FALSE.

   CALL POST_RECEIVES(4)

   ! Finite differences for mass and momentum equations for the second half of the time step

   COMPUTE_FINITE_DIFFERENCES_2: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE COMPUTE_FINITE_DIFFERENCES_2
      CALL OPEN_AND_CLOSE(T(NM),NM)   
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_FINITE_DIFFERENCES_2
      CALL COMPUTE_VELOCITY_FLUX(T(NM),NM,1)
      IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) THEN
         IF (FLUX_LIMITER>=0) THEN
            CALL SCALARF(NM) ! leave this here
            CALL DENSITY_TVD(NM)
         ELSE
            CALL MASS_FINITE_DIFFERENCES(NM)
            CALL DENSITY(NM)
         ENDIF
      ENDIF
   ENDDO COMPUTE_FINITE_DIFFERENCES_2

   ! Exchange min and max density needed by BAROCLINIC correction
   
   IF (BAROCLINIC) THEN
      RHO_UPPER_LOCAL = RHO_UPPER_GLOBAL
      RHO_LOWER_LOCAL = RHO_LOWER_GLOBAL
      CALL MPI_ALLREDUCE(RHO_UPPER_LOCAL,RHO_UPPER_GLOBAL,1,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,IERR)
      CALL MPI_ALLREDUCE(RHO_LOWER_LOCAL,RHO_LOWER_GLOBAL,1,MPI_DOUBLE_PRECISION,MPI_MIN,MPI_COMM_WORLD,IERR)
   ENDIF

   ! Exchange density and mass species

   CALL MESH_EXCHANGE(4)
   IF (PRESSURE_CORRECTION) CALL POST_RECEIVES(5)  

   ! Apply mass and species boundary conditions, update radiation, particles, and re-compute divergence

   COMPUTE_DIVERGENCE_2: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE COMPUTE_DIVERGENCE_2
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_DIVERGENCE_2
      CALL COMPUTE_VELOCITY_FLUX(T(NM),NM,2)
      CALL UPDATE_PARTICLES(T(NM),NM)
      IF(.NOT. WIND_ONLY) CALL RAISED_VEG_MASS_ENERGY_TRANSFER(T(NM),NM)
      IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) THEN
         IF (N_REACTIONS > 0) CALL COMBUSTION (NM)
         IF (SOOT_DEPOSITION) CALL CALC_SOOT_DEPOSITION(NM)
         CALL WALL_BC(T(NM),NM)
         CALL COMPUTE_RADIATION(T(NM),NM)
      ENDIF
      CALL DIVERGENCE_PART_1(T(NM),NM)
   ENDDO COMPUTE_DIVERGENCE_2

   IF (N_ZONE>0) CALL EXCHANGE_DIVERGENCE_INFO

   COMPUTE_PRESSURE_2: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE COMPUTE_PRESSURE_2
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_PRESSURE_2
      CALL DIVERGENCE_PART_2(NM)
      CALL PRESSURE_SOLVER(T(NM),NM)
      CALL EVAC_PRESSURE_LOOP(NM)
   ENDDO COMPUTE_PRESSURE_2

   ! Optional pressure correction

   IF (PRESSURE_CORRECTION) THEN
      CALL MESH_EXCHANGE(5)
      CALL CORRECT_PRESSURE
   ENDIF

   ! Set up the last big exchange of info

   CALL EVAC_MESH_EXCHANGE(T_EVAC,T_EVAC_SAVE,I_EVAC,ICYC,EXCHANGE_EVACUATION,0)
   CALL POST_RECEIVES(6) 

   ! Correct the velocity

   CORRECT_VELOCITY_LOOP: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE CORRECT_VELOCITY_LOOP
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE CORRECT_VELOCITY_LOOP
      CALL VELOCITY_CORRECTOR(NM)
      IF (DIAGNOSTICS) CALL CHECK_DIVERGENCE(NM)
   ENDDO CORRECT_VELOCITY_LOOP

   IF (CHECK_VOLUME_FLOW) CALL COMPUTE_VOLUME_FLOW
 
   ! Exchange velocity and pressure at interpolated boundaries
 
   CALL MESH_EXCHANGE(6)


   ! Force normal components of velocity to match at interpolated boundaries

   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID .OR. .NOT.ACTIVE_MESH(NM)) CYCLE
      IF (EVACUATION_ONLY(NM)) CYCLE
      CALL MATCH_VELOCITY(NM)
   ENDDO

   ! Apply velocity boundary conditions, and update values of HRR, DEVC, etc.

   VELOCITY_BC_LOOP_2: DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID)    CYCLE VELOCITY_BC_LOOP_2
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE VELOCITY_BC_LOOP_2
      CALL VELOCITY_BC(T(NM),NM)
      CALL UPDATE_GLOBAL_OUTPUTS(T(NM),NM)
   ENDDO VELOCITY_BC_LOOP_2

   ! Exchange EVAC information among meshes

   CALL EVAC_EXCHANGE

   ! Check for dumping end of timestep outputs

   CALL UPDATE_CONTROLS(T)
   DO NM=1,NMESHES
      IF (PROCESS(NM)==MYID .AND. ACTIVE_MESH(NM)) CALL DUMP_MESH_OUTPUTS(T(NM),NM) 
   ENDDO
   CALL DUMP_GLOBAL_OUTPUTS(T_MIN)
 
   ! Dump out diagnostics

   IF (DIAGNOSTICS .OR. T_MIN>=T_END) THEN
      CALL WRITE_STRINGS
      CALL EXCHANGE_DIAGNOSTICS
      IF (MYID==0) CALL WRITE_DIAGNOSTICS(T)
   ENDIF
 
   ! Flush output file buffers

   IF (T_MIN>=FLUSH_CLOCK .AND. FLUSH_FILE_BUFFERS) THEN
      IF (MYID==0) CALL FLUSH_GLOBAL_BUFFERS
      IF (MYID==EVAC_PROCESS) CALL FLUSH_EVACUATION_BUFFERS
      DO NM=1,NMESHES
         IF (PROCESS(NM)==MYID) CALL FLUSH_LOCAL_BUFFERS(NM)
      ENDDO
      FLUSH_CLOCK = FLUSH_CLOCK + DT_FLUSH
   ENDIF

   ! Stop the run

   IF (T_MIN>=T_END .OR. PROCESS_STOP_STATUS/=NO_STOP) EXIT MAIN_LOOP
 
   ! Diagnostic print out

   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      IF (MOD(ICYC,3) ==0 .AND. TIMING) THEN
         MPI_WALL_TIME = MPI_WTIME() - MPI_WALL_TIME_START
         WRITE(LU_ERR,'(A,I3,A,I6,A,F12.4)')  ' Mesh ',NM,' ends Iteration ',ICYC,' at ', MPI_WALL_TIME
      ENDIF
   ENDDO
 
ENDDO MAIN_LOOP
 
!***********************************************************************************************************************************
!                                                     END OF TIME STEPPING LOOP
!***********************************************************************************************************************************
 
! Gather up timings for the meshes

DO NM=1,NMESHES
   IF (PROCESS(NM)==MYID) TUSED(1,NM) = SECOND() - TUSED(1,NM)
ENDDO
CALL MPI_GATHERV(TUSED(1,DISPLS(MYID)+1),COUNTS_TIMERS(MYID),MPI_DOUBLE_PRECISION,TUSED, &
                 COUNTS_TIMERS,DISPLS_TIMERS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
IF (MYID==0) CALL TIMINGS

! Stop the calculation

CALL END_FDS

! The list of subroutines called from the main program follows

CONTAINS


SUBROUTINE END_FDS

! End the calculation gracefully, even if there is an error

CHARACTER(100) :: MESSAGE

CALL MPI_REDUCE(PROCESS_STOP_STATUS,STOP_STATUS,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,IERR)
CALL MPI_FINALIZE(IERR)

IF (MYID==0) THEN
   SELECT CASE(STOP_STATUS)
      CASE(NO_STOP)
         WRITE(MESSAGE,'(A)') 'STOP: FDS completed successfully'
         IF (STATUS_FILES) CLOSE(LU_NOTREADY,STATUS='DELETE')
      CASE(INSTABILITY_STOP) 
         WRITE(MESSAGE,'(A)') 'STOP: Numerical Instability'
      CASE(USER_STOP) 
         WRITE(MESSAGE,'(A)') 'STOP: FDS stopped by user'
      CASE(SETUP_STOP) 
         WRITE(MESSAGE,'(A)') 'STOP: FDS was improperly set-up'
   END SELECT
   CALL SHUTDOWN(MESSAGE)
ENDIF

STOP

END SUBROUTINE END_FDS
 
 
 
 
SUBROUTINE EXCHANGE_DIVERGENCE_INFO

! Exchange information mesh to mesh needed for divergence integrals
! First, sum DSUM, PSUM and USUM over all meshes controlled by the active process, then reduce over all processes

INTEGER :: IPZ,IOPZ,IOPZ2
REAL(EB) :: TNOW

TNOW = SECOND()

CONNECTED_ZONES_LOCAL = .FALSE.

DO IPZ=1,N_ZONE
   DSUM_ALL_LOCAL(IPZ) = 0._EB
   PSUM_ALL_LOCAL(IPZ) = 0._EB
   USUM_ALL_LOCAL(IPZ) = 0._EB
   ASUM_ALL_LOCAL(IPZ,:) = 0._EB
   FDS_LEAK_AREA_LOCAL(IPZ,:) = 0._EB
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      IF(EVACUATION_ONLY(NM)) CYCLE
      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)
      ASUM_ALL_LOCAL(IPZ,:) = ASUM_ALL_LOCAL(IPZ,:) + ASUM(IPZ,:,NM)
      FDS_LEAK_AREA_LOCAL(IPZ,:) = FDS_LEAK_AREA_LOCAL(IPZ,:) + FDS_LEAK_AREA(IPZ,:,NM)
      DO IOPZ=0,N_ZONE
         IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) CONNECTED_ZONES_LOCAL(IPZ,IOPZ) = .TRUE.
      ENDDO
   ENDDO
ENDDO

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(ASUM_ALL_LOCAL(1,0),ASUM_ALL(1,0),N_ZONE*(N_ZONE+1),MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
CALL MPI_ALLREDUCE(FDS_LEAK_AREA_LOCAL(0,0),FDS_LEAK_AREA_GLOBAL(0,0),(N_ZONE+1)**2,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)

DO IPZ=1,N_ZONE
   DO NM=1,NMESHES
      IF(EVACUATION_ONLY(NM)) CYCLE
      DSUM(IPZ,NM) = DSUM_ALL(IPZ)
      PSUM(IPZ,NM) = PSUM_ALL(IPZ)
      USUM(IPZ,NM) = USUM_ALL(IPZ)
      ASUM(IPZ,:,NM) = ASUM_ALL(IPZ,:)
      DO IOPZ=0,N_ZONE
         IF (FDS_LEAK_AREA_GLOBAL(IPZ,IOPZ)>0._EB) &
            FDS_LEAK_AREA_RATIO(IPZ,IOPZ,NM) = FDS_LEAK_AREA(IPZ,IOPZ,NM)/FDS_LEAK_AREA_GLOBAL(IPZ,IOPZ)
      ENDDO
      CONNECTED_ZONES(IPZ,:,NM) = CONNECTED_ZONES_GLOBAL(IPZ,:)
   ENDDO
ENDDO

! Connect zones to others which are not directly connected

DO NM=1,NMESHES
   IF(EVACUATION_ONLY(NM)) CYCLE 
   DO IPZ=1,N_ZONE
      DO IOPZ=1,N_ZONE
         IF (IOPZ==IPZ) CYCLE
         IF (CONNECTED_ZONES(IPZ,IOPZ,NM)) THEN
            DO IOPZ2=0,N_ZONE
               IF (IOPZ==IOPZ2) CYCLE
               IF (CONNECTED_ZONES(IOPZ,IOPZ2,NM)) CONNECTED_ZONES(IPZ,IOPZ2,NM) = .TRUE.
               IF (CONNECTED_ZONES(IOPZ,IOPZ2,NM)) CONNECTED_ZONES(IOPZ2,IPZ,NM) = .TRUE.
            ENDDO
         ENDIF
      ENDDO
   ENDDO
ENDDO

TUSED(2,:)=TUSED(2,:) + SECOND() - TNOW
END SUBROUTINE EXCHANGE_DIVERGENCE_INFO

 
SUBROUTINE INITIALIZE_MESH_EXCHANGE(NM)
 
! Create arrays by which info is to exchanged across meshes
 
INTEGER IMIN,IMAX,JMIN,JMAX,KMIN,KMAX,NOM,IOR,IW
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M2,M
LOGICAL FOUND
 
M=>MESHES(NM)
IF (.NOT.EVACUATION_ONLY(NM)) ALLOCATE(MESHES(NM)%OMESH(NMESHES))
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (NOM==NM .AND. PERIODIC_TEST==0) CYCLE OTHER_MESH_LOOP
 
   M2=>MESHES(NOM)
   IMIN=0 
   IMAX=M2%IBP1
   JMIN=0 
   JMAX=M2%JBP1
   KMIN=0 
   KMAX=M2%KBP1
   NIC(NOM,NM) = 0
   FOUND = .FALSE.


   IF (EVACUATION_ONLY(NOM)) CYCLE OTHER_MESH_LOOP ! Issue 257 bug fix

   SEARCH_LOOP: DO IW=1,M%NEWC
      IF (M%IJKW(9,IW)/=NOM) CYCLE SEARCH_LOOP
      NIC(NOM,NM) = NIC(NOM,NM) + 1
      FOUND = .TRUE.
      IOR = M%IJKW(4,IW)
      NOT_PERIODIC: IF (PERIODIC_TEST==0) THEN
         SELECT CASE(IOR)
            CASE( 1) 
               IMIN=MAX(IMIN,M%IJKW(10,IW)-1)
            CASE(-1) 
               IMAX=MIN(IMAX,M%IJKW(13,IW))
            CASE( 2) 
               JMIN=MAX(JMIN,M%IJKW(11,IW)-1)
            CASE(-2) 
               JMAX=MIN(JMAX,M%IJKW(14,IW))
            CASE( 3) 
               KMIN=MAX(KMIN,M%IJKW(12,IW)-1)
            CASE(-3) 
               KMAX=MIN(KMAX,M%IJKW(15,IW))
         END SELECT
      ENDIF NOT_PERIODIC
   ENDDO SEARCH_LOOP
 
   IF ( M2%XS>=M%XS .AND. M2%XF<=M%XF .AND. M2%YS>=M%YS .AND. M2%YF<=M%YF .AND. M2%ZS>=M%ZS .AND. M2%ZF<=M%ZF ) FOUND = .TRUE.
 
   IF (.NOT.FOUND) CYCLE OTHER_MESH_LOOP
 
   I_MIN(NOM,NM) = IMIN
   I_MAX(NOM,NM) = IMAX
   J_MIN(NOM,NM) = JMIN
   J_MAX(NOM,NM) = JMAX
   K_MIN(NOM,NM) = KMIN
   K_MAX(NOM,NM) = KMAX
 
   ALLOCATE(M%OMESH(NOM)% RHO(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%RHOS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%RHO  = RHOA
   M%OMESH(NOM)%RHOS = RHOA
   ALLOCATE(M%OMESH(NOM)% D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%DS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%D  = 0._EB
   M%OMESH(NOM)%DS = 0._EB
   ALLOCATE(M%OMESH(NOM)%  MU(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%MU = 0._EB
   ALLOCATE(M%OMESH(NOM)%   H(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%H = 0._EB
   ALLOCATE(M%OMESH(NOM)%   U(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%  US(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%U  = U0
   M%OMESH(NOM)%US = U0
   ALLOCATE(M%OMESH(NOM)%   V(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%  VS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%V  = V0
   M%OMESH(NOM)%VS = V0
   ALLOCATE(M%OMESH(NOM)%   W(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%  WS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%W  = W0
   M%OMESH(NOM)%WS = W0
   ALLOCATE(M%OMESH(NOM)% DUDT(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)% DVDT(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)% DWDT(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   M%OMESH(NOM)%DUDT = 0._EB
   M%OMESH(NOM)%DVDT = 0._EB
   M%OMESH(NOM)%DWDT = 0._EB
 
   IF (N_SPECIES>0) THEN
      ALLOCATE(M%OMESH(NOM)%  YY(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,N_SPECIES))
      ALLOCATE(M%OMESH(NOM)% YYS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,N_SPECIES))
      DO N=1,N_SPECIES
      M%OMESH(NOM)%YY(:,:,:,N)  = SPECIES(N)%YY0
      M%OMESH(NOM)%YYS(:,:,:,N) = SPECIES(N)%YY0
      ENDDO
   ENDIF
 
   ! Wall arrays
 
   ALLOCATE(M%OMESH(NOM)%IJKW(15,M2%NEWC))
   ALLOCATE(M%OMESH(NOM)%BOUNDARY_TYPE(0:M2%NEWC))
   ALLOCATE(M%OMESH(NOM)%WALL(0:M2%NEWC))
 
   ! Particle and Droplet Orphan Arrays
 
   IF (DROPLET_FILE) THEN
      M%OMESH(NOM)%N_DROP_ORPHANS = 0
      M%OMESH(NOM)%N_DROP_ORPHANS_DIM = N_DROP_ADOPT_MAX
      ALLOCATE(M%OMESH(NOM)%DROPLET(M%OMESH(NOM)%N_DROP_ORPHANS_DIM), STAT=IZERO)
      CALL ChkMemErr('INIT','DROPLET',IZERO)
   ENDIF
 
ENDDO OTHER_MESH_LOOP
 
END SUBROUTINE INITIALIZE_MESH_EXCHANGE
 
 

SUBROUTINE DOUBLE_CHECK(NM)
 
! Double check exchange pairs
 
INTEGER NOM
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M2,M
 
M=>MESHES(NM)
IF (EVACUATION_ONLY(NM)) RETURN
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (NOM==NM .AND. PERIODIC_TEST==0) CYCLE OTHER_MESH_LOOP
   IF (EVACUATION_ONLY(NOM)) CYCLE OTHER_MESH_LOOP ! Issue 257 bug fix

   IF (NIC(NM,NOM)==0 .AND. NIC(NOM,NM)>0) THEN
      M2=>MESHES(NOM)
      ALLOCATE(M%OMESH(NOM)%IJKW(15,M2%NEWC))
      ALLOCATE(M%OMESH(NOM)%BOUNDARY_TYPE(0:M2%NEWC))
      ALLOCATE(M%OMESH(NOM)%WALL(0:M2%NEWC))
   ENDIF
 
ENDDO OTHER_MESH_LOOP
 
END SUBROUTINE DOUBLE_CHECK
 
 

SUBROUTINE POST_RECEIVES(CODE)

USE GLOBAL_CONSTANTS, ONLY: NRA=>NUMBER_RADIATION_ANGLES,NSB=>NUMBER_SPECTRAL_BANDS
INTEGER, INTENT(IN) :: CODE
INTEGER :: SNODE
REAL(EB) :: TNOW

TNOW = SECOND()

MESH_LOOP: DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID) CYCLE MESH_LOOP
   IF (EVACUATION_ONLY(NM)) CYCLE MESH_LOOP 
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (EVACUATION_ONLY(NOM)) CYCLE OTHER_MESH_LOOP
   SNODE = PROCESS(NOM)
   IF (PROCESS(NM)==SNODE) CYCLE OTHER_MESH_LOOP
   IF (NIC(NM,NOM)==0 .AND. NIC(NOM,NM)==0) CYCLE OTHER_MESH_LOOP
   IF (CODE>0 .AND. (.NOT.ACTIVE_MESH(NOM).OR..NOT.ACTIVE_MESH(NM))) CYCLE OTHER_MESH_LOOP
 
   IF (DEBUG) THEN
      WRITE(LU_ERR,'(I3,A,I3,A,I2)') NM,' posting receives from ',NOM,' code=',code
      IF (CODE==0) WRITE(LU_ERR,'(A,I3,A,I3,A,I5)') 'NIC(',NM,',',NOM,')=',NIC(NM,NOM)
   ENDIF
 
   M =>MESHES(NM)
   M4=>MESHES(NOM)
   M3=>MESHES(NM)%OMESH(NOM)
 
   TAG   = TAGS(NM,NOM,CODE)
 
   INITIALIZATION_IF: IF (CODE==0) THEN
 
      IF (NIC(NM,NOM)>0) THEN
         ALLOCATE(M3%REAL_RECV_PKG1(MAX(NIC(NM,NOM),NIC(NOM,NM))*(4+N_SPECIES)+1))
         ALLOCATE(M3%REAL_RECV_PKG2(MAX(NIC(NM,NOM),NIC(NOM,NM))*(8          )+1))
         ALLOCATE(M3%REAL_RECV_PKG3(MAX(NIC(NM,NOM),NIC(NOM,NM))*(4+N_SPECIES)+1))
         ALLOCATE(M3%REAL_RECV_PKG4(MAX(NIC(NM,NOM),NIC(NOM,NM))*(8          )+1))
         ALLOCATE(M3%REAL_RECV_PKG5((NRA*NSB+1)*NIC(NM,NOM)+1))
         ALLOCATE(M3%REAL_RECV_PKG7(MAX(NIC(NM,NOM),NIC(NOM,NM))*(2          )+1))
      ENDIF
 
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%IJKW(1,1),15*M4%NEWC,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
 
      IF (NIC(NM,NOM)>0 .OR. NIC(NOM,NM)>0) THEN
         ALLOCATE(M3%REAL_RECV_PKG6(13*N_DROP_ADOPT_MAX))
         ALLOCATE(M3%INTG_RECV_PKG1( 3*N_DROP_ADOPT_MAX))
         ALLOCATE(M3%LOGI_RECV_PKG1( 2*N_DROP_ADOPT_MAX))
      ENDIF
 
   ENDIF INITIALIZATION_IF

   ! First posting for density and mass fraction
 
   IF (CODE==1 .AND. NIC(NM,NOM)>0) THEN
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%REAL_RECV_PKG1(1),SIZE(M3%REAL_RECV_PKG1),MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF

   ! First posting for PRESSURE_CORRECTION
 
   IF ((CODE==2 .OR. CODE==5) .AND. NIC(NM,NOM)>0) THEN
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%REAL_RECV_PKG7(1),SIZE(M3%REAL_RECV_PKG7),MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF

   ! First posting for pressure

   IF (CODE==3 .AND. NIC(NM,NOM)>0) THEN
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%REAL_RECV_PKG2(1),SIZE(M3%REAL_RECV_PKG2),MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF


   ! Second posting for density and mass fraction
 
   IF (CODE==4 .AND. NIC(NM,NOM)>0) THEN
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%REAL_RECV_PKG3(1),SIZE(M3%REAL_RECV_PKG3),MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF

   ! BOUNDARY_TYPE

   IF (CODE==0 .OR. CODE==6) THEN
      N_REQ = N_REQ+1
      CALL MPI_IRECV(M3%BOUNDARY_TYPE(0),M4%NEWC+1,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD, REQ(N_REQ),IERR)
   ENDIF

   ! Second posting for pressure and only posting for velocities

   IF (CODE==6 .AND. NIC(NM,NOM)>0) THEN
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%REAL_RECV_PKG4(1),SIZE(M3%REAL_RECV_PKG4),MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF

   ! Radiation

   IF (CODE==6 .AND. EXCHANGE_RADIATION .AND. NIC(NM,NOM)>0) THEN
      IWW = NIC(NM,NOM)
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%REAL_RECV_PKG5(1),(NRA*NSB+1)*IWW+1,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF
 
   ! Droplets
 
   IF (CODE==6 .AND. DROPLET_FILE .AND. (NIC(NM,NOM)>0 .OR. NIC(NOM,NM)>0)) THEN
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%N_DROP_ADOPT,1,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      BUFFER_SIZE=13*N_DROP_ADOPT_MAX
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%REAL_RECV_PKG6(1),BUFFER_SIZE,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      BUFFER_SIZE=3*N_DROP_ADOPT_MAX
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%INTG_RECV_PKG1(1),BUFFER_SIZE,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      BUFFER_SIZE=2*N_DROP_ADOPT_MAX
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M3%LOGI_RECV_PKG1(1),BUFFER_SIZE,MPI_LOGICAL,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF

ENDDO OTHER_MESH_LOOP

ENDDO MESH_LOOP

   ! Receive EVACuation information

DO NOM=1,NMESHES
   SNODE = PROCESS(NOM)
   IF (CODE==6 .AND. EXCHANGE_EVACUATION .AND. MYID==EVAC_PROCESS .AND. .NOT.EVACUATION_ONLY(NOM)) THEN
      M4=>MESHES(NOM)
      TAG = NOM*(EVAC_PROCESS+1)*CODE*10
      IWW = (M4%IBAR+2)*(M4%JBAR+2)*(M4%KBAR+2)
      IF (N_SPECIES>0) THEN
         N_REQ=N_REQ+1
         CALL MPI_IRECV(M4%YY(0,0,0,1),IWW*N_SPECIES,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ENDIF
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M4%RHO(0,0,0),IWW,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M4%RSUM(0,0,0),IWW,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M4%TMP(0,0,0),IWW,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M4%UII(0,0,0),IWW,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M4%CELL_INDEX(0,0,0),IWW,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      IWW = MAXVAL(M4%CELL_INDEX)
      N_REQ=N_REQ+1
      CALL MPI_IRECV(M4%SOLID(0),IWW,MPI_LOGICAL,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF
ENDDO
 
TUSED(11,:)=TUSED(11,:) + SECOND() - TNOW
END SUBROUTINE POST_RECEIVES
 
 
 
SUBROUTINE MESH_EXCHANGE(CODE)
 
! Exchange Information between Meshes
 
USE GLOBAL_CONSTANTS, ONLY: NRA=>NUMBER_RADIATION_ANGLES,NSB=>NUMBER_SPECTRAL_BANDS
REAL(EB) :: TNOW
INTEGER, INTENT(IN) :: CODE
INTEGER NM,II,JJ,KK,LL,NC,N,NN,RNODE,SNODE,IMIN,IMAX,JMIN,JMAX,KMIN,KMAX
INTEGER :: NN1,NN2
 
TNOW = SECOND()
 
SENDING_MESH_LOOP: DO NM=1,NMESHES

   IF (PROCESS(NM)/=MYID) CYCLE SENDING_MESH_LOOP
   IF (EVACUATION_ONLY(NM)) CYCLE SENDING_MESH_LOOP ! Issue 257 bug fix
 
   IF (MOD(ICYC,3)==0 .AND. TIMING .AND. ACTIVE_MESH(NM)) THEN
      MPI_WALL_TIME = MPI_WTIME() - MPI_WALL_TIME_START
      WRITE(LU_ERR,'(A,I3,A,I1,A,F12.4)')  ' Mesh ',NM,' enters Mesh Exchange ',CODE,' at ', MPI_WALL_TIME
   ENDIF

! Information about Mesh NM is packed into SEND packages and shipped out to the other meshes (machines) via MPI
 
RECEIVING_MESH_LOOP: DO NOM=1,NMESHES
 
   SNODE = PROCESS(NOM)
   RNODE = PROCESS(NM)

   IF (EVACUATION_ONLY(NOM)) CYCLE RECEIVING_MESH_LOOP ! Issue 257 bug fix
   IF (NIC(NOM,NM)==0 .AND. NIC(NM,NOM)==0)  CYCLE RECEIVING_MESH_LOOP
 
   IF (CODE>0) THEN
      IF (.NOT.ACTIVE_MESH(NM) .OR. .NOT.ACTIVE_MESH(NOM))  CYCLE RECEIVING_MESH_LOOP
   ENDIF
 
   IF (DEBUG) THEN
      WRITE(LU_ERR,'(I3,A,I3,A,I2,A,I5)') NM,' sending data to ',NOM,' code=',CODE, ' tag=',TAGS(NM,NOM,CODE)
   ENDIF
 
   M =>MESHES(NM)
   M3=>MESHES(NM)%OMESH(NOM)
   M4=>MESHES(NOM)
 
   TAG = TAGS(NM,NOM,CODE)
   IMIN = I_MIN(NOM,NM)
   IMAX = I_MAX(NOM,NM)
   JMIN = J_MIN(NOM,NM)
   JMAX = J_MAX(NOM,NM)
   KMIN = K_MIN(NOM,NM)
   KMAX = K_MAX(NOM,NM)

   ! Initial information exchange
 
   INITIALIZE_SEND_IF: IF (CODE==0) THEN
 
      IF (NIC(NOM,NM)>0 .AND. RNODE/=SNODE) THEN
         ALLOCATE(M3%REAL_SEND_PKG1(MAX(NIC(NOM,NM),NIC(NM,NOM))*(4+N_SPECIES)+1))
         ALLOCATE(M3%REAL_SEND_PKG2(MAX(NIC(NOM,NM),NIC(NM,NOM))*(8          )+1))
         ALLOCATE(M3%REAL_SEND_PKG3(MAX(NIC(NOM,NM),NIC(NM,NOM))*(4+N_SPECIES)+1))
         ALLOCATE(M3%REAL_SEND_PKG4(MAX(NIC(NOM,NM),NIC(NM,NOM))*(8          )+1))
         ALLOCATE(M3%REAL_SEND_PKG5((NRA*NSB+1)*NIC(NOM,NM)+1))
         ALLOCATE(M3%REAL_SEND_PKG7(MAX(NIC(NOM,NM),NIC(NM,NOM))*(2          )+1))
      ENDIF
 
      IF (RNODE/=SNODE) THEN
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M%IJKW(1,1),15*M%NEWC, MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         M%OMESH(NOM)%IJKW = M4%IJKW(:,1:M4%NEWC)
      ENDIF
 
      IF (DROPLET_FILE .AND.  (NIC(NOM,NM)>0 .OR. NIC(NM,NOM)>0) .AND. RNODE/=SNODE) THEN
         ALLOCATE(M3%REAL_SEND_PKG6(13*N_DROP_ADOPT_MAX))
         ALLOCATE(M3%INTG_SEND_PKG1( 3*N_DROP_ADOPT_MAX))
         ALLOCATE(M3%LOGI_SEND_PKG1( 2*N_DROP_ADOPT_MAX))
      ENDIF
 
   ENDIF INITIALIZE_SEND_IF

   ! Exchange of density and species mass fractions following the PREDICTOR update

   IF (CODE==1 .AND. NIC(NOM,NM)>0) THEN
      IF (RNODE/=SNODE) THEN
         LL = 0
         IWW = 0
         PACK_REAL_SEND_PKG1: DO IW=1,M4%NEWC
            IF (M3%IJKW(9,IW)/=NM) CYCLE PACK_REAL_SEND_PKG1
            DO KK=M3%IJKW(12,IW),M3%IJKW(15,IW)
               DO JJ=M3%IJKW(11,IW),M3%IJKW(14,IW)
                  DO II=M3%IJKW(10,IW),M3%IJKW(13,IW)
                     IWW = IWW + 1
                     M3%REAL_SEND_PKG1(LL+1) = REAL(IW,EB)
                     M3%REAL_SEND_PKG1(LL+2) = M%RHOS(II,JJ,KK)
                     M3%REAL_SEND_PKG1(LL+3) = M%MU(II,JJ,KK)
                     M3%REAL_SEND_PKG1(LL+4) = M%D(II,JJ,KK)
                     IF (N_SPECIES>0) M3%REAL_SEND_PKG1(LL+5:LL+4+N_SPECIES) = M%YYS(II,JJ,KK,1:N_SPECIES)
                     LL = LL+4+N_SPECIES
                  ENDDO
               ENDDO
            ENDDO
         ENDDO PACK_REAL_SEND_PKG1
         M3%REAL_SEND_PKG1(IWW*(4+N_SPECIES)+1) = -999.0_EB
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%REAL_SEND_PKG1(1),IWW*(4+N_SPECIES)+1,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         M2=>MESHES(NOM)%OMESH(NM)
         M2%RHOS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)= M%RHOS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%MU(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)  = M%MU(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)   = M%D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         IF (N_SPECIES>0) M2%YYS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,1:N_SPECIES)= M%YYS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,1:N_SPECIES)
      ENDIF
   ENDIF

   ! Exchange velocity/pressure info for PRESSURE_CORRECTION

   IF ((CODE==2 .OR. CODE==5) .AND. NIC(NOM,NM)>0) THEN
      IF (RNODE/=SNODE) THEN
         LL = 0
         IWW = 0
         PACK_REAL_SEND_PKG7: DO IW=1,M4%NEWC
            IF (M3%IJKW(9,IW)/=NM) CYCLE PACK_REAL_SEND_PKG7
            DO KK=M3%IJKW(12,IW),M3%IJKW(15,IW)
               DO JJ=M3%IJKW(11,IW),M3%IJKW(14,IW)
                  DO II=M3%IJKW(10,IW),M3%IJKW(13,IW)
                     IWW = IWW + 1
                     M3%REAL_SEND_PKG7(LL+1) = REAL(IW,EB)
                     SELECT CASE(M3%IJKW(4,IW))
                        CASE( 1)
                           M3%REAL_SEND_PKG7(LL+2) = -M%FVX(II,JJ,KK)-(M%H(II+1,JJ,KK)-M%H(II,JJ,KK))*M%RDXN(II)
                        CASE(-1)
                           M3%REAL_SEND_PKG7(LL+2) = -M%FVX(II-1,JJ,KK)-(M%H(II,JJ,KK)-M%H(II-1,JJ,KK))*M%RDXN(II-1)
                        CASE( 2)
                           M3%REAL_SEND_PKG7(LL+2) = -M%FVY(II,JJ,KK)-(M%H(II,JJ+1,KK)-M%H(II,JJ,KK))*M%RDYN(JJ)
                        CASE(-2)
                           M3%REAL_SEND_PKG7(LL+2) = -M%FVY(II,JJ-1,KK)-(M%H(II,JJ,KK)-M%H(II,JJ-1,KK))*M%RDYN(JJ-1)
                        CASE( 3)
                           M3%REAL_SEND_PKG7(LL+2) = -M%FVZ(II,JJ,KK)-(M%H(II,JJ,KK+1)-M%H(II,JJ,KK))*M%RDZN(KK)
                        CASE(-3)
                           M3%REAL_SEND_PKG7(LL+2) = -M%FVZ(II,JJ,KK-1)-(M%H(II,JJ,KK)-M%H(II,JJ,KK-1))*M%RDZN(KK-1)
                     END SELECT
                     LL = LL+2
                  ENDDO
               ENDDO
            ENDDO
         ENDDO PACK_REAL_SEND_PKG7
         M3%REAL_SEND_PKG7(IWW*2+1) = -999.0_EB
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%REAL_SEND_PKG7(1),IWW*2+1,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
          M2=>MESHES(NOM)%OMESH(NM)
          DO IW=1,M4%NEWC
            IF (M3%IJKW(9,IW)/=NM) CYCLE
            DO KK=M3%IJKW(12,IW),M3%IJKW(15,IW)
               DO JJ=M3%IJKW(11,IW),M3%IJKW(14,IW)
                  DO II=M3%IJKW(10,IW),M3%IJKW(13,IW)
                     SELECT CASE(M3%IJKW(4,IW))
                        CASE( 1)
                           M2%DUDT(II,JJ,KK)   = -M%FVX(II,JJ,KK)-(M%H(II+1,JJ,KK)-M%H(II,JJ,KK))*M%RDXN(II)
                        CASE(-1)
                           M2%DUDT(II-1,JJ,KK) = -M%FVX(II-1,JJ,KK)-(M%H(II,JJ,KK)-M%H(II-1,JJ,KK))*M%RDXN(II-1)
                        CASE( 2)
                           M2%DVDT(II,JJ,KK)   = -M%FVY(II,JJ,KK)-(M%H(II,JJ+1,KK)-M%H(II,JJ,KK))*M%RDYN(JJ)
                        CASE(-2)
                           M2%DVDT(II,JJ-1,KK) = -M%FVY(II,JJ-1,KK)-(M%H(II,JJ,KK)-M%H(II,JJ-1,KK))*M%RDYN(JJ-1)
                        CASE( 3)
                           M2%DWDT(II,JJ,KK)   = -M%FVZ(II,JJ,KK)-(M%H(II,JJ,KK+1)-M%H(II,JJ,KK))*M%RDZN(KK)
                        CASE(-3)
                           M2%DWDT(II,JJ,KK-1) = -M%FVZ(II,JJ,KK-1)-(M%H(II,JJ,KK)-M%H(II,JJ,KK-1))*M%RDZN(KK-1)
                     END SELECT
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
      ENDIF
   ENDIF


   ! Send pressure information at the end of the PREDICTOR stage of the time step
 
   IF (CODE==3 .AND. NIC(NOM,NM)>0) THEN
      IF (RNODE/=SNODE) THEN
         LL = 0
         IWW = 0
         PACK_REAL_SEND_PKG2: DO IW=1,M4%NEWC
            IF (M3%IJKW(9,IW)/=NM) CYCLE PACK_REAL_SEND_PKG2
            DO KK=M3%IJKW(12,IW),M3%IJKW(15,IW)
               DO JJ=M3%IJKW(11,IW),M3%IJKW(14,IW)
                  DO II=M3%IJKW(10,IW),M3%IJKW(13,IW)
                     IWW = IWW + 1
                     M3%REAL_SEND_PKG2(LL+1) = REAL(IW,EB)
                     M3%REAL_SEND_PKG2(LL+2) = M%H(II,JJ,KK)
                     M3%REAL_SEND_PKG2(LL+3) = M%US(II,JJ,KK)
                     M3%REAL_SEND_PKG2(LL+4) = M%VS(II,JJ,KK)
                     M3%REAL_SEND_PKG2(LL+5) = M%WS(II,JJ,KK)
                     M3%REAL_SEND_PKG2(LL+6) = M%US(II-1,JJ,KK)
                     M3%REAL_SEND_PKG2(LL+7) = M%VS(II,JJ-1,KK)
                     M3%REAL_SEND_PKG2(LL+8) = M%WS(II,JJ,KK-1)
                     LL = LL+8
                  ENDDO
               ENDDO
            ENDDO
         ENDDO PACK_REAL_SEND_PKG2
         M3%REAL_SEND_PKG2(IWW*8+1) = -999.0_EB
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%REAL_SEND_PKG2(1),IWW*8+1,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         M2=>MESHES(NOM)%OMESH(NM)
         M2%H(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)  = M%H(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%US(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%US(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%VS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%VS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%WS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%WS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
      ENDIF
   ENDIF

   ! Exchange density and mass fraction following CORRECTOR update
 
   IF (CODE==4 .AND. NIC(NOM,NM)>0) THEN
      IF (RNODE/=SNODE) THEN
         LL = 0
         IWW = 0
         PACK_REAL_SEND_PKG3: DO IW=1,M4%NEWC
            IF (M3%IJKW(9,IW)/=NM) CYCLE PACK_REAL_SEND_PKG3
            DO KK=M3%IJKW(12,IW),M3%IJKW(15,IW)
               DO JJ=M3%IJKW(11,IW),M3%IJKW(14,IW)
                  DO II=M3%IJKW(10,IW),M3%IJKW(13,IW)
                     IWW = IWW + 1
                     M3%REAL_SEND_PKG3(LL+1) = REAL(IW,EB)
                     M3%REAL_SEND_PKG3(LL+2) = M%RHO(II,JJ,KK)
                     M3%REAL_SEND_PKG3(LL+3) = M%MU(II,JJ,KK)
                     M3%REAL_SEND_PKG3(LL+4) = M%DS(II,JJ,KK)
                     IF (N_SPECIES>0) M3%REAL_SEND_PKG3(LL+5:LL+4+N_SPECIES) = M%YY(II,JJ,KK,1:N_SPECIES)
                     LL = LL+4+N_SPECIES
                  ENDDO
                ENDDO
             ENDDO
         ENDDO PACK_REAL_SEND_PKG3
         M3%REAL_SEND_PKG3(IWW*(4+N_SPECIES)+1) = -999.0_EB
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%REAL_SEND_PKG3(1),IWW*(4+N_SPECIES)+1,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         M2=>MESHES(NOM)%OMESH(NM)
         M2%RHO(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%RHO(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%MU(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)  = M%MU(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%DS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)  = M%DS(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         IF (N_SPECIES>0) M2%YY(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,1:N_SPECIES)= M%YY(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX,1:N_SPECIES)
      ENDIF
   ENDIF

   ! Exchange BOUNDARY_TYPE following the CORRECTOR stage of the time step

   IF (CODE==0 .OR. CODE==6) THEN
      IF (RNODE/=SNODE) THEN
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M%BOUNDARY_TYPE(0),M%NEWC+1,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         M2=>MESHES(NOM)%OMESH(NM)
         M2%BOUNDARY_TYPE(0:M%NEWC) = M%BOUNDARY_TYPE(0:M%NEWC)
      ENDIF
   ENDIF

   ! Exchange pressure and velocities following CORRECTOR stage of time step
 
   IF (CODE==6 .AND. NIC(NOM,NM)>0) THEN
      IF (RNODE/=SNODE) THEN
         LL = 0
         IWW = 0
         PACK_REAL_SEND_PKG4: DO IW=1,M4%NEWC
            IF (M3%IJKW(9,IW)/=NM) CYCLE PACK_REAL_SEND_PKG4
            DO KK=M3%IJKW(12,IW),M3%IJKW(15,IW)
               DO JJ=M3%IJKW(11,IW),M3%IJKW(14,IW)
                  DO II=M3%IJKW(10,IW),M3%IJKW(13,IW)
                     IWW = IWW + 1
                     M3%REAL_SEND_PKG4(LL+1) = REAL(IW,EB)
                     M3%REAL_SEND_PKG4(LL+2) = M%H(II,JJ,KK)
                     M3%REAL_SEND_PKG4(LL+3) = M%U(II,JJ,KK)
                     M3%REAL_SEND_PKG4(LL+4) = M%V(II,JJ,KK)
                     M3%REAL_SEND_PKG4(LL+5) = M%W(II,JJ,KK)
                     M3%REAL_SEND_PKG4(LL+6) = M%U(II-1,JJ,KK)
                     M3%REAL_SEND_PKG4(LL+7) = M%V(II,JJ-1,KK)
                     M3%REAL_SEND_PKG4(LL+8) = M%W(II,JJ,KK-1)
                     LL = LL+8
                  ENDDO
                ENDDO
             ENDDO
         ENDDO PACK_REAL_SEND_PKG4
         M3%REAL_SEND_PKG4(IWW*8+1) = -999.0_EB
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%REAL_SEND_PKG4(1),IWW*8+1,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         M2=>MESHES(NOM)%OMESH(NM)
         M2%H(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%H(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%U(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%U(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%V(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%V(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
         M2%W(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%W(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
      ENDIF
   ENDIF

   ! Send out radiation info

   SEND_RADIATION: IF ( CODE==6 .AND. NIC(NOM,NM)>0 .AND. EXCHANGE_RADIATION) THEN
      IF (RNODE/=SNODE) THEN
         IWW=0
         LL =0
         PACK_REAL_SEND_PKG5: DO IW=1,M4%NEWC
            IF (M3%IJKW(9,IW)/=NM) CYCLE PACK_REAL_SEND_PKG5
            IWW = IWW+1
            LL  = LL +1
            M3%REAL_SEND_PKG5(LL) = REAL(IW,EB)
            DO NN2=1,NSB
               DO NN1=1,NRA
                  LL = LL + 1
                  M3%REAL_SEND_PKG5(LL) = M3%WALL(IW)%ILW(NN1,NN2)
               ENDDO
            ENDDO
         ENDDO PACK_REAL_SEND_PKG5
         M3%REAL_SEND_PKG5(LL+1) = -999.0_EB
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%REAL_SEND_PKG5(1),(NRA*NSB+1)*IWW+1,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         DO IW=1,M4%NEWC
            IF (M4%IJKW(9,IW)==NM)  &
               M4%WALL(IW)%ILW(1:NRA,1:NSB) = M3%WALL(IW)%ILW(1:NRA,1:NSB)
         ENDDO
      ENDIF
   ENDIF SEND_RADIATION
 
   ! Get Number of Droplet Orphans (droplets that have left other meshes and are waiting to be picked up)
 
   IF (CODE==6 .AND. DROPLET_FILE) THEN
      IF (RNODE/=SNODE) THEN
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%N_DROP_ORPHANS,1,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ELSE
         M2=>MESHES(NOM)%OMESH(NM)
         M2%N_DROP_ADOPT = MIN(M3%N_DROP_ORPHANS,N_DROP_ADOPT_MAX)
         IF (M4%NLP+M2%N_DROP_ADOPT>M4%NLPDIM) CALL RE_ALLOCATE_DROPLETS(1,NOM,0,N_DROP_ADOPT_MAX)
      ENDIF
   ENDIF
 
   ! Sending/Receiving Droplet Buffer Arrays
 
   IF_SEND_DROPLETS: IF (CODE==6 .AND. DROPLET_FILE) THEN 
 
      IF (SNODE/=RNODE) THEN

         NC = 13
         DO N=1,MIN(M3%N_DROP_ORPHANS,N_DROP_ADOPT_MAX)
            M3%REAL_SEND_PKG6((N-1)*NC+1)  = M3%DROPLET(N)%X
            M3%REAL_SEND_PKG6((N-1)*NC+2)  = M3%DROPLET(N)%Y
            M3%REAL_SEND_PKG6((N-1)*NC+3)  = M3%DROPLET(N)%Z
            M3%REAL_SEND_PKG6((N-1)*NC+4)  = M3%DROPLET(N)%TMP
            M3%REAL_SEND_PKG6((N-1)*NC+5)  = M3%DROPLET(N)%U
            M3%REAL_SEND_PKG6((N-1)*NC+6)  = M3%DROPLET(N)%V
            M3%REAL_SEND_PKG6((N-1)*NC+7)  = M3%DROPLET(N)%W
            M3%REAL_SEND_PKG6((N-1)*NC+8)  = M3%DROPLET(N)%R
            M3%REAL_SEND_PKG6((N-1)*NC+9)  = M3%DROPLET(N)%PWT
            M3%REAL_SEND_PKG6((N-1)*NC+10) = M3%DROPLET(N)%A_X
            M3%REAL_SEND_PKG6((N-1)*NC+11) = M3%DROPLET(N)%A_Y
            M3%REAL_SEND_PKG6((N-1)*NC+12) = M3%DROPLET(N)%A_Z
            M3%REAL_SEND_PKG6((N-1)*NC+13) = M3%DROPLET(N)%T
         ENDDO
         BUFFER_SIZE = NC*MIN(M3%N_DROP_ORPHANS,N_DROP_ADOPT_MAX)
         BUFFER_SIZE = MAX(1,BUFFER_SIZE)
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%REAL_SEND_PKG6(1),BUFFER_SIZE,MPI_DOUBLE_PRECISION,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
    
         NC = 3
         DO N=1,MIN(M3%N_DROP_ORPHANS,N_DROP_ADOPT_MAX)
            M3%INTG_SEND_PKG1((N-1)*NC+1) = M3%DROPLET(N)%IOR
            M3%INTG_SEND_PKG1((N-1)*NC+2) = M3%DROPLET(N)%CLASS
            M3%INTG_SEND_PKG1((N-1)*NC+3) = M3%DROPLET(N)%TAG
         ENDDO
         BUFFER_SIZE = NC*MIN(M3%N_DROP_ORPHANS,N_DROP_ADOPT_MAX)
         BUFFER_SIZE = MAX(1,BUFFER_SIZE)
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%INTG_SEND_PKG1(1),BUFFER_SIZE,MPI_INTEGER,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
    
         NC = 2
         DO N=1,MIN(M3%N_DROP_ORPHANS,N_DROP_ADOPT_MAX)
            M3%LOGI_SEND_PKG1((N-1)*NC+1) = M3%DROPLET(N)%SHOW
            M3%LOGI_SEND_PKG1((N-1)*NC+2) = M3%DROPLET(N)%SPLAT
         ENDDO
         BUFFER_SIZE = NC*MIN(M3%N_DROP_ORPHANS,N_DROP_ADOPT_MAX)
         BUFFER_SIZE = MAX(1,BUFFER_SIZE)
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M3%LOGI_SEND_PKG1(1),BUFFER_SIZE,MPI_LOGICAL,SNODE,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)


      ELSE
    
         M2=>MESHES(NOM)%OMESH(NM)
         IF (M2%N_DROP_ADOPT>0) THEN
            M4%DROPLET(M4%NLP+1:M4%NLP+M2%N_DROP_ADOPT)=  M3%DROPLET(1:M2%N_DROP_ADOPT)
            M4%NLP = M4%NLP + M2%N_DROP_ADOPT
         ENDIF 
 
      ENDIF

   ENDIF IF_SEND_DROPLETS
 
ENDDO RECEIVING_MESH_LOOP

ENDDO SENDING_MESH_LOOP

! Send information needed by EVACuation routine

DO NM=1,NMESHES
   IF (CODE==6 .AND. EXCHANGE_EVACUATION .AND. MYID/=EVAC_PROCESS .AND. PROCESS(NM)==MYID .AND. .NOT.EVACUATION_ONLY(NM)) THEN
      M => MESHES(NM)
      TAG = NM*(EVAC_PROCESS+1)*CODE*10
      IWW = (M%IBAR+2)*(M%JBAR+2)*(M%KBAR+2)
      IF (N_SPECIES>0) THEN
         N_REQ=N_REQ+1
         CALL MPI_ISEND(M%YY(0,0,0,1),IWW*N_SPECIES,MPI_DOUBLE_PRECISION,EVAC_PROCESS,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      ENDIF
      N_REQ=N_REQ+1
      CALL MPI_ISEND(M%RHO(0,0,0),IWW,MPI_DOUBLE_PRECISION,EVAC_PROCESS,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_ISEND(M%RSUM(0,0,0),IWW,MPI_DOUBLE_PRECISION,EVAC_PROCESS,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_ISEND(M%TMP(0,0,0),IWW,MPI_DOUBLE_PRECISION,EVAC_PROCESS,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_ISEND(M%UII(0,0,0),IWW,MPI_DOUBLE_PRECISION,EVAC_PROCESS,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      N_REQ=N_REQ+1
      CALL MPI_ISEND(M%CELL_INDEX(0,0,0),IWW,MPI_INTEGER,EVAC_PROCESS,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
      IWW = MAXVAL(M%CELL_INDEX)
      N_REQ=N_REQ+1
      CALL MPI_ISEND(M%SOLID(0),IWW,MPI_LOGICAL,EVAC_PROCESS,TAG,MPI_COMM_WORLD,REQ(N_REQ),IERR)
   ENDIF
ENDDO


! Information from Mesh NM is RECV'ed by Mesh NOM.  NOM is the receiver, NM is the sender.

CALL MPI_WAITALL(N_REQ,REQ(1:N_REQ),ARRAY_OF_STATUSES,IERR)
 
SEND_MESH_LOOP: DO NOM=1,NMESHES
IF (PROCESS(NOM)/=MYID) CYCLE SEND_MESH_LOOP
IF (EVACUATION_ONLY(NOM)) CYCLE SEND_MESH_LOOP ! Issue 257 bug fix
 
RECV_MESH_LOOP: DO NM=1,NMESHES
 
   SNODE = PROCESS(NOM)
   RNODE = PROCESS(NM)
   IF (EVACUATION_ONLY(NM)) CYCLE RECV_MESH_LOOP ! Issue 257 bug fix
   IF (NIC(NOM,NM)==0 .AND. NIC(NM,NOM)==0) CYCLE RECV_MESH_LOOP
 
   IF (CODE>0) THEN
      IF (.NOT.ACTIVE_MESH(NM) .OR. .NOT.ACTIVE_MESH(NOM))  CYCLE RECV_MESH_LOOP
   ENDIF
 
   IF (DEBUG) THEN
      WRITE(LU_ERR,'(I3,A,I3,A,I2,A,I7)') NOM,' receiving data from ',NM, ' code=',CODE,' tag=',TAGS(NM,NOM,CODE)
   ENDIF
 
   M =>MESHES(NM)
   M2=>MESHES(NOM)%OMESH(NM)
   M4=>MESHES(NOM)
 
   TAG = TAGS(NM,NOM,CODE)

   ! Receive information before the time stepping starts needed for radiation exchange
 
   IF (CODE==0 .AND. RADIATION) THEN
      DO IW=1,M%NEWC
         IF (M2%IJKW(9,IW)==NOM) THEN
            ALLOCATE(M2%WALL(IW)%ILW(NRA,NSB))
            M2%WALL(IW)%ILW = SIGMA*TMPA4*RPI
         ENDIF
      ENDDO
   ENDIF 
 
   ! Unpack densities and species mass fractions following PREDICTOR exchange

   IF (CODE==1 .AND. NIC(NOM,NM)>0 .AND. RNODE/=SNODE) THEN
      LL = 0
      UNPACK_REAL_RECV_PKG1: DO 
         IW = NINT(M2%REAL_RECV_PKG1(LL+1))
         IF (IW==-999) EXIT UNPACK_REAL_RECV_PKG1
         DO KK=M4%IJKW(12,IW),M4%IJKW(15,IW)
            DO JJ=M4%IJKW(11,IW),M4%IJKW(14,IW)
               DO II=M4%IJKW(10,IW),M4%IJKW(13,IW)
                  M2%RHOS(II,JJ,KK)  = M2%REAL_RECV_PKG1(LL+2)
                  M2%MU(II,JJ,KK)    = M2%REAL_RECV_PKG1(LL+3)
                  M2%D(II,JJ,KK)     = M2%REAL_RECV_PKG1(LL+4)
                  IF (N_SPECIES>0) M2%YYS(II,JJ,KK,1:N_SPECIES)= M2%REAL_RECV_PKG1(LL+5:LL+4+N_SPECIES)
                  LL = LL+4+N_SPECIES
               ENDDO
            ENDDO
         ENDDO
      ENDDO UNPACK_REAL_RECV_PKG1
   ENDIF 

   ! Unpack densities and species mass fractions following PREDICTOR exchange

   IF ((CODE==2 .OR. CODE==5) .AND. NIC(NOM,NM)>0 .AND. RNODE/=SNODE) THEN
      LL = 0
      UNPACK_REAL_RECV_PKG7: DO
         IW = NINT(M2%REAL_RECV_PKG7(LL+1))
         IF (IW==-999) EXIT UNPACK_REAL_RECV_PKG7
         DO KK=M4%IJKW(12,IW),M4%IJKW(15,IW)
            DO JJ=M4%IJKW(11,IW),M4%IJKW(14,IW)
               DO II=M4%IJKW(10,IW),M4%IJKW(13,IW)
                  SELECT CASE(M4%IJKW(4,IW))
                     CASE( 1)
                        M2%DUDT(II,JJ,KK)   = M2%REAL_RECV_PKG7(LL+2)
                     CASE(-1)
                        M2%DUDT(II-1,JJ,KK) = M2%REAL_RECV_PKG7(LL+2)
                     CASE( 2)
                        M2%DVDT(II,JJ,KK)   = M2%REAL_RECV_PKG7(LL+2)
                     CASE(-2)
                        M2%DVDT(II,JJ-1,KK) = M2%REAL_RECV_PKG7(LL+2)
                     CASE( 3)
                        M2%DWDT(II,JJ,KK)   = M2%REAL_RECV_PKG7(LL+2)
                     CASE(-3)
                        M2%DWDT(II,JJ,KK-1) = M2%REAL_RECV_PKG7(LL+2)
                  END SELECT
                  LL = LL+2
               ENDDO
            ENDDO
         ENDDO
      ENDDO UNPACK_REAL_RECV_PKG7
   ENDIF

   ! Unpack pressure following PREDICTOR stage of time step
 
   IF (CODE==3 .AND. NIC(NOM,NM)>0 .AND. RNODE/=SNODE) THEN
      LL = 0
      UNPACK_REAL_RECV_PKG2: DO
         IW = NINT(M2%REAL_RECV_PKG2(LL+1))
         IF (IW==-999) EXIT UNPACK_REAL_RECV_PKG2
         DO KK=M4%IJKW(12,IW),M4%IJKW(15,IW)
            DO JJ=M4%IJKW(11,IW),M4%IJKW(14,IW)
               DO II=M4%IJKW(10,IW),M4%IJKW(13,IW)
                  M2%H(II,JJ,KK)    = M2%REAL_RECV_PKG2(LL+2)
                  M2%US(II,JJ,KK)   = M2%REAL_RECV_PKG2(LL+3)
                  M2%VS(II,JJ,KK)   = M2%REAL_RECV_PKG2(LL+4)
                  M2%WS(II,JJ,KK)   = M2%REAL_RECV_PKG2(LL+5)
                  M2%US(II-1,JJ,KK) = M2%REAL_RECV_PKG2(LL+6)
                  M2%VS(II,JJ-1,KK) = M2%REAL_RECV_PKG2(LL+7)
                  M2%WS(II,JJ,KK-1) = M2%REAL_RECV_PKG2(LL+8)
                  LL = LL+8
               ENDDO
            ENDDO
         ENDDO
      ENDDO UNPACK_REAL_RECV_PKG2
   ENDIF 

   ! Unpack density and species mass fractions following CORRECTOR update

   IF (CODE==4 .AND. NIC(NOM,NM)>0 .AND. RNODE/=SNODE) THEN
      LL = 0
      UNPACK_REAL_RECV_PKG3: DO 
         IW = NINT(M2%REAL_RECV_PKG3(LL+1))
         IF (IW==-999) EXIT UNPACK_REAL_RECV_PKG3
         DO KK=M4%IJKW(12,IW),M4%IJKW(15,IW)
            DO JJ=M4%IJKW(11,IW),M4%IJKW(14,IW)
               DO II=M4%IJKW(10,IW),M4%IJKW(13,IW)
                  M2%RHO(II,JJ,KK) = M2%REAL_RECV_PKG3(LL+2)
                  M2%MU(II,JJ,KK)  = M2%REAL_RECV_PKG3(LL+3)
                  M2%DS(II,JJ,KK)  = M2%REAL_RECV_PKG3(LL+4)
                  IF (N_SPECIES>0) M2%YY(II,JJ,KK,1:N_SPECIES)= M2%REAL_RECV_PKG3(LL+5:LL+4+N_SPECIES)
                  LL = LL+4+N_SPECIES
               ENDDO
            ENDDO
         ENDDO
      ENDDO UNPACK_REAL_RECV_PKG3
   ENDIF

   ! Unpack pressure and velocities at the end of the CORRECTOR stage of the time step
 
   IF (CODE==6 .AND. NIC(NOM,NM)>0 .AND. RNODE/=SNODE) THEN
      LL = 0
      UNPACK_REAL_RECV_PKG4: DO
         IW = NINT(M2%REAL_RECV_PKG4(LL+1))
         IF (IW==-999) EXIT UNPACK_REAL_RECV_PKG4
         DO KK=M4%IJKW(12,IW),M4%IJKW(15,IW)
            DO JJ=M4%IJKW(11,IW),M4%IJKW(14,IW)
               DO II=M4%IJKW(10,IW),M4%IJKW(13,IW)
                  M2%H(II,JJ,KK)   = M2%REAL_RECV_PKG4(LL+2)
                  M2%U(II,JJ,KK)   = M2%REAL_RECV_PKG4(LL+3)
                  M2%V(II,JJ,KK)   = M2%REAL_RECV_PKG4(LL+4)
                  M2%W(II,JJ,KK)   = M2%REAL_RECV_PKG4(LL+5)
                  M2%U(II-1,JJ,KK) = M2%REAL_RECV_PKG4(LL+6)
                  M2%V(II,JJ-1,KK) = M2%REAL_RECV_PKG4(LL+7)
                  M2%W(II,JJ,KK-1) = M2%REAL_RECV_PKG4(LL+8)
                  LL = LL+8
               ENDDO
            ENDDO
         ENDDO
      ENDDO UNPACK_REAL_RECV_PKG4
   ENDIF


   ! Unpack radiation information at the end of the CORRECTOR stage of the time step

   RECEIVE_RADIATION: IF ( CODE==6 .AND. NIC(NOM,NM)>0 .AND. EXCHANGE_RADIATION .AND. RNODE/=SNODE) THEN
      LL = 0
      UNPACK_REAL_RECV_PKG5: DO 
         LL  = LL + 1
         IW = NINT(M2%REAL_RECV_PKG5(LL))
         IF (IW==-999) EXIT UNPACK_REAL_RECV_PKG5
         DO NN2=1,NSB
            DO NN1=1,NRA
               LL = LL + 1
               M4%WALL(IW)%ILW(NN1,NN2) = M2%REAL_RECV_PKG5(LL)
            ENDDO
         ENDDO
      ENDDO UNPACK_REAL_RECV_PKG5
   ENDIF RECEIVE_RADIATION
 
   ! Get Number of Droplet Orphans
 
   IF (CODE==6 .AND. DROPLET_FILE .AND. RNODE/=SNODE) THEN
      M2%N_DROP_ADOPT = MIN(M2%N_DROP_ADOPT,N_DROP_ADOPT_MAX)
      IF (M4%NLP+M2%N_DROP_ADOPT>M4%NLPDIM) CALL RE_ALLOCATE_DROPLETS(1,NOM,0,N_DROP_ADOPT_MAX)
   ENDIF
 
   ! Sending/Receiving Droplet Buffer Arrays
 
   IF_RECEIVE_DROPLETS: IF (CODE==6 .AND. DROPLET_FILE .AND. RNODE/=SNODE) THEN 
      IF_DROPLETS_SENT: IF (M2%N_DROP_ADOPT>0) THEN
 
         NC = 13
         DO N=M4%NLP+1,M4%NLP+M2%N_DROP_ADOPT
            NN = N-M4%NLP-1
            M4%DROPLET(N)%X   = M2%REAL_RECV_PKG6((NN)*NC+1) 
            M4%DROPLET(N)%Y   = M2%REAL_RECV_PKG6((NN)*NC+2) 
            M4%DROPLET(N)%Z   = M2%REAL_RECV_PKG6((NN)*NC+3) 
            M4%DROPLET(N)%TMP = M2%REAL_RECV_PKG6((NN)*NC+4) 
            M4%DROPLET(N)%U   = M2%REAL_RECV_PKG6((NN)*NC+5) 
            M4%DROPLET(N)%V   = M2%REAL_RECV_PKG6((NN)*NC+6) 
            M4%DROPLET(N)%W   = M2%REAL_RECV_PKG6((NN)*NC+7) 
            M4%DROPLET(N)%R   = M2%REAL_RECV_PKG6((NN)*NC+8) 
            M4%DROPLET(N)%PWT = M2%REAL_RECV_PKG6((NN)*NC+9) 
            M4%DROPLET(N)%A_X = M2%REAL_RECV_PKG6((NN)*NC+10) 
            M4%DROPLET(N)%A_Y = M2%REAL_RECV_PKG6((NN)*NC+11) 
            M4%DROPLET(N)%A_Z = M2%REAL_RECV_PKG6((NN)*NC+12) 
            M4%DROPLET(N)%T   = M2%REAL_RECV_PKG6((NN)*NC+13) 
         ENDDO
 
         NC = 3
         DO N=M4%NLP+1,M4%NLP+M2%N_DROP_ADOPT
            NN = N-M4%NLP-1
            M4%DROPLET(N)%IOR    = M2%INTG_RECV_PKG1((NN)*NC+1) 
            M4%DROPLET(N)%CLASS  = M2%INTG_RECV_PKG1((NN)*NC+2) 
            M4%DROPLET(N)%TAG    = M2%INTG_RECV_PKG1((NN)*NC+3) 
         ENDDO
 
         NC = 2
         DO N=M4%NLP+1,M4%NLP+M2%N_DROP_ADOPT
            NN = N-M4%NLP-1
            M4%DROPLET(N)%SHOW  = M2%LOGI_RECV_PKG1((NN)*NC+1)  
            M4%DROPLET(N)%SPLAT = M2%LOGI_RECV_PKG1((NN)*NC+2)  
         ENDDO
 
         M4%NLP = M4%NLP + M2%N_DROP_ADOPT
 
      ENDIF IF_DROPLETS_SENT
   ENDIF IF_RECEIVE_DROPLETS
 
ENDDO RECV_MESH_LOOP

IF (MOD(ICYC,3)==0 .AND. TIMING .AND. ACTIVE_MESH(NOM)) THEN
   MPI_WALL_TIME = MPI_WTIME() - MPI_WALL_TIME_START
   WRITE(LU_ERR,'(A,I3,A,I1,A,F12.4)')  ' Mesh ',NOM,' exits  Mesh Exchange ',CODE,' at ', MPI_WALL_TIME
ENDIF

ENDDO SEND_MESH_LOOP
 
TUSED(11,:)=TUSED(11,:) + SECOND() - TNOW
END SUBROUTINE MESH_EXCHANGE
 
 

SUBROUTINE WRITE_STRINGS
 
! Write character strings out to the .smv file
 
INTEGER :: N,NOM,N_STRINGS_DUM
CHARACTER(80), ALLOCATABLE, DIMENSION(:) :: STRING_DUM
 
! All meshes send their STRINGs to node 0
 
DO NM=1,NMESHES
   IF (PROCESS(NM)==MYID .AND. MYID>0) THEN
      CALL MPI_SEND(MESHES(NM)%N_STRINGS,1,MPI_INTEGER,0,1,MPI_COMM_WORLD,IERR)
      IF (MESHES(NM)%N_STRINGS>0) CALL MPI_SEND(MESHES(NM)%STRING(1),MESHES(NM)%N_STRINGS*80,MPI_CHARACTER,0,NM, &
                                                MPI_COMM_WORLD,IERR)
   ENDIF
ENDDO
 
! Node 0 receives the STRINGs and writes them to the .smv file
 
IF (MYID==0) THEN
   DO N=1,MESHES(1)%N_STRINGS
      WRITE(LU_SMV,'(A)') TRIM(MESHES(1)%STRING(N))
   ENDDO
   OTHER_MESH_LOOP: DO NOM=2,NMESHES 
      IF (PROCESS(NOM)>0) THEN
         CALL MPI_RECV(N_STRINGS_DUM,1,MPI_INTEGER,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
         IF (N_STRINGS_DUM>0) THEN
            ALLOCATE(STRING_DUM(N_STRINGS_DUM))
            CALL MPI_RECV(STRING_DUM(1),N_STRINGS_DUM*80,MPI_CHARACTER,PROCESS(NOM),NOM,MPI_COMM_WORLD,STATUS,IERR)
         ENDIF
      ELSE
         N_STRINGS_DUM = MESHES(NOM)%N_STRINGS
         IF (N_STRINGS_DUM>0) THEN
            ALLOCATE(STRING_DUM(N_STRINGS_DUM))
            STRING_DUM(1:N_STRINGS_DUM) = MESHES(NOM)%STRING(1:N_STRINGS_DUM)
         ENDIF
      ENDIF
      DO N=1,N_STRINGS_DUM
         WRITE(LU_SMV,'(A)') TRIM(STRING_DUM(N))
      ENDDO
      IF (ALLOCATED(STRING_DUM)) DEALLOCATE(STRING_DUM)
   ENDDO OTHER_MESH_LOOP
ENDIF
 
! All STRING arrays are zeroed out
 
DO NM=1,NMESHES
   IF (PROCESS(NM)==MYID) MESHES(NM)%N_STRINGS = 0
ENDDO
 
END SUBROUTINE WRITE_STRINGS
 
 

SUBROUTINE EXCHANGE_DIAGNOSTICS
 
INTEGER  :: NOM,NECYC,CNT
REAL(EB) :: T_SUM, TNOW
 
TNOW = SECOND()

DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID) CYCLE
   T_SUM = 0.
   SUM_LOOP: DO I=2,N_TIMERS_DIM
      T_SUM = T_SUM + TUSED(I,NM)
   ENDDO SUM_LOOP
ENDDO
 
DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID) CYCLE
   NECYC          = MAX(1,NTCYC(NM)-NCYC(NM))
   T_PER_STEP(NM) = (T_SUM-T_ACCUM(NM))/REAL(NECYC,EB)
   T_ACCUM(NM)    = T_SUM
   NCYC(NM)       = NTCYC(NM)
ENDDO
 
DISP = DISPLS(MYID)+1
CNT  = COUNTS(MYID)
CALL MPI_GATHERV(T(DISP),         CNT,MPI_DOUBLE_PRECISION,T,         COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
CALL MPI_GATHERV(T_ACCUM(DISP),   CNT,MPI_DOUBLE_PRECISION,T_ACCUM,   COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
CALL MPI_GATHERV(T_PER_STEP(DISP),CNT,MPI_DOUBLE_PRECISION,T_PER_STEP,COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
CALL MPI_GATHERV(NTCYC(DISP),     CNT,MPI_INTEGER,         NTCYC,     COUNTS,DISPLS,MPI_INTEGER,         0,MPI_COMM_WORLD,IERR)
CALL MPI_GATHERV(HRR(DISP),       CNT,MPI_DOUBLE_PRECISION,HRR,       COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
CALL MPI_GATHERV(RHRR(DISP),      CNT,MPI_DOUBLE_PRECISION,RHRR,      COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
 
! All nodes greater than 0 send various values to node 0

DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID .OR. MYID==0) CYCLE
   CALL MPI_SEND(MESHES(NM)%DT,                      1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%CFL,                     1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%DIVMX,                   1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%DIVMN,                   1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%RESMAX,                  1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%POIS_PTB,                1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%POIS_ERR,                1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%VN,                      1,MPI_DOUBLE_PRECISION, 0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%ICFL,                    1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%JCFL,                    1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%KCFL,                    1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%IMX,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%JMX,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%KMX,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%IMN,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%JMN,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%KMN,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%IRM,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%JRM,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%KRM,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%I_VN,                    1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%J_VN,                    1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%K_VN,                    1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%NLP,                     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
   CALL MPI_SEND(MESHES(NM)%PRESSURE_ITERATIONS,     1,MPI_INTEGER,          0,1,MPI_COMM_WORLD,IERR)
ENDDO

! Node 0 receives various values from all other nodes


DO NOM=1,NMESHES
   IF (PROCESS(NOM)==0 .OR. MYID/=0) CYCLE
   CALL MPI_RECV(MESHES(NOM)%DT,                      1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%CFL,                     1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%DIVMX,                   1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%DIVMN,                   1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%RESMAX,                  1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%POIS_PTB,                1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%POIS_ERR,                1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%VN,                      1,MPI_DOUBLE_PRECISION,PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%ICFL,                    1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%JCFL,                    1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%KCFL,                    1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%IMX,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%JMX,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%KMX,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%IMN,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%JMN,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%KMN,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%IRM,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%JRM,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%KRM,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%I_VN,                    1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%J_VN,                    1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%K_VN,                    1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%NLP,                     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
   CALL MPI_RECV(MESHES(NOM)%PRESSURE_ITERATIONS,     1,MPI_INTEGER,         PROCESS(NOM),1,MPI_COMM_WORLD,STATUS,IERR)
ENDDO
 
TUSED(11,:) = TUSED(11,:) + SECOND() - TNOW
END SUBROUTINE EXCHANGE_DIAGNOSTICS


SUBROUTINE CORRECT_PRESSURE

USE MATH_FUNCTIONS, ONLY : GAUSSJ
REAL(EB) :: A(NCGC,NCGC),B(NCGC),TNOW
INTEGER :: IERROR,NM
LOGICAL :: RECOMPUTE_A_GLOBAL

TNOW = SECOND()

! Determine if any mesh's RECOMPUTE_A=.TRUE. and if yes, make it .TRUE. for all

CALL MPI_ALLREDUCE(RECOMPUTE_A,RECOMPUTE_A_GLOBAL,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,IERR)
RECOMPUTE_A = RECOMPUTE_A_GLOBAL

! Initializations

IF (RECOMPUTE_A) A     = 0._EB
B     = 0._EB

! Set up linear system of equations for coarse pressure solver

MESH_LOOP_2: DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID) CYCLE MESH_LOOP_2
   CALL COMPUTE_A_B(A,B,NM)
ENDDO MESH_LOOP_2

IF (RECOMPUTE_A) THEN
   A = TRANSPOSE(A)
   CALL MPI_GATHERV(A(1,DISPLS_PC(MYID)+1),COUNTS2D_PC(MYID),MPI_DOUBLE_PRECISION,A,COUNTS2D_PC,DISPLS2D_PC, &
                    MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
ENDIF

CALL MPI_GATHERV(B(DISPLS_PC(MYID)+1),COUNTS_PC(MYID),MPI_DOUBLE_PRECISION,B,COUNTS_PC,DISPLS_PC, &
                 MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)

IF (MYID==0) THEN
   IF (RECOMPUTE_A) THEN
      A = TRANSPOSE(A)
      CALL GAUSSJ(A,NCGC,NCGC,B,1,1,IERROR)
      AINV = A               ! Store inverse of A matrix
   ELSE
      B = MATMUL(AINV,B)
   ENDIF
ENDIF

CALL MPI_BCAST(B,NCGC,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)

MESH_LOOP_3: DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID) CYCLE MESH_LOOP_3
   CALL COMPUTE_CORRECTION_PRESSURE(B,NM)
ENDDO MESH_LOOP_3

TUSED(5,:) = TUSED(5,:) + SECOND() - TNOW
END SUBROUTINE CORRECT_PRESSURE



SUBROUTINE COMPUTE_VOLUME_FLOW

REAL(EB), DIMENSION(NMESHES,NMESHES) :: VDOT,VDOT_LOC
REAL(EB) :: ERROR
TYPE (MESH_TYPE), POINTER :: M
INTEGER :: NM,II,JJ,KK,IOR,IW,NOM

VDOT     = 0._EB
VDOT_LOC = 0._EB

DO NM=1,NMESHES
   IF (PROCESS(NM)/=MYID) CYCLE
   M=>MESHES(NM)
   DO IW=1,M%NEWC
      IF (M%BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) THEN
         NOM = M%IJKW(9,IW)
         II  = M%IJKW(1,IW)
         JJ  = M%IJKW(2,IW)
         KK  = M%IJKW(3,IW)
         IOR = M%IJKW(4,IW)
         SELECT CASE(IOR)
            CASE( 1)
               VDOT_LOC(NM,NOM) = VDOT_LOC(NM,NOM) + M%R(II)*M%DY(JJ)*M%DZ(KK)*M%U(0,JJ,KK)
            CASE(-1)
               VDOT_LOC(NM,NOM) = VDOT_LOC(NM,NOM) + M%R(II-1)*M%DY(JJ)*M%DZ(KK)*M%U(M%IBAR,JJ,KK)
            CASE( 2)
               VDOT_LOC(NM,NOM) = VDOT_LOC(NM,NOM) + M%DX(II)*M%DZ(KK)*M%V(II,0,KK)
            CASE(-2)
               VDOT_LOC(NM,NOM) = VDOT_LOC(NM,NOM) + M%DX(II)*M%DZ(KK)*M%V(II,M%JBAR,KK)
            CASE( 3)
               VDOT_LOC(NM,NOM) = VDOT_LOC(NM,NOM) + M%RC(II)*M%DX(II)*M%DY(JJ)*M%W(II,JJ,0)
            CASE(-3)
               VDOT_LOC(NM,NOM) = VDOT_LOC(NM,NOM) + M%RC(II)*M%DX(II)*M%DY(JJ)*M%W(II,JJ,M%KBAR)
         END SELECT
      ENDIF
   ENDDO
ENDDO

CALL MPI_ALLREDUCE(VDOT_LOC(1,1),VDOT(1,1),NMESHES**2,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)

IF (MYID==0) THEN
   DO NM=1,NMESHES
      DO NOM=1,NMESHES
         ERROR = 2._EB*ABS(VDOT(NM,NOM)-VDOT(NOM,NM))/(ABS(VDOT(NM,NOM)+VDOT(NOM,NM))+1.E-10_EB)
         IF (NM<NOM .AND. ERROR>1.E-5_EB) THEN
            WRITE(LU_ERR,'(A,I3,A,I3,A,E12.6)') 'Volume Flow Error, Meshes ',NM,' and ',NOM,' = ',ERROR
         ENDIF
      ENDDO
   ENDDO
ENDIF

END SUBROUTINE COMPUTE_VOLUME_FLOW



SUBROUTINE DUMP_GLOBAL_OUTPUTS(T)

! Dump HRR data to CHID_hrr.csv, MASS data to CHID_mass.csv, DEVICE data to _devc.csv

REAL(EB) :: T,TNOW
INTEGER :: N,CNT
INTEGER :: NM

TNOW = SECOND()

! Dump out HRR info  after first "gathering" data to node 0

DISP = DISPLS(MYID)+1
CNT  = COUNTS(MYID)

IF_DUMP_HRR: IF (T>=HRR_CLOCK) THEN
   CALL MPI_ALLGATHERV(HRR_TIME_INTERVAL(DISP),CNT,MPI_DOUBLE_PRECISION,HRR_TIME_INTERVAL,COUNTS,DISPLS, &
                       MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
   IF (MINVAL(HRR_TIME_INTERVAL,MASK=.NOT.EVACUATION_ONLY)>0._EB) THEN
      CALL MPI_GATHERV(HRR_SUM(DISP), CNT,MPI_DOUBLE_PRECISION,HRR_SUM, COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHERV(RHRR_SUM(DISP),CNT,MPI_DOUBLE_PRECISION,RHRR_SUM,COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHERV(CHRR_SUM(DISP),CNT,MPI_DOUBLE_PRECISION,CHRR_SUM,COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHERV(FHRR_SUM(DISP),CNT,MPI_DOUBLE_PRECISION,FHRR_SUM,COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      CALL MPI_GATHERV(MLR_SUM(DISP), CNT,MPI_DOUBLE_PRECISION,MLR_SUM, COUNTS,DISPLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      IF (MYID==0) CALL DUMP_HRR(T)
      HRR_CLOCK = HRR_CLOCK + DT_HRR
      HRR_SUM   = 0._EB
      RHRR_SUM  = 0._EB
      CHRR_SUM  = 0._EB
      FHRR_SUM  = 0._EB
      MLR_SUM   = 0._EB
      HRR_TIME_INTERVAL = 0._EB
   ENDIF
ENDIF IF_DUMP_HRR

! Spectral energy output

IF (NMESHES==1 .AND. PERIODIC_TEST==2 .AND. T>=SPEC_CLOCK) CALL SPECTRAL_OUTPUT(1)

! Dump out vegetation data
IF (N_TREES_OUT > 0 .AND. T > VEG_CLOCK) THEN
      COUNTS_TREE(:) = COUNTS(:)*N_TREES_OUT*4
      DISPLS_TREE(:) = DISPLS(:)*N_TREES_OUT*4
      CALL MPI_GATHERV(TREE_OUTPUT_DATA(1,1,DISPLS(MYID)+1),COUNTS_TREE(MYID),MPI_DOUBLE_PRECISION, &
                       TREE_OUTPUT_DATA,COUNTS_TREE,DISPLS_TREE,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,IERR)
      IF (MYID==0) CALL DUMP_VEG(T)
      VEG_CLOCK = VEG_CLOCK + DT_VEG
ENDIF

! Dump out Evac info

IF (MYID==EVAC_PROCESS) CALL EVAC_CSV(T)

! Dump out Mass info after first "gathering" data to node 0


IF_DUMP_MASS: IF (T>=MINT_CLOCK) THEN
   CALL MPI_ALLGATHERV(MINT_TIME_INTERVAL(DISP),CNT,MPI_DOUBLE_PRECISION,MINT_TIME_INTERVAL,COUNTS,DISPLS, &
                       MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)
   IF (MINVAL(MINT_TIME_INTERVAL,MASK=.NOT.EVACUATION_ONLY)>0.) THEN
      CALL MPI_GATHERV(MINT_SUM(0,DISP),COUNTS_MASS(MYID),MPI_DOUBLE_PRECISION,MINT_SUM,COUNTS_MASS,DISPLS_MASS, &
                       MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD, IERR)
      IF (MYID==0) CALL DUMP_MASS(T)
      MINT_CLOCK    = MINT_CLOCK + DT_MASS
      MINT_SUM   = 0._EB
      MINT_TIME_INTERVAL = 0._EB
   ENDIF
ENDIF IF_DUMP_MASS

! Exchange DEVICE parameters among meshes and dump out DEVICE info after first "gathering" data to node 0
 
IF (N_DEVC>0) THEN
  
   ! Exchange the CURRENT_STATE of each DEViCe

   STATE_LOC(1:N_DEVC) = .FALSE.  ! _LOC is a temporary array that holds the STATE value for the devices on each node
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      DO N=1,N_DEVC
         IF (DEVICE(N)%MESH==NM) STATE_LOC(N) = DEVICE(N)%CURRENT_STATE 
      ENDDO
   ENDDO
   CALL MPI_ALLREDUCE(STATE_LOC(1),STATE_GLB(1),N_DEVC,MPI_LOGICAL,MPI_LXOR,MPI_COMM_WORLD,IERR)
   DEVICE(1:N_DEVC)%CURRENT_STATE = STATE_GLB(1:N_DEVC)

   ! Exchange the PRIOR_STATE of each DEViCe

   STATE_LOC(1:N_DEVC) = .FALSE.  ! _LOC is a temporary array that holds the STATE value for the devices on each node
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      DO N=1,N_DEVC
         IF (DEVICE(N)%MESH==NM) STATE_LOC(N) = DEVICE(N)%PRIOR_STATE
      ENDDO
   ENDDO
   CALL MPI_ALLREDUCE(STATE_LOC(1),STATE_GLB(1),N_DEVC,MPI_LOGICAL,MPI_LXOR,MPI_COMM_WORLD,IERR)
   DEVICE(1:N_DEVC)%PRIOR_STATE = STATE_GLB(1:N_DEVC)

   ! Exchange the INSTANT_VALUE of each DEViCe

   TC_LOC(1:N_DEVC) = 0._EB 
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      DO N=1,N_DEVC
         IF (DEVICE(N)%MESH==NM) TC_LOC(N) = DEVICE(N)%INSTANT_VALUE
      ENDDO
   ENDDO
   CALL MPI_ALLREDUCE(TC_LOC(1),TC_GLB(1),N_DEVC,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
   DEVICE(1:N_DEVC)%INSTANT_VALUE = TC_GLB(1:N_DEVC)

   ! Exchange the T_CHANGE of each DEViCe

   TC_LOC(1:N_DEVC) = 0._EB
   DO NM=1,NMESHES
      IF (PROCESS(NM)/=MYID) CYCLE
      DO N=1,N_DEVC
         IF (DEVICE(N)%MESH==NM) TC_LOC(N) = DEVICE(N)%T_CHANGE
      ENDDO
   ENDDO
   CALL MPI_ALLREDUCE(TC_LOC(1),TC_GLB(1),N_DEVC,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)
   DEVICE(1:N_DEVC)%T_CHANGE = TC_GLB(1:N_DEVC)

ENDIF

! Exchange information about Devices that is only needed at print-out time

IF_DUMP_DEVC: IF (T>=DEVC_CLOCK .AND. N_DEVC>0) THEN

   ! Exchange the current COUNT of each DEViCe

   TI_LOC(1:N_DEVC) = DEVICE(1:N_DEVC)%TIME_INTERVAL
   CALL MPI_ALLREDUCE(TI_LOC(1),TI_GLB(1),N_DEVC,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IERR)

   ! Get the current VALUEs of all DEViCes into DEVICE(:)%VALUE on node 0

   IF (MINVAL(TI_GLB)>0._EB) THEN
      TC_LOC(1:N_DEVC) = DEVICE(1:N_DEVC)%VALUE 
      CALL MPI_REDUCE(TC_LOC(1),TC_GLB(1),N_DEVC,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,IERR)
      IF (MYID==0) THEN
         DEVICE(1:N_DEVC)%VALUE         = TC_GLB(1:N_DEVC)
         DEVICE(1:N_DEVC)%TIME_INTERVAL = TI_GLB(1:N_DEVC)
         CALL DUMP_DEVICES(T)
      ENDIF
      DEVC_CLOCK = DEVC_CLOCK + DT_DEVC
      DEVICE(1:N_DEVC)%VALUE = 0._EB
      DEVICE(1:N_DEVC)%TIME_INTERVAL = 0._EB
   ENDIF

ENDIF IF_DUMP_DEVC

! Dump CONTROL info. No gathering required as CONTROL is updated on all meshes

IF (T>=CTRL_CLOCK .AND. N_CTRL>0) THEN
   IF (MYID==0) CALL DUMP_CONTROLS(T)
   CTRL_CLOCK = CTRL_CLOCK + DT_CTRL
ENDIF

TUSED(7,:) = TUSED(7,:) + SECOND() - TNOW
END SUBROUTINE DUMP_GLOBAL_OUTPUTS


SUBROUTINE INITIALIZE_EVAC
 
! Initialize evacuation meshes
 
DO NM=1,NMESHES
   IF (MYID==EVAC_PROCESS .AND. .NOT.EVACUATION_ONLY(NM)) THEN
      M=>MESHES(NM)
      !EVAC: SOLID, CELL_INDEX, OBST_INDEX_C, OBSTRUCTION are allocated in READ_OBST for evac process.
      IF (N_SPECIES>0) THEN
         ALLOCATE(M%YY(0:M%IBP1,0:M%JBP1,0:M%KBP1,N_SPECIES),STAT=IZERO)
         CALL ChkMemErr('MAIN','Evac YY',IZERO)
         M%YY=0._EB
      ENDIF
      ALLOCATE(M%RHO(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('MAIN','Evac RHO',IZERO)
      M%RHO=RHOA
      ALLOCATE(M%RSUM(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('MAIN','Evac RSUM',IZERO)
      M%RSUM=RSUM0
      ALLOCATE(M%TMP(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('MAIN','Evac TMP',IZERO)
      M%TMP=TMPA
      ALLOCATE(M%UII(0:M%IBP1,0:M%JBP1,0:M%KBP1),STAT=IZERO)
      CALL ChkMemErr('MAIN','Evac UII',IZERO)
      M%UII=4._EB*SIGMA*TMPA4
   ENDIF
   IF (PROCESS(NM)/=MYID) CYCLE
   IF (EVACUATION_GRID(NM)) PART_CLOCK(NM) = T_EVAC + DT_PART
   IF (MYID/=EVAC_PROCESS) CYCLE
   IF (ANY(EVACUATION_GRID)) CALL INITIALIZE_EVACUATION(NM,MESH_STOP_STATUS(NM))
   IF (EVACUATION_GRID(NM)) CALL DUMP_EVAC(T_EVAC,NM)
ENDDO
IF (ANY(EVACUATION_GRID)) ICYC = -EVAC_TIME_ITERATIONS

END SUBROUTINE INITIALIZE_EVAC

SUBROUTINE INIT_EVAC_DUMPS
 
! Initialize evacuation dumps
 
T_EVAC  = - EVAC_DT_FLOWFIELD*EVAC_TIME_ITERATIONS + T_BEGIN
T_EVAC_SAVE = T_EVAC
IF (ANY(EVACUATION_GRID) .AND. MYID==EVAC_PROCESS) CALL INITIALIZE_EVAC_DUMPS

END SUBROUTINE INIT_EVAC_DUMPS

SUBROUTINE EVAC_CSV(T)
REAL(EB), INTENT(IN) :: T
 
! Dump out Evac info
IF (T>=EVAC_CLOCK .AND. ANY(EVACUATION_GRID)) THEN
   CALL DUMP_EVAC_CSV(T)
   EVAC_CLOCK = EVAC_CLOCK + DT_HRR
ENDIF

END SUBROUTINE EVAC_CSV

SUBROUTINE EVAC_EXCHANGE
IMPLICIT NONE
LOGICAL EXCHANGE_EVACUATION
INTEGER NM
 
! Fire mesh information ==> Evac meshes
 
IF (.NOT.ANY(EVACUATION_GRID)) RETURN
! IF (ANY(EVACUATION_GRID)) CALL EVAC_MESH_EXCHANGE(T_EVAC,T_EVAC_SAVE,I_EVAC,ICYC,EXCHANGE_EVACUATION,2)
IF (MYID/=EVAC_PROCESS) CALL EVAC_MESH_EXCHANGE(T_EVAC,T_EVAC_SAVE,I_EVAC,ICYC,EXCHANGE_EVACUATION,1)
IF (MYID==EVAC_PROCESS) CALL EVAC_MESH_EXCHANGE(T_EVAC,T_EVAC_SAVE,I_EVAC,ICYC,EXCHANGE_EVACUATION,2)

! Update evacuation devices

DO NM=1,NMESHES
   IF (ACTIVE_MESH(NM)) CYCLE
   IF (.NOT.EVACUATION_GRID(NM)) CYCLE
   IF (MYID/=MAX(0,EVAC_PROCESS)) CYCLE
   CALL UPDATE_GLOBAL_OUTPUTS(T(NM),NM)      
ENDDO

END SUBROUTINE EVAC_EXCHANGE

SUBROUTINE EVAC_PRESSURE_LOOP(NM)
 
! Evacuation flow field calculation
 
INTEGER, INTENT(IN) :: NM
INTEGER :: N
 
IF (EVACUATION_ONLY(NM)) THEN
   PRESSURE_ITERATION_LOOP: DO N=1,EVAC_PRESSURE_ITERATIONS
      CALL NO_FLUX
      IF (EVACUATION_ONLY(NM)) MESHES(NM)%FVZ = 0._EB
      CALL PRESSURE_SOLVER(T(NM),NM)
   ENDDO PRESSURE_ITERATION_LOOP
END IF

END SUBROUTINE EVAC_PRESSURE_LOOP


SUBROUTINE EVAC_MAIN_LOOP
 
! Call evacuation routine and adjust time steps for evac meshes
 
REAL(EB) :: T_FIRE, FIRE_DT
 
IF (.NOT.ANY(EVACUATION_GRID)) RETURN
 
IF (ANY(EVACUATION_ONLY).AND.(ICYC <= 0)) THEN
   ACTIVE_MESH = .FALSE.
END IF
EVAC_DT = EVAC_DT_STEADY_STATE
IF (ICYC < 1) EVAC_DT = EVAC_DT_FLOWFIELD
T_FIRE = T_EVAC + EVAC_DT
IF (ICYC > 0) THEN
   IF (.NOT.ALL(EVACUATION_ONLY)) THEN
      T_FIRE = MINVAL(T,MASK= (.NOT.EVACUATION_ONLY).AND.ACTIVE_MESH)
      DO NM=1,NMESHES
         IF (PROCESS(NM)==MYID) DT_NEXT_SYNC(NM) = MESHES(NM)%DT_NEXT
      ENDDO

      CALL MPI_ALLGATHERV(DT_NEXT_SYNC(DISPLS(MYID)+1),COUNTS(MYID),MPI_DOUBLE_PRECISION,DT_NEXT_SYNC,COUNTS,DISPLS, &
           MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)

     FIRE_DT = MINVAL(DT_NEXT_SYNC,MASK= (.NOT.EVACUATION_ONLY).AND.ACTIVE_MESH)
     T_FIRE = T_FIRE + FIRE_DT
   ENDIF
ENDIF
EVAC_TIME_STEP_LOOP: DO WHILE (T_EVAC < T_FIRE)
   T_EVAC = T_EVAC + EVAC_DT
   IF (MYID==EVAC_PROCESS) CALL PREPARE_TO_EVACUATE(ICYC)
   DO NM=1,NMESHES
      IF (EVACUATION_ONLY(NM)) THEN
         ACTIVE_MESH(NM) = .FALSE.
         CHANGE_TIME_STEP(NM) = .FALSE.
         MESHES(NM)%DT     = EVAC_DT
         MESHES(NM)%DT_NEXT = EVAC_DT
         T(NM)  = T_EVAC
         IF (ICYC <= 1 .And. .Not. BTEST(I_EVAC,2) ) THEN
            IF (ICYC <= 0) ACTIVE_MESH(NM) = .TRUE.
            IF (ICYC <= 0) T(NM) = T_EVAC + EVAC_DT_FLOWFIELD*EVAC_TIME_ITERATIONS - EVAC_DT
         ENDIF
         IF (EVACUATION_GRID(NM) ) THEN
            IF (PROCESS(NM)==MYID) CALL EVACUATE_HUMANS(T_EVAC,NM,ICYC)
            IF (T_EVAC >= PART_CLOCK(NM)) THEN
               IF (PROCESS(NM)==MYID) CALL DUMP_EVAC(T_EVAC,NM)
               DO
                  PART_CLOCK(NM) = PART_CLOCK(NM) + DT_PART
                  IF (PART_CLOCK(NM) >= T_EVAC) EXIT
               ENDDO
            ENDIF
         ENDIF
      ENDIF
   ENDDO
   IF (ICYC < 1) EXIT EVAC_TIME_STEP_LOOP
ENDDO EVAC_TIME_STEP_LOOP

END SUBROUTINE EVAC_MAIN_LOOP


SUBROUTINE GET_REVISION_NUMBER(REV_NUMBER,REV_DATE)
USE isodefs, ONLY : GET_REV_smvv
USE POIS, ONLY : GET_REV_pois
USE VEGE, ONLY : GET_REV_vege
USE COMP_FUNCTIONS, ONLY : GET_REV_func
USE MESH_POINTERS, ONLY : GET_REV_mesh
USE RADCALV, ONLY : GET_REV_irad
USE DCDFLIB, ONLY : GET_REV_ieva
INTEGER,INTENT(INOUT) :: REV_NUMBER
CHARACTER(255),INTENT(INOUT) :: REV_DATE
INTEGER :: MODULE_REV
CHARACTER(255) :: MODULE_DATE

CALL GET_REV_cons(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_ctrl(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_devc(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_divg(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_dump(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_evac(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_fire(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_func(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_ieva(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_init(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_irad(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_mass(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_mesh(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_part(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_pois(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_prec(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_pres(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_radi(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_read(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_smvv(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_type(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_vege(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_velo(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF
CALL GET_REV_wall(MODULE_REV,MODULE_DATE)
IF (MODULE_REV > REV_NUMBER) THEN
   REV_NUMBER = MODULE_REV
   WRITE(REV_DATE,'(A)') MODULE_DATE
ENDIF

END SUBROUTINE GET_REVISION_NUMBER

END PROGRAM FDS

Show details Hide details

Change log

r5401 by mcgratta on Dec 28, 2009 Diff

FDS Source: Issue 939. Add iteration loop
in pressure solver to better enforce no-
flux BC.

Go to:

Project members, sign in to write a code review

Older revisions

r5359 by mcgratta on Dec 21, 2009 Diff

FDS Source: Issue 966. Use time
integration for DEVC and HRR output
rather than pure averages.

r5252 by ru...@nist.gov on Dec 07, 2009 Diff

In main, added arrays COUNTS_TREE and
DISPLS_TREE for use in dumping veg
data

In subroutine DUMP_GLOBAL_OUTPUTS
...

r5210 by mcgratta on Dec 03, 2009 Diff

FDS Source: Update version number to
5.4.3.

All revisions of this file

File info

Size: 98508 bytes, 2714 lines

View raw file

File properties

svn:eol-style: LF
svn:keywords: Date Revision Author URL Id