!##############################################################
      PROGRAM PCOL
!##############################################################
!     AUTHOR: A. BADARUDIN, UM, 2011.
!
!     TECH LOG:
!     14/07/11 - Inlet, Outlet.
!                search 'INLET','OUTLET','PROFILE'.
!     15/07/11 - Grid stretching 'STRETCH', profile 'PROFILE'.
!     17/07/11 - LUD 'LUD'.
!     14/08/11 - Temp Eqn 'T'.
!     KOMEGA - H:\UMResearch\simple_structured\YinJen\YJLee_
!     06062011\sduct_stretchedgrid\3d_simple_nonortho_duct_turb
!     17/06/16 - Bypassing a block of cells. 'BLK'.
!     02/08/16 - Input file for students 'STUDIN', no need to
!                open source file for input.
!     22/06/16 - k-epsilon turbulence model. 'TURB'
!     15/03/16 - No MPI (serial).
!     19/03/16 - LUD update 'LUD'.
!
!     TESTED: //, CAD/M Server. See cgbpar.out,cgbpar.xls
!     TESTER:
!     COMPILE AND LINK:
!     gfortran -o cgbpar cgbpar.f90 global.f90
!
! BLK:
!##############################################################
!     PROGRAM PCOL
!##############################################################
!     TNIBK=2
!     NIMBK=TNIBK+1
!     XBK2X(1)=1
!     XBK2X(2)=5
!     XBK2X)NIMBK)=NIM
! BLK:
!     DO IBK=1,NIMBK
!       XBK(1)=X(XBK2X(1))
!     END DO
! BLK:
!     DO KBK=1,NKMBK
!     DO IBK=1,NIMBK
!     DO JBK=1,NJMBK
!       IJKBK=(KBK-1)*NIJBK+(IBK-1)*NJBK+J-1
!       BKFLG(IJKBK)=.TRUE.
!     END DO
!     END DO
!     END DO
!     BKFLG(0)=.FALSE.
! BLK:
!###########################################################
!     SUBROUTINE CALCUVW
!###########################################################
!     DO KBK=1,NKMBK
!     DO IBK=1,NIMBK
!     DO JBK=1,NJMBK
!     IJKBK=(KBK-1)*NIJBK+(IBK-1)*NJBK+J-1
!     IF (BKFLG(IJKBK).EQ.TRUE) THEN
!     DO K=BZS(IJKBK),BZE(IJKBK)
!       DZ=Z(K)-Z(K-1)
!       DO I=BIS(IJKBK),BIE(IJKBK)
!         DX=X(I)-X(I-1)
!         DO J=BYS(IJKBK),BYE(IJKBK)
!           DY=Y(J)-Y(J-1)
!           IJK=LK(K)+LI(I)+J
!
!.....CELL-FACE PRESSURE, CELL-CENTER GRADIENT, SOURCE
! 3D
!         IF (IBK.LT.NIMBK) THEN
!           IF (BKFLG(IJKBK+TNIBK).EQ..FALSE..AND.&
!             &I.EQ.BIE(IJKBK)) THEN
!             IK=(K-1)*NI+I
!             PE=PEBK(IBK,IK)!PW=PWBK(IBK,IK) FOR WEST BNDY PRESS.
!           ELSE
!             PE=P(IJK+NJ)*FX(I)+P(IJK)*(1.-FX(I))
!           END IF
!         ELSE
!           PE=P(IJK+NJ)*FX(I)+P(IJK)*(1.-FX(I))
!         END IF
!         END DO
!       END DO
!     END DO
!     END IF
!     END DO
!     END DO
!     END DO
! BLK:
!#############################################################
!     SUBROUTINE PBOUND(FI,FIS1,FIB2)
!#############################################################
!     DO KBK=1,NKMBK
!     DO IBK=1,NIMBK
!     DO JBK=1,NJMBK
!     IJKBK=(KBK-1)*NIJBK+(IBK-1)*NJBK+J-1
!     IF (BKFLG(IJKBK).EQ.TRUE) THEN
!     DO K=BZS(IJKBK),BZE(IJKBK)
!       DO I=BIS(IJKBK),BIE(IJKBK)
!       IF (IBK.LT.NIMBK) THEN
!         IF (BKFLG(IJKBK-1).EQ..FALSE...AND.&
!             &I.EQ.BIE(IJKBK)) THEN
!           IK =(K-1)*NI+I
!           IJK=LK(K)+LI(I)+BYS(1)-1
!           FISBK(IBK,IK)=FI(IJK+1)+(FI(IJK+1)-FI(IJK+2))*FY(BYS(1))
!           !FINBK(IBK,IK)=... FOR NORTH BNDY FI
!         ELSE
!           IJK=LK(K)+LI(I)+1
!           FI(IJK)=FI(IJK+1)+(FI(IJK+1)-FI(IJK+2))*FY(2)
!         END IF
!       ELSE
!         IJK=LK(K)+LI(I)+1
!         FI(IJK)=FI(IJK+1)+(FI(IJK+1)-FI(IJK+2))*FY(2)
!       END IF
!       END DO
!     END DO
!     END IF
!     END DO
!     END DO
!     END DO
!
!     EXTEND TO OTHER RELEVANT SUBROUTINES, LOOPS, COORD DIRS ETC
!---------------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
! MPI
      INTEGER :: ITER
      DOUBLE PRECISION :: SOURCE
! TURB
      INTEGER :: IFI
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION (:) :: FI
! STRETCH
      INTEGER :: ICTR,JCTR,BYS1M1
!--------------------------------------------------------------
!
!.....MPI global comm
!
!
!.....MPI virtual topology
!
!.....I/O FILE NAMES. SUNIT-.res , TUNIT-.out .
! STUDIN
      FILIN='cgbpar.inp'
      OPEN (UNIT=5,FILE=FILIN)
      REWIND 5
! SYNTMON
      SUNIT=3
      TUNIT=2
! PLOT3D
      JUNIT=31
      KUNIT=51
!
!-----------------------------------------------------------
!.....READ INPUT DATA FROM UNIT 5
!-----------------------------------------------------------
! STUDIN
! TURB
      READ(5,6) TITLE
    6 FORMAT(A80)
      READ(5,*) LREAD,LWRITE,LLUD
      READ(5,*) MAXIT,IMON,JMON,KMON
      READ(5,*) IPR,JPR,KPR,SORMAX,SLARGE
      READ(5,*) DENSIT,VISC,PRM
      READ(5,*) TH,TC
      READ(5,*) UIN,VIN,WIN,PIN,TIN
      READ(5,*) TEIN,EDIN,ULID
      READ(5,*) (LCAL(I),I=1,NPHI)
      READ(5,*) (URF(I),I=1,NPHI)
      READ(5,*) (SOR(I),I=1,NPHI)
      READ(5,*) (NSW(I),I=1,NPHI)
      READ(5,*) (GDS(I),I=1,NPHI)
      READ(5,*) NI,NJ,NK
      READ(5,*) XMIN,YMIN,ZMIN
      READ(5,*) XMAX,YMAX,ZMAX
      READ(5,*) BYS1,BZE1
      READ(5,*) DFAC
! NANO
!     CPN is the specific heat at constant pressure for base fluid
!     (kJ/kgK). CPS is the specific heat at constant pressure for
!     solid nanoparticles (kJ/kgK). DENS is the desity of solid
!     nanoparticles (kg/m^3). PHIS is the volume fraction of solid
!     nanoparticles (nondimensional).
!
      READ(5,*) CPN,CPS,DENS,PHIS
! NANO
!     FK is the thermal conductivity of fluid (W/mK). SK is the
!     thermal conductivity of solid nanoparticles (W/mK).
!
      READ(5,*) FK,SKN
! STUDIN
      PRINT *,'TITLE:',TITLE
      PRINT *,'LREAD,LWRITE,LLUD:',LREAD,LWRITE,LLUD
      PRINT *,'MAXIT,IMON,JMON,KMON:',MAXIT,IMON,JMON,KMON
      PRINT *,'IPR,JPR,KPR,SORMAX,SLARGE:',IPR,JPR,KPR,SORMAX,SLARGE
      PRINT *,'DENSIT,VISC,PRM:',DENSIT,VISC,PRM
      PRINT *,'TH,TC:',TH,TC
      PRINT *,'UIN,VIN,WIN,PIN,TIN:',UIN,VIN,WIN
      PRINT *,'PIN,TIN:',PIN,TIN
      PRINT *,'TEIN,EDIN,ULID:',TEIN,EDIN,ULID
      PRINT *,'(LCAL(I),I=1,4):',(LCAL(I),I=1,4)
      PRINT *,'(LCAL(I),I=5,NPHI):',(LCAL(I),I=5,NPHI)
      PRINT *,'(URF(I),I=1,4):',(URF(I),I=1,4)
      PRINT *,'(URF(I),I=5,NPHI):',(URF(I),I=5,NPHI)
      PRINT *,'(SOR(I),I=1,4):',(SOR(I),I=1,4)
      PRINT *,'(SOR(I),I=5,NPHI):',(SOR(I),I=5,NPHI)
      PRINT *,'(NSW(I),I=1,4):',(NSW(I),I=1,4)
      PRINT *,'(NSW(I),I=5,NPHI):',(NSW(I),I=5,NPHI)
      PRINT *,'(GDS(I),I=1,4):',(GDS(I),I=1,4)
      PRINT *,'(GDS(I),I=5,NPHI):',(GDS(I),I=5,NPHI)
      PRINT *,'NI,NJ,NK:',NI,NJ,NK
      PRINT *,'XMIN,YMIN,ZMIN:',XMIN,YMIN,ZMIN
      PRINT *,'XMAX,YMAX,ZMAX:',XMAX,YMAX,ZMAX
      PRINT *,'BYS1,BZE1:',BYS1,BZE1
      PRINT *,'DFAC:',DFAC
      PRINT *
!
!.....INPUT AND BOUNDARY DATA, INITIALIZATION, OUTPUT TITLE,
!     ETC.
!
!     If LREAD is set true, results from previous run are read
!     before starting computation; if LWRITE is set true, the
!     results of calculation are written onto a file for post-
!     processing or continuation in a later run; if LLUD is set
!     true, 2nd order upwind scheme is activated when GDS(IU)
!     is greater than SMALL.
! STUDIN
!     LREAD=.FALSE.
!     LWRITE=.TRUE.
!     LLUD=.FALSE.
!
!.....SET SOME CONTROL VARIABLES
!
!     MAXIT is the maximum number of outer iterations to be performed
!     (e.g. one can run 10 outer iterations with LTEST=TRUE to check
!     if everything is OK, then resume calculation with LTEST=FALSE);
!     IMON and JMON are the I and J index of the monitoring location
!     (variable values at this location are printed after every outer
!     iteration); IPR and JPR are the I and J index of the node at
!     which the pressure is kept fixed (usually zero; reference
!     location); SORMAX is the level of the residual norm at which
!     outer iterations are stoped (converged solution); SLARGE is the
!     level of the residual norm at which iterations are stoped because
!     of divergence;
! STUDIN
!     MAXIT=3000
!     IMON=2
!     JMON=2
!     KMON=2
!     IPR=2
!     JPR=2
!     KPR=2
!     SORMAX=1.E-3
!     SLARGE=1.E3
!
!     DENSIT is the fluid density (here assumed constant); VISC is
!     the fluid dynamic viscosity (here assumed constant); PRANL is
!     the fluid Prandtl number; TH, TC, are the hot wall, cold wall
!     temperature, respectively (Degrees Celcius).
! STUDIN
!     DENSIT=100.
! TURB
!     VISC=1.E-3
! TURB - CHANGE TO PRANL
      PRANL=PRM
      PRR=1./PRANL
! T
! STUDIN
!     TH=100.
!     TC=0.
!
!     UIN, VIN, PIN, and TIN are the values of U, V, P, and T used
!     to initialize fields (usually zero, or some mean values);
!     ULID is the lid velocity for the lid-driven flow;
! STUDIN
!     UIN=0.
!     VIN=0.
!     WIN=0.
!     PIN=0.
!     TIN=0.
! TURB
!     TEIN=1.E-5
!     EDIN=1.E-7
! TURB
!     ULID=1.0
! T
      IU=1
      IV=2
      IW=3
      IP=4
      IEN=5
! TURB
      CALL MODDAT
      ZERO=0.
!
!.....FLUID,SOLID,NANOFLUID PARAMETERS.
! NANO
!     DNA1 IS RATIO OF DYNAMIC VISCOSITY OF NANOFLUID
!     COMPARED TO BASE FLUID.
!     DNA2 IS RATIO OF DENSITY OF NANOFLUID COMPARED TO
!     BASE FLUID.
!     DNA3 IS RATIO OF PRANDTL NUMBER OF NANOFLUID
!     COMPARED TO BASE FLUID.
!     RATK IS RATIO OF CONDUCTIVITY OF NANOFLUID COMPARED
!     TO BASE FLUID.
!
      DNA2=(1.-PHIS)*DENSIT+PHIS*DENS
      RATK=(SKN+2.*FK-2.*PHIS*(FK-SKN))/(SKN+2.*FK+PHIS*(FK-SKN))
      DNA3=(1.-PHIS)+PHIS*((DENS*CPS)/(DENSIT*CPN))
      DNA3=RATK/DNA3
!
!     LCAL(I) defines which equations are to be solved (I defines
!     variable as follows: 1 -> U, 2 -> V, 3 -> W, 4 -> P',
!     5 -> T).
! STUDIN
!     LCAL(IU)=.TRUE.
!     LCAL(IV)=.TRUE.
!     LCAL(IW)=.TRUE.
!     LCAL(IP)=.TRUE.
!     LCAL(IEN)=.FALSE.
!     LCAL(ITE)=.TRUE.
!     LCAL(IED)=.TRUE.
!     LCAL(IVIS)=.TRUE.
!
!     URF(I) is the under-relaxation factor for the Ith variable.
! STUDIN
!     URF(IU)=0.1
!     URF(IV)=0.1
!     URF(IW)=0.1
!     URF(IP)=0.025
!     URF(IEN)=0.1
!     URF(ITE)=0.1!0.7
!     URF(IED)=0.05!0.1!0.7
!     URF(IVIS)=1.
!
!     SOR(I) is the required ratio of reduction of the residual norm
!     during inner iterations for Ith variable before they are stopped
!     (e.g. value 0.2 means the residual norm should be reduced by a
!     factor of 5; this is usually sufficient for inner iterations
!     before updating the matrix).
! STUDIN
!     SOR(IU)=0.2
!     SOR(IV)=0.2
!     SOR(IW)=0.2
!     SOR(IP)=0.2
!     SOR(IEN)=0.2
!     SOR(ITE)=0.01
!     SOR(IED)=0.01
!     SOR(IVIS)=0.
!
!     NSW(I) is the maximum allowed number of inner iterations for the
!     Ith variable (for U, V, and T, one inner iteration by SIP is
!     sufficiant; for P', 5 to 10 may be required to satisfy the
!     convergencee criterion).
! STUDIN
!     NSW(IU)=5
!     NSW(IV)=5
!     NSW(IW)=5
!     NSW(IP)=20
!     NSW(IEN)=5
!     NSW(ITE)=100
!     NSW(IED)=100!10
!     NSW(IVIS)=0
!
!     GDS(I) is the blending factor for UDS and CDS in the equation for
!     Ith variable (convective terms; value 1.0 means CDS (second
!     order), 0.0 means UDS (first order), any value between 0.0 and
!     1.0 can be used). The value 1.0 is recomended, except for coarse
!     grids, in case convergence problems are encountered.
! STUDIN
!     GDS(IU)=0.
!     GDS(IV)=0.
!     GDS(IW)=0.
!     GDS(IP)=0.
!     GDS(IEN)=0.
!     GDS(ITE)=0.
!     GDS(IED)=0.
!     GDS(IVIS)=0.
! OUTLET
      SMALL=1.E-30
!
!.....READ GRID DATA
! BLK
! STUDIN
!     NI=22
!     NJ=22
!     NK=602
      NIJ=NI*NJ
      NIK=NI*NK
!
      ALLOCATE (LI(1:NI))
      ALLOCATE (LK(1:NK))
!
      DO I=1,NI
        LI(I)=(I-1)*NJ
      END DO
      DO K=1,NK
        LK(K)=(K-1)*NIJ
      END DO
!
      NIM=NI-1
      NJM=NJ-1
      NKM=NK-1
      NIJK=NI*NJ*NK
!
      ALLOCATE (X(1:NI),XC(1:NI))
      ALLOCATE (Y(1:NJ),YC(1:NJ))
      ALLOCATE (Z(1:NK),ZC(1:NK))
!
!.....Memory allocation
!
      X=0.
! STUDIN
!     XMIN=0.
!     XMAX=0.02
! 3D
      Y=0.
! STUDIN
!     YMIN=0.
!     YMAX=0.02
! 3D
! BLK
      Z=0.
! STUDIN
!     ZMIN=0.
!     ZMAX=0.6
! STUDIN
! STRETCH
!     DFAC=1.2
! PLOT3D
      MACH=REAL(ULID/343.2)
      ALPHA=REAL(0.)
      REYNOLDS=REAL(DENSIT*ULID*(YMAX-YMIN)/VISC)
      TIME=REAL(0.)
!
!.....GRID GENERATION - XDIR. SEE CFD WIKI FOR MORE DETAILS.
!     SEE http://www.cfd-online.com/Wiki/Structured_mesh_generation
! STRETCH
      X(1)=0.
! STRETCH
      DX=1./DBLE((NI-2)/2)
      DO I=2,(NI-2)/2+1,1
        X(I)=X(I-1)+DX
      END DO
! STRETCH
      DO I=2,(NI-2)/2+1,1
        X(I)=1.+(TANH(DFAC*(X(I)-1.)/2.))/(TANH(DFAC/2.))
      END DO
! STRETCH
      DO I=2,(NI-2)/2+1,1
        X(I)=X(I)*((XMAX-XMIN)/2.)+XMIN
      END DO
! STRETCH
      X(NIM)=XMAX
! STRETCH
      ICTR=2
      X(1)=XMIN
      DO I=NIM-1,NI/2+1,-1
        X(I)=X(I+1)-(X(ICTR)-X(ICTR-1))
        ICTR=ICTR+1
      END DO
! STRETCH
      X(1)=XMIN
      X(NI)=X(NIM)
!
!.....GRID GENERATION - YDIR BLK 2.
!     SEE http://www.cfd-online.com/Wiki/Structured_mesh_generation
! STRETCH
      BYS1M1=BYS1-1
      Y(1)=0.
! STRETCH
      DY=1./DBLE((BYS1M1-1)/2)
      DO J=2,(BYS1M1-1)/2+1,1
        Y(J)=Y(J-1)+DY
      END DO
! STRETCH
      DO J=2,(BYS1M1-1)/2+1,1
        Y(J)=1.+(TANH(DFAC*(Y(J)-1.)/2.))/(TANH(DFAC/2.))
      END DO
! STRETCH
      Y(BYS1M1)=YMIN+(YMAX-YMIN)*((DBLE(BYS1M1)-1.)/(DBLE(NJ)-2.))
      DO J=2,(BYS1M1-1)/2+1,1
        Y(J)=Y(J)*((Y(BYS1M1)-YMIN)/2.)+YMIN
      END DO
! STRETCH
      JCTR=2
      Y(1)=YMIN
      DO J=BYS1M1-1,(BYS1M1+1)/2+1,-1
        Y(J)=Y(J+1)-(Y(JCTR)-Y(JCTR-1))
        JCTR=JCTR+1
      END DO
! STRETCH
      Y(1)=YMIN
!
!.....GRID GENERATION - YDIR BLKS 1 AND 3.
!     SEE http://www.cfd-online.com/Wiki/Structured_mesh_generation
! STRETCH
      Y(BYS1M1)=0.
! STRETCH
      DY=1./DBLE((NJ-BYS1M1-1)/2)
      DO J=BYS1M1+1,BYS1M1+((NJ-BYS1M1-1)/2),1
        Y(J)=Y(J-1)+DY
      END DO
! STRETCH
      DO J=BYS1M1+1,BYS1M1+((NJ-BYS1M1-1)/2),1
        Y(J)=1.+(TANH(DFAC*(Y(J)-1.)/2.))/(TANH(DFAC/2.))
      END DO
! STRETCH
      Y(BYS1M1)=YMIN+(YMAX-YMIN)*((DBLE(BYS1M1)-1.)/(DBLE(NJ)-2.))
      DO J=BYS1M1+1,BYS1M1+((NJ-BYS1M1-1)/2),1
        Y(J)=Y(J)*((YMAX-Y(BYS1M1))/2.)+Y(BYS1M1)
      END DO
! STRETCH
      Y(NJM)=YMAX
! STRETCH
      JCTR=BYS1
      Y(BYS1M1)=YMIN+(YMAX-YMIN)*((DBLE(BYS1M1)-1.)/(DBLE(NJ)-2.))
      DO J=NJM-1,BYS1M1+((NJ-BYS1M1-1)/2)+1,-1
        Y(J)=Y(J+1)-(Y(JCTR)-Y(JCTR-1))
        JCTR=JCTR+1
      END DO
! STRETCH
      Y(NJ)=Y(NJM)
!
!     Y(1)=YMIN
!     DY=(YMAX-YMIN)/DBLE(NJ-2)
!     DO J=2,NJM,1
!       Y(J)=Y(J-1)+DY
!     END DO
!     Y(NJ)=Y(NJM)
!
!     PRINT *
!     DO J=1,NJ,1
!       PRINT *,'J,Y(J):',J,Y(J)
!     END DO
!     PRINT *
!
!.....GRID GENERATION - ZDIR BLKS 2 AND 3.
!     SEE http://www.cfd-online.com/Wiki/Structured_mesh_generation
! STRETCH
      Z(BZE1)=0.
! STRETCH
      DZ=1./DBLE(NK-BZE1-1)
      DO K=BZE1+1,BZE1+(NK-BZE1-1),1
        Z(K)=Z(K-1)+DZ
      END DO
! STRETCH
      DO K=BZE1+1,BZE1+(NK-BZE1-1),1
        Z(K)=1.+(TANH(DFAC*(Z(K)-1.)/2.))/(TANH(DFAC/2.))
      END DO
!
      Z(BZE1)=ZMIN+(ZMAX-ZMIN)*((DBLE(BZE1)-1.)/(DBLE(NK)-2.))
      DO K=BZE1+1,BZE1+(NK-BZE1-1),1
        Z(K)=Z(K)*(ZMAX-Z(BZE1))+Z(BZE1)
      END DO
!
      DZ=(ZMAX-ZMIN)/DBLE(NK-2)
      DO K=BZE1+1,BZE1+(NK-BZE1-1),1
        Z(K)=Z(K-1)+DZ
      END DO
!
      Z(NKM)=ZMAX
      Z(NK) =ZMAX
!
!.....GRID GENERATION - ZDIR BLK 1.
!     SEE http://www.cfd-online.com/Wiki/Structured_mesh_generation
! STRETCH
      Z(1)=1.
! STRETCH
      DZ=1./DBLE(BZE1-1)
      DO K=2,BZE1,1
        Z(K)=Z(K-1)-DZ
      END DO
! STRETCH
      DO K=2,BZE1,1
        Z(K)=ABS((TANH(DFAC*(Z(K)-1.)/2.))/(TANH(DFAC/2.)))
      END DO
! STRETCH
      Z(BZE1)=ZMIN+(ZMAX-ZMIN)*((DBLE(BZE1)-1.)/(DBLE(NK)-2.))
      DO K=2,BZE1-1,1
        Z(K)=Z(K)*(Z(BZE1)-ZMIN)+ZMIN
      END DO
! STRETCH
      Z(1)=ZMIN
!
!.....X- COORDINATES OF CV-CENTERS
!
      XC=0.
      DO I=2,NIM
        XC(I)=0.5*(X(I)+X(I-1))
      END DO
      XC(1)=X(1)
      XC(NI)=X(NIM)
!
!.....Y- COORDINATES OF CV-CENTERS
!
      YC=0.
      DO J=2,NJM
        YC(J)=0.5*(Y(J)+Y(J-1))
      END DO
      YC(1)=Y(1)
      YC(NJ)=Y(NJM)
!
!.....Z- COORDINATES OF CV-CENTERS
! 3D
      ZC=0.
      DO K=2,NKM
        ZC(K)=0.5*(Z(K)+Z(K-1))
      END DO
      ZC(1)=Z(1)
      ZC(NK)=Z(NKM)
!
!.....INTERPOLATION FACTORS (see Sect. 4.4.2, Eq. (4.14))
! 3D
      ALLOCATE (FX(1:NIM))
      ALLOCATE (FY(1:NJM))
      ALLOCATE (FZ(1:NKM))
!
      FX=0.
      DO I=1,NIM
        FX(I)=(X(I)-XC(I))/(XC(I+1)-XC(I))
      END DO
!
      FY=0.
      DO J=1,NJM
        FY(J)=(Y(J)-YC(J))/(YC(J+1)-YC(J))
      END DO
! 3D
      FZ=0.
      DO K=1,NKM
        FZ(K)=(Z(K)-ZC(K))/(ZC(K+1)-ZC(K))
      END DO
!
!---------------------------------------------------
!.....BOUNDARY AND INITIAL CONDITIONS
!---------------------------------------------------
!
!.....NORTH WALL VELOCITY (FOR LID-DRIVEN CAVITY)
! 3D
! T
      ALLOCATE (U(1:NIJK))
      ALLOCATE (V(1:NIJK))
      ALLOCATE (W(1:NIJK))
      ALLOCATE (UO(1:NIJK))
      ALLOCATE (VO(1:NIJK))
      ALLOCATE (WO(1:NIJK))
      ALLOCATE (SU(1:NIJK))
      ALLOCATE (SV(1:NIJK))
      ALLOCATE (SW(1:NIJK))
      ALLOCATE (P(1:NIJK))
      ALLOCATE (PP(1:NIJK))
      ALLOCATE (PO(1:NIJK))
      ALLOCATE (F1(1:NIJK))
      ALLOCATE (F2(1:NIJK))
      ALLOCATE (F3(1:NIJK))
      ALLOCATE (DPX(1:NIJK))
      ALLOCATE (DPY(1:NIJK))
      ALLOCATE (DPZ(1:NIJK))
      ALLOCATE (AE(1:NIJK))
      ALLOCATE (AW(1:NIJK))
      ALLOCATE (AN(1:NIJK))
      ALLOCATE (AS(1:NIJK))
      ALLOCATE (AT(1:NIJK))
      ALLOCATE (AB(1:NIJK))
      ALLOCATE (AP(1:NIJK))
      ALLOCATE (APU(1:NIJK))
      ALLOCATE (APV(1:NIJK))
      ALLOCATE (APW(1:NIJK))
      ALLOCATE (T(1:NIJK))
      ALLOCATE (TO(1:NIJK))
! TURB
      ALLOCATE (ED(1:NIJK))
      ALLOCATE (GEN(1:NIJK))
      ALLOCATE (TE(1:NIJK))
      ALLOCATE (VIS(1:NIJK))
      ALLOCATE (VISW(1:NIJK))
      ALLOCATE (VOL(1:NIJK))
      ALLOCATE (YPL(1:NIJK))
!
!.....LOCAL VARS
! TURB
      ALLOCATE (FI(1:NIJK))
!
!.....ZERO OUT ALL THE ABOVE
!
      U=0.;V=0.;W=0.
      UO=0.;VO=0.;WO=0.
      SU=0.;SV=0.;SW=0.
      P=0.;PP=0.;PO=0.
      F1=0.;F2=0.;F3=0.
      DPX=0.;DPY=0.;DPZ=0.
      AE=0.;AW=0.;AN=0.;AS=0.;AT=0.;AB=0.;AP=0.
      APU=0.;APV=0.;APW=0.
      T=0.;TO=0.
! TURB
      ED=0.
      GEN=0.
      TE=0.
      VIS=0.
      VISW=0.
      VOL=0.
      YPL=0.
!
!.....INITIALIZE LOCAL VARS
! TURB
      FI=0.
!
!.....BLK
! BLK
      ALLOCATE (BYS(0:3))
      ALLOCATE (BYE(0:3))
      ALLOCATE (BZS(0:3))
      ALLOCATE (BZE(0:3))
! BLK
      ALLOCATE (US1(1:NIK))
      ALLOCATE (VS1(1:NIK))
      ALLOCATE (WS1(1:NIK))
      ALLOCATE (PS1(1:NIK))
      ALLOCATE (PPS1(1:NIK))
      ALLOCATE (TS1(1:NIK))
! BLK
      ALLOCATE (UB2(1:NIJ))
      ALLOCATE (VB2(1:NIJ))
      ALLOCATE (WB2(1:NIJ))
      ALLOCATE (PB2(1:NIJ))
      ALLOCATE (PPB2(1:NIJ))
      ALLOCATE (TB2(1:NIJ))
! BLK
      BYS=0;BYE=0;BZS=0;BZE=0
! BLK
      US1=0.;VS1=0.;WS1=0.;PS1=0.;PPS1=0.;TS1=0.
      UB2=0.;VB2=0.;WB2=0.;PB2=0.;PPB2=0.;TB2=0.
!
!.....INLET OUTLET
! OUTLET
      ALLOCATE (FMO(1:NIJ))
      ALLOCATE (FMI(1:NIJ))
! OUTLET
      FMO=0.
      FMI=0.
!
!.....BLK
! BLK
! STUDIN
!     BYS(1)=6;       BYE(1)=NJM
      BYS(1)=BYS1;    BYE(1)=NJM
      BYS(2)=2;       BYE(2)=BYS(1)-1
      BYS(3)=BYE(2)+1;BYE(3)=NJM
!     BZS(1)=2;       BZE(1)=201
      BZS(1)=2;       BZE(1)=BZE1
      BZS(2)=BZE(1)+1;BZE(2)=NKM
      BZS(3)=BZS(2);  BZE(3)=NKM
!
!.....INTERNAL CELL VOLUME
! TURB
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DZ=Z(K)-Z(K-1)
        DO I=2,NIM
          DX=X(I)-X(I-1)
          DO J=BYS(B),BYE(B)
            DY=Y(J)-Y(J-1)
            IJK=LK(K)+LI(I)+J
            VOL(IJK)=DX*DY*DZ
          END DO
        END DO
      END DO
      END DO
!
!---------------------------------------------------
!.....TEMPERATURE, TE, ED, VIS BOUNDARIES BLK1
!---------------------------------------------------
!
!.....SOUTH BOUNDARY (ISOTHERMAL WALL)
! BLK1
      DO K=2,BZE(1)
        DO I=2,NIM
          IK =(K-1)*NI+I
          IJK=LK(K)+LI(I)+BYS(1)
          TS1(IK)=TH!T(IJK-1)=TH
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....NORTH BOUNDARY (ISOTHERMAL WALL)
! BLK1
! NANO
      DO K=2,BZE(1)!NKM
        DO I=2,NIM
          IJK=LK(K)+LI(I)+NJM
          T(IJK+1)=TH
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....WEST BOUNDARY (ADIABATIC WALL)
!
      DO K=BZS(1),BZE(1)
        DO J=BYS(1),BYE(1)
          IJK=LK(K)+LI(1)+J
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....EAST BOUNDARY (ADIABATIC WALL)
!
      DO K=BZS(1),BZE(1)
        DO J=BYS(1),BYE(1)
          IJK=LK(K)+LI(NI)+J
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....BOTTOM BOUNDARY (ISOTHERMAL INLET)
! BLK1
! INLET
! TURB
! NANO NO CHANGES
      DO I=2,NIM
      DO J=BYS(1),NJM
        IJK=LK(2)+LI(I)+J
        T(IJK-NIJ)=TH
        TE(IJK-NIJ)=TEIN
        ED(IJK-NIJ)=EDIN
        VIS(IJK-NIJ)=VISC
        VISW(IJK-NIJ)=VIS(IJK-NIJ)
      END DO
      END DO
!
!.....TOP IS NOT A BOUNDARY
!
!---------------------------------------------------
!.....TEMPERATURE, TE, ED, VIS  BOUNDARIES BLK2
!---------------------------------------------------
! BLK
!.....SOUTH BOUNDARY (ISOTHERMAL WALL)
!
      DO K=BZS(2),NKM
        DO I=2,NIM
          IJK=LK(K)+LI(I)+2
          T(IJK-1)=TC
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....NORTH IS NOT A BOUNDARY
!
!.....WEST BOUNDARY (ADIABATIC WALL)
!
      DO K=BZS(2),BZE(2)
        DO J=BYS(2),BYE(2)
          IJK=LK(K)+LI(1)+J
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....EAST BOUNDARY (ADIABATIC WALL)
!
      DO K=BZS(2),BZE(2)
        DO J=BYS(2),BYE(2)
          IJK=LK(K)+LI(NI)+J
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....BOTTOM BOUNDARY (ISOTHERMAL WALL)
! WALL
! BLK
      DO I=2,NIM
        DO J=2,BYE(2)
          IJ =(I-1)*NJ+J
          IJK=LK(BZS(2))+LI(I)+J
          TB2(IJ)=TC!T(IJK-NIJ)=TC
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....TOP BOUNDARY (OUTLET, SEE SUBROUTINES BCT, CALCSC, MODVIS)
!
      DO I=2,NIM
      DO J=BYS(2),BYE(2)
        IJK=LK(NKM)+LI(I)+J
        VIS(IJK+NIJ)=VISC
        VISW(IJK+NIJ)=VIS(IJK+NIJ)
      END DO
      END DO
!
!---------------------------------------------------
!.....TEMPERATURE, TE, ED, VIS  BOUNDARIES BLK3
!---------------------------------------------------
! BLK
!.....SOUTH IS NOT A BOUNDARY
!
!.....NORTH BOUNDARY (ISOTHERMAL WALL)
! BLK
! NANO NO CHANGES
      DO K=BZS(3),NKM
        DO I=2,NIM
          IJK=LK(K)+LI(I)+NJM
          T(IJK+1)=TH
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....WEST BOUNDARY (ADIABATIC WALLS, SEE SUBROUTINE BCT)
!
      DO K=BZS(3),BZE(3)
        DO J=BYS(3),BYE(3)
          IJK=LK(K)+LI(1)+J
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....EAST BOUNDARY (ADIABATIC WALLS, SEE SUBROUTINE BCT)
!
      DO K=BZS(3),BZE(3)
        DO J=BYS(3),BYE(3)
          IJK=LK(K)+LI(NI)+J
          VIS(IJK)=VISC
          VISW(IJK)=VIS(IJK)
        END DO
      END DO
!
!.....BOTTOM IS NOT A BOUNDARY
!
!.....TOP BOUNDARY (OUTLET, SEE SUBROUTINES BCT, CALCSC, MODVIS)
!
      DO I=2,NIM
      DO J=BYS(3),BYE(3)
        IJK=LK(NKM)+LI(I)+J
        VIS(IJK+NIJ)=VISC
        VISW(IJK+NIJ)=VIS(IJK+NIJ)
      END DO
      END DO
!
!.....INLET BOUNDARY CONDITIONS
! INLET
      CALL BCIN
!
!.....INITIAL VARIABLE VALUES (INITIAL CONDITIONS)
! 3D
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            U(IJK)=UIN
            V(IJK)=VIN
            W(IJK)=WIN
            T(IJK)=TIN
            P(IJK)=PIN
! TURB
            TE(IJK)=TEIN
            ED(IJK)=EDIN
!
            UO(IJK)=UIN
            VO(IJK)=VIN
            WO(IJK)=WIN
            TO(IJK)=TIN
! INLET
! TURB
            VIS(IJK)=VISC
!
          END DO
        END DO
      END DO
!
      END DO
!
!------------------------------------------------------
!.....INITIAL OUTPUT - PRINTOUT OF FLOW PARAMETERS
!------------------------------------------------------
!
      OPEN (UNIT=TUNIT,FILE='cgbpar.out',ACCESS='APPEND'&
      &,STATUS='NEW')
      WRITE(TUNIT,601) 'Flow over a backstep',DENSIT,VISC
  601 FORMAT(1H1,//,10X,A80,/,10X,50('*'),/,10X,&
     &        ' FLUID DENSITY    :  ',1P1E10.4,/,10X,&
     &        ' DYNAMIC VISCOSITY:  ',1P1E10.4)
      IF(ULID.NE.0.) THEN
        WRITE(TUNIT,*) '          MAX. INL VELOCITY:  ',ULID
      ENDIF
      WRITE(TUNIT,*) '  '
      WRITE(TUNIT,*) '          UNDERRELAXATION  FACTORS'
      WRITE(TUNIT,*) '          ========================'
      WRITE(TUNIT,*) '          U-VELOCITY  :  ',URF(IU)
      WRITE(TUNIT,*) '          V-VELOCITY  :  ',URF(IV)
      WRITE(TUNIT,*) '          W-VELOCITY  :  ',URF(IW)
      WRITE(TUNIT,*) '          PRESSURE    :  ',URF(IP)
      WRITE(TUNIT,*) '  '
      WRITE(TUNIT,*) '          SPATIAL BLENDING FACTORS &
      &(CDS-UDS)'
      WRITE(TUNIT,*) '&
      &          =================================='
      WRITE(TUNIT,*) '          U-VELOCITY  :  ',GDS(IU)
      WRITE(TUNIT,*) '          V-VELOCITY  :  ',GDS(IV)
      WRITE(TUNIT,*) '          W-VELOCITY  :  ',GDS(IW)
      WRITE(TUNIT,*) '  '
      WRITE(TUNIT,*) '  '
      WRITE(TUNIT,*) '  '
      WRITE(TUNIT,*) '     *****************************'
      WRITE(TUNIT,*) '  '
      WRITE(TUNIT,600) IMON,JMON,KMON
!
!.....READ RESULTS OF PREVIOUS TIME STEP (IF CONTINUATION)
! BLK
! T
! TURB
      IF (LREAD) THEN
        OPEN (UNIT=SUNIT,FILE='cgbpar.res',STATUS='UNKNOWN',&
        &FORM='UNFORMATTED')
          READ(SUNIT) (U(IJK),IJK=1,NIJK),&
                     &(V(IJK),IJK=1,NIJK),&
                     &(W(IJK),IJK=1,NIJK),&
                     &(P(IJK),IJK=1,NIJK),&
                     &(T(IJK),IJK=1,NIJK),&
                     &(TE(IJK),IJK=1,NIJK),&
                     &(ED(IJK),IJK=1,NIJK),&
                     &(VIS(IJK),IJK=1,NIJK),&
                     &(F1(IJK),IJK=1,NIJK),&
                     &(F2(IJK),IJK=1,NIJK),&
                     &(F3(IJK),IJK=1,NIJK)
        CLOSE(SUNIT)
      ENDIF
!
!==============================================
!.....TIME LOOP
!==============================================
!
!.....DEFINE MONITORING LOCATION (NODE WITH I=IMON, J=JMON,
!     K=KMON)
!
      IJKMON=LK(KMON)+LI(IMON)+JMON
!
!.....SET STARTING ITERTION NUMBER FOR LUD
! LUD
      LUST=MAXIT-3
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!.....OUTER ITERATIONS (SIMPLE RELAXATIONS)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      DO ITER=1,MAXIT
! LUD
        IF (LLUD.AND.GDS(IU).LT.SMALL) THEN
          IF (ITER.GT.LUST) THEN
            LFAC=1.
          ELSE
            LFAC=0.
          END IF
        END IF
!
        IF(LCAL(IU)) CALL CALCUVW
!
        IF(LCAL(IP)) CALL CALCP
! T
! TURB
        IF (LCAL(IEN)) THEN
          FI=T
          IFI=IEN
          CALL CALCSC(FI,IFI)
          T=FI
        END IF
! TURB
        IF (LCAL(ITE)) THEN
          FI=TE
          IFI=ITE
          CALL CALCSC(FI,IFI)
          TE=FI
        END IF
! TURB
        IF (LCAL(IED)) THEN
          FI=ED
          IFI=IED
          CALL CALCSC(FI,IFI)
          ED=FI
        END IF
! TURB
        IF (LCAL(IVIS)) THEN
          CALL MODVIS
        END IF
!
        PRINT *,ITER,' / ',MAXIT
!
!.....CHECK CONVERGENCE OF OUTER ITERATIONS
! T
! TURB
        WRITE(TUNIT,606) ITER,RESOR(IU),RESOR(IW),RESOR(IP),&
        &RESOR(IEN),RESOR(ITE),RESOR(IED),U(IJKMON),W(IJKMON),&
        &P(IJKMON),T(IJKMON),TE(IJKMON),ED(IJKMON)
! TURB
!       SOURCE=MAX(RESOR(IU),RESOR(IV),RESOR(IW),RESOR(IP),&
!       &RESOR(IEN),RESOR(ITE))
        SOURCE=RESOR(IP)
!
        IF(SOURCE.GT.SLARGE) GO TO 510
        IF(SOURCE.LT.SORMAX) GO TO 250
!
      END DO
!
!     CLOSE(TUNIT)
!
  250 CONTINUE
!
!.....CONVERGED: IF UNSTEADY FLOW, PRINT AND SAVE NPRTth
!     SOLUTION
!
      CALL P3D
!
      WRITE(TUNIT,*)
      WRITE(TUNIT,*) 'IJK   U   V   W   P   TE   ED'
      DO J=NJ,1,-1
      IJK=LK(202)+LI(11)+J
      WRITE(TUNIT,*) IJK,U(IJK),V(IJK),W(IJK),P(IJK),TE(IJK),ED(IJK)
      END DO
!
      CLOSE(TUNIT)
!
!==============================================================
!.....ALL TIME STEPS DONE; SAVE LAST SOLUTION FOR CONTINUATION
!==============================================================
! BLK
! T
! TURB
      IF (LWRITE) THEN
        OPEN (UNIT=SUNIT,FILE='cgbpar.res',STATUS='UNKNOWN',&
        &FORM='UNFORMATTED')
          WRITE(SUNIT) (U(IJK),IJK=1,NIJK),&
                      &(V(IJK),IJK=1,NIJK),&
                      &(W(IJK),IJK=1,NIJK),&
                      &(P(IJK),IJK=1,NIJK),&
                      &(T(IJK),IJK=1,NIJK),&
                      &(TE(IJK),IJK=1,NIJK),&
                      &(ED(IJK),IJK=1,NIJK),&
                      &(VIS(IJK),IJK=1,NIJK),&
                      &(F1(IJK),IJK=1,NIJK),&
                      &(F2(IJK),IJK=1,NIJK),&
                      &(F3(IJK),IJK=1,NIJK)
        CLOSE(SUNIT)
      ENDIF
!
      STOP
!
!==============================================================
!......MESSAGE FOR DIVERGENCE
!==============================================================
!
  510 PRINT *,'  *** TERMINATED - OUTER ITERATIONS&
      & DIVERGING ***'
!
      WRITE(TUNIT,*)
      WRITE(TUNIT,*) 'IJK   U   V   W   P   TE   ED'
      DO J=NJ,1,-1
      IJK=LK(202)+LI(11)+J
      WRITE(TUNIT,*) IJK,U(IJK),V(IJK),W(IJK),P(IJK),TE(IJK),ED(IJK)
      END DO
!
      CLOSE(TUNIT)
!
      STOP
!
!==============================================================
!......FORMAT SPECIFICATIONS
!==============================================================
! TURB
  600 FORMAT(1X,'ITER.',3X,'I--------------------ABSOLUTE&
      & RESIDUAL SOURCE SUMS-------------------I',3X,&
      &'I---------------FIELD VALUES AT MONITORING LOCATION&
      & (',I3,',',I3,',',I3,&
      &')-------------I',/,2X,'NO.',9X,'U',11X,'W',11X,'MASS',11X,'T',&
      &11X,'TE',9X,'ED',16X,'U',11X,'W',11X,'P',11X,'T',11X,'TE',&
      &11X,'ED',/)
! T
! TURB
  606 FORMAT(1X,I4,2X,1P6E12.4,5X,1P6E12.4)
!
      END PROGRAM PCOL
!
!###########################################################
      SUBROUTINE CALCUVW
!###########################################################
!-----------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
!
      DOUBLE PRECISION :: URFU,URFV,URFW
      DOUBLE PRECISION :: PE,PW,PN,PS,PT,PB
      DOUBLE PRECISION :: FUUDS,FVUDS,FWUDS,FUCDS,FVCDS,FWCDS
! LUD
      DOUBLE PRECISION :: ALIFUE=0.,ALIFUN=0.,ALIFUT=0.
      DOUBLE PRECISION :: ALIFVE=0.,ALIFVN=0.,ALIFVT=0.
      DOUBLE PRECISION :: ALIFWE=0.,ALIFWN=0.,ALIFWT=0.
! LUD
      DOUBLE PRECISION :: FULUD=0.,FVLUD=0.,FWLUD=0.
      DOUBLE PRECISION :: SAIMUE=0.,SAIMUN=0.,SAIMUT=0.
      DOUBLE PRECISION :: SAIMVE=0.,SAIMVN=0.,SAIMVT=0.
      DOUBLE PRECISION :: SAIMWE=0.,SAIMWN=0.,SAIMWT=0.
      DOUBLE PRECISION :: SAIPUE=0.,SAIPUN=0.,SAIPUT=0.
      DOUBLE PRECISION :: SAIPVE=0.,SAIPVN=0.,SAIPVT=0.
      DOUBLE PRECISION :: SAIPWE=0.,SAIPWN=0.,SAIPWT=0.
!-----------------------------------------------------------
! T
!
!.....RECIPROCAL VALUES OF UNDER-RELAXATION FACTORS FOR U AND V
! 3D
      URFU=1./URF(IU)
      URFV=1./URF(IV)
      URFW=1./URF(IW)
!
!.....SET BOUNDARY PRESSURE (LINEAR EXTRAPOLATION FROM INSIDE)
!
! MPI
      CALL PBOUND(P,PS1,PB2)
!
!.....INITIALIZE TEMPORARILY STORED VARIABLES
! 3D
      DO IJK=1,NIJK
        SU(IJK)=0.
        SV(IJK)=0.
        SW(IJK)=0.
        APU(IJK)=0.
        APV(IJK)=0.
        APW(IJK)=0.
      END DO
!
!==========================================================
!.....FLUXES THROUGH INTERNAL EAST CV FACES
!=========================================================
!     F1(IJK) is the mass flux through the east face (outward
!     normal directed to E); FX(I) is the ratio of distance
!     from P to cell face, to distance from P to E; IJK
!     denotes node P and IJKE node E.
!     Contribution of convective and diffusive fluxes from
!     east face to AE(P), AW(E), and source terms at both
!     P and E are calculated; contributions to AP(P) and
!     AP(E) come through the sum of neighbor coefficients
!     and are not explicitly calculated.
!----------------------------------------------------------
! 3D
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
       DO I=2,NIM-1
!
!.....INTERPOLATION FACTORS, DISTANCE FROM P TO E (SAME FOR ALL J)
!
        FXE =FX(I)
        FXP =1.-FXE
        DXPE=XC(I+1)-XC(I)
! 3D
        DO J=BYS(B),BYE(B)
!
          IJK =LK(K)+LI(I)+J
          IJKE=IJK+NJ
!
!.....CELL FACE AREA S = DY*DZ
! 3D
          S=(Y(J)-Y(J-1))*(Z(K)-Z(K-1))
!
!.....CALCULATE VISCOSITY
! TURB
          VISI=VIS(IJKE)*FX(I)+VIS(IJK)*(1.-FX(I))
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX (SAME FOR U AND V)
!
!         D=VISC*S/DXPE
! TURB
          D=VISI*S/DXPE
! NANO
          DNA1=VISI/((1.-PHIS)**2.5)
          D=DNA1*S/DXPE
!
!.....EXPLICIT CONVECTIVE FLUXES FOR UDS AND CDS
! 3D
!
          IF (F1(IJK).LE.0.) THEN
            CE=F1(IJK)
          ELSE
            CE=0.
          END IF
          IF (F1(IJK).GT.0.) THEN
            CP=F1(IJK)
          ELSE
            CP=0.
          END IF
!
!.....LUD
!
          FULUD=0.;FVLUD=0.;FWLUD=0.
          IF (I.GT.2.AND.I.LT.NIM-1) THEN
          IF (J.GT.2.AND.J.LT.NJM-1) THEN
          IF (K.GT.2.AND.K.LT.NKM-1) THEN
!
!.....EQS 5.83, VERSTEEG & MALALASEKERA. ALIFUE OF CELL IJ IS ALIFUW OF
!     CELL IJE, SO, NO NEED TO REPEAT.
! LUD
            IF (F1(IJK).LT.0.) THEN
              ALIFUE=0.
              ALIFVE=0.
              ALIFWE=0.
              SAIMUE=(U(IJK+NJ+NJ)-U(IJK+NJ))/((U(IJK+NJ)-U(IJK))+SMALL)
              SAIMVE=(V(IJK+NJ+NJ)-V(IJK+NJ))/((V(IJK+NJ)-V(IJK))+SMALL)
              SAIMWE=(W(IJK+NJ+NJ)-W(IJK+NJ))/((W(IJK+NJ)-W(IJK))+SMALL)
            END IF
! LUD
            IF (F1(IJK).GE.0.) THEN
              ALIFUE=1.
              ALIFVE=1.
              ALIFWE=1.
              SAIPUE=(U(IJK)-U(IJK-NJ))/((U(IJK+NJ)-U(IJK))+SMALL)
              SAIPVE=(V(IJK)-V(IJK-NJ))/((V(IJK+NJ)-V(IJK))+SMALL)
              SAIPWE=(W(IJK)-W(IJK-NJ))/((W(IJK+NJ)-W(IJK))+SMALL)
            END IF
! LUD2
            T1=(1.-ALIFUE)*SAIMUE-ALIFUE*SAIPUE
            T2=U(IJK+NJ)-U(IJK)
            FULUD=0.5*F1(IJK)*T1*T2
! LUD2
            T1=(1.-ALIFVE)*SAIMVE-ALIFVE*SAIPVE
            T2=V(IJK+NJ)-V(IJK)
            FVLUD=0.5*F1(IJK)*T1*T2
! LUD2
            T1=(1.-ALIFWE)*SAIMWE-ALIFWE*SAIPWE
            T2=W(IJK+NJ)-W(IJK)
            FWLUD=0.5*F1(IJK)*T1*T2
!
          END IF
          END IF
          END IF
!
!.....COEFFICIENTS AE(P) AND AW(E) DUE TO UDS
! 3D
          AE(IJK) = CE-D
          AW(IJKE)=-CP-D
! LUD
! 3D
          FUUDS=CP*U(IJK)+CE*U(IJKE)
          FVUDS=CP*V(IJK)+CE*V(IJKE)
          FWUDS=CP*W(IJK)+CE*W(IJKE)
          FUCDS=F1(IJK)*(U(IJKE)*FXE+U(IJK)*FXP)
          FVCDS=F1(IJK)*(V(IJKE)*FXE+V(IJK)*FXP)
          FWCDS=F1(IJK)*(W(IJKE)*FXE+W(IJK)*FXP)
!
!.....SOURCE TERM CONTRIBUTIONS AT P AND E DUE TO DEFERRED
!     CORRECTION
! 3D
! LUD
          SU(IJK) =SU(IJK) +GDS(IU)*(FUUDS-FUCDS)+LFAC*FULUD
          SU(IJKE)=SU(IJKE)-GDS(IU)*(FUUDS-FUCDS)-LFAC*FULUD
          SV(IJK) =SV(IJK) +GDS(IU)*(FVUDS-FVCDS)+LFAC*FVLUD
          SV(IJKE)=SV(IJKE)-GDS(IU)*(FVUDS-FVCDS)-LFAC*FVLUD
          SW(IJK) =SW(IJK) +GDS(IU)*(FWUDS-FWCDS)+LFAC*FWLUD
          SW(IJKE)=SW(IJKE)-GDS(IU)*(FWUDS-FWCDS)-LFAC*FWLUD
! LUD
!
!.....LUD
!
!
        END DO
       END DO
      END DO
!
      END DO
!
!=========================================================
!.....FLUXES THROUGH INTERNAL NORTH CV FACES
!=========================================================
!     F2(IJK) is the mass flux through the north face (outward
!     normal directed to N); FY(J) is the ratio of distance
!     from P to cell face, to distance from P to N; IJK
!     denotes node P and IJKN node N.
!     Contribution of convective and diffusive fluxes from
!     north face to AN(P), AS(N), and source terms at both
!     P and N are calculated; contributions to AP(P) and
!     AP(N) come through the sum of neighbor coefficients
!     and are not explicitly calculated.
!----------------------------------------------------------
! 3D
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
       IF (BYE(B).EQ.NJM) THEN
         JE=BYE(B)-1
       ELSE
         JE=BYE(B)
       END IF
!
       DO J=BYS(B),JE
!
!.....INTERPOLATION FACTORS, DISTANCE FROM P TO N (SAME FOR ALL !     J)
!
        FYN =FY(J)
        FYP =1.-FYN
        DYPN=YC(J+1)-YC(J)
! 3D
        DO I=2,NIM
          IJK =LK(K)+LI(I)+J
          IJKN=IJK+1
!
!.....CELL FACE AREA S = DX*DZ
! 3D
          S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
!
!.....CALCULATE VISCOSITY
! TURB
          VISI=VIS(IJKN)*FY(J)+VIS(IJK)*(1.-FY(J))
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX (SAME FOR U AND !     V)
!
!         D=VISC*S/DYPN
! TURB
          D=VISI*S/DYPN
! NANO
          DNA1=VISI/((1.-PHIS)**2.5)
          D=DNA1*S/DYPN
!
!.....EXPLICIT CONVECTIVE FLUXES FOR UDS AND CDS
!
          IF (F2(IJK).LE.0.) THEN
            CN=F2(IJK)
          ELSE
            CN=0.
          END IF
          IF (F2(IJK).GT.0.) THEN
            CP=F2(IJK)
          ELSE
            CP=0.
          END IF
!
!.....LUD
! LUD
          FULUD=0.;FVLUD=0.;FWLUD=0.
          IF (I.GT.2.AND.I.LT.NIM-1) THEN
          IF (J.GT.2.AND.J.LT.NJM-1) THEN
          IF (K.GT.2.AND.K.LT.NKM-1) THEN
!
!.....EQS 5.83, VERSTEEG & MALALASEKERA. ALIFUE OF CELL IJ IS ALIFUW
!     OF CELL IJN, SO, NO NEED TO REPEAT.
! LUD
            IF (F2(IJK).LT.0.) THEN
              ALIFUN=0.
              ALIFVN=0.
              ALIFWN=0.
              SAIMUN=(U(IJK+1+1)-U(IJK+1))/((U(IJK+1)-U(IJK))+SMALL)
              SAIMVN=(V(IJK+1+1)-V(IJK+1))/((V(IJK+1)-V(IJK))+SMALL)
              SAIMWN=(W(IJK+1+1)-W(IJK+1))/((W(IJK+1)-W(IJK))+SMALL)
            END IF
! LUD
            IF (F2(IJK).GE.0.) THEN
              ALIFUN=1.
              ALIFVN=1.
              ALIFWN=1.
              SAIPUN=(U(IJK)-U(IJK-1))/((U(IJK+1)-U(IJK))+SMALL)
              SAIPVN=(V(IJK)-V(IJK-1))/((V(IJK+1)-V(IJK))+SMALL)
              SAIPWN=(W(IJK)-W(IJK-1))/((W(IJK+1)-W(IJK))+SMALL)
            END IF
! LUD2
            T1=(1.-ALIFUN)*SAIMUN-ALIFUN*SAIPUN
            T2=U(IJK+1)-U(IJK)
            FULUD=0.5*F2(IJK)*T1*T2
! LUD2
            T1=(1.-ALIFVN)*SAIMVN-ALIFVN*SAIPVN
            T2=V(IJK+1)-V(IJK)
            FVLUD=0.5*F2(IJK)*T1*T2
! LUD2
            T1=(1.-ALIFWN)*SAIMWN-ALIFWN*SAIPWN
            T2=W(IJK+1)-W(IJK)
            FWLUD=0.5*F2(IJK)*T1*T2
!
          END IF
          END IF
          END IF
!
!.....COEFFICIENTS AN(P) AND AS(N) DUE TO UDS
! 3D
          AN(IJK) = CN-D
          AS(IJKN)=-CP-D
! LUD
! 3D
          FUUDS=CP*U(IJK)+CN*U(IJKN)
          FVUDS=CP*V(IJK)+CN*V(IJKN)
          FWUDS=CP*W(IJK)+CN*W(IJKN)
          FUCDS=F2(IJK)*(U(IJKN)*FYN+U(IJK)*FYP)
          FVCDS=F2(IJK)*(V(IJKN)*FYN+V(IJK)*FYP)
          FWCDS=F2(IJK)*(W(IJKN)*FYN+W(IJK)*FYP)
!
!.....SOURCE TERM CONTRIBUTIONS AT P AND N DUE TO DEFERRED
!     CORRECTION
! 3D
! LUD
          SU(IJK) =SU(IJK) +GDS(IU)*(FUUDS-FUCDS)+LFAC*FULUD
          SU(IJKN)=SU(IJKN)-GDS(IU)*(FUUDS-FUCDS)-LFAC*FULUD
          SV(IJK) =SV(IJK) +GDS(IU)*(FVUDS-FVCDS)+LFAC*FVLUD
          SV(IJKN)=SV(IJKN)-GDS(IU)*(FVUDS-FVCDS)-LFAC*FVLUD
          SW(IJK) =SW(IJK) +GDS(IU)*(FWUDS-FWCDS)+LFAC*FWLUD
          SW(IJKN)=SW(IJKN)-GDS(IU)*(FWUDS-FWCDS)-LFAC*FWLUD
        END DO
       END DO
      END DO
!
      END DO
! 3D
!=========================================================
!.....FLUXES THROUGH INTERNAL TOP CV FACES
!=========================================================
!     F3(IJK) is the mass flux through the top face (outward
!     normal directed to T); FZ(K) is the ratio of distance
!     from P to cell face, to distance from P to T; IJK
!     denotes node P and IJKT node T.
!     Contribution of convective and diffusive fluxes from
!     top face to AT(P), AB(T), and source terms at both
!     P and T are calculated; contributions to AP(P) and
!     AP(T) come through the sum of neighbor coefficients
!     and are not explicitly calculated.
!----------------------------------------------------------
! MPI for U(IJK+NIJ)
!
! MPI
!
! MPI
!
! MPI
! 3D
! BLK
      DO B=1,3
!
      IF (BZE(B).EQ.NKM) THEN
        KE=BZE(B)-1
      ELSE
        KE=BZE(B)
      END IF
!
      DO K=BZS(B),KE
!
       FZT=FZ(K)
       FZP=1.-FZT
       DZPT=ZC(K+1)-ZC(K)
!
       DO J=BYS(B),BYE(B)
!
        DO I=2,NIM
          IJK =LK(K)+LI(I)+J
          IJKT=IJK+NIJ
!
!.....CELL FACE AREA S = DX*RN*1
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
!
!.....CALCULATE VISCOSITY
! TURB
          VISI=VIS(IJKT)*FZ(K)+VIS(IJK)*(1.-FZ(K))
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX (SAME FOR U AND V)
!
!         D=VISC*S/DZPT
! TURB
          D=VISI*S/DZPT
! NANO
          DNA1=VISI/((1.-PHIS)**2.5)
          D=DNA1*S/DZPT
!
!.....EXPLICIT CONVECTIVE FLUXES FOR UDS AND CDS
!
          IF (F3(IJK).LE.0.) THEN
            CT=F3(IJK)
          ELSE
            CT=0.
          END IF
          IF (F3(IJK).GT.0.) THEN
            CP=F3(IJK)
          ELSE
            CP=0.
          END IF
!
!.....LUD
!
      FULUD=0.;FVLUD=0.;FWLUD=0.
      IF (I.GT.2.AND.I.LT.NIM-1) THEN
      IF (J.GT.2.AND.J.LT.NJM-1) THEN
      IF (K.GT.2.AND.K.LT.NKM-1) THEN
!
!.....EQS 5.83, VERSTEEG & MALALASEKERA. ALIFUE OF CELL IJ IS ALIFUW
!     OF CELL IJN, SO, NO NEED TO REPEAT.
! LUD
        IF (F3(IJK).LT.0.) THEN
          ALIFUT=0.
          ALIFVT=0.
          ALIFWT=0.
          SAIMUT=(U(IJK+NIJ+NIJ)-U(IJK+NIJ))/((U(IJK+NIJ)-U(IJK))+SMALL)
          SAIMVT=(V(IJK+NIJ+NIJ)-V(IJK+NIJ))/((V(IJK+NIJ)-V(IJK))+SMALL)
          SAIMWT=(W(IJK+NIJ+NIJ)-W(IJK+NIJ))/((W(IJK+NIJ)-W(IJK))+SMALL)
        END IF
! LUD
        IF (F3(IJK).GE.0.) THEN
          ALIFUT=1.
          ALIFVT=1.
          ALIFWT=1.
          SAIPUT=(U(IJK)-U(IJK-NIJ))/((U(IJK+NIJ)-U(IJK))+SMALL)
          SAIPVT=(V(IJK)-V(IJK-NIJ))/((V(IJK+NIJ)-V(IJK))+SMALL)
          SAIPWT=(W(IJK)-W(IJK-NIJ))/((W(IJK+NIJ)-W(IJK))+SMALL)
        END IF
! LUD
        T1=(1.-ALIFUT)*SAIMUT-ALIFUT*SAIPUT
        T2=U(IJK+NIJ)-U(IJK)
        FULUD=0.5*F3(IJK)*T1*T2
! LUD
        T1=(1.-ALIFVT)*SAIMVT-ALIFVT*SAIPVT
        T2=V(IJK+NIJ)-V(IJK)
        FVLUD=0.5*F3(IJK)*T1*T2
! LUD
        T1=(1.-ALIFWT)*SAIMWT-ALIFWT*SAIPWT
        T2=W(IJK+NIJ)-W(IJK)
        FWLUD=0.5*F3(IJK)*T1*T2
! LUD
      END IF
      END IF
      END IF
!
!.....COEFFICIENTS AE(P) AND AW(E) DUE TO UDS
!
          AT(IJK) = CT-D
          AB(IJKT)=-CP-D
!
          FUUDS=CP*U(IJK)+CT*U(IJKT)
          FVUDS=CP*V(IJK)+CT*V(IJKT)
          FWUDS=CP*W(IJK)+CT*W(IJKT)
          FUCDS=F3(IJK)*(U(IJKT)*FZT+U(IJK)*FZP)
          FVCDS=F3(IJK)*(V(IJKT)*FZT+V(IJK)*FZP)
          FWCDS=F3(IJK)*(W(IJKT)*FZT+W(IJK)*FZP)
!
!.....SOURCE TERM CONTRIBUTIONS AT P AND N DUE TO DEFERRED
!     CORRECTION
! 3D
! LUD
          SU(IJK) =SU(IJK) +GDS(IU)*(FUUDS-FUCDS)+LFAC*FULUD
          SU(IJKT)=SU(IJKT)-GDS(IU)*(FUUDS-FUCDS)-LFAC*FULUD
          SV(IJK) =SV(IJK) +GDS(IU)*(FVUDS-FVCDS)+LFAC*FVLUD
          SV(IJKT)=SV(IJKT)-GDS(IU)*(FVUDS-FVCDS)-LFAC*FVLUD
          SW(IJK) =SW(IJK) +GDS(IU)*(FWUDS-FWCDS)+LFAC*FWLUD
          SW(IJKT)=SW(IJKT)-GDS(IU)*(FWUDS-FWCDS)-LFAC*FWLUD
!
        END DO
       END DO
      END DO
!
      END DO
!
! MPI for AB(IJK+NIJ)
!
!=============================================================
!.....VOLUME INTEGRALS (SOURCE TERMS)
!=============================================================
!     Cell-face pressure calculated using linear interpolation;
!     cell volume is VOL, RP is the radius at node P; DX and DY
!     are the width and height of the cell. Contribution to AP
!     coefficient from volume integrals is stored temporarily
!     in arrays APU and APV for U and V, respectively; these
!     arrays are later used to store 1/AP, which is needed in
!     the pressure-correction equation.
!--------------------------------------------------------------
! BLK
      DO B=1,3
!
! MPI for P(IJK-NIJ)
!
      DO K=BZS(B),BZE(B)
        DZ=Z(K)-Z(K-1)
!
        DO I=2,NIM
          DX=X(I)-X(I-1)
!
          DO J=BYS(B),BYE(B)
            DY=Y(J)-Y(J-1)
!
            IJK=LK(K)+LI(I)+J
!
!...... CELL-FACE PRESSURE, CELL-CENTER GRADIENT, SOURCE
! 3D
            PE=P(IJK+NJ)*FX(I)+P(IJK)*(1.-FX(I))
            PW=P(IJK)*FX(I-1)+P(IJK-NJ)*(1.-FX(I-1))
!
            PN=P(IJK+1)*FY(J)+P(IJK)*(1.-FY(J))
! BLK
            IF (B.EQ.1.AND.J.EQ.BYS(B)) THEN
              IK=(K-1)*NI+I
              PS=PS1(IK)
            ELSE
              PS=P(IJK)*FY(J-1)+P(IJK-1)*(1.-FY(J-1))
            END IF
!
            PT=P(IJK+NIJ)*FZ(K)+P(IJK)*(1.-FZ(K))
! BLK
            IF (B.EQ.2.AND.K.EQ.BZS(B)) THEN
              IJ =(I-1)*NJ+J
              PB=PB2(IJ)
            ELSE
              PB=P(IJK)*FZ(K-1)+P(IJK-NIJ)*(1.-FZ(K-1))
            END IF
!
            DPX(IJK)=(PE-PW)/DX
            DPY(IJK)=(PN-PS)/DY
            DPZ(IJK)=(PT-PB)/DZ
! TURB
            SU(IJK)=SU(IJK)-DPX(IJK)*VOL(IJK)
            SV(IJK)=SV(IJK)-DPY(IJK)*VOL(IJK)
            SW(IJK)=SW(IJK)-DPZ(IJK)*VOL(IJK)
!
          END DO
        END DO
      END DO
!
      END DO
!
!=============================================================
!.....PROBLEM MODIFICATIONS - BOUNDARY CONDITIONS
!=============================================================
! 3D
      CALL BCUVW
!
!=============================================================
!.....UNDER-RELAXATION, SOLVING EQUATION SYSTEM FOR U-VELOCITY
!=============================================================
! 3D
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
      DO I=2,NIM
!
      DO J=BYS(B),BYE(B)
!
        IJK=LK(K)+LI(I)+J
        AP(IJK)=(-AE(IJK)-AW(IJK)-AN(IJK)-AS(IJK)-AT(IJK)-&
               &AB(IJK)+APU(IJK))*URFU
        SU(IJK)=SU(IJK)+(1.-URF(IU))*AP(IJK)*U(IJK)
        APU(IJK)=1./AP(IJK)
      END DO
      END DO
      END DO
!
      END DO
!
! MPI for U(IJK+NIJ)
!
      CALL CGSTAB(U,IU)
!
!=============================================================
!.....UNDER-RELAXATION, SOLVING EQUATION SYSTEM FOR V-VELOCITY
!=============================================================
! 3D
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
      DO I=2,NIM
!
      DO J=BYS(B),BYE(B)
!
        IJK=LK(K)+LI(I)+J
        AP(IJK)=(-AE(IJK)-AW(IJK)-AN(IJK)-AS(IJK)-AT(IJK)-&
               &AB(IJK)+APV(IJK))*URFV
        SU(IJK)=SV(IJK)+(1.-URF(IV))*AP(IJK)*V(IJK)
        APV(IJK)=1./AP(IJK)
      END DO
      END DO
      END DO
!
      END DO
!
! MPI for V(IJK+NIJ)
!
      CALL CGSTAB(V,IV)
! 3D
!=============================================================
!.....UNDER-RELAXATION, SOLVING EQUATION SYSTEM FOR W-VELOCITY
!=============================================================
! 3D
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
      DO I=2,NIM
!
      DO J=BYS(B),BYE(B)
!
        IJK=LK(K)+LI(I)+J
        AP(IJK)=(-AE(IJK)-AW(IJK)-AN(IJK)-AS(IJK)-AT(IJK)-&
               &AB(IJK)+APW(IJK))*URFW
        SU(IJK)=SW(IJK)+(1.-URF(IW))*AP(IJK)*W(IJK)
        APW(IJK)=1./AP(IJK)
      END DO
      END DO
      END DO
!
      END DO
! 3D
! MPI for W(IJK+NIJ)
!
      CALL CGSTAB(W,IW)
!
      END SUBROUTINE CALCUVW
!
!##############################################################
      SUBROUTINE CALCP
!##############################################################
!     This routine assembles and solves the pressure-correction
!     equation. Cell-face values of velocity components, used
!     to calculate the mass fluxes, are obtained by linear
!     interpolation and then corrected by adding a term
!     proportional to the third derivative of pressure and
!     squared grid spacing, as described in Chap. 7, Sect.
!     7.5.3.
!--------------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
!
      INTEGER :: IJPREF
!
      DOUBLE PRECISION :: DPXE,DPYN,DPZT
      DOUBLE PRECISION :: DPXEL,DPYNL,DPZTL
      DOUBLE PRECISION :: VOLE,VOLN,VOLT
      DOUBLE PRECISION :: SUM,PPO
      DOUBLE PRECISION :: UEL,VNL,WTL
      DOUBLE PRECISION :: APUE,APVN,APWT
!     DOUBLE PRECISION :: UE,VN,WT
      DOUBLE PRECISION :: PPE,PPW,PPN,PPS,PPT,PPB
!--------------------------------------------------------------
! T
!
!.....SET INDICES, INITIALIZE FIELDS
! INLET OUTLET
      SUM=0.
!
!.....INITIALIZATION OF TEMPORARILY STORED VARIABLES
! INLET OUTLET
! 3D
      DO IJK=1,NIJK
        SU(IJK)=0.
        AP(IJK)=0.
      END DO
!
!============================================================
!.....EAST CV FACES (S - AREA, VOLE - VOLUME BETWEEN P AND E)
!============================================================
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
       DO I=2,NIM-1
!
!.....INTERPOLATION FACTORS, DISTANCE FROM P TO E (SAME FOR ALL J)
!
        DXPE=XC(I+1)-XC(I)
        FXE=FX(I)
        FXP=1.-FXE
!
        DO J=BYS(B),BYE(B)
!
          IJK=LK(K)+LI(I)+J
          IJKE=IJK+NJ
!
          S=(Y(J)-Y(J-1))*(Z(K)-Z(K-1))
          VOLE=DXPE*S
          D=DENSIT*S
! NANO
          D=DNA2*S
!
!.....INTERPOLATED CELL FACE QUANTITIES (PRESSURE GRAD., U AND
!     1/AP). Note: pressure gradient is interpolated midway
!     between P and E, since the gradient calculated at cell
!     face is second order accurate at that location; the
!     velocity is interpolated linearly, to achieve second order
!     accuracy at cell face center.
!
          DPXEL=0.5*(DPX(IJKE)+DPX(IJK))
          UEL=U(IJKE)*FXE+U(IJK)*FXP
          APUE=APU(IJKE)*FXE+APU(IJK)*FXP
!
!.....CELL FACE GRADIENT, VELOCITY AND MASS FLUX
!
          DPXE=(P(IJKE)-P(IJK))/DXPE
          UE=UEL-APUE*VOLE*(DPXE-DPXEL)
          F1(IJK)=D*UE
!
!.....COEFFICIENTS OF P' EQUATION, AE(P) AND AW(E)
!
          AE(IJK)=-D*APUE*S
          AW(IJKE)=AE(IJK)
!
        END DO
       END DO
      END DO
!
      END DO
!
!=============================================================
!.....NORTH CV FACES (S - AREA, VOLN - VOLUME BETWEEN P AND N)
!=============================================================
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
       IF (BYE(B).EQ.NJM) THEN
         JE=BYE(B)-1
       ELSE
         JE=BYE(B)
       END IF
!
       DO J=BYS(B),JE
!
        DYPN=YC(J+1)-YC(J)
        FYN=FY(J)
        FYP=1.-FYN
!
        DO I=2,NIM
          IJK=LK(K)+LI(I)+J
          IJKN=IJK+1
!
          S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
          VOLN=S*DYPN
          D=DENSIT*S
! NANO
          D=DNA2*S
!
!.....INTERPOLATED CELL-FACE QUANTITIES (PRESSURE GRAD., U AND
!     1/AP)
!
          DPYNL=0.5*(DPY(IJKN)+DPY(IJK))
          VNL=V(IJKN)*FYN+V(IJK)*FYP
          APVN=APV(IJKN)*FYN+APV(IJK)*FYP
!
!.....CELL-FACE GRADIENT, VELOCITY AND MASS FLUX
!
          DPYN=(P(IJKN)-P(IJK))/DYPN
          VN=VNL-APVN*VOLN*(DPYN-DPYNL)
          F2(IJK)=D*VN
!
!.....COEFFICIENTS OF P' EQUATION, AN(P) AND AS(N)
!
          AN(IJK)=-D*APVN*S
          AS(IJKN)=AN(IJK)
!
        END DO
       END DO
      END DO
!
      END DO
!
!=============================================================
!.....TOP CV FACES (S - AREA, VOLT - VOLUME BETWEEN P AND T)
!=============================================================
! 3D
!=========================================================
!.....FLUXES THROUGH INTERNAL TOP CV FACES
!=========================================================
!     F3(IJK) is the mass flux through the top face (outward
!     normal directed to T); FZ(K) is the ratio of distance
!     from P to cell face, to distance from P to T; IJK
!     denotes node P and IJKT node T.
!     Contribution of convective and diffusive fluxes from
!     top face to AT(P), AB(T), and source terms at both
!     P and T are calculated; contributions to AP(P) and
!     AP(T) come through the sum of neighbor coefficients
!     and are not explicitly calculated.
!----------------------------------------------------------
! MPI for DPZ(IJK-NIJ)
!
! MPI
!
! MPI
! BLK
      DO B=1,3
!
      IF (BZE(B).EQ.NKM) THEN
        KE=BZE(B)-1
      ELSE
        KE=BZE(B)
      END IF
!
      DO K=BZS(B),KE
!
       DZPT=ZC(K+1)-ZC(K)
       FZT=FZ(K)
       FZP=1.-FZT
!
       DO J=BYS(B),BYE(B)
!
        DO I=2,NIM
          IJK =LK(K)+LI(I)+J
          IJKT=IJK+NIJ
!
!.....CELL FACE AREA S = DX*RN*1
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
          VOLT=S*DZPT
          D=DENSIT*S
! NANO
          D=DNA2*S
!
!.....INTERPOLATED CELL-FACE QUANTITIES (PRESSURE GRAD., U AND
!     1/AP)
!
          DPZTL=0.5*(DPZ(IJKT)+DPZ(IJK))
          WTL=W(IJKT)*FZT+W(IJK)*FZP
          APWT=APW(IJKT)*FZT+APW(IJK)*FZP
!
!.....CELL-FACE GRADIENT, VELOCITY AND MASS FLUX
!
          DPZT=(P(IJKT)-P(IJK))/DZPT
          WT=WTL-APWT*VOLT*(DPZT-DPZTL)
          F3(IJK)=D*WT
!
!.....COEFFICIENTS OF P' EQUATION, AN(P) AND AS(N)
!
          AT(IJK)=-D*APWT*S
          AB(IJKT)=AT(IJK)
!
        END DO
       END DO
      END DO
!
      END DO
!
! MPI for AB(IJK+NIJ)
!
!===============================================================
!.....BOUNDARY CONDITIONS: PRESCRIBED MASS FLUXES, ZERO
!     CORRECTION (EQUIVALENT TO ZERO NORMAL GRADIENT FOR P';
!     COEFFICIENT FOR THE BOUNDARY NODE IS ZERO, NO SPECIAL
!     TREATMENT REQUIRED).
!===============================================================
!
!==============================================================
!.....BOUNDARY FLUX
!==============================================================
!
!.....INLET BOUNDARIES (MASS FLUXES PRESCRIBED IN ROUTINE
!     'BCIN')
! BLK
      DO I=2,NIM
        DO J=BYS(1),NJM
          IJ =LI(I)+J
          IJK=LK(2)+LI(I)+J
          SU(IJK)=SU(IJK)+FMI(IJ)
        END DO
      END DO
!
!.....EXTRAPOLATED VELOCITY AT OUTLET BOUNDARY, OUTLET MASS
!     FLUXES
! OUTLET
! BLK23
      FLOWO=0.
!
      DO I=2,NIM
        DO J=2,NJM
!
          IJ =LI(I)  +J
          IJK=LK(NKM)+LI(I)+J
!
          U(IJK+NIJ)=U(IJK)
          V(IJK+NIJ)=V(IJK)
          W(IJK+NIJ)=W(IJK)
! NANO
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
          FMO(IJ)=DNA2*S*W(IJK+NIJ)
          FLOWO=FLOWO+FMO(IJ)
!
        END DO
      END DO
!
!.....CORRECT MASS FLUX TO SATISFY GLOBAL MASS CONSERVATION &
!     ADD TO SOURCE
! OUTLET
! BLK
      FAC=FLOMAS/(FLOWO+SMALL)
!
      DO I=2,NIM
        DO J=2,NJM
!
          IJ =LI(I)  +J
          IJK=LK(NKM)+LI(I)+J
!
          FMO(IJ)   =FMO(IJ)   *FAC
!
          U(IJK+NIJ)=U(IJK+NIJ)*FAC
          V(IJK+NIJ)=V(IJK+NIJ)*FAC
          W(IJK+NIJ)=W(IJK+NIJ)*FAC
!
          SU(IJK)=SU(IJK)-FMO(IJ)
!
        END DO
      END DO
!
!===============================================================
!.....SORCE TERM AND COEFFICIENT OF NODE P
!===============================================================
!
! MPI for F3(IJK-NIJ)
!
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
      DO I=2,NIM
!
      DO J=BYS(B),BYE(B)
!
        IJK=LK(K)+LI(I)+J
        AP(IJK)=AP(IJK)-AE(IJK)-AW(IJK)-AN(IJK)-AS(IJK)-&
               &AT(IJK)-AB(IJK)
        SU(IJK)=SU(IJK)+F1(IJK-NJ)-F1(IJK)+F2(IJK-1)-F2(IJK)+&
               &F3(IJK-NIJ)-F3(IJK)
        SUM=SUM+SU(IJK)
        PP(IJK)=0.
      END DO
      END DO
      END DO
!
      END DO

!
!===============================================================
!.....SOLVE EQUATIONS SYSTEM FOR P' AND APPLY CORRECTIONS
!===============================================================
!
! MPI for PP(IJK+NIJ)
!
!     PRINT *,'         SUM, FLOMAS, FLOWO = ',SUM,FLOMAS,FLOWO
!     PRINT *,'         IJKMAX, U, V, W, TE, ED = ',&
!     &IJKMAX,U(IJKMAX),V(IJKMAX),W(IJKMAX),TE(IJKMAX),ED(IJKMAX)
!
      CALL CGS(PP,IP)
!
!.....CALCULATE PRESSURE CORRECTION AT BOUNDARIES
!
! MPI for PP(IJK+NIJ)
!
      CALL PBOUND(PP,PPS1,PPB2)
!
!.....VALUE OF P' AT REFERENCE LOCATION TO BE SUBTRACTED FROM
!     ALL P'
!
      IJPREF=LK(KPR)+LI(IPR)+JPR
      IF(IJPREF.GT.(0).AND.IJPREF.LT.(NIJK+1)) THEN
        PPO=PP(IJPREF)
      END IF
!
!.....CORRECT EAST MASS FLUXES
!
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
      DO I=2,NIM-1
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(I)+J
        F1(IJK)=F1(IJK)+AE(IJK)*(PP(IJK+NJ)-PP(IJK))
      END DO
      END DO
      END DO
!
      END DO

!
!.....CORRECT NORTH MASS FLUXES
!
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
      DO I=2,NIM
!
      IF (BYE(B).EQ.NJM) THEN
        JE=BYE(B)-1
      ELSE
        JE=BYE(B)
      END IF
!
      DO J=BYS(B),JE
        IJK=LK(K)+LI(I)+J
        F2(IJK)=F2(IJK)+AN(IJK)*(PP(IJK+1)-PP(IJK))
      END DO
      END DO
      END DO
!
      END DO
!
!.....CORRECT TOP MASS FLUXES
! MPI for PP(IJK+NIJ)
! BLK
      DO B=1,3
!
      IF (BZE(B).EQ.NKM) THEN
        KE=BZE(B)-1
      ELSE
        KE=BZE(B)
      END IF
!
      DO K=BZS(B),KE
!
      DO I=2,NIM
!
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(I)+J
        F3(IJK)=F3(IJK)+AT(IJK)*(PP(IJK+NIJ)-PP(IJK))
      END DO
      END DO
      END DO
!
      END DO
!
!.....CORRECT PRESSURE AND VELOCITIES AT CELL CENTER
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
        DZ=Z(K)-Z(K-1)
!
        DO I=2,NIM
          DX=X(I)-X(I-1)
!
          DO J=BYS(B),BYE(B)
            DY=Y(J)-Y(J-1)
!
            IJK=LK(K)+LI(I)+J
!
            PPE=PP(IJK+NJ) *FX(I)  +PP(IJK)    *(1.-FX(I))
            PPW=PP(IJK)    *FX(I-1)+PP(IJK-NJ) *(1.-FX(I-1))
            PPN=PP(IJK+1)  *FY(J)  +PP(IJK)    *(1.-FY(J))
! BLK
            IF (B.EQ.1.AND.J.EQ.BYS(B)) THEN
              IK=(K-1)*NI+I
              PPS=PPS1(IK)
            ELSE
              PPS=PP(IJK)*FY(J-1)+PP(IJK-1)*(1.-FY(J-1))
            END IF
!
            PPT=PP(IJK+NIJ)*FZ(K)  +PP(IJK)    *(1.-FZ(K))
! BLK
            IF (B.EQ.2.AND.K.EQ.BZS(B)) THEN
              IJ =(I-1)*NJ+J
              PPB=PPB2(IJ)
            ELSE
              PPB=PP(IJK)*FZ(K-1)+PP(IJK-NIJ)*(1.-FZ(K-1))
            END IF
!
            U(IJK)=U(IJK)-(PPE-PPW)*DY*DZ*APU(IJK)
            V(IJK)=V(IJK)-(PPN-PPS)*DX*DZ*APV(IJK)
            W(IJK)=W(IJK)-(PPT-PPB)*DX*DY*APW(IJK)
            P(IJK)=P(IJK)+URF(IP)*(PP(IJK)-PPO)
!
          END DO
        END DO
      END DO
!
      END DO
!
      END SUBROUTINE CALCP
! T
! TURB
!###########################################################
      SUBROUTINE CALCSC(FI,IFI)
!###########################################################
      USE GLOBAL
      IMPLICIT NONE
!
      DOUBLE PRECISION, DIMENSION(1:NIJK), INTENT(INOUT) :: FI
      INTEGER, INTENT(IN) :: IFI
!
      DOUBLE PRECISION :: URFI
      DOUBLE PRECISION :: FUDS,FCDS
! LUD
      DOUBLE PRECISION :: ALFFIE=0.,ALFFIN=0.,ALFFIT=0.
! LUD
      DOUBLE PRECISION :: FFILUD=0.
      DOUBLE PRECISION :: SAIMFIE=0.,SAIMFIN=0.,SAIMFIT=0.
      DOUBLE PRECISION :: SAIPFIE=0.,SAIPFIN=0.,SAIPFIT=0.
!----------------------------------------------------------
!
!.....INITIALIZATION OF TEMPORARILY STORED VARIABLES
!
      DO IJK=1,NIJK
        SU(IJK)=0.
        AP(IJK)=0.
      END DO
!
      URFI=1./URF(IFI)
!
!==========================================================
!.....FLUXES THROUGH INTERNAL EAST CV-FACES
!==========================================================
! 3D
! BLK
      DO B=1,3!MODE-A
!
      DO K=BZS(B),BZE(B)!MODE-A
!
!     DO K=2,NKM,1!MODE-B
!
       DO I=2,NIM-1,1
!
!.....INTERPOLATION FACTORS, DISTANCE FROM P TO E (SAME FOR ALL J)
!
        FXE =FX(I)
        FXP =1.-FXE
        DXPE=XC(I+1)-XC(I)
! 3D
        DO J=BYS(B),BYE(B)!MODE-A
!
!        DO J=2,NJM!MODE-B
!
          IJK=LK(K)+LI(I)+J
          IJKE=IJK+NJ
!
!.....CELL FACE AREA S = DY*RE*1
!
          S=(Y(J)-Y(J-1))*(Z(K)-Z(K-1))
!
!.....INTERPOLATED CELL FACE VALUES
! TURB
          VISI=VIS(IJKE)*FX(I)+VIS(IJK)*(1.-FX(I))-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
          IF(IFI.EQ.IEN) DCOEF=(VISC+VISI/SIGT)/PRANL
          IF(IFI.EQ.ITE) DCOEF=VISC+VISI/SIGTE
          IF(IFI.EQ.IED) DCOEF=VISC+VISI/SIGED
! TURB
          D=DCOEF*S/DXPE
! NANO
          DNA1=DCOEF/((1.-PHIS)**2.5)
          D=DNA1*PRR*DNA3*S/DXPE
!
!.....EXPLICIT CONVECTIVE FLUX FOR UDS AND CDS
!
          IF (F1(IJK).LT.0.) THEN
            CE=F1(IJK)
          ELSE
            CE=0.
          END IF
          IF (F1(IJK).GT.0.) THEN
            CP=F1(IJK)
          ELSE
            CP=0.
          END IF
!
!.....LUD
! LUD
          FFILUD=0.
          IF (I.GT.2.AND.I.LT.NIM-1) THEN
          IF (J.GT.2.AND.J.LT.NJM-1) THEN
          IF (K.GT.2.AND.K.LT.NKM-1) THEN
!
!.....EQS 5.83, VERSTEEG & MALALASEKERA. ALIFUE OF CELL IJ IS
!     ALIFUW OF CELL IJE, SO, NO NEED TO REPEAT.
! LUD
            IF (F1(IJK).LT.0.) THEN
              ALFFIE=0.
              SAIMFIE=(FI(IJK+NJ+NJ)-FI(IJK+NJ))/&
                    &((FI(IJK+NJ)-FI(IJK))+SMALL)
            END IF
! LUD
            IF (F1(IJK).GE.0.) THEN
              ALFFIE=1.
              SAIPFIE=(FI(IJK)-FI(IJK-NJ))/((FI(IJK+NJ)-FI(IJK))+SMALL)
            END IF
! LUD2
            T1=(1.-ALFFIE)*SAIMFIE-ALFFIE*SAIPFIE
            T2=FI(IJK+NJ)-FI(IJK)
            FFILUD=0.5*F1(IJK)*T1*T2
! LUD
          END IF
          END IF
          END IF
! TURB
          FUDS=CP*FI(IJK)+CE*FI(IJKE)
          FCDS=F1(IJK)*(FI(IJKE)*FXE+FI(IJK)*FXP)
!
!.....COEFFICIENTS AE(P) AND AW(E) DUE TO UDS
!
          AE(IJK) = CE-D
          AW(IJKE)=-CP-D
!
!.....SOURCE TERM CONTRIBUTIONS AT P AND E DUE TO DEFERRED
!     CORRECTION
! TURB
! LUD
          SU(IJK) =SU(IJK) +GDS(IFI)*(FUDS-FCDS)+LFAC*FFILUD
          SU(IJKE)=SU(IJKE)-GDS(IFI)*(FUDS-FCDS)-LFAC*FFILUD
        END DO
       END DO
      END DO
!
      END DO!MODE-A
!
!=========================================================
!.....FLUXES THROUGH INTERNAL NORTH CV FACES
!=========================================================
! 3D
! BLK
      DO B=1,3!MODE-A
!
      DO K=BZS(B),BZE(B)!MODE-A
!
!     DO K=2,NKM!MODE-B
!
      IF (BYE(B).EQ.NJM) THEN!MODE-A
        JE=BYE(B)-1!MODE-A
      ELSE!MODE-A
        JE=BYE(B)!MODE-A
      END IF!MODE-A
!
      DO J=BYS(B),JE!MODE-A
!
!     DO J=2,NJM-1,1!MODE-B
!
!.....INTERPOLATION FACTORS, DISTANCE FROM P TO N (SAME FOR ALL J)
!
         FYN =FY(J)
         FYP =1.-FYN
         DYPN=YC(J+1)-YC(J)
!
         DO I=2,NIM,1
          IJK =LK(K)+LI(I)+J
          IJKN=IJK+1
!
!.....CELL FACE AREA S = DX*RN*1
!
          S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
!
!.....INTERPOLATED CELL FACE VALUES
! TURB
          VISI=VIS(IJKN)*FY(J)+VIS(IJK)*(1.-FY(J))-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
          IF(IFI.EQ.IEN) DCOEF=(VISC+VISI/SIGT)/PRANL
          IF(IFI.EQ.ITE) DCOEF=VISC+VISI/SIGTE
          IF(IFI.EQ.IED) DCOEF=VISC+VISI/SIGED
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX (SAME FOR U AND V)
! TURB
          D=DCOEF*S/DYPN
! NANO
          DNA1=DCOEF/((1.-PHIS)**2.5)
          D=DNA1*PRR*DNA3*S/DYPN
!
!.....EXPLICIT CONVECTIVE FLUXES FOR UDS AND CDS
!
          IF (F2(IJK).LT.0.) THEN
            CN=F2(IJK)
          ELSE
            CN=0.
          END IF
          IF (F2(IJK).GT.0.) THEN
            CP=F2(IJK)
          ELSE
            CP=0.
          END IF
!
!.....LUD
! LUD
          FFILUD=0.
          IF (I.GT.2.AND.I.LT.NIM-1) THEN
          IF (J.GT.2.AND.J.LT.NJM-1) THEN
          IF (K.GT.2.AND.K.LT.NKM-1) THEN
!
!.....EQS 5.83, VERSTEEG & MALALASEKERA. ALIFUE OF CELL IJ IS
!     ALIFUW OF CELL IJN, SO, NO NEED TO REPEAT.
! LUD
            IF (F2(IJK).LT.0.) THEN
              ALFFIN=0.
              SAIMFIN=(FI(IJK+1+1)-FI(IJK+1))/&
                    &((FI(IJK+1)-FI(IJK))+SMALL)
            END IF
! LUD
            IF (F2(IJK).GE.0.) THEN
              ALFFIN=1.
              SAIPFIN=(FI(IJK)-FI(IJK-1))/((FI(IJK+1)-FI(IJK))+SMALL)
            END IF
! LUD2
            T1=(1.-ALFFIN)*SAIMFIN-ALFFIN*SAIPFIN
            T2=FI(IJK+1)-FI(IJK)
            FFILUD=0.5*F2(IJK)*T1*T2
! LUD
          END IF
          END IF
          END IF
! TURB
          FUDS=CP*FI(IJK)+CN*FI(IJKN)
          FCDS=F2(IJK)*(FI(IJKN)*FYN+FI(IJK)*FYP)
!
!.....COEFFICIENTS AE(P) AND AW(E) DUE TO UDS
!
          AN(IJK) = CN-D
          AS(IJKN)=-CP-D
!
!.....SOURCE TERM CONTRIBUTIONS AT P AND E DUE TO DEFERRED
!     CORRECTION
! TURB
! LUD
          SU(IJK) =SU(IJK) +GDS(IFI)*(FUDS-FCDS)+LFAC*FFILUD
          SU(IJKN)=SU(IJKN)-GDS(IFI)*(FUDS-FCDS)-LFAC*FFILUD
!
        END DO
       END DO
      END DO
!
      END DO!MODE-A
!
!=========================================================
!.....FLUXES THROUGH INTERNAL TOP CV FACES
!=========================================================
!     F3(IJK) is the mass flux through the top face (outward
!     normal directed to T); FZ(K) is the ratio of distance
!     from P to cell face, to distance from P to T; IJK
!     denotes node P and IJKT node T.
!     Contribution of convective and diffusive fluxes from
!     top face to AT(P), AB(T), and source terms at both
!     P and T are calculated; contributions to AP(P) and
!     AP(T) come through the sum of neighbor coefficients
!     and are not explicitly calculated.
!----------------------------------------------------------
! 3D
! BLK
      DO B=1,3!MODE-A
!
      IF (BZE(B).EQ.NKM) THEN!MODE-A
        KE=BZE(B)-1!MODE-A
      ELSE!MODE-A
        KE=BZE(B)!MODE-A
      END IF!MODE-A
!
      DO K=BZS(B),KE!MODE-A
!
!     DO K=2,NKM-1!MODE-B
!
       FZT=FZ(K)
       FZP=1.-FZT
       DZPT=ZC(K+1)-ZC(K)
!
       DO J=BYS(B),BYE(B)!MODE-A
!
!      DO J=2,NJM!MODE-B
!
        DO I=2,NIM
          IJK =LK(K)+LI(I)+J
          IJKT=IJK+NIJ
!
!.....CELL FACE AREA S = DX*RN*1
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
!
!.....INTERPOLATED CELL FACE VALUES
! TURB
          VISI=VIS(IJKT)*FZ(K)+VIS(IJK)*(1.-FZ(K))-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
          IF(IFI.EQ.IEN) DCOEF=(VISC+VISI/SIGT)/PRANL
          IF(IFI.EQ.ITE) DCOEF=VISC+VISI/SIGTE
          IF(IFI.EQ.IED) DCOEF=VISC+VISI/SIGED
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX (SAME FOR U AND V)
! TURB
          D=DCOEF*S/DZPT
! NANO
          DNA1=DCOEF/((1.-PHIS)**2.5)
          D=DNA1*PRR*DNA3*S/DZPT
!
!.....EXPLICIT CONVECTIVE FLUXES FOR UDS AND CDS
!
          IF (F3(IJK).LT.0.) THEN
            CT=F3(IJK)
          ELSE
            CT=0.
          END IF
          IF (F3(IJK).GT.0.) THEN
            CP=F3(IJK)
          ELSE
            CP=0.
          END IF
!
!.....LUD
! LUD
          FFILUD=0.
          IF (I.GT.2.AND.I.LT.NIM-1) THEN
          IF (J.GT.2.AND.J.LT.NJM-1) THEN
          IF (K.GT.2.AND.K.LT.NKM-1) THEN
!
!.....EQS 5.83, VERSTEEG & MALALASEKERA. ALIFUE OF CELL IJ IS
!     ALIFUW OF CELL !IJN, SO, NO NEED TO REPEAT.
! LUD
            IF (F3(IJK).LT.0.) THEN
              ALFFIT=0.
              SAIMFIT=(FI(IJK+NIJ+NIJ)-FI(IJK+NIJ))/&
                    &((FI(IJK+NIJ)-FI(IJK))+SMALL)
            END IF
! LUD
            IF (F3(IJK).GE.0.) THEN
              ALFFIT=1.
              SAIPFIT=(FI(IJK)-FI(IJK-NIJ))/&
                    &((FI(IJK+NIJ)-FI(IJK))+SMALL)
            END IF
! LUD
            T1=(1.-ALFFIT)*SAIMFIT-ALFFIT*SAIPFIT
            T2=FI(IJK+NIJ)-FI(IJK)
            FFILUD=0.5*F3(IJK)*T1*T2
! LUD
          END IF
          END IF
          END IF
! TURB
          FUDS=CP*FI(IJK)+CT*FI(IJKT)
          FCDS=F3(IJK)*(FI(IJKT)*FZT+FI(IJK)*FZP)
!
!.....COEFFICIENTS AE(P) AND AW(E) DUE TO UDS
!
          AT(IJK) = CT-D
          AB(IJKT)=-CP-D
!
!.....SOURCE TERM CONTRIBUTIONS AT P AND E DUE TO DEFERRED
!     CORRECTION
! TURB
! LUD
          SU(IJK) =SU(IJK) +GDS(IFI)*(FUDS-FCDS)+LFAC*FFILUD
          SU(IJKT)=SU(IJKT)-GDS(IFI)*(FUDS-FCDS)-LFAC*FFILUD
!
        END DO
       END DO
      END DO
!
      END DO!MODE-A
!
!=============================================================
!.....VOLUME INTEGRALS (SOURCE TERMS)
!=============================================================
!
!
!.....BOTTOM BOUNDARY (ISOTHERMAL INLET)
! INLET
! BLK1
      DO I=2,NIM
        DO J=BYS(1),NJM
!
          IJ =LI(I)+J
          IJK=LK(2)+LI(I)+J
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
!         D=VISC*PRR*S/(ZC(2)-ZC(1))
!
!.....INTERPOLATED CELL FACE VALUES. VIS(IJK-NIJ)
!     VALUES SET IN MAIN PROGRAM ABOVE
! TURB
          VISI=VIS(IJK-NIJ)-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
          IF(IFI.EQ.IEN) DCOEF=(VISC+VISI/SIGT)/PRANL
          IF(IFI.EQ.ITE) DCOEF=VISC+VISI/SIGTE
          IF(IFI.EQ.IED) DCOEF=VISC+VISI/SIGED
!
          D=DCOEF*S/(ZC(2)-ZC(1))
! NANO
          DNA1=DCOEF/((1.-PHIS)**2.5)
          D=DNA1*DNA3*S/(ZC(2)-ZC(1))
! NANO
          IF(IFI.EQ.IEN) D=DNA1*DNA3*S/(ZC(2)-ZC(1))
          IF(IFI.EQ.ITE) D=DNA1*S/(ZC(2)-ZC(1))
          IF(IFI.EQ.IED) D=DNA1*S/(ZC(2)-ZC(1))
!
          IF (FMI(IJ).GT.0.) THEN
            CP=FMI(IJ)
          ELSE
            CP=0.
          END IF
!
          CB=-D-CP
! TURB
          AP(IJK)=AP(IJK)-CB
          SU(IJK)=SU(IJK)-CB*FI(IJK-NIJ)
!
        END DO
      END DO
!
!.....OUTLET (TOP) BOUNDARIES (CONSTANT GRADIENT BETWEEN
!     BOUNDARY & CV-CENTER ASSUMED). SHEAR FORCE IN X,Y-DIR,
!     NORMAL FORCE IN Z-DIR).
! 3D
! OUTLET
! BLK23
      DO I=2,NIM
        DO J=2,NJM
          IJ =LI(I)+J
          IJK=LK(NKM)+LI(I)+J
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
!         D=VISC*PRR*S/(ZC(NK)-ZC(NKM))
!
!.....INTERPOLATED CELL FACE VALUES. VIS(IJK+NIJ)
!     VALUES SET IN MAIN PROGRAM ABOVE (INITIAL
!     VARIABLE VALUES)
! TURB
          VISI=VIS(IJK+NIJ)-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
          IF(IFI.EQ.IEN) DCOEF=(VISC+VISI/SIGT)/PRANL
          IF(IFI.EQ.ITE) DCOEF=VISC+VISI/SIGTE
          IF(IFI.EQ.IED) DCOEF=VISC+VISI/SIGED
!
          D=DCOEF*S/(ZC(NK)-ZC(NKM))
!
          IF (FMO(IJ).LE.0.) THEN
            CT=FMO(IJ)
          ELSE
            CT=0.
          END IF
!
          CB=-D+CT
! TURB
          AP(IJK)=AP(IJK)-CB
          SU(IJK)=SU(IJK)-CB*FI(IJK+NIJ)
!
        END DO
      END DO
!
!=============================================================
!.....PROBLEM MODIFICATIONS - BOUNDARY CONDITIONS
!=============================================================
! TURB
!.....WALL BOUNDARY CONDITIONS AND SOURCES FOR TEMPERATURE
!
      IF(IFI.EQ.IEN) CALL BCT
! TURB
!.....WALL BOUNDARY CONDITIONS AND SOURCES FOR TURB. KIN. ENRGY
!
      IF(IFI.EQ.ITE) CALL KINE
! TURB
!.....WALL BOUNDARY CONDITIONS AND SOURCES FOR DISSIPATION RATE
!
      IF(IFI.EQ.IED) CALL DISE
!
!==============================================================
!.....UNDER-RELAXATION, SOLVING EQUATION SYSTEM FOR TEMPERATURE
!==============================================================
! 3D
! BLK
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
!
      DO I=2,NIM
!
      DO J=BYS(B),BYE(B)
!
          IJK=LK(K)+LI(I)+J
          AP(IJK)=(AP(IJK)-AW(IJK)-AE(IJK)-AN(IJK)-AS(IJK)-&
                 &AT(IJK)-AB(IJK))*URFI
          SU(IJK)=SU(IJK)+(1.-URF(IFI))*AP(IJK)*FI(IJK)
      END DO
      END DO
      END DO
!
      END DO
!
      CALL CGSTAB(FI,IFI)
!
!.....TOP BOUNDARY (OUTLET)
! OUTLET
! BLK23
      DO I=2,NIM
        DO J=2,NJM
          IJK=LK(NK)+LI(I)+J
          FI(IJK)=FI(IJK-NIJ)
        END DO
      END DO
!
      END SUBROUTINE CALCSC
! MPI
!###############################################################
      SUBROUTINE CGSTAB(FI,IFI)
!###############################################################
!    This routine incorporates the CGSTAB solver for seven-
!    diagonal, non-symmetric coefficient matrices (suitable for
!    convection/diffusion problems). See Sect. 5.3.7 for
!    details. Array index IJK calculated from indices I, J, and
!    K according to Table 3.1.
!---------------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
! MPI
!
! 3D
! MPI
      DOUBLE PRECISION, DIMENSION(1:NIJK), INTENT(INOUT) :: FI
      INTEGER, INTENT(IN) :: IFI
! 3D
! MPI
      INTEGER :: L
      DOUBLE PRECISION :: ALF,RES0,BET,RESL,RSM
      DOUBLE PRECISION :: BETO,GAM,OMG,UKRESO,VRES,VV
! 3D
      DOUBLE PRECISION :: DD(1:NIJK)
      DOUBLE PRECISION :: PK(1:NIJK),ZK(1:NIJK),RES(1:NIJK)
      DOUBLE PRECISION :: RESO(1:NIJK),UK(1:NIJK),VK(1:NIJK)
!---------------------------------------------------------------
!     RESMAX=1.E-10
!
!.....CALCULATE INITIAL RESIDUAL VECTOR
! MPI
      RES0=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            RES(IJK)=SU(IJK) - AP(IJK)*FI(IJK)-&
                 &AE(IJK)*FI(IJK+NJ) - AW(IJK)*FI(IJK-NJ) -&
                 &AN(IJK)*FI(IJK+1)  - AS(IJK)*FI(IJK-1)  -&
                 &AT(IJK)*FI(IJK+NIJ)- AB(IJK)*FI(IJK-NIJ)
            RES0=RES0+ABS(RES(IJK))
          END DO
        END DO
      END DO
      END DO
!
!.....CALCULATE ELEMENTS OF PRECONDITIONING MATRIX DIAGONAL
!
!      DO K=2,NKM
!        DO I=2,NIM
!          DO J=2,NJM
!            IJK=LK(K)+LI(I)+J
!            DD(IJK)=1./(AP(IJK)-AW(IJK)*DD(IJK-NJ) *AE(IJK-NJ)-&
!                               &AS(IJK)*DD(IJK-1)  *AN(IJK-1)-&
!                               &AB(IJK)*DD(IJK-NIJ)*AT(IJK-NIJ))
!          END DO
!        END DO
!      END DO
!
!.....INITIALIZE WORKING ARRAYS AND CONSTANTS
! MPI
      DO K=2,NKM
        DO I=2,NIM
          DO J=2,NJM
            IJK=LK(K)+LI(I)+J
            RESO(IJK)=RES(IJK)
            PK(IJK)=0.
            UK(IJK)=0.
            ZK(IJK)=0.
            VK(IJK)=0.
          END DO
        END DO
      END DO
      ALF=1.
      BETO=1.
      GAM=1.
!
!.....START INNER ITERATIONS
!
      DO L=1,NSW(IFI)
!
!..... CALCULATE BETA AND OMEGA
! MPI
      BET=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            BET=BET+RES(IJK)*RESO(IJK)
          END DO
        END DO
      END DO
      END DO
!
! MPI. NEED TO SUM UP BET
!
      OMG=BET*GAM/(ALF*BETO+1.E-30)
      BETO=BET
!
!..... CALCULATE PK
! MPI
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            PK(IJK)=RES(IJK)+OMG*(PK(IJK)-ALF*UK(IJK))
          END DO
        END DO
      END DO
      END DO
!
!.....SOLVE (M ZK = PK) - FORWARD SUBSTITUTION
!
!      DO K=2,NKM
!        DO I=2,NIM
!          DO J=2,NJM
!            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=(PK(IJK)-AW(IJK)*ZK(IJK-NJ)-&
!                            &AS(IJK)*ZK(IJK-1)-&
!                            &AB(IJK)*ZK(IJK-NI*NJ))*DD(IJK)
!          END DO
!        END DO
!      END DO
!
!      DO K=2,NKM
!        DO I=2,NIM
!          DO J=2,NJM
!            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=ZK(IJK)/(DD(IJK)+1.E-30)
!          END DO
!        END DO
!      END DO
!
!.....BACKWARD SUBSTITUTION
! MPI
! BLK
      DO B=1,3
      DO K=BZE(B),BZS(B),-1
        DO I=NIM,2,-1
          DO J=BYE(B),BYS(B),-1
            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=(ZK(IJK)-AE(IJK)*ZK(IJK+NJ)-&
!                            &AN(IJK)*ZK(IJK+1)-&
!                            &AT(IJK)*ZK(IJK+NIJ))*DD(IJK)
            ZK(IJK)=RESO(IJK)
          END DO
        END DO
      END DO
      END DO
! DEBUG
!.....CALCULATE UK = A.PK
! MPI for ZK(IJK+NIJ)
!
! MPI for ZK(IJK-NIJ)
!
!.....WEST BOUNDARY (PERI, SHEAR FORCE IN Y-DIR, DU/DX=0)
! MPI
! PERI CLEARED FOR PREX
!
!.....EAST BOUNDARY (PERI, SHEAR FORCE IN Y-DIR, DU/DX=0)
! MPI
! PERI CLEARED FOR PREX
!
! MPI
!
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            UK(IJK)=AP(IJK)*ZK(IJK)   +AE(IJK)*ZK(IJK+NJ)+&
                   &AW(IJK)*ZK(IJK-NJ)+AN(IJK)*ZK(IJK+1)+&
                   &AS(IJK)*ZK(IJK-1) +AT(IJK)*ZK(IJK+NIJ)+&
                   &AB(IJK)*ZK(IJK-NIJ)
          END DO
        END DO
      END DO
      END DO
!
!.....CALCULATE SCALAR PRODUCT UK.RESO AND GAMMA
! MPI
      UKRESO=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            UKRESO=UKRESO+UK(IJK)*RESO(IJK)
          END DO
        END DO
      END DO
      END DO
!
! MPI. NEED TO SUM UP UKRESO
!
      GAM=BET/(UKRESO+1.E-30)
!
!     PRINT *,'L,IFI,GAM,BET,UKRESO=',L,IFI,GAM,BET,UKRESO
!
!
!.....UPDATE (FI) AND CALCULATE W (OVERWRITE RES - IT IS RES-
!     UPDATE)
! MPI
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            FI(IJK)=FI(IJK)+GAM*ZK(IJK)
            RES(IJK)=RES(IJK)-GAM*UK(IJK)
          END DO
        END DO
      END DO
      END DO
!
!.....SOLVE (M Y = W); Y OVERWRITES ZK; FORWARD SUBSTITUTION
!
!      DO K=2,NKM
!        DO I=2,NIM
!          DO J=2,NJM
!            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=(RES(IJK)-AW(IJK)*ZK(IJK-NJ)-&
!                             &AS(IJK)*ZK(IJK-1)-&
!                             &AB(IJK)*ZK(IJK-NIJ))*DD(IJK)
!           END DO
!         END DO
!      END DO
!
!      DO K=2,NKM
!        DO I=2,NIM
!          DO J=2,NJM
!            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=ZK(IJK)/(DD(IJK)+1.E-30)
!          END DO
!        END DO
!      END DO
!
!.....BACKWARD SUBSTITUTION
! MPI
! BLK
      DO B=1,3
      DO K=BZE(B),BZS(B),-1
        DO I=NIM,2,-1
          DO J=BYE(B),BYS(B),-1
            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=(ZK(IJK)-AE(IJK)*ZK(IJK+NJ)-&
!                            &AN(IJK)*ZK(IJK+1)-&
!                            &AT(IJK)*ZK(IJK+NIJ))*DD(IJK)
            ZK(IJK)=RES(IJK)
          END DO
        END DO
      END DO
      END DO
!
!.....CALCULATE V = A.Y (VK = A.ZK)
! MPI for ZK(IJK+NIJ)
!
! MPI for ZK(IJK-NIJ)
!
! MPI
!
!.....WEST BOUNDARY (PERI, SHEAR FORCE IN Y-DIR, DU/DX=0)
! MPI
! PERI CLEARED FOR PREX
!
!.....EAST BOUNDARY (PERI, SHEAR FORCE IN Y-DIR, DU/DX=0)
! MPI
! PERI CLEARED FOR PREX
!
! MPI
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            VK(IJK)=AP(IJK)*ZK(IJK)    +AE(IJK)*ZK(IJK+NJ)+&
                   &AW(IJK)*ZK(IJK-NJ) +AN(IJK)*ZK(IJK+1)+&
                   &AS(IJK)*ZK(IJK-1)  +AT(IJK)*ZK(IJK+NIJ)+&
                   &AB(IJK)*ZK(IJK-NIJ)
          END DO
        END DO
      END DO
      END DO
!
!.....CALCULATE ALPHA (ALF)
! MPI
      VRES=0.
      VV=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            VRES=VRES+VK(IJK)*RES(IJK)
            VV=VV+VK(IJK)*VK(IJK)
          END DO
        END DO
      END DO
      END DO
!
! MPI. NEED TO SUM UP UKRESO
!
      ALF=VRES/(VV+1.E-30)
!
!.....UPDATE VARIABLE (FI) AND RESIDUAL (RES) VECTORS
! MPI
      RESL=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            FI(IJK)=FI(IJK)+ALF*ZK(IJK)
            RES(IJK)=RES(IJK)-ALF*VK(IJK)
            RESL=RESL+ABS(RES(IJK))
          END DO
        END DO
      END DO
      END DO
!
!....CHECK CONVERGENCE
!MPI
      IF (L.EQ.1) THEN
        RESOR(IFI)=RES0
      END IF
!
!MPI
!
      RSM=RESL/(RESOR(IFI)+1.E-20)
      IF (RSM.LT.SOR(IFI)) RETURN
!
!.....END OF ITERATION LOOP
!
      END DO
!
      END SUBROUTINE CGSTAB
!
!###############################################################
      SUBROUTINE CGS(FI,IFI)
!###############################################################
!    This routine incorporates the Incomplete Cholesky
!    preconditioned Conjugate Gradient solver for symmetric
!    matrices in 3D problems with seven-diagonal matrix
!    structure (see Sect. 5.3.6). Array index IJK converted from
!    3D indices I, J, and K according to Table 3.1. NS is the
!    number of inner iterations (sweeps).
!---------------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
! MPI
! 3D
! MPI
      DOUBLE PRECISION, DIMENSION(1:NIJK), INTENT(INOUT) :: FI
      INTEGER, INTENT(IN) :: IFI
! 3D
! MPI
      INTEGER :: L
      DOUBLE PRECISION :: ALF,RES0,S0,SK,BET
      DOUBLE PRECISION :: RSM,PKAPK,RESL
! 3D
      DOUBLE PRECISION :: DD(1:NIJK)
      DOUBLE PRECISION :: PK(1:NIJK),ZK(1:NIJK),RES(1:NIJK)
!---------------------------------------------------------------
!
!.....INITALIZE WORKING ARRAYS
! MPI
      DO K=1,NK
      DO I=1,NI
      DO J=1,NJ
        IJK=LK(K)+LI(I)+J
        PK(IJK)=0.
        ZK(IJK)=0.
!       DD(IJK)=0.
        RES(IJK)=0.
      END DO
      END DO
      END DO
!
!.....CALCULATE INITIAL RESIDUAL VECTOR AND THE NORM
! MPI
      RES0=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            RES(IJK)=SU(IJK)-AE(IJK)*FI(IJK+NJ)-&
                    &AW(IJK)*FI(IJK-NJ) -AN(IJK)*FI(IJK+1)-&
                    &AS(IJK)*FI(IJK-1)  -AT(IJK)*FI(IJK+NIJ)-&
                    &AB(IJK)*FI(IJK-NIJ)-AP(IJK)*FI(IJK)
            RES0=RES0+ABS(RES(IJK))
          END DO
        END DO
      END DO
      END DO
!
!.....CALCULATE ELEMENTS OF DIAGONAL PRECONDITIONING MATRIX
!
!     DO K=2,NKM
!       DO I=2,NIM
!         DO J=2,NJM
!           IJK=LK(K)+LI(I)+J
!            DD(IJK)=1./(AP(IJK)-AW(IJK)**2*DD(IJK-NJ)-&
!                               &AS(IJK)**2*DD(IJK-1)-&
!                               &AB(IJK)**2*DD(IJK-NI*NJ))
!         END DO
!       END DO
!     END DO
!
      S0=1.E20
!
!....START INNER ITERATIONS
!
      DO L=1,NSW(IFI)
!
!.....SOLVE FOR ZK(IJK) -- FORWARD SUBSTITUTION
!MPI
!      DO K=2,NKM
!        DO I=2,NIM
!          DO J=2,NJM
!            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=(RES(IJK)-AW(IJK)*ZK(IJK-NJ)-&
!                   &AS(IJK)*ZK(IJK-1)-&
!                   &AB(IJK)*ZK(IJK-NI*NJ))*DD(IJK)
!          END DO
!        END DO
!      END DO
!
!      DO K=2,NKM
!        DO I=2,NIM
!          DO J=2,NJM
!            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=ZK(IJK)/(DD(IJK)+1.E-30)
!          END DO
!        END DO
!      END DO
!
!.....BACKWARD SUBSTITUTION; CALCULATE SCALAR PRODUCT SK
! MPI
      SK=0.
! BLK
      DO B=1,3
      DO K=BZE(B),BZS(B),-1
        DO I=NIM,2,-1
          DO J=BYE(B),BYS(B),-1
            IJK=LK(K)+LI(I)+J
!            ZK(IJK)=(ZK(IJK)-AE(IJK)*ZK(IJK+NJ)-&
!                   &AN(IJK)*ZK(IJK+1)-&
!                   &AT(IJK)*ZK(IJK+NI*NJ))*DD(IJK)
            ZK(IJK)=RES(IJK)
            SK=SK+RES(IJK)*ZK(IJK)
          END DO
        END DO
      END DO
      END DO
!
!.....CALCULATE BETA AND NEW SEARCH VECTOR PK
! MPI
      IF (L.EQ.1) THEN
        PK=RES
      ELSE
        BET=SK/S0
! BLK
      DO B=1,3
        DO K=BZS(B),BZE(B)
          DO I=2,NIM
            DO J=BYS(B),BYE(B)
              IJK=LK(K)+LI(I)+J
              PK(IJK)=ZK(IJK)+BET*PK(IJK)
            END DO
          END DO
        END DO
      END DO
      END IF
!
!.....CALCULATE SCALAR PRODUCT (PK.A PK) AND ALPHA (OVERWRITE
!     ZK)
! MPI
!
! MPI for PK(IJK+NIJ)
!
! MPI for PK(IJK-NIJ)
!
! MPI
!
!.....WEST BOUNDARY (PERI, SHEAR FORCE IN Y-DIR, DU/DX=0)
! MPI
! PERI
!
!.....EAST BOUNDARY (PERI, SHEAR FORCE IN Y-DIR, DU/DX=0)
! MPI
! PERI
!
! MPI
      PKAPK=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            ZK(IJK)=AP(IJK)*PK(IJK)   +AE(IJK)*PK(IJK+NJ)+&
                   &AW(IJK)*PK(IJK-NJ)+AN(IJK)*PK(IJK+1)+&
                   &AS(IJK)*PK(IJK-1) +AT(IJK)*PK(IJK+NIJ)+&
                   &AB(IJK)*PK(IJK-NIJ)
            PKAPK=PKAPK+PK(IJK)*ZK(IJK)
          END DO
        END DO
      END DO
      END DO
! MPI
!
! MPI
      ALF=SK/PKAPK
!
!.....CALCULATE VARIABLE CORRECTION, NEW RESIDUAL VECTOR, AND
!     NORM
! MPI
      RESL=0.
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
        DO I=2,NIM
          DO J=BYS(B),BYE(B)
            IJK=LK(K)+LI(I)+J
            FI(IJK)=FI(IJK)+ALF*PK(IJK)
            RES(IJK)=RES(IJK)-ALF*ZK(IJK)
            RESL=RESL+ABS(RES(IJK))
          END DO
        END DO
      END DO
      END DO
! MPI
      S0=SK
!
!....CHECK CONVERGENCE
! MPI
      IF (L.EQ.1) THEN
        RESOR(IFI)=RES0
      END IF
!
! MPI
      RSM=RESL/(RESOR(IFI)+1.E-30)
      IF (RSM.LT.SOR(IFI)) RETURN
!
!.....END OF ITERATION LOOP
!
      END DO
!
      END SUBROUTINE CGS
!
!#############################################################
      SUBROUTINE PBOUND(FI,FIS1,FIB2)
!#############################################################
!     This routine calculates boundary values of pressure or
!     pressure-correction by extrapolating (linearly) from
!     inside.
!-------------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
!
      DOUBLE PRECISION, DIMENSION(1:NIJK), INTENT(INOUT) :: FI
      DOUBLE PRECISION, DIMENSION(1:NIK), INTENT(INOUT) :: FIS1
      DOUBLE PRECISION, DIMENSION(1:NIJ), INTENT(INOUT) :: FIB2
      INTEGER :: NJ2,NIJ2
!-------------------------------------------------------------
!
!.....SOUTH BOUNDARIES
! BLK1
      DO K=2,BZE(1)
      DO I=2,NIM
        IK =(K-1)*NI+I
        IJK=LK(K)+LI(I)+BYS(1)-1
        FIS1(IK)=FI(IJK+1)+(FI(IJK+1)-FI(IJK+2))*FY(BYS(1))
      END DO
      END DO
! BLK2
      DO K=BZS(2),NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+1
        FI(IJK)=FI(IJK+1)+(FI(IJK+1)-FI(IJK+2))*FY(2)
      END DO
      END DO
!
!.....NORTH BOUNDARIES
! BLK12
      DO K=2,NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+NJ
        FI(IJK)=FI(IJK-1)+(FI(IJK-1)-FI(IJK-2))*(1.-FY(NJM-1))
      END DO
      END DO
!
!.....WEST AND EAST BOUNDARIES
! BLK123
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
      NJ2=2*NJ
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(1)+J
        FI(IJK)=FI(IJK+NJ)+(FI(IJK+NJ)-FI(IJK+NJ2))*FX(2)
        IJK=LK(K)+LI(NI)+J
        FI(IJK)=FI(IJK-NJ)+(FI(IJK-NJ)-FI(IJK-NJ2))*&
               &(1.-FX(NIM-1))
      END DO
      END DO
!
      END DO
!
!.....BOTTOM BOUNDARIES
! BLK1
      DO I=2,NIM
      NIJ2=2*NIJ
      DO J=BYS(1),NJM
        IJK=LK(1)+LI(I)+J
        FI(IJK)=FI(IJK+NIJ)+(FI(IJK+NIJ)-FI(IJK+NIJ2))*FZ(2)
      END DO
      END DO
! BLK2
      DO I=2,NIM
      NIJ2=2*NIJ
      DO J=2,BYE(2)
        IJ =(I-1)*NJ+J
        IJK=LK(BZS(2)-1)+LI(I)+J
        FIB2(IJ)=FI(IJK+NIJ)+(FI(IJK+NIJ)-FI(IJK+NIJ2))*FZ(BZS(2))
      END DO
      END DO
!
!.....TOP BOUNDARY
! BLK
      DO I=2,NIM
      NIJ2=2*NIJ
      DO J=2,NJM
        IJK=LK(NK)+LI(I)+J
        FI(IJK)=FI(IJK-NIJ)+(FI(IJK-NIJ)-FI(IJK-NIJ2))*&
               &(1.-FZ(NKM-1))
      END DO
      END DO
!
      END SUBROUTINE PBOUND
!
!#############################################################
      SUBROUTINE BCUVW
!#############################################################
!     In this routine, boundary conditions for U and V equations
!     are implemented, i.e. fluxes through boundary cell faces
!     are approximated. Here, the boundaries encountered in
!     cavity flows are considered; inlet and outlet boundaries
!     require different treatment, see Sect. 7.7.
!-------------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
!----------------------------------------------------------
!
!.....SOUTH BOUNDARY (WALL; SHEAR FORCE IN X,Z-DIR, DV/DY=0)
! 3D
! BLK1
      DO K=2,BZE(1)
      DO I=2,NIM
        IK =(K-1)*NI+I
        IJK=LK(K)+LI(I)+BYS(1)
!
        S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
! TURB
        VISS=VISC
        IF(LCAL(ITE).AND.YPL(IJK).GT.CTRANS) VISS=VISW(IJK)
        D=VISS*S/(YC(BYS(1))-Y(BYS(1)-1))
! NANO
        DNA1=VISS/((1.-PHIS)**2.5)
        D=DNA1*S/(YC(BYS(1))-Y(BYS(1)-1))
!
        APU(IJK)=APU(IJK)+D
        SU(IJK) =SU(IJK)+D*US1(IK)!U(IJK-1)
!
        APW(IJK)=APW(IJK)+D
        SW(IJK) =SW(IJK)+D*WS1(IK)!W(IJK-1)
      END DO
      END DO
! BLK2
      DO K=BZS(2),NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+2
!
        S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
! TURB
        VISS=VISC
        IF(LCAL(ITE).AND.YPL(IJK).GT.CTRANS) VISS=VISW(IJK)
        D=VISS*S/(YC(2)-YC(1))
! NANO
        DNA1=VISS/((1.-PHIS)**2.5)
        D=DNA1*S/(YC(2)-YC(1))
!
        APU(IJK)=APU(IJK)+D
        SU(IJK) =SU(IJK)+D*U(IJK-1)
!
        APW(IJK)=APW(IJK)+D
        SW(IJK) =SW(IJK)+D*W(IJK-1)
      END DO
      END DO
!
!.....NORTH BOUNDARY (WALL, SHEAR FORCE IN X,Z-DIR, DV/DY=0)
! 3D
! BLK13
      DO K=2,NKM
      DO I=2,NIM
        IJK =LK(K)+LI(I)+NJM
!
        S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
! TURB
        VISS=VISC
        IF(LCAL(ITE).AND.YPL(IJK).GT.CTRANS) VISS=VISW(IJK)
        D=VISS*S/(YC(NJ)-YC(NJM))
! NANO
        DNA1=VISS/((1.-PHIS)**2.5)
        D=DNA1*S/(YC(NJ)-YC(NJM))
!
        APU(IJK)=APU(IJK)+D
        SU(IJK) =SU(IJK) +D*U(IJK+1)
!
        APW(IJK)=APW(IJK)+D
        SW(IJK) =SW(IJK) +D*W(IJK+1)
      END DO
      END DO
!
!.....WEST BOUNDARY (WALL, SHEAR FORCE IN Y,Z-DIR, DU/DX=0)
! 3D
! BLK123
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
        DO J=BYS(B),BYE(B)
          IJK=LK(K)+LI(2)+J
!
          S=(Y(J)-Y(J-1))*(Z(K)-Z(K-1))
! TURB
          VISS=VISC
          IF(LCAL(ITE).AND.YPL(IJK).GT.CTRANS) VISS=VISW(IJK)
          D=VISS*S/(XC(2)-XC(1))
! NANO
          DNA1=VISS/((1.-PHIS)**2.5)
          D=DNA1*S/(XC(2)-XC(1))
!
          APV(IJK)=APV(IJK)+D
          SV(IJK) =SV(IJK) +D*V(IJK-NJ)
!
          APW(IJK)=APW(IJK)+D
          SW(IJK) =SW(IJK) +D*W(IJK-NJ)
!
!.....EAST BOUNDARY (WALL, SHEAR FORCE IN Y,Z-DIR, DU/DX=0)
! 3D
! BLK123
          IJK=LK(K)+LI(NIM)+J
!
          S=(Y(J)-Y(J-1))*(Z(K)-Z(K-1))
! TURB
          VISS=VISC
          IF(LCAL(ITE).AND.YPL(IJK).GT.CTRANS) VISS=VISW(IJK)
          D=VISS*S/(XC(NI)-XC(NIM))
! NANO
          DNA1=VISS/((1.-PHIS)**2.5)
          D=DNA1*S/(XC(NI)-XC(NIM))
!
          APV(IJK)=APV(IJK)+D
          SV(IJK) =SV(IJK) +D*V(IJK+NJ)
!
          APW(IJK)=APW(IJK)+D
          SW(IJK) =SW(IJK) +D*W(IJK+NJ)
        END DO
      END DO
!
      END DO
!
!.....BOTTOM BOUNDARY (INLET, SHEAR FORCE IN X,Y-DIR, NORMAL
!     FORCE IN Z-DIR)
! 3D
! INLET
! BLK1
      DO I=2,NIM
        DO J=BYS(1),NJM
!
          IJ =LI(I)+J
          IJK=LK(2)+LI(I)+J
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
! TURB
          VISS=VISC
          IF(LCAL(ITE).AND.YPL(IJK).GT.CTRANS) VISS=VISW(IJK)
          D=VISS*S/(ZC(2)-ZC(1))
! NANO
          DNA1=VISS/((1.-PHIS)**2.5)
          D=DNA1*S/(ZC(2)-ZC(1))
!
          IF (FMI(IJ).GT.0.) THEN
            CP=FMI(IJ)
          ELSE
            CP=0.
          END IF
!
          CB=-D-CP
!
          APU(IJK)=APU(IJK)-CB
          SU(IJK) =SU(IJK) -CB*U(IJK-NIJ)
!
          APV(IJK)=APV(IJK)-CB
          SV(IJK) =SV(IJK) -CB*V(IJK-NIJ)
!
          APW(IJK)=APW(IJK)-CB
          SW(IJK) =SW(IJK) -CB*W(IJK-NIJ)
!
        END DO
      END DO
!
!.....BOTTOM BOUNDARY (WALL, SHEAR FORCE IN X,Y-DIR, DW/DZ=0)
! BLK2
      DO I=2,NIM
        DO J=2,BYE(2)
!
          IJ =(I-1)*NJ+J
          IJK=LK(BZS(2))+LI(I)+J
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
! TURB
          VISS=VISC
          IF(LCAL(ITE).AND.YPL(IJK).GT.CTRANS) VISS=VISW(IJK)
          D=VISS*S/(ZC(BZS(2))-Z(BZS(2)-1))
! NANO
          DNA1=VISS/((1.-PHIS)**2.5)
          D=DNA1*S/(ZC(BZS(2))-Z(BZS(2)-1))
!
          APU(IJK)=APU(IJK)+D
          SU(IJK) =SU(IJK) +D*UB2(IJ)
!
          APV(IJK)=APV(IJK)+D
          SV(IJK) =SV(IJK) +D*VB2(IJ)
!
        END DO
      END DO
!
!.....OUTLET (TOP) BOUNDARIES (CONSTANT GRADIENT BETWEEN
!     BOUNDARY & CV-CENTER ASSUMED). SHEAR FORCE IN X,Y-DIR,
!     NORMAL FORCE IN Z-DIR).
! 3D
! OUTLET
! BLK23
      DO I=2,NIM
        DO J=2,NJM
          IJ =LI(I)+J
          IJK=LK(NKM)+LI(I)+J
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
          D=VIS(IJK+NIJ)*S/(ZC(NK)-ZC(NKM))
!
          IF (FMO(IJ).LE.0.) THEN
            CT=FMO(IJ)
          ELSE
            CT=0.
          END IF
!
          CB=-D+CT
!
          APU(IJK)=APU(IJK)-CB
          SU(IJK) =SU(IJK) -CB*U(IJK+NIJ)
!
          APV(IJK)=APV(IJK)-CB
          SV(IJK) =SV(IJK) -CB*V(IJK+NIJ)
!
          APW(IJK)=APW(IJK)-CB
          SW(IJK)=SW(IJK)  -CB*W(IJK+NIJ)
        END DO
      END DO
!
      END SUBROUTINE BCUVW
! T
!#############################################################
      SUBROUTINE BCT
!#############################################################
!     In this routine, boundary conditions for the temperature
!     equation are implemented, i.e. heat fluxes through the
!     boundary cell faces are calculated. Here, specified wall
!     temperature and adiabatic wall (zero heat flux) are
!     considered; treatment at symmetry planes is the same as
!     for an adiabatic wall, but inlet and outlet require
!     different treatment, see Sect. 7.7.
!-------------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
!-------------------------------------------------------------
!
!.....SOUTH BOUNDARY (ISOTHERMAL WALL)
! BLK1
      DO K=2,BZE(1)
      DO I=2,NIM
        IK =(K-1)*NI+I
        IJK=LK(K)+LI(I)+BYS(1)
!
        S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
!
!.....INTERPOLATED CELL FACE VALUES
! TURB
        VISI=VISW(IJK)-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
        DCOEF=(VISC+VISI/SIGT)/PRANL
        D=DCOEF*S/(YC(BYS(1))-Y(BYS(1)-1))
! NANO
        DNA1=DCOEF/((1.-PHIS)**2.5)
        D=DNA1*DNA3*S/(YC(BYS(1))-Y(BYS(1)-1))
!
        AP(IJK)=AP(IJK)+D
        SU(IJK)=SU(IJK)+D*TS1(IK)!T(IJK-1)
      END DO
      END DO
! BLK2
      DO K=BZS(2),NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+2
!
        S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
!
!.....INTERPOLATED CELL FACE VALUES
! TURB
        VISI=VISW(IJK)-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
        DCOEF=(VISC+VISI/SIGT)/PRANL
        D=DCOEF*S/(YC(2)-YC(1))
! NANO
        DNA1=DCOEF/((1.-PHIS)**2.5)
        D=DNA1*DNA3*S/(YC(2)-YC(1))
!
        AP(IJK)=AP(IJK)+D
        SU(IJK)=SU(IJK)+D*T(IJK-1)
      END DO
      END DO
!
!.....NORTH BOUNDARY (ISOTHERMAL WALL)
!
! BLK13
      DO K=2,NKM
        DO I=2,NIM
          IJK=LK(K)+LI(I)+NJM
!
          S=(X(I)-X(I-1))*(Z(K)-Z(K-1))
!
!.....INTERPOLATED CELL FACE VALUES
! TURB
          VISI=VISW(IJK)-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
          DCOEF=(VISC+VISI/SIGT)/PRANL
          D=DCOEF*S/(YC(NJ)-YC(NJM))
! NANO
          DNA1=DCOEF/((1.-PHIS)**2.5)
          D=DNA1*DNA3*S/(YC(NJ)-YC(NJM))
!
          AP(IJK)=AP(IJK)+D
          SU(IJK)=SU(IJK)+D*T(IJK+1)
        END DO
      END DO
!
!.....WEST BOUNDARY (ADIABATIC WALL)
! BLK123
      DO B=1,3
!
      DO K=BZS(B),BZE(B)
        DO J=BYS(B),BYE(B)
          IJK=LK(K)+LI(1)+J
          T(IJK)=T(IJK+NJ)
!
!.....EAST BOUNDARY (ADIABATIC WALL)
! BLK123
          IJK=LK(K)+LI(NI)+J
          T(IJK)=T(IJK-NJ)
        END DO
      END DO
!
      END DO
!
!.....BOTTOM BOUNDARY (ISOTHERMAL INLET), SEE CALCSC SUBROUTINE
! BLK1
!
!.....BOTTOM BOUNDARY (ISOTHERMAL WALL)
! BLK2
      DO I=2,NIM
        DO J=2,BYE(2)
!
          IJ =(I-1)*NJ+J
          IJK=LK(BZS(2))+LI(I)+J
!
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
!
!.....INTERPOLATED CELL FACE VALUES
! TURB
          VISI=VISW(IJK)-VISC
!
!.....COEFFICIENT RESULTING FROM DIFFUSIVE FLUX
! TURB
          DCOEF=(VISC+VISI/SIGT)/PRANL
          D=DCOEF*S/(ZC(BZS(2))-Z(BZS(2)-1))
! NANO
          DNA1=DCOEF/((1.-PHIS)**2.5)
          D=DNA1*DNA3*S/(ZC(BZS(2))-Z(BZS(2)-1))
!
          AP(IJK)=AP(IJK)+D
          SU(IJK)=SU(IJK)+D*TB2(IJ)!T(IJK-NIJ)
        END DO
      END DO
!
!.....TOP BOUNDARY (ZERO GRAD OUTLET), SEE CALCSC SUBROUTINE
! BLK23
!
      END SUBROUTINE BCT
! TURB
!###############################################################
      SUBROUTINE KINE
!###############################################################
!     This routine assembles the source terms (volume integrals)
!     and applies wall boundary conditions for the turbulent
!     kinetic energy equation.
!===============================================================
      USE GLOBAL
      IMPLICIT NONE
!
      DOUBLE PRECISION :: DUX,DUY,DUZ
      DOUBLE PRECISION :: DVX,DVY,DVZ
      DOUBLE PRECISION :: DWX,DWY,DWZ
!
      DOUBLE PRECISION :: G11,G12,G13
      DOUBLE PRECISION :: G21,G22,G23
      DOUBLE PRECISION :: G31,G32,G33
!
      DOUBLE PRECISION :: XTW,YTW,ZTW
!--------------------------------------------------------------
!
!.....VOLUMETRIC SOURCES OF THE MODELLED TURBULENT KINETIC ENERGY
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DZ=Z(K)-Z(K-1)
      DO I=2,NIM
        DX=X(I)-X(I-1)
        DO J=BYS(B),BYE(B)
          DY=Y(J)-Y(J-1)
          IJK=LK(K)+LI(I)+J
!
!.....CELL-FACE U,V,W, CELL-CENTER GRADIENT
!
          UE=U(IJK+NJ)*FX(I)+U(IJK)*(1.-FX(I))
          UW=U(IJK)*FX(I-1)+U(IJK-NJ)*(1.-FX(I-1))
          UN=U(IJK+1)*FY(J)+U(IJK)*(1.-FY(J))
          US=U(IJK)*FY(J-1)+U(IJK-1)*(1.-FY(J-1))
          UT=U(IJK+NIJ)*FZ(K)+U(IJK)*(1.-FZ(K))
          UB=U(IJK)*FZ(K-1)+U(IJK-NIJ)*(1.-FZ(K-1))

          VE=V(IJK+NJ)*FX(I)+V(IJK)*(1.-FX(I))
          VW=V(IJK)*FX(I-1)+V(IJK-NJ)*(1.-FX(I-1))
          VN=V(IJK+1)*FY(J)+V(IJK)*(1.-FY(J))
          VS=V(IJK)*FY(J-1)+V(IJK-1)*(1.-FY(J-1))
          VT=V(IJK+NIJ)*FZ(K)+V(IJK)*(1.-FZ(K))
          VB=V(IJK)*FZ(K-1)+V(IJK-NIJ)*(1.-FZ(K-1))
!
          WE=W(IJK+NJ)*FX(I)+W(IJK)*(1.-FX(I))
          WW=W(IJK)*FX(I-1)+W(IJK-NJ)*(1.-FX(I-1))
          WN=W(IJK+1)*FY(J)+W(IJK)*(1.-FY(J))
          WS=W(IJK)*FY(J-1)+W(IJK-1)*(1.-FY(J-1))
          WT=W(IJK+NIJ)*FZ(K)+W(IJK)*(1.-FZ(K))
          WB=W(IJK)*FZ(K-1)+W(IJK-NIJ)*(1.-FZ(K-1))
!
          DUX=(UE-UW)/DX
          DUY=(UN-US)/DY
          DUZ=(UT-UB)/DZ
!
          DVX=(VE-VW)/DX
          DVY=(VN-VS)/DY
          DVZ=(VT-VB)/DZ
!
          DWX=(WE-WW)/DX
          DWY=(WN-WS)/DY
          DWZ=(WT-WB)/DZ
!
!         EQ. (9.28) FERZIGER AND PERIC 2ND ED.
!
          G11=2.*DUX**2
          G12=(DUY+DVX)*DUY
          G13=(DUZ+DWX)*DUZ
          G21=(DVX+DUY)*DVX
          G22=2.*DVY**2
          G23=(DVZ+DWY)*DVZ
          G31=(DWX+DUZ)*DWX
          G32=(DWY+DVZ)*DWY
          G33=2.*DWZ**2
          GEN(IJK)=(VIS(IJK)-VISC)*(G11+G12+G13+G21+G22+G23+G31+G32+G33)
!
!         EQ. (9.28) FERZIGER AND PERIC 2ND ED.
!
          SU(IJK)=SU(IJK)+GEN(IJK)*VOL(IJK)
!
!         SECOND LAST PARAGRAPH, PP. 280, FERZIGER AND PERIC 2ND ED.
!
          AP(IJK)=AP(IJK)+DENSIT*VOL(IJK)*ED(IJK)/(TE(IJK)+SMALL)
!
        END DO
      END DO
      END DO
      END DO
!
!.....SOUTH BOUNDARY (WALL)
! BLK1
      DO K=2,BZE(1)
      DO I=2,NIM
        IK =(K-1)*NI+I
        IJK=LK(K)+LI(I)+BYS(1)
!
        SU(IJK)=SU(IJK)-GEN(IJK)*VOL(IJK)
        VISS=VISC
        IF(LCAL(ITE).AND.(YPL(IJK)).GT.CTRANS) VISS=VISW(IJK)
!
!       UNIT TANGENT VECTOR ON A WALL BASED ON VELOCITY COMPONENTS
!
        XTW=U(IJK)/(SQRT(U(IJK)**2+W(IJK)**2)+SMALL)
        ZTW=W(IJK)/(SQRT(U(IJK)**2+W(IJK)**2)+SMALL)
!
!       EQNS. (8.74) AND (8.76) FERZIGER AND PERIC 2ND ED.
!
! NANO
!
        DNA1=VISS/((1.-PHIS)**2.5)
        TAU=DNA1*((U(IJK)-US1(IK))*XTW+(W(IJK)-WS1(IK))*ZTW)/&
           &(YC(BYS(1))-Y(BYS(1)-1))
!
!       EQNS. (9.39) AND (9.38) FERZIGER AND PERIC 2ND ED.
!
        GEN(IJK)=ABS(TAU)*CMU25*SQRT(MAX(ZERO,TE(IJK)))/&
                &((YC(BYS(1))-Y(BYS(1)-1))*CAPPA)
        SU(IJK)=SU(IJK)+GEN(IJK)*VOL(IJK)
      END DO
      END DO
!
!.....SOUTH BOUNDARY (WALL)
! BLK2
      DO K=BZS(2),NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+2
!
        SU(IJK)=SU(IJK)-GEN(IJK)*VOL(IJK)
        VISS=VISC
        IF(LCAL(ITE).AND.(YPL(IJK)).GT.CTRANS) VISS=VISW(IJK)
!
!       UNIT TANGENT VECTOR ON A WALL BASED ON VELOCITY COMPONENTS
!
        XTW=U(IJK)/(SQRT(U(IJK)**2+W(IJK)**2)+SMALL)
        ZTW=W(IJK)/(SQRT(U(IJK)**2+W(IJK)**2)+SMALL)
!
!       EQNS. (8.74) AND (8.76) FERZIGER AND PERIC 2ND ED.
!
! NANO
!
        DNA1=VISS/((1.-PHIS)**2.5)
        TAU=DNA1*((U(IJK)-U(IJK-1))*XTW&
                &+(W(IJK)-W(IJK-1))*ZTW)/(YC(2)-YC(1))
!
!       EQNS. (9.39) AND (9.38) FERZIGER AND PERIC 2ND ED.
!
        GEN(IJK)=ABS(TAU)*CMU25*SQRT(MAX(ZERO,TE(IJK)))/&
                &((YC(2)-YC(1))*CAPPA)
        SU(IJK)=SU(IJK)+GEN(IJK)*VOL(IJK)
      END DO
      END DO
!
!.....NORTH BOUNDARY (WALL)
! BLK13
      DO K=2,NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+NJM
!
        SU(IJK)=SU(IJK)-GEN(IJK)*VOL(IJK)
        VISS=VISC
        IF(LCAL(ITE).AND.(YPL(IJK)).GT.CTRANS) VISS=VISW(IJK)
!
!     UNIT TANGENT VECTOR ON A WALL BASED ON VELOCITY COMPONENTS
!
        XTW=U(IJK)/(SQRT(U(IJK)**2+W(IJK)**2)+SMALL)
        ZTW=W(IJK)/(SQRT(U(IJK)**2+W(IJK)**2)+SMALL)
!
!     EQNS. (8.74) AND (8.76) FERZIGER AND PERIC 2ND ED.
!
! NANO
!
        DNA1=VISS/((1.-PHIS)**2.5)
        TAU=DNA1*((U(IJK+1)-U(IJK))*XTW&
                &+(W(IJK+1)-W(IJK))*ZTW)/(YC(NJ)-YC(NJM))
!
!     EQNS. (9.39) AND (9.38) FERZIGER AND PERIC 2ND ED.
!
        GEN(IJK)=ABS(TAU)*CMU25*SQRT(MAX(ZERO,TE(IJK)))/&
                &((YC(NJ)-YC(NJM))*CAPPA)
        SU(IJK)=SU(IJK)+GEN(IJK)*VOL(IJK)
      END DO
      END DO
!
!.....WEST BOUNDARY (WALL)
! BLK123
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(2)+J
!
        SU(IJK)=SU(IJK)-GEN(IJK)*VOL(IJK)
        VISS=VISC
        IF(LCAL(ITE).AND.(YPL(IJK)).GT.CTRANS) VISS=VISW(IJK)
!
! UNIT TANGENT VECTOR ON A WALL BASED ON VELOCITY COMPONENTS
!
        YTW=V(IJK)/(SQRT(V(IJK)**2+W(IJK)**2)+SMALL)
        ZTW=W(IJK)/(SQRT(V(IJK)**2+W(IJK)**2)+SMALL)
!
! EQNS. (9.39) AND (9.38) FERZIGER AND PERIC 2ND ED.
!
! NANO
!
        DNA1=VISS/((1.-PHIS)**2.5)
        TAU=DNA1*((V(IJK)-V(IJK-NJ))*YTW&
                &+(W(IJK)-W(IJK-NJ))*ZTW)/(XC(2)-XC(1))
        GEN(IJK)=ABS(TAU)*CMU25*SQRT(MAX(ZERO,TE(IJK)))/&
                &((XC(2)-XC(1))*CAPPA)
        SU(IJK)=SU(IJK)+GEN(IJK)*VOL(IJK)
      END DO
      END DO
      END DO
!
!.....EAST BOUNDARY (WALL)
! BLK123
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(NIM)+J
!
        SU(IJK)=SU(IJK)-GEN(IJK)*VOL(IJK)
        VISS=VISC
        IF(LCAL(ITE).AND.(YPL(IJK)).GT.CTRANS) VISS=VISW(IJK)
!
! UNIT TANGENT VECTOR ON A WALL BASED ON VELOCITY COMPONENTS
!
        YTW=V(IJK)/(SQRT(V(IJK)**2+W(IJK)**2)+SMALL)
        ZTW=W(IJK)/(SQRT(V(IJK)**2+W(IJK)**2)+SMALL)
!
! EQNS. (9.39) AND (9.38) FERZIGER AND PERIC 2ND ED.
!
! NANO
!
        DNA1=VISS/((1.-PHIS)**2.5)
        TAU=DNA1*((V(IJK+NJ)-V(IJK))*YTW&
                &+(W(IJK+NJ)-W(IJK))*ZTW)/(XC(NI)-XC(NIM))
        GEN(IJK)=ABS(TAU)*CMU25*SQRT(MAX(ZERO,TE(IJK)))/&
                &((XC(NI)-XC(NIM))*CAPPA)
        SU(IJK)=SU(IJK)+GEN(IJK)*VOL(IJK)
      END DO
      END DO
      END DO
!
!.....BOTTOM BOUNDARY (INLET), SEE MAIN PROGRAM AND CALCSC
!     SUBROUTINE
! BLK1
!
!.....BOTTOM BOUNDARY (WALL)
! BLK2
      DO I=2,NIM
      DO J=2,BYE(2)
!
        IJ =(I-1)*NJ+J
        IJK=LK(BZS(2))+LI(I)+J
!
        SU(IJK)=SU(IJK)-GEN(IJK)*VOL(IJK)
        VISS=VISC
        IF(LCAL(ITE).AND.(YPL(IJK)).GT.CTRANS) VISS=VISW(IJK)
!
! UNIT TANGENT VECTOR ON A WALL BASED ON VELOCITY COMPONENTS
!
        XTW=U(IJK)/(SQRT(U(IJK)**2+V(IJK)**2)+SMALL)
        YTW=V(IJK)/(SQRT(U(IJK)**2+V(IJK)**2)+SMALL)
!
!       EQNS. (8.74) AND (8.76) FERZIGER AND PERIC 2ND ED.
!
! NANO
!
        DNA1=VISS/((1.-PHIS)**2.5)
        TAU=DNA1*((U(IJK)-UB2(IJ))*XTW&
                &+(V(IJK)-VB2(IJ))*YTW)/(ZC(BZS(2))-Z(BZS(2)-1))
!
! EQNS. (9.39) AND (9.38) FERZIGER AND PERIC 2ND ED.
!
        GEN(IJK)=ABS(TAU)*CMU25*SQRT(MAX(ZERO,TE(IJK)))/&
                &((ZC(BZS(2))-Z(BZS(2)-1))*CAPPA)
        SU(IJK)=SU(IJK)+GEN(IJK)*VOL(IJK)
      END DO
      END DO
!
!.....TOP BOUNDARY (OUTLET), SEE CALCSC SUBROUTINE
! BLK23
      RETURN
      END SUBROUTINE KINE
! TURB
!###############################################################
      SUBROUTINE DISE
!###############################################################
!     This routine assembles the source terms (volume integrals)
!     and applies wall boundary conditions for the dissipation
!     rate equation.
!===============================================================
      USE GLOBAL
      IMPLICIT NONE
!
      DOUBLE PRECISION :: TMP
!--------------------------------------------------------------
!
!.....VOLUMETRIC SOURCES FOR THE MODELLED DISSIPATION RATE
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DZ=Z(K)-Z(K-1)
      DO I=2,NIM
        DX=X(I)-X(I-1)
        DO J=BYS(B),BYE(B)
          DY=Y(J)-Y(J-1)
          IJK=LK(K)+LI(I)+J
!
! EQ. (9.30) FERZIGER AND PERIC 2ND ED.
!
          TMP=ED(IJK)*VOL(IJK)/(TE(IJK)+SMALL)
          SU(IJK)=SU(IJK)+CE1*TMP*GEN(IJK)
          AP(IJK)=AP(IJK)+CE2*TMP*DENSIT
        END DO
      END DO
      END DO
      END DO
!
!.....WALL BOUNDARIES APPROXIMATED WITH WALL FUNCTIONS.
!     FOR CORRECT VALUES OF DISSIPATION ALL COEFFICIENTS HAVE
!     TO BE ZERO, SU EQUAL THE DISSIPATION, AND AP = 1.
!
!.....SOUTH BOUNDARY (WALL)
! BLK1
      DO K=2,BZE(1)
      DO I=2,NIM
        IK =(K-1)*NI+I
        IJK=LK(K)+LI(I)+BYS(1)
!
!     EQ. (9.40) FERZIGER AND PERIC 2ND ED.
!
        ED(IJK)=CMU75*(MAX(ZERO,TE(IJK)))**1.5/&
               &(CAPPA*(YC(BYS(1))-Y(BYS(1)-1)))
        SU(IJK)=ED(IJK)
        AP(IJK)=1.
        AS(IJK)=0.
        AN(IJK)=0.
        AB(IJK)=0.
        AT(IJK)=0.
        AW(IJK)=0.
        AE(IJK)=0.
      END DO
      END DO
!
!.....SOUTH BOUNDARY (WALL)
! BLK2
      DO K=BZS(2),NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+2
!
!     EQ. (9.40) FERZIGER AND PERIC 2ND ED.
!
        ED(IJK)=CMU75*(MAX(ZERO,TE(IJK)))**1.5/(CAPPA*(YC(2)-YC(1)))
        SU(IJK)=ED(IJK)
        AP(IJK)=1.
        AS(IJK)=0.
        AN(IJK)=0.
        AB(IJK)=0.
        AT(IJK)=0.
        AW(IJK)=0.
        AE(IJK)=0.
      END DO
      END DO
!
!.....NORTH BOUNDARY (WALL)
! BLK13
      DO K=2,NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+NJM
!
!     EQ. (9.40) FERZIGER AND PERIC 2ND ED.
!
        ED(IJK)=CMU75*(MAX(ZERO,TE(IJK)))**1.5/(CAPPA*(YC(NJ)-YC(NJM)))
        SU(IJK)=ED(IJK)
        AP(IJK)=1.
        AS(IJK)=0.
        AN(IJK)=0.
        AB(IJK)=0.
        AT(IJK)=0.
        AW(IJK)=0.
        AE(IJK)=0.
      END DO
      END DO
!
!.....WEST BOUNDARY (WALL)
! BLK123
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(2)+J
!
!     EQ. (9.40) FERZIGER AND PERIC 2ND ED.
!
        ED(IJK)=CMU75*(MAX(ZERO,TE(IJK)))**1.5/(CAPPA*(XC(2)-XC(1)))
        SU(IJK)=ED(IJK)
        AP(IJK)=1.
        AS(IJK)=0.
        AN(IJK)=0.
        AB(IJK)=0.
        AT(IJK)=0.
        AW(IJK)=0.
        AE(IJK)=0.
      END DO
      END DO
      END DO
!
!.....EAST BOUNDARY (WALL)
! BLK123
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(NIM)+J
!
!     EQ. (9.40) FERZIGER AND PERIC 2ND ED.
!
        ED(IJK)=CMU75*(MAX(ZERO,TE(IJK)))**1.5/(CAPPA*(XC(NI)-XC(NIM)))
        SU(IJK)=ED(IJK)
        AP(IJK)=1.
        AS(IJK)=0.
        AN(IJK)=0.
        AB(IJK)=0.
        AT(IJK)=0.
        AW(IJK)=0.
        AE(IJK)=0.
      END DO
      END DO
      END DO
!
!.....BOTTOM BOUNDARY (INLET), SEE MAIN PROGRAM AND CALCSC
!     SUBROUTINE
! BLK1
!
!.....BOTTOM BOUNDARY (WALL)
! BLK2
      DO I=2,NIM
      DO J=2,BYE(2)
        IJK=LK(BZS(2))+LI(I)+J
!
!     EQ. (9.40) FERZIGER AND PERIC 2ND ED.
!
        ED(IJK)=CMU75*(MAX(ZERO,TE(IJK)))**1.5/&
               &(CAPPA*(ZC(BZS(2))-Z(BZS(2)-1)))
        SU(IJK)=ED(IJK)
        AP(IJK)=1.
        AS(IJK)=0.
        AN(IJK)=0.
        AB(IJK)=0.
        AT(IJK)=0.
        AW(IJK)=0.
        AE(IJK)=0.
      END DO
      END DO
!
!.....TOP BOUNDARY (OUTLET), SEE CALCSC SUBROUTINE
! BLK23
!
      RETURN
      END SUBROUTINE DISE
! TURB
!#############################################################
      SUBROUTINE MODVIS
!#############################################################
!     This routine calculates the eddy viscosity and the
!     dimensionless distance from the wall; also, an effective
!     wall viscosity is calculated so that the shear stress
!     can be computed using the same approximation as for
!     laminar flows (see Eqs. (9.37), (8.74) and (8.73)).
!=============================================================
      USE GLOBAL
      IMPLICIT NONE
!
      DOUBLE PRECISION :: CK,VISCW,VISOLD
!--------------------------------------------------------------
!
!.....EDDY VISCOSITY,DIMENSIONLESS WALL DISTANCE AND LOCAL REYNOLDS-
!     NUMBER
! BLK
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DZ=Z(K)-Z(K-1)
      DO I=2,NIM
        DX=X(I)-X(I-1)
        DO J=BYS(B),BYE(B)
          DY=Y(J)-Y(J-1)
          IJK=LK(K)+LI(I)+J
!
          VISOLD=VIS(IJK)
!
! EQ. (9.31) FERZIGER AND PERIC 2ND ED.
!
          VIS(IJK)=(VISC+CMU*DENSIT*TE(IJK)**2/&
                 &(ED(IJK)+SMALL))*URF(IVIS)&
                 &+VISOLD*(1.-URF(IVIS))
!
!         VIS(IJK)=VISC+0.001*(VIS(IJK)-VISC)
        END DO
      END DO
      END DO
      END DO
!
!.....SOUTH BOUNDARY (WALL)
! BLK1
      DO K=2,BZE(1)
      DO I=2,NIM
        IK =(K-1)*NI+I
        IJK=LK(K)+LI(I)+BYS(1)
!
! EQ. (9.36) FERZIGER AND PERIC 2ND ED.
!
        CK=CMU25*SQRT(MAX(ZERO,TE(IJK)))
!
! EQ. (9.35) FERZIGER AND PERIC 2ND ED.
!
        YPL(IJK)=DENSIT*CK*(YC(BYS(1))-Y(BYS(1)-1))/VISC
!
! EQNS. (9.34),(9.35),(9.36),(9.37) AND FOOTNOTE PP283
! FERZIGER AND PERIC 2ND ED.
!
        VISCW=YPL(IJK)*VISC*CAPPA/LOG(ELOG*YPL(IJK))
        VISW(IJK)=MAX(VISC,VISCW)
        VIS(IJK-1)=VISW(IJK)
!
      END DO
      END DO
!
!.....SOUTH BOUNDARY (WALL)
! BLK2
      DO K=BZS(2),NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+2
!
! EQ. (9.36) FERZIGER AND PERIC 2ND ED.
!
        CK=CMU25*SQRT(MAX(ZERO,TE(IJK)))
!
! EQ. (9.35) FERZIGER AND PERIC 2ND ED.
!
        YPL(IJK)=DENSIT*CK*(YC(2)-YC(1))/VISC
!
! EQNS. (9.34),(9.35),(9.36),(9.37) AND FOOTNOTE PP283
! FERZIGER AND PERIC 2ND ED.
!
        VISCW=YPL(IJK)*VISC*CAPPA/LOG(ELOG*YPL(IJK))
        VISW(IJK)=MAX(VISC,VISCW)
        VIS(IJK-1)=VISW(IJK)
!
      END DO
      END DO
!
!.....NORTH BOUNDARY (WALL)
! BLK13
      DO K=2,NKM
      DO I=2,NIM
        IJK=LK(K)+LI(I)+NJM
!
! EQ. (9.36) FERZIGER AND PERIC 2ND ED.
!
        CK=CMU25*SQRT(MAX(ZERO,TE(IJK)))
!
! EQ. (9.35) FERZIGER AND PERIC 2ND ED.
!
        YPL(IJK)=DENSIT*CK*(YC(NJ)-YC(NJM))/VISC
!
! EQNS. (9.34),(9.35),(9.36),(9.37) AND FOOTNOTE PP283
! FERZIGER AND PERIC 2ND ED.
!
        VISCW=YPL(IJK)*VISC*CAPPA/LOG(ELOG*YPL(IJK))
        VISW(IJK)=MAX(VISC,VISCW)
        VIS(IJK+1)=VISW(IJK)
!
      END DO
      END DO
!
!.....WEST BOUNDARY (WALL)
! BLK123
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(2)+J
!
! EQ. (9.36) FERZIGER AND PERIC 2ND ED.
!
        CK=CMU25*SQRT(MAX(ZERO,TE(IJK)))
!
! EQ. (9.35) FERZIGER AND PERIC 2ND ED.
!
        YPL(IJK)=DENSIT*CK*(XC(2)-XC(1))/VISC
!
! EQNS. (9.34),(9.35),(9.36),(9.37) AND FOOTNOTE PP283
! FERZIGER AND PERIC 2ND ED.
!
        VISCW=YPL(IJK)*VISC*CAPPA/LOG(ELOG*YPL(IJK))
        VISW(IJK)=MAX(VISC,VISCW)
        VIS(IJK-NJ)=VISW(IJK)
!
      END DO
      END DO
      END DO
!
!.....EAST BOUNDARY (WALL)
! BLK123
      DO B=1,3
      DO K=BZS(B),BZE(B)
      DO J=BYS(B),BYE(B)
        IJK=LK(K)+LI(NIM)+J
!
! EQ. (9.36) FERZIGER AND PERIC 2ND ED.
!
        CK=CMU25*SQRT(MAX(ZERO,TE(IJK)))
!
! EQ. (9.35) FERZIGER AND PERIC 2ND ED.
!
        YPL(IJK)=DENSIT*CK*(XC(NI)-XC(NIM))/VISC
!
! EQNS. (9.34),(9.35),(9.36),(9.37) AND FOOTNOTE PP283
! FERZIGER AND PERIC 2ND ED.
!
        VISCW=YPL(IJK)*VISC*CAPPA/LOG(ELOG*YPL(IJK))
        VISW(IJK)=MAX(VISC,VISCW)
        VIS(IJK+NJ)=VISW(IJK)
!
      END DO
      END DO
      END DO
!
!.....BOTTOM BOUNDARY (INLET), SEE CALCSC SUBROUTINE
! BLK1
!
!.....BOTTOM BOUNDARY (WALL)
! BLK2
      DO I=2,NIM
      DO J=2,BYE(2)
!
        IJ =(I-1)*NJ+J
        IJK=LK(BZS(2))+LI(I)+J
!
! EQ. (9.36) FERZIGER AND PERIC 2ND ED.
!
        CK=CMU25*SQRT(MAX(ZERO,TE(IJK)))
!
! EQ. (9.35) FERZIGER AND PERIC 2ND ED.
!
        YPL(IJK)=DENSIT*CK*(ZC(BZS(2))-Z(BZS(2)-1))/VISC
!
! EQNS. (9.34),(9.35),(9.36),(9.37) AND FOOTNOTE PP283
! FERZIGER AND PERIC 2ND ED.
!
        VISCW=YPL(IJK)*VISC*CAPPA/LOG(ELOG*YPL(IJK))
        VISW(IJK)=MAX(VISC,VISCW)
        VIS(IJK-NIJ)=VISW(IJK)
      END DO
      END DO
!
!.....TOP BOUNDARY (ZERO GRAD OUTLET), SEE CALCSC SUBROUTINE
! BLK23
! 3D
! OUTLET
! BLK23
      DO I=2,NIM
      DO J=2,NJM
        IJ =LI(I)+J
        IJK=LK(NKM)+LI(I)+J
        VIS(IJK+NIJ)=VIS(IJK)
      END DO
      END DO
!
      VISMAX=0.
      DO K=2,NKM
      DO I=2,NIM
      DO J=2,NJM
        IJK=LK(K)+LI(I)+J
        IF (VIS(IJK).GT.VISMAX) THEN
          VISMAX=VIS(IJK)
          IJKMAX=IJK
        END IF
      END DO
      END DO
      END DO
!
      RETURN
      END SUBROUTINE MODVIS
! TURB
!########################################################
      SUBROUTINE MODDAT
!########################################################
!     In this routine some constants are assigned values.
!========================================================
      USE GLOBAL
      IMPLICIT NONE
!--------------------------------------------------------------
      ITE =6
      IED =7
      IVIS=8
!
!.....CONSTANTS OF STANDARD  K-EPSILON MODEL
!
      SIGT  =0.9
      SIGTE =1.0
      SIGED =1.3
      CE1   =1.44
      CE2   =1.92
      CMU   =0.09
      CMU25=SQRT(SQRT(CMU))
      CMU75=CMU25**3
      CTRANS=11.63
      CAPPA =0.41
      ELOG  =8.342
!
      RETURN
      END SUBROUTINE MODDAT
! ORIENT
!###########################################################
      SUBROUTINE BCIN
!###########################################################
!     Setting boundary values for each time step
!-----------------------------------------------------------
      USE GLOBAL
      IMPLICIT NONE
!-----------------------------------------------------------
!
! INLET
!
      FLOMAS=0.
! BLK1
! NANO
      DO I=2,NIM
        DO J=BYS(1),NJM,1
          IJ =(I-1)*NJ+J
          IJK=LK(1)+LI(I)+J
          W(IJK)=ULID
          S=(X(I)-X(I-1))*(Y(J)-Y(J-1))
          FMI(IJ)=DNA2*S*W(IJK)
          FLOMAS=FLOMAS+FMI(IJ)
        END DO
      END DO
!
! ORIENT
!
      END SUBROUTINE BCIN
!
!#############################################################
      SUBROUTINE P3D
!#############################################################
!    This routine plots results in PLOT3D format
!-------------------------------------------------------------
! MPI
      USE GLOBAL
      IMPLICIT NONE
!
      INTEGER :: M,N
      CHARACTER*20 :: FILPLOT
!-------------------------------------------------------------
!
!.....POSTPROCESSING. PLOT3D GRID FILES.
!
      IDIM=NI
      JDIM=NJ
      KDIM=NK
      NVAR=4
!
      ALLOCATE(RXC(1:IDIM,1:JDIM,1:KDIM))
      ALLOCATE(RYC(1:IDIM,1:JDIM,1:KDIM))
      ALLOCATE(RZC(1:IDIM,1:JDIM,1:KDIM))
!
      DO I=1,IDIM,1
         DO J=1,JDIM,1
           DO K=1,KDIM,1
             RXC(I,J,K)=REAL(XC(I))
             RYC(I,J,K)=REAL(YC(J))
             RZC(I,J,K)=REAL(ZC(K))
           END DO
         END DO
       END DO
!
      WRITE(FILPLOT,'(A5,E10.4,A4)') 'pcol_',TIME,'.xyz'
      OPEN (UNIT=JUNIT,FORM='FORMATTED',FILE=FILPLOT)
      REWIND JUNIT
      WRITE(JUNIT,*) IDIM,JDIM,KDIM
      WRITE(JUNIT,*) &
      &(((RXC(I,J,K),I=1,IDIM),J=1,JDIM),K=1,KDIM),&
      &(((RYC(I,J,K),I=1,IDIM),J=1,JDIM),K=1,KDIM),&
      &(((RZC(I,J,K),I=1,IDIM),J=1,JDIM),K=1,KDIM)
      WRITE(JUNIT,*)
      CLOSE(JUNIT)
!
!.....POSTPROCESSING. PLOT3D FUNCTION FILES.
!
      ALLOCATE(RU(1:NIJK),RV(1:NIJK),RW(1:NIJK),RP(1:NIJK))
!
      RU=0.
      RV=0.
      RW=0.
      RP=0.
      DO IJK=1,NIJK
        RU(IJK)=REAL(U(IJK))
        RV(IJK)=REAL(V(IJK))
        RW(IJK)=REAL(W(IJK))
        RP(IJK)=REAL(P(IJK))
      END DO
!
      ALLOCATE(RFLDVAR(-1*IDIM*JDIM+1:IDIM*JDIM*(KDIM+&
      &1),1:NVAR))
!
       DO N=1,NVAR
         DO IJK=1,NIJK
           IF (N.EQ.1) RFLDVAR(IJK,N)=RU(IJK)
           IF (N.EQ.2) RFLDVAR(IJK,N)=RV(IJK)
           IF (N.EQ.3) RFLDVAR(IJK,N)=RW(IJK)
           IF (N.EQ.4) RFLDVAR(IJK,N)=RP(IJK)
         END DO
       END DO
!
       WRITE(FILPLOT,'(A5,E10.4,A2)') 'pcol_',TIME,'.q'
       OPEN (UNIT=KUNIT,FORM='FORMATTED',FILE=FILPLOT)
       REWIND KUNIT
       WRITE(KUNIT,*) IDIM,JDIM,KDIM
       WRITE(KUNIT,*) REAL(ULID/343.2), REAL(0.),&
       &REAL(1000.0),REAL(TIME)
       WRITE(KUNIT,*) (((DENSIT,I=1,IDIM),J=1,JDIM),K=1,KDIM)
       WRITE(KUNIT,*) (((DENSIT*RFLDVAR(LK(K)+LI(I)+J,1),&
       &I=1,IDIM),J=1,JDIM),K=1,KDIM)
       WRITE(KUNIT,*) (((DENSIT*RFLDVAR(LK(K)+LI(I)+J,2),&
       &I=1,IDIM),J=1,JDIM),K=1,KDIM)
       WRITE(KUNIT,*) (((DENSIT*RFLDVAR(LK(K)+LI(I)+J,3),&
       &I=1,IDIM),J=1,JDIM),K=1,KDIM)
       WRITE(KUNIT,*) (((RFLDVAR(LK(K)+LI(I)+J,4),I=1,IDIM),&
       &J=1,JDIM),K=1,KDIM)
       WRITE(KUNIT,*)
       CLOSE(KUNIT)
!
      END SUBROUTINE P3D