      subroutine heatt(jkk,lpri,lun11,lfpi,
     $       t,r,cfrac,delr,xee,xpx,abel,
     $       epi,ncn2,bremsa,
     $       idat1,rdat1,kdat1,nptrs,
     $       npar,npnxt,npfi,npfirst,nplin,nlsvn,
     $       npconi2,ncsvn,rlev,ilev,nlev,klev,
     $       zrems,zremso,elumab,elumabo,elum,elumo,
     $       rcem,rccemis,opakc,opakscatt,cemab,fline,flinel,
     $       brcems,cmp1,cmp2,httot,cltot,hmctot,
     $             cllines,clcont,htcomp,clcomp,clbrems)
c
c     this routine calculates heating and cooling.
c     author:  T. Kallman
c
      implicit none
c
      include './PARAM'
c
      integer nptmpdim
      parameter (nptmpdim=200000)
c
c     master data
c     global xstar data
c     master data
      integer idat1(nidat1)
      real rdat1(nrdat1)
      integer nptrs(nptt,ndat2)
      character*1 kdat1(nkdat1)
c     pointers to master data
      integer npar(ndat2),npnxt(ndat2),
     $      npfirst(ntyp)
      integer npfi(ntyp,nni)
c     pointers to line data
      integer nplin(nnnl)
c     pointers to line data
      integer npconi2(ndat2)
c     line emissivities
      real rcem(2,nnnl)
c     energy bins
      real epi(ncn)
c     continuum flux
      real bremsa(ncn)
c     continuum emissivities
      real rccemis(2,ncn),brcems(ncn)
c     continuum opacities
      real opakc(ncn),opakscatt(ncn)
c     state variables
      real r,t,xpx,delr,delrl
c     heating-cooling variables
c     input parameters
      real xee
      integer ncn2,lpri,lun11
      integer nlsvn,ncsvn
      real rdat(nptmpdim)
      real fline(ncn),flinel(ncn)
c     level populations
      real abel(nl)
      real cemab(2,nnml)
      character*1 kblnk
      real zrems(4,ncn),zremso(4,ncn)
      real elum(3,nnnl),elumo(3,nnnl)
      real elumab(2,nnml),elumabo(2,nnml)
      real xeltp,cllines,clcont,cmp1,cmp2,cltot,
     $     hmctot,htcomp,clcomp,clbrems,etst,ekt,epiio,
     $     epii,fac,fpr2,optpp,optp2,tmpc1,tautmp,cfrac,
     $     tmpc2,tmpho,xnx,httot,tmpc,tmph,tmpco,hmctmp
      real elin,ener,eth,r19,ergsev,tmpscat
      character*1 kdat(nptmpdim),kdati(nptmpdim)
      character*1 klev(100,nd)
      real tmp2,tmp2o
      integer lskp
      real rlev(10,nd)
      integer ilev(10,nd)
      integer idat(nptmpdim)
      integer idest1,lfpi,jk,klel,kl,
     $     klion,mlleltp,mllel,mlel,mlion,mt2,nilin,nnz,numcon,
     $     nblin
      integer nlev,nlevmx,mltype,ml,mllz,jkk,ltyp,
     $     lrtyp,lcon,nrdt,nidt,nkdt,lk,kkkl,
     $     lprisv,mm,nbinc,mlpar,mlm
c
c
      data kblnk/' '/
      data ergsev/1.602197e-12/
c
c
      xnx=xpx*xee
      lprisv=lpri
      if (lpri.ge.1) lpri=2
      if (lpri.gt.1) write (lun11,*)'in heatt',httot,cltot,delr,r
      if (lpri.gt.1) write (lun11,*)ncsvn
      numcon=ncn2
      r19=r*(1.e-19)
      fpr2=12.56*r19*r19
c
c     comment these out to implement scattering
      clbrems=0.
      lskp=1
      tmp2=0.
      do kl=1,numcon
        tmp2o=tmp2
        tmp2 =brcems(kl)
        if ( kl.ge.2 ) clbrems=clbrems+(tmp2+tmp2o)
     &                *(epi(kl)-epi(kl-lskp))*ergsev/2.
        enddo
c
      delrl=max(delr,1.e-24*r)
      if (lfpi.le.2) then
        httot=0.
        cltot=0.
        epii=epi(1)
        hmctmp=0.
        tmpc=0.
        tmph=0.     
        do kl=1,numcon
          optpp=opakc(kl)
          optp2=max(1.e-34,optpp)
c         for outward only
          epiio=epii
          epii=epi(kl)
          tmpho=tmph
          tautmp=optp2*delrl
          fac=delrl
          if (tautmp.gt.0.01)
     $      fac=(1.-exp(-tautmp))/(1.e-24+optp2)
          tmph=bremsa(kl)*optp2
          tmpco=tmpc
c         these are all statements needed to implement scattering
c          tmph=bremsa(kl)*optp2*fac*fpr2
c          tmpc1=rccemis(1,kl)+brcems(kl)*(1.-cfrac)/2.
c          tmpc2=rccemis(2,kl)+brcems(kl)*(1.+cfrac)/2.
c          tmpc=(tmpc1+tmpc2+fline(kl))*fac*fpr2
c          hmctmp=-(-tmph+tmpc)
c          zrems(1,kl)=zremso(1,kl)-hmctmp
          tmpc1=rccemis(1,kl)
          tmpc2=rccemis(2,kl)
          tmpc=tmpc1+tmpc2
          zrems(1,kl)=max(0.,zremso(1,kl)
     $                +(-tmph+tmpc+fline(kl))*fac*fpr2)
c         this is slightly dangerous
c          flinel(kl)=flinel(kl)*fac
          zrems(2,kl)=zremso(2,kl)+tmpc1*fac*fpr2
          zrems(3,kl)=zremso(3,kl)+tmpc2*fac*fpr2
          tmpscat=opakscatt(kl)*bremsa(kl)
          zrems(4,kl)=zremso(4,kl)+tmpscat*fac*fpr2
          if (kl.gt.1) then
            httot=httot+(tmph+tmpho)
     $       *(epii-epiio)*ergsev/2.
            cltot=cltot+(tmpc+tmpco)
     $       *(epii-epiio)*ergsev/2.
            endif
          if (lpri.ge.1) write (lun11,9009)kl,epii,optpp,
     $     bremsa(kl),tmph,tmpc,fline(kl),flinel(kl),httot,cltot
     $       ,rccemis(1,kl)+rccemis(2,kl),hmctot
 9009     format (1x,i6,14(1pe12.4))
          enddo
        if (lpri.ge.1)
     $   write (lun11,*)'continuum heating, cooling:',
     $       httot,cltot
        clcont=cltot
        do jkk=1,nlsvn
          jk=jkk
          ml=nplin(jk)
          mlm=ml-1
          call dread(ltyp,lrtyp,lcon,
     $      nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $      idat1,rdat1,kdat1,nptrs,0,lun11)
          elin=abs(rdat(1))
          if ((elin.lt.1.e+8).and.(elin.gt.1.)
     $      .and.(ml.ne.0)) then
            nilin=npar(ml)
            ergsev=1.602197e-12
            ener=ergsev*(12398.54)/amax1(elin,1.e-24)
            etst=ener/ergsev
            tmpc1=rcem(1,jk)
            tmpc2=rcem(2,jk)
            tmpc=tmpc1+tmpc2
c           comment this out to implement scattering
            cltot=cltot+tmpc
            nblin=nbinc(ener,epi,ncn2)
            optp2=opakc(nblin)
            tmph=elumo(1,jk)*optp2/fpr2
            fac=delrl
            elum(1,jk)=max(0.,elumo(1,jk)+(-tmph+tmpc1)*fac*fpr2)
            elum(2,jk)=max(0.,elumo(2,jk)+(-tmph+tmpc2)*fac*fpr2)
            if (lpri.ge.1)
     $       write (lun11,*)jk,nilin,etst,
     $        rcem(1,jk),rcem(2,jk),delrl,fpr2,cltot,
     $        elumo(2,jk),elum(2,jk)
            endif
          enddo
c
c       calculate rrc luminosities
C       First look for element data (jk is element index)
        if (lpri.ge.3) write (lun11,*)'rrc print:'
        klel=11
        mlel=npfirst(klel)
        jk=0
        do while (mlel.ne.0)
          jk=jk+1
          mt2=mlel-1
          call dread(ltyp,lrtyp,lcon,
     $      nrdt,rdat,nidt,idat,nkdt,kdat,mt2,
     $      idat1,rdat1,kdat1,nptrs,0,lun11)
          mllel=idat(nidt)
          xeltp=rdat(1)
          xeltp=abel(mllel)
          nnz=idat(1)
          if (lpri.ge.3)
     $      write (lun11,*)'element:',jk,mlel,mllel,nnz,
     $                  (kdat(mm),mm=1,nkdt)
C         ignore if the abundance is small
          if (xeltp.lt.1.e-10) then
              jkk=jkk+nnz
            else
c             now step thru ions (jkk is ion index)
              klion=12
              mlion=npfirst(klion)
              jkk=0
              kl=0
              do while ((mlion.ne.0).and.(kl.lt.nnz))
                jkk=jkk+1
C               retrieve ion name from kdati
                mlm=mlion-1
                call dread(ltyp,lrtyp,lcon,
     $            nrdt,rdat,nidt,idat,nkdt,kdati,mlm,
     $            idat1,rdat1,kdat1,nptrs,0,lun11)
C               if not accessing the same element, skip to the next element
                mlleltp=idat(nidt-1)
                if (mlleltp.eq.mllel) then
                  kl=kl+1
                  if (lpri.ge.3)
     $              write (lun11,*)'  ion:',kl,jkk,mlion,mlleltp,
     $                          (kdati(mm),mm=1,nkdt)
c                 now find level data
c                 step thru types
                  nlevmx=0
                  mltype=13
                  ml=npfi(mltype,jkk)
                  mllz=0
                  if (ml.ne.0) mllz=npar(ml)
c                 step thru records of this type
                  mlpar=0
                  if (ml.ne.0) mlpar=npar(ml)
                  do while ((ml.ne.0).and.(mlpar.eq.mllz))
                    mlm=ml-1
                    call dread(ltyp,lrtyp,lcon,
     $                nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $                idat1,rdat1,kdat1,nptrs,0,lun11)
                    nlev=idat(nidt-1)
                    nlevmx=max(nlevmx,nlev)
                    if ((nlev.gt.0).and.(nlev.le.nd)) then
                      do  lk=1,nrdt
                        rlev(lk,nlev)=rdat(lk)
                        enddo
                      do lk=1,nidt
                        ilev(lk,nlev)=idat(lk)
                        enddo
                      do lk=1,nkdt
                        klev(lk,nlev)=kdat(lk)
                        enddo
                      do lk=nkdt+1,20
                        klev(lk,nlev)=kblnk
                        enddo
                      endif
                    ml=npnxt(ml)
                    if (ml.ne.0) mlpar=npar(ml)
                    enddo
                  nlev=nlevmx
                  mltype=7
                  ml=npfi(mltype,jkk)
                  mllz=0
                  if (ml.ne.0) mllz=npar(ml)
                  mlpar=0
                  if (ml.ne.0) mlpar=npar(ml)
                  do while ((ml.ne.0).and.(mlpar.eq.mllz))
c                   step thru records of this type
                    mlm=ml-1
                    call dread(ltyp,lrtyp,lcon,
     $                nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $                idat1,rdat1,kdat1,nptrs,0,lun11)
                    kkkl=npconi2(ml)
                    idest1=idat(nidt-1)
                    if ((kkkl.gt.0).and.(kkkl.le.ndat2)
     $                .and.(cemab(1,kkkl).gt.1.e-36)) then
                      eth=rlev(4,idest1)-rlev(1,idest1)
                      tmpc=(cemab(1,kkkl)+cemab(2,kkkl))
                      optp2=0.
                      fac=delrl
                      tmph=elumabo(1,kkkl)*optp2/fpr2
                      elumab(1,kkkl)=max(0.,elumabo(1,kkkl)
     $                            +(-tmph+tmpc)*fac*fpr2/2.)
                      elumab(2,kkkl)=max(0.,elumabo(2,kkkl)
     $                            +(-tmph+tmpc)*fac*fpr2/2.)
                      if (lpri.ge.3)
     $                  write (lun11,981)kkkl,eth,idest1,
     $                    cemab(1,kkkl),cemab(2,kkkl),
     $                    elumabo(1,kkkl),elumab(1,kkkl)
 981                    format (1x,i6,1pe11.3,i6,6(1pe11.3))
                      endif
                    ml=npnxt(ml)
                    if (ml.ne.0) mlpar=npar(ml)
                    enddo
                  endif
C               Go to next ion
                mlion=npnxt(mlion)
                enddo
            endif
          mlel=npnxt(mlel)
C         Go to next element
          enddo
       cllines=cltot-clcont
       if (lpri.ge.1)
     $   write (lun11,*)'line cooling',cltot,cllines
        endif
c
c
      fac=delrl
      ekt = t*(0.861707)
      htcomp = cmp1*xnx*ergsev
      clcomp = ekt*cmp2*xnx*ergsev
c      these statements needed to implement scattering
c      httot=htcomp*fac*fpr2+httot
c      cltot=clcomp*fac*fpr2+cltot
      httot=htcomp+httot
      cltot=clcomp+cltot+clbrems
      hmctot=(httot-cltot)
     $   /((abs(httot)+abs(cltot))/2.)
      if (lpri.ge.1) write (lun11,9953)htcomp,clcomp,cmp1,cmp2,
     $   clbrems,httot,cltot,hmctot
 9953 format (1h , ' compton heating, cooling=',8e12.4)
      lpri=lprisv
c

      return
      end
