c
c###################################################################
c
      subroutine verner(pt, pcf, kncf, pfit, kermsg)
c
c     Verner, Yakovlev, Band, Trzhaskovskaya, 
c     At. Data. Nucl. Data Tables 55 2233 (1993)
c
c     this is a subroutine to calculate photoionization cross-sections
c     in cm[2] as function of photon energy in ev for an electron
c
c     sigma_{ph}(k) = S_0*F(y), y =k/k_0,                                 (1)
c     F(y) = [(y-1)^2+(y_w)^2][y^{0.5P-xl-5.5}][(1+sqrt(y/y_a))^{-P}]     (2)
c     
c     pt = k         photon energy (eV)
c
c     pcf is the coefficient data array, where
c     pcf(1) =  E_{th}       Ionization threshold energy in eV
c     pcf(2) =  k_{max}      Maximum photon energy in eV to 
                              which the fitting was performed                                           
c     pcf(3) =  k_0          Fit parameter used in Eq.(1) in eV
c     pcf(4) =  S_0          Fit parameter used in Eq.(1) in Mb
c     pcf(5) =  P,           Dimensionless fit parameters used in Eq.(2)
c     pcf(6) =  y_a,         Dimensionless fit parameters used in Eq.(2)
c     pcf(7) =  y_w          Dimensionless fit parameters used in Eq.(2)
c     pcf(8) =  del          Relative root-mean-square error of fitting in percent
c     pcf(9) =  xl           Angular momentum quantum number

c     kermsg = blank if no errors
c
c     pfit = photo ionization cross-sections in cm[2]
c
c     written by H. K. Chung , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
      double precision pt, pcf, pfit, q
      dimension pcf(9)
      character*(*) kermsg
      data ry/13.58/
c
      kermsg = ' '

      eth=pcf(1)
      xmax=pcf(2)
      x0=pcf(3)
      S0=pcf(4)
      p=pcf(5)
      ya=pcf(6)
      yw=pcf(7)
      del=pcf(8)
      xl=pcf(9)
c

      if(pt.le.eth)then
c         goto 900
         pt=eth
      endif

      y=pt/x0
      q=-(5.5+xl-0.5*p)
      a1=y**q
      a2=(y-1.)**2 + yw**2
      a3=(1.+sqrt(y/ya))**(-p)
      fy=a1*a2*a3
c
      pfit = S0 * 1.e-18 * fy

      return

c 900  pfit=1.e-30

      return
c
      end