    MODULE nag_pde_ell_mg_ex02_mod

      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Default Accessibility ..
      PUBLIC
      ! .. Intrinsic Functions ..
      INTRINSIC KIND
      ! .. Parameters ..
      INTEGER, PARAMETER :: wp = KIND(1.0D0)

    CONTAINS

      SUBROUTINE pde_coeff(x,y,coeff,i_comm,r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC COS, SIN
        ! .. Scalar Arguments ..
        REAL (wp), INTENT (IN) :: x, y
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), INTENT (OUT) :: coeff(7)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:)
        ! .. Executable Statements ..
        coeff(1:6) = r_comm(1:6)

        ! PSI = (-ALPHA-GAMMA+PHI)*SIN(X)*SIN(Y) + BETA*COS(X)*COS(Y) +
        ! .      + DELTA*COS(X)*SIN(Y) + EPSLON*SIN(X)*COS(Y)

        coeff(7) = (-coeff(1)-coeff(3)+coeff(6))*SIN(x)*SIN(y) + &
         coeff(2)*COS(x)*COS(y) + coeff(4)*COS(x)*SIN(y) + &
         coeff(5)*SIN(x)*COS(y)

      END SUBROUTINE pde_coeff

      SUBROUTINE bound_cond(x,y,a,b,c,bnd,i_comm,r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC SIN
        ! .. Parameters ..
        INTEGER, PARAMETER :: bottom = 0, left = 3, right = 1, top = 2
        REAL (wp), PARAMETER :: one = 1.0_wp
        REAL (wp), PARAMETER :: zero = 0.0_wp
        ! .. Scalar Arguments ..
        INTEGER, INTENT (IN) :: bnd
        REAL (wp), INTENT (OUT) :: a, b, c
        REAL (wp), INTENT (IN) :: x, y
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:)
        ! .. Executable Statements ..

        IF (bnd==top .OR. bnd==right) THEN
          ! Solution prescribed
          a = one
          b = zero
          c = SIN(x)*SIN(y)
        ELSE IF (bnd==bottom) THEN
          ! Derivative prescribed
          a = zero
          b = one
          c = -SIN(x)
        ELSE IF (bnd==left) THEN
          ! Derivative prescribed
          a = zero
          b = one
          c = -SIN(y)
        END IF

      END SUBROUTINE bound_cond

    END MODULE nag_pde_ell_mg_ex02_mod
    PROGRAM nag_pde_ell_mg_ex02

      ! Example Program Text for nag_pde_ell_mg
      ! NAG fl90, Release 4. NAG Copyright 2000.

      ! .. Use Statements ..
      USE nag_examples_io, ONLY : nag_std_out
      USE nag_pde_ell_mg, ONLY : nag_pde_ell_mg_sol, nag_pde_ell_rect
      USE nag_pde_ell_mg_ex02_mod, ONLY : pde_coeff, bound_cond, wp
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC EPSILON, MAX, REAL, SIN, SQRT
      ! .. Parameters ..
      INTEGER, PARAMETER :: levels = 3
      INTEGER, PARAMETER :: nx = 2**levels + 1
      INTEGER, PARAMETER :: ny = nx
      INTEGER, PARAMETER :: nxy = nx*ny
      REAL (wp), PARAMETER :: fifty = 50.0_wp
      REAL (wp), PARAMETER :: one = 1.0_wp
      REAL (wp), PARAMETER :: zero = 0.0_wp
      ! .. Local Scalars ..
      INTEGER :: i, j, k1, k2, max_iter, num_iter
      REAL (wp) :: acc, hx, hy, rms_err, x_max, x_min, y_max, y_min
      LOGICAL :: upwind_diff
      ! .. Local Arrays ..
      REAL (wp) :: a(nxy,7), rhs(nxy), r_comm(6), u(nxy), x(nxy), y(nxy)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) 'Example Program Results for nag_pde_ell_mg_sol'
      WRITE (nag_std_out,*) ' '
      ! r_comm(1:6) contains the coefficient alpha, beta, gamma, delta,
      ! epsilon and phi appearing in the example PDE.
      ! They are stored for the use in subroutine pde_coeff
      r_comm(1) = one
      r_comm(2) = zero
      r_comm(3) = one
      r_comm(4) = fifty
      r_comm(5) = fifty
      r_comm(6) = zero

      x_min = zero
      x_max = one
      y_min = zero
      y_max = one
      hx = (x_max-x_min)/REAL(nx-1,kind=wp)
      hy = (y_max-y_min)/REAL(ny-1,kind=wp)

      y(1:nx) = y_min
      x(1:nx) = x_min + hx*(/ (REAL(i-1,kind=wp),i=1,nx) /)
      x(nx) = x_max

      y(nxy+1-nx:nxy) = y_max
      x(nxy+1-nx:nxy) = x(1:nx)

      DO j = 2, ny - 1
        k2 = j*nx
        k1 = k2 + 1 - nx
        y(k1:k2) = y_min + REAL(j-1,kind=wp)*hy
        x(k1:k2) = x(1:nx)
      END DO

      ! Discretize the equations :
      ! Set-up operator, right-hand side

      CALL nag_pde_ell_rect(pde_coeff,bound_cond,nx,ny,x_min,x_max,y_min, &
       y_max,a,rhs,r_comm=r_comm)

      ! Set-up initial guess

      u(1:nxy) = zero

      ! Solve the equation

      acc = MAX(1.0E-6_wp,SQRT(nx*EPSILON(zero)))
      max_iter = 50
      CALL nag_pde_ell_mg_sol(nx,ny,a,rhs,u,acc=acc,max_iter=max_iter, &
       num_iter=num_iter)

      ! Print out the Solution

      WRITE (nag_std_out,*) ' '
      WRITE (nag_std_out,*) ' Exact solution above computed solution'
      WRITE (nag_std_out,*) ' '
      WRITE (nag_std_out,999) '  I/J', (i,i=1,nx)
      rms_err = zero
      DO j = ny, 1, -1
        WRITE (nag_std_out,*) ' '
        WRITE (nag_std_out,998) j, (SIN(x(i+(j-1)*nx))*SIN(y(i+(j-1)*nx)),i=1, &
         nx)
        WRITE (nag_std_out,998) j, (u(i+(j-1)*nx),i=1,nx)
        DO i = 1, nx
          rms_err = rms_err + (SIN(x(i+(j-1)*nx))*SIN(y(i+ &
           (j-1)*nx))-u(i+(j-1)*nx))**2
        END DO
      END DO
      rms_err = SQRT(rms_err/REAL(nx*ny,kind=wp))
      WRITE (nag_std_out,*) ' '
      WRITE (nag_std_out,997) ' Number of Iteration ', num_iter
      WRITE (nag_std_out,996) ' RMS Error ', rms_err

      ! Now discretize and solve the equations using upwinding differences
      ! Set-up operator, right-hand side

      upwind_diff = .TRUE.
      CALL nag_pde_ell_rect(pde_coeff,bound_cond,nx,ny,x_min,x_max,y_min, &
       y_max,a,rhs,r_comm=r_comm,upwind_diff=upwind_diff)

      ! Set-up initial guess

      u(1:nxy) = zero

      ! Solve the equation

      CALL nag_pde_ell_mg_sol(nx,ny,a,rhs,u,acc=acc,max_iter=max_iter, &
       num_iter=num_iter)

      ! Print out the Solution

      WRITE (nag_std_out,*) ' '
      WRITE (nag_std_out,*) ' Exact solution above computed solution'
      WRITE (nag_std_out,*) ' '
      WRITE (nag_std_out,999) '  I/J', (i,i=1,nx)
      rms_err = zero
      DO j = ny, 1, -1
        WRITE (nag_std_out,*) ' '
        WRITE (nag_std_out,998) j, (SIN(x(i+(j-1)*nx))*SIN(y(i+(j-1)*nx)),i=1, &
         nx)
        WRITE (nag_std_out,998) j, (u(i+(j-1)*nx),i=1,nx)
        DO i = 1, nx
          rms_err = rms_err + (SIN(x(i+(j-1)*nx))*SIN(y(i+ &
           (j-1)*nx))-u(i+(j-1)*nx))**2
        END DO
      END DO
      rms_err = SQRT(rms_err/REAL(nx*ny,kind=wp))
      WRITE (nag_std_out,*) ' '
      WRITE (nag_std_out,997) ' Number of Iteration ', num_iter
      WRITE (nag_std_out,996) ' RMS Error ', rms_err

999   FORMAT (1X,A,10I7:/(6X,10I7))
998   FORMAT (1X,I3,2X,10F7.3:/(6X,10F7.3))
997   FORMAT (1X,A,I3)
996   FORMAT (1X,A,1P,E10.2)

    END PROGRAM nag_pde_ell_mg_ex02
