*
*-------------------------------------------------------------------------
*
      subroutine zzto(rsi,zmi,wmi,zgi,tqmi,hmi,alszi,ipm,irm,sig,esig)
      implicit real*8(a-h,o-z)
      character*1 ofl,otop,oqed,out,oscheme
*
      parameter (ndimmx=4)
      parameter(nbf=2,nt=1,lwa=(nbf*(3*nbf+13))/2,lwat=(nt*(3*nt+13))/2)
*
      common/zztout/out
      common/zztafsz/g2z
      common/zztoqed/oqed
      common/zzted/eps,delta
      common/zztosch/oscheme
      common/zztofl/ofl,otop
      common/zztiaz/alzr,alzi,alzir,alzii
      common/zztcmplx/sw(2),sz(2),rsw(2),rsz(2)
      common/zztmass/em,rmm,tm,uqm,dqm,cqm,sqm,bqm
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
      common/zztstf/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta
      common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2
*
      dimension vk(ndimmx)
      dimension zs(2),zl(2)
      dimension st(nt),fvect(nt),wat(lwat)
      dimension fvecw(nbf),waw(lwa),waz(lwa),fvecz(nbf)
      dimension pggfz(2),pgglqz(2),p3qz(2),fzz(2),fwz(2)
*
      external s14aaf
      external g05cbf,zztoborn
      external d01gcf,zztoregion
      external c05nbf,c02ajf,x02ajf,zztofsw,zztofsz,zztofst
*
      wm= wmi
      zm= zmi
      hm= hmi
      zg= zgi
      alsz= alszi
      tqm= tqmi
      rs= rsi
      s= rs*rs
      delta= 0.d0
      eps= 1.d-30
      gf= 1.16639d-5/sqrt(2.d0)
      em= 0.51099906d-3
      pi= 3.141592653589793238462643d0
      pis= pi*pi
      pih= 0.5d0*pi
      tm= 1.7771d0
      bqm= 4.7d0
      cqm= 1.55d0
      alphai= 137.0359895d0
      alpha= 1.d0/alphai
      rmm= 0.10565839d0
      uqm= 0.005d0
      dqm= 0.010d0
      sqm= 0.15d0
      dmy= 0.d0
*
      em2= em*em
      rmm2= rmm*rmm
      tm2= tm*tm
      uqm2= uqm*uqm
      dqm2= dqm*dqm
      sqm2= sqm*sqm
      cqm2= cqm*cqm
      bqm2= bqm*bqm
      dmy2= dmy*dmy
*
      wm2= wm*wm
      zm2= zm*zm
      hm2= hm*hm
      alwi= 128.07d0
      alw= 1.d0/alwi
      if(oscheme.eq.'a') then
       sth2= 0.5d0*pi*alw/gf/wm2
       g2= 4.d0*pi*alw/sth2
      else if(oscheme.eq.'g') then
       sth2= 1.d0-wm2/zm2
       g2= 8.d0*gf*wm2
      endif
      cth2= 1.d0-sth2
      cth4= cth2*cth2
      cth8= cth4*cth4
      g4= g2*g2
      g8= g4*g4
*
      aexp= alpha/pi
      rln= log(s/em2)
      beta= 2.d0*aexp*(rln-1.d0)
      eta= beta
      eob= 1.d0
      feta= eta/16.d0
      hbet= beta/2.d0
      heta= eta/2.d0
      hbeti= 2.d0/beta
      omhb= 1.d0-hbet
      ophb= 1.d0+hbet
      ge= 0.5772156649d0
      arge= hbet*(3.d0/4.d0-ge)
      ifg= 1
      gamb= s14aaf(ophb,ifg)
      if(ifg.ne.0) then
         stop
      else 
         d0gl= exp(arge)/gamb
      endif
* 
      if(ofl.eq.'y') then
       wmc= wm2*wm
       wg= 1.5d0*gf*wmc/pi
       if(otop.eq.'f') then
        tqm2= tqm*tqm
       else if(otop.eq.'d') then
        gmt= 150.d0
        st(1)= gmt*gmt
        tol= sqrt(x02ajf())
        ifail= 1
        call c05nbf(zztofst,nt,st,fvect,tol,wat,lwat,ifail)
        tqm2= st(1)
        tqm= sqrt(tqm2)
       endif
*
       pzr= -zm2
       pzi= 0.d0
       call zztopoleg(pzr,pzi,pggfz,pgglqz,pggnpz)
       call zztopolez0(pggf0,fw0)
       xal= alphai-0.25d0*(pggfz(1)-pggf0+pggnpz)/pi
       yal= -0.25d0*(pggfz(2)+pgglqz(2))/pi
       alzir= xal
       alzii= yal
       alzr= xal/(xal*xal+yal*yal)
       alzi= -yal/(xal*xal+yal*yal)
*
       apis= 16.d0*pis
       p2zr= -zm2
       p2zi= 0.d0
       call zztopole(p2zr,p2zi,p3qz,fzz,fwz)
       bx= 1.d0/8.d0/gf/zm2+(fw0-fzz(1))/zm2/apis
       xih= -p3qz(2)/apis
       agc= 4.d0*pi*alzr
       bgc= -1.d0-8.d0*pi*xih*alzi
       cgc= -4.d0*pi*xih*xih*alzr+bx
       ifail= 0
       call c02ajf(agc,bgc,cgc,zs,zl,ifail)
       g2z= zs(1)
*
       sw(1)= wm2
       sw(2)= -wm*wg
       tol= sqrt(x02ajf())
       ifail= 1
       call c05nbf(zztofsw,nbf,sw,fvecw,tol,waw,lwa,ifail)
       wmu= sqrt(sw(1))
       wlg= -sw(2)/wmu
       do l=1,2
        rsw(l)= sw(l)/s
       enddo
*
       dszi= 1.d-1
       szi0= zm*zg
       izt= 1
    2  sz(1)= zm2
       sz(2)= -szi0*(1.d0-(izt-1)*dszi)
       tol= sqrt(x02ajf())
       ifail= 1
       call c05nbf(zztofsz,nbf,sz,fvecz,tol,waz,lwa,ifail)
       if(ifail.eq.0) then
        zmu= sqrt(sz(1))
        zlg= -sz(2)/zmu
        do l=1,2
         rsz(l)= sz(l)/s
        enddo
       else
        if(izt.lt.20) then
         izt= izt+1
         go to 2
        else
         print*,' Failure in searching for Z pole'
        endif
       endif
       if(out.eq.'y') then
        print 3,wmu,wlg
        print 4,zmu,zlg
        print 5,tqm
       endif
    3 format(/' W pole     ',2e20.5)
    4 format(/' Z pole     ',2e20.5)
    5 format(/' Derived top mass (GeV)   =',e20.5)
      endif
      tqm2= tqm*tqm
*
      if(oqed.eq.'n') then
       ndim= 2
      else if(oqed.eq.'y') then
       ndim= 4
      endif
      npts= ipm
      nrand= irm
      itrans= 0 
*
      ifail= 0 
      call g05cbf(0)
      call d01gcf(ndim,zztoborn,zztoregion,npts,vk,nrand,
     #            itrans,sig,esig,ifail) 
      print 1,ifail
    1 format(/' On exit IFAIL = ',i1,/)
*
      return
      end
*
*----------------------------------------------------------------
*
      real*8 function zztoborn(ndim,x)
      implicit real*8(a-h,o-z)
      character*4 ochan
      character*1 ofl,otop,onu,oqedf,oqcdf
*
      common/zztn/onu
      common/zztch/ochan
      common/zztk/ik,icall
      common/zztqedfsr/oqedf
      common/zztqedfsr/oqcdf
      common/zztofl/ofl,otop
      common/zztcmplx/sw(2),sz(2),rsw(2),rsz(2)
      common/zztmass/em,rmm,tm,uqm,dqm,cqm,sqm,bqm
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
      common/zztstf/alpha,hbet,hbeti,omhb,eob,d0gl,eta,feta,beta
*
      dimension x(ndim),fmu(9),scal(9),oim(9),cqms(9),sqms(9),
     #          af(9),bf(9),acqfr(9),acqfi(9),bcqfr(9),bcqfi(9),
     #          asqfr(9),asqfi(9),bsqfr(9),bsqfi(9),aszz(9,9),
     #          wh(9,9),qfs(9),pm(9)
      dimension pggfzm(2),pgglqzm(2),fw(2),p3qw(2),fzw(2),
     #          p3qsw(2),p3qsz(2),fzsz(2),fwsz(2),pggfsz(2),
     #          pgglqsz(2),p3qcq(2),p3qsq(2),fzcq(2),fzsq(2),
     #          fwcq(2),fwsq(2),pggfcq(2),pggfsq(2),pgglqcq(2),
     #          pgglqsq(2)
*
      icall= icall+1
*
      s= rs*rs
      ccw= gf*wm2/pis
      apis= 16.d0*pis
      rz2= pis/6.d0
*
      do i=1,9
       do j=1,9
        aszz(i,j)= 0.d0
       enddo
      enddo
*
      qfs(1)= 0.d0
      qfs(2)= 1.d0
      qfs(3)= 1.d0
      qfs(4)= 1.d0
      qfs(5)= 4.d0/9.d0
      qfs(6)= 1.d0/9.d0
      qfs(7)= 4.d0/9.d0
      qfs(8)= 1.d0/9.d0
      qfs(9)= 1.d0/9.d0
*
      pm(1)= 0.d0
      pm(2)= em
      pm(3)= rmm
      pm(4)= tm
      pm(5)= uqm
      pm(6)= dqm
      pm(7)= cqm
      pm(8)= sqm
      pm(9)= bqm
*
      fmu(1)= 3.d0
      fmu(2)= 1.d0
      fmu(3)= 1.d0
      fmu(4)= 1.d0
      do i=5,9
       fmu(i)= 3.d0
      enddo
*
      do i=1,9
       do j=1,9
        wh(i,j)= 0.d0
       enddo
      enddo
      if(onu.eq.'y') then
       wh(1,1)= 9.d0/2.d0
      else if(onu.eq.'n') then
       wh(1,1)= 0.d0
      endif
      do j=2,9
       wh(1,j)= 3.d0*fmu(j)
      enddo
      do i=2,9
       do j=i,9
        wh(i,j)= fmu(i)*fmu(j)
       enddo
      enddo
*
      do i=1,9
       if(i.lt.7) then
        scal(i)= 0.d0
       else if(i.eq.7) then
        scal(i)= 2.d0*cqm
       else if(i.eq.8) then
        scal(i)= 0.d0
       else if(i.eq.9) then
        scal(i)= 2.d0*bqm
       endif
      enddo     
      em2= em*em
      rmm2= rmm*rmm
      tm2= tm*tm
      uqm2= uqm*uqm
      dqm2= dqm*dqm
      sqm2= sqm*sqm
      cqm2= cqm*cqm
      bqm2= bqm*bqm
      rzm2= zm2/s
*
      oim(1)= 0.d0
      oim(2)= 4.d0*em2/s
      oim(3)= 4.d0*rmm2/s
      oim(4)= 4.d0*tm2/s
      oim(5)= 4.d0*uqm2/s
      oim(6)= 4.d0*dqm2/s
      cqms(1)= 0.d0
      cqms(2)= em2/s
      cqms(3)= rmm2/s
      cqms(4)= tm2/s
      cqms(5)= uqm2/s
      cqms(6)= dqm2/s
      sqms(1)= 0.d0
      sqms(2)= em2/s
      sqms(3)= rmm2/s
      sqms(4)= tm2/s
      sqms(5)= uqm2/s
      sqms(6)= dqm2/s
*
      scali= scal(7)
      call zztocorrqcd(scali,zm,alsz,bqm2,cqm2,rbm2,rcm2,rsm2)
      oim(7)= 4.d0*rcm2/s
      oim(8)= 4.d0*sqm2/s
      cqms(8)= sqm2/s
      sqms(8)= sqm2/s
      scali= scal(9)
      call zztocorrqcd(scali,zm,alsz,bqm2,cqm2,rbm2,rcm2,rsm2)
      oim(9)= 4.d0*rbm2/s
*
      af(1)= 0.5d0
      bf(1)= 0.5d0
      af(2)= -0.5d0+2.d0*sth2
      bf(2)= -0.5d0
      af(3)= -0.5d0+2.d0*sth2
      bf(3)= -0.5d0
      af(4)= -0.5d0+2.d0*sth2
      bf(4)= -0.5d0
      af(5)= 0.5d0-4.d0/3.d0*sth2
      bf(5)= 0.5d0
      af(6)= -0.5d0+2.d0/3.d0*sth2
      bf(6)= -0.5d0
      af(7)= 0.5d0-4.d0/3.d0*sth2
      bf(7)= 0.5d0
      af(8)= -0.5d0+2.d0/3.d0*sth2
      bf(8)= -0.5d0
      af(9)= -0.5d0+2.d0/3.d0*sth2
      bf(9)= -0.5d0
*
      if(ndim.eq.2) then
       xcqx= x(1)
       xsqx= x(2)
      else if(ndim.eq.4) then
       uvx= x(1)
       vvx= x(2)
       xcqx= x(3)
       xsqx= x(4)
      endif
*
      do i= 1,9
       do j= i,9
*
        if(ndim.eq.2) then
         ueps= 0.d0
         uv= 1.d0
         uvs= uv*uv
         ujc= 1.d0
         veps= 0.d0
         vv= 1.d0
         vjc= 1.d0
        else if(ndim.eq.4) then
         ul= (sqrt(oim(i))+sqrt(oim(j)))*
     #       (sqrt(oim(i))+sqrt(oim(j)))
         omul= 1.d0-ul
         ueps= omul*(1.d0-uvx)**hbeti
         uv= 1.d0-ueps
         ujc= omul**hbet
         uvs= uv*uv
*
         vvl= (sqrt(oim(i))+sqrt(oim(j)))*
     #        (sqrt(oim(i))+sqrt(oim(j)))
         vvll= vvl
*
         vvu1= uv
         vvu= uv
         uvl= uv-vvll
*
         vkf= (uv-vvu)/uvl
         if(vkf.eq.0.d0) then
          if(vvx.eq.1.d0) then
           veps= 0.d0
          else
           veps= uvl*(1.d0-vvx)**hbeti
          endif
          avkf= 1.d0
         else
          veps= uvl*(1.d0-(1.d0-vkf**hbet)*vvx)**hbeti
          avkf= 1.d0-vkf**hbet
         endif
         vv= uv-veps
         vjc0= 1.d0-vvll/uv
         vjc= vjc0**hbet*avkf
        endif
*
        vvs= vv*vv
        xm= uv
        xp= vv/uv  
        sh= vv*s
        xmop= xm/xp
        opxp= 1.d0+xp
        opxm= 1.d0+xm
        omxp= veps/uv
        omxm= ueps
        if(ndim.eq.2) then
         stfp= 1.d0
         stfm= 1.d0
        else if(ndim.eq.4) 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
*
       xcql= oim(i)/vv
       xcqu= (1.d0-sqrt(oim(j)/vv))*(1.d0-sqrt(oim(j)/vv))
       if(ofl.eq.'n') then
        zmas= vv*xcql-rzm2
        zmbs= vv*xcqu-rzm2
       else if(ofl.eq.'y') then
        zmas= vv*xcql-rsz(1)
        zmbs= vv*xcqu-rsz(1)
       endif
       szg= zg/zm
       szgs= szg*szg
       rzm= zm/rs
       rzg= zg/rs
       szg= zg/zm
       szgs= szg*szg
       rzmg= rzm*rzg
       opszgs= 1.d0+szgs
       s0z= rzm2/opszgs
       srzis= rsz(2)*rsz(2)
       aszi= abs(rsz(2))
       if(ofl.eq.'y') then
         vzmg= aszi*vv
       else
         vzmg= rzmg*vv
       endif
       if(ofl.eq.'n') then
        atma= (zmas+szgs*xcql)/rzmg
        atmb= (zmbs+szgs*xcqu)/rzmg
       else if(ofl.eq.'y') then
        atma= zmas/aszi
        atmb= zmbs/aszi
       endif
       if(atma.gt.1.d0.and.atmb.gt.1.d0) then
        atma= 1.d0/atma
        atma= atan(atma)
        zmat= (pih-atma)/vzmg
        atmb= 1.d0/atmb
        atmb= atan(atmb)
        zmbt= (pih-atmb)/vzmg
        xcqjc= (-atmb+atma)/vzmg
       else if(atma.gt.1.d0.and.atmb.lt.-1.d0) then
        atma= 1.d0/atma
        atma= atan(atma)
        zmat= (pih-atma)/vzmg
        atmb= -1.d0/atmb
        atmb= atan(atmb)
        zmbt= (-pih+atmb)/vzmg
        xcqjc= (-pi+atmb+atma)/vzmg
       else if(atma.gt.1.d0.and.abs(atmb).lt.1.d0) then
        atma= 1.d0/atma
        atma= atan(atma)
        zmat= (pih-atma)/vzmg
        atmb= atan(atmb)
        zmbt= atmb/vzmg
        xcqjc= (-pih+atmb+atma)/vzmg
       else if(atma.lt.-1.d0.and.atmb.gt.1.d0) then
        atma= -1.d0/atma
        atma= atan(atma)
        zmat= (-pih+atma)/vzmg
        atmb= 1.d0/atmb
        atmb= atan(atmb)
        zmbt= (pih-atmb)/vzmg
        xcqjc= (pi-atmb-atma)/vzmg
       else if(atma.lt.-1.d0.and.atmb.lt.-1.d0) then
        atma= -1.d0/atma
        atma= atan(atma)
        zmat= (-pih+atma)/vzmg
        atmb= -1.d0/atmb
        atmb= atan(atmb)
        zmbt= (-pih+atmb)/vzmg
        xcqjc= (atmb-atma)/vzmg
       else if(atma.lt.-1.d0.and.abs(atmb).lt.1.d0) then
        atma= -1.d0/atma
        atma= atan(atma)
        zmat= (-pih+atma)/vzmg
        atmb= atan(atmb)
        zmbt= atmb/vzmg
        xcqjc= (pih+atmb-atma)/vzmg
       else if(abs(atma).lt.1.d0.and.atmb.gt.1.d0) then
        atma= atan(atma)
        zmat= atma/vzmg
        atmb= 1.d0/atmb
        atmb= atan(atmb)
        zmbt= (pih-atmb)/vzmg
        xcqjc= (pih-atmb-atma)/vzmg
       else if(abs(atma).lt.1.d0.and.atmb.lt.-1.d0) then
        atma= atan(atma)
        zmat= atma/vzmg
        atmb= -1.d0/atmb
        atmb= atan(atmb)
        zmbt= (-pih+atmb)/vzmg
        xcqjc= (-pih+atmb-atma)/vzmg
       else if(abs(atma).lt.1.d0.and.abs(atmb).lt.1.d0) then
        atma= atan(atma)
        zmat= atma/vzmg
        atmb= atan(atmb)
        zmbt= atmb/vzmg
        xcqjc= (atmb-atma)/vzmg
       endif
       zmv= xcqjc*xcqx+zmat
       iftn= 1
       atnm= vzmg*zmv
       if(ofl.eq.'n') then
        xcq= s0z/vv*(1.d0+szg*s07aaf(atnm,iftn))
       else if(ofl.eq.'y') then
        xcq= (rsz(1)+aszi*s07aaf(atnm,iftn))/vv
       endif
       xcqjac= vv*xcqjc
       xcqs= xcq*xcq
       xcqc= xcqs*xcq
*
       xsql= oim(j)/vv
       if(xcq.lt.0.d0) then
        zztoborn= 0.d0
        ik= ik+1
        return
       endif
       xsqu= (1.d0-sqrt(xcq))*(1.d0-sqrt(xcq))
       if(ofl.eq.'n') then
        zpas= vv*xsql-rzm2
        zpbs= vv*xsqu-rzm2
       else if(ofl.eq.'y') then
        zpas= vv*xsql-rsz(1)
        zpbs= vv*xsqu-rsz(1)
       endif
       if(ofl.eq.'n') then
        atpa= (zpas+szgs*xsql)/rzmg
        atpb= (zpbs+szgs*xsqu)/rzmg
       else if(ofl.eq.'y') then
        atpa= zpas/aszi
        atpb= zpbs/aszi
       endif
       if(atpa.gt.1.d0.and.atpb.gt.1.d0) then
        atpa= 1.d0/atpa
        atpa= atan(atpa)
        zpat= (pih-atpa)/vzmg
        atpb= 1.d0/atpb
        atpb= atan(atpb)
        zpbt= (pih-atpb)/vzmg
        xsqjc= (-atpb+atpa)/vzmg
       else if(atpa.gt.1.d0.and.atpb.lt.-1.d0) then
        atpa= 1.d0/atpa
        atpa= atan(atpa)
        zpat= (pih-atpa)/vzmg
        atpb= -1.d0/atpb
        atpb= atan(atpb)
        zpbt= (-pih+atpb)/vzmg
        xsqjc= (-pi+atpb+atpa)/vzmg
       else if(atpa.gt.1.d0.and.abs(atpb).lt.1.d0) then
        atpa= 1.d0/atpa
        atpa= atan(atpa)
        zpat= (pih-atpa)/vzmg
        atpb= atan(atpb)
        zpbt= atpb/vzmg
        xsqjc= (-pih+atpb+atpa)/vzmg
       else if(atpa.lt.-1.d0.and.atpb.gt.1.d0) then
        atpa= -1.d0/atpa
        atpa= atan(atpa)
        zpat= (-pih+atpa)/vzmg
        atpb= 1.d0/atpb
        atpb= atan(atpb)
        zpbt= (pih-atpb)/vzmg
        xsqjc= (pi-atpb-atpa)/vzmg
       else if(atpa.lt.-1.d0.and.atpb.lt.-1.d0) then
        atpa= -1.d0/atpa
        atpa= atan(atpa)
        zpat= (-pih+atpa)/vzmg
        atpb= -1.d0/atpb
        atpb= atan(atpb)
        zpbt= (-pih+atpb)/vzmg
        xsqjc= (atpb-atpa)/vzmg
       else if(atpa.lt.-1.d0.and.abs(atpb).lt.1.d0) then
        atpa= -1.d0/atpa
        atpa= atan(atpa)
        zpat= (-pih+atpa)/vzmg
        atpb= atan(atpb)
        zpbt= atpb/vzmg
        xsqjc= (pih+atpb-atpa)/vzmg
       else if(abs(atpa).lt.1.d0.and.atpb.gt.1.d0) then
        atpa= atan(atpa)
        zpat= atpa/vzmg
        atpb= 1.d0/atpb
        atpb= atan(atpb)
        zpbt= (pih-atpb)/vzmg
        xsqjc= (pih-atpb-atpa)/vzmg
       else if(abs(atpa).lt.1.d0.and.atpb.lt.-1.d0) then
        atpa= atan(atpa)
        zpat= atpa/vzmg
        atpb= -1.d0/atpb
        atpb= atan(atpb)
        zpbt= (-pih+atpb)/vzmg
        xsqjc= (-pih+atpb-atpa)/vzmg
       else if(abs(atpa).lt.1.d0.and.abs(atpb).lt.1.d0) then
        atpa= atan(atpa)
        zpat= atpa/vzmg
        atpb= atan(atpb)
        zpbt= atpb/vzmg
        xsqjc= (atpb-atpa)/vzmg
       endif
       zpv= xsqjc*xsqx+zpat
       iftn= 1
       atnp= vzmg*zpv
       if(ofl.eq.'n') then
        xsq= s0z/vv*(1.d0+szg*s07aaf(atnp,iftn))
       else if(ofl.eq.'y') then
        xsq= (rsz(1)+aszi*s07aaf(atnp,iftn))/vv
       endif
       if(xsq.lt.0.d0) then
        zztoborn= 0.d0
        ik= ik+1
        return
       endif
       xsqjac= vv*xsqjc
       xsqs= xsq*xsq
       xsqc= xsqs*xsq
*
       if(i.eq.7) then
        scali= sqrt(xcq*vv*s)
        call zztocorrqcd(scali,zm,alsz,bqm2,cqm2,rbm2,rcm2,rsm2)
        cqms(7)= rcm2/s
       else if(i.eq.9) then
        scali= sqrt(xcq*vv*s)
        call zztocorrqcd(scali,zm,alsz,bqm2,cqm2,rbm2,rcm2,rsm2)
        cqms(9)= rbm2/s
       endif
       if(j.eq.7) then
        scalj= sqrt(xsq*vv*s)
        call zztocorrqcd(scalj,zm,alsz,bqm2,cqm2,rbm2,rcm2,rsm2)
        sqms(7)= rcm2/s
       else if(j.eq.9) then
        scalj= sqrt(xsq*vv*s)
        call zztocorrqcd(scalj,zm,alsz,bqm2,cqm2,rbm2,rcm2,rsm2)
        sqms(9)= rbm2/s
       endif
*
       if(oqedf.eq.'y') then
        call zztopolez0(pggf0,fw0)
        p2cqr= xcq*vv*s
        p2cqi= 0.d0
        rp2cqr= -sqrt(p2cqr)
        call zztopoleg(p2cqr,p2cqi,pggfcq,pgglqcq,pggnpcq)
        if(abs(rp2cqr).gt.40.d0) then
         call zztohadr5(rp2cqr,derh,ederh)
         derh= 0.25d0*alpha/pi*pggnpcq
         derl= 0.25d0*alpha/pi*(pggfcq(1)-pggf0)
         der= derl+derh
         alphacq= alpha/(1.d0-der)
        else
         arecq= alphai-0.25d0*(pggfcq(1)-pggf0+pggnpcq)/pi 
         alphacq= 1.d0/arecq
        endif
        p2sqr= xsq*vv*s
        p2sqi= 0.d0
        rp2sqr= -sqrt(p2sqr)
        call zztopoleg(p2sqr,p2sqi,pggfsq,pgglqsq,pggnpsq)
        if(abs(rp2sqr).gt.40.d0) then
         call zztohadr5(rp2sqr,derh,ederh)
         derh= 0.25d0*alpha/pi*pggnpsq
         derl= 0.25d0*alpha/pi*(pggfsq(1)-pggf0)
         der= derl+derh
         alphasq= alpha/(1.d0-der)
        else
         aresq= alphai-0.25d0*(pggfsq(1)-pggf0+pggnpsq)/pi 
         alphasq= 1.d0/aresq
        endif
*
        zim= -1.d-30
        zcq= oim(i)
        zsq= oim(j)
        zcqs= zcq*zcq
        zsqs= zsq*zsq
        omzcq= 1.d0-zcq
        omzcqs= omzcq*omzcq
        omzsq= 1.d0-zsq
        omzsqs= omzsq*omzsq
        if(i.eq.1) then
         factcq= 0.d0
        else
         factcq= -2.d0*omzcqs+4.d0*((zcq+0.5d0*zcqs+2.d0*
     #           log(omzcq))*log(vv/cqms(i))+zcq*(1.d0+0.5d0*zcq)*
     #           log(zcq)+2.d0*rz2-2.d0*zztoli2(omzcq,zim)-
     #           2.d0*log(omzcq)+5.d0/4.d0-3.d0*zcq-0.25d0*zcqs)
        endif
        if(j.eq.1) then
         factsq= 0.d0
        else
         factsq= -2.d0*omzsqs+4.d0*((zsq+0.5d0*zsqs+2.d0*
     #           log(omzsq))*log(vv/sqms(j))+zsq*(1.d0+0.5d0*zsq)*
     #           log(zsq)+2.d0*rz2-2.d0*zztoli2(omzsq,zim)-
     #           2.d0*log(omzsq)+5.d0/4.d0-3.d0*zsq-0.25d0*zsqs)
        endif
        qedcq= 1.d0+alphacq/4.d0/pi*qfs(i)*factcq
        qedsq= 1.d0+alphasq/4.d0/pi*qfs(j)*factsq
       endif
*
       if(oqcdf.eq.'y') then
        if(i.le.4) then
         qcdcq= 0.d0
        else
         nf= 5
         scalcq= sqrt(xcq*vv*s)
         if(i.eq.7.or.i.eq.9) then
          alscq= zztoralphas(zm,scalcq,alsz,nf)
         else
          if(scalcq.gt.1.d0) then
           alscq= zztoralphas(zm,scalcq,alsz,nf)
          else
           scalcq= 1.d0
           alscq= zztoralphas(zm,scalcq,alsz,nf)
          endif
         endif
         vcqs= 1.d0-4.d0*pm(i)*pm(i)/xcq/vv/s
         if(vcqs.lt.0.d0) then
          rvcq= 0.d0
          racq= 0.d0
         else
          vcq= sqrt(vcqs)
          vcqc= vcqs*vcq
          omvcq= 1.d0-vcq
          opvcq= 1.d0+vcq
          opvcqs= opvcq*opvcq
          opvcqc= opvcqs*opvcq
          rcq= omvcq/opvcq
          rcqs= rcq*rcq
          rim= -1.d-30
          if(vcq.lt.0.19d0) then
           bkvcq= 3.d0/4.d0*pi/alscq*(2.d0/3.d0*pi*alscq/vcq+
     #            (1.d0-16.d0/3.d0*alscq/pi)-1.d0)
           bkacq= 3.d0/4.d0*pi/alscq*(2.d0/3.d0*pi*alscq/vcq+
     #            (1.d0-8.d0/3.d0*alscq/pi)-1.d0)
          else if(vcq.gt.0.19d0.and.vcq.lt.0.91d0) then
           avcq= (1.d0+vcqs)*(zztoli2(rcqs,rim)+2.d0*
     #           zztoli2(rcq,rim)+log(rcq)*log(opvcqc/8.d0/vcqs))+
     #           3.d0*vcq*log((1.d0-vcqs)/4.d0/vcq)-vcq*log(vcq)
           pvvcq= 33.d0/24.d0+vcqs*(22.d0/24.d0-7.d0/24.d0*vcqs)   
           pavcq= 21.d0/32.d0+vcqs*(59.d0/32.d0-vcqs*(19.d0/32.d0+
     #            3.d0/32.d0*vcqs))   
           qvvcq= vcq*(5.d0/4.d0-3.d0/4.d0*vcqs)
           qavcq= vcq*(-21.d0/16.d0+vcqs*(30.d0/16.d0+
     #            3.d0/16.d0*vcqs))
           bkvcq= 1.d0/vcq*(avcq-pvvcq/(1.d0-vcqs/3.d0)*log(rcq)+
     #            qvvcq/(1.d0-vcqs/3.d0))
           bkacq= 1.d0/vcq*(avcq-pavcq/vcqs*log(rcq)+qavcq/vcqs)
          else if(vcq.gt.0.91d0) then
           if(i.ne.7.and.i.ne.9) then
            pmi= pm(i)
            pmis= pmi*pmi/xcq/vv/s
            pmiq= pmis*pmis
            bkvcq= 3.d0/4.d0*pi/alscq*(2.d0/vcq/(3.d0-vcqs)*(
     #             1.d0-6.d0*pmiq+alscq/pi*(1.d0+12.d0*pmis+
     #             pmiq*(10.d0-24.d0*log(pmis))))-1.d0)
            bkacq= 3.d0/4.d0*pi/alscq*(1.d0/vcqc*(1.d0-6.d0*pmis+
     #             6.d0*pmiq+alscq/pi*(1.d0-pmis*(6.d0+12.d0*
     #             log(pmis))+pmiq*(-22.d0+24.d0*log(pmis))))-1.d0)
           else 
            pmis= pm(i)*pm(i)*(1.d0-alscq/pi*(8.d0/3.d0-2.d0)*
     #            log(pm(i)*pm(i)/xcq/vv/s))/xcq/vv/s
            pmiq= pmis*pmis
            bkvcq= 3.d0/4.d0*pi/alscq*(2.d0/vcq/(3.d0-vcqs)*(
     #             1.d0-6.d0*pmiq+alscq/pi*(1.d0+12.d0*pmis-22.d0*
     #             pmiq))-1.d0)
            bkacq= 3.d0/4.d0*pi/alscq*(1.d0/vcqc*(1.d0-6.d0*pmis+
     #             6.d0*pmiq+alscq/pi*(1.d0-22.d0*pmis+
     #             10.d0*pmiq))-1.d0)
           endif
          endif
          rvcq= 4.d0/3.d0*alscq/pi*bkvcq
          racq= 4.d0/3.d0*alscq/pi*bkacq
         endif
        endif
*
        if(j.le.4) then
         qcdsq= 0.d0
        else
         nf= 5
         scalsq= sqrt(xsq*vv*s)
         if(j.eq.7.or.j.eq.9) then
          alssq= zztoralphas(zm,scalsq,alsz,nf)
         else
          if(scalsq.gt.1.d0) then
           alssq= zztoralphas(zm,scalsq,alsz,nf)
          else
           scalsq= 1.d0
           alssq= zztoralphas(zm,scalsq,alsz,nf)
          endif
         endif
         vsqs= 1.d0-4.d0*pm(i)*pm(i)/xsq/vv/s
         if(vsqs.lt.0.d0) then
          rvsq= 0.d0
          rasq= 0.d0
         else
          vsq= sqrt(vsqs)
          vsqs= vsq*vsq
          vsqc= vsqs*vsq
          omvsq= 1.d0-vsq
          opvsq= 1.d0+vsq
          opvsqs= opvsq*opvsq
          opvsqc= opvsqs*opvsq
          rsq= omvsq/opvsq
          rsqs= rsq*rsq
          rim= -1.d-30
          if(vsq.lt.0.19d0) then
           bkvsq= 3.d0/4.d0*pi/alssq*(2.d0/3.d0*pi*alssq/vsq+
     #            (1.d0-16.d0/3.d0*alssq/pi)-1.d0)
           bkasq= 3.d0/4.d0*pi/alssq*(2.d0/3.d0*pi*alssq/vsq+
     #            (1.d0-8.d0/3.d0*alssq/pi)-1.d0)
          else if(vsq.gt.0.19d0.and.vsq.lt.0.91d0) then
           avsq= (1.d0+vsqs)*(zztoli2(rsqs,rim)+2.d0*
     #           zztoli2(rsq,rim)+log(rsq)*log(opvsqc/8.d0/vsqs))+
     #           3.d0*vsq*log((1.d0-vsqs)/4.d0/vsq)-vsq*log(vsq)
           pvvsq= 33.d0/24.d0+vsqs*(22.d0/24.d0-7.d0/24.d0*vsqs)   
           pavsq= 21.d0/32.d0+vsqs*(59.d0/32.d0-vsqs*(19.d0/32.d0+
     #            3.d0/32.d0*vsqs))   
           qvvsq= vsq*(5.d0/4.d0-3.d0/4.d0*vsqs)
           qavsq= vsq*(-21.d0/16.d0+vsqs*(30.d0/16.d0+
     #            3.d0/16.d0*vsqs))
           bkvsq= 1.d0/vsq*(avsq-pvvsq/(1.d0-vsqs/3.d0)*log(rsq)+
     #            qvvsq/(1.d0-vsqs/3.d0))
           bkasq= 1.d0/vsq*(avsq-pavsq/vsqs*log(rsq)+qavsq/vsqs)
          else if(vsq.gt.0.91d0) then
           if(j.ne.7.and.j.ne.9) then
            pmj= pm(j)
            pmjs= pmj*pmj/xsq/vv/s
            pmjq= pmjs*pmjs
            bkvsq= 3.d0/4.d0*pi/alssq*(2.d0/vsq/(3.d0-vsqs)*(
     #             1.d0-6.d0*pmjq+alssq/pi*(1.d0+12.d0*pmjs+
     #             pmjq*(10.d0-24.d0*log(pmjs))))-1.d0)
            bkasq= 3.d0/4.d0*pi/alssq*(1.d0/vsqc*(1.d0-6.d0*pmjs+
     #             6.d0*pmjq+alssq/pi*(1.d0-pmjs*(6.d0+12.d0*
     #             log(pmjs))+pmjq*(-22.d0+24.d0*log(pmjs))))-1.d0)
           else 
            pmjs= pm(j)*pm(j)*(1.d0-alssq/pi*(8.d0/3.d0-2.d0)*
     #            log(pm(j)*pm(j)/xsq/vv/s))/xsq/vv/s
            pmjq= pmjs*pmjs
            bkvsq= 3.d0/4.d0*pi/alssq*(2.d0/vsq/(3.d0-vsqs)*(
     #             1.d0-6.d0*pmjq+alssq/pi*(1.d0+12.d0*pmjs-22.d0*
     #             pmjq))-1.d0)
            bkasq= 3.d0/4.d0*pi/alssq*(1.d0/vsqc*(1.d0-6.d0*pmjs+
     #             6.d0*pmjq+alssq/pi*(1.d0-22.d0*pmjs+
     #             10.d0*pmjq))-1.d0)
           endif
          endif
          rvsq= 4.d0/3.d0*alssq/pi*bkvsq
          rasq= 4.d0/3.d0*alssq/pi*bkasq
         endif
        endif
       endif
*
       if(ofl.eq.'y') then
        call zztopolez0(pggf0,fw0)
*
        p2r= -zm2
        p2i= 0.d0
        call zztopoleg(p2r,p2i,pggfzm,pgglqzm,pggnpzm)
        arezm= alphai-0.25d0*(pggfzm(1)-pggf0+pggnpzm)/pi
        aizm= -0.25d0*(pggfzm(2)+pgglqzm(2))/pi
*
        p2r= -wm2
        p2i= 0.d0
        call zztopole(p2r,p2i,p3qw,fzw,fw)
*
        p2r= -sw(1)
        p2i= -sw(2)
        call zztopole(p2r,p2i,p3qsw,fzsw,fwsw)
        xgsw0= 1.d0/8.d0/gf/wm2
        xgswr= 1.d0+0.5d0*ccw*((fw0-fw(1))/wm2+
     #         p3qw(1)-p3qsw(1))
        xgswi= -0.5d0*ccw*p3qsw(2)
        agrsw= xgsw0*xgswr
        agisw= xgsw0*xgswi
        agmsw= agrsw*agrsw+agisw*agisw
        grsw= agrsw/agmsw
        gisw= -agisw/agmsw
*
        p2r= -sz(1)
        p2i= -sz(2)
        call zztopole(p2r,p2i,p3qsz,fzsz,fwsz)
        xgsz0= 1.d0/8.d0/gf/wm2
        xgszr= 1.d0+0.5d0*ccw*((fw0-fw(1))/wm2+
     #         p3qw(1)-p3qsz(1))
        xgszi= -0.5d0*ccw*p3qsz(2)
        agrsz= xgsz0*xgszr
        agisz= xgsz0*xgszi
        agmsz= agrsz*agrsz+agisz*agisz
        grsz= agrsz/agmsz
        gisz= -agisz/agmsz
        call zztopoleg(p2r,p2i,pggfsz,pgglqsz,pggnpsz)
        aresz= arezm+0.25d0*(pggfzm(1)+pgglqzm(1)-
     #         pggfsz(1)-pgglqsz(1))/pi
        aisz= aizm+0.25d0*(pggfzm(2)+pgglqzm(2)-
     #        pggfsz(2)-pgglqsz(2))/pi
        ermz= aresz*aresz+aisz*aisz
        ersz= 4.d0*pi*aresz/ermz
        eisz= -4.d0*pi*aisz/ermz
        gmsz= grsz*grsz+gisz*gisz
        strrsz= ersz*agrsz-eisz*agisz
        strisz= ersz*agisz+eisz*agrsz
        ctrrsz= 1.d0-strrsz
        ctrisz= -strisz
*
        p2r= -xcq*sh
        p2i= 0.d0
        call zztopole(p2r,p2i,p3qcq,fzcq,fwcq)
        xgcq0= 1.d0/8.d0/gf/wm2
        xgcqr= 1.d0+0.5d0*ccw*((fw0-fw(1))/wm2+
     #         p3qw(1)-p3qcq(1))
        xgcqi= -0.5d0*ccw*p3qcq(2)
        agrcq= xgcq0*xgcqr
        agicq= xgcq0*xgcqi
        gmcq= agrcq*agrcq+agicq*agicq
        grcq= agrcq/gmcq
        gicq= -agicq/gmcq
        call zztopoleg(p2r,p2i,pggfcq,pgglqcq,pggnpcq)
        arecq= arezm+0.25d0*(pggfzm(1)+pgglqzm(1)-
     #         pggfcq(1)-pgglqcq(1))/pi
        aicq= aizm+0.25d0*(pggfzm(2)+pgglqzm(2)-
     #        pggfcq(2)-pgglqcq(2))/pi
        emcq= arecq*arecq+aicq*aicq
        ercq= 4.d0*pi*arecq/emcq
        eicq= -4.d0*pi*aicq/emcq
        strrcq= ercq*agrcq-eicq*agicq
        stricq= ercq*agicq+eicq*agrcq
        tstrrcq= 2.d0*strrcq
        tstricq= 2.d0*stricq
        ctrrcq= 1.d0-strrcq
        ctricq= -stricq
        ctrmcq= ctrrcq*ctrrcq+ctricq*ctricq
        ctrrcqi= ctrrcq/ctrmcq
        ctricqi= -ctricq/ctrmcq
        ttrrcq= strrcq*ctrrcqi-stricq*ctricqi
        ttricq= strrcq*ctricqi+stricq*ctrrcqi
        prcqr= xcq*sh-sz(1)
        prcqi= -sz(2)
        prcqm= prcqr*prcqr+prcqi*prcqi
        prcqri= prcqr/prcqm
        prcqii= -prcqi/prcqm
        ratgrcq= grcq*agrsz-gicq*agisz
        ratgicq= grcq*agisz+gicq*agrsz
        ratcrcq= ctrrsz*ctrrcqi-ctrisz*ctricqi
        ratcicq= ctrrsz*ctricqi+ctrisz*ctrrcqi
        ratgcrcq= grcq*ctrrcqi-gicq*ctricqi
        ratgcicq= grcq*ctricqi+gicq*ctrrcqi
        ratrcq= 1.d0-ratgrcq*ratcrcq+ratgicq*ratcicq
        raticq= -ratgrcq*ratcicq-ratgicq*ratcrcq
        arhrcq= xcq*sh-sz(1)+ratrcq*sz(1)-raticq*sz(2)+
     #          (ratgcrcq*(fzcq(1)-fzsz(1))-ratgcicq*(fzcq(2)-
     #          fzsz(2)))/apis
        arhicq= -sz(2)+ratrcq*sz(2)+raticq*sz(1)+(ratgcrcq*(fzcq(2)-
     #           fzsz(2))+ratgcicq*(fzcq(1)-fzsz(1)))/apis
        brhrcq= arhrcq*prcqri-arhicq*prcqii
        brhicq= arhrcq*prcqii+arhicq*prcqri
        brhmcq= brhrcq*brhrcq+brhicq*brhicq
*
        p2r= -xsq*sh
        p2i= 0.d0
        call zztopole(p2r,p2i,p3qsq,fzsq,fwsq)
        xgsq0= 1.d0/8.d0/gf/wm2
        xgsqr= 1.d0+0.5d0*ccw*((fw0-fw(1))/wm2+
     #         p3qw(1)-p3qsq(1))
        xgsqi= -0.5d0*ccw*p3qsq(2)
        agrsq= xgsq0*xgsqr
        agisq= xgsq0*xgsqi
        gmsq= agrsq*agrsq+agisq*agisq
        grsq= agrsq/gmsq
        gisq= -agisq/gmsq
        call zztopoleg(p2r,p2i,pggfsq,pgglqsq,pggnpsq)
        aresq= arezm+0.25d0*(pggfzm(1)+pgglqzm(1)-
     #         pggfsq(1)-pgglqsq(1))/pi
        aisq= aizm+0.25d0*(pggfzm(2)+pgglqzm(2)-
     #        pggfsq(2)-pgglqsq(2))/pi
        emsq= aresq*aresq+aisq*aisq
        ersq= 4.d0*pi*aresq/emsq
        eisq= -4.d0*pi*aisq/emsq
        strrsq= ersq*agrsq-eisq*agisq
        strisq= ersq*agisq+eisq*agrsq
        tstrrsq= 2.d0*strrsq
        tstrisq= 2.d0*strisq
        ctrrsq= 1.d0-strrsq
        ctrisq= -strisq
        ctrmsq= ctrrsq*ctrrsq+ctrisq*ctrisq
        ctrrsqi= ctrrsq/ctrmsq
        ctrisqi= -ctrisq/ctrmsq
        ttrrsq= strrsq*ctrrsqi-strisq*ctrisqi
        ttrisq= strrsq*ctrisqi+strisq*ctrrsqi
        prsqr= xsq*sh-sz(1)
        prsqi= -sz(2)
        prsqm= prsqr*prsqr+prsqi*prsqi
        prsqri= prsqr/prsqm
        prsqii= -prsqi/prsqm
        ratgrsq= grsq*agrsz-gisq*agisz
        ratgisq= grsq*agisz+gisq*agrsz
        ratcrsq= ctrrsz*ctrrsqi-ctrisz*ctrisqi
        ratcisq= ctrrsz*ctrisqi+ctrisz*ctrrsqi
        ratgcrsq= grsq*ctrrsqi-gisq*ctrisqi
        ratgcisq= grsq*ctrisqi+gisq*ctrrsqi
        ratrsq= 1.d0-ratgrsq*ratcrsq+ratgisq*ratcisq
        ratisq= -ratgrsq*ratcisq-ratgisq*ratcrsq
        arhrsq= xsq*sh-sz(1)+ratrsq*sz(1)-ratisq*sz(2)+
     #          (ratgcrsq*(fzsq(1)-fzsz(1))-ratgcisq*(fzsq(2)-
     #          fzsz(2)))/apis
        arhisq= -sz(2)+ratrsq*sz(2)+ratisq*sz(1)+(ratgcrsq*(fzsq(2)-
     #           fzsz(2))+ratgcisq*(fzsq(1)-fzsz(1)))/apis
        brhrsq= arhrsq*prsqri-arhisq*prsqii
        brhisq= arhrsq*prsqii+arhisq*prsqri
        brhmsq= brhrsq*brhrsq+brhisq*brhisq
       endif
*
       rlks= 1.d0+xcqs+xsqs-2.d0*(xcq+xsq+xcq*xsq)
       if(rlks.lt.0.d0) then
        zztoborn= 0.d0
        ik= ik+1
        return
       endif
       rlk= sqrt(rlks)
       tauu= 0.5d0*(-1.d0+xcq+xsq+rlk)
       taul= 0.5d0*(-1.d0+xcq+xsq-rlk)
       rltau= log(tauu/taul)
       taui= -1.d0/tauu+1.d0/taul
       rlups= -log(taul/tauu)
       upsi= 1.d0/taul-1.d0/tauu
       taus= 0.5d0*(tauu*tauu-taul*taul)
       upss= -0.5d0*(tauu*tauu-taul*taul)+(xcq+xsq-1.d0)*(tauu-taul)
       one= tauu-taul
*
       sqq= 1.d0-xcq-xsq
       sqq= 1.d0/sqq
*
       if(ofl.eq.'n') then
        ae= -0.5d0+2.d0*sth2
        be= -0.5d0
        abesqp= ae*ae+be*be
        abesqm= ae*ae-be*be
        rabesq= ae*be
        abecqp= ae*ae+be*be
        abecqm= ae*ae-be*be
        rabecq= ae*be
       else if(ofl.eq.'y') then
        aesqr= -0.5d0+2.d0*strrsq
        aesqi= 2.d0*strisq
        besqr= -0.5d0
        besqi= 0.d0
        aecqr= -0.5d0+2.d0*strrcq
        aecqi= 2.d0*stricq
        becqr= -0.5d0
        becqi= 0.d0
        abesqp= aesqr*aesqr+aesqi*aesqi+
     #          besqr*besqr+besqi*besqi
        abesqm= aesqr*aesqr+aesqi*aesqi-
     #          besqr*besqr-besqi*besqi
        rabesq= aesqr*besqr-aesqi*besqi
        abecqp= aecqr*aecqr+aecqi*aecqi+
     #          becqr*becqr+becqi*becqi
        abecqm= aecqr*aecqr+aecqi*aecqi-
     #          becqr*becqr-becqi*becqi
        rabecq= aecqr*becqr-aecqi*becqi
       endif
       if(ofl.eq.'n') then
        abcqp= af(i)*af(i)+bf(i)*bf(i)
        abcqm= af(i)*af(i)-bf(i)*bf(i)
        absqp= af(j)*af(j)+bf(j)*bf(j)
        absqm= af(j)*af(j)-bf(j)*bf(j)
        acqs= af(i)*af(i)
        bcqs= bf(i)*bf(i)
        asqs= af(j)*af(j)
        bsqs= bf(j)*bf(j)
       else if(ofl.eq.'y') then
        acqfr(1)= 0.5d0
        bcqfr(1)= 0.5d0
        acqfr(2)= -0.5d0+2.d0*strrcq
        bcqfr(2)= -0.5d0
        acqfr(3)= -0.5d0+2.d0*strrcq
        bcqfr(3)= -0.5d0
        acqfr(4)= -0.5d0+2.d0*strrcq
        bcqfr(4)= -0.5d0
        acqfr(5)= 0.5d0-4.d0/3.d0*strrcq
        bcqfr(5)= 0.5d0
        acqfr(6)= -0.5d0+2.d0/3.d0*strrcq
        bcqfr(6)= -0.5d0
        acqfr(7)= 0.5d0-4.d0/3.d0*strrcq
        bcqfr(7)= 0.5d0
        acqfr(8)= -0.5d0+2.d0/3.d0*strrcq
        bcqfr(8)= -0.5d0
        acqfr(9)= -0.5d0+2.d0/3.d0*strrcq
        bcqfr(9)= -0.5d0
        acqfi(1)= 0.d0
        bcqfi(1)= 0.d0
        acqfi(2)= 2.d0*stricq
        bcqfi(2)= 0.d0
        acqfi(3)= 2.d0*stricq
        bcqfi(3)= 0.d0
        acqfi(4)= 2.d0*stricq
        bcqfi(4)= 0.d0
        acqfi(5)= -4.d0/3.d0*stricq
        bcqfi(5)= 0.d0
        acqfi(6)= 2.d0/3.d0*stricq
        bcqfi(6)= 0.d0
        acqfi(7)= -4.d0/3.d0*stricq
        bcqfi(7)= 0.d0
        acqfi(8)= 2.d0/3.d0*stricq
        bcqfi(8)= 0.d0
        acqfi(9)= 2.d0/3.d0*stricq
        bcqfi(9)= 0.d0
*
        asqfr(1)= 0.5d0
        bsqfr(1)= 0.5d0
        asqfr(2)= -0.5d0+2.d0*strrsq
        bsqfr(2)= -0.5d0
        asqfr(3)= -0.5d0+2.d0*strrsq
        bsqfr(3)= -0.5d0
        asqfr(4)= -0.5d0+2.d0*strrsq
        bsqfr(4)= -0.5d0
        asqfr(5)= 0.5d0-4.d0/3.d0*strrsq
        bsqfr(5)= 0.5d0
        asqfr(6)= -0.5d0+2.d0/3.d0*strrsq
        bsqfr(6)= -0.5d0
        asqfr(7)= 0.5d0-4.d0/3.d0*strrsq
        bsqfr(7)= 0.5d0
        asqfr(8)= -0.5d0+2.d0/3.d0*strrsq
        bsqfr(8)= -0.5d0
        asqfr(9)= -0.5d0+2.d0/3.d0*strrsq
        bsqfr(9)= -0.5d0
        asqfi(1)= 0.d0
        bsqfi(1)= 0.d0
        asqfi(2)= 2.d0*strisq
        bsqfi(2)= 0.d0
        asqfi(3)= 2.d0*strisq
        bsqfi(3)= 0.d0
        asqfi(4)= 2.d0*strisq
        bsqfi(4)= 0.d0
        asqfi(5)= -4.d0/3.d0*strisq
        bsqfi(5)= 0.d0
        asqfi(6)= 2.d0/3.d0*strisq
        bsqfi(6)= 0.d0
        asqfi(7)= -4.d0/3.d0*strisq
        bsqfi(7)= 0.d0
        asqfi(8)= 2.d0/3.d0*strisq
        bsqfi(8)= 0.d0
        asqfi(9)= 2.d0/3.d0*strisq
        bsqfi(9)= 0.d0
*
        abcqp= acqfr(i)*acqfr(i)+bcqfr(i)*bcqfr(i)+
     #         acqfi(i)*acqfi(i)+bcqfi(i)*bcqfi(i) 
        abcqm= acqfr(i)*acqfr(i)-bcqfr(i)*bcqfr(i)+
     #         acqfi(i)*acqfi(i)-bcqfi(i)*bcqfi(i) 
        absqp= asqfr(j)*asqfr(j)+bsqfr(j)*bsqfr(j)+
     #         asqfi(j)*asqfi(j)+bsqfi(j)*bsqfi(j) 
        absqm= asqfr(j)*asqfr(j)-bsqfr(j)*bsqfr(j)+
     #         asqfi(j)*asqfi(j)-bsqfi(j)*bsqfi(j) 
        acqs= acqfr(i)*acqfr(i)+acqfi(i)*acqfi(i)
        bcqs= bcqfr(i)*bcqfr(i)+bcqfi(i)*bcqfi(i)
        asqs= asqfr(j)*asqfr(j)+asqfi(j)*asqfi(j)
        bsqs= bsqfr(j)*bsqfr(j)+bsqfi(j)*bsqfi(j)
       endif
*
       betcqs= 1.d0-4.d0*cqms(i)/xcq/vv
       if(betcqs.lt.0.d0) then
        zztoborn= 0.d0
        return
       endif
       betcq= sqrt(betcqs)
       betsqs= 1.d0-4.d0*sqms(j)/xsq/vv
       if(betsqs.lt.0.d0) then
        zztoborn= 0.d0
        return
       endif
       betsq= sqrt(betsqs)
       bet1= 1.d0-betcq-betsq+betcq*betsq
       bet2= betcq-1.d0
       bet3= betsq-1.d0
       betvcq= (3.d0-betcqs)/2.d0
       betvsq= (3.d0-betsqs)/2.d0
*
       xscq= xsq*xcq
       xscqs= xsqs*xcqs
       xscs= xsq*xcqs
       xcss= xsqs*xcq
       xscc= xsqc*xcq
       xcsc= xsq*xcqc
       xcmb1= xscq+xscs+xcss
       xcmb2= xcsc+2.d0*xscqs+xscc
       xcmb3= -xscq+1.d0/3.d0*(xsq-xsqs)
       xcmb4= -xscq-1.d0/3.d0*(xsq-xsqs)
*
       abcsqp= abcqp*absqp
       abcsqpm= abcqp*absqm
       abcsqmp= abcqm*absqp
       abcsqm= abcqm*absqm
       abecsqp= abecqp*abesqp
       rabescq= rabesq*rabecq
       betvcsq= betvcq*betvsq
       acsqs= acqs*asqs
       bcsqs= bcqs*bsqs
       abcsqs= acqs*bsqs
       bacsqs= bcqs*asqs
*
       szz= bet1*sqq*xcmb2*abcsqp*abecsqp*rltau*(-4./9.)
       szz= szz+bet1*sqq*xcmb2*abcsqp*abecsqp*rlups*(-4./9.)
       szz= szz+bet1*sqq*xcmb2*abcsqp*rabescq*rltau*(-16./9.)
       szz= szz+bet1*sqq*xcmb2*abcsqp*rabescq*rlups*(-16./9.)
       szz= szz+bet1*sqq*xcmb2*abcsqpm*abecsqp*rltau*(4./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqpm*abecsqp*rlups*(4./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqpm*rabescq*rltau*(16./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqpm*rabescq*rlups*(16./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqmp*abecsqp*rltau*(4./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqmp*abecsqp*rlups*(4./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqmp*rabescq*rltau*(16./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqmp*rabescq*rlups*(16./3.)
       szz= szz+bet1*sqq*xcmb2*abcsqm*abecsqp*rltau*(-4.)
       szz= szz+bet1*sqq*xcmb2*abcsqm*abecsqp*rlups*(-4.)
       szz= szz+bet1*sqq*xcmb2*abcsqm*rabescq*rltau*(-16.)
       szz= szz+bet1*sqq*xcmb2*abcsqm*rabescq*rlups*(-16.)
       szz= szz+bet1*xscqs*abcsqp*abecsqp*taui*(-2./9.)
       szz= szz+bet1*xscqs*abcsqp*abecsqp*upsi*(-2./9.)
       szz= szz+bet1*xscqs*abcsqp*rabescq*taui*(-8./9.)
       szz= szz+bet1*xscqs*abcsqp*rabescq*upsi*(-8./9.)
       szz= szz+bet1*xscqs*abcsqpm*abecsqp*taui*(2./3.)
       szz= szz+bet1*xscqs*abcsqpm*abecsqp*upsi*(2./3.)
       szz= szz+bet1*xscqs*abcsqpm*rabescq*taui*(8./3.)
       szz= szz+bet1*xscqs*abcsqpm*rabescq*upsi*(8./3.)
       szz= szz+bet1*xscqs*abcsqmp*abecsqp*taui*(2./3.)
       szz= szz+bet1*xscqs*abcsqmp*abecsqp*upsi*(2./3.)
       szz= szz+bet1*xscqs*abcsqmp*rabescq*taui*(8./3.)
       szz= szz+bet1*xscqs*abcsqmp*rabescq*upsi*(8./3.)
       szz= szz+bet1*xscqs*abcsqm*abecsqp*taui*(-2.)
       szz= szz+bet1*xscqs*abcsqm*abecsqp*upsi*(-2.)
       szz= szz+bet1*xscqs*abcsqm*rabescq*taui*(-8.)
       szz= szz+bet1*xscqs*abcsqm*rabescq*upsi*(-8.)
       szz= szz+bet1*xcmb1*abcsqp*abecsqp*rltau*(-2./9.)
       szz= szz+bet1*xcmb1*abcsqp*abecsqp*rlups*(-2./9.)
       szz= szz+bet1*xcmb1*abcsqp*rabescq*rltau*(-8./9.)
       szz= szz+bet1*xcmb1*abcsqp*rabescq*rlups*(-8./9.)
       szz= szz+bet1*xcmb1*abcsqpm*abecsqp*rltau*(2./3.)
       szz= szz+bet1*xcmb1*abcsqpm*abecsqp*rlups*(2./3.)
       szz= szz+bet1*xcmb1*abcsqpm*rabescq*rltau*(8./3.)
       szz= szz+bet1*xcmb1*abcsqpm*rabescq*rlups*(8./3.)
       szz= szz+bet1*xcmb1*abcsqmp*abecsqp*rltau*(2./3.)
       szz= szz+bet1*xcmb1*abcsqmp*abecsqp*rlups*(2./3.)
       szz= szz+bet1*xcmb1*abcsqmp*rabescq*rltau*(8./3.)
       szz= szz+bet1*xcmb1*abcsqmp*rabescq*rlups*(8./3.)
       szz= szz+bet1*xcmb1*abcsqm*abecsqp*rltau*(-2.)
       szz= szz+bet1*xcmb1*abcsqm*abecsqp*rlups*(-2.)
       szz= szz+bet1*xcmb1*abcsqm*rabescq*rltau*(-8.)
       szz= szz+bet1*xcmb1*abcsqm*rabescq*rlups*(-8.)
       szz= szz+bet1*xcmb3*abcsqp*abecsqp*taus*(1./3.)
       szz= szz+bet1*xcmb3*abcsqp*abecsqp*upss*(1./3.)
       szz= szz+bet1*xcmb3*abcsqp*abecsqp*one*(2./3.)
       szz= szz+bet1*xcmb3*abcsqp*rabescq*taus*(4./3.)
       szz= szz+bet1*xcmb3*abcsqp*rabescq*upss*(4./3.)
       szz= szz+bet1*xcmb3*abcsqp*rabescq*one*(8./3.)
       szz= szz+bet1*xcmb3*abcsqpm*abecsqp*taus*(-1.)
       szz= szz+bet1*xcmb3*abcsqpm*abecsqp*upss*(-1.)
       szz= szz+bet1*xcmb3*abcsqpm*abecsqp*one*(-2.)
       szz= szz+bet1*xcmb3*abcsqpm*rabescq*taus*(-4.)
       szz= szz+bet1*xcmb3*abcsqpm*rabescq*upss*(-4.)
       szz= szz+bet1*xcmb3*abcsqpm*rabescq*one*(-8.)
       szz= szz+bet1*xcmb3*abcsqmp*abecsqp*one*(-2./3.)
       szz= szz+bet1*xcmb3*abcsqmp*rabescq*one*(-8./3.)
       szz= szz+bet1*xcmb3*abcsqm*abecsqp*one*(2.)
       szz= szz+bet1*xcmb3*abcsqm*rabescq*one*(8.)
       szz= szz+bet1*xcmb4*abcsqp*abecsqp*taus*(-1./3.)
       szz= szz+bet1*xcmb4*abcsqp*abecsqp*upss*(-1./3.)
       szz= szz+bet1*xcmb4*abcsqp*rabescq*taus*(-4./3.)
       szz= szz+bet1*xcmb4*abcsqp*rabescq*upss*(-4./3.)
       szz= szz+bet1*xcmb4*abcsqpm*abecsqp*taus*(1)
       szz= szz+bet1*xcmb4*abcsqpm*abecsqp*upss*(1)
       szz= szz+bet1*xcmb4*abcsqpm*rabescq*taus*(4)
       szz= szz+bet1*xcmb4*abcsqpm*rabescq*upss*(4)
       szz= szz+bet1*xcmb4*abcsqmp*abecsqp*one*(-2./3.)
       szz= szz+bet1*xcmb4*abcsqmp*rabescq*one*(-8./3.)
       szz= szz+bet1*xcmb4*abcsqm*abecsqp*one*(2.)
       szz= szz+bet1*xcmb4*abcsqm*rabescq*one*(8.)
       szz= szz+bet1*abcsqp*abecsqp*taus*(2./9.*xsqs)
       szz= szz+bet1*abcsqp*abecsqp*upss*(2./9.*xsqs)
       szz= szz+bet1*abcsqp*rabescq*taus*(8./9.*xsqs)
       szz= szz+bet1*abcsqp*rabescq*upss*(8./9.*xsqs)
       szz= szz+bet1*abcsqpm*abecsqp*taus*(-2./3.*xsqs)
       szz= szz+bet1*abcsqpm*abecsqp*upss*(-2./3.*xsqs)
       szz= szz+bet1*abcsqpm*rabescq*taus*(-8./3.*xsqs)
       szz= szz+bet1*abcsqpm*rabescq*upss*(-8./3.*xsqs)
       szz= szz+bet2*sqq*xcmb2*abcsqp*abecsqp*rltau*(-16./9.)
       szz= szz+bet2*sqq*xcmb2*abcsqp*abecsqp*rlups*(-16./9.)
       szz= szz+bet2*sqq*xcmb2*abcsqp*rabescq*rltau*(-64./9.)
       szz= szz+bet2*sqq*xcmb2*abcsqp*rabescq*rlups*(-64./9.)
       szz= szz+bet2*sqq*xcmb2*abcsqmp*abecsqp*rltau*(16./3.)
       szz= szz+bet2*sqq*xcmb2*abcsqmp*abecsqp*rlups*(16./3.)
       szz= szz+bet2*sqq*xcmb2*abcsqmp*rabescq*rltau*(64./3.)
       szz= szz+bet2*sqq*xcmb2*abcsqmp*rabescq*rlups*(64./3.)
       szz= szz+bet2*xscqs*abcsqp*abecsqp*taui*(-8./9.)
       szz= szz+bet2*xscqs*abcsqp*abecsqp*upsi*(-8./9.)
       szz= szz+bet2*xscqs*abcsqp*rabescq*taui*(-32./9.)
       szz= szz+bet2*xscqs*abcsqp*rabescq*upsi*(-32./9.)
       szz= szz+bet2*xscqs*abcsqmp*abecsqp*taui*(8./3.)
       szz= szz+bet2*xscqs*abcsqmp*abecsqp*upsi*(8./3.)
       szz= szz+bet2*xscqs*abcsqmp*rabescq*taui*(32./3.)
       szz= szz+bet2*xscqs*abcsqmp*rabescq*upsi*(32./3.)
       szz= szz+bet2*xcmb1*abcsqp*abecsqp*rltau*(-8./9.)
       szz= szz+bet2*xcmb1*abcsqp*abecsqp*rlups*(-8./9.)
       szz= szz+bet2*xcmb1*abcsqp*rabescq*rltau*(-32./9.)
       szz= szz+bet2*xcmb1*abcsqp*rabescq*rlups*(-32./9.)
       szz= szz+bet2*xcmb1*abcsqmp*abecsqp*rltau*(8./3.)
       szz= szz+bet2*xcmb1*abcsqmp*abecsqp*rlups*(8./3.)
       szz= szz+bet2*xcmb1*abcsqmp*rabescq*rltau*(32./3.)
       szz= szz+bet2*xcmb1*abcsqmp*rabescq*rlups*(32./3.)
       szz= szz+bet2*xcmb3*abcsqp*abecsqp*taus*(4./3.)
       szz= szz+bet2*xcmb3*abcsqp*abecsqp*upss*(4./3.)
       szz= szz+bet2*xcmb3*abcsqp*abecsqp*one*(8./3.)
       szz= szz+bet2*xcmb3*abcsqp*rabescq*taus*(16./3.)
       szz= szz+bet2*xcmb3*abcsqp*rabescq*upss*(16./3.)
       szz= szz+bet2*xcmb3*abcsqp*rabescq*one*(32./3.)
       szz= szz+bet2*xcmb3*abcsqmp*abecsqp*one*(-8./3.)
       szz= szz+bet2*xcmb3*abcsqmp*rabescq*one*(-32./3.)
       szz= szz+bet2*xcmb4*abcsqp*abecsqp*taus*(-4./3.)
       szz= szz+bet2*xcmb4*abcsqp*abecsqp*upss*(-4./3.)
       szz= szz+bet2*xcmb4*abcsqp*rabescq*taus*(-16./3.)
       szz= szz+bet2*xcmb4*abcsqp*rabescq*upss*(-16./3.)
       szz= szz+bet2*xcmb4*abcsqmp*abecsqp*one*(-8./3.)
       szz= szz+bet2*xcmb4*abcsqmp*rabescq*one*(-32./3.)
       szz= szz+bet2*abcsqp*abecsqp*taus*(8./9.*xsqs)
       szz= szz+bet2*abcsqp*abecsqp*upss*(8./9.*xsqs)
       szz= szz+bet2*abcsqp*rabescq*taus*(32./9.*xsqs)
       szz= szz+bet2*abcsqp*rabescq*upss*(32./9.*xsqs)
       szz= szz+bet3*sqq*xcmb2*abcsqp*abecsqp*rltau*(-16./9.)
       szz= szz+bet3*sqq*xcmb2*abcsqp*abecsqp*rlups*(-16./9.)
       szz= szz+bet3*sqq*xcmb2*abcsqp*rabescq*rltau*(-64./9.)
       szz= szz+bet3*sqq*xcmb2*abcsqp*rabescq*rlups*(-64./9.)
       szz= szz+bet3*sqq*xcmb2*abcsqpm*abecsqp*rltau*(16./3.)
       szz= szz+bet3*sqq*xcmb2*abcsqpm*abecsqp*rlups*(16./3.)
       szz= szz+bet3*sqq*xcmb2*abcsqpm*rabescq*rltau*(64./3.)
       szz= szz+bet3*sqq*xcmb2*abcsqpm*rabescq*rlups*(64./3.)
       szz= szz+bet3*xscqs*abcsqp*abecsqp*taui*(-8./9.)
       szz= szz+bet3*xscqs*abcsqp*abecsqp*upsi*(-8./9.)
       szz= szz+bet3*xscqs*abcsqp*rabescq*taui*(-32./9.)
       szz= szz+bet3*xscqs*abcsqp*rabescq*upsi*(-32./9.)
       szz= szz+bet3*xscqs*abcsqpm*abecsqp*taui*(8./3.)
       szz= szz+bet3*xscqs*abcsqpm*abecsqp*upsi*(8./3.)
       szz= szz+bet3*xscqs*abcsqpm*rabescq*taui*(32./3.)
       szz= szz+bet3*xscqs*abcsqpm*rabescq*upsi*(32./3.)
       szz= szz+bet3*xcmb1*abcsqp*abecsqp*rltau*(-8./9.)
       szz= szz+bet3*xcmb1*abcsqp*abecsqp*rlups*(-8./9.)
       szz= szz+bet3*xcmb1*abcsqp*rabescq*rltau*(-32./9.)
       szz= szz+bet3*xcmb1*abcsqp*rabescq*rlups*(-32./9.)
       szz= szz+bet3*xcmb1*abcsqpm*abecsqp*rltau*(8./3.)
       szz= szz+bet3*xcmb1*abcsqpm*abecsqp*rlups*(8./3.)
       szz= szz+bet3*xcmb1*abcsqpm*rabescq*rltau*(32./3.)
       szz= szz+bet3*xcmb1*abcsqpm*rabescq*rlups*(32./3.)
       szz= szz+bet3*xcmb3*abcsqp*abecsqp*taus*(-2./3.)
       szz= szz+bet3*xcmb3*abcsqp*abecsqp*upss*(-2./3.)
       szz= szz+bet3*xcmb3*abcsqp*rabescq*taus*(-8./3.)
       szz= szz+bet3*xcmb3*abcsqp*rabescq*upss*(-8./3.)
       szz= szz+bet3*xcmb3*abcsqpm*abecsqp*taus*(2.)
       szz= szz+bet3*xcmb3*abcsqpm*abecsqp*upss*(2.)
       szz= szz+bet3*xcmb3*abcsqpm*rabescq*taus*(8.)
       szz= szz+bet3*xcmb3*abcsqpm*rabescq*upss*(8.)
       szz= szz+bet3*xcmb4*abcsqp*abecsqp*taus*(2./3.)
       szz= szz+bet3*xcmb4*abcsqp*abecsqp*upss*(2./3.)
       szz= szz+bet3*xcmb4*abcsqp*abecsqp*one*(4./3.)
       szz= szz+bet3*xcmb4*abcsqp*rabescq*taus*(8./3.)
       szz= szz+bet3*xcmb4*abcsqp*rabescq*upss*(8./3.)
       szz= szz+bet3*xcmb4*abcsqp*rabescq*one*(16./3.)
       szz= szz+bet3*xcmb4*abcsqpm*abecsqp*taus*(-2.)
       szz= szz+bet3*xcmb4*abcsqpm*abecsqp*upss*(-2.)
       szz= szz+bet3*xcmb4*abcsqpm*abecsqp*one*(-4.)
       szz= szz+bet3*xcmb4*abcsqpm*rabescq*taus*(-8.)
       szz= szz+bet3*xcmb4*abcsqpm*rabescq*upss*(-8.)
       szz= szz+bet3*xcmb4*abcsqpm*rabescq*one*(-16.)
       szz= szz+bet3*abcsqp*abecsqp*taus*(-4./9.*xsqs)
       szz= szz+bet3*abcsqp*abecsqp*upss*(-4./9.*xsqs)
       szz= szz+bet3*abcsqp*rabescq*taus*(-16./9.*xsqs)
       szz= szz+bet3*abcsqp*rabescq*upss*(-16./9.*xsqs)
       szz= szz+bet3*abcsqpm*abecsqp*taus*(4./3.*xsqs)
       szz= szz+bet3*abcsqpm*abecsqp*upss*(4./3.*xsqs)
       szz= szz+bet3*abcsqpm*rabescq*taus*(16./3.*xsqs)
       szz= szz+bet3*abcsqpm*rabescq*upss*(16./3.*xsqs)
       szz= szz+sqq*xcmb2*abcsqp*abecsqp*rltau*(-64./9.)
       szz= szz+sqq*xcmb2*abcsqp*abecsqp*rlups*(-64./9.)
       szz= szz+sqq*xcmb2*abcsqp*rabescq*rltau*(-256./9.)
       szz= szz+sqq*xcmb2*abcsqp*rabescq*rlups*(-256./9.)
       szz= szz+xscqs*abcsqp*abecsqp*taui*(-32./9.)
       szz= szz+xscqs*abcsqp*abecsqp*upsi*(-32./9.)
       szz= szz+xscqs*abcsqp*rabescq*taui*(-128./9.)
       szz= szz+xscqs*abcsqp*rabescq*upsi*(-128./9.)
       szz= szz+xcmb1*abcsqp*abecsqp*rltau*(-32./9.)
       szz= szz+xcmb1*abcsqp*abecsqp*rlups*(-32./9.)
       szz= szz+xcmb1*abcsqp*rabescq*rltau*(-128./9.)
       szz= szz+xcmb1*abcsqp*rabescq*rlups*(-128./9.)
       szz= szz+xcmb3*abcsqp*abecsqp*taus*(-8./3.)
       szz= szz+xcmb3*abcsqp*abecsqp*upss*(-8./3.)
       szz= szz+xcmb3*abcsqp*rabescq*taus*(-32./3.)
       szz= szz+xcmb3*abcsqp*rabescq*upss*(-32./3.)
       szz= szz+xcmb4*abcsqp*abecsqp*taus*(8./3.)
       szz= szz+xcmb4*abcsqp*abecsqp*upss*(8./3.)
       szz= szz+xcmb4*abcsqp*abecsqp*one*(16./3.)
       szz= szz+xcmb4*abcsqp*rabescq*taus*(32./3.)
       szz= szz+xcmb4*abcsqp*rabescq*upss*(32./3.)
       szz= szz+xcmb4*abcsqp*rabescq*one*(64./3.)
       szz= szz+abcsqp*abecsqp*taus*(-16./9.*xsqs)
       szz= szz+abcsqp*abecsqp*upss*(-16./9.*xsqs)
       szz= szz+abcsqp*rabescq*taus*(-64./9.*xsqs)
       szz= szz+abcsqp*rabescq*upss*(-64./9.*xsqs)
       if(oqcdf.eq.'y'.and.i.gt.4) then
*
      qcdcq= betvcsq*xscqs*abecsqp*taui*rvcq*acsqs*(-32./9.)
      qcdcq= qcdcq+betvcsq*xscqs*abecsqp*upsi*rvcq*acsqs
     # *(-32./9.)
      qcdcq= qcdcq+betvcsq*xscqs*rabescq*taui*rvcq*acsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcsq*xscqs*rabescq*upsi*rvcq*acsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcsq*xcmb1*abecsqp*rltau*rvcq*acsqs
     # *(-32./9.)
      qcdcq= qcdcq+betvcsq*xcmb1*abecsqp*rlups*rvcq*acsqs
     # *(-32./9.)
      qcdcq= qcdcq+betvcsq*xcmb1*rabescq*rltau*rvcq*acsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcsq*xcmb1*rabescq*rlups*rvcq*acsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcsq*xcmb2*abecsqp*rltau*rvcq*sqq*acqs*
     #  asqs*(-64./9.)
      qcdcq= qcdcq+betvcsq*xcmb2*abecsqp*rlups*rvcq*sqq*acqs*
     #  asqs*(-64./9.)
      qcdcq= qcdcq+betvcsq*xcmb2*rabescq*rltau*rvcq*sqq*acqs*
     #  asqs*(-256./9.)
      qcdcq= qcdcq+betvcsq*xcmb2*rabescq*rlups*rvcq*sqq*acqs*
     #  asqs*(-256./9.)
      qcdcq= qcdcq+betvcsq*xcmb3*abecsqp*one*rvcq*acsqs
     # *(32./9.)
      qcdcq= qcdcq+betvcsq*xcmb3*rabescq*one*rvcq*acsqs
     # *(128./9.)
      qcdcq= qcdcq+betvcsq*xcmb4*abecsqp*one*rvcq*acsqs
     # *(32./9.)
      qcdcq= qcdcq+betvcsq*xcmb4*rabescq*one*rvcq*acsqs
     # *(128./9.)
      qcdcq= qcdcq+betvcq*xscqs*abecsqp*taui*rvcq*abcsqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+betvcq*xscqs*abecsqp*upsi*rvcq*abcsqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+betvcq*xscqs*rabescq*taui*rvcq*abcsqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcq*xscqs*rabescq*upsi*rvcq*abcsqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcq*xcmb1*abecsqp*rltau*rvcq*abcsqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+betvcq*xcmb1*abecsqp*rlups*rvcq*abcsqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+betvcq*xcmb1*rabescq*rltau*rvcq*abcsqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcq*xcmb1*rabescq*rlups*rvcq*abcsqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+betvcq*xcmb2*abecsqp*rltau*rvcq*sqq*abcsqs*
     #  betsqs*(-64./9.)
      qcdcq= qcdcq+betvcq*xcmb2*abecsqp*rlups*rvcq*sqq*abcsqs*
     #  betsqs*(-64./9.)
      qcdcq= qcdcq+betvcq*xcmb2*rabescq*rltau*rvcq*sqq*abcsqs*
     #  betsqs*(-256./9.)
      qcdcq= qcdcq+betvcq*xcmb2*rabescq*rlups*rvcq*sqq*abcsqs*
     #  betsqs*(-256./9.)
      qcdcq= qcdcq+betvcq*xcmb3*abecsqp*one*rvcq*abcsqs*betsqs
     # *(32./9.)
      qcdcq= qcdcq+betvcq*xcmb3*rabescq*one*rvcq*abcsqs*betsqs
     # *(128./9.)
      qcdcq= qcdcq+betvcq*xcmb4*abecsqp*one*rvcq*abcsqs*betsqs
     # *(32./9.)
      qcdcq= qcdcq+betvcq*xcmb4*rabescq*one*rvcq*abcsqs*betsqs
     # *(128./9.)
      qcdcq= qcdcq+betvsq*xscqs*abecsqp*taui*racq*bacsqs*betcqs
     # *(-32./9.)
      qcdcq= qcdcq+betvsq*xscqs*abecsqp*upsi*racq*bacsqs*betcqs
     # *(-32./9.)
      qcdcq= qcdcq+betvsq*xscqs*rabescq*taui*racq*bacsqs*betcqs
     # *(-128./9.)
      qcdcq= qcdcq+betvsq*xscqs*rabescq*upsi*racq*bacsqs*betcqs
     # *(-128./9.)
      qcdcq= qcdcq+betvsq*xcmb1*abecsqp*rltau*racq*bacsqs*betcqs
     # *(-32./9.)
      qcdcq= qcdcq+betvsq*xcmb1*abecsqp*rlups*racq*bacsqs*betcqs
     # *(-32./9.)
      qcdcq= qcdcq+betvsq*xcmb1*rabescq*rltau*racq*bacsqs*betcqs
     # *(-128./9.)
      qcdcq= qcdcq+betvsq*xcmb1*rabescq*rlups*racq*bacsqs*betcqs
     # *(-128./9.)
      qcdcq= qcdcq+betvsq*xcmb2*abecsqp*rltau*racq*sqq*bacsqs*
     #  betcqs*(-64./9.)
      qcdcq= qcdcq+betvsq*xcmb2*abecsqp*rlups*racq*sqq*bacsqs*
     #  betcqs*(-64./9.)
      qcdcq= qcdcq+betvsq*xcmb2*rabescq*rltau*racq*sqq*bacsqs*
     #  betcqs*(-256./9.)
      qcdcq= qcdcq+betvsq*xcmb2*rabescq*rlups*racq*sqq*bacsqs*
     #  betcqs*(-256./9.)
      qcdcq= qcdcq+betvsq*xcmb3*abecsqp*one*racq*bacsqs*betcqs
     # *(32./9.)
      qcdcq= qcdcq+betvsq*xcmb3*rabescq*one*racq*bacsqs*betcqs
     # *(128./9.)
      qcdcq= qcdcq+betvsq*xcmb4*abecsqp*one*racq*bacsqs*betcqs
     # *(32./9.)
      qcdcq= qcdcq+betvsq*xcmb4*rabescq*one*racq*bacsqs*betcqs
     # *(128./9.)
      qcdcq= qcdcq+xscqs*abecsqp*taui*racq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+xscqs*abecsqp*upsi*racq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+xscqs*rabescq*taui*racq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+xscqs*rabescq*upsi*racq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+xcmb1*abecsqp*rltau*racq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+xcmb1*abecsqp*rlups*racq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdcq= qcdcq+xcmb1*rabescq*rltau*racq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+xcmb1*rabescq*rlups*racq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdcq= qcdcq+xcmb2*abecsqp*rltau*racq*sqq*bcsqs*betcqs*
     #  betsqs*(-64./9.)
      qcdcq= qcdcq+xcmb2*abecsqp*rlups*racq*sqq*bcsqs*betcqs*
     #  betsqs*(-64./9.)
      qcdcq= qcdcq+xcmb2*rabescq*rltau*racq*sqq*bcsqs*betcqs*
     #  betsqs*(-256./9.)
      qcdcq= qcdcq+xcmb2*rabescq*rlups*racq*sqq*bcsqs*betcqs*
     #  betsqs*(-256./9.)
      qcdcq= qcdcq+xcmb3*abecsqp*one*racq*bcsqs*betcqs*betsqs
     # *(32./9.)
      qcdcq= qcdcq+xcmb3*rabescq*one*racq*bcsqs*betcqs*betsqs
     # *(128./9.)
      qcdcq= qcdcq+xcmb4*abecsqp*one*racq*bcsqs*betcqs*betsqs
     # *(32./9.)
      qcdcq= qcdcq+xcmb4*rabescq*one*racq*bcsqs*betcqs*betsqs
     # *(128./9.)
*
       else
        qcdcq= 0.d0
       endif
       if(oqcdf.eq.'y'.and.j.gt.4) then
*
      qcdsq= betvcsq*xscqs*abecsqp*taui*rvsq*acsqs*(-32./9.)
      qcdsq= qcdsq+betvcsq*xscqs*abecsqp*upsi*rvsq*acsqs
     # *(-32./9.)
      qcdsq= qcdsq+betvcsq*xscqs*rabescq*taui*rvsq*acsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcsq*xscqs*rabescq*upsi*rvsq*acsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcsq*xcmb1*abecsqp*rltau*rvsq*acsqs
     # *(-32./9.)
      qcdsq= qcdsq+betvcsq*xcmb1*abecsqp*rlups*rvsq*acsqs
     # *(-32./9.)
      qcdsq= qcdsq+betvcsq*xcmb1*rabescq*rltau*rvsq*acsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcsq*xcmb1*rabescq*rlups*rvsq*acsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcsq*xcmb2*abecsqp*rltau*rvsq*sqq*acqs*
     #  asqs*(-64./9.)
      qcdsq= qcdsq+betvcsq*xcmb2*abecsqp*rlups*rvsq*sqq*acqs*
     #  asqs*(-64./9.)
      qcdsq= qcdsq+betvcsq*xcmb2*rabescq*rltau*rvsq*sqq*acqs*
     #  asqs*(-256./9.)
      qcdsq= qcdsq+betvcsq*xcmb2*rabescq*rlups*rvsq*sqq*acqs*
     #  asqs*(-256./9.)
      qcdsq= qcdsq+betvcsq*xcmb3*abecsqp*one*rvsq*acsqs
     # *(32./9.)
      qcdsq= qcdsq+betvcsq*xcmb3*rabescq*one*rvsq*acsqs
     # *(128./9.)
      qcdsq= qcdsq+betvcsq*xcmb4*abecsqp*one*rvsq*acsqs
     # *(32./9.)
      qcdsq= qcdsq+betvcsq*xcmb4*rabescq*one*rvsq*acsqs
     # *(128./9.)
      qcdsq= qcdsq+betvcq*xscqs*abecsqp*taui*rasq*abcsqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+betvcq*xscqs*abecsqp*upsi*rasq*abcsqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+betvcq*xscqs*rabescq*taui*rasq*abcsqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcq*xscqs*rabescq*upsi*rasq*abcsqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcq*xcmb1*abecsqp*rltau*rasq*abcsqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+betvcq*xcmb1*abecsqp*rlups*rasq*abcsqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+betvcq*xcmb1*rabescq*rltau*rasq*abcsqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcq*xcmb1*rabescq*rlups*rasq*abcsqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+betvcq*xcmb2*abecsqp*rltau*rasq*sqq*abcsqs*
     #  betsqs*(-64./9.)
      qcdsq= qcdsq+betvcq*xcmb2*abecsqp*rlups*rasq*sqq*abcsqs*
     #  betsqs*(-64./9.)
      qcdsq= qcdsq+betvcq*xcmb2*rabescq*rltau*rasq*sqq*abcsqs*
     #  betsqs*(-256./9.)
      qcdsq= qcdsq+betvcq*xcmb2*rabescq*rlups*rasq*sqq*abcsqs*
     #  betsqs*(-256./9.)
      qcdsq= qcdsq+betvcq*xcmb3*abecsqp*one*rasq*abcsqs*betsqs
     # *(32./9.)
      qcdsq= qcdsq+betvcq*xcmb3*rabescq*one*rasq*abcsqs*betsqs
     # *(128./9.)
      qcdsq= qcdsq+betvcq*xcmb4*abecsqp*one*rasq*abcsqs*betsqs
     # *(32./9.)
      qcdsq= qcdsq+betvcq*xcmb4*rabescq*one*rasq*abcsqs*betsqs
     # *(128./9.)
      qcdsq= qcdsq+betvsq*xscqs*abecsqp*taui*rvsq*bacsqs*betcqs
     # *(-32./9.)
      qcdsq= qcdsq+betvsq*xscqs*abecsqp*upsi*rvsq*bacsqs*betcqs
     # *(-32./9.)
      qcdsq= qcdsq+betvsq*xscqs*rabescq*taui*rvsq*bacsqs*betcqs
     # *(-128./9.)
      qcdsq= qcdsq+betvsq*xscqs*rabescq*upsi*rvsq*bacsqs*betcqs
     # *(-128./9.)
      qcdsq= qcdsq+betvsq*xcmb1*abecsqp*rltau*rvsq*bacsqs*betcqs
     # *(-32./9.)
      qcdsq= qcdsq+betvsq*xcmb1*abecsqp*rlups*rvsq*bacsqs*betcqs
     # *(-32./9.)
      qcdsq= qcdsq+betvsq*xcmb1*rabescq*rltau*rvsq*bacsqs*betcqs
     # *(-128./9.)
      qcdsq= qcdsq+betvsq*xcmb1*rabescq*rlups*rvsq*bacsqs*betcqs
     # *(-128./9.)
      qcdsq= qcdsq+betvsq*xcmb2*abecsqp*rltau*rvsq*sqq*bacsqs*
     #  betcqs*(-64./9.)
      qcdsq= qcdsq+betvsq*xcmb2*abecsqp*rlups*rvsq*sqq*bacsqs*
     #  betcqs*(-64./9.)
      qcdsq= qcdsq+betvsq*xcmb2*rabescq*rltau*rvsq*sqq*bacsqs*
     #  betcqs*(-256./9.)
      qcdsq= qcdsq+betvsq*xcmb2*rabescq*rlups*rvsq*sqq*bacsqs*
     #  betcqs*(-256./9.)
      qcdsq= qcdsq+betvsq*xcmb3*abecsqp*one*rvsq*bacsqs*betcqs
     # *(32./9.)
      qcdsq= qcdsq+betvsq*xcmb3*rabescq*one*rvsq*bacsqs*betcqs
     # *(128./9.)
      qcdsq= qcdsq+betvsq*xcmb4*abecsqp*one*rvsq*bacsqs*betcqs
     # *(32./9.)
      qcdsq= qcdsq+betvsq*xcmb4*rabescq*one*rvsq*bacsqs*betcqs
     # *(128./9.)
      qcdsq= qcdsq+xscqs*abecsqp*taui*rasq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+xscqs*abecsqp*upsi*rasq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+xscqs*rabescq*taui*rasq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+xscqs*rabescq*upsi*rasq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+xcmb1*abecsqp*rltau*rasq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+xcmb1*abecsqp*rlups*rasq*bcsqs*betcqs*betsqs
     # *(-32./9.)
      qcdsq= qcdsq+xcmb1*rabescq*rltau*rasq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+xcmb1*rabescq*rlups*rasq*bcsqs*betcqs*betsqs
     # *(-128./9.)
      qcdsq= qcdsq+xcmb2*abecsqp*rltau*rasq*sqq*bcsqs*betcqs*
     #  betsqs*(-64./9.)
      qcdsq= qcdsq+xcmb2*abecsqp*rlups*rasq*sqq*bcsqs*betcqs*
     #  betsqs*(-64./9.)
      qcdsq= qcdsq+xcmb2*rabescq*rltau*rasq*sqq*bcsqs*betcqs*
     #  betsqs*(-256./9.)
      qcdsq= qcdsq+xcmb2*rabescq*rlups*rasq*sqq*bcsqs*betcqs*
     #  betsqs*(-256./9.)
      qcdsq= qcdsq+xcmb3*abecsqp*one*rasq*bcsqs*betcqs*betsqs
     # *(32./9.)
      qcdsq= qcdsq+xcmb3*rabescq*one*rasq*bcsqs*betcqs*betsqs
     # *(128./9.)
      qcdsq= qcdsq+xcmb4*abecsqp*one*rasq*bcsqs*betcqs*betsqs
     # *(32./9.)
      qcdsq= qcdsq+xcmb4*rabescq*one*rasq*bcsqs*betcqs*betsqs
     # *(128./9.)
*
       else
        qcdsq= 0.d0
       endif
       aszz(i,j)= wh(i,j)*betcq*betsq*(szz+qcdcq+qcdsq)
*
       if(i.eq.j.and.i.ne.1) then
        aszz(i,j)= aszz(i,j)/2.d0
       endif
*
       enddo
      enddo
*
      szzt= 0.d0
      if(ochan.eq.'tota') then
       do i=1,9
        do j=i,9
         szzt= szzt+aszz(i,j)
        enddo
       enddo
      else if(ochan.eq.'qqqq') then
       do i=5,9
        do j=i,9
         szzt= szzt+aszz(i,j)
        enddo
       enddo
      else if(ochan.eq.'llll') then
       do i=2,4
        do j=i,4
         szzt= szzt+aszz(i,j)
        enddo
       enddo
      else if(ochan.eq.'qqnn') then
       do i=1,1
        do j=5,9
         szzt= szzt+aszz(i,j)
        enddo
       enddo
      else if(ochan.eq.'qqll') then
       do i=2,4
        do j=5,9
         szzt= szzt+aszz(i,j)
        enddo
       enddo
      else if(ochan.eq.'llnn') then
       do i=1,1
        do j=2,4
         szzt= szzt+aszz(i,j)
        enddo
       enddo
      else if(ochan.eq.'nnnn') then
       do i=1,1
        do j=1,1
         szzt= szzt+aszz(i,j)
        enddo
       enddo
      endif
*
      if(ofl.eq.'n') then
       factc= g8/cth8
      else if(ofl.eq.'y') then
       factc= (grcq*grcq+gicq*gicq)*(grsq*grsq+gisq*gisq)
       factc= factc/(ctrrcq*ctrrcq+ctricq*ctricq)/
     #       (ctrrsq*ctrrsq+ctrisq*ctrisq)
       factc= factc/brhmcq/brhmsq
      endif
      if(oqedf.eq.'y') then
       factc= factc*qedcq*qedsq
      endif
      fact= factc*vjc*ujc*xcqjac*xsqjac*stf*vv/s
      factn= 1.d0/256.d0*pis*pi/2.d0/4.d0/2.d0/
     #       (2.d0*pi)**8*0.38937966d9
      zztoborn= fact*factn*szzt
*
      return
      end