! 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
    real*8,dimension(0:MSIZE+1,0:NSIZE) :: V
    integer :: i
    return

contains

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.1e-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 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 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

