      subroutine phint53new(stmpp,etmpp,ntmp,ethi,pirt,rrrt,piht,rrcl,
     $ abund1,abund2,ptmp1,ptmp2,xpx,opakab,
     $ opakc,rccemis,lpri,epi,ncn2,bremsa,t,swrat,xnx,crit,
     $ lfast,lun11)
c
c
c     this routine does the integration over the spectrum as required by
c     photo.
c     uses power law piecewise analytic integrals, assuming \sigma~e^{-3}
c        for recombination integrals.
c     special version for photemis
c
      implicit none
c
      integer ntmp
c
      include './PARAM'
c
      real  bremsa(ncn),epi(ncn),stmpp(ntmp),etmpp(ntmp)
      real rccemis(2,ncn),opakc(ncn)
      real sgbar(ncn)
      integer lpri,ncn2,lfast,lun11
      real ethi,pirt,rrrt,piht,rrcl,abund1,abund2,ptmp1,ptmp2,xpx,
     $     opakab,t,swrat,xnx,crit
      real eth,ergsev,bk,tm,bktm,ener,sgtmp,epii,sgtp,optmp,
     $     sumr,sumh,sumi,sumc,tempi,
     $     bremtmp,tempr,sumho,exptst,
     $     atmp2,rnist,tsq,ethsht,
     $     wwir,wwih,bbnurjp,e2t,e2,e1,e1o,bremtmpp,
     $     epiip,e2o,tempc,tempcp,exptmpp,rctmp2,rctmp1,
     $     s2,s2to,s2o,s2t,sgtpp,sum,tempip,temphp,temprp,temph,
     $     wwirp,wwicp,wwic,wwii,wwihp,wwiip,enermx
      integer lprisv,numcon2,nphint,nb1,kl,jk,nbn,klmax,nbinc
c
      logical done
c
      data ergsev/1.602197e-12/
      data bk/1.38062e-16/
c
c
      eth=ethi
      lprisv=lpri
c
      tm=t*1.e4
      bktm=bk*tm/ergsev
      tsq = sqrt(t)
      rnist=(2.61e-21)*swrat/t/tsq
c
      if (lpri.ge.1)
     $  write (lun11,*)'in phint53:',
     $      eth,xnx,swrat,t,
     $      stmpp(1),etmpp(1),ntmp,lfast,ptmp1,ptmp2
     $      ,abund1,abund2,rnist,crit
      ethsht=eth/bktm
      ethsht=amax1(ethsht,0.)
c
      numcon2=max(2,ncn2/50)
      nphint=ncn2-numcon2
c
      sumr = 0.
      sumh = 0.
      sumho = 0.
      sumc = 0.
      sumi=0.
      ener=ethi+etmpp(1)*(13.598)
      nb1=nbinc(ener,epi,ncn2)
      do while ((epi(nb1).lt.ener).and.(nb1.lt.nphint))
        nb1=nb1+1
        enddo
      nb1=nb1-1
      if (lpri.ge.1) write (lun11,*)'nb1=',nb1,ener,nphint
      if (nb1.ge.nphint) return
c      if (epi(kl+1).lt.ethi) return
      tempr=0.

c     step through cross section, map onto kl grid
      jk=1
      enermx=ethi+etmpp(ntmp)*(13.598)
      nbn=nbinc(enermx,epi,ncn2)
      nbn=max(nbn,min(nb1+1,ncn2-1))
      ener=ethi+etmpp(jk)*(13.598)
      sgtmp=stmpp(1)
      sgbar(max(1,nb1-1))=0.
      kl=nb1
      jk=1
      e1=epi(kl)
      e2=ethi+etmpp(jk)*(13.598)
      s2=stmpp(jk)
      if (e1.lt.e2) then
        kl=kl+1
        e1=epi(kl)
        endif
      e1o=e2
      sum=0.
      done=.false.
      do while (.not.done)
c       step through where jk grid is finer
        if (lpri.ne.0) write (lun11,*)'mapping loop:',jk,kl,e2,e1
        do while ((e2.lt.e1).and.(jk.lt.(ntmp-1)))
          jk=jk+1
          e2o=e2
          s2o=s2
          e2=ethi+etmpp(jk)*(13.598)
          s2=stmpp(jk)
          sum=sum+(s2+s2o)*(e2-e2o)/2.
          if (lpri.ne.0) write (lun11,*)'jk loop',jk,e2,s2,e2o,s2o,sum
          enddo
c       kl bin exceeds jk bin, subtract off extra
        sum=sum-(s2+s2o)*(e2-e2o)/2.
c       now interpolate to find value at kl bin
        e2t=e1
        if (e2-e2o.gt.1.e-8) then
             s2t=s2o+(s2-s2o)*(e2t-e2o)/(e2-e2o+1.e-24)
           else
             s2t=s2o
           endif
        s2to=s2t
c       now update sum at kl bin
        sum=sum+(s2t+s2o)*(e2t-e2o)/2.
c       save
        if (abs(e1-e1o).gt.1.e-34) then
            sgbar(kl)=sum/(e1-e1o)
          else
            sgbar(kl)=0.
          endif
        if (lpri.ne.0) write (lun11,*)'saving:',kl,e1,sgbar(kl),sum,s2t
        e1o=e1
c       increment kl
        kl=kl+1
        e1=epi(kl)
c       step through where kl grid is finer
        do while (e1.lt.e2)
          e2t=e1
          if (e2-e2o.gt.1.e-8) then
              s2t=s2o+(s2-s2o)*(e2t-e2o)/(e2-e2o)
            else
              s2t=s2o
            endif
          s2to=s2t
          sum=(s2t+s2to)*(e1-e1o)/2.
          sgbar(kl)=sum/(e1-e1o)
          if (lpri.ne.0) write (lun11,*)'kl loop',kl,e1,
     $                                   sgbar(kl),sum,s2t
          e1o=e1
          kl=kl+1
          e1=epi(kl)
          enddo
c       update sum for remaining bit
        sum=(s2+s2t)*(e2-e2t)/2.
        if (lpri.ne.0) write (lun11,*)'testing for done:',kl,nphint,
     $                                                    jk,ntmp
        if ((kl.gt.nphint-1).or.(jk.ge.ntmp-1))
     $       done=.true.
        enddo
      klmax=kl-1
c
c
c     preliminary setup
      sgtpp=sgbar(nb1)
      bremtmpp=bremsa(nb1)/(12.56)
      epiip=epi(nb1)
      temprp=(12.56)*sgtpp*bremtmpp/epiip
      temphp=temprp*epiip
      exptst=(epiip-eth)/bktm
      exptmpp=exp(-exptst)
      bbnurjp=epiip*epiip*epiip*(1.571e+22)*2.
      tempip=rnist*(bremtmpp+bbnurjp*(ptmp1+ptmp2))
     $  *sgtpp*exptmpp/epiip
      tempcp=tempip*epiip
c
      kl=nb1
      epii=epi(kl)
      if (lpri.ne.0) write (lun11,*)'kl=',kl,klmax,sumh,sumho
      rctmp1=0.
      rctmp2=0.
      do while (kl.lt.klmax)
c
c       the basics
        sgtmp=max(0.,sgbar(kl))
        sgtp=sgtmp
        sgtpp=sgbar(kl+1)
        bremtmp=bremsa(kl)/(12.56)
        bremtmpp=bremsa(kl+1)/(12.56)
        epii=epi(kl)
        epiip=epi(kl+1)
c
c       pi rate
        tempr=temprp
        temprp=(12.56)*sgtpp*bremtmpp/epiip
        wwir=(epiip-epii)/2.
        wwirp=wwir
        sumr = sumr + (tempr*wwir+temprp*wwirp)
c
c       heat
        temph=temphp
        temphp=temprp*epiip
        wwih=wwir
        wwihp=wwih
        sumho=sumh
        sumh = sumh + (temph*wwih+temphp*wwihp)
c
c       rec
        exptst=(epiip-eth)/bktm
        if (exptst.lt.20.) then
          bremtmpp=bremsa(kl+1)/(25.3)
          exptmpp=exp(-exptst)
          bbnurjp=epiip*epiip*epiip*(1.571e+22)*2.
          tempi=tempip
          tempip=rnist*(bremtmpp+bbnurjp*(ptmp1+ptmp2))
     $        *sgtpp*exptmpp/epiip
          wwii=wwir
          wwiip=wwir
          sumi = sumi + (tempi*wwii+tempip*wwiip)
c
c         cool
          tempc=tempcp
          tempcp=tempip*epiip
          wwic=wwir
          wwicp=wwir
          sumc = sumc+tempc*wwic+tempcp*wwicp
c
          atmp2=tempc
          rctmp1=abund2*atmp2*ptmp1*xpx*xnx
          rccemis(1,kl)=rccemis(1,kl)+rctmp1
          rctmp2=abund2*atmp2*ptmp2*xpx*xnx
          rccemis(2,kl)=rccemis(2,kl)+rctmp2
c
          endif
c
c       emiss and opac
c       the emission must be fudged to get the right cooling with a
c         trapezoid integration.
        optmp=abund1*xpx*sgtp
        if (kl.le.(nb1+1)) opakab=optmp
        opakc(kl)=opakc(kl)+optmp
c
c       print
        if (lpri.ge.1)
     $      write (lun11,901)kl,epii,sgtp,bremtmp,
     $         tempr,temprp,wwir,wwirp,sumr,
     $         temph,temphp,wwih,wwihp,sumh,
     $         tempi,tempip,wwii,wwiip,sumi,
     $         tempc,tempcp,wwic,wwicp,sumc
     $         ,rctmp1,rctmp2,exptst
 901       format(1x,'found something',i6,25(1pe11.3))
        if (sgtp.lt.0.) stop 'phint error'
c
        kl=kl+1
        enddo
c
      pirt = pirt + sumr
      rrrt = rrrt + xnx*sumi
      piht = piht + sumh*ergsev
      rrcl = rrcl + xnx*sumc*ergsev

      if (lpri.ge.1) write (lun11,*)'in phint53:',eth,pirt,rrrt
     $         ,piht,rrcl
      lpri=lprisv
c
      do kl=nb1,nbn
        sgbar(kl)=0.
        enddo
c
c
      return
      end
