!   D03RAF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: bndry1, bndry2,                      &
                                          compute_wkspace_lens, monit1,        &
                                          monit2, monit_dummy, pdedef1,        &
                                          pdedef2, pdeiv1, pdeiv2,             &
                                          print_statistics
!     .. Parameters ..
      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, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: xmax = one
      Real (Kind=nag_wp), Parameter, Public :: ymax = one
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: xmin = zero
      Real (Kind=nag_wp), Parameter, Public :: ymin = zero
      Integer, Parameter, Public       :: itrace = 0, nin = 5, nout = 6,       &
                                          npde1 = 1, npde2 = 2
    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)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: activ_energy = 20.0_nag_wp
        Real (Kind=nag_wp), Parameter  :: diffusion = 0.1_nag_wp
        Real (Kind=nag_wp), Parameter  :: heat_release = 1.0_nag_wp
        Real (Kind=nag_wp), Parameter  :: reaction_rate = 5.0_nag_wp
!       .. 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)             :: damkohler
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        damkohler = reaction_rate*exp(activ_energy)/                           &
          (heat_release*activ_energy)
        !$Omp Do
        Do i = 1, npts
          res(i,1:npde) = ut(i,1:npde) - diffusion*(uxx(i,1:npde)+uyy(i,1:npde &
            )) - damkohler*(1.0E0_nag_wp+heat_release-u(i,1:npde))*exp(        &
            -activ_energy/u(i,1:npde))
        End Do
        !$Omp End Do
        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, k, level, npts
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sum
!       .. Executable Statements ..
        If (tlast) Then

!         Print solution

          level = nlev
          Write (nout,99999) level, t
          Write (nout,99998)
          npts = ngpts(level)
          ipsol = lsol(level)
          k = sum(ngpts(1:nlev-1))
          Do i = 1, npts, 4
            Write (nout,99997) xpts(k+i), ypts(k+i), sol(ipsol+i)
          End Do

          Write (nout,*)

        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 (1X,1P,E11.4,2(1X,1P,E11.3))
      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)             :: b1, fourpi
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sin
!       .. Executable Statements ..
        fourpi = 4.0_nag_wp*x01aaf(fourpi)
        Do i = 1, npts
          b1 = 1.0E0_nag_wp + alpha*x(i)*y(i) + beta*sin(fourpi*x(i))*sin(     &
            fourpi*y(i))
          res(i,1) = ut(i,1) - (uxx(i,1)+uyy(i,1)) -                           &
            u(i,1)*(b1-u(i,1)-0.5E-6_nag_wp*u(i,2))
          res(i,2) = -0.05E0_nag_wp*(uxx(i,2)+uyy(i,2)) -                      &
            u(i,2)*(-b1+1.0E4_nag_wp*u(i,1)-u(i,2))
        End Do
        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, k, level, npts
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sum
!       .. Executable Statements ..
        If (tlast) Then

!         Print solution

          level = nlev
          Write (nout,99999) level, t
          Write (nout,99998)
          npts = ngpts(level)
          ipsol = lsol(level)
          k = sum(ngpts(1:nlev-1))

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

          Write (nout,*)

        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 (1P,2(1X,E11.3),2X,E11.3,2X,E11.3)
      End Subroutine monit2
      Subroutine monit_dummy(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)
!       .. Executable Statements ..
        Return
      End Subroutine monit_dummy
      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
!       .. Local Arrays ..
        Integer                        :: istats(4)
!       .. 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))')                                             &
          '             Total  number (rounded) of     ',                      &
          '          Residual  Jacobian    Newton   Lin sys',                  &
          '             evals     evals     iters     iters', ' At level '

        Do j = 1, maxlev
          If (iwk(j+2)/=0) Then
            istats(1:4) = iwk(j+2:j+2+3*maxlev:maxlev)
            Call round_statistics(istats)
            Write (nout,99997) j, istats(1:4)
          End If
        End Do

        Write (nout,'(1X,3(/,A))') '             Maximum  number of   ',       &
          '              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
      Subroutine round_statistics(istat)

!       .. Array Arguments ..
        Integer, Intent (Inout)        :: istat(4)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: lt
        Integer                        :: i, k
!       .. Intrinsic Procedures ..
        Intrinsic                      :: int, log, real
!       .. Executable Statements ..
        lt = log(10.0_nag_wp)
        Do i = 1, 4
          k = int(log(real(istat(i),kind=nag_wp))/lt)
          k = 10**k
          istat(i) = k*((istat(i)+k/2)/k)
        End Do
      End Subroutine round_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 d03rafe_mod, Only: bndry1, compute_wkspace_lens, itrace, monit1,   &
                               monit_dummy, nin, npde1, one, pdedef1, pdeiv1,  &
                               print_statistics, xmax, xmin, ymax, ymin, zero
        Use nag_library, Only: d03raf, nag_wp
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: tols, tolt, tout, ts
        Integer                        :: i, 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(:)
        Integer                        :: opti(4)
        Logical, Allocatable           :: lwk(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: 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

        dt(1:3) = (/0.1E-2_nag_wp,zero,zero/)
        twant(1:2) = (/0.24_nag_wp,0.25_nag_wp/)
        ts = zero

!       Specify that we are starting the integration in time (ind = 0
!       normally).
!       Note: we have parallelized the loop in the function pdedef1 using
!       OpenMP so set alternative value of ind to indicate that this can be
!       run in parallel if we are using a multithreaded implementation.
!       Either option is OK for serial NAG Library implementations from
!       Mark 25 onwards.
        ind = 10

        nx = 41
        ny = 41
        tols = 0.5_nag_wp
        tolt = 0.01_nag_wp
        opti(1) = 6
        opti(2:4) = 0
        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 i = 1, 2
          tout = twant(i)

!         ifail: behaviour on error exit
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
          If (i==1) Then
!           Use monit_dummy to avoid output first time around
            Call d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt,   &
              pdedef1,bndry1,pdeiv1,monit_dummy,opti,optr,rwk,lenrwk,iwk,      &
              leniwk,lwk,lenlwk,itrace,ind,ifail)
          Else
            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)
          End If

          Call print_statistics(ts,iwk,maxlev)

        End Do

        Return
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use d03rafe_mod, Only: bndry2, compute_wkspace_lens, itrace, monit2,   &
                               monit_dummy, nin, npde2, one, pdedef2, pdeiv2,  &
                               print_statistics, xmax, xmin, ymax, ymin, zero
        Use nag_library, Only: d03raf, nag_wp
!       .. Parameters ..
        Integer, Parameter             :: opti(4) = (/4,0,0,0/)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: tols, tolt, tout, ts
        Integer                        :: i, 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

!       Specify that we are starting the integration in time (ind = 0
!       normally).
!       Note: In this second example we have not added OpenMP directives to
!       parallelize the loop in the function pdedef2. Thus the alternative
!       ind=10 must not be specified here, as this will not function correctly
!       if a multithreaded implementation of the NAG Library is used. Adding
!       OpenMP to pdedef2, that would enable ind=10 to be used here safely, is
!       left as an exercise for the reader.
        ind = 0

        nx = 11
        ny = 11
        tols = 0.1_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 i = 1, 2
          tout = twant(i)

          ifail = 0
          If (i==1) Then
!           Use monit_dummy to avoid output first time around
            Call d03raf(npde,ts,tout,dt,xmin,xmax,ymin,ymax,nx,ny,tols,tolt,   &
              pdedef2,bndry2,pdeiv2,monit_dummy,opti,optr,rwk,lenrwk,iwk,      &
              leniwk,lwk,lenlwk,itrace,ind,ifail)
          Else
            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)
          End If

          Call print_statistics(ts,iwk,maxlev)

        End Do

        Return
      End Subroutine ex2
    End Program d03rafe