******Integrable function.
      REAL*4 FUNCTION fun(x)
      IMPLICIT NONE
      REAL*4 x
      fun = x**2
      END


******Reading parameters from a file.
      SUBROUTINE InputParams
      IMPLICIT NONE
      REAL*4 a, b, dx
      INTEGER n_steps
      COMMON/params/a, b, dx, n_steps
      OPEN(1, FILE = 'input.txt')
      READ(1, *)a, b, n_steps
      dx = (b - a)/n_steps
      CLOSE(1)
      END


******Creating a grid.
      SUBROUTINE CreateGrid(grid, n_steps)
      IMPLICIT NONE
      REAL*4 grid
      INTEGER n_steps, i
      DIMENSION grid(n_steps + 1)
      REAL*4 a, b, dx, point
      COMMON/params/a, b, dx
      n_steps = n_steps + 1
      point = a
      grid(1) = point
      DO i = 2, n_steps
        point = point + dx
        grid(i) = point
      END DO
      END


******Method of trapeziums.
      REAL*4 FUNCTION Trap(grid, n_steps)
      IMPLICIT NONE
      REAL*4 grid
      INTEGER n_steps, i
      DIMENSION grid(n_steps)
      REAL*4 a, b, dx, sum, fun
      COMMON/params/a, b, dx
      sum = 0.0
      DO i = 2, n_steps - 1
        sum = sum + fun(grid(i))
      END DO
      Trap = (sum + (fun(grid(1)) + fun(grid(n_steps)))/2.0)*dx
      END


******Gauss2 method.
      REAL*4 FUNCTION Gauss2(grid, n_steps)
      IMPLICIT NONE
      REAL*4 grid, q, x
      INTEGER n_steps, i, j
      DIMENSION grid(n_steps), x(2)
      REAL*4 sum, sum_int, a, b, h, x1, fun
      COMMON/params/a, b
      DATA x / -0.5773503, 0.5773503 /
      q = 1.0
      h = (b-a)/(n_steps - 1)
      DO i = 1,2
        sum_int = 0.0
        DO j = 2, n_steps
            x1 = (grid(j-1)+grid(j)+x(i)*h)/2
            sum_int=sum_int+h*fun(x1)
        END DO
        sum = sum + q * sum_int
      END DO
      Gauss2 = sum/2.0
      END


******Main function.
      PROGRAM LAB4
      IMPLICIT NONE
      REAL*4 array
      DIMENSION array(1000000)
      REAL*4 a, b, dx, ans, Trap, Gauss2
      INTEGER n_steps
      COMMON/params/a, b, dx, n_steps

      CALL InputParams

      CALL CreateGrid(array(1), n_steps)

      ans = Trap(array(1), n_steps)
      PRINT *,'Method of trapeziums:'
      PRINT *, ans

      ans = Gauss2(array(1), n_steps)
      PRINT *,'Gauss2 method:'
      PRINT *, ans

      PAUSE
      END