!   D03RAF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d03rafe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter        :: activ_energy = 20.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: alpha = 50.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: beta = 300.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: diffusion = 0.1_nag_wp
      Real (Kind=nag_wp), Parameter        :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: reaction_rate = 5.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: xmax = one
      Real (Kind=nag_wp), Parameter        :: ymax = one
      Real (Kind=nag_wp), Parameter        :: zero = 0.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: heat_release = one
      Real (Kind=nag_wp), Parameter        :: xmin = zero
      Real (Kind=nag_wp), Parameter        :: ymin = zero
      Integer, Parameter                   :: itrace = 0, nin = 5, nout = 6,   &
                                              npde1 = 1, npde2 = 2
!     .. Local Scalars ..
      Real (Kind=nag_wp), Save             :: damkohler
      Integer, Save                        :: iout
    Contains
      Subroutine pdeiv1(npts,npde,t,x,y,u)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: u(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: x(npts), y(npts)
!       .. Executable Statements ..
        u(1:npts,1:npde) = one
        Return
      End Subroutine pdeiv1
      Subroutine pdedef1(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: res(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npts,npde), ut(npts,npde),   &
                                                ux(npts,npde), uxx(npts,npde), &
                                                uxy(npts,npde), uy(npts,npde), &
                                                uyy(npts,npde), x(npts), y(npts)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        res(1:npts,1:npde) = ut(1:npts,1:npde) - diffusion*(uxx(1:npts,1:npde) &
          +uyy(1:npts,1:npde)) - damkohler*(one+heat_release-u(1:npts,1:npde)) &
          *exp(-activ_energy/u(1:npts,1:npde))
        Return
      End Subroutine pdedef1
      Subroutine bndry1(npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res)

!       .. Use Statements ..
        Use nag_library, Only: x02ajf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: nbpts, npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: res(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npts,npde), ut(npts,npde),   &
                                                ux(npts,npde), uy(npts,npde),  &
                                                x(npts), y(npts)
        Integer, Intent (In)                 :: lbnd(nbpts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: tol
        Integer                              :: i, j
!       .. Intrinsic Procedures ..
        Intrinsic                            :: abs
!       .. Executable Statements ..
        tol = 10._nag_wp*x02ajf()

        Do i = 1, nbpts
          j = lbnd(i)

          If (abs(x(j))<=tol) Then
            res(j,1:npde) = ux(j,1:npde)
          Else If (abs(x(j)-one)<=tol) Then
            res(j,1:npde) = u(j,1:npde) - one
          Else If (abs(y(j))<=tol) Then
            res(j,1:npde) = uy(j,1:npde)
          Else If (abs(y(j)-one)<=tol) Then
            res(j,1:npde) = u(j,1:npde) - one
          End If
        End Do

        Return
      End Subroutine bndry1
      Subroutine monit1(npde,t,dt,dtnew,tlast,nlev,ngpts,xpts,ypts,lsol,sol, &
        ierr)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: dt, dtnew, t
        Integer, Intent (Inout)              :: ierr
        Integer, Intent (In)                 :: nlev, npde
        Logical, Intent (In)                 :: tlast
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: sol(*), xpts(*), ypts(*)
        Integer, Intent (In)                 :: lsol(nlev), ngpts(nlev)
!       .. Local Scalars ..
        Integer                              :: i, ipsol, level, npts
!       .. Executable Statements ..
        If (tlast) Then

!         Print solution

          If (iout==2) Then
            level = nlev
            Write (nout,99999) level, t
            Write (nout,99998)
            npts = ngpts(level)
            ipsol = lsol(level)

            Do i = 1, npts, 4
              Write (nout,99997) xpts(i), ypts(i), sol(ipsol+i)
            End Do

            Write (nout,*)
          End If

        End If

        Return
99999   Format (1X,'Solution at every 4th grid point in level',I10, &
          ' at time ',F8.4,':')
99998   Format (1X/7X,'x',10X,'y',8X,'approx u'/)
99997   Format (3(1X,E11.4))
      End Subroutine monit1
      Subroutine pdeiv2(npts,npde,t,x,y,u)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: u(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: x(npts), y(npts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: fourpi
!       .. Intrinsic Procedures ..
        Intrinsic                            :: sin
!       .. Executable Statements ..
        fourpi = 4.0_nag_wp*x01aaf(fourpi)
        u(1:npts,1) = 10.0_nag_wp + (16.0_nag_wp*x(1:npts)*(one-x(1:npts))*y(1 &
          :npts)*(one-y(1:npts)))**2
        u(1:npts,2) = -one - alpha*x(1:npts)*y(1:npts) - &
          beta*sin(fourpi*x(1:npts))*sin(fourpi*y(1:npts)) + &
          1.0E4_nag_wp*u(1:npts,1)
        Return
      End Subroutine pdeiv2
      Subroutine pdedef2(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: res(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npts,npde), ut(npts,npde),   &
                                                ux(npts,npde), uxx(npts,npde), &
                                                uxy(npts,npde), uy(npts,npde), &
                                                uyy(npts,npde), x(npts), y(npts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: fourpi
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable      :: b1(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: sin
!       .. Executable Statements ..
        fourpi = 4.0_nag_wp*x01aaf(fourpi)
        Allocate (b1(npts))
        b1(1:npts) = one + alpha*x(1:npts)*y(1:npts) + &
          beta*sin(fourpi*x(1:npts))*sin(fourpi*y(1:npts))
        res(1:npts,1) = ut(1:npts,1) - (uxx(1:npts,1)+uyy(1:npts,1)) - &
          u(1:npts,1)*(b1(1:npts)-u(1:npts,1)-0.5E-6_nag_wp*u(1:npts,2))
        res(1:npts,2) = -0.05_nag_wp*(uxx(1:npts,2)+uyy(1:npts,2)) - &
          u(1:npts,2)*(-b1(1:npts)+1.0E4_nag_wp*u(1:npts,1)-u(1:npts,2))
        Deallocate (b1)
        Return
      End Subroutine pdedef2
      Subroutine bndry2(npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res)

!       .. Use Statements ..
        Use nag_library, Only: x02ajf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: nbpts, npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: res(npts,npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npts,npde), ut(npts,npde),   &
                                                ux(npts,npde), uy(npts,npde),  &
                                                x(npts), y(npts)
        Integer, Intent (In)                 :: lbnd(nbpts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: tol
        Integer                              :: i, j
!       .. Intrinsic Procedures ..
        Intrinsic                            :: abs
!       .. Executable Statements ..
        tol = 10.0_nag_wp*x02ajf()

        Do i = 1, nbpts
          j = lbnd(i)

          If (abs(x(j))<=tol .Or. abs(x(j)-one)<=tol) Then
            res(j,1:npde) = ux(j,1:npde)
          Else If (abs(y(j))<=tol .Or. abs(y(j)-one)<=tol) Then
            res(j,1:npde) = uy(j,1:npde)
          End If
        End Do

        Return
      End Subroutine bndry2
      Subroutine monit2(npde,t,dt,dtnew,tlast,nlev,ngpts,xpts,ypts,lsol,sol, &
        ierr)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: dt, dtnew, t
        Integer, Intent (Inout)              :: ierr
        Integer, Intent (In)                 :: nlev, npde
        Logical, Intent (In)                 :: tlast
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: sol(*), xpts(*), ypts(*)
        Integer, Intent (In)                 :: lsol(nlev), ngpts(nlev)
!       .. Local Scalars ..
        Integer                              :: i, ipsol, level, npts
!       .. Executable Statements ..
        If (tlast) Then

!         Print solution

          If (iout==2) Then
            level = nlev
            Write (nout,99999) level, t
            Write (nout,99998)
            npts = ngpts(level)
            ipsol = lsol(level)

            Do i = 1, npts, 2
              Write (nout,99997) xpts(i), ypts(i), sol(ipsol+i), &
                sol(ipsol+npts+i)
            End Do

            Write (nout,*)
          End If

        End If

        Return
99999   Format (1X,'Solution at every 2nd grid point in level',I10, &
          ' at time ',F8.4,':')
99998   Format (1X/7X,'x',10X,'y',9X,'approx c1',3X,'approx c2'/)
99997   Format (2(1X,E11.4),2X,E11.4,2X,E11.4)
      End Subroutine monit2
      Subroutine compute_wkspace_lens(maxlev,npde,maxpts,lenrwk,leniwk,lenlwk)

!       Returns suitable workspace lengths for the two problems
!       being solved, based on trial-and-error.

!       .. Scalar Arguments ..
        Integer, Intent (Out)                :: leniwk, lenlwk, lenrwk
        Integer, Intent (In)                 :: maxlev, maxpts, npde
!       .. Executable Statements ..
        lenrwk = 2*maxpts*npde*(5*maxlev+18*npde+9) + 2*maxpts
        leniwk = 2*maxpts*(14+5*maxlev) + 7*maxlev + 2
        lenlwk = 2*maxpts + 400
        Return
      End Subroutine compute_wkspace_lens
      Subroutine print_statistics(ts,iwk,maxlev)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: ts
        Integer, Intent (In)                 :: maxlev
!       .. Array Arguments ..
        Integer, Intent (In)                 :: iwk(6*maxlev+2)
!       .. Local Scalars ..
        Integer                              :: i, j
!       .. Executable Statements ..
        Write (nout,'(1X,A)') 'Statistics:'
        Write (nout,99999) 'Time = ', ts
        Write (nout,99998) 'Total number of accepted timesteps =', iwk(1)
        Write (nout,99998) 'Total number of rejected timesteps =', iwk(2)
        Write (nout,'(1X,4(/,A))') &
          '             T o t a l  n u m b e r  o f    ', &
          '          Residual  Jacobian    Newton   Lin sys', &
          '             evals     evals     iters     iters', ' At level '

        Do j = 1, maxlev
          If (iwk(j+2)/=0) Then
            Write (nout,99997) j, (iwk(j+2+i*maxlev),i=0,3)
          End If
        End Do

        Write (nout,'(1X,3(/,A))') &
          '             M a x i m u m  n u m b e r  o f', &
          '              Newton iters    Lin sys iters ', ' At level '

        Do j = 1, maxlev
          If (iwk(j+2)/=0) Then
            Write (nout,99996) j, (iwk(j+2+i*maxlev),i=4,5)
          End If
        End Do
        Write (nout,*)

        Return
99999   Format (1X,A,F8.4)
99998   Format (1X,A,I5)
99997   Format (I8,4I10)
99996   Format (I8,2I14)
      End Subroutine print_statistics
    End Module d03rafe_mod
    Program d03rafe

!     D03RAF Example Main Program

!     .. Use Statements ..
      Use d03rafe_mod, Only: nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Executable Statements ..
      Write (nout,*) 'D03RAF Example Program Results'

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use nag_library, Only: d03raf, nag_wp
        Use d03rafe_mod, Only: activ_energy, bndry1, compute_wkspace_lens,     &
                               damkohler, heat_release, iout, itrace, monit1,  &
                               nin, npde1, one, pdedef1, pdeiv1,               &
                               print_statistics, reaction_rate, xmax, xmin,    &
                               ymax, ymin, zero
!       .. Parameters ..
        Integer, Parameter                   :: opti(4) = 0
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: tols, tolt, tout, ts
        Integer                              :: ifail, ind, leniwk, lenlwk,    &
                                                lenrwk, maxlev, npde, npts,    &
                                                nx, ny
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: dt(3), twant(2)
        Real (Kind=nag_wp), Allocatable      :: optr(:,:), rwk(:)
        Integer, Allocatable                 :: iwk(:)
        Logical, Allocatable                 :: lwk(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp, max
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 1'
        Write (nout,*)
!       Skip heading in data file
        Read (nin,*)
        Read (nin,*) npts

        npde = npde1
        damkohler = reaction_rate*exp(activ_energy)/ &
          (heat_release*activ_energy)

        dt(1:3) = (/0.1E-2_nag_wp,zero,zero/)
        twant(1:2) = (/0.24_nag_wp,0.25_nag_wp/)
        ts = zero
        ind = 0
        nx = 21
        ny = 21
        tols = 0.5_nag_wp
        tolt = 0.01_nag_wp
        maxlev = max(opti(1),3)

        Call compute_wkspace_lens(maxlev,npde,npts,lenrwk,leniwk,lenlwk)

        Allocate (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde))
        optr(1:3,1:npde) = one

        Do iout = 1, 2
          tout = twant(iout)

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
          ifail = 0
          Call d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt, &
            pdedef1,bndry1,pdeiv1,monit1,opti,optr,rwk,lenrwk,iwk,leniwk,lwk, &
            lenlwk,itrace,ind,ifail)

          Call print_statistics(ts,iwk,maxlev)

        End Do

        Return
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use nag_library, Only: d03raf, nag_wp
        Use d03rafe_mod, Only: bndry2, compute_wkspace_lens, iout, itrace,     &
                               monit2, nin, npde2, one, pdedef2, pdeiv2,       &
                               print_statistics, xmax, xmin, ymax, ymin, zero
!       .. Parameters ..
        Integer, Parameter                   :: opti(4) = (/4,0,0,0/)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: tols, tolt, tout, ts
        Integer                              :: ifail, ind, leniwk, lenlwk,    &
                                                lenrwk, maxlev, npde, npts,    &
                                                nx, ny
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: dt(3), twant(2)
        Real (Kind=nag_wp), Allocatable      :: optr(:,:), rwk(:)
        Integer, Allocatable                 :: iwk(:)
        Logical, Allocatable                 :: lwk(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: max
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 2'
        Write (nout,*)
        Read (nin,*)
        Read (nin,*) npts

        npde = npde2
        dt(1:3) = (/0.5E-3_nag_wp,1.0E-6_nag_wp,zero/)
        twant(1:2) = (/0.01_nag_wp,0.025_nag_wp/)
        ts = zero
        ind = 0
        nx = 11
        ny = 11
        tols = 0.075_nag_wp
        tolt = 0.1_nag_wp
        maxlev = max(opti(1),3)

        Call compute_wkspace_lens(maxlev,npde,npts,lenrwk,leniwk,lenlwk)

        Allocate (rwk(lenrwk),iwk(leniwk),lwk(lenlwk),optr(3,npde))
        optr(1,1:npde) = (/250.0_nag_wp,1.5E6_nag_wp/)
        optr(2:3,1:npde) = one

        Do iout = 1, 2
          tout = twant(iout)

          ifail = 0
          Call d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt, &
            pdedef2,bndry2,pdeiv2,monit2,opti,optr,rwk,lenrwk,iwk,leniwk,lwk, &
            lenlwk,itrace,ind,ifail)

          Call print_statistics(ts,iwk,maxlev)

        End Do

        Return
      End Subroutine ex2
    End Program d03rafe
