      subroutine func2(jkk,kl,ilimh,lpriz,lun11,lfpi,vturbi,
     $       t,trad,r,xee,xpx,xh1,xh0,cfrac,
     $       epi,ncn2,bremsa,bremsint,
     $       tau0,tauc,opakc,
     $       idat1,rdat1,kdat1,nptrs,np2,
     $       npar,npnxt,npfi,npfirst,
     $       nplin,nplini,nlsvn,npcon,npconi,npilev,npilevi,
     $       npconi2,ncsvn,rates,vsav,idrates,
     $       rniss,rlev,ilev,nmat,
     $          nlpt,iltp,nlev,klev,ajis,cjis,
     $         rrrtot)
c
c     this routine calculates rates affecting level populations
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),nplini(ndat2)
c     pointers to line data
      integer npcon(nnml),npconi2(ndat2),npconi(ndat2)
      integer npilev(nd,nni),npilevi(nnml)
c     energy bins
      real epi(ncn)
c     continuum flux
      real bremsa(ncn),bremsint(ncn)
c     continuum emissivities
      real rccemis(2,ncn)
      real fline(ncn),flinel(ncn)
c     continuum opacities
      real opakc(ncn),opakscatt(ncn)
c     element abundances
c     the saved rates
      real rates(4,ndat2)
      integer idrates(2,ndat2)
      real vsav(4,ndat2)
c     state variables
      real r,t,xpx
c     heating-cooling variables
c     input parameters
      real trad
      real vturbi,xee
      integer ncn2,lpri,lun11,lfpi,np2
      integer nlsvn,ncsvn
      real rlev(10,nd)
      integer ilev(10,nd),nlpt(nd),iltp(nd)
      character*1 kdat(nptmpdim),klev(100,nd)
      character*48 kdesc2
      real rdat(nptmpdim)
      integer idat(nptmpdim)
      real tsq,ans1,ans2,xh1,xh0,cfrac,rrrtot
      real abund1,abund2,ptmp1,ptmp2,ans3,ans4,opakb1
      integer idest1,idest2,idest3,idest4
c     line optical depths
      real tau0(2,nnnl)
c     continuum flux
      real tauc(2,nnml)
c
      character*1 kblnk
      real rniss(nnml)
      real ajis(nds,nds),cjis(nds,nds)
      real rrrt3,tau1,tau2,airtmp,rrrtot2,e1,e2,pescl,pescv,ptmp
      integer nlev,lpriz,nlevmx,mltype,ml,mllz,mlrdesc,lpriu,
     $        lforce,llo,lup,ilimh,jkk,kl,nmat,ml1,ltyp,
     $        lrtyp,lcon,nrdt,nidt,nkdt,lk,ll,kkkl,jkkl,mlpar,mlm
c
      data kblnk/' '/
c
c
c      if (lpri.ne.0)
c     $  write (lun11,*)'in func2, inputs:',t,
c     $         xee,xpx,xnx
c      if (lpri.ne.0) write (lun11,*)'in func2:',ndo
c
c
      lpri=lpriz
      tsq=sqrt(t)
c     zero temporaries
      rrrt3=rrrtot
      rrrtot=0.
      rrrtot2=0.
      nlev=0
      do ml1=1,nmat
        rniss(ml1)=0.
        do ll=1,100
          klev(ll,ml1)=kblnk
          enddo
        enddo
      nmat=0
c     now find level data
c     step thru types
      nlevmx=0
      mltype=13
      ml=npfi(mltype,jkk)
      mllz=npar(ml)
      mlpar=npar(ml)
c     step thru records of this type
      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
          nlpt(nlev)=ml
          iltp(nlev)=ltyp
c         if (lpri.ne.0) write (lun11,9101)
c     $     ml,nlev,ltyp,lrtyp,(rdat(mm),mm=1,4),idat(1),idat(2),
c     $     idat(3),(kdat(mm),mm=1,8)
c 9101     format (1x,'level quantities:',4i6,4(1pe12.5),3i6,8a1)
          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
          endif
        ml=npnxt(ml)
        if (ml.ne.0) mlpar=npar(ml)
        enddo
      nlev=nlevmx
c
c     now find the rates affecting this ion
      mltype=0
      do while (mltype.lt.ntyp)
        mltype=mltype+1
        mlrdesc=mltype
        if (((mlrdesc.lt.10).or.(mlrdesc.gt.13))
     $    .and.(mlrdesc.gt.0)) then
          ml=npfi(mltype,jkk)
          if (ml.gt.0) then
            mllz=npar(ml)
c           step thru records of this type
            mlpar=npar(ml)
            do while ((ml.ne.0).and.(mlpar.eq.mllz))
              if (nptrs(3,ml).eq.mlrdesc) then
                mlm=ml-1
                call dread(ltyp,lrtyp,lcon,
     $            nrdt,rdat,nidt,idat,nkdt,kdat,mlm,
     $            idat1,rdat1,kdat1,nptrs,0,lun11)
c               calculate rates
                lpriu=lpri
                lforce=1
                abund1=0.
                abund2=0.
                ptmp2=1.
                ptmp1=0.
                if ((lrtyp.eq.7).or.((lrtyp.eq.1).and.(ltyp.ne.53))) 
     $            then
                  kkkl=npconi2(ml)
                  if ((kkkl.gt.0).and.(kkkl.le.ndat2)) then
                    tau1=tauc(1,kkkl)
                    tau2=tauc(2,kkkl)
                    ptmp1=pescv(tau1)*(1.-cfrac)
                    ptmp2=pescv(tau2)*(1.-cfrac)
     $                  +2.*pescv(tau1+tau2)*cfrac
                    ptmp=(ptmp1+ptmp2)
                    call ucalc(ltyp,lrtyp,ml,lcon,jkk,vturbi,
     $                nrdt,rdat,nidt,idat,nkdt,kdat,ans1,ans2,
     $                ans3,ans4,idest1,idest2,idest3,idest4,
     $                abund1,abund2,ptmp1,ptmp2,xpx,opakb1,
     $                opakc,opakscatt,rccemis,fline,lpriu,kdesc2,
     $                r,t,trad,tsq,xee,xh1,xh0,
     $                epi,ncn2,bremsa,bremsint,
     $                rlev,ilev,nlpt,iltp,nlev,klev,lfpi,lun11,
     $                idat1,rdat1,kdat1,nptrs,np2,
     $                npar,npnxt,npfi,npfirst,
     $                nplin,nplini,nlsvn,npcon,npconi,npilev,npilevi,
     $                npconi2,ncsvn,rates,vsav,idrates,lforce)
                    if (lrtyp.eq.1) ans2=0.
                    llo=idest1
                    lup=idest2
                    if (kl.eq.ilimh) lup=min(nlev,lup)
                    if ((idest1.gt.0).and.(idest2.gt.0).and.
     $                (llo.le.nd).and.(lup.le.nd)) then
                      airtmp=ans2
                      rrrt3=rrrt3-airtmp*max(0.,1.-ptmp)
                      ajis(llo,lup)=ajis(llo,lup)+airtmp*ptmp
                      ajis(llo,llo)=ajis(llo,llo)-ans1
                      ajis(lup,lup)=ajis(lup,lup)-airtmp*ptmp
                      ajis(lup,llo)=ajis(lup,llo)+ans1
                      cjis(llo,lup)=cjis(llo,lup)-ans3*xpx
                      cjis(lup,llo)=cjis(lup,llo)+ans4*ptmp*xpx
                      nmat=max(nmat,llo)
                      nmat=max(nmat,lup)
                      rrrtot=rrrtot+airtmp*ptmp
                      if (idest2.le.nlev)
     $                  rrrtot2=rrrtot2+airtmp*ptmp
                      if ((lpri.ge.1))
     $                  write (lun11,9004)jkk,lrtyp,ltyp,idest1,
     $                    idest2,llo,lup,ml,ans1,ans2,ans3,ans4,
     $                    ajis(llo,lup),ajis(lup,llo),
     $                    ajis(llo,llo),ajis(lup,lup),ptmp
9004                    format(1x,8i6,' level',9(1pe10.3),
     $                     7(1pe10.3))
                      endif
                     endif
                   endif
                 if ((lrtyp.eq.5).or.(lrtyp.eq.40)) then
                   call ucalc(ltyp,lrtyp,ml,lcon,jkk,vturbi,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,ans1,ans2,
     $               ans3,ans4,idest1,idest2,idest3,idest4,
     $               abund1,abund2,ptmp1,ptmp2,xpx,opakb1,
     $               opakc,opakscatt,rccemis,fline,lpriu,kdesc2,
     $               r,t,trad,tsq,xee,xh1,xh0,
     $               epi,ncn2,bremsa,bremsint,
     $               rlev,ilev,nlpt,iltp,nlev,klev,lfpi,lun11,
     $               idat1,rdat1,kdat1,nptrs,np2,
     $               npar,npnxt,npfi,npfirst,
     $               nplin,nplini,nlsvn,npcon,npconi,npilev,npilevi,
     $               npconi2,ncsvn,rates,vsav,idrates,lforce)
                   llo=idest1
                   lup=idest2
                   if (kl.eq.ilimh) lup=min(nlev,lup)
                   if ((llo.ne.0).and.(lup.ne.0).and.
     $               (llo.le.nd).and.(lup.le.nd)) then
                     ajis(llo,lup)=ajis(llo,lup)+ans2
                     ajis(llo,llo)=ajis(llo,llo)-ans1
                     ajis(lup,lup)=ajis(lup,lup)-ans2
                     ajis(lup,llo)=ajis(lup,llo)+ans1
                     nmat=max(nmat,llo)
                     nmat=max(nmat,lup)
                     if (lpri.ge.1)
     $                 write (lun11,9004)jkk,lrtyp,ltyp,idest1,idest2,
     $                   llo,lup,ml,ans1,ans2,ans3,ans4,
     $                   ajis(llo,lup),ajis(lup,llo),
     $                   ajis(llo,llo),ajis(lup,lup)
                     endif
                   endif
                 if ((lrtyp.eq.3).or.(lrtyp.eq.23)) then
                   call ucalc(ltyp,lrtyp,ml,lcon,jkk,vturbi,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,ans1,ans2,
     $               ans3,ans4,idest1,idest2,idest3,idest4,
     $               abund1,abund2,ptmp1,ptmp2,xpx,opakb1,
     $               opakc,opakscatt,rccemis,fline,lpriu,kdesc2,
     $               r,t,trad,tsq,xee,xh1,xh0,
     $               epi,ncn2,bremsa,bremsint,
     $               rlev,ilev,nlpt,iltp,nlev,klev,lfpi,lun11,
     $               idat1,rdat1,kdat1,nptrs,np2,
     $               npar,npnxt,npfi,npfirst,
     $               nplin,nplini,nlsvn,npcon,npconi,npilev,npilevi,
     $               npconi2,ncsvn,rates,vsav,idrates,lforce)
                   if ((idest1.gt.0).and.(idest2.gt.0).and.
     $               (idest1.le.nlev).and.(idest2.le.nlev).and.
     $               (idest1.le.nd).and.(idest2.le.nd))
     $               then
                     e1=rlev(1,idest1)
                     e2=rlev(1,idest2)
                     if ((e1/(1.e-24+e2)-1.).lt.0.01) then
                         lup=idest2
                         llo=idest1
                       else
                         lup=idest1
                         llo=idest2
                       endif
                     if (kl.eq.ilimh) lup=min(nlev,lup)
                     ajis(llo,lup)=ajis(llo,lup)+ans2
                     ajis(llo,llo)=ajis(llo,llo)-ans1
                     ajis(lup,lup)=ajis(lup,lup)-ans2
                     ajis(lup,llo)=ajis(lup,llo)+ans1
                     nmat=max(nmat,llo)
                     nmat=max(nmat,lup)
                     if (lpri.ge.1)
     $                 write (lun11,9004)jkk,lrtyp,ltyp,idest1,
     $                   idest2,llo,lup,ml,ans1,ans2,ans3,ans4,
     $                   ajis(llo,lup),ajis(lup,llo),
     $                   ajis(llo,llo),ajis(lup,lup)
                     endif
                   endif
                 if ((lrtyp.eq.4).or.(lrtyp.eq.14).or.(lrtyp.eq.9))
     $             then
                   call ucalc(ltyp,lrtyp,ml,lcon,jkk,vturbi,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,ans1,ans2,
     $               ans3,ans4,idest1,idest2,idest3,idest4,
     $               abund1,abund2,ptmp1,ptmp2,xpx,opakb1,
     $               opakc,opakscatt,rccemis,fline,lpriu,kdesc2,
     $               r,t,trad,tsq,xee,xh1,xh0,
     $               epi,ncn2,bremsa,bremsint,
     $               rlev,ilev,nlpt,iltp,nlev,klev,lfpi,lun11,
     $               idat1,rdat1,kdat1,nptrs,np2,
     $               npar,npnxt,npfi,npfirst,
     $               nplin,nplini,nlsvn,npcon,npconi,npilev,npilevi,
     $               npconi2,ncsvn,rates,vsav,idrates,lforce)
                   if ((idest1.gt.0).and.(idest2.gt.0).and.
     $               (idest1.le.nlev).and.(idest2.le.nlev).and.
     $               (idest2.le.nd).and.(idest1.le.nd)) then
                     e1=rlev(1,idest1)
                     e2=rlev(1,idest2)
                     if (e1.lt.e2) then
                         lup=idest2
                         llo=idest1
                       else
                         lup=idest1
                         llo=idest2
                       endif
                     if (kl.eq.ilimh) lup=min(nlev,lup)
                     jkkl=nplini(ml)
                     if ((jkkl.lt.nnnl).and.(jkkl.gt.0)) then
                       tau1=tau0(1,jkkl)
                       tau2=tau0(2,jkkl)
                       ptmp1=pescl(tau1)*(1.-cfrac)
                       ptmp2=pescl(tau2)*(1.-cfrac)
     $                        +2.*pescl(tau1+tau2)*cfrac
                       ans1=ans1*(ptmp1+ptmp2)
                       ans4=ans4*(ptmp1+ptmp2)
c                      these terms represent radiative exctiation.
c                       must turn them off to attain LTE.
                       ajis(lup,llo)=ajis(lup,llo)+ans2
                       ajis(llo,llo)=ajis(llo,llo)-ans2
                       ajis(lup,lup)=ajis(lup,lup)-ans1
                       ajis(llo,lup)=ajis(llo,lup)+ans1
                       cjis(llo,lup)=cjis(llo,lup)-ans3*xpx
                       cjis(lup,llo)=cjis(lup,llo)+ans4*xpx
                       nmat=max(nmat,llo)
                       nmat=max(nmat,lup)
                       if ((lpri.ge.1))
     $                 write (lun11,9009)jkk,lrtyp,ltyp,idest1,
     $                   idest2,llo,lup,ml,ans1,ans2,ans3,ans4,
     $                   ajis(lup,llo),ajis(llo,lup),
     $                   ajis(llo,llo),ajis(lup,lup),
     $                   tau1,tau2,jkkl,ptmp1,ptmp2
 9009                  format (1x,8i6,' level',10(1pe10.3),
     $                       i6,2(1pe10.3))
                       endif
                     endif
                   endif
                 if (lrtyp.eq.42) then
                   call ucalc(ltyp,lrtyp,ml,lcon,jkk,vturbi,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,ans1,ans2,
     $               ans3,ans4,idest1,idest2,idest3,idest4,
     $               abund1,abund2,ptmp1,ptmp2,xpx,opakb1,
     $               opakc,opakscatt,rccemis,fline,lpriu,kdesc2,
     $               r,t,trad,tsq,xee,xh1,xh0,
     $               epi,ncn2,bremsa,bremsint,
     $               rlev,ilev,nlpt,iltp,nlev,klev,lfpi,lun11,
     $               idat1,rdat1,kdat1,nptrs,np2,
     $               npar,npnxt,npfi,npfirst,
     $               nplin,nplini,nlsvn,npcon,npconi,npilev,
     $               npilevi,npconi2,ncsvn,rates,vsav,idrates,
     $               lforce)
                   if ((idest1.gt.0).and.(idest2.gt.0).and.
     $               (idest1.le.nlev).and.(idest2.le.nlev).and.
     $               (idest2.le.nd).and.(idest1.le.nd))
     $               then
                     e1=rlev(1,idest1)
                     e2=rlev(1,idest2)
                     if (e1.lt.e2) then
                         lup=idest2
                         llo=idest1
                       else
                         lup=idest1
                         llo=idest2
                       endif
                     ajis(lup,llo)=ajis(lup,llo)+ans1
                     ajis(llo,llo)=ajis(llo,llo)-ans1
                     cjis(llo,lup)=cjis(llo,lup)-ans3*xpx
                     nmat=max(nmat,llo)
                     nmat=max(nmat,lup)
                     if ((lpri.ge.1))
     $                 write (lun11,9004)jkk,lrtyp,ltyp,idest1,idest2,
     $                 llo,lup,ml,ans1,ans2,ans3,ans4,
     $                 ajis(lup,llo),ajis(llo,lup),
     $                 ajis(llo,llo),ajis(lup,lup)
                     endif
                   endif
                 if (lrtyp.eq.41) then
                   call ucalc(ltyp,lrtyp,ml,lcon,jkk,vturbi,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,ans1,ans2,
     $               ans3,ans4,idest1,idest2,idest3,idest4,
     $               abund1,abund2,ptmp1,ptmp2,xpx,opakb1,
     $               opakc,opakscatt,rccemis,fline,lpriu,kdesc2,
     $               r,t,trad,tsq,xee,xh1,xh0,
     $               epi,ncn2,bremsa,bremsint,
     $               rlev,ilev,nlpt,iltp,nlev,klev,lfpi,lun11,
     $               idat1,rdat1,kdat1,nptrs,np2,
     $               npar,npnxt,npfi,npfirst,
     $               nplin,nplini,nlsvn,npcon,npconi,npilev,
     $               npilevi,npconi2,ncsvn,rates,vsav,idrates,
     $               lforce)
                   if ((idest1.gt.0).and.(idest2.gt.0).and.
     $               (idest1.le.nlev).and.
     $               (idest2.le.nd).and.(idest1.le.nd))
     $               then
                     lup=idest2
                     llo=idest1
                     ajis(lup,llo)=ajis(lup,llo)+ans1
                     ajis(llo,llo)=ajis(llo,llo)-ans1
                     nmat=max(nmat,llo)
                     nmat=max(nmat,lup)
                     if ((lpri.ge.1))
     $                 write (lun11,9004)jkk,lrtyp,ltyp,idest1,idest2,
     $                   llo,lup,ml,ans1,ans2,ans3,ans4,
     $                   ajis(lup,llo),ajis(llo,lup),
     $                   ajis(llo,llo),ajis(lup,lup)
                     endif
                   endif
                 endif
               ml=npnxt(ml)
               if (ml.ne.0) mlpar=npar(ml)
               enddo
            endif
          endif
        enddo
c
      if (lpri.ne.0)
     $   write (lun11,*)'rrrtt=',rrrtot,rrrtot2
c
c
      return
      end
