Merge remote-tracking branch 'github/master'

# Conflicts:
#	Manuals/FDS_User_Guide/FDS_User_Guide.tex

Author: drjfloyd

Source/cons.f90

@@ -448,6 +448,7 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: OUTPUT_CHEM_IT=.FALSE.
 LOGICAL :: REAC_SOURCE_CHECK=.FALSE.
 LOGICAL :: COMPUTE_ADIABATIC_FLAME_TEMPERATURE=.FALSE.              !< Report adiabatic flame temperature per REAC in LU_OUTPUT
+LOGICAL :: VARIABLE_CFT=.FALSE.                                     !< Experimental critical flame temp concept
 
 REAL(EB) :: RSUM0                                     !< Initial specific gas constant, \f$ R \sum_i Z_{i,0}/W_i \f$
 
@@ -611,8 +612,7 @@ MODULE GLOBAL_CONSTANTS
 INTEGER :: LEVEL_SET_MODE=0               !< Indicator of the type of level set calculation to be done
 LOGICAL :: LEVEL_SET_COUPLED_FIRE=.TRUE.  !< Indicator for fire and wind level set coupling
 LOGICAL :: LEVEL_SET_COUPLED_WIND=.TRUE.  !< Indicator for fire and wind level set coupling
-LOGICAL :: LEVEL_SET_ELLIPSE=.TRUE.       !< Indicator of Richards elliptical level set formulation
-LOGICAL :: LSET_TAN2
+LOGICAL :: LEVEL_SET_ELLIPSE=.TRUE.       !< Placeholder for future level set spread formulations
 
 ! Parameters for Terrain and Wind simulation needs
 

Source/data.f90

@@ -34,7 +34,7 @@ MODULE OUTPUT_DATA
               INTEGRATED=.FALSE.,SOLID_PHASE=.FALSE.,GAS_PHASE=.TRUE.,INTEGRATED_PARTICLES=.FALSE.,&
               MASS_FRACTION=.FALSE.,SPEC_ID_REQUIRED=.FALSE.,PART_ID_REQUIRED=.FALSE.,INSIDE_SOLID=.FALSE.,&
               MATL_ID_REQUIRED=.FALSE.,TIME_AVERAGED=.TRUE.,DEVC_APPROPRIATE=.TRUE.,&
-              DUCT_ID_REQUIRED=.FALSE.,NODE_ID_REQUIRED=.FALSE.,HVAC=.FALSE.,QUANTITY2_REQUIRED=.FALSE.,&
+              DUCT_ID_REQUIRED=.FALSE.,NODE_ID_REQUIRED=.FALSE.,HVAC=.FALSE.,&
               CELL_L_REQUIRED=.FALSE.,HVAC_SMV=.FALSE.,PROF_APPROPRIATE=.FALSE., ELEM_ID_REQUIRED=.FALSE.
 END TYPE OUTPUT_QUANTITY_TYPE
 

Source/dump.f90

@@ -996,7 +996,7 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_
 IF (SUPPRESSION .AND. .NOT.EXTINCT) THEN
    SELECT CASE(EXTINCT_MOD)
       CASE(EXTINCTION_1); CALL EXTINCT_1(EXTINCT,ZZ_0,TMP_IN)
-      CASE(EXTINCTION_2); CALL EXTINCT_2(EXTINCT,ZZ_0,ZZ_MIXED,TMP_IN)
+      CASE(EXTINCTION_2); CALL EXTINCT_2(EXTINCT,ZZ_0,ZZ_MIXED,TMP_IN,TAU_MIX)
    END SELECT
 ENDIF
 
@@ -1417,10 +1417,10 @@ END SUBROUTINE EXTINCT_1
 !> \param ZZ_IN Array of lumped species mass fractions in the mixed part of the grid cell at the end of the time step
 !> \param TMP_IN Initial temperature of the grid cell
 
-SUBROUTINE EXTINCT_2(EXTINCT,ZZ_0,ZZ_IN,TMP_IN)
+SUBROUTINE EXTINCT_2(EXTINCT,ZZ_0,ZZ_IN,TMP_IN,TAU_MIX)
 
 USE PHYSICAL_FUNCTIONS, ONLY: GET_ENTHALPY
-REAL(EB),INTENT(IN) :: TMP_IN,ZZ_IN(1:N_TRACKED_SPECIES),ZZ_0(1:N_TRACKED_SPECIES)
+REAL(EB),INTENT(IN) :: TMP_IN,ZZ_IN(1:N_TRACKED_SPECIES),ZZ_0(1:N_TRACKED_SPECIES),TAU_MIX
 LOGICAL, INTENT(INOUT) :: EXTINCT
 REAL(EB) :: ZZ_HAT_0(1:N_TRACKED_SPECIES),ZZ_HAT(1:N_TRACKED_SPECIES),H_0,H_CRIT,PHI_TILDE,CFT
 INTEGER :: NS,NR
@@ -1443,7 +1443,11 @@ SUBROUTINE EXTINCT_2(EXTINCT,ZZ_0,ZZ_IN,TMP_IN)
    RETURN
 ENDIF
 
-CFT = SUM_CFT/SUM_ZZ
+IF (VARIABLE_CFT) THEN
+   CFT = 1522.1_EB*(1000._EB*TAU_MIX)**(-0.038_EB)  ! Experimental mechanism, U of Ghent
+ELSE
+   CFT = SUM_CFT/SUM_ZZ
+ENDIF
 
 ! Compute the modified cell equivalence ratio
 

Source/fire.f90

@@ -731,7 +731,7 @@ SUBROUTINE SET_CUTCELLS_3D
 REAL(EB), DIMENSION(MAX_DIM) :: PLNORMAL
 INTEGER,  DIMENSION(MAX_DIM) :: INDX1
 REAL(EB) :: X1PLN, X3RAY
-LOGICAL :: TRI_ONPLANE_ONLY, RAYTRACE_X2_ONLY
+LOGICAL :: TRI_ONPLANE_ONLY, RAYTRACE_X2_ONLY, MESH_BOUNDARY_FACE
 LOGICAL, SAVE :: FIRST_CALL = .TRUE.
 INTEGER :: NCUTFACE_IAXIS, NCUTFACE_JAXIS, NCUTFACE_KAXIS, ICE1, ICF1, NFACE, IERR, &
            NCUTEDGE_IBCC, NCUTEDGE_IBCF
@@ -1401,7 +1401,6 @@ SUBROUTINE SET_CUTCELLS_3D
       ENDIF
    ENDDO
 
-
    IF (ONE_CC_PER_CARTESIAN_CELL) THEN
       ! Here Block all cells that have volume less (or equal) than the first largest cell found.
       DO ICC1=1,MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH
@@ -1540,6 +1539,51 @@ SUBROUTINE SET_CUTCELLS_3D
          K=K+1
       ENDDO
    ENDDO
+
+   ! Block cells trapped against a mesh boundary whose only same-mesh gas opening is a small cut-face.
+   DO ICC1=1,MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH
+      CC=>MESHES(NM)%CUT_CELL(ICC1)
+      DO J=1,CC%NCELL
+         IF(CC%NOADVANCE(J)/=NOT_BLOCKED) CYCLE
+         SUM_FACE=0; SUM_CCELL=0; ICF1=0
+         DO I=1,CC%CCELEM(1,J)
+            IFACE  = CC%CCELEM(I+1,J)
+            SELECT CASE(CC%FACE_LIST(1,IFACE))
+            CASE(CC_FTYPE_CFGAS,CC_FTYPE_RCGAS)
+               SIDE   = CC%FACE_LIST(2,IFACE)
+               X1AXIS = CC%FACE_LIST(3,IFACE)
+               IIF = CC%IJK(IAXIS) + ADDI(SIDE,X1AXIS)
+               JJF = CC%IJK(JAXIS) + ADDJ(SIDE,X1AXIS)
+               KKF = CC%IJK(KAXIS) + ADDK(SIDE,X1AXIS)
+               MESH_BOUNDARY_FACE = .FALSE.
+               SELECT CASE(X1AXIS)
+               CASE(IAXIS); MESH_BOUNDARY_FACE = IIF==0 .OR. IIF==M%IBAR
+               CASE(JAXIS); MESH_BOUNDARY_FACE = JJF==0 .OR. JJF==M%JBAR
+               CASE(KAXIS); MESH_BOUNDARY_FACE = KKF==0 .OR. KKF==M%KBAR
+               END SELECT
+               IF(MESH_BOUNDARY_FACE) CYCLE
+               SELECT CASE(CC%FACE_LIST(1,IFACE))
+               CASE(CC_FTYPE_CFGAS)
+                  SUM_FACE = SUM_FACE+1
+                  ICF1     = CC%FACE_LIST(4,IFACE)
+               CASE(CC_FTYPE_RCGAS)
+                  SUM_CCELL = SUM_CCELL+1
+               END SELECT
+            END SELECT
+         ENDDO
+         IF(SUM_FACE/=1 .OR. SUM_CCELL>0) CYCLE
+         CF => M%CUT_FACE(ICF1)
+         IIF=CF%IJK(IAXIS); JJF=CF%IJK(JAXIS); KKF=CF%IJK(KAXIS)
+         SELECT CASE(CF%IJK(KAXIS+1))
+         CASE(IAXIS); ACRT = DY(JJF)*DZ(KKF)
+         CASE(JAXIS); ACRT = DZ(KKF)*DX(IIF)
+         CASE(KAXIS); ACRT = DX(IIF)*DY(JJF)
+         END SELECT
+         IF(SUM(CF%AREA(1:CF%NFACE))/ACRT>=CCVOL_LINK) CYCLE
+         CC%NOADVANCE(J)=BLOCKED_CAVITY_CELL
+         K=K+1
+      ENDDO
+   ENDDO
    IF (K>0) THEN
       CALL BLOCK_SMALL_UNLINKED_CUTCELLS(NM,SUM_CCELL)
       IF(SUM_CCELL>0) THEN ! Rebuild incidences and cell linking information:

Source/geom.f90

@@ -332,8 +332,7 @@ MODULE MESH_VARIABLES
    INTEGER, ALLOCATABLE, DIMENSION(:,:) :: LS_KLO_TERRAIN,LS_KHI_TERRAIN,K_LS,LS_SURF_INDEX
    INTEGER :: N_TERRAIN_SLCF=0
 
-   REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: FLUX0_LS,FLUX1_LS,PHI_LS,PHI1_LS,ROS_BACKU, &
-                                            ROS_HEAD,ROS_FLANK,WIND_EXP, &
+   REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: FLUX0_LS,FLUX1_LS,PHI_LS,PHI1_LS,ROS_HEAD, &
                                             SR_X_LS,SR_Y_LS,U_LS,V_LS,Z_LS,DZTDX,DZTDY,MAG_ZT, &
                                             PHI_WS,UMF,THETA_ELPS,PHI_S_X,PHI_S_Y,PHI_W,LS_WORK1,LS_WORK2
 
@@ -476,8 +475,7 @@ MODULE MESH_POINTERS
 INTEGER, POINTER, DIMENSION(:,:,:) :: K_AGL_SLICE
 INTEGER, POINTER, DIMENSION(:,:) :: LS_KLO_TERRAIN,LS_KHI_TERRAIN,K_LS,LS_SURF_INDEX
 INTEGER, POINTER :: N_TERRAIN_SLCF
-REAL(EB), POINTER, DIMENSION(:,:) :: FLUX0_LS,FLUX1_LS,PHI_LS,PHI1_LS,ROS_BACKU, &
-                                     ROS_HEAD,ROS_FLANK,WIND_EXP, &
+REAL(EB), POINTER, DIMENSION(:,:) :: FLUX0_LS,FLUX1_LS,PHI_LS,PHI1_LS,ROS_HEAD, &
                                      SR_X_LS,SR_Y_LS,U_LS,V_LS,Z_LS,DZTDX,DZTDY,MAG_ZT, &
                                      PHI_WS,UMF,THETA_ELPS,PHI_S_X,PHI_S_Y,PHI_W,LS_WORK1,LS_WORK2
 
@@ -868,10 +866,7 @@ SUBROUTINE POINT_TO_MESH(NM)
 FLUX1_LS => M%FLUX1_LS
 PHI_LS => M%PHI_LS
 PHI1_LS => M%PHI1_LS
-ROS_BACKU => M%ROS_BACKU
 ROS_HEAD => M%ROS_HEAD
-ROS_FLANK => M%ROS_FLANK
-WIND_EXP => M%WIND_EXP
 SR_X_LS => M%SR_X_LS
 SR_Y_LS => M%SR_Y_LS
 U_LS => M%U_LS

Source/mesh.f90

@@ -4522,7 +4522,7 @@ SUBROUTINE READ_COMB
                 EXTINCTION_MODEL,FINITE_RATE_MIN_TEMP, FIXED_MIX_TIME,FREE_BURN_TEMPERATURE, &
                 FUEL_ID_FOR_AFT,INITIAL_UNMIXED_FRACTION, MAX_CHEMISTRY_SUBSTEPS, MAX_EQUIV_RATIO, MIN_EQUIV_RATIO, &
                 N_FIXED_CHEMISTRY_SUBSTEPS, ODE_MIN_ATOL,ODE_REL_ERROR,ODE_SOLVER,SUPPRESSION,TAU_CHEM, &
-                TAU_FLAME,TURBULENT_FLAME_SPEED,RAMP_ZETA_0,USE_MIXED_ZN_AFT_TMP,ZZ_MIN_GLOBAL
+                TAU_FLAME,TURBULENT_FLAME_SPEED,RAMP_ZETA_0,USE_MIXED_ZN_AFT_TMP,VARIABLE_CFT,ZZ_MIN_GLOBAL
 
 ODE_SOLVER    = 'null'
 ODE_REL_ERROR = -1._EB
@@ -7893,11 +7893,10 @@ SUBROUTINE READ_SURF(QUICK_READ)
 INTEGER :: NPPC,N,IOS,NL,NN,NNN,NNNN,N_LIST,LEAK_PATH(2),DUCT_PATH(2),RGB(3),NR,IL
 INTEGER ::  N_LAYER_CELLS_MAX(MAX_LAYERS),VEG_LSET_FUEL_INDEX,INDEX_LIST(MAX_MATERIALS**2)
 INTEGER :: CHILD_LAYER(MAX_LAYERS,N_MATL),CHILD_SURF(N_MATL)
-REAL(EB) :: VEG_LSET_IGNITE_TIME,VEG_LSET_QCON,VEG_LSET_ROS_HEAD,VEG_LSET_ROS_FLANK,VEG_LSET_ROS_BACK, &
-            VEG_LSET_WIND_EXP,VEG_LSET_BETA,VEG_LSET_HT,VEG_LSET_SIGMA,VEG_LSET_ROS_00, &
+REAL(EB) :: VEG_LSET_IGNITE_TIME,VEG_LSET_QCON,VEG_LSET_BETA,VEG_LSET_HT,VEG_LSET_SIGMA,VEG_LSET_ROS_00, &
             VEG_LSET_M1,VEG_LSET_M10,VEG_LSET_M100,VEG_LSET_MLW,VEG_LSET_MLH,VEG_LSET_SURF_LOAD,VEG_LSET_FIREBASE_TIME,&
             VEG_LSET_CHAR_FRACTION,VEG_LSET_WIND_HEIGHT,VEL_PART,INIT_PER_AREA,TIME_STEP_FACTOR
-LOGICAL :: DEFAULT,VEG_LSET_SPREAD,VEG_LSET_TAN2,TGA_ANALYSIS,VARIABLE_THICKNESS,HT3D,THERM_THICK,VEG_LSET_ROS_FIXED,SKIP_INRAD
+LOGICAL :: DEFAULT,VEG_LSET_SPREAD,TGA_ANALYSIS,VARIABLE_THICKNESS,HT3D,THERM_THICK,VEG_LSET_ROS_FIXED,SKIP_INRAD
 LOGICAL, ALLOCATABLE, DIMENSION(:) :: DUPLICATE
 ! Ember generating variables
 REAL(EB) :: EMBER_GENERATION_HEIGHT(2),EMBER_IGNITION_POWER_MEAN,EMBER_IGNITION_POWER_SIGMA,EMBER_TRACKING_RATIO,EMBER_YIELD
@@ -7933,8 +7932,8 @@ SUBROUTINE READ_SURF(QUICK_READ)
                 TIME_STEP_FACTOR,TMP_BACK,TMP_FRONT,TMP_FRONT_INITIAL,TMP_GAS_BACK,TMP_GAS_FRONT,TMP_INNER,TRANSPARENCY,&
                 VEG_LSET_BETA,VEG_LSET_CHAR_FRACTION,VEG_LSET_FIREBASE_TIME,VEG_LSET_FUEL_INDEX,VEG_LSET_HT,VEG_LSET_IGNITE_TIME,&
                 VEG_LSET_M1,VEG_LSET_M10,VEG_LSET_M100,VEG_LSET_MLW,VEG_LSET_MLH,VEG_LSET_QCON,&
-                VEG_LSET_ROS_00,VEG_LSET_ROS_BACK,VEG_LSET_ROS_FLANK,VEG_LSET_ROS_HEAD,VEG_LSET_ROS_FIXED,VEG_LSET_SIGMA,&
-                VEG_LSET_SURF_LOAD,VEG_LSET_TAN2,VEG_LSET_WIND_EXP,VEG_LSET_WIND_RAMP,VEG_LSET_WIND_HEIGHT,&
+                VEG_LSET_ROS_00,VEG_LSET_ROS_FIXED,VEG_LSET_SIGMA,&
+                VEG_LSET_SURF_LOAD,VEG_LSET_WIND_RAMP,VEG_LSET_WIND_HEIGHT,&
                 VEL,VEL_BULK,VEL_GRAD,VEL_PART,VEL_T,VOLUME_FLOW,WIDTH,XYZ,Z0,Z_0
 
 ! Count the SURF lines in the input file
@@ -8209,27 +8208,20 @@ SUBROUTINE READ_SURF(QUICK_READ)
    ! Level set vegetation fire spread specific
 
    VEG_LSET_SPREAD = .FALSE.
-   IF (VEG_LSET_IGNITE_TIME < 1.E6_EB .OR. VEG_LSET_FUEL_INDEX>0 .OR. &
-       VEG_LSET_ROS_00>0._EB .OR. VEG_LSET_ROS_HEAD>0._EB) VEG_LSET_SPREAD = .TRUE.
+   IF (VEG_LSET_IGNITE_TIME < 1.E6_EB .OR. VEG_LSET_FUEL_INDEX>0 .OR. VEG_LSET_ROS_00>0._EB) VEG_LSET_SPREAD = .TRUE.
    IF (VEG_LSET_SPREAD .AND. LEVEL_SET_MODE==0) THEN
       WRITE(MESSAGE,'(A,A,A)') 'ERROR(305): SURF ',TRIM(ID),' indicates a level set simulation, but LEVEL_SET_MODE not set on MISC.'
       CALL SHUTDOWN(MESSAGE) ; RETURN
    ENDIF
    IF (VEG_LSET_FUEL_INDEX>0 .AND. LEVEL_SET_COUPLED_FIRE) HRRPUA = 1._EB  ! HRRPUA to be set properly later
    IF (VEG_LSET_ROS_00    >0 .AND. LEVEL_SET_COUPLED_FIRE) HRRPUA = 1._EB
-   IF (VEG_LSET_ROS_HEAD  >0 .AND. LEVEL_SET_COUPLED_FIRE) HRRPUA = 1._EB
 
    SF%VEG_LSET_SPREAD       = VEG_LSET_SPREAD
    SF%VEG_LSET_ROS_00       = VEG_LSET_ROS_00       ! no-wind, no-slope RoS (m/s), Rothermel model
-   SF%VEG_LSET_ROS_HEAD     = VEG_LSET_ROS_HEAD     ! head fire rate of spread (m/s), McArthur model
-   SF%VEG_LSET_ROS_FLANK    = VEG_LSET_ROS_FLANK    ! flank fire rate of spread, McArthur model
-   SF%VEG_LSET_ROS_BACK     = VEG_LSET_ROS_BACK     ! back fire rate of spread, McArthur model
    SF%VEG_LSET_ROS_FIXED    = VEG_LSET_ROS_FIXED    ! ROS is a fixed value and does not change with slope or wind
-   SF%VEG_LSET_WIND_EXP     = VEG_LSET_WIND_EXP     ! exponent on wind cosine in ROS formula
    SF%VEG_LSET_SIGMA        = VEG_LSET_SIGMA * 0.01 ! SAV for Farsite emulation in LSET converted to 1/cm
    SF%VEG_LSET_HT           = VEG_LSET_HT
    SF%VEG_LSET_BETA         = VEG_LSET_BETA
-   SF%VEG_LSET_TAN2         = VEG_LSET_TAN2
    SF%VEG_LSET_IGNITE_T     = VEG_LSET_IGNITE_TIME
    SF%VEG_LSET_QCON         =-VEG_LSET_QCON*1000._EB ! convert from kW/m^2 to W/m^2
    SF%VEG_LSET_M1           = VEG_LSET_M1
@@ -9184,12 +9176,7 @@ SUBROUTINE SET_SURF_DEFAULTS
 
 VEG_LSET_IGNITE_TIME    = 1.E9_EB
 VEG_LSET_ROS_00         = 0.0_EB
-VEG_LSET_ROS_HEAD       = 0.0_EB
-VEG_LSET_ROS_FLANK      = 0.0_EB
-VEG_LSET_ROS_BACK       = 0.0_EB
 VEG_LSET_ROS_FIXED      = .FALSE.
-VEG_LSET_WIND_EXP       = 1.0_EB
-VEG_LSET_TAN2           = .FALSE.
 VEG_LSET_HT             = -1.0_EB
 VEG_LSET_BETA           = -1.0_EB
 VEG_LSET_SIGMA          = -1.0_EB
@@ -16618,11 +16605,6 @@ SUBROUTINE GET_QUANTITY_INDEX(SMOKEVIEW_LABEL,SMOKEVIEW_BAR_LABEL,OUTPUT_INDEX,O
 
       OUTPUT_INDEX = ND
 
-      IF (OUTPUT_QUANTITY(ND)%QUANTITY2_REQUIRED .AND. OUTPUT2_INDEX==0) THEN
-         WRITE(MESSAGE,'(3A)')  'ERROR(1015): Output QUANTITY ',TRIM(QUANTITY),' requires a QUANTITY2.'
-         CALL SHUTDOWN(MESSAGE) ; RETURN
-      ENDIF
-
       IF (OUTPUT_QUANTITY(ND)%SPEC_ID_REQUIRED .AND. (Y_INDEX<1 .AND. Z_INDEX<0)) THEN
          IF (SPEC_ID=='null') THEN
             WRITE(MESSAGE,'(3A)')  'ERROR(1016): Output QUANTITY ',TRIM(QUANTITY),' requires a SPEC_ID.'

Source/read.f90

@@ -1013,9 +1013,8 @@ MODULE TYPES
 
    ! Level Set Firespread
 
-   LOGICAL :: VEG_LSET_SPREAD,VEG_LSET_TAN2,VEG_LSET_ROS_FIXED
-   REAL(EB) :: VEG_LSET_IGNITE_T,VEG_LSET_ROS_HEAD,VEG_LSET_ROS_00,VEG_LSET_QCON,VEG_LSET_ROS_FLANK,VEG_LSET_ROS_BACK, &
-               VEG_LSET_WIND_EXP,VEG_LSET_SIGMA,VEG_LSET_HT,VEG_LSET_BETA,&
+   LOGICAL :: VEG_LSET_SPREAD,VEG_LSET_ROS_FIXED
+   REAL(EB) :: VEG_LSET_IGNITE_T,VEG_LSET_ROS_00,VEG_LSET_QCON,VEG_LSET_SIGMA,VEG_LSET_HT,VEG_LSET_BETA,&
                VEG_LSET_M1,VEG_LSET_M10,VEG_LSET_M100,VEG_LSET_MLW,VEG_LSET_MLH,VEG_LSET_SURF_LOAD,VEG_LSET_FIREBASE_TIME, &
                VEG_LSET_CHAR_FRACTION,VEG_LSET_WIND_HEIGHT, &
                B_ROTH,C_ROTH,BETA_ROTH