    MODULE nag_pde_parab_1d_ex06_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) :: eps
        ! .. Executable Statements ..
        eps = r_comm(2)

        finish = 1
        q(1) = -x*du_dx(1)*v(1)*vdot(1)
        r(1) = eps*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(:)
        ! .. Local Scalars ..
        REAL (wp) :: eps, inveps
        ! .. Executable Statements ..
        eps = r_comm(2)
        inveps = 1.0_wp/eps

        finish = 1

        IF (bound==0) THEN
          beta(1) = 1.0_wp
          gamma(1) = -v(1)*EXP(inveps*t)
        ELSE
          beta(1) = 1.0_wp
          gamma(1) = -v(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) :: eps, inveps, t_start
        ! .. Executable Statements ..
        t_start = r_comm(1)
        npts_pde = SIZE(x_pde)
        eps = r_comm(2)
        inveps = 1.0_wp/eps

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

      END SUBROUTINE init_value

      SUBROUTINE remesh_fun(t,x_pde,u,r,fmon,i_comm,r_comm)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC ABS, MAX, MIN, SIZE
        ! .. Scalar Arguments ..
        REAL (wp), INTENT (IN) :: t
        ! .. Array Arguments ..
        INTEGER, OPTIONAL, INTENT (IN) :: i_comm(:)
        REAL (wp), INTENT (OUT) :: fmon(:)
        REAL (wp), INTENT (IN) :: r(:,:), u(:,:), x_pde(:)
        REAL (wp), OPTIONAL, INTENT (IN) :: r_comm(:)
        ! .. Local Scalars ..
        INTEGER :: i, k, l, npts_pde
        REAL (wp) :: drdx, h
        ! .. Executable Statements ..
        npts_pde = SIZE(x_pde)

        DO i = 1, npts_pde - 1
          k = MAX(1,i-1)
          l = MIN(npts_pde,i+1)
          h = 0.5_wp*(x_pde(l)-x_pde(k))
          drdx = (r(1,i+1)-r(1,i))/h
          fmon(i) = ABS(drdx)
        END DO
        fmon(npts_pde) = fmon(npts_pde-1)

      END SUBROUTINE remesh_fun

      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(:)
        ! .. Local Scalars ..
        REAL (wp) :: eps
        ! .. Executable Statements ..
        eps = r_comm(2)

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

      END SUBROUTINE ode_coef

      FUNCTION exact_sol(x,t,eps)
        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC EXP, SIZE
        ! .. Scalar Arguments ..
        REAL (wp), INTENT (IN) :: eps, t
        ! .. Array Arguments ..
        REAL (wp), INTENT (IN) :: x(:)
        ! .. Function Return Value ..
        REAL (wp) :: exact_sol(SIZE(x))
        ! .. Local Scalars ..
        INTEGER :: n
        ! .. Executable Statements ..
        n = SIZE(x)

        exact_sol(1:n) = 1.0_wp - x(1:n)
        exact_sol = EXP(exact_sol*(t/eps)) - 1.0_wp

      END FUNCTION exact_sol

    END MODULE nag_pde_parab_1d_ex06_mod
    PROGRAM nag_pde_parab_1d_ex06

      ! 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_cntrl_wp => &
       nag_pde_parab_1d_cntrl_dp, nag_pde_parab_1d_cntrl_init, &
       nag_pde_parab_1d_fd, nag_pde_interp_1d_fd
      USE nag_pde_parab_1d_ex06_mod, ONLY : pde_coef, bound_cond, init_value, &
       ode_coef, remesh_fun, exact_sol, wp
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC EPSILON, MAX, PRECISION, REAL, SIZE, SQRT, SUM
      ! .. Parameters ..
      INTEGER, PARAMETER :: itype = 1, ninterp = 5, npde = 1, npts_ode = 1, &
       npts_pde = 201, num_ode = 1
      INTEGER, PARAMETER :: neqn = npde*npts_pde + num_ode
      REAL (wp), PARAMETER :: x_interp(ninterp) = (/ 0.3_wp, 0.4_wp, 0.5_wp, &
       0.6_wp, 0.7_wp/)
      ! .. Local Scalars ..
      INTEGER :: i, it, num_jac_eval, num_resid_eval, num_time_iter, &
       num_time_step, n_remesh
      REAL (wp) :: acc, bound_int, dx_mesh, eps, l2norm, l2norm0, t_end, &
       t_start
      LOGICAL :: first_call, remesh
      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(2), u(neqn), u_interp(npde,ninterp,itype), &
       x_ode(npts_ode), x_pde(npts_pde)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) &
       ' Example Program Results for nag_pde_parab_1d_ex06 '

      acc = MAX(SQRT(EPSILON(acc)),0.5E-4_wp)

      ! Set initial conditions

      remesh = .TRUE.
      n_remesh = 3
      dx_mesh = 0.5_wp
      first_call = .TRUE.

      rel_loc_tol = acc
      abs_loc_tol = acc

      WRITE (nag_std_out,998) acc, npts_pde
      WRITE (nag_std_out,993) n_remesh

      ! Initialize the structure control and set required control parameters

      CALL nag_pde_parab_1d_cntrl_init(control)

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

      bound_int = 2.0_wp/REAL(npts_pde-1,kind=wp)
      IF (PRECISION(bound_int)<=6) bound_int = 8.0_wp/REAL(npts_pde-1,kind=wp)

      DO it = 1, 5

        SELECT CASE (it)
        CASE (1)
          eps = 1.0_wp
        CASE (2)
          eps = 0.1_wp
        CASE (3)
          eps = 0.08_wp
        CASE (4)
          eps = 0.06_wp
        CASE (5)
          eps = 0.04_wp
        END SELECT

        t_start = 1.0E-4_wp

        ! 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
        r_comm(1) = t_start
        r_comm(2) = eps

        WRITE (nag_std_out,992) eps
        WRITE (nag_std_out,*)

        t_end = 0.2_wp

        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,rel_loc_tol=rel_loc_tol,abs_loc_tol=abs_loc_tol, &
         remesh=remesh,n_remesh=n_remesh,dx_mesh=dx_mesh,control=control, &
         remesh_fun=remesh_fun,num_time_step=num_time_step, &
         bound_int=bound_int,num_resid_eval=num_resid_eval, &
         num_jac_eval=num_jac_eval,num_time_iter=num_time_iter,r_comm=r_comm)

        ! Interpolate at required spatial points

        CALL nag_pde_interp_1d_fd(u,x_pde,x_interp,u_interp)

        exact_u(1:ninterp) = exact_sol(x_interp(1:ninterp),t_start,eps)

        WRITE (nag_std_out,999) t_start
        WRITE (nag_std_out,996) x_interp(1:ninterp)

        WRITE (nag_std_out,995) u_interp(1,1:ninterp,1)
        WRITE (nag_std_out,994) exact_u(1:ninterp)

        exact_u(1:SIZE(x_pde)) = exact_sol(x_pde,t_start,eps)

        l2norm = SUM((u(1:npde*npts_pde)-exact_u(1:npde*npts_pde))**2)
        l2norm = SQRT(l2norm)

        l2norm0 = SUM((exact_u(1:npde*npts_pde))**2)
        l2norm0 = SQRT(l2norm0)

        l2norm = l2norm/l2norm0

        WRITE (nag_std_out,990) 'L2norm relative error = ', l2norm

        ! Print integration statistics

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

      END DO

999   FORMAT (' T = ',F6.3)
998   FORMAT (//'  Accuracy requirement =',E10.3,' Number of points = ',I3/)
997   FORMAT (' Number of integration steps in time = ',I6/' Number o', &
       'f function evaluations = ',I6/' Number of Jacobian eval','uations =', &
       I6/' Number of iterations = ',I6/)
996   FORMAT (1X,'X           ',5F9.4)
995   FORMAT (1X,'Approx sol. ',5F9.4)
994   FORMAT (1X,'Exact  sol. ',5F9.4/)
993   FORMAT (2X,'Remeshing every',I3,' time steps'/)
992   FORMAT (2X,'E =',F8.3)
991   FORMAT (1X,F9.4,1X,F9.4)
990   FORMAT (1X,A,1P,E8.1)

    END PROGRAM nag_pde_parab_1d_ex06
