subroutine bkhsgo(freq,et,d,b,na,a,ss,nmax,iz,nsh,sg)
c     =====================================================
c
c
c subroutine to calculate K and L shell photoionization cross-sections
c       -based on Tim Kallman's bkhsgo subroutine from XSTAR, modified
c        by Omer Blaes 5-7-98.
c       -na.ne.2 bug corrected on 2/24/00 by Omer Blaes
c
c freq is photon frequency in Hz (note that this subroutine immediately
c converts it into eV)
c
c et is threshold energy in eV
c
c iz is the ionization stage of the species being photoionized (1=neutral etc.)
c
c ss is the iz'th element of the array ss(nmax) in Tim's original version
c
c sg is returned as the contribution to the photoionization cross-section
c in cm^2 due to whatever process is being considered.
c
c     this routine does the work in computing cross sections by the
c     method of barfield, et. al.
c
c
c
      INCLUDE 'IMPLIC.FOR'
      dimension b(5),a(11,5)
c
      data sigth/1.e-34/
c
      tmp1 = 0.
      jj = 1
      epii = 4.1357e-15*freq
      sg=0.
      if ( epii.gt.et ) then
        xx = epii*(1.e-3) - d
        if ( xx.le.0. ) return
        do nna=1,na
          if ( xx.ge.b(nna) ) jj = jj + 1
        enddo
        if ( jj.le.na ) then
        if(xx.lt.0.) xx=0.
          yy = log10(xx)
          tmp = 0.
          do kkk = 1,11
            kk = 12 - kkk
            tmp = a(kk,jj) + yy*tmp
          end do
          if(tmp.lt.-50.) tmp=-50.
          if(tmp.gt.24.) tmp=24.
          sgtmp = 10.**(tmp-24.)
        nelec=nmax+1-iz
        if(nelec.gt.nsh) nelec=nsh
          enelec = float(nelec)
          tmp1o = tmp1
          tmp1=sgtmp*ss
        if(tmp1.lt.sigth*enelec) tmp1=sigth*enelec
          if ( epii.ge.5.e+4 ) then
           if(tmp1.gt.tmp1o) tmp1 = tmp1o
          end if
          sg = sg + tmp1
        endif
      endif
      return
      end