AIP(I,J)=AIM(I+1,J)-FLOW
      IF(J.EQ.M2) GOTO 2160
      FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
      FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)*
     + RHO(I-1,J))
      GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+
     + 1.0E-30)*XCVI(I)
      GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)*
     + GAM(I-1,J)+1.E-30)*XCVIP(I-1)
      DIFF=RMN(J+1)*2.*(GM+GMM)
      GOTO 2170
2160  FL=XCVI(I)*V(I,M1)*RHO(I,M1)
      FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1)
      DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1)
2170  FLOW=RMN(J+1)*(FL+FLM)
      CALL DIFLOW
      AJM(I,J+1)=ACOF+AMAX1(0.,FLOW)
      AJP(I,J)=AJM(I,J+1)-FLOW
      VOL=YCVR(J)*XCVS(I)
      APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1))/(XCVS(I)*DT)
      AP(I,J)=AP(I,J)-APT
      CON(I,J)=CON(I,J)+APT*U(I,J)
      AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJM(I,J)+AJP(I,J))/
     + RELAX(NF)
      CON(I,J)=CON(I,J)*VOL+AP(I,J)*U(I,J)*REL
      DU(I,J)=VOL/(XDIF(I)*SX(J))
      CON(I,J)=CON(I,J)+DU(I,J)*(P(I-1,J)-P(I,J))
      DU(I,J)=DU(I,J)/AP(I,J)
2130  CONTINUE
      CALL SOLVE
2110  CONTINUE
*-----------------   COEFFICIENTS FOR THE V EQUATION   -----------------
      NF=2
      IF(.NOT.LSOLVE(NF)) GOTO 2210
      IST=2
      JST=3
      CALL USER(6)
      REL=1.-RELAX(NF)
      DO 2220 I=2,L2
      AREA=R(1)*XCV(I)
      FLOW=AREA*V(I,2)*RHO(I,1)
      DIFF=AREA*GAM(I,1)/YCV(2)
      CALL DIFLOW
2220  AJM(I,3)=ACOF+AMAX1(0.,FLOW)
      DO 2230 J=3,M2
      FL=ARXJ(J)*U(2,J)*RHO(1,J)
      FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1)
      FLOW=FL+FLM
      DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J))
      CALL DIFLOW
      AIM(2,J)=ACOF+AMAX1(0.,FLOW)
      DO 2230 I=2,L2
      IF(I.EQ.L2) GOTO 2240
      FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
      FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)*
     + RHO(I,J-1))
      GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*
     + GAM(I,J)+1.E-30)*ARXJ(J)
      GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)*
     + GAM(I,J-1)+1.E-30)*ARXJP(J-1)
      DIFF=2.*(GM+GMM)/SXMN(J)
      GOTO 2250
2240  FL=ARXJ(J)*U(L1,J)*RHO(L1,J)
      FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1)
      DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J))
2250  FLOW=FL+FLM
      CALL DIFLOW
      AIM(I+1,J)=ACOF+AMAX1(0.,FLOW)
      AIP(I,J)=AIM(I+1,J)-FLOW
      IF(J.EQ.M2) GOTO 2260
      AREA=R(J)*XCV(I)
      FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J)
      FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))*RMN(J+1)
      FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I)
      DIFF=AREA*GAM(I,J)/YCV(J)
      GOTO 2270
2260  AREA=R(M1)*XCV(I)
      FLOW=AREA*V(I,M1)*RHO(I,M1)
      DIFF=AREA*GAM(I,M1)/YCV(M2)
2270  CALL DIFLOW
      AJM(I,J+1)=ACOF+AMAX1(0.,FLOW)
      AJP(I,J)=AJM(I,J+1)-FLOW
      VOL=YCVRS(J)*XCV(I)
      SXT=SX(J)
      IF(J.EQ.M2) SXT=SX(M1)
      SXB=SX(J-1)
      IF(J.EQ.3) SXB=SX(1)
      APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*
     + RHO(I,J-1)*0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT)
      AP(I,J)=AP(I,J)-APT
      CON(I,J)=CON(I,J)+APT*V(I,J)
      AP(I,J)=(-AP(I,J)*VOL+AIM(I,J)+AIP(I,J)+AJM(I,J)+AJP(I,J))/
     + RELAX(NF)
      CON(I,J)=CON(I,J)*VOL+AP(I,J)*V(I,J)*REL
      DV(I,J)=VOL/YDIF(J)
      CON(I,J)=CON(I,J)+DV(I,J)*(P(I,J-1)-P(I,J))
      DV(I,J)=DV(I,J)/AP(I,J)
2230  CONTINUE
      CALL SOLVE
2210  CONTINUE
*---------  COEFFICIENT FOR THE PRESSURE CORRECTION EQUATION  ----------
      NF=3
      IF(.NOT.LSOLVE(NF)) GOTO 2310
      IST=2
      JST=2
      CALL USER(6)
      SMAX=0.
      SSUM=0.
      DO 2320 J=2,M2
      DO 2320 I=2,L2
      VOL=YCVR(J)*XCV(I)
2320  CON(I,J)=CON(I,J)*VOL
      DO 2330 I=2,L2
      ARHO=R(1)*XCV(I)*RHO(I,1)
      CON(I,2)=CON(I,2)+ARHO*V(I,2)
2330  AJM(I,2)=0.
      DO 2340 J=2,M2
      ARHO=ARX(J)*RHO(1,J)
      CON(2,J)=CON(2,J)+ARHO*U(2,J)
      AIM(2,J)=0.
      DO 2340 I=2,L2
      IF(I.EQ.L2) GOTO 2350
      ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
      FLOW=ARHO*U(I+1,J)
      CON(I,J)=CON(I,J)-FLOW
      CON(I+1,J)=CON(I+1,J)+FLOW
      AIP(I,J)=ARHO*DU(I+1,J)
      AIM(I+1,J)=AIP(I,J)
      GOTO 2360
2350  ARHO=ARX(J)*RHO(L1,J)
      CON(I,J)=CON(I,J)-ARHO*U(L1,J)
      AIP(I,J)=0.
2360  IF(J.EQ.M2) GOTO 2370
      ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
      FLOW=ARHO*V(I,J+1)
      CON(I,J)=CON(I,J)-FLOW
      CON(I,J+1)=CON(I,J+1)+FLOW
      AJP(I,J)=ARHO*DV(I,J+1)
      AJM(I,J+1)=AJP(I,J)
      GOTO 2380
2370  ARHO=RMN(M1)*XCV(I)*RHO(I,M1)
      CON(I,J)=CON(I,J)-ARHO*V(I,M1)
      AJP(I,J)=0.
2380  AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)
      PC(I,J)=0.
      SMAX=AMAX1(SMAX,ABS(CON(I,J)))
      SSUM=SSUM+CON(I,J)
2340  CONTINUE
      CALL SOLVE
*---------  COME HERE TO CORRECT THE PRESSURE AND VELOCITIES  ----------
      DO 2390 J=2,M2
      DO 2390 I=2,L2
      P(I,J)=P(I,J)+PC(I,J)*RELAX(NP)
      IF(I.NE.2)U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J))
      IF(J.NE.2)V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J))
2390  CONTINUE
2310  CONTINUE
*-----------------  COEFFICIENTS FOR OTHER EQUATIONS  ------------------
      IST=2
      JST=2
      DO 2410 N=4,NFMAX
      NF=N
      IF(.NOT.LSOLVE(NF)) GOTO 2410
      CALL USER(6)
      REL=1.-RELAX(NF)
      DO 2420 I=2,L2
      AREA=R(1)*XCV(I)
      FLOW=AREA*V(I,2)*RHO(I,1)
      DIFF=AREA*GAM(I,1)/YDIF(2)
      CALL DIFLOW
2420  AJM(I,2)=ACOF+AMAX1(0.,FLOW)
      DO 2430 J=2,M2
      FLOW=ARX(J)*U(2,J)*RHO(1,J)
      DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J))
      CALL DIFLOW
      AIM(2,J)=ACOF+AMAX1(0.,FLOW)
      DO 2430 I=2,L2
      IF(I.EQ.L2) GOTO 2440
      FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
      DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+
     + XCV(I+1)*GAM(I,J)+1.E-30)*SX(J))
      GOTO 2450
2440  FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
      DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J))
2450  CALL DIFLOW
      AIM(I+1,J)=ACOF+AMAX1(0.,FLOW)
      AIP(I,J)=AIM(I+1,J)-FLOW
      AREA=RMN(J+1)*XCV(I)
      IF(J.EQ.M2) GOTO 2460
      FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
      DIFF=AREA*2.*GAM(I,J+1)*GAM(I,J)/(YCV(J)*GAM(I,J+1)+
     + YCV(J+1)*GAM(I,J)+1.E-30)
      GOTO 2470
2460  FLOW=AREA*V(I,M1)*RHO(I,M1)
      DIFF=AREA*GAM(I,M1)/YDIF(M1)
2470  CALL DIFLOW
      AJM(I,J+1)=ACOF+AMAX1(0.,FLOW)
      AJP(I,J)=AJM(I,J+1)-FLOW
      VOL=YCVR(J)*XCV(I)
      APT=RHO(I,J)/DT
      AP(I,J)=AP(I,J)-APT
      CON(I,J)=CON(I,J)+APT*F(I,J,NF)
      AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))/
     + RELAX(NF)
      CON(I,J)=CON(I,J)*VOL+AP(I,J)*F(I,J,NF)*REL
2430  CONTINUE
      CALL SOLVE
2410  CONTINUE
      TIME=TIME+DT
      ITER=ITER+1
      IF(ITER.EQ.LAST) LSTOP=.TRUE.
      RETURN
      END
*=======================================================================
      SUBROUTINE SUPPLY(K)
*-----------------------------------------------------------------------
$INCLUDE:'SIMPLE.INC'
*-----------------------------------------------------------------------
110   FORMAT(1X,26(1H*),3X,A10,3X,26(1H*))
120   FORMAT(1X,'I =',I8,6I10)
130   FORMAT(1X,'J')
140   FORMAT(1X,I2,2X,1P7E10.2)
150   FORMAT(1X,' ')
160   FORMAT(1X,'I =',I8,6I10)
170   FORMAT(1X,'X = ',1P7E10.2)
180   FORMAT(1X,'XU =',1P7E10.2)
190   FORMAT(1X,'TH =',1P7E10.2)
200   FORMAT(1X,'J =',I8,6I10)
210   FORMAT(1X,'Y = ',1P7E10.2)
220   FORMAT(1X,'YV =',1P7E10.2)
*-----------------------------------------------------------------------
      GOTO (1000,2000) K
*-----------------------------------------------------------------------
*     ENTRY UGRID
1000  XU(2)=0.
      DX=XL/FLOAT(L1-2)
      DO 1010 I=3,L1
1010  XU(I)=XU(I-1)+DX
      YV(2)=0.
      DY=YL/FLOAT(M1-2)
      DO 1020 J=3,M1
1020  YV(J)=YV(J-1)+DY
      RETURN
*-----------------------------------------------------------------------
*     ENTRY PRINT
2000  IF(.NOT.LPRINT(3)) GOTO 2010
*------------------   CALCULATE THE STREAM FUNCTION   ------------------
      F(2,2,3)=0.
      DO 2020 I=2,L1
      IF(I.NE.2) F(I,2,3)=F(I-1,2,3)-RHO(I-1,1)*V(I-1,2)*R(1)*XCV(I-1)
      DO 2020 J=3,M1
      RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1)
2020  F(I,J,3)=F(I,J-1,3)+RHOM*U(I,J-1)*ARX(J-1)
2010  CONTINUE
      IF(.NOT.LPRINT(NP)) GOTO 2030
*-----------  CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION  -----------
      DO 2040 J=2,M2
      P(1,J)=(P(2,J)*XCVS(3)-P(3,J)*XDIF(2))/XDIF(3)
2040  P(L1,J)=(P(L2,J)*XCVS(L2)-P(L3,J)*XDIF(L1))/XDIF(L2)
      DO 2050 I=2,L2
      P(I,1)=(P(I,2)*YCVS(3)-P(I,3)*YDIF(2))/YDIF(3)
2050  P(I,M1)=(P(I,M2)*YCVS(M2)-P(I,M3)*YDIF(M1))/YDIF(M2)
      P(1,1)=P(2,1)+P(1,2)-P(2,2)
      P(L1,1)=P(L2,1)+P(L1,2)-P(L2,2)
      P(1,M1)=P(2,M1)+P(1,M2)-P(2,M2)
      P(L1,M1)=P(L2,M1)+P(L1,M2)-P(L2,M2)
      PREF=P(IPREF,JPREF)
      DO 2060 J=1,M1
      DO 2060 I=1,L1
2060  P(I,J)=P(I,J)-PREF
2030  CONTINUE
*-----------------------------------------------------------------------
      WRITE(10,150)
      IEND=1
2080  IF(IEND.EQ.L1) GOTO 2070
      IBEG=IEND+1
      IEND=IEND+7
      IEND=MIN0(IEND,L1)
      WRITE(10,160) (I,I=IBEG,IEND)
      WRITE(10,180) (XU(I),I=IBEG,IEND)
      GOTO 2080
2070  WRITE(10,150)
      JEND=1
2100  IF(JEND.EQ.M1) GOTO 2090
      JBEG=JEND+1
      JEND=JEND+7
      JEND=MIN0(JEND,M1)
      WRITE(10,200) (J,J=JBEG,JEND)
      WRITE(10,220) (YV(J),J=JBEG,JEND)
      GOTO 2100
2090  CONTINUE
*-----------------------------------------------------------------------
      WRITE(10,150)
      IEND=0
2140  IF(IEND.EQ.L1) GOTO 2110
      IBEG=IEND+1
      IEND=IEND+7
      IEND=MIN0(IEND,L1)
      WRITE(10,160) (I,I=IBEG,IEND)
      IF(MODE.EQ.3) GOTO 2120
      WRITE(10,170) (X(I),I=IBEG,IEND)
      GOTO 2130
2120  WRITE(10,190) (X(I),I=IBEG,IEND)
2130  GOTO 2140
2110  WRITE(10,150)
      JEND=0
2160  IF(JEND.EQ.M1) GOTO 2150
      JBEG=JEND+1
      JEND=JEND+7
      JEND=MIN0(JEND,M1)
      WRITE(10,200) (J,J=JBEG,JEND)
      WRITE(10,210) (Y(J),J=JBEG,JEND)
      GOTO 2160
2150  CONTINUE
*-----------------------------------------------------------------------
      DO 2170 N=1,NGAM
      NF=N
      IF(.NOT.LPRINT(NF)) GOTO 2170
      WRITE(10,150)
      WRITE(10,110) TITLE(NF)
      IFST=1
      JFST=1
      IF(NF.EQ.1.OR.NF.EQ.3) IFST=2
      IF(NF.EQ.2.OR.NF.EQ.3) JFST=2
      IBEG=IFST-7
2190  CONTINUE
      IBEG=IBEG+7
      IEND=IBEG+6
      IEND=MIN0(IEND,L1)
      WRITE(10,150)
      WRITE(10,120) (I,I=IBEG,IEND)
      WRITE(10,130)
      JFL=JFST+M1
      DO 2180 JJ=JFST,M1
      J=JFL-JJ
      WRITE(10,140) J,(F(I,J,NF),I=IBEG,IEND)
2180  CONTINUE
      IF(IEND.LT.L1) GOTO 2190
2170  CONTINUE
      RETURN
      END