Program d03eafe

!     D03EAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: d03eaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: alpha, alpha0, alpsav, p, q
      Integer                          :: i, ifail, j, ldc, m, n, n1, n1p1,    &
                                          np1, np4, npts
      Logical                          :: dorm, ext, stage1
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: c(:,:), phi(:), phid(:), x(:), y(:)
      Integer, Allocatable             :: icint(:)
!     .. Executable Statements ..
      Write (nout,*) 'D03EAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) m, n
      n1 = 2*(n+m) - 2
      n1p1 = n1 + 1
      np1 = n + 1
      np4 = n + 4
      ldc = n + 1
      Allocate (c(ldc,np4),phi(n),phid(n),x(n1p1),y(n1p1),icint(np1))
      Read (nin,*) ext, dorm
      If (.Not. dorm) Then
        Read (nin,*) alpha0
        alpha = alpha0
      End If
      Read (nin,*)(x(i),y(i),i=1,n1+1)
      Read (nin,*)(phi(i),phid(i),i=1,n)

      stage1 = .True.

!     ifail: behaviour on error exit
!            =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call d03eaf(stage1,ext,dorm,n,p,q,x,y,n1p1,phi,phid,alpha,c,ldc,np4,     &
        icint,np1,ifail)

      alpsav = alpha
      Write (nout,*)
      If (.Not. dorm) Then
        If (ext) Then
          Write (nout,*) 'Exterior Neumann problem'
          Write (nout,*)
          Write (nout,99998) 'C=', alpha0
        Else
          Write (nout,*) 'Interior Neumann problem'
          Write (nout,*)
          Write (nout,99999) 'Total integral =', alpha0
        End If
      Else If (ext) Then
        Write (nout,*)
        Write (nout,*) 'Exterior problem'
        Write (nout,*)
        Write (nout,99999) 'Computed C =', alpsav
      End If
      j = 2
      Write (nout,*)
      Write (nout,*) 'Nodes'
      Write (nout,99996) 'X', 'Y', 'PHI', 'PHID'
      Do i = 1, n
        If (x(j)==x(j-1) .And. y(j)==y(j-1)) Then
          j = j + 2
        End If
        Write (nout,99997) x(j), y(j), phi(i), phid(i)
        j = j + 2
      End Do
      stage1 = .False.
      Write (nout,*)
      Write (nout,*) 'Selected points'
      Write (nout,*) '           X              Y              PHI'
      Read (nin,*) npts
      Do i = 1, npts
        Read (nin,*) p, q
        alpha = alpsav

        ifail = 0
        Call d03eaf(stage1,ext,dorm,n,p,q,x,y,n1p1,phi,phid,alpha,c,ldc,np4,   &
          icint,np1,ifail)

        Write (nout,99997) p, q, alpha
      End Do

99999 Format (1X,A,F15.2)
99998 Format (1X,A,E15.4)
99997 Format (1X,4F15.2)
99996 Format (1X,A12,A15,A17,A16)
    End Program d03eafe