MODULE MASS ! Compute the mass equation differences  USE PRECISION_PARAMETERSUSE MESH_POINTERSIMPLICIT NONEPRIVATECHARACTER(255), PARAMETER :: massid='$Id: mass.f90 709 2007-09-28 17:47:43Z mcgratta $'CHARACTER(255), PARAMETER :: massrev='$Revision: 709 $'CHARACTER(255), PARAMETER :: massdate='$Date: 2007-09-28 13:47:43 -0400 (Fri, 28 Sep 2007) $'REAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYPREAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP,DPPUBLIC MASS_FINITE_DIFFERENCES,DENSITY,GET_REV_mass  CONTAINS SUBROUTINE MASS_FINITE_DIFFERENCES(NM)USE COMP_FUNCTIONS, ONLY: SECONDUSE GLOBAL_CONSTANTS, ONLY: N_SPECIES,ISOTHERMAL,NULL_BOUNDARY,POROUS_BOUNDARY,PREDICTOR,CORRECTOR,EVACUATION_ONLYINTEGER, INTENT(IN) :: NMREAL(EB) :: FXYZ,PMDT,UDRHODN,TNOWINTEGER  :: I,J,K,N,II,JJ,KK,IIG,JJG,KKG,IW,IORREAL(EB), POINTER, DIMENSION(:) :: UWPREAL(EB), POINTER, DIMENSION(:,:,:) :: UDRHODX,VDRHODY,WDRHODZ,EPSX,EPSY,EPSZ IF (EVACUATION_ONLY(NM)) RETURNIF (SOLID_PHASE_ONLY) RETURNTNOW=SECOND()CALL POINT_TO_MESH(NM) IF (PREDICTOR) THEN   UU => U   VV => V   WW => W   DP => D   RHOP => RHO   UWP  => UW   PMDT = DTELSE   UU => US   VV => VS   WW => WS   DP => DS   RHOP => RHOS   UWP  => UWS   PMDT = -DTENDIF! Define local CFL numbers EPSX => WORK1EPSY => WORK2EPSZ => WORK3DO K=0,KBAR   DO J=0,JBAR      DO I=0,IBAR         EPSX(I,J,K) = PMDT*UU(I,J,K)*RDXN(I)         EPSY(I,J,K) = PMDT*VV(I,J,K)*RDYN(J)         EPSZ(I,J,K) = PMDT*WW(I,J,K)*RDZN(K)      ENDDO   ENDDOENDDO ! Compute spatial differences for density equation NOT_ISOTHERMAL_IF: IF (.NOT.ISOTHERMAL) THEN    UDRHODX => WORK4   VDRHODY => WORK5   WDRHODZ => WORK6      DO K=0,KBAR      DO J=0,JBAR         DO I=0,IBAR            UDRHODX(I,J,K) = UU(I,J,K)*(RHOP(I+1,J,K)-RHOP(I,J,K))*RDXN(I)            VDRHODY(I,J,K) = VV(I,J,K)*(RHOP(I,J+1,K)-RHOP(I,J,K))*RDYN(J)            WDRHODZ(I,J,K) = WW(I,J,K)*(RHOP(I,J,K+1)-RHOP(I,J,K))*RDZN(K)         ENDDO      ENDDO   ENDDO    WLOOP: DO IW=1,NWC      IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WLOOP      II  = IJKW(1,IW)       IIG = IJKW(6,IW)      JJ  = IJKW(2,IW)       JJG = IJKW(7,IW)      KK  = IJKW(3,IW)       KKG = IJKW(8,IW)      IOR = IJKW(4,IW)      UDRHODN = UWP(IW)*(RHO_W(IW)-RHOP(IIG,JJG,KKG))*RDN(IW)      SELECT CASE(IOR)         CASE( 1)            UDRHODX(II,JJ,KK)   = UDRHODN         CASE(-1)             UDRHODX(II-1,JJ,KK) = UDRHODN         CASE( 2)             VDRHODY(II,JJ,KK)   = UDRHODN         CASE(-2)             VDRHODY(II,JJ-1,KK) = UDRHODN         CASE( 3)             WDRHODZ(II,JJ,KK)   = UDRHODN         CASE(-3)             WDRHODZ(II,JJ,KK-1) = UDRHODN      END SELECT   ENDDO WLOOP      DO K=1,KBAR      DO J=1,JBAR         DO I=1,IBAR            IF (SOLID(CELL_INDEX(I,J,K))) CYCLE            FXYZ   = .5_EB*(UDRHODX(I,J,K)  *(1._EB-EPSX(I,J,K))   +  &                            UDRHODX(I-1,J,K)*(1._EB+EPSX(I-1,J,K)) +  &                            VDRHODY(I,J,K)  *(1._EB-EPSY(I,J,K))   +  &                            VDRHODY(I,J-1,K)*(1._EB+EPSY(I,J-1,K)) +  &                            WDRHODZ(I,J,K)  *(1._EB-EPSZ(I,J,K))   +  &                            WDRHODZ(I,J,K-1)*(1._EB+EPSZ(I,J,K-1)) )            FRHO(I,J,K) = FXYZ + RHOP(I,J,K)*DP(I,J,K)         ENDDO      ENDDO   ENDDO ENDIF NOT_ISOTHERMAL_IF ! Compute the species equation differences IF (N_SPECIES > 0) THEN   IF (PREDICTOR) YYP => YY   IF (CORRECTOR) YYP => YYS   UDRHODX => WORK4   VDRHODY => WORK5   WDRHODZ => WORK6ENDIF SPECIES_LOOP: DO N=1,N_SPECIES    DO K=0,KBAR      DO J=0,JBAR         DO I=0,IBAR            UDRHODX(I,J,K) = UU(I,J,K)*( RHOP(I+1,J,K)*YYP(I+1,J,K,N) - RHOP(I,J,K)*YYP(I,J,K,N) )*RDXN(I)            VDRHODY(I,J,K) = VV(I,J,K)*( RHOP(I,J+1,K)*YYP(I,J+1,K,N) - RHOP(I,J,K)*YYP(I,J,K,N) )*RDYN(J)            WDRHODZ(I,J,K) = WW(I,J,K)*( RHOP(I,J,K+1)*YYP(I,J,K+1,N) - RHOP(I,J,K)*YYP(I,J,K,N) )*RDZN(K)         ENDDO      ENDDO   ENDDO    ! Correct U d(RHO*Y)/dx etc. on boundaries    WLOOP2: DO IW=1,NWC      IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WLOOP2      II  = IJKW(1,IW)       IIG = IJKW(6,IW)      JJ  = IJKW(2,IW)       JJG = IJKW(7,IW)      KK  = IJKW(3,IW)       KKG = IJKW(8,IW)      IOR = IJKW(4,IW)      UDRHODN = UWP(IW)*( RHO_W(IW)*YY_W(IW,N) - RHOP(IIG,JJG,KKG)*YYP(IIG,JJG,KKG,N) )*RDN(IW)      SELECT CASE(IOR)         CASE( 1)            UDRHODX(II,JJ,KK)   = UDRHODN         CASE(-1)            UDRHODX(II-1,JJ,KK) = UDRHODN         CASE( 2)            VDRHODY(II,JJ,KK)   = UDRHODN         CASE(-2)             VDRHODY(II,JJ-1,KK) = UDRHODN         CASE( 3)             WDRHODZ(II,JJ,KK)   = UDRHODN         CASE(-3)             WDRHODZ(II,JJ,KK-1) = UDRHODN      END SELECT   ENDDO WLOOP2   ! Sum up the convective and diffusive terms in the transport equation and store in DEL_RHO_D_DEL_Y    DO K=1,KBAR      DO J=1,JBAR         DO I=1,IBAR            FXYZ   = .5_EB*(UDRHODX(I,J,K)  *(1._EB-EPSX(I,J,K))   +  &                            UDRHODX(I-1,J,K)*(1._EB+EPSX(I-1,J,K)) +  &                            VDRHODY(I,J,K)  *(1._EB-EPSY(I,J,K))   +  &                            VDRHODY(I,J-1,K)*(1._EB+EPSY(I,J-1,K)) +  &                            WDRHODZ(I,J,K)  *(1._EB-EPSZ(I,J,K))   +  &                            WDRHODZ(I,J,K-1)*(1._EB+EPSZ(I,J,K-1)) )             DEL_RHO_D_DEL_Y(I,J,K,N) = -DEL_RHO_D_DEL_Y(I,J,K,N) + FXYZ + RHOP(I,J,K)*YYP(I,J,K,N)*DP(I,J,K)          ENDDO      ENDDO   ENDDO ENDDO SPECIES_LOOP TUSED(3,NM)=TUSED(3,NM)+SECOND()-TNOWEND SUBROUTINE MASS_FINITE_DIFFERENCES  SUBROUTINE DENSITY(NM)! Update the density and species mass fractionsUSE COMP_FUNCTIONS, ONLY: SECOND USE PHYSICAL_FUNCTIONS, ONLY : GET_MOLECULAR_WEIGHT2USE GLOBAL_CONSTANTS, ONLY: N_SPECIES,CO_PRODUCTION,I_PROG_F,I_PROG_CO,I_FUEL,TMPMAX,TMPMIN,EVACUATION_ONLY,PREDICTOR,CORRECTOR, &                            CHANGE_TIME_STEP,ISOTHERMAL,TMPA,N_SPEC_DILUENTS, N_ZONE,MIXTURE_FRACTION_SPECIES, &                            GAS_SPECIES, MIXTURE_FRACTION, R0 REAL(EB) :: WFAC,DTRATIO,OMDTRATIO,Z_2,TNOWINTEGER  :: I,J,K,NINTEGER, INTENT(IN) :: NMREAL(EB), POINTER, DIMENSION(:,:,:) :: R_SUM_DILUENTS IF (EVACUATION_ONLY(NM)) RETURNIF (SOLID_PHASE_ONLY) RETURNTNOW=SECOND()CALL POINT_TO_MESH(NM)PREDICTOR_STEP: SELECT CASE (PREDICTOR)CASE(.TRUE.) PREDICTOR_STEP   IF (.NOT.CHANGE_TIME_STEP(NM)) THEN      DO N=1,N_SPECIES         DO K=1,KBAR            DO J=1,JBAR               DO I=1,IBAR                  YYS(I,J,K,N) = RHO(I,J,K)*YY(I,J,K,N) - DT*DEL_RHO_D_DEL_Y(I,J,K,N)               ENDDO             ENDDO         ENDDO      ENDDO   ELSE      DTRATIO   = DT/DTOLD      OMDTRATIO = 1._EB - DTRATIO      DO N=1,N_SPECIES         DO K=1,KBAR            DO J=1,JBAR               DO I=1,IBAR                  YYS(I,J,K,N) = OMDTRATIO*RHO(I,J,K) *YY(I,J,K,N) +  DTRATIO*RHOS(I,J,K)*YYS(I,J,K,N)               ENDDO            ENDDO         ENDDO      ENDDO   ENDIF   ! Predict the density at the next time step (RHOS or RHO^*)   IF (.NOT.ISOTHERMAL) THEN      DO K=1,KBAR         DO J=1,JBAR            DO I=1,IBAR               RHOS(I,J,K) = RHO(I,J,K)-DT*FRHO(I,J,K)            ENDDO          ENDDO      ENDDO   ELSE      FORALL (I=0:IBP1,J=0:JBP1,K=0:KBP1) RHOS(I,J,K) = PBAR_S(K,PRESSURE_ZONE(I,J,K))/(TMPA*SPECIES(0)%RCON)      DO N=1,N_SPECIES         WFAC = 1._EB - SPECIES(N)%RCON/SPECIES(0)%RCON         DO K=1,KBAR            DO J=1,JBAR               DO I=1,IBAR                  RHOS(I,J,K) = RHOS(I,J,K) + WFAC*YYS(I,J,K,N)               ENDDO             ENDDO         ENDDO      ENDDO   ENDIF    ! Correct densities above or below clip limits   CALL CHECK_DENSITY   ! Extract mass fraction from RHO * YY   DO N=1,N_SPECIES      DO K=1,KBAR         DO J=1,JBAR            DO I=1,IBAR               YYS(I,J,K,N) = YYS(I,J,K,N)/RHOS(I,J,K)            ENDDO          ENDDO      ENDDO   ENDDO   ! Correct mass fractions above or below clip limits   CALL CHECK_MASS_FRACTION   ! Predict background pressure at next time step   DO I=1,N_ZONE      PBAR_S(:,I) = PBAR(:,I) + D_PBAR_DT(I)*DT   ENDDO   ! Compute mixture fraction and diluent sums: Y_SUM=Sum(Y_i), Z_SUM=Sum(Z_i)   IF (MIXTURE_FRACTION) THEN      Z_SUM  =  0._EB      Y_SUM  =  0._EB      IF (N_SPEC_DILUENTS > 0) THEN         R_SUM_DILUENTS => WORK4         R_SUM_DILUENTS =  0._EB      ENDIF      DO N=1,N_SPECIES         IF (SPECIES(N)%MODE==MIXTURE_FRACTION_SPECIES) Z_SUM = Z_SUM + YYS(:,:,:,N)         IF (SPECIES(N)%MODE==GAS_SPECIES) THEN            Y_SUM = Y_SUM + YYS(:,:,:,N)            R_SUM_DILUENTS(:,:,:) = R_SUM_DILUENTS(:,:,:) + SPECIES(N)%RCON*YYS(:,:,:,N)         ENDIF      ENDDO   ENDIF   ! Compute molecular weight term RSUM=R0*SUM(Y_i/M_i)    IF (N_SPECIES>0 .AND. .NOT.MIXTURE_FRACTION) THEN      RSUM = SPECIES(0)%RCON      DO N=1,N_SPECIES         WFAC = SPECIES(N)%RCON - SPECIES(0)%RCON         RSUM(:,:,:) = RSUM(:,:,:) + WFAC*YYS(:,:,:,N)      ENDDO      IF (ISOTHERMAL) FORALL (I=0:IBP1,J=0:JBP1,K=0:KBP1) RHOS(I,J,K) = PBAR_S(K,PRESSURE_ZONE(I,J,K))/(TMPA*RSUM(I,J,K))   ENDIF   IF (MIXTURE_FRACTION) THEN      DO K=1,KBAR         DO J=1,JBAR            DO I=1,IBAR               IF (CO_PRODUCTION) THEN                  Z_2 = YYS(I,J,K,I_PROG_CO)               ElSE                  Z_2 = 0._EB               ENDIF               CALL GET_MOLECULAR_WEIGHT2(YYS(I,J,K,I_FUEL),Z_2,YYS(I,J,K,I_PROG_F),Y_SUM(I,J,K),RSUM(I,J,K))               RSUM(I,J,K) = R0/RSUM(I,J,K)            ENDDO         ENDDO      ENDDO      IF (N_SPEC_DILUENTS > 0) RSUM = RSUM*(1._EB-Y_SUM) + R_SUM_DILUENTS   ENDIF   ! Extract predicted temperature at next time step from Equation of State   IF (.NOT.ISOTHERMAL) THEN      IF (N_SPECIES==0) THEN         FORALL (I=0:IBP1,J=0:JBP1,K=0:KBP1) TMP(I,J,K) = PBAR_S(K,PRESSURE_ZONE(I,J,K))/(SPECIES(0)%RCON*RHOS(I,J,K))      ELSE         FORALL (I=0:IBP1,J=0:JBP1,K=0:KBP1) TMP(I,J,K) = PBAR_S(K,PRESSURE_ZONE(I,J,K))/(RSUM(I,J,K)*RHOS(I,J,K))      ENDIF      TMP = MAX(TMPMIN,MIN(TMPMAX,TMP))   ENDIF! The CORRECTOR step   CASE(.FALSE.) PREDICTOR_STEP   ! Correct species mass fraction at next time step (YY here actually means YY*RHO)   DO N=1,N_SPECIES      DO K=1,KBAR         DO J=1,JBAR            DO I=1,IBAR               YY(I,J,K,N) = .5_EB*(RHO(I,J,K)*YY(I,J,K,N) + RHOS(I,J,K)*YYS(I,J,K,N) - DT*DEL_RHO_D_DEL_Y(I,J,K,N) )             ENDDO         ENDDO      ENDDO   ENDDO   ! Correct density at next time step   IF (.NOT.ISOTHERMAL) THEN      DO K=1,KBAR         DO J=1,JBAR            DO I=1,IBAR               RHO(I,J,K) = .5_EB*(RHO(I,J,K)+RHOS(I,J,K)-DT*FRHO(I,J,K))            ENDDO         ENDDO      ENDDO   ELSE      FORALL (I=0:IBP1,J=0:JBP1,K=0:KBP1) RHO(I,J,K) = PBAR(K,PRESSURE_ZONE(I,J,K))/(SPECIES(0)%RCON*TMPA)      DO N=1,N_SPECIES         WFAC = 1._EB - SPECIES(N)%RCON/SPECIES(0)%RCON         FORALL (I=1:IBAR,J=1:JBAR,K=1:KBAR) RHO(I,J,K) = RHO(I,J,K) + WFAC*YY(I,J,K,N)      ENDDO   ENDIF   ! Correct densities above or below clip limits   CALL CHECK_DENSITY    ! Extract Y_n from rho*Y_n   DO N=1,N_SPECIES      DO K=1,KBAR         DO J=1,JBAR            DO I=1,IBAR               YY(I,J,K,N) = YY(I,J,K,N)/RHO(I,J,K)            ENDDO          ENDDO      ENDDO   ENDDO   ! Correct mass fractions above or below clip limits   CALL CHECK_MASS_FRACTION   ! Correct background pressure   DO I=1,N_ZONE      PBAR(:,I) = .5_EB*(PBAR(:,I) + PBAR_S(:,I) + D_PBAR_S_DT(I)*DT)   ENDDO    ! Compute mixture fraction and diluent sums: Y_SUM=Sum(Y_i), Z_SUM=Sum(Z_i)   IF (MIXTURE_FRACTION) THEN      Z_SUM  =  0._EB      Y_SUM  =  0._EB      IF (N_SPEC_DILUENTS > 0) THEN         R_SUM_DILUENTS => WORK4         R_SUM_DILUENTS =  0._EB         ENDIF      DO N=1,N_SPECIES         IF (SPECIES(N)%MODE==MIXTURE_FRACTION_SPECIES) Z_SUM = Z_SUM + YY(:,:,:,N)         IF (SPECIES(N)%MODE==GAS_SPECIES) THEN            Y_SUM = Y_SUM + YY(:,:,:,N)            R_SUM_DILUENTS(:,:,:) = R_SUM_DILUENTS(:,:,:) + SPECIES(N)%RCON*YY(:,:,:,N)         ENDIF      ENDDO   ENDIF   ! Compute molecular weight term RSUM=R0*SUM(Y_i/M_i)    IF (N_SPECIES>0 .AND. .NOT. MIXTURE_FRACTION) THEN      RSUM = SPECIES(0)%RCON      DO N=1,N_SPECIES         WFAC = SPECIES(N)%RCON - SPECIES(0)%RCON         RSUM(:,:,:) = RSUM(:,:,:) + WFAC*YY(:,:,:,N)      ENDDO      IF (ISOTHERMAL) FORALL (I=0:IBP1,J=0:JBP1,K=0:KBP1) RHO(I,J,K) = PBAR(K,PRESSURE_ZONE(I,J,K))/(TMPA*RSUM(I,J,K))   ENDIF   IF (MIXTURE_FRACTION) THEN      DO K=1,KBAR         DO J=1,JBAR            DO I=1,IBAR               IF (CO_PRODUCTION) THEN                  Z_2 = YY(I,J,K,I_PROG_CO)               ElSE                  Z_2 = 0._EB               ENDIF               CALL GET_MOLECULAR_WEIGHT2(YY(I,J,K,I_FUEL),Z_2,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),RSUM(I,J,K))               RSUM(I,J,K) = R0/RSUM(I,J,K)            ENDDO