    MODULE nag_pde_parab_1d_ex07_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(:)
        ! .. Executable Statements ..

        finish = 1
        q(1) = -x*du_dx(1)*v(1)*vdot(1)
        r(1) = du_dx(1)
        p(1,1) = v(1)*v(1)

      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
        ! .. 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
          gamma(1) = -v(1)*vdot(1)
        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 EXP, 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) :: t_start
        ! .. Executable Statements ..
        t_start = r_comm(1)
        npts_pde = SIZE(x_pde)

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

      END SUBROUTINE init_value

      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 exact_sol(t,x_pde,exact_u)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC EXP, SIZE
        ! .. Scalar Arguments ..
        REAL (wp), INTENT (IN) :: t
        ! .. Array Arguments ..
        REAL (wp), INTENT (OUT) :: exact_u(:)
        REAL (wp), INTENT (IN) :: x_pde(:)
        ! .. Local Scalars ..
        INTEGER :: npts_pde
        ! .. Executable Statements ..
        npts_pde = SIZE(x_pde)

        exact_u(1:npts_pde) = EXP(t*(1.0_wp-x_pde(1:npts_pde))) - 1.0_wp

      END SUBROUTINE exact_sol

    END MODULE nag_pde_parab_1d_ex07_mod
    PROGRAM nag_pde_parab_1d_ex07

      ! 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_fd, nag_deallocate
      USE nag_pde_parab_1d_ex07_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 :: npde = 1, npts_ode = 1, npts_pde = 21, 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(npde*npts_pde), &
       rel_loc_tol(neqn), r_comm(1), u(neqn), x_ode(npts_ode), x_pde(npts_pde)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) &
       ' Example Program Results for nag_pde_parab_1d_ex07 '

      ! Set initial conditions and accuracy

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

      WRITE (nag_std_out,997) acc, npts_pde

      ! Set spatial mesh points

      DO i = 1, npts_pde
        x_pde(i) = REAL((i-1),kind=wp)/REAL((npts_pde-1),kind=wp)
      END DO

      x_ode(1) = 1.0_wp

      WRITE (nag_std_out,999) x_pde(1), x_pde(5), x_pde(9), x_pde(13), &
       x_pde(21)

      ! Initialize the structure control and set required control parameters

      CALL nag_pde_parab_1d_cntrl_init(control)

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

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

        CALL nag_pde_parab_1d_fd(pde_coef,bound_cond,init_value,t_start,t_end, &
         x_pde,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.

        CALL exact_sol(t_end,x_pde,exact_u)

        WRITE (nag_std_out,998) t_start
        WRITE (nag_std_out,995) u(1), u(5), u(9), u(13), u(21), u(22)
        WRITE (nag_std_out,994) exact_u(1), exact_u(5), exact_u(9), &
         exact_u(13), exact_u(21), t_start
      END DO

      ! Free structure comm_ode  allocated by NAG fl90

      CALL nag_deallocate(comm_ode)

      ! Print integration statistics

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

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

    END PROGRAM nag_pde_parab_1d_ex07
