FDS Source: Issues firemodels#13516 and firemodels#13562. Move INTERN…

…AL_RADIATION to ONE_D

Source/dump.f90

@@ -2970,7 +2970,7 @@ SUBROUTINE INITIALIZE_DIAGNOSTIC_FILE(DT)
    WRITE(LU_OUTPUT,'(A,ES10.3)')  '                                       500 K: ', ML%H(500)*0.001_EB
    WRITE(LU_OUTPUT,'(A,ES10.3)')  '                                       800 K: ', ML%H(800)*0.001_EB
 
-   IF (ML%KAPPA_S<5.0E4_EB) THEN
+   IF (ML%KAPPA_S<4.9E4_EB) THEN
       WRITE(LU_OUTPUT,'(A,F8.2)') '     Absorption coefficient (1/m) ',ML%KAPPA_S
    ENDIF
 

Source/func.f90

@@ -1820,6 +1820,8 @@ SUBROUTINE PACK_BOUNDARY_ONE_D(NM,IC,RC,LC,OS,OD_INDEX,UNPACK_IT,COUNT_ONLY)
    LC=LC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%LOGICALS(LC) , ONE_D%HT3D_LAYER(NL) , UNPACK_IT)
 ENDDO
 
+LC=LC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%LOGICALS(LC) , ONE_D%INTERNAL_RADIATION , UNPACK_IT)
+
 I1 = RC+1 ; RC = I1 + ONE_D%N_CELLS_MAX - 1
 IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%RHO_C_S(1:RC-I1+1) , UNPACK_IT)
 
@@ -1855,7 +1857,6 @@ SUBROUTINE REALLOCATE_BOUNDARY_ONE_D(ONE_D)
 IF (ALLOCATED(ONE_D%N_LAYER_CELLS_MAX)) DEALLOCATE(ONE_D%N_LAYER_CELLS_MAX)   ; ALLOCATE(ONE_D%N_LAYER_CELLS_MAX(ONE_D%N_LAYERS))
 IF (ALLOCATED(ONE_D%RAMP_IHS_INDEX)) DEALLOCATE(ONE_D%RAMP_IHS_INDEX)         ; ALLOCATE(ONE_D%RAMP_IHS_INDEX(ONE_D%N_LAYERS))
 IF (ALLOCATED(ONE_D%MATL_INDEX)) DEALLOCATE(ONE_D%MATL_INDEX)                 ; ALLOCATE(ONE_D%MATL_INDEX(ONE_D%N_MATL))
-
 IF (ALLOCATED(ONE_D%M_DOT_S_PP)) DEALLOCATE(ONE_D%M_DOT_S_PP)                 ; ALLOCATE(ONE_D%M_DOT_S_PP(ONE_D%N_MATL))
 IF (ALLOCATED(ONE_D%X)) DEALLOCATE(ONE_D%X)                                   ; ALLOCATE(ONE_D%X(0:ONE_D%N_CELLS_MAX))
 IF (ALLOCATED(ONE_D%DX_OLD)) DEALLOCATE(ONE_D%DX_OLD)                         ; ALLOCATE(ONE_D%DX_OLD(ONE_D%N_CELLS_OLD))
@@ -1892,6 +1893,7 @@ END SUBROUTINE REALLOCATE_BOUNDARY_ONE_D
 
 SUBROUTINE INITIALIZE_BOUNDARY_ONE_D(NM,OD_INDEX,SURF_INDEX)
 
+USE GLOBAL_CONSTANTS, ONLY: RADIATION
 INTEGER, INTENT(IN) :: NM,OD_INDEX,SURF_INDEX
 TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
 TYPE(SURFACE_TYPE), POINTER :: SF
@@ -1952,9 +1954,11 @@ SUBROUTINE INITIALIZE_BOUNDARY_ONE_D(NM,OD_INDEX,SURF_INDEX)
 ONE_D%SMALLEST_CELL_SIZE(1:SF%N_LAYERS) = SF%SMALLEST_CELL_SIZE(1:SF%N_LAYERS)
 ONE_D%PART_MASS = 0._EB
 ONE_D%PART_ENTHALPY = 0._EB
-DO NN=1,SF%N_MATL
+DO NN=1,ONE_D%N_MATL
    ONE_D%MATL_COMP(NN)%RHO(0:ONE_D%N_CELLS_INI+1) = SF%RHO_0(0:ONE_D%N_CELLS_INI+1,NN)
+   IF (RADIATION .AND. MATERIAL(ONE_D%MATL_INDEX(NN))%KAPPA_S<4.9E4_EB) ONE_D%INTERNAL_RADIATION = .TRUE.
 ENDDO
+IF (RADIATION .AND. ANY(SF%KAPPA_S(1:SF%N_LAYERS)>0._EB)) ONE_D%INTERNAL_RADIATION = .TRUE.
 ONE_D%BACK_MESH = NM
 ONE_D%BACK_SURF = SURF_INDEX
 ONE_D%PYROLYSIS_MODEL = SF%PYROLYSIS_MODEL

Source/init.f90

@@ -4139,6 +4139,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
       DO NL=1,ONE_D%N_LAYERS
          ONE_D%MATL_COMP(NN)%MASS_FRACTION(NL) = MATL_MASS_FRACTION(NL,NN)
       ENDDO
+      IF (MATERIAL(ONE_D%MATL_INDEX(NN))%KAPPA_S<4.9E4_EB) ONE_D%INTERNAL_RADIATION = .TRUE.
    ENDDO
    DEALLOCATE(ONE_D%HEAT_SOURCE)       ; ALLOCATE(ONE_D%HEAT_SOURCE(ONE_D%N_LAYERS))       ; ONE_D%HEAT_SOURCE = 0._EB
    DEALLOCATE(ONE_D%RAMP_IHS_INDEX)    ; ALLOCATE(ONE_D%RAMP_IHS_INDEX(ONE_D%N_LAYERS))    ; ONE_D%RAMP_IHS_INDEX = 0._EB

Source/part.f90

@@ -781,8 +781,9 @@ SUBROUTINE PARTICLE_FACE_INSERT(WALL_INDEX,CFACE_INDEX)
          ENDIF
 
          ! Ember flag to be used for outputs
-         IF (INSERT_TYPE==1 .AND. &
-            (ANY(SF%EMBER_GENERATION_HEIGHT>=0._EB) .OR. SF%EMBER_YIELD>0._EB)) LP%EMBER=.TRUE.
+
+         IF (INSERT_TYPE==1 .AND. (ANY(SF%EMBER_GENERATION_HEIGHT>=0._EB) .OR. SF%EMBER_YIELD>0._EB)) LP%EMBER=.TRUE.
+
          ! Assign particle position on the cell face
 
          CALL RANDOM_NUMBER(RN)
@@ -871,7 +872,7 @@ SUBROUTINE PARTICLE_FACE_INSERT(WALL_INDEX,CFACE_INDEX)
 
          CALL INITIALIZE_SINGLE_PARTICLE
 
-         IF (INSERT_TYPE==2) MESHES(NM)%BOUNDARY_ONE_D(LP%OD_INDEX)%TMP = TMP_PART
+         IF (INSERT_TYPE==2 .AND. LPC%SOLID_PARTICLE) MESHES(NM)%BOUNDARY_ONE_D(LP%OD_INDEX)%TMP = TMP_PART
 
          IF (.NOT.LPC%MASSLESS_TRACER .AND. .NOT.LPC%MASSLESS_TARGET) MASS_SUM = MASS_SUM + LP%PWT*LP%MASS
 

Source/read.f90

@@ -6857,6 +6857,7 @@ SUBROUTINE READ_MATL
 CHARACTER(LABEL_LENGTH), DIMENSION(MAX_MATERIALS,MAX_REACTIONS) :: MATL_ID
 CHARACTER(LABEL_LENGTH), ALLOCATABLE, DIMENSION(:) :: SEARCH_PHRASE,MATL_NAME_RESERVED
 TYPE(MATERIAL_TYPE), POINTER :: ML2=>NULL()
+TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
 INTEGER :: N,NN,NNN,IOS,NR,N_REACTIONS,N_MATL_RESERVED,N_MATL_READ
 NAMELIST /MATL/ A,ABSORPTION_COEFFICIENT,ADJUST_H,ALLOW_SHRINKING,ALLOW_SWELLING,BOILING_TEMPERATURE,&
                 CONDUCTIVITY,CONDUCTIVITY_RAMP,DENSITY,E,EMISSIVITY,FYI,&
@@ -7140,10 +7141,11 @@ SUBROUTINE READ_MATL
       DO NN=1,MAX_MATERIALS
          IF (ML%RESIDUE_MATL_NAME(NN,NR)/='null') ML%N_RESIDUE(NR) = ML%N_RESIDUE(NR) + 1
       ENDDO
-      DO NN=1,MAX_LPC
+      DO NN=1,N_LAGRANGIAN_CLASSES
+         LPC => LAGRANGIAN_PARTICLE_CLASS(NN)
          IF (PART_ID(NN,NR)/='null') THEN
             ML%N_LPC(NR) = ML%N_LPC(NR) + 1
-            LAGRANGIAN_PARTICLE_CLASS(NN)%INCLUDE_BOUNDARY_ONE_D_TYPE = .TRUE.
+            IF (LPC%SOLID_PARTICLE) LPC%INCLUDE_BOUNDARY_ONE_D_TYPE = .TRUE.
          ENDIF
       ENDDO
    ENDDO
@@ -7890,24 +7892,25 @@ SUBROUTINE READ_SURF(QUICK_READ)
 
    IF (ANY(MOISTURE_FRACTION>TWO_EPSILON_EB) .OR. &
        ANY(MASS_PER_VOLUME>TWO_EPSILON_EB)   .OR. &
-       ANY(SURFACE_VOLUME_RATIO>TWO_EPSILON_EB)) THEN
+       ANY(SF%SURFACE_VOLUME_RATIO>TWO_EPSILON_EB)) THEN
 
       ! Determine convective heat transfer coefficient based on element, not surface geometry
 
-      IF (ANY(MASS_PER_VOLUME>0._EB) .AND. SURFACE_VOLUME_RATIO(1)>TWO_EPSILON_EB) &
-         CONVECTION_LENGTH_SCALE = 4._EB/SURFACE_VOLUME_RATIO(1)
+      IF (ANY(MASS_PER_VOLUME>0._EB) .AND. SF%SURFACE_VOLUME_RATIO(1)>TWO_EPSILON_EB) &
+         CONVECTION_LENGTH_SCALE = 4._EB/SF%SURFACE_VOLUME_RATIO(1)
 
       ! Loop over layers and make adjustments to specified densities and moisture content
 
       LAYER_LOOP_2: DO NL=1,MAX_LAYERS
 
          ! Convert SURFACE_VOLUME_RATIO into a THICKNESS
 
-         IF (SURFACE_VOLUME_RATIO(NL)>TWO_EPSILON_EB .AND. THICKNESS(NL)<0._EB .AND. .NOT.HT3D .AND. .NOT.VARIABLE_THICKNESS) THEN
+         IF (SF%SURFACE_VOLUME_RATIO(NL)>TWO_EPSILON_EB .AND. THICKNESS(NL)<0._EB .AND. &
+             .NOT.HT3D .AND. .NOT.VARIABLE_THICKNESS) THEN
             SELECT CASE(GEOMETRY)
-               CASE('CARTESIAN')   ; THICKNESS(NL) = 1._EB/SURFACE_VOLUME_RATIO(NL)
-               CASE('CYLINDRICAL') ; THICKNESS(NL) = 2._EB/SURFACE_VOLUME_RATIO(NL)
-               CASE('SPHERICAL')   ; THICKNESS(NL) = 3._EB/SURFACE_VOLUME_RATIO(NL)
+               CASE('CARTESIAN')   ; THICKNESS(NL) = 1._EB/SF%SURFACE_VOLUME_RATIO(NL)
+               CASE('CYLINDRICAL') ; THICKNESS(NL) = 2._EB/SF%SURFACE_VOLUME_RATIO(NL)
+               CASE('SPHERICAL')   ; THICKNESS(NL) = 3._EB/SF%SURFACE_VOLUME_RATIO(NL)
             END SELECT
          ENDIF
 
@@ -7945,7 +7948,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
                SUM_D = SUM_D + MATL_MASS_FRACTION(NL,NN)/MATERIAL(NNN)%RHO_S
             ENDDO
             SF%PACKING_RATIO(NL) = MASS_PER_VOLUME(NL)*SUM_D
-            SF%KAPPA_S(NL) = SHAPE_FACTOR*SF%PACKING_RATIO(NL)*SURFACE_VOLUME_RATIO(NL)
+            SF%KAPPA_S(NL) = SF%SHAPE_FACTOR*SF%PACKING_RATIO(NL)*SF%SURFACE_VOLUME_RATIO(NL)
             SF%DENSITY_ADJUST_FACTOR(NL,:) = SF%PACKING_RATIO(NL)*SF%DENSITY_ADJUST_FACTOR(NL,:)
             EMISSIVITY      = 1._EB
             EMISSIVITY_BACK = 1._EB
@@ -9105,31 +9108,16 @@ SUBROUTINE PROC_SURF_2
 
    ! Determine if surface has internal radiation
 
-   SF%INTERNAL_RADIATION = .FALSE.
    IF (RADIATION) THEN
       DO NL=1,SF%N_LAYERS
-         IF (SF%KAPPA_S(NL)>0._EB) SF%INTERNAL_RADIATION = .TRUE.
-         DO NN =1,SF%N_LAYER_MATL(NL)
-            ML => MATERIAL(SF%LAYER_MATL_INDEX(NL,NN))
-            IF (ML%KAPPA_S<5.0E4_EB) SF%INTERNAL_RADIATION = .TRUE.
-         ENDDO
-      ENDDO
-   ENDIF
-
-   ! Internal radiation only allowed for Cartesian geometry
-
-   IF (SF%INTERNAL_RADIATION .AND. .NOT.SF%GEOMETRY==SURF_CARTESIAN) THEN
-      WRITE(MESSAGE,'(3A)') 'ERROR(355): SURF ',TRIM(SF%ID),' not Cartesian and cannot have a MATL with an ABSORPTION_COEFFICIENT.'
-      CALL SHUTDOWN(MESSAGE) ; RETURN
-   ENDIF
-
-   ! In case of internal radiation, do not allow zero-emissivity
-
-   IF (SF%INTERNAL_RADIATION) THEN
-      DO NL=1,SF%N_LAYERS
-         DO NN =1,SF%N_LAYER_MATL(NL)
+         DO NN=1,SF%N_LAYER_MATL(NL)
             ML => MATERIAL(SF%LAYER_MATL_INDEX(NL,NN))
-            IF (ML%EMISSIVITY == 0._EB) THEN
+            IF (ML%KAPPA_S<4.9E4_EB .AND. .NOT.SF%GEOMETRY==SURF_CARTESIAN) THEN
+               WRITE(MESSAGE,'(3A)') 'ERROR(355): SURF ',TRIM(SF%ID),&
+                                     ' not Cartesian and cannot have a MATL with an ABSORPTION_COEFFICIENT.'
+               CALL SHUTDOWN(MESSAGE) ; RETURN
+            ENDIF
+            IF (ML%KAPPA_S<4.9E4_EB .AND. ML%EMISSIVITY==0._EB) THEN
                WRITE(MESSAGE,'(5A)') 'ERROR(356): SURF ',TRIM(SF%ID),' zero emissivity of MATL ',&
                                      TRIM(MATL_NAME(SF%LAYER_MATL_INDEX(NL,NN))),' inconsistent with internal radiation.'
                CALL SHUTDOWN(MESSAGE) ; RETURN

Source/type.f90

@@ -252,6 +252,7 @@ MODULE TYPES
    INTEGER :: PYROLYSIS_MODEL=0 !< Indicator of a pyrolysis model used in depth
 
    LOGICAL, ALLOCATABLE, DIMENSION(:) :: HT3D_LAYER           !< (1:ONE_D\%N_LAYERS) Indicator that layer in 3D
+   LOGICAL :: INTERNAL_RADIATION=.FALSE.                      !< Indicator that internal radiation transport done in solid
 
 END TYPE BOUNDARY_ONE_D_TYPE
 
@@ -961,7 +962,7 @@ MODULE TYPES
    CHARACTER(LABEL_LENGTH), ALLOCATABLE, DIMENSION(:) :: MATL_NAME
    CHARACTER(LABEL_LENGTH), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_ID
    REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION
-   LOGICAL :: BURN_AWAY,ADIABATIC,INTERNAL_RADIATION,USER_DEFINED=.TRUE., &
+   LOGICAL :: BURN_AWAY,ADIABATIC,USER_DEFINED=.TRUE., &
               FREE_SLIP=.FALSE.,NO_SLIP=.FALSE.,SPECIFIED_NORMAL_VELOCITY=.FALSE.,SPECIFIED_TANGENTIAL_VELOCITY=.FALSE., &
               SPECIFIED_NORMAL_GRADIENT=.FALSE.,CONVERT_VOLUME_TO_MASS=.FALSE.,&
               BOUNDARY_FUEL_MODEL=.FALSE.,SET_H=.FALSE.,DIRICHLET_FRONT=.FALSE.,DIRICHLET_BACK=.FALSE.,BLOWING=.FALSE.

Source/wall.f90

@@ -1949,7 +1949,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 ! Set mass and energy fluxes to zero prior to time sub-iteration
 
 B1%Q_CON_F = 0._EB
-IF (SF%INTERNAL_RADIATION) Q_RAD_OUT_OLD = B1%Q_RAD_OUT
+IF (ONE_D%INTERNAL_RADIATION) Q_RAD_OUT_OLD = B1%Q_RAD_OUT
 B1%Q_RAD_OUT = 0._EB
 
 IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN
@@ -2208,7 +2208,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 
    ! Calculate internal radiation for Cartesian geometry only
 
-   IF (SF%INTERNAL_RADIATION) THEN
+   IF (ONE_D%INTERNAL_RADIATION) THEN
       DO I=1,NWP
          IF (SF%KAPPA_S(LAYER_INDEX(I))<0._EB) THEN
             VOLSUM = 0._EB
@@ -2249,7 +2249,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
    ! Explicitly update the temperature field and adjust time step if the change in temperature exceeds DELTA_TMP_MAX
 
    IF (ICYC>WALL_INCREMENT) THEN
-      IF (SF%INTERNAL_RADIATION) THEN
+      IF (ONE_D%INTERNAL_RADIATION) THEN
          Q_NET_F = Q_CON_F
          Q_NET_B = Q_CON_B
       ELSE
@@ -2278,7 +2278,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 
    ! Store the mass and energy fluxes from this time sub-iteration
 
-   IF (SF%INTERNAL_RADIATION) THEN
+   IF (ONE_D%INTERNAL_RADIATION) THEN
       Q_RAD_OUT_OLD = Q_RAD_OUT
       B1%Q_RAD_OUT = B1%Q_RAD_OUT + Q_RAD_OUT*DT_BC_SUB/DT_BC
    ENDIF
@@ -2737,7 +2737,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
    IF (SF%DIRICHLET_FRONT) THEN
       RFACF2 = -1._EB
       QDXKF  = 2._EB*B1%TMP_F
-   ELSEIF (.NOT.RADIATION .OR. SF%INTERNAL_RADIATION .OR. ISOLATED_THIN_WALL) THEN
+   ELSEIF (.NOT.RADIATION .OR. ONE_D%INTERNAL_RADIATION .OR. ISOLATED_THIN_WALL) THEN
       RFACF  = 0.5_EB*HTCF
       RFACF2 = (KODXF-RFACF)/(KODXF+RFACF)
       QDXKF  = (HTCF*B1%TMP_G  + Q_LIQUID_F)/(KODXF+RFACF)
@@ -2750,7 +2750,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
    IF (SF%DIRICHLET_BACK) THEN
       RFACB2 = -1._EB
       QDXKB  = 2._EB*B1%TMP_B
-   ELSEIF (.NOT.RADIATION .OR. SF%INTERNAL_RADIATION .OR. ISOLATED_THIN_WALL_BACK) THEN
+   ELSEIF (.NOT.RADIATION .OR. ONE_D%INTERNAL_RADIATION .OR. ISOLATED_THIN_WALL_BACK) THEN
       RFACB  = 0.5_EB*HTCB
       RFACB2 = (KODXB-RFACB)/(KODXB+RFACB)
       QDXKB  = (HTCB*TMP_GAS_BACK + Q_LIQUID_B)/(KODXB+RFACB)
@@ -2791,8 +2791,8 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
    B1%TMP_B  = 0.5_EB*(ONE_D%TMP(NWP)+ONE_D%TMP(NWP+1))
 
    B1%Q_CON_F = B1%Q_CON_F - 0.5_EB*HTCF*DT_BC_SUB*B1%TMP_F
-   IF (RADIATION .AND. .NOT.SF%INTERNAL_RADIATION) B1%Q_RAD_OUT = B1%Q_RAD_OUT + &
-                                                   DT_BC_SUB*(B1%TMP_F_OLD**4+2._EB*B1%TMP_F_OLD**3*(B1%TMP_F-B1%TMP_F_OLD))
+   IF (RADIATION .AND. .NOT.ONE_D%INTERNAL_RADIATION) B1%Q_RAD_OUT = B1%Q_RAD_OUT + &
+                                                      DT_BC_SUB*(B1%TMP_F_OLD**4+2._EB*B1%TMP_F_OLD**3*(B1%TMP_F-B1%TMP_F_OLD))
 
    ! Clipping for excessively high or low temperatures
 
@@ -2826,7 +2826,7 @@ SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,
 IF (ALLOCATED(RHO_DOT) .AND. ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) DEALLOCATE(RHO_DOT)
 
 B1%Q_CON_F = B1%Q_CON_F / DT_BC
-IF (RADIATION .AND. .NOT. SF%INTERNAL_RADIATION) B1%Q_RAD_OUT = B1%Q_RAD_OUT / DT_BC * SIGMA * B1%EMISSIVITY
+IF (RADIATION .AND. .NOT. ONE_D%INTERNAL_RADIATION) B1%Q_RAD_OUT = B1%Q_RAD_OUT / DT_BC * SIGMA * B1%EMISSIVITY
 IF (.NOT. SF%BOUNDARY_FUEL_MODEL) B1%HEAT_TRANS_COEF = HTCF
 
 ! If any gas massflux or particle mass flux is non-zero or the surface temperature exceeds the ignition temperature,