SUBROUTINE xskdconv (ear, ne, param, ifl, photar, photer)

      implicit none

      INTEGER ne, ifl
      REAL ear(0:ne), param(*), photar(ne), photer(*)

      CHARACTER*255 filenm

      integer MAXAR
      parameter (MAXAR = 1000)

      real eparam(9)
      real ephotar(MAXAR),p(MAXAR),ephoter(MAXAR),e(0:MAXAR),pcon(MAXAR)
      real emin, emax, elo, ehi, z, de
      real e1, e2
      logical firstcall,ccheck
      INTEGER i, nex, ix, jx
      
      save firstcall
      data firstcall/.true./


c This model convolves the input spectrum with
c the line profile from a Keplerian Accretion disk around
c both Schwarzschild & Kerr (arbitrary spin) black holes.
c
c It requires the use of an extensive set of transfer functions, which must be
c downloaded from the XSPEC website. These transfer functions should be placed
c in $LHEASOFT/../spectral/xspec/manager/XSKDLINE/
c
c It is not called directly from XSPEC, rather, it is called by the files
c xskdcN. Here, N is either 0 or 1 and denotes different modes of operation for
c the model. The convolution is performed by performing multiple calls to
c xskdline.f, which generates the convolution kernel, which is then used to
c determine the blurring of each energy bin. For further information how each 
c of the modes works, consult xskdcN.f, xskdline.f as appropriate.
c
c Finally, the user should note that chatter>=12 will case the model to dump
c the current line parameters to the screen

c=============

      if (firstcall) then
         write(*,*) 'this is a convolution model by Chris Done'
         write(*,*) 'USE THE EXTEND COMMAND'
         firstcall=.false.
      end if

c Extend energy array
      nex=MAXAR
      emin = log10(0.5 * ear(1))
      emax = log10(2.0 * ear(ne))
      de = (emax - emin) / real(nex)
      do ix = 0, nex
         e(ix) = 10.0**(emin + ix * de)
      end do
      
c     rebin old spectrum on new energy array
      i = 1
      do ix = 1, nex
         p(ix) = 0.0
         pcon(ix)=0.0
         emin = e(ix - 1)
         emax = e(ix)

c        only do if overlap on grid - otherwise keep p(ix)=0
         if (emin.gt.ear(0).and.emin.lt.ear(ne-1)) then
            do while(ear(i) < emin)
               i = i + 1
            end do

            do while(ear(i - 1) < emax)
               e1 = max(ear(i - 1), emin)
               e2 = min(ear(i), emax)
               p(ix) = p(ix) + photar(i) * (e2 - e1)/(ear(i)-ear(i-1))
               i = i + 1
            end do
            i = i - 1
         end if
      end do


c     set line parameters
      eparam(2)=param(1)
      eparam(3)=param(2)
      eparam(4)=param(3)
      eparam(5)=param(4)
      eparam(6)=param(5)
      eparam(7)=param(6)
      eparam(8)=param(7)
      eparam(9)=param(8)

c  convolution with kdline
      do ix=1,nex,1
         ccheck=.FALSE.
         If (ix.EQ.1) ccheck=.TRUE.
         eparam(1)=0.5*(e(ix)+e(ix-1))
         call xskdline (ccheck, e, nex, eparam, ifl,ephotar,ephoter)
         do jx=1,nex,1
            pcon(jx)=pcon(jx)+p(ix)*ephotar(jx)
         end do
         param(1)=eparam(2)
         param(2)=eparam(3)
         param(3)=eparam(4)
         param(4)=eparam(5)
         param(5)=eparam(6)
         param(6)=eparam(7)
         param(7)=eparam(8)
      end do

c rebinning
      ix = 1
      do i = 1, ne
         photar(i) = 0.0
         emin = ear(i - 1)
         emax = ear(i)
         do while(e(ix) < emin)
           ix = ix + 1
         end do

         do while(e(ix - 1) < emax)
           e1 = max(e(ix - 1), emin)
           e2 = min(e(ix), emax)
           photar(i) = photar(i) + pcon(ix) * (e2 - e1)/(e(ix)-e(ix-1))
           ix = ix + 1
         end do
         ix = ix - 1

      end do

      RETURN
      END