subroutine neu(u2,t2)
*
* 1.  Photo-, Plasma, Pair, and Recombination neutrino emissivities are
*     calculated following
*     Beaudet, Petrosian, & Salpeter 1967, ApJ, 150, 979.
* 2.  Bremsstrahlung neutrino emissivity is calculated following
*     Festa & Ruderman 1969, Phys. Rev., 180, 1227.
* 3.  See also:
*     Cazzola, De Zotti, & Saggion 1971, Phys. Rev., D3, 1722
*     Flowers 1973, ApJ, 180, 911.
* 4.  Note: According to more recent studies, these overestimate the real 
*     neutrino rates by factors of 2 or more.
*
      implicit double precision(a-h,o-z)

      common/neut/en(5),fn
      common/temp/s1,r1,s2,r2,b2,p2,tx2,e2,xc2,xo2,fca2,f2,q2,w2,c
      common/je/je
      common/c/az(3),zz3(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x
     1,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi

      data a1,a2,a3,a4/5.167e-1,-3.104e0,-2.801e0,7.577e0/
      data a10,a11,a12,b11,b12,b13,c1/2.320e-7,8.449e-8,1.787e-8,
     12.581e-2,1.734e-2,6.990e-4,5.6457e-1/
      data a20,a21,a22,b21,b22,b23,c2/4.886e10,7.580e10,6.023e10,
     16.290e-3,7.483e-3,3.061e-4,1.5654e00/
      data a30,a31,a32,b31,b32,b33,c3/6.002e19,2.084e20,1.872e21,
     19.383e-1,-4.141e-01,5.829e-2,5.5924e00/

      z=zz3(je)
      if ( je .eq. 1 ) then
         ce=1.
      else
         ce=2.
         if ( je .eq. 3 ) then
            azx = 1. / ( xc2/12. + xo2/16. )
            z   = azx / 2.
         endif
      endif
      rho=10.**u2/ce 
*
* x is lambda and y is xi from BPS
*
      x=10.**(t2-9.77307)
      y=rho**(1./3)/(1000*x)
      yyy=y*y*y
      xx=x*x
      zz=z*z
      xyy=x*y*y
      do 1 i=1,5
        en(i)=0
1     continue
      fn=0
*
*     plasma neutrinos (dominant at most luminosities)
*
      if(x.lt.1e-3) go to 8
      if(y.gt.60) go to 5
      fpl=dexp(-c1*y)*(a10+y*(a11+y*a12))/(yyy+(b11+(b12+b13/x)/x)/x)
      en(1)=rho*rho*fpl/ce
*
*     photoneutrinos (negligible for us)
*
      if(y.gt.15) go to 3
      fph=dexp(-c2*y)*(a20+y*(a21+y*a22))/(yyy+(b21+(b22+b23/x)/x)/x)
      en(2)=xx*xx*x*fph/ce
*
*     pair neutrinos (noticable at high luminosities)
*
    3 if(x.lt.4e-2.or.y.gt.5) go to 4
      fpa=dexp(-c3*y)*(a30+y*(a31+y*a32))/(yyy+(b31+(b32+b33/x)/x)/x)
      g=1-xx*(13.04-xx*(133.5+xx*(1534.+xx*918.6)))
      en(3)=g*fpa/(rho*ce*dexp(2.d0/x))
*
*     recombination neutrinos (a small contributor)
*
4     if(y.gt.10) go to 5
      u=24.280*xyy*ce**(2./3) 
      if(u.gt.30) go to 5
      v=zz/(3.7772e04*x)
      frc=xx/dexp(u+v)
      en(4)=1.850e-4*zz*zz*rho*frc/ce
*
*     gas bremsstrahlung neutrinos (dominant at low L)
*
5     u=u2-dlog10(ce)
      v=a1*u
      fbr=4.775e-1*atan(1.707e0*(u-t2)+4.566e0)+5.250e0
      w=a2+v+(t2-8)*fbr-dsqrt(v*v+a3*u+a4)
      en(5)=(2/ce)*z*10.**w

      do 7 i=1,5
         fn=fn+en(i)
7     continue
8     return
      end 

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