      PROGRAM flow
      IMPLICIT NONE
      INTEGER           :: i
      INTEGER           :: j
      INTEGER,PARAMETER :: Nx = 70
      INTEGER,PARAMETER :: Ny = 20
      INTEGER,PARAMETER :: ix = 25!Position of the obstacle
      INTEGER,PARAMETER :: iT = 9!x dimension of the obstacle
      INTEGER,PARAMETER :: jW = 5!y dimension of the obstacle
      INTEGER,PARAMETER :: M  = 200!No. of iterations after which the convergence criterion is checked
      INTEGER           :: im,jm!iteration counter
      REAL*16           :: psi(0:Nx+1,0:Ny+1)
      REAL*16           :: zeta(0:Nx+1,0:Ny+1)
      REAL*16           :: psiold(0:Nx+1,0:Ny+1)
      REAL*16           :: R = 0.0001q0!lattice Reynolds number
      REAL*16           :: error = 10.q0!error for the convergence criterion
      REAL*16,PARAMETER :: precisao = 1.q-8!requested precision for the convergence criterion
      REAL*16,PARAMETER :: w = 0.3q0!relaxation parameter
      LOGICAL           :: readfromfile
      OPEN(100,file='flow.dat',status='unknown')
      write(100,*)'Program flow.f95'
      write(100,*)' '
      write(100,*)'Solves the stationary flow past a rectangular obstacle via the iterative method on a dicretized lattice'
      write(100,*)'Lattice size Nx =',Nx
      write(100,*)'Lattice size Ny =',Ny
      write(100,*)'Obstacle width 2W =',2*jW
      write(100,*)'Obstacle depth T =',iT
      write(100,*)'Obstacle distance from the upstream =',ix
      write(100,*)'Lattice Reynolds number = v0 / viscosity =',R
      write(100,*)'Physical Reynolds number = 2W v0 / viscosity =',R/real(2*jW)
      write(100,*)' '
      write(100,*)'No. of iterations after which interactions are checked =',M
      write(100,*)'Relaxation parameter =',w
!!INITIAL CONDITION
      readfromfile = .TRUE.
!      readfromfile = .FALSE.
      if(readfromfile)then
         OPEN(152,file='stream-field.dat',status='unknown')
         OPEN(153,file='vortex-field.dat',status='unknown')
         do i=0,Nx+1
            do j=0,Ny+1
               read(152,*)im,im,psi(i,j)
               read(153,*)im,im,zeta(i,j)
            enddo
         enddo 
         CLOSE(152)
         CLOSE(153)
      else
         do j=1,Ny
            psi(:,j) = real(j,16)!initial condition
         enddo
         zeta(:,:)    = 0.q0!initial condition and boundary condition
      endif
      psi(:,0)     = 0.q0!Boundary AB and DE
      zeta(:,0)    = 0.q0!Boundary AB and DE
      psi(:,Ny+1)  = 2.q0 + psi(:,Ny-1)!Boundary GH
      zeta(:,Ny)   = 0.q0!Boundary GH
      psi(0,:)     = psi(2,:)!Boundary AH
      psi(Nx+1,:)  = psi(Nx-1,:)!Boundary FG
      zeta(Nx+1,:) = zeta(Nx-1,:)!Boundary FG
      do j=0,jW
         psi(ix,j)     = 0.q0 !Boundary BC
         zeta(ix,j)    = -2.q0 * psi(ix-1,j) !Boundary BC
         psi(ix+iT,j)  = 0.q0 !Boundary DE
         zeta(ix+iT,j) = -2.q0 * psi(ix+1,j) !Boundary DE
      enddo
      do i=ix+1,ix+iT-1
         psi(i,jW)  = 0.q0 !Boundary CD
         zeta(i,jW) = -2.q0 * psi(i,jW+1) !Boundary CD
      enddo
!!ITERATE THE METHOD
      R           = R / 16.q0!divide by 16 for latter convenience
      psiold(:,:) = psi(:,:)
      jm          = 0
      do while(error.gt.precisao.and.jm*M.lt.60000)
         do im=1,M!perform the iterations
            CALL updatepsi(Nx,Ny,ix,iT,jW,w,psi,zeta)
            CALL updateBCDE(Nx,Ny,ix,iT,jW,psi,zeta)
            CALL updatezeta(Nx,Ny,ix,iT,jW,w,R,psi,zeta)
         enddo
         CALL checkerror(Nx,Ny,ix,iT,jW,psi,psiold,error)
         jm = jm + 1
      enddo
!!EXPORT DATA
      CALL exportdata(Nx,Ny,ix,iT,jW,psi,zeta)
      write(100,*)' '
      write(100,*)'Precision required for stop =',precisao
      write(100,*)'error =',error
      write(100,*)'No. of iterations required =',M*jm
      write(100,*)'Done.'
      CLOSE(100)
      write(6,*)'Done.'
      STOP
      END PROGRAM flow
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!END OF PROGRAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE updatepsi(Nx,Ny,ix,iT,jW,w,psi,zeta)
      IMPLICIT NONE
      INTEGER,INTENT(IN)    :: Nx
      INTEGER,INTENT(IN)    :: Ny
      INTEGER,INTENT(IN)    :: ix
      INTEGER,INTENT(IN)    :: iT
      INTEGER,INTENT(IN)    :: jW
      REAL*16,INTENT(IN)    :: w
      REAL*16,INTENT(INOUT) :: psi(0:Nx+1,0:Ny+1)
      REAL*16,INTENT(IN)    :: zeta(0:Nx+1,0:Ny+1)
      INTEGER               :: i!x-counter
      INTEGER               :: j!y-counter
      REAL*16               :: aux
!!!!!!UPDATE psi ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!!
      do i=1,ix-1
         do j=1,Ny
            aux = psi(i+1,j) + psi(i-1,j) + psi(i,j+1) + psi(i,j-1) + zeta(i,j)
            aux = aux / 4.q0
            psi(i,j) = psi(i,j)*(1.q0 - w) + w*aux!apply the update
         enddo
      enddo
!!!!!!UPDATE psi ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!
      do i=ix,ix+iT
         do j=jW+1,Ny
            aux = psi(i+1,j) + psi(i-1,j) + psi(i,j+1) + psi(i,j-1) + zeta(i,j)
            aux = aux / 4.q0
            psi(i,j) = psi(i,j)*(1.q0 - w) + w*aux!apply the update
         enddo
      enddo
!!!!!!UPDATE psi ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!!!
      do i=ix+iT+1,Nx
         do j=1,Ny
            aux = psi(i+1,j) + psi(i-1,j) + psi(i,j+1) + psi(i,j-1) + zeta(i,j)
            aux = aux / 4.q0
            psi(i,j) = psi(i,j)*(1.q0 - w) + w*aux!apply the update
         enddo
      enddo
!!!!!!BOUNDARY CONDITIONS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      psi(0,:)     = psi(2,:)!update boundary AH
      psi(:,Ny+1)  = 2.q0 + psi(:,Ny-1)!update boundary GH
      psi(Nx+1,:)  = psi(Nx-1,:)!update boundary FG
      RETURN
      END SUBROUTINE updatepsi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE updatezeta(Nx,Ny,ix,iT,jW,w,R,psi,zeta)
      IMPLICIT NONE
      INTEGER,INTENT(IN)    :: Nx
      INTEGER,INTENT(IN)    :: Ny
      INTEGER,INTENT(IN)    :: ix
      INTEGER,INTENT(IN)    :: iT
      INTEGER,INTENT(IN)    :: jW
      REAL*16,INTENT(IN)    :: w
      REAL*16,INTENT(IN)    :: R
      REAL*16,INTENT(IN)    :: psi(0:Nx+1,0:Ny+1)
      REAL*16,INTENT(INOUT) :: zeta(0:Nx+1,0:Ny+1)
      INTEGER               :: i!x-counter
      INTEGER               :: j!y-counter
      REAL*16               :: aux
!!!!!!UPDATE zeta ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!
      do i=1,ix-1
         do j=1,Ny-1
            aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1))
            aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j))
            aux = R*aux
            aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0
            zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update
         enddo
      enddo
!!!!!!UPDATE zeta ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!
!here, we need to be careful because of the ambiguity of the boundary condition at the corners C and D
      do i=ix+1,ix+iT-1!the corners C and D do not enter in this sum
         do j=jW+1,Ny-1
            aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1))
            aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j))
            aux = R*aux
            aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0
            zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update
         enddo
      enddo
      i=ix
      do j=jW+2,Ny-1!the corners C and D do not enter in this sum
         aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1))
         aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j))
         aux = R*aux
         aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0
         zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update
      enddo
      i=ix+iT
      do j=jW+2,Ny-1!the corners C and D do not enter in this sum
         aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1))
         aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j))
         aux = R*aux
         aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0
         zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update
      enddo
      i = ix
      j = jW+1!the updtate involves corner C
      aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) + 2.q0*psi(i,j) )!here, we used the value of zeta(corner C) = -2 psi(i,j+1)
      aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j))
      aux = R*aux
      aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0
      zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update to the site just besides corner C
      i = ix+iT
      j = jW+1!the updtate involves corner D
      aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) + 2.q0*psi(i,j) )!here, we used the value of zeta(corner D) = -2 psi(i,j+1)
      aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j))
      aux = R*aux
      aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0
      zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update to the site just besides corner D
!!!!!!UPDATE zeta ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!!
      do i=ix+iT+1,Nx
         do j=1,Ny-1
            aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1))
            aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j))
            aux = R*aux
            aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0
            zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update
         enddo
      enddo
!!!!!!BOUNDARY CONDITIONS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      zeta(Nx+1,:)  = zeta(Nx-1,:)!update boundary FG
      RETURN
      END SUBROUTINE updatezeta
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE updateBCDE(Nx,Ny,ix,iT,jW,psi,zeta)
      IMPLICIT NONE
      INTEGER,INTENT(IN)    :: Nx
      INTEGER,INTENT(IN)    :: Ny
      INTEGER,INTENT(IN)    :: ix
      INTEGER,INTENT(IN)    :: iT
      INTEGER,INTENT(IN)    :: jW
      REAL*16,INTENT(IN)    :: psi(0:Nx+1,0:Ny+1)
      REAL*16,INTENT(INOUT) :: zeta(0:Nx+1,0:Ny+1)
      INTEGER               :: i!x-counter
      INTEGER               :: j!y-counter
!!!!!!UPDATE zeta ON BOUNDARIES BC AND DE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      do j=1,jW
         zeta(ix,j)    = -2.q0 * psi(ix-1,j) 
         zeta(ix+iT,j) = -2.q0 * psi(ix+1,j)
      enddo
!!!!!!UPDATE zeta ON BOUNDARY CD!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      do i=ix+1,ix+iT-1
         zeta(i,jW) = -2.q0 * psi(i,jW+1)
      enddo
      RETURN
      END SUBROUTINE updateBCDE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE exportdata(Nx,Ny,ix,iT,jW,psi,zeta)
      IMPLICIT NONE
      INTEGER,INTENT(IN)    :: Nx
      INTEGER,INTENT(IN)    :: Ny
      INTEGER,INTENT(IN)    :: ix
      INTEGER,INTENT(IN)    :: iT
      INTEGER,INTENT(IN)    :: jW
      REAL*16,INTENT(IN)    :: psi(0:Nx+1,0:Ny+1)
      REAL*16,INTENT(IN)    :: zeta(0:Nx+1,0:Ny+1)
      REAL*16               :: vx,vy,dpts,aux
      INTEGER               :: i!x-counter
      INTEGER               :: j!y-counter
      INTEGER               :: k,npts
      npts = 20
      dpts = real(npts,16)
      OPEN(150,file='velocity-field.dat',status='unknown')
      OPEN(151,file='velocity-field2.dat',status='unknown')
!!!!!!EXPORT DATA ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!
      do i=1,ix-1
         do j=1,Ny
            vx = ( psi(i,j+1) - psi(i,j-1) ) / 2.q0
            vy = ( psi(i-1,j) - psi(i+1,j) ) / 2.q0
            write(151,*)i,j,vx,vy
            do k=0,npts
               aux = real(k,16)/dpts
               write(150,*)real(i,16)+vx*aux,real(j,16)+vy*aux
            enddo
            do k=0,npts/2
               aux = real(k,16)/dpts*2.q0
               write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)-0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)+0.2q0*vx*aux
               write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)+0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)-0.2q0*vx*aux
            enddo
         enddo
      enddo
!!!!!!EXPORT DATA ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!
      do i=ix,ix+iT
         do j=jW+1,Ny
            vx = ( psi(i,j+1) - psi(i,j-1) ) / 2.q0
            vy = ( psi(i-1,j) - psi(i+1,j) ) / 2.q0
            write(151,*)i,j,vx,vy
            do k=0,npts
               aux = real(k,16)/dpts
               write(150,*)real(i,16)+vx*aux,real(j,16)+vy*aux
            enddo
            do k=0,npts/2
               aux = real(k,16)/dpts*2.q0
               write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)-0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)+0.2q0*vx*aux
               write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)+0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)-0.2q0*vx*aux
            enddo
         enddo
      enddo
!!!!!!EXPORT DATA ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!!
      do i=ix+iT+1,Nx
         do j=1,Ny
            vx = ( psi(i,j+1) - psi(i,j-1) ) / 2.q0
            vy = ( psi(i-1,j) - psi(i+1,j) ) / 2.q0
            write(151,*)i,j,vx,vy
            do k=0,npts
               aux = real(k,16)/dpts
               write(150,*)real(i,16)+vx*aux,real(j,16)+vy*aux
            enddo
            do k=0,npts/2
               aux = real(k,16)/dpts*2.q0
               write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)-0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)+0.2q0*vx*aux
               write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)+0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)-0.2q0*vx*aux
            enddo
         enddo
      enddo
      CLOSE(150)
      CLOSE(151)
      OPEN(152,file='stream-field.dat',status='unknown')
      OPEN(153,file='vortex-field.dat',status='unknown')
      do i=0,Nx+1
         do j=0,Ny+1
            write(152,*)i,j,psi(i,j)
            write(153,*)i,j,zeta(i,j)
         enddo
      enddo
      CLOSE(152)
      CLOSE(153)
      RETURN
      END SUBROUTINE exportdata
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE checkerror(Nx,Ny,ix,iT,jW,psi,psiold,error)
      IMPLICIT NONE
      INTEGER,INTENT(IN)    :: Nx
      INTEGER,INTENT(IN)    :: Ny
      INTEGER,INTENT(IN)    :: ix
      INTEGER,INTENT(IN)    :: iT
      INTEGER,INTENT(IN)    :: jW
      REAL*16,INTENT(IN)    :: psi(0:Nx+1,0:Ny+1)
      REAL*16,INTENT(INOUT) :: psiold(0:Nx+1,0:Ny+1)
      REAL*16,INTENT(OUT)   :: error
      REAL*16               :: vx,vy,dpts,aux
      INTEGER               :: i!x-counter
      INTEGER               :: j!y-counter
      INTEGER               :: k,npts
      error = 0.q0
!!!!!!EXPORT DATA ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!
      do i=1,ix-1
         do j=1,Ny
            error = error + abs( psi(i,j) - psiold(i,j) )
         enddo
      enddo
!!!!!!EXPORT DATA ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!
      do i=ix,ix+iT
         do j=jW+1,Ny
            error = error + abs( psi(i,j) - psiold(i,j) )
         enddo
      enddo
!!!!!!EXPORT DATA ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!!
      do i=ix+iT+1,Nx
         do j=1,Ny
            error = error + abs( psi(i,j) - psiold(i,j) )
         enddo
      enddo
      error       = error / real(Nx*Ny-it*jw,16)
      psiold(:,:) = psi(:,:)
      RETURN
      END SUBROUTINE checkerror
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
