[求助]版主和高手们请帮忙看看我的程序错在哪？？！！

wlh

我的问题是：一维浅水方程
底坡项是：b(x)=0.2-0.05*(x-10)^2.当8<x<12时，其他b(x)=0
初值条件是：给定入流单宽流量hu=0.18m^2/s.出流处z=h=0.33m
我写的程序如下：请您看看怎么办啊。谢谢!!!

C        shallow equation use DG method
      program main
      implicit double precision (a-h,o-z)
      parameter (g=9.8d0,nx=200)
      double precision u0(0:2,0:nx+1,3),u1(0:2,0:nx+1,3)
      double precision u2(0:2,0:nx+1,3),right(0:2,0:nx+1,2)
      double precision x(0:nx+1),flux(0:nx+1,2),rint(0:2,0:nx+1,2)
c     u(l,i,m) l==coeficient 0,1,2 ** i==pionts 0,1,2,....,nx+1 ** m==1(h),2(h*u)
      
      write(*,*) ';Input T:';
      read(*,*) supt
      
      call initial(x,u0,nx,dx)
      t=0.                     
      do while(t.lt.supt)
        dt=getdt(u0,nx,dx,g)
        if(t+dt.gt.supt) dt=supt-t
        t=t+dt
        write(*,*) t
        call limiter(u0,nx,dx,g)
            !pause  ';2   ';
        call calright(u0,nx,flux,rint,right,g)  
        do  1 i=1,nx
          do 1 m=1,2
            u1(0,i,m)=u0(0,i,m)+dt*right(0,i,m)/dx
            u1(1,i,m)=u0(1,i,m)+dt*right(1,i,m)*3./dx
            u1(2,i,m)=u0(2,i,m)+dt*right(2,i,m)*45.d0/4./dx
    1   continue
        do  2 i=1,nx
    bminus=u0(0,i,3)+u0(1,i,3)+2.d0/3.*u0(2,i,3)
          bplus=u0(0,i+1,3)-u0(1,i+1,3)+2.d0/3.*u0(2,i+1,3)
          do 2 m=1,2
            u1(0,i,m)=u1(0,i,m)+0.25*dt*(bplus+bminus)*(u1(0,i,m)+
     %         u0(0,i,m))/dx
            u1(1,i,m)=u1(1,i,m)+0.25*dt*(bplus+bminus)*(u1(1,i,m)+
     %         u0(1,i,m))*3./dx
            u1(2,i,m)=u1(2,i,m)+0.25*dt*(bplus+bminus)*(u1(2,i,m)+
     %         u0(2,i,m))*45.d0/4./dx
    2   continue         
        call boundary(u1,nx)
             !  pause ';3';
        call limiter(u1,nx,dx,g)                                      
        call calright(u1,nx,flux,rint,right,g)
        do  3 i=1,nx
          do 3 m=1,2
            u2(0,i,m)=0.75*u0(0,i,m)+0.25*(u1(0,i,m)+dt*right(0,i,m)/dx)
            u2(1,i,m)=0.75*u0(1,i,m)
     *               +0.25*(u1(1,i,m)+dt*right(1,i,m)*3./dx)
            u2(2,i,m)=0.75*u0(2,i,m)
     *               +0.25*(u1(2,i,m)+dt*right(2,i,m)*45.d0/4./dx)  
    3   continue
       do  4 i=1,nx
    bminus=u1(0,i,3)+u1(1,i,3)+2.d0/3.*u1(2,i,3)
          bplus=u1(0,i+1,3)-u1(1,i+1,3)+2.d0/3.*u1(2,i+1,3)
          do 4 m=1,2
            u2(0,i,m)=u2(0,i,m)+0.25*dt*(bplus+bminus)*(u1(0,i,m)+
     %         u2(0,i,m))/dx
            u2(1,i,m)=u2(1,i,m)+0.25*dt*(bplus+bminus)*(u1(1,i,m)+
     %         u2(1,i,m))*3./dx
            u2(2,i,m)=u2(2,i,m)+0.25*dt*(bplus+bminus)*(u1(2,i,m)+
     %         u2(2,i,m))*45.d0/4./dx
    4   continue
        call boundary(u2,nx)
        
             !        pause';4';
        call limiter(u2,nx,dx,g)   
        call calright(u2,nx,flux,rint,right,g)
        do  5 i=1,nx
          do 5 m=1,2
            u0(0,i,m)=1.d0/3.*u0(0,i,m)
     *               +2.d0/3.*(u2(0,i,m)+dt*right(0,i,m)/dx)
            u0(1,i,m)=1.d0/3.*u0(1,i,m)
     *               +2.d0/3.*(u2(1,i,m)+dt*right(1,i,m)*3./dx)
            u0(2,i,m)=1.d0/3.*u0(2,i,m)
     *               +2.d0/3.*(u2(2,i,m)+dt*right(2,i,m)*45.d0/4./dx)  
    5   continue
         do  6 i=1,nx
    bminus=u2(0,i,3)+u2(1,i,3)+2.d0/3.*u2(2,i,3)
          bplus=u2(0,i+1,3)-u2(1,i+1,3)+2.d0/3.*u2(2,i+1,3)
          do 6 m=1,2
            u0(0,i,m)=u0(0,i,m)+0.25*dt*(bplus+bminus)*(u2(0,i,m)+
     %         u0(0,i,m))/dx
            u0(1,i,m)=u0(1,i,m)+0.25*dt*(bplus+bminus)*(u2(1,i,m)+
     %         u0(1,i,m))*3./dx
            u0(2,i,m)=u0(2,i,m)+0.25*dt*(bplus+bminus)*(u2(2,i,m)+
     %         u0(2,i,m))*45.d0/4./dx
    6   continue
        call boundary(u0,nx)
              !               pause';5';
      enddo      
      
      write(*,*) ';T=';,t
      call output(x,u0,nx)
      end
      subroutine boundary(u0,nx)
      implicit double precision (a-h,o-z)
      double precision u0(0:2,0:nx+1,3)
      
        h=?          !! 上游水深咋办？
        u=0.18
        u0(0,0,1)=h
        u0(0,0,2)=u  
        h=0.33
        u=?          !!!!!下游流量咋办啊
        u0(0,nx+1,1)=h
        u0(0,nx+1,2)=u
        do l=1,2
          do m=1,2
            u0(l,0,m)=0.0
            u0(l,nx+1,m)=0.0
          enddo
        enddo  
      end
      
      subroutine initial(x,u0,nx,dx)
      implicit double precision (a-h,o-z)
      double precision u0(0:2,0:nx+1,3)
      double precision x(0:nx+1)
      
      xl=0
      xr=25.
      dx=(xr-xl)/nx
      do i=0,nx+1
        x(i)=i*dx
do  l=0,2
  if(x(i).gt.8.and.x(i).lt.12)then
  u0(l,i,3)=0.2-0.05*(x(i)-10)**2
  else
  u0(l,i,3)=0
  endif
enddo
      enddo  
      do 10 i=0,nx+1  
          h=0.33-u0(0,i,3)   !!!!!这里的设置对不对
    u=0
   
        u0(0,i,1)=h
        u0(0,i,2)=u

        do 10 l=1,2
          do 10 m=1,2
            u0(l,i,m)=0.0
   10   continue
      end
      
      double precision function getdt(u,nx,dx,g)
      implicit double precision (a-h,o-z)
      double precision u(0:2,0:nx+1,2)
      a=0.0
      do i=1,nx
        velocity=u(0,i,2)
        c=sqrt(g*u(0,i,1))
        temp=abs(velocity)+c
        if(temp.gt.a) a=temp
      enddo
      getdt=dx/a*0.08
      end                     
      
C     minmod function
      double precision function rminmod(x,y,z,dx)
      implicit double precision (a-h,o-z)
      rm=0.0
      if(abs(x).le.rm*dx*dx) then
        rminmod=x
      else if(x.gt.0.and.y.gt.0.and.z.gt.0) then
        rminmod=min(x,y,z)
      else if(x.lt.0.and.y.lt.0.and.z.lt.0) then
        rminmod=max(x,y,z)
      else
        rminmod=0
      endif        
      end
C          v=Ru
      subroutine product(r,u,v)
      implicit double precision (a-h,o-z)
      double precision r(2,2),u(2),v(2)
      do 50 m=1,2
        v(m)=0.
        do 50 m1=1,2
          v(m)=v(m)+r(m,m1)*u(m1)
   50 continue
      end               
C     limiter  characterized piecewise  
      subroutine limiter(u,nx,dx,g)
      implicit double precision (a-h,o-z)
      double precision u(0:2,0:nx+1,3),evecr(2,2),evecl(2,2)
      double precision v_1(2),v0(2),vx(2),up(2)
      do 20 i=1,nx            
        c=sqrt(u(0,i,1)/g)
        evecl(1,1)=1.
        evecl(1,2)=c
        evecl(2,1)=1.
        evecl(2,2)=-c
        evecr(1,1)=0.5
        evecr(1,2)=0.5
        evecr(2,1)=0.5/c
        evecr(2,2)=-0.5/c
        do mm=1,2
          up(mm)=u(0,i,mm)-u(0,i-1,mm)
        enddo  
        call product(evecl,up,v_1)
        do mm=1,2
          up(mm)=u(0,i+1,mm)-u(0,i,mm)
        enddo  
        call product(evecl,up,v0)
        do mm=1,2
          up(mm)=u(1,i,mm)
        enddo  
        call product(evecl,up,vx)
        
        do mm=1,2
          vx(mm)=rminmod(vx(mm),v_1(mm),v0(mm),dx)
        enddo
        call product(evecr,vx,up)
        do mm=1,2
          if(abs(u(1,i,mm)-up(mm)).gt.1.e-7) then
            u(1,i,mm)=up(mm)
            u(2,i,mm)=0.
          endif  
        enddo      
   20 continue
      call boundary(u,nx)
      end  
C         calculate the numerical flux
      subroutine calflux(u,nx,flux,g)
      implicit double precision (a-h,o-z)
      double precision u(0:2,0:nx+1,3),flux(0:nx+1,3)
      do i=0,nx
        hminus=u(0,i,1)+u(1,i,1)+2.d0/3.*u(2,i,1)
        hplus=u(0,i+1,1)-u(1,i+1,1)+2.d0/3.*u(2,i+1,1)   
        vminus=u(0,i,2)+u(1,i,2)+2.d0/3.*u(2,i,2)
        vplus=u(0,i+1,2)-u(1,i+1,2)+2.d0/3.*u(2,i+1,2)
  bminus=u(0,i,3)+u(1,i,3)+2.d0/3.*u(2,i,3)
        bplus=u(0,i+1,3)-u(1,i+1,3)+2.d0/3.*u(2,i+1,3)
  v1=abs(u(0,i,2))
        v2=abs(u(0,i+1,2))
        c1=sqrt(g*u(0,i,1))
        c2=sqrt(g*u(0,i+1,1))
        alpha=max(v1+c1,v2+c2)
        fplus=hplus*vplus
        fminus=hminus*vminus
        flux(i,1)=0.5*(fplus+fminus-alpha*(hplus-hminus))
        fplus=0.5*vplus*vplus+g*(hplus+bplus)!!!!!!把这种底坡加在数值
        fminus=0.5*vminus*vminus+g*(hminus+bminus)!!!!!通量中行不？？
        flux(i,2)=0.5*(fplus+fminus-alpha*(vplus-vminus))
      enddo
      end  
C     calculate the integrate term
      subroutine calint(u,nx,rint,g)
      implicit double precision (a-h,o-z)
      double precision u(0:2,0:nx+1,3),rint(0:2,0:nx+1,2)
      gd_1=-sqrt(0.6)
      gd0=0.0
      gd1=sqrt(0.6)
      do i=1,nx
        do m=1,2
          rint(0,i,m)=0.
        enddo            
        h_1=u(0,i,1)+u(1,i,1)*gd_1+u(2,i,1)*(gd_1**2-1.d0/3.)
        h0=u(0,i,1)+u(1,i,1)*gd0+u(2,i,1)*(gd0**2-1.d0/3.)
        h1=u(0,i,1)+u(1,i,1)*gd1+u(2,i,1)*(gd1**2-1.d0/3.)
        v_1=u(0,i,2)+u(1,i,2)*gd_1+u(2,i,2)*(gd_1**2-1.d0/3.)
        v0=u(0,i,2)+u(1,i,2)*gd0+u(2,i,2)*(gd0**2-1.d0/3.)
        v1=u(0,i,2)+u(1,i,2)*gd1+u(2,i,2)*(gd1**2-1.d0/3.)
        hv_1=h_1*v_1
        hv0=h0*v0
        hv1=h1*v1
        rint(1,i,1)=5.d0/9.*hv_1+8.d0/9.*hv0+5.d0/9.*hv1
        rint(1,i,2)=5.d0/9.*(0.5*v_1**2+g*h_1)+8.d0/9.*(0.5*v0**2
     %           +g*h0)+5.d0/9.*(0.5*v1**2+g*h1)
        rint(2,i,1)=2.*(5.d0/9.*hv_1*gd_1
     %                 +8.d0/9.*hv0*gd0+5.d0/9.*hv1*gd1)
        rint(2,i,2)=2.*(5.d0/9.*(0.5*v_1**2+g*h_1)*gd_1+8.d0/9.
     %    *(0.5*v0**2+g*h0)*gd0+5.d0/9.*(0.5*v1**2+g*h1)*gd1)
      enddo
      end        
C     calculate the right term
      subroutine calright(u,nx,flux,rint,right,g)
      implicit double precision (a-h,o-z)
      double precision u(0:2,0:nx+1,3),right(0:2,0:nx+1,2)
      double precision flux(0:nx+1,2),rint(0:2,0:nx+1,2)
      
      call calflux(u,nx,flux,g)
      call calint(u,nx,rint,g)
      do 30 i=1,nx                              
        do 30 m=1,2
          right(0,i,m)=-(flux(i,m)-flux(i-1,m))
          right(1,i,m)=-(flux(i,m)+flux(i-1,m))+rint(1,i,m)
          right(2,i,m)=-2.d0/3.*(flux(i,m)-flux(i-1,m))+rint(2,i,m)
   30 continue      
      end
C        output x,u
      subroutine output(x,u,nx)
      implicit double precision (a-h,o-z)
      double precision x(0:nx+1),u(0:2,0:nx+1,3),z(0:nx+1)
      open(unit=1,file=';hight.plt';)
      open(unit=2,file=';veloc.plt';)
      open(unit=3,file=';b.plt';)
      do i=1,nx
        write(1,101) x(i),u(0,i,1)+u(1,i,3)
        write(2,101) x(i),u(0,i,2)
  write(3,101) x(i),u(1,i,3)
      enddo
  101 format(1x,2e20.10)
      close(1)
close(2)
close(3)
      end