      subroutine stpcut(ldirt,lpri,lun11,vturbi,
     $       idat1,rdat1,kdat1,nptrs,np2,
     $       npar,npnxt,npfi,npfirst,
     $       nplin,nplini,nlsvn,npcon,npconi,npilev,npilevi,
     $       npconi2,ncsvn,
     $      epi,ncn2,opakc,oplin,opakab,delr,t,
     $      dpthc,tau0,tauc,eliml,elimh)
c
c     this routine updates.  calculates depths, etc.
c     author:  T. Kallman
c
      implicit none
c
      include './PARAM'
      integer nbtpp
      parameter (nbtpp=10000)
      integer nptmpdim
      parameter (nptmpdim=200000)
c
c     master data
      integer idat1(nidat1),
     $      nptrs(nptt,ndat2)
      real rdat1(nrdat1)
      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 npcon(nnml),npconi2(ndat2),npconi(ndat2)
      integer npilev(nd,nni),npilevi(nnml)
      integer nplin(nnnl),nplini(nnnl)
      real epi(ncn),opakc(ncn),dpthc(2,ncn)
      real tauc(2,nnml)
      real opakab(nnml)
c     line opacities
      real oplin(nnnl)
c     line optical depths
      real tau0(2,nnnl)
      real optmp(ncn)
c     temporaries for dread
      real rdat(nptmpdim)
      integer idat(nptmpdim)
      integer ldon(2)
c
      character*1 kdat(nptmpdim)
      real vturbi,delr,t,eliml,elimh,
     $  dpcrit,bbb,optpp,delea,aatmp,elin,etmp,vth,
     $  vturb,deleturb,deleth,dele,aasmall,
     $  deleused,deleepi,delet,deletpp,e00,dpthmx,
     $  dpthtmp,e0,tst,opsum,optmpo,profile,optp2,
     $  tmpopmx,tmpopo,etptst
      integer ldirt,lpri,lun11,
     $  ncn2,llk,lnn,ln,ml,lup,nilin,nelin,
     $  nbtmp,iion,nitmp,ndtmp,mllz,
     $  iltmp,nlsvn,nlsv,i,ij,lind,
     $  lcon,ldir,mlm,ml1m,ltyp,lrtyp,ml1,ml2,ml1min,
     $  ml1max,mlc,mloff,mlmin,mlmax,ncut,nidt,
     $  nkdt,nrdt,np2,ncsvn
      integer nbinc,mlpar
      real voigte
c     arrays containing printout info
      integer ilsv(nnnl)
      real ewsv(nnnl)
      real oplin2(nnnl)
c
c     temporary grid for use in calculating profile
      real etpp(nbtpp),optpp2(nbtpp)
c
c      real*8 tmpew,tmpewo,tmpop,tmpe,sum,sume
      real tmpew,tmpewo,tmpop,tmpe,sum,sume
c
c     Not used
      integer javi
c
      javi=np2
      javi=npfirst(1)
      javi=nplini(1)
      javi=npcon(1)
      javi=npconi(1)
      javi=npilev(1,1)
      javi=npilevi(1)
      javi=npconi2(1)
      npconi2(1)=javi
c
      if (lpri.ne.0)
     $ write (lun11,*)'in stpcut:',delr,nlsvn
c
      do llk=1,ncn2
        optmp(llk)=0.
        enddo
c
c      go to 9000
c
      bbb=vturbi
      dpcrit=1.e-9
      nlsv=0
      do  lnn=1,nlsvn
        ln=lnn
        ml=nplin(ln)
        oplin2(ln)=0.
        call dread(ltyp,lrtyp,lcon,
     $          nrdt,rdat,nidt,idat,nkdt,kdat,ml-1,
     $          idat1,rdat1,kdat1,nptrs,0,lun11)
c       get damping parameter for fe 2-3 uta
        delea=0.
        if (ltyp.eq.82) delea=rdat(6)*(4.14e-15)
c       get damping parameter for type 50 data
        if (ltyp.eq.50) delea=rdat(3)*(4.14e-15)
        elin=abs(rdat(1))
        lup=idat(2)
        nilin=npar(ml)
        ml=nilin
        call dread(ltyp,lrtyp,lcon,
     $       nrdt,rdat,nidt,idat,nkdt,kdat,nilin-1,
     $       idat1,rdat1,kdat1,nptrs,0,lun11)
        nilin=idat(3)
c       get nuclear mass       
        nelin=npar(ml)
        ml=nelin
        call dread(ltyp,lrtyp,lcon,
     $          nrdt,rdat,nidt,idat,nkdt,kdat,ml-1,
     $          idat1,rdat1,kdat1,nptrs,0,lun11)
        aatmp=rdat(2)
        optpp=oplin(ln)
c     $      *deleth/max(deleth,deleturb)
c     $      *(vth/3.e+5)
        if (lpri.ne.0)
     $   write (lun11,*)ln,elin,aatmp,nilin,nelin,lup,
     $                  optpp,delr
        if ((ln.gt.0).and.(ln.le.nnnl) 
     $    .and.(elin.gt.eliml).and.(elin.lt.elimh) 
     $    .and.(aatmp.gt.0.1).and.(optpp*delr.gt.dpcrit)
     $    .and.(nilin.gt.0).and.(nilin.le.nni))
     $    then

c         line parameters
          etmp=12398.54/elin
          nbtmp=nbinc(etmp,epi,ncn2)
          if (lpri.ne.0)
     $     write (lun11,*)ln,elin,nilin,lup

c         find associated type 86 data
          iion=0
          idat(nidt)=0
          if (lpri.ne.0)
     $     write (lun11,*)'searching for ion'
          do while ((idat(nidt).ne.nilin).and.(iion.lt.nni))
            iion=iion+1
            nitmp=npfi(13,iion)
            call dread(ltyp,lrtyp,lcon,
     $            nrdt,rdat,nidt,idat,nkdt,kdat,nitmp-1,
     $           idat1,rdat1,kdat1,nptrs,0,lun11)
            if (lpri.gt.1)
     $        write (lun11,*)iion,idat(nidt),nilin,nitmp
            enddo
          if (lpri.gt.0)
     $        write (lun11,*)iion,idat(nidt),nilin,nitmp
c         get damping parameter for iron auger damped lines
          ndtmp=npfi(41,iion)
          if (ndtmp.ne.0) then
            if (lpri.ne.0)
     $       write (lun11,*)'  found ion',lup,ndtmp
            if (ndtmp.gt.0) then
              mllz=npar(ndtmp)
              call dread(ltyp,lrtyp,lcon,
     $          nrdt,rdat,nidt,idat,nkdt,kdat,ndtmp-1,
     $          idat1,rdat1,kdat1,nptrs,0,lun11)
              iltmp=idat(2)
              mlpar=mllz
              do while ((ndtmp.ne.0).and.(lup.ne.iltmp)
     $          .and.(mlpar.eq.mllz)) 
                call dread(ltyp,lrtyp,lcon,
     $            nrdt,rdat,nidt,idat,nkdt,kdat,ndtmp-1,
     $            idat1,rdat1,kdat1,nptrs,0,lun11)
                iltmp=idat(2)
                if (lpri.gt.1)
     $           write (lun11,*)'   ',nidt,iltmp,ndtmp
                ndtmp=npnxt(ndtmp)     
                mlpar=0
                if (ndtmp.ne.0) mlpar=npar(ndtmp)
                enddo
              if (lpri.gt.0)
     $           write (lun11,*)'   ',nidt,iltmp,ndtmp
              endif
            if (lup.eq.iltmp) delea=rdat(3)*(4.14e-15)
            endif
c
c         cheat for narrow line plot
c         delea=0.        
c
c         thermal width quantities
          vth=(1.2e+1)*sqrt(t/aatmp)
          vturb=max(bbb,vth)
          e0=etmp
          deleturb=e0*(vturb/3.e+5)
          deleth=e0*(vth/3.e+5)
          dele=deleth+deleturb
          aasmall=delea/(1.e-24+dele)/12.56
          ml1=nbtmp
c         NB this statement forces line center to lie on grid
c          etmp=epi(ml1)
          if (lpri.ne.0)
     $     write (lun11,*)'optpp,dpcrit*opakc(ml1),ml1,deleth,delea:',
     $                    optpp,dpcrit*opakc(ml1),ml1,deleth,delea
          delet=(epi(ml1)-etmp)/dele
          if (optpp*delr.gt.dpcrit)  then
c
            if (lpri.ne.0) write (lun11,*)'calculating profile',
     $                     ml1,epi(ml1),epi(ml1+1),optpp
c
c           calculate profile on temporary grid
            e00=epi(ml1)
            deleepi=epi(ml1+1)-epi(ml1)
            deletpp=dele/8.
            ncut=int(deleepi/deletpp)
            ncut=max(ncut,2)
            deleused=deleepi/float(ncut)
            if (lpri.ne.0) write (lun11,*)'ncut=',ncut,deleused,deletpp,
     $                                    deleepi
c           apply binned opacity
            delet=(e00-etmp)/dele
            if (aasmall.gt.1.e-6) then
                 profile=voigte(abs(delet),aasmall)
              else
                 profile=exp(-delet*delet)/1.772
              endif 
            ml2=nbtpp/2
            etpp(ml2)=e00
            optpp2(ml2)=optpp*profile
            mlc=0
            ldir=1
            ldon(1)=0
            ldon(2)=0
            mlmin=nbtpp
            mlmax=1
            ml1min=ncn+1
            ml1max=0
            do while ((ldon(1)*ldon(2).eq.0).and.(mlc.lt.nbtpp/2)) 
              mlc=mlc+1
              do ij=1,2
                ldir=-ldir
                if (ldon(ij).ne.1) then
                  mlm=ml2+ldir*mlc
                  etptst=e00+float(ldir*mlc)*deleused
                  if ((mlm.le.nbtpp).and.(mlm.ge.1)
     $             .and.(etptst.gt.0.).and.(etptst.lt.epi(ncn2))) then
                    mloff=abs(mlm-ml2)
                    if (mod(mloff,ncut).eq.0) then
                      ml1m=ml1+ldir*(mlc/ncut)
                      ml1min=min(ml1m,ml1min)
                      ml1max=max(ml1m,ml1max)
                      endif
                    mlmin=min(mlm,mlmin)
                    mlmax=max(mlm,mlmax)
                    etpp(mlm)=e00+float(ldir*mlc)*deleused
                    delet=(etpp(mlm)-etmp)/dele
                    if (aasmall.gt.1.e-9) then
                       profile=voigte(abs(delet),aasmall)
                      else
                        profile=exp(-delet*delet)/1.772
                      endif 
                    optpp2(mlm)=optpp*profile
                    tst=optpp2(mlm)*delr
                    if (lpri.ne.0) write (lun11,*) mlm,etpp(mlm),ij,
     $               deleused,delet,mlmin,mlmax,ml1,ml1min,ml1max,
     $               mlc,mloff,mod(mloff,ncut),profile,optpp2(mlm),tst
                    endif
                  if (((tst.lt.dpcrit)
     $                 .or.(mlm.le.1).or.(mlm.ge.nbtpp)
     $                 .or.(etptst.le.0.).or.(etptst.ge.epi(ncn2))
     $                 .or.(mlc.gt.nbtpp).or.(abs(delet).gt.50.)) 
     $                 .and.(ml1min.lt.ml1-2).and.(ml1max.gt.ml1+2)
     $                 .and.(ml1min.ge.1).and.(ml1max.le.ncn))
     $                  ldon(ij)=1
                  endif
                enddo                 
              enddo
            sum=0.
            opsum=0.
            tmpop=0.
            tmpopmx=0.
            sume=0.
            ml1m=ml1min
            if (lpri.ne.0) write (lun11,*)'renormalizing profile'
            nlsv=nlsv+1
            ilsv(nlsv)=ln
            ewsv(nlsv)=0.
            tmpew=0.
            mlmin=max(mlmin,2)        
            mlmax=min(mlmax,nbtpp)        
            ml1m=nbinc(etpp(mlmin),epi,ncn2)
            do mlm=mlmin,mlmax
              mloff=mlm-ml2
              if (etpp(mlm).gt.epi(ml1m)) then                  
                  tmpopo=tmpop
                  tmpop=optpp2(mlm)*delr
                  tmpopmx=max(tmpopmx,tmpop)
                  tmpewo=tmpew
                  tmpew=exp(-tmpop)
                  tmpe=abs(etpp(mlm)-etpp(mlm-1))
                  if (mlm.gt.mlmin+1) then
                    sum=sum+(tmpew+tmpewo)*tmpe/2.
                    sume=sume+tmpe
                    opsum=opsum+(tmpop+tmpopo)*tmpe/2.
                    endif
                  optmpo=optmp(ml1m)
                  if (sume.gt.1.d-34) then
                    if (tmpopmx.gt.1.e-4) then
                        optp2=-log(max(1.d-99,sum)/sume)/delr
                      else
                        optp2=opsum/sume/delr
                      endif
                    do while (etpp(mlm).gt.epi(ml1m)) 
                      optmp(ml1m)=optmp(ml1m)+optp2
c                      if (ml1m.eq.nbtmp+1) then
c                        oplin2(ln)=optp2
c                        endif
                      ml1m=ml1m+1
                      enddo
                    endif
                  if (lpri.ne.0) write (lun11,*)mlm,mloff,ml1m,
     $               epi(ml1m),epi(ml1m+1),etpp(mlm),optmp(ml1m),
     $               optmpo,optpp2(mlm),sum,sume,ewsv(nlsv)
                  ewsv(nlsv)=ewsv(nlsv)+sume-sum
                  tmpopmx=0.
                  opsum=0.
                  sum=0.
                  sume=0.
                else
                  tmpopo=tmpop
                  tmpop=optpp2(mlm)*delr
                  tmpopmx=max(tmpopmx,tmpop)
                  tmpewo=tmpew
                  tmpew=exp(-tmpop)
                  tmpe=abs(etpp(mlm)-etpp(mlm-1))
                  sum=sum+(tmpew+tmpewo)*tmpe/2.
                  sume=sume+tmpe
                  opsum=opsum+(tmpop+tmpopo)*tmpe/2.
c                  if (lpri.ne.0) write (lun11,*)mlm,mloff,ml1m,
c     $                epi(ml1m),etpp(mlm),optpp2(mlm),sum,sume
                endif
              enddo
            endif      
          endif
        enddo
c
 9000 continue
c
c     calculate continuum depths
      if (lpri.ne.0) write (lun11,*)'calculating depths in stpcuta'
      lind=1
      if (ldirt.gt.0) lind=2
      dpthmx=0.
      optpp=0.
      do  i = 1,ncn2
c         opakc(i)=opakc(i)+optmp(i)
         if (lpri.ne.0) write (lun11,*)i,epi(i),opakc(i),optmp(i)
         optpp=min(optmp(i),1.e+3/(1.e-24+delr))
         dpthtmp=(opakc(i)+optpp)*delr
         dpthc(lind,i) = dpthc(lind,i) + dpthtmp
         if (dpthtmp.gt.dpthmx) then
           dpthmx=dpthtmp
           endif
         enddo
c
c     calculate line depths
      do  i = 1,nlsvn
         tau0(lind,i) = tau0(lind,i) + oplin(i)*delr
         if (lpri.ne.0)
     $    write (lun11,*)i,lind,oplin(i),delr,tau0(lind,i)
         enddo
c
c     calculate level depths
c      write (lun11,*)'in stpcut:',delr,lind
      do i = 1,ncsvn
         tauc(lind,i) = tauc(lind,i) + opakab(i)*delr
c         write (lun11,*)i,opakab(i),tauc(lind,i)
         enddo
c
c
      return
      end
