C=================================================================C
C  2-Dimensional Cavity Flow by Stream function and Vorticity     C
C       with FDM and SOR                                          C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         C
C     後 保範（東京工芸大学）                                     C
C=================================================================C
      PARAMETER(NDX=300, NDY=300)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION PHI(0:NDX,0:NDY), OMG(0:NDX,0:NDY)
      DIMENSION U(0:NDX,0:NDY), V(0:NDX,0:NDY), OMGW(0:NDX,0:NDY)
      CHARACTER*11 FLOW1
      DATA FLOW1/'FLOW1-?.txt'/
C  Initial Data
      ID = 0
      CALL XCLOCK(T1,3)
      WRITE(6,*) 'Type In  NX,Re'
      READ(5,*) NX,Re
      WRITE(6,*) 'DT,NT(Total NO.),MT(Output Interval)'
      READ(5,*) DT,NT,MT
      if(NX.gt.NDX) NX = NDX
      NY = NX
      WRITE(6,100) NX,NY,Re,DT,NT,MT
C  Initial Constant
      T  = 0.0
      DH = NX
      H  = 1.0/DH
      H2 = H**2
      D2 = DH**2
C  Set Initial Condition
      do j=0,NY
        do i=0,NX
          OMG(i,j) = 0.0
          PHI(i,j) = 0.0
        end do
      end do
C  SOR Parameter(ALP) 
      ALP = 1.0 + DLOG(NX*1.0D0)/DLOG(NDX*1.2D0)
      EPS = 1.0e-4
C  Main Loop
      CALL XCLOCK(T1,5)
      NITER = 0
      do k=1,NT
        T = T + DT
C   Set Boundary Condition
        do i=0,NX
          OMG(i,0)  = -2.0*PHI(i,1)*D2
          OMG(i,NY) = -2.0*(PHI(i,NY-1) + H)*D2
        end do
        do j=1,NY-1
          OMG(0,j)  = -2.0*PHI(1,j)*D2
          OMG(NX,j) = -2.0*PHI(NX-1,j)*D2
        end do
C   Compute Omega
        OVER = 0.0
        do j=1,NY-1
          do i=1,NX-1
            OMGW(i,j) = OMG(i,j) + DT*D2*( (
     1               - (PHI(i,j+1)-PHI(i,j-1))*(OMG(i+1,j)-OMG(i-1,j)) 
     2               + (PHI(i+1,j)-PHI(i-1,j))*(OMG(i,j+1)-OMG(i,j-1))
     3               ) / 4.0 + (OMG(i,j-1)+OMG(i-1,j)-4.0*OMG(i,j)
     4               + OMG(i+1,j)+OMG(i,j+1)) / Re )
            OVER = OVER + DABS(OMG(i,j))
          end do
        end do
C   Copy back
        do j=1,NY-1
          do i=1,NX-1
            OMG(i,j) = OMGW(i,j)
          end do
        end do
C   Check Divergent
        if(OVER.gt.1.0D10) then
          WRITE(6,*) '** Stop for over flow computation **'
          WRITE(6,*) '  You have to give smaller DT'
          STOP
        end if
C   Compute Phi
        CALL PHISOR(PHI,OMG,NX,NY,NDX,H2,ALP,EPS,ITER)
        NITER = NITER + ITER
C   Output U,V,Phi
        if(MOD(k,MT).eq.0) then
          FLOW1(7:7) = CHAR(ID+ICHAR('0'))
          OPEN(UNIT=1,FILE=FLOW1)
          WRITE(1,100) NX,NY,Re,DT,NT,MT
          WRITE(1,150) T,k,NITER
          WRITE(6,150) T,k,NITER
          NITER = 0
          CALL OUTUVP(PHI,U,V,NX,NY,NDX)
          CLOSE(UNIT=1)
          if(ID.lt.9) ID = ID + 1
        end if
      end do
      CALL XCLOCK(T2,5)
      WRITE(6,200) NX,NT,T2-T1 
C
      STOP
  100 FORMAT('# Cavity Flow Numerical Analysis',/'# NX,NY=',2I5,
     1'  Re=',F6.0,'  DT=',F6.2,'  NT,MT=',2I5)
  150 FORMAT('#Time=',F7.3,'  Step=',I5,'  SOR Loops=',I6) 
  200 FORMAT(1H ,' NX,NT,Time(s)=',2I5,F9.2) 
      END
C=================================================================C
      SUBROUTINE PHISOR(PHI,OMG,NX,NY,NDX,H2,ALP,EPS,ITER)
C=================================================================C
C  Solve Ax=b by SOR with 2 dimensional FDM                       C
C    Given Omega ( Acceleration factor )                          C
C-----------------------------------------------------------------C
C    PHI(0:NDX,0:NY) R*8, I/O, Phi Value                          C
C    OMG(0:NDX,0:NY) R*8, In,  Omega Value                        C
C    NX              I*4, In,  Grid Numbers on X-axis             C
C    NY              I*4, In,  Grid Numbers on Y-axis             C
C    NDX             I*4, In,  First Array Size of PHI,OMG        C
C    H2              R*8, In,  H2=H**2                            C
C    ALP             R*8, In,  SOR Acceleration factor            C
C    EPS             R*8, In,  if ||r||/||b|| <= EPS --> return   C
C    ITER            I*4, Out, Number of Iteration                C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION PHI(0:NDX,0:NY), OMG(0:NDX,0:NY)
C  Get 2-Norm B=D2*OMG
      BN = 0.0
      do j=1,NY-1
        do i=1,NX-1
          BN = BN + (H2*OMG(i,j))**2
        end do
      end do
C   Main Loop
      do k=1,NX*NY
        RN = 0.0
        do j=1,NY-1
          do i=1,NX-1
            R = (H2*OMG(i,j) + PHI(i,j-1) + PHI(i-1,j) + PHI(i+1,j)
     1          + PHI(i,j+1) )/4.0 - PHI(i,j)
            PHI(i,j) = PHI(i,j) + ALP*R
            RN       = RN + (R*4.0)**2
          end do
        end do
C    if(ERR <= EPS) return
        ERR = DSQRT(RN/BN)
        if(ERR.le.EPS) goto 100
      end do
  100 CONTINUE
      ITER = k
C
      RETURN
      END
C=================================================================C
      SUBROUTINE OUTUVP(PHI,U,V,NX,NY,NDX)
C=================================================================C
C  Compute U,V and Output U,V,Phi                                 C
C    U=d(PHI)/dy, V=-d(PHI)/dx                                    C
C-----------------------------------------------------------------C
C    PHI(0:NDX,0:NY) R*8, In,  Phi Value                          C
C    U(0:NDX,0:NY)   R*8, Out, U Value                            C
C    V(0:NDX,0:NY)   R*8, Out, V Value                            C
C    NX              I*4, In,  Grid Numbers on X-axis             C
C    NY              I*4, In,  Grid Numbers on Y-axis             C
C    NDX             I*4, In,  First Array Size of PHI,OMG        C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION PHI(0:NDX,0:NY), U(0:NDX,0:NY), V(0:NDX,0:NY)
C   Initial 
      XV = NX/2.0D0
      YV = NY/2.0D0
C   Compute U,V
      do j=1,NY-1
        do i=1,NX-1
          U(i,j) = (PHI(i,j+1) - PHI(i,j-1))*YV
          V(i,j) = (PHI(i-1,j) - PHI(i+1,j))*XV
        end do
      end do
C    Boundary on Y=0,1
      do i=0,NX
        U(i,NY) = 1.0
        V(i,NY) = 0.0
        U(i,0)  = 0.0
        V(i,0)  = 0.0
      end do
C    Boundary on X=0,1
      do j=1,NY-1
        U(0,j)  = 0.0
        V(0,j)  = 0.0
        U(NX,j) = 0.0
        V(NX,j) = 0.0
      end do
C   Output
      DX = 1.0D0/NX
      DY = 1.0D0/NY
      do j=0,NY
        Y = j*DY
        do i=0,NX
          X = i*DX
          WRITE(1,100) X,Y,PHI(i,j),U(i,j),V(i,j)
        end do
      end do
C
      RETURN
  100 FORMAT(1H ,5F9.5)
      END
C=================================================================C
      SUBROUTINE XCLOCK(CPU,ID)
C=================================================================C
C   CPU time Subroutine                                           C
C     CPU     R*8  Out, CPU Time                                  C
C     ID      I*4  In,  Dummy ( Same Hitachi FORTRAN XCLOCK )     C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         C
C=================================================================C
      REAL*8 CPU
      INTEGER*2 I1, I2, I3, I4
C
      IF(ID.GE.1) THEN
        CALL GETTIM(I1,I2,I3,I4)
        CPU = ( I1*60.0 + I2 )*60.0 + I3 + I4*0.01
      END IF
C
      RETURN
      END