*
*-----WTOXSSW------------------------------------------------------------
*
      real*8 function wtoxssw(ndim,x)
      implicit real*8(a-h,o-z)
*
      common/wtkount/ik
      common/wtistrf/isf
      common/wtipt/ifz(51)
      common/wtfsw/tfactsw
      common/wtnpr/ipr,ipr0
      common/wtcutsw/tcs,s0cut
      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/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
*
      dimension x(ndim),ww(5)
*
      if(ipr.eq.4) then
         ipro= 1
      else if(ipr.eq.5) then
         ipro= 2
      else if(ipr.eq.39) then
         ipro= 3
      endif
      do ii=1,5
       ww(ii)= 0.d0
      enddo
      em2= em*em
      em4= em2*em2
      em6= em4*em2
      s2= s*s
      s3= s2*s
      rem= em/sqrt(s)
      rwms= rwm2
      rwmq= rwms*rwms
      rwgs= rwg*rwg
      rwgm= rwg*rwm
      rwgms= rwgs*rwms
      rwmpgs= rwms+rwgs
      s0= s0cut
*
      if(ndim.eq.4) then
       ymx= x(1)
       ypx= x(2)
       zx= x(3)
       xx= x(4)
      else if(ndim.eq.6) then
       uvx= x(1)
       vvx= x(2)
       ymx= x(3)
       ypx= x(4)
       zx= x(5)
       xx= x(6)
      endif
*
      iz= 1
      ik= ik+1
*
      if(ndim.eq.4) 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 2
       endif
       vvl= 0.25d0*(rem+sqrt(rem*rem+4.d0*s0))*
     #      (rem+sqrt(rem*rem+4.d0*s0))
       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 2
       else
          vjc= vjc0**hbet
       endif
      endif
*
      if(vv.lt.0.d0) then
         iz= 0
         ifz(2)= ifz(2)+1
         go to 2
      endif
      svv= sqrt(vv)
      xm= uv
      xp= vv/uv  
      xmop= xm/xp
      sh= vv*s
      shs= sh*sh
      ssh= sqrt(sh)
      vvs= vv*vv
      vvi= 1.d0/vv
      vvis= vvi*vvi
      vvic= vvis*vvi
      vviq= vvic*vvi
*
      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
*
*-----ym is initialized
*
      yml= (4.d0*s0+tcs)/(4.d0+tcs)/vv
      ymu= 1.d0-em/ssh
      if(yml.gt.ymu) then
         iz= 0
         ifz(3)= ifz(3)+1
         go to 2
      endif
*
      ymjc= ymu-yml
      ym= ymjc*ymx+yml
      yms= ym*ym
      ymvv= ym*vv
      ymvvs= ymvv*ymvv
      ymvvc= ymvvs*ymvv
      ymvvq= ymvvc*ymvv
*      omym= (1.d0-s0/vv)*(1.d0-ymx)+em/ssh*ymx
      omym= 1.d0-ym
      omyms= omym*omym
      omymc= omyms*omym
*
*-----yp is initialized
*
      ypl= s0/ymvv+zv
      ypu= 1.d0
      if(ypl.gt.ypu) then
         iz= 0
         ifz(4)= ifz(4)+1
         go to 2
      endif
      ypjc= ypu-ypl
      yp= ypjc*ypx+ypl
      ypi= 1.d0/yp
      ypis= ypi*ypi
      ypic= ypis*ypi
      ypiq= ypic*ypi
      ypiv= ypiq*ypi
      ypivi= ypiv*ypi
      omyp= 1.d0-yp
      omyps= omyp*omyp
      omypc= omyps*omyp
*
      qh0s= em2*(xm*ym+omxm)*ym/omym
      qhcs= qh0s+0.25d0*xm*xm*omym*tcs*s
*
      z0= qh0s/ym/sh
      zc0= qhcs/ym/sh
      zc1= 1.d0-s0/ymvv
      zc2= yp-s0/ymvv
      zc= dmin1(zc0,zc1,zc2)
      if(zc.eq.zc0) then
       dzc= 0.25d0*xm*xm*omym*tcs*s/ym/sh
      else
       dzc= zc-z0
      endif
*
      if(z0.gt.zc) then
         iz= 0
         ifz(5)= ifz(5)+1
         go to 2
      endif
*
*-----Loop for integration starts here
*     z is initialized
*
      do ii=1,3
*
      if(ii.eq.1) then
       zjac= dzc/zc/z0/ym/sh
       zv= z0*zc/(dzc*zx+z0)
      else if(ii.eq.2) then
       zjac= log(zc/z0)
       zv= exp(zjac*zx+log(z0))
      else if(ii.eq.3) then
       zjac= dzc
       zv= zjac*zx+z0
      else if(ii.eq.4) then
       zjac= dzc/zc/z0/(ym*sh)**2
       zv= z0*zc/(dzc*zx+z0)
      else if(ii.eq.5) then
       zjac= log(zc/z0)/ym/sh
       zv= exp(zjac*zx+log(z0))
      endif
*
*-----x is initialized
*
      if(ipro.eq.2) then
       ums= uqm*uqm/s
       dms= dqm*dqm/s
       rum= sqrt(ums)
       rdm= sqrt(dms)
       tha= (rum-rdm)*(rum-rdm)
       thn= (rum+rdm)*(rum+rdm)
      else if(ipro.eq.1) then
       dms= rmm*rmm/s
      else if(ipro.eq.3) then
       dms= tm*tm/s
      endif
*
      xl1= s0/ymvv+zv
      xl2= (1.d0-0.5d0/zv)*yp
      xl3= thn/ymvv+zv
      xl= dmax1(xl1,xl2,xl3)
      xu= yp
      if(xl.gt.xu) then
         iz= 0
         ifz(6)= ifz(6)+1
         go to 1
      endif
*
      atma= ((xl-zv)*ymvv-rwms)/rwgm
      atmb= ((xu-zv)*ymvv-rwms)/rwgm
*
      if(atma.gt.1.d0.and.atmb.gt.1.d0) then
         atma= 1.d0/atma
         atma= atan(atma)
         zmat= (pih-atma)/rwgm
         atmb= 1.d0/atmb
         atmb= atan(atmb)
         zmbt= (pih-atmb)/rwgm
         xjac0= (-atmb+atma)/rwgm
      else if(atma.gt.1.d0.and.atmb.lt.-1.d0) then
         atma= 1.d0/atma
         atma= atan(atma)
         zmat= (pih-atma)/rwgm
         atmb= -1.d0/atmb
         atmb= atan(atmb)
         zmbt= (-pih+atmb)/rwgm
         xjac0= (-pi+atmb+atma)/rwgm
      else if(atma.gt.1.d0.and.abs(atmb).lt.1.d0) then
         atma= 1.d0/atma
         atma= atan(atma)
         zmat= (pih-atma)/rwgm
         atmb= atan(atmb)
         zmbt= atmb/rwgm
         xjac0= (-pih+atmb+atma)/rwgm
      else if(atma.lt.-1.d0.and.atmb.gt.1.d0) then
         atma= -1.d0/atma
         atma= atan(atma)
         zmat= (-pih+atma)/rwgm
         atmb= 1.d0/atmb
         atmb= atan(atmb)
         zmbt= (pih-atmb)/rwgm
         xjac0= (pi-atmb-atma)/rwgm
      else if(atma.lt.-1.d0.and.atmb.lt.-1.d0) then
         atma= -1.d0/atma
         atma= atan(atma)
         zmat= (-pih+atma)/rwgm
         atmb= -1.d0/atmb
         atmb= atan(atmb)
         zmbt= (-pih+atmb)/rwgm
         xjac0= (atmb-atma)/rwgm
      else if(atma.lt.-1.d0.and.abs(atmb).lt.1.d0) then
         atma= -1.d0/atma
         atma= atan(atma)
         zmat= (-pih+atma)/rwgm
         atmb= atan(atmb)
         zmbt= atmb/rwgm
         xjac0= (pih+atmb-atma)/rwgm
      else if(abs(atma).lt.1.d0.and.atmb.gt.1.d0) then
         atma= atan(atma)
         zmat= atma/rwgm
         atmb= 1.d0/atmb
         atmb= atan(atmb)
         zmbt= (pih-atmb)/rwgm
         xjac0= (pih-atmb-atma)/rwgm
      else if(abs(atma).lt.1.d0.and.atmb.lt.-1.d0) then
         atma= atan(atma)
         zmat= atma/rwgm
         atmb= -1.d0/atmb
         atmb= atan(atmb)
         zmbt= (-pih+atmb)/rwgm
         xjac0= (-pih+atmb-atma)/rwgm
      else if(abs(atma).lt.1.d0.and.abs(atmb).lt.1.d0) then
         atma= atan(atma)
         zmat= atma/rwgm
         atmb= atan(atmb)
         zmbt= atmb/rwgm
         xjac0= (atmb-atma)/rwgm
      endif
*
      xjac= xjac0/ymvv
      argx= rwgm*(xjac0*xx+zmat)
      iftn= 1
      s07= s07aaf(argx,iftn)
      if(iftn.ne.0) then
       print*,' call to s07aaf, IFAIL = ',iftn
      endif
      xv= rwms/ymvv+zv+rwgm/ymvv*s07
      xvs= xv*xv
      xvc= xvs*xv
      xvq= xvc*xv
      xvv= xvq*xv
      xvvi= xvv*xv
*
      pwtm2= ((yp-zv)*ymvv-rwgm*s07)**2+rwms*rwgs
*
      qms= zv*ymvv
      qmss= qms*qms
      qps= (yp-xv)*ymvv
      qpm= -0.5d0*yp*ymvv
      qpms= qpm*qpm
      qra= qps/qpm
      qras= qra*qra
      ophqra= 1.d0+0.5d0*qra
      opqra= 1.d0+qra
      tpqra= 2.d0+qra
      ophqras= ophqra*ophqra
      bqs0= qps+2.d0*qpm
      bqs= bqs0*(1.d0+qms/bqs0)
      qex= qms/qpm
      qexs= qex*qex
      qexx= qms*qex
      qt= qps*qms/qpms
      grm2= qpms*(1.d0-qt)
      grm= sqrt(grm2)
      opgrm= 1.d0+grm/qpm
      omgrm= 1.d0-grm/qpm
      ophqrae= ophqra+0.25d0*qras
*
      if(ipro.eq.2) then
       umds= ums-dms
       dmus= dms-ums
       upds= ums+dms
*
       if(ums.ne.0.d0.and.dms.ne.0.d0) then
        if(qt.lt.1.d-4) then
         denu= vvs*(bqs0+qms-umds)/(4.d0*ums*qpms-
     #         bqs0*qps*qms-qps*qmss-2.d0*qps*qms*
     #         upds+2.d0*qpm*qms*umds)
         dend= vvs*(bqs0+qms+umds)/(4.d0*dms*qpms-
     #         bqs0*qps*qms-qps*qmss-2.d0*qps*qms*
     #         upds-2.d0*qpm*qms*umds)
         denu= (1.d0+2.d0*bqs*ums/(bqs+umds)**2)*denu
         dend= (1.d0+2.d0*bqs*dms/(bqs-umds)**2)*dend
*
         unum= qpm*tpqra+2.d0*dms-ums+
     #         qms*(1.d0-0.5d0*qra*ophqra)-0.25d0*qra*
     #         qexx*ophqrae-0.5d0*ums*
     #         qex+0.5d0*dms*(1.d0-qra)*qex
         uden= -ums+0.5d0*qms*qra*ophqra+0.25d0*qra*
     #         qexx*ophqrae-0.5d0*ums*
     #         qex+0.5d0*dms*opqra*qex
         rlncu= log(unum/uden)
         dnum= qpm*tpqra+2.d0*ums-dms+
     #         qms*(1.d0-0.5d0*qra*ophqra)-0.25d0*qra*
     #         qexx*ophqrae-0.5d0*dms*
     #         qex+0.5d0*ums*(1.d0-qra)*qex
         dden= -dms+0.5d0*qms*qra*ophqra+0.25d0*qra*
     #         qexx*ophqrae-0.5d0*dms*
     #         qex+0.5d0*ums*opqra*qex
         rlncd= log(dnum/dden)
        else
         denu= vvs*(bqs-umds)/(4.d0*ums*qpms-qms*qpm*(
     #         qra*tpqra*qpm+2.d0*opqra*dmus+qra*qms))
         dend= vvs*(bqs+umds)/(4.d0*dms*qpms-qms*qpm*(
     #         qra*tpqra*qpm+2.d0*opqra*umds+qra*qms))
         unum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+
     #         dmus/qpm)+dms*opgrm-ums
         uden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+
     #         dmus/qpm)+dms*omgrm-ums
         dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+
     #         umds/qpm)+ums*opgrm-dms
         dden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+
     #         umds/qpm)+ums*omgrm-dms
         rlncu= -log(unum/uden)
         rlncd= -log(dnum/dden)
        endif
       else 
        if(qt.lt.1.d-4) then
         denu= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms))
         dend= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms))
         unum= qpm*tpqra+qms*(1.d0-0.5d0*qra*ophqra)-
     #         0.25d0*qexx*qra*ophqras
         uden= 0.5d0*qms*qra*(ophqra+0.5d0*qex*ophqras+
     #         qra/8.d0*qexs*(opqra*ophqra))
         if((unum/uden).le.0.d0) then
          iz= 0
          ifz(6)= ifz(6)+1
          go to 1
         endif
         rlncu= log(unum/uden)
         rlncd= rlncd
        else
         denu= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms))
         dend= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms))
         unum= ophqra*qpm+grm*ophqra+0.5d0*qms*opgrm
         uden= ophqra*(qpm-grm)+0.5d0*qms*omgrm
         dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*opgrm
         dden= ophqra*(qpm-grm)+0.5d0*qms*omgrm
         if((unum/uden).le.0.d0) then
          iz= 0
          ifz(6)= ifz(6)+1
          go to 1
         endif
         rlncu= -log(unum/uden)
         rlncd= -log(dnum/dden)
        endif
       endif
*
      else if(ipro.eq.1.or.ipro.eq.3) then
       umds= -dms
       dmus= dms
       upds= dms
*
       if(qt.lt.1.d-4) then
        dend= vvs*(bqs0+qms+umds)/(4.d0*dms*qpms-
     #        bqs0*qps*qms-qps*qmss-2.d0*qps*qms*
     #        upds-2.d0*qpm*qms*umds)
        dend= (1.d0+2.d0*bqs*dms/(bqs-umds)**2)*dend
*
        dnum= qpm*tpqra-dms+qms*(1.d0-0.5d0*qra*ophqra)-
     #        0.25d0*qra*qexx*ophqrae-0.5d0*dms*
     #        qex+0.5d0*ums*(1.d0-qra)*qex
        dden= -dms+0.5d0*qms*qra*ophqra+0.25d0*qra*
     #        qexx*ophqrae-0.5d0*dms*qex+0.5d0*ums*opqra*qex
        rlncd= log(dnum/dden)
       else
        dend= vvs*(bqs+umds)/(4.d0*dms*qpms-qms*qpm*(
     #        qra*tpqra*qpm+2.d0*opqra*umds+qra*qms))
        dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+
     #        umds/qpm)-dms
        dden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+
     #        umds/qpm)-dms
        rlncd= -log(dnum/dden)
       endif
      endif
*
      bd= -xp*xm+(xp+xm)*em2/s
      qyjc= 1.d0-em4/s2*(1.d0+(xp+xm*omxm)/bd)+em6/s3/bd*(xp+xm) 
*
      if(ii.eq.1) then
       psf= em2*(2.d0*ymvvs+omxm*(2.d0*ym*vvs-ymvvs))
      else if(ii.eq.2) then
       psf= -ymvvs+2.d0*vvs*(ym-1.d0)
      else if(ii.ge.3) then
       psf= ymvvs
      endif
*
      tjac0= ujc*vjc*ymjc*ypjc*zjac*xjac/qyjc
      tjac= tjac0*stf*psf/sh
      rmu2= -bqs
      rlk= (rmu2-tha)*(rmu2-thn)
      tjac= tjac*sqrt(rlk)/rmu2
*
      if(ipro.eq.2.and.ii.le.2) then
*
      wcomu1= ypic*rwms*(-8./9.*xvq)
      wcomu1= wcomu1+ypis*rwms*(16./9.*xvc)
      wcomu1= wcomu1+ypi*rwms*(-4./3.*xvs)
      wcomu1= wcomu1+rwms*(4./9.*xv)
      wcomu1= wcomu1+ymvv*ypic*(4./9.*xvv)
      wcomu1= wcomu1+ymvv*ypis*(-8./9.*xvq)
      wcomu1= wcomu1+ymvv*ypi*(2./3.*xvc)
      wcomu1= wcomu1+ymvv*(-2./9.*xvs)
      wcomu1= ypi*ymvv*wcomu1
*
      wcomu2= ypic*(4./9.*xvc)
      wcomu2= wcomu2+ypis*(-8./9.*xvs)
      wcomu2= wcomu2+ypi*(2./3.*xv)
      wcomu2= wcomu2-2./9.
      wcomu2= ypi*rwms*rwmpgs*wcomu2
*
      wcomu= wcomu1+wcomu2
*
      wcomd1= ypiq*rwms*(-2./9.*xvq)
      wcomd1= wcomd1+ypic*rwms*xvc*(4./9.+4./9.*xv)
      wcomd1= wcomd1+ypis*rwms*xvs*(-1./3.-8./9.*xv-2./9.*xvs)
      wcomd1= wcomd1+ypi*rwms*xv*(1./9.+2./3.*xv+4./9.*xvs)
      wcomd1= wcomd1+rwms*omyp*(-1./9.*xv)
      wcomd1= wcomd1+rwms*xv*(-1./9.-1./3.*xv)
      wcomd1= wcomd1+ymvv*ypiq*(1./9.*xvv)
      wcomd1= wcomd1+ymvv*ypic*xvq*(-2./9.-2./9.*xv)
      wcomd1= wcomd1+ymvv*ypis*xvc*(1./6.+4./9.*xv+1./9.*xvs)
      wcomd1= wcomd1+ymvv*ypi*xvs*(-1./18.-1./3.*xv-2./9.*xvs)
      wcomd1= wcomd1+ymvv*omyp*(1./18.*xvs)
      wcomd1= wcomd1+ymvv*xvs*(1./18.+1./6.*xv)
      wcomd1= ymvv*wcomd1
*
      wcomd2= ypiq*(1./9.*xvc)
      wcomd2= wcomd2+ypic*xvs*(-2./9.-2./9.*xv)
      wcomd2= wcomd2+ypis*xv*(1./6.+4./9.*xv+1./9.*xvs)
      wcomd2= wcomd2+ypi*(-1./18.-1./3.*xv-2./9.*xvs)
      wcomd2= wcomd2+omyp*(1./18.)
      wcomd2= wcomd2+(1./18.+1./6.*xv)
      wcomd2= rwms*rwmpgs*wcomd2
*
      wcomd= wcomd1+wcomd2
*
      wcom01= ypiq*rwms*(14./3.*xvq)
      wcom01= wcom01+ypic*rwms*(-28./3.*xvc-10./3.*xvq)
      wcom01= wcom01+ypis*rwms*xvs*(62./9.+20./3.*xv+11./9.*xvs)
      wcom01= wcom01+ypi*rwms*xv*(-20./9.-50./9.*
     #  xv-22./9.*xvs-2./3.*xvc)
      wcom01= wcom01+rwms*omyp*xv*(11./18.+xv)
      wcom01= wcom01+rwms*omyps*(1./3.*xv)
      wcom01= wcom01+rwms*xv*(23./18.+3./2.*xv+4./3.*xvs)
      wcom01= wcom01+ymvv*ypiq*(-7./3.*xvv)
      wcom01= wcom01+ymvv*ypic*xvq*(5.+5./3.*xv)
      wcom01= wcom01+ymvv*ypis*xvc*(-71./18.-11./3.*xv-11./18.*xvs)
      wcom01= wcom01+ymvv*ypi*xvs*(29./18.+59./18.*xv
     #  +14./9.*xvs+1./3.*xvc)
      wcom01= wcom01+ymvv*omyp*xv*(-1./6.-5./36.*xv-7./6.*xvs)
      wcom01= wcom01+ymvv*omyps*(-2./3.*xvs)
      wcom01= wcom01+ymvv*omypc*(-1./6.*xv)
      wcom01= wcom01+ymvv*xvs*(-41./36.-3./4.*xv-xvs)
      wcom01= ymvv*wcom01
*
      wcom02= ypiq*(-7./3.*xvc)
      wcom02= wcom02+ypic*xvs*(13./3.+5./3.*xv)
      wcom02= wcom02+ypis*xv*(-53./18.-3*xv-11./18.*xvs)
      wcom02= wcom02+ypi*(11./18.+41./18.*xv+8./9.*xvs+1./3.*xvc)
      wcom02= wcom02+omyp*(-5./36.-1./6.*xv)
      wcom02= wcom02+(-5./36.-3./4.*xv-1./3.*xvs)
      wcom02= wcom02*rwms*rwmpgs
*
      wcom0= wcom01+wcom02
*
      denum= denu*ums/vv
      dendm= dend*dms/vv
*
      wcomum= ypis/vv*rwms*rwmpgs*(4./9.*xvs)
      wcomum= wcomum+ypi/vv*rwms*rwmpgs*(-8./9.*xv)
      wcomum= wcomum+1./vv*rwms*rwmpgs*(4./9.)
      wcomum= wcomum+ymvv*ypis/vv*rwms*(-8./9.*xvc)
      wcomum= wcomum+ymvv*ypi/vv*rwms*(16./9.*xvs)
      wcomum= wcomum+ymvv/vv*rwms*(-8./9.*xv)
      wcomum= wcomum+ymvvs*ypis/vv*(4./9.*xvq)
      wcomum= wcomum+ymvvs*ypi/vv*(-8./9.*xvc)
      wcomum= wcomum+ymvvs/vv*(4./9.*xvs)
*
      wcomdm= ypis/vv*rwms*rwmpgs*(1./9.*xvs)
      wcomdm= wcomdm+ypi/vv*rwms*rwmpgs*xv*(-2./9.-2./9.*xv)
      wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyp*(2./9.*xv)
      wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyps*(1./9.)
      wcomdm= wcomdm+1./vv*rwms*rwmpgs*xv*(2./9.+1./9.*xv)
      wcomdm= wcomdm+ymvv*ypis/vv*rwms*(-2./9.*xvc)
      wcomdm= wcomdm+ymvv*ypi/vv*rwms*xvs*(4./9.+4./9.*xv)
      wcomdm= wcomdm+ymvv/vv*rwms*omyp*(-4./9.*xvs)
      wcomdm= wcomdm+ymvv/vv*rwms*omyps*(-2./9.*xv)
      wcomdm= wcomdm+ymvv/vv*rwms*xvs*(-4./9.-2./9.*xv)
      wcomdm= wcomdm+ymvvs*ypis/vv*(1./9.*xvq)
      wcomdm= wcomdm+ymvvs*ypi/vv*xvc*(-2./9.-2./9.*xv)
      wcomdm= wcomdm+ymvvs/vv*omyp*(2./9.*xvc)
      wcomdm= wcomdm+ymvvs/vv*omyps*(1./9.*xvs)
      wcomdm= wcomdm+ymvvs/vv*xvc*(2./9.+1./9.*xv)
*
      wcom= wcomu*rlncu+wcomd*rlncd+wcom0+
     #      ymvv*(wcomum*denum+wcomdm*dendm) 
*
      else if(ipro.eq.2.and.ii.ge.3) then
*
      wcomu1= ypic*rwms*(-8./9.*xvq)
      wcomu1= wcomu1+ypis*rwms*(16./9.*xvc)
      wcomu1= wcomu1+ypi*rwms*(-4./3.*xvs)
      wcomu1= wcomu1+rwms*(4./9.*xv)
      wcomu1= wcomu1+ymvv*ypic*(4./9.*xvv)
      wcomu1= wcomu1+ymvv*ypis*(-8./9.*xvq)
      wcomu1= wcomu1+ymvv*ypi*(2./3.*xvc)
      wcomu1= wcomu1+ymvv*(-2./9.*xvs)
      wcomu1= ypi*ymvv*wcomu1
*
      wcomu2= ypic*(4./9.*xvc)
      wcomu2= wcomu2+ypis*(-8./9.*xvs)
      wcomu2= wcomu2+ypi*(2./3.*xv)
      wcomu2= wcomu2-2./9.
      wcomu2= ypi*rwms*rwmpgs*wcomu2
*
      wcomu= wcomu1+wcomu2
*
      wcomd1= ypiq*rwms*(-2./9.*xvq)
      wcomd1= wcomd1+ypic*rwms*xvc*(4./9.+4./9.*xv)
      wcomd1= wcomd1+ypis*rwms*xvs*(-1./3.-8./9.*xv-2./9.*xvs)
      wcomd1= wcomd1+ypi*rwms*xv*(1./9.+2./3.*xv+4./9.*xvs)
      wcomd1= wcomd1+rwms*omyp*(-1./9.*xv)
      wcomd1= wcomd1+rwms*xv*(-1./9.-1./3.*xv)
      wcomd1= wcomd1+ymvv*ypiq*(1./9.*xvv)
      wcomd1= wcomd1+ymvv*ypic*xvq*(-2./9.-2./9.*xv)
      wcomd1= wcomd1+ymvv*ypis*xvc*(1./6.+4./9.*xv+1./9.*xvs)
      wcomd1= wcomd1+ymvv*ypi*xvs*(-1./18.-1./3.*xv-2./9.*xvs)
      wcomd1= wcomd1+ymvv*omyp*(1./18.*xvs)
      wcomd1= wcomd1+ymvv*xvs*(1./18.+1./6.*xv)
      wcomd1= ymvv*wcomd1
*
      wcomd2= ypiq*(1./9.*xvc)
      wcomd2= wcomd2+ypic*xvs*(-2./9.-2./9.*xv)
      wcomd2= wcomd2+ypis*xv*(1./6.+4./9.*xv+1./9.*xvs)
      wcomd2= wcomd2+ypi*(-1./18.-1./3.*xv-2./9.*xvs)
      wcomd2= wcomd2+omyp*(1./18.)
      wcomd2= wcomd2+(1./18.+1./6.*xv)
      wcomd2= rwms*rwmpgs*wcomd2
*
      wcomd= wcomd1+wcomd2
*
      wcom01= ypiq*rwms*(14./3.*xvq)
      wcom01= wcom01+ypic*rwms*(-28./3.*xvc-10./3.*xvq)
      wcom01= wcom01+ypis*rwms*xvs*(62./9.+20./3.*xv+11./9.*xvs)
      wcom01= wcom01+ypi*rwms*xv*(-20./9.-50./9.*
     #  xv-22./9.*xvs-2./3.*xvc)
      wcom01= wcom01+rwms*omyp*xv*(11./18.+xv)
      wcom01= wcom01+rwms*omyps*(1./3.*xv)
      wcom01= wcom01+rwms*xv*(23./18.+3./2.*xv+4./3.*xvs)
      wcom01= wcom01+ymvv*ypiq*(-7./3.*xvv)
      wcom01= wcom01+ymvv*ypic*xvq*(5.+5./3.*xv)
      wcom01= wcom01+ymvv*ypis*xvc*(-71./18.-11./3.*xv-11./18.*xvs)
      wcom01= wcom01+ymvv*ypi*xvs*(29./18.+59./18.*xv
     #  +14./9.*xvs+1./3.*xvc)
      wcom01= wcom01+ymvv*omyp*xv*(-1./6.-5./36.*xv-7./6.*xvs)
      wcom01= wcom01+ymvv*omyps*(-2./3.*xvs)
      wcom01= wcom01+ymvv*omypc*(-1./6.*xv)
      wcom01= wcom01+ymvv*xvs*(-41./36.-3./4.*xv-xvs)
      wcom01= ymvv*wcom01
*
      wcom02= ypiq*(-7./3.*xvc)
      wcom02= wcom02+ypic*xvs*(13./3.+5./3.*xv)
      wcom02= wcom02+ypis*xv*(-53./18.-3*xv-11./18.*xvs)
      wcom02= wcom02+ypi*(11./18.+41./18.*xv+8./9.*xvs+1./3.*xvc)
      wcom02= wcom02+omyp*(-5./36.-1./6.*xv)
      wcom02= wcom02+(-5./36.-3./4.*xv-1./3.*xvs)
      wcom02= wcom02*rwms*rwmpgs
*
      wcom0= wcom01+wcom02
*
      denum= denu*ums/vv
      dendm= dend*dms/vv
*
      wcomum= ypis/vv*rwms*rwmpgs*(4./9.*xvs)
      wcomum= wcomum+ypi/vv*rwms*rwmpgs*(-8./9.*xv)
      wcomum= wcomum+1./vv*rwms*rwmpgs*(4./9.)
      wcomum= wcomum+ymvv*ypis/vv*rwms*(-8./9.*xvc)
      wcomum= wcomum+ymvv*ypi/vv*rwms*(16./9.*xvs)
      wcomum= wcomum+ymvv/vv*rwms*(-8./9.*xv)
      wcomum= wcomum+ymvvs*ypis/vv*(4./9.*xvq)
      wcomum= wcomum+ymvvs*ypi/vv*(-8./9.*xvc)
      wcomum= wcomum+ymvvs/vv*(4./9.*xvs)
*
      wcomdm= ypis/vv*rwms*rwmpgs*(1./9.*xvs)
      wcomdm= wcomdm+ypi/vv*rwms*rwmpgs*xv*(-2./9.-2./9.*xv)
      wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyp*(2./9.*xv)
      wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyps*(1./9.)
      wcomdm= wcomdm+1./vv*rwms*rwmpgs*xv*(2./9.+1./9.*xv)
      wcomdm= wcomdm+ymvv*ypis/vv*rwms*(-2./9.*xvc)
      wcomdm= wcomdm+ymvv*ypi/vv*rwms*xvs*(4./9.+4./9.*xv)
      wcomdm= wcomdm+ymvv/vv*rwms*omyp*(-4./9.*xvs)
      wcomdm= wcomdm+ymvv/vv*rwms*omyps*(-2./9.*xv)
      wcomdm= wcomdm+ymvv/vv*rwms*xvs*(-4./9.-2./9.*xv)
      wcomdm= wcomdm+ymvvs*ypis/vv*(1./9.*xvq)
      wcomdm= wcomdm+ymvvs*ypi/vv*xvc*(-2./9.-2./9.*xv)
      wcomdm= wcomdm+ymvvs/vv*omyp*(2./9.*xvc)
      wcomdm= wcomdm+ymvvs/vv*omyps*(1./9.*xvs)
      wcomdm= wcomdm+ymvvs/vv*xvc*(2./9.+1./9.*xv)
*
      wcom= wcomu*rlncu+wcomd*rlncd+wcom0+
     #      ymvv*(wcomum*denum+wcomdm*dendm) 
*
      dw1= ymvv*ypivi*rlncu*vvis*rwms*(56./3.*xvv)
      dw1= dw1+ymvv*ypivi*rlncu*rwms*(64./3.*xvv)
      dw1= dw1+ymvv*ypivi*rlncd*vvis*rwms*(16./3.*xvv)
      dw1= dw1+ymvv*ypivi*rlncd*rwms*(8*xvv)
      dw1= dw1+ymvv*ypivi*vvis*rwms*(-500./9.*xvv)
      dw1= dw1+ymvv*ypivi*rwms*(-232./3.*xvv)
      dw1= dw1+ymvv*ypiv*rlncu*vvis*rwms*(-48*xvq-136./9.*xvv)
      dw1= dw1+ymvv*ypiv*rlncu*rwms*(-64*xvq-40./3.*xvv)
      dw1= dw1+ymvv*ypiv*rlncd*vvis*rwms*(-40./3.*
     #  xvq-104./9.*xvv)
      dw1= dw1+ymvv*ypiv*rlncd*rwms*(-24.*xvq-16.*xvv)
      dw1= dw1+ymvv*ypiv*vvis*rwms*(400./3.*xvq+680./9.*xvv)
      dw1= dw1+ymvv*ypiv*rwms*(232.*xvq+260./3.*xvv)
      dw1= dw1+ymvv*ypiq*rlncu*vvis*rwms*(362./9.*
     #  xvc+364./9.*xvq)
      dw1= dw1+ymvv*ypiq*rlncu*rwms*(226./3.*xvc+116./3.*xvq)
      dw1= dw1+ymvv*ypiq*rlncd*vvis*rwms*(95./9.*
     #  xvc+278./9.*xvq+8*xvv)
      dw1= dw1+ymvv*ypiq*rlncd*rwms*(85./3.*xvc+142./
     #  3.*xvq+28./3.*xvv)
      dw1= dw1+ymvv*ypiq*vvis*rwms*(-101*xvc-
     #  1754./9.*xvq-28*xvv)
      dw1= dw1+ymvv*ypiq*rwms*(-2411./9.*xvc-758./
     #  3.*xvq-254./9.*xvv)
      dw1= dw1+ymvv*ypic*rlncu*vvis*rwms*(-100./9.
     #  *xvs-328./9.*xvc)
      dw1= dw1+ymvv*ypic*rlncu*rwms*(-376./9.*xvs-380./9.*xvc)
      dw1= dw1+ymvv*ypic*rlncd*vvis*rwms*(-22./9.*
     #  xvs-250./9.*xvc-212./9.*xvq-16./9.*xvv)
      dw1= dw1+ymvv*ypic*rlncd*rwms*(-142./9.*xvs-
     #  496./9.*xvc-80./3.*xvq-4./3.*xvv)
      dw1= dw1+ymvv*ypic*vvis*rwms*(199./9.*xvs+
     #  505./3.*xvc+242./3.*xvq+40./9.*xvv)
      dw1= dw1+ymvv*ypic*rwms*(1285./9.*xvs+2519./9.
     #  *xvc+236./3.*xvq+10./3.*xvv)
      dw1= dw1+ymvv*ypis*rlncu*vvis*rwms*(2./9.*xv+12*xvs)
      dw1= dw1+ymvv*ypis*rlncu*rwms*(86./9.*xv+184./9.*xvs)
      dw1= dw1+ymvv*ypis*rlncd*vvis*rwms*(-1./9.*
     #  xv+82./9.*xvs+223./9.*xvc+6*xvq)
      dw1= dw1+ymvv*ypis*rlncd*rwms*(32./9.*xv+274./9.
     #  *xvs+265./9.*xvc+10./3.*xvq)
      dw1= dw1+ymvv*ypis*vvis*rwms*(29./18.*xv-941.
     #  /18.*xvs-497./6.*xvc-15*xvq-2./3.*xvv)
      dw1= dw1+ymvv*ypis*rwms*(-593./18.*xv-1270./9.
     #  *xvs-169./2.*xvc-25./3.*xvq)
      dw1= dw1+ymvv*ypi*rlncu*vvis*rwms*(-8./9.*xv)
      dw1= dw1+ymvv*ypi*rlncu*rwms*(-4./9.-32./9.*xv)
      dw1= dw1+ymvv*ypi*rlncd*vvis*rwms*(-2./3.*
     #  xv-98./9.*xvs-68./9.*xvc)
      dw1= dw1+ymvv*ypi*rlncd*rwms*(-1./9.-62./9.*
     #  xv-46./3.*xvs-8./3.*xvc)
      dw1= dw1+ymvv*ypi*vvis*rwms*(-1./3.+25./9.
     #  *xv+319./9.*xvs+173./9.*xvc+8./3.*xvq)
      dw1= dw1+ymvv*ypi*rwms*(23./9.+277./9.*xv+383.
     #  /9.*xvs+29./3.*xvc-2./3.*xvq)
      dw1= dw1+ymvv*rlncd*vvis*rwms*omyp*(8./9.*xv)
      dw1= dw1+ymvv*rlncd*vvis*rwms*(7./9.*xv+38./9.*xvs)
      dw1= dw1+ymvv*rlncd*rwms*omyp*(1./9.)
      dw1= dw1+ymvv*rlncd*rwms*(1./9.+10./3.*xv+2./3.*xvs)
      dw1= dw1+ymvv*vvis*rwms*omyp*(-11./9.*xv-8./3.*xvs)
      dw1= dw1+ymvv*vvis*rwms*omyps*(-2./3.*xv)
      dw1= dw1+ymvv*vvis*rwms*(1./3.-61./18.*xv-77./
     #  9.*xvs-4*xvc)
      dw1= dw1+ymvv*rwms*omyp*(-11./18.-7./3.*xv+2./3.*xvs)
      dw1= dw1+ymvv*rwms*omyps*(-1./3.+1./3.*xv)
      dw1= dw1+ymvv*rwms*(-13./9.-73./9.*xv-20./3.*xvs+xvc)
      dw1= dw1+ymvvs*ypivi*rlncu*vvis*(-28./3.*xvvi)
      dw1= dw1+ymvvs*ypivi*rlncu*(-32./3.*xvvi)
      dw1= dw1+ymvvs*ypivi*rlncd*vvis*(-8./3.*xvvi)
      dw1= dw1+ymvvs*ypivi*rlncd*(-4*xvvi)
      dw1= dw1+ymvvs*ypivi*vvis*(250./9.*xvvi)
      dw1= dw1+ymvvs*ypivi*(116./3.*xvvi)
      dw1= dw1+ymvvs*ypiv*rlncu*vvis*(24.*xvv+68./9.*xvvi)
      dw1= dw1+ymvvs*ypiv*rlncu*(104./3.*xvv+20./3.*xvvi)
      dw1= dw1+ymvvs*ypiv*rlncd*vvis*(20./3.*xvv+52./9.*xvvi)
      dw1= dw1+ymvvs*ypiv*rlncd*(40./3.*xvv+8*xvvi)
      dw1= dw1+ymvvs*ypiv*vvis*(-200./3.*xvv-340./9.*xvvi)
      dw1= dw1+ymvvs*ypiv*(-386./3.*xvv-130./3.*xvvi)
      dw1= dw1+ymvvs*ypiq*rlncu*vvis*(-178./9.*xvq-188./9.*xvv)
      dw1= dw1+ymvvs*ypiq*rlncu*(-400./9.*xvq-64./3.*xvv)
      dw1= dw1+ymvvs*ypiq*rlncd*vvis*(-46./9.*
     #  xvq-142./9.*xvv-4*xvvi)
      dw1= dw1+ymvvs*ypiq*rlncd*(-157./9.*xvq-80./
     #  3.*xvv-14./3.*xvvi)
      dw1= dw1+ymvvs*ypiq*vvis*(146./3.*xvq+910./9.*xvv+14*xvvi)
      dw1= dw1+ymvvs*ypiq*(1495./9.*xvq+428./3.*
     #  xvv+127./9.*xvvi)
      dw1= dw1+ymvvs*ypic*rlncu*vvis*(44./9.*xvc+176./9.*xvq)
      dw1= dw1+ymvvs*ypic*rlncu*(244./9.*xvc+232./9.*xvq)
      dw1= dw1+ymvvs*ypic*rlncd*vvis*(8./9.*xvc+
     #  131./9.*xvq+112./9.*xvv+8./9.*xvvi)
      dw1= dw1+ymvvs*ypic*rlncd*(97./9.*xvc+313./9.*
     #  xvq+46./3.*xvv+2./3.*xvvi)
      dw1= dw1+ymvvs*ypic*vvis*(-68./9.*xvc-91
     #  *xvq-130./3.*xvv-20./9.*xvvi)
      dw1= dw1+ymvvs*ypic*(-908./9.*xvc-1612./9.*
     #  xvq-137./3.*xvv-5./3.*xvvi)
      dw1= dw1+ymvvs*ypis*rlncu*vvis*(2./9.*xvs-20./3.*xvc)
      dw1= dw1+ymvvs*ypis*rlncu*(-64./9.*xvs-128./9.*xvc)
      dw1= dw1+ymvvs*ypis*rlncd*vvis*(2./9.*xvs-
     #  44./9.*xvc-128./9.*xvq-10./3.*xvv)
      dw1= dw1+ymvvs*ypis*rlncd*(-25./9.*xvs-64./3.
     #  *xvc-58./3.*xvq-2*xvv)
      dw1= dw1+ymvvs*ypis*vvis*(-23./9.*xvs+
     #  523./18.*xvc+49*xvq+25./3.*xvv+1./3.*xvvi)
      dw1= dw1+ymvvs*ypis*(28*xvs+935./9.*xvc+511./9.*xvq+5*xvv)
      dw1= dw1+ymvvs*ypi*rlncu*vvis*(4./9.*xvs)
      dw1= dw1+ymvvs*ypi*rlncu*(4./9.*xv+28./9.*xvs)
      dw1= dw1+ymvvs*ypi*rlncd*vvis*(1./3.*xvs+
     #  64./9.*xvc+43./9.*xvq)
      dw1= dw1+ymvvs*ypi*rlncd*(1./9.*xv+49./9.*xvs
     #  +101./9.*xvc+2*xvq)
      dw1= dw1+ymvvs*ypi*vvis*(1./3.*xv-19./18.*
     #  xvs-215./9.*xvc-221./18.*xvq-5./3.*xvv)
      dw1= dw1+ymvvs*ypi*(-32./9.*xv-163./6.*xvs
     #  -203./6.*xvc-20./3.*xvq+1./3.*xvv)
      dw1= dw1+ymvvs*rlncd*vvis*omyp*(-7./9.*xvs)
      dw1= dw1+ymvvs*rlncd*vvis*(-5./9.*xvs-28./9.*xvc)
      dw1= dw1+ymvvs*rlncd*omyp*(-1./9.*xv)
      dw1= dw1+ymvvs*rlncd*(-1./9.*xv-8./3.*xvs-2./3.*xvc)
      dw1= dw1+ymvvs*vvis*omyp*(-1./24.+5./12.*xv-
     #  7./18.*xvs+10./3.*xvc)
      dw1= dw1+ymvvs*vvis*omyps*(-2./3.*xv+5./3.*xvs)
      dw1= dw1+ymvvs*vvis*omypc*(1./3.*xv)
      dw1= dw1+ymvvs*vvis*(-5./12.*xv+109./36.*xvs
     #  +49./9.*xvc+10./3.*xvq)
      dw1= dw1+ymvvs*omyp*(5./24.+1./36.*xv+25./6.*xvs-5./6.*xvc)
      dw1= dw1+ymvvs*omyps*(3./2.*xv-1./2.*xvs)
      dw1= dw1+ymvvs*omypc*(1./6.-1./6.*xv)
      dw1= dw1+ymvvs*(85./36.*xv+251./36.*xvs+13./2.*xvc-5./6.*xvq)
      dw1= dw1+ypivi*rlncu*vvis*rwms*rwmpgs*(-28./3.*xvq)
      dw1= dw1+ypivi*rlncu*rwms*rwmpgs*(-32./3.*xvq)
      dw1= dw1+ypivi*rlncd*vvis*rwms*rwmpgs*(-8./3.*xvq)
      dw1= dw1+ypivi*rlncd*rwms*rwmpgs*(-4*xvq)
      dw1= dw1+ypivi*vvis*rwms*rwmpgs*(250./9.*xvq)
      dw1= dw1+ypivi*rwms*rwmpgs*(116./3.*xvq)
      dw1= dw1+ypiv*rlncu*vvis*rwms*rwmpgs*(24.*xvc+68./9.*xvq)
      dw1= dw1+ypiv*rlncu*rwms*rwmpgs*(88./3.*xvc+20./3.*xvq)
      dw1= dw1+ypiv*rlncd*vvis*rwms*rwmpgs*(20./3.*xvc+52./9.*xvq)
      dw1= dw1+ypiv*rlncd*rwms*rwmpgs*(32./3.*xvc+8*xvq)
      dw1= dw1+ypiv*vvis*rwms*rwmpgs*(-200./3.*
     #  xvc-340./9.*xvq)
      dw1= dw1+ypiv*rwms*rwmpgs*(-310./3.*xvc-130./3.*xvq)
      dw1= dw1+ypiq*rlncu*vvis*rwms*rwmpgs*(-184./
     #  9.*xvs-176./9.*xvc)
      dw1= dw1+ypiq*rlncu*rwms*rwmpgs*(-278./9.*xvs-52./3.*xvc)
      dw1= dw1+ypiq*rlncd*vvis*rwms*rwmpgs*(-49./9.
     #  *xvs-136./9.*xvc-4*xvq)
      dw1= dw1+ypiq*rlncd*rwms*rwmpgs*(-98./9.*xvs
     #  -62./3.*xvc-14./3.*xvq)
      dw1= dw1+ypiq*vvis*rwms*rwmpgs*(157./3.*xvs
     #  +844./9.*xvc+14*xvq)
      dw1= dw1+ypiq*rwms*rwmpgs*(916./9.*xvs+110*xvc+127./9.*xvq)
      dw1= dw1+ypic*rlncu*vvis*rwms*rwmpgs*(56./9.*
     #  xv+152./9.*xvs)
      dw1= dw1+ypic*rlncu*rwms*rwmpgs*(44./3.*xv+148./9.*xvs)
      dw1= dw1+ypic*rlncd*vvis*rwms*rwmpgs*(14./9.*
     #  xv+119./9.*xvs+100./9.*xvc+8./9.*xvq)
      dw1= dw1+ypic*rlncd*rwms*rwmpgs*(5*xv+61./3.*
     #  xvs+34./3.*xvc+2./3.*xvq)
      dw1= dw1+ypic*vvis*rwms*rwmpgs*(-131./9.*xv
     #  -232./3.*xvs-112./3.*xvc-20./9.*xvq)
      dw1= dw1+ypic*rwms*rwmpgs*(-383./9.*xv-907./9.
     #  *xvs-33*xvc-5./3.*xvq)
      dw1= dw1+ypis*rlncu*vvis*rwms*rwmpgs*(-4./9.-16./3.*xv)
      dw1= dw1+ypis*rlncu*rwms*rwmpgs*(-22./9.-56./9.*xv)
      dw1= dw1+ypis*rlncd*vvis*rwms*rwmpgs*(-1./9.
     #  -38./9.*xv-95./9.*xvs-8./3.*xvc)
      dw1= dw1+ypis*rlncd*rwms*rwmpgs*(-7./9.-82./9.
     #  *xv-91./9.*xvs-4./3.*xvc)
      dw1= dw1+ypis*vvis*rwms*rwmpgs*(17./18.+209./
     #  9.*xv+203./6.*xvs+20./3.*xvc+1./3.*xvq)
      dw1= dw1+ypis*rwms*rwmpgs*(107./18.+341./9.*xv
     #  +499./18.*xvs+10./3.*xvc)
      dw1= dw1+ypi*rlncu*vvis*rwms*rwmpgs*(4./9.)
      dw1= dw1+ypi*rlncu*rwms*rwmpgs*(4./9.)
      dw1= dw1+ypi*rlncd*vvis*rwms*rwmpgs*(1./3.+
     #  34./9.*xv+25./9.*xvs)
      dw1= dw1+ypi*rlncd*rwms*rwmpgs*(13./9.+37./9.*xv+2./3.*xvs)
      dw1= dw1+ypi*vvis*rwms*rwmpgs*(-31./18.-
     #  104./9.*xv-125./18.*xvs-xvc)
      dw1= dw1+ypi*rwms*rwmpgs*(-83./18.-169./18.*
     #  xv-3*xvs+1./3.*xvc)
      dw1= dw1+rlncd*vvis*rwms*rwmpgs*omyp*(-1./9.)
      dw1= dw1+rlncd*vvis*rwms*rwmpgs*(-2./9.-10./9.*xv)
      dw1= dw1+rlncd*rwms*rwmpgs*(-2./3.)
      dw1= dw1+vvis*rwms*rwmpgs*omyp*(5./18.+1./3.*xv)
      dw1= dw1+vvis*rwms*rwmpgs*(7./9.+22./9.*xv+xvs)
      dw1= dw1+rwms*rwmpgs*omyp*(1./6.-1./6.*xv)
      dw1= dw1+rwms*rwmpgs*(11./9.+5./6.*xv-1./6.*xvs)
      dw1= dw1+denu*ymvvs*ypiq*vviq*rwms*rwmpgs*(-4./9.*xvq)
      dw1= dw1+denu*ymvvs*ypic*vviq*rwms*rwmpgs*(4./3.*xvc)
      dw1= dw1+denu*ymvvs*ypis*vviq*rwms*rwmpgs*(-4./3.*xvs)
      dw1= dw1+denu*ymvvs*ypis*vvis*rwms*rwmpgs*(-2./9.*xvs)
      dw1= dw1+denu*ymvvs*ypi*vviq*rwms*rwmpgs*(4./9.*xv)
      dw1= dw1+denu*ymvvs*ypi*vvis*rwms*rwmpgs*(4./9.*xv)
      dw1= dw1+denu*ymvvs*vvis*rwms*rwmpgs*(-2./9.)
      dw1= dw1+denu*ymvvc*ypiq*vviq*rwms*(8./9.*xvv)
      dw1= dw1+denu*ymvvc*ypic*vviq*rwms*(-8./3.*xvq)
      dw1= dw1+denu*ymvvc*ypis*vviq*rwms*(8./3.*xvc)
      dw1= dw1+denu*ymvvc*ypis*vvis*rwms*(4./9.*xvc)
      dw1= dw1+denu*ymvvc*ypi*vviq*rwms*(-8./9.*xvs)
      dw1= dw1+denu*ymvvc*ypi*vvis*rwms*(-8./9.*xvs)
      dw1= dw1+denu*ymvvc*vvis*rwms*(4./9.*xv)
      dw1= dw1+denu*ymvvq*ypiq*vviq*(-4./9.*xvvi)
      dw1= dw1+denu*ymvvq*ypic*vviq*(4./3.*xvv)
      dw1= dw1+denu*ymvvq*ypis*vviq*(-4./3.*xvq)
      dw1= dw1+denu*ymvvq*ypis*vvis*(-2./9.*xvq)
      dw1= dw1+denu*ymvvq*ypi*vviq*(4./9.*xvc)
      dw1= dw1+denu*ymvvq*ypi*vvis*(4./9.*xvc)
      dw1= dw1+denu*ymvvq*vvis*(-2./9.*xvs)
      dw1= dw1+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-1./9.*xvq)
      dw1= dw1+dend*ymvvs*ypic*vviq*rwms*rwmpgs
     #     *(1./3.*xvc+2./9.*xvq)
      dw1= dw1+dend*ymvvs*ypis*vviq*rwms*rwmpgs
     #     *(-1./3.*xvs-2./3.*xvc-1./9.*xvq)
      dw1= dw1+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./18.*xvs)
      dw1= dw1+dend*ymvvs*ypi*vviq*rwms*rwmpgs
     #     *(1./9.*xv+2./3.*xvs+1./3.*xvc)
      dw1= dw1+dend*ymvvs*ypi*vvis*rwms*rwmpgs
     #     *(1./9.*xv+1./9.*xvs)
      dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-1./9.*xv)
      dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*(-1./9.
     #  *xv-1./3.*xvs)
      dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-1./9.*xv)
      dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./18.)
      dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*(-1./9.*xv-1./18.*xvs)
      dw1= dw1+dend*ymvvc*ypiq*vviq*rwms*(2./9.*xvv)
      dw1= dw1+dend*ymvvc*ypic*vviq*rwms*(-2./3.*xvq-4./9.*xvv)
      dw1= dw1+dend*ymvvc*ypis*vviq*rwms*(2./3.*
     #  xvc+4./3.*xvq+2./9.*xvv)
      dw1= dw1+dend*ymvvc*ypis*vvis*rwms*(1./9.*xvc)
      dw1= dw1+dend*ymvvc*ypi*vviq*rwms*(-2./
     #  9.*xvs-4./3.*xvc-2./3.*xvq)
      dw1= dw1+dend*ymvvc*ypi*vvis*rwms*(-2./9.*xvs-2./9.*xvc)
      dw1= dw1+dend*ymvvc*vviq*rwms*omyp*(2./9.*xvs)
      dw1= dw1+dend*ymvvc*vviq*rwms*(2./9.*xvs+2./3.*xvc)
      dw1= dw1+dend*ymvvc*vvis*rwms*omyp*(2./9.*xvs)
      dw1= dw1+dend*ymvvc*vvis*rwms*omyps*(1./9.*xv)
      dw1= dw1+dend*ymvvc*vvis*rwms*(2./9.*xvs+1./9.*xvc)
      dw1= dw1+dend*ymvvq*ypiq*vviq*(-1./9.*xvvi)
      dw1= dw1+dend*ymvvq*ypic*vviq*(1./3.*xvv+2./9.*xvvi)
      dw1= dw1+dend*ymvvq*ypis*vviq*(-1./3.*
     #  xvq-2./3.*xvv-1./9.*xvvi)
      dw1= dw1+dend*ymvvq*ypis*vvis*(-1./18.*xvq)
      dw1= dw1+dend*ymvvq*ypi*vviq*(1./9.*xvc
     #  +2./3.*xvq+1./3.*xvv)
      dw1= dw1+dend*ymvvq*ypi*vvis*(1./9.*xvc+1./9.*xvq)
      dw1= dw1+dend*ymvvq*vviq*omyp*(-1./9.*xvc)
      dw1= dw1+dend*ymvvq*vviq*(-1./9.*xvc-1./3.*xvq)
      dw1= dw1+dend*ymvvq*vvis*omyp*(-1./9.*xvc)
      dw1= dw1+dend*ymvvq*vvis*omyps*(-1./18.*xvs)
      dw1= dw1+dend*ymvvq*vvis*(-1./9.*xvc-1./18.*xvq)
*
      dw2= ymvv*ypivi*rlncu*vvis*rwms*(56*xvv)
      dw2= dw2+ymvv*ypivi*rlncu*rwms*(64./3.*xvv)
      dw2= dw2+ymvv*ypivi*rlncd*vvis*rwms*(16.*xvv)
      dw2= dw2+ymvv*ypivi*rlncd*rwms*(8*xvv)
      dw2= dw2+ymvv*ypivi*vvis*rwms*(-500./3.*xvv)
      dw2= dw2+ymvv*ypivi*rwms*(-232./3.*xvv)
      dw2= dw2+ymvv*ypiv*rlncu*vvis*rwms*(-144*xvq-136./3.*xvv)
      dw2= dw2+ymvv*ypiv*rlncu*rwms*(-64*xvq-40./3.*xvv)
      dw2= dw2+ymvv*ypiv*rlncd*vvis*rwms*(-40*xvq-104./3.*xvv)
      dw2= dw2+ymvv*ypiv*rlncd*rwms*(-24.*xvq-16.*xvv)
      dw2= dw2+ymvv*ypiv*vvis*rwms*(400.*xvq+680./3.*xvv)
      dw2= dw2+ymvv*ypiv*rwms*(232.*xvq+260./3.*xvv)
      dw2= dw2+ymvv*ypiq*rlncu*vvis*rwms*(362./3.*
     #  xvc+364./3.*xvq)
      dw2= dw2+ymvv*ypiq*rlncu*rwms*(226./3.*xvc+116./3.*xvq)
      dw2= dw2+ymvv*ypiq*rlncd*vvis*rwms*(95./3.*
     #  xvc+278./3.*xvq+24.*xvv)
      dw2= dw2+ymvv*ypiq*rlncd*rwms*(85./3.*xvc+142./
     #  3.*xvq+28./3.*xvv)
      dw2= dw2+ymvv*ypiq*vvis*rwms*(-303*xvc-
     #  1754./3.*xvq-84*xvv)
      dw2= dw2+ymvv*ypiq*rwms*(-2411./9.*xvc-758./
     #  3.*xvq-254./9.*xvv)
      dw2= dw2+ymvv*ypic*rlncu*vvis*rwms*(-100./3.
     #  *xvs-328./3.*xvc)
      dw2= dw2+ymvv*ypic*rlncu*rwms*(-376./9.*xvs-380./9.*xvc)
      dw2= dw2+ymvv*ypic*rlncd*vvis*rwms*(-22./3.*
     #  xvs-250./3.*xvc-212./3.*xvq-16./3.*xvv)
      dw2= dw2+ymvv*ypic*rlncd*rwms*(-142./9.*xvs-
     #  496./9.*xvc-80./3.*xvq-4./3.*xvv)
      dw2= dw2+ymvv*ypic*vvis*rwms*(199./3.*xvs+
     #  505*xvc+242.*xvq+40./3.*xvv)
      dw2= dw2+ymvv*ypic*rwms*(1285./9.*xvs+2519./9.
     #  *xvc+236./3.*xvq+10./3.*xvv)
      dw2= dw2+ymvv*ypis*rlncu*vvis*rwms*(2./3.*xv+36*xvs)
      dw2= dw2+ymvv*ypis*rlncu*rwms*(86./9.*xv+184./9.*xvs)
      dw2= dw2+ymvv*ypis*rlncd*vvis*rwms*(-1./3.*
     #  xv+82./3.*xvs+223./3.*xvc+18*xvq)
      dw2= dw2+ymvv*ypis*rlncd*rwms*(32./9.*xv+274./9.
     #  *xvs+265./9.*xvc+10./3.*xvq)
      dw2= dw2+ymvv*ypis*vvis*rwms*(29./6.*xv-941./
     #  6.*xvs-497./2.*xvc-45*xvq-2*xvv)
      dw2= dw2+ymvv*ypis*rwms*(-593./18.*xv-1270./9.
     #  *xvs-169./2.*xvc-25./3.*xvq)
      dw2= dw2+ymvv*ypi*rlncu*vvis*rwms*(-8./3.*xv)
      dw2= dw2+ymvv*ypi*rlncu*rwms*(-4./9.-32./9.*xv)
      dw2= dw2+ymvv*ypi*rlncd*vvis*rwms*(-2*xv-
     #  98./3.*xvs-68./3.*xvc)
      dw2= dw2+ymvv*ypi*rlncd*rwms*(-1./9.-62./9.*
     #  xv-46./3.*xvs-8./3.*xvc)
      dw2= dw2+ymvv*ypi*vvis*rwms*(-1+25./3.*xv
     #  +319./3.*xvs+173./3.*xvc+8*xvq)
      dw2= dw2+ymvv*ypi*rwms*(23./9.+277./9.*xv+383.
     #  /9.*xvs+29./3.*xvc-2./3.*xvq)
      dw2= dw2+ymvv*rlncd*vvis*rwms*omyp*(8./3.*xv)
      dw2= dw2+ymvv*rlncd*vvis*rwms*(7./3.*xv+38./3.*xvs)
      dw2= dw2+ymvv*rlncd*rwms*omyp*(1./9.)
      dw2= dw2+ymvv*rlncd*rwms*(1./9.+10./3.*xv+2./3.*xvs)
      dw2= dw2+ymvv*vvis*rwms*omyp*(-11./3.*xv-8*xvs)
      dw2= dw2+ymvv*vvis*rwms*omyps*(-2*xv)
      dw2= dw2+ymvv*vvis*rwms*(1.-61./6.*xv-77./3.*xvs-12*xvc)
      dw2= dw2+ymvv*rwms*omyp*(-11./18.-7./3.*xv+2./3.*xvs)
      dw2= dw2+ymvv*rwms*omyps*(-1./3.+1./3.*xv)
      dw2= dw2+ymvv*rwms*(-13./9.-73./9.*xv-20./3.*xvs+xvc)
      dw2= dw2+ymvvs*ypivi*rlncu*vvis*(-28*xvvi)
      dw2= dw2+ymvvs*ypivi*rlncu*(-32./3.*xvvi)
      dw2= dw2+ymvvs*ypivi*rlncd*vvis*(-8*xvvi)
      dw2= dw2+ymvvs*ypivi*rlncd*(-4*xvvi)
      dw2= dw2+ymvvs*ypivi*vvis*(250./3.*xvvi)
      dw2= dw2+ymvvs*ypivi*(116./3.*xvvi)
      dw2= dw2+ymvvs*ypiv*rlncu*vvis*(72*xvv+68./3.*xvvi)
      dw2= dw2+ymvvs*ypiv*rlncu*(104./3.*xvv+20./3.*xvvi)
      dw2= dw2+ymvvs*ypiv*rlncd*vvis*(20*xvv+52./3.*xvvi)
      dw2= dw2+ymvvs*ypiv*rlncd*(40./3.*xvv+8*xvvi)
      dw2= dw2+ymvvs*ypiv*vvis*(-200*xvv-340./3.*xvvi)
      dw2= dw2+ymvvs*ypiv*(-386./3.*xvv-130./3.*xvvi)
      dw2= dw2+ymvvs*ypiq*rlncu*vvis*(-178./3.*xvq-188./3.*xvv)
      dw2= dw2+ymvvs*ypiq*rlncu*(-400./9.*xvq-64./3.*xvv)
      dw2= dw2+ymvvs*ypiq*rlncd*vvis*(-46./3.*
     #  xvq-142./3.*xvv-12*xvvi)
      dw2= dw2+ymvvs*ypiq*rlncd*(-157./9.*xvq-80./
     #  3.*xvv-14./3.*xvvi)
      dw2= dw2+ymvvs*ypiq*vvis*(146*xvq+910./3.*xvv+42.*xvvi)
      dw2= dw2+ymvvs*ypiq*(1495./9.*xvq+428./3.*xvv+127./9.*xvvi)
      dw2= dw2+ymvvs*ypic*rlncu*vvis*(44./3.*xvc+176./3.*xvq)
      dw2= dw2+ymvvs*ypic*rlncu*(244./9.*xvc+232./9.*xvq)
      dw2= dw2+ymvvs*ypic*rlncd*vvis*(8./3.*xvc+
     #  131./3.*xvq+112./3.*xvv+8./3.*xvvi)
      dw2= dw2+ymvvs*ypic*rlncd*(97./9.*xvc+313./9.*
     #  xvq+46./3.*xvv+2./3.*xvvi)
      dw2= dw2+ymvvs*ypic*vvis*(-68./3.*xvc-
     #  273*xvq-130*xvv-20./3.*xvvi)
      dw2= dw2+ymvvs*ypic*(-908./9.*xvc-1612./9.*
     #  xvq-137./3.*xvv-5./3.*xvvi)
      dw2= dw2+ymvvs*ypis*rlncu*vvis*(2./3.*xvs-20*xvc)
      dw2= dw2+ymvvs*ypis*rlncu*(-64./9.*xvs-128./9.*xvc)
      dw2= dw2+ymvvs*ypis*rlncd*vvis*(2./3.*xvs-
     #  44./3.*xvc-128./3.*xvq-10*xvv)
      dw2= dw2+ymvvs*ypis*rlncd*(-25./9.*xvs-64./3.
     #  *xvc-58./3.*xvq-2*xvv)
      dw2= dw2+ymvvs*ypis*vvis*(-23./3.*xvs+
     #  523./6.*xvc+147*xvq+25*xvv+xvvi)
      dw2= dw2+ymvvs*ypis*(28*xvs+935./9.*xvc+511./9.*xvq+5*xvv)
      dw2= dw2+ymvvs*ypi*rlncu*vvis*(4./3.*xvs)
      dw2= dw2+ymvvs*ypi*rlncu*(4./9.*xv+28./9.*xvs)
      dw2= dw2+ymvvs*ypi*rlncd*vvis*(xvs+64./3.*xvc+43./3.*xvq)
      dw2= dw2+ymvvs*ypi*rlncd*(1./9.*xv+49./9.*xvs
     #  +101./9.*xvc+2*xvq)
      dw2= dw2+ymvvs*ypi*vvis*(xv-19./6.*xvs-
     #  215./3.*xvc-221./6.*xvq-5*xvv)
      dw2= dw2+ymvvs*ypi*(-32./9.*xv-163./6.*xvs
     #  -203./6.*xvc-20./3.*xvq+1./3.*xvv)
      dw2= dw2+ymvvs*rlncd*vvis*omyp*(-7./3.*xvs)
      dw2= dw2+ymvvs*rlncd*vvis*(-5./3.*xvs-28./3.*xvc)
      dw2= dw2+ymvvs*rlncd*omyp*(-1./9.*xv)
      dw2= dw2+ymvvs*rlncd*(-1./9.*xv-8./3.*xvs-2./3.*xvc)
      dw2= dw2+ymvvs*vvis*omyp*(-1./8.+5./4.*xv-7./6.*xvs+10*xvc)
      dw2= dw2+ymvvs*vvis*omyps*(-2*xv+5*xvs)
      dw2= dw2+ymvvs*vvis*omypc*(xv)
      dw2= dw2+ymvvs*vvis*(-5./4.*xv+109./12.*xvs
     #  +49./3.*xvc+10*xvq)
      dw2= dw2+ymvvs*omyp*(5./24.+1./36.*xv+25./6.*xvs-5./6.*xvc)
      dw2= dw2+ymvvs*omyps*(3./2.*xv-1./2.*xvs)
      dw2= dw2+ymvvs*omypc*(1./6.-1./6.*xv)
      dw2= dw2+ymvvs*(85./36.*xv+251./36.*xvs+13./2.*xvc-5./6.*xvq)
      dw2= dw2+ypivi*rlncu*vvis*rwms*rwmpgs*(-28*xvq)
      dw2= dw2+ypivi*rlncu*rwms*rwmpgs*(-32./3.*xvq)
      dw2= dw2+ypivi*rlncd*vvis*rwms*rwmpgs*(-8*xvq)
      dw2= dw2+ypivi*rlncd*rwms*rwmpgs*(-4.*xvq)
      dw2= dw2+ypivi*vvis*rwms*rwmpgs*(250./3.*xvq)
      dw2= dw2+ypivi*rwms*rwmpgs*(116./3.*xvq)
      dw2= dw2+ypiv*rlncu*vvis*rwms*rwmpgs*(72*xvc+68./3.*xvq)
      dw2= dw2+ypiv*rlncu*rwms*rwmpgs*(88./3.*xvc+20./3.*xvq)
      dw2= dw2+ypiv*rlncd*vvis*rwms*rwmpgs*(20*xvc+52./3.*xvq)
      dw2= dw2+ypiv*rlncd*rwms*rwmpgs*(32./3.*xvc+8*xvq)
      dw2= dw2+ypiv*vvis*rwms*rwmpgs*(-200*xvc-340./3.*xvq)
      dw2= dw2+ypiv*rwms*rwmpgs*(-310./3.*xvc-130./3.*xvq)
      dw2= dw2+ypiq*rlncu*vvis*rwms*rwmpgs*(-184./
     #  3.*xvs-176./3.*xvc)
      dw2= dw2+ypiq*rlncu*rwms*rwmpgs*(-278./9.*xvs-52./3.*xvc)
      dw2= dw2+ypiq*rlncd*vvis*rwms*rwmpgs*(-49./3.
     #  *xvs-136./3.*xvc-12*xvq)
      dw2= dw2+ypiq*rlncd*rwms*rwmpgs*(-98./9.*xvs
     #  -62./3.*xvc-14./3.*xvq)
      dw2= dw2+ypiq*vvis*rwms*rwmpgs*(157*xvs+844./3.*xvc+42.*xvq)
      dw2= dw2+ypiq*rwms*rwmpgs*(916./9.*xvs+110*xvc+127./9.*xvq)
      dw2= dw2+ypic*rlncu*vvis*rwms*rwmpgs*(56./3.*xv+152./3.*xvs)
      dw2= dw2+ypic*rlncu*rwms*rwmpgs*(44./3.*xv+148./9.*xvs)
      dw2= dw2+ypic*rlncd*vvis*rwms*rwmpgs*(14./3.*
     #  xv+119./3.*xvs+100./3.*xvc+8./3.*xvq)
      dw2= dw2+ypic*rlncd*rwms*rwmpgs*(5*xv+61./3.*
     #  xvs+34./3.*xvc+2./3.*xvq)
      dw2= dw2+ypic*vvis*rwms*rwmpgs*(-131./3.*xv
     #  -232.*xvs-112*xvc-20./3.*xvq)
      dw2= dw2+ypic*rwms*rwmpgs*(-383./9.*xv-907./9.
     #  *xvs-33*xvc-5./3.*xvq)
      dw2= dw2+ypis*rlncu*vvis*rwms*rwmpgs*(-4./3.-16.*xv)
      dw2= dw2+ypis*rlncu*rwms*rwmpgs*(-22./9.-56./9.*xv)
      dw2= dw2+ypis*rlncd*vvis*rwms*rwmpgs*(-1./3.
     #  -38./3.*xv-95./3.*xvs-8*xvc)
      dw2= dw2+ypis*rlncd*rwms*rwmpgs*(-7./9.-82./9.
     #  *xv-91./9.*xvs-4./3.*xvc)
      dw2= dw2+ypis*vvis*rwms*rwmpgs*(17./6.+209./
     #  3.*xv+203./2.*xvs+20*xvc+xvq)
      dw2= dw2+ypis*rwms*rwmpgs*(107./18.+341./9.*xv
     #  +499./18.*xvs+10./3.*xvc)
      dw2= dw2+ypi*rlncu*vvis*rwms*rwmpgs*(4./3.)
      dw2= dw2+ypi*rlncu*rwms*rwmpgs*(4./9.)
      dw2= dw2+ypi*rlncd*vvis*rwms*rwmpgs*(1+34./3.*xv+25./3.*xvs)
      dw2= dw2+ypi*rlncd*rwms*rwmpgs*(13./9.+37./9.*xv+2./3.*xvs)
      dw2= dw2+ypi*vvis*rwms*rwmpgs*(-31./6.-
     #  104./3.*xv-125./6.*xvs-3*xvc)
      dw2= dw2+ypi*rwms*rwmpgs*(-83./18.-169./18.*xv-3*xvs+1./3.*xvc)
      dw2= dw2+rlncd*vvis*rwms*rwmpgs*omyp*(-1./3.)
      dw2= dw2+rlncd*vvis*rwms*rwmpgs*(-2./3.-10./3.*xv)
      dw2= dw2+rlncd*rwms*rwmpgs*(-2./3.)
      dw2= dw2+vvis*rwms*rwmpgs*omyp*(5./6.+xv)
      dw2= dw2+vvis*rwms*rwmpgs*(7./3.+22./3.*xv+3*xvs)
      dw2= dw2+rwms*rwmpgs*omyp*(1./6.-1./6.*xv)
      dw2= dw2+rwms*rwmpgs*(11./9.+5./6.*xv-1./6.*xvs)
      dw2= dw2+denu*ymvvs*ypiq*vviq*rwms*rwmpgs*(-4./3.*xvq)
      dw2= dw2+denu*ymvvs*ypic*vviq*rwms*rwmpgs*(4.*xvc)
      dw2= dw2+denu*ymvvs*ypis*vviq*rwms*rwmpgs*(-4.*xvs)
      dw2= dw2+denu*ymvvs*ypis*vvis*rwms*rwmpgs*(-2./9.*xvs)
      dw2= dw2+denu*ymvvs*ypi*vviq*rwms*rwmpgs*(4./3.*xv)
      dw2= dw2+denu*ymvvs*ypi*vvis*rwms*rwmpgs*(4./9.*xv)
      dw2= dw2+denu*ymvvs*vvis*rwms*rwmpgs*(-2./9.)
      dw2= dw2+denu*ymvvc*ypiq*vviq*rwms*(8./3.*xvv)
      dw2= dw2+denu*ymvvc*ypic*vviq*rwms*(-8*xvq)
      dw2= dw2+denu*ymvvc*ypis*vviq*rwms*(8*xvc)
      dw2= dw2+denu*ymvvc*ypis*vvis*rwms*(4./9.*xvc)
      dw2= dw2+denu*ymvvc*ypi*vviq*rwms*(-8./3.*xvs)
      dw2= dw2+denu*ymvvc*ypi*vvis*rwms*(-8./9.*xvs)
      dw2= dw2+denu*ymvvc*vvis*rwms*(4./9.*xv)
      dw2= dw2+denu*ymvvq*ypiq*vviq*(-4./3.*xvvi)
      dw2= dw2+denu*ymvvq*ypic*vviq*(4.*xvv)
      dw2= dw2+denu*ymvvq*ypis*vviq*(-4.*xvq)
      dw2= dw2+denu*ymvvq*ypis*vvis*(-2./9.*xvq)
      dw2= dw2+denu*ymvvq*ypi*vviq*(4./3.*xvc)
      dw2= dw2+denu*ymvvq*ypi*vvis*(4./9.*xvc)
      dw2= dw2+denu*ymvvq*vvis*(-2./9.*xvs)
      dw2= dw2+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-1./3.*xvq)
      dw2= dw2+dend*ymvvs*ypic*vviq*rwms*rwmpgs*(xvc+2./3.*xvq)
      dw2= dw2+dend*ymvvs*ypis*vviq*rwms*rwmpgs*(-xvs-2*xvc-1./3.*xvq)
      dw2= dw2+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./18.*xvs)
      dw2= dw2+dend*ymvvs*ypi*vviq*rwms*rwmpgs*(1./3.*xv+2*xvs+xvc)
      dw2= dw2+dend*ymvvs*ypi*vvis*rwms*rwmpgs*(1./9.*xv+1./9.*xvs)
      dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-1./3.*xv)
      dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*(-1./3.*xv-xvs)
      dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-1./9.*xv)
      dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./18.)
      dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*(-1./9.*xv-1./18.*xvs)
      dw2= dw2+dend*ymvvc*ypiq*vviq*rwms*(2./3.*xvv)
      dw2= dw2+dend*ymvvc*ypic*vviq*rwms*(-2*xvq-4./3.*xvv)
      dw2= dw2+dend*ymvvc*ypis*vviq*rwms*(2*xvc+4.*xvq+2./3.*xvv)
      dw2= dw2+dend*ymvvc*ypis*vvis*rwms*(1./9.*xvc)
      dw2= dw2+dend*ymvvc*ypi*vviq*rwms*(-2./3.*xvs-4.*xvc-2*xvq)
      dw2= dw2+dend*ymvvc*ypi*vvis*rwms*(-2./9.*xvs-2./9.*xvc)
      dw2= dw2+dend*ymvvc*vviq*rwms*omyp*(2./3.*xvs)
      dw2= dw2+dend*ymvvc*vviq*rwms*(2./3.*xvs+2*xvc)
      dw2= dw2+dend*ymvvc*vvis*rwms*omyp*(2./9.*xvs)
      dw2= dw2+dend*ymvvc*vvis*rwms*omyps*(1./9.*xv)
      dw2= dw2+dend*ymvvc*vvis*rwms*(2./9.*xvs+1./9.*xvc)
      dw2= dw2+dend*ymvvq*ypiq*vviq*(-1./3.*xvvi)
      dw2= dw2+dend*ymvvq*ypic*vviq*(xvv+2./3.*xvvi)
      dw2= dw2+dend*ymvvq*ypis*vviq*(-xvq-2*xvv-1./3.*xvvi)
      dw2= dw2+dend*ymvvq*ypis*vvis*(-1./18.*xvq)
      dw2= dw2+dend*ymvvq*ypi*vviq*(1./3.*xvc+2*xvq+xvv)
      dw2= dw2+dend*ymvvq*ypi*vvis*(1./9.*xvc+1./9.*xvq)
      dw2= dw2+dend*ymvvq*vviq*omyp*(-1./3.*xvc)
      dw2= dw2+dend*ymvvq*vviq*(-1./3.*xvc-xvq)
      dw2= dw2+dend*ymvvq*vvis*omyp*(-1./9.*xvc)
      dw2= dw2+dend*ymvvq*vvis*omyps*(-1./18.*xvs)
      dw2= dw2+dend*ymvvq*vvis*(-1./9.*xvc-1./18.*xvq)
*
      if(ii.eq.3) then
       wcom= wcom*(2.d0-ym)**2/yms-dw1+2.d0/ym*(1.d0-1.d0/ym)*dw2
      else if(ii.eq.4) then
       wcom= 0.5d0*em4*omxm*omxm*(dw1-dw2)     
      else if(ii.eq.5) then
       wcom= em2*((2.d0+0.5d0*omxm)*dw1+omxm*(2.d0/ym-
     #       3.d0/2.d0)*dw2)
      endif
*
      else if((ipro.eq.1.or.ipro.eq.3).and.ii.le.2) then
*
      wcomd1= ypiq*rwms*(-2*xvq)
      wcomd1= wcomd1+ypic*rwms*(4*xvc+4*xvq)
      wcomd1= wcomd1+ypis*rwms*(-3*xvs-8*xvc-2*xvq)
      wcomd1= wcomd1+ypi*rwms*(xv+6*xvs+4*xvc)
      wcomd1= wcomd1+rwms*omyp*(-xv)
      wcomd1= wcomd1+rwms*(-xv-3*xvs)
      wcomd1= wcomd1+ymvv*ypiq*(xvv)
      wcomd1= wcomd1+ymvv*ypic*(-2*xvq-2*xvv)
      wcomd1= wcomd1+ymvv*ypis*(3./2.*xvc+4*xvq+xvv)
      wcomd1= wcomd1+ymvv*ypi*(-1./2.*xvs-3*xvc-2*xvq)
      wcomd1= wcomd1+ymvv*omyp*(1./2.*xvs)
      wcomd1= wcomd1+ymvv*(1./2.*xvs+3./2.*xvc)
      wcomd1= ymvv*wcomd1
*
      wcomd2= ypiq*(xvc)
      wcomd2= wcomd2+ypic*(-2*xvs-2*xvc)
      wcomd2= wcomd2+ypis*(3./2.*xv+4*xvs+xvc)
      wcomd2= wcomd2+ypi*(-1./2.-3*xv-2*xvs)
      wcomd2= wcomd2+omyp*(1./2.)
      wcomd2= wcomd2+(1./2.+3./2.*xv)
      wcomd2= rwms*rwmpgs*wcomd2
*
      wcom01= ypiq*rwms*(22./3.*xvq)
      wcom01= wcom01+ypic*rwms*(-44./3.*xvc-14*xvq)
      wcom01= wcom01+ypis*rwms*(10*xvs+28*xvc+7*xvq)
      wcom01= wcom01+ypi*rwms*(-8./3.*xv-20*xvs-14*xvc-2./3.*xvq)
      wcom01= wcom01+rwms*omyp*(17./6.*xv+xvs)
      wcom01= wcom01+rwms*omyps*(1./3.*xv)
      wcom01= wcom01+rwms*(17./6.*xv+19./2.*xvs+4./3.*xvc)
      wcom01= wcom01+ymvv*ypiq*(-11./3.*xvv)
      wcom01= wcom01+ymvv*ypic*(23./3.*xvq+7*xvv)
      wcom01= wcom01+ymvv*ypis*(-11./2.*xvc-15*xvq-7./2.*xvv)
      wcom01= wcom01+ymvv*ypi*(11./6.*xvs+23./2.*xvc+8*xvq+1./3.*xvv)
      wcom01= wcom01+ymvv*omyp*(-1./6.*xv-19./12.*xvs-7./6.*xvc)
      wcom01= wcom01+ymvv*omyps*(-2./3.*xvs)
      wcom01= wcom01+ymvv*omypc*(-1./6.*xv)
      wcom01= wcom01+ymvv*(-23./12.*xvs-23./4.*xvc-xvq)
      wcom01= ymvv*wcom01
*
      wcom02= ypiq*(-11./3.*xvc)
      wcom02= wcom02+ypic*(7*xvs+7*xvc)
      wcom02= wcom02+ypis*(-9./2.*xv-13*xvs-7./2.*xvc)
      wcom02= wcom02+ypi*(5./6.+17./2.*xv+6*xvs+1./3.*xvc)
      wcom02= wcom02+omyp*(-11./12.-1./6.*xv)
      wcom02= wcom02+(-11./12.-15./4.*xv-1./3.*xvs)
      wcom02= rwms*rwmpgs*wcom02
*
      wcom0= wcom01+wcom02
*
      wcomdm= ypis*vvi*rwms*rwmpgs*xvs
      wcomdm= wcomdm+ypi*vvi*rwms*rwmpgs*(-2*xv-2*xvs)
      wcomdm= wcomdm+vvi*rwms*rwmpgs*omyp*(2*xv)
      wcomdm= wcomdm+vvi*rwms*rwmpgs*omyps
      wcomdm= wcomdm+vvi*rwms*rwmpgs*(2*xv+xvs)
      wcomdm= wcomdm+ymvv*ypis*vvi*rwms*(-2*xvc)
      wcomdm= wcomdm+ymvv*ypi*vvi*rwms*(4*xvs+4*xvc)
      wcomdm= wcomdm+ymvv*vvi*rwms*omyp*(-4*xvs)
      wcomdm= wcomdm+ymvv*vvi*rwms*omyps*(-2*xv)
      wcomdm= wcomdm+ymvv*vvi*rwms*(-4*xvs-2*xvc)
      wcomdm= wcomdm+ymvvs*ypis*vvi*xvq
      wcomdm= wcomdm+ymvvs*ypi*vvi*(-2*xvc-2*xvq)
      wcomdm= wcomdm+ymvvs*vvi*omyp*(2*xvc)
      wcomdm= wcomdm+ymvvs*vvi*omyps*xvs
      wcomdm= wcomdm+ymvvs*vvi*(2*xvc+xvq)
      dendm= dend*dms/vv
      wcomdm= ymvv*dendm*wcomdm
*
      wcom= wcomd*rlncd+wcom0+wcomdm
*
      else if((ipro.eq.1.or.ipro.eq.3).and.ii.ge.3) then
*
      wcomd1= ypiq*rwms*(-2*xvq)
      wcomd1= wcomd1+ypic*rwms*(4*xvc+4*xvq)
      wcomd1= wcomd1+ypis*rwms*(-3*xvs-8*xvc-2*xvq)
      wcomd1= wcomd1+ypi*rwms*(xv+6*xvs+4*xvc)
      wcomd1= wcomd1+rwms*omyp*(-xv)
      wcomd1= wcomd1+rwms*(-xv-3*xvs)
      wcomd1= wcomd1+ymvv*ypiq*(xvv)
      wcomd1= wcomd1+ymvv*ypic*(-2*xvq-2*xvv)
      wcomd1= wcomd1+ymvv*ypis*(3./2.*xvc+4*xvq+xvv)
      wcomd1= wcomd1+ymvv*ypi*(-1./2.*xvs-3*xvc-2*xvq)
      wcomd1= wcomd1+ymvv*omyp*(1./2.*xvs)
      wcomd1= wcomd1+ymvv*(1./2.*xvs+3./2.*xvc)
      wcomd1= ymvv*wcomd1
*
      wcomd2= ypiq*(xvc)
      wcomd2= wcomd2+ypic*(-2*xvs-2*xvc)
      wcomd2= wcomd2+ypis*(3./2.*xv+4*xvs+xvc)
      wcomd2= wcomd2+ypi*(-1./2.-3*xv-2*xvs)
      wcomd2= wcomd2+omyp*(1./2.)
      wcomd2= wcomd2+(1./2.+3./2.*xv)
      wcomd2= rwms*rwmpgs*wcomd2
*
      wcom01= ypiq*rwms*(22./3.*xvq)
      wcom01= wcom01+ypic*rwms*(-44./3.*xvc-14*xvq)
      wcom01= wcom01+ypis*rwms*(10*xvs+28*xvc+7*xvq)
      wcom01= wcom01+ypi*rwms*(-8./3.*xv-20*xvs-14*xvc-2./3.*xvq)
      wcom01= wcom01+rwms*omyp*(17./6.*xv+xvs)
      wcom01= wcom01+rwms*omyps*(1./3.*xv)
      wcom01= wcom01+rwms*(17./6.*xv+19./2.*xvs+4./3.*xvc)
      wcom01= wcom01+ymvv*ypiq*(-11./3.*xvv)
      wcom01= wcom01+ymvv*ypic*(23./3.*xvq+7*xvv)
      wcom01= wcom01+ymvv*ypis*(-11./2.*xvc-15*xvq-7./2.*xvv)
      wcom01= wcom01+ymvv*ypi*(11./6.*xvs+23./2.*xvc+8*xvq+1./3.*xvv)
      wcom01= wcom01+ymvv*omyp*(-1./6.*xv-19./12.*xvs-7./6.*xvc)
      wcom01= wcom01+ymvv*omyps*(-2./3.*xvs)
      wcom01= wcom01+ymvv*omypc*(-1./6.*xv)
      wcom01= wcom01+ymvv*(-23./12.*xvs-23./4.*xvc-xvq)
      wcom01= ymvv*wcom01
*
      wcom02= ypiq*(-11./3.*xvc)
      wcom02= wcom02+ypic*(7*xvs+7*xvc)
      wcom02= wcom02+ypis*(-9./2.*xv-13*xvs-7./2.*xvc)
      wcom02= wcom02+ypi*(5./6.+17./2.*xv+6*xvs+1./3.*xvc)
      wcom02= wcom02+omyp*(-11./12.-1./6.*xv)
      wcom02= wcom02+(-11./12.-15./4.*xv-1./3.*xvs)
      wcom02= rwms*rwmpgs*wcom02
*
      wcom0= wcom01+wcom02
*
      wcomdm= ypis*vvi*rwms*rwmpgs*xvs
      wcomdm= wcomdm+ypi*vvi*rwms*rwmpgs*(-2*xv-2*xvs)
      wcomdm= wcomdm+vvi*rwms*rwmpgs*omyp*(2*xv)
      wcomdm= wcomdm+vvi*rwms*rwmpgs*omyps
      wcomdm= wcomdm+vvi*rwms*rwmpgs*(2*xv+xvs)
      wcomdm= wcomdm+ymvv*ypis*vvi*rwms*(-2*xvc)
      wcomdm= wcomdm+ymvv*ypi*vvi*rwms*(4*xvs+4*xvc)
      wcomdm= wcomdm+ymvv*vvi*rwms*omyp*(-4*xvs)
      wcomdm= wcomdm+ymvv*vvi*rwms*omyps*(-2*xv)
      wcomdm= wcomdm+ymvv*vvi*rwms*(-4*xvs-2*xvc)
      wcomdm= wcomdm+ymvvs*ypis*vvi*xvq
      wcomdm= wcomdm+ymvvs*ypi*vvi*(-2*xvc-2*xvq)
      wcomdm= wcomdm+ymvvs*vvi*omyp*(2*xvc)
      wcomdm= wcomdm+ymvvs*vvi*omyps*xvs
      wcomdm= wcomdm+ymvvs*vvi*(2*xvc+xvq)
      dendm= dend*dms/vv
      wcomdm= ymvv*dendm*wcomdm
*
      wcom= wcomd*rlncd+wcom0+wcomdm
*
      dw1= ymvv*ypivi*rlncd*vvis*rwms*(40*xvv)
      dw1= dw1+ymvv*ypivi*rlncd*rwms*(40*xvv)
      dw1= dw1+ymvv*ypivi*vvis*rwms*(-260./3.*xvv)
      dw1= dw1+ymvv*ypivi*rwms*(-296./3.*xvv)
      dw1= dw1+ymvv*ypiv*rlncd*vvis*rwms*(-104*xvq-80*xvv)
      dw1= dw1+ymvv*ypiv*rlncd*rwms*(-120*xvq-72*xvv)
      dw1= dw1+ymvv*ypiv*vvis*rwms*(640./3.*xvq+520./3.*xvv)
      dw1= dw1+ymvv*ypiv*rwms*(296*xvq+172*xvv)
      dw1= dw1+ymvv*ypiq*rlncd*vvis*rwms*(89*xvc
     #    +218*xvq+48*xvv)
      dw1= dw1+ymvv*ypiq*rlncd*rwms*(141*xvc+214*xvq+36*xvv)
      dw1= dw1+ymvv*ypiq*vvis*rwms*(-511./3.*xvc
     #    -1366./3.*xvq-108*xvv)
      dw1= dw1+ymvv*ypiq*rwms*(-339*xvc-1526./3.*xvq-86*xvv)
      dw1= dw1+ymvv*ypic*rlncd*vvis*rwms*(-26*
     #    xvs-202*xvc-140*xvq-8*xvv)
      dw1= dw1+ymvv*ypic*rlncd*rwms*(-78*xvs-252*
     #    xvc-104*xvq-4*xvv)
      dw1= dw1+ymvv*ypic*vvis*rwms*(133./3.*xvs+
     #    1217./3.*xvc+310*xvq+64./3.*xvv)
      dw1= dw1+ymvv*ypic*rwms*(535./3.*xvs+1753./3.*
     #    xvc+244*xvq+34./3.*xvv)
      dw1= dw1+ymvv*ypis*rlncd*vvis*rwms*(xv+70*
     #    xvs+145*xvc+26*xvq)
      dw1= dw1+ymvv*ypis*rlncd*rwms*(18*xv+142*xvs
     #    +119*xvc+10*xvq)
      dw1= dw1+ymvv*ypis*vvis*rwms*(-1./6.*xv-
     #    805./6.*xvs-635./2.*xvc-209./3.*xvq-2./3.*xvv)
      dw1= dw1+ymvv*ypis*rwms*(-81./2.*xv-316*xvs
     #    -545./2.*xvc-27*xvq)
      dw1= dw1+ymvv*ypi*rlncd*vvis*rwms*(-6*xv-62*xvs-32*xvc)
      dw1= dw1+ymvv*ypi*rlncd*rwms*(-1-34*xv-66*xvs-8*xvc)
      dw1= dw1+ymvv*ypi*vvis*rwms*(-1./3.+29./3.
     #    *xv+135*xvs+257./3.*xvc+8./3.*xvq)
      dw1= dw1+ymvv*ypi*rwms*(3+229./3.*xv+147*
     #    xvs+25*xvc-2./3.*xvq)
      dw1= dw1+ymvv*rlncd*vvis*rwms*omyp*(4*xv)
      dw1= dw1+ymvv*rlncd*vvis*rwms*(5*xv+18*xvs)
      dw1= dw1+ymvv*rlncd*rwms*omyp
      dw1= dw1+ymvv*rlncd*rwms*(1+16*xv+2*xvs)
      dw1= dw1+ymvv*vvis*rwms*omyp*(1./3.-9*xv-8./3.*xvs)
      dw1= dw1+ymvv*vvis*rwms*omyps*(-2./3.*xv)
      dw1= dw1+ymvv*vvis*rwms*(1./3.-19./2.*xv-45*xvs-4*xvc)
      dw1= dw1+ymvv*rwms*omyp*(-19./6.-3*xv+2./3.*xvs)
      dw1= dw1+ymvv*rwms*omyps*(-1./3.+1./3.*xv)
      dw1= dw1+ymvv*rwms*(-10./3.-34*xv-12*xvs+xvc)
      dw1= dw1+ymvvs*ypivi*rlncd*vvis*(-20*xvvi)
      dw1= dw1+ymvvs*ypivi*rlncd*(-20*xvvi)
      dw1= dw1+ymvvs*ypivi*vvis*(130./3.*xvvi)
      dw1= dw1+ymvvs*ypivi*(148./3.*xvvi)
      dw1= dw1+ymvvs*ypiv*rlncd*vvis*(52*xvv+40*xvvi)
      dw1= dw1+ymvvs*ypiv*rlncd*(64*xvv+36*xvvi)
      dw1= dw1+ymvvs*ypiv*vvis*(-320./3.*xvv-260./3.*xvvi)
      dw1= dw1+ymvvs*ypiv*(-482./3.*xvv-86*xvvi)
      dw1= dw1+ymvvs*ypiq*rlncd*vvis*(-44*xvq-110*xvv-24*xvvi)
      dw1= dw1+ymvvs*ypiq*rlncd*(-81*xvq-116*xvv-18*xvvi)
      dw1= dw1+ymvvs*ypiq*vvis*(250./3.*xvq+694./
     #    3.*xvv+54*xvvi)
      dw1= dw1+ymvvs*ypiq*(203*xvq+844./3.*xvv+43*xvvi)
      dw1= dw1+ymvvs*ypic*rlncd*vvis*(12*xvc+103
     #    *xvq+72*xvv+4*xvvi)
      dw1= dw1+ymvvs*ypic*rlncd*(49*xvc+149*xvq+58*xvv+2*xvvi)
      dw1= dw1+ymvvs*ypic*vvis*(-56./3.*xvc-
     #    629./3.*xvq-162*xvv-32./3.*xvvi)
      dw1= dw1+ymvvs*ypic*(-364./3.*xvc-1084./3.*
     #    xvq-139*xvv-17./3.*xvvi)
      dw1= dw1+ymvvs*ypis*rlncd*vvis*(-36*xvc-78*xvq-14*xvv)
      dw1= dw1+ymvvs*ypis*rlncd*(-13*xvs-92*xvc-74*xvq-6*xvv)
      dw1= dw1+ymvvs*ypis*vvis*(-5./3.*xvs+419.
     #    /6.*xvc+177*xvq+115./3.*xvv+1./3.*xvvi)
      dw1= dw1+ymvvs*ypis*(100./3.*xvs+223*xvc+
     #    177*xvq+49./3.*xvv)
      dw1= dw1+ymvvs*ypi*rlncd*vvis*(3*xvs+36*xvc+19*xvq)
      dw1= dw1+ymvvs*ypi*rlncd*(xv+25*xvs+45*xvc+6*xvq)
      dw1= dw1+ymvvs*ypi*vvis*(1./3.*xv-25./6.*
     #    xvs-83*xvc-319./6.*xvq-5./3.*xvv)
      dw1= dw1+ymvvs*ypi*(-4*xv-391./6.*xvs-221.
     #    /2.*xvc-18*xvq+1./3.*xvv)
      dw1= dw1+ymvvs*rlncd*vvis*omyp*(-3*xvs)
      dw1= dw1+ymvvs*rlncd*vvis*(-3*xvs-12*xvc)
      dw1= dw1+ymvvs*rlncd*omyp*(-xv)
      dw1= dw1+ymvvs*rlncd*(-xv-12*xvs-2*xvc)
      dw1= dw1+ymvvs*vvis*omyp*(-1./24.+1./12.*xv+
     #    35./6.*xvs+10./3.*xvc)
      dw1= dw1+ymvvs*vvis*omyps*(-2./3.*xv+5./3.*xvs)
      dw1= dw1+ymvvs*vvis*omypc*(1./3.*xv)
      dw1= dw1+ymvvs*vvis*(-5./12.*xv+25./4.*xvs
     #    +31*xvc+10./3.*xvq)
      dw1= dw1+ymvvs*omyp*(5./24.+13./4.*xv+29./6.*xvs
     #    -5./6.*xvc)
      dw1= dw1+ymvvs*omyps*(3./2.*xv-1./2.*xvs)
      dw1= dw1+ymvvs*omypc*(1./6.-1./6.*xv)
      dw1= dw1+ymvvs*(17./4.*xv+365./12.*xvs+67./6.*
     #    xvc-5./6.*xvq)
      dw1= dw1+ypivi*rlncd*vvis*rwms*rwmpgs*(-20*xvq)
      dw1= dw1+ypivi*rlncd*rwms*rwmpgs*(-20*xvq)
      dw1= dw1+ypivi*vvis*rwms*rwmpgs*(130./3.*xvq)
      dw1= dw1+ypivi*rwms*rwmpgs*(148./3.*xvq)
      dw1= dw1+ypiv*rlncd*vvis*rwms*rwmpgs*(52*xvc+40*xvq)
      dw1= dw1+ypiv*rlncd*rwms*rwmpgs*(56*xvc+36*xvq)
      dw1= dw1+ypiv*vvis*rwms*rwmpgs*(-320./3.*
     #    xvc-260./3.*xvq)
      dw1= dw1+ypiv*rwms*rwmpgs*(-406./3.*xvc-86*xvq)
      dw1= dw1+ypiq*rlncd*vvis*rwms*rwmpgs*(-45*
     #    xvs-108*xvc-24*xvq)
      dw1= dw1+ypiq*rlncd*rwms*rwmpgs*(-60*xvs-98*
     #    xvc-18*xvq)
      dw1= dw1+ypiq*vvis*rwms*rwmpgs*(87*xvs+224*xvc+54*xvq)
      dw1= dw1+ypiq*rwms*rwmpgs*(136*xvs+682./3.*xvc+43*xvq)
      dw1= dw1+ypic*rlncd*vvis*rwms*rwmpgs*(14*xv+
     #    99*xvs+68*xvc+4*xvq)
      dw1= dw1+ypic*rlncd*rwms*rwmpgs*(29*xv+103*xvs
     #    +46*xvc+2*xvq)
      dw1= dw1+ypic*vvis*rwms*rwmpgs*(-77./3.*xv
     #    -196*xvs-148*xvc-32./3.*xvq)
      dw1= dw1+ypic*rwms*rwmpgs*(-173./3.*xv-223*
     #    xvs-105*xvc-17./3.*xvq)
      dw1= dw1+ypis*rlncd*vvis*rwms*rwmpgs*(-1-
     #    34*xv-67*xvs-12*xvc)
      dw1= dw1+ypis*rlncd*rwms*rwmpgs*(-5-50*xv-45*xvs-4*xvc)
      dw1= dw1+ypis*vvis*rwms*rwmpgs*(11./6.+193./
     #    3.*xv+281./2.*xvs+94./3.*xvc+1./3.*xvq)
      dw1= dw1+ypis*rwms*rwmpgs*(49./6.+95*xv+191./2.
     #    *xvs+32./3.*xvc)
      dw1= dw1+ypi*rlncd*vvis*rwms*rwmpgs*(3+26*xv+13*xvs)
      dw1= dw1+ypi*rlncd*rwms*rwmpgs*(9+21*xv+2*xvs)
      dw1= dw1+ypi*vvis*rwms*rwmpgs*(-11./2.-52
     #    *xv-65./2.*xvs-xvc)
      dw1= dw1+ypi*rwms*rwmpgs*(-85./6.-77./2.*xv
     #    -7*xvs+1./3.*xvc)
      dw1= dw1+rlncd*vvis*rwms*rwmpgs*omyp*(-1)
      dw1= dw1+rlncd*vvis*rwms*rwmpgs*(-2-6*xv)
      dw1= dw1+rlncd*rwms*rwmpgs*(-4)
      dw1= dw1+vvis*rwms*rwmpgs*omyp*(11./6.+1./3.*xv)
      dw1= dw1+vvis*rwms*rwmpgs*(11./3.+40./3.*xv+xvs)
      dw1= dw1+rwms*rwmpgs*omyp*(1./6.-1./6.*xv)
      dw1= dw1+rwms*rwmpgs*(17./3.+3./2.*xv-1./6.*xvs)
      dw1= dw1+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-xvq)
      dw1= dw1+dend*ymvvs*ypic*vviq*rwms*rwmpgs*(3*xvc+2*xvq)
      dw1= dw1+dend*ymvvs*ypis*vviq*rwms*rwmpgs*
     # (-3*xvs-6*xvc-xvq)
      dw1= dw1+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./2.*xvs)
      dw1= dw1+dend*ymvvs*ypi*vviq*rwms*rwmpgs*
     # (xv+6*xvs+3*xvc)
      dw1= dw1+dend*ymvvs*ypi*vvis*rwms*rwmpgs*(xv+xvs)
      dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-xv)
      dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*(-xv-3*xvs)
      dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-xv)
      dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./2.)
      dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*(-xv-1./2.*xvs)
      dw1= dw1+dend*ymvvc*ypiq*vviq*rwms*(2*xvv)
      dw1= dw1+dend*ymvvc*ypic*vviq*rwms*(-6*xvq-4*xvv)
      dw1= dw1+dend*ymvvc*ypis*vviq*rwms*(6*xvc+12*xvq+2*xvv)
      dw1= dw1+dend*ymvvc*ypis*vvis*rwms*xvc
      dw1= dw1+dend*ymvvc*ypi*vviq*rwms*(-2*xvs-12*xvc-6*xvq)
      dw1= dw1+dend*ymvvc*ypi*vvis*rwms*(-2*xvs-2*xvc)
      dw1= dw1+dend*ymvvc*vviq*rwms*omyp*(2*xvs)
      dw1= dw1+dend*ymvvc*vviq*rwms*(2*xvs+6*xvc)
      dw1= dw1+dend*ymvvc*vvis*rwms*omyp*(2*xvs)
      dw1= dw1+dend*ymvvc*vvis*rwms*omyps*(xv)
      dw1= dw1+dend*ymvvc*vvis*rwms*(2*xvs+xvc)
      dw1= dw1+dend*ymvvq*ypiq*vviq*(-xvvi)
      dw1= dw1+dend*ymvvq*ypic*vviq*(3*xvv+2*xvvi)
      dw1= dw1+dend*ymvvq*ypis*vviq*(-3*xvq-6*xvv-xvvi)
      dw1= dw1+dend*ymvvq*ypis*vvis*(-1./2.*xvq)
      dw1= dw1+dend*ymvvq*ypi*vviq*(xvc+6*xvq+3*xvv)
      dw1= dw1+dend*ymvvq*ypi*vvis*(xvc+xvq)
      dw1= dw1+dend*ymvvq*vviq*omyp*(-xvc)
      dw1= dw1+dend*ymvvq*vviq*(-xvc-3*xvq)
      dw1= dw1+dend*ymvvq*vvis*omyp*(-xvc)
      dw1= dw1+dend*ymvvq*vvis*omyps*(-1./2.*xvs)
      dw1= dw1+dend*ymvvq*vvis*(-xvc-1./2.*xvq)

      dw2= ymvv*ypivi*rlncd*vvis*rwms*(120*xvv)
      dw2= dw2+ymvv*ypivi*rlncd*rwms*(40*xvv)
      dw2= dw2+ymvv*ypivi*vvis*rwms*(-260*xvv)
      dw2= dw2+ymvv*ypivi*rwms*(-296./3.*xvv)
      dw2= dw2+ymvv*ypiv*rlncd*vvis*rwms*(-312*xvq-240*xvv)
      dw2= dw2+ymvv*ypiv*rlncd*rwms*(-120*xvq-72*xvv)
      dw2= dw2+ymvv*ypiv*vvis*rwms*(640*xvq+520*xvv)
      dw2= dw2+ymvv*ypiv*rwms*(296*xvq+172*xvv)
      dw2= dw2+ymvv*ypiq*rlncd*vvis*rwms*(267*xvc
     #    +654*xvq+144*xvv)
      dw2= dw2+ymvv*ypiq*rlncd*rwms*(141*xvc+214*xvq+36*xvv)
      dw2= dw2+ymvv*ypiq*vvis*rwms*(-511*xvc-
     #    1366*xvq-324*xvv)
      dw2= dw2+ymvv*ypiq*rwms*(-339*xvc-1526./3.*xvq-86*xvv)
      dw2= dw2+ymvv*ypic*rlncd*vvis*rwms*(-78*
     #    xvs-606*xvc-420*xvq-24*xvv)
      dw2= dw2+ymvv*ypic*rlncd*rwms*(-78*xvs-252*
     #    xvc-104*xvq-4*xvv)
      dw2= dw2+ymvv*ypic*vvis*rwms*(133*xvs+1217
     #    *xvc+930*xvq+64*xvv)
      dw2= dw2+ymvv*ypic*rwms*(535./3.*xvs+1753./3.*
     #    xvc+244*xvq+34./3.*xvv)
      dw2= dw2+ymvv*ypis*rlncd*vvis*rwms*(3*xv+210
     #    *xvs+435*xvc+78*xvq)
      dw2= dw2+ymvv*ypis*rlncd*rwms*(18*xv+142*xvs
     #    +119*xvc+10*xvq)
      dw2= dw2+ymvv*ypis*vvis*rwms*(-1./2.*xv-
     #    805./2.*xvs-1905./2.*xvc-209*xvq-2*xvv)
      dw2= dw2+ymvv*ypis*rwms*(-81./2.*xv-316*xvs
     #    -545./2.*xvc-27*xvq)
      dw2= dw2+ymvv*ypi*rlncd*vvis*rwms*(-18*xv-186*xvs-96*xvc)
      dw2= dw2+ymvv*ypi*rlncd*rwms*(-1-34*xv-66*xvs-8*xvc)
      dw2= dw2+ymvv*ypi*vvis*rwms*(-1+29*xv+
     #    405*xvs+257*xvc+8*xvq)
      dw2= dw2+ymvv*ypi*rwms*(3+229./3.*xv+147*
     #    xvs+25*xvc-2./3.*xvq)
      dw2= dw2+ymvv*rlncd*vvis*rwms*omyp*(12*xv)
      dw2= dw2+ymvv*rlncd*vvis*rwms*(15*xv+54*xvs)
      dw2= dw2+ymvv*rlncd*rwms*omyp
      dw2= dw2+ymvv*rlncd*rwms*(1+16*xv+2*xvs)
      dw2= dw2+ymvv*vvis*rwms*omyp*(1-27*xv-8*xvs)
      dw2= dw2+ymvv*vvis*rwms*omyps*(-2*xv)
      dw2= dw2+ymvv*vvis*rwms*(1-57./2.*xv-135*xvs-12*xvc)
      dw2= dw2+ymvv*rwms*omyp*(-19./6.-3*xv+2./3.*xvs)
      dw2= dw2+ymvv*rwms*omyps*(-1./3.+1./3.*xv)
      dw2= dw2+ymvv*rwms*(-10./3.-34*xv-12*xvs+xvc)
      dw2= dw2+ymvvs*ypivi*rlncd*vvis*(-60*xvvi)
      dw2= dw2+ymvvs*ypivi*rlncd*(-20*xvvi)
      dw2= dw2+ymvvs*ypivi*vvis*(130*xvvi)
      dw2= dw2+ymvvs*ypivi*(148./3.*xvvi)
      dw2= dw2+ymvvs*ypiv*rlncd*vvis*(156*xvv+120*xvvi)
      dw2= dw2+ymvvs*ypiv*rlncd*(64*xvv+36*xvvi)
      dw2= dw2+ymvvs*ypiv*vvis*(-320*xvv-260*xvvi)
      dw2= dw2+ymvvs*ypiv*(-482./3.*xvv-86*xvvi)
      dw2= dw2+ymvvs*ypiq*rlncd*vvis*(-132*xvq
     #    -330*xvv-72*xvvi)
      dw2= dw2+ymvvs*ypiq*rlncd*(-81*xvq-116*xvv-18*xvvi)
      dw2= dw2+ymvvs*ypiq*vvis*(250*xvq+694*
     #    xvv+162*xvvi)
      dw2= dw2+ymvvs*ypiq*(203*xvq+844./3.*xvv+43*xvvi)
      dw2= dw2+ymvvs*ypic*rlncd*vvis*(36*xvc+309
     #    *xvq+216*xvv+12*xvvi)
      dw2= dw2+ymvvs*ypic*rlncd*(49*xvc+149*xvq+58*xvv+2*xvvi)
      dw2= dw2+ymvvs*ypic*vvis*(-56*xvc-629*
     #    xvq-486*xvv-32*xvvi)
      dw2= dw2+ymvvs*ypic*(-364./3.*xvc-1084./3.*
     #    xvq-139*xvv-17./3.*xvvi)
      dw2= dw2+ymvvs*ypis*rlncd*vvis*(-108*xvc
     #    -234*xvq-42*xvv)
      dw2= dw2+ymvvs*ypis*rlncd*(-13*xvs-92*xvc-74*xvq-6*xvv)
      dw2= dw2+ymvvs*ypis*vvis*(-5*xvs+419./2.
     #    *xvc+531*xvq+115*xvv+xvvi)
      dw2= dw2+ymvvs*ypis*(100./3.*xvs+223*xvc+
     #    177*xvq+49./3.*xvv)
      dw2= dw2+ymvvs*ypi*rlncd*vvis*(9*xvs+108*xvc+57*xvq)
      dw2= dw2+ymvvs*ypi*rlncd*(xv+25*xvs+45*xvc+6*xvq)
      dw2= dw2+ymvvs*ypi*vvis*(xv-25./2.*xvs-
     #    249*xvc-319./2.*xvq-5*xvv)
      dw2= dw2+ymvvs*ypi*(-4*xv-391./6.*xvs-221.
     #    /2.*xvc-18*xvq+1./3.*xvv)
      dw2= dw2+ymvvs*rlncd*vvis*omyp*(-9*xvs)
      dw2= dw2+ymvvs*rlncd*vvis*(-9*xvs-36*xvc)
      dw2= dw2+ymvvs*rlncd*omyp*(-xv)
      dw2= dw2+ymvvs*rlncd*(-xv-12*xvs-2*xvc)
      dw2= dw2+ymvvs*vvis*omyp*(-1./8.+1./4.*xv+35.
     #    /2.*xvs+10*xvc)
      dw2= dw2+ymvvs*vvis*omyps*(-2*xv+5*xvs)
      dw2= dw2+ymvvs*vvis*omypc*(xv)
      dw2= dw2+ymvvs*vvis*(-5./4.*xv+75./4.*xvs+
     #    93*xvc+10*xvq)
      dw2= dw2+ymvvs*omyp*(5./24.+13./4.*xv+29./6.*xvs-5./6.*xvc)
      dw2= dw2+ymvvs*omyps*(3./2.*xv-1./2.*xvs)
      dw2= dw2+ymvvs*omypc*(1./6.-1./6.*xv)
      dw2= dw2+ymvvs*(17./4.*xv+365./12.*xvs+67./6.*
     #    xvc-5./6.*xvq)
      dw2= dw2+ypivi*rlncd*vvis*rwms*rwmpgs*(-60*xvq)
      dw2= dw2+ypivi*rlncd*rwms*rwmpgs*(-20*xvq)
      dw2= dw2+ypivi*vvis*rwms*rwmpgs*(130*xvq)
      dw2= dw2+ypivi*rwms*rwmpgs*(148./3.*xvq)
      dw2= dw2+ypiv*rlncd*vvis*rwms*rwmpgs*(156*xvc+120*xvq)
      dw2= dw2+ypiv*rlncd*rwms*rwmpgs*(56*xvc+36*xvq)
      dw2= dw2+ypiv*vvis*rwms*rwmpgs*(-320*xvc-260*xvq)
      dw2= dw2+ypiv*rwms*rwmpgs*(-406./3.*xvc-86*xvq)
      dw2= dw2+ypiq*rlncd*vvis*rwms*rwmpgs*(-135*xvs-324*xvc-72*xvq)
      dw2= dw2+ypiq*rlncd*rwms*rwmpgs*(-60*xvs-98*xvc-18*xvq)
      dw2= dw2+ypiq*vvis*rwms*rwmpgs*(261*xvs+
     #    672*xvc+162*xvq)
      dw2= dw2+ypiq*rwms*rwmpgs*(136*xvs+682./3.*xvc+43*xvq)
      dw2= dw2+ypic*rlncd*vvis*rwms*rwmpgs*(42*xv+
     #    297*xvs+204*xvc+12*xvq)
      dw2= dw2+ypic*rlncd*rwms*rwmpgs*(29*xv+103*xvs+46*xvc+2*xvq)
      dw2= dw2+ypic*vvis*rwms*rwmpgs*(-77*xv-588
     #    *xvs-444*xvc-32*xvq)
      dw2= dw2+ypic*rwms*rwmpgs*(-173./3.*xv-223*
     #    xvs-105*xvc-17./3.*xvq)
      dw2= dw2+ypis*rlncd*vvis*rwms*rwmpgs*(-3-
     #    102*xv-201*xvs-36*xvc)
      dw2= dw2+ypis*rlncd*rwms*rwmpgs*(-5-50*xv-45*xvs-4*xvc)
      dw2= dw2+ypis*vvis*rwms*rwmpgs*(11./2.+193*
     #    xv+843./2.*xvs+94*xvc+xvq)
      dw2= dw2+ypis*rwms*rwmpgs*(49./6.+95*xv+191./2.
     #    *xvs+32./3.*xvc)
      dw2= dw2+ypi*rlncd*vvis*rwms*rwmpgs*(9+78*xv+39*xvs)
      dw2= dw2+ypi*rlncd*rwms*rwmpgs*(9+21*xv+2*xvs)
      dw2= dw2+ypi*vvis*rwms*rwmpgs*(-33./2.-
     #    156*xv-195./2.*xvs-3*xvc)
      dw2= dw2+ypi*rwms*rwmpgs*(-85./6.-77./2.*xv
     #    -7*xvs+1./3.*xvc)
      dw2= dw2+rlncd*vvis*rwms*rwmpgs*omyp*(-3)
      dw2= dw2+rlncd*vvis*rwms*rwmpgs*(-6-18*xv)
      dw2= dw2+rlncd*rwms*rwmpgs*(-4)
      dw2= dw2+vvis*rwms*rwmpgs*omyp*(11./2.+xv)
      dw2= dw2+vvis*rwms*rwmpgs*(11+40*xv+3*xvs)
      dw2= dw2+rwms*rwmpgs*omyp*(1./6.-1./6.*xv)
      dw2= dw2+rwms*rwmpgs*(17./3.+3./2.*xv-1./6.*xvs)
      dw2= dw2+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-3*xvq)
      dw2= dw2+dend*ymvvs*ypic*vviq*rwms*rwmpgs*(9*xvc+6*xvq)
      dw2= dw2+dend*ymvvs*ypis*vviq*rwms*rwmpgs*(-9*xvs-18*xvc-3*xvq)
      dw2= dw2+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./2.*xvs)
      dw2= dw2+dend*ymvvs*ypi*vviq*rwms*rwmpgs*(3*xv+18*xvs+9*xvc)
      dw2= dw2+dend*ymvvs*ypi*vvis*rwms*rwmpgs*(xv+xvs)
      dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-3*xv)
      dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*(-3*xv-9*xvs)
      dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-xv)
      dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./2.)
      dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*(-xv-1./2.*xvs)
      dw2= dw2+dend*ymvvc*ypiq*vviq*rwms*(6*xvv)
      dw2= dw2+dend*ymvvc*ypic*vviq*rwms*(-18*xvq-12*xvv)
      dw2= dw2+dend*ymvvc*ypis*vviq*rwms*(18*xvc+36*xvq+6*xvv)
      dw2= dw2+dend*ymvvc*ypis*vvis*rwms*(xvc)
      dw2= dw2+dend*ymvvc*ypi*vviq*rwms*(-6*xvs-36*xvc-18*xvq)
      dw2= dw2+dend*ymvvc*ypi*vvis*rwms*(-2*xvs-2*xvc)
      dw2= dw2+dend*ymvvc*vviq*rwms*omyp*(6*xvs)
      dw2= dw2+dend*ymvvc*vviq*rwms*(6*xvs+18*xvc)
      dw2= dw2+dend*ymvvc*vvis*rwms*omyp*(2*xvs)
      dw2= dw2+dend*ymvvc*vvis*rwms*omyps*(xv)
      dw2= dw2+dend*ymvvc*vvis*rwms*(2*xvs+xvc)
      dw2= dw2+dend*ymvvq*ypiq*vviq*(-3*xvvi)
      dw2= dw2+dend*ymvvq*ypic*vviq*(9*xvv+6*xvvi)
      dw2= dw2+dend*ymvvq*ypis*vviq*(-9*xvq-18*xvv-3*xvvi)
      dw2= dw2+dend*ymvvq*ypis*vvis*(-1./2.*xvq)
      dw2= dw2+dend*ymvvq*ypi*vviq*(3*xvc+18*xvq+9*xvv)
      dw2= dw2+dend*ymvvq*ypi*vvis*(xvc+xvq)
      dw2= dw2+dend*ymvvq*vviq*omyp*(-3*xvc)
      dw2= dw2+dend*ymvvq*vviq*(-3*xvc-9*xvq)
      dw2= dw2+dend*ymvvq*vvis*omyp*(-xvc)
      dw2= dw2+dend*ymvvq*vvis*omyps*(-1./2.*xvs)
      dw2= dw2+dend*ymvvq*vvis*(-xvc-1./2.*xvq)
*
      if(ii.eq.3) then
       wcom= wcom*(2.d0-ym)**2/yms-dw1+2.d0/ym*(1.d0-1.d0/ym)*dw2
      else if(ii.eq.4) then
       wcom= 0.5d0*em4*omxm*omxm*(dw1-dw2)     
      else if(ii.eq.5) then
       wcom= em2*((2.d0+0.5d0*omxm)*dw1+omxm*(2.d0/ym-
     #       3.d0/2.d0)*dw2)
      endif
*
      endif
*
      www= wcom/pwtm2
*
    1 if(iz.eq.0) then
         ww(ii)= 0.d0
         iz= 1
      else
         ww(ii)= www*tjac*tfactsw
      endif
*
      enddo
*
    2 if(iz.eq.0) then
       wwt= 0.d0
      else
       wwt= ww(1)+ww(2)+ww(3)+ww(4)+ww(5)
      endif
*
      if(wwt.lt.0.d0) then
       iz= 0
       ifz(7)= ifz(7)+1
       wtoxssw= 0.d0
      else
       wtoxssw= wwt
      endif
*
      return
      end      
*
*-------------------------------------------------------------------------
*
      real*8 function wtoxsswc(ndim,x)
      implicit real*8(a-h,o-z)
      character*1 opeak
*
      common/wtps/opeak
      common/wtkount/ik
      common/wtistrf/isf
      common/wtmult/inm,jnm
      common/wtipt/ifz(51)
      common/wtfsw/tfactsw
      common/wtnpr/ipr,ipr0
      common/wticuts/iac(4)
      common/wtcutsw/tcs,s0cut
      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/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)
*
      dimension vu(4),vut(4)
      dimension vd(4),vdt(4)
      dimension vn(4),vnt(4)
      dimension x(ndim),ww(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(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
      psf= 0.d0
      wwp= 0.d0
      tjac= 0.d0
      do kk=1,2
        ww(kk)= 0.d0
      enddo
*
      em2= em*em
      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)
*
      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
      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
      qhcs= qh0s+0.25d0*xm*xm*omy*tcs*s
*
      x0= qh0s/y/sh
      xc0= qhcs/y/sh
      xc1= omy/y
      xc= dmin1(xc0,xc1)
*      dxc= 0.25d0*xm*xm*omy*tcs/yvv
      if(xc.eq.xc0) then
       dxc= 0.25d0*xm*xm*omy*tcs/yvv
      else
       dxc= xc-x0
      endif
*
      if(x0.gt.xc) then
         iz= 0
         ifz(5)= ifz(5)+1
         go to 1
      endif
*
      ii= inm
      if(ii.eq.1) then
       xjac= dxc/xc/x0/y/sh
       xv= x0*xc/(dxc*xx+x0)
      else if(ii.eq.2) then
       xjac= log(xc/x0)
       xv= exp(xjac*xx+log(x0))
      else if(ii.eq.3) then
       xjac= dxc
       xv= xjac*xx+x0
      endif
      omxv= 1.d0-xv
      opxv= 1.d0+xv
      opxvs= opxv*opxv
      xvs= xv*xv
*
      if(ii.eq.1) then
       psf= em2*(2.d0+omxm*(2.d0/y-1.d0))
      else if(ii.eq.2) then
       psf= -1.d0+2.d0/y*(1.d0-1.d0/y)
      else if(ii.eq.3) then
       psf= 1.d0
      endif
*
*-----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-rwm2
      zmbs= zmb-rwm2
      atma= zmas/rwmg
      atmb= zmbs/rwmg
*
      if(opeak.eq.'y') then
       if(atma.gt.1.d0.and.atmb.gt.1.d0) then
          atma= 1.d0/atma
          atma= atan(atma)
          zmat= (pih-atma)/rwgm
          atmb= 1.d0/atmb
          atmb= atan(atmb)
          zmbt= (pih-atmb)/rwgm
          x2jc0= (-atmb+atma)/rwgm
       else if(atma.gt.1.d0.and.atmb.lt.-1.d0) then
          atma= 1.d0/atma
          atma= atan(atma)
          zmat= (pih-atma)/rwgm
          atmb= -1.d0/atmb
          atmb= atan(atmb)
          zmbt= (-pih+atmb)/rwgm
          x2jc0= (-pi+atmb+atma)/rwgm
       else if(atma.gt.1.d0.and.abs(atmb).lt.1.d0) then
          atma= 1.d0/atma
          atma= atan(atma)
          zmat= (pih-atma)/rwgm
          atmb= atan(atmb)
          zmbt= atmb/rwgm
          x2jc0= (-pih+atmb+atma)/rwgm
       else if(atma.lt.-1.d0.and.atmb.gt.1.d0) then
          atma= -1.d0/atma
          atma= atan(atma)
          zmat= (-pih+atma)/rwgm
          atmb= 1.d0/atmb
          atmb= atan(atmb)
          zmbt= (pih-atmb)/rwgm
          x2jc0= (pi-atmb-atma)/rwgm
       else if(atma.lt.-1.d0.and.atmb.lt.-1.d0) then
          atma= -1.d0/atma
          atma= atan(atma)
          zmat= (-pih+atma)/rwgm
          atmb= -1.d0/atmb
          atmb= atan(atmb)
          zmbt= (-pih+atmb)/rwgm
          x2jc0= (atmb-atma)/rwgm
       else if(atma.lt.-1.d0.and.abs(atmb).lt.1.d0) then
          atma= -1.d0/atma
          atma= atan(atma)
          zmat= (-pih+atma)/rwgm
          atmb= atan(atmb)
          zmbt= atmb/rwgm
          x2jc0= (pih+atmb-atma)/rwgm
       else if(abs(atma).lt.1.d0.and.atmb.gt.1.d0) then
          atma= atan(atma)
          zmat= atma/rwgm
          atmb= 1.d0/atmb
          atmb= atan(atmb)
          zmbt= (pih-atmb)/rwgm
          x2jc0= (pih-atmb-atma)/rwgm
       else if(abs(atma).lt.1.d0.and.atmb.lt.-1.d0) then
          atma= atan(atma)
          zmat= atma/rwgm
          atmb= -1.d0/atmb
          atmb= atan(atmb)
          zmbt= (-pih+atmb)/rwgm
          x2jc0= (-pih+atmb-atma)/rwgm
       else if(abs(atma).lt.1.d0.and.abs(atmb).lt.1.d0) then
          atma= atan(atma)
          zmat= atma/rwgm
          atmb= atan(atmb)
          zmbt= atmb/rwgm
          x2jc0= (atmb-atma)/rwgm
       endif
*
       x2jc= x2jc0/yvv
       argx2= rwgm*(x2jc0*x2x+zmat)
       iftn= 1
       s07= s07aaf(argx2,iftn)
       x2= (rwms+rwgm*s07)/yvv
       pmjc= 1.d0
       pwsm2= rwgms*(1.d0+s07*s07)
*
      else if(opeak.eq.'n') then         
       smjc0= zmb-zma
       x2= (smjc0*x2x+zma)/yvv
       pmjc= 1.d0/((yvv*x2-rwm2)**2+rwmgs)
       x2jc= smjc0/yvv
       pwsm2= 1.d0/pmjc
      endif
      omx2= 1.d0-x2
      omx2s= omx2*omx2
      x2s= x2*x2
      x2c= x2s*x2
      x2q= x2c*x2
*
*-----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(10)= ifz(10)+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(12)= ifz(12)+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(13)= ifz(13)+1
         go to 1
       endif
       if(jj.ge.4.and.(ups-dmsh).gt.0.d0) then
         iz= 0
         ifz(13)= ifz(13)+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(13)= ifz(13)+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(14)= ifz(14)+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(15)= ifz(15)+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(18)= ifz(18)+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(19)= ifz(19)+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(20)= ifz(20)+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(21)= ifz(21)+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)
      pwtm2= ((1.d0+xv-x2+z)*yvv+rwms)**2+rwms*rwgs
*
      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(22)= ifz(22)+1
       go to 1
      endif
*
*-----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(23)= ifz(23)+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(24)= ifz(24)+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(25)= ifz(25)+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(26)= ifz(26)+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(28)= ifz(28)+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(29)= ifz(29)+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(30)= ifz(30)+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(31)= ifz(31)+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(32)= ifz(32)+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(33)= ifz(33)+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(34)= ifz(34)+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(35)= ifz(35)+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(36)= ifz(36)+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(37)= ifz(37)+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(38)= ifz(38)+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(39)= ifz(39)+1
       go to 2
      endif
      vels= xv*(omy-xv*y)
      vel= sqrt(vels)
      if(abs(vel).gt.1.d0) then
       iz= 0
       ifz(40)= ifz(40)+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(41)= ifz(41)+1
        go to 2
       endif
       if(eut.lt.eutmn) then
        iz= 0
        ifz(42)= ifz(42)+1
        go to 2
       endif
       if(edt.lt.edtmn) then
        iz= 0
        ifz(43)= ifz(43)+1
        go to 2
       endif
       if(ctut.lt.asa(3).or.ctut.gt.bsa(3)) then
        iz= 0
        ifz(44)= ifz(44)+1
        go to 2
       endif
       if(ctdt.lt.asa(4).or.ctdt.gt.bsa(4)) then
        iz= 0
        ifz(45)= ifz(45)+1
        go to 2
       endif
      else if(ipro.eq.1.or.ipro.eq.3) then
       if(edt.lt.edtmn) then
        iz= 0
        ifz(46)= ifz(46)+1
        go to 2
       endif
       if(ctdt.lt.asa(4).or.ctdt.gt.bsa(4)) then
        iz= 0
        ifz(47)= ifz(47)+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(47)= ifz(47)+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
*
      if(ipro.eq.2.and.ii.le.2) then
       if(jj.eq.1) then
*
      wcom= pwsm2*y*zxi*(2./9.*vvs*tau*x2-2./9.*
     #    vvs*tau+4./9.*vvs*x1*x2-4./9.*vvs*x1+2./9.*vvs)
      wcom= wcom+pwsm2*y*(-2./9.*vvs*x1)
      wcom= wcom+pwtm2*y*(-vvs*tau*x1+vvs*tau)
      wcom= wcom+y*rwms*rwmpgs*(4./3.*vvs*x1
     #    -4./3.*vvs*x1s)
      wcom= wcom+ys*rwms*(2*vvc*tau*x1*x2+
     #    vvc*tau*x1*z-2./3.*vvc*tau*x1+2./3.*vvc*tau*x2-2*
     #    vvc*tau*z-2*vvc*tau+vvc*taus*x2-vvc*taus-
     #    10./3.*vvc*x1*x2-1./3.*vvc*x1*z+4./3.*vvc*x1+14./
     #    3.*vvc*x1s*x2-4./3.*vvc*x1s*z-3*vvc*x1s+
     #    vvc*x2-2./3.*vvc)
      wcom= wcom+yc*(vvq*tau*x1*x2*z+2./3.*
     #    vvq*tau*x1*x2-2*vvq*tau*x1*x2s+2*vvq*tau*x1*z+
     #    vvq*tau*x1*zs+vvq*tau*x1+2*vvq*tau*x2*z+2*vvq
     #    *tau*x2-2./3.*vvq*tau*x2s-2*vvq*tau*z-vvq*tau*
     #    zs-vvq*tau+vvq*taus*x2*z+vvq*taus*x2-
     #    vvq*taus*x2s-5./3.*vvq*x1*x2*z-4./3.*vvq*x1*x2
     #    +2*vvq*x1*x2s+10./3.*vvq*x1s*x2*z+3*vvq*x1s
     #    *x2-10./3.*vvq*x1s*x2s+vvq*x2*z+2./3.*vvq*x2
     #    -vvq*x2s)
*
       else if(jj.eq.2) then
*
      wcom= pwsm2*y*zxi*(-4./9.*vvs*x1*x2+4./9.*
     #    vvs*x1s*x2+2./9.*vvs*x2)
      wcom= wcom+pwsm2*y*(-2./9.*vvs*x1-4./9.*
     #    vvs*x1s+4./9.*vvs)
      wcom= wcom+y*rwms*rwmpgs*(-4./3.*vvs*x1
     #    +8./3.*vvs*x1s-4./3.*vvs*x1c)
      wcom= wcom+ys*rwms*(4./3.*vvc*x1*x2-2./
     #    3.*vvc*x1*z-2./3.*vvc*x1-4*vvc*x1s*x2+4./3.*
     #    vvc*x1s*z+4./3.*vvc*x1s+8./3.*vvc*x1c*x2-4./
     #    3.*vvc*x1c*z-4./3.*vvc*x1c+2./3.*vvc*x2)
      wcom= wcom+yc*(2./3.*vvq*x1*x2*z+2./3.*
     #    vvq*x1*x2-4./3.*vvq*x1s*x2*z-4./3.*vvq*x1s*x2
     #    +4./3.*vvq*x1s*x2s+4./3.*vvq*x1c*x2*z+4./3.*
     #    vvq*x1c*x2-4./3.*vvq*x1c*x2s-2./3.*vvq*x2s)
*
       else if(jj.eq.3) then
*
      wcom= pwsm2*rums*(16./9.*vv*x1-8./9.*vv*x1s-8./9.*vv)
*
       else if(jj.eq.4) then
*
      wcom= pwsm2*y*zxi*(2./9.*vvs*ups*x2-2./9.*
     #    vvs*ups-4./9.*vvs*x1*x2+4./9.*vvs*x1-2./9.*vvs)
      wcom= wcom+pwsm2*y*(2./9.*vvs*x1+4./9.*
     #    vvs*x2+1./9.*vvs*z-4./9.*vvs)
      wcom= wcom+y*rwms*rwmpgs*(-1./3.*vvs*
     #    ups*x2+1./3.*vvs*ups+2./3.*vvs*x1*x2+1./3.*vvs*x1
     #    *z+1./3.*vvs*x1-2./3.*vvs*x1s-2./3.*vvs*x2*z-
     #    2./3.*vvs*x2+1./3.*vvs)
      wcom= wcom+ys*rwms*(1./3.*vvc*ups*x1
     #    -1./3.*vvc*ups*x2*z-2./3.*vvc*ups*x2+2./3.*vvc*
     #    ups*x2s-1./3.*vvc*ups*z-1./3.*vvc*ups-2./3.*vvc
     #    *x1*x2-4./3.*vvc*x1*x2s+2*vvc*x1*z+1./3.*vvc*x1
     #    *zs+2*vvc*x1+4./3.*vvc*x1s*x2-2./3.*vvc*x1s
     #    *z-vvc*x1s-2./3.*vvc*x2*z-2./3.*vvc*x2*zs-2.
     #    /3.*vvc*x2+4./3.*vvc*x2s*z+4./3.*vvc*x2s-4./3.
     #    *vvc*z-2./3.*vvc*zs-vvc)
      wcom= wcom+yc*(-1./3.*vvq*ups*x1*x2+
     #    1./3.*vvq*ups*x2*z+1./3.*vvq*ups*x2+1./3.*vvq*ups*
     #    x2s*z+1./3.*vvq*ups*x2s-1./3.*vvq*ups*x2c-2*
     #    vvq*x1*x2*z-1./3.*vvq*x1*x2*zs-2*vvq*x1*x2-1./3.
     #    *vvq*x1*x2s*z+1./3.*vvq*x1*x2s+2./3.*vvq*x1*
     #    x2c+2./3.*vvq*x1s*x2*z+vvq*x1s*x2-2./3.*vvq
     #    *x1s*x2s+4./3.*vvq*x2*z+2./3.*vvq*x2*zs+vvq
     #    *x2+2./3.*vvq*x2s*z+2./3.*vvq*x2s*zs+1./3.*
     #    vvq*x2s-2./3.*vvq*x2c*z-2./3.*vvq*x2c)
*
       else if(jj.eq.5) then
*
      wcom= pwsm2*y*zxi*(-4./9.*vvs*x1*x2+4./9.*
     #    vvs*x1s*x2+2./9.*vvs*x2)
      wcom= wcom+pwsm2*y*(-4./9.*vvs*x1*x2+1./
     #    9.*vvs*x1*z+8./9.*vvs*x1-4./9.*vvs*x1s+2./9.*
     #    vvs*x2*z+2./9.*vvs*x2-1./9.*vvs*z+1./9.*vvs*
     #    zs-4./9.*vvs)
      wcom= wcom+y*rwms*rwmpgs*(2./3.*vvs*x1*x2*
     #    z+4./3.*vvs*x1*x2+vvs*x1*z+1./3.*vvs*x1*zs+2./
     #    3.*vvs*x1-2./3.*vvs*x1s*x2-2./3.*vvs*x1s*z-4./
     #    3.*vvs*x1s+2./3.*vvs*x1c-2./3.*vvs*x2*z-1./3.
     #    *vvs*x2*zs-2./3.*vvs*x2-1./3.*vvs*z-1./3.*vvs
     #    *zs)
      wcom= wcom+ys*rwms*(-2./3.*vvc*x1*x2*z
     #    -4./3.*vvc*x1*x2-4./3.*vvc*x1*x2s*z-8./3.*vvc*
     #    x1*x2s+10./3.*vvc*x1*z+5./3.*vvc*x1*zs+1./3.*
     #    vvc*x1*zc+2*vvc*x1+2./3.*vvc*x1s*x2*z+8./3.*
     #    vvc*x1s*x2+4./3.*vvc*x1s*x2s-8./3.*vvc*x1s*
     #    z-2./3.*vvc*x1s*zs-2*vvc*x1s-4./3.*vvc*
     #    x1c*x2+2./3.*vvc*x1c*z+2./3.*vvc*x1c-1./3.*
     #    vvc*x2*zc+4./3.*vvc*x2s*z+2./3.*vvc*x2s*zs
     #    +4./3.*vvc*x2s-4./3.*vvc*z-vvc*zs-1./3.*
     #    vvc*zc-2./3.*vvc)
      wcom= wcom+yc*(-10./3.*vvq*x1*x2*z-5./
     #    3.*vvq*x1*x2*zs-1./3.*vvq*x1*x2*zc-2*vvq*x1*x2
     #    -1./3.*vvq*x1*x2s*z-1./3.*vvq*x1*x2s*zs+2./3.
     #    *vvq*x1*x2s+2./3.*vvq*x1*x2c*z+4./3.*vvq*x1*
     #    x2c+8./3.*vvq*x1s*x2*z+2./3.*vvq*x1s*x2*zs+
     #    2*vvq*x1s*x2-4./3.*vvq*x1s*x2s-2./3.*vvq*
     #    x1s*x2c-2./3.*vvq*x1c*x2*z-2./3.*vvq*x1c*x2
     #    +2./3.*vvq*x1c*x2s+4./3.*vvq*x2*z+vvq*x2*zs
     #    +1./3.*vvq*x2*zc+2./3.*vvq*x2+1./3.*vvq*x2s*z
     #    +1./3.*vvq*x2s*zs+1./3.*vvq*x2s*zc-2./3.*
     #    vvq*x2c*z-1./3.*vvq*x2c*zs-2./3.*vvq*x2c)
*
       else if(jj.eq.6) then
*
      wcom= pwsm2*rdms*(4./9.*vv*x1*z+4./9.*vv*x1-2./9.
     #    *vv*x1s-4./9.*vv*z-2./9.*vv*zs-2./9.*vv)
*
       endif
      else if(ipro.eq.2.and.ii.eq.3) then
       if(jj.eq.1) then
*
      wcom= pwsm2*y*zxi*(2./9.*vvs*tau*x2-2./9.*
     #    vvs*tau+4./9.*vvs*x1*x2-4./9.*vvs*x1+2./9.*vvs)
      wcom= wcom+pwsm2*y*(-2./9.*vvs*x1)
      wcom= wcom+pwtm2*y*(-vvs*tau*x1+vvs*tau)
      wcom= wcom+y*rwms*rwmpgs*(4./3.*vvs*x1-4./3.*vvs*x1s)
      wcom= wcom+ys*rwms*(2*vvc*tau*x1*x2+
     #    vvc*tau*x1*z-2./3.*vvc*tau*x1+2./3.*vvc*tau*x2-2*
     #    vvc*tau*z-2*vvc*tau+vvc*taus*x2-vvc*taus-
     #    10./3.*vvc*x1*x2-1./3.*vvc*x1*z+4./3.*vvc*x1+14./
     #    3.*vvc*x1s*x2-4./3.*vvc*x1s*z-3*vvc*x1s+
     #    vvc*x2-2./3.*vvc)
      wcom= wcom+yc*(vvq*tau*x1*x2*z+2./3.*
     #    vvq*tau*x1*x2-2*vvq*tau*x1*x2s+2*vvq*tau*x1*z+
     #    vvq*tau*x1*zs+vvq*tau*x1+2*vvq*tau*x2*z+2*vvq
     #    *tau*x2-2./3.*vvq*tau*x2s-2*vvq*tau*z-vvq*tau*
     #    zs-vvq*tau+vvq*taus*x2*z+vvq*taus*x2-
     #    vvq*taus*x2s-5./3.*vvq*x1*x2*z-4./3.*vvq*x1*x2
     #    +2*vvq*x1*x2s+10./3.*vvq*x1s*x2*z+3*vvq*x1s
     #    *x2-10./3.*vvq*x1s*x2s+vvq*x2*z+2./3.*vvq*x2
     #    -vvq*x2s)
      dw1= pwsm2*y*zxi*(4./9.*vvs*x1-2./9.*
     #    vvs-16./9.*tau*x1*x2-8./3.*tau*x1+8./3.*tau*x1s+8.
     #    /9.*tau*x2+4./9.*tau+8./9.*taus*x1-4./9.*taus*x2-
     #    4./9.*taus+20./9.*x1*x2+4./3.*x1-20./9.*x1s*x2-4*
     #    x1s+8./3.*x1c-4./9.*x2)
      dw1= dw1+pwsm2*y*(16./9.+4./3.*tau*x1-
     #    8./9.*tau-44./9.*x1+32./9.*x1s)
      dw1= dw1+pwtm2*y*(-2*vvs*x1+vvs*x1s+vvs)
      dw1= dw1+ys*rwms*(-16./3.*vv*tau*x1*
     #    x2+4./3.*vv*tau*x1*z+16./3.*vv*tau*x1+8./3.*vv*tau*x2
     #    -8./3.*vv*tau-4./3.*vv*taus*x2+4./3.*vv*taus+20./
     #    3.*vv*x1*x2-4./3.*vv*x1*z-20./3.*vv*x1-20./3.*vv*x1s*
     #    x2+8./3.*vv*x1s*z+20./3.*vv*x1s-4./3.*vv*x2+4./3.
     #    *vv+2*vvc*tau*x1-3*vvc*tau-2*vvc*x1*x2+2*vvc*
     #    x1*z+5./3.*vvc*x1-4./3.*vvc*x1s+2*vvc*x2-2*
     #    vvc*z-2*vvc)
      dw1= dw1+yc*(-16./3.*vvs*tau*x1*x2*z
     #    -37./3.*vvs*tau*x1*x2+28./3.*vvs*tau*x1*x2s+vvs
     #    *tau*x1*z+5./2.*vvs*tau*x1+11./3.*vvs*tau*x2-8./3.*
     #    vvs*tau*x2s-1./2.*vvs*tau-2*vvs*taus*x2*z-13./
     #    3.*vvs*taus*x2+10./3.*vvs*taus*x2s+vvs*taus
     #    +4./3.*vvs*x1*x2*z+23./3.*vvs*x1*x2-20./3.*vvs*x1
     #    *x2s-1./2.*vvs*x1-14./3.*vvs*x1s*x2*z-32./3.*
     #    vvs*x1s*x2+26./3.*vvs*x1s*x2s+vvs*x1s*z+3.
     #    /2.*vvs*x1s-4./3.*vvs*x2+4./3.*vvs*x2s+3*
     #    vvq*tau*x1*z+5./2.*vvq*tau*x1+3*vvq*tau*x2-2*
     #    vvq*tau*z-5./2.*vvq*tau+vvq*taus*x2-4*vvq*x1*
     #    x2*z-11./3.*vvq*x1*x2+2*vvq*x1*x2s+3*vvq*x1*z
     #    +2*vvq*x1*zs+3./2.*vvq*x1+10./3.*vvq*x1s*x2
     #    -1./2.*vvq*x1s+4*vvq*x2*z+3*vvq*x2-2*vvq*
     #    x2s-2*vvq*z-vvq*zs-vvq)
      dw2= pwsm2*y*zxi*(4./9.*vvs*x1-2./9.*
     #    vvs-16./3.*tau*x1*x2-8*tau*x1+8*tau*x1s+8./3.*tau
     #    *x2+4./3.*tau+8./3.*taus*x1-4./3.*taus*x2-4./3.*
     #    taus+20./3.*x1*x2+4*x1-20./3.*x1s*x2-12*x1s+8
     #    *x1c-4./3.*x2)
      dw2= dw2+pwsm2*y*(16./3.+4*tau*x1-8./3.
     #    *tau-44./3.*x1+32./3.*x1s)
      dw2= dw2+pwtm2*y*(-2*vvs*x1+vvs*x1s+vvs)
      dw2= dw2+ys*rwms*(-16*vv*tau*x1*x2
     #    +4*vv*tau*x1*z+16*vv*tau*x1+8*vv*tau*x2-8*vv*tau-4*
     #    vv*taus*x2+4*vv*taus+20*vv*x1*x2-4*vv*x1*z-20*vv*
     #    x1-20*vv*x1s*x2+8*vv*x1s*z+20*vv*x1s-4*vv*x2+
     #    4*vv+2*vvc*tau*x1-3*vvc*tau-2*vvc*x1*x2+2*vvc
     #    *x1*z+5./3.*vvc*x1-4./3.*vvc*x1s+2*vvc*x2-2*
     #    vvc*z-2*vvc)
      dw2= dw2+yc*(-16*vvs*tau*x1*x2*z-
     #    37*vvs*tau*x1*x2+28*vvs*tau*x1*x2s+3*vvs*tau*x1*z
     #    +15./2.*vvs*tau*x1+11*vvs*tau*x2-8*vvs*tau*x2s
     #    -3./2.*vvs*tau-6*vvs*taus*x2*z-13*vvs*taus*x2
     #    +10*vvs*taus*x2s+3*vvs*taus+4*vvs*x1*x2*z
     #    +23*vvs*x1*x2-20*vvs*x1*x2s-3./2.*vvs*x1-14*
     #    vvs*x1s*x2*z-32*vvs*x1s*x2+26*vvs*x1s*x2s
     #    +3*vvs*x1s*z+9./2.*vvs*x1s-4*vvs*x2+4*vvs
     #    *x2s+3*vvq*tau*x1*z+5./2.*vvq*tau*x1+3*vvq*tau*
     #    x2-2*vvq*tau*z-5./2.*vvq*tau+vvq*taus*x2-4*
     #    vvq*x1*x2*z-11./3.*vvq*x1*x2+2*vvq*x1*x2s+3*
     #    vvq*x1*z+2*vvq*x1*zs+3./2.*vvq*x1+10./3.*vvq*
     #    x1s*x2-1./2.*vvq*x1s+4*vvq*x2*z+3*vvq*x2-2*
     #    vvq*x2s-2*vvq*z-vvq*zs-vvq)
*
       else if(jj.eq.2) then
*
      wcom= pwsm2*y*zxi*(-4./9.*vvs*x1*x2+4./9.*
     #    vvs*x1s*x2+2./9.*vvs*x2)
      wcom= wcom+pwsm2*y*(-2./9.*vvs*x1-4./9.*
     #    vvs*x1s+4./9.*vvs)
      wcom= wcom+y*rwms*rwmpgs*(-4./3.*vvs*x1
     #    +8./3.*vvs*x1s-4./3.*vvs*x1c)
      wcom= wcom+ys*rwms*(4./3.*vvc*x1*x2-2./
     #    3.*vvc*x1*z-2./3.*vvc*x1-4*vvc*x1s*x2+4./3.*
     #    vvc*x1s*z+4./3.*vvc*x1s+8./3.*vvc*x1c*x2-4./
     #    3.*vvc*x1c*z-4./3.*vvc*x1c+2./3.*vvc*x2)
      wcom= wcom+yc*(2./3.*vvq*x1*x2*z+2./3.*
     #    vvq*x1*x2-4./3.*vvq*x1s*x2*z-4./3.*vvq*x1s*x2
     #    +4./3.*vvq*x1s*x2s+4./3.*vvq*x1c*x2*z+4./3.*
     #    vvq*x1c*x2-4./3.*vvq*x1c*x2s-2./3.*vvq*x2s)
      dw1= pwsm2*y*zxi*(-4./9.*vvs*x1*x2-8./9.
     #    *vvs*x1+8./9.*vvs*x1s+4./9.*vvs*x2+2./9.*vvs
     #    -4./9.*x1*x2+4./3.*x1s*x2+8./9.*x1s-8./9.*x1c*
     #    x2-16./9.*x1c+8./9.*x1q)
      dw1= dw1+pwsm2*y*(-8./9.+40./9.*x1-20./
     #    3.*x1s+28./9.*x1c)
      dw1= dw1+ys*rwms*(-4./3.*vv*x1*x2+4./
     #    3.*vv*x1+4*vv*x1s*x2-4./3.*vv*x1s*z-4*vv*x1s-8./
     #    3.*vv*x1c*x2+4./3.*vv*x1c*z+8./3.*vv*x1c-4./3.*
     #    vvc*x1*x2+4./3.*vvc*x1s-4./3.*vvc*x1c+4./3.*
     #    vvc*x2-2./3.*vvc)
      dw1= dw1+yc*(-4./3.*vvs*x1*x2+4./3.*
     #    vvs*x1*x2s+4./3.*vvs*x1s*x2*z+4*vvs*x1s*x2-
     #    4*vvs*x1s*x2s-4./3.*vvs*x1c*x2*z-8./3.*vvs*
     #    x1c*x2+8./3.*vvs*x1c*x2s+4./3.*vvq*x1*x2s-4.
     #    /3.*vvq*x1s*x2+4./3.*vvq*x1c*x2+2./3.*vvq*x2-
     #    4./3.*vvq*x2s)
      dw2= pwsm2*y*zxi*(-4./9.*vvs*x1*x2-8./9.
     #    *vvs*x1+8./9.*vvs*x1s+4./9.*vvs*x2+2./9.*vvs
     #    -4./3.*x1*x2+4*x1s*x2+8./3.*x1s-8./3.*x1c*x2-
     #    16./3.*x1c+8./3.*x1q)
      dw2= dw2+pwsm2*y*(-8./3.+40./3.*x1-20*
     #    x1s+28./3.*x1c)
      dw2= dw2+ys*rwms*(-4*vv*x1*x2+4*vv*
     #    x1+12*vv*x1s*x2-4*vv*x1s*z-12*vv*x1s-8*vv*x1c
     #    *x2+4*vv*x1c*z+8*vv*x1c-4./3.*vvc*x1*x2+4./3.*
     #    vvc*x1s-4./3.*vvc*x1c+4./3.*vvc*x2-2./3.*
     #    vvc)
      dw2= dw2+yc*(-4*vvs*x1*x2+4*vvs*x1
     #    *x2s+4*vvs*x1s*x2*z+12*vvs*x1s*x2-12*vvs*
     #    x1s*x2s-4*vvs*x1c*x2*z-8*vvs*x1c*x2+8*vvs
     #    *x1c*x2s+4./3.*vvq*x1*x2s-4./3.*vvq*x1s*x2+
     #    4./3.*vvq*x1c*x2+2./3.*vvq*x2-4./3.*vvq*x2s)
*
       else if(jj.eq.3) then
*
      wcom= pwsm2*rums*(16./9.*vv*x1-8./9.*vv*x1s-8./9.*vv)
      dw1= pwsm2*y*(-8./9.*vvs*x1+4./9.*vvs*
     #    x1s+4./9.*vvs-8./9.*x1+8./3.*x1s-8./3.*x1c+8./9.*x1q)
      dw1= dw1+pwsm2*rums*(-8./9.*vvi*x1
     #    +16./9.*vvi*x1s-8./9.*vvi*x1c-4./9.*vv*x1+4./9.*vv)
      dw2= pwsm2*y*(-8./9.*vvs*x1+4./9.*vvs*
     #    x1s+4./9.*vvs-8./3.*x1+8*x1s-8*x1c+8./3.*x1q)
      dw2= dw2+pwsm2*rums*(-8./3.*vvi*x1
     #    +16./3.*vvi*x1s-8./3.*vvi*x1c-4./9.*vv*x1+4./9.*vv)
*
       else if(jj.eq.4) then
*
      wcom= pwsm2*y*zxi*(2./9.*vvs*ups*x2-2./9.*
     #    vvs*ups-4./9.*vvs*x1*x2+4./9.*vvs*x1-2./9.*vvs)
      wcom= wcom+pwsm2*y*(2./9.*vvs*x1+4./9.*
     #    vvs*x2+1./9.*vvs*z-4./9.*vvs)
      wcom= wcom+y*rwms*rwmpgs*(-1./3.*vvs*
     #    ups*x2+1./3.*vvs*ups+2./3.*vvs*x1*x2+1./3.*vvs*x1
     #    *z+1./3.*vvs*x1-2./3.*vvs*x1s-2./3.*vvs*x2*z-
     #    2./3.*vvs*x2+1./3.*vvs)
      wcom= wcom+ys*rwms*(1./3.*vvc*ups*x1
     #    -1./3.*vvc*ups*x2*z-2./3.*vvc*ups*x2+2./3.*vvc*
     #    ups*x2s-1./3.*vvc*ups*z-1./3.*vvc*ups-2./3.*vvc
     #    *x1*x2-4./3.*vvc*x1*x2s+2*vvc*x1*z+1./3.*vvc*x1
     #    *zs+2*vvc*x1+4./3.*vvc*x1s*x2-2./3.*vvc*x1s
     #    *z-vvc*x1s-2./3.*vvc*x2*z-2./3.*vvc*x2*zs-2.
     #    /3.*vvc*x2+4./3.*vvc*x2s*z+4./3.*vvc*x2s-4./3.
     #    *vvc*z-2./3.*vvc*zs-vvc)
      wcom= wcom+yc*(-1./3.*vvq*ups*x1*x2+
     #    1./3.*vvq*ups*x2*z+1./3.*vvq*ups*x2+1./3.*vvq*ups*
     #    x2s*z+1./3.*vvq*ups*x2s-1./3.*vvq*ups*x2c-2*
     #    vvq*x1*x2*z-1./3.*vvq*x1*x2*zs-2*vvq*x1*x2-1./3.
     #    *vvq*x1*x2s*z+1./3.*vvq*x1*x2s+2./3.*vvq*x1*
     #    x2c+2./3.*vvq*x1s*x2*z+vvq*x1s*x2-2./3.*vvq
     #    *x1s*x2s+4./3.*vvq*x2*z+2./3.*vvq*x2*zs+vvq
     #    *x2+2./3.*vvq*x2s*z+2./3.*vvq*x2s*zs+1./3.*
     #    vvq*x2s-2./3.*vvq*x2c*z-2./3.*vvq*x2c)
      dw1= pwsm2*y*zxi*(-4./9.*vvs*x1+2./9.*
     #    vvs-16./9.*ups*x1*x2-8./3.*ups*x1+8./3.*ups*x1s+8.
     #    /9.*ups*x2+4./9.*ups-8./9.*upss*x1+4./9.*upss*x2+
     #    4./9.*upss-20./9.*x1*x2-4./3.*x1+20./9.*x1s*x2+4*
     #    x1s-8./3.*x1c+4./9.*x2)
      dw1= dw1+pwsm2*y*(20./9.+1./9.*vvs-
     #    22./9.*ups*x1+14./9.*ups*x2+4./3.*ups-40./9.*x1*x2-20.
     #    /9.*x1*z-22./3.*x1+46./9.*x1s+16./9.*x2*z+34./9.*x2
     #    +2./9.*x2s+4./3.*z)
      dw1= dw1+y*rwms*rwmpgs*(2./3.*vvs*x1
     #    -4./3.*vvs*x2+1./3.*vvs)
      dw1= dw1+ys*rwms*(8./3.*vv*ups*x1*x2
     #    -2./3.*vv*ups*x1*z-8./3.*vv*ups*x1-2./3.*vv*ups*x2*z-
     #    vv*ups*x2-4./3.*vv*ups*x2s+2*vv*ups*z+7./3.*vv*ups-
     #    2./3.*vv*upss*x2+2./3.*vv*upss+8./3.*vv*x1*x2*z+3*
     #    vv*x1*x2+8./3.*vv*x1*x2s-19./3.*vv*x1*z-4./3.*vv*x1*
     #    zs-6*vv*x1-10./3.*vv*x1s*x2+4./3.*vv*x1s*z+10./
     #    3.*vv*x1s+2./3.*vv*x2*z+2./3.*vv*x2*zs+1./3.*vv*x2
     #    -8./3.*vv*x2s*z-8./3.*vv*x2s+14./3.*vv*z+2*vv*
     #    zs+8./3.*vv-1./3.*vvc*ups*x2-1./3.*vvc*ups-2./3.
     #    *vvc*x1*x2+vvc*x1*z+8./3.*vvc*x1-2./3.*vvc*
     #    x1s-2*vvc*x2*z-4./3.*vvc*x2+8./3.*vvc*x2s-5.
     #    /3.*vvc*z-2*vvc)
      dw1= dw1+yc*(2./3.*vvs*ups*x1*x2*z+8.
     #    /3.*vvs*ups*x1*x2-8./3.*vvs*ups*x1*x2s-2*vvs*ups*
     #    x2*z-7./3.*vvs*ups*x2+2./3.*vvs*ups*x2s*z+vvs*
     #    ups*x2s+4./3.*vvs*ups*x2c-2./3.*vvs*upss*x2+2.
     #    /3.*vvs*upss*x2s+19./3.*vvs*x1*x2*z+4./3.*vvs*
     #    x1*x2*zs+6*vvs*x1*x2-8./3.*vvs*x1*x2s*z-3*vvs
     #    *x1*x2s-8./3.*vvs*x1*x2c-4./3.*vvs*x1s*x2*z-
     #    10./3.*vvs*x1s*x2+10./3.*vvs*x1s*x2s-14./3.*
     #    vvs*x2*z-2*vvs*x2*zs-8./3.*vvs*x2-2./3.*vvs*
     #    x2s*z-2./3.*vvs*x2s*zs-1./3.*vvs*x2s+8./3.*
     #    vvs*x2c*z+8./3.*vvs*x2c+1./3.*vvq*ups*x2+1./3.
     #    *vvq*ups*x2s-vvq*x1*x2*z-8./3.*vvq*x1*x2+2./3.*
     #    vvq*x1s*x2+5./3.*vvq*x2*z+2*vvq*x2+2*vvq*
     #    x2s*z+vvq*x2s-4./3.*vvq*x2c)
      dw2= pwsm2*y*zxi*(-4./9.*vvs*x1+2./9.*
     #    vvs-16./3.*ups*x1*x2-8*ups*x1+8*ups*x1s+8./3.*ups
     #    *x2+4./3.*ups-8./3.*upss*x1+4./3.*upss*x2+4./3.*
     #    upss-20./3.*x1*x2-4*x1+20./3.*x1s*x2+12*x1s-8
     #    *x1c+4./3.*x2)
      dw2= dw2+pwsm2*y*(20./3.+1./9.*vvs-
     #    22./3.*ups*x1+14./3.*ups*x2+4*ups-40./3.*x1*x2-20./3.
     #    *x1*z-22*x1+46./3.*x1s+16./3.*x2*z+34./3.*x2+2./3.
     #    *x2s+4*z)
      dw2= dw2+y*rwms*rwmpgs*(2./3.*vvs*x1
     #    -4./3.*vvs*x2+1./3.*vvs)
      dw2= dw2+ys*rwms*(8*vv*ups*x1*x2-2*
     #    vv*ups*x1*z-8*vv*ups*x1-2*vv*ups*x2*z-3*vv*ups*x2-4*
     #    vv*ups*x2s+6*vv*ups*z+7*vv*ups-2*vv*upss*x2+2*vv*
     #    upss+8*vv*x1*x2*z+9*vv*x1*x2+8*vv*x1*x2s-19*vv*x1
     #    *z-4*vv*x1*zs-18*vv*x1-10*vv*x1s*x2+4*vv*x1s*z
     #    +10*vv*x1s+2*vv*x2*z+2*vv*x2*zs+vv*x2-8*vv*
     #    x2s*z-8*vv*x2s+14*vv*z+6*vv*zs+8*vv-1./3.*
     #    vvc*ups*x2-1./3.*vvc*ups-2./3.*vvc*x1*x2+vvc*x1
     #    *z+8./3.*vvc*x1-2./3.*vvc*x1s-2*vvc*x2*z-4./3.
     #    *vvc*x2+8./3.*vvc*x2s-5./3.*vvc*z-2*vvc)
      dw2= dw2+yc*(2*vvs*ups*x1*x2*z+8*
     #    vvs*ups*x1*x2-8*vvs*ups*x1*x2s-6*vvs*ups*x2*z-7
     #    *vvs*ups*x2+2*vvs*ups*x2s*z+3*vvs*ups*x2s+4*
     #    vvs*ups*x2c-2*vvs*upss*x2+2*vvs*upss*x2s+
     #    19*vvs*x1*x2*z+4*vvs*x1*x2*zs+18*vvs*x1*x2-8*
     #    vvs*x1*x2s*z-9*vvs*x1*x2s-8*vvs*x1*x2c-4*
     #    vvs*x1s*x2*z-10*vvs*x1s*x2+10*vvs*x1s*x2s
     #    -14*vvs*x2*z-6*vvs*x2*zs-8*vvs*x2-2*vvs*
     #    x2s*z-2*vvs*x2s*zs-vvs*x2s+8*vvs*x2c*z
     #    +8*vvs*x2c+1./3.*vvq*ups*x2+1./3.*vvq*ups*x2s
     #    -vvq*x1*x2*z-8./3.*vvq*x1*x2+2./3.*vvq*x1s*x2
     #    +5./3.*vvq*x2*z+2*vvq*x2+2*vvq*x2s*z+vvq*
     #    x2s-4./3.*vvq*x2c)
*
       else if(jj.eq.5) then
*
      wcom= pwsm2*y*zxi*(-4./9.*vvs*x1*x2+4./9.*
     #    vvs*x1s*x2+2./9.*vvs*x2)
      wcom= wcom+pwsm2*y*(-4./9.*vvs*x1*x2+1./
     #    9.*vvs*x1*z+8./9.*vvs*x1-4./9.*vvs*x1s+2./9.*
     #    vvs*x2*z+2./9.*vvs*x2-1./9.*vvs*z+1./9.*vvs*
     #    zs-4./9.*vvs)
      wcom= wcom+y*rwms*rwmpgs*(2./3.*vvs*x1*x2*
     #    z+4./3.*vvs*x1*x2+vvs*x1*z+1./3.*vvs*x1*zs+2./
     #    3.*vvs*x1-2./3.*vvs*x1s*x2-2./3.*vvs*x1s*z-4./
     #    3.*vvs*x1s+2./3.*vvs*x1c-2./3.*vvs*x2*z-1./3.
     #    *vvs*x2*zs-2./3.*vvs*x2-1./3.*vvs*z-1./3.*vvs*zs)
      wcom= wcom+ys*rwms*(-2./3.*vvc*x1*x2*z
     #    -4./3.*vvc*x1*x2-4./3.*vvc*x1*x2s*z-8./3.*vvc*
     #    x1*x2s+10./3.*vvc*x1*z+5./3.*vvc*x1*zs+1./3.*
     #    vvc*x1*zc+2*vvc*x1+2./3.*vvc*x1s*x2*z+8./3.*
     #    vvc*x1s*x2+4./3.*vvc*x1s*x2s-8./3.*vvc*x1s*
     #    z-2./3.*vvc*x1s*zs-2*vvc*x1s-4./3.*vvc*
     #    x1c*x2+2./3.*vvc*x1c*z+2./3.*vvc*x1c-1./3.*
     #    vvc*x2*zc+4./3.*vvc*x2s*z+2./3.*vvc*x2s*zs
     #    +4./3.*vvc*x2s-4./3.*vvc*z-vvc*zs-1./3.*
     #    vvc*zc-2./3.*vvc)
      wcom= wcom+yc*(-10./3.*vvq*x1*x2*z-5./
     #    3.*vvq*x1*x2*zs-1./3.*vvq*x1*x2*zc-2*vvq*x1*x2
     #    -1./3.*vvq*x1*x2s*z-1./3.*vvq*x1*x2s*zs+2./3.
     #    *vvq*x1*x2s+2./3.*vvq*x1*x2c*z+4./3.*vvq*x1*
     #    x2c+8./3.*vvq*x1s*x2*z+2./3.*vvq*x1s*x2*zs+
     #    2*vvq*x1s*x2-4./3.*vvq*x1s*x2s-2./3.*vvq*
     #    x1s*x2c-2./3.*vvq*x1c*x2*z-2./3.*vvq*x1c*x2
     #    +2./3.*vvq*x1c*x2s+4./3.*vvq*x2*z+vvq*x2*zs
     #    +1./3.*vvq*x2*zc+2./3.*vvq*x2+1./3.*vvq*x2s*z
     #    +1./3.*vvq*x2s*zs+1./3.*vvq*x2s*zc-2./3.*
     #    vvq*x2c*z-1./3.*vvq*x2c*zs-2./3.*vvq*x2c)
      dw1= pwsm2*y*zxi*(-4./9.*vvs*x1*x2-8./9.
     #    *vvs*x1+8./9.*vvs*x1s+4./9.*vvs*x2+2./9.*vvs
     #    -4./9.*x1*x2+4./3.*x1s*x2+8./9.*x1s-8./9.*x1c*
     #    x2-16./9.*x1c+8./9.*x1q)
      dw1= dw1+pwsm2*y*(8./9.-2./9.*vvs*x1+4.
     #    /9.*vvs*x2+1./3.*vvs*z-28./9.*x1*x2*z-52./9.*x1*x2
     #    -4./9.*x1*x2s-4*x1*z-2./3.*x1*zs-14./3.*x1+10./
     #    3.*x1s*x2+8./3.*x1s*z+20./3.*x1s-26./9.*x1c+8.
     #    /3.*x2*z+2./3.*x2*zs+22./9.*x2+4./9.*x2s*z+4./9.*
     #    x2s+4./3.*z+4./9.*zs)
      dw1= dw1+y*rwms*rwmpgs*(2*vvs*x1*x2+2./
     #    3.*vvs*x1*z+4./3.*vvs*x1-4./3.*vvs*x1s-4./3.*
     #    vvs*x2*z-2*vvs*x2-1./3.*vvs*z)
      dw1= dw1+ys*rwms*(1./3.*vv*x1*x2*z+8./
     #    3.*vv*x1*x2s*z+8./3.*vv*x1*x2s-19./3.*vv*x1*z-11./3.
     #    *vv*x1*zs-2./3.*vv*x1*zc-10./3.*vv*x1-2*vv*x1s*x2
     #    *z-2*vv*x1s*x2-4./3.*vv*x1s*x2s+13./3.*vv*x1s*z
     #    +4./3.*vv*x1s*zs+11./3.*vv*x1s+4./3.*vv*x1c*x2
     #    -2./3.*vv*x1c*z-4./3.*vv*x1c+5./3.*vv*x2*z+5./3.*
     #    vv*x2*zs+2./3.*vv*x2*zc+2./3.*vv*x2-8./3.*vv*x2s*
     #    z-4./3.*vv*x2s*zs-4./3.*vv*x2s+8./3.*vv*z+7./3.
     #    *vv*zs+2./3.*vv*zc+vv+4./3.*vvc*x1*x2*z-4./3.*
     #    vvc*x1*x2-4*vvc*x1*x2s+5*vvc*x1*z+vvc*x1*zs
     #    +14./3.*vvc*x1+2*vvc*x1s*x2-2*vvc*x1s*z-10./
     #    3.*vvc*x1s+2./3.*vvc*x1c-2*vvc*x2*z-5./3.*
     #    vvc*x2*zs-2./3.*vvc*x2+8./3.*vvc*x2s*z+4*
     #    vvc*x2s-3*vvc*z-4./3.*vvc*zs-2*vvc)
      dw1= dw1+yc*(19./3.*vvs*x1*x2*z+11./3.
     #    *vvs*x1*x2*zs+2./3.*vvs*x1*x2*zc+10./3.*vvs*x1*
     #    x2-1./3.*vvs*x1*x2s*z-8./3.*vvs*x1*x2c*z-8./3.*
     #    vvs*x1*x2c-13./3.*vvs*x1s*x2*z-4./3.*vvs*x1s*
     #    x2*zs-11./3.*vvs*x1s*x2+2*vvs*x1s*x2s*z+2*
     #    vvs*x1s*x2s+4./3.*vvs*x1s*x2c+2./3.*vvs*
     #    x1c*x2*z+4./3.*vvs*x1c*x2-4./3.*vvs*x1c*x2s
     #    -8./3.*vvs*x2*z-7./3.*vvs*x2*zs-2./3.*vvs*x2*
     #    zc-vvs*x2-5./3.*vvs*x2s*z-5./3.*vvs*x2s*
     #    zs-2./3.*vvs*x2s*zc-2./3.*vvs*x2s+8./3.*
     #    vvs*x2c*z+4./3.*vvs*x2c*zs+4./3.*vvs*x2c-
     #    5*vvq*x1*x2*z-vvq*x1*x2*zs-14./3.*vvq*x1*x2-2*
     #    vvq*x1*x2s*z+2*vvq*x1*x2c+2*vvq*x1s*x2*z+10.
     #    /3.*vvq*x1s*x2-2./3.*vvq*x1s*x2s-2./3.*vvq*
     #    x1c*x2+3*vvq*x2*z+4./3.*vvq*x2*zs+2*vvq*x2+
     #    7./3.*vvq*x2s*z)
      dw1= dw1+yc*(5./3.*vvq*x2s*zs+2./3.
     #    *vvq*x2s-4./3.*vvq*x2c*z-2*vvq*x2c)
      dw2= pwsm2*y*zxi*(-4./9.*vvs*x1*x2-8./9.
     #    *vvs*x1+8./9.*vvs*x1s+4./9.*vvs*x2+2./9.*vvs
     #    -4./3.*x1*x2+4*x1s*x2+8./3.*x1s-8./3.*x1c*x2-
     #    16./3.*x1c+8./3.*x1q)
      dw2= dw2+pwsm2*y*(8./3.-2./9.*vvs*x1+4.
     #    /9.*vvs*x2+1./3.*vvs*z-28./3.*x1*x2*z-52./3.*x1*x2
     #    -4./3.*x1*x2s-12*x1*z-2*x1*zs-14*x1+10*x1s*x2
     #    +8*x1s*z+20*x1s-26./3.*x1c+8*x2*z+2*x2*zs
     #    +22./3.*x2+4./3.*x2s*z+4./3.*x2s+4*z+4./3.*zs)
      dw2= dw2+y*rwms*rwmpgs*(2*vvs*x1*x2+2./
     #    3.*vvs*x1*z+4./3.*vvs*x1-4./3.*vvs*x1s-4./3.*
     #    vvs*x2*z-2*vvs*x2-1./3.*vvs*z)
      dw2= dw2+ys*rwms*(vv*x1*x2*z+8*vv*x1*
     #    x2s*z+8*vv*x1*x2s-19*vv*x1*z-11*vv*x1*zs-2*vv*
     #    x1*zc-10*vv*x1-6*vv*x1s*x2*z-6*vv*x1s*x2-4*vv*
     #    x1s*x2s+13*vv*x1s*z+4*vv*x1s*zs+11*vv*x1s
     #    +4*vv*x1c*x2-2*vv*x1c*z-4*vv*x1c+5*vv*x2*z+5*
     #    vv*x2*zs+2*vv*x2*zc+2*vv*x2-8*vv*x2s*z-4*vv*
     #    x2s*zs-4*vv*x2s+8*vv*z+7*vv*zs+2*vv*zc+3*
     #    vv+4./3.*vvc*x1*x2*z-4./3.*vvc*x1*x2-4*vvc*x1*
     #    x2s+5*vvc*x1*z+vvc*x1*zs+14./3.*vvc*x1+2*
     #    vvc*x1s*x2-2*vvc*x1s*z-10./3.*vvc*x1s+2./3.
     #    *vvc*x1c-2*vvc*x2*z-5./3.*vvc*x2*zs-2./3.*
     #    vvc*x2+8./3.*vvc*x2s*z+4*vvc*x2s-3*vvc*z-
     #    4./3.*vvc*zs-2*vvc)
      dw2= dw2+yc*(19*vvs*x1*x2*z+11*vvs*
     #    x1*x2*zs+2*vvs*x1*x2*zc+10*vvs*x1*x2-vvs*x1*
     #    x2s*z-8*vvs*x1*x2c*z-8*vvs*x1*x2c-13*vvs*
     #    x1s*x2*z-4*vvs*x1s*x2*zs-11*vvs*x1s*x2+6*
     #    vvs*x1s*x2s*z+6*vvs*x1s*x2s+4*vvs*x1s*
     #    x2c+2*vvs*x1c*x2*z+4*vvs*x1c*x2-4*vvs*x1c
     #    *x2s-8*vvs*x2*z-7*vvs*x2*zs-2*vvs*x2*zc-3
     #    *vvs*x2-5*vvs*x2s*z-5*vvs*x2s*zs-2*vvs*
     #    x2s*zc-2*vvs*x2s+8*vvs*x2c*z+4*vvs*x2c*
     #    zs+4*vvs*x2c-5*vvq*x1*x2*z-vvq*x1*x2*zs-
     #    14./3.*vvq*x1*x2-2*vvq*x1*x2s*z+2*vvq*x1*x2c+
     #    2*vvq*x1s*x2*z+10./3.*vvq*x1s*x2-2./3.*vvq*
     #    x1s*x2s-2./3.*vvq*x1c*x2+3*vvq*x2*z+4./3.*
     #    vvq*x2*zs+2*vvq*x2+7./3.*vvq*x2s*z+5./3.*
     #    vvq*x2s*zs+2./3.*vvq*x2s-4./3.*vvq*x2c*z-
     #    2*vvq*x2c)
*
       else if(jj.eq.6) then
*
      wcom= pwsm2*rdms*(4./9.*vv*x1*z+4./9.*vv*x1-2./9.
     #    *vv*x1s-4./9.*vv*z-2./9.*vv*zs-2./9.*vv)
      dw1= pwsm2*y*(-2./9.*vvs*x1*z-2./9.*vvs*x1
     #    +1./9.*vvs*x1s+2./9.*vvs*z+1./9.*vvs*zs+1./
     #    9.*vvs-4./3.*x1*x2*z-4./9.*x1*x2*zs-8./9.*x1*x2-4.
     #    /9.*x1*x2s*z-4./9.*x1*x2s-4./9.*x1*z-2./9.*x1*zs
     #    -2./9.*x1+8./9.*x1s*x2*z+10./9.*x1s*x2+2./9.*
     #    x1s*x2s+8./9.*x1s*z+2./9.*x1s*zs+2./3.*x1s
     #    -4./9.*x1c*x2-4./9.*x1c*z-2./3.*x1c+2./9.*x1q
     #    +4./9.*x2*z+2./9.*x2*zs+2./9.*x2+4./9.*x2s*z+2./
     #    9.*x2s*zs+2./9.*x2s)
      dw1= dw1+pwsm2*rdms*(-4./3.*vvi*x1*
     #    x2*z-14./9.*vvi*x1*x2-2./9.*vvi*x1*x2s-14./
     #    9.*vvi*x1*z-4./9.*vvi*x1*zs-10./9.*vvi*
     #    x1+8./9.*vvi*x1s*x2+10./9.*vvi*x1s*z+14./
     #    9.*vvi*x1s-2./3.*vvi*x1c+10./9.*vvi*x2
     #    *z+4./9.*vvi*x2*zs+2./3.*vvi*x2+2./9.*
     #    vvi*x2s*z+2./9.*vvi*x2s+4./9.*vvi*z+
     #    2./9.*vvi*zs+2./9.*vvi+1./3.*vv*x1-1./3.*vv
     #    *z-1./3.*vv)
      dw2= pwsm2*y*(-2./9.*vvs*x1*z-2./9.*vvs*x1
     #    +1./9.*vvs*x1s+2./9.*vvs*z+1./9.*vvs*zs+1./
     #    9.*vvs-4*x1*x2*z-4./3.*x1*x2*zs-8./3.*x1*x2-4./3.
     #    *x1*x2s*z-4./3.*x1*x2s-4./3.*x1*z-2./3.*x1*zs-2.
     #    /3.*x1+8./3.*x1s*x2*z+10./3.*x1s*x2+2./3.*x1s*
     #    x2s+8./3.*x1s*z+2./3.*x1s*zs+2*x1s-4./3.*
     #    x1c*x2-4./3.*x1c*z-2*x1c+2./3.*x1q+4./3.*x2*z
     #    +2./3.*x2*zs+2./3.*x2+4./3.*x2s*z+2./3.*x2s*
     #    zs+2./3.*x2s)
      dw2= dw2+pwsm2*rdms*(-4*vvi*x1*x2*z
     #    -14./3.*vvi*x1*x2-2./3.*vvi*x1*x2s-14./3.*
     #    vvi*x1*z-4./3.*vvi*x1*zs-10./3.*vvi*x1
     #    +8./3.*vvi*x1s*x2+10./3.*vvi*x1s*z+14./3.
     #    *vvi*x1s-2*vvi*x1c+10./3.*vvi*x2*z+4.
     #    /3.*vvi*x2*zs+2*vvi*x2+2./3.*vvi*x2s*z
     #    +2./3.*vvi*x2s+4./3.*vvi*z+2./3.*vvi*
     #    zs+2./3.*vvi+1./3.*vv*x1-1./3.*vv*z-1./3.*vv)
*
       endif
*
       wcom= wcom*(2.d0-y)**2/ys-dw1+2.d0/y*(1.d0-1.d0/y)*dw2
*
      else if((ipro.eq.1.or.ipro.eq.3).and.ii.le.2) then
       if(jj.eq.1) then
*
      wcom= pwtm2*y*(-vvs*tau*x1+vvs*tau)
      wcom= wcom+ys*rwms*(2*vvc*tau*x1*x2+
     #    vvc*tau*x1*z-2*vvc*tau*z-2*vvc*tau+vvc*taus*
     #    x2-vvc*taus-2*vvc*x1*x2-vvc*x1*z+2*vvc*
     #    x1s*x2-vvc*x1s+vvc*x2)
      wcom= wcom+yc*(vvq*tau*x1*x2*z-2*vvq
     #    *tau*x1*x2s+2*vvq*tau*x1*z+vvq*tau*x1*zs+vvq*
     #    tau*x1+2*vvq*tau*x2*z+2*vvq*tau*x2-2*vvq*tau*z-
     #    vvq*tau*zs-vvq*tau+vvq*taus*x2*z+vvq*taus
     #    *x2-vvq*taus*x2s-vvq*x1*x2*z+2*vvq*x1*x2s
     #    +2*vvq*x1s*x2*z+vvq*x1s*x2-2*vvq*x1s*x2s
     #    +vvq*x2*z-vvq*x2s)
*
       else if(jj.eq.2) then
        wcom= 0.d0
*
       else if(jj.eq.3) then
        wcom= 0.d0
*
       else if(jj.eq.4) then
*
      wcom= pwsm2*y*(vvs*z)
      wcom= wcom+y*rwms*rwmpgs*(-vvs*ups*x2
     #    +vvs*ups+2*vvs*x1*x2+vvs*x1*z+vvs*x1-2*
     #    vvs*x1s-2*vvs*x2*z-2*vvs*x2+vvs)
      wcom= wcom+ys*rwms*(vvc*ups*x1-vvc
     #    *ups*x2*z-2*vvc*ups*x2+2*vvc*ups*x2s-vvc*ups*z
     #    -vvc*ups-2*vvc*x1*x2-4*vvc*x1*x2s+6*vvc*x1*
     #    z+vvc*x1*zs+6*vvc*x1+4*vvc*x1s*x2-2*vvc*
     #    x1s*z-3*vvc*x1s-2*vvc*x2*z-2*vvc*x2*zs-2*
     #    vvc*x2+4*vvc*x2s*z+4*vvc*x2s-4*vvc*z-2*
     #    vvc*zs-3*vvc)
      wcom= wcom+yc*(-vvq*ups*x1*x2+vvq*
     #    ups*x2*z+vvq*ups*x2+vvq*ups*x2s*z+vvq*ups*x2s
     #    -vvq*ups*x2c-6*vvq*x1*x2*z-vvq*x1*x2*zs-6*
     #    vvq*x1*x2-vvq*x1*x2s*z+vvq*x1*x2s+2*vvq*x1*
     #    x2c+2*vvq*x1s*x2*z+3*vvq*x1s*x2-2*vvq*x1s
     #    *x2s+4*vvq*x2*z+2*vvq*x2*zs+3*vvq*x2+2*
     #    vvq*x2s*z+2*vvq*x2s*zs+vvq*x2s-2*vvq*
     #    x2c*z-2*vvq*x2c)
*
       else if(jj.eq.5) then
*
      wcom= pwsm2*y*(-vvs*x1*z+vvs*z+vvs*zs)
      wcom= wcom+y*rwms*rwmpgs*(2*vvs*x1*x2*z+
     #    4*vvs*x1*x2+3*vvs*x1*z+vvs*x1*zs+2*vvs*x1-2
     #    *vvs*x1s*x2-2*vvs*x1s*z-4*vvs*x1s+2*vvs*
     #    x1c-2*vvs*x2*z-vvs*x2*zs-2*vvs*x2-vvs*z
     #    -vvs*zs)
      wcom= wcom+ys*rwms*(-2*vvc*x1*x2*z-4
     #    *vvc*x1*x2-4*vvc*x1*x2s*z-8*vvc*x1*x2s+10*
     #    vvc*x1*z+5*vvc*x1*zs+vvc*x1*zc+6*vvc*x1+2
     #    *vvc*x1s*x2*z+8*vvc*x1s*x2+4*vvc*x1s*x2s-
     #    8*vvc*x1s*z-2*vvc*x1s*zs-6*vvc*x1s-4*
     #    vvc*x1c*x2+2*vvc*x1c*z+2*vvc*x1c-vvc*x2*
     #    zc+4*vvc*x2s*z+2*vvc*x2s*zs+4*vvc*x2s
     #    -4*vvc*z-3*vvc*zs-vvc*zc-2*vvc)
      wcom= wcom+yc*(-10*vvq*x1*x2*z-5*vvq
     #    *x1*x2*zs-vvq*x1*x2*zc-6*vvq*x1*x2-vvq*x1*
     #    x2s*z-vvq*x1*x2s*zs+2*vvq*x1*x2s+2*vvq*x1
     #    *x2c*z+4*vvq*x1*x2c+8*vvq*x1s*x2*z+2*vvq*
     #    x1s*x2*zs+6*vvq*x1s*x2-4*vvq*x1s*x2s-2*
     #    vvq*x1s*x2c-2*vvq*x1c*x2*z-2*vvq*x1c*x2+2
     #    *vvq*x1c*x2s+4*vvq*x2*z+3*vvq*x2*zs+vvq*
     #    x2*zc+2*vvq*x2+vvq*x2s*z+vvq*x2s*zs+
     #    vvq*x2s*zc-2*vvq*x2c*z-vvq*x2c*zs-2*
     #    vvq*x2c)
*
       else if(jj.eq.6) then
*
      wcom= pwsm2*rdms*(4*vv*x1*z+4*vv*x1-2*vv*x1s-
     #    4*vv*z-2*vv*zs-2*vv)
*
       endif
      else if((ipro.eq.1.or.ipro.eq.3).and.ii.eq.3) then
       if(jj.eq.1) then
*
      wcom= pwtm2*y*(-vvs*tau*x1+vvs*tau)
      wcom= wcom+ys*rwms*(2*vvc*tau*x1*x2+
     #    vvc*tau*x1*z-2*vvc*tau*z-2*vvc*tau+vvc*taus*
     #    x2-vvc*taus-2*vvc*x1*x2-vvc*x1*z+2*vvc*
     #    x1s*x2-vvc*x1s+vvc*x2)
      wcom= wcom+yc*(vvq*tau*x1*x2*z-2*vvq
     #    *tau*x1*x2s+2*vvq*tau*x1*z+vvq*tau*x1*zs+vvq*
     #    tau*x1+2*vvq*tau*x2*z+2*vvq*tau*x2-2*vvq*tau*z-
     #    vvq*tau*zs-vvq*tau+vvq*taus*x2*z+vvq*taus
     #    *x2-vvq*taus*x2s-vvq*x1*x2*z+2*vvq*x1*x2s
     #    +2*vvq*x1s*x2*z+vvq*x1s*x2-2*vvq*x1s*x2s
     #    +vvq*x2*z-vvq*x2s)
      dw1= pwtm2*y*(-2*vvs*x1+vvs*x1s+vvs)
      dw1= dw1+ys*rwms*(2*vvc*tau*x1-3*
     #    vvc*tau-2*vvc*x1*x2+2*vvc*x1*z+vvc*x1+2*vvc
     #    *x2-2*vvc*z-2*vvc)
      dw1= dw1+yc*(-4*vvs*tau*x1*x2*z-7*
     #    vvs*tau*x1*x2+4*vvs*tau*x1*x2s+vvs*tau*x1*z+5./
     #    2.*vvs*tau*x1+vvs*tau*x2-1./2.*vvs*tau-2*vvs*
     #    taus*x2*z-3*vvs*taus*x2+2*vvs*taus*x2s+
     #    vvs*taus+vvs*x1*x2-1./2.*vvs*x1-2*vvs*x1s*
     #    x2*z-4*vvs*x1s*x2+2*vvs*x1s*x2s+vvs*x1s*z
     #    +3./2.*vvs*x1s+3*vvq*tau*x1*z+5./2.*vvq*tau*x1
     #    +3*vvq*tau*x2-2*vvq*tau*z-5./2.*vvq*tau+vvq*
     #    taus*x2-4*vvq*x1*x2*z-3*vvq*x1*x2+2*vvq*x1*
     #    x2s+3*vvq*x1*z+2*vvq*x1*zs+3./2.*vvq*x1+2*
     #    vvq*x1s*x2-1./2.*vvq*x1s+4*vvq*x2*z+3*vvq*
     #    x2-2*vvq*x2s-2*vvq*z-vvq*zs-vvq)
      dw2= pwtm2*y*(-2*vvs*x1+vvs*x1s+vvs)
      dw2= dw2+ys*rwms*(2*vvc*tau*x1-3*
     #    vvc*tau-2*vvc*x1*x2+2*vvc*x1*z+vvc*x1+2*vvc
     #    *x2-2*vvc*z-2*vvc)
      dw2= dw2+yc*(-12*vvs*tau*x1*x2*z-
     #    21*vvs*tau*x1*x2+12*vvs*tau*x1*x2s+3*vvs*tau*x1*z
     #    +15./2.*vvs*tau*x1+3*vvs*tau*x2-3./2.*vvs*tau-6
     #    *vvs*taus*x2*z-9*vvs*taus*x2+6*vvs*taus*x2s
     #    +3*vvs*taus+3*vvs*x1*x2-3./2.*vvs*x1-6*vvs*
     #    x1s*x2*z-12*vvs*x1s*x2+6*vvs*x1s*x2s+3*
     #    vvs*x1s*z+9./2.*vvs*x1s+3*vvq*tau*x1*z+5./2.*
     #    vvq*tau*x1+3*vvq*tau*x2-2*vvq*tau*z-5./2.*vvq*
     #    tau+vvq*taus*x2-4*vvq*x1*x2*z-3*vvq*x1*x2+2*
     #    vvq*x1*x2s+3*vvq*x1*z+2*vvq*x1*zs+3./2.*vvq
     #    *x1+2*vvq*x1s*x2-1./2.*vvq*x1s+4*vvq*x2*z+3
     #    *vvq*x2-2*vvq*x2s-2*vvq*z-vvq*zs-vvq)
*
       else if(jj.eq.2) then
        wcom= 0.d0
        dw1= 0.d0
        dw2= 0.d0
*
       else if(jj.eq.3) then
        wcom= 0.d0
        dw1= 0.d0
        dw2= 0.d0
*
       else if(jj.eq.4) then
*
      wcom= pwsm2*y*(vvs*z)
      wcom= wcom+y*rwms*rwmpgs*(-vvs*ups*x2
     #    +vvs*ups+2*vvs*x1*x2+vvs*x1*z+vvs*x1-2*
     #    vvs*x1s-2*vvs*x2*z-2*vvs*x2+vvs)
      wcom= wcom+ys*rwms*(vvc*ups*x1-vvc
     #    *ups*x2*z-2*vvc*ups*x2+2*vvc*ups*x2s-vvc*ups*z
     #    -vvc*ups-2*vvc*x1*x2-4*vvc*x1*x2s+6*vvc*x1*
     #    z+vvc*x1*zs+6*vvc*x1+4*vvc*x1s*x2-2*vvc*
     #    x1s*z-3*vvc*x1s-2*vvc*x2*z-2*vvc*x2*zs-2*
     #    vvc*x2+4*vvc*x2s*z+4*vvc*x2s-4*vvc*z-2*
     #    vvc*zs-3*vvc)
      wcom= wcom+yc*(-vvq*ups*x1*x2+vvq*
     #    ups*x2*z+vvq*ups*x2+vvq*ups*x2s*z+vvq*ups*x2s
     #    -vvq*ups*x2c-6*vvq*x1*x2*z-vvq*x1*x2*zs-6*
     #    vvq*x1*x2-vvq*x1*x2s*z+vvq*x1*x2s+2*vvq*x1*
     #    x2c+2*vvq*x1s*x2*z+3*vvq*x1s*x2-2*vvq*x1s
     #    *x2s+4*vvq*x2*z+2*vvq*x2*zs+3*vvq*x2+2*
     #    vvq*x2s*z+2*vvq*x2s*zs+vvq*x2s-2*vvq*
     #    x2c*z-2*vvq*x2c)
*
      dw1= pwsm2*y*(vvs-2*ups*x1+2*ups*x2-8*x1*
     #    x2-4*x1*z-6*x1+6*x1s+4*x2*z+6*x2+2*x2s)
      dw1= dw1+y*rwms*rwmpgs*(2*vvs*x1-4*
     #    vvs*x2+vvs)
      dw1= dw1+ys*rwms*(8*vv*ups*x1*x2-2*
     #    vv*ups*x1*z-8*vv*ups*x1-2*vv*ups*x2*z-3*vv*ups*x2-4*
     #    vv*ups*x2s+6*vv*ups*z+7*vv*ups-2*vv*upss*x2+2*vv*
     #    upss+8*vv*x1*x2*z+9*vv*x1*x2+8*vv*x1*x2s-19*vv*x1
     #    *z-4*vv*x1*zs-18*vv*x1-10*vv*x1s*x2+4*vv*x1s*z
     #    +10*vv*x1s+2*vv*x2*z+2*vv*x2*zs+vv*x2-8*vv*
     #    x2s*z-8*vv*x2s+14*vv*z+6*vv*zs+8*vv-vvc*ups
     #    *x2-vvc*ups-2*vvc*x1*x2+3*vvc*x1*z+8*vvc*x1
     #    -2*vvc*x1s-6*vvc*x2*z-4*vvc*x2+8*vvc*x2s
     #    -5*vvc*z-6*vvc)
      dw1= dw1+yc*(2*vvs*ups*x1*x2*z+8*
     #    vvs*ups*x1*x2-8*vvs*ups*x1*x2s-6*vvs*ups*x2*z-7
     #    *vvs*ups*x2+2*vvs*ups*x2s*z+3*vvs*ups*x2s+4*
     #    vvs*ups*x2c-2*vvs*upss*x2+2*vvs*upss*x2s+
     #    19*vvs*x1*x2*z+4*vvs*x1*x2*zs+18*vvs*x1*x2-8*
     #    vvs*x1*x2s*z-9*vvs*x1*x2s-8*vvs*x1*x2c-4*
     #    vvs*x1s*x2*z-10*vvs*x1s*x2+10*vvs*x1s*x2s
     #    -14*vvs*x2*z-6*vvs*x2*zs-8*vvs*x2-2*vvs*
     #    x2s*z-2*vvs*x2s*zs-vvs*x2s+8*vvs*x2c*z
     #    +8*vvs*x2c+vvq*ups*x2+vvq*ups*x2s-3*vvq*
     #    x1*x2*z-8*vvq*x1*x2+2*vvq*x1s*x2+5*vvq*x2*z+6
     #    *vvq*x2+6*vvq*x2s*z+3*vvq*x2s-4*vvq*x2c)
      dw2= pwsm2*y*(vvs-6*ups*x1+6*ups*x2-24*x1*
     #    x2-12*x1*z-18*x1+18*x1s+12*x2*z+18*x2+6*x2s)
      dw2= dw2+y*rwms*rwmpgs*(2*vvs*x1-4*
     #    vvs*x2+vvs)
      dw2= dw2+ys*rwms*(24*vv*ups*x1*x2-6*
     #    vv*ups*x1*z-24*vv*ups*x1-6*vv*ups*x2*z-9*vv*ups*x2-12
     #    *vv*ups*x2s+18*vv*ups*z+21*vv*ups-6*vv*upss*x2+6*
     #    vv*upss+24*vv*x1*x2*z+27*vv*x1*x2+24*vv*x1*x2s-57
     #    *vv*x1*z-12*vv*x1*zs-54*vv*x1-30*vv*x1s*x2+12*vv*
     #    x1s*z+30*vv*x1s+6*vv*x2*z+6*vv*x2*zs+3*vv*x2-
     #    24*vv*x2s*z-24*vv*x2s+42*vv*z+18*vv*zs+24*vv-
     #    vvc*ups*x2-vvc*ups-2*vvc*x1*x2+3*vvc*x1*z+8*
     #    vvc*x1-2*vvc*x1s-6*vvc*x2*z-4*vvc*x2+8*
     #    vvc*x2s-5*vvc*z-6*vvc)
      dw2= dw2+yc*(6*vvs*ups*x1*x2*z+24*
     #    vvs*ups*x1*x2-24*vvs*ups*x1*x2s-18*vvs*ups*x2*z
     #    -21*vvs*ups*x2+6*vvs*ups*x2s*z+9*vvs*ups*x2s
     #    +12*vvs*ups*x2c-6*vvs*upss*x2+6*vvs*upss*
     #    x2s+57*vvs*x1*x2*z+12*vvs*x1*x2*zs+54*vvs*x1*
     #    x2-24*vvs*x1*x2s*z-27*vvs*x1*x2s-24*vvs*x1*
     #    x2c-12*vvs*x1s*x2*z-30*vvs*x1s*x2+30*vvs*
     #    x1s*x2s-42*vvs*x2*z-18*vvs*x2*zs-24*vvs*x2
     #    -6*vvs*x2s*z-6*vvs*x2s*zs-3*vvs*x2s+24*
     #    vvs*x2c*z+24*vvs*x2c+vvq*ups*x2+vvq*ups*
     #    x2s-3*vvq*x1*x2*z-8*vvq*x1*x2+2*vvq*x1s*x2+
     #    5*vvq*x2*z+6*vvq*x2+6*vvq*x2s*z+3*vvq*x2s
     #    -4*vvq*x2c)
*
       else if(jj.eq.5) then
*
      wcom= pwsm2*y*(-vvs*x1*z+vvs*z+vvs*zs)
      wcom= wcom+y*rwms*rwmpgs*(2*vvs*x1*x2*z+
     #    4*vvs*x1*x2+3*vvs*x1*z+vvs*x1*zs+2*vvs*x1-2
     #    *vvs*x1s*x2-2*vvs*x1s*z-4*vvs*x1s+2*vvs*
     #    x1c-2*vvs*x2*z-vvs*x2*zs-2*vvs*x2-vvs*z
     #    -vvs*zs)
      wcom= wcom+ys*rwms*(-2*vvc*x1*x2*z-4
     #    *vvc*x1*x2-4*vvc*x1*x2s*z-8*vvc*x1*x2s+10*
     #    vvc*x1*z+5*vvc*x1*zs+vvc*x1*zc+6*vvc*x1+2
     #    *vvc*x1s*x2*z+8*vvc*x1s*x2+4*vvc*x1s*x2s-
     #    8*vvc*x1s*z-2*vvc*x1s*zs-6*vvc*x1s-4*
     #    vvc*x1c*x2+2*vvc*x1c*z+2*vvc*x1c-vvc*x2*
     #    zc+4*vvc*x2s*z+2*vvc*x2s*zs+4*vvc*x2s
     #    -4*vvc*z-3*vvc*zs-vvc*zc-2*vvc)
      wcom= wcom+yc*(-10*vvq*x1*x2*z-5*vvq
     #    *x1*x2*zs-vvq*x1*x2*zc-6*vvq*x1*x2-vvq*x1*
     #    x2s*z-vvq*x1*x2s*zs+2*vvq*x1*x2s+2*vvq*x1
     #    *x2c*z+4*vvq*x1*x2c+8*vvq*x1s*x2*z+2*vvq*
     #    x1s*x2*zs+6*vvq*x1s*x2-4*vvq*x1s*x2s-2*
     #    vvq*x1s*x2c-2*vvq*x1c*x2*z-2*vvq*x1c*x2+2
     #    *vvq*x1c*x2s+4*vvq*x2*z+3*vvq*x2*zs+vvq*
     #    x2*zc+2*vvq*x2+vvq*x2s*z+vvq*x2s*zs+
     #    vvq*x2s*zc-2*vvq*x2c*z-vvq*x2c*zs-2*
     #    vvq*x2c)
      dw1= pwsm2*y*(-2*vvs*x1+3*vvs*z+2*vvs
     #    -12*x1*x2*z-16*x1*x2-4*x1*x2s-8*x1*z-2*x1*zs-
     #    6*x1+10*x1s*x2+8*x1s*z+12*x1s-6*x1c+8*x2*z
     #    +2*x2*zs+6*x2+4*x2s*z+4*x2s)
      dw1= dw1+y*rwms*rwmpgs*(6*vvs*x1*x2+2*
     #    vvs*x1*z+4*vvs*x1-4*vvs*x1s-4*vvs*x2*z-6*
     #    vvs*x2-vvs*z)
      dw1= dw1+ys*rwms*(vv*x1*x2*z+8*vv*x1*
     #    x2s*z+8*vv*x1*x2s-19*vv*x1*z-11*vv*x1*zs-2*vv*
     #    x1*zc-10*vv*x1-6*vv*x1s*x2*z-6*vv*x1s*x2-4*vv*
     #    x1s*x2s+13*vv*x1s*z+4*vv*x1s*zs+11*vv*x1s
     #    +4*vv*x1c*x2-2*vv*x1c*z-4*vv*x1c+5*vv*x2*z+5*
     #    vv*x2*zs+2*vv*x2*zc+2*vv*x2-8*vv*x2s*z-4*vv*
     #    x2s*zs-4*vv*x2s+8*vv*z+7*vv*zs+2*vv*zc+3*
     #    vv+4*vvc*x1*x2*z-4*vvc*x1*x2-12*vvc*x1*x2s+15
     #    *vvc*x1*z+3*vvc*x1*zs+14*vvc*x1+6*vvc*x1s*
     #    x2-6*vvc*x1s*z-10*vvc*x1s+2*vvc*x1c-6*
     #    vvc*x2*z-5*vvc*x2*zs-2*vvc*x2+8*vvc*x2s*z
     #    +12*vvc*x2s-9*vvc*z-4*vvc*zs-6*vvc)
      dw1= dw1+yc*(19*vvs*x1*x2*z+11*vvs*
     #    x1*x2*zs+2*vvs*x1*x2*zc+10*vvs*x1*x2-vvs*x1*
     #    x2s*z-8*vvs*x1*x2c*z-8*vvs*x1*x2c-13*vvs*
     #    x1s*x2*z-4*vvs*x1s*x2*zs-11*vvs*x1s*x2+6*
     #    vvs*x1s*x2s*z+6*vvs*x1s*x2s+4*vvs*x1s*
     #    x2c+2*vvs*x1c*x2*z+4*vvs*x1c*x2-4*vvs*x1c
     #    *x2s-8*vvs*x2*z-7*vvs*x2*zs-2*vvs*x2*zc-3
     #    *vvs*x2-5*vvs*x2s*z-5*vvs*x2s*zs-2*vvs*
     #    x2s*zc-2*vvs*x2s+8*vvs*x2c*z+4*vvs*x2c*
     #    zs+4*vvs*x2c-15*vvq*x1*x2*z-3*vvq*x1*x2*zs
     #    -14*vvq*x1*x2-6*vvq*x1*x2s*z+6*vvq*x1*x2c+6
     #    *vvq*x1s*x2*z+10*vvq*x1s*x2-2*vvq*x1s*x2s
     #    -2*vvq*x1c*x2+9*vvq*x2*z+4*vvq*x2*zs+6*
     #    vvq*x2+7*vvq*x2s*z+5*vvq*x2s*zs+2*vvq*
     #    x2s-4*vvq*x2c*z-6*vvq*x2c)
      dw2= pwsm2*y*(-2*vvs*x1+3*vvs*z+2*vvs
     #    -36*x1*x2*z-48*x1*x2-12*x1*x2s-24*x1*z-6*x1*zs
     #    -18*x1+30*x1s*x2+24*x1s*z+36*x1s-18*x1c+
     #    24*x2*z+6*x2*zs+18*x2+12*x2s*z+12*x2s)
      dw2= dw2+y*rwms*rwmpgs*(6*vvs*x1*x2+2*
     #    vvs*x1*z+4*vvs*x1-4*vvs*x1s-4*vvs*x2*z-6*
     #    vvs*x2-vvs*z)
      dw2= dw2+ys*rwms*(3*vv*x1*x2*z+24*vv*
     #    x1*x2s*z+24*vv*x1*x2s-57*vv*x1*z-33*vv*x1*zs-6*
     #    vv*x1*zc-30*vv*x1-18*vv*x1s*x2*z-18*vv*x1s*x2-
     #    12*vv*x1s*x2s+39*vv*x1s*z+12*vv*x1s*zs+33*vv*
     #    x1s+12*vv*x1c*x2-6*vv*x1c*z-12*vv*x1c+15*vv*
     #    x2*z+15*vv*x2*zs+6*vv*x2*zc+6*vv*x2-24*vv*x2s*z
     #    -12*vv*x2s*zs-12*vv*x2s+24*vv*z+21*vv*zs+6*
     #    vv*zc+9*vv+4*vvc*x1*x2*z-4*vvc*x1*x2-12*vvc*
     #    x1*x2s+15*vvc*x1*z+3*vvc*x1*zs+14*vvc*x1+6*
     #    vvc*x1s*x2-6*vvc*x1s*z-10*vvc*x1s+2*vvc*
     #    x1c-6*vvc*x2*z-5*vvc*x2*zs-2*vvc*x2+8*vvc
     #    *x2s*z+12*vvc*x2s-9*vvc*z-4*vvc*zs-6*
     #    vvc)
      dw2= dw2+yc*(57*vvs*x1*x2*z+33*vvs*
     #    x1*x2*zs+6*vvs*x1*x2*zc+30*vvs*x1*x2-3*vvs*x1
     #    *x2s*z-24*vvs*x1*x2c*z-24*vvs*x1*x2c-39*vvs
     #    *x1s*x2*z-12*vvs*x1s*x2*zs-33*vvs*x1s*x2+18
     #    *vvs*x1s*x2s*z+18*vvs*x1s*x2s+12*vvs*x1s*
     #    x2c+6*vvs*x1c*x2*z+12*vvs*x1c*x2-12*vvs*
     #    x1c*x2s-24*vvs*x2*z-21*vvs*x2*zs-6*vvs*x2*
     #    zc-9*vvs*x2-15*vvs*x2s*z-15*vvs*x2s*zs-
     #    6*vvs*x2s*zc-6*vvs*x2s+24*vvs*x2c*z+12*
     #    vvs*x2c*zs+12*vvs*x2c-15*vvq*x1*x2*z-3*
     #    vvq*x1*x2*zs-14*vvq*x1*x2-6*vvq*x1*x2s*z+6*
     #    vvq*x1*x2c+6*vvq*x1s*x2*z+10*vvq*x1s*x2-2*
     #    vvq*x1s*x2s-2*vvq*x1c*x2+9*vvq*x2*z+4*vvq
     #    *x2*zs+6*vvq*x2+7*vvq*x2s*z+5*vvq*x2s*zs
     #    +2*vvq*x2s-4*vvq*x2c*z-6*vvq*x2c)
*
       else if(jj.eq.6) then
*
      wcom= pwsm2*rdms*(4*vv*x1*z+4*vv*x1-2*vv*x1s-
     #    4*vv*z-2*vv*zs-2*vv)
      dw1= pwsm2*y*(-2*vvs*x1*z-2*vvs*x1+vvs
     #    *x1s+2*vvs*z+vvs*zs+vvs-12*x1*x2*z-4*x1*
     #    x2*zs-8*x1*x2-4*x1*x2s*z-4*x1*x2s-4*x1*z-2*x1
     #    *zs-2*x1+8*x1s*x2*z+10*x1s*x2+2*x1s*x2s+8
     #    *x1s*z+2*x1s*zs+6*x1s-4*x1c*x2-4*x1c*z-
     #    6*x1c+2*x1q+4*x2*z+2*x2*zs+2*x2+4*x2s*z+2
     #    *x2s*zs+2*x2s)
      dw1= dw1+pwsm2*rdms*(-12*vvi*x1*x2*
     #    z-14*vvi*x1*x2-2*vvi*x1*x2s-14*vvi*x1*
     #    z-4*vvi*x1*zs-10*vvi*x1+8*vvi*x1s*x2
     #    +10*vvi*x1s*z+14*vvi*x1s-6*vvi*x1c
     #    +10*vvi*x2*z+4*vvi*x2*zs+6*vvi*x2+2*
     #    vvi*x2s*z+2*vvi*x2s+4*vvi*z+2*
     #    vvi*zs+2*vvi+3*vv*x1-3*vv*z-3*vv)
      dw2= pwsm2*y*(-2*vvs*x1*z-2*vvs*x1+vvs
     #    *x1s+2*vvs*z+vvs*zs+vvs-36*x1*x2*z-12*x1*
     #    x2*zs-24*x1*x2-12*x1*x2s*z-12*x1*x2s-12*x1*z-
     #    6*x1*zs-6*x1+24*x1s*x2*z+30*x1s*x2+6*x1s*
     #    x2s+24*x1s*z+6*x1s*zs+18*x1s-12*x1c*x2-
     #    12*x1c*z-18*x1c+6*x1q+12*x2*z+6*x2*zs+6*x2
     #    +12*x2s*z+6*x2s*zs+6*x2s)
      dw2= dw2+pwsm2*rdms*(-36*vvi*x1*x2*
     #    z-42*vvi*x1*x2-6*vvi*x1*x2s-42*vvi*x1*
     #    z-12*vvi*x1*zs-30*vvi*x1+24*vvi*x1s*
     #    x2+30*vvi*x1s*z+42*vvi*x1s-18*vvi*
     #    x1c+30*vvi*x2*z+12*vvi*x2*zs+18*vvi*
     #    x2+6*vvi*x2s*z+6*vvi*x2s+12*vvi*z+
     #    6*vvi*zs+6*vvi+3*vv*x1-3*vv*z-3*vv)
*
       endif
       wcom= wcom*(2.d0-y)**2/ys-dw1+2.d0/y*(1.d0-1.d0/y)*dw2
*
      endif
      www= y*wcom/pwtm2
*
    2 if(iz.eq.0) then
         ww(kk)= 0.d0
         iz= 1
         inull= inull+1
      else
         ww(kk)= www*tjac*tfactsw*stf*psf/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
       wtoxsswc= 0.d0
      else
       wtoxsswc= wwtt
      endif
*
      return
      end