c
       subroutine neutrino(dl,tl,x,y,c,o,z,
     1                     eplas,ephot,epair,erec,ebrem,etot)
       implicit double precision (a-h,o-z)
       double precision mue
c
c  Neutrino rates from Itoh et al. (1996), Ap.J. Supp 102, 411
c    (plama, photo, pair, bremmstrahlung)
c  and Beaudet, Petrosian, and Salpeter (1967), Ap. J. 150, 979.
c    (recomb.)
c
c  Standardized for ISUEVO by SDK: December 1990
c  original BPS updated to Itoh: September 1996
c
       cln=dlog(10.d0)
c       t=10.d0**tl
       t=dexp(cln*tl)
c       rho=10.d0**dl
       rho=dexp(cln*dl)
       mue=2.d0/(1.d0+X)
       zbar=x+2.d0*y+6.d0*c+8.d0*o
       if (t.lt.1.d7.or.rho.lt.1.d3) then
         eplas=0.d0
         ephot=0.d0
         epair=0.d0
         erec=0.d0
         ebrem=0.d0
         etot=0.d0
         return
       endif
c
       flam=t/5.9302d9
       cv=0.5d0+2.0d0*0.232d0
       cvp=1.d0-cv
       ca=0.5d0
       cap=1.d0-ca
       temp=1.018d0*dexp(0.6666667*dlog(rho/mue))
       tf=5.9302d9*(dexp(0.5*dlog(1.d0+temp))-1)
       call PLASMA(flam,rho,mue,t,cv,cvp,ca,cap,qplas)
       call BREMM(t,rho,mue,x,y,c,o,zbar,cv,cvp,ca,cap,qbrem)
       call PAIR(flam,rho,mue,cv,cvp,ca,cap,qpair)
       call PHOTO(flam,rho,t,mue,cv,cvp,ca,cap,qphot)
       call RECOMB(flam,t,rho,zbar,qrec)
       eplas=qplas/rho
       ephot=qphot/rho
       epair=qpair/rho
       erec=qrec/rho
       ebrem=qbrem/rho
       etot=eplas+ephot+epair+erec+ebrem
       return
       end
c
       subroutine PLASMA(flam,rho,mue,t,cv,cvp,ca,cap,qplas)
c from Itoh et al. 1996, checks based on plots in that paper and tables
c in the CDROM version
       implicit double precision (a-h,o-z)
       double precision mue
       top=1.1095d11*rho/mue 
       temp=dexp(0.6666667d0*dlog(1.019d-6*rho/mue))
       bot=t*t*dexp(0.5d0*dlog(1+temp))
       gam=dexp(0.5d0*dlog(top/bot))
       ft=2.4d0+0.6d0*dexp(0.5d0*dlog(gam))+0.51d0*gam
     1         +1.25d0*dexp(1.5d0*dlog(gam))
       fl=(8.6d0*gam*gam+1.35d0*dexp(3.5d0*dlog(gam)))/
     1    (225d0-17*gam+gam*gam)
       x=0.16666667*(17.5+dlog10(2*rho/mue)-3*dlog10(t))
       y=0.16666667*(-24.5+dlog10(2*rho/mue)+3*dlog10(t))
       n=2
       if (x*x.gt.0.49d0.or.y.lt.0) then
         fxy=1.d0
       else
         fxy=1.05d0+(0.39-1.25*x-0.35d0*sin(4.5d0*x)
     1               -0.3d0*dexp(-(4.5d0*x+0.9d0)**2))
     2         *exp(-(min(0.d0,y-1.6d0+1.25*x)/(0.57d0-0.25d0*x))**2)
       endif
       qv=3.00d21*flam**9*gam**6*dexp(-gam)*(ft+fl)*fxy
       qplas=(cv**2+n*cvp**2)*qv
       return
       end
c
       subroutine BREMM(t,rho,mue,x,y,c,o,z,cv,cvp,ca,cap,qbrem)
c Itoh et al 1996; checked with their figures 6 and 7.  For carbon,
c doesn't do well at high density regime, but okay elsewhere
c requires separate computation routines for He, C,and O
       implicit double precision (a-h,o-z)
c
       double precision mue
c
       cln=2.302585093d0
       t8=t/1.d8
       rome=rho/mue
       temp=1.018d0*dexp(0.6666667*dlog(rho/1d6/mue))
       tf=5.9302d9*(dexp(0.5*dlog(1.d0+temp))-1.d0)
c
       if ((y.gt.0.5d0.and.t.gt.0.01d0*tf).or.(t.gt.0.3d0*tf)) then
c weakly degenerate electrons
         a0=23.5d0
         a1=6.83d4
         a2=7.81d8
         a3=230.d0
         a4=6.75d5
         a5=7.66d9
         b1=1.47d0
         b2=0.0329
         b3=7.75d5*dexp(1.5d0*dlog(t8))
     1             +247d0*dexp(3.85d0*dlog(t8))
         b4=4.07d0+0.0240d0*dexp(1.4d0*dlog(t8))
         b5=4.57d-5*dexp(-0.11d0*dlog(t8))
c
         bot1=a3+a4/(t8*t8)+a5/(dexp(5.d0*dlog(t8)))
         bot1=bot1*(1.d0+1.d-9*rome)
         bot2=b3/rome+b4+b5*dexp(0.656d0*dlog(rome))
         ggas=1.d0/bot1 + 1.d0/bot2
         eta=rome/(7.05d6*dexp(1.5d0*dlog(t8))+5.12d4*t8*t8*t8)
c
         bot1=a0+a1/t8/t8+a2*dexp(-5d0*dlog(t8)) 
         bot2=1.d0+b1/eta+b2/eta/eta
         fgas=1.d0/bot1+1.26d0*(1.d0+1.d0/eta)/bot2
c
         n=2
c zbar is the mean value (mass-weighted) of the nuclear charge Z
         zbar=x+2.d0*y+6.d0*c+8.d0*o
         qgas=0.5738d0*(2.d0*zbar)*dexp(6.d0*dlog(t8))*rho
         temp1=0.5*((cv*cv+ca*ca)+n*(cvp*cvp+cap*cap))*fgas
         temp2=0.5*((cv*cv-ca*ca)+n*(cvp*cvp-cap*cap))*ggas
         qgas=qgas*(temp1-temp2)
         qbrem=qgas
         return
       else
c
c strongly degenerate electrons (only valid for 1<gamma<160
c
         fliqy=0.d0
         fliqc=0.d0
         fliqo=0.d0
         gliqy=0.d0
         gliqc=0.d0
         gliqo=0.d0
c
         if (y.gt.0) call FGLIQY(t,rho,mue,y,fliqy,gliqy)
         if (c.gt.0) call FGLIQC(t,rho,mue,c,fliqc,gliqc)
         if (o.gt.0) call FGLIQO(t,rho,mue,1.d0-x-y-c,fliqo,gliqo)
c
         n=2
         temp=0.5738d0*dexp(6.d0*dlog(t8))*rho
         facf=0.5*((cv*cv+ca*ca)+n*(cvp*cvp+cap*cap))
         facg=0.5*((cv*cv-ca*ca)+n*(cvp*cvp-cap*cap))
         qbrem=temp*(y*1.d0*(facf*fliqy-facg*gliqy)
     1              +c*3.d0*(facf*fliqc-facg*gliqc)
     2              +o*4.d0*(facf*fliqo-facg*gliqo))
       endif
       return
       end
c
       subroutine PHOTO(flam,rho,t,mue,cv,cvp,ca,cap,qphot)
c From Itoh et al. 1996, results checked with tables from 
c Itoh et al.(1989 Ap.J. 339, 354),
       implicit double precision (a-h,o-z)
       dimension a(0:2),b(3)
       dimension c78(0:6,0:2),d78(1:5,0:2)
       dimension c89(0:6,0:2),d89(1:5,0:2)
       data b/6.290d-3,7.483d-3,3.061d-4/
       data c78/1.008d11,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,
     1   8.156d10,9.728d8,-3.806d9,-4.384d9,-5.774d9,-5.249d9,-5.153d9,
     2   1.067d11,-9.782d9,-7.193d9,-6.936d9,-6.893d9,-7.041d9,-7.193d9/ 
       data d78/0.d0,0.d0,0.d0,0.d0,0.d0,
     1   -1.879d10,-9.667d9,-5.602d9,-3.370d9,-1.825d9,
     2   -2.919d10,-1.185d10,-7.270d9,-4.222d9,-1.560d9/
       data c89/9.889d10,-4.524d8,-6.088d6,4.269d7,5.172d7,4.910d7,
     *         4.388d7,
     1   1.813d11,-7.556d9,-3.304d9,-1.031d9,-1.764d9,-1.851d9,-1.928d9,
     2   9.750d10,3.484e10,5.199d9,-1.695d9,-2.865d9,-3.395d9,-3.418d9/
       data d89/-1.135d8,1.256d8,5.149d7,3.436d7,1.005d7,
     1           1.652d9,-3.119d9,-1.839d9,-1.458d9,-8.956d8,
     2          -1.548d10,-9.338d9,-5.899d9,-3.035d9,-1.598d9/
c
       double precision mue
       cln=2.302585093d0
c       xi=(((rho/mue)/1e9)**(1./3.))/flam
       xi=(dexp(0.3333333d0*dlog(rho/mue/1d9)))/flam
       if (t.lt.1.d8) then
         tau=dlog10(t/1.d7)
         c=0.5654d0+tau
       elseif (t.gt.1.d8.and.t.lt.1.d9) then
         tau=dlog10(t/1.d8)
         c=1.5654d0
       endif
       do 10 i=0,2
         sum=0.d0
         do 20 j=1,5
           arg=5.235988d0*j*tau
           if (t.ge.1.d7.and.t.lt.1.d8) then
             term=c78(j,i)*cos(arg)+d78(j,i)*sin(arg)
           else
             term=c89(j,i)*cos(arg)+d89(j,i)*sin(arg)
           endif
           sum=sum+term
20       continue
         if (t.ge.1.d7.and.t.lt.1.d8) then
           a(i)=0.5d0*c78(0,i)+sum+0.5d0*c78(6,i)*cos(31.41593d0*tau)
         else
           a(i)=0.5d0*c89(0,i)+sum+0.5d0*c89(6,i)*cos(31.41593d0*tau)
         endif
10     continue
       call CALC(a,b,c,xi,flam,fphoto)
       n=2
       temp=1.875d8*flam+1.653d8*flam**2+8.499d8*flam**3
     1                                  -1.604d8*flam**4
       smq=0.666d0*dexp(-2.066d0*dlog(1.d0+2.045d0*flam))
     1                                           /(1.d0+rho/mue/temp)
       temp1=0.5*((cv**2+ca**2)+n*(cvp**2+cap**2))
       temp2=((cv**2-ca**2)+n*(cvp**2-cap**2))/
     1       ((cv**2+ca**2)+n*(cvp**2+cap**2))
       qphot=temp1*(1.d0-temp2*smq)*(rho/mue)*flam**5*fphoto
       return
       end
c
       subroutine PAIR(flam,rho,mue,cv,cvp,ca,cap,qpair)
c from Itoh et al. 1987, 1996; checked with Itoh et al. 1996, fig3
       implicit double precision (a-h,o-z)
       dimension a(0:2),b(3)
       data a/6.002d19,2.084d20,1.872d21/
       data b/9.383d-1,-4.141d-1,5.829d-2/
c
       double precision mue
       xi=(dexp(1.d0/3.d0*dlog(rho/mue/1e9)))/flam
       c=5.5924d0
       rflam=dexp(0.5*dlog(flam))
       flam2=flam*flam
       flam3=flam2*flam
       flam4=flam2*flam2
       flam6=flam2*flam4
       flam8=flam2*flam6
       call CALC(a,b,c,xi,flam,fpair)
       g=1.d0-13.04d0*flam2+133.5d0*flam4
     1               +1534d0*flam6+918.6d0*flam8
       gef=g*dexp(-2.d0/flam)*fpair
       n=2
       temp=1.d0/(10.748d0*flam2+0.3967d0*rflam+1.0050d0)
       smq=temp*(dexp(-0.3d0*dlog(1.d0+(rho/mue)
     1                /(7.692d7*flam3+9.715d6*rflam))))
       cvmm=(cv**2-ca**2)+n*(cvp**2-cap**2)
       cvpp=(cv**2+ca**2)+n*(cvp**2+cap**2)
       qpair=0.5d0*cvpp*(1.d0+(cvmm/cvpp)*smq)*gef
       return
       end
c
       subroutine CALC(a,b,c,xi,flam,f)
       implicit double precision (a-h,o-z)
       dimension a(0:2),b(3)
       top=(a(0)+a(1)*xi+a(2)*(xi**2))*dexp(-c*xi)	
       bottom=(xi**3)+(b(1)/flam)+(b(2)/(flam**2))+(b(3)/(flam**3))
       f=top/bottom
       return
       end
c
       subroutine RECOMB(flam,t,rho,z,qrec)
       implicit double precision (a-h,o-z)
       zeta=1.579d5*z*z/t
       temp=1.850d-4*(z**6)*((2.d0*z)**(-2))*rho*rho*flam*flam
       qrec=dexp(-zeta-(2.428d-5*(rho**(0.66667d0))/flam))*temp
c       qrec=dexp(-zeta-(2.428d-5*(dexp(2.d0/3.d0*rho))/flam))*temp
       return
       end
c
       subroutine FGLIQY(t,rho,mue,y,fliqy,gliqy)
c Itoh et al 1996
       implicit double precision (a-h,o-z)
c fitting coefficients for helum
       dimension ay(0:5),by(1:4),ey(0:5),fy(1:4)
       double precision cy,dy,gy,hy
       double precision iy(0:5),jy(1:4),ky,ly,py(0:5),qy(1:4),ry,sy
       double precision aly(0:3),bty(0:3)
       data ay/ 0.09037d0,-0.03009d0,-0.00564d0,-0.00544d0,
     1         -0.00290d0,-0.00224d0/
       data by/-0.02148d0,-0.00817d0,-0.00300d0,-0.00170d0/
       data cy/ 0.00671d0/
       data dy/ 0.28130d0/
       data ey/-0.02006d0, 0.01790d0,-0.00783d0,-0.00021d0,
     1          0.00024d0,-0.00014d0/
       data fy/ 0.00538d0,-0.00175d0,-0.00346d0,-0.00031d0/
       data gy/-0.02199d0/
       data hy/ 0.17300d0/
c
       data iy/ 0.00192d0,-0.00301d0,-0.00073d0, 0.00182d0,
     1          0.00037d0, 0.00116d0/
       data jy/ 0.01706d0,-0.00753d0, 0.00066d0,-0.00060d0/
       data ky/-0.01021d0/
       data ly/ 0.06417d0/
       data py/-0.01112d0, 0.00603d0,-0.00149d0, 0.00047d0,
     1          0.00040d0, 0.00028d0/
       data qy/ 0.00422d0,-0.00009d0,-0.00066d0,-0.00003d0/
       data ry/-0.00561d0/
       data sy/ 0.03522d0/
c
       data aly/-0.07980d0, 0.17057d0, 1.51980d0,-0.61058d0/
       data bty/-0.05881d0, 0.00165d0, 1.82700d0,-0.76993d0/
c
       double precision mue
c
       cln=2.302585093d0
       t8=t/1.d8
       rome=rho/mue
       temp=1.018d0*dexp(0.6666667*dlog(rome/1d6))
       tf=5.9302d9*(dexp(0.5*dlog(1.d0+temp))-1)
c
c strongly degenerate electrons (only valid for 1<gamma<160
c
       gamma=9.100d-1/t8*dexp(0.3333d0*dlog(y*rho/1.d6/4.d0))
       u=6.283185d-1*(dlog10(rho)-3.d0)
       sumacos=0.d0
       sumecos=0.d0
       sumicos=0.d0
       sumpcos=0.d0
       do 10 m=1,5
         arg=1.d0*m*u
         sumacos=sumacos+ay(m)*cos(arg)
         sumecos=sumecos+ey(m)*cos(arg)
         sumicos=sumicos+iy(m)*cos(arg)
         sumpcos=sumpcos+py(m)*cos(arg)
10     continue
       sumbsin=0.d0
       sumfsin=0.d0
       sumjsin=0.d0
       sumqsin=0.d0
       do 20 m=1,4
         arg=1.d0*m*u
         sumbsin=sumbsin+by(m)*sin(arg)
         sumfsin=sumfsin+fy(m)*sin(arg)
         sumjsin=sumjsin+jy(m)*sin(arg)
         sumqsin=sumqsin+qy(m)*sin(arg)
20     continue
       v=0.d0
       w=0.d0
       do 30 m=0,3
         tempg=dexp(-m/3.d0*dlog(gamma))
         v=v+aly(m)*tempg
         w=w+bty(m)*tempg
30     continue
       fu1=ay(0)/2.d0+sumacos+sumbsin+cy*u+dy
       fu160=ey(0)/2.d0+sumecos+sumfsin+gy*u+hy
       gu1=iy(0)/2.d0+sumicos+sumjsin+ky*u+ly
       gu160=py(0)/2.d0+sumpcos+sumqsin+ry*u+sy
       fliqy=v*fu1+(1.d0-v)*fu160
       gliqy=w*gu1+(1.d0-w)*gu160
c
       return
       end
c
       subroutine FGLIQC(t,rho,mue,c,fliqc,gliqc)
c Itoh et al 1996
       implicit double precision (a-h,o-z)
c fitting coefficients for carbon
       dimension ac(0:5),bc(1:4),ec(0:5),fc(1:4)
       double precision cc,dc,gc,hc
       double precision ic(0:5),jc(1:4),kc,lc,pc(0:5),qc(1:4),rc,sc
       double precision alc(0:3),btc(0:3)
       data ac/ 0.17946d0,-0.05821d0,-0.01089d0,-0.01147d0,
     1         -0.00656d0,-0.00519d0/
       data bc/-0.04969d0,-0.01584d0,-0.00504d0,-0.00281d0/
       data cc/ 0.00945d0/
       data dc/ 0.34529d0/
       data ec/ 0.06781d0,-0.00944d0,-0.01289d0,-0.00589d0,
     1         -0.00404d0,-0.00330d0/
       data fc/-0.02213d0,-0.01136d0,-0.00467d0,-0.00131d0/
       data gc/-0.02342d0/
       data hc/ 0.24819d0/
c
       data ic/ 0.00766d0,-0.00710d0,-0.00028d0, 0.00232d0,
     1          0.00044d0, 0.00158d0/
       data jc/ 0.02300d0,-0.01078d0, 0.00118d0,-0.00089d0/
       data kc/-0.01259d0/
       data lc/ 0.07917d0/
       data pc/-0.00769d0, 0.00356d0,-0.00184d0, 0.00146d0,
     1          0.00031d0, 0.00069d0/
       data qc/ 0.01052d0,-0.00354d0,-0.00014d0,-0.00018d0/
       data rc/-0.00829d0/
       data sc/ 0.05211d0/
c
       data alc/-0.05483d0,-0.01946d0, 1.86310d0,-0.78873d0/
       data btc/-0.06711d0, 0.06859d0, 1.74360d0,-0.74498d0/
c
       double precision mue
c
       cln=2.302585093d0
       t8=t/1.d8
       rome=rho/mue
       temp=1.018d0*dexp(0.6666667*dlog(rome/1d6))
       tf=5.9302d9*(dexp(0.5*dlog(1.d0+temp))-1)
c
c strongly degenerate electrons (only valid for 1<gamma<160
c
       gamma=8.190d0/t8*dexp(0.3333d0*dlog(c*rho/1.d6/12.d0))
       u=6.283185d-1*(dlog10(rho)-3.d0)
       sumacos=0.d0
       sumecos=0.d0
       sumicos=0.d0
       sumpcos=0.d0
       do 10 m=1,5
         sumacos=sumacos+ac(m)*cos(m*u)
         sumecos=sumecos+ec(m)*cos(m*u)
         sumicos=sumicos+ic(m)*cos(m*u)
         sumpcos=sumpcos+pc(m)*cos(m*u)
10     continue
       sumbsin=0.d0
       sumfsin=0.d0
       sumjsin=0.d0
       sumqsin=0.d0
       do 20 m=1,4
         sumbsin=sumbsin+bc(m)*sin(m*u)
         sumfsin=sumfsin+fc(m)*sin(m*u)
         sumjsin=sumjsin+jc(m)*sin(m*u)
         sumqsin=sumqsin+qc(m)*sin(m*u)
20     continue
       v=0.d0
       w=0.d0
       do 30 m=0,3
         tempg=dexp(-m/3.d0*dlog(gamma))
         v=v+alc(m)*tempg
         w=w+btc(m)*tempg
30     continue
       fu1=ac(0)/2.d0+sumacos+sumbsin+cc*u+dc
       fu160=ec(0)/2.d0+sumecos+sumfsin+gc*u+hc
       gu1=ic(0)/2.d0+sumicos+sumjsin+kc*u+lc
       gu160=pc(0)/2.d0+sumpcos+sumqsin+rc*u+sc
       fliqc=v*fu1+(1.d0-v)*fu160
       gliqc=w*gu1+(1.d0-w)*gu160
c
       return
       end
c
c
       subroutine FGLIQO(t,rho,mue,o,fliqo,gliqo)
c Itoh et al 1996
       implicit double precision (a-h,o-z)
c fitting coefficients for oxygen
       dimension ao(0:5),bo(1:4),eo(0:5),fo(1:4)
       double precision co,do,go,ho
       double precision io(0:5),jo(1:4),ko,lo,po(0:5),qo(1:4),ro,so
       double precision alo(0:3),bto(0:3)
c
       data ao/ 0.20933d0,-0.06740d0,-0.01293d0,-0.01352d0,
     1         -0.00776d0,-0.00613d0/
       data bo/-0.05950d0,-0.01837d0,-0.00567d0,-0.00310d0/
       data co/ 0.00952d0/
       data do/ 0.36029d0/
       data eo/ 0.09304d0,-0.01656d0,-0.01489d0,-0.00778d0,
     1         -0.00520d0,-0.00418d0/
       data fo/-0.03076d0,-0.00522d0,-0.00161d0,-0.02513d0/
       data go/-0.02513d0/
       data ho/ 0.27480d0/
c
       data io/ 0.00951d0,-0.00838d0,-0.00011d0, 0.00244d0,
     1          0.00046d0, 0.00168d0/
       data jo/ 0.02455d0,-0.01167d0, 0.00132d0,-0.00097d0/
       data ko/-0.01314d0/
       data lo/ 0.08263d0/
       data po/-0.00700d0, 0.00295d0,-0.00184d0, 0.00166d0,
     1          0.00032d0, 0.00082d0/
       data qo/ 0.01231d0,-0.00445d0, 0.00002d0,-0.00026d0/
       data ro/-0.00921d0/
       data so/ 0.05786d0/
c
       data alo/-0.06597d0, 0.06048d0, 1.74860d0,-0.74308d0/
       data bto/-0.07123d0, 0.10865d0, 1.70150d0,-0.73653d0/
c
       double precision mue
c
       cln=2.302585093d0
       t8=t/1.d8
       rome=rho/mue
       temp=1.018d0*dexp(0.6666667*dlog(rome/1d6))
       tf=5.9302d9*(dexp(0.5*dlog(1.d0+temp))-1)
c
c strongly degenerate electrons (only valid for 1<gamma<160
c
       gamma=1.456d1/t8*dexp(0.3333d0*dlog(o*rho/1.d6/16.d0))
       u=6.283185d-1*(dlog10(rho)-3.d0)
       sumacos=0.d0
       sumecos=0.d0
       sumicos=0.d0
       sumpcos=0.d0
       do 10 m=1,5
         sumacos=sumacos+ao(m)*cos(m*u)
         sumecos=sumecos+eo(m)*cos(m*u)
         sumicos=sumicos+io(m)*cos(m*u)
         sumpcos=sumpcos+po(m)*cos(m*u)
10     continue
       sumbsin=0.d0
       sumfsin=0.d0
       sumjsin=0.d0
       sumqsin=0.d0
       do 20 m=1,4
         sumbsin=sumbsin+bo(m)*sin(m*u)
         sumfsin=sumfsin+fo(m)*sin(m*u)
         sumjsin=sumjsin+jo(m)*sin(m*u)
         sumqsin=sumqsin+qo(m)*sin(m*u)
20     continue
       v=0.d0
       w=0.d0
       do 30 m=0,3
         tempg=dexp(-m/3.d0*dlog(gamma))
         v=v+alo(m)*tempg
         w=w+bto(m)*tempg
30     continue
       fu1=ao(0)/2.d0+sumacos+sumbsin+co*u+do
       fu160=eo(0)/2.d0+sumecos+sumfsin+go*u+ho
       gu1=io(0)/2.d0+sumicos+sumjsin+ko*u+lo
       gu160=po(0)/2.d0+sumpcos+sumqsin+ro*u+so
       fliqo=v*fu1+(1.d0-v)*fu160
       gliqo=w*gu1+(1.d0-w)*gu160
c
       return
       end
c