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

    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
!dec$ attributes offload : mic :: pdeiv1
       SUBROUTINE pdeiv1(npts,npde,t,x,y,u)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: pdeiv1
       END SUBROUTINE pdeiv1
!dec$ attributes offload : mic :: pdedef1
       SUBROUTINE pdedef1(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: pdedef1
       END SUBROUTINE pdedef1
!dec$ attributes offload : mic :: bndry1
       SUBROUTINE bndry1(npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res)

!         .. Use Statements ..
          USE nag_library, ONLY : x02ajf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: bndry1
       END SUBROUTINE bndry1
!dec$ attributes offload : mic :: monit1
       SUBROUTINE monit1(npde,t,dt,dtnew,tlast,nlev,ngpts,xpts,ypts,lsol,sol, &
             ierr)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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))
!dec$ attributes offload : mic :: monit1
       END SUBROUTINE monit1
!dec$ attributes offload : mic :: pdeiv2
       SUBROUTINE pdeiv2(npts,npde,t,x,y,u)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: pdeiv2
       END SUBROUTINE pdeiv2
!dec$ attributes offload : mic :: pdedef2
       SUBROUTINE pdedef2(npts,npde,t,x,y,u,ut,ux,uy,uxx,uxy,uyy,res)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: pdedef2
       END SUBROUTINE pdedef2
!dec$ attributes offload : mic :: bndry2
       SUBROUTINE bndry2(npts,npde,t,x,y,u,ut,ux,uy,nbpts,lbnd,res)

!         .. Use Statements ..
          USE nag_library, ONLY : x02ajf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: bndry2
       END SUBROUTINE bndry2
!dec$ attributes offload : mic :: monit2
       SUBROUTINE monit2(npde,t,dt,dtnew,tlast,nlev,ngpts,xpts,ypts,lsol,sol, &
             ierr)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)
!dec$ attributes offload : mic :: monit2
       END SUBROUTINE monit2
!dec$ attributes offload : mic :: compute_wkspace_lens
       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.

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: compute_wkspace_lens
       END SUBROUTINE compute_wkspace_lens
!dec$ attributes offload : mic :: print_statistics
       SUBROUTINE print_statistics(ts,iwk,maxlev)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)
!dec$ attributes offload : mic :: print_statistics
       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
!dec$ attributes offload : mic :: ex1
       SUBROUTINE ex1

!         .. Use Statements ..
          USE nag_library, ONLY : d03raf
          USE d03rafe_mod, ONLY : activ_energy, bndry1, compute_wkspace_lens,  &
                                  damkohler, heat_release, iout, itrace,       &
                                  monit1, nag_wp, nin, npde1, one, pdedef1,    &
                                  pdeiv1, print_statistics, reaction_rate,     &
                                  xmax, xmin, ymax, ymin, zero
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: ex1
       END SUBROUTINE ex1
!dec$ attributes offload : mic :: ex2
       SUBROUTINE ex2

!         .. Use Statements ..
          USE nag_library, ONLY : d03raf
          USE d03rafe_mod, ONLY : bndry2, compute_wkspace_lens, iout, itrace,  &
                                  monit2, nag_wp, nin, npde2, one, pdedef2,    &
                                  pdeiv2, print_statistics, xmax, xmin, ymax,  &
                                  ymin, zero
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: ex2
       END SUBROUTINE ex2
    END PROGRAM d03rafe
