! Integrate Oxygen Diffusion Model found in the book
! Crank, John. Free and Moving Boundary Problems,
! Oxford Univ. Press, (1984), p. 19-20 and p. 261-262.
      use feynman_kac_int
      use hqsval_int
      use mp_types
      use norm_int
      use umach_int
      implicit none

      integer :: nout
      real(kind(1e0)), allocatable :: xgrid(:), tgrid(:), y(:,:),&
                                      yprime(:,:), f(:,:), s(:)
      real(kind(1e0)) :: dx, x_min, x_max, zero=0.0e0, one=1.0e0
      real(kind(1e0)) :: atol(1), rtol(1)
      type(s_fcn_data) :: fk_ox2
      integer :: i, nint, n, nunit, ntimes = 8
      integer, parameter :: ndeg = 6, nlbc = 1, nrbc = 2
      real(kind(1e0)), external :: fkcoef_ox2, fkinitcond_ox2
      external fkbc_ox2, fkforce_ox2
      
      call umach(2,nout)

! Define number of equally spaced intervals for elements
      nint = 100
! Allocate the space needed for the integration process
      n = 3*(nint+1)
      allocate(xgrid(nint+1), y(n,0:ntimes), yprime(n,0:ntimes),&
               tgrid(ntimes), f(1,ntimes), s(ntimes))

! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
      allocate(fk_ox2 % rdata (1), fk_ox2 % idata (1))

      atol(1) = 0.5e-2
      rtol(1) = 0.5e-2
      fk_ox2 % rdata(1) = atol(1)
      fk_ox2 % idata(1) = ndeg

! Define interval endpoints
      x_min = zero
      x_max = one
! Define interval widths
      dx = (x_max-x_min)/real(nint)
      xgrid(1) = x_min
      xgrid(nint+1) = x_max
! Define grid points of interval
      do i=2,nint
        xgrid(i) = xgrid(i-1) + dx
      end do
! Define time integration output points
! These correspond to published values in Crank's book, p. 261-262 
      tgrid = (/0.04e0,0.06e0,0.10e0,0.12e0,0.14e0,&
                0.16e0,0.18e0,0.185e0/)
      call feynman_kac (xgrid, tgrid, nlbc, nrbc, fkcoef_ox2,&
                        fkinitcond_ox2, fkbc_ox2, y, yprime,&
                        ATOL = atol, RTOL = rtol,&
                        FKFORCE = fkforce_ox2, FCN_DATA = fk_ox2)
! Summarize output at the left end
      do i=1,ntimes
        f(:,i)= hqsval ((/zero/), xgrid, y(:,i))
      end do
! Check differences of evaluation and published left end values
      f(1,:) = f(1,:) - (/2.743e-01, 2.236e-01, 1.432e-01,&
           1.091e-01, 7.786e-02, 4.883e-02, 2.179e-02, 1.534e-02/)
      write(nout,*) "Oxygen Depletion Model, from Crank's "//&
                    "Book, p. 261-262,"
      write(nout,*) "'Free and Moving Boundary Value Problems'"
      if (norm(f(1,:)) < ntimes * atol(1)) then
        write(nout,*) "FEYNMAN_KAC Example 6 - Fixed Sealed "//&
                      "Surface Values are correct"
      else
        write(nout,*) "FEYNMAN_KAC Example 6 - does not agree with"//&
                      " published left end values"
      end if
! Define known position of free boundary at the time points
      s = (/0.9992e0,0.9918e0,0.9350e0,0.8792e0,&
            0.7989e0,0.6834e0,0.5011e0,0.4334e0/)
                       
! Evaluate and verify solution is small near free boundary -
      do i=1,ntimes
        f(:,i) = hqsval ((/s(i)/), xgrid, y(:,i))
      end do

      if (norm(f(1,:)) < ntimes * atol(1)) then
        write(nout,*) "FEYNMAN_KAC Example 6 - Free Boundary "//&
                   "Position Values are correct"
      else
        write(nout,*) "FEYNMAN_KAC Example 6 - does not agree "//&
                   "with published free boundary values"
      end if
      end

      function fkcoef_ox2 (x, tx, iflag, fk_ox2) result(value)
      use mp_types
      implicit none
! Coefficient valuation routine for the Oxygen Diffusion
! Model found in Crank's book, p. 20
! Input/Ouput variables
      integer, intent(inout) :: iflag
      real(kind(1e0)), intent(in) :: x, tx
      type(s_fcn_data), optional :: fk_ox2
      real(kind(1e0)) :: value

! Local variables
      real(kind(1e0)) :: zero = 0.0e0, two = 2.0e0         

      select case (iflag)
        case (1) ! Factor DSigma/Dx(x,t)
          value = zero
        case (2) !  Factor Sigma(x,t)
          value = sqrt(two)
        case (3) !  Factor Mu (x,t)
          value = zero
        case (4) !  Factor Kappa (x,t)
          value = zero
      end select
! Signal no dependence on tx=t=time for any coefficient.
      iflag = 0
      return
      end function fkcoef_ox2

      function fkinitcond_ox2(x, fk_ox2) result(value)
      use mp_types
      implicit none
      real(kind(1e0)), intent(in) :: x
      type (s_fcn_data), optional :: fk_ox2
      real(kind(1e0)) :: value

      real(kind(1e0)) :: half = 0.5e0, one = 1.0e0

      value = half * (one - x)**2
      return           
      end function fkinitcond_ox2

      subroutine fkbc_ox2 (tx, iflag, bccoefs, fk_ox2)
      use mp_types
      implicit none
! Evaluation routine for Oxygen Diffusion Model 
! boundary conditions.  
! Input/Ouput variables
      real(kind(1e0)), intent(in) :: tx
      integer, intent(inout) :: iflag
      real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
      type (s_fcn_data), optional :: fk_ox2

! Local variables
      real(kind(1e0)) :: zero = 0.0e0, one = 1.0e0, atol

      atol = fk_ox2 % rdata(1)
      select case (iflag)
        case (1) ! Left Boundary Condition,  at X_min=0
! There is a rapid blending of the boundary condition to achieve
! a zero derivative value at the left end.
! The initial data has the derivative with value one.
! This boundary condition essentially abruptly changes that
! derivative value to zero.
! Returning iflag=1 signals time dependence.  This is important
! for this problem.
          bccoefs(1,1:4) = (/0.0e0, one, 0.0e0, exp(tx/atol**2)/)
          return                     
        case (2) ! Right Boundary Condition, at X_max=1
          bccoefs(1,1:4) = (/one, 0.0e0, 0.0e0, 0.0e0/)
          bccoefs(2,1:4) = (/0.0e0, one, 0.0e0, 0.0e0/)
      end select 
      iflag = 0 ! Signal no dependence on tx=time.
      end subroutine fkbc_ox2

      subroutine fkforce_ox2 (interval, t, hx, y, xlocal, qw, u,&
                              phi, dphi, fk_ox2)
! Evaluation routine for Oxygen Diffusion model forcing function.
      use mp_types
      implicit none
      integer, parameter :: local = 6
      integer :: i, j, l, mu, ndeg
      integer, intent(in) :: interval 

      real(kind(1e0)), intent(in) :: y(:), t, hx, qw(:),&
                                     xlocal(:), u(:,:)
      real(kind(1e0)), intent(out) :: phi(:), dphi(:,:)
      type (s_fcn_data), optional :: fk_ox2

      real(kind(1e0)) :: yl(local), bf(local)
      real(kind(1e0)) :: value, zero = 0.0e0, one = 1.0e0, rt

      yl = y(3*interval-2:3*interval+3)
      phi = zero

      value = fk_ox2 % rdata(1)
      ndeg = fk_ox2 % idata(1)

      mu = 2
      do j=1,local
        do l=1,ndeg
          bf(1:3) = u(l,1:3)
          bf(4:6) = u(l,7:9) 
          rt = dot_product(yl,bf)
          rt = one - (value/(rt + value))**mu              
          phi(j) = phi(j) + qw(l) * bf(j) * RT
        end do
      end do
      phi = phi * hx
! This is the local derivative matrix for the forcing term -        
      dphi = zero
      do j=1,local
        do i = 1,local
          do l=1,ndeg
            bf(1:3) = u(l,1:3)
            bf(4:6) = u(l,7:9)
            rt = dot_product(yl,bf)
            rt = one/(rt + value)
            dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j) *&
                        rt**(mu+1)
          end do
        end do
      end do
      dphi = mu * dphi * hx * value**mu
      return                                
      end subroutine fkforce_ox2

