! Compute Constant Elasticity of Variance Model for Vanilla Call
      use feynman_kac_int
      use hqsval_int
      use mp_types
      use umach_int

      implicit none

! The set of strike prices
      real(kind(1e0)) :: ks(3)=(/15.0e0,20.0e0,25.0e0/)
! The set of sigma values
      real(kind(1e0)) :: sigma(3) = (/0.2e0, 0.3e0, 0.4e0/)
! The set of model diffusion powers
      real(kind(1e0)) :: alpha(3) = (/2.0e0,1.0e0,0.0e0/)
! Time values for the options
      integer, parameter :: nt = 3
      real(kind(1e0)) :: time(nt)=(/1.e0/12., 4.e0/12., 7.e0/12./)
! Values of the underlying where evaluation are made
      integer, parameter :: nv = 3, nlbc = 3, nrbc = 3
      real(kind(1e0)) :: xs(nv) = (/19.0e0,20.0e0,21.0e0/)
! Value of the interest rate and continuous dividend
      real(kind(1e0)) :: r = 0.05e0, dividend = 0.0e0
! Values of the min and max underlying values modeled
      real(kind(1e0)) :: x_min = 0.0e0, x_max = 60.0e0
 
! Define parameters for the integration step.
      integer, parameter :: nx = 121, nint = nx-1, n = 3*nx
      real(kind(1e0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt),&
                         dx, f(nv,nt)
      type(s_fcn_data) fcn_data
      integer :: nout
      real(kind(1e0)), external :: fkcoef, fkinitcond               
      external fkbc      

      integer :: i,i1,i2,i3,j
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
      allocate(fcn_data % rdata (6))
! Define an equally-spaced grid of points for the underlying price
      dx = (x_max-x_min)/real(nint)
      xgrid(1) = x_min
      xgrid(nx) = x_max
      do i = 2,nx-1
        xgrid(i) = xgrid(i-1) + dx
      end do

      call umach(2, nout)
      write(nout,'(T05,A)') "Constant Elasticity of Variance Model "//&
                            "for Vanilla Call"
      write(nout,'(T10,"Interest Rate ", F7.3, T38,"Continuous '//&
                   'Dividend ", F7.3  )') r, dividend
      write(nout,'(T10,"Minimum and Maximum Prices of Underlying ",'//&
                   '2F7.2)') x_min, x_max
      write(nout,'(T10,"Number of equally spaced spline knots ",I4,'//&
                   '/T10,"Number of unknowns ",I4)')&
                   nx-1,n
      write(nout,'(/T10,"Time in Years Prior to Expiration ",2X,'//&
                   '3F7.4)') time
      write(nout,'(T10,"Option valued at Underlying Prices  ",'//&
                   '3F7.2)') xs

      do i1 = 1,3       ! Loop over power
        do i2=1,3     ! Loop over volatility
          do i3=1,3  ! Loop over strike price
! Pass data through into evaluation routines.       
            fcn_data % rdata =& 
               (/ks(i3),x_max,sigma(i2),alpha(i1),r,dividend/)
            call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef,&
                              fkinitcond, fkbc, y, yprime,&
                              FCN_DATA = fcn_data)
! Evaluate solution at vector of points XS(:), at each time value
! prior to expiration.                        
            do i=1,nt
              f(:,i) = hqsval (xs, xgrid, y(:,i))
            end do 
            write(nout,'(/T05,"Strike=",F5.2," Sigma=", F5.2,'//&
              '" Alpha=", F5.2,/(T25," Call Option Values ",'//&
              'X,3F7.4))') ks(I3),sigma(I2),&
              alpha(i1),(f(i,:),i=1,nv)
          end do !i3 - Strike price loop
        end do !i2 - Sigma loop
      end do !i1 - Alpha loop
      end

! These functions and routines define the coefficients, payoff
! and boundary conditions.
      function fkcoef (X, TX, iflag, fcn_data)
        use mp_types
        implicit none
        real(kind(1e0)), intent(in) :: X, TX
        integer, intent(inout) :: iflag
         type(s_fcn_data), optional :: fcn_data
         real(kind(1e0)) :: fkcoef
         real(kind(1e0)) :: sigma, interest_rate, alpha, dividend,&
                            zero = 0.0e0, half = 0.5e0
                   
        sigma = fcn_data % rdata(3)
        alpha = fcn_data % rdata(4)
        interest_rate = fcn_data % rdata(5)
        dividend = fcn_data % rdata(6)
        select case (iflag)
          case (1)
! The coefficient derivative d(sigma)/dx
            fkcoef = half*alpha*sigma*x**(alpha*half-1.0e0)
! The coefficient sigma(x)				
          case (2)
            fkcoef = sigma*x**(alpha*half) 
          case (3)
! The coefficient mu(x)
            fkcoef = (interest_rate - dividend) * x
          case (4)
! The coefficient kappa(x)
            fkcoef = interest_rate
          end select
! Note that there is no time dependence
        iflag = 0
        return 
      end function fkcoef

      function fkinitcond(x, fcn_data)
        use mp_types
        implicit none
        real(kind(1e0)), intent(in) :: x
        type (s_fcn_data), optional :: fcn_data
        real(kind(1e0)) :: fkinitcond
        real(kind(1e0)) :: zero = 0.0e0
        real(kind(1e0)) :: strike_price

        strike_price = fcn_data % rdata(1)
! The payoff function
        fkinitcond = max(x - strike_price, zero)
        return
      end function fkinitcond

      subroutine fkbc (tx, iflag, bccoefs, fcn_data)
        use mp_types
        implicit none
        real(kind(1e0)), intent(in) :: tx
        integer, intent(inout) :: iflag
        real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
        type (s_fcn_data), optional :: fcn_data
        real(kind(1e0)) :: x_max, df, interest_rate, strike_price
        
        strike_price = fcn_data % rdata(1)
        x_max = fcn_data % rdata(2)
        interest_rate = fcn_data % rdata(5)
        select case (iflag)
          case (1)
            bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, 0.0e0/)
            bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 0.0e0/)
            bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
! Note no time dependence at left end  
            iflag = 0
          case (2)
            df = exp(interest_rate * tx)
            bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0,&
                             x_max - df*strike_price/)   
            bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 1.0e0/)
            bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
        end select
      end subroutine fkbc
