C=================================================================C
C  2-Dimensional Cavity Flow by velocity and pressure             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 U(0:NDX,0:NDY), V(0:NDX,0:NDY), P(0:NDX,0:NDY)
      DIMENSION UW(0:NDX,0:NDY), VW(0:NDX,0:NDY)
      DIMENSION B(0:NDX,0:NDY)
      CHARACTER*12 FLOW2
      DATA FLOW2/'FLOW2-??.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
      D5 = DH/2.0
      H  = 1.0/DH
      H2 = H**2
      D2 = DH**2
C  Set Initial Condition
      do j=0,NY
        do i=0,NX
          U(i,j) = 0.0
          V(i,j) = 0.0
          P(i,j) = 0.0
        end do
      end do
C    U=1.0 on L1
      do i=0,NX
        U(i,NY) = 1.0
      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   Compute U,V
        OVER = 0.0
        do j=1,NY-1
          do i=1,NX-1
            UW(i,j) = U(i,j) + DT*( ( U(i-1,j)**2 - U(i+1,j)**2
     1             + U(i,j-1)*V(i,j-1) - U(i,j+1)*V(i,j+1)
     2             + P(i-1,j) - P(i+1,j) )*D5 - ( 4.0*U(i,j) - U(i,j-1)
     3             - U(i-1,j) - U(i+1,j) - U(i,j+1) )*D2/Re )
            VW(i,j) = V(i,j) + DT*( ( V(i,j-1)**2 - V(i,j+1)**2
     1             + U(i-1,j)*V(i-1,j) - U(i+1,j)*V(i+1,j)
     2             + P(i,j-1) - P(i,j+1) )*D5 - ( 4.0*V(i,j) - V(i,j-1)
     3             - V(i-1,j) - V(i+1,j) - V(i,j+1) )*D2/Re )
            OVER   = OVER + abs(U(i,j)) + abs(V(i,j))
          end do
        end do
C   Copy back U,V
        do j=1,NY-1
          do i=1,NX-1
            U(i,j) = UW(i,j)
            V(i,j) = VW(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 P
        CALL PSOR(P,U,V,B,NX,NY,NDX,DT,H,ALP,EPS,ITER)
        NITER = NITER + ITER
C   Output U,V,Phi
        if(MOD(k,MT).eq.0) then
          NI1 = ID/10
          NI2 = ID - NI1*10 
          FLOW2(7:7) = CHAR(NI1+ICHAR('0'))
          FLOW2(8:8) = CHAR(NI2+ICHAR('0'))
          OPEN(UNIT=1,FILE=FLOW2)
          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(P,U,V,NX,NY,NDX)
          CLOSE(UNIT=1)
          if(ID.lt.99) 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 PSOR(P,U,V,B,NX,NY,NDX,DT,H,ALP,EPS,ITER)
C=================================================================C
C  Solve Ax=b by SOR with 2 dimensional FDM                       C
C    Get P with U,V ( Acceleration factor )                       C
C-----------------------------------------------------------------C
C    P(0:NDX,0:NY)   R*8, I/O, P vector                           C
C    U(0:NDX,0:NY)   R*8, In,  U vector                           C
C    V(0:NDX,0:NY)   R*8, In,  V vector                           C
C    B(0:NDX,0:NDY)  R*8, Wk,  Work                               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    DT              R*8, In,  Delta T                            C
C    H               R*8, In,  H=1.0/NX                           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 P(0:NDX,0:NY), U(0:NDX,0:NY), V(0:NDX,0:NY)
      DIMENSION B(0:NDX,0:NY)
C  Set B and BN (B norm) 
      BN = 0.0
      do j=0,NY
        do i=0,NX
C    Boundary on L1
          if(j.ne.NY) then
            U0P = U(i,j+1)
            V0P = V(i,j+1)
          else
            U0P = 1.0
            V0P = 0.0
          end if
C    Boundary on L2
          if(i.ne.0) then
            UM0 = U(i-1,j)
            VM0 = V(i-1,j)
          else
            UM0 = U(1,j)
            VM0 = -V(1,j)
          end if
C    Boundary on L3
          if(j.ne.0) then
            U0M = U(i,j-1)
            V0M = V(i,j-1)
          else
            U0M = -U(i,1) 
            V0M = V(i,1)
          end if
C    Boundary on L4
          if(i.ne.NX) then
            UP0 = U(i+1,j)
            VP0 = V(i+1,j)
          else
            UP0 = U(NX-1,j)
            VP0 = -V(NX-1,j)
          end if
C    Computation
          B(i,j) = ( (UP0 - UM0)**2 + (V0P - V0M)**2 )/4.0
     1           + (VP0 - VM0)*(U0P - U0M)/2.0
     2           - (UP0 - UM0 + V0P - V0M)*H/(2.0*DT)
          BN = BN + B(i,j)**2
        end do
      end do
C   Main Loop
      do k=1,NX*NY
        RN = 0.0
        do j=0,NY
          JP = j + 1
          JM = j - 1
C    Boundary on L1,L3
          if(j.eq.NY) JP = NY - 1
          if(j.eq.0)  JM = 1
          do i=0,NX
            IP = i + 1
            IM = i - 1
C    Boundary on L2,L4
            if(i.eq.NX) IP = NX - 1
            if(i.eq.0)  IM = 1
C    Computation
            R  = (B(i,j)+P(i,JM)+P(IM,j)+P(IP,j)+P(i,JP))/4.0 - P(i,j)
            P(i,j) = P(i,j) + ALP*R
C    Reset on P1
            if(i.eq.0 .and. j.eq.0) then
              P(i,j) = 0.0
              R = 0.0
            end if
            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(P,U,V,NX,NY,NDX)
C=================================================================C
C  Compute U,V and Output P,U,V                                   C
C-----------------------------------------------------------------C
C    P(0:NDX,0:NY)   R*8, In,  P vector                           C
C    U(0:NDX,0:NY)   R*8, Out, U vector                           C
C    V(0:NDX,0:NY)   R*8, Out, V vector                           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 P(0:NDX,0:NY), U(0:NDX,0:NY), V(0:NDX,0:NY)
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,P(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