请高手帮我看个一维欧拉方程Roe格式的程序

昭明

这个程序用于求解一维欧拉方程，用Roe格式写成，检查了很多遍，死活搞不定，请高人指点，不胜感激。
!        1d-euler
        parameter(Nx=50,gama=1.4,eps=0.0001)
        real,dimension(3,-Nx:Nx):: q,qt,f,qn
        real,dimension(-Nx:Nx)::rou,u,E,p,h
        real,dimension(3,-Nx+1:Nx)::flux,r1,r2,r3
        dt=0.01
        dx=0.1

        open(1,file='tec.plt')
!        write(1,*)'filetype=solution'
        write(1,*) 'variables= "x","r","u","p"'
!----initial

        do i=-Nx,0
        rou(i)=1.;        u(i)=0.;        p(i)=1.
        enddo
        do i=1,Nx
        rou(i)=0.125        ;        u(i)=0.        ;        p(i)=0.1
        enddo

        do i=-Nx,Nx
        q(1,i)=rou(i)
        q(2,i)=rou(i)*u(i)
        q(3,i)= p(i)/(gama-1.) + rou(i)*u(i)**2./2.
        enddo

!-----roe---
        do k=1,20 !时间循环，时间一长就出错。

        do i=-Nx+1,Nx-1
        qn(:,i)=q(:,i)
        enddo

        do i=-Nx,Nx
       
        rou(i)=q(1,i)
        u(i)=q(2,i)/q(1,i)
        E(i)=q(3,i)/q(1,i)
        p(i)=(gama-1.)*rou(i)*(E(i)-0.5*u(i)**2.) !p=(r-1.)*rou*(E-.5*u**2)
        h(i)=E(i)+p(i)/rou(i)
        h(i)=0.5*u(i)**2.+gama/(gama-1.)*p(i)/rou(i)

        f(1,i)=q(2,i)
        f(2,i)=q(2,i)**2./q(1,i)+p(i)
        f(3,i)=q(2,i)*q(3,i)/q(1,i) + p(i)*q(2,i)/q(1,i)
       
!        f(1,i)=rou(i)*u(i)
!        f(2,i)=rou(i)*u(i)**2.+p(i)
!        f(3,i)=u(i)*p(i)*gama/(gama-1.) + 0.5*rou(i)*u(i)**3.
        if(p(i)/rou(i)<0)then
            write(*,*) "error, N=", i
                endif
        enddo

!print*,q(1,0)
        do i=-Nx+1,Nx
!各变量的roe平均
        t1=sqrt(q(1,i))        ;        t2=sqrt(q(1,i-1))
        ra=t1*t2 !密度
        ua=( u(i)*t1+u(i-1)*t2 )/(t1+t2) !速度
        ha=( t1*h(i)+t2*h(i-1) )/( t1+t2 ) !焓
        ca=sqrt( (gama-1.)*(ha-0.5*ua**2.) ) !声速
        pa=(ra*ha-ra*ua**2./2) * (gama-1.)/gama !压强
        pa=ra*ca**2./gama
               
        lad1=ua-ca                       
        lad2=ua                       
        lad3=ua+ca
        if(abs(lad1)<=eps) lad1=(lad1**2.+eps**2.)/2./eps !熵修正
        if(abs(lad2)<=eps) lad2=(lad2**2.+eps**2.)/2./eps
        if(abs(lad3)<=eps) lad3=(lad3**2.+eps**2.)/2./eps

        r1(1,i)=1.                        ;        r2(1,i)=1.                        ;        r3(1,i)=1.  !右行特征向量
        r1(2,i)=ua-ca                ;        r2(2,i)=ua                        ;        r3(2,i)=ua+ca
        r1(3,i)=ha-ua*ca        ;        r2(3,i)=0.5*ua**2.        ;        r3(3,i)=ha+ua*ca
       
        dr=rou(i)-rou(i-1)        ;        du=u(i)-u(i-1)  ;        dp=p(i)-p(i-1)

        a1=(dp-pa*ca*du)/2./ca**2.        ;        a2=dr-dp/ca**2. ;        a3=(dp+pa*ca*du)/2./ca**2.
!        a3=(dp+ra*ca*du)/2./ca**2.
       
        flux(:,i)=0.5*( f(:,i-1)+f(:,i) ) - 0.5*( abs(lad1)*a1*r1(:,i) + abs(lad2)*a2*r2(:,i) + abs(lad3)*a3*r3(:,i) )

        enddo
!---comput the n+1
        !R-K
        do i=-Nx+1,Nx-1
        q(:,i)=q(:,i)-( flux(:,i+1)-flux(:,i) )*dt/dx
        enddo
        do i=-Nx+1,Nx-1
        q(:,i)=0.75*qn(:,i)+0.25*q(:,i)-0.25*( flux(:,i+1)-flux(:,i) )*dt/dx
        enddo
        do i=-Nx+1,Nx-1
        q(:,i)=1.*qn(:,i)/3.+2.*q(:,i)/3-2./3.*( flux(:,i+1)-flux(:,i) )*dt/dx
        enddo
!--------boundary condition
!        q(:,-Nx)=q(:,-Nx+1)        ;        q(:,Nx)=q(:,Nx-1)

        print*,k
!        write(1,*) 'ZONE T="STP:',k*dt,'sec", STRANDID=1,SOLUTIONTIME=',k
!write(1,*) 'zone f=point,i=',2*nx+1,'   strandid=1,solutiontime=',k
        do i=-Nx,Nx
        write(1,*) i*dx,q(1,i),q(2,i)/q(1,i),(gama-1.)*(q(3,i)-0.5*q(2,i)**2./q(1,i))
        enddo



        enddo

        end

shirazbj

学习一下