    MODULE nag_pde_parab_1d_ex04_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 ..

        finish = 1
        eps = r_comm(1)

        q(1) = u(1)*du_dx(1)
        r(1) = eps*du_dx(1)
        p(1,1) = 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
        ! .. 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) :: a, b, c, eps, ue, x
        ! .. Executable Statements ..

        eps = r_comm(1)
        finish = 1

        beta(1) = 0.0_wp
        IF (bound==0) THEN
          x = 0.0_wp
          a = (x-0.25_wp-0.75_wp*t)/(4.0_wp*eps)
          b = (0.9_wp*x-0.325_wp-0.495_wp*t)/(2.0_wp*eps)
          IF (a>0.0_wp .AND. a>b) THEN
            a = EXP(-a)
            c = (0.8_wp*x-0.4_wp-0.24_wp*t)/(4.0_wp*eps)
            c = EXP(c)
            ue = (0.5_wp+0.1_wp*c+a)/(1.0_wp+c+a)
          ELSE IF (b>0.0_wp .AND. b>=a) THEN
            b = EXP(-b)
            c = (-0.8_wp*x+0.4_wp+0.24_wp*t)/(4.0_wp*eps)
            c = EXP(c)
            ue = (0.1_wp+0.5_wp*c+b)/(1.0_wp+c+b)
          ELSE
            a = EXP(a)
            b = EXP(b)
            ue = (1.0_wp+0.5_wp*a+0.1_wp*b)/(1.0_wp+a+b)
          END IF
        ELSE
          x = 1.0_wp
          a = (x-0.25_wp-0.75_wp*t)/(4.0_wp*eps)
          b = (0.9_wp*x-0.325_wp-0.495_wp*t)/(2.0_wp*eps)

          IF (a>0.0_wp .AND. a>b) THEN
            a = EXP(-a)
            c = (0.8_wp*x-0.4_wp-0.24_wp*t)/(4.0_wp*eps)
            c = EXP(c)
            ue = (0.5_wp+0.1_wp*c+a)/(1.0_wp+c+a)
          ELSE IF (b>0.0_wp .AND. b>=a) THEN
            b = EXP(-b)
            c = (-0.8_wp*x+0.4_wp+0.24_wp*t)/(4.0_wp*eps)
            c = EXP(c)
            ue = (0.1_wp+0.5_wp*c+b)/(1.0_wp+c+b)
          ELSE
            a = EXP(a)
            b = EXP(b)
            ue = (1.0_wp+0.5_wp*a+0.1_wp*b)/(1.0_wp+a+b)
          END IF
        END IF
        gamma(1) = u(1) - ue

      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 :: i, npts_pde
        REAL (wp) :: a, b, c, eps, t
        ! .. Executable Statements ..
        eps = r_comm(1)
        npts_pde = SIZE(x_pde)

        t = 0.0_wp
        DO i = 1, npts_pde
          a = (x_pde(i)-0.25_wp-0.75_wp*t)/(4.0_wp*eps)
          b = (0.9_wp*x_pde(i)-0.325_wp-0.495_wp*t)/(2.0_wp*eps)
          IF (a>0.0_wp .AND. a>b) THEN
            a = EXP(-a)
            c = (0.8_wp*x_pde(i)-0.4_wp-0.24_wp*t)/(4.0_wp*eps)
            c = EXP(c)
            u(1,i) = (0.5_wp+0.1_wp*c+a)/(1.0_wp+c+a)
          ELSE IF (b>0.0_wp .AND. b>=a) THEN
            b = EXP(-b)
            c = (-0.8_wp*x_pde(i)+0.4_wp+0.24_wp*t)/(4.0_wp*eps)
            c = EXP(c)
            u(1,i) = (0.1_wp+0.5_wp*c+b)/(1.0_wp+c+b)
          ELSE
            a = EXP(a)
            b = EXP(b)
            u(1,i) = (1.0_wp+0.5_wp*a+0.1_wp*b)/(1.0_wp+a+b)
          END IF
        END DO

      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

      FUNCTION exact_sol(t,x,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 :: i, n
        REAL (wp) :: a, b, c
        ! .. Executable Statements ..
        n = SIZE(x)

        DO i = 1, n
          a = (x(i)-0.25_wp-0.75_wp*t)/(4.0_wp*eps)
          b = (0.9_wp*x(i)-0.325_wp-0.495_wp*t)/(2.0_wp*eps)
          c = (0.8_wp*x(i)-0.4_wp-0.24_wp*t)/(4.0_wp*eps)
          IF (a>0.0_wp .AND. a>b) THEN
            a = EXP(-a)
            b = EXP(c)
            c = (0.5_wp+0.1_wp*b+a)
          ELSE IF (b>0.0_wp .AND. b>=a) THEN
            b = EXP(-b)
            a = EXP(-c)
            c = (0.1_wp+0.5_wp*a+b)
          ELSE
            a = EXP(a)
            b = EXP(b)
            c = (1.0_wp+0.5_wp*a+0.1_wp*b)
          END IF
          exact_sol(i) = c/(1.0_wp+a+b)
        END DO

      END FUNCTION exact_sol

    END MODULE nag_pde_parab_1d_ex04_mod
    PROGRAM nag_pde_parab_1d_ex04

      ! 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_pde_interp_1d_fd, nag_deallocate
      USE nag_pde_parab_1d_ex04_mod, ONLY : pde_coef, bound_cond, init_value, &
       remesh_fun, exact_sol, wp
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC REAL
      ! .. Parameters ..
      INTEGER, PARAMETER :: itype = 1, ninterp = 5, npde = 1, npts_pde = 61
      INTEGER, PARAMETER :: neqn = npde*npts_pde
      ! .. Local Scalars ..
      INTEGER :: i, it, num_jac_eval, num_resid_eval, num_time_iter, &
       num_time_step, n_remesh
      REAL (wp) :: acc, dx_mesh, eps, t_end, t_start
      LOGICAL :: first_call, remesh
      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(ninterp), rel_loc_tol(neqn), &
       r_comm(1), u(neqn), u_interp(npde,ninterp,itype), x_interp(ninterp), &
       x_pde(npts_pde)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) &
       ' Example Program Results for nag_pde_parab_1d_ex04 '

      ! Set initial conditions and accuracy

      acc = 0.5E-4_wp
      eps = 0.005_wp
      t_start = 0.0_wp
      remesh = .TRUE.
      n_remesh = 3
      dx_mesh = 0.5_wp
      first_call = .TRUE.
      r_comm(1) = eps
      rel_loc_tol = acc
      abs_loc_tol = acc

      WRITE (nag_std_out,998) 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

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

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

        first_call = .FALSE.

        SELECT CASE (it)
        CASE (1)
          x_interp = (/ 0.3_wp, 0.4_wp, 0.5_wp, 0.6_wp, 0.7_wp/)
        CASE (2)
          x_interp = (/ 0.4_wp, 0.5_wp, 0.6_wp, 0.7_wp, 0.8_wp/)
        CASE (3)
          x_interp = (/ 0.6_wp, 0.65_wp, 0.7_wp, 0.75_wp, 0.8_wp/)
        CASE (4)
          x_interp = (/ 0.7_wp, 0.75_wp, 0.8_wp, 0.85_wp, 0.9_wp/)
        CASE (5)
          x_interp = (/ 0.8_wp, 0.85_wp, 0.9_wp, 0.95_wp, 1.0_wp/)
        END SELECT

        ! Interpolate at required spatial points

        CALL nag_pde_interp_1d_fd(u,x_pde,x_interp,u_interp)

        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)

        ! Exact solution

        exact_u = exact_sol(t_start,x_interp,eps)

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

      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 (' 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)

    END PROGRAM nag_pde_parab_1d_ex04
