MODULE VELO! Module computes the velocity flux terms, baroclinic torque correction terms, and performs the CFL CheckUSE PRECISION_PARAMETERSUSE GLOBAL_CONSTANTSUSE MESH_POINTERSIMPLICIT NONEPRIVATECHARACTER(255), PARAMETER :: veloid='$Id: velo.f90 645 2007-09-19 20:15:22Z drjfloyd $'CHARACTER(255), PARAMETER :: velorev='$Revision: 645 $'CHARACTER(255), PARAMETER :: velodate='$Date: 2007-09-19 16:15:22 -0400 (Wed, 19 Sep 2007) $'PUBLIC COMPUTE_VELOCITY_FLUX,VELOCITY_PREDICTOR,VELOCITY_CORRECTOR,NO_FLUX,GET_REV_veloPRIVATE VELOCITY_FLUX,VELOCITY_FLUX_ISOTHERMAL,VELOCITY_FLUX_CYLINDRICAL CONTAINS SUBROUTINE COMPUTE_VELOCITY_FLUX(T,NM)USE COMP_FUNCTIONS, ONLY: SECONDREAL(EB), INTENT(IN) :: TREAL(EB) :: TNOWINTEGER, INTENT(IN) :: NMIF (SOLID_PHASE_ONLY) RETURNTNOW = SECOND()IF (DNS .AND. ISOTHERMAL .AND. N_SPECIES==0) THEN   CALL VELOCITY_FLUX_ISOTHERMAL(NM)ELSE   IF (.NOT.CYLINDRICAL) CALL VELOCITY_FLUX(T,NM)   IF (     CYLINDRICAL) CALL VELOCITY_FLUX_CYLINDRICAL(T,NM)ENDIFTUSED(4,NM) = TUSED(4,NM) + SECOND() - TNOWEND SUBROUTINE COMPUTE_VELOCITY_FLUXSUBROUTINE VELOCITY_FLUX(T,NM)! Compute convective and diffusive termsUSE PHYSICAL_FUNCTIONS, ONLY: GET_MU2USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMPINTEGER, INTENT(IN) :: NMREAL(EB) :: T,MUX,MUY,MUZ,UP,UM,VP,VM,WP,WM,VTRM, &            DTXYDY,DTXZDZ,DTYZDZ,DTXYDX,DTXZDX,DTYZDY, &            DUDX,DVDY,DWDZ,DUDY,DUDZ,DVDX,DVDZ,DWDX,DWDY, &            VOMZ,WOMY,UOMY,VOMX,UOMZ,WOMX,PMDT,MPDT, &            S2,C,CDXDYDZTT,AH,RRHO,GX,GY,GZ,Z_2, &            TXXP,TXXM,TYYP,TYYM,TZZP,TZZM,DTXXDX,DTYYDY,DTZZDZ, &            EPSUP,EPSUM,EPSVP,EPSVM,EPSWP,EPSWM,MU_SUM,DUMMY=0._EBINTEGER :: II,JJ,KK,I,J,K,IW,IIG,JJG,KKG,ITMP,N,IE,IECREAL(EB), POINTER, DIMENSION(:,:,:) :: TXY,TXZ,TYZ,OMX,OMY,OMZ, UU,VV,WW,RHOP,DPREAL(EB), POINTER, DIMENSION(:,:,:,:) :: YYP CALL POINT_TO_MESH(NM) IF (PREDICTOR) THEN   UU => U   VV => V   WW => W   DP => D     RHOP => RHO   IF (N_SPECIES > 0) YYP => YYELSE   UU => US   VV => VS   WW => WS   DP => DS   RHOP => RHOS   IF (N_SPECIES > 0) YYP => YYSENDIFTXY => WORK1TXZ => WORK2TYZ => WORK3OMX => WORK4OMY => WORK5OMZ => WORK6 CALC_MU: IF (PREDICTOR) THEN        LES_VS_DNS: IF (LES) THEN   ! Smagorinsky model (LES)      C = CSMAG**2      IF (EVACUATION_ONLY(NM)) C = (0.9_EB)**2      KLOOP: DO K=1,KBAR         JLOOP: DO J=1,JBAR            ILOOP: DO I=1,IBAR               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE ILOOP               IF (.NOT.TWO_D) THEN                  CDXDYDZTT = C*(DX(I)*DY(J)*DZ(K))**TWTH               ELSE                  CDXDYDZTT = C*DX(I)*DZ(K)               ENDIF               DUDX = RDX(I)*(UU(I,J,K)-UU(I-1,J,K))               DVDY = RDY(J)*(VV(I,J,K)-VV(I,J-1,K))               DWDZ = RDZ(K)*(WW(I,J,K)-WW(I,J,K-1))               DUDY = 0.25_EB*RDY(J)*(UU(I,J+1,K)-UU(I,J-1,K)+UU(I-1,J+1,K)-UU(I-1,J-1,K))               DVDX = 0.25_EB*RDX(I)*(VV(I+1,J,K)-VV(I-1,J,K)+VV(I+1,J-1,K)-VV(I-1,J-1,K))               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))               DWDY = 0.25_EB*RDY(J)*(WW(I,J+1,K)-WW(I,J-1,K)+WW(I,J+1,K-1)-WW(I,J-1,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))                DVDZ = 0.25_EB*RDZ(K)*(VV(I,J,K+1)-VV(I,J,K-1)+VV(I,J-1,K+1)-VV(I,J-1,K-1))                S2   = 2._EB*(DUDX*DUDX + DVDY*DVDY + DWDZ*DWDZ ) + (DUDY+DVDX)**2 + (DUDZ+DWDX)**2 + &                      (DVDZ+DWDY)**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 JLOOP      ENDDO KLOOP    ELSE LES_VS_DNS ! DNS viscosity       MIXTURE_FRACTION_DNS: IF (.NOT.MIXTURE_FRACTION) THEN         DO K=1,KBAR            DO J=1,JBAR               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         ENDDO      ELSE MIXTURE_FRACTION_DNS         Z_2 = 0._EB         DO K=1,KBAR            DO J=1,JBAR               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),Z_2,YY(I,J,K,I_PROG_F),Y_SUM(I,J,K),MU(I,J,K),ITMP)                                    ENDIF                  MU_SUM = MU(I,J,K)*(1-Y_SUM(I,J,K))                  DO N=1,N_SPECIES                     IF(SPECIES(N)%MODE/=MIXTURE_FRACTION_SPECIES) MU_SUM = MU_SUM + YYP(I,J,K,N)*(SPECIES(N)%MU(ITMP)-MU(I,J,K))                  ENDDO                  MU(I,J,K) = MU_SUM               ENDDO IVLOOP2            ENDDO         ENDDO      ENDIF MIXTURE_FRACTION_DNS    ENDIF LES_VS_DNS ! Mirror viscosity into solids    WALL_LOOP: DO IW=1,NWC      IF (BOUNDARY_TYPE(IW)==NULL_BOUNDARY .OR. BOUNDARY_TYPE(IW)==POROUS_BOUNDARY) CYCLE WALL_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)      MU(II,JJ,KK) = MU(IIG,JJG,KKG)   ENDDO WALL_LOOP       MU(0,0:JBP1,0)       = MU(   1,0:JBP1,1)   MU(IBP1,0:JBP1,0)    = MU(IBAR,0:JBP1,1)   MU(IBP1,0:JBP1,KBP1) = MU(IBAR,0:JBP1,KBAR)   MU(0,0:JBP1,KBP1)    = MU(   1,0:JBP1,KBAR)   MU(0:IBP1,0,0)       = MU(0:IBP1,   1,1)   MU(0:IBP1,JBP1,0)    = MU(0:IBP1,JBAR,1)   MU(0:IBP1,JBP1,KBP1) = MU(0:IBP1,JBAR,KBAR)   MU(0:IBP1,0,KBP1)    = MU(0:IBP1,   1,KBAR)   MU(0,0,0:KBP1)       = MU(   1,   1,0:KBP1)   MU(IBP1,0,0:KBP1)    = MU(IBAR,   1,0:KBP1)   MU(IBP1,JBP1,0:KBP1) = MU(IBAR,JBAR,0:KBP1)   MU(0,JBP1,0:KBP1)    = MU(   1,JBAR,0:KBP1) ENDIF CALC_MU ! Compute vorticity and stress tensor componentsDO K=0,KBAR   DO J=0,JBAR      DO I=0,IBAR         DUDY = RDYN(J)*(UU(I,J+1,K)-UU(I,J,K))         DVDX = RDXN(I)*(VV(I+1,J,K)-VV(I,J,K))         DUDZ = RDZN(K)*(UU(I,J,K+1)-UU(I,J,K))         DWDX = RDXN(I)*(WW(I+1,J,K)-WW(I,J,K))         DVDZ = RDZN(K)*(VV(I,J,K+1)-VV(I,J,K))         DWDY = RDYN(J)*(WW(I,J+1,K)-WW(I,J,K))         OMX(I,J,K) = DWDY - DVDZ         OMY(I,J,K) = DUDZ - DWDX         OMZ(I,J,K) = DVDX - DUDY         MUX = 0.25_EB*(MU(I,J+1,K)+MU(I,J,K)+MU(I,J,K+1)+MU(I,J+1,K+1))         MUY = 0.25_EB*(MU(I+1,J,K)+MU(I,J,K)+MU(I,J,K+1)+MU(I+1,J,K+1))         MUZ = 0.25_EB*(MU(I+1,J,K)+MU(I,J,K)+MU(I,J+1,K)+MU(I+1,J+1,K))         TXY(I,J,K) = MUZ*(DVDX + DUDY)         TXZ(I,J,K) = MUY*(DUDZ + DWDX)         TYZ(I,J,K) = MUX*(DVDZ + DWDY)      ENDDO   ENDDOENDDO! Correct vorticity and stress tensor components at solid edgesEDGE_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(1)         OMX(II,JJ,KK) = OME_E(IE)         TYZ(II,JJ,KK) = TAU_E(IE)      CASE(2)         OMY(II,JJ,KK) = OME_E(IE)         TXZ(II,JJ,KK) = TAU_E(IE)      CASE(3)         OMZ(II,JJ,KK) = OME_E(IE)         TXY(II,JJ,KK) = TAU_E(IE)   END SELECTENDDO EDGE_LOOP! Compute gravity components GX  = EVALUATE_RAMP(T,DUMMY,I_RAMP_GX)*GVEC(1)GY  = EVALUATE_RAMP(T,DUMMY,I_RAMP_GY)*GVEC(2)GZ  = 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 x-direction flux term FVX DO K=1,KBAR   DO J=1,JBAR      DO I=0,IBAR         WP    = WW(I,J,K)   + WW(I+1,J,K)         WM    = WW(I,J,K-1) + WW(I+1,J,K-1)         VP    = VV(I,J,K)   + VV(I+1,J,K)         VM    = VV(I,J-1,K) + VV(I+1,J-1,K)         EPSWP = 1._EB + WP*MPDT*RDZN(K)         EPSWM = 1._EB + WM*PMDT*RDZN(K-1)         EPSVP = 1._EB + VP*MPDT*RDYN(J)         EPSVM = 1._EB + VM*PMDT*RDYN(J-1)         WOMY  = EPSWP*WP*OMY(I,J,K) + EPSWM*WM*OMY(I,J,K-1)         VOMZ  = EPSVP*VP*OMZ(I,J,K) + EPSVM*VM*OMZ(I,J-1,K)         RRHO  = 2._EB/(RHOP(I,J,K)+RHOP(I+1,J,K))         AH    = RHO_0(K)*RRHO - 1._EB            DVDY  = (VV(I+1,J,K)-VV(I+1,J-1,K))*RDY(J)         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*(DVDY+DWDZ) )         DVDY  = (VV(I,J,K)-VV(I,J-1,K))*RDY(J)         DWDZ  = (WW(I,J,K)-WW(I,J,K-1))*RDZ(K)         TXXM  = MU(I,J,K)  *( FOTH*DP(I,J,K)   - 2._EB*(DVDY+DWDZ) )         DTXXDX= RDXN(I)*(TXXP      -TXXM)         DTXYDY= RDY(J) *(TXY(I,J,K)-TXY(I,J-1,K))         DTXZDZ= RDZ(K) *(TXZ(I,J,K)-TXZ(I,J,K-1))         VTRM  = RRHO*(DTXXDX + DTXYDY + DTXZDZ)         FVX(I,J,K) = 0.25_EB*(WOMY - VOMZ) + GX*AH - VTRM       ENDDO    ENDDO   ENDDO    ! Compute y-direction flux term FVYDO K=1,KBAR   DO J=0,JBAR      DO I=1,IBAR         UP    = UU(I,J,K)   + UU(I,J+1,K)         UM    = UU(I-1,J,K) + UU(I-1,J+1,K)         WP    = WW(I,J,K)   + WW(I,J+1,K)         WM    = WW(I,J,K-1) + WW(I,J+1,K-1)         EPSUP = 1._EB + UP*MPDT*RDXN(I)         EPSUM = 1._EB + UM*PMDT*RDXN(I-1)         EPSWP = 1._EB + WP*MPDT*RDZN(K)         EPSWM = 1._EB + WM*PMDT*RDZN(K-1)         WOMX  = EPSWP*WP*OMX(I,J,K) + EPSWM*WM*OMX(I,J,K-1)         UOMZ  = EPSUP*UP*OMZ(I,J,K) + EPSUM*UM*OMZ(I-1,J,K)         RRHO  = 2._EB/(RHOP(I,J,K)+RHOP(I,J+1,K))         AH    = RHO_0(K)*RRHO - 1._EB         DUDX  = (UU(I,J+1,K)-UU(I-1,J+1,K))*RDX(I)         DWDZ  = (WW(I,J+1,K)-WW(I,J+1,K-1))*RDZ(K)         TYYP  = MU(I,J+1,K)*( FOTH*DP(I,J+1,K) - 2._EB*(DUDX+DWDZ) )         DUDX  = (UU(I,J,K)-UU(I-1,J,K))*RDX(I)         DWDZ  = (WW(I,J,K)-WW(I,J,K-1))*RDZ(K)         TYYM  = MU(I,J,K)  *( FOTH*DP(I,J,K)   - 2._EB*(DUDX+DWDZ) )         DTXYDX= RDX(I) *(TXY(I,J,K)-TXY(I-1,J,K))         DTYYDY= RDYN(J)*(TYYP      -TYYM)         DTYZDZ= RDZ(K) *(TYZ(I,J,K)-TYZ(I,J,K-1))         VTRM  = RRHO*(DTXYDX + DTYYDY + DTYZDZ)         FVY(I,J,K) = 0.25_EB*(UOMZ - WOMX) + GY*AH - VTRM       ENDDO   ENDDO   ENDDO    ! Compute z-direction flux term FVZ DO K=0,KBAR   DO J=1,JBAR      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)         VP    = VV(I,J,K)   + VV(I,J,K+1)         VM    = VV(I,J-1,K) + VV(I,J-1,K+1)         EPSUP = 1._EB + UP*MPDT*RDXN(I)         EPSUM = 1._EB + UM*PMDT*RDXN(I-1)         EPSVP = 1._EB + VP*MPDT*RDYN(J)         EPSVM = 1._EB + VM*PMDT*RDYN(J-1)         UOMY  = EPSUP*UP*OMY(I,J,K) + EPSUM*UM*OMY(I-1,J,K)         VOMX  = EPSVP*VP*OMX(I,J,K) + EPSVM*VM*OMX(I,J-1,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  = (UU(I,J,K+1)-UU(I-1,J,K+1))*RDX(I)         DVDY  = (VV(I,J,K+1)-VV(I,J-1,K+1))*RDY(J)         TZZP  = MU(I,J,K+1)*( FOTH*DP(I,J,K+1) - 2._EB*(DUDX+DVDY) )         DUDX  = (UU(I,J,K)-UU(I-1,J,K))*RDX(I)         DVDY  = (VV(I,J,K)-VV(I,J-1,K))*RDY(J)         TZZM  = MU(I,J,K)  *( FOTH*DP(I,J,K)   - 2._EB*(DUDX+DVDY) )         DTXZDX= RDX(I) *(TXZ(I,J,K)-TXZ(I-1,J,K))         DTYZDY= RDY(J) *(TYZ(I,J,K)-TYZ(I,J-1,K))         DTZZDZ= RDZN(K)*(TZZP      -TZZM)         VTRM  = RRHO*(DTXZDX + DTYZDY + DTZZDZ)         FVZ(I,J,K) = 0.25_EB*(VOMX - UOMY) + GZ*AH - VTRM             ENDDO   ENDDO   ENDDO    ! Baroclinic torque correction IF (BAROCLINIC) CALL BAROCLINIC_CORRECTION ! Adjust FVX, FVY and FVZ at solid, internal obstructions for no flux CALL NO_FLUX END SUBROUTINE VELOCITY_FLUX  SUBROUTINE VELOCITY_FLUX_ISOTHERMAL(NM) ! Compute the velocity flux at cell edges (ISOTHERMAL DNS ONLY) REAL(EB) :: UP,UM,VP,VM,WP,WM,VTRM,VOMZ,WOMY,UOMY,VOMX,UOMZ,WOMX, &            DVDZ,DVDX,DWDY,DWDX,DUDZ,DUDY,PMDT,MPDT, &            EPSUP,EPSUM,EPSVP,EPSVM,EPSWP,EPSWMINTEGER :: I,J,K,II,JJ,KK,IE,IECREAL(EB), POINTER, DIMENSION(:,:,:) :: OMX,OMY,OMZ,UU,VV,WWINTEGER, INTENT(IN) :: NM CALL POINT_TO_MESH(NM) IF (PREDICTOR) THEN   UU => U   VV => V   WW => WELSE   UU => US   VV => VS   WW => WSENDIF OMX => WORK4OMY => WORK5OMZ => WORK6 ! Compute vorticity and stress tensor components DO K=0,KBAR   DO J=0,JBAR      DO I=0,IBAR         DUDY = RDYN(J)*(UU(I,J+1,K)-UU(I,J,K))         DVDX = RDXN(I)*(VV(I+1,J,K)-VV(I,J,K))         DUDZ = RDZN(K)*(UU(I,J,K+1)-UU(I,J,K))         DWDX = RDXN(I)*(WW(I+1,J,K)-WW(I,J,K))         DVDZ = RDZN(K)*(VV(I,J,K+1)-VV(I,J,K))         DWDY = RDYN(J)*(WW(I,J+1,K)-WW(I,J,K))         OMX(I,J,K) = DWDY - DVDZ         OMY(I,J,K) = DUDZ - DWDX         OMZ(I,J,K) = DVDX - DUDY      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(1)         OMX(II,JJ,KK) = OME_E(IE)      CASE(2)         OMY(II,JJ,KK) = OME_E(IE)      CASE(3)         OMZ(II,JJ,KK) = OME_E(IE)   END SELECTENDDO EDGE_LOOP ! 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 x-direction flux term FVX DO K=1,KBAR   DO J=1,JBAR      DO I=0,IBAR         WP    = WW(I,J,K)   + WW(I+1,J,K)         WM    = WW(I,J,K-1) + WW(I+1,J,K-1)         VP    = VV(I,J,K)   + VV(I+1,J,K)         VM    = VV(I,J-1,K) + VV(I+1,J-1,K)         EPSWP = 1._EB + WP*MPDT*RDZN(K)         EPSWM = 1._EB + WM*PMDT*RDZN(K-1)         EPSVP = 1._EB + VP*MPDT*RDYN(J)         EPSVM = 1._EB + VM*PMDT*RDYN(J-1)         WOMY  = EPSWP*WP*OMY(I,J,K) + EPSWM*WM*OMY(I,J,K-1)         VOMZ  = EPSVP*VP*OMZ(I,J,K) + EPSVM*VM*OMZ(I,J-1,K)         VTRM  = RREDZ(K)*(OMY(I,J,K)-OMY(I,J,K-1)) - RREDY(J)*(OMZ(I,J,K)-OMZ(I,J-1,K))         FVX(I,J,K) = 0.25_EB*(WOMY - VOMZ) - VTRM      ENDDO   ENDDOENDDO! Compute y-direction flux term FVY DO K=1,KBAR   DO J=0,JBAR      DO I=1,IBAR         UP    = UU(I,J,K)   + UU(I,J+1,K)         UM    = UU(I-1,J,K) + UU(I-1,J+1,K)         WP    = WW(I,J,K)   + WW(I,J+1,K)         WM    = WW(I,J,K-1) + WW(I,J+1,K-1)         EPSUP = 1._EB + UP*MPDT*RDXN(I)         EPSUM = 1._EB + UM*PMDT*RDXN(I-1)         EPSWP = 1._EB + WP*MPDT*RDZN(K)         EPSWM = 1._EB + WM*PMDT*RDZN(K-1)         WOMX  = EPSWP*WP*OMX(I,J,K) + EPSWM*WM*OMX(I,J,K-1)         UOMZ  = EPSUP*UP*OMZ(I,J,K) + EPSUM*UM*OMZ(I-1,J,K)         VTRM  = RREDX(I)*(OMZ(I,J,K)-OMZ(I-1,J,K)) - RREDZ(K)*(OMX(I,J,K)-OMX(I,J,K-1))         FVY(I,J,K) = 0.25_EB*(UOMZ - WOMX) - VTRM      ENDDO   ENDDOENDDO ! Compute z-direction flux term FVZ DO K=0,KBAR   DO J=1,JBAR      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)         VP    = VV(I,J,K)   + VV(I,J,K+1)         VM    = VV(I,J-1,K) + VV(I,J-1,K+1)         EPSUP = 1._EB + UP*MPDT*RDXN(I)         EPSUM = 1._EB + UM*PMDT*RDXN(I-1)         EPSVP = 1._EB + VP*MPDT*RDYN(J)         EPSVM = 1._EB + VM*PMDT*RDYN(J-1)         UOMY  = EPSUP*UP*OMY(I,J,K) + EPSUM*UM*OMY(I-1,J,K)         VOMX  = EPSVP*VP*OMX(I,J,K) + EPSVM*VM*OMX(I,J-1,K)         VTRM  = RREDY(J)*(OMX(I,J,K)-OMX(I,J-1,K)) - RREDX(I)*(OMY(I,J,K)-OMY(I-1,J,K))         FVZ(I,J,K) = 0.25_EB*(VOMX - UOMY) - VTRM      ENDDO   ENDDO