ENDDO ! Adjust FVX, FVY and FVZ at solid, internal obstructions for no flux CALL NO_FLUX END SUBROUTINE VELOCITY_FLUX_ISOTHERMAL  SUBROUTINE VELOCITY_FLUX_CYLINDRICAL(T,NM)USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP USE PHYSICAL_FUNCTIONS, ONLY: GET_MU2! Compute convective and diffusive terms for 2D axisymmetric REAL(EB) :: T,DMUDXINTEGER :: I0INTEGER, INTENT(IN) :: NMREAL(EB) :: MUY,UP,UM,WP,WM,VTRM, &            DTXZDZ,DTXZDX, &            DUDX,DWDZ,DUDZ,DWDX,WOMY,UOMY,PMDT,MPDT, &            S2,C,CDXDYDZTT,AH,RRHO,GX,GZ,&            TXXP,TXXM,TZZP,TZZM,DTXXDX,DTZZDZ, &            EPSUP,EPSUM,EPSWP,EPSWM,MU_SUM,DUMMY=0._EBINTEGER :: II,JJ,KK,I,J,K,IW,IIG,JJG,KKG,ITMP,N,IE,IECREAL(EB), POINTER, DIMENSION(:,:,:) :: TXZ,OMY,UU,WW,RHOP,DPREAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP CALL POINT_TO_MESH(NM) IF (PREDICTOR) THEN   UU => U   WW => W   DP => D     RHOP => RHO   IF (N_SPECIES>0) YYP => YYELSE   UU => US   WW => WS   DP => DS   RHOP => RHOS   IF (N_SPECIES>0) YYP => YYSENDIF TXZ => WORK2OMY => WORK5 CALC_MU: IF (PREDICTOR) THEN        LES_VS_DNS: IF (LES) THEN   ! Smagorinsky model (LES)      C = CSMAG**2      J = 1      KLOOP: DO K=1,KBAR         ILOOP: DO I=1,IBAR            IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ILOOP            CDXDYDZTT = C*DX(I)*DZ(K)            DUDX = RDX(I)*(UU(I,J,K)-UU(I-1,J,K))            DWDZ = RDZ(K)*(WW(I,J,K)-WW(I,J,K-1))            DWDX = 0.25_EB*RDX(I)*(WW(I+1,J,K)-WW(I-1,J,K)+WW(I+1,J,K-1)-WW(I-1,J,K-1))            DUDZ = 0.25_EB*RDZ(K)*(UU(I,J,K+1)-UU(I,J,K-1)+UU(I-1,J,K+1)-UU(I-1,J,K-1))            S2   = 2._EB*(DUDX*DUDX     +  DWDZ*DWDZ ) +  (DUDZ+DWDX)**2 -  TWTH*DP(I,J,K)**2            S2   = MAX(0._EB,S2)            ITMP = 0.1_EB*TMP(I,J,K)            MU(I,J,K) = MAX(SPECIES(0)%MU(ITMP), RHOP(I,J,K)*CDXDYDZTT*SQRT(S2))         ENDDO ILOOP      ENDDO KLOOP   ELSE LES_VS_DNS           ! DNS viscosity      MIXTURE_FRACTION_DNS: IF (.NOT.MIXTURE_FRACTION) THEN         J = 1         DO K=1,KBAR            IVLOOP: DO I=1,IBAR               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE IVLOOP               ITMP = 0.1_EB*TMP(I,J,K)               MU_SUM = SPECIES(0)%MU(ITMP)               DO N=1,N_SPECIES                  MU_SUM = MU_SUM + YYP(I,J,K,N)*(SPECIES(N)%MU(ITMP)-SPECIES(0)%MU(ITMP))               ENDDO               MU(I,J,K) = MU_SUM            ENDDO IVLOOP         ENDDO      ELSE MIXTURE_FRACTION_DNS         J = 1         DO K=1,KBAR            IVLOOP2: DO I=1,IBAR               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE IVLOOP2               ITMP = 0.1_EB*TMP(I,J,K)               IF(CO_PRODUCTION) THEN                  CALL GET_MU2(YY(I,J,K,I_FUEL),YY(I,J,K,I_PROG_CO),YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),MU(I,J,K),ITMP)               ELSE                  CALL GET_MU2(YY(I,J,K,I_FUEL),0._EB,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),MU(I,J,K),ITMP)                                 ENDIF            ENDDO IVLOOP2         ENDDO      ENDIF MIXTURE_FRACTION_DNS    ENDIF LES_VS_DNS ! Mirror viscosity into solids    DO IW=1,NWC      II  = IJKW(1,IW)      JJ  = IJKW(2,IW)      KK  = IJKW(3,IW)      IIG = IJKW(6,IW)      JJG = IJKW(7,IW)      KKG = IJKW(8,IW)      MU(II,JJ,KK) = MU(IIG,JJG,KKG)   ENDDO ENDIF CALC_MU ! Compute vorticity and stress tensor components DO K=0,KBAR   DO J=0,JBAR      DO I=0,IBAR         DUDZ = RDZN(K)*(UU(I,J,K+1)-UU(I,J,K))         DWDX = RDXN(I)*(WW(I+1,J,K)-WW(I,J,K))         OMY(I,J,K) = DUDZ - DWDX         MUY = 0.25_EB*(MU(I+1,J,K)+MU(I,J,K)+MU(I,J,K+1)+MU(I+1,J,K+1))         TXZ(I,J,K) = MUY*(DUDZ + DWDX)      ENDDO   ENDDOENDDO ! Correct vorticity and stress tensor components at solid edges EDGE_LOOP: DO IE=1,N_EDGES   II  = IJKE(1,IE)   JJ  = IJKE(2,IE)   KK  = IJKE(3,IE)   IEC = IJKE(4,IE)   SELECT CASE(IEC)      CASE(2)         OMY(II,JJ,KK) = OME_E(IE)         TXZ(II,JJ,KK) = TAU_E(IE)   END SELECTENDDO EDGE_LOOP ! Compute gravity components GX  = 0._EBGZ  = EVALUATE_RAMP(T,DUMMY,I_RAMP_GZ)*GVEC(3) ! Upwind/Downwind bias factors IF (PREDICTOR) THEN   PMDT =  0.5_EB*DT   MPDT = -0.5_EB*DTELSE   PMDT = -0.5_EB*DT   MPDT =  0.5_EB*DTENDIF ! Compute r-direction flux term FVX IF (XS==0) THEN    I0 = 1ELSE   I0 = 0ENDIF J = 1DO K= 1,KBAR   DO I=I0,IBAR      WP    = WW(I,J,K)   + WW(I+1,J,K)      WM    = WW(I,J,K-1) + WW(I+1,J,K-1)      EPSWP = 1._EB + WP*MPDT*RDZN(K)      EPSWM = 1._EB + WM*PMDT*RDZN(K-1)      WOMY  = EPSWP*WP*OMY(I,J,K) + EPSWM*WM*OMY(I,J,K-1)      RRHO  = 2._EB/(RHOP(I,J,K)+RHOP(I+1,J,K))      AH    = RHO_0(K)*RRHO - 1._EB         DWDZ  = (WW(I+1,J,K)-WW(I+1,J,K-1))*RDZ(K)      TXXP  = MU(I+1,J,K)*( FOTH*DP(I+1,J,K) - 2._EB*DWDZ )      DWDZ  = (WW(I,J,K)-WW(I,J,K-1))*RDZ(K)      TXXM  = MU(I,J,K)  *( FOTH*DP(I,J,K) -2._EB*DWDZ )      DTXXDX= RDXN(I)*(TXXP      -TXXM)      DTXZDZ= RDZ(K) *(TXZ(I,J,K)-TXZ(I,J,K-1))      DMUDX = (MU(I+1,J,K)-MU(I,J,K))*RDXN(I)      VTRM  = RRHO*( DTXXDX + DTXZDZ - 2._EB*UU(I,J,K)*DMUDX/R(I) )       FVX(I,J,K) = 0.25_EB*WOMY + GX*AH - VTRM    ENDDOENDDO    ! Compute z-direction flux term FVZ J = 1DO K=0,KBAR   DO I=1,IBAR      UP    = UU(I,J,K)   + UU(I,J,K+1)      UM    = UU(I-1,J,K) + UU(I-1,J,K+1)      EPSUP = 1._EB + UP*MPDT*RDXN(I)      EPSUM = 1._EB + UM*PMDT*RDXN(I-1)      UOMY  = EPSUP*UP*OMY(I,J,K) + EPSUM*UM*OMY(I-1,J,K)      RRHO  = 2._EB/(RHOP(I,J,K)+RHOP(I,J,K+1))      AH    = 0.5_EB*(RHO_0(K)+RHO_0(K+1))*RRHO - 1._EB      DUDX  = (R(I)*UU(I,J,K+1)-R(I-1)*UU(I-1,J,K+1))*RDX(I)*RRN(I)      TZZP  = MU(I,J,K+1)*( FOTH*DP(I,J,K+1) - 2._EB*DUDX )      DUDX  = (R(I)*UU(I,J,K)-R(I-1)*UU(I-1,J,K))*RDX(I)*RRN(I)      TZZM  = MU(I,J,K)  *( FOTH*DP(I,J,K)   - 2._EB*DUDX )      DTXZDX= RDX(I) *(R(I)*TXZ(I,J,K)-R(I-1)*TXZ(I-1,J,K))*RRN(I)      DTZZDZ= RDZN(K)*(TZZP      -TZZM)      VTRM  = RRHO*(DTXZDX + DTZZDZ)      FVZ(I,J,K) = -0.25_EB*UOMY + GZ*AH - VTRM    ENDDOENDDO    ! Baroclinic torque correction terms IF (BAROCLINIC) CALL BAROCLINIC_CORRECTION ! Adjust FVX and FVZ at solid, internal obstructions for no flux CALL NO_FLUX END SUBROUTINE VELOCITY_FLUX_CYLINDRICAL  SUBROUTINE NO_FLUX! Set FVX,FVY,FVZ inside internal blockages to maintain no fluxUSE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP REAL(EB) :: RFODTINTEGER  :: IC2,IC1,N,I,J,K,IW,II,JJ,KK,IORREAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WWTYPE (OBSTRUCTION_TYPE), POINTER :: OB IF (PREDICTOR) THEN   UU => U   VV => V   WW => WELSE   UU => US   VV => VS   WW => WSENDIF RFODT = RELAXATION_FACTOR/DT ! Drive velocity components inside solid obstructions towards zero OBST_LOOP: DO N=1,N_OBST   OB=>OBSTRUCTION(N)   DO K=OB%K1+1,OB%K2      DO J=OB%J1+1,OB%J2         LOOP1: DO I=OB%I1  ,OB%I2            IC1 = CELL_INDEX(I,J,K)            IC2 = CELL_INDEX(I+1,J,K)            IF (SOLID(IC1) .AND. SOLID(IC2)) FVX(I,J,K) = -RDXN(I)*(H(I+1,J,K)-H(I,J,K)) + RFODT*UU(I,J,K)         ENDDO LOOP1      ENDDO    ENDDO    DO K=OB%K1+1,OB%K2      DO J=OB%J1  ,OB%J2         LOOP2: DO I=OB%I1+1,OB%I2            IC1 = CELL_INDEX(I,J,K)            IC2 = CELL_INDEX(I,J+1,K)            IF (SOLID(IC1) .AND. SOLID(IC2)) FVY(I,J,K) = -RDYN(J)*(H(I,J+1,K)-H(I,J,K)) + RFODT*VV(I,J,K)         ENDDO LOOP2      ENDDO    ENDDO    DO K=OB%K1  ,OB%K2      DO J=OB%J1+1,OB%J2         LOOP3: DO I=OB%I1+1,OB%I2            IC1 = CELL_INDEX(I,J,K)            IC2 = CELL_INDEX(I,J,K+1)            IF (SOLID(IC1) .AND. SOLID(IC2)) FVZ(I,J,K) = -RDZN(K)*(H(I,J,K+1)-H(I,J,K)) + RFODT*WW(I,J,K)         ENDDO LOOP3      ENDDO    ENDDO ENDDO OBST_LOOP ! Add normal velocity to FVX, etc. for surface cells WALL_LOOP: DO IW=1,NWC !!! IF (.NOT.OBSTRUCTION(OBST_INDEX_W(IW))%SAWTOOTH) CYCLE WALL_LOOP  ! Allow velocity flux through surface of no SAWTOOTH obsts.IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==OPEN_BOUNDARY) CYCLE WALL_LOOPIF (BOUNDARY_TYPE(IW)==INTERPOLATED_BOUNDARY) CYCLE WALL_LOOPII  = IJKW(1,IW)JJ  = IJKW(2,IW)KK  = IJKW(3,IW)IOR = IJKW(4,IW)SELECT CASE (BOUNDARY_TYPE(IW))   CASE (SOLID_BOUNDARY,POROUS_BOUNDARY)      SELECT CASE(IOR)         CASE( 1)             FVX(II,JJ,KK)   = -RDXN(II)  *(H(II+1,JJ,KK)-H(II,JJ,KK)) + RFODT*(UU(II,JJ,KK)   + UWS(IW))         CASE(-1)             FVX(II-1,JJ,KK) = -RDXN(II-1)*(H(II,JJ,KK)-H(II-1,JJ,KK)) + RFODT*(UU(II-1,JJ,KK) - UWS(IW))         CASE( 2)             FVY(II,JJ,KK)   = -RDYN(JJ)  *(H(II,JJ+1,KK)-H(II,JJ,KK)) + RFODT*(VV(II,JJ,KK)   + UWS(IW))         CASE(-2)            FVY(II,JJ-1,KK) = -RDYN(JJ-1)*(H(II,JJ,KK)-H(II,JJ-1,KK)) + RFODT*(VV(II,JJ-1,KK) - UWS(IW))         CASE( 3)             FVZ(II,JJ,KK)   = -RDZN(KK)  *(H(II,JJ,KK+1)-H(II,JJ,KK)) + RFODT*(WW(II,JJ,KK)   + UWS(IW))         CASE(-3)             FVZ(II,JJ,KK-1) = -RDZN(KK-1)*(H(II,JJ,KK)-H(II,JJ,KK-1)) + RFODT*(WW(II,JJ,KK-1) - UWS(IW))      END SELECT   CASE (MIRROR_BOUNDARY)      SELECT CASE(IOR)         CASE( 1)            FVX(II  ,JJ,KK) = 0._EB         CASE(-1)            FVX(II-1,JJ,KK) = 0._EB         CASE( 2)            FVY(II  ,JJ,KK) = 0._EB         CASE(-2)            FVY(II,JJ-1,KK) = 0._EB         CASE( 3)            FVZ(II  ,JJ,KK) = 0._EB         CASE(-3)            FVZ(II,JJ,KK-1) = 0._EB      END SELECTEND SELECT  ENDDO WALL_LOOP END SUBROUTINE NO_FLUX  SUBROUTINE VELOCITY_PREDICTOR(T,NM,STOP_STATUS)USE COMP_FUNCTIONS, ONLY: SECOND ! Estimates the velocity components at the next time step REAL(EB), INTENT(IN) :: TREAL(EB) :: TNOW,RHSINTEGER  :: STOP_STATUS,I,J,KINTEGER, INTENT(IN) :: NMIF (SOLID_PHASE_ONLY) RETURN TNOW=SECOND() CALL POINT_TO_MESH(NM)DO K=1,KBAR   DO J=1,JBAR      DO I=0,IBAR         RHS = -FVX(I,J,K) - RDXN(I)*(H(I+1,J,K)-H(I,J,K))         US(I,J,K) = U(I,J,K) + DT*RHS      ENDDO    ENDDO ENDDO DO K=1,KBAR   DO J=0,JBAR      DO I=1,IBAR         RHS = -FVY(I,J,K) - RDYN(J)*(H(I,J+1,K)-H(I,J,K))         VS(I,J,K) = V(I,J,K) + DT*RHS      ENDDO    ENDDO ENDDO  DO K=0,KBAR   DO J=1,JBAR      DO I=1,IBAR         RHS = -FVZ(I,J,K) - RDZN(K)*(H(I,J,K+1)-H(I,J,K))         WS(I,J,K) = W(I,J,K) + DT*RHS      ENDDO    ENDDO ENDDO ! Check the stability criteria DTOLD = DTCALL CHECK_STABILITY(NM) IF (DT < DTINT*1.E-4) THEN  ! The time step has gotten too small. Kill the job.   STOP_STATUS = INSTABILITY_STOP   RETURNENDIF IF (.NOT.CHANGE_TIME_STEP(NM)) CALL VELOCITY_BC(T,NM) TUSED(4,NM)=TUSED(4,NM)+SECOND()-TNOWEND SUBROUTINE VELOCITY_PREDICTOR  SUBROUTINE VELOCITY_CORRECTOR(T,NM)USE COMP_FUNCTIONS, ONLY: SECOND! Correct the velocity componentsREAL(EB), INTENT(IN) :: T REAL(EB) :: TNOW,RHSINTEGER  :: I,J,KINTEGER, INTENT(IN) :: NM IF (SOLID_PHASE_ONLY) RETURNTNOW=SECOND() CALL POINT_TO_MESH(NM)DO K=1,KBAR   DO J=1,JBAR      DO I=0,IBAR         RHS = -FVX(I,J,K) - RDXN(I)*(H(I+1,J,K)-H(I,J,K))         U(I,J,K) = .5_EB*(U(I,J,K) + US(I,J,K) + DT*RHS)      ENDDO    ENDDO ENDDO DO K=1,KBAR   DO J=0,JBAR      DO I=1,IBAR         RHS = -FVY(I,J,K) - RDYN(J)*(H(I,J+1,K)-H(I,J,K))         V(I,J,K) = .5_EB*(V(I,J,K) + VS(I,J,K) + DT*RHS)      ENDDO    ENDDO ENDDO  DO K=0,KBAR   DO J=1,JBAR      DO I=1,IBAR         RHS = -FVZ(I,J,K) - RDZN(K)*(H(I,J,K+1)-H(I,J,K))         W(I,J,K) = .5_EB*(W(I,J,K) + WS(I,J,K) + DT*RHS)      ENDDO    ENDDO ENDDO CALL VELOCITY_BC(T,NM) TUSED(4,NM)=TUSED(4,NM)+SECOND()-TNOWEND SUBROUTINE VELOCITY_CORRECTOR  SUBROUTINE VELOCITY_BC(T,NM)! Assert tangential velocity boundary conditionsUSE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP REAL(EB) :: BC,MUA,T,FVT,UP,UM,VP,VM,WP,WM,DUDY,DUDZ,DVDX,DVDZ,DWDX,DWDY,TSIINTEGER  :: I,J,K,IBC,NOM1,NOM2,IIO1,IIO2,JJO1,JJO2,KKO1,KKO2,NM,IE,II,JJ,KK,IEC,IOR,IWM,IWP,ICMM,ICMP,ICPM,ICPP,ICREAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,U_Y,U_Z,V_X,V_Z,W_X,W_YTYPE (SURFACE_TYPE), POINTER :: SF! Point to the appropriate velocity fieldIF (PREDICTOR) THEN   UU => US   VV => VS   WW => WSELSE   UU => U   VV => V   WW => WENDIF! Set the boundary velocity place holder to some large negative numberIF (CORRECTOR) THEN   UVW_GHOST = -1.E6_EB   U_Y => WORK1   U_Z => WORK2   V_X => WORK3   V_Z => WORK4   W_X => WORK5   W_Y => WORK6   U_Y = -1.E6_EB   U_Z = -1.E6_EB   V_X = -1.E6_EB   V_Z = -1.E6_EB   W_X = -1.E6_EB   W_Y = -1.E6_EBENDIF! Loop over all cell edges and determine the appropriate velocity BCsEDGE_LOOP: DO IE=1,N_EDGES   II   = IJKE( 1,IE)   JJ   = IJKE( 2,IE)   KK   = IJKE( 3,IE)   IEC  = IJKE( 4,IE)   IBC  = IJKE( 5,IE)   IOR  = IJKE( 6,IE)   NOM1 = IJKE( 7,IE)   IIO1 = IJKE( 8,IE)   JJO1 = IJKE( 9,IE)   KKO1 = IJKE(10,IE)   NOM2 = IJKE(11,IE)