1        program diff
  2        integer xmax,ymax
  3  c     Размер сетки по оси X
  4        parameter (xmax=4000)
  5  c     Размер сетки по оси Y
  6        parameter (ymax=6000)
  7        real*8 f0(xmax,ymax), f1(xmax,ymax)
  8        real*8 df
  9        integer x,y,n
 10  c     Задаем начальные значения
 11        call initial(f0,xmax,ymax)
 12  c     Задаем граничные значения
 13        call boundary(f0,xmax,ymax)
 14  c     Начинаем основной итерационный цикл
 15        df=0.0
 16        n=0
 17      1 continue
 18  c     При смене итераций меняем местами массивы
 19        if ( ((n/2)*2) .EQ. n ) then
 20  c       Обрабатываем f0 для каждой четной итерации
 21          call iter(f0,f1,xmax,ymax,df)
 22        else
 23  c       Обрабатываем f1 для каждой нечетной итерации
 24          call iter(f1,f0,xmax,ymax,df)
 25        endif
 26  c     Увеличиваем номер цикла на единицу
 27        n=n+1
 28  c     Выводим на экран разницу между итерациями
 29        write(*,*) 'Diff:',df
 30  c     Если разница больше заданной, то решение не найдено
 31  c     и мы идем на следующий цикл
 32        if (df > 1e-2) goto 1
 33  c     Печатаем количество итераций, потребовавшихся
 34  c     для нахождения решения
 35        write(*,*) 'Iteration count:', n
 36        stop
 37        end
 38  
 39  c     Подпрограмма вычисления искомой функции
 40        subroutine iter(f0,f1,xmax,ymax,df)
 41        integer xmax,ymax
 42        real*8 f0(xmax,ymax), f1(xmax,ymax)
 43        real*8 dt,dx,dy
 44        real*8 df
 45        integer x,y,n
 46        dt=0.01
 47        dx=0.5
 48        dy=0.5
 49  c     Вычисляем функцию в ячейках сетки
 50        do x=2,xmax-1
 51          do y=2,ymax-1
 52          f1(x,y)=f0(x,y)+dt*(
 53       x    (f0(x+1,y)-2*f0(x,y)+f0(x-1,y))/(dx*dx)
 54       x    +
 55       x    (f0(x,y+1)-2*f0(x,y)+f0(x,y-1))/(dy*dy)
 56       x    )
 57          end do
 58        end do
 59  c     Копируем границу для следующей итерации
 60        do x=1,xmax
 61          f1(x,1)=f0(x,1)
 62          f1(x,ymax)=f0(x,ymax)
 63        end do
 64        do y=1,ymax
 65          f1(1,y)=f0(1,y)
 66          f1(xmax,y)=f0(xmax,y)
 67        end do
 68  c     Находим максимальную дельту
 69        df=0.0
 70        do x=1,xmax
 71          do y=1,ymax
 72            if ( df < abs(f0(x,y)-f1(x,y)) ) then
 73                 df = abs(f0(x,y)-f1(x,y))
 74            endif
 75          end do
 76        end do
 77        return
 78        end
 79        
 80  c     Подпрограмма задания начальных значений
 81        subroutine initial(f0,xmax,ymax)
 82        integer xmax,ymax,x,y
 83        real*8 f0(xmax,ymax)
 84  c     Задаем начальные значения массива
 85        do x=1,xmax
 86          do y=1,ymax
 87            f0(x,y)=0.0
 88          end do
 89        end do
 90        return
 91        end
 92  
 93  c     Подпрограмма задания граничных условий
 94        subroutine boundary(f0,xmax,ymax)
 95        integer xmax,ymax,x,y
 96        real*8 f0(xmax,ymax)
 97        common myid,np
 98  c     Задаем граничные значения массива
 99        do y=1,ymax
100          f0(1,y)=0.0
101          f0(xmax,y)=0.0
102        end do
103        do x=1,xmax
104          f0(x,1)=sin((x*1.0)/(xmax/2))
105   	     f0(x,ymax)=0.0
106        end do
107        f0(xmax/2,ymax)=-5.0
108        return
109        end