PROGRAM d03fafe

!      D03FAF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : d03faf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: dx, dy, dz, error, lambda, pertrb,   &
                                          t, xf, xs, yf, ys, yy, zf, zs, zz
       INTEGER                         :: i, ifail, j, k, l, lbdcnd, ldf,      &
                                          ldf2, lwrk, m, maxlm, mbdcnd, n,     &
                                          nbdcnd
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: bdxf(:,:), bdxs(:,:), bdyf(:,:),     &
                                          bdys(:,:), bdzf(:,:), bdzs(:,:),     &
                                          f(:,:,:), w(:), x(:), y(:), z(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          abs, cos, real, sin
!      .. Executable Statements ..
       WRITE (nout,*) 'D03FAF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) l, m, n, maxlm
       ldf = l + 1
       ldf2 = m + 1
       lwrk = 2*(n+1)*maxlm + 3*l + 3*m + 4*n + 6000
       ALLOCATE (bdxf(ldf2,n+1),bdxs(ldf2,n+1),bdyf(ldf,n+1),bdys(ldf,n+1), &
          bdzf(ldf,m+1),bdzs(ldf,m+1),f(ldf,ldf2,n+1),w(lwrk),x(l+1),y(m+1), &
          z(n+1))
       READ (nin,*) lambda
       READ (nin,*) xs, xf
       READ (nin,*) ys, yf
       READ (nin,*) zs, zf
       READ (nin,*) lbdcnd, mbdcnd, nbdcnd

!      Define the grid points for later use.

       dx = (xf-xs)/real(l,kind=nag_wp)
       DO i = 1, l + 1
          x(i) = xs + real(i-1,kind=nag_wp)*dx
       END DO
       dy = (yf-ys)/real(m,kind=nag_wp)
       DO j = 1, m + 1
          y(j) = ys + real(j-1,kind=nag_wp)*dy
       END DO
       dz = (zf-zs)/real(n,kind=nag_wp)
       DO k = 1, n + 1
          z(k) = zs + real(k-1,kind=nag_wp)*dz
       END DO

!      Define the array of derivative boundary values.

       DO j = 1, m + 1
          yy = sin(y(j))
          bdzf(1:l+1,j) = -yy*x(1:l+1)**4
       END DO

!      Note that for this example all other boundary arrays are
!      dummy variables.

!      We define the function boundary values in the F array.

       f(1,1:m+1,1:n+1) = 0.0_nag_wp
       DO k = 1, n + 1
          zz = cos(z(k))
          DO j = 1, m + 1
             f(l+1,j,k) = zz*sin(y(j))
          END DO
       END DO
       f(1:l+1,1:m+1,1) = -bdzf(1:l+1,1:m+1)

!      Define the values of the right hand side of the Helmholtz
!      equation.

       DO k = 2, n + 1
          zz = 4.0_nag_wp*cos(z(k))
          DO j = 1, m + 1
             yy = sin(y(j))*zz
             DO i = 2, l
                f(i,j,k) = yy*x(i)**2*(3.0_nag_wp-x(i)**2)
             END DO
          END DO
       END DO

!      Call D03FAF to generate and solve the finite difference equation.
!      ifail: behaviour on error exit   
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft     
       ifail = 0
       CALL d03faf(xs,xf,l,lbdcnd,bdxs,bdxf,ys,yf,m,mbdcnd,bdys,bdyf,zs,zf,n, &
          nbdcnd,bdzs,bdzf,lambda,ldf,ldf2,f,pertrb,w,lwrk,ifail)

!      Compute discretization error.  The exact solution to the
!      problem is

!      U(X,Y,Z) = X**4*SIN(Y)*COS(Z)

       error = 0.0_nag_wp
       DO k = 1, n + 1
          zz = cos(z(k))
          DO j = 1, m + 1
             yy = sin(y(j))*zz
             DO i = 1, l + 1
                t = abs(f(i,j,k)-yy*x(i)**4)
                IF (t>error) error = t
             END DO
          END DO
       END DO
       WRITE (nout,*)
       WRITE (nout,99999) error

99999  FORMAT (1X,'Maximum component of discretization error =',1P,E13.6)
    END PROGRAM d03fafe