subroutine ww(inpts,inrand,wmi,rsi,s0i,thc,tc,alszi,sig,esig)
      implicit real*8(a-h,o-z)
      character*1 oqcd,oral,opeak
*
      parameter (nin=98,nout=99,ndimmx=7,liw=ndimmx+2,
     #           lw=ndimmx*(ndimmx-1)/2+12*ndimmx)
*
      common/wtpk/opeak
      common/wtifz/ifz(10),ifzzf,ik
      common/wtcut/opcut,omcut,tcs,s0
      common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis,
     #              pih,tfact
      common/wtapar/s,rwm,rwm2,rwg,rwmg,swg,swgs,opswgs,s0w,
     #              g2,sth2,isf
*
      dimension vk(ndimmx)
*
      external d01gcf,wtoregion,wwcc20
*
      do i=1,10
          ifz(i)= 0
      enddo
      ifzzf= 0
      ik= 0
      iqcd= 1
*      oqcd= 'y'
      oqcd= 'n'
      oral= 'r'
      iosf= 1
      ios= 2
      opeak= 'y'
*
      cfct= 0.38937966d9
      alphai= 137.0359895d0
      em= 0.51099906d-3
      pi= 3.141592653589793238462643d0
      pis= pi*pi
      pih= 0.5d0*pi
      ge= 0.5772156649d0
      gf= 8.2476172696d-6
      zm= 91.1866d0
      rs= rsi
      s= rs*rs
      cut= cos(thc/180.d0*pi)
      opcut= 1.d0+cut
      omcut= 1.d0-cut
      tcs= tc*tc
      s0= s0i
*
      wm= wmi
      alsz= alszi
*
      nf= 5
      if(iqcd.eq.2) then
         x1= 0.001d0
         x2= 10.0d0
         xacc= 1.d-12
         nf= 5
         qcdl= wtoqcdlam(nf,als,wm,x1,x2,xacc)
      else
         qcdl= 0.d0
      endif
      als= wtoralphas(zm,wm,alsz,nf)
*
      qgfopi= 0.25d0*gf/pi
      wm2= wm*wm
      wmc= wm2*wm
      zm2= zm*zm
      wg0= 1.5d0*gf*wmc/pi
      if(oqcd.eq.'n') then
         wg= wg0
      else
         wg= wg0*(1.d0+2.d0/3.d0*als/pi)
      endif
      wg2= wg*wg
*
      rwm= wm/rs
      rwm2= wm2/s
      rwg= wg/rs
      swg= wg/wm
      swgs= swg*swg
      rwmg= rwm*rwg
      opswgs= 1.d0+swgs
      s0w= rwm2/opswgs
*
      alwi= 128.07d0
*
*-----e.m. quantities are initialized
*
      em2= em*em
      alpha= 1.d0/alphai
      alw= 1.d0/alwi
      aexp= alpha/pi
      rln= log(s/em2)
      beta= 2.d0*aexp*(rln-1.d0)
      isf= iosf
      if(iosf.eq.1) then
         eta= beta
         eob= 1.d0
      else if(iosf.gt.1) then
         eta= 2.d0*aexp*rln
         eob= rln/(rln-1.d0)
      endif
      feta= eta/16.d0
      hbet= beta/2.d0
      heta= eta/2.d0
      hbeti= 2.d0/beta
      omhb= 1.d0-hbet
      ophb= 1.d0+hbet
      if(iosf.lt.3) then 
         arge= hbet*(3.d0/4.d0-ge)
      else if(iosf.eq.3) then
         arge= -hbet*ge+3.d0/4.d0*heta
      endif
      ifg= 1
      gamb= s14aaf(ophb,ifg)
      if(ifg.ne.0) then
         stop
      else 
         d0gl= exp(arge)/gamb
      endif
*
      if(ios.eq.1) then
         sth2= 0.5d0*pi*alw/gf/wm2
         g2= 4.d0*pi*alw/sth2
      else if(ios.eq.2) then
         sth2= 1.d0-wm2/zm2
         g2= 8.d0*gf*wm2
      endif
      g4= g2*g2
      g8= g4*g4
      sth4= sth2*sth2
      tfact= g8*sth4/2.d0/(2.d0*pi)**8*pi**3/2.d0*3.d0*cfct
*
      if(inpts.le.6) then
         jnpts= inpts
      else if(inpts.eq.7) then
         jnpts= 2011*101
         vk(1)= 0.1d+1
         vk(2)= 0.72725d+5
         vk(3)= 0.118296d+6
         vk(4)= 0.107084d+6
         vk(5)= 0.1938d+4
         vk(6)= 0.185127d+6
         vk(7)= 0.14844d+6
      else if(inpts.eq.8) then
         jnpts= 10193*101
         vk(1)= 0.1d+1
         vk(2)= 0.1019431d+7
         vk(3)= 0.35353d+6
         vk(4)= 0.708948d+6
         vk(5)= 0.951714d+6
         vk(6)= 0.197618d+6
         vk(7)= 0.54816d+6
      else if(inpts.eq.9) then
         jnpts= 22807*151
         vk(1)= 0.1d+01
         vk(2)= 0.1619459d+7
         vk(3)= 0.226133d+7
         vk(4)= 0.256381d+7
         vk(5)= 0.230245d+7
         vk(6)= 0.855081d+6
         vk(7)= 0.609193d+6
      else if(inpts.eq.10) then
         jnpts= 32771*181
         vk(1)= 0.1d+01
         vk(2)= 0.4960373d+7
         vk(3)= 0.4851623d+7
         vk(4)= 0.3262017d+7
         vk(5)= 0.722217d+6
         vk(6)= 0.4644124d+7
         vk(7)= 0.3212165d+7
      else if(inpts.eq.11) then
         jnpts= 32771*379
         vk(1)= 0.1d+01
         vk(2)= 0.9646626d+07
         vk(3)= 0.8128723d+07
         vk(4)= 0.9521278d+07
         vk(5)= 0.4720279d+07
         vk(6)= 0.7036407d+07
         vk(7)= 0.1868554d+07
      endif
*
      ndim= 7
*      ndim= 5
      npts= jnpts
      itrans= 0
      nrand= inrand
      ifail= 0
      call d01gcf(ndim,wwcc20,wtoregion,npts,vk,nrand,
     #            itrans,sig,esig,ifail)
*
      ifzt= 0
      do i=1,10
          ifzt= ifzt+ifz(i)
      enddo
*      print 1,ik
*      print 2,1.d2*(1.d0-1.d0*ifzt/ik)
*      print 3,ifzzf
*      do i=1,10
*          print 5,ifz(i)
*      enddo
*      print 4
    1 format(/' # of calls      = ',i10)
    2 format(/' 1-nc/int (%)    = ',f10.3)
    3 format(/' # of neg M^2    = ',i10)
    4 format(/' End of Call DUMP -----------------------------',/)
    5 format(1x,i10)
      return
      end
*
*------------------------------------------------------------------
*
      real*8 function wwcc20(ndim,x)
      implicit real*8(a-h,o-z)
      character*1 opeak
*
      common/wtpk/opeak
      common/wtifz/ifz(10),ifzzf,ik
      common/wtcut/opcut,omcut,tcs,s0
      common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis,
     #              pih,tfact
      common/wtapar/s,rwm,rwm2,rwg,rwmg,swg,swgs,opswgs,s0w,
     #              g2,sth2,isf
*
      dimension x(ndim),zrt1(2),zrt2(2)
*
      external c02ajf,s09aaf
*
*-----order of integration is: uv,vv; y; x2,x1,tau,z
*
      ik= ik+1
      em= 0.51099906d-3
      em2= em*em
      rem= em/sqrt(s)
*
      if(ndim.eq.5) then
       yx= x(1)
       x2x= x(2)
       x1x= x(3)
       taux= x(4)
       zx= x(5)
      else
       uvx= x(1)
       vvx= x(2)
       yx= x(3)
       x2x= x(4)
       x1x= x(5)
       taux= x(6)
       zx= x(7)
      endif
*
      iz= 1
*
      if(ndim.eq.5) then
         ueps= 0.d0
         uv= 1.d0
         uvs= uv*uv
         ujc= 1.d0
         veps= 0.d0
         vv= 1.d0
         vjc= 1.d0
      else
       if(uvx.gt.1.d0) then
          iz= 0
          ifz(1)= ifz(1)+1
          go to 1
       endif
       vvl= (sqrt(s0)+0.5d0*rem)*(sqrt(s0)+0.5d0*rem)
       omul= 1.d0-vvl
       ueps= omul*(1.d0-uvx)**hbeti
       uv= 1.d0-ueps
       ujc= omul**hbet
       uvs= uv*uv
*
       uvl= uv-vvl
       veps= uvl*(1.d0-vvx)**hbeti
       vv= uv-veps
       vjc0= 1.d0-vvl/uv
       if(vjc0.le.0.d0) then
          iz= 0
          ifz(2)= ifz(2)+1
          go to 1
       else
          vjc= vjc0**hbet
       endif
      endif
*
      if(vv.lt.0.d0) then
         iz= 0
         ifz(2)= ifz(2)+1
         go to 1
      endif
      svv= sqrt(vv)
      vvs= vv*vv
      xm= uv
      xp= vv/uv  
      xmop= xm/xp
      sh= vv*s
      ssh= sqrt(sh)
*
      opxp= 1.d0+xp
      opxm= 1.d0+xm
      omxp= veps/uv
      omxm= ueps
      if(isf.eq.0) then
         stfp= 1.d0
         stfm= 1.d0
      else if(isf.gt.0) then
         if(omxp.eq.0) then
            stfp= d0gl
         else
            rcpx= 0.25d0*opxp*opxp
            rcpy= xp
            iflp= 1
            rclp= s21baf(rcpx,rcpy,iflp)
            stfp= d0gl+eob*omxp**omhb*(-0.5d0*opxp+
     #            feta*(-4.d0*opxp*log(omxp)+
     #            3.d0*opxp*log(xp)+4.d0*rclp-5.d0-xp))
         endif
         if(omxm.eq.0) then
            stfm= d0gl
         else
            rcmx= 0.25d0*opxm*opxm
            rcmy= xm
            iflm= 1
            rclm= s21baf(rcmx,rcmy,iflm)
            stfm= d0gl+eob*omxm**omhb*(-0.5d0*opxm+
     #            feta*(-4.d0*opxm*log(omxm)+
     #            3.d0*opxm*log(xm)+4.d0*rclm-5.d0-xm))
         endif
      endif
*
      stf= stfp*stfm
*
*-----y is initialized
*
      yl= s0/vv
      yu= 1.d0-em/ssh
      yjc= yu-yl
      y= yjc*yx+yl
      ys= y*y
      yvv= y*vv
      omy= (1.d0-s0/vv)*(1.d0-yx)+em/ssh*yx
      omys= omy*omy
      omyc= omys*omy
*
*-----x2 is initialized
*
      x2l= s0/yvv
      x2u= 1.d0
      vwmg= rwmg*yvv
      zma= y*vv*x2l
      zmb= y*vv*x2u
      zmas= zma-rwm2
      zmbs= zmb-rwm2
      atma= zmas/rwmg
      atmb= zmbs/rwmg
*
      if(opeak.eq.'y') then
       if(atma.gt.1.d0.and.atmb.gt.1.d0) then
          atma= 1.d0/atma
          atma= atan(atma)
          zmat= (pih-atma)/vwmg
          atmb= 1.d0/atmb
          atmb= atan(atmb)
          zmbt= (pih-atmb)/vwmg
          smjc0= (-atmb+atma)/vwmg
       else if(atma.gt.1.d0.and.atmb.lt.-1.d0) then
          atma= 1.d0/atma
          atma= atan(atma)
          zmat= (pih-atma)/vwmg
          atmb= -1.d0/atmb
          atmb= atan(atmb)
          zmbt= (-pih+atmb)/vwmg
          smjc0= (-pi+atmb+atma)/vwmg
       else if(atma.gt.1.d0.and.abs(atmb).lt.1.d0) then
          atma= 1.d0/atma
          atma= atan(atma)
          zmat= (pih-atma)/vwmg
          atmb= atan(atmb)
          zmbt= atmb/vwmg
          smjc0= (-pih+atmb+atma)/vwmg
       else if(atma.lt.-1.d0.and.atmb.gt.1.d0) then
          atma= -1.d0/atma
          atma= atan(atma)
          zmat= (-pih+atma)/vwmg
          atmb= 1.d0/atmb
          atmb= atan(atmb)
          zmbt= (pih-atmb)/vwmg
          smjc0= (pi-atmb-atma)/vwmg
       else if(atma.lt.-1.d0.and.atmb.lt.-1.d0) then
          atma= -1.d0/atma
          atma= atan(atma)
          zmat= (-pih+atma)/vwmg
          atmb= -1.d0/atmb
          atmb= atan(atmb)
          zmbt= (-pih+atmb)/vwmg
          smjc0= (atmb-atma)/vwmg
       else if(atma.lt.-1.d0.and.abs(atmb).lt.1.d0) then
          atma= -1.d0/atma
          atma= atan(atma)
          zmat= (-pih+atma)/vwmg
          atmb= atan(atmb)
          zmbt= atmb/vwmg
          smjc0= (pih+atmb-atma)/vwmg
       else if(abs(atma).lt.1.d0.and.atmb.gt.1.d0) then
          atma= atan(atma)
          zmat= atma/vwmg
          atmb= 1.d0/atmb
          atmb= atan(atmb)
          zmbt= (pih-atmb)/vwmg
          smjc0= (pih-atmb-atma)/vwmg
       else if(abs(atma).lt.1.d0.and.atmb.lt.-1.d0) then
          atma= atan(atma)
          zmat= atma/vwmg
          atmb= -1.d0/atmb
          atmb= atan(atmb)
          zmbt= (-pih+atmb)/vwmg
          smjc0= (-pih+atmb-atma)/vwmg
       else if(abs(atma).lt.1.d0.and.abs(atmb).lt.1.d0) then
          atma= atan(atma)
          zmat= atma/vwmg
          atmb= atan(atmb)
          zmbt= atmb/vwmg
          smjc0= (atmb-atma)/vwmg
       endif
*
       zmv= smjc0*x2x+zmat
       iftn= 1
       atnm= vwmg*zmv
       x2= (rwm2+rwmg*s07aaf(atnm,iftn))/(y*vv)
       x2jc= y*vv*smjc0
       pmjc= 1.d0
*
      else if(opeak.eq.'n') then         
       smjc0= zmb-zma
       x2= (smjc0*x2x+zma)/(y*vv)
       pmjc= 1.d0/((y*vv*x2-rwm2)**2+rwmgs)
       x2jc= smjc0
      endif
      omx2= 1.d0-x2
*
*-----x1 is initialized
*
      x1l= x2
      x1u= 1.d0
      if(x1l.gt.x1u) then
         iz= 0
         ifz(3)= ifz(3)+1
         go to 1
      endif
      x1jc= x1u-x1l
      x1= x1jc*x1x+x1l
      omx1= omx2*(1.d0-x1x)
      omx1s= omx1*omx1
      x12= x1*x2
      x21s= x2*x1s
*
*-----tau is initialized
*
      taum= -opcut*xm*y*x1/(opcut*xm*y+omcut*xp)
      taup= -omcut*xm*y*x1/(omcut*xm*y+opcut*xp)
      taumm= -x1
      taupp= 0.d0
      taul= dmax1(taum,taumm)
      tauu= dmin1(taup,taupp)
      if(taul.gt.tauu) then
         iz= 0
         ifz(4)= ifz(4)+1
         go to 1
      endif
      taujc= tauu-taul
      tau= taujc*taux+taul
      taus= tau*tau
      optau= 1.d0+tau
*
*-----z is initialized
*
      if((1.d0-x2).eq.0.d0.or.x1.eq.0.d0) then
         iz= 0
         ifz(5)= ifz(5)+1
         go to 1
      else
         ct= 1.d0+2.d0*tau/x1
         cp= 1.d0+2.d0*(x2-x1)/x1/(1.d0-x2)
         if(ct.lt.-1.d0.or.ct.gt.1.d0) then
            iz= 0
            ifz(5)= ifz(5)+1
            go to 1
         else
            st= sqrt(1.d0-ct*ct)
         endif
         if(cp.lt.-1.d0.or.cp.gt.1.d0) then
            iz= 0
            ifz(5)= ifz(5)+1
            go to 1
         else
            sp= sqrt(1.d0-cp*cp)
         endif
      endif
      zca= -4.d0
      zcb= 4.d0*(ct*(2.d0*x2-x1-x12)/x1-1.d0+x2)
      zcc= -((1.d0-x2)*ct-(2.d0*x2-x1-x12)/x1)**2
      ifvv= 1
      call c02ajf(zca,zcb,zcc,zrt1,zrt2,ifvv)
      if(ifvv.ne.0.or.zrt1(2).ne.0.d0) then
         iz= 0
         ifz(6)= ifz(6)+1
         go to 1
      endif
      zm= dmin1(zrt1(1),zrt2(1))
      zp= dmax1(zrt1(1),zrt2(1))
*
      azmm= -(opcut*xm*y*(tau-x2+x1)+omcut*xp*optau)/
     #       (opcut*xm*y+omcut*xp)
      azpp= -(omcut*xm*y*(tau-x2+x1)+opcut*xp*optau)/
     #       (omcut*xm*y+opcut*xp)
      zmm= dmin1(azmm,azpp)
      zpp= dmax1(azmm,azpp)
      zmmm= x2-1.d0
      zppp= 0.d0
      zpppp= x2-x1-tau
*
      za= dmax1(zm,zmm,zmmm)
      zb= dmin1(zp,zpp,zppp,zpppp)
*
      if(za.gt.zb) then
         iz= 0
         ifz(7)= ifz(7)+1
         go to 1
      endif
*
      zcb= 0.5d0*zcb
      sqz= sqrt(zcb*zcb+4.d0*zcc)
      if(za.eq.zm) then
         ta= -pi/4.d0
      else
         ifas= 1
         aa= (zcb-4.d0*za)/sqz
         ta= -0.5d0*s09aaf(aa,ifas)
         if(ifas.ne.0) then
          iz= 0
          ifz(8)= ifz(8)+1
          go to 1
         endif
      endif
      if(zb.eq.zp) then
         tb= pi/4.d0
      else
         ifas= 1
         bb= (zcb-4.d0*zb)/sqz
         tb= -0.5d0*s09aaf(bb,ifas)
         if(ifas.ne.0) then
          iz= 0
          ifz(8)= ifz(8)+1
          go to 1
         endif
      endif
      tl= dmin1(ta,tb)
      tu= dmax1(ta,tb)
      arg= (tu-tl)*zx+tl
      z= 0.25d0*(zcb+sqz*sin(2.d0*arg))
      zs= z*z
      zjc= tu-tl
*
      if(st.ne.0.d0.and.sp.ne.0.d0) then
         cphi= (2.d0*z/omx2+1.d0-ct*cp)/st/sp
         if(cphi.lt.-1.d0.or.cphi.gt.1.d0) then
            iz= 0
            ifz(9)= ifz(9)+1
            go to 1
         endif
      endif
*
      psjc= 0.25d0/x1
*
      qyjc= 1.d0-2.d0*em2/s*(xp+xm)/vv
*
      qh0s= em2*y/(1.d0-y)*(xm*y+1.d0-xm)
      qhcs= qh0s+0.25d0*xm*xm*(1.d0-y)*tcs*s
      by= 1.d0-xm*(1.d0-y)
      bk= 0.25d0*xm*xm*tcs
      psf= -omyc/y/by*bk*s/(omys*bk*s+y*by*em2)*
     #     (1.d0+xm+2.d0*(1.d0-xm)/y)+
     #     (1.d0-2.d0/y+2.d0/(y*y))*log(1.d0+omys*bk*s/
     #     (em2*y*by))
*
      upsi= 1.d0+tau+z
      pwsm2= (x2*yvv-rwm2)*(x2*yvv-rwm2)+rwmg*rwmg
      pwsr= rwm2-x2*yvv
      x2t= x2-z-1.d0
      pwtm2= (x2t*yvv-rwm2)*(x2t*yvv-rwm2)+rwmg*rwmg
      pwtr= rwm2-x2t*yvv
      rwmgs= rwmg*rwmg
      x1s= x1*x1
      x1c= x1s*x1
      x21s= x2*x1s
      x2s= x2*x2
      yc= ys*y
*
      tjac0= ujc*vjc*yjc*x2jc*x1jc*taujc*zjc*psjc/qyjc*pmjc
      tjac= tjac0*stf*psf/sh
*
      w_11= tau*(1.d0+z*x1+2.d0*z-2.d0*x12+x1+zs)
      w_11= w_11+taus*(1.d0+z-x2)
      w_11= w_11+z*x1+z*x1s+2.d0*x12+x1*zs-2.d0*x21s-x2+x1s
*
      w_22= y*pwtm2*tau*(1.d0-x1)
*
      w_33= 4.d0/9.d0*pwsm2*y/tau*(1.d0-x1)
*
      w_44= tau/upsi*(-1.d0-x1)
      w_44= w_44+tau
      w_44= w_44-1/upsi*taus
      w_44= w_44-1.d0/upsi*x1
      w_44= w_44+x1
      w_44= 1.d0/9.d0*w_44
*
      w_12= pwtr*tau*(-2.d0+z*x1-2.d0*z+2.d0*x12)
      w_12= w_12+pwtr*taus*(-1.d0+x2)
      w_12= w_12+pwtr*(-z*x1-2.d0*x12+2.d0*x21s+x2-x1s)
*
      w_13= pwsr/tau*(2./3.*z*x1-4./3.*z*x1s-4./3.*x12+2.
     #    /3.*x1+4./3.*x21s+2./3.*x2-4./3.*x1s)
      w_13= w_13+pwsr*tau*(-2./3.*x1+2./3.*x2)
      w_13= w_13+pwsr*(-2./3.-2./3.*z*x1+4./3.*x12-2./3.*x1s)
*
      w_14= pwsr*tau/upsi*(-1./3.-2./3.*x12+x1-2./3.*x1s)
      w_14= w_14+pwsr*tau*(1./3.+2./3.*z)
      w_14= w_14+pwsr/upsi*taus*(1./3.-1./3.*x1-1./3.*x2)
      w_14= w_14+pwsr/upsi*(2./3.*x12-2./3.*x21s-1./3.*x2)
      w_14= w_14+pwsr*(2./3.*z*x1+1./3.*x1s)
*
      w_23= pwsr*pwtr/tau*(-4./3.*x1+8./3.*x1s-4./3.*x1c)
      w_23= w_23+pwsr*pwtr*(4./3.*x1-4./3.*x1s)
      w_23= w_23+1/tau*rwmgs*(-4./3.*x1+8./3.*x1s-4./3.*x1c)
      w_23= w_23+rwmgs*(4./3.*x1-4./3.*x1s)
*
      w_24= pwsr*pwtr*tau/upsi*(1./3.+2./3.*x12+1./3.*x1
     #    -2./3.*x1s)
      w_24= w_24+pwsr*pwtr*tau*(-2./3.+1./3.*x1)
      w_24= w_24+pwsr*pwtr/upsi*taus*(1./3.-1./3.*x1+1./
     #    3.*x2)
      w_24= w_24+pwsr*pwtr/upsi*(-2./3.*x12+2./3.*x2*
     #    x1s+1./3.*x2+2./3.*x1s-2./3.*x1c)
      w_24= w_24+pwsr*pwtr*(-1./3.*x1)
      w_24= w_24+tau/upsi*rwmgs*(1./3.+2./3.*x12+1./3.
     #    *x1-2./3.*x1s)
      w_24= w_24+tau*rwmgs*(-2./3.+1./3.*x1)
      w_24= w_24+1/upsi*rwmgs*taus*(1./3.-1./3.*x1+1./3.*x2)
      w_24= w_24+1/upsi*rwmgs*(-2./3.*x12+2./3.*x21s
     #    +1./3.*x2+2./3.*x1s-2./3.*x1c)
      w_24= w_24+rwmgs*(-1./3.*x1)
*
      w_34= 1/tau/upsi*(-4./9.*x12+4./9.*x21s+2./9.*x2)
      w_34= w_34+1/tau*(2./9.*x1-4./9.*x1s)
      w_34= w_34+tau/upsi*(-2./9.+2./9.*x1+2./9.*x2)
      w_34= w_34+1/upsi*(2./9.+4./9.*x12-2./3.*x1+4./9.*x1s)
      w_34= w_34-2./9.*x1
*
      w_11= w_11*yc
      w_44= w_44*y*pwsm2
      w_12= w_12*ys
      w_13= w_13*ys
      w_14= w_14*ys
      w_23= w_23*y
      w_24= w_24*y
      w_34= w_34*y*pwsm2
*
*      www= -0.25d0*(w_11+w_22+w_33+w_44+w_12+w_13+w_14+
*     #              w_23+w_24+w_34)/pwtm2
*      www= -0.5d0*(w_11+w_22+w_33+w_44+w_12+w_13+w_14+
*     #              w_23+w_24+w_34)/pwtm2
      www= -0.5d0*(w_11+w_22+w_12)/pwtm2
*
      if((www*tjac*tfact).lt.0.d0) then
       iz= 0
       ifz(10)= ifz(10)+1
       ifzzf= ifzzf+1
      endif
*
    1 if(iz.eq.0) then
         wwcc20= 0.d0
      else
         wwcc20= www*tjac*tfact
      endif
*
      return
      end      
*
*-----wtoregion-----------------------------------------------------
*
      subroutine wtoregion(n,x,j,a,b)
      implicit real*8(a-h,o-z)
*
      dimension x(n)
*
      if(j.gt.1) then
         do i=1,j-1
             if(x(i).lt.0.d0.or.x(i).gt.1.d0) then
                stop
             endif
         enddo
      endif
*
      a= 0.d0
      b= 1.d0
*
      return
      end
*
*-----wtoqcdlam-----------------------------------------------------------------
*-----computes lamba^nf_msbar from alpha_s(rs0)
*
      real*8 function wtoqcdlam(nf,als,rs,x1,x2,xacc)
      implicit real*8(a-h,o-z)
*
      parameter (jmax=50,nout=21)
*
      fmid= wtoqcdscale(nf,als,rs,x2)
      f= wtoqcdscale(nf,als,rs,x1)
      if (f*fmid.ge.0.d0) then
          write(nout,fmt=*) 'root must be bracketed for bisection'
          print*,'root must be bracketed for bisection'
          write(nout,fmt=1) als
          print 1,als
    1     format(/' error detected by wtoqcdlam ',/
     #            ' current value of alpha_s  = ',e20.5)
          stop
      endif
      if (f.lt.0.d0) then
          wtoqcdlam= x1
          dx= x2-x1
      else
          wtoqcdlam= x2
          dx= x1-x2
      endif
      do j=1,jmax
          dx= dx*0.5d0
          xmid= wtoqcdlam+dx
          fmid= wtoqcdscale(nf,als,rs,xmid)
          if (fmid.le.0.d0) wtoqcdlam= xmid
          if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return
      enddo
      pause 'too many bisections'
      end       
*
*-----wtoqcdscale-------------------------------------------------------------
*-----computes lamba^nf_msbar from alpha_s(rs0)
*
      real*8 function wtoqcdscale(nf,als,rs,x)
      implicit real*8(a-h,o-z)
*
      common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis,
     #              pih,tfact
*
      qcdb0= 11.d0-2.d0/3.d0*nf
      qcdb1= 102.d0-38.d0/3.d0*nf
      qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nf+325.d0/27.d0*nf*nf)
      qcda= 2.d0*log(rs/x)
      wtoqcdscale= als-(4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*
     #            log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-
     #            0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)))
*
      return
      end
*
*-----wtoralphas--------------------------------------------------------------
*-----computes alpha_s(rs) for nf flavors given als= alpha_s(wm)
*
      real*8 function wtoralphas(rs0,rs,als,nf)
      implicit real*8(a-h,o-z)
*
      common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta,pi,pis,
     #              pih,tfact
*
      dimension b0(5),b1(5),b2(5)
*
*-----limits for lambda_5 are 1 mev < lambda_5 < 10 gev
*
      bqm= 4.7d0
      cqm= 1.55d0
*
      if(als.eq.0.d0) then
         wtoralphas= 0.d0
      else
         x1= 0.001d0
         x2= 10.0d0
         xacc= 1.d-12
         qcdl= wtoqcdlam(nf,als,rs0,x1,x2,xacc)
         pqcdl5= qcdl
         do i=1,5
             b0(i)= (11.d0-2.d0/3.d0*i)/4.d0
             b1(i)= (102.d0-38.d0/3.d0*i)/16.d0
             b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0-
     #              325.d0/27.d0*i))/64.d0
         enddo
*
         if(rs.lt.bqm) then
            rat= bqm/qcdl
            rl= 2.d0*log(rat)
            rll= log(rl)
            rb= log(b0(5)/b0(4))
            fac= b1(5)/b0(5)-b1(4)/b0(4)
            facp= b2(5)/b0(5)-b2(4)/b0(4)
            fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4)
            rhs= (b0(5)-b0(4))*rl+fac*rll-b1(4)/b0(4)*rb+
     #           b1(5)/b0(5)/b0(5)*fac*rll/rl+1.d0/b0(5)/rl*(
     #           fac2-facp-7.d0/72.d0)
            rhs= rhs/b0(4)
            ratl2= exp(rhs)
            qcdl= qcdl/sqrt(ratl2)      
            pqcdl4= qcdl
            nfe= nf-1
            if(rs.lt.cqm) then
               rat= cqm/qcdl
               rl= 2.d0*log(rat)
               rll= log(rl)
               rb= log(b0(4)/b0(3))
               fac= b1(4)/b0(4)-b1(3)/b0(3)
               facp= b2(4)/b0(4)-b2(3)/b0(3)
               fac2= b1(4)*b1(4)/b0(4)/b0(4)-b1(3)*b1(3)/b0(3)/b0(3)
               rhs= (b0(4)-b0(3))*rl+fac*rll-b1(3)/b0(3)*rb+
     #              b1(4)/b0(4)/b0(4)*fac*rll/rl+1.d0/b0(4)/rl*(
     #              fac2-facp-7.d0/72.d0)
               rhs= rhs/b0(3)
               ratl2= exp(rhs)
               qcdl= qcdl/sqrt(ratl2)      
               pqcdl3= qcdl
               nfe= nf-2
            endif
         else
            nfe= nf
         endif
*
         qcdb0= 11.d0-2.d0/3.d0*nfe
         qcdb1= 102.d0-38.d0/3.d0*nfe
         qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+
     #          325.d0/27.d0*nfe*nfe)
         qcda= 2.d0*log(rs/qcdl)
*
         wtoralphas= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*
     #            log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-
     #            0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0))
*
      endif
*
      return
      end
*