c23456789012345678901234567890
C     Program to compute entrance flow
C     Vort(Io+1,JJ) given by Stokes flow
      DIMENSION U(202,42),V(202,42),VORT(202,42),PSI(202,42)
      DIMENSION VORTN1(202,42),YCORD(42)
      DIMENSION CP(2,102),SUM(102),DYVOR(102),AVEDY(102)
      DIMENSION CPP(2,102)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    1 FORMAT(5F12.9)
C 
C      SET INPUT DATA 
C
      N=5
      II=N*20
      IO=N*10
      JJ=N*10
      AL=4.
      RE=1.
      EPSI=1.E-04
      EVORT=1.E-04
C 
C      COMPUTE DX , BETA , AND F
C 
      AII=FLOAT(II) 
      AJJ=FLOAT(JJ)
      DX=  AL/AII
      BETA=AL*FLOAT(JJ)/AII
      BETA2=BETA**2 
      PI=ACOS(-1.)
      BE2P1=BETA2+1.
      DT=.5/((1.+BETA)/DX+BE2P1*4./RE/DX**2)  
      E=(COS(PI/FLOAT(2*II+1))+BETA2*COS(PI/FLOAT(JJ)))/BE2P1       
      ETA=E**2
      F=2.*(1-SQRT(1.-ETA))/ETA       
C 
C       SET INITIAL CONDITIONS AND BOUNDARY VALUES
C 
      DO 20 I=1,II+1
      DO 10 J=1,JJ+1
      U(I,J)=1.   
      V(I,J)=0. 
      VORT(I,J)=0.
   10 PSI(I,J)=(FLOAT(J)-1.)/AJJ
   20 CONTINUE
      DO 25 I=IO+1,II+1      
      U(I,JJ+1)=0.
      V(I,JJ+1)=0.
      PSI(I,JJ+1)=1.
   25 VORT(I,JJ+1)=3.       
      DO 27 J=1,JJ+1
      Y=(FLOAT(J)-1.)/AJJ       
      V(II+1,J)=0. 
      U(II+1,J)=1.5*(1.-Y**2)
      PSI(II+1,J)=1.5*Y-.5*Y**3      
   27 VORT(II+1,J)=3.*Y    
C 
C      SOLVE FOR VORTICITY AT INTERIOR POINTS 
C 
   30 ISW2=0
      DO 38 I=2,II 
      DO 37 J=2,JJ
      DELSQ=VORT(I+1,J)+VORT(I-1,J)+BETA2*VORT(I,J+1)+BETA2*VORT(I,J- 
     11)-2.*BE2P1*VORT(I,J) 
      IF(U(I,J).GT.0.) GO TO 31       
      CONU=U(I+1,J)*VORT(I+1,J)-U(I,J)*VORT(I,J)
      GO TO 32
   31 CONU=U(I,J)*VORT(I,J)-U(I-1,J)*VORT(I-1,J)
   32 IF(V(I,J).GT.0) GO TO 35
      CONV=V(I,J+1)*VORT(I,J+1)-V(I,J)*VORT(I,J)
      GO TO 36    
   35 CONV=V(I,J)*VORT(I,J)-V(I,J-1)*VORT(I,J-1)
   36 DVORT=DT/DX*(-(CONU+BETA*CONV)+2./RE/DX*DELSQ)      
      VORTN1(I,J)=VORT(I,J)+DVORT     
   37 IF(DVORT.GT.EVORT) ISW2=1       
   38 CONTINUE    
C 
C       UPDATE VORTICITY MATRIX 
C 
      DO 40 I=2,II 
      DO 39 J=2,JJ
   39 VORT(I,J)=VORTN1(I,J) 
   40 CONTINUE
C 
C       SOLVE FOR STREAM FUNCTION 
C 
   41 ISW1=0
      DO 50 I=2,II 
      DO 49 J=2,JJ
      DSTR=PSI(I+1,J)+PSI(I-1,J)+BETA2*PSI(I,J+1)+BETA2*PSI(I,J-1)-2.*
     1BE2P1*PSI(I,J)+VORT(I,J)*DX**2  
      PSI(I,J)=PSI(I,J)+F/2./BE2P1*DSTR 
   49 IF(DSTR.GT.EPSI) ISW1=1 
   50 CONTINUE
      IF(ISW1.GT.0) GO TO 41
C 
C      CALCULATE U AND V VELOCITIES AT INTERIOR POINTS
C 
      DO 60 I=2,II 
      DO 59,J=2,JJ
      U(I,J)=(PSI(I,J+1)-PSI(I,J-1))*BETA/2./DX 
   59 V(I,J)=(PSI(I-1,J)-PSI(I+1,J))/2./DX    
   60 CONTINUE
C 
C      CALCULATE CENTERLINE & STAG STREAMLINE VALUES OF U 
C 
      DO 75 I=2,II 
   75 U(I,1)=PSI(I,2)*BETA/DX 
      DO 80 I=2,IO
   80 U(I,JJ+1)=(1.-PSI(I,JJ))*BETA/DX
C 
C      CALCULATE VORT ON THE WALLS
C 
      DO 90 I=IO+2,II
      VORTEMP=VORT(I,JJ+1)  
      VORT(I,JJ+1)=(1.-PSI(I,JJ))*2*BETA2/DX**2 
      DVORT=VORT(I,JJ+1)-VORTEMP      
   90 IF(DVORT.GT.EVORT) ISW2=1       
      VORT(IO+1,JJ+1)=VORT(IO+1,JJ)+(1-PSI(IO+1,JJ))/DX**2*5*BETA2/4
      IF(ISW2.GT.0) GO TO 30
C 
C      PRINT OUTPUT 
C 
      DO 95 I=1,II+1,2*N
      XCORD=(FLOAT(I)-FLOAT(IO)-1.)*AL/AII
      MN=40+N
      WRITE(MN,*) XCORD
      DO 95 J=1,JJ+1,2*N
      YCORD(J)=(FLOAT(J)-1.)/AJJ
      WRITE(MN,1) YCORD(J),PSI(I,J),VORT(I,J),U(I,J),V(I,J)
   95 CONTINUE
      END