编写的程序，用于解决差分的sor法，不知道到底哪里错了，求高人指点，感激不尽~~~~~~

program example

    implicit none
    real , parameter :: jmax=11,pi=3.1415926
    real , dimension (jmax,jmax):: A,B,C,D,E,F,U
    integer :: i,j
    real :: midl,rjac
    do i=1,jmax
        do j=1,jmax
        A(i,j)=1.0
        B(i,j)=1.0
        C(i,j)=1.0
        D(i,j)=1.0
        E(i,j)=-4.0
        F(i,j)=0.0
        U(i,j)=0.0
        end do
    end do
    midl=jmax/2+1
    F(midl,midl)=2.0
    rjac=cos(pi/jmax)
    call sor(A,B,C,D,E,F,U,jmax,rjac)
    print 100,'SOR solution:'
    100 format(1x,a)
    do i=1,jmax
        print 200,(u(i,j),j=1,jmax)
        200 format(1x,11f6.2)
    end do
    print *,'Test that sulotion satisfies Difference Eqns:'
    do i=2,jmax-1
        do j=2,jmax-1
            f(i,j)=u(i+1,j)+u(i-1,j)+u(i,j+1)&
                   +u(i,j-1)-4.0*u(i,j)
        end do
        print 300,(f(i,j),j=2,jmax-1)
        300 format(7x,11f6.2)
     end do
contains
       subroutine sor(a,b,c,d,e,f,u,jmax,rjac)
            implicit none
            real , intent(in) :: rjac
            integer,intent(in) :: jmax
            real , intent(in),dimension(jmax,jmax) :: a,b,c,d,e,f
            real , intent(inout),dimension(jmax,jmax) :: u
            real , parameter :: maxits=100,eps=1.d-5,zero=0.d0,half=.5d0,&
                       qtr=0.25d0,one=1.d0
            real :: anormf,anorm,omega,resid
            integer :: j,l,n
                anormf=zero
                do j=2,jmax-1
                    do l=2,jmax-1
                        anormf=anormf+abs(f(j,l))
                    end do
                end do
                omega=one
                do n=1,maxits
                    anorm=zero
                    do j=2,jmax-1
                        do l=2,jmax-1
                            if (mod(l+j,2)==mod(n,2)) then
                                resid=a(j,l)*u(j+1,l)+b(j,l)*u(j-1,l)&
                                      +c(j,l)*u(j,l+1)+d(j,l)*u(j,l-1)&
                                      +e(j,l)*u(j,l)-f(j,l)
                                anorm=anorm+abs(resid)
                                u(j,l)=u(j,l)-omega*resid/e(j,l)
                            end if
                        end do
                    end do
                    if(n==1)then
                        omega=one/(one-half*rjac**2)
                    else
                        omega=one/(one-qtr*rjac**2*omega)
                    end if
                    if((n>1).and.(anorm<eps*anormf)) return
                end do
                pause 'maxits exceeded'
            end subroutine sor

end  program example