main.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

PROGRAM FDS  
 
! Fire Dynamics Simulator, Main Program, Single 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 VELO
USE RAD
USE MEMORY_FUNCTIONS
USE COMP_FUNCTIONS, ONLY : SECOND, WALL_CLOCK_TIME, SHUTDOWN
USE DEVICE_VARIABLES
USE WALL_ROUTINES
USE FIRE
USE CONTROL_FUNCTIONS
USE EVAC

IMPLICIT NONE
 
! Miscellaneous declarations
CHARACTER(255), PARAMETER :: mainid='$Id: main.f90 719 2007-10-01 17:09:23Z mcgratta $'
CHARACTER(255), PARAMETER :: mainrev='$Revision: 719 $'
CHARACTER(255), PARAMETER :: maindate='$Date: 2007-10-01 13:09:23 -0400 (Mon, 01 Oct 2007) $'
LOGICAL  :: EX,DIAGNOSTICS
INTEGER  :: MYID=0,LO10,NM,IZERO,REVISION_NUMBER
REAL(EB) :: T_MAX,T_MIN
CHARACTER(255) :: REVISION_DATE
REAL(EB), ALLOCATABLE, DIMENSION(:) :: T,DT_SYNC,DTNEXT_SYNC
INTEGER, ALLOCATABLE, DIMENSION(:) ::  MESH_STOP_STATUS
LOGICAL, ALLOCATABLE, DIMENSION(:) ::  ACTIVE_MESH
INTEGER NOM,IMIN,IMAX,JMIN,JMAX,KMIN,KMAX,IW
INTEGER, PARAMETER :: N_DROP_ADOPT_MAX=10000
TYPE (MESH_TYPE), POINTER :: M,M4
TYPE (OMESH_TYPE), POINTER :: M2,M3

! Start wall clock timing

WALL_CLOCK_START = WALL_CLOCK_TIME()
 
! Assign a compilation date, version number, revision number

WRITE(REVISION_DATE,'(A)') mainrev(INDEX(mainrev,':')+1:LEN_TRIM(mainrev)-2)
READ (REVISION_DATE,'(I5)') REVISION_NUMBER
WRITE(REVISION_DATE,'(A)') maindate
CALL GET_REVISION_NUMBER(REVISION_NUMBER,REVISION_DATE)
SVN_REVISION_NUMBER = REVISION_NUMBER
WRITE(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 used to indicate the major version
SERIAL         = .TRUE.
 
! Read input from CHID.data file (All Nodes)

CALL READ_DATA(MYID)

CALL EVAC_READ_DATA
 
! Open and write to Smokeview file 
 
CALL ASSIGN_FILE_NAMES
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 array (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)
T     = T_BEGIN
MESH_STOP_STATUS = NO_STOP
CALL INITIALIZE_GLOBAL_VARIABLES
IF (RADIATION) CALL INIT_RADIATION
DO NM=1,NMESHES
   CALL INITIALIZE_MESH_VARIABLES(NM)
ENDDO
! Allocate and initialize mesh variable exchange arrays
DO NM=1,NMESHES
   CALL INITIALIZE_MESH_EXCHANGE(NM)
ENDDO
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
   CALL DOUBLE_CHECK(NM)
ENDDO
 
! Potentially read data from a previous calculation 
 
DO NM=1,NMESHES
   IF (RESTART) CALL READ_RESTART(T(NM),NM)
ENDDO
 
! Initialize output files containing global data 
 
CALL INITIALIZE_GLOBAL_DUMPS

CALL INIT_EVAC_DUMPS
 
! Initialize output files that are mesh-specific
 
DO NM=1,NMESHES
   CALL INITIALIZE_MESH_DUMPS(NM)
   CALL INITIALIZE_DROPLETS(NM)
   CALL INITIALIZE_TREES(NM)
   CALL INITIALIZE_EVAC(NM)
ENDDO
 
! Write out character strings to .smv file
 
CALL WRITE_STRINGS

! Initialize Mesh Exchange Arrays (All Nodes)

CALL MESH_EXCHANGE(0)

! Make an initial dump of ambient values

DO NM=1,NMESHES
   CALL UPDATE_OUTPUTS(T(NM),NM)      
   CALL DUMP_MESH_OUTPUTS(T(NM),NM)
ENDDO
CALL UPDATE_CONTROLS(T)
CALL DUMP_GLOBAL_OUTPUTS(T(1))

! Check for changes in VENT or OBSTruction control and device status at t=T_BEGIN

OBST_VENT_LOOP: DO NM=1,NMESHES
   CALL OPEN_AND_CLOSE(T(NM),NM)
ENDDO OBST_VENT_LOOP

! ********************************************************************
!                      MAIN TIMESTEPPING LOOP
! ********************************************************************

MAIN_LOOP: DO  
   ICYC  = ICYC + 1 
   ! Check for program stops

   INQUIRE(FILE=TRIM(CHID)//'.stop',EXIST=EX)
   IF (EX) MESH_STOP_STATUS = USER_STOP
 
   ! 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
 
   IF (SYNCHRONIZE) THEN
      DTNEXT_SYNC(1:NMESHES) = MESHES(1:NMESHES)%DTNEXT
      DO NM=1,NMESHES
         IF (SYNC_TIME_STEP(NM)) THEN
            MESHES(NM)%DTNEXT = MINVAL(DTNEXT_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)%DTNEXT<=T_MAX) ACTIVE_MESH(NM) = .TRUE.
            IF (GLOBAL_STOP_STATUS/=NO_STOP) ACTIVE_MESH(NM) = .TRUE.
         ENDIF
      ENDDO
   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
   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. GLOBAL_STOP_STATUS/=NO_STOP) DIAGNOSTICS = .TRUE.
   
   ! If no meshes are due to be updated, update them all
 
   IF (ALL(.NOT.ACTIVE_MESH)) ACTIVE_MESH = .TRUE.
   CALL EVAC_MAIN_LOOP

!=====================================================================
!  Predictor Step
!=====================================================================

   PREDICTOR = .TRUE.
   CORRECTOR = .FALSE.
   
   COMPUTE_FINITE_DIFFERENCES_1: DO NM=1,NMESHES
      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 (.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=1,NMESHES
         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=1,NMESHES
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_PRESSURE_LOOP
         CALL DIVERGENCE_PART_2(NM)
         CALL PRESSURE_SOLVER(NM)
         CALL EVAC_PRESSURE_LOOP(NM)
      ENDDO COMPUTE_PRESSURE_LOOP
 
      IF (PRESSURE_CORRECTION) CALL CORRECT_PRESSURE
      PREDICT_VELOCITY_LOOP: DO NM=1,NMESHES
         IF (.NOT.ACTIVE_MESH(NM)) CYCLE PREDICT_VELOCITY_LOOP
         CALL VELOCITY_PREDICTOR(T(NM),NM,MESH_STOP_STATUS(NM))
         IF (MESH_STOP_STATUS(NM)==INSTABILITY_STOP) GLOBAL_STOP_STATUS = INSTABILITY_STOP
      ENDDO PREDICT_VELOCITY_LOOP
      IF (GLOBAL_STOP_STATUS/=NO_STOP) EXIT CHANGE_TIME_STEP_LOOP

      IF (SYNCHRONIZE .AND. ANY(CHANGE_TIME_STEP)) THEN
         CHANGE_TIME_STEP = .TRUE.
         DT_SYNC(1:NMESHES) = MESHES(1:NMESHES)%DT
         DTNEXT_SYNC(1:NMESHES) = MESHES(1:NMESHES)%DTNEXT
         DO NM=1,NMESHES
            IF (EVACUATION_ONLY(NM)) CHANGE_TIME_STEP(NM) = .FALSE.
            MESHES(NM)%DTNEXT = MINVAL(DTNEXT_SYNC,MASK=SYNC_TIME_STEP)
            MESHES(NM)%DT     = MINVAL(DT_SYNC,MASK=SYNC_TIME_STEP)
         ENDDO
      ENDIF
 
      IF (.NOT.ANY(CHANGE_TIME_STEP)) EXIT CHANGE_TIME_STEP_LOOP
 
   ENDDO CHANGE_TIME_STEP_LOOP
   CHANGE_TIME_STEP = .FALSE.
   
   DO NM=1,NMESHES
      IF (ACTIVE_MESH(NM)) T(NM) = T(NM) + MESHES(NM)%DT
   ENDDO

   ! Exchange information among meshes

   CALL MESH_EXCHANGE(1)

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  Corrector Step
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   CORRECTOR = .TRUE.
   PREDICTOR = .FALSE.

   COMPUTE_FINITE_DIFFERENCES_2: DO NM=1,NMESHES
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_FINITE_DIFFERENCES_2
      CALL COMPUTE_VELOCITY_FLUX(T(NM),NM)     
      IF (.NOT.ISOTHERMAL .OR. N_SPECIES>0) THEN
         CALL MASS_FINITE_DIFFERENCES(NM)
         CALL DENSITY(NM)
         ! Do combustion, then apply thermal, species and density boundary conditions and solve for radiation
         IF (N_REACTIONS > 0) CALL COMBUSTION (NM)
         CALL WALL_BC(T(NM),NM)
         CALL COMPUTE_RADIATION(NM)
      ENDIF
      CALL UPDATE_PARTICLES(T(NM),NM)
      CALL DIVERGENCE_PART_1(T(NM),NM)
   ENDDO COMPUTE_FINITE_DIFFERENCES_2
     
   CALL EXCHANGE_DIVERGENCE_INFO
   
   COMPUTE_PRESSURE_LOOP_2: DO NM=1,NMESHES
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE COMPUTE_PRESSURE_LOOP_2
      CALL DIVERGENCE_PART_2(NM)
      CALL PRESSURE_SOLVER(NM)
      CALL EVAC_PRESSURE_LOOP(NM)
   ENDDO COMPUTE_PRESSURE_LOOP_2
   
   IF (PRESSURE_CORRECTION) CALL CORRECT_PRESSURE 
     
   CORRECT_VELOCITY_LOOP: DO NM=1,NMESHES
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE CORRECT_VELOCITY_LOOP
      CALL OPEN_AND_CLOSE(T(NM),NM)
      CALL VELOCITY_CORRECTOR(T(NM),NM)
      IF (DIAGNOSTICS) CALL CHECK_DIVERGENCE(NM)
   ENDDO CORRECT_VELOCITY_LOOP

   OUTPUT_LOOP: DO NM=1,NMESHES
      IF (.NOT.ACTIVE_MESH(NM)) CYCLE OUTPUT_LOOP
      CALL UPDATE_OUTPUTS(T(NM),NM)      
      CALL DUMP_MESH_OUTPUTS(T(NM),NM)
   ENDDO OUTPUT_LOOP


!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Exchange information among meshes
   CALL MESH_EXCHANGE(2)
   CALL EVAC_EXCHANGE
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   ! Write character strings out to the .smv file
 
   CALL WRITE_STRINGS

   ! Exchange info for diagnostic print out
 
   IF (DIAGNOSTICS) CALL EXCHANGE_DIAGNOSTICS

   ! Dump global quantities like HRR, MASS, and DEViCes. 
   
   CALL DUMP_GLOBAL_OUTPUTS(MINVAL(T))
   CALL UPDATE_CONTROLS(T)
 
   ! Dump out diagnostics
 
   IF (DIAGNOSTICS) CALL WRITE_DIAGNOSTICS(T)
 
   ! Stop the run
   
   IF (T_MIN>=T_END .OR. GLOBAL_STOP_STATUS/=NO_STOP) EXIT MAIN_LOOP
 
   ! Flush Buffers 
   
   IF (MOD(ICYC,10)==0 .AND. FLUSH_FILE_BUFFERS) THEN
      CALL FLUSH_GLOBAL_BUFFERS
      DO NM=1,NMESHES
         CALL FLUSH_LOCAL_BUFFERS(NM)
      ENDDO
   ENDIF

ENDDO MAIN_LOOP
 
!***********************************************************************
!                          END OF TIMESTEP
!***********************************************************************
 
TUSED(1,1:NMESHES) = SECOND() - TUSED(1,1:NMESHES)
 
CALL TIMINGS
 
SELECT CASE(GLOBAL_STOP_STATUS)
   CASE(NO_STOP)
      CALL SHUTDOWN('STOP: FDS completed successfully')
   CASE(INSTABILITY_STOP)
      CALL SHUTDOWN('STOP: Numerical Instability')
   CASE(USER_STOP)
      CALL SHUTDOWN('STOP: FDS stopped by user')
END SELECT
 
 
CONTAINS

 
SUBROUTINE EXCHANGE_DIVERGENCE_INFO

! Exchange information mesh to mesh used to compute global pressure integrals

INTEGER :: IPZ
REAL(EB) :: DSUM_ALL,PSUM_ALL,USUM_ALL

DO IPZ=1,N_ZONE
   DSUM_ALL = 0._EB
   PSUM_ALL = 0._EB
   USUM_ALL = 0._EB
   DO NM=1,NMESHES
      DSUM_ALL = DSUM_ALL + DSUM(IPZ,NM)
      PSUM_ALL = PSUM_ALL + PSUM(IPZ,NM)
      USUM_ALL = USUM_ALL + USUM(IPZ,NM)
   ENDDO
   DSUM(IPZ,1:NMESHES) = DSUM_ALL
   PSUM(IPZ,1:NMESHES) = PSUM_ALL
   USUM(IPZ,1:NMESHES) = USUM_ALL
ENDDO

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)
 
ALLOCATE(M%OMESH(NMESHES))
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
   IF (NOM==NM) CYCLE OTHER_MESH_LOOP
 
   M2=>MESHES(NOM)
   IMIN=0
   JMIN=0
   KMIN=0
   IMAX=M2%IBP1
   JMAX=M2%JBP1
   KMAX=M2%KBP1
   NIC(NOM,NM) = 0
   FOUND = .FALSE.
   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)
      SELECT CASE(IOR)
         CASE( 1)
            IMIN=MAX(IMIN,M%IJKW(10,IW)-1)
         CASE(-1) 
            IMAX=MIN(IMAX,M%IJKW(10,IW))
         CASE( 2) 
            JMIN=MAX(JMIN,M%IJKW(11,IW)-1)
         CASE(-2) 
            JMAX=MIN(JMAX,M%IJKW(11,IW))
         CASE( 3) 
            KMIN=MAX(KMIN,M%IJKW(12,IW)-1)
         CASE(-3) 
            KMAX=MIN(KMAX,M%IJKW(12,IW))
      END SELECT
   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)% TMP(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%   H(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%   U(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%   V(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)%   W(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)% FVX(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)% FVY(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   ALLOCATE(M%OMESH(NOM)% FVZ(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
   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))
   ENDIF
 
   ! Wall arrays
 
   ALLOCATE(M%OMESH(NOM)%BOUNDARY_TYPE(0:M2%NEWC))
   M%OMESH(NOM)%BOUNDARY_TYPE(0:M2%NEWC) = M2%BOUNDARY_TYPE(0:M2%NEWC)
   ALLOCATE(M%OMESH(NOM)%IJKW(12,M2%NEWC))
   M%OMESH(NOM)%IJKW(1:12,1:M2%NEWC) = M2%IJKW(1:12,1: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 = 1000
      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)
 
OTHER_MESH_LOOP: DO NOM=1,NMESHES
   IF (NOM==NM) CYCLE OTHER_MESH_LOOP
   IF (NIC(NM,NOM)==0 .AND. NIC(NOM,NM)>0) THEN
      M2=>MESHES(NOM)
      ALLOCATE(M%OMESH(NOM)%IJKW(12,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 MESH_EXCHANGE(CODE)

! Exchange Information between Meshes

USE RADCONS, ONLY :NSB,NRA 
REAL(EB) :: TNOW 
INTEGER, INTENT(IN) :: CODE
INTEGER :: NM
 
TNOW = SECOND()
 
MESH_LOOP: DO NM=1,NMESHES
   OTHER_MESH_LOOP: DO NOM=1,NMESHES
 
      IF (CODE==0 .AND. NIC(NOM,NM)<1 .AND. NIC(NM,NOM)>0 .AND. I_MIN(NOM,NM)<0 .AND. RADIATION) THEN
         M =>MESHES(NM)
         M2=>MESHES(NOM)%OMESH(NM)
         DO IW=1,M%NEWC
            IF (M%IJKW(9,IW)==NOM) THEN
               ALLOCATE(M2%WALL(IW)%ILW(NRA,NSB))
               M2%WALL(IW)%ILW = SIGMA*TMPA4*RPI
            ENDIF
         ENDDO