! Algorithm for solving Navier-Stokes Equations from
! Numerical Analysis for Applied Mathematics, Science
! and Engineering, Chapter 9 by Donald Greenspan and
! Vincent Casulli, Addison-Wesley, 1988

program main
    implicit none
    integer,parameter :: MSIZE=10
    integer,parameter :: NSIZE=10
    real*8,dimension(1:MSIZE,1:NSIZE) :: P
    real*8,dimension(0:MSIZE,0:NSIZE+1) :: U,UR
    real*8,dimension(0:MSIZE+1,0:NSIZE) :: V,VR
    real*8 :: x,y,dx,dy,dt,omega,eps,nu
    integer :: i,j
    nu=0.1d0
    dx=1d0/MSIZE
    dy=1d0/NSIZE
    dt=0.01d0
    omega=1d0
    eps=0.001d0
    do i=1,MSIZE-1
        do j=1,NSIZE
            x=i*dx
            y=(j-0.5d0)*dy
            U(i,j)=u0(x,y)
        end do
    end do
    do i=1,MSIZE
        do j=1,NSIZE-1
            x=(i-0.5d0)*dx
            y=j*dy
            V(i,j)=v0(x,y)
        end do
    end do
    do i=1,MSIZE
        do j=1,NSIZE
            P(i,j)=0
        end do
    end do
    call dobndry(U,V)
    call velprint(U)
    print '("------")'
    call dostep
    return

contains

function Fu(i,j)
    implicit none
    integer :: i,j
    real*8 :: Fu,vbar,r
    vbar=(v(i,j)+v(i,j-1)+v(i+1,j-1)+v(i+1,j))/4
    r=(u(i+1,j)-2*u(i,j)+u(i-1,j))/(dx*dx)
    r=r+(u(i,j+1)-2*u(i,j)+u(i,j-1))/(dy*dy)
    r=r*nu
    r=r-u(i,j)*(u(i+1,j)-u(i-1,j))/(2*dx)
    r=r-vbar*(u(i,j+1)-u(i,j-1))/(2*dy)
    r=u(i,j)+dt*r
end function

function Gv(i,j)
    implicit none
    integer :: i,j
    real*8 :: Gv,ubar,r
    ubar=(u(i,j)+u(i,j+1)+u(i-1,j+1)+u(i-1,j))/4
    r=(v(i+1,j)-2*v(i,j)+v(i-1,j))/(dx*dx)
    r=r+(v(i,j+1)-2*v(i,j)+v(i,j-1))/(dy*dy)
    r=r*nu
    r=r-ubar*(v(i+1,j)-v(i-1,j))/(2*dx)
    r=r-v(i,j)*(v(i,j+1)-v(i,j-1))/(2*dy)
    r=v(i,j)+dt*r
end function

function u0(x,y)
    implicit none
    real*8 :: x,y,u0
    u0=0
    return
end function

function v0(x,y)
    implicit none
    real*8 :: x,y,v0
    v0=0
    return
end function

function ugamma(x,y)
    implicit none
    real*8 x,y,ugamma
    if (abs(y-1).le.1d-5) then
        ugamma=1
    else
        ugamma=0
    end if
end function

function vgamma(x,y)
    implicit none
    real*8 :: x,y,vgamma
    vgamma=0
    return
end function

subroutine dobndry(u,v)
    implicit none
    real*8,dimension(0:MSIZE,0:NSIZE+1) :: u
    real*8,dimension(0:MSIZE+1,0:NSIZE) :: v
    integer :: i,j
    real*8 :: x,y
    do i=1,MSIZE
        x=(i-0.5d0)*dx
        v(i,0)=vgamma(x,0d0)
        v(i,NSIZE)=vgamma(x,1d0)
    end do
    do j=1,NSIZE
        y=(j-0.5d0)*dy
        u(0,j)=ugamma(0d0,y)
        u(MSIZE,j)=ugamma(1d0,y)
    end do
    do i=1,MSIZE-1
        x=i*dx
        u(i,NSIZE+1)=2*ugamma(x,1d0)-u(i,NSIZE)
        u(i,0)=2*ugamma(x,0d0)-u(i,1)
    end do
    do j=1,NSIZE-1
        y=j*dy
        v(MSIZE+1,j)=2*vgamma(1d0,y)-u(MSIZE,j)
        v(0,j)=2*vgamma(0d0,y)-u(1,j)
    end do
end subroutine

subroutine dostep
    implicit none
    integer :: i,j,r
    real*8 :: dp
    do i=1,MSIZE-1
        do j=1,NSIZE
            UR(i,j)=Fu(i,j)-dt*(P(i+1,j)-P(i,j))/dx
        end do
    end do
    do i=1,MSIZE
        do j=1,NSIZE-1
            VR(i,j)=Gv(i,j)-dt*(P(i,j+1)-P(i,j))/dy
        end do
    end do
    do r=1,126
        call dobndry(UR,VR)
        do i=1,MSIZE
            do j=1,NSIZE
                dp=(UR(i,j)-UR(i-1,j))/dx
                dp=dp+(VR(i,j)-VR(i,j-1))/dy
                dp=-omega*dp/(2*dt*(1/(dx*dx)+1/(dy*dy)))
                P(i,j)=P(i,j)+dp
                UR(i,j)=UR(i,j)+dt/dx*dp
                UR(i-1,j)=UR(i-1,j)-dt/dx*dp
                VR(i,j)=VR(i,j)+dt/dy*dp
                VR(i,j-1)=VR(i,j-1)-dt/dy*dp
            end do
        end do
    end do
    call velprint(UR)
end subroutine

subroutine sor(w,x,b)
    real*8 :: w
    real*8,dimension(1:) :: x,b
    integer :: i,n
    n=size(x)
    x(1)=x(1)-w/2*(b(1)+2*x(1)-x(2));
    do i=2,n-1
        x(i)=x(i)-w/2*(b(i)-x(i-1)+2*x(i)-x(i+1))
    end do
    x(n)=x(n)-w/2*(b(n)-x(n-1)+2*x(n))
end subroutine

subroutine multA(y,x)
    real*8,dimension(1:) :: x,y
    integer :: i,n
    n=size(x)
    y(1)=-2*x(1)+x(2)
    do i=2,n-1
        y(i)=x(i-1)-2*x(i)+x(i+1)
    end do
    y(n)=x(n-1)-2*x(n)
end subroutine

subroutine velprint(x)
    real*8,dimension(0:,0:) :: x
    integer :: M,N
    M=size(x(:,0))
    N=size(x(0,:))
    print '("M=",I0,"N=",I0)',M,N
    do j=N-1,0,-1
        print '(*(F5.3," "))',(x(i,j),i=0,M-1)
    end do
end subroutine

subroutine dvecprint(x)
    real*8,dimension(1:) :: x
    integer :: N,i
    N=size(x)
    do i=1,N
        print '(G0," ")',x(i)
    end do
end subroutine

end program
