subroutine gnutrino(u2,t2)
*
*  Neutrino rates from Itoh, Adachi, Nakagawa, Kohyama, & Munkata, 1989, 
*  ApJ, 339, 354. (+Errata 1990, ApJ, 360, 741)
*    (Plasma, Photo, Pair)
* Recombination neutrinos from Kohyama, Itoh, Obama, & Mutoh 1993, ApJ,
*  415, 267
*  Bremsstrahlung from:
* 1. Liquid Metal: Itoh & Kohyama, 1983, ApJ, 275, 858.
* 2. Crystal Lattice: Itoh etal., 1984, ApJ, 279, 413.
* 3. Low T Quantum Corr.: Itoh etal., 1984, ApJ, 280, 787.
* 4. Phonon Contr.: Itoh etal., 1984, ApJ, 285, 304.
* 5. Partially Degen. e: Munkata, Kohyama & Itoh, 1987, ApJ, 316, 708.
* Prepared for WDXD by PAB: Oct. 1990
* Last messed with on 08 17 94
*
       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/comp/amhyhe,amheca,x(4)
      common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x
     1,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi
       common/rhomooe/rhomooe
       common/phys/flam,gamma,degtemp,tfermi
       double precision mue
       data gamc,gamo,gamhe,gamh/3577200.,5778200.,573264.,227500./
*
* initialize
*
       t=10.d0**t2
       rho=10.d0**u2
       do 10 i=1,5
         en(i) = 0.0
10     continue
       fn = 0.0
       qbrem = 0.0
       qplas = 0.0
       qphot = 0.0
       qpair = 0.0
       qrec  = 0.0
*
* Escape clause (enter when T < 10^7K or Rho < 1 g/cc)
* With this, I don't need an escape clause for flam
*
       if(t2.lt.7.0 .or. u2.lt.0.0)then
         return
       endif
*
* compute gamma and mue for mixture
*
       if(je.eq.1)then
*
* assume pure H for now. (Note: for H, I shouldn't be here anyway!)
*
         ce = 1.
         mue = 1.
         gamma = gamh * 10.**(u2/3. - t2)
         zbar = 1.0
         aavg = 1.0
        else
         if(je.eq.2)then
*
* pure helium
*
           omue = 1./(x(2)/2.)
           gamma = gamhe * 10.**(u2/3. - t2)
           zbar = 2.0
	   aavg = 4.0
          else
*
* helium, carbon, oxygen, and iron (the latter for tests)
*
           gammahe = gamhe * 10.**(u2/3. - t2)
           gammac  = gamc  * 10.**(u2/3. - t2)
           gammao  = gamo  * 10.**(u2/3. - t2)
           gamma = xc2*gammac + xo2*gammao
           omue = xc2*6./12. + xo2*8./16.
           zbar = xc2*6.0 + xo2*8.0
	   aavg = xc2*12.0 + xo2*16.0
         endif
       endif
*
* zbar is the mean charge on a nucleus, for bremm and pair.  
* flam was x in neu 
*
       flam=t/5.9302d+9
       mue = 1.d0 / omue
       rhomooe = rho/mue
       rhomue23 = (rhomooe)**(2./3.)
       sqfac = 1. + 1.018*rhomue23
       tfermi = 5.9302d+9*(dsqrt(sqfac) - 1.)
       degtemp = t/tfermi
       y = (rhomooe)**(1./3.)/(1000.*flam)
       s2theta = 0.217
       cv=0.5d0+2.0d0*s2theta
       cvp=1.d0-cv
       ca=0.5d0
       cap=1.d0-ca
       if(y .le. 100.)then
         call plasmon(flam,rho,t2,mue,cv,cvp,ca,cap,qplas)
       endif
       if(y .le. 100.)then
         call photon(flam,t2,rho,mue,cv,cvp,ca,cap,qphot)
       endif
       if(flam.ge.0.10)then
         call pairn(flam,rho,mue,cv,cvp,ca,cap,qpair)
       endif
       call recombine(aavg,zbar,t,rho,cv,cvp,ca,cap,qrec)       
       if(u2.gt.4.0 .and. degtemp.gt.0.3)then
         call brempde(u2,t2,gamma,qbrempde)
         qbrem = qbrempde
         if(qbrem.lt.0.0)then
           qbrem = 0.0
         endif
*
* skip over other neutrino rates if we came here.
*
         go to 50
       endif
       if(u2.ge.1.0)then
         if(gamma.gt.1.0 .and. gamma.lt.171.0)then
           if(x(2).gt.0.00001)then
             call bremliqh(u2,t2,gamma,qbremlh)
            else 
             qbremlh = 0.0
           endif
           if(x(3).gt.0.00001)then
             call bremliqc(u2,t2,gamma,qbremlc)
            else 
             qbremlc = 0.0
           endif
           if(x(4).gt.0.00001)then
             call bremliqo(u2,t2,gamma,qbremlo)
            else 
             qbremlo = 0.0
           endif
           qbrem = qbremlh*x(2) + qbremlc*x(3) + qbremlo*x(4)
         elseif(gamma.ge.171.0)then
           if(x(3).gt.0.00001)then
             call bremsolc(u2,t2,gamma,qbremsc)
            else 
             qbremsc = 0.0
           endif
           if(x(4).gt.0.00001)then
             call bremsolo(u2,t2,gamma,qbremso)
            else 
             qbremso = 0.0
           endif
           qbrem = qbremsc*x(3) + qbremso*x(4)
           if(qbrem.lt.0.0)then
             qbrem = 0.0
           endif
         endif
       endif
*
* convert to lumonisities
*
50     eplas=qplas/rho
       en(1) = eplas
       ephot=qphot/rho
       en(2) = ephot
       epair=qpair/rho
       en(3) = epair
       erec=qrec/rho
       en(4) = erec
*
* output of Bremsstrahlung routines already is divided by rho
*
       ebrem=qbrem
       en(5) = ebrem
       etot=eplas+ephot+epair+erec+ebrem
       fn=etot
       return
       end

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