    MODULE nag_pde_parab_1d_ex01_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
        ! .. Scalar Arguments ..
        INTEGER, INTENT (OUT) :: finish
        REAL (wp), INTENT (IN) :: t, x
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), INTENT (IN) :: du_dx(:), u(:)
        REAL (wp), INTENT (OUT) :: p(:,:), q(:), r(:)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:), v(:), vdot(:)
        ! .. Local Scalars ..
        REAL (wp) :: alpha
        ! .. Executable Statements ..
        alpha = r_comm(1)

        finish = 1

        q(1) = 4.0_wp*alpha*(u(2)+x*du_dx(2))
        q(2) = 0.0_wp
        r(1) = x*du_dx(1)
        r(2) = du_dx(2) - u(1)*u(2)
        p(1,1:2) = 0.0_wp
        p(2,1) = 0.0_wp
        p(2,2) = 1.0_wp - x*x

      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

        IF (bound==0) THEN
          beta(1:2) = (/ 0.0_wp, 1.0_wp/)
          gamma(1) = u(1)
          gamma(2) = -u(1)*u(2)
        ELSE
          beta(1:2) = (/ 1.0_wp, 0.0_wp/)
          gamma(1:2) = -u(1:2)
        END IF

      END SUBROUTINE bound_cond

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

        u(1,1:npts_pde) = 2.0_wp*alpha*x_pde(1:npts_pde)
        u(2,1:npts_pde) = 1.0_wp

      END SUBROUTINE init_value

    END MODULE nag_pde_parab_1d_ex01_mod
    PROGRAM nag_pde_parab_1d_ex01

      ! 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_fd, nag_pde_interp_1d_fd, nag_deallocate
      USE nag_pde_parab_1d_ex01_mod, ONLY : pde_coef, bound_cond, init_value, &
       wp
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC EPSILON, MAX, REAL, SIN, SQRT
      ! .. Parameters ..
      INTEGER, PARAMETER :: itype = 1, ninterp = 6, npde = 2, npts_pde = 20
      REAL (wp), PARAMETER :: x_interp(ninterp) = (/ 0.0_wp, 0.4_wp, 0.6_wp, &
       0.8_wp, 0.9_wp, 1.0_wp/)
      ! .. Local Scalars ..
      INTEGER :: i, it, num_jac_eval, num_resid_eval, num_time_iter, &
       num_time_step
      REAL (wp) :: acc, alpha, hx, pi, piby2, t_end, t_start
      LOGICAL :: first_call
      CHARACTER (1) :: coord_sys
      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_pde), rel_loc_tol(npde*npts_pde), &
       r_comm(1), u(npde*npts_pde), u_interp(npde,ninterp,itype), &
       x_pde(npts_pde)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) &
       ' Example Program Results for nag_pde_parab_1d_ex01 '

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

      ! Set initial conditions

      alpha = 1.0_wp
      acc = MAX(SQRT(EPSILON(alpha)),1.0E-4_wp)
      hx = piby2/REAL((npts_pde-1),kind=wp)
      t_start = 0.0_wp
      t_end = 0.1E-4_wp
      first_call = .TRUE.
      coord_sys = 'P'

      r_comm(1) = alpha
      rel_loc_tol = acc
      abs_loc_tol = acc

      ! Set spatial mesh points

      x_pde(1) = 0.0_wp
      x_pde(npts_pde) = 1.0_wp

      DO i = 2, npts_pde - 1
        x_pde(i) = SIN(hx*REAL((i-1),kind=wp))
      END DO

      WRITE (nag_std_out,999) acc, alpha
      WRITE (nag_std_out,998) (x_interp(i),i=1,6)

      ! 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_fd(pde_coef,bound_cond,init_value,t_start,t_end, &
         x_pde,u,first_call=first_call,comm_ode=comm_ode,coord_sys=coord_sys, &
         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_fd(u,x_pde,x_interp,u_interp, &
         coord_sys=coord_sys)

        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,997) num_time_step, num_resid_eval, num_jac_eval, &
       num_time_iter

999   FORMAT (//' Accuracy requirement  = ',E12.5/' Parameter ALPHA =', &
       '       ',E12.3/)
998   FORMAT ('   T  /  X  ',6F8.4/)
997   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/)
996   FORMAT (1X,F6.4,' U(1)',6F8.4)
995   FORMAT (8X,'U(2)',6F8.4/)

    END PROGRAM nag_pde_parab_1d_ex01
