MODULE WALL_ROUTINES ! Compute the wall boundary conditions USE PRECISION_PARAMETERSUSE GLOBAL_CONSTANTSUSE MESH_POINTERS IMPLICIT NONEPRIVATECHARACTER(255), PARAMETER :: wallid='$Id: wall.f90 721 2007-10-01 18:33:34Z drjfloyd $'CHARACTER(255), PARAMETER :: wallrev='$Revision: 721 $'CHARACTER(255), PARAMETER :: walldate='$Date: 2007-10-01 14:33:34 -0400 (Mon, 01 Oct 2007) $'PUBLIC WALL_BC,GET_REV_wall CONTAINSSUBROUTINE WALL_BC(T,NM)USE COMP_FUNCTIONS, ONLY: SECONDREAL(EB) :: TNOWREAL(EB), INTENT(IN) :: TINTEGER, INTENT(IN) :: NMIF (EVACUATION_ONLY(NM)) RETURNTNOW=SECOND()CALL POINT_TO_MESH(NM)IF (.NOT. ISOTHERMAL) CALL THERMAL_BC(T)CALL SPECIES_BC(T)CALL DENSITY_BCTUSED(6,NM)=TUSED(6,NM)+SECOND()-TNOWCONTAINS SUBROUTINE THERMAL_BC(T)! Thermal boundary conditions for adiabatic, fixed temperature, fixed flux and interpolated boundaries.! One dimensional heat transfer and pyrolysis is done in PYROLYSIS, which is called at the end of this routine.USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP REAL(EB) :: DT_BC,T,TSI,TMP_G,DTMP,TMP_OTHER,CP_TERM,RHOWAL,RAMP_FACTOR,QNET,FDERIV,TMP_EXTERIORINTEGER  :: IOR,II,JJ,KK,IBC,IIG,JJG,KKG, IWREAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOPTYPE (SURFACE_TYPE), POINTER :: SFTYPE (VENTS_TYPE), POINTER :: VT IF (PREDICTOR) THEN   UU => U   VV => V   WW => W   RHOP => RHOSELSE   UU => US   VV => VS   WW => WS   RHOP => RHOENDIF ! Loop through all boundary cells and apply heat transfer method, except for thermally-thick cells HEAT_FLUX_LOOP: DO IW=1,NWC   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE HEAT_FLUX_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)   IOR = IJKW(4,IW)   IBC = IJKW(5,IW)   SF  => SURFACE(IBC)    METHOD_OF_HEAT_TRANSFER: SELECT CASE(SF%THERMAL_BC_INDEX)       CASE (NO_CONVECTION) METHOD_OF_HEAT_TRANSFER         TMP_F(IW) = TMP(IIG,JJG,KKG)         TMP_W(IW) = TMP_F(IW)      CASE (INFLOW_OUTFLOW) METHOD_OF_HEAT_TRANSFER           TMP_F(IW) = TMP(IIG,JJG,KKG)         TMP_EXTERIOR = TMP_0(KK)         IF (VENT_INDEX(IW)>0) THEN            VT => VENTS(VENT_INDEX(IW))            IF (VT%TMP_EXTERIOR>0._EB) TMP_EXTERIOR = VT%TMP_EXTERIOR         ENDIF         SELECT CASE(IOR)            CASE( 1)                IF (UU(II,JJ,KK)>=0._EB)   TMP_F(IW) = TMP_EXTERIOR            CASE(-1)                IF (UU(II-1,JJ,KK)<=0._EB) TMP_F(IW) = TMP_EXTERIOR            CASE( 2)                IF (VV(II,JJ,KK)>=0._EB)   TMP_F(IW) = TMP_EXTERIOR            CASE(-2)                IF (VV(II,JJ-1,KK)<=0._EB) TMP_F(IW) = TMP_EXTERIOR            CASE( 3)                IF (WW(II,JJ,KK)>=0._EB)   TMP_F(IW) = TMP_EXTERIOR            CASE(-3)                IF (WW(II,JJ,KK-1)<=0._EB) TMP_F(IW) = TMP_EXTERIOR            END SELECT         TMP(II,JJ,KK) = TMP_F(IW)         TMP_W(IW) = TMP_F(IW)       CASE (SPECIFIED_TEMPERATURE) METHOD_OF_HEAT_TRANSFER         TMP_G = TMP(IIG,JJG,KKG)         IF (TW(IW)==T_BEGIN) THEN            TSI = T         ELSE            TSI = T - TW(IW)         ENDIF         TMP_F(IW) = TMP_0(KK) + EVALUATE_RAMP(TSI,SF%TAU(TIME_TEMP),SF%RAMP_INDEX(TIME_TEMP))*(SF%TMP_FRONT-TMP_0(KK))         DTMP = TMP_G - TMP_F(IW)         HEAT_TRANS_COEF(IW) = HEAT_TRANSFER_COEFFICIENT(IW,IIG,JJG,KKG,IOR,TMP_G,DTMP)         QCONF(IW) = HEAT_TRANS_COEF(IW)*DTMP         RHOWAL    = 0.5_EB*(RHOP(IIG,JJG,KKG)+RHO_W(IW))         CP_TERM   = MAX(0._EB,-CP_GAMMA*UW(IW)*RHOWAL)         TMP_W(IW) = ( (RDN(IW)*KW(IW)-0.5_EB*CP_TERM)*TMP_G + CP_TERM*TMP_F(IW)-QCONF(IW) )/(0.5_EB*CP_TERM+RDN(IW)*KW(IW))         TMP_W(IW) = MAX(TMPMIN,TMP_W(IW))      CASE (ADIABATIC_INDEX) METHOD_OF_HEAT_TRANSFER         TMP_G = TMP(IIG,JJG,KKG)         TMP_OTHER = TMP_F(IW)         ADLOOP: DO            DTMP = TMP_G - TMP_OTHER            HEAT_TRANS_COEF(IW) = HEAT_TRANSFER_COEFFICIENT(IW,IIG,JJG,KKG,IOR,TMP_G,DTMP)            IF (RADIATION) THEN               QNET = HEAT_TRANS_COEF(IW)*DTMP + QRADIN(IW) - E_WALL(IW) * SIGMA * TMP_OTHER ** 4               FDERIV = -HEAT_TRANS_COEF(IW) -  4._EB * E_WALL(IW) * SIGMA * TMP_OTHER ** 3            ELSE               QNET = HEAT_TRANS_COEF(IW)*DTMP               FDERIV = -HEAT_TRANS_COEF(IW)            ENDIF            IF (FDERIV /= 0._EB) TMP_OTHER = TMP_OTHER - QNET / FDERIV            IF (ABS(TMP_OTHER - TMP_F(IW)) / TMP_F(IW) < 0.0001) THEN               TMP_F(IW) = TMP_OTHER               EXIT ADLOOP            ELSE               TMP_F(IW) = TMP_OTHER               CYCLE ADLOOP            ENDIF         ENDDO ADLOOP         DTMP = TMP_G - TMP_F(IW)         HEAT_TRANS_COEF(IW) = HEAT_TRANSFER_COEFFICIENT(IW,IIG,JJG,KKG,IOR,TMP_G,DTMP)         QCONF(IW) = HEAT_TRANS_COEF(IW)*DTMP         RHOWAL    = 0.5_EB*(RHOP(IIG,JJG,KKG)+RHO_W(IW))         CP_TERM   = MAX(0._EB,-CP_GAMMA*UW(IW)*RHOWAL)         TMP_W(IW) = ( (RDN(IW)*KW(IW)-0.5_EB*CP_TERM)*TMP_G + CP_TERM*TMP_F(IW)-QCONF(IW) )/(0.5_EB*CP_TERM+RDN(IW)*KW(IW))         TMP_W(IW) = MAX(TMPMIN,TMP_W(IW))      CASE (SPECIFIED_HEAT_FLUX) METHOD_OF_HEAT_TRANSFER         IF (TW(IW)==T_BEGIN) THEN            TSI = T         ELSE            TSI = T - TW(IW)         ENDIF         RAMP_FACTOR = EVALUATE_RAMP(TSI,SF%TAU(TIME_HEAT),SF%RAMP_INDEX(TIME_HEAT))         TMP_F(IW) = TMPA + RAMP_FACTOR*(SF%TMP_FRONT-TMPA)         QCONF(IW) = -RAMP_FACTOR*SF%CONVECTIVE_HEAT_FLUX*AREA_ADJUST(IW)         RHOWAL    = 0.5_EB*(RHOP(IIG,JJG,KKG)+RHO_W(IW))         TMP_G     = TMP(IIG,JJG,KKG)         CP_TERM   = MAX(0._EB,-CP_GAMMA*UW(IW)*RHOWAL)         TMP_W(IW) = ( (RDN(IW)*KW(IW)-0.5_EB*CP_TERM)*TMP_G + CP_TERM*TMP_F(IW)-QCONF(IW) )/(0.5_EB*CP_TERM+RDN(IW)*KW(IW))         TMP_W(IW) = MAX(TMPMIN,TMP_W(IW))       CASE (INTERPOLATED_BC) METHOD_OF_HEAT_TRANSFER          TMP_G = TMP(IIG,JJG,KKG)         TMP_OTHER =OMESH(IJKW(9,IW))%TMP(IJKW(10,IW),IJKW(11,IW),IJKW(12,IW))         IF (CELL_VOLUME_RATIO(IW)<0.5 .OR.  CELL_VOLUME_RATIO(IW)>2.0) THEN            TMP_W(IW) = TMP_G            SELECT CASE(IOR)               CASE( 1)                   IF (UU(II,JJ,KK)>=0._EB)   TMP_W(IW) = TMP_OTHER               CASE(-1)                   IF (UU(II-1,JJ,KK)<=0._EB) TMP_W(IW) = TMP_OTHER               CASE( 2)                   IF (VV(II,JJ,KK)>=0._EB)   TMP_W(IW) = TMP_OTHER               CASE(-2)                   IF (VV(II,JJ-1,KK)<=0._EB) TMP_W(IW) = TMP_OTHER               CASE( 3)                   IF (WW(II,JJ,KK)>=0._EB)   TMP_W(IW) = TMP_OTHER               CASE(-3)                   IF (WW(II,JJ,KK-1)<=0._EB) TMP_W(IW) = TMP_OTHER            END SELECT         ELSE            TMP_W(IW) = TMP_OTHER         ENDIF         TMP_F(IW) = TMP_W(IW)         QCONF(IW) = KW(IW)*(TMP_G-TMP_W(IW))*RDN(IW)         TMP(II,JJ,KK) = TMP_W(IW)    END SELECT METHOD_OF_HEAT_TRANSFER ! Record wall temperature in the ghost cell    IF (SOLID(CELL_INDEX(II,JJ,KK))) TMP(II,JJ,KK) = MAX(100._EB,MIN(4900._EB,TMP_W(IW))) ENDDO HEAT_FLUX_LOOP ! For thermally-thick boundary conditions, call the routine PYROLYSIS IF (CORRECTOR) THEN   WALL_COUNTER = WALL_COUNTER + 1   IF (WALL_COUNTER==WALL_INCREMENT) THEN      DT_BC    = T - BC_CLOCK      BC_CLOCK = T      CALL PYROLYSIS(T,DT_BC)      WALL_COUNTER = 0   ENDIFENDIF END SUBROUTINE THERMAL_BC  SUBROUTINE SPECIES_BC(T)USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP REAL(EB) T,YY_WALL,YY_G,DENOM,YY_OTHER(MAX_SPECIES),YY_EXTERIOR(MAX_SPECIES),RHO_G,UN,DD,EPSB,MFT,TSIINTEGER IBC,IIG,JJG,KKG,IOR,IWB,IW, II, JJ, KK, NTYPE (SURFACE_TYPE), POINTER :: SFTYPE (VENTS_TYPE), POINTER :: VTREAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYPREAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP IF (PREDICTOR) THEN   UU => U   VV => V   WW => W   RHOP => RHOS   IF (N_SPECIES > 0) YYP => YYSELSE   UU => US   VV => VS   WW => WS   RHOP => RHO   IF (N_SPECIES > 0) YYP => YYENDIF  ! Loop through the wall cells, apply mass boundary conditionsWALL_CELL_LOOP: DO IW=1,NWC   IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_CELL_LOOP   IBC = IJKW(5,IW)   SF  => SURFACE(IBC)   IF (N_SPECIES==0 .AND. .NOT. SF%SPECIES_BC_INDEX==SPECIFIED_MASS_FLUX) CYCLE WALL_CELL_LOOP   IF (N_SPECIES==0 .AND. SF%SPECIES_BC_INDEX==SPECIFIED_MASS_FLUX .AND. SF%MASS_FLUX(0) == 0._EB) CYCLE WALL_CELL_LOOP   II  = IJKW(1,IW)   JJ  = IJKW(2,IW)   KK  = IJKW(3,IW)   IOR = IJKW(4,IW)   IIG = IJKW(6,IW)   JJG = IJKW(7,IW)   KKG = IJKW(8,IW)      METHOD_OF_MASS_TRANSFER: SELECT CASE(SF%SPECIES_BC_INDEX)       CASE (NO_MASS_FLUX) METHOD_OF_MASS_TRANSFER          IF (.NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) YY_W(IW,1:N_SPECIES) = YYP(IIG,JJG,KKG,1:N_SPECIES)         IF (BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) THEN            YY_EXTERIOR(1:N_SPECIES) = SPECIES(1:N_SPECIES)%YY0            VT => VENTS(VENT_INDEX(IW))            DO N=1,N_SPECIES               IF (VT%MASS_FRACTION(N)>-1._EB) YY_EXTERIOR(N) = VT%MASS_FRACTION(N)            ENDDO            SELECT CASE(IOR)               CASE( 1)                   IF (UU(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)               CASE(-1)                  IF (UU(II-1,JJ,KK)<0._EB) YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)               CASE( 2)                   IF (VV(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)               CASE(-2)                  IF (VV(II,JJ-1,KK)<0._EB) YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)               CASE( 3)                  IF (WW(II,JJ,KK)>0._EB)   YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)               CASE(-3)                  IF (WW(II,JJ,KK-1)<0._EB) YY_W(IW,1:N_SPECIES)=YY_EXTERIOR(1:N_SPECIES)            END SELECT            YYP(II,JJ,KK,1:N_SPECIES)=YY_W(IW,1:N_SPECIES)         ENDIF         IF ( SF%LEAK_PATH(1)>-1 .AND. UWS(IW)<0._EB) THEN            RHO_G = RHOP(IIG,JJG,KKG)            UN = -UWS(IW)            IF (PREDICTOR) EPSB = -.5_EB*UN**2*DT*RDN(IW)            IF (CORRECTOR) EPSB =  .5_EB*UN**2*DT*RDN(IW)            SPECIES_LOOP_1: DO N=1,N_SPECIES               DD    = RHODW(IW,N)*RDN(IW)               YY_G  = YYP(IIG,JJG,KKG,N)               DENOM = DD + (.5_EB*UN+EPSB)*RHO_W(IW)               YY_W(IW,N) = ( MASSFLUX(IW,N) + YY_G*(DD + (EPSB-.5_EB*UN)*RHO_G) ) / DENOM            ENDDO SPECIES_LOOP_1         ENDIF       CASE (SPECIFIED_MASS_FRACTION) METHOD_OF_MASS_TRANSFER         IF (TW(IW)==T_BEGIN) THEN            IF (PREDICTOR) TSI = T + DT            IF (CORRECTOR) TSI = T         ELSE            IF (PREDICTOR) TSI = T + DT - TW(IW)            IF (CORRECTOR) TSI = T      - TW(IW)         ENDIF         DO N=1,N_SPECIES            IF (UWS(IW) <= 0._EB) THEN               YY_WALL = SPECIES(N)%YY0 + EVALUATE_RAMP(TSI,SF%TAU(N),SF%RAMP_INDEX(N))*(SF%MASS_FRACTION(N)-SPECIES(N)%YY0)            ELSE               YY_WALL = YYP(IIG,JJG,KKG,N)            ENDIF            IF (DNS) YY_W(IW,N) = 2._EB*YY_WALL - YYP(IIG,JJG,KKG,N)            IF (LES) YY_W(IW,N) =       YY_WALL          ENDDO       CASE (SPECIFIED_MASS_FLUX) METHOD_OF_MASS_TRANSFER         ! If the current time is before the "activation" time, TW, apply simple BCs and get out         IF (T < TW(IW) .AND. N_SPECIES > 0) THEN            IF (.NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) YY_W(IW,1:N_SPECIES)      = YYP(IIG,JJG,KKG,1:N_SPECIES)            IF (SOLID(CELL_INDEX(II,JJ,KK)))         YYP(II,JJ,KK,1:N_SPECIES) = YY_W(IW,1:N_SPECIES)            IF (PREDICTOR) UWS(IW)   = 0._EB             MASSFLUX(IW,1:N_SPECIES) = 0._EB            ACTUAL_BURN_RATE(IW)     = 0._EB            CYCLE WALL_CELL_LOOP         ENDIF         MFT = 0._EB         ! If the user has specified the burning rate, evaluate the ramp and other related parameters         SUM_MASSFLUX_LOOP: DO N=0,N_SPECIES            IF (SF%MASS_FLUX(N) > 0._EB) THEN  ! Use user-specified ramp-up of mass flux               IF (TW(IW)==T_BEGIN) THEN                  TSI = T               ELSE                  TSI = T - TW(IW)               ENDIF               MASSFLUX(IW,N) = EVALUATE_RAMP(TSI,SF%TAU(N),SF%RAMP_INDEX(N))*SF%MASS_FLUX(N)               IF (N==I_FUEL) THEN                  IF (EW(IW)>0._EB) MASSFLUX(IW,N) = MASSFLUX(IW,N)*EXP(-EW(IW))                  ACTUAL_BURN_RATE(IW) = MASSFLUX(IW,N)               ENDIF            ENDIF            MASSFLUX(IW,N) = MASSFLUX(IW,N)*AREA_ADJUST(IW)            IF (N==I_FUEL) ACTUAL_BURN_RATE(IW) = ACTUAL_BURN_RATE(IW)*AREA_ADJUST(IW)            MFT = MFT + MASSFLUX(IW,N)         ENDDO SUM_MASSFLUX_LOOP          ! Add total fuel consumed to various summing arrays         CONSUME_FUEL: IF (CORRECTOR .AND. SF%THERMALLY_THICK .AND. I_FUEL /= 0) THEN              IF (SF%SURFACE_DENSITY>0._EB .AND. SF%BACKING==EXPOSED) THEN               IWB = WALL_INDEX_BACK(IW)               IF (BOUNDARY_TYPE(IWB)==SOLID_BOUNDARY) MASS_LOSS(IWB) = MASS_LOSS(IWB) + MASSFLUX(IW,I_FUEL)*DT            ENDIF            MASS_LOSS(IW) = MASS_LOSS(IW) + MASSFLUX(IW,I_FUEL)*DT            OBSTRUCTION(OBST_INDEX_W(IW))%MASS = OBSTRUCTION(OBST_INDEX_W(IW))%MASS - MASSFLUX(IW,I_FUEL)*DT*AW(IW)         ENDIF CONSUME_FUEL         ! Compute the ghost cell value of the species to get the right mass flux          RHO_G = RHOP(IIG,JJG,KKG)         UN    = 2._EB*MFT/(RHO_W(IW)+RHO_G)         IF (PREDICTOR) UWS(IW) = -UN         IF (PREDICTOR) EPSB = -.5_EB*UN**2*DT*RDN(IW)         IF (CORRECTOR) EPSB =  .5_EB*UN**2*DT*RDN(IW)         SPECIES_LOOP: DO N=1,N_SPECIES            DD    = RHODW(IW,N)*RDN(IW)            YY_G  = YYP(IIG,JJG,KKG,N)            DENOM = DD + (.5_EB*UN+EPSB)*RHO_W(IW)            YY_W(IW,N) = ( MASSFLUX(IW,N) + YY_G*(DD + (EPSB-.5_EB*UN)*RHO_G) ) / DENOM         ENDDO SPECIES_LOOP      CASE (INTERPOLATED_BC) METHOD_OF_MASS_TRANSFER          YY_OTHER(1:N_SPECIES) = OMESH(IJKW(9,IW))%YY(IJKW(10,IW),IJKW(11,IW),IJKW(12,IW),1:N_SPECIES)         IF (CELL_VOLUME_RATIO(IW)<0.5_EB .OR. CELL_VOLUME_RATIO(IW)>2._EB) THEN            IF (.NOT.SOLID(CELL_INDEX(IIG,JJG,KKG))) YY_W(IW,1:N_SPECIES) = YYP(IIG,JJG,KKG,1:N_SPECIES)             SELECT CASE(IOR)               CASE( 1)