subroutine pairn(flam,rho,mue,cv,cvp,ca,cap,qpair)
*
* Pair production neutrinos
* Munkata etal version is unchanged except for log T > 10, which
* we better not hit, although I've included the appropriate coefficients.
*
       implicit double precision (a-h,o-z)

       dimension a(0:2),b(3)
       double precision mue
       xi=(((rho/mue)/1e9)**(1./3.))/flam
       a(0)=6.002e19
       a(1)=2.084e20
       a(2)=1.872e21
       n = 2
*
* pick out appropriate b's and c's according to temp
* flam = 1.6863 is T = 10^10 K
*
       if(flam.lt.1.6863)then
         b(1)=9.383e-1
         b(2)=-4.141e-1
         b(3)=5.829e-2
         c=5.5924
        else
         b(1)= 1.2383
         b(2)=-8.141e-1
         b(3)=0.0
         c=4.9924
       endif
       call qcalc(a,b,c,xi,flam,f)
       g=1.d0-13.04d0*(flam**2)+133.5d0*(flam**4)+1534d0*(flam**6)
     1       +918.6d0*(flam**8)
       gef=g*dexp(-2.d0/flam)*f
       temp=1.d0/(10.7480d0*flam**2+0.3967d0*flam**0.5d0+1.0050d0)
       smq=temp*(1.d0+(rho/mue)/
     1     (7.692d7*flam**3+9.715d6*flam**0.5d0))**(-0.3)
       cvm=(cv**2-ca**2)+n*(cvp**2-cap**2)
       cvp=(cv**2+ca**2)+n*(cvp**2+cap**2)
       qpair=0.5d0*cvp*(1.d0+(cvm/cvp)*smq)*gef
       return
       end

************************************************************************