Résoudre l'équation d'advection avec terme source par volumes finis

[turbulent]

Bonjour tout le monde;
J'ai conçu un code en fortran pour résoudre l'équation d'advection avec terme source par la méthode de volume finis : Ct + F(C)x = G(C)
le code n'affiche aucune erreur, mais il ne tourne pas, aidez s'il vous plait à détecter l'erreur.

      implicit none
      real ros,rof,lambda,beta,
     .k0,dt,dx,dh,l,d,mu,nc 
      real q(2000,110)
	real phi(2000,110)
	real c(2000,110)
	real cm(2000,110)
	real m(2000,110)
      real F(2000,110)	
	real G(2000,110)	
c   	real a
	real x(110)
	real xc(110)
	integer i,j,a
      
	l=0.03d0
	nc=32
	mu=0.003005d0
	d=0.007978846d0
	
	
	dt=0.5d0
	dx=l/(nc-2)
	dh=20.d0
	ros=2671
	rof=1000
	k0=0.00003d0
	lambda=0.003d0
		

	
	do  j=2,nc-1
	phi(1,j)=0.33d0

	c(1,j)=0.d0
	m(1,j)=0.d0
	
	enddo	

	do  i=1,2000
	
	phi(i,1)=1.d0
      phi(i,nc)=1.d0
	c(i,1)=0.d0
c	phi(i,22

	c(i,nc)=c(i,nc-1)
      c(i,nc+1)=c(i,nc)
c      c(i,nc+1)=0.1d0
	phi(i,nc+1)=1.d0
C	q(i,nc+1)=q(i,nc)
	q(i,nc+1)=0.5
C	c(i,nc+2)=c(i,nc)
C	phi(i,nc+2)=1.d0
c	q(i,nc+2)=q(i,nc)
C	q(i,nc+2)=0.5
	enddo	
	
	do i=2,2000
	do j=2,nc 
	x(j)=j*dx
c	xc(j)=dx/2+(j-1)*dx
      do a=j,j+1
	xc(a)=dx/2+(a-1)*dx
c      x(j)=j*dx
c	x(j-0.5)=x(j-1)+dx/2
      q(i-1,j) =k0*((ros/rof-1)*c(i-1,j)+1)*((1-phi(1,j))**2)*(phi
     .(i-1,j)**3)*dh/((1+2.5*c(i-1,j))*(phi(1,j)**3)*(1-phi(i-1,j))**2) 
      	
	phi(i,j)=phi(i-1,j)+dt* lambda *(1-phi(i-1,j)) * q(i-1,j)
      if (phi(i,j).gt.0.33005d0) then
      phi(i,j)=0.33005d0
  	endif
c      c(i,j) = 1-phi(i-1,j)*(1-c(i-1,j))/phi(i,j) - q(i-1,j)
c     .*dt *(c(i-1,j)-c(i-1,j-1))/(phi(i,j) * dx) 
      
c	a=dx/2,l-dx/2 
c	if (x(j).gt.xc(a))      then
	if (x(j).gt.xc(a))      then
c	F(i,xc(a))=q(i,j)*c(i,j)/phi(i,j) 
	F(i,x(j))=q(i,j)*c(i,j)/phi(i,j) 
	else
c	F(i,xc(a))=q(i,j-1)*c(i,j-1)/phi(i,j-1)
	F(i,x(j))=q(i,j-1)*c(i,j-1)/phi(i,j-1)
	endif
c	G(i,j)=(1-(c(i,j)-dt*(F(i,xc(a+1))-F(i,xc(a))))/dx)/phi(i,j)+q(i,j)
c     .*(c(i,j)-dt*(F(i,xc(a+1))-F(i,xc(a)))/dx)*(1/phi(i,j+1)-1/phi(i,j)
c     .)/dx

	G(i,j)=(1-(c(i,j)-dt*(F(i,x(j+1))-F(i,x(j))))/dx)/phi(i,j)+q(i,j)
     .*(c(i,j)-dt*(F(i,x(j+1))-F(i,x(j)))/dx)*(1/phi(i,j+1)-1/phi(i,j)
     .)/dx 
	
	c(i+1,j) =c(i,j)-dt*(F(i,xc(j+1))-F(i,xc(j)))/dx+dt*G(i,j)  
c	c(i+1,j) =c(i,j)-dt*(F(i,xc(a+1))-F(i,xc(a)))/dx+dt*G(i,j)  
			
	cm(i,j)=1000*ros*c(i,j)/((1-c(i,j))*rof+ros*c(i,j))
      
	m(i,j)=m(i-1,j)+1000*(ros*c(i-1,j)*q(i-1,j)*dt)
	enddo
	enddo
	enddo
	
	open(17,file='VF.txt', status='unknown')
	write(17,10)(i*dt,phi(i,nc),cm(i,nc),m(i,nc),i=1,2000)
10    format(/(4G16.6))
	end

-----

[Jack]

Je t'invite à (re)lire les règles participatives de ce forum:
http://forums.futura-sciences.com/programmation-langages-algorithmique/441632-regles-participatives-nouveau-forum.html

Comme c'est ta première intervention, j'ai ajouté les balises code

PS : houla, l'indentation est à revoir complètement. Un problème d'éditeur?