subroutine allardt(xl,t,hneutr,hcharg,prof)
c     ===========================================
c
c     quasi-molecular opacity for Lyman alpha, with T-dependent
c     profile
c
c     Input:  xl:  wavelength in [A]
c             hneutr:  neutral H particle density [cm-3]
c             hcharg: ionized H particle density [cm-3]
c     Output: prof:  Lyman alpha line profile, normalized to 1.0e8
c             if integrated over A;
c             It then renormalized by multiplying by
c             8.853e-29*lambda_0^2*f_ij
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      parameter (NXMAX=1400,NNMAX=5,NTAMAX=6)
      parameter (xnorma=8.8528e-29*1215.6*1215.6*0.41618)
      common /calphatd/xlalpd(NXMAX,NTAMAX),plalpd(NXMAX,NNMAX,NTAMAX),
     *       stnead(ntamax),stnchd(ntamax),
     *       vneuad(ntamax),vchaad(ntamax),
     *       talpd(ntamax),nxalpd(ntamax),ntalpd
c
      prof=0.
c
c     find the two partial tables close to actual T
c
      it0=0
      do it=1,ntalpd
         it0=it
         if(t.lt.talpd(it)) then
            it0=it-1
            go to 10
         end if
      end do
   10 continue
      if(it0.eq.0) then
         it0=1
         go to 20
      end if
      if(it0.ge.ntalpd) then
         it0=ntalpd
         go to 20
      end if
      go to 30
   20 continue
c
      if(xl.lt.xlalpd(1,it0).or.xl.gt.xlalpd(nxalpd(it0),it0)) return
      vn1=hneutr/stnead(it0)
      vn2=hcharg/stnchd(it0)
      vns=vn1*vneuad(it0)+vn2*vchaad(it0)
      vn11=vn1*vn1
      vn22=vn2*vn2
      vn12=vn1*vn2
      xnorm=1.0/(1.0+vns+0.5*vns*vns)
c
      jl=0
      ju=nxalpd(it0)+1
 110  if(ju-jl.gt.1) then
         jm=(ju+jl)/2
         if(xl.gt.xlalpd(jm,it0)) then
            jl=jm
         else
            ju=jm
         endif
         go to 110
      endif
      j=jl
c
      if(j.eq.0) j=1
      if(j.eq.nxalpd(it0)) j=j-1
      a1=(xl-xlalpd(j,it0))/(xlalpd(j+1,it0)-xlalpd(j,it0))
      p1=  vn1*((1.0-a1)*plalpd(j,1,it0)+a1*plalpd(j+1,1,it0))
      p11=vn11*((1.0-a1)*plalpd(j,2,it0)+a1*plalpd(j+1,2,it0))
      p2=  vn2*((1.0-a1)*plalpd(j,3,it0)+a1*plalpd(j+1,3,it0))
      p22=vn22*((1.0-a1)*plalpd(j,4,it0)+a1*plalpd(j+1,4,it0))
      p12=vn12*((1.0-a1)*plalpd(j,5,it0)+a1*plalpd(j+1,5,it0))
      prof=(p1+p2+p11+p22+p12)*xnorm*xnorma
      return
c
   30 continue
c
c     interpolate in the tables for different T
c
c     the lower T
c
      if(xl.lt.xlalpd(1,it0).or.xl.gt.xlalpd(nxalpd(it0),it0)) return
      vn1=hneutr/stnead(it0)
      vn2=hcharg/stnchd(it0)
      vns=vn1*vneuad(it0)+vn2*vchaad(it0)
      vn11=vn1*vn1
      vn22=vn2*vn2
      vn12=vn1*vn2
      xnorm=1.0/(1.0+vns+0.5*vns*vns)
      jl=0
      ju=nxalpd(it0)+1
 120  if(ju-jl.gt.1) then
         jm=(ju+jl)/2
         if(xl.gt.xlalpd(jm,it0)) then
            jl=jm
         else
            ju=jm
         endif
         go to 120
      endif
      j=jl
c
      if(j.eq.0) j=1
      if(j.eq.nxalpd(it0)) j=j-1
      a1=(xl-xlalpd(j,it0))/(xlalpd(j+1,it0)-xlalpd(j,it0))
      p1=  vn1*((1.0-a1)*plalpd(j,1,it0)+a1*plalpd(j+1,1,it0))
      p11=vn11*((1.0-a1)*plalpd(j,2,it0)+a1*plalpd(j+1,2,it0))
      p2=  vn2*((1.0-a1)*plalpd(j,3,it0)+a1*plalpd(j+1,3,it0))
      p22=vn22*((1.0-a1)*plalpd(j,4,it0)+a1*plalpd(j+1,4,it0))
      p12=vn12*((1.0-a1)*plalpd(j,5,it0)+a1*plalpd(j+1,5,it0))
      prof0=(p1+p2+p11+p22+p12)*xnorm*xnorma
c
c     the higher T
c
      it0=it0+1
      if(xl.lt.xlalpd(1,it0).or.xl.gt.xlalpd(nxalpd(it0),it0)) return
      vn1=hneutr/stnead(it0)
      vn2=hcharg/stnchd(it0)
      vns=vn1*vneuad(it0)+vn2*vchaad(it0)
      vn11=vn1*vn1
      vn22=vn2*vn2
      vn12=vn1*vn2
      xnorm=1.0/(1.0+vns+0.5*vns*vns)
      jl=0
      ju=nxalpd(it0)+1
 130  if(ju-jl.gt.1) then
         jm=(ju+jl)/2
         if(xl.gt.xlalpd(jm,it0)) then
            jl=jm
         else
            ju=jm
         endif
         go to 130
      endif
      j=jl
c
      if(j.eq.0) j=1
      if(j.eq.nxalpd(it0)) j=j-1
      a1=(xl-xlalpd(j,it0))/(xlalpd(j+1,it0)-xlalpd(j,it0))
      p1=  vn1*((1.0-a1)*plalpd(j,1,it0)+a1*plalpd(j+1,1,it0))
      p11=vn11*((1.0-a1)*plalpd(j,2,it0)+a1*plalpd(j+1,2,it0))
      p2=  vn2*((1.0-a1)*plalpd(j,3,it0)+a1*plalpd(j+1,3,it0))
      p22=vn22*((1.0-a1)*plalpd(j,4,it0)+a1*plalpd(j+1,4,it0))
      p12=vn12*((1.0-a1)*plalpd(j,5,it0)+a1*plalpd(j+1,5,it0))
      prof1=(p1+p2+p11+p22+p12)*xnorm*xnorma
c
c     final profile coefficient
c
      dt=talpd(it0)-talpd(it0-1)
      prof=(prof0*(talpd(it0)-t)+prof1*(t-talpd(it0-1)))/dt
c
      return
      end