T%C2(T%NOC(1)+1,1) = M%XF-M%XS      CASE(2)         T%C1(T%NOC(2)+1,2) = M%YF-M%YS         T%C2(T%NOC(2)+1,2) = M%YF-M%YS      CASE(3)         T%C1(T%NOC(3)+1,3) = M%ZF-M%ZS         T%C2(T%NOC(3)+1,3) = M%ZF-M%ZS      END SELECTC      DO N=1,T%NOC(IC)+1      T%C3(N,IC) = (T%C2(N,IC)-T%C2(N-1,IC))/(T%C1(N,IC)-T%C1(N-1,IC))      ENDDOC      ENDIF C      ENDDO ICLOOPC      DEALLOCATE(A)      DEALLOCATE(XX)      DEALLOCATE(ND)CC Set up grid stretching arraysC      ALLOCATE(M%R(0:M%IBAR),STAT=IZERO)      CALL ChkMemErr('READ','R',IZERO)      ALLOCATE(M%RC(0:M%IBAR+1),STAT=IZERO)      CALL ChkMemErr('READ','RC',IZERO) ; M%RC = 1.      ALLOCATE(M%RRN(0:M%IBP1),STAT=IZERO)      CALL ChkMemErr('READ','RRN',IZERO) ; M%RRN = 1.      ALLOCATE(M%X(0:M%IBAR),STAT=IZERO)      CALL ChkMemErr('READ','X',IZERO)      ALLOCATE(M%XC(0:M%IBP1),STAT=IZERO)      CALL ChkMemErr('READ','XC',IZERO)      ALLOCATE(M%HX(0:M%IBP1),STAT=IZERO)      CALL ChkMemErr('READ','HX',IZERO)      ALLOCATE(M%DX(0:M%IBP1),STAT=IZERO)      CALL ChkMemErr('READ','DX',IZERO)      ALLOCATE(M%RDX(0:M%IBP1),STAT=IZERO)      CALL ChkMemErr('READ','RDX',IZERO)      ALLOCATE(M%DXN(0:M%IBAR),STAT=IZERO)      CALL ChkMemErr('READ','DXN',IZERO)      ALLOCATE(M%RDXN(0:M%IBAR),STAT=IZERO)      CALL ChkMemErr('READ','RDXN',IZERO)      ALLOCATE(M%Y(0:M%JBAR),STAT=IZERO)      CALL ChkMemErr('READ','Y',IZERO)      ALLOCATE(M%YC(0:M%JBP1),STAT=IZERO)      CALL ChkMemErr('READ','YC',IZERO)      ALLOCATE(M%HY(0:M%JBP1),STAT=IZERO)      CALL ChkMemErr('READ','HY',IZERO)      ALLOCATE(M%DY(0:M%JBP1),STAT=IZERO)      CALL ChkMemErr('READ','DY',IZERO)      ALLOCATE(M%RDY(0:M%JBP1),STAT=IZERO)      CALL ChkMemErr('READ','RDY',IZERO)      ALLOCATE(M%DYN(0:M%JBAR),STAT=IZERO)      CALL ChkMemErr('READ','DYN',IZERO)      ALLOCATE(M%RDYN(0:M%JBAR),STAT=IZERO)      CALL ChkMemErr('READ','RDYN',IZERO)      ALLOCATE(M%Z(0:M%KBAR),STAT=IZERO)      CALL ChkMemErr('READ','Z',IZERO)      ALLOCATE(M%ZC(0:M%KBP1),STAT=IZERO)      CALL ChkMemErr('READ','ZC',IZERO)      ALLOCATE(M%HZ(0:M%KBP1),STAT=IZERO)      CALL ChkMemErr('READ','HZ',IZERO)      ALLOCATE(M%DZ(0:M%KBP1),STAT=IZERO)      CALL ChkMemErr('READ','DZ',IZERO)      ALLOCATE(M%RDZ(0:M%KBP1),STAT=IZERO)      CALL ChkMemErr('READ','RDZ',IZERO)      ALLOCATE(M%DZN(0:M%KBAR),STAT=IZERO)      CALL ChkMemErr('READ','DZN',IZERO)      ALLOCATE(M%RDZN(0:M%KBAR),STAT=IZERO)      CALL ChkMemErr('READ','RDZN',IZERO)CC Define X grid stretching termsC      M%DXMIN = 1000.      DO I=1,M%IBAR      XI    = (REAL(I,EB)-.5)*M%DXI      M%HX(I) = GP(XI,1,NM)      M%DX(I) = M%HX(I)*M%DXI      M%DXMIN = MIN(M%DXMIN,M%DX(I))      IF (M%HX(I).LE.0.) THEN         WRITE(MESSAGE,'(A,I2)')      .       'ERROR: x transformation not monotonic, mesh ',NM         CALL SHUTDOWN(MESSAGE)         ENDIF      M%RDX(I) = 1./M%DX(I)      ENDDOC      M%HX(0)    = M%HX(1)      M%HX(M%IBP1) = M%HX(M%IBAR)      M%DX(0)    = M%DX(1)      M%DX(M%IBP1) = M%DX(M%IBAR)      M%RDX(0)    = 1./M%DX(1)      M%RDX(M%IBP1) = 1./M%DX(M%IBAR)C      DO I=0,M%IBAR      XI     = I*M%DXI      M%X(I) = M%XS + G(XI,1,NM)      IF (CYLINDRICAL) THEN ; M%R(I) = M%X(I)                       ELSE ; M%R(I) = 1. ; ENDIF      M%DXN(I)  = 0.5*(M%DX(I)+M%DX(I+1))      M%RDXN(I) = 1./M%DXN(I)      ENDDO      M%X(0)      = M%XS      M%X(M%IBAR) = M%XFC      DO I=1,M%IBAR      M%XC(I) = 0.5*(M%X(I)+M%X(I-1))      ENDDO      M%XC(0)      = M%XS - 0.5*M%DX(0)      M%XC(M%IBP1) = M%XF + 0.5*M%DX(M%IBP1)C      IF (CYLINDRICAL) THEN        DO I=1,M%IBAR      M%RRN(I) = 2./(M%R(I)+M%R(I-1))      M%RC(I)  = 0.5*(M%R(I)+M%R(I-1))      ENDDO      M%RRN(0)    = M%RRN(1)      M%RRN(M%IBP1) = M%RRN(M%IBAR)      ENDIFCC Define Y grid stretching termsC      M%DYMIN = 1000.      DO J=1,M%JBAR      ETA   = (REAL(J,EB)-.5)*M%DETA      M%HY(J) = GP(ETA,2,NM)      M%DY(J) = M%HY(J)*M%DETA      M%DYMIN = MIN(M%DYMIN,M%DY(J))      IF (M%HY(J).LE.0.) THEN         WRITE(MESSAGE,'(A,I2)')      .       'ERROR: y transformation not monotonic, mesh ',NM         CALL SHUTDOWN(MESSAGE)         ENDIF      M%RDY(J) = 1./M%DY(J)      ENDDOC      M%HY(0)    = M%HY(1)      M%HY(M%JBP1) = M%HY(M%JBAR)      M%DY(0)    = M%DY(1)      M%DY(M%JBP1) = M%DY(M%JBAR)      M%RDY(0)    = 1./M%DY(1)      M%RDY(M%JBP1) = 1./M%DY(M%JBAR)C      DO J=0,M%JBAR      ETA     = J*M%DETA      M%Y(J)    = M%YS + G(ETA,2,NM)      M%DYN(J)  = 0.5*(M%DY(J)+M%DY(J+1))      M%RDYN(J) = 1./M%DYN(J)      ENDDOC      M%Y(0)      = M%YS      M%Y(M%JBAR) = M%YFC      DO J=1,M%JBAR      M%YC(J) = 0.5*(M%Y(J)+M%Y(J-1))      ENDDO      M%YC(0)      = M%YS - 0.5*M%DY(0)      M%YC(M%JBP1) = M%YF + 0.5*M%DY(M%JBP1)CC Define Z grid stretching termsC      M%DZMIN = 1000.      DO K=1,M%KBAR      ZETA  = (REAL(K,EB)-.5)*M%DZETA      M%HZ(K) = GP(ZETA,3,NM)      M%DZ(K) = M%HZ(K)*M%DZETA      M%DZMIN = MIN(M%DZMIN,M%DZ(K))      IF (M%HZ(K).LE.0.) THEN         WRITE(MESSAGE,'(A,I2)')      .       'ERROR: z transformation not monotonic, mesh ',NM         CALL SHUTDOWN(MESSAGE)         ENDIF      M%RDZ(K) = 1./M%DZ(K)      ENDDOC      M%HZ(0)    = M%HZ(1)      M%HZ(M%KBP1) = M%HZ(M%KBAR)      M%DZ(0)    = M%DZ(1)      M%DZ(M%KBP1) = M%DZ(M%KBAR)      M%RDZ(0)    = 1./M%DZ(1)      M%RDZ(M%KBP1) = 1./M%DZ(M%KBAR)C      DO K=0,M%KBAR      ZETA      = K*M%DZETA      M%Z(K)    = M%ZS + G(ZETA,3,NM)      M%DZN(K)  = 0.5*(M%DZ(K)+M%DZ(K+1))      M%RDZN(K) = 1./M%DZN(K)      ENDDOC      M%Z(0)      = M%ZS      M%Z(M%KBAR) = M%ZFC      DO K=1,M%KBAR      M%ZC(K) = 0.5*(M%Z(K)+M%Z(K-1))      ENDDO      M%ZC(0)      = M%ZS - 0.5*M%DZ(0)      M%ZC(M%KBP1) = M%ZF + 0.5*M%DZ(M%KBP1)CC Set up arrays that will return coordinate positionsC      NIPX   = 100*M%IBAR      NIPY   = 100*M%JBAR      NIPZ   = 100*M%KBAR      NIPXS  = NINT(NIPX*M%DX(0)/(M%XF-M%XS))      NIPXF  = NINT(NIPX*M%DX(M%IBP1)/(M%XF-M%XS))      NIPYS  = NINT(NIPY*M%DY(0)/(M%YF-M%YS))      NIPYF  = NINT(NIPY*M%DY(M%JBP1)/(M%YF-M%YS))      NIPZS  = NINT(NIPZ*M%DZ(0)/(M%ZF-M%ZS))      NIPZF  = NINT(NIPZ*M%DZ(M%KBP1)/(M%ZF-M%ZS))      M%RDXINT = REAL(NIPX,EB)/(M%XF-M%XS)      M%RDYINT = REAL(NIPY,EB)/(M%YF-M%YS)      M%RDZINT = REAL(NIPZ,EB)/(M%ZF-M%ZS)C      ALLOCATE(M%CELLSI(-NIPXS:NIPX+NIPXF),STAT=IZERO)      CALL ChkMemErr('READ','CELLSI',IZERO)      ALLOCATE(M%CELLSJ(-NIPYS:NIPY+NIPYF),STAT=IZERO)      CALL ChkMemErr('READ','CELLSJ',IZERO)      ALLOCATE(M%CELLSK(-NIPZS:NIPZ+NIPZF),STAT=IZERO)      CALL ChkMemErr('READ','CELLSK',IZERO)C      DO I=-NIPXS,NIPX+NIPXF      M%CELLSI(I) = GINV(REAL(I,EB)/M%RDXINT,1,NM)*M%RDXI      M%CELLSI(I) = MAX(M%CELLSI(I),-0.9_EB)      M%CELLSI(I) = MIN(M%CELLSI(I),REAL(M%IBAR)+0.9_EB)      ENDDO      DO J=-NIPYS,NIPY+NIPYF      M%CELLSJ(J) = GINV(REAL(J,EB)/M%RDYINT,2,NM)*M%RDETA      M%CELLSJ(J) = MAX(M%CELLSJ(J),-0.9_EB)      M%CELLSJ(J) = MIN(M%CELLSJ(J),REAL(M%JBAR)+0.9_EB)      ENDDO      DO K=-NIPZS,NIPZ+NIPZF      M%CELLSK(K) = GINV(REAL(K,EB)/M%RDZINT,3,NM)*M%RDZETA      M%CELLSK(K) = MAX(M%CELLSK(K),-0.9_EB)      M%CELLSK(K) = MIN(M%CELLSK(K),REAL(M%KBAR)+0.9_EB)      ENDDOC      ENDDO MESH_LOOPCC      CONTAINSC      INTEGER FUNCTION IFAC(II,N)      INTEGER II,N      IFAC = 1      DO I=II-N+1,II      IFAC = IFAC*I      ENDDO      END FUNCTION IFAC      END SUBROUTINE READ_TRANCC      SUBROUTINE READ_TIMEC      REAL(EB) :: DT,VEL_CHAR      INTEGER :: NM      NAMELIST /TIME/ DT,TWFIN,FYI,WALL_INCREMENT,SYNCHRONIZE,     .                EVAC_DT_FLOWFIELD,EVAC_DT_STEADY_STATE      TYPE (MESH_TYPE), POINTER :: MC      DT             =-1.      TWFIN          = 1.0      WALL_INCREMENT = 2      SET_UP         = .FALSE.      SYNCHRONIZE    = .FALSE.      EVAC_DT_FLOWFIELD = 0.01      EVAC_DT_STEADY_STATE = 0.05C      TIME_LOOP: DO      CALL CHECKREAD('TIME',LU5,IOS) ; IF (IOS.EQ.1) EXIT TIME_LOOP      READ(LU5,TIME,END=21,ERR=22,IOSTAT=IOS)   22 IF (IOS.GT.0) CALL SHUTDOWN('ERROR: Problem with TIME line')      ENDDO TIME_LOOP   21 REWIND(LU5)C      IF (TWFIN.LE.0.) SET_UP = .TRUE.C      IF (SYNCHRONIZE) SYNC_TIME_STEP = .TRUE.      IF (ANY(SYNC_TIME_STEP)) SYNCHRONIZE = .TRUE.C      MESH_LOOP: DO NM=1,NMESHES      M=>MESH(NM)      IF (DT.GT.0.) THEN         M%DT = DT      ELSE         VEL_CHAR = 0.2*SQRT(10.*(M%ZF-M%ZS))         M%DT     = (M%DXMIN*M%DYMIN*M%DZMIN)**(1./3.)/VEL_CHAR      ENDIF      IF (EVACUATION_ONLY(NM)) THEN         SYNC_TIME_STEP(NM) = .FALSE.         M%DT = EVAC_DT_FLOWFIELD         ENDIF      ENDDO MESH_LOOPC      END SUBROUTINE READ_TIMECC      SUBROUTINE READ_MISCC      REAL(EB) MAX_OVER_PRESSURE,DTSAM_PART,GAUGE_TEMPERATURE      NAMELIST /MISC/ PR,SC,TMPA,TMPO,GVEC,RF,NFRAMES,FYI,     .                CSMAG,RAMP_G,BAROCLINIC,     .                DT0DZ,ISOTHERMAL,INCOMPRESSIBLE,     .                PINF,DATABASE,SURF_DEFAULT,EVAC_SURF_DEFAULT,     .                DENSITY,DATABASE_DIRECTORY,RENDER_FILE,     .                C_FORCED,C_VERTICAL,C_HORIZONTAL,RESTART,     .                DTCORE,BACKGROUND_SPECIES,MW,LES,DNS,     .                VISCOSITY,THERMAL_CONDUCTIVITY,NOISE,     .                RADIATION,GAMMA,BNDF_DEFAULT,REACTION,     .                MAX_OVER_PRESSURE,AUTOMATIC_Z,U0,V0,W0,HUMIDITY,     .                POROUS_FLOOR,SUPPRESSION,CHECK_POISSON,     .                TEXTURE_ORIGIN,NSTRATA,SMOKE3D,     .                THICKEN_OBSTRUCTIONS,LEAK_AREA,     .                DROP_VERTICAL_VELOCITY,     .                DROP_HORIZONTAL_VELOCITY,RHO_SOOT,     .                DTSAM_PART,NPPS,DTPAR,DTSPAR,DEBUG,TIMING,     .                GAUGE_TEMPERATURE,CHARACTERISTIC_VELOCITY,     .                FLUSH_FILE_BUFFERS,MAXIMUM_DROPLETS,     .                EVAC_PRESSURE_ITERATIONS,EVAC_TIME_ITERATIONS,     .                TEXTURE_DIRECTORYCC Physical constantsC      R0      = 8314.3    ! Universal Gas Constant (J/kmol/K)      R1      = 1.9862E-3 ! Universal Gas Constant (kcal/mol/K)      TMPA    = 20.       ! Ambient temperature (C)      TMPO    = -10000.       GRAV    = 9.81      ! Acceleration of gravity (m/s**2)      GAMMA   = 1.4       ! Heat capacity ratio for air      PINF    = 101325.   ! Ambient pressure (Pa)      TMPM    = 273.15    ! Melting temperature of water (K)      SIGMA   = 5.67E-8   ! Stefan-Boltzmann constant (W/m**2/K**4)      HUMIDITY= -1.       ! Relative Humidity      RHO_SOOT= 1850.     ! Density of soot particle (kg/m3)CC Empirical constantsC      C_VERTICAL   = 1.31 ! Vertical free convection (Holman, Table 7-2)      C_HORIZONTAL = 1.52 ! Horizontal free convection       C_FORCED     = 0.037 ! Forced convection coefficientCC Miscellanious parametersC      PI      = 4.*ATAN(1.0_EB)      RPI     = 1./PI      TWOPI   = 2.*PI      PIO2    = PI/2.      ONTH    = 1./3.      THFO    = 3./4.      ONSI    = 1./6.      TWTH    = 2./3.      FOTH    = 4./3.      RFPI    = 1./(4.*PI)CC Numerical ParametersC      U0 = 0.      V0 = 0.      W0 = 0.      BACKGROUND_SPECIES = 'AIR'      VISCOSITY = -1.      THERMAL_CONDUCTIVITY = -1.      MU_USER = -1.      K_USER  = -1.      MW      = 0.         DENSITY = -1.      RESTART = .FALSE.      RADIATION      = .TRUE.      SUPPRESSION    = .TRUE.      CHECK_POISSON  = .FALSE.      BAROCLINIC     = .FALSE.      NOISE          = .TRUE.      ISOTHERMAL     = .FALSE.        INCOMPRESSIBLE = .FALSE.      BNDF_DEFAULT   = .TRUE.      LES            = .TRUE.      DNS            = .FALSE.      COMBUSTION_MODEL = 1      AUTOMATIC_Z    = .TRUE.      MAX_OVER_PRESSURE = 1000000.      POROUS_FLOOR = .TRUE.      TEXTURE_ORIGIN(1) = 0.      TEXTURE_ORIGIN(2) = 0.      TEXTURE_ORIGIN(3) = 0.      LEAK_AREA = 0.      GAUGE_TEMPERATURE = -1.      CHARACTERISTIC_VELOCITY = 1.      MAXIMUM_DROPLETS = 500000CC EVACuation parametersC      EVAC_PRESSURE_ITERATIONS = 50      EVAC_TIME_ITERATIONS     = 50CC LES parametersC      CSMAG                = 0.20     ! Smagorinsky constant      PR                   = -1.0     ! Turbulent Prandtl number      SC                   = -1.0     ! Turbulent Schmidt numberCC MiscC      DATABASE             = 'null'      DATABASE_DIRECTORY   = 'null'      TEXTURE_DIRECTORY    = 'null'      RAMP_G               = 'null'      SURF_DEFAULT         = 'INERT'      EVAC_SURF_DEFAULT    = 'INERT'      RENDER_FILE          = 'null'      REACTION             = 'null'      DTCORE               = 1000000.      GVEC(1)              = 0.        ! x-component of gravity       GVEC(2)              = 0.        ! y-component of gravity       GVEC(3)              = -GRAV     ! z-component of gravity       DT0DZ                = -1000000.      NFRAMES              = 1000      ! Number of output data sets      DTSAM_PART           = -1.      NPPS                 = 100000    ! Number Particles Per Set      DTPAR                = 0.05      ! Particle Insertion Interval      DTSPAR               = 0.05      ! Droplet Insertion Interval      RF                   = 1.00      ! Relaxation factor for no-flux      DROPLET_FILE         = .FALSE.      NSTRATA              = 7         ! Number bins for drop dist.      IF (.NOT.TWO_D) SMOKE3D = .TRUE.      IF (     TWO_D) SMOKE3D = .FALSE.