Merge pull request #14041 from mcgratta/pyro_depth

FDS Source: Fix problem with BURN_AWAY and VARIABLE_THICKNESS at boundaries

Author: mcgratta

Source/data.f90

@@ -1509,6 +1509,10 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
 OUTPUT_QUANTITY(-47)%INSIDE_SOLID = .TRUE.
 OUTPUT_QUANTITY(-47)%BNDF_APPROPRIATE = .FALSE.
 
+OUTPUT_QUANTITY(-48)%NAME= 'PYROLYSIS DEPTH'
+OUTPUT_QUANTITY(-48)%UNITS= 'm'
+OUTPUT_QUANTITY(-48)%SHORT_NAME= 'p-depth'
+
 OUTPUT_QUANTITY(-51)%NAME = 'ENTHALPY FLUX WALL'
 OUTPUT_QUANTITY(-51)%UNITS= 'kW/m2'
 OUTPUT_QUANTITY(-51)%SHORT_NAME = 'hf'

Source/dump.f90

@@ -8929,6 +8929,13 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL_IN
          SOLID_PHASE_OUTPUT = 0.5_EB*( ONE_D%X(I_DEPTH-1) + ONE_D%X(I_DEPTH) )
       ENDIF
 
+   CASE(48) ! PYROLYSIS DEPTH
+      IF (SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN
+         SOLID_PHASE_OUTPUT = ONE_D%PYROLYSIS_DEPTH
+      ELSE
+         SOLID_PHASE_OUTPUT = 0._EB
+      ENDIF
+
    CASE(51)  ! ENTHALPY FLUX WALL
       ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
       CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,B1%TMP_F)
@@ -9793,7 +9800,7 @@ SUBROUTINE DUMP_PROF(T,NM)
       IF (NWP==0) CYCLE PROF_LOOP
       CALL GET_WALL_NODE_WEIGHTS(NWP,ONE_D%N_LAYERS,ONE_D%N_LAYER_CELLS,ONE_D%LAYER_THICKNESS,SF%GEOMETRY, &
          ONE_D%X(0:NWP),SF%LAYER_DIVIDE,DX_S(1:NWP),RDX_S(0:NWP+1),RDXN_S(0:NWP),DX_WGT_S(0:NWP),DXF,DXB,LAYER_INDEX,MF_FRAC,&
-         SF%INNER_RADIUS)
+         SF%INNER_RADIUS,ONE_D%PYROLYSIS_DEPTH)
    ELSE
       NWP = SF%N_CELLS_INI
       IF (NWP==0) CYCLE PROF_LOOP

Source/func.f90

@@ -659,6 +659,7 @@ MODULE COMP_OPERATORS
    MODULE PROCEDURE EQUATE_INTEGERS
    MODULE PROCEDURE EQUATE_INTEGER_VECTORS
    MODULE PROCEDURE EQUATE_LOGICALS
+   MODULE PROCEDURE EQUATE_LOGICAL_VECTORS
 END INTERFACE
 
 CONTAINS
@@ -715,6 +716,16 @@ SUBROUTINE EQUATE_LOGICALS(A,B,SWAP)
 ENDIF
 END SUBROUTINE EQUATE_LOGICALS
 
+SUBROUTINE EQUATE_LOGICAL_VECTORS(A,B,SWAP)
+LOGICAL, INTENT(INOUT), DIMENSION(:) :: A,B
+LOGICAL, INTENT(IN) :: SWAP
+IF (SWAP) THEN
+   B = A
+ELSE
+   A = B
+ENDIF
+END SUBROUTINE EQUATE_LOGICAL_VECTORS
+
 END MODULE COMP_OPERATORS
 
 
@@ -1786,56 +1797,39 @@ SUBROUTINE PACK_BOUNDARY_ONE_D(NM,IC,RC,LC,OS,OD_INDEX,UNPACK_IT,COUNT_ONLY)
    IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC) , ONE_D%RAMP_IHS_INDEX(NL) , UNPACK_IT)
 ENDDO
 
-I1 = RC+1 ; RC = I1 + ONE_D%N_MATL - 1
-IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%M_DOT_S_PP(1:RC-I1+1) , UNPACK_IT)
-
-I1 = RC+1 ; RC = I1 + ONE_D%N_CELLS_MAX
-IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%X(0:RC-I1) , UNPACK_IT)
-
-I1 = RC+1 ; RC = I1 + ONE_D%N_CELLS_OLD - 1
-IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%DX_OLD(1:RC-I1+1) , 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%TMP(0:RC-I1) , 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%DELTA_TMP(0:RC-I1) , UNPACK_IT)
-
-DO NL=1,ONE_D%N_LAYERS
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%LAYER_THICKNESS(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%LAYER_THICKNESS_OLD(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%MINIMUM_LAYER_THICKNESS(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%MIN_DIFFUSIVITY(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%DDSUM(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%SMALLEST_CELL_SIZE(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%STRETCH_FACTOR(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%HEAT_SOURCE(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%CELL_SIZE_FACTOR(NL) , UNPACK_IT)
-   RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC) , ONE_D%CELL_SIZE(NL) , UNPACK_IT)
-   LC=LC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%LOGICALS(LC) , ONE_D%HT3D_LAYER(NL) , UNPACK_IT)
+DO NN=1,ONE_D%N_MATL
+IC=IC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%INTEGERS(IC) , ONE_D%MATL_INDEX(NN) , 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)
-
-I1 = RC+1 ; RC = I1 + ONE_D%N_CELLS_MAX + 1
-IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%K_S(0:RC-I1) , UNPACK_IT)
-
-I1 = RC+1 ; RC = I1 + ONE_D%N_LPC - 1
-IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%PART_MASS(1:RC-I1+1) , UNPACK_IT)
-
-I1 = RC+1 ; RC = I1 + ONE_D%N_LPC - 1
-IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%PART_ENTHALPY(1:RC-I1+1) , UNPACK_IT)
+RC=RC+1                             ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC)   ,ONE_D%PYROLYSIS_DEPTH        , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_MATL-1      ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%M_DOT_S_PP(1:RC-I1+1)  , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_CELLS_MAX   ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%X(0:RC-I1)             , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_CELLS_OLD-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%DX_OLD(1:RC-I1+1)      , 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%TMP(0:RC-I1)           , 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%DELTA_TMP(0:RC-I1)     , 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)
+I1=RC+1 ; RC=I1+ONE_D%N_CELLS_MAX+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%K_S(0:RC-I1)           , UNPACK_IT)
+
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%LAYER_THICKNESS(1:RC-I1+1)    , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%LAYER_THICKNESS_OLD(1:RC-I1+1), UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%MIN_LAYER_THICKNESS(1:RC-I1+1), UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%MIN_DIFFUSIVITY(1:RC-I1+1)    , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%DDSUM(1:RC-I1+1)              , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%SMALLEST_CELL_SIZE(1:RC-I1+1) , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%STRETCH_FACTOR(1:RC-I1+1)     , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%HEAT_SOURCE(1:RC-I1+1)        , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%CELL_SIZE_FACTOR(1:RC-I1+1)   , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%CELL_SIZE(1:RC-I1+1)          , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LPC-1    ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%PART_MASS(1:RC-I1+1)          , UNPACK_IT)
+I1=RC+1 ; RC=I1+ONE_D%N_LPC-1    ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC),ONE_D%PART_ENTHALPY(1:RC-I1+1)      , UNPACK_IT)
+
+LC=LC+1                          ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%LOGICALS(LC)    , ONE_D%INTERNAL_RADIATION    , UNPACK_IT)
+I1=LC+1 ; LC=I1+ONE_D%N_LAYERS-1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%LOGICALS(I1:LC) , ONE_D%HT3D_LAYER(1:LC-I1+1) , UNPACK_IT)
 
 DO NN=1,ONE_D%N_MATL
-   I1 = RC+1 ; RC = I1 + ONE_D%N_LAYERS - 1
+   I1=RC+1 ; RC=I1+ONE_D%N_LAYERS-1
    IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%MATL_COMP(NN)%MASS_FRACTION(1:RC-I1+1) , UNPACK_IT)
-ENDDO
-
-DO NN=1,ONE_D%N_MATL
-   I1 = RC+1 ; RC = I1 + (ONE_D%N_CELLS_MAX+2) - 1
+   I1=RC+1 ; RC=I1+ONE_D%N_CELLS_MAX+1
    IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(I1:RC) , ONE_D%MATL_COMP(NN)%RHO(0:RC-I1) , UNPACK_IT)
 ENDDO
 
@@ -1860,8 +1854,8 @@ SUBROUTINE REALLOCATE_BOUNDARY_ONE_D(ONE_D)
 IF (ALLOCATED(ONE_D%LAYER_THICKNESS)) DEALLOCATE(ONE_D%LAYER_THICKNESS)       ; ALLOCATE(ONE_D%LAYER_THICKNESS(ONE_D%N_LAYERS))
 IF (ALLOCATED(ONE_D%LAYER_THICKNESS_OLD)) DEALLOCATE(ONE_D%LAYER_THICKNESS_OLD) 
    ALLOCATE(ONE_D%LAYER_THICKNESS_OLD(ONE_D%N_LAYERS))
-IF (ALLOCATED(ONE_D%MINIMUM_LAYER_THICKNESS)) DEALLOCATE(ONE_D%MINIMUM_LAYER_THICKNESS)       
-   ALLOCATE(ONE_D%MINIMUM_LAYER_THICKNESS(ONE_D%N_LAYERS))
+IF (ALLOCATED(ONE_D%MIN_LAYER_THICKNESS)) DEALLOCATE(ONE_D%MIN_LAYER_THICKNESS)       
+   ALLOCATE(ONE_D%MIN_LAYER_THICKNESS(ONE_D%N_LAYERS))
 IF (ALLOCATED(ONE_D%HT3D_LAYER)) DEALLOCATE(ONE_D%HT3D_LAYER)                 ; ALLOCATE(ONE_D%HT3D_LAYER(ONE_D%N_LAYERS))
 IF (ALLOCATED(ONE_D%MIN_DIFFUSIVITY)) DEALLOCATE(ONE_D%MIN_DIFFUSIVITY)       ; ALLOCATE(ONE_D%MIN_DIFFUSIVITY(ONE_D%N_LAYERS))
 IF (ALLOCATED(ONE_D%RHO_C_S)) DEALLOCATE(ONE_D%RHO_C_S)                       ; ALLOCATE(ONE_D%RHO_C_S(ONE_D%N_CELLS_MAX))
@@ -1936,7 +1930,7 @@ SUBROUTINE INITIALIZE_BOUNDARY_ONE_D(NM,OD_INDEX,SURF_INDEX)
 ENDIF
 ONE_D%DELTA_TMP = 0._EB
 ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS) = SF%LAYER_THICKNESS(1:SF%N_LAYERS)
-ONE_D%MINIMUM_LAYER_THICKNESS(1:ONE_D%N_LAYERS) = SF%MINIMUM_LAYER_THICKNESS(1:SF%N_LAYERS)
+ONE_D%MIN_LAYER_THICKNESS(1:ONE_D%N_LAYERS) = SF%MIN_LAYER_THICKNESS(1:SF%N_LAYERS)
 ONE_D%HT3D_LAYER(1:ONE_D%N_LAYERS) = SF%HT3D_LAYER(1:SF%N_LAYERS)
 ONE_D%MIN_DIFFUSIVITY(1:ONE_D%N_LAYERS) = SF%MIN_DIFFUSIVITY(1:SF%N_LAYERS)
 ONE_D%STRETCH_FACTOR(1:ONE_D%N_LAYERS)  = SF%STRETCH_FACTOR(1:SF%N_LAYERS)
@@ -2855,20 +2849,21 @@ END SUBROUTINE GET_WALL_NODE_COORDINATES
 !> \param LAYER_INDEX Array of indices indicating the layer to which each interior cell belongs
 !> \param MF_FRAC Array containing the fraction of each cells mass that is assigned to the front surface
 !> \param INNER_RADIUS Inner radius of hollow cylinder or sphere (m)
+!> \param X_DIVIDE Depth at which pyrolyzates move to back side (m)
 
 SUBROUTINE GET_WALL_NODE_WEIGHTS(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
-         LAYER_THICKNESS,GEOMETRY,X_S,LAYER_DIVIDE,DX,RDX,RDXN,DX_WGT,DXF,DXB,LAYER_INDEX,MF_FRAC,INNER_RADIUS)
+         LAYER_THICKNESS,GEOMETRY,X_S,LAYER_DIVIDE,DX,RDX,RDXN,DX_WGT,DXF,DXB,LAYER_INDEX,MF_FRAC,INNER_RADIUS,X_DIVIDE)
 
 ! Get the wall internal coordinates
 
 INTEGER, INTENT(IN)     :: N_CELLS, N_LAYERS, N_LAYER_CELLS(N_LAYERS),GEOMETRY
 REAL(EB), INTENT(IN)    :: X_S(0:N_CELLS),LAYER_THICKNESS(1:N_LAYERS),LAYER_DIVIDE,INNER_RADIUS
 INTEGER, INTENT(OUT)    :: LAYER_INDEX(0:N_CELLS+1)
 REAL(EB), INTENT(OUT)   :: DX(1:N_CELLS),RDX(0:N_CELLS+1),RDXN(0:N_CELLS),DX_WGT(0:N_CELLS),DXF,DXB, &
-                           MF_FRAC(1:N_CELLS)
+                           MF_FRAC(1:N_CELLS),X_DIVIDE
 
 INTEGER :: I, II, NL, I_GRAD
-REAL(EB) :: R, THICKNESS, X_DIVIDE
+REAL(EB) :: R, THICKNESS
 
    THICKNESS = SUM(LAYER_THICKNESS)
 
@@ -2955,6 +2950,7 @@ SUBROUTINE GET_WALL_NODE_WEIGHTS(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
    IF (LAYER_DIVIDE >= REAL(N_LAYERS,EB)) THEN
 
       MF_FRAC = 1._EB
+      X_DIVIDE = THICKNESS
 
    ELSE
 
@@ -6172,4 +6168,4 @@ SUBROUTINE ACCUMULATE_STRING(STRING_SIZE,MYSTR,ACCSTR,ACCSTR_T_LEN,ACCSTR_USE_LE
 END SUBROUTINE ACCUMULATE_STRING
 
 END MODULE MISC_FUNCTIONS
-   
+   

Source/init.f90

@@ -1652,7 +1652,7 @@ END SUBROUTINE ADJUST_HT3D_WALL_CELLS
 
 SUBROUTINE REALLOCATE_ONE_D_ARRAYS(NM,WALL_CELL,THIN_WALL_CELL)
 
-USE GEOMETRY_FUNCTIONS, ONLY: GET_N_LAYER_CELLS,GET_WALL_NODE_COORDINATES,GET_WALL_NODE_WEIGHTS
+USE GEOMETRY_FUNCTIONS, ONLY: GET_N_LAYER_CELLS,GET_WALL_NODE_COORDINATES
 USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_REAL_ARRAY,REALLOCATE_INTEGER_ARRAY,PACK_WALL,PACK_THIN_WALL
 INTEGER, INTENT(IN) :: NM
 INTEGER, INTENT(IN), OPTIONAL :: WALL_CELL,THIN_WALL_CELL
@@ -3804,7 +3804,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
            N_MATL_OBST,N_LAYERS,N_MATLS,IIF,JJF,KKF
 INTEGER, DIMENSION(MAX_MATERIALS) :: MATL_INDEX_OBST,MATL_INDEX
 REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: MATL_MASS_FRACTION_OBST,MATL_MASS_FRACTION
-REAL(EB), DIMENSION(0:MAX_LAYERS) :: LAYER_THICKNESS,MINIMUM_LAYER_THICKNESS
+REAL(EB), DIMENSION(0:MAX_LAYERS) :: LAYER_THICKNESS,MIN_LAYER_THICKNESS
 REAL(EB), DIMENSION(MAX_LAYERS) :: LAYER_THICKNESS_OBST,HEAT_SOURCE,HEAT_SOURCE_OBST,&
                                    STRETCH_FACTOR,STRETCH_FACTOR_OBST,CELL_SIZE,CELL_SIZE_OBST,&
                                    CELL_SIZE_FACTOR,CELL_SIZE_FACTOR_OBST
@@ -4044,7 +4044,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
    MATL_INDEX(1:N_MATLS) = MATL_INDEX_OBST(1:N_MATLS)  ! MATL_INDEX_OBST is taken from the OBSTs that make up the solid
    MATL_MASS_FRACTION = 0._EB
    LAYER_THICKNESS = 0._EB
-   MINIMUM_LAYER_THICKNESS = 0._EB
+   MIN_LAYER_THICKNESS = 0._EB
    HT3D_LAYER = .FALSE.
    FRONT_LINING_THICKNESS = 0._EB
    BACK_LINING_THICKNESS = 0._EB
@@ -4077,7 +4077,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
       IF (.NOT.SF%LINING) EXIT
       N_LAYERS = N_LAYERS + 1
       LAYER_THICKNESS(N_LAYERS)  = SF%LAYER_THICKNESS(N_LAYERS) 
-      MINIMUM_LAYER_THICKNESS(N_LAYERS)  = SF%MINIMUM_LAYER_THICKNESS(N_LAYERS) 
+      MIN_LAYER_THICKNESS(N_LAYERS)  = SF%MIN_LAYER_THICKNESS(N_LAYERS) 
       HT3D_LAYER(N_LAYERS)       = .FALSE.
       HEAT_SOURCE(N_LAYERS)      = SF%HEAT_SOURCE(N_LAYERS)
       RAMP_IHS_INDEX(N_LAYERS)   = SF%RAMP_IHS_INDEX(N_LAYERS)
@@ -4103,7 +4103,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
    DO NL=1,N_LAYERS_OBST
       N_LAYERS = N_LAYERS + 1
       LAYER_THICKNESS(N_LAYERS)  = LAYER_THICKNESS_OBST(NL)
-      MINIMUM_LAYER_THICKNESS(N_LAYERS) = SF%MINIMUM_LAYER_THICKNESS(1)
+      MIN_LAYER_THICKNESS(N_LAYERS) = SF%MIN_LAYER_THICKNESS(1)
       HT3D_LAYER(N_LAYERS)       = .TRUE.
       HEAT_SOURCE(N_LAYERS)      = HEAT_SOURCE_OBST(NL)
       RAMP_IHS_INDEX(N_LAYERS)   = RAMP_IHS_INDEX_OBST(NL)
@@ -4124,7 +4124,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
       IF (.NOT.SF_BACK%LINING) EXIT
       N_LAYERS = N_LAYERS + 1
       LAYER_THICKNESS(N_LAYERS)  = SF_BACK%LAYER_THICKNESS(SF_BACK%N_LAYERS-NL+1)
-      MINIMUM_LAYER_THICKNESS(N_LAYERS)  = SF_BACK%MINIMUM_LAYER_THICKNESS(SF_BACK%N_LAYERS-NL+1)
+      MIN_LAYER_THICKNESS(N_LAYERS)  = SF_BACK%MIN_LAYER_THICKNESS(SF_BACK%N_LAYERS-NL+1)
       HT3D_LAYER(N_LAYERS)       = .FALSE.
       HEAT_SOURCE(N_LAYERS)      = SF_BACK%HEAT_SOURCE(SF_BACK%N_LAYERS-NL+1)
       RAMP_IHS_INDEX(N_LAYERS)   = SF_BACK%RAMP_IHS_INDEX(SF_BACK%N_LAYERS-NL+1)
@@ -4147,10 +4147,10 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
    DEALLOCATE(ONE_D%MATL_COMP)       ; ALLOCATE(ONE_D%MATL_COMP(ONE_D%N_MATL))
    DEALLOCATE(ONE_D%MATL_INDEX)      ; ALLOCATE(ONE_D%MATL_INDEX(ONE_D%N_MATL))
    DEALLOCATE(ONE_D%LAYER_THICKNESS) ; ALLOCATE(ONE_D%LAYER_THICKNESS(ONE_D%N_LAYERS))
-   DEALLOCATE(ONE_D%MINIMUM_LAYER_THICKNESS) ; ALLOCATE(ONE_D%MINIMUM_LAYER_THICKNESS(ONE_D%N_LAYERS))
+   DEALLOCATE(ONE_D%MIN_LAYER_THICKNESS) ; ALLOCATE(ONE_D%MIN_LAYER_THICKNESS(ONE_D%N_LAYERS))
    DEALLOCATE(ONE_D%HT3D_LAYER)      ; ALLOCATE(ONE_D%HT3D_LAYER(ONE_D%N_LAYERS))
    ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS) = LAYER_THICKNESS(1:ONE_D%N_LAYERS)
-   ONE_D%MINIMUM_LAYER_THICKNESS(1:ONE_D%N_LAYERS) = MINIMUM_LAYER_THICKNESS(1:ONE_D%N_LAYERS)
+   ONE_D%MIN_LAYER_THICKNESS(1:ONE_D%N_LAYERS) = MIN_LAYER_THICKNESS(1:ONE_D%N_LAYERS)
    ONE_D%HT3D_LAYER(1:ONE_D%N_LAYERS) = HT3D_LAYER(1:ONE_D%N_LAYERS)
    DO NN=1,ONE_D%N_MATL
       ALLOCATE(ONE_D%MATL_COMP(NN)%MASS_FRACTION(ONE_D%N_LAYERS))
@@ -4960,7 +4960,7 @@ SUBROUTINE REASSIGN_WALL_CELLS(T,NM)
 
 SUBROUTINE GET_BOUNDARY_TYPE
 
-INTEGER :: IOR,IIG,JJG,KKG,IW_OLD,IERR,PRESSURE_BC_TYPE,ICG_OLD,II
+INTEGER :: IOR,IIG,JJG,KKG,ICO,IW_OLD,IERR,PRESSURE_BC_TYPE,ICG_OLD,II
 TYPE (BOUNDARY_PROP1_TYPE), POINTER :: B1,B1_OLD
 TYPE (BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D_OLD
 TYPE (WALL_TYPE), POINTER :: WC_OLD
@@ -5036,42 +5036,66 @@ SUBROUTINE GET_BOUNDARY_TYPE
 ! Special cases 2: HT3D solid shifts the position of the burned away surface to the exposed surface position.
 
 SF => SURFACE(WC%SURF_INDEX)
+
 IF (REMOVE .AND. ( (SF%THERMAL_BC_INDEX==THERMALLY_THICK.AND.(SF%VARIABLE_THICKNESS.OR.SF%HT_DIM>1)) &
                   .OR. SF%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED ) ) THEN
+
    BC => MESHES(NM)%BOUNDARY_COORD(WC%BC_INDEX)
    IIG = BC%IIG
    JJG = BC%JJG
    KKG = BC%KKG
    IOR = BC%IOR
    ICG_OLD  = 0
    SELECT CASE(IOR)
-      CASE(-1) ; IF (IIG>1)    ICG_OLD = CELL_INDEX(IIG-1,JJG,KKG)
-      CASE( 1) ; IF (IIG<IBAR) ICG_OLD = CELL_INDEX(IIG+1,JJG,KKG)
-      CASE(-2) ; IF (JJG>1)    ICG_OLD = CELL_INDEX(IIG,JJG-1,KKG)
-      CASE( 2) ; IF (JJG<JBAR) ICG_OLD = CELL_INDEX(IIG,JJG+1,KKG)
-      CASE(-3) ; IF (KKG>1)    ICG_OLD = CELL_INDEX(IIG,JJG,KKG-1)
-      CASE( 3) ; IF (KKG<KBAR) ICG_OLD = CELL_INDEX(IIG,JJG,KKG+1)
+      CASE(-1) ; ICG_OLD = CELL_INDEX(IIG-1,JJG,KKG)
+      CASE( 1) ; ICG_OLD = CELL_INDEX(IIG+1,JJG,KKG)
+      CASE(-2) ; ICG_OLD = CELL_INDEX(IIG,JJG-1,KKG)
+      CASE( 2) ; ICG_OLD = CELL_INDEX(IIG,JJG+1,KKG)
+      CASE(-3) ; ICG_OLD = CELL_INDEX(IIG,JJG,KKG-1)
+      CASE( 3) ; ICG_OLD = CELL_INDEX(IIG,JJG,KKG+1)
    END SELECT
-   IW_OLD = CELL(ICG_OLD)%WALL_INDEX(-IOR)
-   IF (IW_OLD>0) THEN
-      WC_OLD => MESHES(NM)%WALL(IW_OLD)
+
+   IF (MESHES(NM)%CELL(ICG_OLD)%EXTERIOR) THEN
+      EWC => MESHES(NM)%EXTERNAL_WALL(CELL(ICG)%WALL_INDEX(IOR))
+      NOM = EWC%NOM
+      IF (NOM>0) THEN
+         IIO = EWC%IIO_MIN
+         JJO = EWC%JJO_MIN
+         KKO = EWC%KKO_MIN
+         ICG_OLD = MESHES(NOM)%CELL_INDEX(IIO,JJO,KKO)
+         IW_OLD = MESHES(NOM)%CELL(ICG_OLD)%WALL_INDEX(-IOR)
+      ELSE
+         IW_OLD = 0
+      ENDIF
+   ELSE
+      NOM = NM
+      IW_OLD = CELL(ICG_OLD)%WALL_INDEX(-IOR)
+   ENDIF
+
+   SWAP: IF (IW_OLD>0) THEN
+      WC_OLD => MESHES(NOM)%WALL(IW_OLD)
+      IF (WC_OLD%OD_INDEX==0) EXIT SWAP
       IF (SF%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED) THEN
          B1     => MESHES(NM)%BOUNDARY_PROP1(WC%B1_INDEX)
-         B1_OLD => MESHES(NM)%BOUNDARY_PROP1(WC_OLD%B1_INDEX)
+         B1_OLD => MESHES(NOM)%BOUNDARY_PROP1(WC_OLD%B1_INDEX)
          IF (WC_OLD%SURF_INDEX==WC%SURF_INDEX) B1%T_IGN = B1_OLD%T_IGN
-      ELSEIF (.NOT.CELL(ICG_OLD)%SOLID .AND. .NOT.CELL(ICG)%SOLID .AND. CELL(IC)%SOLID .AND. &
+      ELSEIF (.NOT.MESHES(NOM)%CELL(ICG_OLD)%SOLID .AND. .NOT.MESHES(NM)%CELL(ICG)%SOLID .AND. MESHES(NM)%CELL(IC)%SOLID .AND. &
               SUM(BOUNDARY_ONE_D(WC_OLD%OD_INDEX)%N_LAYER_CELLS(:))>0) THEN
-         WC%OD_INDEX = WC_OLD%OD_INDEX
+         IF (NOM/=NM) THEN
+            MESHES(NM)%BOUNDARY_ONE_D(WC%OD_INDEX) = MESHES(NOM)%BOUNDARY_ONE_D(WC_OLD%OD_INDEX)
+         ELSE
+            WC%OD_INDEX = WC_OLD%OD_INDEX
+         ENDIF
          WC%BOUNDARY_TYPE = SOLID_BOUNDARY
-         ONE_D_OLD => MESHES(NM)%BOUNDARY_ONE_D(WC_OLD%OD_INDEX)
+         ONE_D_OLD => MESHES(NOM)%BOUNDARY_ONE_D(WC_OLD%OD_INDEX)
          IF (ONE_D_OLD%BACK_MESH>0 .AND. ONE_D_OLD%BACK_MESH/=NM) THEN
             OS => OMESH(ONE_D_OLD%BACK_MESH)%WALL_SEND_BUFFER
             DO II=1,OS%N_ITEMS
                IF (OS%ITEM_INDEX(II)==IW_OLD) OS%ITEM_INDEX(II) = IW
             ENDDO
          ENDIF
       ENDIF
-   ENDIF
+   ENDIF SWAP
 ENDIF
 
 END SUBROUTINE GET_BOUNDARY_TYPE