$DEBUG

*----------------------------MAIN PROGRAM-----------------------------------
****************************************************************************
      LOGICAL LSTOP
      COMMON/CNTL/LSTOP
****************************************************************************
*     CALL SETUP0
      OPEN(08,FILE='RESULT')
      CALL GRID
      CALL SETUP1
      CALL START
   10 CALL DENSE
      CALL BOUND
      CALL OUTPUT
      IF(.NOT.LSTOP) GO TO 15
C*     CLOSE(08)
      STOP
   15 CALL SETUP2
      GO TO 10
      END
*---------------------------------------------------------------------------
      SUBROUTINE DIFLOW
****************************************************************************
      COMMON/COEF/FLOW,DIFF,ACOF
****************************************************************************
      ACOF=DIFF
      IF(FLOW .LT. 1.E-35)RETURN
      TEMP=DIFF-ABS(FLOW)*0.1
      ACOF=0.
      IF(TEMP .LE. 0. ) RETURN
      TEMP=TEMP/DIFF
      ACOF=DIFF*TEMP**5
      RETURN
      END
*--------------------------------------------------------------------------
      SUBROUTINE SOLVE
****************************************************************************
      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(JD),QT(JD)
      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

*****************************************************************************
      ISTF=IST-1
      JSTF=JST-1
      KSTF=KST-1
      IT1=L2+IST
      JT1=M2+JST
      KT1=N2+KST
****************************************************************************
      DO 999 NT=1,NTIMES(NF)
      DO 999 N=NF,NF
*-----------------------------------------------------------------------
      DO 90 K=KST,N2
      DO 90 J=JST,M2
      PT(ISTF)=0.
      QT(ISTF)=F(ISTF,J,K,N)
      DO 70 I=IST,L2
      DENOM=AP(I,J,K)-PT(I-1)*AIM(I,J,K)
      PT(I)=AIP(I,J,K)/DENOM
      TEMP=CON(I,J,K)+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
     &+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
      QT(I)=(TEMP+AIM(I,J,K)*QT(I-1))/DENOM
   70 CONTINUE
      DO 80 II=IST,L2
      I=IT1-II
   80 F(I,J,K,N)=F(I+1,J,K,N)*PT(I)+QT(I)
   90 CONTINUE
*-----------------------------------------------------------------------
      DO 190 KK=KST,N2
      K=KT1-KK
      DO 190 JJ=JST,M2
      J=JT1-JJ
      IF(J.EQ.M2.AND.K.EQ.N2) GOTO 190
      PT(ISTF)=0.
      QT(ISTF)=F(ISTF,J,K,N)
      DO 170 I=IST,L2
      DENOM=AP(I,J,K)-PT(I-1)*AIM(I,J,K)
      PT(I)=AIP(I,J,K)/DENOM
      TEMP=CON(I,J,K)+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
     &+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
      QT(I)=(TEMP+AIM(I,J,K)*QT(I-1))/DENOM
  170 CONTINUE
      DO 180 II=IST,L2
      I=IT1-II
  180 F(I,J,K,N)=F(I+1,J,K,N)*PT(I)+QT(I)
  190 CONTINUE
*-----------------------------------------------------------------------
      DO 290 K=KST,N2
      DO 290 I=IST,L2
      PT(JSTF)=0.
      QT(JSTF)=F(I,JSTF,K,N)
      DO 270 J=JST,M2
      DENOM=AP(I,J,K)-PT(J-1)*AJM(I,J,K)
      PT(J)=AJP(I,J,K)/DENOM
      TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
     &+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
      QT(J)=(TEMP+AJM(I,J,K)*QT(J-1))/DENOM
  270 CONTINUE
      DO 280 JJ=JST,M2
      J=JT1-JJ
  280 F(I,J,K,N)=F(I,J+1,K,N)*PT(J)+QT(J)
  290 CONTINUE
*-----------------------------------------------------------------------
      DO 390 KK=KST,N2
      K=KT1-KK
      DO 390 II=IST,L2
      I=IT1-II
      IF(I.EQ.L2.AND.K.EQ.N2) GOTO 390
      PT(JSTF)=0.
      QT(JSTF)=F(I,JSTF,K,N)
      DO 370 J=JST,M2
      DENOM=AP(I,J,K)-PT(J-1)*AJM(I,J,K)
      PT(J)=AJP(I,J,K)/DENOM
      TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
     &+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
      QT(J)=(TEMP+AJM(I,J,K)*QT(J-1))/DENOM
  370 CONTINUE
      DO 380 JJ=JST,M2
      J=JT1-JJ
  380 F(I,J,K,N)=F(I,J+1,K,N)*PT(J)+QT(J)
  390 CONTINUE
*-----------------------------------------------------------------------
      DO 490 J=JST,M2
      DO 490 I=IST,L2
      PT(KSTF)=0.
      QT(KSTF)=F(I,J,KSTF,N)
      DO 470 K=KST,N2
      DENOM=AP(I,J,K)-PT(K-1)*AKM(I,J,K)

	IF(DENOM.EQ.0.) THEN 
	WRITE(*,998) DENOM,PT(K-1),AKM(I,J,K),AP(I,J,K),NF,I,J,K
      WRITE(2,998) DENOM,PT(K-1),AKM(I,J,K),AP(I,J,K),NF,I,J,K
998	FORMAT(1X,'DENOM=',F15.5,4X,'PT(K-1)=',
     & F15.5,4X,'AKM(I,J,K)=',F15.5,4X,'AP(I,J,K)=',F15.5,4X,'NF=',I3,
	& 'I=',I3,2X,'J=',I3,2X,'K=',I3)

	ENDIF

      PT(K)=AKP(I,J,K)/DENOM

      IF(DENOM.EQ.0.) WRITE(*,*) 'OK'

      TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
     &+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
      QT(K)=(TEMP+AKM(I,J,K)*QT(K-1))/DENOM
  470 CONTINUE
      DO 480 KK=KST,N2
      K=KT1-KK
  480 F(I,J,K,N)=F(I,J,K+1,N)*PT(K)+QT(K)
  490 CONTINUE
*-----------------------------------------------------------------------
      DO 590 JJ=JST,M2
      J=JT1-JJ
      DO 590 II=IST,L2
      I=IT1-II
      IF(I.EQ.L2.AND.J.EQ.M2) GOTO 590
      PT(KSTF)=0.
      QT(KSTF)=F(I,J,KSTF,N)
      DO 570 K=KST,N2
      DENOM=AP(I,J,K)-PT(K-1)*AKM(I,J,K)
      PT(K)=AKP(I,J,K)/DENOM
      TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
     &+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
      QT(K)=(TEMP+AKM(I,J,K)*QT(K-1))/DENOM
  570 CONTINUE
      DO 580 KK=KST,N2
      K=KT1-KK
  580 F(I,J,K,N)=F(I,J,K+1,N)*PT(K)+QT(K)
  590 CONTINUE
************************************************************************
  999 CONTINUE
      DO 600 K=2,N2
      DO 600 J=2,M2
      DO 600 I=2,L2
      CON(I,J,K)=0.
      AP(I,J,K)=0.
  600 CONTINUE
      RETURN
      END
      SUBROUTINE SETUP
************************************************************************
      DOUBLE PRECISION TITLE
      PARAMETER (ID=25,JD=55,KD=25)
      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/CNTL/LSTOP
      COMMON/SORC/SMAX,SSUM
      COMMON/COEF/FLOW,DIFF,ACOF
      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))
************************************************************************
    1 FORMAT(//15X,'COMPUTATION IN CARTISIAN COORDINATES(SS)')
    2 FORMAT(//15X,'COMPUTATION FOR AXISYMMETRICAL SITUATION')
    3 FORMAT(//15X,' COMPUTATION IN POLAR COORDINATES  ')
    4 FORMAT(1X,14X,40(1H*),//)
*
      DATA NFMAX,NP,NRHO,NGAM/10,11,12,13/
      DATA LSTOP,LSOLVE,LPRINT/24*.FALSE./
      DATA LBLK/10*.FALSE./
      DATA MODE,LAST,TIME,ITER/1,5,0.,0/
      DATA RELAX,NTIMES/13*1.,10*1/
      DATA DT,IPREF,JPREF,KPREF,RHOCON/1.E+10,1,1,1,1./
*-----------------------------------------------------------------------
      ENTRY SETUP1
      L2=L1-1
      L3=L2-1
      M2=M1-1
      M3=M2-1
      N2=N1-1
      N3=N2-1
      X(1)=XU(2)
      DO 5 I=2,L2
  5   X(I)=0.5*(XU(I+1)+XU(I))
      X(L1)=XU(L1)
      Y(1)=YV(2)
      DO 10 J=2,M2
   10 Y(J)=0.5*(YV(J+1)+YV(J))
      Y(M1)=YV(M1)
      Z(1)=ZW(2)
      DO 12 K=2,N2
   12 Z(K)=0.5*(ZW(K+1)+ZW(K))
      Z(N1)=ZW(N1)
      DO 15 I=2,L1
   15 XDIF(I)=X(I)-X(I-1)
      DO 18 I=2,L2
   18 XCV(I)=XU(I+1)-XU(I)
      DO 20 I=3,L2
   20 XCVS(I)=XDIF(I)
      XCVS(3)=XCVS(3)+XDIF(2)
      XCVS(L2)=XCVS(L2)+XDIF(L1)
      DO 22 I=3,L3
      XCVI(I)=0.5*XCV(I)
   22 XCVIP(I)=XCVI(I)
      XCVIP(2)=XCV(2)
      XCVI(L2)=XCV(L2)
      DO 35 J=2,M1
   35 YDIF(J)=Y(J)-Y(J-1)
      DO 40 J=2,M2
   40 YCV(J)=YV(J+1)-YV(J)
      DO 45 J=3,M2
   45 YCVS(J)=YDIF(J)
      YCVS(3)=YCVS(3)+YDIF(2)
      YCVS(M2)=YCVS(M2)+YDIF(M1)
      DO 48 J=3,M3
      YCVJ(J)=0.5*YCV(J)
   48 YCVJP(J)=YCVJ(J)
      YCVJP(2)=YCV(2)
      YCVJ(M2)=YCV(M2)
      DO 50 K=2,N1
   50 ZDIF(K)=Z(K)-Z(K-1)
      DO 52 K=2,N2
   52 ZCV(K)=ZW(K+1)-ZW(K)
      DO 55 K=3,N2
   55 ZCVS(K)=ZDIF(K)
      ZCVS(3)=ZCVS(3)+ZDIF(2)
      ZCVS(N2)=ZCVS(N2)+ZDIF(N1)
      DO 58 K=3,N3
      ZCVK(K)=0.5*ZCV(K)
   58 ZCVKP(K)=ZCVK(K)
      ZCVKP(2)=ZCV(2)
      ZCVK(N2)=ZCV(N2)
      DO 60 K=2,N2
      DO 62 J=2,M2
   62 AX(J,K)=YCV(J)*ZCV(K)
      DO 64 I=2,L2
   64 AY(I,K)=XCV(I)*ZCV(K)
   60 CONTINUE
      DO 66 J=2,M2
      DO 66 I=2,L2
   66 AZ(I,J)=XCV(I)*YCV(J)
      DO 85 I=3,L2
      FX(I)=0.5*XCV(I-1)/XDIF(I)
   85 FXM(I)=1.-FX(I)
      FX(2)=0.
      FXM(2)=1.
      FX(L1)=1.
      FXM(L1)=0.
      DO 90 J=3,M2
      FY(J)=0.5*YCV(J-1)/YDIF(J)
   90 FYM(J)=1.-FY(J)
      FY(2)=0.
      FYM(2)=1.
      FY(M1)=1.
      FYM(M1)=0.
      DO 92 K=3,N2
      FZ(K)=0.5*ZCV(K-1)/ZDIF(K)
   92 FZM(K)=1.-FZ(K)
      FZ(2)=0.
      FZM(2)=1.
      FZ(N1)=1.
      FZM(N1)=0.
*---CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE----
      DO 95 K=1,N1
      DO 95 J=1,M1
      DO 95 I=1,L1
      PC(I,J,K)=0.
      U(I,J,K)=0.
      V(I,J,K)=0.
      W(I,J,K)=0.
      CON(I,J,K)=0.
      AP(I,J,K)=0.
      RHO(I,J,K)=RHOCON
      P(I,J,K)=0.