FLOW=ARHO*W(I,J,K+1)
      CON(I,J,K)=CON(I,J,K)-FLOW
      CON(I,J,K+1)=CON(I,J,K+1)+FLOW
      AKP(I,J,K)=ARHO*DW(I,J,K+1)
      AKM(I,J,K+1)=AKP(I,J,K)
      GO TO 409
  408 ARHO=AZ(I,J)*RHO(I,J,N1)
      CON(I,J,K)=CON(I,J,K)-ARHO*W(I,J,N1)
      AKP(I,J,K)=0.
  409 AP(I,J,K)=AIP(I,J,K)+AIM(I,J,K)+AJP(I,J,K)+AJM(I,J,K)+AKP(I,J,K)
     &+AKM(I,J,K)
      PC(I,J,K)=0.
      SMAX=AMAX1(SMAX,ABS(CON(I,J,K)))
      SSUM=SSUM+CON(I,J,K)
  403 CONTINUE
      CALL SOLVE
*---COMEE HERE TO CORRECT THE PRESSURE AND VELOCITIES
      DO 501 K=2,N2
      DO 501 J=2,M2
      DO 501 I=2,L2
      P(I,J,K)=P(I,J,K)+PC(I,J,K)*RELAX(NP)
      IF(I .NE. 2) U(I,J,K)=U(I,J,K)+DU(I,J,K)*(PC(I-1,J,K)-PC(I,J,K))
      IF(J .NE. 2) V(I,J,K)=V(I,J,K)+DV(I,J,K)*(PC(I,J-1,K)-PC(I,J,K))
      IF(K .NE. 2) W(I,J,K)=W(I,J,K)+DW(I,J,K)*(PC(I,J,K-1)-PC(I,J,K))
  501 CONTINUE
  500 CONTINUE
*---COEFFICIENTS FOR OTHER EQUATIONS----
      IST=2
      JST=2
      KST=2
      DO 600 N=5,NFMAX
      NF=N
      IF(.NOT. LSOLVE(NF)) GO TO 600
      CALL GAMSOR
      REL=1.-RELAX(NF)
      DO 599 I=2,L2
      DO 601 K=2,N2
      FLOW=AY(I,K)*V(I,2,K)*RHO(I,1,K)
      DIFF=AY(I,K)*GAM(I,1,K)/YDIF(2)
      CALL DIFLOW
  601 AJM(I,2,K)=ACOF+AMAX1(0.,FLOW)
      DO 602 J=2,M2
      FLOW=AZ(I,J)*RHO(I,J,1)*W(I,J,2)
      DIFF=AZ(I,J)*GAM(I,J,1)/ZDIF(2)
      CALL DIFLOW
  602 AKM(I,J,2)=ACOF+AMAX1(0.,FLOW)
  599 CONTINUE
      DO 603 K=2,N2
      DO 603 J=2,M2
      FLOW=AX(J,K)*U(2,J,K)*RHO(1,J,K)
      DIFF=AX(J,K)*GAM(1,J,K)/XDIF(2)
      CALL DIFLOW
      AIM(2,J,K)=ACOF+AMAX1(0.,FLOW)
      DO 603 I=2,L2
      IF(I .EQ. L2) GO TO 604
      FLOW=AX(J,K)*U(I+1,J,K)*(FX(I+1)*RHO(I+1,J,K)+FXM(I+1)*RHO(I,J,K))
      DIFF=AX(J,K)*2.*GAM(I,J,K)*GAM(I+1,J,K)/(XCV(I)*GAM(I+1,J,K)+
     &   XCV(I+1)*GAM(I,J,K)+1.0E-30)
      GO TO 605
  604 FLOW=AX(J,K)*U(L1,J,K)*RHO(L1,J,K)
      DIFF=AX(J,K)*GAM(L1,J,K)/XDIF(L1)
  605 CALL DIFLOW
      AIM(I+1,J,K)=ACOF+AMAX1(0.,FLOW)
      AIP(I,J,K)=AIM(I+1,J,K)-FLOW
      IF(J .EQ. M2) GO TO 606
      FLOW=AY(I,K)*V(I,J+1,K)*(FY(J+1)*RHO(I,J+1,K)+FYM(J+1)*RHO(I,J,K))
      DIFF=AY(I,K)*2.*GAM(I,J,K)*GAM(I,J+1,K)/(YCV(J)*GAM(I,J+1,K)+
     &   YCV(J+1)*GAM(I,J,K)+1.0E-30)
      GO TO 607
  606 FLOW=AY(I,K)*V(I,M1,K)*RHO(I,M1,K)
      DIFF=AY(I,K)*GAM(I,M1,K)/YDIF(M1)
  607 CALL DIFLOW
      AJM(I,J+1,K)=ACOF+AMAX1(0.,FLOW)
      AJP(I,J,K)=AJM(I,J+1,K)-FLOW
      IF(K .EQ. N2) GO TO 608
      FLOW=AZ(I,J)*W(I,J,K+1)*(FZ(K+1)*RHO(I,J,K+1)+FZM(K+1)*RHO(I,J,K))
      DIFF=AZ(I,J)*2.*GAM(I,J,K)*GAM(I,J,K+1)/(ZCV(K)*GAM(I,J,K+1)+
     &   ZCV(K+1)*GAM(I,J,K)+1.0E-30)
      GO TO 609
  608 FLOW=AZ(I,J)*W(I,J,N1)*RHO(I,J,N1)
      DIFF=AZ(I,J)*GAM(I,J,N1)/ZDIF(N1)
  609 CALL DIFLOW
      AKM(I,J,K+1)=ACOF+AMAX1(0.,FLOW)
      AKP(I,J,K)=AKM(I,J,K+1)-FLOW
      VOL=AX(J,K)*XCV(I)
      APT=RHO(I,J,K)/DT
      AP(I,J,K)=AP(I,J,K)-APT
      CON(I,J,K)=CON(I,J,K)+APT*F(I,J,K,NF)
      AP(I,J,K)=(-AP(I,J,K)*VOL+AIP(I,J,K)+AIM(I,J,K)+AJP(I,J,K)+
     &AJM(I,J,K)+AKP(I,J,K)+AKM(I,J,K))/RELAX(NF)
      CON(I,J,K)=CON(I,J,K)*VOL+REL*AP(I,J,K)*F(I,J,K,NF)
  603 CONTINUE
      CALL SOLVE
  600 CONTINUE
      ITER=ITER+1
      IF(ITER.GE.LAST) LSTOP=.TRUE.
      RETURN
      END
*
      SUBROUTINE SUPPLY
************************************************************************
      PARAMETER (ID=25,JD=55,KD=25)
      DOUBLE PRECISION TITLE
      LOGICAL LSOLVE,LPRINT,LBLK,LSTOP
      COMMON F(ID,JD,KD,10),P(ID,JD,KD),RHO(ID,JD,KD),GAM(ID,JD,KD),
     &CON(ID,JD,KD),AP(ID,JD,KD),AKP(ID,JD,KD),AKM(ID,JD,KD),
     &AIP(ID,JD,KD),AIM(ID,JD,KD),AJP(ID,JD,KD),AJM(ID,JD,KD),
     &X(ID),XU(ID),XDIF(ID),XCV(ID),XCVS(ID),XCVI(ID),XCVIP(ID),
     &Y(JD),YV(JD),YDIF(JD),YCV(JD),YCVS(JD),YCVJ(JD),YCVJP(JD),
     &Z(KD),ZW(KD),ZDIF(KD),ZCV(KD),ZCVS(KD),ZCVK(KD),ZCVKP(KD),
     &AX(JD,KD),AY(ID,KD),AZ(ID,JD)
      COMMON DU(ID,JD,KD),DV(ID,JD,KD),DW(ID,JD,KD),
     & FX(ID),FXM(ID),FY(JD),FYM(JD),FZ(KD),FZM(KD),PT(ID),QT(ID)
      COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,N1,N2,N3,
     & IST,JST,KST,ITER,LAST,TITLE(13),RELAX(13),TIME,DT,XL,YL,ZL,
     & IPREF,JPREF,KPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),
     & RHOCON
	COMMON/DIA/DIAMETER
      COMMON/CNTL/LSTOP
      COMMON/SORC/SMAX,SSUM
      COMMON/COEF/FLOW,DIFF,ACOF
      COMMON/STEP/MSTEP1,MSTEP2,MDSTEP
      DIMENSION U(ID,JD,KD),V(ID,JD,KD),W(ID,JD,KD),PC(ID,JD,KD),
     &T(ID,JD,KD)
      DIMENSION XPLUS(ID,JD,KD),AKE(ID,JD,KD)
     & ,DIS(ID,JD,KD),AMUT(ID,JD,KD),GEN(ID,JD,KD)
      EQUIVALENCE (F(1,1,1,1),U(1,1,1)),(F(1,1,1,2),V(1,1,1)),
     &(F(1,1,1,3),W(1,1,1)),(F(1,1,1,4),PC(1,1,1)),(F(1,1,1,5),T(1,1,1))
     &,(F(1,1,1,6),AKE(1,1,1)),(F(1,1,1,7),DIS(1,1,1))
     & ,(F(1,1,1,8),AMUT(1,1,1)),(F(1,1,1,9),GEN(1,1,1))
      DATA MSTEP1,MSTEP2,MDSTEP/1,1,1/
************************************************************************
   10 FORMAT(1X,26(1H*),3X,A10,3X,26(1H*))
   20 FORMAT(1X,4H I =,I5,19I6)
   21 FORMAT(1X,4H K =,I5,19I6)
   30 FORMAT(1X,1HJ)
   31 FORMAT(1X,1HK)
   40 FORMAT(1X,I2,3X,1P17E7.0)
   41 FORMAT(1X,I2,3X,20F6.3)
   42 FORMAT(1X,I2,3X,20F6.2)
   50 FORMAT(1X,1H )
   51 FORMAT(1X,' I =',2X,7(I4,5X))
   52 FORMAT(1X,' X =',1P7E9.2)
   53 FORMAT(1X,'TH =',1P7E9.2)
   54 FORMAT(1X,'J =',2X,7(I4,5X))
   55 FORMAT(1X,'Y =',1P7E9.2)
   56 FORMAT(1X,'K=',2X,7(I4,5X))
   57 FORMAT(1X,'Z=',1P7E9.2)
   58 FORMAT(1X,26(1H-),3X,'K=',I3,3X,26(1H-))
   59 FORMAT(1X,26(1H-),3X,'J=',I3,3X,26(1H-))
   60 FORMAT(1X,26(1H-),3X,'I=',I3,3X,26(1H-))
************************************************************************
*
      ENTRY UGRID
      XU(2)=0.
      DX=XL/FLOAT(L1-2)
      DO 1 I=3,L1
    1 XU(I)=XU(I-1)+DX
      YV(2)=0.
      DY=YL/FLOAT(M1-2)
      DO 2 J=3,M1
    2 YV(J)=YV(J-1)+DY
      ZW(2)=0.
      DZ=ZL/FLOAT(N1-2)
      DO 3 K=3,N1
    3 ZW(K)=ZW(K-1)+DZ
      RETURN

      ENTRY MYGRID
	RAD=DIAMETER/2.
      XU(2)=0.
	NX=(L1-2)/2
	XU(L1)=XL
	RX=(XL/2-NX*RAD)/(RAD*NX*(NX-1)/2.)
      DO 1001 I=3,(L1+2)/2
	XU(I)=XU(I-1)+RAD+RX*(I-3)*RAD
	IX=L1-I+2
1001	XU(IX)=XU(IX+1)-RAD-RX*(I-3)*RAD

      YV(2)=0.
      DY=YL/FLOAT(M1-2)
      DO 1002 J=3,M1
 1002 YV(J)=YV(J-1)+DY

C      YV(2)=0.
C	NY=(M1-2)/2
C	YV(M1)=YL
C	RY=(YL/2-NY*RAD)/(RAD*NY*(NY-1)/2.)
C      DO 1002 J=3,(M1+2)/2
C	YV(J)=YV(J-1)+RAD+RY*(J-3)*RAD
C	JY=M1-J+2
C1002	YV(JY)=YV(JY+1)-RAD-RY*(J-3)*RAD

      ZW(2)=0.
	NZ=(N1-2)/2
	ZW(N1)=ZL
	RZ=(ZL/2-NZ*RAD)/(RAD*NZ*(NZ-1)/2.)
      DO 1003 K=3,(N1+2)/2
	ZW(K)=ZW(K-1)+RAD+RZ*(K-3)*RAD
	KZ=N1-K+2
1003	ZW(KZ)=ZW(KZ+1)-RAD-RZ*(K-3)*RAD
      RETURN




************************************************************************
       ENTRY PRINT
      WRITE (8,50)
      IEND=0
  301 IF(IEND .EQ. L1) GO TO 310
      IBEG=IEND+1
      IEND=IEND+7
      IEND=MIN0(IEND,L1)
      WRITE (8,50)
      WRITE(8,51) (I,I=IBEG,IEND)
      WRITE(8,52) (XU(I),I=IBEG,IEND)
  303 GO TO 301
  310 JEND=0
      WRITE(8,50)
  311 IF(JEND .EQ. M1) GO TO 320
      JBEG=JEND+1
      JEND=JEND+7
      JEND=MIN0(JEND,M1)
      WRITE(8,50)
      WRITE(8,54) (J,J=JBEG,JEND)
      WRITE(8,55) (YV(J),J=JBEG,JEND)
      GO TO 311
  320 KEND=0
      WRITE(8,50)
  321 IF(KEND.EQ.N1) GOTO 330
      KBEG=KEND+1
      KEND=KEND+7
      KEND=MIN0(KEND,N1)
      WRITE(8,50)
      WRITE(8,56) (K,K=KBEG,KEND)
      WRITE(8,57) (ZW(K),K=KBEG,KEND)
      GO TO 321
  330 CONTINUE

      ENTRY SECTION
      IF(MODE.NE.1) GOTO 119
      DO 999 K=MSTEP1,MSTEP2,MDSTEP
      WRITE(8,58) K
      DO 999 N=1,NGAM
      NF=N
      IF(.NOT. LPRINT(NF)) GO TO 999
      WRITE(8,50)
      WRITE(8,10) TITLE(NF)
      IFST=1
      JFST=1
      IF(NF .EQ. 1 ) IFST=2
      IF(NF .EQ. 2 ) JFST=2
      IBEG=IFST-20
  110 CONTINUE
      IBEG=IBEG+20
      IEND=IBEG+19
      IEND=MIN0(IEND,L1)
      WRITE(8,50)
      WRITE(8,20) (I,I=IBEG,IEND)
      WRITE(8,30)
      JFL=JFST+M1
      DO 115 JJ=JFST,M1
      J=JFL-JJ
      IF(NF.LE.3) WRITE(8,41) J,(F(I,J,K,NF),I=IBEG,IEND)
      IF(NF.GT.3.AND.NF.NE.5) WRITE(8,40)
     &J,(F(I,J,K,NF),I=IBEG,IEND)
      IF(NF.EQ.5) WRITE(8,42) J,(F(I,J,K,NF),I=IBEG,IEND)
  115 CONTINUE
      IF(IEND .LT. L1) GO TO 110
  999 CONTINUE

  119 IF(MODE.NE.2) GOTO 129
      DO 998 J=MSTEP1,MSTEP2,MDSTEP
      WRITE(8,59) J
      DO 998 N=1,NGAM
      NF=N
      IF(.NOT. LPRINT(NF)) GO TO 998
      WRITE(8,50)
      WRITE(8,10) TITLE(NF)
      IFST=1
      KFST=1
      IF(NF .EQ. 1) IFST=2
      IF(NF .EQ. 3) KFST=2
      IBEG=IFST-20
  120 CONTINUE
      IBEG=IBEG+20
      IEND=IBEG+19
      IEND=MIN0(IEND,L1)
      WRITE(8,50)
      WRITE(8,20) (I,I=IBEG,IEND)
      WRITE(8,31)
      KFL=KFST+N1
      DO 160 KK=KFST,N1
      K=KFL-KK
      IF(NF.LE.3) WRITE(8,41) K,(F(I,J,K,NF),I=IBEG,IEND)
      IF(NF.EQ.4) WRITE(8,40) K,(F(I,J,K,NF),I=IBEG,IEND)
      IF(NF.EQ.5) WRITE(8,42) K,(F(I,J,K,NF),I=IBEG,IEND)
  160 CONTINUE
      IF(IEND .LT. L1) GO TO 120
  998 CONTINUE

  129 IF(MODE.NE.3) RETURN
      DO 997 I=MSTEP1,MSTEP2,MDSTEP
      WRITE(8,60) I
      DO 997 N=1,NGAM
      NF=N
      IF(.NOT. LPRINT(NF)) GO TO 997
      WRITE(8,50)
      WRITE(8,10) TITLE(NF)
      KFST=1
      JFST=1
      IF(NF .EQ. 3 ) KFST=2
      IF(NF .EQ. 2 ) JFST=2
      KBEG=KFST-20
  130 CONTINUE
      KBEG=KBEG+20
      KEND=KBEG+19
      KEND=MIN0(KEND,N1)
      WRITE(8,50)
      WRITE(8,21) (K,K=KBEG,KEND)
      WRITE(8,30)
      JFL=JFST+M1
      DO 135 JJ=JFST,M1
      J=JFL-JJ
      IF(NF.LE.3) WRITE(8,41) J,(F(I,J,K,NF),K=KBEG,KEND)
      IF(NF.EQ.4) WRITE(8,40) J,(F(I,J,K,NF),K=KBEG,KEND)
      IF(NF.EQ.5) WRITE(8,42) J,(F(I,J,K,NF),K=KBEG,KEND)
  135 CONTINUE
      IF(KEND .LT. N1) GO TO 130
  997 CONTINUE
      RETURN
      END