PROGRAM FDS  ! Fire Dynamics Simulator, Main Program, Multiple CPU version.USE PRECISION_PARAMETERSUSE MESH_VARIABLESUSE GLOBAL_CONSTANTSUSE TRANUSE DUMPUSE READ_INPUTUSE INITUSE DIVGUSE PRESUSE MASSUSE PARTUSE VELOUSE RADUSE MEMORY_FUNCTIONSUSE COMP_FUNCTIONS, ONLY : SECOND, WALL_CLOCK_TIME, SHUTDOWNUSE MATH_FUNCTIONS, ONLY : GAUSSJUSE DEVICE_VARIABLESUSE WALL_ROUTINESUSE FIREUSE RADCONSUSE CONTROL_FUNCTIONS!EVAC:USE EVACIMPLICIT NONE! Miscellaneous declarationsCHARACTER(255), PARAMETER :: mainmpiid='$Id: main_mpi.f90 719 2007-10-01 17:09:23Z mcgratta $'CHARACTER(255), PARAMETER :: mainmpirev='$Revision: 719 $'CHARACTER(255), PARAMETER :: mainmpidate='$Date: 2007-10-01 13:09:23 -0400 (Mon, 01 Oct 2007) $'LOGICAL  :: EX,DIAGNOSTICS,EXCHANGE_RADIATION=.TRUE.INTEGER  :: LO10,NM,IZERO,DATE_TIME(8),NN,REVISION_NUMBERCHARACTER(10) :: BIG_BEN(3)CHARACTER(255) :: REVISION_DATEREAL(EB) :: T_MAX,T_MINREAL(EB), ALLOCATABLE, DIMENSION(:) ::  T,TC_GLB,TC_LOC,DT_SYNC, DTNEXT_SYNC,DSUM_ALL,PSUM_ALL,USUM_ALLINTEGER, ALLOCATABLE, DIMENSION(:) ::  MESH_STOP_STATUS,COUNT_GLB,COUNT_LOCLOGICAL, ALLOCATABLE, DIMENSION(:) ::  ACTIVE_MESH,STATE_GLB,STATE_LOCINTEGER NOM,IWW,IWINTEGER, PARAMETER :: N_DROP_ADOPT_MAX=10000TYPE (MESH_TYPE), POINTER :: M,M4TYPE (OMESH_TYPE), POINTER :: M2,M3 ! MPI stuff! INCLUDE '/usr/local/include/mpif.h'  ! Uncomment this line only if the compiler needs help finding mpif.hINCLUDE 'mpif.h'INTEGER :: N,MYID=0,NUMPROCS=1,I,IERR,STATUS(MPI_STATUS_SIZE)INTEGER :: RNODE,BUFFER_SIZE,TAG,PNAMELENINTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: TAGSINTEGER, ALLOCATABLE, DIMENSION(:) :: REQ,PREQINTEGER, ALLOCATABLE, DIMENSION(:,:) :: ARRAY_OF_STATUSESINTEGER, ALLOCATABLE, DIMENSION(:,:) :: ARRAY_OF_STATUSES2INTEGER :: N_REQ,N_PREQCHARACTER(MPI_MAX_PROCESSOR_NAME) PNAME ! 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) WRITE(LU_ERR,'(A,I2,A,I2,A,A)') 'Mesh ',MYID+1,' of ', NUMPROCS,' is alive on ',PNAME(1:PNAMELEN)! Start wall clock timingWALL_CLOCK_START = WALL_CLOCK_TIME() ! Assign a compilation date (All Nodes)WRITE(REVISION_DATE,'(A)') mainmpirev(INDEX(mainmpirev,':')+1:LEN_TRIM(mainmpirev)-2)READ (REVISION_DATE,'(I5)') REVISION_NUMBERWRITE(REVISION_DATE,'(A)') mainmpidateCALL GET_REVISION_NUMBER(REVISION_NUMBER,REVISION_DATE)SVN_REVISION_NUMBER = REVISION_NUMBERWRITE(COMPILE_DATE,'(A)') REVISION_DATE(INDEX(REVISION_DATE,'(')+1:INDEX(REVISION_DATE,')')-1)WRITE(VERSION_STRING,'(A)') '5.0.0'VERSION_NUMBER = 5.0  ! Just use to indicate the major versionPARALLEL       = .TRUE. ! Read input from CHID.data file (All Nodes) CALL READ_DATA(MYID) IF (NMESHES/=NUMPROCS) CALL SHUTDOWN('ERROR: Number of meshes not equal to '// 'number of threads') ! Read input for EVACUATION routine ! IF (ANY(EVACUATION_GRID)) CALL READ_EVAC ! Open and write to Smokeview file (Master Node Only) CALL ASSIGN_FILE_NAMESIF (MYID==0) CALL WRITE_SMOKEVIEW_FILE ! 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)ALLOCATE(T(NMESHES),STAT=IZERO)CALL ChkMemErr('MAIN','T',IZERO)ALLOCATE(DT_SYNC(NMESHES),STAT=IZERO)CALL ChkMemErr('MAIN','DT_SYNC',IZERO)ALLOCATE(DTNEXT_SYNC(NMESHES),STAT=IZERO)CALL ChkMemErr('MAIN','DTNEXT_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 meshesALLOCATE(COUNT_LOC(N_DEVC),STAT=IZERO)CALL ChkMemErr('MAIN','COUNT_LOC',IZERO) ALLOCATE(COUNT_GLB(N_DEVC),STAT=IZERO)CALL ChkMemErr('MAIN','COUNT_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 meshesIF (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)ENDIF! Start the clockT     = T_BEGINMESH_STOP_STATUS = NO_STOP ! Create unique tags for all mesh exchanges ALLOCATE(REQ(NMESHES*NMESHES*10)) REQ = MPI_REQUEST_NULLALLOCATE(PREQ(NMESHES*NMESHES*10)) PREQ = MPI_REQUEST_NULLALLOCATE(ARRAY_OF_STATUSES(MPI_STATUS_SIZE,NMESHES*NMESHES*10))ALLOCATE(ARRAY_OF_STATUSES2(MPI_STATUS_SIZE,NMESHES*NMESHES*10))ALLOCATE(TAGS(NMESHES,NMESHES,0:2))TAG = 0DO NM=1,NMESHES   DO NOM=NM,NMESHES      TAG = TAG+1      TAGS(NM,NOM,0) = TAG      TAGS(NOM,NM,0) = TAG   ENDDOENDDOTAGS(:,:,1) = TAGS(:,:,0) + 1000TAGS(:,:,2) = TAGS(:,:,0) + 2000 ! Initialize global parameters (All Nodes) CALL INITIALIZE_GLOBAL_VARIABLESCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! Initialize radiation (All Nodes) IF (RADIATION) CALL INIT_RADIATIONCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! Allocate and initialize mesh-specific variables DO NM=MYID+1,NMESHES,NUMPROCS   CALL INITIALIZE_MESH_VARIABLES(NM)ENDDOCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! Allocate and initialize mesh variable exchange arrays DO NM=MYID+1,NMESHES,NUMPROCSCALL INITIALIZE_MESH_EXCHANGE(NM)ENDDOCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) CALL MPI_ALLGATHER(I_MIN(1,MYID+1),NMESHES,MPI_INTEGER,I_MIN, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)CALL MPI_ALLGATHER(I_MAX(1,MYID+1),NMESHES,MPI_INTEGER,I_MAX, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)CALL MPI_ALLGATHER(J_MIN(1,MYID+1),NMESHES,MPI_INTEGER,J_MIN, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)CALL MPI_ALLGATHER(J_MAX(1,MYID+1),NMESHES,MPI_INTEGER,J_MAX, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)CALL MPI_ALLGATHER(K_MIN(1,MYID+1),NMESHES,MPI_INTEGER,K_MIN, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)CALL MPI_ALLGATHER(K_MAX(1,MYID+1),NMESHES,MPI_INTEGER,K_MAX, NMESHES,MPI_INTEGER,MPI_COMM_WORLD,IERR)CALL MPI_ALLGATHER(NIC(1,MYID+1),  NMESHES,MPI_INTEGER,NIC,   NMESHES,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=MYID+1,NMESHES,NUMPROCS   CALL DOUBLE_CHECK(NM)ENDDOCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! Potentially read data from a previous calculation  DO NM=MYID+1,NMESHES,NUMPROCS   IF (RESTART) CALL READ_RESTART(T(NM),NM)ENDDOCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! Initialize output files containing global data (Master Node Only) IF (MYID==0) CALL INITIALIZE_GLOBAL_DUMPSCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! Initialize output files that are mesh-specific DO NM=MYID+1,NMESHES,NUMPROCS   CALL INITIALIZE_MESH_DUMPS(NM)   CALL INITIALIZE_DROPLETS(NM)   CALL INITIALIZE_TREES(NM)!     IF (ANY(EVACUATION_GRID)) CALL INITIALIZE_EVACUATION(NM)   CALL POST_RECEIVES(NM,0)ENDDOCALL MPI_BARRIER(MPI_COMM_WORLD, IERR) ! Write out character strings to .smv file CALL WRITE_STRINGS ! Initialize Mesh Exchange Arrays (All Nodes) CALL MESH_EXCHANGE(0)CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)! Make an initial dump of ambient valuesDO NM=MYID+1,NMESHES,NUMPROCS   CALL UPDATE_OUTPUTS(T(NM),NM)         CALL DUMP_MESH_OUTPUTS(T(NM),NM)ENDDOCALL MPI_ALLGATHER(T(MYID+1),1,MPI_DOUBLE_PRECISION,T,1, MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)CALL MPI_BARRIER(MPI_COMM_WORLD, IERR)CALL UPDATE_CONTROLS(T)CALL DUMP_GLOBAL_OUTPUTS(T(1)) ! Check for changes in VENT or OBSTruction control and device status at t=T_BEGINOBST_VENT_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS   CALL OPEN_AND_CLOSE(T(NM),NM)  ENDDO OBST_VENT_LOOP!***********************************************************************************************************************************!                                                           MAIN TIMESTEPPING LOOP!*********************************************************************************************************************************** MAIN_LOOP: DO      ICYC  = ICYC + 1     IF (MOD(ICYC,3)==0 .AND. TIMING) THEN      CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)      WRITE(0,'(A,I2,A,I6,A,I2,A,I3.3)')  ' Thread ',MYID+1,' starts iteration',ICYC,' at ', DATE_TIME(7),'.',DATE_TIME(8)   ENDIF    EXCHANGE_RADIATION = .FALSE.   IF (RADIATION) THEN      IF (MOD(ICYC,ANGLE_INCREMENT*TIME_STEP_INCREMENT)==0) EXCHANGE_RADIATION = .TRUE.   ENDIF    ! Synchronize clocks      CALL MPI_ALLGATHER(T(MYID+1),1,MPI_DOUBLE_PRECISION,T,1, MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)   ! Check for program stops    INQUIRE(FILE=TRIM(CHID)//'.stop',EXIST=EX)   IF (EX) MESH_STOP_STATUS = USER_STOP   CALL MPI_ALLGATHER(MESH_STOP_STATUS(MYID+1),1,MPI_INTEGER,MESH_STOP_STATUS,1,MPI_INTEGER,MPI_COMM_WORLD,IERR)    ! Figure out fastest and slowest meshes    T_MAX = -1000000._EB   T_MIN =  1000000._EB   DO NM=1,NMESHES      T_MIN = MIN(T(NM),T_MIN)      T_MAX = MAX(T(NM),T_MAX)      IF (MESH_STOP_STATUS(NM)/=NO_STOP) GLOBAL_STOP_STATUS = MESH_STOP_STATUS(NM)   ENDDO    ! Determine time step    IF (SYNCHRONIZE) THEN      DT_SYNC(MYID+1) = MESHES(MYID+1)%DTNEXT      CALL MPI_ALLGATHER(DT_SYNC(MYID+1),1,MPI_DOUBLE_PRECISION, DT_SYNC,1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR)      IF (SYNC_TIME_STEP(MYID+1)) THEN         MESHES(MYID+1)%DTNEXT = MINVAL(DT_SYNC,MASK=SYNC_TIME_STEP)         T(MYID+1) = MINVAL(T,MASK=SYNC_TIME_STEP)         ACTIVE_MESH(MYID+1) = .TRUE.      ELSE         ACTIVE_MESH(MYID+1) = .FALSE.         IF (T(MYID+1)+MESHES(MYID+1)%DTNEXT <= T_MAX)  ACTIVE_MESH(MYID+1) = .TRUE.         IF (GLOBAL_STOP_STATUS/=NO_STOP) ACTIVE_MESH(MYID+1) = .TRUE.      ENDIF   ELSE      ACTIVE_MESH = .FALSE.      DO NM=1,NMESHES         IF (T(NM)+MESHES(NM)%DTNEXT <= T_MAX) ACTIVE_MESH(NM) = .TRUE.         IF (GLOBAL_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(ICYC,EB))   IF (MOD(ICYC,10**LO10)==0 .OR. MOD(ICYC,100)==0 .OR. T_MIN>=T_END .OR. GLOBAL_STOP_STATUS/=NO_STOP) DIAGNOSTICS = .TRUE.    ! Give every processor the full ACTIVE_MESH array    CALL MPI_ALLGATHER(ACTIVE_MESH(MYID+1), 1, MPI_LOGICAL,  ACTIVE_MESH,1, 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. !     Do not do EVACuation if past the max iteration criteria!!     IF (ANY(EVACUATION_GRID)) THEN!       EVAC_DT = 1000000.!       DO NM=1,NMESHES!         IF (.NOT.EVACUATION_ONLY(NM))!    .         EVAC_DT = MIN(EVAC_DT,MESHES(NM)%DTNEXT)!       ENDDO!       DO NM=1,NMESHES!         IF (EVACUATION_ONLY(NM)) THEN!           IF (ICYC > EVAC_TIME_ITERATIONS) THEN!             ACTIVE_MESH(NM) = .FALSE.!             EVAC_DT = MIN(EVAC_DT, EVAC_DT_STEADY_STATE)!             MESHES(NM)%DT     = EVAC_DT!             T(NM)           = T(NM) + MESHES(NM)%DT!             IF (EVACUATION_GRID(NM) ) THEN!               CALL EVACUATE_HUMANS(T(NM),NM)!               IF (T(NM)>=PART_CLOCK(NM)) THEN!                 CALL DUMP_EVAC(PART_CLOCK(NM),NM)!                 DO!                   PART_CLOCK(NM) = PART_CLOCK(NM) + WPAR!                   IF (PART_CLOCK(NM)>=T(NM)) EXIT!                 ENDDO!               ENDIF!             ENDIF!           ELSE!             ACTIVE_MESH(NM) = .TRUE.!             EVAC_DT = MIN(EVAC_DT, EVAC_DT_FLOWFIELD)!           ENDIF!         ENDIF!       ENDDO!       DO NM=1,NMESHES!         IF (EVACUATION_ONLY(NM)) MESHES(NM)%DTNEXT = EVAC_DT!       ENDDO!     ENDIF!   PREDICTOR = .TRUE.   CORRECTOR = .FALSE.   ! Diagnostic calls   IF (DEBUG) WRITE(0,*) 'Cycle ',ICYC,' Mesh ',MYID+1, ' starting',ACTIVE_MESH(MYID+1)   IF (MOD(ICYC,3)==0 .AND. TIMING) THEN      CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)      IF (ACTIVE_MESH(MYID+1)) WRITE(0,'(A,I2,A,I2,A,I3.3)') ' Thread ',MYID+1,' is active at ', DATE_TIME(7),'.',DATE_TIME(8)   ENDIF    ! Begin the PREDICTOR step    COMPUTE_FINITE_DIFFERENCES_1: DO NM=MYID+1,NMESHES,NUMPROCS      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_FINITE_DIFFERENCES_1      MESHES(NM)%DT = MESHES(NM)%DTNEXT      NTCYC(NM)   = NTCYC(NM) + 1       CALL INSERT_DROPLETS_AND_PARTICLES(T(NM),NM)      CALL COMPUTE_VELOCITY_FLUX(T(NM),NM)      CALL UPDATE_PARTICLES(T(NM),NM) !     IF (EVACUATION_ONLY(NM)) CALL EVACUATE_HUMANS(T(NM),NM)       IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) CALL MASS_FINITE_DIFFERENCES(NM)    ENDDO COMPUTE_FINITE_DIFFERENCES_1    CHANGE_TIME_STEP_LOOP: DO       COMPUTE_DIVERGENCE_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS         IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_DIVERGENCE_LOOP         IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) THEN            CALL DENSITY(NM)            CALL WALL_BC(T(NM),NM)         ENDIF         CALL DIVERGENCE_PART_1(T(NM),NM)      ENDDO COMPUTE_DIVERGENCE_LOOP      CALL EXCHANGE_DIVERGENCE_INFO      COMPUTE_PRESSURE_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS         IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_PRESSURE_LOOP         CALL DIVERGENCE_PART_2(NM)         CALL PRESSURE_SOLVER(NM)      ENDDO COMPUTE_PRESSURE_LOOP      IF (PRESSURE_CORRECTION) CALL CORRECT_PRESSURE(1) !     IF (ANY(EVACUATION_GRID) .AND. EVACUATION_ONLY(NM)) THEN!        PRESSURE_ITERATION_LOOP: DO N=1,EVAC_PRESSURE_ITERATIONS!           CALL NO_FLUX!           CALL PRESSURE_SOLVER(NM)!        ENDDO PRESSURE_ITERATION_LOOP!     ENDIF       PREDICT_VELOCITY_LOOP: DO NM=MYID+1,NMESHES,NUMPROCS         IF (.NOT.ACTIVE_MESH(NM)) CYCLE PREDICT_VELOCITY_LOOP         CALL VELOCITY_PREDICTOR(T(NM),NM,MESH_STOP_STATUS(NM))      ENDDO PREDICT_VELOCITY_LOOP       IF (SYNCHRONIZE) THEN         NM = MYID+1         CALL MPI_ALLGATHER(CHANGE_TIME_STEP(NM),1,MPI_LOGICAL, CHANGE_TIME_STEP,1,MPI_LOGICAL,MPI_COMM_WORLD,IERR)         CALL MPI_ALLGATHER(MESH_STOP_STATUS(NM),1,MPI_INTEGER,MESH_STOP_STATUS,1, MPI_INTEGER,MPI_COMM_WORLD,IERR)