*=======================================================================
*     FORTRAN CODE FOR SIMPLE METHOD (BY S.V. PATANKAR)
*-----------------------------------------------------------------------
      LOGICAL LSTOP
      COMMON/CNTL/LSTOP
*-----------------------------------------------------------------------
      CALL USER(1)
      CALL SETUP(1)
      CALL USER(2)
100   CALL USER(3)
      CALL USER(4)
      CALL USER(5)
      IF(LSTOP) STOP
      CALL SETUP(2)
      GOTO 100
      END
*=======================================================================
      SUBROUTINE DIFLOW
*-----------------------------------------------------------------------
      COMMON/COEF/FLOW,DIFF,ACOF
*-----------------------------------------------------------------------
      ACOF=DIFF
      IF(FLOW.EQ.0.) 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
*-----------------------------------------------------------------------
$INCLUDE:'SIMPLE.INC'
*-----------------------------------------------------------------------
      ISTF=IST-1
      JSTF=JST-1
      IT1=L2+IST
      IT2=L3+IST
      JT1=M2+JST
      JT2=M3+JST
*-----------------------------------------------------------------------
      DO 100 NT=1,NTIMES(NF)
      DO 100 N=NF,NF
*-----------------------------------------------------------------------
      IF(.NOT.LBLK(NF)) GOTO 110
      PT(ISTF)=0.
      QT(ISTF)=0.
      DO 120 I=IST,L2
      BL=0.
      BLP=0.
      BLM=0.
      BLC=0.
      DO 130 J=JST,M2
      BL=BL+AP(I,J)
      IF(J.NE.M2) BL=BL-AJP(I,J)
      IF(J.NE.JST) BL=BL-AJM(I,J)
      BLP=BLP+AIP(I,J)
      BLM=BLM+AIM(I,J)
      BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)+
     + AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)
130   CONTINUE
      DENOM=BL-PT(I-1)*BLM
      IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E35
      PT(I)=BLP/DENOM
      QT(I)=(BLC+BLM*QT(I-1))/DENOM
120   CONTINUE
      BL=0.
      DO 140 II=IST,L2
      I=IT1-II
      BL=BL*PT(I)+QT(I)
      DO 140 J=JST,M2
140   F(I,J,N)=F(I,J,N)+BL
*-----------------------------------------------------------------------
      PT(JSTF)=0.
      QT(JSTF)=0.
      DO 150 J=JST,M2
      BL=0.
      BLP=0.
      BLM=0.
      BLC=0.
      DO 160 I=IST,L2
      BL=BL+AP(I,J)
      IF(I.NE.L2) BL=BL-AIP(I,J)
      IF(I.NE.IST) BL=BL-AIM(I,J)
      BLP=BLP+AJP(I,J)
      BLM=BLM+AJM(I,J)
      BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)+
     + AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)
160   CONTINUE
      DENOM=BL-PT(J-1)*BLM
      IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E+35
      PT(J)=BLP/DENOM
      QT(J)=(BLC+BLM*QT(J-1))/DENOM
150   CONTINUE
      BL=0.
      DO 170 JJ=JST,M2
      J=JT1-JJ
      BL=BL*PT(J)+QT(J)
      DO 170 I=IST,L2
170   F(I,J,N)=F(I,J,N)+BL
110   CONTINUE
*-----------------------------------------------------------------------
      DO 180 J=JST,M2
      PT(ISTF)=0.
      QT(ISTF)=F(ISTF,J,N)
      DO 190 I=IST,L2
      DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
      PT(I)=AIP(I,J)/DENOM
      TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)
      QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
190   CONTINUE
      DO 200 II=IST,L2
      I=IT1-II
200   F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)
180   CONTINUE
*-----------------------------------------------------------------------
      DO 210 JJ=JST,M3
      J=JT2-JJ
      PT(ISTF)=0.
      QT(ISTF)=F(ISTF,J,N)
      DO 220 I=IST,L2
      DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
      PT(I)=AIP(I,J)/DENOM
      TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)
      QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
220   CONTINUE
      DO 230 II=IST,L2
      I=IT1-II
230   F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)
210   CONTINUE
*-----------------------------------------------------------------------
      DO 240 I=IST,L2
      PT(JSTF)=0.
      QT(JSTF)=F(I,JSTF,N)
      DO 250 J=JST,M2
      DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
      PT(J)=AJP(I,J)/DENOM
      TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
      QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
250   CONTINUE
      DO 260 JJ=JST,M2
      J=JT1-JJ
260   F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)
240   CONTINUE
*-----------------------------------------------------------------------
      DO 270 II=IST,L3
      I=IT2-II
      PT(JSTF)=0.
      QT(JSTF)=F(I,JSTF,N)
      DO 280 J=JST,M2
      DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
      PT(J)=AJP(I,J)/DENOM
      TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
      QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
280   CONTINUE
      DO 290 JJ=JST,M2
      J=JT1-JJ
290   F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)
270   CONTINUE
*-----------------------------------------------------------------------
100   CONTINUE
      DO 300 J=2,M2
      DO 300 I=2,L2
      CON(I,J)=0.
      AP(I,J)=0.
300   CONTINUE
      RETURN
      END
*=======================================================================
      SUBROUTINE SETUP(K)
*-----------------------------------------------------------------------
$INCLUDE:'SIMPLE.INC'
      COMMON/CNTL/LSTOP
      COMMON/SORC/SMAX,SSUM
      COMMON/COEF/FLOW,DIFF,ACOF
*-----------------------------------------------------------------------
10    FORMAT(/15X,' COMPUTATION IN CARTESIAN COORDINATES')
20    FORMAT(/15X,'COMPUTATION FOR AXISYMMETRIC SITUATION')
30    FORMAT(/15X,'   COMPUTATION IN POLAR COORDINATES')
40    FORMAT(14X,40(1H*),/)
*-----------------------------------------------------------------------
      GOTO (1000,2000) K
*-----------------------------------------------------------------------
*     ENTRY SETUP1
1000  NFMAX=10
      NP=11
      NRHO=12
      NGAM=13
      LSTOP=.FALSE.
      DO 1010 I=1,10
      LSOLVE(I)=.FALSE.
      NTIMES(I)=1
1010  LBLK(I)=.TRUE.
      DO 1020 I=1,13
      LPRINT(I)=.FALSE.
1020  RELAX(I)=1.
      LAST=5
      TIME=0.
      ITER=0
      DT=1.E+10
      IPREF=1
      JPREF=1
      RHOCON=1.
*-----------------------------------------------------------------------
      L2=L1-1
      L3=L2-1
      M2=M1-1
      M3=M2-1
      X(1)=XU(2)
      DO 1030 I=2,L2
1030  X(I)=0.5*(XU(I)+XU(I+1))
      X(L1)=XU(L1)
      Y(1)=YV(2)
      DO 1040 J=2,M2
1040  Y(J)=0.5*(YV(J+1)+YV(J))
      Y(M1)=YV(M1)
      DO 1050 I=2,L1
1050  XDIF(I)=X(I)-X(I-1)
      DO 1060 I=2,L2
1060  XCV(I)=XU(I+1)-XU(I)
      DO 1070 I=3,L2
1070  XCVS(I)=XDIF(I)
      XCVS(L2)=XCVS(L2)+XDIF(L1)
      XCVS(3)=XCVS(3)+XDIF(2)
      DO 1080 I=3,L3
      XCVI(I)=0.5*XCV(I)
1080  XCVIP(I)=XCVI(I)
      XCVIP(2)=XCV(2)
      XCVI(L2)=XCV(L2)
      DO 1090 J=2,M1
1090  YDIF(J)=Y(J)-Y(J-1)
      DO 1100 J=2,M2
1100  YCV(J)=YV(J+1)-YV(J)
      DO 1110 J=3,M2
1110  YCVS(J)=YDIF(J)
      YCVS(3)=YCVS(3)+YDIF(2)
      YCVS(M2)=YCVS(M2)+YDIF(M1)
      IF(MODE.NE.1) GOTO 1120
      DO 1130 J=1,M1
      RMN(J)=1.0
1130  R(J)=1.0
      GOTO 1140
1120  DO 1150 J=2,M1
1150  R(J)=R(J-1)+YDIF(J)
      RMN(2)=R(1)
      DO 1160 J=3,M2
1160  RMN(J)=RMN(J-1)+YCV(J-1)
      RMN(M1)=R(M1)
1140  CONTINUE 
      DO 1170 J=1,M1
      SX(J)=1.0
      SXMN(J)=1.0
      IF(MODE.NE.3) GOTO 1170
      SX(J)=R(J)
      IF(J.NE.1) SXMN(J)=RMN(J)
1170  CONTINUE
      DO 1180 J=2,M2
      YCVR(J)=R(J)*YCV(J)
      ARX(J)=YCVR(J)
      IF(MODE.NE.3) GOTO 1180
      ARX(J)=YCV(J)
1180  CONTINUE
      DO 1190 J=4,M3
1190  YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J)
      YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3)
      YCVRS(M2)=.5*(R(M1)+R(M3))*YCVS(M2)
      IF(MODE.NE.2) GOTO 1200
      DO 1210 J=3,M3
      ARXJ(J)=.25*(1.+RMN(J)/R(J))*ARX(J)
1210  ARXJP(J)=ARX(J)-ARXJ(J)
      GOTO 1220
1200  DO 1230 J=3,M3
      ARXJ(J)=.5*ARX(J)
1230  ARXJP(J)=ARXJ(J)
1220  ARXJP(2)=ARX(2)
      ARXJ(M2)=ARX(M2)
      DO 1240 J=3,M3
      FV(J)=ARXJP(J)/ARX(J)
1240  FVP(J)=1.0-FV(J)
      DO 1250 I=3,L2
      FX(I)=.5*XCV(I-1)/XDIF(I)
1250  FXM(I)=1.-FX(I)
      FX(2)=0.0
      FXM(2)=1.
      FX(L1)=1.
      FXM(L1)=0.
      DO 1260 J=3,M2
      FY(J)=.5*YCV(J-1)/YDIF(J)
1260  FYM(J)=1.-FY(J)
      FY(2)=0.
      FYM(2)=1.0
      FY(M1)=1.0
      FYM(M1)=0.
*----------  CON,AP,U,V,RHO,PC,P ARRAYS ARE INITIALIZED HERE  ----------
      DO 1270 J=1,M1
      DO 1270 I=1,L1
      PC(I,J)=0.
      U(I,J)=0.
      V(I,J)=0.
      CON(I,J)=0.
      AP(I,J)=0.
      RHO(I,J)=RHOCON
      P(I,J)=0.
1270  CONTINUE
      IF(MODE.EQ.1) WRITE(10,10)
      IF(MODE.EQ.2) WRITE(10,20)
      IF(MODE.EQ.3) WRITE(10,30)
      WRITE(10,40)
      RETURN
*-----------------------------------------------------------------------
*     ENTRY SETUP2
*-----------------   COEFFICIENTS FOR THE U EQUATION   -----------------
2000  NF=1
      IF(.NOT.LSOLVE(NF)) GOTO 2110
      IST=3
      JST=2
      CALL USER(6)
      REL=1.-RELAX(NF)
      DO 2120 I=3,L2
      FL=XCVI(I)*V(I,2)*RHO(I,1)
      FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1)
      FLOW=R(1)*(FL+FLM)
      DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2)
      CALL DIFLOW
2120  AJM(I,2)=ACOF+AMAX1(0.,FLOW)
      DO 2130 J=2,M2
      FLOW=ARX(J)*U(2,J)*RHO(1,J)
      DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J))
      CALL DIFLOW
      AIM(3,J)=ACOF+AMAX1(0.,FLOW)
      DO 2130 I=3,L2
      IF(I.EQ.L2) GOTO 2140
      FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J))
      FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
      FLOW=ARX(J)*0.5*(FL+FLP)
      DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J))
      GOTO 2150
2140  FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
      DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J))
2150  CALL DIFLOW
      AIM(I+1,J)=ACOF+AMAX1(0.,FLOW)