!   D03PCF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d03pcfe_mod

!      D03PCF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: nin = 5, nout = 6, npde = 2
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: alpha
    CONTAINS
       SUBROUTINE pdedef(npde,t,x,u,ux,p,q,r,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t, x
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: npde
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: p(npde,npde), q(npde), r(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde)
!         .. Executable Statements ..
          q(1) = 4.0_nag_wp*alpha*(u(2)+x*ux(2))
          q(2) = 0.0_nag_wp
          r(1) = x*ux(1)
          r(2) = ux(2) - u(1)*u(2)
          p(1,1) = 0.0_nag_wp
          p(1,2) = 0.0_nag_wp
          p(2,1) = 0.0_nag_wp
          p(2,2) = 1.0_nag_wp - x*x
          RETURN
       END SUBROUTINE pdedef
       SUBROUTINE bndary(npde,t,u,ux,ibnd,beta,gamma,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, npde
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: beta(npde), gamma(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde)
!         .. Executable Statements ..
          IF (ibnd==0) THEN
             beta(1) = 0.0_nag_wp
             beta(2) = 1.0_nag_wp
             gamma(1) = u(1)
             gamma(2) = -u(1)*u(2)
          ELSE
             beta(1) = 1.0_nag_wp
             beta(2) = 0.0_nag_wp
             gamma(1) = -u(1)
             gamma(2) = u(2)
          END IF
          RETURN
       END SUBROUTINE bndary
       SUBROUTINE uinit(u,x,npts)

!         Routine for PDE initial conditon

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: npts
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(2,npts)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
!         .. Local Scalars ..
          INTEGER                             :: i
!         .. Executable Statements ..
          DO i = 1, npts
             u(1,i) = 2.0_nag_wp*alpha*x(i)
             u(2,i) = 1.0_nag_wp
          END DO
          RETURN
       END SUBROUTINE uinit
    END MODULE d03pcfe_mod
    PROGRAM d03pcfe

!      D03PCF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03pcf, d03pzf, nag_wp, x01aaf
       USE d03pcfe_mod, ONLY : alpha, bndary, nin, nout, npde, pdedef, uinit
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: acc, hx, pi, piby2, tout, ts
       INTEGER                             :: i, ifail, ind, intpts, it,       &
                                              itask, itrace, itype, lisave,    &
                                              lrsave, m, neqn, npts, nwk
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: rsave(:), u(:,:), uout(:,:,:),   &
                                              x(:), xout(:)
       INTEGER, ALLOCATABLE                :: isave(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real, sin
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PCF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) intpts, npts, itype
       neqn = npde*npts
       lisave = neqn + 24
       nwk = (10+6*npde)*neqn
       lrsave = nwk + (21+3*npde)*npde + 7*npts + 54
       ALLOCATE (rsave(lrsave),u(npde,npts),uout(npde,intpts,itype),x(npts), &
          xout(intpts),isave(lisave))

       READ (nin,*) xout(1:intpts)
       READ (nin,*) acc, alpha
       READ (nin,*) m, itrace
       ind = 0
       itask = 1

!      Set spatial mesh points

       piby2 = 0.5_nag_wp*x01aaf(pi)
       hx = piby2/real(npts-1,kind=nag_wp)
       x(1) = 0.0_nag_wp
       x(npts) = 1.0_nag_wp
       DO i = 2, npts - 1
          x(i) = sin(hx*real(i-1,kind=nag_wp))
       END DO

!      Set initial conditions
       READ (nin,*) ts, tout

!      Set the initial values
       CALL uinit(u,x,npts)

       DO it = 1, 5
          tout = 10.0_nag_wp*tout

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
          ifail = 0
          CALL d03pcf(npde,m,ts,tout,pdedef,bndary,u,npts,x,acc,rsave,lrsave, &
             isave,lisave,itask,itrace,ind,ifail)

          IF (it==1) THEN
             WRITE (nout,99999) acc, alpha
             WRITE (nout,99998) xout(1:6)
          END IF

!         Interpolate at required spatial points

          ifail = 0
          CALL d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)

          WRITE (nout,99996) tout, uout(1,1:intpts,1)
          WRITE (nout,99995) uout(2,1:intpts,1)
       END DO

!      Print integration statistics

       WRITE (nout,99997) isave(1), isave(2), isave(3), isave(5)


99999  FORMAT (//' Accuracy requirement  = ',E12.5/' Parameter ALPHA =', &
          '       ',E12.3/)
99998  FORMAT ('   T  /  X   ',6F8.4/)
99997  FORMAT (' Number of integration steps in time                  ', &
          I4/' Number of residual evaluations of resulting ODE system', &
          I4/' Number of Jacobian evaluations                        ', &
          I4/' Number of iterations of nonlinear solver              ',I4)
99996  FORMAT (1X,F7.4,' U(1)',6F8.4)
99995  FORMAT (9X,'U(2)',6F8.4/)
    END PROGRAM d03pcfe