IF (UU(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES) = YY_OTHER(1:N_SPECIES)               CASE(-1)                  IF (UU(II-1,JJ,KK)<0._EB) YY_W(IW,1:N_SPECIES) = YY_OTHER(1:N_SPECIES)               CASE( 2)                  IF (VV(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES) = YY_OTHER(1:N_SPECIES)               CASE(-2)                  IF (VV(II,JJ-1,KK)<0._EB) YY_W(IW,1:N_SPECIES) = YY_OTHER(1:N_SPECIES)               CASE( 3)                  IF (WW(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES) = YY_OTHER(1:N_SPECIES)               CASE(-3)                  IF (WW(II,JJ,KK-1)<0._EB) YY_W(IW,1:N_SPECIES) = YY_OTHER(1:N_SPECIES)            END SELECT         ELSE            YY_W(IW,1:N_SPECIES) = YY_OTHER(1:N_SPECIES)         ENDIF         YYP(II,JJ,KK,1:N_SPECIES) = YY_W(IW,1:N_SPECIES)   END SELECT METHOD_OF_MASS_TRANSFER   ! Only set species mass fraction in the ghost cell if it is solid       IF (N_SPECIES > 0 .AND. SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) &      YYP(II,JJ,KK,1:N_SPECIES) = YY_W(IW,1:N_SPECIES)ENDDO WALL_CELL_LOOPEND SUBROUTINE SPECIES_BC  SUBROUTINE DENSITY_BC ! Compute density at wall from wall temperatures and mass fractions USE PHYSICAL_FUNCTIONS, ONLY : GET_MOLECULAR_WEIGHT2REAL(EB) :: R_SUM_DILUENTS_W,Y_SUM_W,Z_SUM_W, Z_2, WFACINTEGER  :: IBC,IIG,JJG,KKG,IW,II,JJ,KK, NREAL(EB), POINTER, DIMENSION(:,:) :: PBAR_PTYPE (SURFACE_TYPE), POINTER :: SFREAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP IF (PREDICTOR) THEN    PBAR_P => PBAR_S   RHOP => RHOSELSE    PBAR_P => PBAR    RHOP => RHOENDIFWALL_CELL_LOOP: DO IW=1,NWC    IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_CELL_LOOP   II  = IJKW(1,IW)   JJ  = IJKW(2,IW)   KK  = IJKW(3,IW)   IIG = IJKW(6,IW)   JJG = IJKW(7,IW)   KKG = IJKW(8,IW)   IBC = IJKW(5,IW)   SF => SURFACE(IBC)! Determine ghost cell value of RSUM=R0*Sum(Y_i/M_i) for non-mixture fraction case   IF (.NOT.MIXTURE_FRACTION .AND. N_SPECIES>0) THEN      RSUM_W(IW) = SPECIES(0)%RCON      DO N=1,N_SPECIES         WFAC = SPECIES(N)%RCON - SPECIES(0)%RCON         RSUM_W(IW) = RSUM_W(IW) + WFAC*YY_W(IW,N)      ENDDO   ENDIF ! Determine ghost cell value of RSUM=R0*Sum(Y_i/M_i) for mixture fraction case   IF (MIXTURE_FRACTION) THEN      IF (CO_PRODUCTION) THEN         Z_2 = YY_W(IW,I_PROG_CO)      ELSE         Z_2 = 0._EB      ENDIF      Z_SUM_W  =  0._EB      Y_SUM_W  =  0._EB      IF (N_SPEC_DILUENTS > 0) R_SUM_DILUENTS_W = 0._EB      DO N=1,N_SPECIES         IF (SPECIES(N)%MODE==MIXTURE_FRACTION_SPECIES) THEN            Z_SUM_W = Z_SUM_W + YY_W(IW,N)         ELSEIF (SPECIES(N)%MODE==GAS_SPECIES) THEN            Y_SUM_W = Y_SUM_W + YY_W(IW,N)            R_SUM_DILUENTS_W = R_SUM_DILUENTS_W + SPECIES(N)%RCON*YY_W(IW,N)         ENDIF      ENDDO      CALL GET_MOLECULAR_WEIGHT2(YY_W(IW,I_FUEL),Z_2,YY_W(IW,I_PROG_F),Y_SUM_W,RSUM_W(IW))      RSUM_W(IW) = R0/RSUM_W(IW)      IF (N_SPEC_DILUENTS > 0) RSUM_W(IW) = RSUM_W(IW)*(1._EB-Y_SUM_W) + R_SUM_DILUENTS_W   ENDIF ! Compute ghost cell density   IF (N_SPECIES==0) THEN      RHO_W(IW) = PBAR_P(KK,PRESSURE_ZONE_WALL(IW))/(SPECIES(0)%RCON*TMP_W(IW))   ELSE      RHO_W(IW) = PBAR_P(KK,PRESSURE_ZONE_WALL(IW))/(RSUM_W(IW)*TMP_W(IW))   ENDIF! Actually set the ghost cell value of density in the ghost cell if it is a solid wall   IF ( (SOLID(CELL_INDEX(II,JJ,KK)) .AND. .NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) .OR.  &         BOUNDARY_TYPE(IW)==OPEN_BOUNDARY .OR. BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) THEN      RHOP(II,JJ,KK) = RHO_W(IW)      IF (N_SPECIES>0) RSUM(II,JJ,KK) = RSUM_W(IW)   ENDIF ENDDO WALL_CELL_LOOPEND SUBROUTINE DENSITY_BC  SUBROUTINE PYROLYSIS(T,DT_BC)! Loop through all the boundary cells that require a 1-D heat transfer calcUSE PHYSICAL_FUNCTIONS, ONLY: GET_MOLECULAR_WEIGHT2USE GEOMETRY_FUNCTIONSUSE MATH_FUNCTIONS, ONLY: EVALUATE_RAMPREAL(EB) :: DTMP,QNETF,QDXKF,QDXKB,RR,TMP_G,T,RFACF,RFACB,RFACF2,RFACB2, &            PPCLAUS,PPSURF,TMP_G_B,RHOWAL,CP_TERM,DT_BC, &            DXKF,DXKB,REACTION_RATE,QCONB,DELTA_RHO_S,QRADINB,RFLUX_UP,RFLUX_DOWN,E_WALLB, &            HVRG,Z_2,Y_MF_G,Y_MF_W,YY_S,Y_SUM_W, RSUM_W, RSUM_G, RSUM_S, MFLUX, MFLUX_W, VOLSUM, &            DXF, DXB,HTCB,Q_WATER_F,Q_WATER_B,TMP_F_OLD, DX_GRID, RHO_S0, H_R,DT2_BC,TOLERANCE,C_S_ADJUST_UNITSINTEGER :: IBC,IIG,JJG,KKG,IIB,JJB,KKB,IWB,NWP,I,J,ITMP,NR,NNN,NL,II,JJ,KK,IW,IOR,N,I_OBSTREAL(EB) :: SMALLEST_CELL_SIZE(MAX_LAYERS),THICKNESS,LAYER_THICKNESS_NEW(MAX_LAYERS)REAL(EB),ALLOCATABLE,DIMENSION(:) :: TMP_W_NEWINTEGER  :: N_LAYER_CELLS_NEW(MAX_LAYERS), NWP_NEW,I_GRAD,STEPCOUNTLOGICAL :: POINT_SHRINK, RECOMPUTE,ITERATETYPE (WALL_TYPE), POINTER :: WCTYPE (SURFACE_TYPE), POINTER :: SFTYPE (MATERIAL_TYPE), POINTER :: ML ! Special adjustment of specific heat for steady state applicationsC_S_ADJUST_UNITS = 1000._EB/TIME_SHRINK_FACTOR ! Loop through the thermally-thick wall cellsWALL_CELL_LOOP: DO IW=1,NWC    IF (BOUNDARY_TYPE(IW)/=SOLID_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_CELL_LOOP   IBC = IJKW(5,IW)   SF  => SURFACE(IBC)   IF (SF%THERMAL_BC_INDEX /= THERMALLY_THICK) CYCLE WALL_CELL_LOOP   IF (SUM(WMPUA(IW,:))>0._EB .AND. T>TW(IW)) EW(IW) = EW(IW) + SF%E_COEFFICIENT*SUM(WMPUA(IW,:))*DT_BC   II  = IJKW(1,IW)   JJ  = IJKW(2,IW)   KK  = IJKW(3,IW)   WC  => WALL(IW)   IOR = IJKW(4,IW)   IIG = IJKW(6,IW)   JJG = IJKW(7,IW)   KKG = IJKW(8,IW)   IF (WC%BURNAWAY) CYCLE WALL_CELL_LOOP   SELECT CASE(SF%GEOMETRY)   CASE(SURF_CARTESIAN)      I_GRAD = 0   CASE(SURF_CYLINDRICAL)      I_GRAD = 1   CASE(SURF_SPHERICAL)      I_GRAD = 2   END SELECT   ! Compute convective heat flux at the surface    TMP_G = TMP(IIG,JJG,KKG)   IF (ASSUMED_GAS_TEMPERATURE > 0._EB) TMP_G = ASSUMED_GAS_TEMPERATURE  ! This is just for diagnostic calcs   TMP_F_OLD = TMP_F(IW)   DTMP = TMP_G - TMP_F_OLD   HEAT_TRANS_COEF(IW) = HEAT_TRANSFER_COEFFICIENT(IW,IIG,JJG,KKG,IOR,TMP_G,DTMP)   QCONF(IW) = HEAT_TRANS_COEF(IW)*DTMP    ! Get heat losses from convection and radiation out of back of surface     E_WALLB = MATERIAL(SF%LAYER_MATL_INDEX(SF%N_LAYERS,1))%EMISSIVITY            SELECT CASE(SF%BACKING)      CASE(VOID)  ! Non-insulated backing to an ambient void         DTMP = SF%TMP_BACK - TMP_B(IW)         HTCB = HEAT_TRANSFER_COEFFICIENT(IW,-1,-1,-1,IOR,SF%TMP_BACK,DTMP)         QRADINB   =  E_WALLB*SIGMA*SF%TMP_BACK**4         Q_WATER_B = 0._EB         TMP_G_B = SF%TMP_BACK               CASE(INSULATED)          HTCB      = 0._EB         QRADINB   = 0._EB         E_WALLB   = 0._EB         Q_WATER_B = 0._EB         TMP_G_B   = TMPA                        CASE(EXPOSED)           IWB = WALL_INDEX_BACK(IW)         Q_WATER_B = 0._EB         IF (BOUNDARY_TYPE(IWB)==SOLID_BOUNDARY) THEN            IIB = IJKW(6,IWB)            JJB = IJKW(7,IWB)            KKB = IJKW(8,IWB)            TMP_G_B  = TMP(IIB,JJB,KKB)            DTMP = TMP_G_B - TMP_B(IW)            HTCB = HEAT_TRANSFER_COEFFICIENT(IWB,IIB,JJB,KKB,IOR,TMP_G_B,DTMP)                        HEAT_TRANS_COEF(IWB) = HTCB            QRADINB  = QRADIN(IWB)            IF (NLP>0) Q_WATER_B = -SUM(WCPUA(IWB,:))         ELSE            TMP_G_B  = TMPA            DTMP = TMP_G_B - TMP_B(IW)            HTCB = HEAT_TRANSFER_COEFFICIENT(IW,-1,-1,-1,IOR,TMP_G_B,DTMP)            QRADINB  =  E_WALLB*SIGMA*TMPA4         ENDIF   END SELECT    ! Take away energy flux due to water evaporation    IF (NLP>0) THEN      Q_WATER_F  = -SUM(WCPUA(IW,:))   ELSE      Q_WATER_F  = 0._EB   ENDIF   ! Compute grid for shrinking wall nodes   COMPUTE_GRID: IF (SF%SHRINK) THEN      NWP = SUM(WC%N_LAYER_CELLS)      THICKNESS = SUM(WC%LAYER_THICKNESS)      CALL GET_WALL_NODE_WEIGHTS(NWP,SF%N_LAYERS,WC%N_LAYER_CELLS,THICKNESS,SF%GEOMETRY, &         WC%X_S(0:NWP),DX_S(1:NWP),RDX_S(0:NWP+1),RDXN_S(0:NWP),DX_WGT_S(0:NWP),DXF,DXB,LAYER_INDEX)   ELSE COMPUTE_GRID      NWP         = SF%N_CELLS      DXF         = SF%DXF      DXB         = SF%DXB      DX_S(1:NWP)          = SF%DX(1:NWP)      RDX_S(0:NWP+1)       = SF%RDX(0:NWP+1)      RDXN_S(0:NWP)        = SF%RDXN(0:NWP)      DX_WGT_S(0:NWP)      = SF%DX_WGT(0:NWP)      LAYER_INDEX(0:NWP+1) = SF%LAYER_INDEX(0:NWP+1)   ENDIF COMPUTE_GRID    ! Calculate reaction rates based on the solid phase reactions    Q_S                  = 0._EB   PYROLYSIS_MATERIAL_IF: IF (SF%PYROLYSIS_MODEL==PYROLYSIS_MATERIAL) THEN   MFLUX                = MASSFLUX(IW,I_FUEL)   MASSFLUX(IW,I_FUEL)  = 0._EB   ACTUAL_BURN_RATE(IW) = 0._EB   IF (I_WATER>0) THEN      MFLUX_W           = MASSFLUX(IW,I_WATER)      MASSFLUX(IW,I_WATER)  = 0._EB   ENDIF   POINT_SHRINK         = .FALSE.   WC%SHRINKING         = .FALSE.   IF (SF%SHRINK) X_S_NEW(0:NWP)  = WC%X_S(0:NWP)      POINT_LOOP1: DO I=1,NWP      IF (SF%SHRINK) THEN      POINT_SHRINK = .TRUE.      MATERIAL_LOOP1a: DO N=1,SF%N_MATL         IF (WC%RHO_S(I,N).EQ. 0._EB) CYCLE MATERIAL_LOOP1a         ML  => MATERIAL(SF%MATL_INDEX(N))         IF (ML%N_REACTIONS.EQ.0 .OR. ANY(ML%NU_RESIDUE.GT.0._EB)) THEN            POINT_SHRINK = .FALSE.            EXIT MATERIAL_LOOP1a         ENDIF      ENDDO MATERIAL_LOOP1a      ENDIF               RHO_S0 = SF%LAYER_DENSITY(LAYER_INDEX(I))      VOLSUM = 0._EB      MATERIAL_LOOP1b: DO N=1,SF%N_MATL         ML  => MATERIAL(SF%MATL_INDEX(N))         IF (WC%RHO_S(I,N) .EQ. 0._EB) CYCLE MATERIAL_LOOP1b         IF (ML%PYROLYSIS_MODEL==PYROLYSIS_LIQUID) THEN            VOLSUM = 1._EB            CYCLE MATERIAL_LOOP1b         ENDIF         DO J=1,ML%N_REACTIONS            ! Reaction rate in 1/s            REACTION_RATE = ML%A(J)*(WC%RHO_S(I,N)/RHO_S0)**ML%N_S(J)*EXP(-ML%E(J)/(R0*WC%TMP_S(I)))            ! power term            DTMP = WC%TMP_S(I)-ML%TMP_THR(J)            IF (ML%N_T(J)/=0._EB) THEN               IF (DTMP > 0._EB) THEN                  REACTION_RATE = REACTION_RATE * DTMP**ML%N_T(J)               ELSE                  REACTION_RATE = 0._EB               ENDIF            ELSE ! threshold               IF (DTMP < 0._EB) REACTION_RATE = 0._EB            ENDIF            ! Reaction rate in kg/(m3s)            REACTION_RATE = RHO_S0 * REACTION_RATE            ! Limit reaction rate            REACTION_RATE = MIN(REACTION_RATE , WC%RHO_S(I,N)/DT_BC)            MASSFLUX(IW,I_FUEL)  = MASSFLUX(IW,I_FUEL) +  &               ML%ADJUST_BURN_RATE*ML%NU_FUEL(J)*REACTION_RATE/RDX_S(I)/SF%THICKNESS**I_GRAD            ACTUAL_BURN_RATE(IW) = ACTUAL_BURN_RATE(IW) +  &               ML%NU_FUEL(J)*REACTION_RATE/RDX_S(I)/SF%THICKNESS**I_GRAD            IF (I_WATER>0) MASSFLUX(IW,I_WATER) = MASSFLUX(IW,I_WATER) + &               ML%NU_WATER(J)*REACTION_RATE/RDX_S(I)/SF%THICKNESS**I_GRAD            ITMP = MIN(2000,MAX(0,NINT(WC%TMP_S(I)-TMPA)))            Q_S(I) = Q_S(I) - REACTION_RATE * ML%Q_ARRAY(J,ITMP)            WC%RHO_S(I,N) = WC%RHO_S(I,N) - DT_BC*REACTION_RATE            IF (ML%NU_RESIDUE(J) .GT. 0._EB ) THEN               NNN = SF%RESIDUE_INDEX(N,J)               WC%RHO_S(I,NNN) = WC%RHO_S(I,NNN) + ML%NU_RESIDUE(J)*DT_BC*REACTION_RATE            ENDIF         ENDDO         VOLSUM = VOLSUM + WC%RHO_S(I,N)/ML%RHO_S      ENDDO MATERIAL_LOOP1b      VOLSUM = MIN(VOLSUM,1._EB)      IF (SF%SHRINK) X_S_NEW(I) = X_S_NEW(I-1)+(WC%X_S(I)-WC%X_S(I-1))*VOLSUM      IF (POINT_SHRINK) THEN         IF (VOLSUM<1.0_EB) THEN            WC%SHRINKING=.TRUE.            IF (VOLSUM>0.0_EB) THEN               MATERIAL_LOOP1c: DO N=1,SF%N_MATL                  WC%RHO_S(I,N) = WC%RHO_S(I,N)/VOLSUM               ENDDO MATERIAL_LOOP1c            ENDIF         ENDIF      ENDIF   ENDDO POINT_LOOP1   ! If the fuel or water massflux is non-zero, set the ignition time   IF (TW(IW)>T) THEN      IF (MASSFLUX(IW,I_FUEL) > 0._EB) TW(IW) = T      IF (I_WATER > 0) THEN         IF (MASSFLUX(IW,I_WATER) > 0._EB) TW(IW) = T      ENDIF   ENDIF   ! Special reactions: LIQUID   ! Liquid evaporation can only take place on the surface (1st cell)   POINT_SHRINK = .FALSE.   IF (SF%SHRINK) THEN      POINT_SHRINK = .TRUE.      MATERIAL_LOOP2a: DO N=1,SF%N_MATL         ML  => MATERIAL(SF%MATL_INDEX(N))         IF (ML%PYROLYSIS_MODEL/=PYROLYSIS_LIQUID .AND. WC%RHO_S(1,N)>0._EB) POINT_SHRINK = .FALSE.      ENDDO MATERIAL_LOOP2a      THICKNESS = WC%LAYER_THICKNESS(LAYER_INDEX(1))   ELSE      THICKNESS = SF%LAYER_THICKNESS(LAYER_INDEX(1))   ENDIF   ! Estimate the previous value of liquid mass fluxes. The possibility of multiple liquids not taken into account.   MFLUX = MAX(0._EB,MFLUX - MASSFLUX(IW,I_FUEL))   IF (I_WATER>0) MFLUX_W = MAX(0._EB,MFLUX_W - MASSFLUX(IW,I_WATER))   MATERIAL_LOOP2: DO N=1,SF%N_MATL      ML  => MATERIAL(SF%MATL_INDEX(N))      IF (ML%PYROLYSIS_MODEL/=PYROLYSIS_LIQUID) CYCLE MATERIAL_LOOP2      IF (WC%RHO_S(1,N)==0._EB) CYCLE MATERIAL_LOOP2           IF (ML%NU_FUEL(1)>0._EB) THEN         MFLUX = MFLUX/ML%NU_FUEL(1)      ELSEIF (ML%NU_WATER(1)>0._EB) THEN         IF (I_WATER>0) MFLUX = MFLUX_W/ML%NU_WATER(1)      ENDIF      ! gas phase       Y_MF_G = YY(IIG,JJG,KKG,I_FUEL)      IF (CO_PRODUCTION) THEN         Z_2 = YY(IIG,JJG,KKG,I_PROG_CO)      ELSE         Z_2 = 0._EB      ENDIF      CALL GET_MOLECULAR_WEIGHT2(Y_MF_G,Z_2,YY(IIG,JJG,KKG,I_PROG_F),  Y_SUM(IIG,JJG,KKG),RSUM_G)      ! wall values      Y_MF_W = YY_W(IW,I_FUEL)      IF (CO_PRODUCTION) THEN         Z_2 = YY_W(IW,I_PROG_CO)      ELSE         Z_2 = 0._EB