FDS Source: change wall node and layer removal logic (Issue #14902)

Source/func.f90

@@ -1840,7 +1840,7 @@ REAL(EB) FUNCTION GET_PARTICLE_ENTHALPY(I_LPC,TMP_S)
 USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D_UNIFORM
 REAL(EB), INTENT(IN) :: TMP_S
 REAL(EB) :: RHO,RHO_H,VOL,DTMP,H_S,THICKNESS
-INTEGER :: I,N,ITMP,I_GRAD
+INTEGER :: I,N,ITMP
 INTEGER, INTENT(IN) :: I_LPC
 TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
 TYPE(SURFACE_TYPE), POINTER :: SF
@@ -1852,21 +1852,16 @@ REAL(EB) FUNCTION GET_PARTICLE_ENTHALPY(I_LPC,TMP_S)
    CALL INTERPOLATE1D_UNIFORM(LBOUND(SPECIES(LPC%Y_INDEX)%H_L,1),SPECIES(LPC%Y_INDEX)%H_L,TMP_S,GET_PARTICLE_ENTHALPY)
 ELSE
    SF=>SURFACE(LPC%SURF_INDEX)
-   SELECT CASE(SF%GEOMETRY)
-      CASE(SURF_CARTESIAN)                          ; I_GRAD = 1
-      CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2
-      CASE(SURF_SPHERICAL)                          ; I_GRAD = 3
-   END SELECT
    RHO_H = 0._EB
    RHO = 0._EB
    ITMP = MIN(I_MAX_TEMP-1,INT(TMP_S))
    DTMP = TMP_S-REAL(ITMP,EB)
    THICKNESS = SUM(SF%LAYER_THICKNESS)
    DO I=1,SUM(SF%N_LAYER_CELLS)
       IF (SF%GEOMETRY==SURF_INNER_CYLINDRICAL) THEN
-         VOL = (SF%INNER_RADIUS+SF%X_S(I))**I_GRAD - (SF%INNER_RADIUS+SF%X_S(I-1))**I_GRAD
+         VOL = (SF%INNER_RADIUS+SF%X_S(I))**SF%I_GRAD - (SF%INNER_RADIUS+SF%X_S(I-1))**SF%I_GRAD
       ELSE
-         VOL = (THICKNESS+SF%INNER_RADIUS-SF%X_S(I-1))**I_GRAD - (THICKNESS+SF%INNER_RADIUS-SF%X_S(I))**I_GRAD
+         VOL = (THICKNESS+SF%INNER_RADIUS-SF%X_S(I-1))**SF%I_GRAD - (THICKNESS+SF%INNER_RADIUS-SF%X_S(I))**SF%I_GRAD
       ENDIF
       MATL_REMESH: DO N=1,SF%N_MATL
          IF (SF%RHO_0(I,N)<=TWO_EPSILON_EB) CYCLE MATL_REMESH
@@ -2112,7 +2107,7 @@ REAL(EB) FUNCTION SURFACE_DENSITY(MODE,SF,ONE_D,MATL_INDEX)
 
 INTEGER, INTENT(IN) :: MODE
 INTEGER, INTENT(IN), OPTIONAL :: MATL_INDEX
-INTEGER :: I_GRAD,NWP,II2,N,ITMP
+INTEGER :: NWP,II2,N,ITMP
 REAL(EB) :: WGT,R_S(0:NWP_MAX),EPUM,DTMP
 TYPE(BOUNDARY_ONE_D_TYPE), INTENT(IN), POINTER :: ONE_D
 TYPE(SURFACE_TYPE), INTENT(IN), POINTER :: SF
@@ -2124,12 +2119,6 @@ REAL(EB) FUNCTION SURFACE_DENSITY(MODE,SF,ONE_D,MATL_INDEX)
 
 ELSE THERMALLY_THICK_IF
 
-   SELECT CASE(SF%GEOMETRY)
-      CASE(SURF_CARTESIAN)                           ; I_GRAD = 1
-      CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL)  ; I_GRAD = 2
-      CASE(SURF_SPHERICAL)                           ; I_GRAD = 3
-   END SELECT
-
    NWP = SUM(ONE_D%N_LAYER_CELLS)
    IF (SF%GEOMETRY==SURF_INNER_CYLINDRICAL) THEN
       R_S(0:NWP) = SF%INNER_RADIUS + ONE_D%X(0:NWP)
@@ -2140,8 +2129,9 @@ REAL(EB) FUNCTION SURFACE_DENSITY(MODE,SF,ONE_D,MATL_INDEX)
    SURFACE_DENSITY = 0._EB
    NUMBER_WALL_POINTS_LOOP: DO II2=1,NWP
       AREA_VOLUME_SELECT: SELECT CASE(MODE)
-         CASE(0,2); WGT = ABS(R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(REAL(I_GRAD,EB)*(SF%INNER_RADIUS+SF%THICKNESS)**(I_GRAD-1))
-         CASE(1,3); WGT = ABS(R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(SF%INNER_RADIUS+SF%THICKNESS)**I_GRAD
+         CASE(0,2); WGT = ABS(R_S(II2-1)**SF%I_GRAD-R_S(II2)**SF%I_GRAD)/ &
+                          (REAL(SF%I_GRAD,EB)*(SF%INNER_RADIUS+SF%THICKNESS)**(SF%I_GRAD-1))
+         CASE(1,3); WGT = ABS(R_S(II2-1)**SF%I_GRAD-R_S(II2)**SF%I_GRAD)/(SF%INNER_RADIUS+SF%THICKNESS)**SF%I_GRAD
       END SELECT AREA_VOLUME_SELECT
 
       EPUM = 1._EB ! energy per unit mass
@@ -4063,6 +4053,7 @@ SUBROUTINE PACK_BOUNDARY_ONE_D(NM,IC,RC,LC,OS,OD_INDEX,UNPACK_IT,COUNT_ONLY,CHEC
 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_LAYER_MASS(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)
@@ -4106,6 +4097,7 @@ SUBROUTINE REALLOCATE_BOUNDARY_ONE_D(ONE_D)
    ALLOCATE(ONE_D%LAYER_THICKNESS_OLD(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%MIN_LAYER_MASS)) DEALLOCATE(ONE_D%MIN_LAYER_MASS)         ; ALLOCATE(ONE_D%MIN_LAYER_MASS(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))
@@ -4132,12 +4124,12 @@ END SUBROUTINE REALLOCATE_BOUNDARY_ONE_D
 
 SUBROUTINE INITIALIZE_BOUNDARY_ONE_D(NM,OD_INDEX,SURF_INDEX)
 
-USE GLOBAL_CONSTANTS, ONLY: RADIATION
+USE GLOBAL_CONSTANTS, ONLY: RADIATION,NWP_MAX
 INTEGER, INTENT(IN) :: NM,OD_INDEX,SURF_INDEX
 TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
 TYPE(SURFACE_TYPE), POINTER :: SF
 INTEGER :: NN,I,II
-REAL(EB) :: RAMP_POSITION,TF,TB
+REAL(EB) :: RAMP_POSITION,TF,TB,R(0:NWP_MAX)
 
 ONE_D => MESHES(NM)%BOUNDARY_ONE_D(OD_INDEX)
 SF => SURFACE(SURF_INDEX)
@@ -4153,6 +4145,12 @@ SUBROUTINE INITIALIZE_BOUNDARY_ONE_D(NM,OD_INDEX,SURF_INDEX)
 DO I=1,MIN(ONE_D%N_CELLS_MAX,ONE_D%N_CELLS_INI)
    ONE_D%DX_OLD(I) = ONE_D%X(I)-ONE_D%X(I-1)
 ENDDO
+R(0:ONE_D%N_CELLS_INI) = SF%THICKNESS + SF%INNER_RADIUS - ONE_D%X(0:ONE_D%N_CELLS_INI)
+I = 0
+DO NN = 1,ONE_D%N_LAYERS
+   ONE_D%MIN_LAYER_MASS(NN) = SF%MIN_LAYER_MASS(NN)*(R(I)**SF%I_GRAD - R(I+ONE_D%N_LAYER_CELLS(NN))**SF%I_GRAD)
+   I = I + ONE_D%N_LAYER_CELLS(NN)
+ENDDO
 IF (ALLOCATED(ONE_D%LAYER_THICKNESS_OLD)) ONE_D%LAYER_THICKNESS_OLD(1:ONE_D%N_LAYERS) = SF%LAYER_THICKNESS(1:SF%N_LAYERS)
 IF (SF%RAMP_T_I_INDEX > 0) THEN
    !NOTE: Replicating EVALUATE_RAMP since MODULE MATH_FUNCTIONS uses the MODULE which contains this routine
@@ -4310,6 +4308,7 @@ SUBROUTINE PACK_BOUNDARY_PROP1(NM,IC,RC,LC,OS,B1_INDEX,UNPACK_IT,COUNT_ONLY,SURF
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RHO_G,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%RDN,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%K_G,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%M_DOT_LAYER_PP,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%Q_DOT_G_PP,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%Q_DOT_O2_PP,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),B1%Q_CONDENSE,UNPACK_IT)
@@ -5254,24 +5253,14 @@ END SUBROUTINE GET_WALL_NODE_WEIGHTS
 !> \param X_S_NEW Array of interior cell edge positions after shrinkage or swelling (m)
 !> \param INT_WGT Array of weighting factors for new arrangement of interior cells
 
-SUBROUTINE GET_INTERPOLATION_WEIGHTS(GEOMETRY,NWP,NWP_NEW,INNER_RADIUS,X_S,X_S_NEW,INT_WGT)
+SUBROUTINE GET_INTERPOLATION_WEIGHTS(GEOMETRY,I_GRAD,NWP,NWP_NEW,INNER_RADIUS,X_S,X_S_NEW,INT_WGT)
 
-INTEGER, INTENT(IN)  :: GEOMETRY,NWP,NWP_NEW
+INTEGER, INTENT(IN)  :: GEOMETRY,I_GRAD,NWP,NWP_NEW
 REAL(EB), INTENT(IN) :: X_S(0:NWP), X_S_NEW(0:NWP_NEW), INNER_RADIUS
 REAL(EB), INTENT(OUT) :: INT_WGT(NWP_NEW,NWP)
 
 REAL(EB) XUP,XLOW,XUP_NEW,XLOW_NEW,VOL_NEW,VOL,THICKNESS
-INTEGER I_NEW, I_OLD, I_GRAD
-
-
-SELECT CASE(GEOMETRY)
-   CASE(SURF_CARTESIAN)
-      I_GRAD = 1
-   CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL)
-      I_GRAD = 2
-   CASE(SURF_SPHERICAL)
-      I_GRAD = 3
-END SELECT
+INTEGER I_NEW, I_OLD
 
 I_OLD = 1
 INT_WGT = 0._EB

Source/init.f90

@@ -4188,6 +4188,7 @@ SUBROUTINE FIND_WALL_BACK_INDEX(NM,IW)
    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%MIN_LAYER_THICKNESS) ; ALLOCATE(ONE_D%MIN_LAYER_THICKNESS(ONE_D%N_LAYERS))
+   DEALLOCATE(ONE_D%MIN_LAYER_MASS) ; ALLOCATE(ONE_D%MIN_LAYER_MASS(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%MIN_LAYER_THICKNESS(1:ONE_D%N_LAYERS) = MIN_LAYER_THICKNESS(1:ONE_D%N_LAYERS)

Source/read.f90

@@ -7430,7 +7430,7 @@ SUBROUTINE PROC_MATL
 
 USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM,LINEAR_SYSTEM_SOLVE
 USE PHYSICAL_FUNCTIONS, ONLY: GET_TMP_REF
-INTEGER :: I,N,NN,NL,NLPC,NS,NS2,NR,NR2,Z_INDEX(N_TRACKED_SPECIES,MAX_REACTIONS),IERR,ITMP,I_GRAD,&
+INTEGER :: I,N,NN,NL,NLPC,NS,NS2,NR,NR2,Z_INDEX(N_TRACKED_SPECIES,MAX_REACTIONS),IERR,ITMP,&
            MATL_MATRIX_SIZE,REAC_COUNTER,TEMP_COUNTER,TEMP_MATL(N_MATL,MAX_REACTIONS)
 REAL(EB) :: ANS,H_ADJUST,NU_INERT,H_R_CALC(0:I_MAX_TEMP),SUM_NU(N_MATL,MAX_REACTIONS),DTMP,THICKNESS,VOL,X1
 INTEGER, ALLOCATABLE, DIMENSION(:) :: MATL_MATRIX_POINTER
@@ -7648,11 +7648,6 @@ SUBROUTINE PROC_MATL
                MATL_COEF_VECTOR(REAC_COUNTER) = MATL_COEF_VECTOR(REAC_COUNTER) + ML%NU_LPC(NLPC,NR)*ANS*ML%TMP_REF(NR)
             ELSE LIQUID_IF
                SF=>SURFACE(LPC%SURF_INDEX)
-               SELECT CASE(SF%GEOMETRY)
-                  CASE(SURF_CARTESIAN)                          ; I_GRAD = 1
-                  CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2
-                  CASE(SURF_SPHERICAL)                          ; I_GRAD = 3
-               END SELECT
                ALLOCATE(RHO_H(SF%N_MATL))
                RHO_H = 0._EB
                ALLOCATE(RHO(SF%N_MATL))
@@ -7662,7 +7657,7 @@ SUBROUTINE PROC_MATL
                THICKNESS = SUM(SF%LAYER_THICKNESS)
                X1 = THICKNESS+SF%INNER_RADIUS
                DO I=1,SF%N_LAYERS
-                  VOL = X1**I_GRAD-(X1 - SF%LAYER_THICKNESS(I))**I_GRAD
+                  VOL = X1**SF%I_GRAD-(X1 - SF%LAYER_THICKNESS(I))**SF%I_GRAD
                   X1 = X1 - SF%LAYER_THICKNESS(I)
                   MATL_RHO: DO NN=1,SF%N_MATL
                      IF (SF%MATL_MASS_FRACTION(I,NN)<=TWO_EPSILON_EB) CYCLE MATL_RHO
@@ -7815,7 +7810,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
             REFERENCE_HEAT_FLUX(MAX_QDOTPP_REF),REFERENCE_HEAT_FLUX_TIME_INTERVAL,MINIMUM_SCALING_HEAT_FLUX,&
             MAXIMUM_SCALING_HEAT_FLUX,REFERENCE_THICKNESS(MAX_QDOTPP_REF),&
             AREA_MULTIPLIER,Z_0,PARTICLE_EXTRACTION_VELOCITY,NEAR_WALL_EDDY_VISCOSITY,&
-            DELAMINATION_TMP(MAX_LAYERS),DELAMINATION_DENSITY(MAX_LAYERS)
+            DELAMINATION_TMP(MAX_LAYERS),DELAMINATION_DENSITY(MAX_LAYERS),MINIMUM_LAYER_MASS_FRACTION(MAX_LAYERS)
 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)
@@ -7843,7 +7838,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
                 LEAK_PATH,LEAK_PATH_ID,LENGTH,MASS_FLUX,MASS_FLUX_TOTAL,MASS_FLUX_VAR,MASS_FRACTION,&
                 MASS_TRANSFER_COEFFICIENT, &
                 MATL_ID,MATL_MASS_FRACTION,MASS_PER_VOLUME,MCC_CONVERSION_FACTOR,MINIMUM_BURNOUT_TIME,&
-                MINIMUM_LAYER_THICKNESS,MLRPUA,MOISTURE_CONTENT,MOISTURE_FRACTION,&
+                MINIMUM_LAYER_MASS_FRACTION,MINIMUM_LAYER_THICKNESS,MLRPUA,MOISTURE_CONTENT,MOISTURE_FRACTION,&
                 N_LAYER_CELLS_MAX,NEAR_WALL_EDDY_VISCOSITY,NEAR_WALL_TURBULENCE_MODEL,NET_HEAT_FLUX,&
                 NO_SLIP,NODE_ID,NPPC,NUSSELT_C0,NUSSELT_C1,NUSSELT_C2,NUSSELT_M,&
                 PARTICLE_EXTRACTION_VELOCITY,PARTICLE_MASS_FLUX,PARTICLE_SURFACE_DENSITY,PART_ID,&
@@ -8387,19 +8382,23 @@ SUBROUTINE READ_SURF(QUICK_READ)
    SELECT CASE(GEOMETRY)
       CASE('CARTESIAN')
          SF%GEOMETRY       = SURF_CARTESIAN
+         SF%I_GRAD = 1
          IF (SF%WIDTH>0._EB)                 SF%BACKING = INSULATED
       CASE('CYLINDRICAL')
          SF%GEOMETRY       = SURF_CYLINDRICAL
          IF (SF%INNER_RADIUS<TWO_EPSILON_EB) SF%BACKING = INSULATED
+         SF%I_GRAD = 2
       CASE('INNER CYLINDRICAL')
          SF%GEOMETRY       = SURF_INNER_CYLINDRICAL
          IF (SF%INNER_RADIUS<TWO_EPSILON_EB) THEN
             WRITE(MESSAGE,'(A,A,A)') 'ERROR(312): SURF ',TRIM(SF%ID),' needs an INNER_RADIUS.'
             CALL SHUTDOWN(MESSAGE) ; RETURN
          ENDIF
+         SF%I_GRAD = 2
       CASE('SPHERICAL')
          SF%GEOMETRY       = SURF_SPHERICAL
          IF (SF%INNER_RADIUS<TWO_EPSILON_EB) SF%BACKING = INSULATED
+         SF%I_GRAD = 3
       CASE DEFAULT
          WRITE(MESSAGE,'(A,A,A)') 'ERROR(312): SURF ',TRIM(SF%ID),' GEOMETRY not recognized.'
          CALL SHUTDOWN(MESSAGE) ; RETURN
@@ -8684,6 +8683,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
    ALLOCATE(SF%HEAT_SOURCE(MAX_LAYERS))             ; SF%HEAT_SOURCE = 0._EB
    ALLOCATE(SF%RAMP_IHS_INDEX(MAX_LAYERS))          ; SF%RAMP_IHS_INDEX = 0
    ALLOCATE(SF%MIN_LAYER_THICKNESS(MAX_LAYERS))     ; SF%MIN_LAYER_THICKNESS = 0._EB
+   ALLOCATE(SF%MIN_LAYER_MASS(MAX_LAYERS))          ; SF%MIN_LAYER_MASS = 1._EB
    ALLOCATE(SF%SWELL_RATIO(MAX_LAYERS))             ; SF%SWELL_RATIO = 1._EB
 
    COUNT_LAYERS: DO NL=1,MAX_LAYERS
@@ -8719,7 +8719,7 @@ SUBROUTINE READ_SURF(QUICK_READ)
       IF (SF%LAYER_DENSITY(NL) > 0._EB) SF%LAYER_DENSITY(NL) = 1./SF%LAYER_DENSITY(NL)
       SF%THICKNESS = SF%THICKNESS + SF%LAYER_THICKNESS(NL)
    ENDDO COUNT_LAYERS
-
+   SF%MIN_LAYER_MASS = SF%LAYER_DENSITY*MINIMUM_LAYER_MASS_FRACTION
    ! Define mass flux division point if the user does not specify. For all but
    ! surfaces with exposed backing, all pyrolyzed mass migrates to the front surface.
 
@@ -9014,6 +9014,7 @@ SUBROUTINE SET_SURF_DEFAULTS
 MINIMUM_SCALING_HEAT_FLUX = 0._EB
 MAXIMUM_SCALING_HEAT_FLUX = 1.E9_EB
 MINIMUM_BURNOUT_TIME    = 1.E6_EB
+MINIMUM_LAYER_MASS_FRACTION = 1.E-12_EB
 MINIMUM_LAYER_THICKNESS = -1.E-6_EB  ! The absolute value is the default, m
 MLRPUA                  = 0._EB
 MOISTURE_CONTENT        = 0._EB
@@ -9210,6 +9211,7 @@ SUBROUTINE PROC_SURF_1
                      CALL SHUTDOWN(MESSAGE) ; RETURN
                   ENDIF
                ENDIF
+            CASE(SURF_SPHERICAL)
          END SELECT
       ENDIF
    ENDDO
@@ -9251,7 +9253,7 @@ SUBROUTINE PROC_SURF_2
 LOGICAL :: INDEX_LIST(MAX_LPC)
 REAL(EB) :: R_L(0:MAX_LAYERS), FUEL_MF,HRRPUA,E,MF_SUM,HR_SUM
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: Q_REF_INT
-INTEGER  :: I_GRAD,I_CONE_RAMP,E_PTR
+INTEGER  :: I_CONE_RAMP,E_PTR
 INTEGER  :: PROCESSED(N_RAMP)
 LOGICAL :: BURNING,BLOWING,SUCKING
 TYPE(RAMPS_TYPE),POINTER :: RP,RP_E2T,RP_INT,RP_QREF
@@ -9265,12 +9267,6 @@ SUBROUTINE PROC_SURF_2
    SF => SURFACE(N)
    IF (SF%THERMAL_BC_INDEX==THERMALLY_THICK) SF%INCLUDE_BOUNDARY_ONE_D_TYPE = .TRUE.
 
-   SELECT CASE(SF%GEOMETRY)
-      CASE(SURF_CARTESIAN)                           ; I_GRAD = 1
-      CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL)  ; I_GRAD = 2
-      CASE(SURF_SPHERICAL)                           ; I_GRAD = 3
-   END SELECT
-
    ! Particle Information
 
    SF%PART_INDEX = 0
@@ -9560,7 +9556,7 @@ SUBROUTINE PROC_SURF_2
    DO NL=1,SF%N_LAYERS
       R_L(NL) = R_L(NL-1)-SF%LAYER_THICKNESS(NL)
       SF%SURFACE_DENSITY = SF%SURFACE_DENSITY + SF%LAYER_DENSITY(NL) * &
-         (R_L(NL-1)**I_GRAD-R_L(NL)**I_GRAD)/(REAL(I_GRAD,EB)*SF%THICKNESS**(I_GRAD-1))
+         (R_L(NL-1)**SF%I_GRAD-R_L(NL)**SF%I_GRAD)/(REAL(SF%I_GRAD,EB)*SF%THICKNESS**(SF%I_GRAD-1))
    ENDDO
 
    IF ((ABS(SF%SURFACE_DENSITY) <= TWO_EPSILON_EB) .AND. SF%BURN_AWAY) THEN

Source/type.f90

@@ -241,6 +241,7 @@ MODULE TYPES
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: LAYER_THICKNESS     !< (1:N_LAYERS) Thickness of layer (m)
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: LAYER_THICKNESS_OLD !< (1:N_LAYERS) Thickness of layer (m) at last remesh
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_LAYER_THICKNESS !< (1:N_LAYERS) Minimum thickness of layer (m)
+   REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_LAYER_MASS      !< (1:N_LAYERS) Layer mass for removal (m)
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_DIFFUSIVITY     !< (1:N_LAYERS) Min diffusivity of all matls in layer (m2/s)
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: DDSUM               !< Scaling factor to get minimum cell size
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: SMALLEST_CELL_SIZE  !< Minimum cell size (m)
@@ -321,6 +322,7 @@ MODULE TYPES
    REAL(EB) :: RHO_G                 !< Gas density in near wall cell (kg/m3)
    REAL(EB) :: RDN=1._EB             !< \f$ 1/ \delta n \f$ at the surface (1/m)
    REAL(EB) :: K_G=0.025_EB          !< Thermal conductivity of gas in adjacent gas phase cell near wall
+   REAL(EB) :: M_DOT_LAYER_PP=0._EB  !< Mass of a layer removed in kg used to update OB%MASS
    REAL(EB) :: Q_DOT_G_PP=0._EB      !< Heat release rate per unit area (W/m2)
    REAL(EB) :: Q_DOT_O2_PP=0._EB     !< Heat release rate per unit area (W/m2) due to oxygen consumption
    REAL(EB) :: Q_CONDENSE=0._EB      !< Heat release rate per unit area (W/m2) due to gas condensation
@@ -922,6 +924,7 @@ MODULE TYPES
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: MASS_FRACTION,MASS_FLUX,DDSUM,SMALLEST_CELL_SIZE,SWELL_RATIO
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: REFERENCE_HEAT_FLUX !< Reference flux for the flux scaling pyrolysis model (kW/m2)
    REAL(EB), ALLOCATABLE, DIMENSION(:) :: SPYRO_TH_FACTOR     !< Thickness scaling factor for the flux scaling pyrolysis model
+   REAL(EB), ALLOCATABLE, DIMENSION(:) :: MIN_LAYER_MASS      !< Initial mass in each layer
 
    TYPE(RAMP_ID_TYPE),  ALLOCATABLE, DIMENSION(:) :: RAMP
    INTEGER, DIMENSION(3) :: RGB
@@ -938,6 +941,7 @@ MODULE TYPES
    INTEGER :: N_CONE_CURVES=0 !< Total number of time vs. heat release rate curves specified for the S-pyro model
    INTEGER :: N_MATL          !< Total number of unique materials over all layers
    INTEGER :: NODE_INDEX=0    !< Ductnode index
+   INTEGER :: I_GRAD = 1      !< Geometry factor for volume
    INTEGER, DIMENSION(30) :: ONE_D_REALS_ARRAY_SIZE=0,ONE_D_INTEGERS_ARRAY_SIZE=0,ONE_D_LOGICALS_ARRAY_SIZE=0
    INTEGER, ALLOCATABLE, DIMENSION(:) :: N_LAYER_CELLS      !< Number of wall cells per material layer
    INTEGER, ALLOCATABLE, DIMENSION(:) :: LAYER_INDEX        !< The layer each wall cell belongs to