From: Felipe Abrão on
Please, could somebody explain me what is happening. I simply
allocated a matrix called "a", but figured out that the program is
inserting some strange values to some of its elements (like
1.057229231618489E-307 or 6.941570408926931E+228), even before I set
any value to it! Shouldn't all the values be zero before any desired
value setting? Below is my code, but you don't have to read it all.
You can check the problem just after line 65, where I wrote a command
in order to read the values of matrix "a". The same thing is happening
with matrixs "b" and "c". In advance I thank you all!

PS: I use Fortran 90 in Windows XP.


! EXERCICIO.f90
!
! FUNCTIONS:
! EXERCICIO - Entry point of console application.
!

!
****************************************************************************
!
! PROGRAM: EXERCICIO
!
! PURPOSE: Entry point for the console application.
!
!
****************************************************************************

program EXERCICIO

use IMSL

implicit none

! Variables

integer(4)::i,j,k,l,m,Nx,Ny
integer(4),parameter::ipath=1
integer(4),dimension(6)::iparam

integer(4),allocatable,dimension(:)::irow_u,jcol_u,irow_v,jcol_v,irow_p,jcol_p
real(8)::dx,dy,D_w,D_e,D_s,D_n,D_EE,D_F,D_G,D_H,D_I,D_J,D_K,D_L
real(8),parameter::Re=1.0d2,gama=1.0d0/
Re,u_E=1.0d0,u_G=0.0d0,u_H=0.0d0,v_I=0.0d0,v_J=0.0d0,&

v_K=0.0d0,v_L=0.0d0,error=1.0d-6,alpha_u=5.0d-1,alpha_v=5.0d-1,alpha_p=5.0d-1
!real(8),dimension(N*N)::theta_old
real(8),dimension(5)::rparam

real(8),allocatable,dimension(:)::u,v,p,theta,S,S_v,S_p,X,Y,u_F,AA,BB,CC
!real(8),dimension(N,N)::theta_fake

real(8),allocatable,dimension(:,:)::a,a_inv,b,c,p_fake,p_old,u_fake,u_old,v_fake,v_old,a_P,&

a_W,a_E,a_S,a_N,b_P,b_W,b_E,b_S,b_N,c_P,c_W,c_E,c_S,c_N,F_w_u,F_e_u,F_s_u,F_n_u,S_fake,&
S_v_fake,S_p_fake,F_w_v,F_e_v,F_s_v,F_n_v

! Body of EXERCICIO_1

iparam(1)=0

Nx=5
Ny=5




allocate(a(Ny*Nx,Ny*Nx),a_inv(Ny*Nx,Ny*Nx),b(Ny*Nx,Ny*Nx),c(Ny*Nx,Ny*Nx))

allocate(a_P(Ny,Nx),a_W(Ny,Nx),a_E(Ny,Nx),a_S(Ny,Nx),a_N(Ny,Nx),b_P(Ny,Nx),b_W(Ny,Nx),&

b_E(Ny,Nx),b_S(Ny,Nx),b_N(Ny,Nx),c_P(Ny,Nx),c_W(Ny,Nx),c_E(Ny,Nx),c_S(Ny,Nx),c_N(Ny,Nx),&

p_fake(Ny,Nx),p_old(Ny,Nx),u_fake(Ny,Nx),u_old(Ny,Nx),v_fake(Ny,Nx),v_old(Ny,Nx),F_w_u(Ny,Nx),&

F_e_u(Ny,Nx),F_s_u(Ny,Nx),F_n_u(Ny,Nx),S_fake(Ny,Nx),S_v_fake(Ny,Nx),S_p_fake(Ny,Nx),&
F_w_v(Ny,Nx),F_e_v(Ny,Nx),F_s_v(Ny,Nx),F_n_v(Ny,Nx))

allocate(u(Ny*Nx),v(Ny*Nx),p(Ny*Nx),theta(Ny*Nx),S(Ny*Nx),S_v(Ny*Nx),S_p(Ny*Nx))
allocate(Y(Ny),u_F(Ny))
allocate(X(Nx))


open(7,file='pre.dat')

do i=1,Ny*Nx
do j=1,Ny*Nx
write(7,*) i,j,a(i,j)
end do
end do



dx=1.0d0/dble(Nx)
dy=1.0d0/dble(Ny)
D_w=gama/dx
D_e=gama/dx
D_s=gama/dy
D_n=gama/dy
D_EE=2.0d0*D_w
D_F=2.0d0*D_e
D_G=2.0d0*D_s
D_H=2.0d0*D_n
D_I=D_EE
D_J=D_F
D_K=D_G
D_L=D_H

p_fake=5.0d-1
u_fake=5.0d-1
v_fake=5.0d-1

! do i=1,Ny
! do j=1,Nx
! u_fake(i,j)=dble(i+j)
! v_fake(i,j)=dble(i+j)
! p_fake(i,j)=dble(i+j)
! end do
! end do



theta=5.0d-1
!theta_old=6.0d-1

do i=1,Nx
X(i)=1.0d0/dble((Nx-1))*dble((i-1))
end do

do i=1,Ny
Y(i)=1.0d0/dble((Ny-1))*dble((i-1))
u_F(i)=4.0d0*(Y(i)-Y(i)**2.0d0)
end do

p_old=1.0d0
u_old=1.0d0
v_old=1.0d0

! do while(maxval(dabs(p_fake-p_old)).gt.error)

!do while((maxval(dabs(p_fake-p_old)).gt.error).or.(maxval(dabs(u_fake-
u_old)).gt.error).or.&
!(maxval(dabs(v_fake-v_old)).gt.error))

!do while((maxval(dabs(p_fake(2:Ny-1,2:Nx-1)-
p_old(2:Ny-1,2:Nx-1))).gt.error).or.(maxval(dabs(&
!u_fake(2:Ny-1,2:Nx-1)-u_old(2:Ny-1,2:Nx-1))).gt.error).or.
(maxval(dabs(v_fake(2:Ny-1,2:Nx-1)-&
!v_old(2:Ny-1,2:Nx-1))).gt.error))



! write(*,*) 'diferença de pressão =',maxval(dabs(p_fake-p_old))
! write(*,*) 'diferença de velocidade u =',maxval(dabs(u_fake-u_old))
! write(*,*) 'diferença de velocidade v =',maxval(dabs(v_fake-v_old))


! p_old=p_fake
! u_old=u_fake
! v_old=v_fake


do i=1,Ny
do j=1,Nx
if (j.eq.1) then
F_w_u(2:Ny-1,j)=1.0d0
else
F_w_u(i,j)=1.0d0/2.0d0*(u_fake(i,j)+u_fake(i,j-1))
end if
if (i.eq.1) then
F_w_v(i,j)=0.0d0 !Tenho dúvida.
else
F_w_v(i,j)=1.0d0/2.0d0*(u_fake(i,j)+u_fake(i-1,j))
end if
if (j.eq.Nx) then
F_e_u(i,j)=4.0d0*(Y(i)-Y(i)**2.0d0)
else
F_e_u(i,j)=1.0d0/2.0d0*(u_fake(i,j+1)+u_fake(i,j))
end if
if (i.eq.1) then
F_e_v(i,j)=0.0d0 !Tenho dúvida.
else if (j.eq.Nx) then
F_e_v(i,j)=4.0d0*(Y(i)-Y(i)**2.0d0) !Tenho dúvida.
else
F_e_v(i,j)=1.0d0/2.0d0*(u_fake(i,j+1)+u_fake(i-1,j+1))
end if
if (i.eq.1) then
F_s_u(i,j)=0.0d0
else if (j.eq.1) then
F_s_u(i,j)=0.0d0
else
F_s_u(i,j)=1.0d0/2.0d0*(v_fake(i,j)+v_fake(i,j-1))
end if
if (i.eq.1) then
F_s_v(i,j)=0.0d0
else
F_s_v(i,j)=1.0d0/2.0d0*(v_fake(i-1,j)+v_fake(i,j))
end if
if (i.eq.Ny) then
F_n_u(i,j)=0.0d0
else if (j.eq.1) then
F_n_u(i,j)=0.0d0
else
F_n_u(i,j)=1.0d0/2.0d0*(v_fake(i+1,j)+v_fake(i+1,j-1))
end if
if (i.eq.Ny) then
F_n_v(i,j)=0.0d0
else
F_n_v(i,j)=1.0d0/2.0d0*(v_fake(i,j)+v_fake(i+1,j))
end if
end do
end do







!interior

do i=3,Ny-1
do j=3,Nx-1
!coeficientes de u
a_P(i,j)=(F_e_u(i,j)*dy+F_n_u(i,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_u
a_W(i,j)=F_w_u(i,j)*dy+F_e_u(i,j)*1.0d0/8.0d0*dy+D_w*dy
a_E(i,j)=D_e*dy
a_S(i,j)=F_s_u(i,j)*dx+F_n_u(i,j)*1.0d0/8.0d0*dx+D_s*dx
a_N(i,j)=D_n*dx

S_fake(i,j)=(3.0d0/8.0d0*u_fake(i,j)*dy-2.0d0/8.0d0*u_fake(i,j-1)*dy-1.0d0/8.0d0&
*u_fake(i,j-2)*dy)*F_w_u(i,j)+(-3.0d0/8.0d0*u_fake(i,j+1)*dy
+2.0d0/8.0d0*&
u_fake(i,j)*dy)*F_e_u(i,j)
+(3.0d0/8.0d0*u_fake(i,j)*dx-2.0d0/8.0d0*u_fake(i-1,j)&
*dx-1.0d0/8.0d0*u_fake(i-2,j)*dx)*F_s_u(i,j)+(-3.0d0/8.0d0*u_fake(i
+1,j)*dx+&
2.0d0/8.0d0*u_fake(i,j)*dx)*F_n_u(i,j)+(p_fake(i,j-1)-
p_fake(i,j))*dy+((1.0d0-&
alpha_u)*a_P(i,j))*u_old(i,j)

!coeficientes de v
b_P(i,j)=(F_e_v(i,j)*dy+F_n_v(i,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_v
b_W(i,j)=F_w_v(i,j)*dy+F_e_v(i,j)*1.0d0/8.0d0*dy+D_w*dy
b_E(i,j)=D_e*dy
b_S(i,j)=F_s_v(i,j)*dx+F_n_v(i,j)*1.0d0/8.0d0*dx+D_s*dx
b_N(i,j)=D_n*dx

S_v_fake(i,j)=(3.0d0/8.0d0*v_fake(i,j)*dy-2.0d0/8.0d0*v_fake(i,j-1)*dy-1.0d0/
&
8.0d0*v_fake(i,j-2)*dy)*F_w_v(i,j)+(-3.0d0/8.0d0*v_fake(i,j+1)*dy
+2.0d0/8.0d0*&
v_fake(i,j)*dy)*F_e_v(i,j)
+(3.0d0/8.0d0*v_fake(i,j)*dx-2.0d0/8.0d0*v_fake(i-1,j)&
*dx-1.0d0/8.0d0*v_fake(i-2,j)*dx)*F_s_v(i,j)+(-3.0d0/8.0d0*v_fake(i
+1,j)*dx+&
2.0d0/8.0d0*v_fake(i,j)*dx)*F_n_v(i,j)+(p_fake(i-1,j)-
p_fake(i,j))*dx+((1.0d0-&
alpha_v)*b_P(i,j))*v_old(i,j)

!coeficientes de p
c_W(i,j)=dy/a_P(i,j)*dy
c_E(i,j)=dy/a_E(i,j)*alpha_u*dy
c_S(i,j)=dx/b_P(i,j)*dx
c_N(i,j)=dx/b_N(i,j)*alpha_v*dx
c_P(i,j)=c_W(i,j)+c_E(i,j)+c_S(i,j)+c_N(i,j)
S_p_fake(i,j)=u_fake(i,j)*dy-u_fake(i,j+1)*dy+v_fake(i,j)*dx-
v_fake(i+1,j)*dx

end do
end do

!coluna 2 - com exceção da linha 2

do i=3,Ny-1
!coeficientes de u
a_P(i,2)=(F_e_u(i,2)*dy+F_n_u(i,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_u
a_W(i,2)=F_w_u(i,2)*7.0d0/8.0d0*dy+F_e_u(i,2)*1.0d0/8.0d0*dy+D_w*dy
a_E(i,2)=D_e*dy
a_S(i,2)=F_s_u(i,2)*dx+F_n_u(i,2)*1.0d0/8.0d0*dx+D_s*dx
a_N(i,2)=D_n*dx
S_fake(i,2)=(3.0d0/8.0d0*u_fake(i,2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(i,
2)+(-3.0d0/8.0d0&
*u_fake(i,3)*dy+2.0d0/8.0d0*u_fake(i,2)*dy)*F_e_u(i,
2)+(3.0d0/8.0d0*u_fake(i,2)*dx-&
2.0d0/8.0d0*u_fake(i-1,2)*dx-1.0d0/8.0d0*u_fake(i-2,2)*dx)*F_s_u(i,
2)+(-3.0d0/8.0d0*&
u_fake(i+1,2)*dx+2.0d0/8.0d0*u_fake(i,2)*dx)*F_n_u(i,2)+(p_fake(i,
1)-p_fake(i,2))*dy&
+((1.0d0-alpha_u)*a_P(i,2))*u_old(i,2)

!coeficientes de v
b_P(i,2)=(F_e_v(i,2)*dy+F_n_v(i,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_v
b_W(i,2)=F_w_v(i,2)*7.0d0/8.0d0*dy+F_e_v(i,2)*1.0d0/8.0d0*dy+D_w*dy
b_E(i,2)=D_e*dy
b_S(i,2)=F_s_v(i,2)*dx+F_n_v(i,2)*1.0d0/8.0d0*dx+D_s*dx
b_N(i,2)=D_n*dx
S_v_fake(i,2)=(3.0d0/8.0d0*v_fake(i,
2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(i,2)+(-3.0d0/&
8.0d0*v_fake(i,3)*dy+2.0d0/8.0d0*v_fake(i,2)*dy)*F_e_v(i,
2)+(3.0d0/8.0d0*v_fake(i,2)&

*dx-2.0d0/8.0d0*v_fake(i-1,2)*dx-1.0d0/8.0d0*v_fake(i-2,2)*dx)*F_s_v(i,
2)+(-3.0d0/&
8.0d0*v_fake(i+1,2)*dx+2.0d0/8.0d0*v_fake(i,2)*dx)*F_n_v(i,2)+
(p_fake(i-1,2)-&
p_fake(i,2))*dx+((1.0d0-alpha_v)*b_P(i,2))*v_old(i,2)

!coeficientes de p
c_W(i,2)=dy/a_P(i,2)*dy
c_E(i,2)=dy/a_E(i,2)*alpha_u*dy
c_S(i,2)=dx/b_P(i,2)*dx
c_N(i,2)=dx/b_N(i,2)*alpha_v*dx
c_P(i,2)=c_W(i,2)+c_E(i,2)+c_S(i,2)+c_N(i,2)
S_p_fake(i,2)=u_fake(i,2)*dy-u_fake(i,3)*dy+v_fake(i,2)*dx-v_fake(i
+1,2)*dx

end do

!linha 2 - com exceção da coluna 2

do j=3,Nx-1
!coeficientes de u
a_P(2,j)=(F_e_u(2,j)*dy+F_n_u(2,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_u
a_W(2,j)=F_w_u(2,j)*dy+F_e_u(2,j)*1.0d0/8.0d0*dy+D_w*dy
a_E(2,j)=D_e*dy
a_S(2,j)=F_s_u(2,j)*7.0d0/8.0d0*dx+F_n_u(2,j)*1.0d0/8.0d0*dx+D_s*dx
a_N(2,j)=D_n*dx

S_fake(2,j)=(3.0d0/8.0d0*u_fake(2,j)*dy-2.0d0/8.0d0*u_fake(2,j-1)*dy-1.0d0/8.0d0*&
u_fake(2,j-2)*dy)*F_w_u(2,j)+(-3.0d0/8.0d0*u_fake(2,j+1)*dy
+2.0d0/8.0d0*u_fake(2,j)*&
dy)*F_e_u(2,j)
+(3.0d0/8.0d0*u_fake(2,j)*dx-2.0d0/8.0d0*u_G*dx)*F_s_u(2,j)+(-3.0d0/&
8.0d0*u_fake(3,j)*dx+2.0d0/8.0d0*u_fake(2,j)*dx)*F_n_u(2,j)+
(p_fake(2,j-1)-&
p_fake(2,j))*dy+((1.0d0-alpha_u)*a_P(2,j))*u_old(2,j)

!coeficientes de v
b_P(2,j)=(F_e_v(2,j)*dy+F_n_v(2,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_v
b_W(2,j)=F_w_v(2,j)*dy+F_e_v(2,j)*1.0d0/8.0d0*dy+D_w*dy
b_E(2,j)=D_e*dy
b_S(2,j)=F_s_v(2,j)*7.0d0/8.0d0*dx+F_n_v(2,j)*1.0d0/8.0d0*dx+D_s*dx
b_N(2,j)=D_n*dx

S_v_fake(2,j)=(3.0d0/8.0d0*v_fake(2,j)*dy-2.0d0/8.0d0*v_fake(2,j-1)*dy-1.0d0/8.0d0*&
v_fake(2,j-2)*dy)*F_w_v(2,j)+(-3.0d0/8.0d0*v_fake(2,j+1)*dy
+2.0d0/8.0d0*v_fake(2,j)*&
dy)*F_e_v(2,j)
+(3.0d0/8.0d0*v_fake(2,j)*dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,j)+(-3.0d0/&
8.0d0*v_fake(3,j)*dx+2.0d0/8.0d0*v_fake(2,j)*dx)*F_n_v(2,j)+
(p_fake(1,j)-&
p_fake(2,j))*dx+((1.0d0-alpha_v)*b_P(2,j))*v_old(2,j)

!coeficientes de p
c_W(2,j)=dy/a_P(2,j)*dy
c_E(2,j)=dy/a_E(2,j)*alpha_u*dy
c_S(2,j)=dx/b_P(2,j)*dx
c_N(2,j)=dx/b_N(2,j)*alpha_v*dx
c_P(2,j)=c_W(2,j)+c_E(2,j)+c_S(2,j)+c_N(2,j)
S_p_fake(2,j)=u_fake(2,j)*dy-u_fake(2,j+1)*dy+v_fake(2,j)*dx-
v_fake(3,j)*dx

end do

!coluna 2 - linha 2

!coeficientes de u
a_P(2,2)=(F_e_u(2,2)*dy+F_n_u(2,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_u
a_W(2,2)=F_w_u(2,2)*7.0d0/8.0d0*dy+F_e_u(2,2)*1.0d0/8.0d0*dy
+D_w*dy
a_E(2,2)=D_e*dy
a_S(2,2)=F_s_u(2,2)*7.0d0/8.0d0*dx+F_n_u(2,2)*1.0d0/8.0d0*dx
+D_s*dx
a_N(2,2)=D_n*dx

S_fake(2,2)=(3.0d0/8.0d0*u_fake(2,2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(2,2)+
(-3.0d0/&
8.0d0*u_fake(2,3)*dy
+2.0d0/8.0d0*u_fake(2,2)*dy)*F_e_u(2,2)+(3.0d0/8.0d0*&
u_fake(2,2)*dx-2.0d0/8.0d0*u_G*dx)*F_s_u(2,2)+
(-3.0d0/8.0d0*u_fake(3,2)*dx+2.0d0&
/8.0d0*u_fake(2,2)*dx)*F_n_u(2,2)+(p_fake(2,1)-p_fake(2,2))*dy+
((1.0d0-alpha_u)*&
a_P(2,2))*u_old(2,2)

!coeficientes de v
b_P(2,2)=(F_e_v(2,2)*dy+F_n_v(2,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
alpha_v
b_W(2,2)=F_w_v(2,2)*7.0d0/8.0d0*dy+F_e_v(2,2)*1.0d0/8.0d0*dy
+D_w*dy
b_E(2,2)=D_e*dy
b_S(2,2)=F_s_v(2,2)*7.0d0/8.0d0*dx+F_n_v(2,2)*1.0d0/8.0d0*dx
+D_s*dx
b_N(2,2)=D_n*dx

S_v_fake(2,2)=(3.0d0/8.0d0*v_fake(2,2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(2,2)+
(-3.0d0&
/8.0d0*v_fake(2,3)*dy
+2.0d0/8.0d0*v_fake(2,2)*dy)*F_e_v(2,2)+(3.0d0/8.0d0*&
v_fake(2,2)*dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,2)+
(-3.0d0/8.0d0*v_fake(3,2)*dx+2.0d0&
/8.0d0*v_fake(2,2)*dx)*F_n_v(2,2)+(p_fake(1,2)-p_fake(2,2))*dx+
((1.0d0-alpha_v)*&
b_P(2,2))*v_old(2,2)

!coeficientes de p
c_W(2,2)=dy/a_P(2,2)*dy
c_E(2,2)=dy/a_E(2,2)*alpha_u*dy
c_S(2,2)=dx/b_P(2,2)*dx
c_N(2,2)=dx/b_N(2,2)*alpha_v*dx
c_P(2,2)=c_W(2,2)+c_E(2,2)+c_S(2,2)+c_N(2,2)
S_p_fake(2,2)=u_fake(2,2)*dy-u_fake(2,3)*dy+v_fake(2,2)*dx-
v_fake(3,2)*dx


!fronteira W - com exceção da linha 2

do i=3,Ny-1
!coeficientes de u
a_P(i,1)=(F_e_u(i,1)*7.0d0/8.0d0*dy+F_n_u(i,1)*dx+D_e*dy+D_s*dx
+D_n*dx+D_EE*9.0d0/&
3.0d0*dy)/alpha_u
a_W(i,1)=1.0d-20 !0.0d0
a_E(i,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
a_S(i,1)=F_s_u(i,1)*dx+F_n_u(i,1)*1.0d0/8.0d0*dx+D_s*dx
a_N(i,1)=D_n*dx
S_fake(i,1)=(u_E*dy)*F_w_u(i,1)+(-3.0d0/8.0d0*u_fake(i,2)*dy
+2.0d0/8.0d0*u_E*dy)*&
F_e_u(i,1)+(3.0d0/8.0d0*u_fake(i,
1)*dx-2.0d0/8.0d0*u_fake(i-1,1)*dx-1.0d0/8.0d0*&
u_fake(i-2,1)*dx)*F_s_u(i,1)+(-3.0d0/8.0d0*u_fake(i+1,1)*dx
+2.0d0/8.0d0*u_fake(i,1)*&
dx)*F_n_u(i,1)+D_EE*8.0d0/3.0d0*u_E*dy+((1.0d0-alpha_u)*a_P(i,
1))*u_old(i,1) !Falta o termo de pressão.


! a_P(i,1)=1.0d0
! a_W(i,1)=1.0d-20
! a_E(i,1)=1.0d-20
! a_S(i,1)=1.0d-20
! a_N(i,1)=1.0d-20
! S_fake(i,1)=u_E






!coeficientes de v
b_P(i,1)=(F_e_v(i,1)*7.0d0/8.0d0*dy+F_n_v(i,1)*dx+D_e*dy+D_s*dx
+D_n*dx+D_I*9.0d0/&
3.0d0*dy)/alpha_v
b_W(i,1)=1.0d-20 !0.0d0
b_E(i,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
b_S(i,1)=F_s_v(i,1)*dx+F_n_v(i,1)*1.0d0/8.0d0*dx+D_s*dx
b_N(i,1)=D_n*dx
S_v_fake(i,1)=(v_I*dy)*F_w_v(i,1)+(-3.0d0/8.0d0*v_fake(i,2)*dy
+2.0d0/8.0d0*v_I*dy)*&
F_e_v(i,1)+(3.0d0/8.0d0*v_fake(i,
1)*dx-2.0d0/8.0d0*v_fake(i-1,1)*dx-1.0d0/8.0d0*&
v_fake(i-2,1)*dx)*F_s_v(i,1)+(-3.0d0/8.0d0*v_fake(i+1,1)*dx
+2.0d0/8.0d0*v_fake(i,1)*&
dx)*F_n_v(i,1)+(p_fake(i-1,1)-p_fake(i,1))*dx+D_I*8.0d0/3.0d0*v_I*dy
+((1.0d0-alpha_v&
)*b_P(i,1))*v_old(i,1)

! b_P(i,1)=1.0d0
! b_W(i,1)=1.0d-20
! b_E(i,1)=1.0d-20
! b_S(i,1)=1.0d-20
! b_N(i,1)=1.0d-20
! S_v_fake(i,1)=v_I




!coeficientes de p

c_W(i,1)=dy/a_P(i,1)*dy
c_E(i,1)=dy/a_E(i,1)*alpha_u*dy
c_S(i,1)=dx/b_P(i,1)*dx
c_N(i,1)=dx/b_N(i,1)*alpha_v*dx
c_P(i,1)=c_W(i,1)+c_E(i,1)+c_S(i,1)+c_N(i,1)
S_p_fake(i,1)=u_fake(i,1)*dy-u_fake(i,2)*dy+v_fake(i,1)*dx-v_fake(i
+1,1)*dx

end do

!fronteira W - linha 2

!coeficientes de u
a_P(2,1)=(F_e_u(2,1)*7.0d0/8.0d0*dy+F_n_u(2,1)*dx+D_e*dy+D_s*dx
+D_n*dx+D_EE*9.0d0/&
3.0d0*dy)/alpha_u
a_W(2,1)=1.0d-20 !0.0d0
a_E(2,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
a_S(2,1)=F_s_u(2,1)*7.0d0/8.0d0*dx+F_n_u(2,1)*1.0d0/8.0d0*dx+D_s*dx
a_N(2,1)=D_n*dx
S_fake(2,1)=(u_E*dy)*F_w_u(2,1)+(-3.0d0/8.0d0*u_fake(2,2)*dy
+2.0d0/8.0d0*u_E*dy)*&

F_e_u(2,1)+(3.0d0/8.0d0*u_fake(2,1)*dx-2.0d0/8.0d0*u_G*dx)*F_s_u(2,1)+
(-3.0d0/8.0d0*&
u_fake(3,1)*dx
+2.0d0/8.0d0*u_fake(2,1)*dx)*F_n_u(2,1)+D_EE*8.0d0/3.0d0*u_E*dy+((&
1.0d0-alpha_u)*a_P(2,1))*u_old(2,1) !Falta o termo de pressão.


! a_P(2,1)=1.0d0
! a_W(2,1)=1.0d-20
! a_E(2,1)=1.0d-20
! a_S(2,1)=1.0d-20
! a_N(2,1)=1.0d-20
! S_fake(2,1)=u_E




!coeficientes de v
b_P(2,1)=(F_e_v(2,1)*7.0d0/8.0d0*dy+F_n_v(2,1)*dx+D_e*dy+D_s*dx
+D_n*dx+D_I*9.0d0/&
3.0d0*dy)/alpha_v
b_W(2,1)=1.0d-20 !0.0d0
b_E(2,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
b_S(2,1)=F_s_v(2,1)*7.0d0/8.0d0*dx+F_n_v(2,1)*1.0d0/8.0d0*dx+D_s*dx
b_N(2,1)=D_n*dx
S_v_fake(2,1)=(v_I*dy)*F_w_v(2,1)+(-3.0d0/8.0d0*v_fake(2,2)*dy
+2.0d0/8.0d0*v_I*dy)*&

F_e_v(2,1)+(3.0d0/8.0d0*v_fake(2,1)*dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,1)+
(-3.0d0/8.0d0*&
v_fake(3,1)*dx+2.0d0/8.0d0*v_fake(2,1)*dx)*F_n_v(2,1)+(p_fake(1,1)-
p_fake(2,1))*dx+&
D_I*8.0d0/3.0d0*v_I*dy+((1.0d0-alpha_v)*b_P(2,1))*v_old(2,1)

! b_P(2,1)=1.0d0
! b_W(2,1)=1.0d-20
! b_E(2,1)=1.0d-20
! b_S(2,1)=1.0d-20
! b_N(2,1)=1.0d-20
! S_v_fake(2,1)=v_I




!coeficientes de p
c_W(2,1)=dy/a_P(2,1)*dy
c_E(2,1)=dy/a_E(2,1)*alpha_u*dy
c_S(2,1)=dx/b_P(2,1)*dx
c_N(2,1)=dx/b_N(2,1)*alpha_v*dx
c_P(2,1)=c_W(2,1)+c_E(2,1)+c_S(2,1)+c_N(2,1)
S_p_fake(2,1)=u_fake(2,1)*dy-u_fake(2,2)*dy+v_fake(2,1)*dx-
v_fake(3,1)*dx



!fronteira E - com exceção da linha 2

do i=3,Ny-1
!coeficientes de u
a_P(i,Nx)=(F_n_u(i,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_F*9.0d0/3.0d0*dy)/
alpha_u
a_W(i,Nx)=F_w_u(i,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
a_E(i,Nx)=1.0d-20 !0.0d0
a_S(i,Nx)=F_s_u(i,Nx)*dx+F_n_u(i,Nx)*1.0d0/8.0d0*dx+D_s*dx
a_N(i,Nx)=D_n*dx

S_fake(i,Nx)=(3.0d0/8.0d0*u_fake(i,Nx)*dy-2.0d0/8.0d0*u_fake(i,Nx-1)*dy-1..0d0/8.0d0*&
u_fake(i,Nx-2)*dy)*F_w_u(i,Nx)+(-u_F(i)*dy)*F_e_u(i,Nx)
+(3.0d0/8.0d0*u_fake(i,Nx)*dx&
-2.0d0/8.0d0*u_fake(i-1,Nx)*dx-1.0d0/8.0d0*u_fake(i-2,Nx)*dx)*F_s_u(i,Nx)
+(-3.0d0/&
8.0d0*u_fake(i+1,Nx)*dx+2.0d0/8.0d0*u_fake(i,Nx)*dx)*F_n_u(i,Nx)+
(p_fake(i,Nx-1)-&
p_fake(i,Nx))*dy+D_F*8.0d0/3.0d0*u_F(i)*dy+((1.0d0-
alpha_u)*a_P(i,Nx))*u_old(i,Nx)


! a_P(i,Nx)=1.0d0
! a_W(i,Nx)=1.0d-20
! a_E(i,Nx)=1.0d-20
! a_S(i,Nx)=1.0d-20
! a_N(i,Nx)=1.0d-20
! S_fake(i,Nx)=u_F(i)





!coeficientes de v
b_P(i,Nx)=(F_n_v(i,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_J*9.0d0/3.0d0*dy)/
alpha_v
b_W(i,Nx)=F_w_v(i,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
b_E(i,Nx)=1.0d-20 !0.0d0
b_S(i,Nx)=F_s_v(i,Nx)*dx+F_n_v(i,Nx)*1.0d0/8.0d0*dx+D_s*dx
b_N(i,Nx)=D_n*dx

S_v_fake(i,Nx)=(3.0d0/8.0d0*v_fake(i,Nx)*dy-2.0d0/8.0d0*v_fake(i,Nx-1)*dy-1.0d0/
&
8.0d0*v_fake(i,Nx-2)*dy)*F_w_v(i,Nx)+(-v_J*dy)*F_e_v(i,Nx)
+(3.0d0/8.0d0*v_fake(i,Nx)&

*dx-2.0d0/8.0d0*v_fake(i-1,Nx)*dx-1.0d0/8.0d0*v_fake(i-2,Nx)*dx)*F_s_v(i,Nx)
+(-3.0d0&
/8.0d0*v_fake(i+1,Nx)*dx+2.0d0/8.0d0*v_fake(i,Nx)*dx)*F_n_v(i,Nx)+
(p_fake(i-1,Nx)-&
p_fake(i,Nx))*dx+D_J*8.0d0/3.0d0*v_J*dy+((1.0d0-
alpha_v)*b_P(i,Nx))*v_old(i,Nx)


! b_P(i,Nx)=1.0d0
! b_W(i,Nx)=1.0d-20
! b_E(i,Nx)=1.0d-20
! b_S(i,Nx)=1.0d-20
! b_N(i,Nx)=1.0d-20
! S_v_fake(i,Nx)=v_J




!coeficientes de p
c_W(i,Nx)=dy/a_P(i,Nx)*dy
c_E(i,Nx)=dy/a_E(i,Nx)*alpha_u*dy
c_S(i,Nx)=dx/b_P(i,Nx)*dx
c_N(i,Nx)=dx/b_N(i,Nx)*alpha_v*dx
c_P(i,Nx)=c_W(i,Nx)+c_E(i,Nx)+c_S(i,Nx)+c_N(i,Nx)
S_p_fake(i,Nx)=v_fake(i,Nx)*dx-v_fake(i+1,Nx)*dx !Faltam termos de
u.

end do

!fronteira E - linha 2

!coeficientes de u
a_P(2,Nx)=(F_n_u(2,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_F*9.0d0/3.0d0*dy)/
alpha_u
a_W(2,Nx)=F_w_u(2,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
a_E(2,Nx)=1.0d-20 !0.0d0
a_S(2,Nx)=F_s_u(2,Nx)*7.0d0/8.0d0*dx+F_n_u(2,Nx)*1.0d0/8.0d0*dx
+D_s*dx
a_N(2,Nx)=D_n*dx

S_fake(2,Nx)=(3.0d0/8.0d0*u_fake(2,Nx)*dy-2.0d0/8.0d0*u_fake(2,Nx-1)*dy-1..0d0/8.0d0*&
u_fake(2,Nx-2)*dy)*F_w_u(2,Nx)+(-u_F(2)*dy)*F_e_u(2,Nx)
+(3.0d0/8.0d0*u_fake(2,Nx)*dx&
-2.0d0/8.0d0*u_G*dx)*F_s_u(2,Nx)+(-3.0d0/8.0d0*u_fake(3,Nx)*dx
+2.0d0/8.0d0*&
u_fake(2,Nx)*dx)*F_n_u(2,Nx)+(p_fake(2,Nx-1)-p_fake(2,Nx))*dy
+D_F*8.0d0/3.0d0*u_F(2)&
*dy+((1.0d0-alpha_u)*a_P(2,Nx))*u_old(2,Nx)


! a_P(2,Nx)=1.0d0
! a_W(2,Nx)=1.0d-20
! a_E(2,Nx)=1.0d-20
! a_S(2,Nx)=1.0d-20
! a_N(2,Nx)=1.0d-20
! S_fake(2,Nx)=u_F(2)



!coeficientes de v
b_P(2,Nx)=(F_n_v(2,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_J*9.0d0/3.0d0*dy)/
alpha_v
b_W(2,Nx)=F_w_v(2,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
b_E(2,Nx)=1.0d-20 !0.0d0
b_S(2,Nx)=F_s_v(2,Nx)*7.0d0/8.0d0*dx+F_n_v(2,Nx)*1.0d0/8.0d0*dx
+D_s*dx
b_N(2,Nx)=D_n*dx

S_v_fake(2,Nx)=(3.0d0/8.0d0*v_fake(2,Nx)*dy-2.0d0/8.0d0*v_fake(2,Nx-1)*dy-1.0d0/
&
8.0d0*v_fake(2,Nx-2)*dy)*F_w_v(2,Nx)+(-v_J*dy)*F_e_v(2,Nx)
+(3.0d0/8.0d0*v_fake(2,Nx)&
*dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,Nx)+(-3.0d0/8.0d0*v_fake(3,Nx)*dx
+2.0d0/8.0d0*&
v_fake(2,Nx)*dx)*F_n_v(2,Nx)+(p_fake(1,Nx)-p_fake(2,Nx))*dx
+D_J*8.0d0/3.0d0*v_J*dy+(&
(1.0d0-alpha_v)*b_P(2,Nx))*v_old(2,Nx)


! b_P(2,Nx)=1.0d0
! b_W(2,Nx)=1.0d-20
! b_E(2,Nx)=1.0d-20
! b_S(2,Nx)=1.0d-20
! b_N(2,Nx)=1.0d-20
! S_v_fake(2,Nx)=v_J

!coeficientes de p
c_W(2,Nx)=dy/a_P(2,Nx)*dy
c_E(2,Nx)=dy/a_E(2,Nx)*alpha_u*dy
c_S(2,Nx)=dx/b_P(2,Nx)*dx
c_N(2,Nx)=dx/b_N(2,Nx)*alpha_v*dx
c_P(2,Nx)=c_W(2,Nx)+c_E(2,Nx)+c_S(2,Nx)+c_N(2,Nx)
S_p_fake(2,Nx)=v_fake(2,Nx)*dx-v_fake(3,Nx)*dx !Faltam termos de u.


!fronteira S - com exceção da coluna 2

do j=3,Nx-1
!coeficientes de u
a_P(1,j)=(F_e_u(1,j)*dy+F_n_u(1,j)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy
+D_n*dx+D_G*9.0d0/&
3.0d0*dx)/alpha_u
a_W(1,j)=F_w_u(1,j)*dy+F_e_u(1,j)*1.0d0/8.0d0*dy+D_w*dy
a_E(1,j)=D_e*dy
a_S(1,j)=1.0d-20 !0.0d0
a_N(1,j)=D_n*dx+D_G*1.0d0/3.0d0*dx

S_fake(1,j)=(3.0d0/8.0d0*u_fake(1,j)*dy-2.0d0/8.0d0*u_fake(1,j-1)*dy-1.0d0/8.0d0*&
u_fake(1,j-2)*dy)*F_w_u(1,j)+(-3.0d0/8.0d0*u_fake(1,j+1)*dy
+2.0d0/8.0d0*u_fake(1,j)*&
dy)*F_e_u(1,j)+(u_G*dx)*F_s_u(1,j)+(-3.0d0/8.0d0*u_fake(2,j)*dx
+2.0d0/8.0d0*u_G*dx)*&
F_n_u(1,j)+(p_fake(1,j-1)-p_fake(1,j))*dy+D_G*8.0d0/3.0d0*u_G*dx+
((1.0d0-alpha_u)*&
a_P(1,j))*u_old(1,j)


! a_P(1,j)=1.0d0
! a_W(1,j)=1.0d-20
! a_E(1,j)=1.0d-20
! a_S(1,j)=1.0d-20
! a_N(1,j)=1.0d-20
! S_fake(1,j)=u_G






!coeficientes de v
b_P(1,j)=(F_e_v(1,j)*dy+F_n_v(1,j)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy
+D_n*dx+D_K*9.0d0/&
3.0d0*dx)/alpha_v
b_W(1,j)=F_w_v(1,j)*dy+F_e_v(1,j)*1.0d0/8.0d0*dy+D_w*dy
b_E(1,j)=D_e*dy
b_S(1,j)=1.0d-20 !0.0d0
b_N(1,j)=D_n*dx+D_K*1.0d0/3.0d0*dx

S_v_fake(1,j)=(3.0d0/8.0d0*v_fake(1,j)*dy-2.0d0/8.0d0*v_fake(1,j-1)*dy-1.0d0/8.0d0*&
v_fake(1,j-2)*dy)*F_w_v(1,j)+(-3.0d0/8.0d0*v_fake(1,j+1)*dy
+2.0d0/8.0d0*v_fake(1,j)*&
dy)*F_e_v(1,j)+(v_K*dx)*F_s_v(1,j)+(-3.0d0/8.0d0*v_fake(2,j)*dx
+2.0d0/8.0d0*v_K*dx)*&
F_n_v(1,j)+D_K*8.0d0/3.0d0*v_K*dx+((1.0d0-
alpha_v)*b_P(1,j))*v_old(1,j) !Falta o termo de pressão.

! b_P(1,j)=1.0d0
! b_W(1,j)=1.0d-20
! b_E(1,j)=1.0d-20
! b_S(1,j)=1.0d-20
! b_N(1,j)=1.0d-20
! S_v_fake(1,j)=v_K





!coeficientes de p
c_W(1,j)=dy/a_P(1,j)*dy
c_E(1,j)=dy/a_E(1,j)*alpha_u*dy
c_S(1,j)=dx/b_P(1,j)*dx
c_N(1,j)=dx/b_N(1,j)*alpha_v*dx
c_P(1,j)=c_W(1,j)+c_E(1,j)+c_S(1,j)+c_N(1,j)
S_p_fake(1,j)=u_fake(1,j)*dy-u_fake(1,j+1)*dy+v_fake(1,j)*dx-
v_fake(2,j)*dx

end do

!fronteira S - coluna 2

!coeficientes de u
a_P(1,2)=(F_e_u(1,2)*dy+F_n_u(1,2)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy+D_n*dx
+D_G*9.0d0/3.0d0*dx)/&
alpha_u
a_W(1,2)=F_w_u(1,2)*7.0d0/8.0d0*dy+F_e_u(1,2)*1.0d0/8.0d0*dy+D_w*dy
a_E(1,2)=D_e*dy
a_S(1,2)=1.0d-20 !0.0d0
a_N(1,2)=D_n*dx+D_G*1.0d0/3.0d0*dx

S_fake(1,2)=(3.0d0/8.0d0*u_fake(1,2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(1,2)+
(-3.0d0/8.0d0*&
u_fake(1,3)*dy+2.0d0/8.0d0*u_fake(1,2)*dy)*F_e_u(1,2)+
(u_G*dx)*F_s_u(1,2)+(-3.0d0/8.0d0*&
u_fake(2,2)*dx+2.0d0/8.0d0*u_G*dx)*F_n_u(1,2)+(p_fake(1,1)-
p_fake(1,2))*dy+D_G*8.0d0/3.0d0*&
u_G*dx+((1.0d0-alpha_u)*a_P(1,2))*u_old(1,2)

! a_P(1,2)=1.0d0
! a_W(1,2)=1.0d-20
! a_E(1,2)=1.0d-20
! a_S(1,2)=1.0d-20
! a_N(1,2)=1.0d-20
! S_fake(1,2)=u_G








!coeficientes de v
b_P(1,2)=(F_e_v(1,2)*dy+F_n_v(1,2)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy+D_n*dx
+D_K*9.0d0/3.0d0*dx)/&
alpha_v
b_W(1,2)=F_w_v(1,2)*7.0d0/8.0d0*dy+F_e_v(1,2)*1.0d0/8.0d0*dy+D_w*dy
b_E(1,2)=D_e*dy
b_S(1,2)=1.0d-20 !0.0d0
b_N(1,2)=D_n*dx+D_K*1.0d0/3.0d0*dx

S_v_fake(1,2)=(3.0d0/8.0d0*v_fake(1,2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(1,2)+
(-3.0d0/8.0d0*&
v_fake(1,3)*dy+2.0d0/8.0d0*v_fake(1,2)*dy)*F_e_v(1,2)+
(v_K*dx)*F_s_v(1,2)+(-3.0d0/8.0d0*&
v_fake(2,2)*dx+2.0d0/8.0d0*v_K*dx)*F_n_v(1,2)+D_K*8.0d0/3.0d0*v_K*dx+
((1.0d0-alpha_v)*&
b_P(1,2))*v_old(1,2) !Falta o termo de pressão.

! b_P(1,2)=1.0d0
! b_W(1,2)=1.0d-20
! b_E(1,2)=1.0d-20
! b_S(1,2)=1.0d-20
! b_N(1,2)=1.0d-20
! S_v_fake(1,2)=v_K





!coeficientes de p
c_W(1,2)=dy/a_P(1,2)*dy
c_E(1,2)=dy/a_E(1,2)*alpha_u*dy
c_S(1,2)=dx/b_P(1,2)*dx
c_N(1,2)=dx/b_N(1,2)*alpha_v*dx
c_P(1,2)=c_W(1,2)+c_E(1,2)+c_S(1,2)+c_N(1,2)
S_p_fake(1,2)=u_fake(1,2)*dy-u_fake(1,3)*dy+v_fake(1,2)*dx-
v_fake(2,2)*dx


!fronteira N - com exceção da coluna 2

do j=3,Nx-1
!coeficientes de u
a_P(Ny,j)=(F_e_u(Ny,j)*dy+D_w*dy+D_e*dy+D_s*dx+D_H*9.0d0/3.0d0*dx)/
alpha_u
a_W(Ny,j)=F_w_u(Ny,j)*dy+F_e_u(Ny,j)*1.0d0/8.0d0*dy+D_w*dy
a_E(Ny,j)=D_e*dy
a_S(Ny,j)=F_s_u(Ny,j)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
a_N(Ny,j)=1.0d-20 !0.0d0

S_fake(Ny,j)=(3.0d0/8.0d0*u_fake(Ny,j)*dy-2.0d0/8.0d0*u_fake(Ny,j-1)*dy-1..0d0/8.0d0*&
u_fake(Ny,j-2)*dy)*F_w_u(Ny,j)+(-3.0d0/8.0d0*u_fake(Ny,j+1)*dy
+2.0d0/8.0d0*&
u_fake(Ny,j)*dy)*F_e_u(Ny,j)
+(3.0d0/8.0d0*u_fake(Ny,j)*dx-2.0d0/8.0d0*u_fake(Ny-1,j)&
*dx-1.0d0/8.0d0*u_fake(Ny-2,j)*dx)*F_s_u(Ny,j)+(-u_H*dx)*F_n_u(Ny,j)
+(p_fake(Ny,j-1)&
-p_fake(Ny,j))*dy+D_H*8.0d0/3.0d0*u_H*dx+((1.0d0-
alpha_u)*a_P(Ny,j))*u_old(Ny,j)


! a_P(Ny,j)=1.0d0
! a_W(Ny,j)=1.0d-20
! a_E(Ny,j)=1.0d-20
! a_S(Ny,j)=1.0d-20
! a_N(Ny,j)=1.0d-20
! S_fake(Ny,j)=u_H




!coeficientes de v
b_P(Ny,j)=(F_e_v(Ny,j)*dy+D_w*dy+D_e*dy+D_s*dx+D_L*9.0d0/3.0d0*dx)/
alpha_v
b_W(Ny,j)=F_w_v(Ny,j)*dy+F_e_v(Ny,j)*1.0d0/8.0d0*dy+D_w*dy
b_E(Ny,j)=D_e*dy
b_S(Ny,j)=F_s_v(Ny,j)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
b_N(Ny,j)=1.0d-20 !0.0d0

S_v_fake(Ny,j)=(3.0d0/8.0d0*v_fake(Ny,j)*dy-2.0d0/8.0d0*v_fake(Ny,j-1)*dy-1.0d0/
&
8.0d0*v_fake(Ny,j-2)*dy)*F_w_v(Ny,j)+(-3.0d0/8.0d0*v_fake(Ny,j+1)*dy
+2.0d0/8.0d0*&
v_fake(Ny,j)*dy)*F_e_v(Ny,j)
+(3.0d0/8.0d0*v_fake(Ny,j)*dx-2.0d0/8.0d0*v_fake(Ny-1,j)&
*dx-1.0d0/8.0d0*v_fake(Ny-2,j)*dx)*F_s_v(Ny,j)+(-v_L*dx)*F_n_v(Ny,j)
+(p_fake(Ny-1,j)&
-p_fake(Ny,j))*dx+D_L*8.0d0/3.0d0*v_L*dx+((1.0d0-
alpha_v)*b_P(Ny,j))*v_old(Ny,j)


! b_P(Ny,j)=1.0d0
! b_W(Ny,j)=1.0d-20
! b_E(Ny,j)=1.0d-20
! b_S(Ny,j)=1.0d-20
! b_N(Ny,j)=1.0d-20
! S_v_fake(Ny,j)=v_L





!coeficientes de p
c_W(Ny,j)=dy/a_P(Ny,j)*dy
c_E(Ny,j)=dy/a_E(Ny,j)*alpha_u*dy
c_S(Ny,j)=dx/b_P(Ny,j)*dx
c_N(Ny,j)=dx/b_N(Ny,j)*alpha_v*dx
c_P(Ny,j)=c_W(Ny,j)+c_E(Ny,j)+c_S(Ny,j)+c_N(Ny,j)
S_p_fake(Ny,j)=u_fake(Ny,j)*dy-u_fake(Ny,j+1)*dy !Faltam termos de
v.


end do

!fronteira N - coluna 2

!coeficientes de u
a_P(Ny,2)=(F_e_u(Ny,2)*dy+D_w*dy+D_e*dy+D_s*dx+D_H*9.0d0/3.0d0*dx)/
alpha_u
a_W(Ny,2)=F_w_u(Ny,2)*7.0d0/8.0d0*dy+F_e_u(Ny,2)*1.0d0/8.0d0*dy
+D_w*dy
a_E(Ny,2)=D_e*dy
a_S(Ny,2)=F_s_u(Ny,2)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
a_N(Ny,2)=1.0d-20 !0.0d0
S_fake(Ny,2)=(3.0d0/8.0d0*u_fake(Ny,
2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(Ny,2)+(-3.0d0/&
8.0d0*u_fake(Ny,3)*dy+2.0d0/8.0d0*u_fake(Ny,2)*dy)*F_e_u(Ny,
2)+(3.0d0/8.0d0*&
u_fake(Ny,
2)*dx-2.0d0/8.0d0*u_fake(Ny-1,2)*dx-1.0d0/8.0d0*u_fake(Ny-2,2)*dx)*&
F_s_u(Ny,2)+(-u_H*dx)*F_n_u(Ny,2)+(p_fake(Ny,1)-p_fake(Ny,2))*dy
+D_H*8.0d0/3.0d0*u_H&
*dx+((1.0d0-alpha_u)*a_P(Ny,2))*u_old(Ny,2)

! a_P(Ny,2)=1.0d0
! a_W(Ny,2)=1.0d-20
! a_E(Ny,2)=1.0d-20
! a_S(Ny,2)=1.0d-20
! a_N(Ny,2)=1.0d-20
! S_fake(Ny,2)=u_H




!coeficientes de v
b_P(Ny,2)=(F_e_v(Ny,2)*dy+D_w*dy+D_e*dy+D_s*dx+D_L*9.0d0/3.0d0*dx)/
alpha_v
b_W(Ny,2)=F_w_v(Ny,2)*7.0d0/8.0d0*dy+F_e_v(Ny,2)*1.0d0/8.0d0*dy
+D_w*dy
b_E(Ny,2)=D_e*dy
b_S(Ny,2)=F_s_v(Ny,2)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
b_N(Ny,2)=1.0d-20 !0.0d0
S_v_fake(Ny,2)=(3.0d0/8.0d0*v_fake(Ny,
2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(Ny,2)+(-3.0d0/&
8.0d0*v_fake(Ny,3)*dy+2.0d0/8.0d0*v_fake(Ny,2)*dy)*F_e_v(Ny,
2)+(3.0d0/8.0d0*&
v_fake(Ny,
2)*dx-2.0d0/8.0d0*v_fake(Ny-1,2)*dx-1.0d0/8.0d0*v_fake(Ny-2,2)*dx)*&
F_s_v(Ny,2)+(-v_L*dx)*F_n_v(Ny,2)+(p_fake(Ny-1,2)-p_fake(Ny,2))*dx
+D_L*8.0d0/3.0d0*&
v_L*dx+((1.0d0-alpha_v)*b_P(Ny,2))*v_old(Ny,2)

! b_P(Ny,2)=1.0d0
! b_W(Ny,2)=1.0d-20
! b_E(Ny,2)=1.0d-20
! b_S(Ny,2)=1.0d-20
! b_N(Ny,2)=1.0d-20
! S_v_fake(Ny,2)=v_L







!coeficientes de p
c_W(Ny,2)=dy/a_P(Ny,2)*dy
c_E(Ny,2)=dy/a_E(Ny,2)*alpha_u*dy
c_S(Ny,2)=dx/b_P(Ny,2)*dx
c_N(Ny,2)=dx/b_N(Ny,2)*alpha_v*dx
c_P(Ny,2)=c_W(Ny,2)+c_E(Ny,2)+c_S(Ny,2)+c_N(Ny,2)
S_p_fake(Ny,2)=u_fake(Ny,2)*dy-u_fake(Ny,3)*dy !Faltam termos de v.


!fronteira SW

!coeficientes de u
a_P(1,1)=(F_e_u(1,1)*7.0d0/8.0d0*dy+F_n_u(1,1)*7.0d0/8.0d0*dx+D_e*dy
+D_n*dx+D_EE*&
9.0d0/3.0d0*dy+D_G*9.0d0/3.0d0*dx)/alpha_u
a_W(1,1)=1.0d-20 !0.0d0
a_E(1,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
a_S(1,1)=1.0d-20 !0.0d0
a_N(1,1)=D_n*dx+D_G*1.0d0/3.0d0*dx
S_fake(1,1)=(u_E*dy)*F_w_u(1,1)+(-3.0d0/8.0d0*u_fake(1,2)*dy
+2.0d0/8.0d0*u_E*dy)*&
F_e_u(1,1)+(u_G*dx)*F_s_u(1,1)+(-3.0d0/8.0d0*u_fake(2,1)*dx
+2.0d0/8.0d0*u_G*dx)*&
F_n_u(1,1)+D_EE*8.0d0/3.0d0*u_E*dy+D_G*8.0d0/3.0d0*u_G*dx+((1.0d0-
alpha_u)*a_P(1,1))&
*u_old(1,1) !Falta o termo de pressão.


! a_P(1,1)=1.0d0
! a_W(1,1)=1.0d-20
! a_E(1,1)=1.0d-20
! a_S(1,1)=1.0d-20
! a_N(1,1)=1.0d-20
! S_fake(1,1)=1.0d-30








!coeficientes de v
b_P(1,1)=(F_e_v(1,1)*7.0d0/8.0d0*dy+F_n_v(1,1)*7.0d0/8.0d0*dx+D_e*dy
+D_n*dx+D_I*&
9.0d0/3.0d0*dy+D_K*9.0d0/3.0d0*dx)/alpha_v
b_W(1,1)=1.0d-20 !0.0d0
b_E(1,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
b_S(1,1)=1.0d-20 !0.0d0
b_N(1,1)=D_n*dx+D_K*1.0d0/3.0d0*dx
S_v_fake(1,1)=(v_I*dy)*F_w_v(1,1)+(-3.0d0/8.0d0*v_fake(1,2)*dy
+2.0d0/8.0d0*v_I*dy)*&
F_e_v(1,1)+(v_K*dx)*F_s_v(1,1)+(-3.0d0/8.0d0*v_fake(2,1)*dx
+2.0d0/8.0d0*v_K*dx)*&
F_n_v(1,1)+D_I*8.0d0/3.0d0*v_I*dy+D_K*8.0d0/3.0d0*v_K*dx+((1.0d0-
alpha_v)*b_P(1,1))&
*v_old(1,1) !Falta o termo de pressão.

! b_P(1,1)=1.0d0
! b_W(1,1)=1.0d-20
! b_E(1,1)=1.0d-20
! b_S(1,1)=1.0d-20
! b_N(1,1)=1.0d-20
! S_v_fake(1,1)=1.0d-30




!coeficientes de p
c_W(1,1)=dy/a_P(1,1)*dy
c_E(1,1)=dy/a_E(1,1)*alpha_u*dy
c_S(1,1)=dx/b_P(1,1)*dx
c_N(1,1)=dx/b_N(1,1)*alpha_v*dx
c_P(1,1)=c_W(1,1)+c_E(1,1)+c_S(1,1)+c_N(1,1)
S_p_fake(1,1)=u_fake(1,1)*dy-u_fake(1,2)*dy+v_fake(1,1)*dx-
v_fake(2,1)*dx


!fronteira NW

!coeficientes de u
a_P(Ny,1)=(F_e_u(Ny,1)*7.0d0/8.0d0*dy+D_e*dy+D_s*dx
+D_EE*9.0d0/3.0d0*dy+D_H*9.0d0/&
3.0d0*dx)/alpha_u
a_W(Ny,1)=1.0d-20 !0.0d0
a_E(Ny,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
a_S(Ny,1)=F_s_u(Ny,1)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
a_N(Ny,1)=1.0d-20 !0.0d0
S_fake(Ny,1)=(u_E*dy)*F_w_u(Ny,1)+(-3.0d0/8.0d0*u_fake(Ny,2)*dy
+2.0d0/8.0d0*u_E*dy)*&
F_e_u(Ny,1)+(3.0d0/8.0d0*u_fake(Ny,
1)*dx-2.0d0/8.0d0*u_fake(Ny-1,1)*dx-1.0d0/8.0d0*&
u_fake(Ny-2,1)*dx)*F_s_u(Ny,1)+(-u_H*dx)*F_n_u(Ny,
1)+D_EE*8.0d0/3.0d0*u_E*dy+D_H*&
8.0d0/3.0d0*u_H*dx+((1.0d0-alpha_u)*a_P(Ny,1))*u_old(Ny,1) ! Falta
o termo de pressão.


! a_P(Ny,1)=1.0d0
! a_W(Ny,1)=1.0d-20
! a_E(Ny,1)=1.0d-20
! a_S(Ny,1)=1.0d-20
! a_N(Ny,1)=1.0d-20
! S_fake(Ny,1)=1.0d-30





!coeficientes de v
b_P(Ny,1)=(F_e_v(Ny,1)*7.0d0/8.0d0*dy+D_e*dy+D_s*dx
+D_I*9.0d0/3.0d0*dy+D_L*9.0d0/&
3.0d0*dx)/alpha_v
b_W(Ny,1)=1.0d-20 !0.0d0
b_E(Ny,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
b_S(Ny,1)=F_s_v(Ny,1)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
b_N(Ny,1)=1.0d-20 !0.0d0
S_v_fake(Ny,1)=(v_I*dy)*F_w_v(Ny,1)+(-3.0d0/8.0d0*v_fake(Ny,2)*dy
+2.0d0/8.0d0*v_I*dy&
)*F_e_v(Ny,1)+(3.0d0/8.0d0*v_fake(Ny,
1)*dx-2.0d0/8.0d0*v_fake(Ny-1,1)*dx-1.0d0/8.0d0&
*v_fake(Ny-2,1)*dx)*F_s_v(Ny,1)+(-v_L*dx)*F_n_v(Ny,1)+
(p_fake(Ny-1,1)-p_fake(Ny,1))*&
dx+D_I*8.0d0/3.0d0*v_I*dy+D_L*8.0d0/3.0d0*v_L*dx+((1.0d0-
alpha_v)*b_P(Ny,1))*&
v_old(Ny,1)


! b_P(Ny,1)=1.0d0
! b_W(Ny,1)=1.0d-20
! b_E(Ny,1)=1.0d-20
! b_S(Ny,1)=1.0d-20
! b_N(Ny,1)=1.0d-20
! S_v_fake(Ny,1)=1.0d-30


!coeficientes de p
c_W(Ny,1)=dy/a_P(Ny,1)*dy
c_E(Ny,1)=dy/a_E(Ny,1)*alpha_u*dy
c_S(Ny,1)=dx/b_P(Ny,1)*dx
c_N(Ny,1)=dx/b_N(Ny,1)*alpha_v*dx
c_P(Ny,1)=c_W(Ny,1)+c_E(Ny,1)+c_S(Ny,1)+c_N(Ny,1)
S_p_fake(Ny,1)=u_fake(Ny,1)*dy-u_fake(Ny,2)*dy !Faltam termos de v.


!fronteira SE

!coeficientes de u
a_P(1,Nx)=(F_n_u(1,Nx)*7.0d0/8.0d0*dx+D_w*dy+D_n*dx
+D_F*9.0d0/3.0d0*dy+D_G*9.0d0/&
3.0d0*dx)/alpha_u
a_W(1,Nx)=F_w_u(1,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
a_E(1,Nx)=1.0d-20 !0.0d0
a_S(1,Nx)=1.0d-20 !0.0d0
a_N(1,Nx)=D_n*dx+D_G*1.0d0/3.0d0*dx

S_fake(1,Nx)=(3.0d0/8.0d0*u_fake(1,Nx)*dy-2.0d0/8.0d0*u_fake(1,Nx-1)*dy-1..0d0/8.0d0*&
u_fake(1,Nx-2)*dy)*F_w_u(1,Nx)+(-u_F(1)*dy)*F_e_u(1,Nx)+
(u_G*dx)*F_s_u(1,Nx)+(-3.0d0&
/8.0d0*u_fake(2,Nx)*dx+2.0d0/8.0d0*u_G*dx)*F_n_u(1,Nx)+
(p_fake(1,Nx-1)-p_fake(1,Nx))&
*dy+D_F*8.0d0/3.0d0*u_F(1)*dy+D_G*8.0d0/3.0d0*u_G*dx+((1.0d0-
alpha_u)*a_P(1,Nx))*&
u_old(1,Nx)

! a_P(1,Nx)=1.0d0
! a_W(1,Nx)=1.0d-20
! a_E(1,Nx)=1.0d-20
! a_S(1,Nx)=1.0d-20
! a_N(1,Nx)=1.0d-20
! S_fake(1,Nx)=1.0d-30








!coeficientes de v
b_P(1,Nx)=(F_n_v(1,Nx)*7.0d0/8.0d0*dx+D_w*dy+D_n*dx
+D_J*9.0d0/3.0d0*dy+D_K*9.0d0/&
3.0d0*dx)/alpha_v
b_W(1,Nx)=F_w_v(1,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
b_E(1,Nx)=1.0d-20 !0.0d0
b_S(1,Nx)=1.0d-20 !0.0d0
b_N(1,Nx)=D_n*dx+D_K*1.0d0/3.0d0*dx

S_v_fake(1,Nx)=(3.0d0/8.0d0*v_fake(1,Nx)*dy-2.0d0/8.0d0*v_fake(1,Nx-1)*dy-1.0d0/
&
8.0d0*v_fake(1,Nx-2)*dy)*F_w_v(1,Nx)+(-v_J*dy)*F_e_v(1,Nx)+
(v_K*dx)*F_s_v(1,Nx)+(&
-3.0d0/8.0d0*v_fake(2,Nx)*dx+2.0d0/8.0d0*v_K*dx)*F_n_v(1,Nx)
+D_J*8.0d0/3.0d0*v_J*dy+&
D_K*8.0d0/3.0d0*v_K*dx+((1.0d0-alpha_v)*b_P(1,Nx))*v_old(1,Nx) !
Falta o termo de pressão.


! b_P(1,Nx)=1.0d0
! b_W(1,Nx)=1.0d-20
! b_E(1,Nx)=1.0d-20
! b_S(1,Nx)=1.0d-20
! b_N(1,Nx)=1.0d-20
! S_v_fake(1,Nx)=1.0d-30




!coeficientes de p
c_W(1,Nx)=dy/a_P(1,Nx)*dy
c_E(1,Nx)=dy/a_E(1,Nx)*alpha_u*dy
c_S(1,Nx)=dx/b_P(1,Nx)*dx
c_N(1,Nx)=dx/b_N(1,Nx)*alpha_v*dx
c_P(1,Nx)=c_W(1,Nx)+c_E(1,Nx)+c_S(1,Nx)+c_N(1,Nx)
S_p_fake(1,Nx)=v_fake(1,Nx)*dx-v_fake(2,Nx)*dx !Faltam termos de u.


!fronteira NE

!coeficientes de u
a_P(Ny,Nx)=(D_w*dy+D_s*dx+D_F*9.0d0/3.0d0*dy+D_H*9.0d0/3.0d0*dx)/
alpha_u
a_W(Ny,Nx)=F_w_u(Ny,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
a_E(Ny,Nx)=1.0d-20 !0.0d0
a_S(Ny,Nx)=F_s_u(Ny,Nx)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
a_N(Ny,Nx)=1.0d-20 !0.0d0

S_fake(Ny,Nx)=(3.0d0/8.0d0*u_fake(Ny,Nx)*dy-2.0d0/8.0d0*u_fake(Ny,Nx-1)*dy-1.0d0/
&
8.0d0*u_fake(Ny,Nx-2)*dy)*F_w_u(Ny,Nx)+(-u_F(Ny)*dy)*F_e_u(Ny,Nx)
+(3.0d0/8.0d0*&

u_fake(Ny,Nx)*dx-2.0d0/8.0d0*u_fake(Ny-1,Nx)*dx-1.0d0/8.0d0*u_fake(Ny-2,Nx)*dx)*&
F_s_u(Ny,Nx)+(-u_H*dx)*F_n_u(Ny,Nx)+(p_fake(Ny,Nx-1)-
p_fake(Ny,Nx))*dy+D_F*8.0d0/&
3.0d0*u_F(Ny)*dy+D_H*8.0d0/3.0d0*u_H*dx+((1.0d0-
alpha_u)*a_P(Ny,Nx))*u_old(Ny,Nx)

! a_P(Ny,Nx)=1.0d0
! a_W(Ny,Nx)=1.0d-20
! a_E(Ny,Nx)=1.0d-20
! a_S(Ny,Nx)=1.0d-20
! a_N(Ny,Nx)=1.0d-20
! S_fake(Ny,Nx)=1.0d-30




!coeficientes de v
b_P(Ny,Nx)=(D_w*dy+D_s*dx+D_J*9.0d0/3.0d0*dy+D_L*9.0d0/3.0d0*dx)/
alpha_v
b_W(Ny,Nx)=F_w_v(Ny,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
b_E(Ny,Nx)=1.0d-20 !0.0d0
b_S(Ny,Nx)=F_s_v(Ny,Nx)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
b_N(Ny,Nx)=1.0d-20 !0.0d0

S_v_fake(Ny,Nx)=(3.0d0/8.0d0*v_fake(Ny,Nx)*dy-2.0d0/8.0d0*v_fake(Ny,Nx-1)*dy-1.0d0/
&
8.0d0*v_fake(Ny,Nx-2)*dy)*F_w_v(Ny,Nx)+(-v_J*dy)*F_e_v(Ny,Nx)
+(3.0d0/8.0d0*&

v_fake(Ny,Nx)*dx-2.0d0/8.0d0*v_fake(Ny-1,Nx)*dx-1.0d0/8.0d0*v_fake(Ny-2,Nx)*dx)*&
F_s_v(Ny,Nx)+(-v_L*dx)*F_n_v(Ny,Nx)+(p_fake(Ny-1,Nx)-
p_fake(Ny,Nx))*dx+D_J*8.0d0/&
3.0d0*v_J*dy+D_L*8.0d0/3.0d0*v_L*dx+((1.0d0-
alpha_v)*b_P(Ny,Nx))*v_old(Ny,Nx)

! b_P(Ny,Nx)=1.0d0
! b_W(Ny,Nx)=1.0d-20
! b_E(Ny,Nx)=1.0d-20
! b_S(Ny,Nx)=1.0d-20
! b_N(Ny,Nx)=1.0d-20
! S_v_fake(Ny,Nx)=1.0d-30





!coeficientes de p
c_W(Ny,Nx)=dy/a_P(Ny,Nx)*dy
c_E(Ny,Nx)=dy/a_E(Ny,Nx)*alpha_u*dy
c_S(Ny,Nx)=dx/b_P(Ny,Nx)*dx
c_N(Ny,Nx)=dx/b_N(Ny,Nx)*alpha_v*dx
c_P(Ny,Nx)=c_W(Ny,Nx)+c_E(Ny,Nx)+c_S(Ny,Nx)+c_N(Ny,Nx)
S_p_fake(Ny,Nx)=0.0d0 !Faltam termos de u e de v.


!Fonte

k=1

do i=1,Ny
do j=1,Nx
S(k)=S_fake(i,j)
S_v(k)=S_v_fake(i,j)
S_p(k)=S_p_fake(i,j)
k=k+1
end do
end do


!Matriz dos Coeficientes

k=1

do i=1,Ny
do j=1,Nx
a(k,k)=a_P(i,j)
b(k,k)=b_P(i,j)
c(k,k)=c_P(i,j)
k=k+1
end do
end do


k=2

do i=1,Ny
if (i.eq.1) then
do j=2,Nx
a(k,k-1)=-a_W(i,j)
b(k,k-1)=-b_W(i,j)
c(k,k-1)=-c_W(i,j)
k=k+1
end do
else
do j=1,Nx
a(k,k-1)=-a_W(i,j)
b(k,k-1)=-b_W(i,j)
c(k,k-1)=-c_W(i,j)
k=k+1
end do
end if
end do

k=2

do i=1,Ny
if (i.ne.Ny) then
do j=1,Nx
a(k-1,k)=-a_E(i,j)
b(k-1,k)=-b_E(i,j)
c(k-1,k)=-c_E(i,j)
k=k+1
end do
else
do j=1,Nx-1
a(k-1,k)=-a_E(i,j)
b(k-1,k)=-b_E(i,j)
c(k-1,k)=-c_E(i,j)
k=k+1
end do
end if
end do

k=Nx+1

do i=2,Ny
do j=1,Nx
a(k,k-Nx)=-a_S(i,j)
b(k,k-Nx)=-b_S(i,j)
c(k,k-Nx)=-c_S(i,j)
k=k+1
end do
end do

k=1

do i=1,Ny-1
do j=1,Nx
a(k,k+Nx)=-a_N(i,j)
b(k,k+Nx)=-b_N(i,j)
c(k,k+Nx)=-c_N(i,j)
k=k+1
end do
end do



k=0
l=0
m=0

do i=1,Ny*Nx
do j=1,Ny*Nx
if (dabs(a(i,j)).ne.0.0d0) then
k=k+1
end if
if (dabs(b(i,j)).ne.0.0d0) then
l=l+1
end if
if (dabs(c(i,j)).ne.0.0d0) then
m=m+1
end if
end do
end do

allocate(AA(k),irow_u(k),jcol_u(k))
allocate(BB(l),irow_v(l),jcol_v(l))
allocate(CC(m),irow_p(m),jcol_p(m))

k=1
l=1
m=1

do i=1,Ny*Nx
do j=1,Ny*Nx
if (dabs(a(i,j)).ne.0.0d0) then
AA(k)=a(i,j)
irow_u(k)=i
jcol_u(k)=j
k=k+1
end if
if (dabs(b(i,j)).ne.0.0d0) then
BB(l)=b(i,j)
irow_v(l)=i
jcol_v(l)=j
l=l+1
end if
if (dabs(c(i,j)).ne.0.0d0) then
CC(m)=c(i,j)
irow_p(m)=i
jcol_p(m)=j
m=m+1
end if
end do
end do

call DL4LXG(iparam,rparam)

iparam(5)=1000000


!call DLINRG(Ny*Nx,a,Ny*Nx,a_inv,Ny*Nx)
!CALL DLSARG (Ny*Nx,a,Ny*Nx,S,ipath,u)
call DLSLXG(Ny*Nx,k-1,AA,irow_u,jcol_u,S,ipath,iparam,rparam,u)
call DLSLXG(Ny*Nx,l-1,BB,irow_v,jcol_v,S_v,ipath,iparam,rparam,v)
call DLSLXG(Ny*Nx,m-1,CC,irow_p,jcol_p,S_p,ipath,iparam,rparam,p)

! k=1


! do i=1,Ny
! do j=1,Nx
! p_fake(i,j)=p(k)
! if (j.ne.1) then
! u_fake(i,j)=u(k)+dy/(a_P(i,j)*alpha_u)*(p_fake(i,j-1)-
p_fake(i,j))
! end if
! if (i.ne.1) then
! v_fake(i,j)=v(k)+dx/(b_P(i,j)*alpha_v)*(p_fake(i-1,j)-
p_fake(i,j))
! end if
! p_fake(i,j)=p_old(i,j)+alpha_p*p(k)
! k=k+1
! end do
! end do


p_old=p_fake
u_old=u_fake
v_old=v_fake


k=1

do i=1,Ny
do j=1,Nx
p_fake(i,j)=p(k)
!p_fake(i,j)=p_old(i,j)+alpha_p*p(k)
!p_fake(i,j)=alpha_p*(k)+(1.0d0-alpha_p)*p_old(i,j)
!if (i.ne.1) then
! if (i.ne.Ny) then
! if (j.ne.1) then
! if (j.ne.Nx) then
!u_fake(i,j)=u(k)+dy/(a_P(i,j)*alpha_u)*(p_fake(i,j-1)-
p_fake(i,j))
!v_fake(i,j)=v(k)+dx/(b_P(i,j)*alpha_v)*(p_fake(i-1,j)-
p_fake(i,j))
! u_fake(i,j)=u(k)+dy/a_P(i,j)*(p_fake(i,j-1)-p_fake(i,j))
! v_fake(i,j)=v(k)+dx/b_P(i,j)*(p_fake(i-1,j)-p_fake(i,j))
! u_fake(i,j)=alpha_u*u(k)+(1.0d0-alpha_u)*u_old(i,j)
! v_fake(i,j)=alpha_v*v(k)+(1.0d0-alpha_v)*v_old(i,j)
! end if
! end if
! end if
!end if
if (j.ne.1) then
u_fake(i,j)=u(k)+dy/(a_P(i,j)*alpha_u)*(p_fake(i,j-1)-p_fake(i,j))
!u_fake(i,j)=u(k)+dy/a_P(i,j)*(p_fake(i,j-1)-p_fake(i,j))
!u_fake(i,j)=alpha_u*u(k)
!u_fake(i,j)=alpha_u*u(k)+(1.0d0-alpha_u)*u_old(i,j)
! write(*,*) u_fake(i,j)
end if
if (i.ne.1) then
v_fake(i,j)=v(k)+dx/(b_P(i,j)*alpha_v)*(p_fake(i-1,j)-p_fake(i,j))
!v_fake(i,j)=v(k)+dx/b_P(i,j)*(p_fake(i-1,j)-p_fake(i,j))
!v_fake(i,j)=alpha_v*v(k)
!v_fake(i,j)=alpha_v*v(k)+(1.0d0-alpha_v)*v_old(i,j)
end if
p_fake(i,j)=p_old(i,j)+alpha_p*p(k)
k=k+1
end do
end do




deallocate(AA,irow_u,jcol_u,BB,irow_v,jcol_v,CC,irow_p,jcol_p)

!u=matmul(a_inv,S)
!theta=matmul(a_inv,S)



! end do

open(1,file='u.dat')
write(1,11) Ny,Nx
11 format('variables=X,Y,u',/,'zone t="time','" i=',i4,' j=',i4,'
f=point')

k=1

do i=1,Ny
do j=1,Nx
write(1,*) Y(i),X(j),u(k)
k=k+1
end do
end do

open(2,file='v.dat')
write(2,12) Ny,Nx
12 format('variables=X,Y,v',/,'zone t="time','" i=',i4,' j=',i4,'
f=point')

k=1

do i=1,Ny
do j=1,Nx
write(2,*) Y(i),X(j),v(k)
k=k+1
end do
end do

open(3,file='p.dat')
write(3,13) Ny,Nx
13 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
f=point')

k=1

do i=1,Ny
do j=1,Nx
write(3,*) Y(i),X(j),p(k)
k=k+1
end do
end do

open(4,file='S.dat')
write(4,14) Ny,Nx
14 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
f=point')

k=1

do i=1,Ny
do j=1,Nx
write(4,*) S_fake(i,j),k,S(k)
k=k+1
end do
end do

open(5,file='a.dat')
!write(5,15) Ny,Nx
!15 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
f=point')


do i=1,Ny*Nx
do j=1,Ny*Nx
write(5,*) i,j,a(i,j)
end do
end do

open(6,file='aa.dat')
!write(6,16) Ny,Nx
!16 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
f=point')


do i=1,Ny
do j=1,Nx
write(6,*) i,j,a_P(i,j)
end do
end do



end program EXERCICIO
From: Gordon Sande on
On 2010-07-05 16:37:27 -0300, Felipe Abr�o <felipeabrao(a)gmail.com> said:

> Please, could somebody explain me what is happening. I simply
> allocated a matrix called "a", but figured out that the program is
> inserting some strange values to some of its elements (like
> 1.057229231618489E-307 or 6.941570408926931E+228), even before I set
> any value to it! Shouldn't all the values be zero before any desired
> value setting?

No!

Before you set the value it has NO value!

The "history" is that back when computers were $1000 per hour and programmers
were paid $3.50 per hour it was more efficient to just use whatever value
was left over from the previous program that was executing. It is called
an undefined value. Fortran starts things with an undefined value. This is
not as silly as it sounds as you might have a subroutine that gets used twice.
What values should be used the second time? If you had been using a
good debugging
compiler, and had set the suitable options to have it check, then you
would have
gotten a diagnostic about an undefined value.

Since you are using Windows you might try the Silverfrost Fortran 90 (it
used to be Salford who developed it for student use) which is free for
personal use. Be sure to set the /undef option!

> Below is my code, but you don't have to read it all.
> You can check the problem just after line 65, where I wrote a command
> in order to read the values of matrix "a". The same thing is happening
> with matrixs "b" and "c". In advance I thank you all!
>
> PS: I use Fortran 90 in Windows XP.
>
>
> ! EXERCICIO.f90
> !
> ! FUNCTIONS:
> ! EXERCICIO - Entry point of console application.
> !
>
> !
> ***************************************************************************
> *
> !
> ! PROGRAM: EXERCICIO
> !
> ! PURPOSE: Entry point for the console application.
> !
> !
> ***************************************************************************
> *
>
> program EXERCICIO
>
> use IMSL
>
> implicit none
>
> ! Variables
>
> integer(4)::i,j,k,l,m,Nx,Ny
> integer(4),parameter::ipath=1
> integer(4),dimension(6)::iparam
>
> integer(4),allocatable,dimension(:)::irow_u,jcol_u,irow_v,jcol_v,irow_p,jco
> l_p
> real(8)::dx,dy,D_w,D_e,D_s,D_n,D_EE,D_F,D_G,D_H,D_I,D_J,D_K,D_L
> real(8),parameter::Re=1.0d2,gama=1.0d0/
> Re,u_E=1.0d0,u_G=0.0d0,u_H=0.0d0,v_I=0.0d0,v_J=0.0d0,&
>
> v_K=0.0d0,v_L=0.0d0,error=1.0d-6,alpha_u=5.0d-1,alpha_v=5.0d-1,al
> pha_p=5.0d-1
> !real(8),dimension(N*N)::theta_old
> real(8),dimension(5)::rparam
>
> real(8),allocatable,dimension(:)::u,v,p,theta,S,S_v,S_p,X,Y,u_F,AA,BB,CC
> !real(8),dimension(N,N)::theta_fake
>
> real(8),allocatable,dimension(:,:)::a,a_inv,b,c,p_fake,p_old,u_fake,u_old,v
> _fake,v_old,a_P,&
>
> a_W,a_E,a_S,a_N,b_P,b_W,b_E,b_S,b_N,c_P,c_W,c_E,c_S,c_N,F_w_u,F_e_u,F_s_u,F
> _n_u,S_fake,&
> S_v_fake,S_p_fake,F_w_v,F_e_v,F_s_v,F_n_v
>
> ! Body of EXERCICIO_1
>
> iparam(1)=0
>
> Nx=5
> Ny=5
>
>
>
>
> allocate(a(Ny*Nx,Ny*Nx),a_inv(Ny*Nx,Ny*Nx),b(Ny*Nx,Ny*Nx),c(Ny*Nx,Ny*Nx))
>
> allocate(a_P(Ny,Nx),a_W(Ny,Nx),a_E(Ny,Nx),a_S(Ny,Nx),a_N(Ny,Nx),b_P(Ny,Nx),
> b_W(Ny,Nx),&
>
> b_E(Ny,Nx),b_S(Ny,Nx),b_N(Ny,Nx),c_P(Ny,Nx),c_W(Ny,Nx),c_E(Ny,Nx),c_S(Ny,Nx
> ),c_N(Ny,Nx),&
>
> p_fake(Ny,Nx),p_old(Ny,Nx),u_fake(Ny,Nx),u_old(Ny,Nx),v_fake(Ny,Nx),v_old(N
> y,Nx),F_w_u(Ny,Nx),&
>
> F_e_u(Ny,Nx),F_s_u(Ny,Nx),F_n_u(Ny,Nx),S_fake(Ny,Nx),S_v_fake(Ny,Nx),S_p_fa
> ke(Ny,Nx),&
> F_w_v(Ny,Nx),F_e_v(Ny,Nx),F_s_v(Ny,Nx),F_n_v(Ny,Nx))
>
> allocate(u(Ny*Nx),v(Ny*Nx),p(Ny*Nx),theta(Ny*Nx),S(Ny*Nx),S_v(Ny*Nx),S_p(Ny
> *Nx))
> allocate(Y(Ny),u_F(Ny))
> allocate(X(Nx))
>
>
> open(7,file='pre.dat')
>
> do i=1,Ny*Nx
> do j=1,Ny*Nx
> write(7,*) i,j,a(i,j)
> end do
> end do
>
>
>
> dx=1.0d0/dble(Nx)
> dy=1.0d0/dble(Ny)
> D_w=gama/dx
> D_e=gama/dx
> D_s=gama/dy
> D_n=gama/dy
> D_EE=2.0d0*D_w
> D_F=2.0d0*D_e
> D_G=2.0d0*D_s
> D_H=2.0d0*D_n
> D_I=D_EE
> D_J=D_F
> D_K=D_G
> D_L=D_H
>
> p_fake=5.0d-1
> u_fake=5.0d-1
> v_fake=5.0d-1
>
> ! do i=1,Ny
> ! do j=1,Nx
> ! u_fake(i,j)=dble(i+j)
> ! v_fake(i,j)=dble(i+j)
> ! p_fake(i,j)=dble(i+j)
> ! end do
> ! end do
>
>
>
> theta=5.0d-1
> !theta_old=6.0d-1
>
> do i=1,Nx
> X(i)=1.0d0/dble((Nx-1))*dble((i-1))
> end do
>
> do i=1,Ny
> Y(i)=1.0d0/dble((Ny-1))*dble((i-1))
> u_F(i)=4.0d0*(Y(i)-Y(i)**2.0d0)
> end do
>
> p_old=1.0d0
> u_old=1.0d0
> v_old=1.0d0
>
> ! do while(maxval(dabs(p_fake-p_old)).gt.error)
>
> !do while((maxval(dabs(p_fake-p_old)).gt.error).or.(maxval(dabs(u_fake-
> u_old)).gt.error).or.&
> !(maxval(dabs(v_fake-v_old)).gt.error))
>
> !do while((maxval(dabs(p_fake(2:Ny-1,2:Nx-1)-
> p_old(2:Ny-1,2:Nx-1))).gt.error).or.(maxval(dabs(&
> !u_fake(2:Ny-1,2:Nx-1)-u_old(2:Ny-1,2:Nx-1))).gt.error).or.
> (maxval(dabs(v_fake(2:Ny-1,2:Nx-1)-&
> !v_old(2:Ny-1,2:Nx-1))).gt.error))
>
>
>
> ! write(*,*) 'diferen�a de press�o =',maxval(dabs(p_fake-p_old))
> ! write(*,*) 'diferen�a de velocidade u =',maxval(dabs(u_fake-u_old))
> ! write(*,*) 'diferen�a de velocidade v =',maxval(dabs(v_fake-v_old))
>
>
> ! p_old=p_fake
> ! u_old=u_fake
> ! v_old=v_fake
>
>
> do i=1,Ny
> do j=1,Nx
> if (j.eq.1) then
> F_w_u(2:Ny-1,j)=1.0d0
> else
> F_w_u(i,j)=1.0d0/2.0d0*(u_fake(i,j)+u_fake(i,j-1))
> end if
> if (i.eq.1) then
> F_w_v(i,j)=0.0d0 !Tenho d�vida.
> else
> F_w_v(i,j)=1.0d0/2.0d0*(u_fake(i,j)+u_fake(i-1,j))
> end if
> if (j.eq.Nx) then
> F_e_u(i,j)=4.0d0*(Y(i)-Y(i)**2.0d0)
> else
> F_e_u(i,j)=1.0d0/2.0d0*(u_fake(i,j+1)+u_fake(i,j))
> end if
> if (i.eq.1) then
> F_e_v(i,j)=0.0d0 !Tenho d�vida.
> else if (j.eq.Nx) then
> F_e_v(i,j)=4.0d0*(Y(i)-Y(i)**2.0d0) !Tenho d�vida.
> else
> F_e_v(i,j)=1.0d0/2.0d0*(u_fake(i,j+1)+u_fake(i-1,j+1))
> end if
> if (i.eq.1) then
> F_s_u(i,j)=0.0d0
> else if (j.eq.1) then
> F_s_u(i,j)=0.0d0
> else
> F_s_u(i,j)=1.0d0/2.0d0*(v_fake(i,j)+v_fake(i,j-1))
> end if
> if (i.eq.1) then
> F_s_v(i,j)=0.0d0
> else
> F_s_v(i,j)=1.0d0/2.0d0*(v_fake(i-1,j)+v_fake(i,j))
> end if
> if (i.eq.Ny) then
> F_n_u(i,j)=0.0d0
> else if (j.eq.1) then
> F_n_u(i,j)=0.0d0
> else
> F_n_u(i,j)=1.0d0/2.0d0*(v_fake(i+1,j)+v_fake(i+1,j-1))
> end if
> if (i.eq.Ny) then
> F_n_v(i,j)=0.0d0
> else
> F_n_v(i,j)=1.0d0/2.0d0*(v_fake(i,j)+v_fake(i+1,j))
> end if
> end do
> end do
>
>
>
>
>
>
>
> !interior
>
> do i=3,Ny-1
> do j=3,Nx-1
> !coeficientes de u
> a_P(i,j)=(F_e_u(i,j)*dy+F_n_u(i,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_u
> a_W(i,j)=F_w_u(i,j)*dy+F_e_u(i,j)*1.0d0/8.0d0*dy+D_w*dy
> a_E(i,j)=D_e*dy
> a_S(i,j)=F_s_u(i,j)*dx+F_n_u(i,j)*1.0d0/8.0d0*dx+D_s*dx
> a_N(i,j)=D_n*dx
>
> S_fake(i,j)=(3.0d0/8.0d0*u_fake(i,j)*dy-2.0d0/8.0d0*u_fake(i,j-1)*dy-1.0d
> 0/8.0d0&
> *u_fake(i,j-2)*dy)*F_w_u(i,j)+(-3.0d0/8.0d0*u_fake(i,j+1)*dy
> +2.0d0/8.0d0*&
> u_fake(i,j)*dy)*F_e_u(i,j)
> +(3.0d0/8.0d0*u_fake(i,j)*dx-2.0d0/8.0d0*u_fake(i-1,j)&
> *dx-1.0d0/8.0d0*u_fake(i-2,j)*dx)*F_s_u(i,j)+(-3.0d0/8.0d0*u_fake(i
> +1,j)*dx+&
> 2.0d0/8.0d0*u_fake(i,j)*dx)*F_n_u(i,j)+(p_fake(i,j-1)-
> p_fake(i,j))*dy+((1.0d0-&
> alpha_u)*a_P(i,j))*u_old(i,j)
>
> !coeficientes de v
> b_P(i,j)=(F_e_v(i,j)*dy+F_n_v(i,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_v
> b_W(i,j)=F_w_v(i,j)*dy+F_e_v(i,j)*1.0d0/8.0d0*dy+D_w*dy
> b_E(i,j)=D_e*dy
> b_S(i,j)=F_s_v(i,j)*dx+F_n_v(i,j)*1.0d0/8.0d0*dx+D_s*dx
> b_N(i,j)=D_n*dx
>
> S_v_fake(i,j)=(3.0d0/8.0d0*v_fake(i,j)*dy-2.0d0/8.0d0*v_fake(i,j-1)*dy-1.
> 0d0/
> &
> 8.0d0*v_fake(i,j-2)*dy)*F_w_v(i,j)+(-3.0d0/8.0d0*v_fake(i,j+1)*dy
> +2.0d0/8.0d0*&
> v_fake(i,j)*dy)*F_e_v(i,j)
> +(3.0d0/8.0d0*v_fake(i,j)*dx-2.0d0/8.0d0*v_fake(i-1,j)&
> *dx-1.0d0/8.0d0*v_fake(i-2,j)*dx)*F_s_v(i,j)+(-3.0d0/8.0d0*v_fake(i
> +1,j)*dx+&
> 2.0d0/8.0d0*v_fake(i,j)*dx)*F_n_v(i,j)+(p_fake(i-1,j)-
> p_fake(i,j))*dx+((1.0d0-&
> alpha_v)*b_P(i,j))*v_old(i,j)
>
> !coeficientes de p
> c_W(i,j)=dy/a_P(i,j)*dy
> c_E(i,j)=dy/a_E(i,j)*alpha_u*dy
> c_S(i,j)=dx/b_P(i,j)*dx
> c_N(i,j)=dx/b_N(i,j)*alpha_v*dx
> c_P(i,j)=c_W(i,j)+c_E(i,j)+c_S(i,j)+c_N(i,j)
> S_p_fake(i,j)=u_fake(i,j)*dy-u_fake(i,j+1)*dy+v_fake(i,j)*dx-
> v_fake(i+1,j)*dx
>
> end do
> end do
>
> !coluna 2 - com exce��o da linha 2
>
> do i=3,Ny-1
> !coeficientes de u
> a_P(i,2)=(F_e_u(i,2)*dy+F_n_u(i,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_u
> a_W(i,2)=F_w_u(i,2)*7.0d0/8.0d0*dy+F_e_u(i,2)*1.0d0/8.0d0*dy+D_w*dy
> a_E(i,2)=D_e*dy
> a_S(i,2)=F_s_u(i,2)*dx+F_n_u(i,2)*1.0d0/8.0d0*dx+D_s*dx
> a_N(i,2)=D_n*dx
> S_fake(i,2)=(3.0d0/8.0d0*u_fake(i,2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(i,
> 2)+(-3.0d0/8.0d0&
> *u_fake(i,3)*dy+2.0d0/8.0d0*u_fake(i,2)*dy)*F_e_u(i,
> 2)+(3.0d0/8.0d0*u_fake(i,2)*dx-&
> 2.0d0/8.0d0*u_fake(i-1,2)*dx-1.0d0/8.0d0*u_fake(i-2,2)*dx)*F_s_u(i,
> 2)+(-3.0d0/8.0d0*&
> u_fake(i+1,2)*dx+2.0d0/8.0d0*u_fake(i,2)*dx)*F_n_u(i,2)+(p_fake(i,
> 1)-p_fake(i,2))*dy&
> +((1.0d0-alpha_u)*a_P(i,2))*u_old(i,2)
>
> !coeficientes de v
> b_P(i,2)=(F_e_v(i,2)*dy+F_n_v(i,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_v
> b_W(i,2)=F_w_v(i,2)*7.0d0/8.0d0*dy+F_e_v(i,2)*1.0d0/8.0d0*dy+D_w*dy
> b_E(i,2)=D_e*dy
> b_S(i,2)=F_s_v(i,2)*dx+F_n_v(i,2)*1.0d0/8.0d0*dx+D_s*dx
> b_N(i,2)=D_n*dx
> S_v_fake(i,2)=(3.0d0/8.0d0*v_fake(i,
> 2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(i,2)+(-3.0d0/&
> 8.0d0*v_fake(i,3)*dy+2.0d0/8.0d0*v_fake(i,2)*dy)*F_e_v(i,
> 2)+(3.0d0/8.0d0*v_fake(i,2)&
>
> *dx-2.0d0/8.0d0*v_fake(i-1,2)*dx-1.0d0/8.0d0*v_fake(i-2,2)*dx)*F_s_v(i,
> 2)+(-3.0d0/&
> 8.0d0*v_fake(i+1,2)*dx+2.0d0/8.0d0*v_fake(i,2)*dx)*F_n_v(i,2)+
> (p_fake(i-1,2)-&
> p_fake(i,2))*dx+((1.0d0-alpha_v)*b_P(i,2))*v_old(i,2)
>
> !coeficientes de p
> c_W(i,2)=dy/a_P(i,2)*dy
> c_E(i,2)=dy/a_E(i,2)*alpha_u*dy
> c_S(i,2)=dx/b_P(i,2)*dx
> c_N(i,2)=dx/b_N(i,2)*alpha_v*dx
> c_P(i,2)=c_W(i,2)+c_E(i,2)+c_S(i,2)+c_N(i,2)
> S_p_fake(i,2)=u_fake(i,2)*dy-u_fake(i,3)*dy+v_fake(i,2)*dx-v_fake(i
> +1,2)*dx
>
> end do
>
> !linha 2 - com exce��o da coluna 2
>
> do j=3,Nx-1
> !coeficientes de u
> a_P(2,j)=(F_e_u(2,j)*dy+F_n_u(2,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_u
> a_W(2,j)=F_w_u(2,j)*dy+F_e_u(2,j)*1.0d0/8.0d0*dy+D_w*dy
> a_E(2,j)=D_e*dy
> a_S(2,j)=F_s_u(2,j)*7.0d0/8.0d0*dx+F_n_u(2,j)*1.0d0/8.0d0*dx+D_s*dx
> a_N(2,j)=D_n*dx
>
> S_fake(2,j)=(3.0d0/8.0d0*u_fake(2,j)*dy-2.0d0/8.0d0*u_fake(2,j-1)*dy-1.0d
> 0/8.0d0*&
> u_fake(2,j-2)*dy)*F_w_u(2,j)+(-3.0d0/8.0d0*u_fake(2,j+1)*dy
> +2.0d0/8.0d0*u_fake(2,j)*&
> dy)*F_e_u(2,j)
> +(3.0d0/8.0d0*u_fake(2,j)*dx-2.0d0/8.0d0*u_G*dx)*F_s_u(2,j)+(-3.0d0/&
> 8.0d0*u_fake(3,j)*dx+2.0d0/8.0d0*u_fake(2,j)*dx)*F_n_u(2,j)+
> (p_fake(2,j-1)-&
> p_fake(2,j))*dy+((1.0d0-alpha_u)*a_P(2,j))*u_old(2,j)
>
> !coeficientes de v
> b_P(2,j)=(F_e_v(2,j)*dy+F_n_v(2,j)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_v
> b_W(2,j)=F_w_v(2,j)*dy+F_e_v(2,j)*1.0d0/8.0d0*dy+D_w*dy
> b_E(2,j)=D_e*dy
> b_S(2,j)=F_s_v(2,j)*7.0d0/8.0d0*dx+F_n_v(2,j)*1.0d0/8.0d0*dx+D_s*dx
> b_N(2,j)=D_n*dx
>
> S_v_fake(2,j)=(3.0d0/8.0d0*v_fake(2,j)*dy-2.0d0/8.0d0*v_fake(2,j-1)*dy-1.
> 0d0/8.0d0*&
> v_fake(2,j-2)*dy)*F_w_v(2,j)+(-3.0d0/8.0d0*v_fake(2,j+1)*dy
> +2.0d0/8.0d0*v_fake(2,j)*&
> dy)*F_e_v(2,j)
> +(3.0d0/8.0d0*v_fake(2,j)*dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,j)+(-3.0d0/&
> 8.0d0*v_fake(3,j)*dx+2.0d0/8.0d0*v_fake(2,j)*dx)*F_n_v(2,j)+
> (p_fake(1,j)-&
> p_fake(2,j))*dx+((1.0d0-alpha_v)*b_P(2,j))*v_old(2,j)
>
> !coeficientes de p
> c_W(2,j)=dy/a_P(2,j)*dy
> c_E(2,j)=dy/a_E(2,j)*alpha_u*dy
> c_S(2,j)=dx/b_P(2,j)*dx
> c_N(2,j)=dx/b_N(2,j)*alpha_v*dx
> c_P(2,j)=c_W(2,j)+c_E(2,j)+c_S(2,j)+c_N(2,j)
> S_p_fake(2,j)=u_fake(2,j)*dy-u_fake(2,j+1)*dy+v_fake(2,j)*dx-
> v_fake(3,j)*dx
>
> end do
>
> !coluna 2 - linha 2
>
> !coeficientes de u
> a_P(2,2)=(F_e_u(2,2)*dy+F_n_u(2,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_u
> a_W(2,2)=F_w_u(2,2)*7.0d0/8.0d0*dy+F_e_u(2,2)*1.0d0/8.0d0*dy
> +D_w*dy
> a_E(2,2)=D_e*dy
> a_S(2,2)=F_s_u(2,2)*7.0d0/8.0d0*dx+F_n_u(2,2)*1.0d0/8.0d0*dx
> +D_s*dx
> a_N(2,2)=D_n*dx
>
> S_fake(2,2)=(3.0d0/8.0d0*u_fake(2,2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(2,2)+
> (-3.0d0/&
> 8.0d0*u_fake(2,3)*dy
> +2.0d0/8.0d0*u_fake(2,2)*dy)*F_e_u(2,2)+(3.0d0/8.0d0*&
> u_fake(2,2)*dx-2.0d0/8.0d0*u_G*dx)*F_s_u(2,2)+
> (-3.0d0/8.0d0*u_fake(3,2)*dx+2.0d0&
> /8.0d0*u_fake(2,2)*dx)*F_n_u(2,2)+(p_fake(2,1)-p_fake(2,2))*dy+
> ((1.0d0-alpha_u)*&
> a_P(2,2))*u_old(2,2)
>
> !coeficientes de v
> b_P(2,2)=(F_e_v(2,2)*dy+F_n_v(2,2)*dx+D_w*dy+D_e*dy+D_s*dx+D_n*dx)/
> alpha_v
> b_W(2,2)=F_w_v(2,2)*7.0d0/8.0d0*dy+F_e_v(2,2)*1.0d0/8.0d0*dy
> +D_w*dy
> b_E(2,2)=D_e*dy
> b_S(2,2)=F_s_v(2,2)*7.0d0/8.0d0*dx+F_n_v(2,2)*1.0d0/8.0d0*dx
> +D_s*dx
> b_N(2,2)=D_n*dx
>
> S_v_fake(2,2)=(3.0d0/8.0d0*v_fake(2,2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(2,2)+
> (-3.0d0&
> /8.0d0*v_fake(2,3)*dy
> +2.0d0/8.0d0*v_fake(2,2)*dy)*F_e_v(2,2)+(3.0d0/8.0d0*&
> v_fake(2,2)*dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,2)+
> (-3.0d0/8.0d0*v_fake(3,2)*dx+2.0d0&
> /8.0d0*v_fake(2,2)*dx)*F_n_v(2,2)+(p_fake(1,2)-p_fake(2,2))*dx+
> ((1.0d0-alpha_v)*&
> b_P(2,2))*v_old(2,2)
>
> !coeficientes de p
> c_W(2,2)=dy/a_P(2,2)*dy
> c_E(2,2)=dy/a_E(2,2)*alpha_u*dy
> c_S(2,2)=dx/b_P(2,2)*dx
> c_N(2,2)=dx/b_N(2,2)*alpha_v*dx
> c_P(2,2)=c_W(2,2)+c_E(2,2)+c_S(2,2)+c_N(2,2)
> S_p_fake(2,2)=u_fake(2,2)*dy-u_fake(2,3)*dy+v_fake(2,2)*dx-
> v_fake(3,2)*dx
>
>
> !fronteira W - com exce��o da linha 2
>
> do i=3,Ny-1
> !coeficientes de u
> a_P(i,1)=(F_e_u(i,1)*7.0d0/8.0d0*dy+F_n_u(i,1)*dx+D_e*dy+D_s*dx
> +D_n*dx+D_EE*9.0d0/&
> 3.0d0*dy)/alpha_u
> a_W(i,1)=1.0d-20 !0.0d0
> a_E(i,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
> a_S(i,1)=F_s_u(i,1)*dx+F_n_u(i,1)*1.0d0/8.0d0*dx+D_s*dx
> a_N(i,1)=D_n*dx
> S_fake(i,1)=(u_E*dy)*F_w_u(i,1)+(-3.0d0/8.0d0*u_fake(i,2)*dy
> +2.0d0/8.0d0*u_E*dy)*&
> F_e_u(i,1)+(3.0d0/8.0d0*u_fake(i,
> 1)*dx-2.0d0/8.0d0*u_fake(i-1,1)*dx-1.0d0/8.0d0*&
> u_fake(i-2,1)*dx)*F_s_u(i,1)+(-3.0d0/8.0d0*u_fake(i+1,1)*dx
> +2.0d0/8.0d0*u_fake(i,1)*&
> dx)*F_n_u(i,1)+D_EE*8.0d0/3.0d0*u_E*dy+((1.0d0-alpha_u)*a_P(i,
> 1))*u_old(i,1) !Falta o termo de press�o.
>
>
> ! a_P(i,1)=1.0d0
> ! a_W(i,1)=1.0d-20
> ! a_E(i,1)=1.0d-20
> ! a_S(i,1)=1.0d-20
> ! a_N(i,1)=1.0d-20
> ! S_fake(i,1)=u_E
>
>
>
>
>
>
> !coeficientes de v
> b_P(i,1)=(F_e_v(i,1)*7.0d0/8.0d0*dy+F_n_v(i,1)*dx+D_e*dy+D_s*dx
> +D_n*dx+D_I*9.0d0/&
> 3.0d0*dy)/alpha_v
> b_W(i,1)=1.0d-20 !0.0d0
> b_E(i,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
> b_S(i,1)=F_s_v(i,1)*dx+F_n_v(i,1)*1.0d0/8.0d0*dx+D_s*dx
> b_N(i,1)=D_n*dx
> S_v_fake(i,1)=(v_I*dy)*F_w_v(i,1)+(-3.0d0/8.0d0*v_fake(i,2)*dy
> +2.0d0/8.0d0*v_I*dy)*&
> F_e_v(i,1)+(3.0d0/8.0d0*v_fake(i,
> 1)*dx-2.0d0/8.0d0*v_fake(i-1,1)*dx-1.0d0/8.0d0*&
> v_fake(i-2,1)*dx)*F_s_v(i,1)+(-3.0d0/8.0d0*v_fake(i+1,1)*dx
> +2.0d0/8.0d0*v_fake(i,1)*&
> dx)*F_n_v(i,1)+(p_fake(i-1,1)-p_fake(i,1))*dx+D_I*8.0d0/3.0d0*v_I*dy
> +((1.0d0-alpha_v&
> )*b_P(i,1))*v_old(i,1)
>
> ! b_P(i,1)=1.0d0
> ! b_W(i,1)=1.0d-20
> ! b_E(i,1)=1.0d-20
> ! b_S(i,1)=1.0d-20
> ! b_N(i,1)=1.0d-20
> ! S_v_fake(i,1)=v_I
>
>
>
>
> !coeficientes de p
>
> c_W(i,1)=dy/a_P(i,1)*dy
> c_E(i,1)=dy/a_E(i,1)*alpha_u*dy
> c_S(i,1)=dx/b_P(i,1)*dx
> c_N(i,1)=dx/b_N(i,1)*alpha_v*dx
> c_P(i,1)=c_W(i,1)+c_E(i,1)+c_S(i,1)+c_N(i,1)
> S_p_fake(i,1)=u_fake(i,1)*dy-u_fake(i,2)*dy+v_fake(i,1)*dx-v_fake(i
> +1,1)*dx
>
> end do
>
> !fronteira W - linha 2
>
> !coeficientes de u
> a_P(2,1)=(F_e_u(2,1)*7.0d0/8.0d0*dy+F_n_u(2,1)*dx+D_e*dy+D_s*dx
> +D_n*dx+D_EE*9.0d0/&
> 3.0d0*dy)/alpha_u
> a_W(2,1)=1.0d-20 !0.0d0
> a_E(2,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
> a_S(2,1)=F_s_u(2,1)*7.0d0/8.0d0*dx+F_n_u(2,1)*1.0d0/8.0d0*dx+D_s*dx
> a_N(2,1)=D_n*dx
> S_fake(2,1)=(u_E*dy)*F_w_u(2,1)+(-3.0d0/8.0d0*u_fake(2,2)*dy
> +2.0d0/8.0d0*u_E*dy)*&
>
> F_e_u(2,1)+(3.0d0/8.0d0*u_fake(2,1)*dx-2.0d0/8.0d0*u_G*dx)*F_s_u(2,1)+
> (-3.0d0/8.0d0*&
> u_fake(3,1)*dx
> +2.0d0/8.0d0*u_fake(2,1)*dx)*F_n_u(2,1)+D_EE*8.0d0/3.0d0*u_E*dy+((&
> 1.0d0-alpha_u)*a_P(2,1))*u_old(2,1) !Falta o termo de press�o.
>
>
> ! a_P(2,1)=1.0d0
> ! a_W(2,1)=1.0d-20
> ! a_E(2,1)=1.0d-20
> ! a_S(2,1)=1.0d-20
> ! a_N(2,1)=1.0d-20
> ! S_fake(2,1)=u_E
>
>
>
>
> !coeficientes de v
> b_P(2,1)=(F_e_v(2,1)*7.0d0/8.0d0*dy+F_n_v(2,1)*dx+D_e*dy+D_s*dx
> +D_n*dx+D_I*9.0d0/&
> 3.0d0*dy)/alpha_v
> b_W(2,1)=1.0d-20 !0.0d0
> b_E(2,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
> b_S(2,1)=F_s_v(2,1)*7.0d0/8.0d0*dx+F_n_v(2,1)*1.0d0/8.0d0*dx+D_s*dx
> b_N(2,1)=D_n*dx
> S_v_fake(2,1)=(v_I*dy)*F_w_v(2,1)+(-3.0d0/8.0d0*v_fake(2,2)*dy
> +2.0d0/8.0d0*v_I*dy)*&
>
> F_e_v(2,1)+(3.0d0/8.0d0*v_fake(2,1)*dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,1)+
> (-3.0d0/8.0d0*&
> v_fake(3,1)*dx+2.0d0/8.0d0*v_fake(2,1)*dx)*F_n_v(2,1)+(p_fake(1,1)-
> p_fake(2,1))*dx+&
> D_I*8.0d0/3.0d0*v_I*dy+((1.0d0-alpha_v)*b_P(2,1))*v_old(2,1)
>
> ! b_P(2,1)=1.0d0
> ! b_W(2,1)=1.0d-20
> ! b_E(2,1)=1.0d-20
> ! b_S(2,1)=1.0d-20
> ! b_N(2,1)=1.0d-20
> ! S_v_fake(2,1)=v_I
>
>
>
>
> !coeficientes de p
> c_W(2,1)=dy/a_P(2,1)*dy
> c_E(2,1)=dy/a_E(2,1)*alpha_u*dy
> c_S(2,1)=dx/b_P(2,1)*dx
> c_N(2,1)=dx/b_N(2,1)*alpha_v*dx
> c_P(2,1)=c_W(2,1)+c_E(2,1)+c_S(2,1)+c_N(2,1)
> S_p_fake(2,1)=u_fake(2,1)*dy-u_fake(2,2)*dy+v_fake(2,1)*dx-
> v_fake(3,1)*dx
>
>
>
> !fronteira E - com exce��o da linha 2
>
> do i=3,Ny-1
> !coeficientes de u
> a_P(i,Nx)=(F_n_u(i,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_F*9.0d0/3.0d0*dy)/
> alpha_u
> a_W(i,Nx)=F_w_u(i,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
> a_E(i,Nx)=1.0d-20 !0.0d0
> a_S(i,Nx)=F_s_u(i,Nx)*dx+F_n_u(i,Nx)*1.0d0/8.0d0*dx+D_s*dx
> a_N(i,Nx)=D_n*dx
>
> S_fake(i,Nx)=(3.0d0/8.0d0*u_fake(i,Nx)*dy-2.0d0/8.0d0*u_fake(i,Nx-1)*dy-1
> .0d0/8.0d0*&
> u_fake(i,Nx-2)*dy)*F_w_u(i,Nx)+(-u_F(i)*dy)*F_e_u(i,Nx)
> +(3.0d0/8.0d0*u_fake(i,Nx)*dx&
> -2.0d0/8.0d0*u_fake(i-1,Nx)*dx-1.0d0/8.0d0*u_fake(i-2,Nx)*dx)*F_s_u(i,Nx
> )
> +(-3.0d0/&
> 8.0d0*u_fake(i+1,Nx)*dx+2.0d0/8.0d0*u_fake(i,Nx)*dx)*F_n_u(i,Nx)+
> (p_fake(i,Nx-1)-&
> p_fake(i,Nx))*dy+D_F*8.0d0/3.0d0*u_F(i)*dy+((1.0d0-
> alpha_u)*a_P(i,Nx))*u_old(i,Nx)
>
>
> ! a_P(i,Nx)=1.0d0
> ! a_W(i,Nx)=1.0d-20
> ! a_E(i,Nx)=1.0d-20
> ! a_S(i,Nx)=1.0d-20
> ! a_N(i,Nx)=1.0d-20
> ! S_fake(i,Nx)=u_F(i)
>
>
>
>
>
> !coeficientes de v
> b_P(i,Nx)=(F_n_v(i,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_J*9.0d0/3.0d0*dy)/
> alpha_v
> b_W(i,Nx)=F_w_v(i,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
> b_E(i,Nx)=1.0d-20 !0.0d0
> b_S(i,Nx)=F_s_v(i,Nx)*dx+F_n_v(i,Nx)*1.0d0/8.0d0*dx+D_s*dx
> b_N(i,Nx)=D_n*dx
>
> S_v_fake(i,Nx)=(3.0d0/8.0d0*v_fake(i,Nx)*dy-2.0d0/8.0d0*v_fake(i,Nx-1)*dy
> -1.0d0/
> &
> 8.0d0*v_fake(i,Nx-2)*dy)*F_w_v(i,Nx)+(-v_J*dy)*F_e_v(i,Nx)
> +(3.0d0/8.0d0*v_fake(i,Nx)&
>
> *dx-2.0d0/8.0d0*v_fake(i-1,Nx)*dx-1.0d0/8.0d0*v_fake(i-2,Nx)*dx)*F_s_v(i,Nx
> )
> +(-3.0d0&
> /8.0d0*v_fake(i+1,Nx)*dx+2.0d0/8.0d0*v_fake(i,Nx)*dx)*F_n_v(i,Nx)+
> (p_fake(i-1,Nx)-&
> p_fake(i,Nx))*dx+D_J*8.0d0/3.0d0*v_J*dy+((1.0d0-
> alpha_v)*b_P(i,Nx))*v_old(i,Nx)
>
>
> ! b_P(i,Nx)=1.0d0
> ! b_W(i,Nx)=1.0d-20
> ! b_E(i,Nx)=1.0d-20
> ! b_S(i,Nx)=1.0d-20
> ! b_N(i,Nx)=1.0d-20
> ! S_v_fake(i,Nx)=v_J
>
>
>
>
> !coeficientes de p
> c_W(i,Nx)=dy/a_P(i,Nx)*dy
> c_E(i,Nx)=dy/a_E(i,Nx)*alpha_u*dy
> c_S(i,Nx)=dx/b_P(i,Nx)*dx
> c_N(i,Nx)=dx/b_N(i,Nx)*alpha_v*dx
> c_P(i,Nx)=c_W(i,Nx)+c_E(i,Nx)+c_S(i,Nx)+c_N(i,Nx)
> S_p_fake(i,Nx)=v_fake(i,Nx)*dx-v_fake(i+1,Nx)*dx !Faltam termos de
> u.
>
> end do
>
> !fronteira E - linha 2
>
> !coeficientes de u
> a_P(2,Nx)=(F_n_u(2,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_F*9.0d0/3.0d0*dy)/
> alpha_u
> a_W(2,Nx)=F_w_u(2,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
> a_E(2,Nx)=1.0d-20 !0.0d0
> a_S(2,Nx)=F_s_u(2,Nx)*7.0d0/8.0d0*dx+F_n_u(2,Nx)*1.0d0/8.0d0*dx
> +D_s*dx
> a_N(2,Nx)=D_n*dx
>
> S_fake(2,Nx)=(3.0d0/8.0d0*u_fake(2,Nx)*dy-2.0d0/8.0d0*u_fake(2,Nx-1)*dy-1
> .0d0/8.0d0*&
> u_fake(2,Nx-2)*dy)*F_w_u(2,Nx)+(-u_F(2)*dy)*F_e_u(2,Nx)
> +(3.0d0/8.0d0*u_fake(2,Nx)*dx&
> -2.0d0/8.0d0*u_G*dx)*F_s_u(2,Nx)+(-3.0d0/8.0d0*u_fake(3,Nx)*dx
> +2.0d0/8.0d0*&
> u_fake(2,Nx)*dx)*F_n_u(2,Nx)+(p_fake(2,Nx-1)-p_fake(2,Nx))*dy
> +D_F*8.0d0/3.0d0*u_F(2)&
> *dy+((1.0d0-alpha_u)*a_P(2,Nx))*u_old(2,Nx)
>
>
> ! a_P(2,Nx)=1.0d0
> ! a_W(2,Nx)=1.0d-20
> ! a_E(2,Nx)=1.0d-20
> ! a_S(2,Nx)=1.0d-20
> ! a_N(2,Nx)=1.0d-20
> ! S_fake(2,Nx)=u_F(2)
>
>
>
> !coeficientes de v
> b_P(2,Nx)=(F_n_v(2,Nx)*dx+D_w*dy+D_s*dx+D_n*dx+D_J*9.0d0/3.0d0*dy)/
> alpha_v
> b_W(2,Nx)=F_w_v(2,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
> b_E(2,Nx)=1.0d-20 !0.0d0
> b_S(2,Nx)=F_s_v(2,Nx)*7.0d0/8.0d0*dx+F_n_v(2,Nx)*1.0d0/8.0d0*dx
> +D_s*dx
> b_N(2,Nx)=D_n*dx
>
> S_v_fake(2,Nx)=(3.0d0/8.0d0*v_fake(2,Nx)*dy-2.0d0/8.0d0*v_fake(2,Nx-1)*dy
> -1.0d0/
> &
> 8.0d0*v_fake(2,Nx-2)*dy)*F_w_v(2,Nx)+(-v_J*dy)*F_e_v(2,Nx)
> +(3.0d0/8.0d0*v_fake(2,Nx)&
> *dx-2.0d0/8.0d0*v_K*dx)*F_s_v(2,Nx)+(-3.0d0/8.0d0*v_fake(3,Nx)*dx
> +2.0d0/8.0d0*&
> v_fake(2,Nx)*dx)*F_n_v(2,Nx)+(p_fake(1,Nx)-p_fake(2,Nx))*dx
> +D_J*8.0d0/3.0d0*v_J*dy+(&
> (1.0d0-alpha_v)*b_P(2,Nx))*v_old(2,Nx)
>
>
> ! b_P(2,Nx)=1.0d0
> ! b_W(2,Nx)=1.0d-20
> ! b_E(2,Nx)=1.0d-20
> ! b_S(2,Nx)=1.0d-20
> ! b_N(2,Nx)=1.0d-20
> ! S_v_fake(2,Nx)=v_J
>
> !coeficientes de p
> c_W(2,Nx)=dy/a_P(2,Nx)*dy
> c_E(2,Nx)=dy/a_E(2,Nx)*alpha_u*dy
> c_S(2,Nx)=dx/b_P(2,Nx)*dx
> c_N(2,Nx)=dx/b_N(2,Nx)*alpha_v*dx
> c_P(2,Nx)=c_W(2,Nx)+c_E(2,Nx)+c_S(2,Nx)+c_N(2,Nx)
> S_p_fake(2,Nx)=v_fake(2,Nx)*dx-v_fake(3,Nx)*dx !Faltam termos de u.
>
>
> !fronteira S - com exce��o da coluna 2
>
> do j=3,Nx-1
> !coeficientes de u
> a_P(1,j)=(F_e_u(1,j)*dy+F_n_u(1,j)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy
> +D_n*dx+D_G*9.0d0/&
> 3.0d0*dx)/alpha_u
> a_W(1,j)=F_w_u(1,j)*dy+F_e_u(1,j)*1.0d0/8.0d0*dy+D_w*dy
> a_E(1,j)=D_e*dy
> a_S(1,j)=1.0d-20 !0.0d0
> a_N(1,j)=D_n*dx+D_G*1.0d0/3.0d0*dx
>
> S_fake(1,j)=(3.0d0/8.0d0*u_fake(1,j)*dy-2.0d0/8.0d0*u_fake(1,j-1)*dy-1.0d
> 0/8.0d0*&
> u_fake(1,j-2)*dy)*F_w_u(1,j)+(-3.0d0/8.0d0*u_fake(1,j+1)*dy
> +2.0d0/8.0d0*u_fake(1,j)*&
> dy)*F_e_u(1,j)+(u_G*dx)*F_s_u(1,j)+(-3.0d0/8.0d0*u_fake(2,j)*dx
> +2.0d0/8.0d0*u_G*dx)*&
> F_n_u(1,j)+(p_fake(1,j-1)-p_fake(1,j))*dy+D_G*8.0d0/3.0d0*u_G*dx+
> ((1.0d0-alpha_u)*&
> a_P(1,j))*u_old(1,j)
>
>
> ! a_P(1,j)=1.0d0
> ! a_W(1,j)=1.0d-20
> ! a_E(1,j)=1.0d-20
> ! a_S(1,j)=1.0d-20
> ! a_N(1,j)=1.0d-20
> ! S_fake(1,j)=u_G
>
>
>
>
>
>
> !coeficientes de v
> b_P(1,j)=(F_e_v(1,j)*dy+F_n_v(1,j)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy
> +D_n*dx+D_K*9.0d0/&
> 3.0d0*dx)/alpha_v
> b_W(1,j)=F_w_v(1,j)*dy+F_e_v(1,j)*1.0d0/8.0d0*dy+D_w*dy
> b_E(1,j)=D_e*dy
> b_S(1,j)=1.0d-20 !0.0d0
> b_N(1,j)=D_n*dx+D_K*1.0d0/3.0d0*dx
>
> S_v_fake(1,j)=(3.0d0/8.0d0*v_fake(1,j)*dy-2.0d0/8.0d0*v_fake(1,j-1)*dy-1.
> 0d0/8.0d0*&
> v_fake(1,j-2)*dy)*F_w_v(1,j)+(-3.0d0/8.0d0*v_fake(1,j+1)*dy
> +2.0d0/8.0d0*v_fake(1,j)*&
> dy)*F_e_v(1,j)+(v_K*dx)*F_s_v(1,j)+(-3.0d0/8.0d0*v_fake(2,j)*dx
> +2.0d0/8.0d0*v_K*dx)*&
> F_n_v(1,j)+D_K*8.0d0/3.0d0*v_K*dx+((1.0d0-
> alpha_v)*b_P(1,j))*v_old(1,j) !Falta o termo de press�o.
>
> ! b_P(1,j)=1.0d0
> ! b_W(1,j)=1.0d-20
> ! b_E(1,j)=1.0d-20
> ! b_S(1,j)=1.0d-20
> ! b_N(1,j)=1.0d-20
> ! S_v_fake(1,j)=v_K
>
>
>
>
>
> !coeficientes de p
> c_W(1,j)=dy/a_P(1,j)*dy
> c_E(1,j)=dy/a_E(1,j)*alpha_u*dy
> c_S(1,j)=dx/b_P(1,j)*dx
> c_N(1,j)=dx/b_N(1,j)*alpha_v*dx
> c_P(1,j)=c_W(1,j)+c_E(1,j)+c_S(1,j)+c_N(1,j)
> S_p_fake(1,j)=u_fake(1,j)*dy-u_fake(1,j+1)*dy+v_fake(1,j)*dx-
> v_fake(2,j)*dx
>
> end do
>
> !fronteira S - coluna 2
>
> !coeficientes de u
> a_P(1,2)=(F_e_u(1,2)*dy+F_n_u(1,2)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy+D_n*dx
> +D_G*9.0d0/3.0d0*dx)/&
> alpha_u
> a_W(1,2)=F_w_u(1,2)*7.0d0/8.0d0*dy+F_e_u(1,2)*1.0d0/8.0d0*dy+D_w*dy
> a_E(1,2)=D_e*dy
> a_S(1,2)=1.0d-20 !0.0d0
> a_N(1,2)=D_n*dx+D_G*1.0d0/3.0d0*dx
>
> S_fake(1,2)=(3.0d0/8.0d0*u_fake(1,2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(1,2)+
> (-3.0d0/8.0d0*&
> u_fake(1,3)*dy+2.0d0/8.0d0*u_fake(1,2)*dy)*F_e_u(1,2)+
> (u_G*dx)*F_s_u(1,2)+(-3.0d0/8.0d0*&
> u_fake(2,2)*dx+2.0d0/8.0d0*u_G*dx)*F_n_u(1,2)+(p_fake(1,1)-
> p_fake(1,2))*dy+D_G*8.0d0/3.0d0*&
> u_G*dx+((1.0d0-alpha_u)*a_P(1,2))*u_old(1,2)
>
> ! a_P(1,2)=1.0d0
> ! a_W(1,2)=1.0d-20
> ! a_E(1,2)=1.0d-20
> ! a_S(1,2)=1.0d-20
> ! a_N(1,2)=1.0d-20
> ! S_fake(1,2)=u_G
>
>
>
>
>
>
>
>
> !coeficientes de v
> b_P(1,2)=(F_e_v(1,2)*dy+F_n_v(1,2)*7.0d0/8.0d0*dx+D_w*dy+D_e*dy+D_n*dx
> +D_K*9.0d0/3.0d0*dx)/&
> alpha_v
> b_W(1,2)=F_w_v(1,2)*7.0d0/8.0d0*dy+F_e_v(1,2)*1.0d0/8.0d0*dy+D_w*dy
> b_E(1,2)=D_e*dy
> b_S(1,2)=1.0d-20 !0.0d0
> b_N(1,2)=D_n*dx+D_K*1.0d0/3.0d0*dx
>
> S_v_fake(1,2)=(3.0d0/8.0d0*v_fake(1,2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(1,2)+
> (-3.0d0/8.0d0*&
> v_fake(1,3)*dy+2.0d0/8.0d0*v_fake(1,2)*dy)*F_e_v(1,2)+
> (v_K*dx)*F_s_v(1,2)+(-3.0d0/8.0d0*&
> v_fake(2,2)*dx+2.0d0/8.0d0*v_K*dx)*F_n_v(1,2)+D_K*8.0d0/3.0d0*v_K*dx+
> ((1.0d0-alpha_v)*&
> b_P(1,2))*v_old(1,2) !Falta o termo de press�o.
>
> ! b_P(1,2)=1.0d0
> ! b_W(1,2)=1.0d-20
> ! b_E(1,2)=1.0d-20
> ! b_S(1,2)=1.0d-20
> ! b_N(1,2)=1.0d-20
> ! S_v_fake(1,2)=v_K
>
>
>
>
>
> !coeficientes de p
> c_W(1,2)=dy/a_P(1,2)*dy
> c_E(1,2)=dy/a_E(1,2)*alpha_u*dy
> c_S(1,2)=dx/b_P(1,2)*dx
> c_N(1,2)=dx/b_N(1,2)*alpha_v*dx
> c_P(1,2)=c_W(1,2)+c_E(1,2)+c_S(1,2)+c_N(1,2)
> S_p_fake(1,2)=u_fake(1,2)*dy-u_fake(1,3)*dy+v_fake(1,2)*dx-
> v_fake(2,2)*dx
>
>
> !fronteira N - com exce��o da coluna 2
>
> do j=3,Nx-1
> !coeficientes de u
> a_P(Ny,j)=(F_e_u(Ny,j)*dy+D_w*dy+D_e*dy+D_s*dx+D_H*9.0d0/3.0d0*dx)/
> alpha_u
> a_W(Ny,j)=F_w_u(Ny,j)*dy+F_e_u(Ny,j)*1.0d0/8.0d0*dy+D_w*dy
> a_E(Ny,j)=D_e*dy
> a_S(Ny,j)=F_s_u(Ny,j)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
> a_N(Ny,j)=1.0d-20 !0.0d0
>
> S_fake(Ny,j)=(3.0d0/8.0d0*u_fake(Ny,j)*dy-2.0d0/8.0d0*u_fake(Ny,j-1)*dy-1
> .0d0/8.0d0*&
> u_fake(Ny,j-2)*dy)*F_w_u(Ny,j)+(-3.0d0/8.0d0*u_fake(Ny,j+1)*dy
> +2.0d0/8.0d0*&
> u_fake(Ny,j)*dy)*F_e_u(Ny,j)
> +(3.0d0/8.0d0*u_fake(Ny,j)*dx-2.0d0/8.0d0*u_fake(Ny-1,j)&
> *dx-1.0d0/8.0d0*u_fake(Ny-2,j)*dx)*F_s_u(Ny,j)+(-u_H*dx)*F_n_u(Ny,j)
> +(p_fake(Ny,j-1)&
> -p_fake(Ny,j))*dy+D_H*8.0d0/3.0d0*u_H*dx+((1.0d0-
> alpha_u)*a_P(Ny,j))*u_old(Ny,j)
>
>
> ! a_P(Ny,j)=1.0d0
> ! a_W(Ny,j)=1.0d-20
> ! a_E(Ny,j)=1.0d-20
> ! a_S(Ny,j)=1.0d-20
> ! a_N(Ny,j)=1.0d-20
> ! S_fake(Ny,j)=u_H
>
>
>
>
> !coeficientes de v
> b_P(Ny,j)=(F_e_v(Ny,j)*dy+D_w*dy+D_e*dy+D_s*dx+D_L*9.0d0/3.0d0*dx)/
> alpha_v
> b_W(Ny,j)=F_w_v(Ny,j)*dy+F_e_v(Ny,j)*1.0d0/8.0d0*dy+D_w*dy
> b_E(Ny,j)=D_e*dy
> b_S(Ny,j)=F_s_v(Ny,j)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
> b_N(Ny,j)=1.0d-20 !0.0d0
>
> S_v_fake(Ny,j)=(3.0d0/8.0d0*v_fake(Ny,j)*dy-2.0d0/8.0d0*v_fake(Ny,j-1)*dy
> -1.0d0/
> &
> 8.0d0*v_fake(Ny,j-2)*dy)*F_w_v(Ny,j)+(-3.0d0/8.0d0*v_fake(Ny,j+1)*dy
> +2.0d0/8.0d0*&
> v_fake(Ny,j)*dy)*F_e_v(Ny,j)
> +(3.0d0/8.0d0*v_fake(Ny,j)*dx-2.0d0/8.0d0*v_fake(Ny-1,j)&
> *dx-1.0d0/8.0d0*v_fake(Ny-2,j)*dx)*F_s_v(Ny,j)+(-v_L*dx)*F_n_v(Ny,j)
> +(p_fake(Ny-1,j)&
> -p_fake(Ny,j))*dx+D_L*8.0d0/3.0d0*v_L*dx+((1.0d0-
> alpha_v)*b_P(Ny,j))*v_old(Ny,j)
>
>
> ! b_P(Ny,j)=1.0d0
> ! b_W(Ny,j)=1.0d-20
> ! b_E(Ny,j)=1.0d-20
> ! b_S(Ny,j)=1.0d-20
> ! b_N(Ny,j)=1.0d-20
> ! S_v_fake(Ny,j)=v_L
>
>
>
>
>
> !coeficientes de p
> c_W(Ny,j)=dy/a_P(Ny,j)*dy
> c_E(Ny,j)=dy/a_E(Ny,j)*alpha_u*dy
> c_S(Ny,j)=dx/b_P(Ny,j)*dx
> c_N(Ny,j)=dx/b_N(Ny,j)*alpha_v*dx
> c_P(Ny,j)=c_W(Ny,j)+c_E(Ny,j)+c_S(Ny,j)+c_N(Ny,j)
> S_p_fake(Ny,j)=u_fake(Ny,j)*dy-u_fake(Ny,j+1)*dy !Faltam termos de
> v.
>
>
> end do
>
> !fronteira N - coluna 2
>
> !coeficientes de u
> a_P(Ny,2)=(F_e_u(Ny,2)*dy+D_w*dy+D_e*dy+D_s*dx+D_H*9.0d0/3.0d0*dx)/
> alpha_u
> a_W(Ny,2)=F_w_u(Ny,2)*7.0d0/8.0d0*dy+F_e_u(Ny,2)*1.0d0/8.0d0*dy
> +D_w*dy
> a_E(Ny,2)=D_e*dy
> a_S(Ny,2)=F_s_u(Ny,2)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
> a_N(Ny,2)=1.0d-20 !0.0d0
> S_fake(Ny,2)=(3.0d0/8.0d0*u_fake(Ny,
> 2)*dy-2.0d0/8.0d0*u_E*dy)*F_w_u(Ny,2)+(-3.0d0/&
> 8.0d0*u_fake(Ny,3)*dy+2.0d0/8.0d0*u_fake(Ny,2)*dy)*F_e_u(Ny,
> 2)+(3.0d0/8.0d0*&
> u_fake(Ny,
> 2)*dx-2.0d0/8.0d0*u_fake(Ny-1,2)*dx-1.0d0/8.0d0*u_fake(Ny-2,2)*dx)*&
> F_s_u(Ny,2)+(-u_H*dx)*F_n_u(Ny,2)+(p_fake(Ny,1)-p_fake(Ny,2))*dy
> +D_H*8.0d0/3.0d0*u_H&
> *dx+((1.0d0-alpha_u)*a_P(Ny,2))*u_old(Ny,2)
>
> ! a_P(Ny,2)=1.0d0
> ! a_W(Ny,2)=1.0d-20
> ! a_E(Ny,2)=1.0d-20
> ! a_S(Ny,2)=1.0d-20
> ! a_N(Ny,2)=1.0d-20
> ! S_fake(Ny,2)=u_H
>
>
>
>
> !coeficientes de v
> b_P(Ny,2)=(F_e_v(Ny,2)*dy+D_w*dy+D_e*dy+D_s*dx+D_L*9.0d0/3.0d0*dx)/
> alpha_v
> b_W(Ny,2)=F_w_v(Ny,2)*7.0d0/8.0d0*dy+F_e_v(Ny,2)*1.0d0/8.0d0*dy
> +D_w*dy
> b_E(Ny,2)=D_e*dy
> b_S(Ny,2)=F_s_v(Ny,2)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
> b_N(Ny,2)=1.0d-20 !0.0d0
> S_v_fake(Ny,2)=(3.0d0/8.0d0*v_fake(Ny,
> 2)*dy-2.0d0/8.0d0*v_I*dy)*F_w_v(Ny,2)+(-3.0d0/&
> 8.0d0*v_fake(Ny,3)*dy+2.0d0/8.0d0*v_fake(Ny,2)*dy)*F_e_v(Ny,
> 2)+(3.0d0/8.0d0*&
> v_fake(Ny,
> 2)*dx-2.0d0/8.0d0*v_fake(Ny-1,2)*dx-1.0d0/8.0d0*v_fake(Ny-2,2)*dx)*&
> F_s_v(Ny,2)+(-v_L*dx)*F_n_v(Ny,2)+(p_fake(Ny-1,2)-p_fake(Ny,2))*dx
> +D_L*8.0d0/3.0d0*&
> v_L*dx+((1.0d0-alpha_v)*b_P(Ny,2))*v_old(Ny,2)
>
> ! b_P(Ny,2)=1.0d0
> ! b_W(Ny,2)=1.0d-20
> ! b_E(Ny,2)=1.0d-20
> ! b_S(Ny,2)=1.0d-20
> ! b_N(Ny,2)=1.0d-20
> ! S_v_fake(Ny,2)=v_L
>
>
>
>
>
>
>
> !coeficientes de p
> c_W(Ny,2)=dy/a_P(Ny,2)*dy
> c_E(Ny,2)=dy/a_E(Ny,2)*alpha_u*dy
> c_S(Ny,2)=dx/b_P(Ny,2)*dx
> c_N(Ny,2)=dx/b_N(Ny,2)*alpha_v*dx
> c_P(Ny,2)=c_W(Ny,2)+c_E(Ny,2)+c_S(Ny,2)+c_N(Ny,2)
> S_p_fake(Ny,2)=u_fake(Ny,2)*dy-u_fake(Ny,3)*dy !Faltam termos de v.
>
>
> !fronteira SW
>
> !coeficientes de u
> a_P(1,1)=(F_e_u(1,1)*7.0d0/8.0d0*dy+F_n_u(1,1)*7.0d0/8.0d0*dx+D_e*dy
> +D_n*dx+D_EE*&
> 9.0d0/3.0d0*dy+D_G*9.0d0/3.0d0*dx)/alpha_u
> a_W(1,1)=1.0d-20 !0.0d0
> a_E(1,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
> a_S(1,1)=1.0d-20 !0.0d0
> a_N(1,1)=D_n*dx+D_G*1.0d0/3.0d0*dx
> S_fake(1,1)=(u_E*dy)*F_w_u(1,1)+(-3.0d0/8.0d0*u_fake(1,2)*dy
> +2.0d0/8.0d0*u_E*dy)*&
> F_e_u(1,1)+(u_G*dx)*F_s_u(1,1)+(-3.0d0/8.0d0*u_fake(2,1)*dx
> +2.0d0/8.0d0*u_G*dx)*&
> F_n_u(1,1)+D_EE*8.0d0/3.0d0*u_E*dy+D_G*8.0d0/3.0d0*u_G*dx+((1.0d0-
> alpha_u)*a_P(1,1))&
> *u_old(1,1) !Falta o termo de press�o.
>
>
> ! a_P(1,1)=1.0d0
> ! a_W(1,1)=1.0d-20
> ! a_E(1,1)=1.0d-20
> ! a_S(1,1)=1.0d-20
> ! a_N(1,1)=1.0d-20
> ! S_fake(1,1)=1.0d-30
>
>
>
>
>
>
>
>
> !coeficientes de v
> b_P(1,1)=(F_e_v(1,1)*7.0d0/8.0d0*dy+F_n_v(1,1)*7.0d0/8.0d0*dx+D_e*dy
> +D_n*dx+D_I*&
> 9.0d0/3.0d0*dy+D_K*9.0d0/3.0d0*dx)/alpha_v
> b_W(1,1)=1.0d-20 !0.0d0
> b_E(1,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
> b_S(1,1)=1.0d-20 !0.0d0
> b_N(1,1)=D_n*dx+D_K*1.0d0/3.0d0*dx
> S_v_fake(1,1)=(v_I*dy)*F_w_v(1,1)+(-3.0d0/8.0d0*v_fake(1,2)*dy
> +2.0d0/8.0d0*v_I*dy)*&
> F_e_v(1,1)+(v_K*dx)*F_s_v(1,1)+(-3.0d0/8.0d0*v_fake(2,1)*dx
> +2.0d0/8.0d0*v_K*dx)*&
> F_n_v(1,1)+D_I*8.0d0/3.0d0*v_I*dy+D_K*8.0d0/3.0d0*v_K*dx+((1.0d0-
> alpha_v)*b_P(1,1))&
> *v_old(1,1) !Falta o termo de press�o.
>
> ! b_P(1,1)=1.0d0
> ! b_W(1,1)=1.0d-20
> ! b_E(1,1)=1.0d-20
> ! b_S(1,1)=1.0d-20
> ! b_N(1,1)=1.0d-20
> ! S_v_fake(1,1)=1.0d-30
>
>
>
>
> !coeficientes de p
> c_W(1,1)=dy/a_P(1,1)*dy
> c_E(1,1)=dy/a_E(1,1)*alpha_u*dy
> c_S(1,1)=dx/b_P(1,1)*dx
> c_N(1,1)=dx/b_N(1,1)*alpha_v*dx
> c_P(1,1)=c_W(1,1)+c_E(1,1)+c_S(1,1)+c_N(1,1)
> S_p_fake(1,1)=u_fake(1,1)*dy-u_fake(1,2)*dy+v_fake(1,1)*dx-
> v_fake(2,1)*dx
>
>
> !fronteira NW
>
> !coeficientes de u
> a_P(Ny,1)=(F_e_u(Ny,1)*7.0d0/8.0d0*dy+D_e*dy+D_s*dx
> +D_EE*9.0d0/3.0d0*dy+D_H*9.0d0/&
> 3.0d0*dx)/alpha_u
> a_W(Ny,1)=1.0d-20 !0.0d0
> a_E(Ny,1)=D_e*dy+D_EE*1.0d0/3.0d0*dy
> a_S(Ny,1)=F_s_u(Ny,1)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
> a_N(Ny,1)=1.0d-20 !0.0d0
> S_fake(Ny,1)=(u_E*dy)*F_w_u(Ny,1)+(-3.0d0/8.0d0*u_fake(Ny,2)*dy
> +2.0d0/8.0d0*u_E*dy)*&
> F_e_u(Ny,1)+(3.0d0/8.0d0*u_fake(Ny,
> 1)*dx-2.0d0/8.0d0*u_fake(Ny-1,1)*dx-1.0d0/8.0d0*&
> u_fake(Ny-2,1)*dx)*F_s_u(Ny,1)+(-u_H*dx)*F_n_u(Ny,
> 1)+D_EE*8.0d0/3.0d0*u_E*dy+D_H*&
> 8.0d0/3.0d0*u_H*dx+((1.0d0-alpha_u)*a_P(Ny,1))*u_old(Ny,1) ! Falta
> o termo de press�o.
>
>
> ! a_P(Ny,1)=1.0d0
> ! a_W(Ny,1)=1.0d-20
> ! a_E(Ny,1)=1.0d-20
> ! a_S(Ny,1)=1.0d-20
> ! a_N(Ny,1)=1.0d-20
> ! S_fake(Ny,1)=1.0d-30
>
>
>
>
>
> !coeficientes de v
> b_P(Ny,1)=(F_e_v(Ny,1)*7.0d0/8.0d0*dy+D_e*dy+D_s*dx
> +D_I*9.0d0/3.0d0*dy+D_L*9.0d0/&
> 3.0d0*dx)/alpha_v
> b_W(Ny,1)=1.0d-20 !0.0d0
> b_E(Ny,1)=D_e*dy+D_I*1.0d0/3.0d0*dy
> b_S(Ny,1)=F_s_v(Ny,1)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
> b_N(Ny,1)=1.0d-20 !0.0d0
> S_v_fake(Ny,1)=(v_I*dy)*F_w_v(Ny,1)+(-3.0d0/8.0d0*v_fake(Ny,2)*dy
> +2.0d0/8.0d0*v_I*dy&
> )*F_e_v(Ny,1)+(3.0d0/8.0d0*v_fake(Ny,
> 1)*dx-2.0d0/8.0d0*v_fake(Ny-1,1)*dx-1.0d0/8.0d0&
> *v_fake(Ny-2,1)*dx)*F_s_v(Ny,1)+(-v_L*dx)*F_n_v(Ny,1)+
> (p_fake(Ny-1,1)-p_fake(Ny,1))*&
> dx+D_I*8.0d0/3.0d0*v_I*dy+D_L*8.0d0/3.0d0*v_L*dx+((1.0d0-
> alpha_v)*b_P(Ny,1))*&
> v_old(Ny,1)
>
>
> ! b_P(Ny,1)=1.0d0
> ! b_W(Ny,1)=1.0d-20
> ! b_E(Ny,1)=1.0d-20
> ! b_S(Ny,1)=1.0d-20
> ! b_N(Ny,1)=1.0d-20
> ! S_v_fake(Ny,1)=1.0d-30
>
>
> !coeficientes de p
> c_W(Ny,1)=dy/a_P(Ny,1)*dy
> c_E(Ny,1)=dy/a_E(Ny,1)*alpha_u*dy
> c_S(Ny,1)=dx/b_P(Ny,1)*dx
> c_N(Ny,1)=dx/b_N(Ny,1)*alpha_v*dx
> c_P(Ny,1)=c_W(Ny,1)+c_E(Ny,1)+c_S(Ny,1)+c_N(Ny,1)
> S_p_fake(Ny,1)=u_fake(Ny,1)*dy-u_fake(Ny,2)*dy !Faltam termos de v.
>
>
> !fronteira SE
>
> !coeficientes de u
> a_P(1,Nx)=(F_n_u(1,Nx)*7.0d0/8.0d0*dx+D_w*dy+D_n*dx
> +D_F*9.0d0/3.0d0*dy+D_G*9.0d0/&
> 3.0d0*dx)/alpha_u
> a_W(1,Nx)=F_w_u(1,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
> a_E(1,Nx)=1.0d-20 !0.0d0
> a_S(1,Nx)=1.0d-20 !0.0d0
> a_N(1,Nx)=D_n*dx+D_G*1.0d0/3.0d0*dx
>
> S_fake(1,Nx)=(3.0d0/8.0d0*u_fake(1,Nx)*dy-2.0d0/8.0d0*u_fake(1,Nx-1)*dy-1
> .0d0/8.0d0*&
> u_fake(1,Nx-2)*dy)*F_w_u(1,Nx)+(-u_F(1)*dy)*F_e_u(1,Nx)+
> (u_G*dx)*F_s_u(1,Nx)+(-3.0d0&
> /8.0d0*u_fake(2,Nx)*dx+2.0d0/8.0d0*u_G*dx)*F_n_u(1,Nx)+
> (p_fake(1,Nx-1)-p_fake(1,Nx))&
> *dy+D_F*8.0d0/3.0d0*u_F(1)*dy+D_G*8.0d0/3.0d0*u_G*dx+((1.0d0-
> alpha_u)*a_P(1,Nx))*&
> u_old(1,Nx)
>
> ! a_P(1,Nx)=1.0d0
> ! a_W(1,Nx)=1.0d-20
> ! a_E(1,Nx)=1.0d-20
> ! a_S(1,Nx)=1.0d-20
> ! a_N(1,Nx)=1.0d-20
> ! S_fake(1,Nx)=1.0d-30
>
>
>
>
>
>
>
>
> !coeficientes de v
> b_P(1,Nx)=(F_n_v(1,Nx)*7.0d0/8.0d0*dx+D_w*dy+D_n*dx
> +D_J*9.0d0/3.0d0*dy+D_K*9.0d0/&
> 3.0d0*dx)/alpha_v
> b_W(1,Nx)=F_w_v(1,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
> b_E(1,Nx)=1.0d-20 !0.0d0
> b_S(1,Nx)=1.0d-20 !0.0d0
> b_N(1,Nx)=D_n*dx+D_K*1.0d0/3.0d0*dx
>
> S_v_fake(1,Nx)=(3.0d0/8.0d0*v_fake(1,Nx)*dy-2.0d0/8.0d0*v_fake(1,Nx-1)*dy
> -1.0d0/
> &
> 8.0d0*v_fake(1,Nx-2)*dy)*F_w_v(1,Nx)+(-v_J*dy)*F_e_v(1,Nx)+
> (v_K*dx)*F_s_v(1,Nx)+(&
> -3.0d0/8.0d0*v_fake(2,Nx)*dx+2.0d0/8.0d0*v_K*dx)*F_n_v(1,Nx)
> +D_J*8.0d0/3.0d0*v_J*dy+&
> D_K*8.0d0/3.0d0*v_K*dx+((1.0d0-alpha_v)*b_P(1,Nx))*v_old(1,Nx) !
> Falta o termo de press�o.
>
>
> ! b_P(1,Nx)=1.0d0
> ! b_W(1,Nx)=1.0d-20
> ! b_E(1,Nx)=1.0d-20
> ! b_S(1,Nx)=1.0d-20
> ! b_N(1,Nx)=1.0d-20
> ! S_v_fake(1,Nx)=1.0d-30
>
>
>
>
> !coeficientes de p
> c_W(1,Nx)=dy/a_P(1,Nx)*dy
> c_E(1,Nx)=dy/a_E(1,Nx)*alpha_u*dy
> c_S(1,Nx)=dx/b_P(1,Nx)*dx
> c_N(1,Nx)=dx/b_N(1,Nx)*alpha_v*dx
> c_P(1,Nx)=c_W(1,Nx)+c_E(1,Nx)+c_S(1,Nx)+c_N(1,Nx)
> S_p_fake(1,Nx)=v_fake(1,Nx)*dx-v_fake(2,Nx)*dx !Faltam termos de u.
>
>
> !fronteira NE
>
> !coeficientes de u
> a_P(Ny,Nx)=(D_w*dy+D_s*dx+D_F*9.0d0/3.0d0*dy+D_H*9.0d0/3.0d0*dx)/
> alpha_u
> a_W(Ny,Nx)=F_w_u(Ny,Nx)*dy+D_w*dy+D_F*1.0d0/3.0d0*dy
> a_E(Ny,Nx)=1.0d-20 !0.0d0
> a_S(Ny,Nx)=F_s_u(Ny,Nx)*dx+D_s*dx+D_H*1.0d0/3.0d0*dx
> a_N(Ny,Nx)=1.0d-20 !0.0d0
>
> S_fake(Ny,Nx)=(3.0d0/8.0d0*u_fake(Ny,Nx)*dy-2.0d0/8.0d0*u_fake(Ny,Nx-1)*d
> y-1.0d0/
> &
> 8.0d0*u_fake(Ny,Nx-2)*dy)*F_w_u(Ny,Nx)+(-u_F(Ny)*dy)*F_e_u(Ny,Nx)
> +(3.0d0/8.0d0*&
>
> u_fake(Ny,Nx)*dx-2.0d0/8.0d0*u_fake(Ny-1,Nx)*dx-1.0d0/8.0d0*u_fake(Ny-2,Nx)
> *dx)*&
> F_s_u(Ny,Nx)+(-u_H*dx)*F_n_u(Ny,Nx)+(p_fake(Ny,Nx-1)-
> p_fake(Ny,Nx))*dy+D_F*8.0d0/&
> 3.0d0*u_F(Ny)*dy+D_H*8.0d0/3.0d0*u_H*dx+((1.0d0-
> alpha_u)*a_P(Ny,Nx))*u_old(Ny,Nx)
>
> ! a_P(Ny,Nx)=1.0d0
> ! a_W(Ny,Nx)=1.0d-20
> ! a_E(Ny,Nx)=1.0d-20
> ! a_S(Ny,Nx)=1.0d-20
> ! a_N(Ny,Nx)=1.0d-20
> ! S_fake(Ny,Nx)=1.0d-30
>
>
>
>
> !coeficientes de v
> b_P(Ny,Nx)=(D_w*dy+D_s*dx+D_J*9.0d0/3.0d0*dy+D_L*9.0d0/3.0d0*dx)/
> alpha_v
> b_W(Ny,Nx)=F_w_v(Ny,Nx)*dy+D_w*dy+D_J*1.0d0/3.0d0*dy
> b_E(Ny,Nx)=1.0d-20 !0.0d0
> b_S(Ny,Nx)=F_s_v(Ny,Nx)*dx+D_s*dx+D_L*1.0d0/3.0d0*dx
> b_N(Ny,Nx)=1.0d-20 !0.0d0
>
> S_v_fake(Ny,Nx)=(3.0d0/8.0d0*v_fake(Ny,Nx)*dy-2.0d0/8.0d0*v_fake(Ny,Nx-1)
> *dy-1.0d0/
> &
> 8.0d0*v_fake(Ny,Nx-2)*dy)*F_w_v(Ny,Nx)+(-v_J*dy)*F_e_v(Ny,Nx)
> +(3.0d0/8.0d0*&
>
> v_fake(Ny,Nx)*dx-2.0d0/8.0d0*v_fake(Ny-1,Nx)*dx-1.0d0/8.0d0*v_fake(Ny-2,Nx)
> *dx)*&
> F_s_v(Ny,Nx)+(-v_L*dx)*F_n_v(Ny,Nx)+(p_fake(Ny-1,Nx)-
> p_fake(Ny,Nx))*dx+D_J*8.0d0/&
> 3.0d0*v_J*dy+D_L*8.0d0/3.0d0*v_L*dx+((1.0d0-
> alpha_v)*b_P(Ny,Nx))*v_old(Ny,Nx)
>
> ! b_P(Ny,Nx)=1.0d0
> ! b_W(Ny,Nx)=1.0d-20
> ! b_E(Ny,Nx)=1.0d-20
> ! b_S(Ny,Nx)=1.0d-20
> ! b_N(Ny,Nx)=1.0d-20
> ! S_v_fake(Ny,Nx)=1.0d-30
>
>
>
>
>
> !coeficientes de p
> c_W(Ny,Nx)=dy/a_P(Ny,Nx)*dy
> c_E(Ny,Nx)=dy/a_E(Ny,Nx)*alpha_u*dy
> c_S(Ny,Nx)=dx/b_P(Ny,Nx)*dx
> c_N(Ny,Nx)=dx/b_N(Ny,Nx)*alpha_v*dx
> c_P(Ny,Nx)=c_W(Ny,Nx)+c_E(Ny,Nx)+c_S(Ny,Nx)+c_N(Ny,Nx)
> S_p_fake(Ny,Nx)=0.0d0 !Faltam termos de u e de v.
>
>
> !Fonte
>
> k=1
>
> do i=1,Ny
> do j=1,Nx
> S(k)=S_fake(i,j)
> S_v(k)=S_v_fake(i,j)
> S_p(k)=S_p_fake(i,j)
> k=k+1
> end do
> end do
>
>
> !Matriz dos Coeficientes
>
> k=1
>
> do i=1,Ny
> do j=1,Nx
> a(k,k)=a_P(i,j)
> b(k,k)=b_P(i,j)
> c(k,k)=c_P(i,j)
> k=k+1
> end do
> end do
>
>
> k=2
>
> do i=1,Ny
> if (i.eq.1) then
> do j=2,Nx
> a(k,k-1)=-a_W(i,j)
> b(k,k-1)=-b_W(i,j)
> c(k,k-1)=-c_W(i,j)
> k=k+1
> end do
> else
> do j=1,Nx
> a(k,k-1)=-a_W(i,j)
> b(k,k-1)=-b_W(i,j)
> c(k,k-1)=-c_W(i,j)
> k=k+1
> end do
> end if
> end do
>
> k=2
>
> do i=1,Ny
> if (i.ne.Ny) then
> do j=1,Nx
> a(k-1,k)=-a_E(i,j)
> b(k-1,k)=-b_E(i,j)
> c(k-1,k)=-c_E(i,j)
> k=k+1
> end do
> else
> do j=1,Nx-1
> a(k-1,k)=-a_E(i,j)
> b(k-1,k)=-b_E(i,j)
> c(k-1,k)=-c_E(i,j)
> k=k+1
> end do
> end if
> end do
>
> k=Nx+1
>
> do i=2,Ny
> do j=1,Nx
> a(k,k-Nx)=-a_S(i,j)
> b(k,k-Nx)=-b_S(i,j)
> c(k,k-Nx)=-c_S(i,j)
> k=k+1
> end do
> end do
>
> k=1
>
> do i=1,Ny-1
> do j=1,Nx
> a(k,k+Nx)=-a_N(i,j)
> b(k,k+Nx)=-b_N(i,j)
> c(k,k+Nx)=-c_N(i,j)
> k=k+1
> end do
> end do
>
>
>
> k=0
> l=0
> m=0
>
> do i=1,Ny*Nx
> do j=1,Ny*Nx
> if (dabs(a(i,j)).ne.0.0d0) then
> k=k+1
> end if
> if (dabs(b(i,j)).ne.0.0d0) then
> l=l+1
> end if
> if (dabs(c(i,j)).ne.0.0d0) then
> m=m+1
> end if
> end do
> end do
>
> allocate(AA(k),irow_u(k),jcol_u(k))
> allocate(BB(l),irow_v(l),jcol_v(l))
> allocate(CC(m),irow_p(m),jcol_p(m))
>
> k=1
> l=1
> m=1
>
> do i=1,Ny*Nx
> do j=1,Ny*Nx
> if (dabs(a(i,j)).ne.0.0d0) then
> AA(k)=a(i,j)
> irow_u(k)=i
> jcol_u(k)=j
> k=k+1
> end if
> if (dabs(b(i,j)).ne.0.0d0) then
> BB(l)=b(i,j)
> irow_v(l)=i
> jcol_v(l)=j
> l=l+1
> end if
> if (dabs(c(i,j)).ne.0.0d0) then
> CC(m)=c(i,j)
> irow_p(m)=i
> jcol_p(m)=j
> m=m+1
> end if
> end do
> end do
>
> call DL4LXG(iparam,rparam)
>
> iparam(5)=1000000
>
>
> !call DLINRG(Ny*Nx,a,Ny*Nx,a_inv,Ny*Nx)
> !CALL DLSARG (Ny*Nx,a,Ny*Nx,S,ipath,u)
> call DLSLXG(Ny*Nx,k-1,AA,irow_u,jcol_u,S,ipath,iparam,rparam,u)
> call DLSLXG(Ny*Nx,l-1,BB,irow_v,jcol_v,S_v,ipath,iparam,rparam,v)
> call DLSLXG(Ny*Nx,m-1,CC,irow_p,jcol_p,S_p,ipath,iparam,rparam,p)
>
> ! k=1
>
>
> ! do i=1,Ny
> ! do j=1,Nx
> ! p_fake(i,j)=p(k)
> ! if (j.ne.1) then
> ! u_fake(i,j)=u(k)+dy/(a_P(i,j)*alpha_u)*(p_fake(i,j-1)-
> p_fake(i,j))
> ! end if
> ! if (i.ne.1) then
> ! v_fake(i,j)=v(k)+dx/(b_P(i,j)*alpha_v)*(p_fake(i-1,j)-
> p_fake(i,j))
> ! end if
> ! p_fake(i,j)=p_old(i,j)+alpha_p*p(k)
> ! k=k+1
> ! end do
> ! end do
>
>
> p_old=p_fake
> u_old=u_fake
> v_old=v_fake
>
>
> k=1
>
> do i=1,Ny
> do j=1,Nx
> p_fake(i,j)=p(k)
> !p_fake(i,j)=p_old(i,j)+alpha_p*p(k)
> !p_fake(i,j)=alpha_p*(k)+(1.0d0-alpha_p)*p_old(i,j)
> !if (i.ne.1) then
> ! if (i.ne.Ny) then
> ! if (j.ne.1) then
> ! if (j.ne.Nx) then
> !u_fake(i,j)=u(k)+dy/(a_P(i,j)*alpha_u)*(p_fake(i,j-1)-
> p_fake(i,j))
> !v_fake(i,j)=v(k)+dx/(b_P(i,j)*alpha_v)*(p_fake(i-1,j)-
> p_fake(i,j))
> ! u_fake(i,j)=u(k)+dy/a_P(i,j)*(p_fake(i,j-1)-p_fake(i,j))
> ! v_fake(i,j)=v(k)+dx/b_P(i,j)*(p_fake(i-1,j)-p_fake(i,j))
> ! u_fake(i,j)=alpha_u*u(k)+(1.0d0-alpha_u)*u_old(i,j)
> ! v_fake(i,j)=alpha_v*v(k)+(1.0d0-alpha_v)*v_old(i,j)
> ! end if
> ! end if
> ! end if
> !end if
> if (j.ne.1) then
> u_fake(i,j)=u(k)+dy/(a_P(i,j)*alpha_u)*(p_fake(i,j-1)-p_fake(i,j))
> !u_fake(i,j)=u(k)+dy/a_P(i,j)*(p_fake(i,j-1)-p_fake(i,j))
> !u_fake(i,j)=alpha_u*u(k)
> !u_fake(i,j)=alpha_u*u(k)+(1.0d0-alpha_u)*u_old(i,j)
> ! write(*,*) u_fake(i,j)
> end if
> if (i.ne.1) then
> v_fake(i,j)=v(k)+dx/(b_P(i,j)*alpha_v)*(p_fake(i-1,j)-p_fake(i,j))
> !v_fake(i,j)=v(k)+dx/b_P(i,j)*(p_fake(i-1,j)-p_fake(i,j))
> !v_fake(i,j)=alpha_v*v(k)
> !v_fake(i,j)=alpha_v*v(k)+(1.0d0-alpha_v)*v_old(i,j)
> end if
> p_fake(i,j)=p_old(i,j)+alpha_p*p(k)
> k=k+1
> end do
> end do
>
>
>
>
> deallocate(AA,irow_u,jcol_u,BB,irow_v,jcol_v,CC,irow_p,jcol_p)
>
> !u=matmul(a_inv,S)
> !theta=matmul(a_inv,S)
>
>
>
> ! end do
>
> open(1,file='u.dat')
> write(1,11) Ny,Nx
> 11 format('variables=X,Y,u',/,'zone t="time','" i=',i4,' j=',i4,'
> f=point')
>
> k=1
>
> do i=1,Ny
> do j=1,Nx
> write(1,*) Y(i),X(j),u(k)
> k=k+1
> end do
> end do
>
> open(2,file='v.dat')
> write(2,12) Ny,Nx
> 12 format('variables=X,Y,v',/,'zone t="time','" i=',i4,' j=',i4,'
> f=point')
>
> k=1
>
> do i=1,Ny
> do j=1,Nx
> write(2,*) Y(i),X(j),v(k)
> k=k+1
> end do
> end do
>
> open(3,file='p.dat')
> write(3,13) Ny,Nx
> 13 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
> f=point')
>
> k=1
>
> do i=1,Ny
> do j=1,Nx
> write(3,*) Y(i),X(j),p(k)
> k=k+1
> end do
> end do
>
> open(4,file='S.dat')
> write(4,14) Ny,Nx
> 14 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
> f=point')
>
> k=1
>
> do i=1,Ny
> do j=1,Nx
> write(4,*) S_fake(i,j),k,S(k)
> k=k+1
> end do
> end do
>
> open(5,file='a.dat')
> !write(5,15) Ny,Nx
> !15 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
> f=point')
>
>
> do i=1,Ny*Nx
> do j=1,Ny*Nx
> write(5,*) i,j,a(i,j)
> end do
> end do
>
> open(6,file='aa.dat')
> !write(6,16) Ny,Nx
> !16 format('variables=X,Y,p',/,'zone t="time','" i=',i4,' j=',i4,'
> f=point')
>
>
> do i=1,Ny
> do j=1,Nx
> write(6,*) i,j,a_P(i,j)
> end do
> end do
>
>
>
> end program EXERCICIO


From: mecej4 on
Felipe Abrão wrote:

> Please, could somebody explain me what is happening. I simply
> allocated a matrix called "a", but figured out that the program is
> inserting some strange values to some of its elements (like
> 1.057229231618489E-307 or 6.941570408926931E+228), even before I set
> any value to it! Shouldn't all the values be zero before any desired
> value setting? Below is my code, but you don't have to read it all.
> You can check the problem just after line 65, where I wrote a command
> in order to read the values of matrix "a". The same thing is happening
> with matrixs "b" and "c". In advance I thank you all!
>
> PS: I use Fortran 90 in Windows XP.
<-- long program listing -->

It is your job to initialize variables before using them.

Some compilers have a switch to initialize variables by default to zero.
This, however, does not work with allocatable variables -- you cannot
initialize variables that have not yet been allocated.

-- mecej4
From: e p chandler on

"Felipe Abr�o" <felipeabrao(a)gmail.com> wrote in message
news:bf5a29d3-df12-47e5-82bf-951daffefbea(a)u7g2000yqm.googlegroups.com...
Please, could somebody explain me what is happening. I simply
allocated a matrix called "a", but figured out that the program is
inserting some strange values to some of its elements (like
1.057229231618489E-307 or 6.941570408926931E+228), even before I set
any value to it! Shouldn't all the values be zero before any desired
value setting? Below is my code, but you don't have to read it all.
You can check the problem just after line 65, where I wrote a command
in order to read the values of matrix "a". The same thing is happening
with matrixs "b" and "c". In advance I thank you all!

PS: I use Fortran 90 in Windows XP.

[really long code snipped]

------>
That is what is supposed to happen. If you don't assign a value, then you
should get whatever happens to be in memory. These small and large numbers
are typical "garbage" values. If you want zero, assign zero to it. The good
news is that you can assign a value to an entire array with one statement.

[rant mode on]
While there are some compilers that will zero memory first and there are
some compiler options that may zero memory first, you should not count on
that taking place. Un-initialized variables are a programming error.
Consider the case of a recent poster in this newsgroup. The poor soul has
inherited a code base of approximately 600,000 lines of code. His code
requires that un-initialized variables be set to zero (along with something
else not quite relevant to this discussion.). His compiler has an option to
do that. But he wants to port his code to a different compiler and possibly
different platforms. He has a monumental task which I am extremely happy is
NOT mine!
[rant mode off]
--- Elliot






From: Louis Krupp on
On 7/5/2010 1:37 PM, Felipe Abr�o wrote:
> Please, could somebody explain me what is happening. I simply
> allocated a matrix called "a", but figured out that the program is
> inserting some strange values to some of its elements (like
> 1.057229231618489E-307 or 6.941570408926931E+228), even before I set
> any value to it! Shouldn't all the values be zero before any desired
> value setting? Below is my code, but you don't have to read it all.

You don't have to post it all.

Write a small program to demonstrate the problem you're having, and post
it. People will thank you for doing that.

Louis