    MODULE nag_pde_parab_1d_ex02_mod

      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Default Accessibility ..
      PUBLIC
      ! .. Intrinsic Functions ..
      INTRINSIC KIND
      ! .. Parameters ..
      INTEGER, PARAMETER :: wp = KIND(1.0D0)

    CONTAINS

      SUBROUTINE pde_coef(t,x,u,du_dx,p,q,r,finish,v,vdot,i_comm,r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC SIZE
        ! .. Scalar Arguments ..
        INTEGER, INTENT (OUT) :: finish
        REAL (wp), INTENT (IN) :: t
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), INTENT (IN) :: du_dx(:,:), u(:,:), x(:)
        REAL (wp), INTENT (OUT) :: p(:,:,:), q(:,:), r(:,:)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:), v(:), vdot(:)
        ! .. Local Scalars ..
        INTEGER :: nptl
        ! .. Executable Statements ..
        nptl = SIZE(x)

        finish = 1

        q(1,1:nptl) = u(2,1:nptl)
        q(2,1:nptl) = u(1,1:nptl)*du_dx(2,1:nptl) - du_dx(1,1:nptl)*u(2,1:nptl &
         )
        r(1,1:nptl) = du_dx(1,1:nptl)
        r(2,1:nptl) = du_dx(2,1:nptl)
        p(1,1,1:nptl) = 0.0_wp
        p(1,2,1:nptl) = 0.0_wp
        p(2,1,1:nptl) = 0.0_wp
        p(2,2,1:nptl) = 1.0_wp

      END SUBROUTINE pde_coef

      SUBROUTINE bound_cond(t,u,du_dx,bound,beta,gamma,finish,v,vdot,i_comm, &
         r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Scalar Arguments ..
        INTEGER, INTENT (IN) :: bound
        INTEGER, INTENT (OUT) :: finish
        REAL (wp), INTENT (IN) :: t
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), INTENT (OUT) :: beta(:), gamma(:)
        REAL (wp), INTENT (IN) :: du_dx(:), u(:)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:), v(:), vdot(:)
        ! .. Executable Statements ..

        finish = 1
        beta(1:2) = (/ 1.0_wp, 0.0_wp/)
        gamma(1) = 0.0_wp

        IF (bound==0) THEN
          gamma(2) = u(1) - 1.0_wp
        ELSE
          gamma(2) = u(1) + 1.0_wp
        END IF

      END SUBROUTINE bound_cond

      SUBROUTINE init_value(x,u,v,i_comm,r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC SIN, SIZE
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:)
        REAL (wp), INTENT (OUT) :: u(:,:)
        REAL (wp), OPTIONAL, INTENT (OUT) :: v(:)
        REAL (wp), INTENT (IN) :: x(:)
        ! .. Local Scalars ..
        INTEGER :: npts_pde
        REAL (wp) :: piby2
        ! .. Executable Statements ..
        piby2 = r_comm(1)
        npts_pde = SIZE(x)

        u(1,1:npts_pde) = -SIN(piby2*x(1:npts_pde))
        u(2,1:npts_pde) = -piby2*piby2*u(1,1:npts_pde)

      END SUBROUTINE init_value

    END MODULE nag_pde_parab_1d_ex02_mod
    PROGRAM nag_pde_parab_1d_ex02

      ! Example Program Text for nag_pde_parab_1d
      ! NAG fl90, Release 4. NAG Copyright 2000.

      ! .. Use Statements ..
      USE nag_examples_io, ONLY : nag_std_out
      USE nag_math_constants, ONLY : nag_pi
      USE nag_pde_parab_1d, ONLY : nag_pde_parab_1d_comm_wp => &
       nag_pde_parab_1d_comm_dp, nag_pde_parab_1d_cntrl_wp => &
       nag_pde_parab_1d_cntrl_dp, nag_pde_parab_1d_cntrl_init, &
       nag_pde_parab_1d_coll, nag_pde_interp_1d_coll, nag_deallocate
      USE nag_pde_parab_1d_ex02_mod, ONLY : pde_coef, bound_cond, init_value, &
       wp
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC REAL
      ! .. Parameters ..
      INTEGER, PARAMETER :: deg_poly = 3, itype = 1, ninterp = 6, npde = 2, &
       npts_bk = 10
      INTEGER, PARAMETER :: nel = npts_bk - 1
      INTEGER, PARAMETER :: npts = nel*deg_poly + 1
      REAL (wp), PARAMETER :: x_interp(ninterp) = (/ -1.0_wp, -0.6_wp, &
       -0.2_wp, 0.2_wp, 0.6_wp, 1.0_wp/)
      ! .. Local Scalars ..
      INTEGER :: i, it, num_jac_eval, num_resid_eval, num_time_iter, &
       num_time_step
      REAL (wp) :: acc, pi, piby2, t_end, t_start, npts_bk_1
      LOGICAL :: first_call
      TYPE (nag_pde_parab_1d_comm_wp) :: comm_ode
      TYPE (nag_pde_parab_1d_cntrl_wp) :: control
      ! .. Local Arrays ..
      REAL (wp) :: abs_loc_tol(npde*npts), rel_loc_tol(npde*npts), r_comm(1), &
       u(npde*npts), u_interp(npde,ninterp,itype), x_bkpts(npts_bk)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) &
       ' Example Program Results for nag_pde_parab_1d_ex02 '

      pi = nag_pi(0.0_wp)
      piby2 = 0.5_wp*pi

      ! Set initial conditions and accuracy

      t_start = 0.0_wp
      acc = 1.0E-4_wp
      t_end = 0.1E-4_wp
      first_call = .TRUE.
      r_comm(1) = piby2
      rel_loc_tol = acc
      abs_loc_tol = acc

      ! Set the break-points

      npts_bk_1 = REAL((npts_bk-1),kind=wp)
      DO i = 1, npts_bk
        x_bkpts(i) = -1.0_wp + 2.0_wp*REAL((i-1),kind=wp)/npts_bk_1
      END DO

      WRITE (nag_std_out,999) deg_poly, nel
      WRITE (nag_std_out,998) acc, npts
      WRITE (nag_std_out,997) x_interp(1:ninterp)

      ! Initialize the structure control and set required control parameters

      CALL nag_pde_parab_1d_cntrl_init(control)

      control%l2_norm = .TRUE.

      DO it = 1, 5
        t_end = 10.0_wp*t_end

        CALL nag_pde_parab_1d_coll(pde_coef,bound_cond,init_value,deg_poly, &
         t_start,t_end,x_bkpts,u,first_call=first_call,comm_ode=comm_ode, &
         rel_loc_tol=rel_loc_tol,abs_loc_tol=abs_loc_tol,control=control, &
         num_time_step=num_time_step,num_resid_eval=num_resid_eval, &
         num_jac_eval=num_jac_eval,num_time_iter=num_time_iter,r_comm=r_comm)

        first_call = .FALSE.

        ! Interpolate at required spatial points

        CALL nag_pde_interp_1d_coll(u,x_bkpts,deg_poly,comm_ode,x_interp, &
         u_interp)

        WRITE (nag_std_out,996) t_end, u_interp(1,1:ninterp,1)
        WRITE (nag_std_out,995) u_interp(2,1:ninterp,1)
      END DO

      ! Free structure comm_ode  allocated by NAG fl90

      CALL nag_deallocate(comm_ode)

      ! Print integration statistics

      WRITE (nag_std_out,994) num_time_step, num_resid_eval, num_jac_eval, &
       num_time_iter

999   FORMAT (' Polynomial degree =',I4,'   No. of elements = ',I4)
998   FORMAT (' Accuracy requirement = ',E9.3,'  Number of points = ',I5/)
997   FORMAT ('  T /   X   ',6F8.4/)
996   FORMAT (1X,F6.4,' U(1)',6F8.4)
995   FORMAT (8X,'U(2)',6F8.4/)
994   FORMAT (' Number of integration steps in time                  ', &
       I4/' Number of residual evaluations of resulting ODE sys','tem', &
       I4/' Number of Jacobian evaluations               ','         ', &
       I4/' Number of iterations of nonlinear solve','r              ',I4/)

    END PROGRAM nag_pde_parab_1d_ex02
