SUBROUTINE FLUXI()

      USE LES3D_DATA
      IMPLICIT REAL*8(A-H,O-Z)

      ALLOCATABLE QS(:), FSI(:,:)
      ALLOCATABLE QDIFFX(:), RMU(:), EKCOEF(:)
      ALLOCATABLE DUDX(:), DUDY(:), DUDZ(:)
      ALLOCATABLE DVDX(:), DVDY(:), DVDZ(:)
      ALLOCATABLE DWDX(:), DWDY(:), DWDZ(:)
      ALLOCATABLE DTDX(:), DCDX(:,:)
      ALLOCATABLE DHDX(:), DKDX(:)

      PARAMETER (  R12I = 1.0D0 / 12.0D0,
     >              R6I = 1.0D0 / 6.0D0,
     >            THIRD = 1.0D0 / 3.0D0,
     >             TWO3 = 2.0D0 / 3.0D0 )

      call timer_start(11)
      ALLOCATE( QS(0:IMAX-1), FSI(0:IMAX-1,ND))
      ALLOCATE( QDIFFX(0:IMAX-1), RMU(0:IMAX-1), EKCOEF(0:IMAX-1))
      ALLOCATE( DUDX(0:IMAX-1), DUDY(0:IMAX-1), DUDZ(0:IMAX-1))
      ALLOCATE( DVDX(0:IMAX-1), DVDY(0:IMAX-1), DVDZ(0:IMAX-1))
      ALLOCATE( DWDX(0:IMAX-1), DWDY(0:IMAX-1), DWDZ(0:IMAX-1))
      ALLOCATE( DTDX(0:IMAX-1), DCDX(0:IMAX-1,NSPECI))
      ALLOCATE( DHDX(0:IMAX-1), DKDX(0:IMAX-1))
      call timer_stop(11)
      
      call timer_start(12)
      
      QS=0D0
      FSI=0D0
      QDIFFX=0D0
      RMU=0D0
      EKCOEF=0D0
      DUDX=0D0
      DUDY=0D0
      DUDZ=0D0
      DVDX=0D0
      DVDY=0D0
      DVDZ=0D0
      DWDX=0D0
      DWDY=0D0
      DWDZ=0D0
      DTDX=0D0
      DCDX=0D0
      DHDX=0D0
      DKDX=0D0
      call timer_stop(12)
      call timer_start(13)
      CALL EXTRAPI()
      call timer_stop(13)

      I1 = 0
      I2 = IMAX - 1

      J1 = 1
      J2 = JMAX - 1

      K1 = 1
      K2 = KMAX - 1

      call timer_start(14)

	print*, 'in FLUXI, outside PARALLEL DO'

C$OMP PARALLELDO NUM_THREADS(MIN(K2-K1+1,31))
C$OMP&SHARED (I1,CMU,GDIFFAC,DXI,RDENG,K2,DZI,KK,J2,CNUK,IBDD,XAREA)     
C$OMP&SHARED(K1,I2,DYI,J1,DTVOL,JMAX,DELTA,IADD,N,DCDX,DVDZ,CPK,FSI)     
C$OMP&SHARED(HFK,DWDY,QAV,WAV,CONC,PAV,EK,DU,VAV,HF,UAV,W,V,U,T,Q,H)     
C$OMP&PRIVATE (RHEV,RHEK,TXX,HK,ABD,CPAV,TEMP,SXX,SGSXX,WAVE,RLMBDA)     
C$OMP&PRIVATE (SGSEX,RHAVE,QSPI,VAVE,TZX,TXZ,RK,SGSXY,ICD,II,QSP,QX)     
C$OMP&PRIVATE (UAVE,SYY,SXZ,NS,RD,SGSXZ,TYX,TXY,EKAVE,IBD,SXY,DIV,SZZ)   
C$OMP&PRIVATE (L,J,I,QS,DVDY,DVDX,EKCOEF,DUDZ,DUDY,DUDX,DHDX,QDIFFX)     
C$OMP&PRIVATE (DKDX,RMU,DTDX,DWDZ,DWDX,K)                                
      DO K = K1,K2

         DO J = J1,J2

            ABD = DBLE(IBDD)
!
! EULER STUFF
!
            DO I = I1,I2

               QS(I) = UAV(I,J,K) * XAREA
            END DO
 
            IF(NSCHEME .EQ. 2) THEN
               DO I = I1,I2
                  II = I + 1 - IADD
                  QSP = U(II,J,K) * XAREA
                  QSPI = (QSP - QS(I)) * DBLE((1 - 2 * IADD))
                  IF(QSPI .GT. 0.0D0) QS(I) = 0.5D0 * (QS(I) + QSP)
               END DO
            END IF

            DO L = 1,5
               DO I = I1,I2

                  FSI(I,L) = QAV(I,J,K,L) * QS(I)
               END DO
            END DO

            IF(ISGSK .EQ. 1) THEN
               DO I = I1,I2
                  FSI(I,7) = QAV(I,J,K,7) * QS(I)
               END DO
            END IF

            IF(ICHEM .GT. 0) THEN
               DO NS = 1, NSPECI
                  DO I = I1,I2
                     FSI(I,7+NS) = QAV(I,J,K,7+NS) * QS(I)
                  END DO
               END DO
            END IF
 
            do kk =1,1
            DO I = I1,I2

               FSI(I,2) = FSI(I,2) + PAV(I,J,K) * XAREA
            END DO

            DO I = I1,I2


               FSI(I,5) = FSI(I,5) + PAV(I,J,K) * QS(I)
            END DO
            enddo
!
! VISCOUS STUFF
!
            IF(NSCHEME .EQ. 2) THEN      ! I-DIRECTION DERIVATIVES
 
               DO I = I1,I2
                  DUDX(I) = DXI * (U(I+1,J,K) - U(I,J,K))
                  DVDX(I) = DXI * (V(I+1,J,K) - V(I,J,K))
                  DWDX(I) = DXI * (W(I+1,J,K) - W(I,J,K))
                  DTDX(I) = DXI * (T(I+1,J,K) - T(I,J,K))
               ENDDO

               DO I = I1, I2 
                  DUDY(I) = DYI * (UAV(I,J+1,K) - UAV(I,J-1,K)) * 0.5D0
                  DVDY(I) = DYI * (VAV(I,J+1,K) - VAV(I,J-1,K)) * 0.5D0
                  DWDY(I) = DYI * (WAV(I,J+1,K) - WAV(I,J-1,K)) * 0.5D0
               ENDDO

               DO I = I1, I2 
                  DUDZ(I) = DZI * (UAV(I,J,K+1) - UAV(I,J,K-1)) * 0.5D0
                  DVDZ(I) = DZI * (VAV(I,J,K+1) - VAV(I,J,K-1)) * 0.5D0
                  DWDZ(I) = DZI * (WAV(I,J,K+1) - WAV(I,J,K-1)) * 0.5D0
               END DO
 
               IF(ISGSK .EQ. 1) THEN
                  DO I = I1,I2
                     DKDX(I) = DXI * (EK(I+1,J,K) - EK(I,J,K))
                     DHDX(I) = DXI * ( H(I+1,J,K) -  H(I,J,K))
                  END DO
               END IF
 
               IF(ICHEM .GT. 0) THEN
                  DO NS = 1,NSPECI
                     DO I = I1,I2
                        DCDX(I,NS) = 
     >                   DXI * (CONC(I+1,J,K,NS) - CONC(I,J,K,NS))
                     END DO
                  END DO
               END IF

            ELSE

               DO I = I1,I2


                  II  =  I + IADD
                  IBD = II - IBDD
                  ICD = II + IBDD
   
                  DUDX(I) = 
     >               DXI * ABD * ((U(IBD,J,K) - U(ICD,J,K)) +
     >                    8.0D0 * (U( II,J,K) - U(IBD,J,K))) * R6I
                  DVDX(I) = 
     >               DXI * ABD * ((V(IBD,J,K) - V(ICD,J,K)) +
     >                    8.0D0 * (V( II,J,K) - V(IBD,J,K))) * R6I
                  DWDX(I) = 
     >               DXI * ABD * ((W(IBD,J,K) - W(ICD,J,K)) +
     >                    8.0D0 * (W( II,J,K) - W(IBD,J,K))) * R6I
                  DTDX(I) = 
     >               DXI * ABD * ((T(IBD,J,K) - T(ICD,J,K)) +
     >                    8.0D0 * (T( II,J,K) - T(IBD,J,K))) * R6I
               END DO
 
               IF(JPERIOD .EQ. 1 .AND.
     >           (J .EQ. 1 .OR. J .EQ. JMAX-1)) THEN
                  DO I = I1,I2
                     DUDY(I) = DYI * (UAV(I,J+1,K) - UAV(I,J-1,K))*0.5D0
                     DVDY(I) = DYI * (VAV(I,J+1,K) - VAV(I,J-1,K))*0.5D0
                     DWDY(I) = DYI * (WAV(I,J+1,K) - WAV(I,J-1,K))*0.5D0
                  END DO
               ELSE
                  DO I = I1,I2


                     DUDY(I) =
     >                    DYI * R12I * (UAV(I,J-2,K) - UAV(I,J+2,K) +
     >                         8.0D0 * (UAV(I,J+1,K) - UAV(I,J-1,K)))
                     DVDY(I) =
     >                    DYI * R12I * (VAV(I,J-2,K) - VAV(I,J+2,K) +
     >                         8.0D0 * (VAV(I,J+1,K) - VAV(I,J-1,K)))
                     DWDY(I) =
     >                    DYI * R12I * (WAV(I,J-2,K) - WAV(I,J+2,K) +
     >                         8.0D0 * (WAV(I,J+1,K) - WAV(I,J-1,K)))
                  END DO
               END IF
 
               DO I = I1,I2


                  DUDZ(I) = DZI * R12I * (UAV(I,J,K-2) - UAV(I,J,K+2) +
     >                           8.0D0 * (UAV(I,J,K+1) - UAV(I,J,K-1)))
                  DVDZ(I) = DZI * R12I * (VAV(I,J,K-2) - VAV(I,J,K+2) +
     >                           8.0D0 * (VAV(I,J,K+1) - VAV(I,J,K-1)))
                  DWDZ(I) = DZI * R12I * (WAV(I,J,K-2) - WAV(I,J,K+2) +
     >                           8.0D0 * (WAV(I,J,K+1) - WAV(I,J,K-1)))
               END DO
 
               IF(ISGSK .EQ. 1) THEN
                  DO I = I1,I2
                     II  =  I + IADD
                     IBD = II - IBDD
                     ICD = II + IBDD
 
                     DKDX(I) = 
     >                 DXI * ABD * ((EK(IBD,J,K) - EK(ICD,J,K)) +
     >                      8.0D0 * (EK(II,J,K)  - EK(IBD,J,K))) * R6I
                     DHDX(I) = 
     >                 DXI * ABD * (( H(IBD,J,K) -  H(ICD,J,K)) +
     >                      8.0D0 * ( H(II,J,K)  -  H(IBD,J,K))) * R6I
                  END DO
               END IF

               IF(ICHEM .GT. 0) THEN
                  DO NS = 1,NSPECI
                     DO I = I1,I2
                        II  =  I + IADD
                        IBD = II - IBDD
                        ICD = II + IBDD
 
                        DCDX(I,NS) = DXI *
     >              ABD * ((CONC(IBD,J,K,NS) - CONC(ICD,J,K,NS)) +
     >            8.0D0 *  (CONC( II,J,K,NS) - CONC(IBD,J,K,NS))) * R6I
                     END DO
                  END DO
               END IF

            END IF

            DO I = I1,I2


               QDIFFX(I) = 0.0D0
            END DO
 
            DO I = I1,I2
               TEMP   = 0.5D0 * (T(I,J,K) + T(I+1,J,K))
               RMU(I) = CMU * TEMP * SQRT(TEMP) / (TEMP + 110.0D0)
            END DO
 
            IF(ICHEM .EQ. 2) THEN
               DO NS = 1,NSPECI
                  DO I = I1,I2
                     TEMP = 0.5D0 * (T(I,J,K) + T(I+1,J,K))
                     RD   = RMU(I) / PRANDLT
                     HK   = RD * (CPK(NS) * TEMP + HFK(NS))
                     QDIFFX(I) = QDIFFX(I) + HK * DCDX(I,NS)
                  END DO
               END DO
            END IF
 
            DO I = I1,I2


               UAVE = 0.5D0 * (U(I,J,K)    + U(I+1,J,K))
               VAVE = 0.5D0 * (V(I,J,K)    + V(I+1,J,K))
               WAVE = 0.5D0 * (W(I,J,K)    + W(I+1,J,K))
               CPAV = 0.5D0 * (HF(I,J,K,2) + HF(I+1,J,K,2))
 
               DIV    =  DUDX(I) + DVDY(I) + DWDZ(I)
               RK     =  CPAV * RMU(I) / PRANDLT
               RLMBDA = -TWO3 * RMU(I)
 
               TXX = -2.0D0 * RMU(I) * DUDX(I) - RLMBDA * DIV
               TXY = -RMU(I) * (DUDY(I) + DVDX(I))
               TXZ = -RMU(I) * (DWDX(I) + DUDZ(I))
               TYX =  TXY
               TZX =  TXZ
 
               FSI(I,2) = FSI(I,2) + TXX * XAREA
               FSI(I,3) = FSI(I,3) + TYX * XAREA
               FSI(I,4) = FSI(I,4) + TZX * XAREA
 
               QX = -RK * DTDX(I) - QDIFFX(I)
 
               FSI(I,5) = FSI(I,5) +
     >         (TXX * UAVE + TXY * VAVE + TXZ * WAVE + QX) * XAREA
            END DO

            IF(ICHEM .EQ. 1) THEN
               DO I = I1,I2
                  RD = RMU(I) * GDIFFAC
                  FSI(I,8) = FSI(I,8) - RD * DCDX(I,1) * XAREA
               END DO
            ELSE IF(ICHEM .EQ. 2) THEN
               DO NS = 1,NSPECI
                  DO I = I1,I2
                     RD = RMU(I) / PRANDLT
                     FSI(I,7+NS) =
     >                     FSI(I,7+NS) - RD * DCDX(I,NS) * XAREA
                  END DO
               END DO
            END IF
 
            IF(ISGSK .EQ. 1) THEN
               DO I = I1,I2
                  RHAVE  = 0.5D0 * (Q(I,J,K,1,N) + Q(I+1,J,K,1,N))
                  EKAVE  = 0.5D0 * (EK(I,J,K)    + EK(I+1,J,K))
!                  CNUAVE = 0.5D0 * (CNU(I,J,K)   + CNU(I+1,J,K))
!                  CENAVE = 0.5D0 * (CEN(I,J,K)   + CEN(I+1,J,K))
 
!                  EKCOEF(I) = RHAVE * CNUAVE * SQRT(EKAVE) * DELTA
                  EKCOEF(I) = RHAVE * CNUK * SQRT(EKAVE) * DELTA
                  RHEV  = 2.0D0 * EKCOEF(I)
!                  RDENG = CENAVE * EKCOEF(I)
                  RHEK  = TWO3 * RHAVE * EKAVE
 
                  SXX = DUDX(I)
                  SYY = DVDY(I)
                  SZZ = DWDZ(I)
                  SXY = 0.5D0 * (DUDY(I) + DVDX(I))
                  SXZ = 0.5D0 * (DUDZ(I) + DWDX(I))
 
                  DIV  = (SXX + SYY + SZZ) * THIRD

                  SGSXX = - RHEV * (SXX - DIV) + RHEK
                  SGSXY = - RHEV * SXY
                  SGSXZ = - RHEV * SXZ
 
                  SGSEX = - RDENG * DHDX(I)
 
                  FSI(I,2) = FSI(I,2) + SGSXX * XAREA
                  FSI(I,3) = FSI(I,3) + SGSXY * XAREA
                  FSI(I,4) = FSI(I,4) + SGSXZ * XAREA
                  FSI(I,5) = FSI(I,5) + SGSEX * XAREA
               END DO
 
               DO I = I1,I2
!                  CENAVE = 0.5D0 * (CEN(I,J,K) + CEN(I+1,J,K))
!                  RDENG  = CENAVE * EKCOEF(I) + RMU(I) / PRANDLT
                  RDENG = EKCOEF(I) + RMU(I) / PRANDLT
                  FSI(I,7) = FSI(I,7) - RDENG * DKDX(I) * XAREA
               END DO
 
               IF(ICHEM .GT. 0) THEN
                  DO NS = 1,NSPECI
                     DO I = I1,I2
                        FSI(I,7+NS) = FSI(I,7+NS) -
     >                       EKCOEF(I) * DCDX(I,NS) * XAREA
                     END DO
                  END DO
               END IF

            END IF
!
! DF/DX = (DT / VOL) * (F_I+1 - F_I)
!
            DO L = 1, 5 
               DO I = I1+1, I2


                  DU(I,J,K,L) = -DTVOL * (FSI(I,L) - FSI(I-1,L))
               END DO
            END DO

            IF(ISGSK .EQ. 1) THEN
               DO I = I1+1, I2
                  DU(I,J,K,7) = -DTVOL * (FSI(I,7) - FSI(I-1,7))
               END DO
            END IF

            IF(ICHEM .GT. 0) THEN
               DO NS = 1, NSPECI
                  DO I = I1+1, I2
                     DU(I,J,K,7+NS) =
     >                     -DTVOL * (FSI(I,7+NS) - FSI(I-1,7+NS))
                  END DO
               END DO
            END IF

         END DO		! DO J

      END DO		! DO K

      DEALLOCATE( QS, FSI,
     >            QDIFFX, RMU, EKCOEF,
     >            DUDX, DUDY, DUDZ,
     >            DVDX, DVDY, DVDZ,
     >            DWDX, DWDY, DWDZ,
     >            DTDX, DCDX,
     >            DHDX, DKDX)

	print*, ' in FLUXI, after DEALLOCATE'

      RETURN
      END
!------------------------------------------------------------------------