C     ************************************************
C     **  TWO-DIMENSIONAL POTENTIAL FLOW SOLVED BY  **
C     **                                            **
C     **  FINITE DIFFERENCE METHOD                  **
C     **                                            **
C     **  IN GENERAL COORDINATE                     **
C     **                                            **
C     **  SOR ITERATION METHOD                      **
C     **                                            **
C     **                                            **
C     **  UNIT 2: RESULTS OUTPUT                    **
C     **                                            **
C     ************************************************
C
C     ******************
C     ** MAIN ROUTINE **
C     ******************
C
C
         CALL GRID
C
         CALL INIT
C
         CALL VTRANS
C
         CALL CLEAR
C
         CALL BOUND
C
         CALL POISS
C
         CALL UVP
C
         CALL OUTPUT
C
C
      STOP
      END
C
C     *****************************
C     ** SUBROUTINE INITIAL DATA **
C     *****************************
C
      SUBROUTINE INIT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /CONV1/ EPS1,OMEGA1,ITER1
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
C
C     *** DIVIDE NUMBER ***
C
         NXM1=NX-1
         NYM1=NY-1
         NXP1=NX+1
         NYP1=NY+1
C
C     *** ITERATION PARAMETER ***
C
         ITER1=1000
C
C     *** CONVERGENCE PARAMETER ***
C
         EPS1=1.D-6
C
C     *** RELAXATION PARAMETER ***
C
         OMEGA1=1.8D0
C
C
      RETURN
      END
C
C     **************************
C     ** SUBROUTINE GRID DATA **
C     **************************
C
      SUBROUTINE GRID
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MX=201,MY=101)
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
      COMMON /COORD/ X(MX,MY),Y(MX,MY)
C
C
      READ(1,500) NX,NY
  500 FORMAT(2I10)
C
      READ(1,502)((X(I,J),Y(I,J),I=1,NX),J=1,NY)
  502 FORMAT(2F15.10)
C
C
      RETURN
      END
C
C     ************************************
C     ** INDEPENDENT VARIABLE TRANSLATE **
C     ************************************
C
      SUBROUTINE VTRANS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MX=201,MY=101)
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
      COMMON /COORD/ X(MX,MY),Y(MX,MY)
      COMMON /COEFF/ C1(MX,MY),C2(MX,MY),C3(MX,MY)
     &              ,C4(MX,MY),C5(MX,MY)
      COMMON /TRANS/ XX(MX,MY),XE(MX,MY),YX(MX,MY),YE(MX,MY)
      COMMON /JACOBI/ YAC(MX,MY)
C
      WRITE(*,*)'VTRANS'
C
      DO 100 J=1,NY
      DO 110 I=1,NX
C
C     *****  DX/DXI *****
C
      IF(I.EQ.1) THEN
C
        XX(I,J)=0.5D0*(-X(3,J)+4.D0*X(2,J)-3.D0*X(1,J))
        GOTO 12
C
      END IF
C
      IF(I.EQ.NX) THEN
C
        XX(I,J)=0.5D0*(X(NX-2,J)-4.D0*X(NX-1,J)+3.D0*X(NX,J))
        GOTO 12
C
      END IF
C
        XX(I,J)=0.5D0*(X(I+1,J)-X(I-1,J))
C
C     *****  DX/DETA  *****
C
  12  CONTINUE
C
      IF(J.EQ.1) THEN
C
        XE(I,J)=0.5D0*(-X(I,3)+4.D0*X(I,2)-3.D0*X(I,1))
        GOTO 13
C
      END IF
C
      IF(J.EQ.NY) THEN
C
        XE(I,J)=0.5D0*( X(I,NY-2)-4.D0*X(I,NY-1)+3.D0*X(I,NY))
        GOTO 13
C
      END IF
C
        XE(I,J)=0.5D0*(X(I,J+1)-X(I,J-1))
C
C     *****  DY/DXI  *****
C
   13 CONTINUE
C
      IF(I.EQ.1) THEN
C
        YX(I,J)=0.5D0*(-Y(3,J)+4.D0*Y(2,J)-3.D0*Y(1,J))
        GOTO 22
C
      END IF
C
      IF(I.EQ.NX) THEN
C
        YX(I,J)=0.5D0*(Y(NX-2,J)-4.D0*Y(NX-1,J)+3.D0*Y(NX,J))
        GOTO 22
C
      END IF
C
        YX(I,J)=0.5D0*(Y(I+1,J)-Y(I-1,J))
C
C     *****  DY/DETA  *****
C
   22 CONTINUE
C
      IF(J.EQ.1) THEN
C
        YE(I,J)=0.5D0*(-Y(I,3)+4.D0*Y(I,2)-3.D0*Y(I,1))
        GOTO 23
C
      END IF
C
      IF(J.EQ.NY) THEN
C
        YE(I,J)=0.5D0*(Y(I,NY-2)-4.D0*Y(I,NY-1)+3.D0*Y(I,NY))
        GOTO 23
C
      END IF
C
        YE(I,J)=0.5D0*(Y(I,J+1)-Y(I,J-1))
C
   23 CONTINUE
C
C     *****  JACOBIAN  *****
C
      YAC(I,J)=XX(I,J)*YE(I,J)-XE(I,J)*YX(I,J)
C
C
         IP1=I+1
         IM1=I-1
         JP1=J+1
         JM1=J-1
C
      IF(I.EQ.1  ) IM1=NXM1
      IF(I.EQ.NX) IP1=2
      IF((J.EQ.1).OR.(J.EQ.NY)) GOTO 110
C
C     *** 2ND ORDER DERIVATIVE ***
C
      XXX=X(IP1,J)-2.D0*X(I,J)+X(IM1,J)
      XXE=0.25D0*(X(IP1,JP1)-X(IM1,JP1)-X(IP1,JM1)+X(IM1,JM1))
      XEE=X(I,JP1)-2.D0*X(I,J)+X(I,JM1)
      YXX=Y(IP1,J)-2.D0*Y(I,J)+Y(IM1,J)
      YXE=0.25D0*(Y(IP1,JP1)-Y(IM1,JP1)-Y(IP1,JM1)+Y(IM1,JM1))
      YEE=Y(I,JP1)-2.D0*Y(I,J)+Y(I,JM1)
C
      ALPHA=XE(I,J)*XE(I,J)+YE(I,J)*YE(I,J)
      BETA =XX(I,J)*XE(I,J)+YX(I,J)*YE(I,J)
      GAMMA=XX(I,J)*XX(I,J)+YX(I,J)*YX(I,J)
C
      YACS=1.D0/YAC(I,J)/YAC(I,J)
      YACC=YACS/YAC(I,J)
C
C     *****  DJ/DXI  *****
C
      YACX=XXX*YE(I,J)+YXE*XX(I,J)-XXE*YX(I,J)-YXX*XE(I,J)
C
C     *****  DJ/DETA  *****
C
      YACE=XXE*YE(I,J)+YEE*XX(I,J)-XEE*YX(I,J)-YXE*XE(I,J)
C
C     *****  COEFICIENT OF LAPLACIAN  *****
C
      C1(I,J)=ALPHA*YACS
      C2(I,J)=-2.D0*BETA*YACS
      C3(I,J)=GAMMA*YACS
      C4(I,J)=(XXE*XE(I,J)+YXE*YE(I,J)-XEE*XX(I,J)
     &        -YEE*YX(I,J))*YACS-(ALPHA*YACX-BETA*YACE)*YACC
      C5(I,J)=(XXE*XX(I,J)+YXE*YX(I,J)-XXX*XE(I,J)
     &        -YXX*YE(I,J))*YACS-(GAMMA*YACE-BETA*YACX)*YACC
C
C
  110 CONTINUE
  100 CONTINUE
C
C
      RETURN
      END
C
C     **********************************
C     ** SUBROUTINE ARRAY INITIALIZED **
C     **********************************
C
      SUBROUTINE CLEAR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MX=201,MY=101)
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
      COMMON /RESULT/ PSI(MX,MY)
      COMMON /COORD/ X(MX,MY),Y(MX,MY)
C
C
C     *****  IMPULSIVE START  *****
C
      DO 100 J=1,NY
      DO 100 I=1,NX
C
         PSI(I,J)=0.D0
C
  100 CONTINUE
C
C
      DO 110 J=1,NY
      DO 110 I=1,NX
C
         PSI(I,J)=Y(I,J)
C
  110 CONTINUE
C
C
C
      RETURN
      END
C
C     *************************************
C     ** SUBSTITUTING BOUNDARY CONDITION **
C     *************************************
C
      SUBROUTINE BOUND
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MX=201,MY=101)
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
      COMMON /COORD/ X(MX,MY),Y(MX,MY)
      COMMON /RESULT/ PSI(MX,MY)
C
C
C     **** UNIFORM FLOW CONDITION ****
C
      DO 100 I=1,NX
C
         PSI(I,NY)=Y(I,J)
         PSI(I,1)=0.D0
C
  100 CONTINUE
C
C
      RETURN
      END
C
C     *****************************************
C     ** POISSON EQUATION OF STREAM FUNCTION **
C     *****************************************
C
      SUBROUTINE POISS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MX=201,MY=101)
      COMMON /CONV1/ EPS1,OMEGA1,ITER1
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
      COMMON /COEFF/ C1(MX,MY),C2(MX,MY),C3(MX,MY)
     &              ,C4(MX,MY),C5(MX,MY)
      COMMON /RESULT/ PSI(MX,MY)
      COMMON /JACOBI/ YAC(MX,MY)
C
C
C     ***** SOR ITERATION LOOP *****
C
      DO 500 LOOP=1,ITER1
C
         ERR=0.D0
C
C     *** CAL. INTERNAL POINT ***
C
      DO 110 J=2,NYM1
C
         JP1=J+1
         JM1=J-1
C
C
      DO 120 II=2,NX
C
         I=II
       IF(MOD(LOOP,2).EQ.0) I=NX+2-I
C
         IP1=I+1
         IM1=I-1
C
      IF(I.EQ.NX) IP1=2
C
C
         PO=PSI(I,J)
C
C     ***** FINITE DIFFERENCE EQUATION FOR POISSON EQ. *****
C
      P1=C1(I,J)*(PSI(IP1,J)+PSI(IM1,J))
      P2=C2(I,J)*(PSI(IP1,JP1)-PSI(IM1,JP1)
     &           -PSI(IP1,JM1)+PSI(IM1,JM1))*0.25D0
      P3=C3(I,J)*(PSI(I,JP1)+PSI(I,JM1))
      P4=C4(I,J)*(PSI(IP1,J)-PSI(IM1,J))*0.5D0
      P5=C5(I,J)*(PSI(I,JP1)-PSI(I,JM1))*0.5D0
C
C     **** RESIDUAL *****
C
        PP=P1+P2+P3+P4+P5
        DP=PP/(2.D0*(C1(I,J)+C3(I,J)))-PO
C
        PSI(I,J)=PO+OMEGA1*DP
C
        ERR=ERR+DP*DP
C
C
  120 CONTINUE
C
C     *** BRANCH CUT CONDITION ***
C
         PSI(1,J)=PSI(NX,J)
C
  110 CONTINUE
C
         WRITE(*,*)'LOOP & ERROR=',LOOP,ERR
         IF(ERR.LT.EPS1) GOTO 99
C
  500 CONTINUE
C
   99 CONTINUE
C
C
      WRITE(*,*)
      WRITE(*,610) LOOP
  610 FORMAT('  POISSON ITERATION=',I5)
C
C
      RETURN
      END
C
C     *****************************************
C     ** POISSON EQUATION OF STREAM FUNCTION **
C     *****************************************
C
      SUBROUTINE UVP
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MX=201,MY=101)
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
      COMMON /COEFF/ C1(MX,MY),C2(MX,MY),C3(MX,MY)
     &              ,C4(MX,MY),C5(MX,MY)
      COMMON /TRANS/ XX(MX,MY),XE(MX,MY),YX(MX,MY),YE(MX,MY)
      COMMON /RESULT/ PSI(MX,MY)
      COMMON /JACOBI/ YAC(MX,MY)
      COMMON /PRIMTV/ U(MX,MY),V(MX,MY),P(MX,MY)
      COMMON /CPRESS/ CP(MX)
C
C
      DO 110 J=2,NYM1
C
         JP1=J+1
         JM1=J-1
C
      DO 100 I=2,NX
C
         IP1=I+1
         IM1=I-1
C
      IF(I.EQ.NX) IP1=2
      IF(I.EQ.1  ) IM1=NXM1
C
C
       U(I,J)=-(PSI(IP1,J)-PSI(IM1,J))*.5D0*XE(I,J)/YAC(I,J)
     &        +(PSI(I,JP1)-PSI(I,JM1))*.5D0*XX(I,J)/YAC(I,J)
C
       V(I,J)=-(PSI(IP1,J)-PSI(IM1,J))*.5D0*YE(I,J)/YAC(I,J)
     &        +(PSI(I,JP1)-PSI(I,JM1))*.5D0*YX(I,J)/YAC(I,J)
C
       P(I,J)=3.D0-.5D0*(U(I,J)**2+V(I,J)**2)
C
  100 CONTINUE
C
         U(1,J)=U(NX,J)
         V(1,J)=V(NX,J)
         P(1,J)=P(NX,J)
C
  110 CONTINUE
C
      DO 120 I=2,NX
C
         IP1=I+1
         IM1=I-1
C
      IF(I.EQ.NX) IP1=2
      IF(I.EQ.1  ) IM1=NXM1
C
       U(I,1)=-(PSI(IP1,1)-PSI(IM1,1))*.5D0*XE(I,1)/YAC(I,1)
     &   +(-PSI(I,3)+4.D0*PSI(I,2)-3.D0*PSI(I,1))*.5D0*XX(I,1)/YAC(I,1)
C
       V(I,1)=-(PSI(IP1,1)-PSI(IM1,1))*.5D0*YE(I,1)/YAC(I,1)
     &   +(-PSI(I,3)+4.D0*PSI(I,2)-3.D0*PSI(I,1))*.5D0*YX(I,1)/YAC(I,1)
C
       P(I,1)=3.D0-.5D0*(U(I,1)**2+V(I,1)**2)
C
  120 CONTINUE
C
        U(1,1)=U(1,2)
        V(1,1)=V(1,2)
C       U(1,1)=U(NX,1)
C       V(1,1)=V(NX,1)
C
       P(1,1)=3.D0-.5D0*(U(1,1)**2+V(1,1)**2)
       P(NX,1)=P(1,1)
C      P(1,1)=P(NX,1)
C
         J=NY
C
      DO 130 I=1,NX
C
       U(I,J)=1.D0
       V(I,J)=0.D0
C
       P(I,J)=3.D0-.5D0*(U(I,J)**2+V(I,J)**2)
C
  130 CONTINUE
C
      DO 140 I=1,NX
C
       CP(I)=1.D0-(U(I,1)**2+V(I,1)**2)
C
  140 CONTINUE
C
C
      RETURN
      END
C
C     ****************************
C     ** OUTPUT RESULT TO FILE  **
C     ****************************
C
      SUBROUTINE OUTPUT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(MX=201,MY=101)
      COMMON /DIVIDE/ NXP1,NYP1,NX,NY,NXM1,NYM1
      COMMON /RESULT/ PSI(MX,MY)
      COMMON /PRIMTV/ U(MX,MY),V(MX,MY),P(MX,MY)
      COMMON /CPRESS/ CP(MX)
C
C
      WRITE(2,600) NX,NY
  600 FORMAT(2I10)
C
      WRITE(2,602)((PSI(I,J),I=1,NX),J=1,NY)
  602 FORMAT(F15.7)
C
      WRITE(3,603) NX,NY
  603 FORMAT(2I10)
C
      WRITE(3,604)((U(I,J),V(I,J),P(I,J),I=1,NX),J=1,NY)
  604 FORMAT(3F15.7)
C
      WRITE(4,610)(CP(I),I=1,NX)
  610 FORMAT(F15.7)
C
C
      RETURN
      END