    MODULE nag_pde_parab_1d_ex05_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) = -x(1:nptl)*du_dx(1,1:nptl)*v(1)*vdot(1)
        r(1,1:nptl) = du_dx(1,1:nptl)
        p(1,1,1:nptl) = v(1)*v(1)

      END SUBROUTINE pde_coef

      SUBROUTINE ode_coef(t,v,vdot,x_ode,ucp,ducp_dx,rcp,ducp_dt,d2ucp_dtdx,f, &
         finish,g_in_f,i_comm,r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Scalar Arguments ..
        INTEGER, INTENT (OUT) :: finish
        REAL (wp), INTENT (IN) :: t
        LOGICAL, OPTIONAL, INTENT (IN) :: g_in_f
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), INTENT (IN) :: d2ucp_dtdx(:,:), ducp_dt(:,:), ducp_dx(:,:), &
         rcp(:,:), ucp(:,:), v(:), vdot(:), x_ode(:)
        REAL (wp), INTENT (OUT) :: f(:)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:)
        ! .. Executable Statements ..

        IF (g_in_f) THEN
          f(1) = vdot(1) - v(1)*ucp(1,1) - ducp_dx(1,1) - 1.0_wp - t
          finish = 1
        ELSE
          f(1) = vdot(1)
          finish = -1
        END IF

      END SUBROUTINE ode_coef

      SUBROUTINE bound_cond(t,u,du_dx,bound,beta,gamma,finish,v,vdot,i_comm, &
         r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC EXP
        ! .. 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) = 1.0_wp
        IF (bound==0) THEN
          gamma(1) = -v(1)*EXP(t)
        ELSE IF (bound==1) THEN
          gamma(1) = -v(1)*vdot(1)
        END IF

      END SUBROUTINE bound_cond

      SUBROUTINE init_value(x,u,v,i_comm,r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC EXP, 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) :: t_start
        ! .. Executable Statements ..
        t_start = r_comm(1)
        npts_pde = SIZE(x)

        v(1) = t_start
        u(1,1:npts_pde) = EXP(t_start*(1.0_wp-x(1:npts_pde))) - 1.0_wp

      END SUBROUTINE init_value

      SUBROUTINE exact_sol(time,x,u)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC EXP, SIZE
        ! .. Scalar Arguments ..
        REAL (wp), INTENT (IN) :: time
        ! .. Array Arguments ..
        REAL (wp), INTENT (OUT) :: u(:)
        REAL (wp), INTENT (IN) :: x(:)
        ! .. Local Scalars ..
        INTEGER :: npts
        ! .. Executable Statements ..
        npts = SIZE(x)

        u(1:npts) = EXP(time*(1.0_wp-x(1:npts))) - 1.0_wp

      END SUBROUTINE exact_sol

    END MODULE nag_pde_parab_1d_ex05_mod
    PROGRAM nag_pde_parab_1d_ex05

      ! 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_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_deallocate
      USE nag_pde_parab_1d_ex05_mod, ONLY : pde_coef, bound_cond, init_value, &
       ode_coef, exact_sol, wp
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC REAL
      ! .. Parameters ..
      INTEGER, PARAMETER :: deg_poly = 2, npde = 1, npts_bk = 11, npts_ode = 1
      INTEGER, PARAMETER :: npts_pde = (npts_bk-1)*deg_poly + 1
      INTEGER, PARAMETER :: num_ode = 1
      INTEGER, PARAMETER :: neqn = npde*npts_pde + num_ode
      ! .. Local Scalars ..
      INTEGER :: i, it, num_jac_eval, num_resid_eval, num_time_iter, &
       num_time_step
      REAL (wp) :: acc, t_end, t_start
      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(neqn), exact_u(npts_bk), rel_loc_tol(neqn), &
       r_comm(1), u(neqn), x_bkpts(npts_bk), x_ode(npts_ode)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) &
       ' Example Program Results for nag_pde_parab_1d_ex05 '

      acc = 1.0E-4_wp

      ! Set initial conditions

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

      ! Set the break-points

      DO i = 1, npts_bk
        x_bkpts(i) = REAL((i-1),kind=wp)/REAL((npts_bk-1),kind=wp)
      END DO

      x_ode(1) = 1.0_wp

      WRITE (nag_std_out,999) deg_poly, (npts_bk-1)
      WRITE (nag_std_out,996) acc, npts_pde
      WRITE (nag_std_out,998) x_bkpts(1), x_bkpts(3), x_bkpts(5), x_bkpts(7), &
       x_bkpts(11)

      ! Initialize the structure control and set required control parameters

      CALL nag_pde_parab_1d_cntrl_init(control)

      control%l2_norm = .TRUE.
      control%matrix_type = 'F'

      DO it = 1, 5
        t_end = 0.1_wp*(2.0_wp**it)

        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,ode_coef=ode_coef, &
         num_ode=num_ode,x_ode=x_ode,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.

        ! Check against the exact solution

        CALL exact_sol(t_end,x_bkpts,exact_u)

        WRITE (nag_std_out,997) t_start
        WRITE (nag_std_out,994) u(1), u(5), u(9), u(13), u(21), u(22)
        WRITE (nag_std_out,993) exact_u(1), exact_u(3), exact_u(5), &
         exact_u(7), exact_u(11), t_start
      END DO

      ! Free structure comm_ode  allocated by NAG fl90

      CALL nag_deallocate(comm_ode)

      ! Print integration statistics

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

999   FORMAT (' Degree of Polynomial =',I4,'   No. of elements =',I4/)
998   FORMAT ('  X        ',5F9.3/)
997   FORMAT (' T = ',F6.3)
996   FORMAT (//'  Simple coupled PDE using BDF '/'  Accuracy require', &
       'ment =',E10.3,' Number of points = ',I4/)
995   FORMAT (' Number of integration steps in time = ',I6/' Number o', &
       'f function evaluations = ',I6/' Number of Jacobian eval','uations =', &
       I6/' Number of iterations = ',I6/)
994   FORMAT (1X,'App.  sol.  ',F7.3,4F9.3,'  ODE sol. =',F8.3)
993   FORMAT (1X,'Exact sol.  ',F7.3,4F9.3,'  ODE sol. =',F8.3/)

    END PROGRAM nag_pde_parab_1d_ex05
