*
*-------------------------------------------------------------------------
*
      real*8 function wtoxsswcfl(ndim,x)
      implicit real*8(a-h,o-z)
      character*1 opeak,ofl,rio,ospe
*
      common/wtai/rio
      common/wtfls/ofl
      common/wtcat/ipro
      common/wtps/opeak
      common/wtkount/ik
      common/wtswfl/ospe
      common/wtistrf/isf
      common/wtipt/ifz(51)
      common/wtfsw/tfactsw
      common/wtmult/inm,jnm
      common/wtnpr/ipr,ipr0
      common/wticuts/iac(4)
      common/wtcutsw/tcs,s0cut
      common/wtcmplx/sw(2),sz(2),rsw(2),rsz(2)
      common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy
      common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi
      common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2,
     #              tqm2,dmy2
      common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,g8,tfact,pih,alw,
     #              eta,feta,beta,g2,tfacth
      common/wtapar/ars,s,rwm,rwm2,rwg,rwmg,swg,swgs,opswgs,sth2,cth2,
     #            hsth2,tsth2,scth2,asth2,tth2,rzm,rzm2,rzg,rzmg,szg,
     #            szgs,opszgs,sth4,cth4,ve,vf,vfp,rbqm2,rszw,rszw2,
     #            s0w,s0z
      common/wtcuts/aim(6),bim(6),ae(4),asa(4),bsa(4),afsa(6),bfsa(6),
     #            ombsa(4),opbsa(4),teq,rae(4),omasa(4),opasa(4),
     #            sg12,cg12,sg13,cg13,sg14,cg14,sg23,cg23,sg24,
     #            cg24,sg34,cg34,sct120,sct130,sct140,sct230,
     #            sct240,sct340,sgam(4),cgam(4)
      common/wtampli/gsm2,gtm2,rl1,rl2,rl3,xv,x1,x2,z,tau,ups,strrx,
     #               ctrrx,ums,dms,vv,y,pwsm2,gstm2,zxi,pwtm2             
      common/wtamplia/rs1wg(14),rs2wg(14),rt1wg(14),rt2wg(14),
     #                rwwg(14,14),rs1wz(14),rs2wz(14),rt1wz(14),
     #                rt2wz(14),rwwz(14,14),rwwgz(14,14)
      common/wtampli2/vvs,vvc,vvq,ys,yc,yq,bxp,bxps,bxpc,bxpi,bxx,o,
     #                os,oc,oq,x12,x12s,x12c,x2z,x2zs,x2zc,x2zq,zi,
     #                zis,tx,txs,ux,uxs,cst,cqt,csti,cqti,vpe,vme,
     #                vpd,vmd,vpu,vmu,vpes,vmes,vpds,vmds,vpus,vmus
*
      dimension co(2)
      dimension vu(4),vut(4)
      dimension vd(4),vdt(4)
      dimension vn(4),vnt(4)
      dimension x(ndim),ww(2)
      dimension cp_11(2),cp_12(2)
      dimension cm_11(2),cm_12(2)
      dimension cph_11(2),cph_12(2)
      dimension cmh_11(2),cmh_12(2)
      dimension c03(2,6),c02(2,4),c01(2,2),c00(2)
      dimension ct3(2,6),ct2(2,4),ct1(2,2),ct0(2)
      dimension cp_21(2),cp_22(2),cp_23(2),cp_24(2)
      dimension cm_21(2),cm_22(2),cm_23(2),cm_24(2)
      dimension cp_31(2),cp_32(2),cp_33(2),cp_34(2),
     #          cp_35(2),cp_36(2)
      dimension cm_31(2),cm_32(2),cm_33(2),cm_34(2),
     #          cm_35(2),cm_36(2)
      dimension ctt3(2,6),ctt2(2,4),ctt1(2,2),ctt0(2)
      dimension cph_21(2),cph_22(2),cph_23(2),cph_24(2)
      dimension cmh_21(2),cmh_22(2),cmh_23(2),cmh_24(2)
      dimension cph_31(2),cph_32(2),cph_33(2),cph_34(2),
     #          cph_35(2),cph_36(2)
      dimension cmh_31(2),cmh_32(2),cmh_33(2),cmh_34(2),
     #          cmh_35(2),cmh_36(2)
      dimension w(10,2),dw(10,2),wg(14,2),wz(14,2)
      dimension p3qw(2),fzw(2),fw(2),pggfzm(2),pgglqzm(2),
     #          pggnpzm(2),p3qsw(2),fzsw(2),fwsw(2),
     #          p3qsz(2),fzsz(2),fwsz(2),pggfsz(2),
     #          pgglqsz(2),pggnpsz(2)
      dimension pggfx(2),pgglqx(2),pggnpx(2),p3qx(2),fzx(2),fwx(2),
     #          p3qx2(2),fzx2(2),fwsx2(2),p3qy(2),fzy(2),fwsy(2)
*
      external s09aaf
*
*-----order of integration is: uv,vv; y; x2,x1,tau/ups,z,phi_u/phi_d
*
      ik= ik+1
      inull= 0
      taujc= 0.d0
      upsjc= 0.d0
      if(ofl.ne.'c'.and.rio.ne.'a'.and.ospe.ne.'y') then
       print*,'I/O parameter error in FL'
       stop
      endif
      co(1)= 1.d0
      co(2)= 0.d0
      swis= sw(2)*sw(2)
      srwis= rsw(2)*rsw(2)
      srzis= rsz(2)*rsz(2)
      aswi= abs(rsw(2))
      wm2= wm*wm
      zm2= zm*zm
      ali= alphai
      call wtopolez0(pggf0,fw0)
*
      p2r= -wm2
      p2i= 0.d0
      call wtopole(p2r,p2i,p3qw,fzw,fw)
*
      p2r= -zm2
      p2i= 0.d0
      call wtopoleg(p2r,p2i,pggfzm,pgglqzm,pggnpzm)
      ccz= gf*zm2/pis
      ccw= gf*wm2/pis
      apis= 16.d0*pis
      arezm= ali-0.25d0*(pggfzm(1)-pggf0+pggnpzm(1))/pi
      aizm= -0.25d0*(pggfzm(2)+pgglqzm(2))/pi
*
*-----running g at sw
*
      p2r= -sw(1)
      p2i= -sw(2)
      call wtopole(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
*
*-----running g at sz
*
      p2r= -sz(1)
      p2i= -sz(2)
      call wtopole(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
*
*-----running sine at sz
*
      call wtopoleg(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
*
      if(ipr.eq.4) then
         ipro= 1
      else if(ipr.eq.5) then
         ipro= 2
      else if(ipr.eq.39) then
         ipro= 3
      endif
      wwt= 0.d0
      wwp= 0.d0
      tjac= 0.d0
      do kk=1,2
        ww(kk)= 0.d0
      enddo
*
      em2= em*em
      em4= em2*em2
      rem2= em2/s
      rem= em/sqrt(s)
      rwms= rwm2
      rwmq= rwms*rwms
      rwgs= rwg*rwg
      rwgm= rwg*rwm
      rwgms= rwgs*rwms
      rwmgs= rwmg*rwmg
      rwmpgs= rwms+rwgs
      if(ipro.eq.2) then
       ums= uqm*uqm/s
       dms= dqm*dqm/s
      else if(ipro.eq.1) then
       ums= 0.d0
       dms= rmm*rmm/s
      else if(ipro.eq.3) then
       ums= 0.d0
       dms= tm*tm/s
      endif
      um= sqrt(ums)
      dm= sqrt(dms)
      s0= s0cut
*
      if(ndim.eq.7) then
       yx= x(1)
       xx= x(2)
       x2x= x(3)
       x1x= x(4)
       angx= x(5)
       zx= x(6)
       phx= x(7)
      else if(ndim.eq.9) then 
       uvx= x(1)
       vvx= x(2)
       yx= x(3)
       xx= x(4)
       x2x= x(5)
       x1x= x(6)
       angx= x(7)
       zx= x(8)
       phx= x(9)
      endif
      omangx= 1.d0-angx
*
      iz= 1
*
      if(ndim.eq.7) 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
       vcut= s0
       vvl= 0.25d0*(rem+sqrt(rem*rem+4.d0*vcut))**2
       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(3)= ifz(3)+1
         go to 1
      endif
      svv= sqrt(vv)
      vvs= vv*vv
      vvc= vvs*vv
      vvq= vvc*vv
      vvi= 1.d0/vv
      rdms= dms*vv
      rums= ums*vv
      xm= uv
      xms= xm*xm
      xp= vv/uv  
      xps= xp*xp
      xmop= xm/xp
      sh= vv*s
      ssh= sqrt(sh)
      vwmg= aswi*vv
*
      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
*
      yl1= s0/vv
      yl2= (um+dm)*(um+dm)/vv
      yl= dmax1(yl1,yl2)
      yu= 1.d0-em/ssh
      if(yl.gt.yu) then
         iz= 0
         ifz(4)= ifz(4)+1
         go to 1
      endif
      yjc= yu-yl
      y= yjc*yx+yl
      ys= y*y
      yc= ys*y
      yq= yc*y
      yvv= y*vv
      yvvs= yvv*yvv
      omy= 1.d0-y
      omys= omy*omy
      omyc= omys*omy
      dmsh= dms/yvv
      umsh= ums/yvv
      dmh= sqrt(dmsh)
      umh= sqrt(umsh)
      umshs= umsh*umsh
      dmshs= dmsh*dmsh
      umshc= umshs*umsh
      dmshc= dmshs*dmsh
      umshq= umshc*umsh
      dmshq= dmshc*dmsh
      umshf= umshq*umsh
*
*-----X is initialized
*
*      qh0s= em2*(xm*y+omxm)*y/omy
      ctc= cos(sqrt(tcs))
      ctcs= ctc*ctc
      betaes= 1.d0-4.d0*rem2
      betae= sqrt(betaes)
      argq= omys-4.d0*(1.d0-betaes*ctcs)/(1.d0+betaes)**2*rem2
      argq= sqrt(argq)
      argq= omy+argq
      chiq= 0.5d0*(1.d0+betaes)/(1.d0+betae*ctc)*argq*sqrt(s)
      qhcs= -2.d0*em2+0.5d0*sqrt(s)/chiq*((1.d0-betae*ctc)*chiq*chiq+
     #      (1.d0+betae*ctc)*em2)
      bdh= -xp*xm+(xp+xm)*rem2
      ah= xm-(1.d0-xm)*rem2
      bh= xm-1.d0
      ch= rem2*(1.d0-xp)/bdh
      dh= -xp/bdh+rem2*(1.d0+xp)/bdh
      eh= (1.d0-xm)/bdh*(xp-rem2)
      qhcs= (ah-rem2*bh*ch/dh)*qhcs+em2*bh/dh*(y-eh)
*
      argq0= omys-4.d0*(1.d0-betaes)/(1.d0+betaes)**2*rem2
      argq0= sqrt(argq0)
      argq0= omy+argq0
      chiq0= 0.5d0*(1.d0+betaes)/(1.d0+betae)*argq0*sqrt(s)
      qh0s= -2.d0*em2+0.5d0*sqrt(s)/chiq0*((1.d0-betae)*chiq0*chiq0+
     #      (1.d0+betae)*em2)
      qh0s= (ah-rem2*bh*ch/dh)*qh0s+em2*bh/dh*(y-eh)

*
      x0= qh0s/y/sh
      xc0= qhcs/y/sh
      xc1= omy/y
      xc= dmin1(xc0,xc1)
      dxc= xc-x0
      x0s= x0*x0
      xcs= xc*xc
*
      if(x0.gt.xc) then
         iz= 0
         ifz(5)= ifz(5)+1
         go to 1
      endif
*
*-----ii=1,4; 1 me^4/X^3, 2 me^2/X^2, 3 1/X, 4 rest
*
      ii= inm
      if(ii.eq.1) then
       xjac= 0.5d0*(1/x0s-1/xcs)
       xv= xjac*xx-0.5d0/x0s
       xv= sqrt(-0.5d0/xv)
       xjac= xjac/(y*sh)**2
      else if(ii.eq.2) then
       xjac= dxc/xc/x0/y/sh
       xv= x0*xc/(dxc*xx+x0)
      else if(ii.eq.3) then
       xjac= log(xc/x0)
       xv= exp(xjac*xx+log(x0))
      else if(ii.eq.4) then
       xjac= dxc
       xv= xjac*xx+x0
      endif
      omxv= 1.d0-xv
      opxv= 1.d0+xv
      opxvs= opxv*opxv
      xvs= xv*xv
      xpo= xv+1.d0
      xpos= xpo*xpo
*
*-----Quantities are initialized in the X-channel
*
      p2r= -xv*y*sh
      p2i= 0.d0
      call wtopole(p2r,p2i,p3qx,fzx,fwx)
      xgx0= 1.d0/8.d0/gf/wm2
      xgxr= 1.d0+0.5d0*ccw*((fw0-fw(1))/wm2+
     #      p3qw(1)-p3qx(1))
      xgxi= -0.5d0*ccw*p3qx(2)
      agrx= xgx0*xgxr
      agix= xgx0*xgxi
      gmx= agrx*agrx+agix*agix
      grx= agrx/gmx
      gix= -agix/gmx
      call wtopoleg(p2r,p2i,pggfx,pgglqx,pggnpx)
      arex= arezm+0.25d0*(pggfzm(1)+pgglqzm(1)-
     #      pggfx(1)-pgglqx(1))/pi
      aix= aizm+0.25d0*(pggfzm(2)+pgglqzm(2)-
     #     pggfx(2)-pgglqx(2))/pi
      emx= arex*arex+aix*aix
      erx= 4.d0*pi*arex/emx
      eix= -4.d0*pi*aix/emx
      strrx= erx*agrx-eix*agix
      strix= erx*agix+eix*agrx
*
      ctrrx= 1.d0-strrx
      ctrix= -strix
      ctrrsz= 1.d0-strrsz
      ctrisz= -strisz
      ctrmx= ctrrx*ctrrx+ctrix*ctrix
      ctrrxi= ctrrx/ctrmx
      ctrixi= -ctrix/ctrmx
      ttrrx= strrx*ctrrxi-strix*ctrixi
      ttrix= strrx*ctrixi+strix*ctrrxi
*
*-----propagator in X-channel
*
      prxr= xv*y*sh-sz(1)
      prxi= -sz(2)
      prxm= prxr*prxr+prxi*prxi
      prxri= prxr/prxm
      prxii= -prxi/prxm
      ratgr= grx*agrsz-gix*agisz
      ratgi= grx*agisz+gix*agrsz
      ratcr= ctrrsz*ctrrxi-ctrisz*ctrixi
      ratci= ctrrsz*ctrixi+ctrisz*ctrrxi
      ratgcr= grx*ctrrxi-gix*ctrixi
      ratgci= grx*ctrixi+gix*ctrrxi
      ratr= 1.d0-ratgr*ratcr+ratgi*ratci
      rati= -ratgr*ratci-ratgi*ratcr
      arhrx= xv*y*sh-sz(1)+ratr*sz(1)-rati*sz(2)+(ratgcr*(fzsz(1)-
     #       fzsz(1))-ratgci*(fzx(2)-fzsz(2)))/apis
      arhix= -sz(2)+ratr*sz(2)+rati*sz(1)+(ratgcr*(fzsz(2)-
     #       fzsz(2))+ratgci*(fzx(1)-fzsz(1)))/apis
      brhrx= arhrx*prxri-arhix*prxii
      brhix= arhrx*prxii+arhix*prxri
      brhmx= brhrx*brhrx+brhix*brhix
      rhrx= brhrx/brhmx
      rhix= -brhix/brhmx
*
*-----flux-funxtions are given
*
      psf_gg_11= 0.5d0*omxm*omxm*em4
      psf_gg_12= -0.5d0*omxm*omxm*em4
      psf_gg_21= em2*(2.d0+0.5d0*omxm)
      psf_gg_22= em2*omxm*(2.d0/y-3.d0/2.d0)
      psf_gg_31= -1.d0
      psf_gg_32= 2.d0*(1.d0-1.d0/y)
      psf_gg_41= 0.d0
      psf_gg_42= 2.d0/ys*(1.d0/xpo+1.d0)/xpo-
     #           2.d0/y*(1.d0/xpo+1.d0)/xpo+1.d0/xpo
*
*-----the following ones are given but not used
*
      verx= 4.d0*strrx-1.d0
      psf_gz_41= -y*verx
      psf_gz_42= verx*(-2.d0/y/xpos+2.d0/xpos+y*(1.d0-1.d0/xpo))
      psf_gz_43= (y-2.d0)/xpo
*
      v2a2rx= verx*verx+1.d0
      psf_zz_41= -v2a2rx*xv*ys
      psf_zz_42= v2a2rx*xv*(ys-2.d0/xpos-ys/xpo-y/xpos)
      psf_zz_43= 2.d0*verx*xv*y*(y-2.d0)/xpo
      psf_zz_44= -0.5d0*xv*ys*v2a2rx
*
*-----x2 is initialized
*
      x2l1= (um+dm)*(um+dm)/yvv
      x2l2= s0/yvv
      x2l= dmax1(x2l1,x2l2)
      x2u= 1.d0
      if(x2l.gt.x2u) then
         iz= 0
         ifz(6)= ifz(6)+1
         go to 1
      endif
      zma= yvv*x2l
      zmb= yvv*x2u
      zmas= zma-rsw(1)
      zmbs= zmb-rsw(1)
      atma= zmas/aswi
      atmb= zmbs/aswi
*
      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
          x2jc0= (-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
          x2jc0= (-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
          x2jc0= (-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
          x2jc0= (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
          x2jc0= (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
          x2jc0= (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
          x2jc0= (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
          x2jc0= (-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
          x2jc0= (atmb-atma)/vwmg
       endif
*
       x2jc= x2jc0/yvv
       argx2= vwmg*(x2jc0*x2x+zmat)
       iftn= 1
       s07= s07aaf(argx2,iftn)
       x2= (rsw(1)+aswi*s07)/yvv
       pmjc= 1.d0
*
      else if(opeak.eq.'n') then         
       smjc0= zmb-zma
       x2= (smjc0*x2x+zma)/yvv
       pmjc= 1.d0/((yvv*x2-rsw(1))**2+srwis)
       x2jc= smjc0/yvv
      endif
      omx2= 1.d0-x2
      omx2s= omx2*omx2
      x2s= x2*x2
      x2c= x2s*x2
      x2q= x2c*x2
*
*-----Quantities are initialized in the x2-channel
*
      p2x2r= -x2*y*sh
      p2x2i= 0.d0
      call wtopole(p2x2r,p2x2i,p3qx2,fzx2,fwsx2)
      xgx20= 1.d0/8.d0/gf/wm2
      xgx2r= 1.d0+0.5d0*ccw*((fw0-fw(1))/wm2+p3qw(1)-p3qx2(1))
      xgx2i= -0.5d0*ccw*p3qx2(2)
      agrx2= xgx20*xgx2r
      agix2= xgx20*xgx2i
      agx2m= agrx2*agrx2+agix2*agix2
      grx2= agrx2/agx2m
      gix2= -agix2/agx2m
*
*-----propagator in x2-channel
*
      px2r= x2*y*sh-sw(1)
      px2ir= px2r/(px2r*px2r+swis)
      px2ii= sw(2)/(px2r*px2r+swis)
      xrhrx2= fwsx2(1)-fwsw(1)+sw(1)*(p3qsw(1)-p3qx2(1))-
     #       sw(2)*(p3qsw(2)-p3qx2(2))
      xrhix2= fwsx2(2)-fwsw(2)+sw(1)*(p3qsw(2)-p3qx2(2))+
     #       sw(2)*(p3qsw(1)-p3qx2(1))
      arhrx2= 1.d0+1.d0/apis*(px2ir*(grx2*xrhrx2-gix2*xrhix2)-
     #        px2ii*(grx2*xrhix2+gix2*xrhrx2))
      arhix2= 1.d0/apis*(px2ii*(grx2*xrhrx2-gix2*xrhix2)+
     #        px2ir*(grx2*xrhix2+gix2*xrhrx2))
      arhx2m= arhrx2*arhrx2+arhix2*arhix2
      rhrx2= arhrx2/arhx2m
      rhix2= -arhix2/arhx2m
*
*-----x1 is initialized
*
      xlk= x2*(x2-2.d0*(dmsh+umsh))+(dmsh-umsh)*(dmsh-umsh)
      xlk= sqrt(xlk)
      x1l1= dmsh+2.d0*umh
      x1l2= x2
      x1l3= dmsh+(umsh-dmsh+x2*(1.d0+umsh-dmsh+x2)-omx2*xlk)/x2/2.d0
      x1l= dmax1(x1l1,x1l2,x1l3)
      x1u1= 1.d0
      x1u2= (1.d0-dmh)**2+x2
      x1u3= 1.d0-dmsh*omx2/x2
      x1u4= dmsh+(umsh-dmsh+x2*(1.d0+umsh-dmsh+x2)+omx2*xlk)/x2/2.d0
      x1u= dmin1(x1u1,x1u2,x1u3,x1u4)
      if(x1l.gt.x1u) then
         iz= 0
         ifz(7)= ifz(7)+1
         go to 1
      endif
      x1jc= x1u-x1l
      x1= x1jc*x1x+x1l
      omx1= 1.d0-x1
      omx1s= omx1*omx1
      x1mo= x1-1.d0
      x1mos= x1mo*x1mo
      x1moc= x1mos*x1mo
      x1moq= x1moc*x1mo
      x1mx2= x1-x2
      x2mx1= x2-x1
      x2mx1s= x2mx1*x2mx1
      x2mx1c= x2mx1s*x2mx1
      x2mx1q= x2mx1c*x2mx1
      x2mx1f= x2mx1q*x2mx1
      x2mx1sx= x2mx1f*x2mx1
      x1s= x1*x1
      x1c= x1s*x1
      x1q= x1c*x1
      x1f= x1q*x1
      edm= 1.d0+x2-x1
      edms= edm*edm
      edmc= edms*edm
      edmq= edmc*edm
      edmf= edmq*edm
      edmsx= edmf*edm
      reu= x1-dmsh
      reus= reu*reu
      reuc= reus*reu
      red= 1.d0+x2-x1+dmsh
      reds= red*red
      redc= reds*red
      betus= 1.d0-4.d0*umsh/reus
      betds= 1.d0-4.d0*dmsh/reds
      if(betus.lt.0.d0.or.betds.lt.0.d0) then
         iz= 0
         ifz(8)= ifz(8)+1
         go to 1
      endif
      betu= sqrt(betus)
      betd= sqrt(betds)
*
      jj= jnm
      if(jj.le.3) then
       psjc= 1.d0/reu/betu/opxv
      else if(jj.ge.4) then
       psjc= 1.d0/red/betd/opxv
      endif
*
*-----tau/upsilon is initialized
*
*-----Here is the tau branch
*
      if(jj.le.3) then
       if(ipro.eq.2) then
        taul0= -x1-xv+umsh+dmsh+0.5d0*(1.d0-betu)*opxv*reu
        tauu0= -x1-xv+umsh+dmsh+0.5d0*(1.d0+betu)*opxv*reu
       else 
        taul0= -reu-xv
        tauu0= xv*(reu-1.d0)
       endif
*
       phta= umsh*dmsh*(-1+x2)
       phta= phta+umsh*(1-x1*x2-x1+x2s)
       phta= phta+umshs
       phta= phta+dmsh*(-x1*x2+x1-x2+x2s)
       phta= phta-x2*x1mo*x2mx1
*
       phtb= xv*umsh*dmsh*(-1-2*x1*x2+2*x2+x2s)
       phtb= phtb+xv*umsh*dmshs*(-1+x2)
       phtb= phtb+xv*umsh*(2-2*x1*x2-x1*x2s-3*x1+x1s*x2+
     #    x1s+2*x2s)
       phtb= phtb+xv*umshs*dmsh
       phtb= phtb+xv*umshs*(2-x1)
       phtb= phtb+xv*dmsh*(-2*x1*x2-2*x1*x2s+2*x1+2*x1s*x2
     #    -x1s-2*x2+3*x2s)
       phtb= phtb+xv*dmshs*(-x1*x2+x1-x2+x2s)
       phtb= phtb+xv*(x1*x2*x1mo*x2mx1-2*x2*x1mo*x2mx1)
       phtb= phtb+umsh*dmsh*(-1+4*x1*x2-2*x1+2*x2-3*x2s)
       phtb= phtb+umsh*dmshs*(1-x2)
       phtb= phtb+umsh*(2*x1*x2+3*x1*x2s+x1-3*x1s*x2-x1s-2*x2s)
       phtb= phtb+umshs*dmsh*(1-2*x2)
       phtb= phtb+umshs*(-2+2*x1*x2+3*x1-2*x2s)
       phtb= phtb+umshc*(-2)
       phtb= phtb+dmsh*(2*x1*x2s-2*x1s*x2+x1s-x2s)
       phtb= phtb+dmshs*(x1*x2-x1+x2-x2s)
       phtb= phtb-x1*x2*x1mo*x2mx1
*
       phtc= xv*umsh*dmsh*(-1+2*x1*x2+4*x1*x2s+2*x1-5*x1s*
     #    x2-2*x2s)
       phtc= phtc+xv*umsh*dmshs*(4*x1*x2-2*x1-2*x2s)
       phtc= phtc+xv*umsh*dmshc*(1-x2)
       phtc= phtc+xv*umsh*(2*x1*x2s+x1-2*x1s*x2-2*x1s*x2s
     #    -2*x1s+2*x1c*x2+x1c)
       phtc= phtc+xv*umshs*dmsh*(-2+2*x1*x2+2*x1-x2s)
       phtc= phtc+xv*umshs*dmshs*(-x2)
       phtc= phtc+xv*umshs*(x1*x2s+2*x1-x1s*x2-2*x1s)
       phtc= phtc+xv*umshc*dmsh*(-1)
       phtc= phtc+xv*umshc*(x1)
       phtc= phtc+xv*dmsh*(4*x1*x2s-3*x1s*x2-3*x1s*x2s+
     #    x1s+3*x1c*x2-x1c-x2s)
       phtc= phtc+xv*dmshs*(3*x1*x2s-x1-3*x1s*x2+2*x1s+
     #    x2-2*x2s)
       phtc= phtc+xv*dmshc*(x1*x2-x1+x2-x2s)
       phtc= phtc+xv*(-x1*x2*x1mo*x2mx1+x1s*x2*x1mo*
     #    x2mx1)
       phtc= phtc+xvs*umsh*dmsh*(-3*x1*x2+x1+2*x2s)
       phtc= phtc+xvs*umsh*dmshs*(-1+x2)
       phtc= phtc+xvs*umsh*(1-2*x1*x2-2*x1*x2s-2*x1+2*x1s
     #    *x2+x1s+2*x2s)
       phtc= phtc+xvs*umshs*dmsh*(x2)
       phtc= phtc+xvs*umshs*(2-x1*x2-2*x1+x2s)
       phtc= phtc+xvs*umshc
       phtc= phtc+xvs*dmsh*(-x1*x2-2*x1*x2s+x1+2*x1s*x2
     #    -x1s-x2+2*x2s)
       phtc= phtc+xvs*dmshs*(-x1*x2+x1-x2+x2s)
       phtc= phtc+xvs*(x1*x2*x1mo*x2mx1-x2*x1mo*x2mx1)
       phtc= phtc+umsh*dmsh*(-x1*x2-2*x1*x2s+x1+2*x1s*x2-
     #    x1s-x2+2*x2s)
       phtc= phtc+umsh*dmshs*(-x1*x2+x1-x2+x2s)
       phtc= phtc+umsh*(-x1*x2-2*x1*x2s+2*x1s*x2+x1s*
     #    x2s-x1c*x2+x2s)
       phtc= phtc+umshs*dmsh*(-3*x1*x2+x1+2*x2s)
       phtc= phtc+umshs*dmshs*(-1+x2)
       phtc= phtc+umshs*(1-2*x1*x2-2*x1*x2s-2*x1+2*x1s*x2
     #    +x1s+2*x2s)
       phtc= phtc+umshc*dmsh*(x2)
       phtc= phtc+umshc*(2-x1*x2-2*x1+x2s)
       phtc= phtc+umshq
*
       phtdisc= umsh*dmsh*opxvs*(-20*x1*x2s+8*x1*x2c+4*x1*x2q
     #    +10*x1s*x2+16*x1s*x2s-8*x1s*x2c+6*x1s*
     #    x2q-14*x1c*x2+2*x1c*x2s-14*x1c*x2c+2*
     #    x1c+2*x1q*x2+8*x1q*x2s-2*x1q+8*x2c-8*x2q)
       phtdisc= phtdisc+umsh*dmshs*opxvs*(10*x1*x2-20*x1*x2s+4*x1*
     #    x2c-6*x1*x2q+10*x1s*x2-8*x1s*x2s+18*x1s*
     #    x2c-8*x1s+6*x1c*x2-12*x1c*x2s+2*x1c-2*
     #    x2s+8*x2c-2*x2q)
       phtdisc= phtdisc+umsh*dmshc*opxvs*(-2*x1*x2+10*x1*x2s-10*x1
     #    *x2c+2*x1-10*x1s*x2+8*x1s*x2s+2*x1s-2*x2+2*x2q)
       phtdisc= phtdisc+umsh*dmshq*opxvs*(4*x1*x2-2*x1*x2s-2*x1+2*
     #    x2-4*x2s+2*x2c)
       phtdisc= phtdisc+umsh*opxvs*(8*x1*x2c+8*x1*x2q-2*x1s*x2s
     #    -16*x1s*x2c-2*x1s*x2q-2*x1c*x2+4*x1c*
     #    x2s+4*x1c*x2c-2*x1c*x2q+4*x1q*x2+4*x1q*
     #    x2c-2*x1f*x2-2*x1f*x2s-4*x2q)
       phtdisc= phtdisc+umshs*dmsh*opxvs*(-8*x1*x2-16*x1*x2s+24*x1
     #    *x2c-2*x1*x2q-10*x1+16*x1s*x2-8*x1s*x2s+6
     #    *x1s*x2c+10*x1s-10*x1c*x2-4*x1c*x2s+2*x1c+8*x2-8*x2q)
       phtdisc= phtdisc+umshs*dmshs*opxvs*(1-20*x1*x2+4*x1*x2s-6
     #    *x1*x2c+10*x1+6*x1s*x2+6*x1s*x2s-6*x1s-8*
     #    x2+20*x2s-8*x2c+x2q)
       phtdisc= phtdisc+umshs*dmshc*opxvs*(-2+2*x1*x2-4*x1*x2s
     #    +2*x1+2*x2c)
       phtdisc= phtdisc+umshs*dmshq*opxvs*(1-2*x2+x2s)
       phtdisc= phtdisc+umshs*opxvs*(8*x1*x2+16*x1*x2s+16*x1*x2c
     #    +8*x1*x2q-16*x1s*x2-12*x1s*x2s-16*x1s*x2c
     #    +x1s*x2q+x1s+4*x1c*x2+4*x1c*x2s-2*x1c
     #    *x2c-2*x1c+4*x1q*x2+x1q*x2s+x1q-8*x2s-8*x2q)
       phtdisc= phtdisc+umshc*dmsh*opxvs*(8+8*x1*x2+4*x1*x2s-20*x1
     #    +6*x1s*x2+2*x1s-8*x2c)
       phtdisc= phtdisc+umshc*dmshs*opxvs*(-2-6*x1*x2+2*x1+8*x2
     #    -2*x2s)
       phtdisc= phtdisc+umshc*dmshc*opxvs*(-2+2*x2)
       phtdisc= phtdisc+umshc*opxvs*(-4+16*x1*x2+16*x1*x2s+8*x1
     #    *x2c+8*x1-16*x1s*x2-2*x1s*x2s-2*x1s-2*
     #    x1c*x2-2*x1c-16*x2s-4*x2q)
       phtdisc= phtdisc+umshq*dmsh*opxvs*(8-2*x1-8*x2)
       phtdisc= phtdisc+umshq*dmshs*opxvs
       phtdisc= phtdisc+umshq*opxvs*(-8+8*x1*x2+8*x1+x1s-8*x2s)
       phtdisc= phtdisc+umshf*opxvs*(-4)
       phtdisc= phtdisc+dmsh*opxvs*(-2*x1*x2q+2*x1s*x2c+6*x1s
     #    *x2q+2*x1c*x2s-10*x1c*x2c-4*x1c*x2q-2*
     #    x1q*x2+2*x1q*x2s+8*x1q*x2c+2*x1f*x2-4*x1f*x2s)
       phtdisc= phtdisc+dmshs*opxvs*(2*x1*x2c-6*x1*x2q-6*x1s*
     #    x2s+6*x1s*x2c+6*x1s*x2q+2*x1c*x2+6*x1c*
     #    x2s-12*x1c*x2c-6*x1q*x2+6*x1q*x2s+x1q+x2q)
       phtdisc= phtdisc+dmshc*opxvs*(2*x1*x2s+2*x1*x2c-4*x1*x2q
     #    +2*x1s*x2-10*x1s*x2s+8*x1s*x2c+6*x1c*x2
     #    -4*x1c*x2s-2*x1c-2*x2c+2*x2q)
       phtdisc= phtdisc+dmshq*opxvs*(-2*x1*x2+4*x1*x2s-2*x1*x2c
     #    -2*x1s*x2+x1s*x2s+x1s+x2s-2*x2c+x2q)
       phtdisc= phtdisc+opxvs*(x2s*x1mos*x2mx1s*x1s)
       if(phtdisc.lt.1.d-15) then
        taul= taul0
        tauu= tauu0
       else
        if(phtb.gt.0.d0) then
         phtq= -0.5d0*(phtb+sqrt(phtdisc))
        else if(phtb.lt.0.d0) then
         phtq= -0.5d0*(phtb-sqrt(phtdisc))
        endif
        ataul= phtq/phta
        atauu= phtc/phtq
        taupl= dmin1(ataul,atauu)
        taupu= dmax1(ataul,atauu)
        taul= dmax1(taul0,taupl)
        tauu= dmin1(tauu0,taupu)
       endif
       if(taul.gt.tauu) then
         iz= 0
         ifz(9)= ifz(9)+1
         go to 1
       endif
*
*-----Here is the ups branch
*
      else if(jj.ge.4) then
       upsl0= -1.d0-x2+x1-xv+0.5d0*(1.d0-betd)*opxv*red
       upsu0= -1.d0-x2+x1-xv+0.5d0*(1.d0+betd)*opxv*red
*
       phta= umsh*dmsh*(-1+x2)
       phta= phta+umsh*(1-x1*x2-x1+x2s)
       phta= phta+umshs
       phta= phta+dmsh*(-x1*x2+x1-x2+x2s)
       phta= phta-x2*x1mo*x2mx1
*
       phtb= xv*umsh*dmsh*(-2+2*x1*x2+2*x2-2*x2s)
       phtb= phtb+xv*umsh*dmshs*(1-x2)
       phtb= phtb+xv*umsh*(1+2*x1*x2s-x1s*x2-x1s-x2+x2s-x2c)
       phtb= phtb+xv*umshs*dmsh*(-1)
       phtb= phtb+xv*umshs*(1+x1-x2)
       phtb= phtb+xv*dmsh*(-2*x1*x2+3*x1*x2s+x1-2*x1s*x2+
     #    x1s-x2+x2s-x2c)
       phtb= phtb+xv*dmshs*(x1*x2-x1+x2-x2s)
       phtb= phtb+xv*(-x1*x2*x1mo*x2mx1-x2*x1mo*x2mx1+
     #    x2s*x1mo*x2mx1)
       phtb= phtb+umsh*dmsh*(-2+2*x1)
       phtb= phtb+umsh*dmshs*(1-x2)
       phtb= phtb+umsh*(1-2*x1*x2-2*x1*x2s-2*x1+x1s*x2+
     #    x1s+x2+x2s+x2c)
       phtb= phtb+umshs*dmsh*(-1)
       phtb= phtb+umshs*(1-x1+x2)
       phtb= phtb+dmsh*(2*x1*x2-x1*x2s+x1-x1s-x2-x2s+x2c)
       phtb= phtb+dmshs*(x1*x2-x1+x2-x2s)
       phtb= phtb+x1*x2*x1mo*x2mx1-x2*x1mo*x2mx1-x2s*x1mo*x2mx1
*
       phtc= xv*umsh*dmsh*(4*x1*x2s-2*x1s*x2-2*x2c)
       phtc= phtc+xv*umsh*dmshs*(x1*x2-x1+x2-x2s)
       phtc= phtc+xv*umsh*(3*x1*x2+3*x1*x2s+3*x1*x2c+x1-3*
     #    x1s*x2-3*x1s*x2s-2*x1s+x1c*x2+x1c-x2-
     #    x2s-x2c-x2q)
       phtc= phtc+xv*umshs*dmsh*x1mx2
       phtc= phtc+xv*umshs*(2*x1*x2+x1-x1s-x2-x2s)
       phtc= phtc+xv*dmsh*(-2*x1*x2+x1*x2s+4*x1*x2c+x1s*
     #    x2-5*x1s*x2s+x1s+2*x1c*x2-x1c+x2s-x2c-x2q)
       phtc= phtc+xv*dmshs*(-2*x1*x2+2*x1*x2s-x1s*x2+
     #    x1s+x2s-x2c)
       phtc= phtc+xv*(-x1*x2*x1mo*x2mx1-2*x1*x2s*x1mo*
     #    x2mx1+x1s*x2*x1mo*x2mx1+x2s*x1mo*x2mx1+
     #    x2c*x1mo*x2mx1)
       phtc= phtc+xvs*umsh*dmsh*(x1*x2-x1+x2-x2s)
       phtc= phtc+xvs*umsh*(x1*x2+2*x1*x2s+x1-x1s*x2-x1s-x2-x2c)
       phtc= phtc+xvs*umshs*x1mx2
       phtc= phtc+xvs*dmsh*(-2*x1*x2+2*x1*x2s-x1s*x2+
     #    x1s+x2s-x2c)
       phtc= phtc+xvs*(-x1*x2*x1mo*x2mx1+x2s*x1mo*x2mx1)
       phtc= phtc+umsh*dmsh*(x1*x2+2*x1*x2s+x1-x1s*x2-x1s
     #    -x2-x2c)
       phtc= phtc+umsh*dmshs*(x1*x2-x1+x2-x2s)
       phtc= phtc+umshs*dmsh*x1mx2
       phtc= phtc+dmsh*(2*x1*x2s+x1*x2c-x1s*x2-2*x1s*
     #    x2s+x1c*x2-x2c)
       phtc= phtc+dmshs*(-2*x1*x2+2*x1*x2s-x1s*x2+x1s+x2s-x2c)
*
       phtdisc= umsh*dmsh*opxvs*(-8*x2mx1*x1mos-6*x2mx1*x1moc+2*
     #    x2mx1s*x1mo-2*x2mx1s*x1mos+4*x2mx1s*x1moc+2*
     #    x2mx1c*x1mo-4*x2mx1c*x1mos-6*x2mx1c*x1moc+2*
     #    x2mx1c-2*x2mx1q*x1mo-10*x2mx1q*x1mos+6*x2mx1q
     #    -2*x2mx1f*x1mo+6*x2mx1f+2*x2mx1sx)
       phtdisc= phtdisc+umsh*dmshs*opxvs*(12*x2mx1*x1mos+6*x2mx1*
     #    x1moc+6*x2mx1s*x1mo-2*x2mx1s*x1mos-6*x2mx1s*
     #    x1moc-4*x2mx1c*x1mo-6*x2mx1c*x1mos-2*x2mx1c
     #    +6*x2mx1q*x1mo+4*x2mx1q+6*x2mx1f)
       phtdisc= phtdisc+umsh*dmshc*opxvs*(-8*x2mx1*x1mos-2*x2mx1*
     #    x1moc-10*x2mx1s*x1mo+2*x2mx1s*x1mos+10*x2mx1c
     #    *x1mo-2*x2mx1c+6*x2mx1q)
       phtdisc= phtdisc+umsh*dmshq*opxvs*(2*x2mx1*x1mos+4*x2mx1s*
     #    x1mo+2*x2mx1c)
       phtdisc= phtdisc+umsh*opxvs*(2*x2mx1*x1mos+2*x2mx1*x1moc-2*
     #    x2mx1s*x1mo+2*x2mx1s*x1mos+2*x2mx1s*x1moc-8*
     #    x2mx1c*x1mo-6*x2mx1c*x1mos-2*x2mx1c*x1moc-12*
     #    x2mx1q*x1mo-10*x2mx1q*x1mos-2*x2mx1q*x1moc-8*
     #    x2mx1f*x1mo-4*x2mx1f*x1mos-2*x2mx1sx*x1mo)
       phtdisc= phtdisc+umshs*dmsh*opxvs*(10*x2mx1*x1mo+8*x2mx1*
     #    x1mos+4*x2mx1s*x1mo-8*x2mx1s*x1mos+2*x2mx1s+2*
     #    x2mx1c*x1mo+4*x2mx1c*x1mos+8*x2mx1c+8*x2mx1q*
     #    x1mo+10*x2mx1q+4*x2mx1f-4*x1mos)
       phtdisc= phtdisc+umshs*dmshs*opxvs*(-6*x2mx1*x1mo-10*x2mx1*
     #    x1mos-2*x2mx1s*x1mo+6*x2mx1s*x1mos-6*x2mx1s
     #    +12*x2mx1c*x1mo+8*x2mx1c+6*x2mx1q+6*x1mos)
       phtdisc= phtdisc+umshs*dmshc*opxvs*(-2*x2mx1*x1mo+4*x2mx1*
     #    x1mos+8*x2mx1s*x1mo+2*x2mx1s+4*x2mx1c-4*x1mos)
       phtdisc= phtdisc+umshs*dmshq*opxvs*(2*x2mx1*x1mo+x2mx1s+
     #    x1mos)
       phtdisc= phtdisc+umshs*opxvs*(-4*x2mx1*x1mo-2*x2mx1*x1mos-
     #    10*x2mx1s*x1mo-6*x2mx1s*x1mos+x2mx1s-6*x2mx1c
     #    *x1mo-2*x2mx1c*x1mos+4*x2mx1c+2*x2mx1q*x1mo+
     #    x2mx1q*x1mos+6*x2mx1q+2*x2mx1f*x1mo+4*x2mx1f
     #    +x2mx1sx+x1mos)
       phtdisc= phtdisc+umshc*dmsh*opxvs*(-4*x2mx1*x1mo-2*x2mx1+6*
     #    x2mx1s*x1mo+4*x2mx1s+6*x2mx1c+6*x1mo)
       phtdisc= phtdisc+umshc*dmshs*opxvs*(6*x2mx1*x1mo-2*x2mx1+6*
     #    x2mx1s-6*x1mo)
       phtdisc= phtdisc+umshc*dmshc*opxvs*(2*x2mx1+2*x1mo)
       phtdisc= phtdisc+umshc*opxvs*(-2*x2mx1*x1mo+2*x2mx1+2*
     #    x2mx1s*x1mo+6*x2mx1s+2*x2mx1c*x1mo+6*x2mx1c+2
     #    *x2mx1q-2*x1mo)
       phtdisc= phtdisc+umshq*dmsh*opxvs*(-2+2*x2mx1)
       phtdisc= phtdisc+umshq*dmshs*opxvs
       phtdisc= phtdisc+umshq*opxvs*(1+2*x2mx1+x2mx1s)
       phtdisc= phtdisc+dmsh*opxvs*(-4*x2mx1s*x1mos-6*x2mx1s*
     #    x1moc-2*x2mx1s*x1moq-2*x2mx1c*x1mo-10*x2mx1c*
     #    x1mos-4*x2mx1c*x1moc+2*x2mx1c*x1moq-6*
     #    x2mx1q*x1mo-8*x2mx1q*x1mos+2*x2mx1q*x1moc-6*
     #    x2mx1f*x1mo-2*x2mx1f*x1mos-2*x2mx1sx*x1mo)
       phtdisc= phtdisc+dmshs*opxvs*(6*x2mx1s*x1mos+6*x2mx1s*
     #    x1moc+x2mx1s*x1moq+6*x2mx1c*x1mo+8*x2mx1c*
     #    x1mos-2*x2mx1c*x1moc+4*x2mx1q*x1mo-6*x2mx1q*
     #    x1mos+x2mx1q-2*x2mx1f*x1mo+2*x2mx1f+x2mx1sx)
       phtdisc= phtdisc+dmshc*opxvs*(-4*x2mx1s*x1mos-2*x2mx1s*
     #    x1moc-6*x2mx1c*x1mo-2*x2mx1c*x1mos+2*x2mx1q*
     #    x1mo-2*x2mx1q+2*x2mx1f)
       phtdisc= phtdisc+dmshq*opxvs*(x2mx1s*x1mos+2*x2mx1c*x1mo
     #    +x2mx1q)
       phtdisc= phtdisc+opxvs*(x2mx1s*x1mos*x2s*edms)
*
       if(phtdisc.lt.1.d-15) then
        upsl= upsl0
        upsu= upsu0
       else
        if(phtb.gt.0.d0) then
         phtq= -0.5d0*(phtb+sqrt(phtdisc))
        else if(phtb.lt.0.d0) then
         phtq= -0.5d0*(phtb-sqrt(phtdisc))
        endif
        aupsl= phtq/phta
        aupsu= phtc/phtq
        upspl= dmin1(aupsl,aupsu)
        upspu= dmax1(aupsl,aupsu)
        upsl= dmax1(upsl0,upspl)
        upsu= dmin1(upsu0,upspu)
       endif
       if(upsl.gt.upsu) then
         iz= 0
         ifz(9)= ifz(9)+1
         go to 1
       endif
      endif
*
*-----Quark propagators are mapped away (for no angular cuts)
*
      if(asa(3).eq.-1.d0.and.asa(4).eq.-1.d0.and.
     #   bsa(3).eq.1.d0.and.bsa(4).eq.1.d0) then
       if(jj.eq.1) then
        taujc= betu*opxv*reu
        tau= taujc*angx+taul0
        taus= tau*tau
       else if(jj.eq.2) then
        taujc= log((umsh-tauu0)/(umsh-taul0))
        tau= umsh-exp(taujc*angx+log(umsh-taul0))
        taus= tau*tau
       else if(jj.eq.3) then
        taujc= 1.d0/(umsh-tauu0)-1.d0/(umsh-taul0)
        tau= umsh-1.d0/(taujc*angx+1.d0/(umsh-taul0))
        taus= tau*tau
       else if(jj.eq.4) then
        upsjc= betd*opxv*red
        ups= upsjc*angx+upsl0
        upss= ups*ups
       else if(jj.eq.5) then
        upsjc= log((dmsh-upsu0)/(dmsh-upsl0))
        ups= dmsh-exp(upsjc*angx+log(dmsh-upsl0))
        upss= ups*ups
       else if(jj.eq.6) then
        upsjc= 1.d0/(dmsh-upsu0)-1.d0/(dmsh-upsl0)
        ups= dmsh-1.d0/(upsjc*angx+1.d0/(dmsh-upsl0))
        upss= ups*ups
       endif
      else
       if(jj.le.3) then
        taujc0= tauu-taul
        tau= taujc0*angx+taul
        taus= tau*tau
       else if(jj.ge.4) then
        upsjc0= upsu-upsl
        ups= upsjc0*angx+upsl
        upss= ups*ups
       endif
       if(jj.le.3.and.(tau-umsh).gt.0.d0) then
         iz= 0
         ifz(10)= ifz(10)+1
         go to 1
       endif
       if(jj.ge.4.and.(ups-dmsh).gt.0.d0) then
         iz= 0
         ifz(10)= ifz(10)+1
         go to 1
       endif
       if(jj.eq.1) then
        taujc= taujc0
       else if(jj.eq.2) then
        taujc= taujc0/(tau-umsh)
       else if(jj.eq.3) then
        taujc= taujc0/(tau-umsh)/(tau-umsh)
       else if(jj.eq.4) then
        upsjc= upsjc0
       else if(jj.eq.5) then
        upsjc= upsjc0/(ups-dmsh)
       else if(jj.eq.6) then
        upsjc= upsjc0/(ups-dmsh)/(ups-dmsh)
       endif
      endif
*
*-----z is initialized
*
      if(omx2.eq.0.d0.or.x1.eq.0.d0) then
         iz= 0
         ifz(11)= ifz(11)+1
         go to 1
      else
         if(jj.le.3) then
          ct= (-omxv-2.d0*(tau+xv-umsh)/reu)/opxv/betu
          cp= (1.d0-2.d0*(x1-x2-umsh)/omx2/reu)/betu
         else if(jj.ge.4) then
          ct= (-omxv-2.d0*(ups+xv-dmsh)/red)/opxv/betd
          cp= (1.d0-2.d0*(omx1+umsh)/omx2/red)/betd
         endif
         if(ct.lt.-1.d0.or.ct.gt.1.d0) then
            iz= 0
            ifz(12)= ifz(12)+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(13)= ifz(13)+1
            go to 1
         else
            sp= sqrt(1.d0-cp*cp)
         endif
      endif
*
      zca= -1.d0
      zcb= -xv-0.5d0*omx2*(omxv+opxv*ct*cp)
      zcc= -xvs+xv*omx2s-xv*omxv*omx2-0.25d0*omx2s*opxvs*
     #     (ct*ct+cp*cp)-xv*opxv*omx2*ct*cp-0.5d0*omxv*opxv*
     #     omx2s*ct*cp
      zdisc= 0.5d0*opxv*omx2*st*sp
      if(zcb.gt.0.d0) then
       zcq= -(zcb+zdisc)
      else if(zcb.lt.0.d0) then
       zcq= -(zcb-zdisc)
      endif
      azl= -zcq
      azu= zcc/zcq
      zl= dmin1(azl,azu)
      zu= dmax1(azl,azu)
      if(jj.le.3) then
       zl1= -tau-x1+x2-xv+umsh+dmsh-0.5d0*(1.d0+betd)*opxv*red
       zu1= -tau-x1+x2-xv+umsh+dmsh-0.5d0*(1.d0-betd)*opxv*red
       zl2= -1.d0+x2-xv
       zu2= -xv*x2
*
       phca= 16.d0*umsh
       phca= phca+8.d0*dmsh*x1
       phca= phca-4.d0*dmshs
       phca= phca-4.d0*x1s
*
       phcb= xv*umsh*dmsh*(8)
       phcb= phcb+xv*umsh*(32-8*x1+16*x2)
       phcb= phcb+xv*dmsh*(-8+8*x1+16*x2)
       phcb= phcb+xv*dmshs*(-8)
       phcb= phcb+xv*(-16*x1*x2-8*x1+16*x2)
       phcb= phcb+umsh*dmsh*(-8*x2)
       phcb= phcb+umsh*(16+8*x1*x2+16*x1-32*x2+16*tau)
       phcb= phcb+umshs*(-16)
       phcb= phcb+dmsh*(8*x1+8*x2*tau-8*x2-8*tau)
       phcb= phcb-8*x1*x2*tau+8*x1*x2-8*x1*tau-8*x1s+16*x2*tau
*
       phcc= xv*umsh*dmsh*(16-24*x2)
       phcc= phcc+xv*umsh*(16+32*x1*x2+8*x1+8*x2*tau-32*x2-16
     #    *x2s+8*tau)
       phcc= phcc+xv*umshs*(-16-8*x2)
       phcc= phcc+xv*dmsh*(16*x1*x2-8*x1+8*x2*tau+8*x2-16*x2s-8*tau)
       phcc= phcc+xv*(-16*x1*x2*tau+16*x1*x2+16*x1*x2s-8*x1
     #    -16*x1s*x2+24*x2*tau+8*x2-16*x2s-8*tau)
       phcc= phcc+xvs*umsh*dmsh*(8)
       phcc= phcc+xvs*umsh*(8+16*x2)
       phcc= phcc+xvs*umshs*(-4)
       phcc= phcc+xvs*dmsh*(-8+16*x2)
       phcc= phcc+xvs*dmshs*(-4)
       phcc= phcc+xvs*(-4-16*x1*x2+16*x2)
       phcc= phcc+umsh*(-8*x1*x2+16*x1-24*x2*tau-16*x2+8*
     #    x2s*tau+8*x2s+16*tau)
       phcc= phcc+umshs*(-16+16*x2-4*x2s)
       phcc= phcc+8*x1*x2*tau+8*x1*x2-8*x1*tau-4*x1s+8*x2*tau
     #    +8*x2*taus-8*x2s*tau-4*x2s*taus-4*x2s-4*taus
*
       phdisc= 256.d0*phta*(tau-taupl)*(tau-taupu)
       if(phdisc.lt.1.d-15) then
        zl3= zl
        zu3= zu
       else
        if(phcb.gt.0.d0) then
         phcq= -0.5d0*(phcb+sqrt(phdisc))
        else if(phcb.lt.0.d0) then
         phcq= -0.5d0*(phcb-sqrt(phdisc))
        endif
        azl3= phcq/phca
        azu3= phcc/phcq
        zl3= dmin1(azl3,azu3)
        zu3= dmax1(azl3,azu3)
       endif
       zln= dmax1(zl,zl1,zl2,zl3)
       zun= dmin1(zu,zu1,zu2,zu3)
      else if(jj.ge.4) then
       zl1= -1.d0+x1-ups-xv-0.5d0*(1.d0+betu)*opxv*reu
       zu1= -1.d0+x1-ups-xv-0.5d0*(1.d0-betu)*opxv*reu
       zl2= -1.d0+x2-xv
       zu2= -xv*x2
*
       phca= dmsh*(8+8*x1-8*x2)
       phca= phca+dmshs*(-4)
       phca= phca-4+8*x1*x2+8*x1-4*x1s-8*x2-4*x2s
*
       phcb= xv*umsh*dmsh*(8)
       phcb= phcb+xv*umsh*(-8-8*x1+8*x2)
       phcb= phcb+xv*dmsh*(16+8*x1-8*x2)
       phcb= phcb+xv*dmshs*(-8)
       phcb= phcb+xv*(-8+16*x1*x2+8*x1-8*x2-16*x2s)
       phcb= phcb+umsh*dmsh*(8)
       phcb= phcb+umsh*(-8+8*x1-8*x2-16*ups)
       phcb= phcb+dmsh*(16-8*x1*x2-8*x2*ups-16*x2+8*x2s+8*ups)
       phcb= phcb+dmshs*(-8+8*x2)
       phcb= phcb-8+8*x1*x2*ups+8*x1*x2+8*x1*ups+16*x1-8*x1s
     #    -8*x2-8*x2s*ups-8*ups
*
       phcc= xv*umsh*dmsh*(16-8*x2)
       phcc= phcc+xv*umsh*(-16+16*x1*x2+8*x1-8*x2*ups-16*
     #    x2s-8*ups)
       phcc= phcc+xv*umshs*(-8)
       phcc= phcc+xv*dmsh*(16-8*x1-8*x2*ups-8*x2+8*ups)
       phcc= phcc+xv*dmshs*(-8+8*x2)
       phcc= phcc+xv*(-8+16*x1*x2*ups+16*x1*x2+16*x1*x2s+8
     #    *x1-16*x1s*x2+8*x2*ups-16*x2s*ups-16*x2s-8*ups)
       phcc= phcc+xvs*umsh*dmsh*(8)
       phcc= phcc+xvs*umsh*(-8)
       phcc= phcc+xvs*umshs*(-4)
       phcc= phcc+xvs*dmsh*(8)
       phcc= phcc+xvs*dmshs*(-4)
       phcc= phcc+xvs*(-4+16*x1*x2-16*x2s)
       phcc= phcc+umsh*dmsh*(8-8*x2)
       phcc= phcc+umsh*(-8+8*x1+8*x2*ups-8*ups)
       phcc= phcc+umshs*(-4)
       phcc= phcc+dmsh*(8+8*x1*x2-8*x1-16*x2*ups-8*x2+8*x2s
     #    *ups+8*ups)
       phcc= phcc+dmshs*(-4+8*x2-4*x2s)
       phcc= phcc-4-8*x1*x2*ups+8*x1*ups+8*x1-4*x1s+8*x2*ups
     #    +8*x2*upss-4*x2s*upss-8*ups-4*upss
*
       phdisc= 256.d0*phta*(ups-upspl)*(ups-upspu)
       if(phdisc.lt.1.d-15) then
        zl3= zl
        zu3= zu
       else
        if(phcb.gt.0.d0) then
         phcq= -0.5d0*(phcb+sqrt(phdisc))
        else if(phcb.lt.0.d0) then
         phcq= -0.5d0*(phcb-sqrt(phdisc))
        endif
        azl3= phcq/phca
        azu3= phcc/phcq
        zl3= dmin1(azl3,azu3)
        zu3= dmax1(azl3,azu3)
       endif
       zln= dmax1(zl,zl1,zl2,zl3)
       zun= dmin1(zu,zu1,zu2,zu3)
      endif
*
      if(zln.gt.zun) then
       iz= 0
       ifz(14)= ifz(14)+1
       go to 1
      endif
*
*-----the overall sqrt is mapped away
*
      sqz= 0.5d0*opxv*omx2*st*sp
      if(sqz.eq.0.d0) then
       iz= 0
       ifz(15)= ifz(15)+1
       go to 1
      endif
      dzl= (zln-zl)/sqz
      dzu= (zu-zun)/sqz
      if(zl.eq.zln) then
       tl= -pi/2.d0
      else
       if(dzl.gt.0.d0.and.dzl.lt.1.d-10) then
        tl= -pi/2.d0+sqrt(2.d0*dzl)*(1.d0-dzl/12.d0)
       else
        ifas= 1
        aa= (zcb-zln)/sqz
        tl= -s09aaf(aa,ifas)
        if(ifas.ne.0) then
         iz= 0
         ifz(16)= ifz(16)+1
          go to 1
        endif
       endif
      endif
      if(zu.eq.zun) then
       tu= pi/2.d0
      else
       if(dzu.gt.0.d0.and.dzu.lt.1.d-10) then
        tu= pi/2.d0-sqrt(2.d0*dzu)*(1.d0-dzu/12.d0)
       else
        ifas= 1
        bb= (zcb-zun)/sqz
        tu= -s09aaf(bb,ifas)
        if(ifas.ne.0) then
         iz= 0
         ifz(17)= ifz(17)+1
         go to 1
        endif
       endif
      endif
*
      zjc= tu-tl
      arg= zjc*zx+tl
      z= zcb+sqz*sin(arg)
      zs= z*z
      zc= zs*z
      zxi= 1.d0/(z+1.d0+2.d0*xv)
*
      bcphi= z+xv+0.5d0*omx2*(opxv*ct*cp+omxv)
      acphi= 0.5d0*omx2*opxv*st*sp
      if((acphi+bcphi).lt.0.d0.or.(acphi-bcphi).lt.0.d0) then
       iz= 0
       ifz(18)= ifz(18)+1
       go to 1
      endif
*
      by= -1.d0-xv+x2-z
*
*-----Quantities are initialized in the Y-channel
*
      p2yr= -by*y*sh
      p2yi= 0.d0
      call wtopole(p2yr,p2yi,p3qy,fzy,fwsy)
      xgy0= 1.d0/8.d0/gf/wm2
      xgyr= 1.d0+0.5d0*ccw*((fw0-fw(1))/wm2+p3qw(1)-p3qy(1))
      xgyi= -0.5d0*ccw*p3qy(2)
      agry= xgy0*xgyr
      agiy= xgy0*xgyi
      agym= agry*agry+agiy*agiy
      gry= agry/agym
      giy= -agiy/agym
*
*-----propagator in Y-channel
*
      pyr= by*y*sh-sw(1)
      pyir= pyr/(pyr*pyr+swis)
      pyii= sw(2)/(pyr*pyr+swis)
      xrhry= fwsy(1)-fwsw(1)+sw(1)*(p3qsw(1)-p3qy(1))-
     #       sw(2)*(p3qsw(2)-p3qy(2))
      xrhiy= fwsy(2)-fwsw(2)+sw(1)*(p3qsw(2)-p3qy(2))+
     #       sw(2)*(p3qsw(1)-p3qy(1))
      arhry= 1.d0+1.d0/apis*(pyir*(gry*xrhry-giy*xrhiy)-
     #        pyii*(gry*xrhiy+giy*xrhry))
      arhiy= 1.d0/apis*(pyii*(gry*xrhry-giy*xrhiy)+
     #        pyir*(gry*xrhiy+giy*xrhry))
      arhym= arhry*arhry+arhiy*arhiy
      rhry= arhry/arhym
      rhiy= -arhiy/arhym
*
*-----vertices are initialized
*
      jfl= 1
      p1s= -by*y*sh
      p2s= -x2*y*sh
      ps= xv*y*sh
      rm12= rnm2
      rm22= rnm2
      rm32= rnm2
      call wtocff(jfl,p1s,p2s,ps,rm12,rm22,rm32,c00,c01,c02,c03)
      if(jfl.eq.-1) then
         iz= 0
         ifz(19)= ifz(19)+1
         go to 1
      endif
*
      jfl= 1
      p1s= -by*y*sh
      p2s= -x2*y*sh
      ps= xv*y*sh
      rm12= rnm2
      rm22= tqm2
      rm32= rnm2
      call wtocff(jfl,p1s,p2s,ps,rm12,rm22,rm32,ct0,ct1,ct2,ct3)
      if(jfl.eq.-1) then
         iz= 0
         ifz(20)= ifz(20)+1
         go to 1
      endif
*
      jfl= 1
      p1s= -by*y*sh
      p2s= -x2*y*sh
      ps= xv*y*sh
      rm12= tqm2
      rm22= rnm2
      rm32= tqm2
      call wtocff(jfl,p1s,p2s,ps,rm12,rm22,rm32,ctt0,ctt1,ctt2,ctt3)
      if(jfl.eq.-1) then
         iz= 0
         ifz(21)= ifz(21)+1
         go to 1
      endif
*
      do k=1,2
       cph_11(k)= ctt1(k,1)+0.5d0*ct1(k,1)-3.d0/2.d0*c01(k,1) 
       cph_12(k)= ctt1(k,2)+0.5d0*ct1(k,2)-3.d0/2.d0*c01(k,2)  
       cph_21(k)= ctt2(k,1)+0.5d0*ct2(k,1)-3.d0/2.d0*c02(k,1)  
       cph_22(k)= ctt2(k,2)+0.5d0*ct2(k,2)-3.d0/2.d0*c02(k,2)   
       cph_23(k)= ctt2(k,3)+0.5d0*ct2(k,3)-3.d0/2.d0*c02(k,3)   
       cph_24(k)= ctt2(k,4)+0.5d0*ct2(k,4)-3.d0/2.d0*c02(k,4)   
       cph_31(k)= ctt3(k,1)+0.5d0*ct3(k,1)-3.d0/2.d0*c03(k,1)      
       cph_32(k)= ctt3(k,2)+0.5d0*ct3(k,2)-3.d0/2.d0*c03(k,2)    
       cph_33(k)= ctt3(k,3)+0.5d0*ct3(k,3)-3.d0/2.d0*c03(k,3)    
       cph_34(k)= ctt3(k,4)+0.5d0*ct3(k,4)-3.d0/2.d0*c03(k,4)    
       cph_35(k)= ctt3(k,5)+0.5d0*ct3(k,5)-3.d0/2.d0*c03(k,5)    
       cph_36(k)= ctt3(k,6)+0.5d0*ct3(k,6)-3.d0/2.d0*c03(k,6)    
*
       cp_11(k)= ctt1(k,1)+ct1(k,1)-2.d0*c01(k,1) 
       cp_12(k)= ctt1(k,2)+ct1(k,2)-2.d0*c01(k,2)  
       cp_21(k)= ctt2(k,1)+ct2(k,1)-2.d0*c02(k,1)  
       cp_22(k)= ctt2(k,2)+ct2(k,2)-2.d0*c02(k,2)   
       cp_23(k)= ctt2(k,3)+ct2(k,3)-2.d0*c02(k,3)   
       cp_24(k)= ctt2(k,4)+ct2(k,4)-2.d0*c02(k,4)   
       cp_31(k)= ctt3(k,1)+ct3(k,1)-2.d0*c03(k,1)      
       cp_32(k)= ctt3(k,2)+ct3(k,2)-2.d0*c03(k,2)    
       cp_33(k)= ctt3(k,3)+ct3(k,3)-2.d0*c03(k,3)    
       cp_34(k)= ctt3(k,4)+ct3(k,4)-2.d0*c03(k,4)    
       cp_35(k)= ctt3(k,5)+ct3(k,5)-2.d0*c03(k,5)    
       cp_36(k)= ctt3(k,6)+ct3(k,6)-2.d0*c03(k,6)    
*
       cmh_11(k)= ctt1(k,1)-0.5d0*ct1(k,1)-1.d0/2.d0*c01(k,1) 
       cmh_12(k)= ctt1(k,2)-0.5d0*ct1(k,2)-1.d0/2.d0*c01(k,2)  
       cmh_21(k)= ctt2(k,1)-0.5d0*ct2(k,1)-1.d0/2.d0*c02(k,1)  
       cmh_22(k)= ctt2(k,2)-0.5d0*ct2(k,2)-1.d0/2.d0*c02(k,2)   
       cmh_23(k)= ctt2(k,3)-0.5d0*ct2(k,3)-1.d0/2.d0*c02(k,3)   
       cmh_24(k)= ctt2(k,4)-0.5d0*ct2(k,4)-1.d0/2.d0*c02(k,4)   
       cmh_31(k)= ctt3(k,1)-0.5d0*ct3(k,1)-1.d0/2.d0*c03(k,1)      
       cmh_32(k)= ctt3(k,2)-0.5d0*ct3(k,2)-1.d0/2.d0*c03(k,2)    
       cmh_33(k)= ctt3(k,3)-0.5d0*ct3(k,3)-1.d0/2.d0*c03(k,3)    
       cmh_34(k)= ctt3(k,4)-0.5d0*ct3(k,4)-1.d0/2.d0*c03(k,4)    
       cmh_35(k)= ctt3(k,5)-0.5d0*ct3(k,5)-1.d0/2.d0*c03(k,5)    
       cmh_36(k)= ctt3(k,6)-0.5d0*ct3(k,6)-1.d0/2.d0*c03(k,6)    
*
       cm_11(k)= ctt1(k,1)-ct1(k,1)
       cm_12(k)= ctt1(k,2)-ct1(k,2)
       cm_21(k)= ctt2(k,1)-ct2(k,1)
       cm_22(k)= ctt2(k,2)-ct2(k,2)
       cm_23(k)= ctt2(k,3)-ct2(k,3)
       cm_24(k)= ctt2(k,4)-ct2(k,4)
       cm_31(k)= ctt3(k,1)-ct3(k,1)
       cm_32(k)= ctt3(k,2)-ct3(k,2)
       cm_33(k)= ctt3(k,3)-ct3(k,3)
       cm_34(k)= ctt3(k,4)-ct3(k,4)
       cm_35(k)= ctt3(k,5)-ct3(k,5)
       cm_36(k)= ctt3(k,6)-ct3(k,6)
      enddo
*
*-----3 massless generations
*
      x0h= y*sh
      xh= xv*y*sh
      x2h= x2*y*sh
      yh= by*y*sh
*
      do k=1,2
       w(1,k)= x2h*(-c03(k,2)-c03(k,3)+2.d0*c03(k,4))
     #            +xh*(c01(k,2)+c02(k,2)-c03(k,3)+c03(k,4))
     #            +yh*(c02(k,1)+c02(k,2)-2.d0*c02(k,3)+
     #            c03(k,1)-2.d0*c03(k,3)+c03(k,4))
     #            -2.d0/3.d0*co(k)+2.d0*c02(k,4)-2.d0*c03(k,5)+
     #            2.d0*c03(k,6)+0.5d0*p3qx2(k)
       w(2,k)= x2h*(-2.d0*c02(k,2)+2.d0*c02(k,3)+
     #            3.d0*c03(k,2)+3.d0*c03(k,3)-6.d0*c03(k,4))
     #            +xh*(-c01(k,2)-3.d0*c02(k,2)+2.d0*c02(k,3)+
     #            3.d0*c03(k,3)-3.d0*c03(k,4))+yh*(-2.d0*
     #            c01(k,1)+2.d0*c01(k,2)-5*c02(k,1)-3.d0*c02(k,2)
     #            +8.d0*c02(k,3)-3.d0*c03(k,1)+6.d0*c03(k,3)-
     #            3.d0*c03(k,4))+2.d0*c02(k,4)+6.d0*c03(k,5)-
     #            6.d0*c03(k,6)
       w(3,k)= x2h*(-2.d0*c02(k,2)+2.d0*c02(k,3)-c03(k,2)+
     #            c03(k,3))+xh*(c01(k,2)+c02(k,2)+2.d0*
     #            c02(k,3)+c03(k,3)+c03(k,4))+yh*(-c02(k,1)+
     #            c02(k,2)-c03(k,1)+c03(k,4))+2.d0*c02(k,4)+
     #            2.d0*c03(k,5)+2.d0*c03(k,6)
       w(4,k)= 8.d0*x0h*(-c02(k,2)+c02(k,3)+c03(k,3)-
     #            c03(k,4))
       w(5,k)= 0.d0
       w(6,k)= 0.d0
       w(7,k)= 0.d0
       w(8,k)= x0h*(-4.d0*c01(k,2)-8.d0*c02(k,2)-4.d0*
     #            c02(k,3)-8.d0*c03(k,4))
       w(9,k)= 2.d0*x2h*(c01(k,2)+c02(k,3)-c03(k,2)+
     #            c03(k,4))+2.d0*xh*(c01(k,2)+c02(k,2)+
     #            c02(k,3)+c03(k,4))+2.d0*yh*(-c01(k,1)-
     #            c02(k,1)+c02(k,2)-c02(k,3)-c03(k,3)+c03(k,4))
     #            -2.d0/3.d0*co(k)+4*c02(k,4)+4*c03(k,6)+
     #            0.5d0*p3qx2(k)
       w(10,k)= 0.d0
      enddo
      do if=1,10
       do k=1,2
        w(if,k)= 3.d0*w(if,k)
       enddo
      enddo
*
*-----t-b subtracted
*
      do k=1,2
       w(1,k)= w(1,k)+
     # x2h*(-0.5d0*cph_32(k)-0.5d0*cph_33(k)+cph_34(k))
     # +0.5d0*xh*(cph_12(k)+cph_22(k)-cph_33(k)+cph_34(k))
     # +0.5d0*yh*(cph_21(k)+cph_22(k)-2.d0*cph_23(k)+cph_31(k)
     # -2.d0*cph_33(k)+cph_34(k))
     # +0.5d0*tqm2*(ctt1(k,1)-ctt1(k,2))
     # +0.5d0*ctt0(k)*tqm2 
     # +cph_24(k)-cph_35(k)+cph_36(k)
       w(2,k)= w(2,k)+
     # x2h*(-cph_22(k)+cph_23(k)+1.5d0*cph_32(k)+1.5d0*cph_33(k)
     # -3.d0*cph_34(k))
     # +0.5d0*xh*(-cph_12(k)-3.d0*cph_22(k)+2.d0*cph_23(k)
     # +3.d0*cph_33(k)-3.d0*cph_34(k))
     # +yh*(-cph_11(k)+cph_12(k)-2.5d0*cph_21(k)-1.5d0*cph_22(k)
     # +4.d0*cph_23(k)-1.5d0*cph_31(k)+3.d0*cph_33(k)-1.5d0*cph_34(k))
     # +0.5d0*tqm2*(-ctt1(k,1)+ctt1(k,2))
     # -0.5d0*ctt0(k)*tqm2
     # +cph_24(k)+3.d0*cph_35(k)-3.d0*cph_36(k)
       w(3,k)= w(3,k)+
     # x2h*(-cph_22(k)+cph_23(k)-0.5d0*cph_32(k)+0.5d0*cph_33(k))
     # +0.5d0*xh*(cph_12(k)+cph_22(k)+2.d0*cph_23(k)
     # +cph_33(k)+cph_34(k))
     # +0.5d0*yh*(-cph_21(k)+cph_22(k)-cph_31(k)+cph_34(k))
     # +0.5d0*tqm2*(ctt1(k,1)+ctt1(k,2))
     # +0.5d0*ctt0(k)*tqm2 
     # +cph_24(k)+cph_35(k)+cph_36(k)
       w(4,k)= w(4,k)+
     # 4.d0*x0h*(-cph_22(k)+cph_23(k)+cph_33(k)-cph_34(k))
       w(5,k)= w(5,k)+
     # x2h*(-cmh_32(k)+cmh_33(k))
     # +xh*(cmh_12(k)+cmh_22(k)+2.d0*cmh_23(k)
     # +cmh_33(k)+cmh_34(k))
     # +yh*(-2.d0*cmh_11(k)+2.d0*cmh_12(k)
     # -3.d0*cmh_21(k)+cmh_22(k)
     # +2.d0*cmh_23(k)-cmh_31(k)+cmh_34(k))
     # +tqm2*(ctt1(k,1)+ctt1(k,2))
     # +ctt0(k)*tqm2 
     # +6.d0*cmh_24(k)+6.d0*cmh_35(k)+6.d0*cmh_36(k)
       w(6,k)= w(6,k)+
     # x2h*(-cmh_32(k)-cmh_33(k)+2.d0*cmh_34(k))
     # +xh*(-cmh_12(k)+cmh_22(k)-2.d0*cmh_23(k)
     # -cmh_33(k)+cmh_34(k))
     # +yh*(cmh_21(k)+cmh_22(k)-2.d0*cmh_23(k)
     # +cmh_31(k)-2.d0*cmh_33(k)+cmh_34(k))
     # +tqm2*(-ctt1(k,1)+ctt1(k,2))
     # -ctt0(k)*tqm2
     # -2.d0*cmh_24(k)-6.d0*cmh_35(k)+6.d0*cmh_36(k)
       w(7,k)= w(7,k)+
     # 4.d0*x0h*(cmh_12(k)+cmh_23(k))
       w(8,k)= w(8,k)+
     # x0h*(-2.d0*cmh_12(k)-4.d0*cph_22(k)
     # -2.d0*cph_23(k)-4.d0*cph_34(k))
       w(9,k)= w(9,k)+
     # x2h*(cph_12(k)+cph_23(k)-cph_32(k)+cph_34(k))
     # +xh*(cph_12(k)+cph_22(k)+cph_23(k)+cph_34(k))
     # +yh*(-cmh_11(k)-cph_21(k)+cph_22(k)-cph_23(k)
     # -cph_33(k)+cph_34(k))
     # -tqm2*ctt1(k,2) 
     # +2.d0*cph_24(k)+2.d0*cph_36(k)
       w(10,k)= w(10,k)+
     # 4.d0*x0h*(cmh_12(k)+cmh_23(k))
      enddo
*
*-----additions for Z
*
      do k=1,2
       do if=1,10
        wg(if,k)= w(if,k) 
       enddo
      enddo
      do k=1,2
       dw(1,k)= 0.125d0*x2h*(cm_32(k)+cm_33(k)-2.d0*cm_34(k))
     # +0.125d0*xh*(-cm_12(k)-cm_22(k)+cm_33(k)-cm_34(k))
     # +0.125d0*yh*(-cm_21(k)-cm_22(k)+2.d0*cm_23(k)-cm_31(k)
     # +2.d0*cm_33(k)-cm_34(k))
     # +0.5d0*tqm2*(-ctt1(k,1)+ctt1(k,2))
     # +0.25d0*(-cm_24(k)+cm_35(k)-cm_36(k))
     # -0.5d0*ctt0(k)*tqm2
       dw(2,k)= 0.125d0*x2h*(2.d0*cm_22(k)-2.d0*cm_23(k)
     # -3.d0*cm_32(k)-3.d0*cm_33(k)+6.d0*cm_34(k))
     # +0.125d0*xh*(cm_12(k)+3.d0*cm_22(k)-2*cm_23(k)
     # -3.d0*cm_33(k)+3.d0*cm_34(k))
     # +0.125d0*yh*(2.d0*cm_11(k)-2.d0*cm_12(k)+5*cm_21(k)
     # +3.d0*cm_22(k)-8*cm_23(k)+3.d0*cm_31(k)-6.d0*cm_33(k)
     # +3.d0*cm_34(k))+0.5d0*tqm2*(ctt1(k,1)-ctt1(k,2))
     # +0.25d0*(-cm_24(k)-3.d0*cm_35(k)+3.d0*cm_36(k))
     # +0.5d0*ctt0(k)*tqm2 
       dw(3,k)=  0.125d0*x2h*(2.d0*cm_22(k)-2.d0*cm_23(k)
     # +cm_32(k)-cm_33(k))+0.125d0*xh*(-cm_12(k)-cm_22(k)
     # -2.d0*cm_23(k)-cm_33(k)-cm_34(k))
     # +0.125d0*yh*(cm_21(k)-cm_22(k)+cm_31(k)-cm_34(k))
     # +0.5d0*tqm2*(-ctt1(k,1)-ctt1(k,2))
     # +0.25d0*(-cm_24(k)-cm_35(k)-cm_36(k))
     # -0.5d0*ctt0(k)*tqm2 
       dw(4,k)= x0h*(cm_22(k)-cm_23(k)-cm_33(k)+cm_34(k))
       dw(5,k)= 0.25d0*x2h*(cp_32(k)-cp_33(k))
     # +0.25d0*xh*(-cp_12(k)-cp_22(k)
     # -2.d0*cp_23(k)-cp_33(k)-cp_34(k))
     # +0.25d0*yh*(2.d0*cp_11(k)-2.d0*cp_12(k)+3.d0*cp_21(k)
     # -cp_22(k)-2.d0*cp_23(k)+cp_31(k)-cp_34(k))
     # +tqm2*(-ctt1(k,1)-ctt1(k,2))
     # +1.5d0*(-cp_24(k)-cp_35(k)-cp_36(k))
     # -ctt0(k)*tqm2 
       dw(6,k)= 0.25d0*x2h*(cp_32(k)+cp_33(k)-2.d0*cp_34(k))
     # +0.25d0*xh*(cp_12(k)-cp_22(k)+2.d0*cp_23(k)
     # +cp_33(k)-cp_34(k))
     # +0.25d0*yh*(-cp_21(k)-cp_22(k)+2.d0*cp_23(k)-cp_31(k)
     # +2.d0*cp_33(k)-cp_34(k))
     # +tqm2*(ctt1(k,1)-ctt1(k,2))
     # +0.5d0*( cp_24(k)+3.d0*cp_35(k)-3.d0*cp_36(k))
     # +ctt0(k)*tqm2 
       dw(7,k)= -x0h*(cp_12(k)+cp_23(k))
       dw(8,k)= 0.5d0*x0h*(cm_12(k)+2.d0*cm_22(k)
     # +cm_23(k)+2.d0*cm_34(k))
       dw(9,k)= 0.25d0*x2h*(-cm_12(k)-cm_23(k)
     # +cm_32(k)-cm_34(k))
     # +0.25d0*xh*(-cm_12(k)-cm_22(k)-cm_23(k)-cm_34(k))
     # +0.25d0*yh*(cm_11(k)+cm_21(k)-cm_22(k)+cm_23(k)
     # +cm_33(k)-cm_34(k))+tqm2*ctt1(k,2) 
     # -0.5d0*(cm_24(k)+cm_36(k))
       dw(10,k)= -x0h*(cp_12(k)+cp_23(k))
      enddo
      do k=1,2
       do if=1,10
        wz(if,k)= w(if,k)+dw(if,k)/ctrrx
       enddo
      enddo
      do k=1,2
       do if=1,10
        wg(if,k)= wg(if,k)/apis
        wz(if,k)= wz(if,k)/apis
       enddo
      enddo
*
*-----phi is initialized
*
      if(jj.le.3) then
       cpu= cos(2.d0*pi*phx)
       spu= sin(2.d0*pi*phx)
      else if(jj.ge.4) then
       cpd= cos(2.d0*pi*phx)
       spd= sin(2.d0*pi*phx)
      endif
*
*-----here one starts with u,d variables and then n variables
*
      if(jj.le.3) then
       eu= x1-dmsh
       ctu= (1.d0-2.d0*(tau+x1+xv-umsh-dmsh)/opxv/reu)/betu
       if(abs(ctu).gt.1.d0) then
        iz= 0
        ifz(22)= ifz(22)+1
        go to 1
       endif
       epsu= 1.d0-abs(ctu)
       if(epsu.lt.1.d-5) then
        stu= sqrt(2.d0*epsu)*(1.d0-epsu/4.d0*(1.d0+epsu/8.d0))
       else
        stu= sqrt(1.d0-ctu*ctu)
       endif
       if(stu.eq.0.d0) then
        iz= 0
        ifz(23)= ifz(23)+1
        go to 1
       endif
*
       ed= 1.d0+x2-x1+dmsh
       ctd= (1.d0+2.d0*(tau+x1-x2+z+xv-umsh-dmsh)/opxv/red)/betd
       if(abs(ctd).gt.1.d0) then
        iz= 0
        ifz(24)= ifz(24)+1
        go to 1
       endif
       epsd= 1.d0-abs(ctd)
       if(epsd.lt.1.d-5) then
        std= sqrt(2.d0*epsd)*(1.d0-epsd/4.d0*(1.d0+epsd/8.d0))
       else
        std= sqrt(1.d0-ctd*ctd)
       endif
       aca= betu*betd*stu*std*cpu
       acb= betu*betd*stu*std*spu
       acca= -1.d0+2.d0*(x2-umsh-dmsh)/reu/red
       acc= betu*betd*ctu*ctd+acca
       discazn= acb*acb-acc*acc+aca*aca
       if(discazn.lt.0.d0) then
        iz= 0
        ifz(25)= ifz(25)+1
        go to 1
       endif
       discazn= sqrt(discazn)
*
       if(acc.eq.aca) then
        if(acb.eq.0.d0) then
         iz= 0
         ifz(26)= ifz(26)+1
         go to 1
        else
         byp= -aca/acb
         bym= -aca/acb
        endif
       else
        if(acb.gt.0.d0) then
         acqn= -(acb+discazn)
        else if(acb.lt.0.d0) then
         acqn= -(acb-discazn)
        endif
        byp= acqn/(acc-aca)
        bym= (acc+aca)/acqn
       endif
*
       cpdp= (1.d0-byp*byp)/(1.d0+byp*byp)
       cpdm= (1.d0-bym*bym)/(1.d0+bym*bym)
       spdp= 2.d0*byp/(1.d0+byp*byp)
       spdm= 2.d0*bym/(1.d0+bym*bym)
*
       en= omx2
       ctn= 1.d0-2.d0*(1.d0-x2+z+xv)/opxv/omx2
       actn= -(betu*eu*ctu+betd*ed*ctd)/en
       if(abs(ctn).gt.1.d0) then
        iz= 0
        ifz(27)= ifz(27)+1
        go to 1
       endif
       epsn= 1.d0-abs(ctn)
       if(epsn.lt.1.d-5) then
        stn= sqrt(2.d0*epsn)*(1.d0-epsn/4.d0*(1.d0+epsn/8.d0))
       else
        stn= sqrt(1.d0-ctn*ctn)
       endif
       cpnp= -(betu*eu*stu*cpu+betd*ed*std*cpdp)/en/stn
       spnp= -(betu*eu*stu*spu+betd*ed*std*spdp)/en/stn
       cpnm= -(betu*eu*stu*cpu+betd*ed*std*cpdm)/en/stn
       spnm= -(betu*eu*stu*spu+betd*ed*std*spdm)/en/stn
       tstn= 1.d0-(omx2s-(betu*eu*ctu+betd*ed*ctd)**2)/(en*stn)**2
*
*-----here one starts with d,u variables and then n variables
*
      else if(jj.ge.4) then
       ed= 1.d0+x2-x1+dmsh
       ctd= (1.d0-2.d0*(1.d0+x2-x1+ups+xv)/opxv/red)/betd
       if(abs(ctd).gt.1.d0) then
        iz= 0
        ifz(22)= ifz(22)+1
        go to 1
       endif
       epsd= 1.d0-abs(ctd)
       if(epsd.lt.1.d-5) then
        std= sqrt(2.d0*epsd)*(1.d0-epsd/4.d0*(1.d0+epsd/8.d0))
       else
        std= sqrt(1.d0-ctd*ctd)
       endif
       if(std.eq.0.d0) then
        iz= 0
        ifz(23)= ifz(23)+1
        go to 1
       endif
*
       eu= x1-dmsh
       ctu= (1.d0+2.d0*(omx1+ups+z+xv)/opxv/reu)/betu
       if(abs(ctu).gt.1.d0) then
        iz= 0
        ifz(24)= ifz(24)+1
        go to 1
       endif
       epsu= 1.d0-abs(ctu)
       if(epsu.lt.1.d-5) then
        stu= sqrt(2.d0*epsu)*(1.d0-epsu/4.d0*(1.d0+epsu/8.d0))
       else
        stu= sqrt(1.d0-ctu*ctu)
       endif
       aca= betu*betd*stu*std*cpd
       acb= betu*betd*stu*std*spd
       acca= -1.d0+2.d0*(x2-umsh-dmsh)/reu/red
       acc= betu*betd*ctu*ctd+acca
       discazn= acb*acb-acc*acc+aca*aca
       if(discazn.lt.0.d0) then
        iz= 0
        ifz(25)= ifz(25)+1
        go to 1
       endif
       discazn= sqrt(discazn)
       if(acb.gt.0.d0) then
        acqn= -(acb+discazn)
       else if(acb.lt.0.d0) then
        acqn= -(acb-discazn)
       endif
       byp= acqn/(acc-aca)
       bym= (acc+aca)/acqn
       cpup= (1.d0-byp*byp)/(1.d0+byp*byp)
       cpum= (1.d0-bym*bym)/(1.d0+bym*bym)
       spup= 2.d0*byp/(1.d0+byp*byp)
       spum= 2.d0*bym/(1.d0+bym*bym)
*
       en= omx2
       ctn= 1.d0-2.d0*(1.d0-x2+z+xv)/opxv/omx2
       if(abs(ctn).gt.1.d0) then
        iz= 0
        ifz(26)= ifz(26)+1
        go to 1
       endif
       epsn= 1.d0-abs(ctn)
       if(epsn.lt.1.d-5) then
        stn= sqrt(2.d0*epsn)*(1.d0-epsn/4.d0*(1.d0+epsn/8.d0))
       else
        stn= sqrt(1.d0-ctn*ctn)
       endif
       cpnp= -(betu*eu*stu*cpup+betd*ed*std*cpd)/en/stn
       cpnm= -(betu*eu*stu*cpum+betd*ed*std*cpd)/en/stn
       spnp= -(betu*eu*stu*spup+betd*ed*std*spd)/en/stn
       spnm= -(betu*eu*stu*spum+betd*ed*std*spd)/en/stn
       tstn= 1.d0-(omx2s-(betu*eu*ctu+betd*ed*ctd)**2)/(en*stn)**2
      endif
*
*-----third loop starts here
*
      do kk=1,2
*
      if(kk.eq.1) then
       if(jj.le.3) then
        if(abs(cpdp).le.1.d0.and.abs(cpnp).le.1.d0.and.
     #     abs(spnp).le.1.d0) then
         cpd= cpdp
         spd= spdp
         cpn= cpnp
         spn= spnp
        else
         iz= 0
         ifz(28)= ifz(28)+1
         go to 2
        endif
       else if(jj.ge.4) then
        if(abs(cpup).le.1.d0.and.abs(cpnp).le.1.d0.and.
     #     abs(spnp).le.1.d0) then
         cpu= cpup
         spu= spup
         cpn= cpnp
         spn= spnp
        else
         iz= 0
         ifz(29)= ifz(29)+1
         go to 2
        endif
       endif
      else if(kk.eq.2) then
       if(jj.le.3) then
        if(abs(cpdm).le.1.d0.and.abs(cpnm).le.1.d0.and.
     #     abs(spnm).le.1.d0) then
         cpd= cpdm
         spd= spdm
         cpn= cpnm
         spn= spnm
        else
         iz= 0
         ifz(30)= ifz(30)+1
         go to 2
        endif
       else if(jj.ge.4) then
        if(abs(cpum).le.1.d0.and.abs(cpnm).le.1.d0.and.
     #     abs(spnm).le.1.d0) then
         cpu= cpum
         spu= spum
         cpn= cpnm
         spn= spnm
        else
         iz= 0
         ifz(31)= ifz(31)+1
         go to 2
        endif
       endif
      endif
*
      vfact= xp+xm*y
      vfacts= vfact*vfact
      velp= 1.d0-2.d0*xm*y/vfact+xv*xm*y/vfact*(-6.d0+2.d0*y+
     #      6.d0*xm*y/vfact-2.d0*xm*ys/vfact)+xvs*xm*y/vfact*(
     #      -8.d0+10.d0*y-2.d0*ys+26.d0*xm*y/vfact-22.d0*xm*ys/vfact+ 
     #      4.d0*xm*yc/vfact-18.d0*xms*ys/vfacts+12.d0*xms*yc/vfacts- 
     #      2.d0*xms*yq/vfacts)
      if(abs(velp).gt.1.d0) then
       iz= 0
       ifz(32)= ifz(32)+1
       go to 2
      endif
      betpl= 2.d0*sqrt(xp*xm*y)/vfact*(1.d0+
     #       xv*(3.d0/2.d0-3.d0*xm*y/vfact+xm*ys/vfact-0.5d0*y)+
     #       xvs*(7.d0/8.d0-17.d0/2.d0*xm*y/vfact+8.d0*xm*ys/vfact- 
     #       3.d0/2.d0*xm*yc/vfact+9.d0*xms*ys/vfacts-6.d0*xms*yc/
     #       vfacts+xms*yq/vfacts-7.d0/4.d0*y+3.d0/8.d0*ys))
*
      if(abs(betpl).gt.1.d0) then
       iz= 0
       ifz(33)= ifz(33)+1
       go to 2
      endif
      betl= 1.d0-0.5d0*xv*omy+xvs/8.d0*(-1.d0+6.d0*y-ys)
      if(abs(betl).gt.1.d0) then
       iz= 0
       ifz(34)= ifz(34)+1
       go to 2
      endif
      vels= xv*(omy-xv*y)
      vel= sqrt(vels)
      if(abs(vel).gt.1.d0) then
       iz= 0
       ifz(35)= ifz(35)+1
       go to 2
      endif
*
      ysn= sqrt(y*sh)
*
      vu(1)= 0.5d0*eu*stu*cpu*ysn
      vu(2)= 0.5d0*eu*stu*spu*ysn
      vu(3)= 0.5d0*eu*ctu*ysn
      vu(4)= 0.5d0*eu*ysn
*
      vd(1)= 0.5d0*ed*std*cpd*ysn
      vd(2)= 0.5d0*ed*std*spd*ysn
      vd(3)= 0.5d0*ed*ctd*ysn
      vd(4)= 0.5d0*ed*ysn
*
      vn(1)= 0.5d0*en*stn*cpn*ysn
      vn(2)= 0.5d0*en*stn*spn*ysn
      vn(3)= 0.5d0*en*ctn*ysn
      vn(4)= 0.5d0*en*ysn
*
      vut(1)= vel*(vu(4)-vu(3))-vu(1)
      vut(2)= vu(2)
      vut(3)= (vel*(1.d0+velp)*vu(1)+(vel*vel-1.d0)*vu(3)-
     #        (vel*vel+velp)*vu(4))/betl/betpl
      vut(4)= (-vel*(1.d0+velp)*vu(1)+velp*(1.d0-vel*vel)*vu(3)+
     #        (1.d0+vel*vel*velp)*vu(4))/betl/betpl
*
      vdt(1)= vel*(vd(4)-vd(3))-vd(1)
      vdt(2)= vd(2)
      vdt(3)= (vel*(1.d0+velp)*vd(1)+(vel*vel-1.d0)*vd(3)-
     #        (vel*vel+velp)*vd(4))/betl/betpl
      vdt(4)= (-vel*(1.d0+velp)*vd(1)+velp*(1.d0-vel*vel)*vd(3)+
     #        (1.d0+vel*vel*velp)*vd(4))/betl/betpl
*
      vnt(1)= vel*(vn(4)-vn(3))-vn(1)
      vnt(2)= vn(2)
      vnt(3)= (vel*(1.d0+velp)*vn(1)+(vel*vel-1.d0)*vn(3)-
     #        (vel*vel+velp)*vn(4))/betl/betpl
      vnt(4)= (-vel*(1.d0+velp)*vn(1)+velp*(1.d0-vel*vel)*vn(3)+
     #        (1.d0+vel*vel*velp)*vn(4))/betl/betpl
*
      eut= vut(4)
      ctut= vut(3)/eut
      edt= vdt(4)
      ctdt= vdt(3)/edt
      ent= vnt(4)
      ctnt= vnt(3)/ent
*
      if(ipro.eq.2) then
       eutmn= dmax1(uqm,ae(3))
       edtmn= dmax1(dqm,ae(4))
      else if(ipro.eq.1) then
       edtmn= dmax1(rmm,ae(4))
      else if(ipro.eq.3) then
       edtmn= dmax1(tm,ae(4))
      endif
*
      if(ipro.eq.2) then
       if(eut.lt.eutmn.and.edt.lt.edtmn) then
        iz= 0
        ifz(36)= ifz(36)+1
        go to 2
       endif
       if(eut.lt.eutmn) then
        iz= 0
        ifz(37)= ifz(37)+1
        go to 2
       endif
       if(edt.lt.edtmn) then
        iz= 0
        ifz(38)= ifz(38)+1
        go to 2
       endif
       if(ctut.lt.asa(3).or.ctut.gt.bsa(3)) then
        iz= 0
        ifz(39)= ifz(39)+1
        go to 2
       endif
       if(ctdt.lt.asa(4).or.ctdt.gt.bsa(4)) then
        iz= 0
        ifz(40)= ifz(40)+1
        go to 2
       endif
      else if(ipro.eq.1.or.ipro.eq.3) then
       if(edt.lt.edtmn) then
        iz= 0
        ifz(36)= ifz(36)+1
        go to 2
       endif
       if(ctdt.lt.asa(4).or.ctdt.gt.bsa(4)) then
        iz= 0
        ifz(37)= ifz(37)+1
        go to 2
       endif
       vuts= vut(1)*vut(1)+vut(2)*vut(2)+vut(3)*vut(3)-
     #       vut(4)*vut(4)
       vdts= vdt(1)*vdt(1)+vdt(2)*vdt(2)+vdt(3)*vdt(3)-
     #       vdt(4)*vdt(4)
       hadim= vut(1)*vdt(1)+vut(2)*vdt(2)+vut(3)*vdt(3)-
     #        vut(4)*vdt(4)
       hadim= 2.d0*hadim+vuts+vdts
       hadim= -hadim
       if(hadim.lt.(s0*s)) then
        iz= 0
        ifz(38)= ifz(38)+1
        go to 2
       endif
      endif
*
      if(jj.le.3) then
       tjac= ujc*vjc*yjc*x2jc*x1jc*xjac*taujc*zjc*psjc*pmjc
      else if(jj.ge.4) then
       tjac= ujc*vjc*yjc*x2jc*x1jc*xjac*upsjc*zjc*psjc*pmjc
      endif
*
*-----ingredients for the amplitude squared
*
      gsm2= grx2*grx2+gix2*gix2
      gtm2= gry*gry+giy*giy
      gstm2= gsm2*gtm2
      pwsm2= ((x2*yvv-rsw(1))*(x2*yvv-rsw(1))+srwis)/arhx2m
      pwtm2= ((by*yvv-rsw(1))*(by*yvv-rsw(1))+srwis)/arhym
      rl3= x2*yvv-rsw(1)
*
      crhr= rhry*rhrx2+rhiy*rhix2
      crhi= rhry*rhix2-rhiy*rhrx2 
      ccr= grx2*gry+gix2*giy
      cci= gix2*gry-grx2*giy 
*
*-----special numeration 11 = 2+3, 12= 2-3, 13= 5+6, 14= 5-6
*
      tmpr1= wg(2,1)+wg(3,1)
      tmpr2= wg(2,1)-wg(3,1)
      tmpr3= wg(5,1)+wg(6,1)
      tmpr4= wg(5,1)-wg(6,1)
      tmpi1= wg(2,2)+wg(3,2)
      tmpi2= wg(2,2)-wg(3,2)
      tmpi3= wg(5,2)+wg(6,2)
      tmpi4= wg(5,2)-wg(6,2)
*
      wg(11,1)= tmpr1
      wg(11,2)= tmpi1
      wg(12,1)= tmpr2
      wg(12,2)= tmpi2
      wg(13,1)= tmpr3
      wg(13,2)= tmpi3
      wg(14,1)= tmpr4
      wg(14,2)= tmpi4
*
      tmpr1= wz(2,1)+wz(3,1)
      tmpr2= wz(2,1)-wz(3,1)
      tmpr3= wz(5,1)+wz(6,1)
      tmpr4= wz(5,1)-wz(6,1)
      tmpi1= wz(2,2)+wz(3,2)
      tmpi2= wz(2,2)-wz(3,2)
      tmpi3= wz(5,2)+wz(6,2)
      tmpi4= wz(5,2)-wz(6,2)
*
      wz(11,1)= tmpr1
      wz(11,2)= tmpi1
      wz(12,1)= tmpr2
      wz(12,2)= tmpi2
      wz(13,1)= tmpr3
      wz(13,2)= tmpi3
      wz(14,1)= tmpr4
      wz(14,2)= tmpi4
*
*-----real1 = Re(Gs*Gtc*Pwt*PwsC)
*
      rl1= by*x2*yvvs*crhr*ccr+by*x2*yvvs*crhi*cci 
     #     -rsw(1)*by*yvv*crhr*ccr-rsw(1)*by*yvv*crhi*cci
     #     -rsw(1)*x2*yvv*crhr*ccr-rsw(1)*x2*yvv*crhi*cci 
     #     +rsw(1)*rsw(1)*crhr*ccr+rsw(1)*rsw(1)*crhi*cci
     #     -rsw(2)*by*yvv*crhr*cci+rsw(2)*by*yvv*crhi*ccr
     #     +rsw(2)*x2*yvv*crhr*cci-rsw(2)*x2*yvv*crhi*ccr 
     #     +rsw(2)*rsw(2)*crhr*ccr+rsw(2)*rsw(2)*crhi*cci
*
*-----real2 = Re(Gs*Gtc*Pwt)
*
      rl2= by*yvv*ccr*rhry-by*yvv*cci*rhiy
     #     -rsw(1)*ccr*rhry+rsw(1)*cci*rhiy
     #     +rsw(2)*ccr*rhiy+rsw(2)*cci*rhry
*
*-----Re(w(n)*Gs*PwsC)
*-----Re(w(n)*Gs)
*
      ccrhsr= grx2*rhrx2+gix2*rhix2
      ccrhsi= grx2*rhix2-gix2*rhrx2
      do if=1,14
       rs1wg(if)= x2*yvv*(wg(if,1)*ccrhsr+wg(if,2)*ccrhsi)
     #             +rsw(1)*(-wg(if,1)*ccrhsr-wg(if,2)*ccrhsi)
     #             +rsw(2)*(wg(if,1)*ccrhsi-wg(if,2)*ccrhsr)
       rs2wg(if)= wg(if,1)*grx2-wg(if,2)*gix2
       rs1wz(if)= x2*yvv*(wz(if,1)*ccrhsr+wz(if,2)*ccrhsi)
     #             +rsw(1)*(-wz(if,1)*ccrhsr-wz(if,2)*ccrhsi)
     #             +rsw(2)*(wz(if,1)*ccrhsi-wz(if,2)*ccrhsr)
       rs2wz(if)= wz(if,1)*grx2-wz(if,2)*gix2
      enddo
*
*-----Re(w(n)*Gt*PwtC)
*-----Re(w(n)*Gt)
*
      ccrhtr= gry*rhry+giy*rhiy
      ccrhti= gry*rhiy-giy*rhry
      do if=1,14
       rt1wg(if)= by*yvv*(wg(if,1)*ccrhtr+wg(if,2)*ccrhti)
     #             +rsw(1)*(-wg(if,1)*ccrhtr-wg(if,2)*ccrhti)
     #             +rsw(2)*(wg(if,1)*ccrhti-wg(if,2)*ccrhtr)
       rt2wg(if)= wg(if,1)*gry-wg(if,2)*giy
       rt1wz(if)= by*yvv*(wz(if,1)*ccrhtr+wz(if,2)*ccrhti)
     #             +rsw(1)*(-wz(if,1)*ccrhtr-wz(if,2)*ccrhti)
     #             +rsw(2)*(wz(if,1)*ccrhti-wz(if,2)*ccrhtr)
       rt2wz(if)= wz(if,1)*gry-wz(if,2)*giy
      enddo
*
*-----Re(w(i)*wC(j))
*
      do if=1,14
       do jf=1,14
        rwwg(if,jf)= wg(if,1)*wg(jf,1)+wg(if,2)*wg(jf,2)
        rwwz(if,jf)= wz(if,1)*wz(jf,1)+wz(if,2)*wz(jf,2)
        rwwgz(if,jf)= wg(if,1)*wz(jf,1)+wg(if,2)*wz(jf,2)
       enddo
      enddo
*
*-----------------------------------------------------------------------
*
      vvs= vv*vv
      vvc= vvs*vv
      vvq= vvc*vv
      ys= y*y
      yc= ys*y
      yq= yc*y
*
      bxp= xv+1.d0
      bxps= bxp*bxp
      bxpc= bxps*bxp
      bxpi= 1.d0/bxp
      bxx= xv/bxps
      o= 1.d0-x1
      os= o*o
      oc= os*o
      oq= oc*o
      x12= x1-x2
      x12s= x12*x12
      x12c= x12s*x12
      x2z= x2-z
      x2zs= x2z*x2z
      x2zc= x2zs*x2z
      x2zq= x2zc*x2z
      zi= 1.d0/z
      zis= zi*zi
      if(jj.le.3) then
       tx= tau+xv+x1
       txs= tx*tx
      else if(jj.ge.4) then
       ux= ups+xv+x1
       uxs= ux*ux
      endif
*
*-----running quantities in the X-channel
*
      cst= ctrrx 
      cqt= cst*cst
      csti=  1.d0/cst
      cqti= 1.d0/cqt
      vpe= -0.5d0+strrx
      vme= strrx+0.25d0
      vpd= -0.5d0+strrx/3.d0
      vmd= strrx/3.d0+0.25d0
      vpu= 0.5d0+2.d0/3.d0*strrx
      vmu= -2.d0/3.d0*strrx-0.25d0
      vpes= vpe*vpe
      vmes= vme*vme
      vpds= vpd*vpd
      vmds= vmd*vmd
      vpus= vpu*vpu
      vmus= vmu*vmu

*
*-----The GG amplitude, proportional to g^2, g^3 and g^4
*
      if(jj.eq.1) then
       w1gg= wto_wgg11s()
       w1gg= w1gg+wto_wgg11c()
       w1gg= w1gg+gstm2*wto_wgg11q()
       w2gg= wto_wgg21s()
       w2gg= w2gg+wto_wgg21c()
       w2gg= w2gg+gstm2*wto_wgg21q()
      else if(jj.eq.2) then
       w1gg= wto_wgg12s()
       w1gg= w1gg+wto_wgg12c()
       w2gg= wto_wgg22s()
       w2gg= w2gg+wto_wgg22c()
      else if(jj.eq.3) then
       w1gg= wto_wgg13s()
       w2gg= wto_wgg23s()
      else if(jj.eq.4) then
       w1gg= wto_wgg14s()
       w1gg= w1gg+wto_wgg14c()
       w2gg= wto_wgg24s()
       w2gg= w2gg+wto_wgg24c()
      else if(jj.eq.5) then
       w1gg= wto_wgg15s()
       w1gg= w1gg+wto_wgg15c()
       w2gg= wto_wgg25s()
       w2gg= w2gg+wto_wgg25c()
      else if(jj.eq.6) then
       w1gg= wto_wgg16s()
       w2gg= wto_wgg26s()
      endif
*
      pbts= ((xv*yvv-rsw(1))*(xv*yvv-rsw(1))+srwis)/arhrx
      pbt= -sqrt(pbts)
*
      if(ii.eq.1) then
       gg= erx*erx*(w1gg*psf_gg_11+w2gg*psf_gg_12)
       wcom= gg
      else if(ii.eq.2) then
       gg= erx*erx*(w1gg*psf_gg_21+w2gg*psf_gg_22)
       wcom= gg
      else if(ii.eq.3) then
       gg= erx*erx*(w1gg*psf_gg_31+w2gg*psf_gg_32)
       wcom= gg
      else if(ii.eq.4) then
       gg= erx*erx*(w1gg*psf_gg_41+w2gg*psf_gg_42)
       wcom= gg
      endif
*
      www= y*wcom/pwtm2*arhx2m
*
    2 if(iz.eq.0) then
         ww(kk)= 0.d0
         iz= 1
         inull= inull+1
      else
         ww(kk)= www*tjac*tfactsw*stf/sh
      endif
*
      enddo
*
    1 if(iz.eq.0) then
         wwp= 0.d0
      else
         if(inull.gt.0) then
          wwp= 0.d0
         else
          wwp= ww(1)+ww(2)
         endif
      endif
*
      wwt= wwp
      wwtt= wwt
*
      if(inull.gt.0) then
       wtoxsswcfl= 0.d0
      else
       wtoxsswcfl= wwtt
      endif
*
      return
      end