      subroutine linopac(lpri,lun11,optpp,rcem1,rcem2,elin,vturbi,t,
     $                   aatmp,delea,epi,ncn2,opakc,opakscatt,rccemis,
     $                   fline,lfast)
c
      implicit none
c
c     this routine puts line opacity into continuum bins
c     author:  T. Kallman
c 
      include './PARAM'
c
      real epi(ncn), opakc(ncn),optmp(ncn), opakscatt(ncn),
     $     rccemis(2,ncn),fline(ncn)
      integer ldon(2)
c
      real vturbi,bbb,elin,aatmp,vth,t,vturb,e0,
     $     deleturb,deleth,dele,aasmall,
     $     optpp,delet,profile,voigte,tst,ergsev,
     $     delea,dpcrit,expo, prftmp,sum,rcem1,rcem2
      integer ml1,nbinc,ncn2,ldir,mlc,mlm,lpri,lun11,
     $     ij,lfast,mlmin,mlmax
c
      data dpcrit/1.e-1/,ergsev/1.602197e-12/
c
c      if (lfast.ge.4) return
c
c     for scattering model, add in line opacity
      bbb=vturbi
      elin=abs(elin)
c     thermal width quantities
      vth=(1.3E+1)*sqrt(t/aatmp)
      vturb=bbb
      e0=(12398.54)/amax1(elin,1.E-24)
      deleturb=e0*(vturb/3.E+5)
      deleth=e0*(vth/3.E+5)
      dele=deleth+deleturb
      aasmall=delea/(1.E-24+dele)/12.56
      ml1=nbinc(e0,epi,ncn2)
      ml1=max(min(ncn-1,ml1),2)
      e0=epi(ml1)
      elin=12398.54/e0
      if (lpri.gt.1) write (lun11,*) 
     &   'optpp,dpcrit*opakc(ml1),ml1,deleth,delea:',
     &    optpp, dpcrit*opakc(ml1),ml1,deleth,delea
c     here is what we do to get the heating right
      if (lfast.gt.2) then
          prftmp=dele*2./(epi(ml1+1)-epi(ml1-1))
          tst=optpp*prftmp
c          opakc(ml1)=opakc(ml1)+tst
          fline(ml1)=fline(ml1)+(rcem1+rcem2)*prftmp/dele/ergsev
          opakscatt(ml1)=opakscatt(ml1)+tst
          if (lpri.gt.1) write (lun11,*)optpp,prftmp,tst,dele
        else
          if (optpp.gt.dpcrit*opakc(ml1)) then
c           apply binned opacity
            delet=(epi(ml1)-e0)/dele
            if (aasmall.gt.1.E-6) then
                profile=voigte(delet,aasmall)
              else
                profile=expo(-delet*delet)/1.772
              endif
            optmp(ml1)=profile
            if (lpri.gt.1) write (lun11,*) ml1, epi(ml1), epi(ml1+1), 
     &                            opakc(ml1), tst
            mlc=0
            ldir=1
            ldon(1)=0
            ldon(2)=0
            mlmin=ml1-1
            mlmax=ml1
            do while ((ldon(1)*ldon(2).eq.0).and.(mlc.lt.10000))
              mlc=mlc+1
              do ij=1, 2
                ldir=-ldir
                if (ldon(ij).ne.1) then
                  mlm=ml1+ldir*mlc
                  if ((mlm.le.ncn2).and.(mlm.ge.1)) then
                    mlmin=min(mlm,mlmin)
                    mlmax=max(mlm,mlmax)
                    delet=(epi(mlm)-e0)/dele
                    if (aasmall.gt.1.E-6) then
                         profile=voigte(delet,aasmall)
                       else
                         profile=expo(-delet*delet)/1.772
                       endif
                    optmp(mlm)=profile
                    tst=optpp*profile
                    if (lpri.gt.1) write (lun11,*) mlm, epi(mlm), ij, 
     &                delet, sqrt(aasmall), profile, tst, opakc(mlm)
                    endif
                  if ((tst.lt.dpcrit*opakc(mlm)).or.(mlm.le.1).or.
     &              (mlm.ge.ncn).or.(mlc.gt.10000).or.
     &              (abs(delet).gt.2000.)) ldon(ij)=1
                  endif
                enddo
              enddo
            sum=0.
            mlmin=max(mlmin,2)
            do mlm=mlmin,mlmax+1
              sum=sum+(optmp(mlm)+optmp(mlm-1))*(epi(mlm)-epi(mlm-1))/2.
              enddo
            do mlm=mlmin,mlmax
              rccemis(1,mlm)=rccemis(1,mlm)+optmp(mlm)*rcem1/sum/ergsev
              rccemis(2,mlm)=rccemis(2,mlm)+optmp(mlm)*rcem2/sum/ergsev
              opakc(mlm)=opakc(mlm)+optmp(mlm)*optpp*dele/sum
              opakscatt(mlm)=opakscatt(mlm)+optmp(mlm)*optpp*dele/sum
              enddo            
            endif
          endif
c
      return
      end
