program main
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real::	dx, dy, dt
 real::	xnue, density, length, height, time, inlet_velocity, outlet_pressure
 real,dimension(0:md,0:nd)::	u, v, p, u_old, v_old 
 real,dimension(0:md)::		xp
 real,dimension(0:nd)::		yp
 real,dimension(0:md,0:nd)::	ap, ae, aw, an, as, bb
 integer::	m, n, istep, istep_max
 integer::	i,j

!-----------------

! ----------------
! set up condistions

call  physical_conditions (xnue, density, length, height, time, inlet_velocity, outlet_pressure)
call  grid_conditions (xp, yp, dx, dy, dt, xnue, density, length, height, time, inlet_velocity, m, n, istep_max)
call  output_grid (xp, yp, m, n)

istep	=0
time	=istep* dt

call  initial_conditions (p, u, v, xp, yp, length, height, inlet_velocity, outlet_pressure, m, n)
call  boundary(p, u, v, xp, yp, length, height, inlet_velocity, outlet_pressure, m, n)

! print initial conditions
! call  output_solution (p, u, v, m, n)

! ----------------
! MAC algorithm start

do istep = 1, istep_max

time	=istep* dt
write(*,*)'--- time_steps= ',istep, ' --  time = ',time

do i = 0, m+1
do j = 0, n+1
 u_old(i,j)	=u(i,j)
 v_old(i,j)	=v(i,j)
end do
end do

call  solve_p (p, u, v, u_old, v_old, xnue, density, height, yp, dx, dy, dt, m, n)
call  solve_u (p, u, v, u_old, v_old, xnue, density, dx, dy, dt, m, n)
call  solve_v (p, u, v, u_old, v_old, xnue, density, dx, dy, dt, m, n)

call  boundary(p, u, v, xp, yp, length, height, inlet_velocity, outlet_pressure, m, n)

!call output_solution (p, u, v, m, n)

end do

call  output_divergent (p, u, v, dx, dy, m, n)

! MAC algorithm end
! ----------------

! print conditions (recall)
! call  physical_conditions (xnue, density, length, height, time, inlet_velocity, outlet_pressure)
! call  grid_conditions (xp, yp, dx, dy, dt, xnue, density, length, height, time, inlet_velocity, m, n, istep_max)

! print solutions 
! call  output_solution (p, u, v, m, n)
call  output_solution_post (p, u, v, m, n)
call  output_paraview (p, xp, yp, m, n, u, v)

stop
end program main
!******************

!  solve variables  

!******************
subroutine  solve_p (p, u, v, u_old, v_old, xnue, density, height, yp, dx, dy, dt, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in)::	dx, dy, dt
 real,intent(in)::	xnue, density, height
 real,intent(inout),dimension(0:md,0:nd)::	u, v, p, u_old, v_old 
 real,intent(in),dimension(0:nd)::		yp
 integer,intent(in)::	m, n

!-----------------
! local variables 
 real::		u_stg, v_stg
 real,dimension(0:md,0:nd)::	ap, ae, aw, an, as, bb
 integer::	i, j

! ---------------
interface
subroutine  solve_matrix (p, ap, ae, aw, an, as, bb, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(inout),dimension(0:md,0:nd)::	p
 real,intent(in),dimension(0:md,0:nd)::		ap, ae, aw, an, as, bb
 integer,intent(in)::	m, n
end subroutine  solve_matrix

subroutine  boundrary_matrix (p, ap, ae, aw, an, as, bb, m, n, height, yp)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in)::	height
 real,intent(in),dimension(0:md,0:nd)::	p
 real,intent(inout),dimension(0:md,0:nd)::	ap, ae, aw, an, as, bb
 real,intent(in),dimension(0:nd)::		yp
 integer,intent(in)::	m, n
end subroutine  boundrary_matrix
end interface

! ----------------
do i = 1, m
do j = 1, n
! ----------------
!   velocity u
! ----------------
! convection_x  (1st upwind scheme)
u(i,j)=u_old(i,j)						&
      -dt*max(u_old(i,j),0.)*(u_old(i,j)-u_old(i-1,j))/dx	&  ! u>0 
      -dt*min(u_old(i,j),0.)*(u_old(i+1,j)-u_old(i,j))/dx	   ! u<0

! convection_y
u(i,j)=u(i,j)							&
      -dt*max(v_old(i,j),0.)*(u_old(i,j)-u_old(i,j-1))/dy	&  ! v>0 
      -dt*min(v_old(i,j),0.)*(u_old(i,j+1)-u_old(i,j))/dy	   ! v<0

! diffusion_x
u(i,j)=u(i,j)							&
      +dt*xnue*(u_old(i+1,j)-2*u_old(i,j)+u_old(i-1,j))/dx/dx	

! diffusion_y
u(i,j)=u(i,j)							&
      +dt*xnue*(u_old(i,j+1)-2*u_old(i,j)+u_old(i,j-1))/dy/dy 	

! ----------------
!   velocity v
! ----------------
! convection_x  (1st upwind scheme)
v(i,j)=v_old(i,j)						&
      -dt*max(u_old(i,j),0.)*(v_old(i,j)-v_old(i-1,j))/dx	&    ! u>0 
      -dt*min(u_old(i,j),0.)*(v_old(i+1,j)-v_old(i,j))/dx		 ! u<0

! convection_y
v(i,j)=v(i,j)							&
      -dt*max(v_old(i,j),0.)*(v_old(i,j)-v_old(i,j-1))/dy	&    ! v>0 
      -dt*min(v_old(i,j),0.)*(v_old(i,j+1)-v_old(i,j))/dy	     ! v<0

! diffusion_x
v(i,j)=v(i,j)							&
      +dt*xnue*(v_old(i+1,j)-2*v_old(i,j)+v_old(i-1,j))/dx/dx	

! diffusion_y
v(i,j)=v(i,j)							&
      +dt*xnue*(v_old(i,j+1)-2*v_old(i,j)+v_old(i,j-1))/dy/dy	
end do
end do

! ----------------
! matrix solution 

do i = 1, m
do j = 1, n
ae(i,j)= dt/density/dx/dx
aw(i,j)= dt/density/dx/dx
an(i,j)= dt/density/dy/dy
as(i,j)= dt/density/dy/dy
ap(i,j)= -ae(i,j)-aw(i,j)-an(i,j)-as(i,j)

bb(i,j)= (u(i+1,j)-u(i-1,j))/dx/2+(v(i,j+1)-v(i,j-1))/dy/2
end do
end do

call boundrary_matrix (p, ap, ae, aw, an, as, bb, m, n, height, yp)

call solve_matrix (p, ap, ae, aw, an, as, bb, m, n)

! ----------------
return
end subroutine solve_p
!******************

!******************
subroutine  solve_matrix (p, ap, ae, aw, an, as, bb, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(inout),dimension(0:md,0:nd)::	p
 real,intent(in),dimension(0:md,0:nd)::		ap, ae, aw, an, as, bb
 integer,intent(in)::	m, n

! local variables
real::		relux_factor, error
real,dimension(0:md, 0:nd)::	p_old
integer::	i, j, iter, iter_max

! ----------------
!   SOR algorithm
! ----------------
iter_max = 300		! SOR max interation steps
relux_factor=1.8	! SOR reluxation factor

do iter = 1, iter_max

error=0.

do i = 1, m
do j = 1, n
 p_old(i,j) = p(i,j)
end do
end do

do i = 1, m
do j = 1, n
p(i,j) = (  bb(i,j)					&
          - ae(i,j)*p_old(i+1,j) -aw(i,j)*p(i-1,j)	&
          - an(i,j)*p_old(i,j+1) -as(i,j)*p(i,j-1) )	&
         /ap(i,j)    * relux_factor			&
        + p_old(i,j) * (1.-relux_factor)

 error = max(error, abs( p(i,j) - p_old(i,j) ) )
end do
end do

!write(*,*)'iteration no.', iter,'  -- error=', error

end do

write(*,*)'SOR iteration no.', iter,'  -- error=', error

! ----------------

return
end subroutine solve_matrix
!******************

!******************
subroutine  boundrary_matrix (p, ap, ae, aw, an, as, bb, m, n, height, yp)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in)::	height
 real,intent(in),dimension(0:md,0:nd)::		p
 real,intent(inout),dimension(0:md,0:nd)::	ap, ae, aw, an, as, bb
 real,intent(in),dimension(0:nd)::		yp
 integer,intent(in)::	m, n

! local variables
integer	i, j

! ----------------
! inlet (dp/x=0 at i=1)
do j= 1, n
 ae(1,j) =ae(1,j)+aw(1,j)
 aw(1,j) =0.
end do

! outlet (p=outlet_pressure at i=m)
do j= 1, n
 bb(m,j) =ap(m,j)*p(m,j)
 ae(m,j) = 0.
 aw(m,j) = 0.
 an(m,j) = 0.
 as(m,j) = 0.
end do

! lower wall (dp/dy=0. at j=1)
do i= 1, m
 an(i,1) =an(i,1)+as(i,1)
 as(i,1) = 0.
end do

! symmetry  (dp/dy=0. at j=n)
do i= 1, m
 as(i,n) =as(i,n)+an(i,n)
 an(i,n) = 0.
end do
! ----------------

return
end subroutine  boundrary_matrix
!******************

!******************
subroutine  solve_u (p, u, v, u_old, v_old, xnue, density, dx, dy, dt, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in)::	dx, dy, dt
 real,intent(in)::	xnue, density
 real,intent(inout),dimension(0:md,0:nd)::	u, v, p, u_old, v_old 
 integer,intent(in)::	m, n

! local variables
integer::	i, j

! ----------------
do i = 1, m
do j = 1, n
! convection_x  (1st upwind scheme)
! (already calculated in solve_p)

! convection_y
! (already calculated in solve_p)

! diffusion_x
! (already calculated in solve_p)

! diffusion_y
! (already calculated in solve_p)

! pressure
u(i,j)=u(i,j) -dt/density*(p(i+1,j)-p(i-1,j))/dx/2

end do
end do

! ----------------
return
end subroutine solve_u
!******************

!******************
subroutine  solve_v (p, u, v, u_old, v_old, xnue, density, dx, dy, dt, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in)::	dx, dy, dt
 real,intent(in)::	xnue, density
 real,intent(inout),dimension(0:md,0:nd)::	u, v, p, u_old, v_old 
 integer,intent(in)::	m, n

! local variables
integer::	i, j

! ----------------
do i = 1, m
do j = 1, n
! convection_x  (1st upwind scheme)
! (already calculated in solve_p)

! convection_y
! (already calculated in solve_p)

! diffusion_x
! (already calculated in solve_p)

! diffusion_y
! (already calculated in solve_p)

! pressure
v(i,j)=v(i,j) -dt/density*(p(i,j+1)-p(i,j-1))/dy/2

end do
end do
! ----------------
return
end subroutine solve_v
!******************

!  conditions  

!******************
subroutine  boundary(p, u, v, xp, yp, length, height, inlet_velocity, outlet_pressure, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in)::	length, height, inlet_velocity, outlet_pressure
 real,intent(inout),dimension(0:md,0:nd)::	u, v, p
 real,intent(in),dimension(0:md)::		xp
 real,intent(in),dimension(0:nd)::		yp
 integer,intent(in)::	m, n

! local variables
integer	i, j

! ----------------
! inlet (u=inlet_velocity, v=0., dp/dx=0 at i=1)
do j= 1, n
 u(1,j) =inlet_velocity
 v(1,j) =0.
 u(0,j) =inlet_velocity 	! dummy
 v(0,j) =0.					! dummy
 p(0,j) =p(2,j)
end do

! outlet (du/dx=0., v=0., p=outlet_pressure at i=m)
do j= 1, n
 u(m+1,j) =u(m-1,j)
 v(m+1,j)=0.
 v(m,j)=0.					! dummy
 p(m,j) =outlet_pressure
 p(m+1,j)=outlet_pressure
end do

! lower wall (u=0., v=0., dp/dy=0. at j=1)
do i= 0, m+1
 u(i,1) =0.
 v(i,1) =0.
 u(i,0) =0.					! dummy
 v(i,0) =0.					! dummy
 p(i,0) =p(i,2)
end do

! symmetry  (du/dy=0., v(i,n)=0., dp/dy=0. at j=n)
do i= 0, m+1
 u(i,n+1) = u(i,n-1)  
 v(i,n) =0.
 v(i,n+1) =0.      			! dummy  
 p(i,n+1) =p(i,n-1)
end do
! ----------------

return
end subroutine boundary
!*****************************

!*****************************
subroutine physical_conditions (xnue, density, length, height, time, inlet_velocity, outlet_pressure)
 implicit none
 real,intent(inout)::	xnue, density, length, height, time, inlet_velocity, outlet_pressure

! local variables
 real::	reynolds_no

! ----------------
xnue		=1.e-4  ! (m2/s) :water
density	=1000.  ! (kg/m3):water

length		=1.0     ! (m)
height		=0.1     ! (m)
time		=100.	  ! (s)
inlet_velocity=0.1     ! (m/s)
outlet_pressure=0.     ! (N/m2)

reynolds_no=height*inlet_velocity/xnue

write(*,*) 
write(*,*) 'xnue ='	, xnue
write(*,*) 'density ='	, density
write(*,*) 'length ='	, length
write(*,*) 'height ='	, height
write(*,*) 'time ='	, time
write(*,*) 'inlet_velocity ='	, inlet_velocity
write(*,*) 'outlet_pressure ='	, outlet_pressure
write(*,*) 'reynolds_no='	, reynolds_no
write(*,*) 
! ----------------

return
end subroutine physical_conditions
!******************

!******************
subroutine  grid_conditions (xp, yp, dx, dy, dt, xnue, density, length, height, time, inlet_velocity, m, n, istep_max)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(inout)::	dx, dy, dt
 real,intent(in)::	xnue, density, length, height, time, inlet_velocity
 real,intent(inout),dimension(0:md)::		xp
 real,intent(inout),dimension(0:nd)::		yp
 integer,intent(inout)::	m, n, istep_max

! local variables
real::	cfl_no, pecret_no
integer::	i, j

! ----------------
m	=10
n	=10
istep_max=1000
dx	=length/real(m)
dy	=height/real(n)
dt	=time/real(istep_max)

cfl_no		= inlet_velocity*dt/dx
pecret_no	= inlet_velocity*dx/xnue

write(*,*)
write(*,*) 'm, n =', m, n
write(*,*) 'dx, dy =', dx, dy
write(*,*) 'dt =', dt
write(*,*) 'cfl_no =', cfl_no
write(*,*) 'pecret_no =', pecret_no
write(*,*) 

xp(1)	=0.
do i = 2, m
 xp(i)	=xp(i-1) + dx
end do
xp(0)	=xp(1) - dx			! dummy
xp(m+1)	=xp(m) + dx			! dummy

yp(1)	=0.
do j = 2, n
 yp(j)	=yp(j-1) + dy
end do
yp(0)	=yp(1) - dy			! dummy
yp(n+1)	=yp(n) + dy			! dummy
! ----------------

return
end subroutine  grid_conditions
!******************

!******************
subroutine  initial_conditions (p, u, v, xp, yp, length, height, inlet_velocity, outlet_pressure, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in)::	length, height, inlet_velocity, outlet_pressure
 real,intent(out),dimension(0:md,0:nd)::	u, v, p 
 real,intent(in),dimension(0:md)::		xp
 real,intent(in),dimension(0:nd)::		yp
 integer,intent(in)::	m, n

! local variables
integer::	i, j

! ----------------
do j = 1, n
do i = 1, m
 u(i,j)=inlet_velocity
 v(i,j)=0.
 p(i,j)=outlet_pressure
end do
end do
! ----------------

return
end subroutine initial_conditions
!******************

! output

!******************
subroutine  output_solution (p, u, v, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in),dimension(0:md,0:nd)::	u, v, p 
 integer,intent(in)::	m, n

! local variables
integer::	i, j

! ----------------
write(*,*)

write(*,*)'velocity u '
do j = 0, n+1
write(*,*) (u(i,j), i=0,m+1)
end do
write(*,*)

write(*,*)'velocity v '
do j = 0, n+1
write(*,*) (v(i,j), i=0,m+1)
end do
write(*,*)

write(*,*)'pressure'
do j = 0, n+1
write(*,*) (p(i,j), i=0,m+1)
end do
write(*,*)
! ----------------

return
end subroutine output_solution
!******************

!******************
subroutine  output_grid (xp, yp, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in),dimension(0:md)::	xp
 real,intent(in),dimension(0:nd)::	yp
 integer,intent(in)::	m, n

! local variables
integer::	i, j

open (60, file='grid.dat', status='replace')

! ----------------
write(60,*)'m, n =', m, n
write(60,*)'grid points ='
write(60,*) (xp(i), i=1,m)
write(60,*) (yp(j), j=1,n)
! ----------------

close (60)

return
end subroutine output_grid
!******************

!******************
subroutine  output_solution_post (p, u, v, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in),dimension(0:md,0:nd)::	u, v, p
 integer,intent(in)::	m, n

! local variables
integer::	i, j
real,dimension(0:md, 0:nd)::	u_cnt, v_cnt

open (61, file='solution_uvp.dat', status='replace')

! ----------------
! interpolation at p-center grid

do i = 1, m
do j = 1, n
 u_cnt(i,j)=u(i,j)
 v_cnt(i,j)=v(i,j)
end do
end do

do j = 1, n
 u_cnt(0,j)=u_cnt(1,j)
 v_cnt(0,j)=v_cnt(1,j)
 u_cnt(m+1,j)=u_cnt(m,j)
 v_cnt(m+1,j)=v_cnt(m,j)
end do

do i = 0, m+1
 u_cnt(i,0)=u_cnt(i,1)
 v_cnt(i,0)=v_cnt(i,1)
 u_cnt(i,n+1)=u_cnt(i,n)
 v_cnt(i,n+1)=v_cnt(i,n)
end do

!-----------------
write(61,*)'m, n =', m, n

write(61,*)'velocity u '
do j = 1, n
write(61,*) (u_cnt(i,j), i=1,m)
end do

write(61,*)'velocity v '
do j = 1, n
write(61,*) (v_cnt(i,j), i=1,m)
end do

write(61,*)'pressure'
do j = 1, n
write(61,*) (p(i,j), i=1,m)
end do

! ----------------
close (61)

return
end subroutine output_solution_post
!******************

!******************
subroutine  output_paraview (p, xp, yp, m, n, u, v)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
! real,intent(in),dimension(0:md,0:nd)::	p
 real,intent(in),dimension(0:md)::	xp
 real,intent(in),dimension(0:nd)::	yp
 real,dimension(0:md, 0:nd)::	u, v, p
 integer,intent(in)::	m, n
 integer::	i, j

open(50,file='output_paraview.vtk',status="unknown",form="formatted",position="rewind")
!open(*,file='solution.vtk',status="replace")

    write(50,"('# vtk DataFile Version 3.0')")
    write(50,"('2D flow')")
    write(50,"('ASCII ')")
    
    write(50,"('DATASET STRUCTURED_GRID')")
    write(50,"('DIMENSIONS ',3(1x,i3))") m, n, 1
    
    write(50,"('POINTS ',i9,' float')") m*n
    do j=1,n
    do i=1,m
      write(50,"(3(f16.4,1x))") xp(i), yp(j), 0.0d0
    enddo
    enddo
    
    write(50,"('POINT_DATA ',i9)") m*n
    
!! velocity vector
    write(50,"('VECTORS velocity float')")
    do j=1,n
    do i=1,m
      write(50,"(3(f16.4,1x))") u(i,j), v(i,j), 0.0d0
    enddo
    enddo
    
!! pressure
    write(50,"('SCALARS pressure float')")
    write(50,"('LOOKUP_TABLE default')")
    do j=1,n
    do i=1,m
!      write(50,"(f16.4)") p(i,j)
      write(50,"(3(f16.4,1x))") p(i,j)
    enddo
    enddo
close(50)

return
end subroutine  output_paraview
!******************

!******************
subroutine  output_divergent (p, u, v, dx, dy, m, n)
 implicit none
 integer,parameter::	md=101, nd=101     ! md, nd > grid size (m,n)
 real,intent(in),dimension(0:md,0:nd)::	u, v, p 
 real,intent(in)::	dx, dy
 integer,intent(in)::	m, n

! local variables
integer::	i, j
real,dimension(0:md,0:nd)::	div

open (62, file='divergent.dat', status='replace')

do i = 1, m
do j = 1, n
div(i,j)= (u(i+1,j)-u(i-1,j))/dx/2+(v(i,j+1)-v(i,j-1))/dy/2
end do
end do

write(62,*)
write(62,*)'divergent velocity'
do j = 1, n
write(62,*) (div(i,j), i=1,m)
end do
write(62,*)

close(62)

end subroutine  output_divergent
!******************
