Sequential program

C=======================================================================
C=======================================================================
C
C       example program poisson:
C               uses Jacobi relaxation to solve the Poisson equation
C
C=======================================================================
C=======================================================================

        program poisson

C=======grid size
        integer NX
        integer NY
        parameter (NX=10)
        parameter (NY=10)

        real H
        parameter (H=0.05)
C       how often to check for convergence
        integer NCHECK
        parameter (NCHECK=10)
C       convergence criterion
        real TOL
        parameter (TOL=0.00001)
C       maximum number of steps
        integer MAXSTEPS
        parameter (MAXSTEPS=1000)

        external F
        external G

C=======grid variables
        real uk(1:NX,1:NY), ukp1(1:NX,1:NY)

        real diff, diffmax

C=======initialization
C       interior points
        do i = 2, NX-1
        do j = 2, NY-1
                uk(i,j) = F(i,j,NX,NY,H)
        enddo
        enddo
C       boundary points
        do j = 1, NY
                uk(1,j) = G(1,j,NX,NY,H)
                uk(NX,j) = G(NX,j,NX,NY,H)
        enddo
        do i = 2, NX-1
                uk(i,1) = G(i,1,NX,NY,H)
                uk(i,NY) = G(i,NY,NX,NY,H)
        enddo

C=======main loop (until convergence)
        diffmax = TOL + 1.0

        do k=1,MAXSTEPS

C=======grid computation (compute new values)
        do i = 2, NX-1
        do j = 2, NY-1
                ukp1(i,j) = 0.25*(H*H*F(i,j,NX,NY,H) 
     -                  + uk(i-1,j) + uk(i,j-1) 
     -                  + uk(i+1,j) + uk(i,j+1) )
        enddo
        enddo
C=======convergence test (global max)
        if (mod(k,NCHECK) .eq. 0) then
                diffmax = 0.0
                do i = 2, NX-1
                do j = 2, NY-1
                        diff = abs(ukp1(i,j) - uk(i,j))
                        if (diff .gt. diffmax) diffmax = diff
                enddo
                enddo
        endif
C=======grid computation (copy new values to old values)
        do i = 2, NX-1
        do j = 2, NY-1
                uk(i,j) = ukp1(i,j)
        enddo
        enddo

C=======convergence check
        if (diffmax .le. TOL) go to 1000

        enddo

C=======sequential output
1000    continue
        print*, 'NX, NY = ', NX, NY
        print*, 'H = ', H
        print*, 'tolerance = ', TOL
        call Fprint
        call Gprint

        if (diffmax .le. TOL) then
                print*, 'convergence occurred in ', k, ' steps'
        else
                print*, 'no convergence in ', MAXSTEPS, ' steps; ',
     -                  'max. difference ', diffmax
        endif

        do i = 1, NX
                if (NX .gt. 10) print*, ' '
                print 9999, (uk(i,j), j = 1, NY)
        enddo
9999    format(10F8.4)
                
        end

C=======================================================================

        real function F(i,j,nx,ny,h)
        integer i, j, nx, ny
        real h

        F = 0.0
        end

        subroutine Fprint
        print*, 'F(i,j) = 0.0'
        end

C=======================================================================

        real function G(i,j,nx,ny,h)
        integer i, j, ny, ny
        real h

        G = (i+j)*h
        end

        subroutine Gprint
        print*, 'G(i,j) = (i+j)*H'
        end

C=======================================================================