      subroutine shock(tshock,taus,beta,tempr,tau,npt)

      real tshock, taus, beta, tempr(*), tau(*)
      integer npt
*
*        Provides temperature and ionization timescale for shock models.
*
*        INPUT:  tshock -  mean temperature at the shock (K)
*                taus   -  ionization parameter nt (s/cm^3)
*                beta   -  electron/mean temperature ratio at the shock
*        OUTPUT: tempr  -  temperatures (K)
*                tau    -  ionization timescale (s/cm^3)
*                npt    -  number of points
*
      real tmin, tmax, taumax, heabun, xmu, xmue, convn, r
      real delf1, tl1, delf, sum, tf, tl, tiz
      integer k, l, i
      real elrat
      external elrat
*        
      tmin=1.0
      tmax=1.0
      npt=2
      tempr(1) = tmin
      tempr(2) = tmax
      tau(1) = 0.0
      tau(2) = 1.0
      taumax = tau(npt)
      if( beta .lt. 1. )then
         heabun = 1.0/10**1.01
         xmu = (1.+4.*heabun)/(2.+3*heabun)
         xmue = (1.+2.*heabun)/(1.+4*heabun)
         convn = 1./(xmu*xmue)
         r = convn*taus/tshock**1.5
         k = 1
         tau(k) = 0.0
         tempr(k) = beta*tmax
         delf1 = 0.01
         tl1 = elrat(r,delf1,beta)*tmax
         do 30 l=2,100
            delf = float(l)*0.01
            tl = elrat(r,delf,beta)*tmax
            if( abs((tl-tempr(k))/tl) .gt. 0.02 )then
               k = k + 1
               tempr(k) = tl1
               tau(k) = delf1
            endif
            delf1 = delf
            tl1 = tl
 30      continue
         k = k + 1
         delf = 1.0
         tempr(k) = elrat(r,delf,beta)*tmin
         tau(k) = taumax
         npt = k
      endif
*
      go to 19
*      tauh = 0.5*tau(npt)
*      call locate(tau,npt,tauh,j)
*      if( (tauh-tau(j))/(tau(j+1)-tau(j)) .lt. 0.5 )then
*         nh = j
*       else
*         nh = j + 1
*      endif
*      if( nh .eq. 1 .or. nh .eq. npt )then
         sum = 0.0
         do 40 i = 2,npt
            sum = sum + 0.5*(tempr(i-1)+tempr(i))*(tau(i)-tau(i-1))
 40      continue
         tiz = sum/tau(npt)
*         tau(1) = 0.0
*         tau(2) = tau(npt)
*         tempr(1) = 2.*tiz - tempr(npt)
*         tempr(2) = tempr(npt)
**         tempr(1) = tempr(2)
*         npt = 2
         tf = tempr(npt)
*         npt = 5
         npt = 11
         do 42 i=1,npt
            tau(i) = float(i-1)/float(npt-1)*taumax
            tempr(i) = tf + 2.*(tiz-tf)*(1.-float(i-1)/float(npt-1))
 42      continue
*         tm = 0.5*( tf + tiz )
*         ti = 2.*tiz - tm
*         tempr(3) = tf
*         tempr(2) = 2.*tm - tf
*         tempr(1) = 2.*ti - tempr(2)
*         tau(2) = 0.5*taumax
*         tau(3) = taumax
*         tau(1) = 0.0
*         npt = 3
*       else
*         tf = tempr(npt)
*         sum = 0.0
*         do 45 i = nh+1,npt
*            sum = sum + 0.5*(tempr(i-1)+tempr(i))*(tau(i)-tau(i-1))
* 45      continue
*         tm = sum/( tau(npt) - tau(nh) )
*         sum = 0.0
*         do 50 i = 2,nh
*            sum = sum + 0.5*(tempr(i-1)+tempr(i))*(tau(i)-tau(i-1))
* 50      continue
*         ti = sum/tau(nh)
*         tempr(3) = tf
*         tempr(2) = 2.*tm - tf
*         tempr(1) = 2.*ti - tempr(2)
*
*         tauh=tau(nh)
*         tau(2) = tauh
*         tau(3) = taumax
*         tau(1) = 0.0
*         npt = 3
*      endif
 19   continue
*        Normalize temperature and ionization timescale to physical units.
      do 20 i = 1,npt
         tempr(i) = tshock*tempr(i)
         tau(i) = taus*tau(i)
 20   continue
      return
      end
*
*
*

c ************************************************************************
      SUBROUTINE doshock(nt, delarr, nionp, ei, ntemp, eisum, t, fout,
     &                   ntempr, tempr, tau, ts, taul, tau0, beta, 
     &                   nzmax, ntmp, ionel, ionstage)

      INTEGER nt, nionp, ntemp, ntempr

      REAL delarr(nt), ei(nionp, nt), eisum(nionp, ntemp)
      REAL tempr(ntempr), tau(ntempr), t(ntemp)
      REAL fout(nionp)

      INTEGER ionel(nionp), ionstage(nionp)

      REAL ts, taul, tau0, beta

      INTEGER nzmax, ntmp


      integer nem, iem, ntt
      parameter ( nem = 400 )
      parameter ( ntt = 5 )


c       the lower and upper limiting values of log temperature,
c               the increment in log temp.
c
      real tll,tlu,dtl
      parameter (tll=4.0e0,tlu=9.0e0,dtl=4.0e-3)

      integer lt, it, i
      integer ion, iion, izone
      real eitot, sumem, sumt, dela
      real taumin, deltau, tauim1, taui, tausiz
      integer npt, n
      real delem
      real tt

          do 198 it = 1,nt
            delarr(it) = 0.0
            do 199 ion = 1,nionp
               ei(ion,it) = 0.0
 199        continue
 198     continue
      if( taul .eq. 0.0 )then
         n=200
         taumin = 1.e8
         deltau = alog10(tau0/taumin)/float(n-1)
         tauim1 = 0.0
       else
         n=20
         deltau = ( tau0 - taul )/float(n)
         tauim1 = taul
      endif         
      do 4 izone=1,n
         if( taul .eq. 0.0 )then
            taui = taumin*10**( float(izone-1)*deltau )
            delem = ( taui - tauim1 )/tau0
          else
            taui = taul + float(izone)*deltau
            delem = ( taui - tauim1 )/( tau0 - taul )
         endif
         tausiz = 0.5*(taui+tauim1)
         tauim1 = taui
         call shock(ts,tausiz,beta,tempr,tau,npt)
         call ionsneqr(tempr,tau,npt,nzmax,nionp,fout,ionel,
     &                 ionstage)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       calculate temperature index
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
            lt=(alog10(tempr(npt))-tll)/dtl+1.5
            if (lt.lt.1)  lt=1
            if (lt.gt.nt) lt=nt
            it = lt
            do 1200 iion=1,nionp
 1200       ei(iion,it)=fout(iion)*delem+ei(iion,it)
            delarr(it) = delarr(it) + delem
 200     continue
 4    continue
*          Rebin temperature zones for spectral calculations.
        ntmp = 0
        do 115 i=1,nt-ntt+1,ntt
           eitot = 0.0
           do 116 it=i,i+ntt-1
              if( delarr(it) .gt. 0.0 ) eitot = 1.0
 116       continue
           if( eitot .gt. 0.0 )then
              ntmp = ntmp + 1
              sumem = 0.0
              sumt = 0.0
              do 117 ion=1,nionp
                 eisum(ion,ntmp) = 0.0
 117          continue
              do 118 it=i,i+ntt-1
                 tt=it-1
                 tt=tll+dtl*tt
                 tt=10.0e0**tt
                 dela = delarr(it)
                 sumem = dela + sumem
                 sumt = dela*tt + sumt
                 do 119 ion=1,nionp
                    eisum(ion,ntmp) = ei(ion,it) + eisum(ion,ntmp)
 119             continue
 118          continue
              t(ntmp) = sumt/sumem
           endif
 115    continue

        RETURN
        END
