*
*----------------------------------------------------------------
*
      subroutine wto(rs,wmi,zmi,zgi,tqmi,hmi,alszi,sig,esig)
      implicit real*8 (a-h,o-z)
      character*1,oxcm,ocoul,oud,oqcd,oglu,opeak,ofl,ofsr,oral,obin,
     #            oprt,opeakn,om,osm,oww,ozz,opeaka,ostop,rio,
     #            oanom,otrans,omh,ord,otop,ostore,omssm,ockm,
     #            omdist,opglu,osmass,oint,oseed
      character*2,ofs,oev
      character*4,otype
*
      parameter (nin=98,nout=99,ndimmx=9,ninv=10,liw=ndimmx+2,
     #           lw=ndimmx*(ndimmx-1)/2+12*ndimmx)
      parameter (itmxmx=100000,im=100,npos=512,ifmax=10000)
      parameter(nbf=2,nt=1,lwa=(nbf*(3*nbf+13))/2,lwat=(nt*(3*nt+13))/2)
*
      common/wtmp/zrm
      common/wtai/rio
      common/wtfb/oev
      common/wtcw/oww
      common/wtcz/ozz
      common/wtmod/om
      common/wtfs/ofs
      common/wtud/oud
      common/wtoi/oint
      common/wtii/iint
      common/wtsfh/ip0
      common/wtkr/ireg
      common/wtcb/obin
      common/wtfls/ofl
      common/wtqcd/als
      common/wtgg/oglu
      common/wtrms/rsm2
      common/wtim/ostop
      common/wtprt/oprt
      common/wtfsr/ofsr
      common/wtopa/delc
      common/wtps/opeak
      common/wthiggs/hm
      common/wtdis/dist
      common/wtsmod/osm
      common/wtafsz/g2z
      common/wtkount/ik
      common/wtckm/ockm
      common/wtickm/ickm
      common/wthqcd/alsh
      common/wtistrf/isf
      common/wtbme/bfact
      common/wthiggo/omh
      common/wtoral/oral
      common/wtlmsb/qcdl
      common/wtqcdz/alsz
      common/wtaqcd/oqcd
      common/wtap/opeaka
      common/wtcqcd/iqcd
      common/wtcac/oanom
      common/wtgnum/itmx
      common/wtsp/psg(4)
      common/wtfsm/fsm(4)
      common/wtpqcd/opglu
      common/wopst/ostore
      common/wtsf/ix0,it0
      common/wtpsn/opeakn
      common/wtomd/omdist
      common/wtseed/oseed
      common/wtcoul/ocoul
      common/wttopop/otop
      common/wtchi/hch(36)
      common/wtspl/prx0(9)
      common/wopstm/osmass
      common/wtshel/otrans
      common/wtipt/ifz(51)
      common/wtmssmo/omssm
      common/wtcurrent/oxcm
      common/wtnf/ifl(npos)
      common/wticuts/iac(4)
      common/wtrm/rbm2,rcm2
      common/wtnpr/ipr,ipr0
      common/wtisa/isaa,isab
      common/wttopt/ios,iosf
      common/wtspln/prx0n(9)
      common/wtmd/arrinv(10)
      common/wtvckm/vckm(3,3)
      common/wtmatx/colf(8,8)
      common/wtparc/xap(ninv)
      common/wtcx/xscmx(npos)
      common/wthx/xshmx(npos)
      common/wtochannel/otype
      common/wtpoints/ipm,irm
      common/wtparh/xaph(ninv)
      common/wtpmx/xmx(npos,9)
      common/wtdis2/distm,distp
      common/wtbac/dg1z,dkg,rlg
      common/wtcx0/xscmx0(npos)
      common/wthx0/xshmx0(npos)
      common/wtrmss/chcm2,chsm2
      common/wttc/itc,itcc,itcn             
      common/wtpmxh/xmxh(npos,9)
      common/wtlb/abp,bbp,abm,bbm
      common/wtstor/stry(npos,ifmax)
      common/wtncc/chf2,chfp2,conc(10)
      common/wtcclr/vupl,vupr,vdpl,vdpr
      common/wtiaz/alzr,alzi,alzir,alzii
      common/wtalqed/arezm,aizm,aresz,aisz
      common/wtcmplx/sw(2),sz(2),rsw(2),rsz(2)
      common/wtmssmi/am,tbeta,rmu,scalm,bat,bab
      common/wtnclr/vel,ver,velr,vfl,vfr,vfpl,vfpr
      common/wtcchannel/chu,chup,chd,chdp,fcuc,fcdc
      common/wtnchannel/chf,chfp,tif,tifp,fcun,fcdn
      common/wthapar/rhm,rhm2,rhg,rhmg,shg,shgs,opshgs
      common/wtee/qch,qch2,vqr,vql,hbe(24),hbo(24),hmp(24)
      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/wtacc/acg1g,aclg,ackg,acg4g,acktg,acltg,acg5g,acg1z,aclz,
     #           ackz,acg4z,acktz,acltz,acg5z
      common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2,
     #              tqm2,dmy2
      common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,g8,tfact,pih,alw,
     #              eta,feta,beta,g2,tfacth
      common/wtacchannel/omchu,opchu,omchup,opchup,omchdp,opchdp,
     #                   omchd,opchd,hchup,hchu,hchdp,hchd
      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/wtsubreg/dsm,usm,dsp,usp,rl(6),rr(6),srl(6),sdsm,sdsp,vvl1,
     #              vvl2,vvl3,ul,omul,suml
      common/wtcuts/aim(6),bim(6),ae(4),asa(4),bsa(4),afsa(6),bfsa(6),
     #            ombsa(4),opbsa(4),teq,rae(4),omasa(4),opasa(4),
     #            sg12,cg12,sg13,cg13,sg14,cg14,sg23,cg23,sg24,
     #            cg24,sg34,cg34,sct120,sct130,sct140,sct230,
     #            sct240,sct340,sgam(4),cgam(4)
      common/wtmssm/ams,shm,shms,bhm,bhms,sbeta,cbeta,salpha,calpha,
     #              sbma,cbma,rshm,rshm2,rshg,rshmg,sshg,sshgs,
     #              opsshgs,rbhm,rbhm2,rbhg,rbhmg,sbhg,sbhgs,opsbhgs,
     #              ram,ram2,rag,ramg,sag,sags,opsags
      common/wtmssmc/chms,chm,rchn,rchm2,rchg,rchmg,schg,schgs,opschgs
*
      dimension xmx0(9)
      dimension sclm(4)
      dimension xmxh0(9)
      dimension idhep(6)
      dimension phep(5,6)
      dimension aphep(5,4)
      dimension vk(ndimmx)
      dimension prx(ndimmx)
      dimension zs(2),zl(2)
      dimension aprx(ndimmx)
      dimension prxa(ninv,itmxmx)
      dimension axstr(100),ir(100)
      dimension vol(npos,4),volt(4)
      dimension exc(im),extst(itmxmx)
      dimension excl(im),qtrf(itmxmx)
      dimension volh(npos,8),volht(8)
      dimension excwp(im),cwpst(itmxmx)
      dimension volhq(npos,12),volhtq(12)
      dimension volhb(npos,20),volhtb(20)
      dimension st(nt),fvect(nt),wat(lwat)
      dimension gbl(npos,9,4),gbu(npos,9,4)
      dimension gblh(npos,9,8),gbuh(npos,9,8)
      dimension gblhq(npos,9,12),gbuhq(npos,9,12)
      dimension gblhb(npos,9,20),gbuhb(npos,9,20)
      dimension exce(im),exetst(itmxmx),exc4j(itmxmx)
      dimension wmx(lw),bl(ndimmx),bu(ndimmx),iw(liw)
      dimension fvecw(nbf),waw(lwa),waz(lwa),fvecz(nbf)
      dimension pggfz(2),pgglqz(2),p3qz(2),fzz(2),fwz(2)
      dimension amps(itmxmx),amps1(itmxmx),amps2(itmxmx)
      dimension siga(4),sigan(4),sigah(8),sigahq(12),sigahb(20)
      dimension xsmx(npos,4),xsmxn(npos,4),xsmxh(npos,8),
     #          xsmxhq(npos,12),xsmxhb(npos,20)
      dimension oxsmx(npos,4),oxsmxn(npos,4),oxsmxh(npos,8),
     #          oxsmxhq(npos,12),oxsmxhb(npos,20)
      dimension exc1(im),extst1(itmxmx),exc2(im),extst2(itmxmx),
     #          exc3(im),exc4(im),exc5(im),exc6(im),exc7(im),
     #          extst3(itmxmx),extst4(itmxmx),extst5(itmxmx),
     #          extst6(itmxmx),excm(im)
*
      external d01gcf,wtoregion,wtoxsc,wtoxsn,wtoxsm43,wtoxsnh19,
     #         wtoxsnh24,wtoxsnh49,wtoxsnh32,wtoxsnh64,wtoxsn32,
     #         wtoxsn64,wtoxsnh26,wtoxsnsh64,wtoxssc
      external s14aaf,s21baf,x04cbf,x04abf,s09abf
      external g05caf,g05cbf,g05ccf,e04jaf
      external c05nbf,c02ajf,x02ajf,wtofsw,wtofsz,wtofst
*
      if(osm.eq.'g'.and.oseed.eq.'y') then
         call g05ccf
      endif
*
      ik= 0
      do i=1,51
          ifz(i)= 0
      enddo
*
      wm= wmi
      zm= zmi
      zg= zgi
      alsz= alszi
*
      nf= 5
      if(iqcd.eq.2) then
         x1= 0.001d0
         x2= 10.0d0
         xacc= 1.d-12
         nf= 5
         qcdl= wtoqcdlam(nf,als,wm,x1,x2,xacc)
      else
         qcdl= 0.d0
      endif
      als= wtoralphas(zm,wm,alsz,nf)
*
      ars= rs
      pis= pi*pi
      em2= em*em
      amm2= rmm*rmm
      tm2= tm*tm
      rnm2= rnm*rnm
      uqm2= uqm*uqm
      dqm2= dqm*dqm
      bqm2= bqm*bqm
      cqm2= cqm*cqm
      sqm2= sqm*sqm
      bqm2= bqm*bqm
      dmy2= dmy*dmy
      hm= hmi
*
      do i=1,4
          rae(i)= ae(i)/rs
      enddo
      teq= rae(1)+rae(2)+rae(3)+rae(4)
      sct120= 2.d0*rae(1)*rae(2)*(1.d0-bfsa(1))
      sct130= 2.d0*rae(1)*rae(3)*(1.d0-bfsa(2))
      sct140= 2.d0*rae(1)*rae(4)*(1.d0-bfsa(3))
      sct230= 2.d0*rae(2)*rae(3)*(1.d0-bfsa(4))
      sct240= 2.d0*rae(2)*rae(4)*(1.d0-bfsa(5))
      sct340= 2.d0*rae(3)*rae(4)*(1.d0-bfsa(6))
*
*-----parameters
*
      alpha= 1.d0/alphai
      alw= 1.d0/alwi
      pih= 0.5d0*pi
      s= rs*rs
      qgfopi= 0.25d0*gf/pi
      wm2= wm*wm
      wmc= wm2*wm
      zm2= zm*zm
      wg0= 1.5d0*gf*wmc/pi
      if(oqcd.eq.'n') then
         wg= wg0
      else
         wg= wg0*(1.d0+2.d0/3.d0*als/pi)
      endif
      wg2= wg*wg
      if(ios.eq.1) then
         sth2= 0.5d0*pi*alw/gf/wm2
         g2= 4.d0*pi*alw/sth2
      else if(ios.eq.2) then
         sth2= 1.d0-wm2/zm2
         g2= 8.d0*gf*wm2
      else if(ios.eq.3) then
         scal= -zm2
         call wtopself(scal,pggf)
         derl= 0.25d0*alpha/pi*pggf
         call wtohadr5(zm,derh,ederh)
         der= derl+derh
         alz= (1.d0-der)/alpha
         args= 1.d0-2.d0*pi/alz/gf/zm2
         sth2= 0.5d0*(1.d0-sqrt(args))
         g2= 8.d0*gf*zm2*(1.d0-sth2)
      endif
      g4= g2*g2
      g8= g4*g4
      cth2= 1.d0-sth2
      tth2= sth2/cth2
      hsth2= 0.5d0*sth2
      tsth2= 2.d0*sth2
      scth2= 0.125d0/cth2
      asth2= 0.5d0-sth2
      sth4= sth2*sth2
      cth4= cth2*cth2
      psg(1)= 1.d0
      psg(2)= -1.d0
      psg(3)= -1.d0
      psg(4)= 1.d0
      if(ofl.eq.'a'.and.oanom.eq.'r') then
         acg1g= 0.d0
         aclg= rlg
         ackg= dkg
         acg4g= 0.d0
         acktg= 0.d0
         acltg= 0.d0
         acg5g= 0.d0
         acg1z= dg1z
         aclz= rlg
         ackz= dg1z-sth2/cth2*dkg
         acg4z= 0.d0
         acktz= 0.d0
         acltz= 0.d0
         acg5z= 0.d0
      endif
*
*-----mass&width parameters
*
      rwm= wm/rs
      rwm2= wm2/s
      rwg= wg/rs
      swg= wg/wm
      swgs= swg*swg
      if(oww.eq.'r'.or.oww.eq.'f') then
         rwmg= rwm*rwg
      else if(oww.eq.'i') then
         rwmg= rwm*(1.d0-rwg*rwg/rwm/rwm)*rwg
      else if(oww.eq.'t') then
         rwmg= (1.d0+swgs)*rwm*rwg
      endif
      opswgs= 1.d0+swgs
      s0w= rwm2/opswgs
      rzm= zm/rs
      rzm2= zm2/s
      rzg= zg/rs
      if(ozz.eq.'r'.or.ozz.eq.'f') then
         rzmg= rzm*rzg
      else if(ozz.eq.'i') then
         rzmg= rzm*(1.d0-rzg*rzg/rzm/rzm)*rzg
      endif
      szg= zg/zm
      szgs= szg*szg
      opszgs= 1.d0+szgs
      rszw= zg/zm
      rszw2= rszw*rszw
      s0z= rzm2/opszgs
*
      tqm= tqmi
      if(ofl.eq.'c') then
         if(rio.eq.'i') then
          tqm2= tqm*tqm
         else if(rio.eq.'a') then
          if(otop.eq.'f') then
             tqm2= tqm*tqm
          else if(otop.eq.'d') then
*
*-----Iteration for m_t starts
*
             gmt= 150.d0
             st(1)= gmt*gmt
             tol= sqrt(x02ajf())
             ifail= 1
             call c05nbf(wtofst,nt,st,fvect,tol,wat,lwat,ifail)
             tqm2= st(1)
             tqm= sqrt(tqm2)
          endif
         endif
*
*-----alpha(zm)
*
         pzr= -zm2
         pzi= 0.d0
         call wtopoleg(pzr,pzi,pggfz,pgglqz,pggnpz)
         call wtopolez0(pggf0,fw0)
         xal= alphai-0.25d0*(pggfz(1)-pggf0+pggnpz)/pi
         yal= -0.25d0*(pggfz(2)+pgglqz(2))/pi
         alzir= xal
         alzii= yal
         alzr= xal/(xal*xal+yal*yal)
         alzi= -yal/(xal*xal+yal*yal)
*
*-----Re(1/g^2(zm))
*
         apis= 16.d0*pis
         p2zr= -zm2
         p2zi= 0.d0
         call wtopole(p2zr,p2zi,p3qz,fzz,fwz)
         bx= 1.d0/8.d0/gf/zm2+(fw0-fzz(1))/zm2/apis
         xih= -p3qz(2)/apis
         agc= 4.d0*pi*alzr
         bgc= -1.d0-8.d0*pi*xih*alzi
         cgc= -4.d0*pi*xih*xih*alzr+bx
         ifail= 0
         call c02ajf(agc,bgc,cgc,zs,zl,ifail)
         g2z= zs(1)
*
*-----Iteration for sw starts
*
         sw(1)= wm2
         sw(2)= -wm*wg
         tol= sqrt(x02ajf())
         ifail= 1
         call c05nbf(wtofsw,nbf,sw,fvecw,tol,waw,lwa,ifail)
         wmu= sqrt(sw(1))
         wlg= -sw(2)/wmu
         do l=1,2
             rsw(l)= sw(l)/s
         enddo
*
*-----Iteration for sz starts
*
         dszi= 1.d-1
         szi0= zm*zg
         izt= 1
 1005    sz(1)= zm2
         sz(2)= -szi0*(1.d0-(izt-1)*dszi)
         tol= sqrt(x02ajf())
         ifail= 1
         call c05nbf(wtofsz,nbf,sz,fvecz,tol,waz,lwa,ifail)
         if(ifail.eq.0) then
            zmu= sqrt(sz(1))
            zlg= -sz(2)/zmu
            do l=1,2
                rsz(l)= sz(l)/s
            enddo
         else
            if(izt.lt.20) then
               izt= izt+1
               go to 1005
            else
               print*,' Failure in searching for Z pole'
            endif
         endif
      endif
      tqm2= tqm*tqm
*
      if((otype.eq.'nc21'.or.otype.eq.'nc25'.or.
     #   otype.eq.'nc33'.or.otype.eq.'nc50'.or.
     #   otype.eq.'nc68').and.omssm.eq.'n') then
         call wtocorrqcd(hm,als,bqm2,cqm2,rbm2,rcm2,rsm2)
         if(oqcd.eq.'y') then
            nf= 5
*            alsh= wtoralphas(wm,hm,als,nf)
*            hg= qgfopi*hm*(3.d0*(rbm2+rcm2)*(1.d0+alsh/pi*(
*     #          17.d0/3.d0+(35.94d0-1.36d0*nf)*alsh/pi))+tm*tm)
*            hgb= qgfopi*hm*(3.d0*rbm2*(1.d0+alsh/pi*(
*     #           17.d0/3.d0+(35.94d0-1.36d0*nf)*alsh/pi)))
*            hgr= qgfopi*hm*(3.d0*rcm2*(1.d0+alsh/pi*(
*     #           17.d0/3.d0+(35.94d0-1.36d0*nf)*alsh/pi))+tm*tm)
*/////////////////////////////////////////////////////////////////////////
            alsh= wtoralphas(wm,hm,als,nf)
            scalh= -hm*hm
            call wtopself(scalh,pggfh)
            derlh= 0.25d0*alpha/pi*pggfh
            call wtohadr5(hm,derhh,ederhh)
            derh= derlh+derhh
            alh= alpha/(1.d0-derh)
            hm2= hm*hm
            hctotc= alsh/pi*(17.d0/3.d0+(35.94d0-1.36d0*nf)*
     #              alsh/pi)
            hcorrb= -6.d0*rbm2/hm2+0.472*alh/pi+alsh/pi*(0.651d0*
     #              alh/pi+5.667d0-40.d0*rbm2/hm2+alsh/pi*(
     #              29.14671d0-87.725d0*rbm2/hm2+41.75760d0*alsh/pi))
            clnt= 2*log(hm/tqm)
            hcorrt= alsh*alsh/pis*(3.111d0-0.667*clnt+rbm2/hm2*(-10.d0+
     #              4.d0*clnt+4.d0/3.d0*log(rbm2/hm2))+hm2/tqm2*(
     #              0.241d0-0.070d0*clnt)+alsh/pi*(
     #              50.474d0-clnt*(8.167d0+1.278*clnt)))+gf*tqm2/
     #              8.d0/pis*(1+alsh/pi*(-4.913d0+alsh/pi*(-72.117d0-
     #              20.945d0*clnt)))
            hctotb= hcorrb+hcorrt
            hg= qgfopi*hm*(3.d0*rcm2*(1.d0+hctotc)+3.d0*rbm2*
     #          (1.d0+hctotb)+tm*tm)
            print*,qgfopi*hm*3.d0*rcm2*(1.d0+hctotc),
     #             qgfopi*hm*3.d0*rbm2*(1.d0+hctotb),
     #             qgfopi*hm*tm*tm
            hgb= qgfopi*hm*3.d0*rbm2*(1.d0+hctotb)
            hgr= qgfopi*hm*(3.d0*rcm2*(1.d0+hctotc)+tm*tm)
*/////////////////////////////////////////////////////////////////////////
         else if(oqcd.eq.'n') then
            hg= qgfopi*hm*(3.d0*(rbm2+rcm2)+tm*tm)
            hgb= qgfopi*hm*3.d0*rbm2
            hgr= qgfopi*hm*(3.d0*rcm2+tm*tm)
         endif
         if(oglu.eq.'y') then
            nf= 5
            alsh= wtoralphas(wm,hm,als,nf)
            hg= hg+gf/36.d0/pi**3*alsh*alsh*hm**3*(1.d0+
     #          (95.d0/4.d0-7.d0/6.d0*nf)*alsh/pi)
            print*,gf/36.d0/pi**3*alsh*alsh*hm**3*(1.d0+
     #          (95.d0/4.d0-7.d0/6.d0*nf)*alsh/pi)
            hgr= hgr+gf/36.d0/pi**3*alsh*alsh*hm**3*(1.d0+
     #          (95.d0/4.d0-7.d0/6.d0*nf)*alsh/pi)
         endif
*         hm2= hm*hm
         hm3= hm2*hm
         if(hm.gt.(2.d0*wm)) then
            xwh= wm2/hm2
            btw= sqrt(1.d0-4.d0*xwh)
            hgw= gf*hm3/8.d0/pi*(1.d0-4.d0*xwh*(1.d0-3.d0*xwh))*
     #           btw
            print*,hgw
            hg= hg+hgw
            hgr= hgr+hgw
         endif
         i
O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O





















































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































bqm2= rbm2/s 
O
      else if(omssm.eq.'y') then
O
         rbqm2= shbm2/s(
O
      endifc

O
*0
O
      chf2= chf*chf(
O
      chfp2= chfp*chfp"
O
      conc(1)= chf*chfp*sth4/16.d0
O
      conc(2)= 0.25d0*chf*tth2
O
      conc(3)= 0.25d0*chfp*tth23#
O
      conc(4)= 0.25d0*chf*chfp*tth2
O

O
      conc(5)= 6.25d-2/cth4 #
O
      conc(6)= chf*chfp2*sth4/16.d0d#
O
      conc(7)= chf2*chfp*sth4/16.d0=
O
      conc(8)= 3.125d-2/cth2
O
      conc(9)= chf*sth2/8.d0
O
      conc(10)= chf*sth2/16.d0

O
* 
O
      ve= -0.5d0+2.d0*sth2
O
      vf= tif-2.d0*chf*sth2 
O
      vfp= tifp-2.d0*chfp*sth2
O
      vel= ve-0.5d0.
O
      vfl= vf+tifn
O
      vfpl= vfp+tifp
O
      ver= ve+0.5d0t
O
      velr= ver*vel0
O
      vfr= vf-tif
O

O
      vfpr= vfp-tifp
O
      vel2= vel*veli
O
      ver2= ver*ver 
O
      vfl2= vfl*vflm
O
      vfr2= vfr*vfr*
O
      vfpl2= vfpl*vfpl
O
      vfpr2= vfpr*vfpr 
O
      vdp= -0.5d0-2.d0*chdp*sth2
O
      vdpl= vdp-0.5d0d
O
      vdpr= vdp+0.5d0g
O
      vup= 0.5d0-2.d0*chup*sth2h
O
      vupl= vup+0.5d0 
O
      vupr= vup-0.5d0z1
O
      if(otype.eq.'nc48'.or.otype.eq.'nc50') then/
O
         vqr= vfrh
O
         vql= vfl(
O
         ver3= ver2*ver*
O
         vel3= vel2*vel 
O
         vqr2= vqr*vqr
O
         vql2= vql*vql

O
*t
O
         hbo(1)= ver*vql
O
         hbo(2)= velrh
O
         hbo(3)= ver2*vel*vqld
O
         hbo(4)= vel*vqr
O
         hbo(5)= velr0
O
         hbo(6)= ver*vel2*vqrh
O
         hbo(7)= ver*vqr
O
         hbo(8)= velr 
O
         hbo(9)= ver2*vel*vqr)
O
         hbo(10)= vel*vql0
O
         hbo(11)= velr
O
         hbo(12)= ver*vel2*vql
O
         hbo(13)= ver*vqrr
O
         hbo(14)= ver2
O
         hbo(15)= ver3*vqr
O
         hbo(16)= vel*vqll
O
         hbo(17)= vel2  
O
         hbo(18)= vel3*vql
O
         hbo(19)= ver*vql 
O
         hbo(20)= ver2 
O

O
         hbo(21)= ver3*vql
O
         hbo(22)= vel*vqra
O
         hbo(23)= vel2 z
O
         hbo(24)= vel3*vqr -

O
*a
O
         hbe(1)= vel*vql  
O
         hbe(2)= velrw
O
         hbe(3)= ver*vel2*vql 
O
         hbe(4)= ver*vqr
O
         hbe(5)= velr 
O
         hbe(6)= ver2*vel*vqr 
O
         hbe(7)= vel*vqr
O
         hbe(8)= velrh
O
         hbe(9)= ver*vel2*vqrg
O
         hbe(10)= ver*vqll
O
         hbe(11)= velr
O
         hbe(12)= ver2*vel*vql
O
         hbe(13)= ver*vqrw
O
         hbe(14)= ver2 =
O
         hbe(15)= ver3*vqr
O
         hbe(16)= vel*vql 
O
         hbe(17)= vel2
O
         hbe(18)= vel3*vql
O
         hbe(19)= ver*vql
O

O
         hbe(20)= ver2
O
         hbe(21)= ver3*vql
O
         hbe(22)= vel*vqr
O

O
         hbe(23)= vel2 z
O
         hbe(24)= vel3*vqr

O
* 
O
         hmp(1)= vel*vql
O
         hmp(2)= ver*vql
O
         hmp(3)= velr*vql2
O
         hmp(4)= ver*vqr
O
         hmp(5)= vel*vqr
O
         hmp(6)= velr*vqr2
O
         hmp(7)= vel*vqr
O
         hmp(8)= ver*vqr r
O
         hmp(9)= velr*vqr2
O
         hmp(10)= ver*vql 
O
         hmp(11)= vel*vql 
O
         hmp(12)= velr*vql2
O

O
         hmp(13)= ver*vqr_
O
         hmp(14)= ver*vqrg
O
         hmp(15)= ver2*vqr2)
O
         hmp(16)= vel*vqll
O
         hmp(17)= vel*vql 
O
         hmp(18)= vel2*vql2l
O
         hmp(19)= ver*vql,
O
         hmp(20)= ver*vql 
O
         hmp(21)= ver2*vql2 
O
         hmp(22)= vel*vqr 
O
         hmp(23)= vel*vqr 
O
         hmp(24)= vel2*vqr2 

O
      else
O
         vqr= 0.d0
O
         vql= 0.d0
O
         do i=1,24
O
             hbe(i)= 0.d0l
O
             hbo(i)= 0.d0 
O
             hmp(i)= 0.d0g
O
         enddo
O
      endif 

O
*y
O
      hch(1)= vfl*ver/16.d02
O
      hch(2)= vfr*vel/16.d0
O

O
      hch(3)= vfr*ver/16.d0 
O
      hch(4)= vfl*vel/16.d0
O

O
      hch(5)= vfpl*ver/16.d0
O
      hch(6)= vfpr*vel/16.d0
O
      hch(7)= vfpr*ver/16.d0
O
      hch(8)= vfpl*vel/16.d0!
O
      hch(9)= ver2*vfl*vfpl/16.d0t"
O
      hch(10)= vel2*vfr*vfpr/16.d0"
O
      hch(11)= ver2*vfl*vfpr/16.d0"
O
      hch(12)= vel2*vfr*vfpl/16.d0"
O
      hch(13)= ver2*vfr*vfpl/16.d0"
O
      hch(14)= vel2*vfl*vfpr/16.d0"
O
      hch(15)= ver2*vfr*vfpr/16.d0"
O
      hch(16)= vel2*vfl*vfpl/16.d0
O
      hch(17)= vfr*vfpl/16.d0a
O
      hch(18)= vfl*vfpr/16.d0
O

O
      hch(19)= vfl*vfpl/16.d0
O

O
      hch(20)= vfr*vfpr/16.d0 "
O
      hch(21)= ver*vfr*vfpl2/16.d0"
O
      hch(22)= vel*vfl*vfpr2/16.d0"
O
      hch(23)= ver*vfl*vfpl2/16.d0"
O
      hch(24)= vel*vfr*vfpr2/16.d0"
O
      hch(25)= ver*vfr*vfpr2/16.d0"
O
      hch(26)= vel*vfl*vfpl2/16.d0"
O
      hch(27)= vel*vfr*vfpl2/16.d0"
O
      hch(28)= ver*vfl*vfpr2/16.d0"
O
      hch(29)= ver*vfr2*vfpl/16.d0"
O
      hch(30)= vel*vfl2*vfpr/16.d0"
O
      hch(31)= ver*vfl2*vfpl/16.d0"
O
      hch(32)= vel*vfr2*vfpr/16.d0"
O
      hch(33)= ver*vfr2*vfpr/16.d0"
O
      hch(34)= vel*vfl2*vfpl/16.d0"
O
      hch(35)= ver*vfl2*vfpr/16.d0"
O
      hch(36)= vel*vfr2*vfpl/16.d0

O
*r
O
      hchup= tsth2*chup=
O
      hchu= tsth2*chu 
O
      hchdp= tsth2*chdp 
O
      hchd= tsth2*chd 

O
* 
O
*-----CC

O
* 
O
      g2f= 8.d0*gf*wm2
O
      cth= sqrt(cth2) 
O
      gca= sqrt(g2)+
O
      gcf= sqrt(g2f)+
O
      gvea= 0.25d0*gca/cth*(4.d0*sth2-1.d0)i
O
      gaea= -0.25d0*gca/cthr 
O
      gwfa= 0.5d0*gca/sqrt(2.d0)
O
      g3va= -gca*cth+
O
      gvef= 0.25d0*gcf/cth*(4.d0*sth2-1.d0).
O
      gaef= -0.25d0*gcf/ctho 
O
      gwff= 0.5d0*gcf/sqrt(2.d0)
O
      g3vf= -gcf*cth

O
*#%
O
*-----e.m. quantities are initializedh

O
*
O

O
      aexp= alpha/pi
O
      rln= log(s/em2)2 
O
      beta= 2.d0*aexp*(rln-1.d0)
O
      isf= iosf 
O
      if(iosf.eq.1) then
O
         eta= beta
O
         eob= 1.d0
O
      else if(iosf.gt.1) then2
O
         eta= 2.d0*aexp*rln 
O
         eob= rln/(rln-1.d0)
O
      endifl
O
      feta= eta/16.d0 
O
      hbet= beta/2.d0d
O
      heta= eta/2.d0
O
      hbeti= 2.d0/beta
O
      omhb= 1.d0-hbetn
O
      ophb= 1.d0+hbet 
O
      if(iosf.lt.3) then c"
O
         arge= hbet*(3.d0/4.d0-ge)
O
      else if(iosf.eq.3) then)&
O
         arge= -hbet*ge+3.d0/4.d0*heta
O
      endif/
O
      ifg= 1
O
      gamb= s14aaf(ophb,ifg)
O
      if(ifg.ne.0) thens
         stop
O

O
      else c
O
         d0gl= exp(arge)/gambw
O
      endif,

O
*f
O
*-----splitting&transformationsh

O
*i
O
      pin= 256.d0*pi**7w
O
      if(oxcm.eq.'c') then,
O
         tfact= g8/128.d0/pin*fcuc*fcdc*cfct
O
      else if(oxcm.eq.'n') thenm1
O
         if((otype.eq.'nc64'.or.otype.eq.'nc68').3!
O
     #      and.ofs.eq.'qq') then %
O
            tfact= g8/128.d0/pin*cfct*
O
         else /
O
            tfact= g8/128.d0/pin*fcun*fcdn*cfct/
O
         endif
O
      else if(oxcm.eq.'m') then "
O
         tfact= g8/128.d0/pin*cfct
O
      endifl
O
      if(otype.eq.'nc68') then-
O
         tfacth= g8/128.d0/pin*fcun*fcdn*cfct3

O
      else
O
         tfacth= 0.d0
O

O
      endif 

O
* 
O
      do i=1,6
O
          rl(i)= aim(i)*aim(i)
O
          rr(i)= bim(i)*bim(i)
O
          srl(i)= sqrt(rl(i))#
O
      enddo 
O
      dsm= rl(1)
O
      usm= rr(1)
O
      dsp= rl(6)
O
      usp= rr(6)

O
*s
O
      sdsm= sqrt(dsm)l
O
      sdsp= sqrt(dsp) 

O
* "
O
*-----limits for u from cuts on IM

O
* 8
O
      vvl1= (sdsm+sdsp)*(sdsm+sdsp)                     ?
O
      vvl2= (srl(2)+srl(5))*(srl(2)+srl(5))                    t?
O
      vvl3= (srl(3)+srl(4))*(srl(3)+srl(4))                    #
O
      ul= dmax1(vvl1,vvl2,vvl3)(
O
      omul= 1.d0-ul#
O
      suml= 0.d0
O
      do j=1,4
O
          suml= suml+rae(j)b
O
      enddo                 

O
*=
O
      if(oxcm.eq.'c') then
O
         if(omdist.eq.'y') then/
O
            ndim= 2/
         else/
O
          if(itc.eq.11) then
O
             if(isf.eq.0) then+
O
                Print*,' Born not allowed '
O

O
                stop
O
             else 
O
                ndim= 7.
O
             endif 
O
          else if(itc.ge.7) then
O
             if(isf.eq.0) then
O
                ndim= 6s
O
             else 
O
                ndim= 8i
O
             endif
O
          else
O
             if(isf.eq.0) then
O
                ndim= 7r
O
             else*
O
                ndim= 9 
O
             endif
O
          endiff
O
         endif.
O
      else if(oxcm.eq.'n'.or.oxcm.eq.'m') then
O
         if(itc.gt.0) then
O
            if(isf.eq.0) then
O

O
               ndim= 6
O
            else  
O
               ndim= 8
O
            endif)
         else 
O
            if(isf.eq.0) thent
O
               ndim= 7
O
            else
O
               ndim= 9
O
            endif)
O
         endif

O
      else
O
         if(isf.eq.0) then
O
            ndim= 7 
         elser
O
            ndim= 9
O

O
         endif
O
      endift

O
*-
O
      if(ndim.eq.2) then
O
         if(inpts.le.6) then
O
            jnpts= inpts!
O
         else if(inpts.eq.7) thend
O
            vk(1)= 1.d0 
O
            vk(2)= 0.141603d6h
O
         endif*
O
      else if(ndim.eq.6.or.ndim.eq.8) then
O
         if(inpts.le.6) then
O
            jnpts= inpts!
O
         else if(inpts.eq.8) then 
O
            jnpts= 10193*101
O
            if(ndim.eq.6) then
O
               vk(1)= 1.d0"
O
               vk(2)= 0.672849d+06!
O
               vk(3)= 0.53093d+05 "
O
               vk(4)= 0.164857d+06"
O
               vk(5)= 0.114815d+06 
O
               vk(6)= 0.3215d+04#
O
            else if(ndim.eq.8) then
O

O
               vk(1)= 1.d0"
O
               vk(2)= 0.757894d+06"
O
               vk(3)= 0.784365d+06"
O
               vk(4)= 0.236855d+06"
O
               vk(5)= 0.347946d+06"
O
               vk(6)= 0.524281d+06"
O
               vk(7)= 0.128976d+06"
O
               vk(8)= 0.805687d+06
O
            endifb!
O
         else if(inpts.eq.9) thena
O
            jnpts= 22807*151
O
            if(ndim.eq.6) then
O
               vk(1)= 1.d0#
O
               vk(2)= 0.2903683d+07 "
O
               vk(3)= 0.278237d+06"
O
               vk(4)= 0.413956d+06#
O
               vk(5)= 0.1366666d+07
O
#
O
               vk(6)= 0.1522064d+07 #
O
            else if(ndim.eq.8) then=
O
               vk(1)= 1.d0#
O
               vk(2)= 0.1464524d+07("
O
               vk(3)= 0.329469d+07#
O
               vk(4)= 0.2417287d+07,!
O
               vk(5)= 0.33812d+05s#
O
               vk(6)= 0.2709542d+07e#
O
               vk(7)= 0.1615901d+07r"
O
               vk(8)= 0.249863d+06
O
            endifr
         else 1
O
            print*,' NPTS choice is Not allowed '.
O
            stop
O
         endif

O
      else
O
         if(inpts.le.6) then
O
            jnpts= inpts!
O
         else if(inpts.eq.7) thene
O
            jnpts= 2011*101=
O
            if(ndim.eq.7) then
O
               vk(1)= 0.1d+1 
O
               vk(2)= 0.72725d+5!
O
               vk(3)= 0.118296d+6d!
O
               vk(4)= 0.107084d+6v
O
               vk(5)= 0.1938d+4/!
O
               vk(6)= 0.185127d+6  
O
               vk(7)= 0.14844d+6#
O
            else if(ndim.eq.9) then1
O
               vk(1)= 0.1d+1 
O
               vk(2)= 0.82871d+5 
O
               vk(3)= 0.13509d+5!
O
               vk(4)= 0.159618d+6
O
 
O
               vk(5)= 0.99403d+5 
O
               vk(6)= 0.53186d+5 
O
               vk(7)= 0.68306d+5 
O
               vk(8)= 0.86067d+5 
O
               vk(9)= 0.12481d+5
O
            endifc!
O
         else if(inpts.eq.8) then 
O
            jnpts= 10193*101
O
            if(ndim.eq.7) then
O
               vk(1)= 0.1d+1"
O
               vk(2)= 0.1019431d+7 
O
               vk(3)= 0.35353d+6!
O
               vk(4)= 0.708948d+6 !
O
               vk(5)= 0.951714d+6x!
O
               vk(6)= 0.197618d+63 
O
               vk(7)= 0.54816d+6#
O
            else if(ndim.eq.9) then3
O
               vk(1)= 0.1d+01(!
O
               vk(2)= 0.755242d+6 !
O
               vk(3)= 0.911407d+6*!
O
               vk(4)= 0.442285d+6s!
O
               vk(5)= 0.850204d+6r!
O
               vk(6)= 0.572366d+6 "
O
               vk(7)= 0.1026802d+7!
O
               vk(8)= 0.892453d+6 
O
            vk(9)= 0.685582d+6
O
            endif)!
O
         else if(inpts.eq.9) thenh
O
            jnpts= 22807*151
O
            if(ndim.eq.7) then
O
               vk(1)= 0.1d+01 "
O
               vk(2)= 0.1619459d+7!
O
               vk(3)= 0.226133d+7i!
O
               vk(4)= 0.256381d+7s!
O
               vk(5)= 0.230245d+7/!
O
               vk(6)= 0.855081d+6h!
O
               vk(7)= 0.609193d+6 #
O
            else if(ndim.eq.9) then1
O
               vk(1)= 0.1d+01 "
O
               vk(2)= 0.2055793d+7!
O
               vk(3)= 0.767734d+6+"
O
               vk(4)= 0.3183104d+7"
O
               vk(5)= 0.2813063d+7"
O
               vk(6)= 0.2463708d+7 
O
               vk(7)= 0.58258d+5"
O
               vk(8)= 0.2817562d+7"
O
               vk(9)= 0.1276513d+7
O
            endifg"
O
         else if(inpts.eq.10) then
O
            jnpts= 32771*181
O
            if(ndim.eq.7) then
O
               vk(1)= 0.1d+01 "
O
               vk(2)= 0.4960373d+7"
O
               vk(3)= 0.4851623d+7"
O
               vk(4)= 0.3262017d+7!
O
               vk(5)= 0.722217d+6("
O
               vk(6)= 0.4644124d+7"
O
               vk(7)= 0.3212165d+7#
O
            else if(ndim.eq.9) then0
O
               vk(1)= 0.1d+01t"
O
               vk(2)= 0.2498015d+7"
O
               vk(3)= 0.4246511d+7"
O
               vk(4)= 0.4724489d+7"
O
               vk(5)= 0.3448063d+7"
O
               vk(6)= 0.1119927d+7"
O
               vk(7)= 0.2141959d+7!
O
               vk(8)= 0.115857d+7#!
O
               vk(9)= 0.287463d+70
O
            endifs"
O
         else if(inpts.eq.11) then
O
            jnpts= 32771*379
O
            if(ndim.eq.7) then
O
               vk(1)= 0.1d+01(#
O
               vk(2)= 0.9646626d+07g#
O
               vk(3)= 0.8128723d+07d#
O
               vk(4)= 0.9521278d+07s#
O
               vk(5)= 0.4720279d+076#
O
               vk(6)= 0.7036407d+07.#
O
               vk(7)= 0.1868554d+07 #
O
            else if(ndim.eq.9) thene
O
               vk(1)= 0.1d+01 #
O
               vk(2)= 0.4365962d+07c"
O
               vk(3)= 0.650771d+07#
O
               vk(4)= 0.6559665d+07 #
O
               vk(5)= 0.9838408d+07,#
O
               vk(6)= 0.2561851d+07t$
O
               vk(7)= 0.11842175d+08"
O
               vk(8)= 0.208211d+06#
O
               vk(9)= 0.6217272d+07=
O
            endifr
O
         endif
O
      endif0

O
*r
O
      ik= 0t
O
      do i=1,51
O

O
          ifz(i)= 0(
O
      enddo2

O
*m 
O
*-----Not available for changing

O
*m
O
      otrans= 'n'q

O
*b'
O
*-----4(8) xs and 4(8) max are computed
O


O
* (
O
*-----begin generation-evaluation branch

O
* 
O
      if(om.eq.'g') then
O
       if(osm.eq.'n') then

O
*s!
O
*-----all CC but CC12, XS is SIGAl

O
*=1
O
         if(oxcm.eq.'c'.and.otype.ne.'cc12') thenh
O
          do i0=1,2s
O
           do j0=1,2
O
            ix0= i0h
O
            it0= j0 
O
            ip= 2*(i0-1)+j0a
O
            pa0= 0.d0 
O
            pb0= 1.d0,
O
            np0= 9)
O
            call g05faf(pa0,pb0,np0,prx0),
O
            do i=1,npos
O

O
                xscmx(i)= 0.d0
O
            enddo,
O
            npts= jnptsc
O
            itrans= 02
O
            nrand= inrand 
O
            ifail= 0
O
            call g05cbf(0)<
O
            call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand,.
O
     #                  itrans,sig,esig,ifail)
O
            siga(ip)= sig.
O
            do i=1,npos=#
O
                xscmx0(i)= xscmx(i)m
O
            enddoc

O
*h 
O
*-----Improving upon the maximum

O
* 
O
            ibound= 0s
O
            do i=1,npos &
O
             call wtomxregion(i,bl,bu)
O
             istop= 0 
O
             do j=1,9)+
O
                 if((xmx(i,j).lt.bl(j)).or.c-
O
     #              (xmx(i,j).gt.bu(j))) then3"
O
                    istop= istop+1
O
                 endif
O
             enddo0
O
             if(istop.gt.0.or.ostop.eq.'s') then#
O
              xsmx(i,ip)= xscmx0(i)
O

O
             else0
O
              irst= 0m(
O
              if(xscmx0(i).eq.0.d0) then 
O
                 xscmx0(i)= 1.d0
O
                 do j=1,9+1
O
                     xmx0(j)= 0.5d0*(bl(j)+bu(j))a
O
                 enddo
O
              else
O
                 do j=1,9#&
O
                     xmx0(j)= xmx(i,j)
O
                 enddo
O
              endift
O
 1200         ifail= 1
O
              ireg= i*5
O
              call e04jaf(ndim,ibound,bl,bu,xmx0,fmx,h.
O
     #                    iw,liw,wmx,lw,ifail)!
O
              if(ifail.eq.2) thenr#
O
                 if(irst.eq.0) thenr
O
                    irst= 1 
O
                    go to 1200
O
                 elseg#
O
                    tmx1= xscmx0(i)
O
(
O
                    tmx2= -fmx*xscmx0(i)0
O
                    xsmx(i,ip)= dmax1(tmx1,tmx2)
O
                 endif&
O
              else if(ifail.gt.5) then&
O
                 xsmx(i,ip)= xscmx0(i)
O
              else+
O
                 xsmx(i,ip)= -fmx*xscmx0(i)=
O
              endifg
O
             endif
O
            enddo 

O
*o
O
            if(ip.eq.1) then%
O
               print 1999,ip,siga(ip) 2
O
               open(nin,file='dummy',status='new')
O
               do i=1,9 +
O
                   write(nin,fmt=*) prx0(i)s
O
               enddo
O
               do i=1,npos4
O
                write(nin,fmt=*) siga(ip),xsmx(i,ip)
O
               enddo
O
               close(nin) 
O
            else%
O
               print 1999,ip,siga(ip)=B
O
               open(nin,file='dummy',status='old',access='append')
O
               do i=1,npos4
O
                write(nin,fmt=*) siga(ip),xsmx(i,ip)
O
               enddo
O
               close(nin) 
O
            endifs
O
           enddo
O
          enddo/1
O
         else if(oxcm.eq.'n'.or.(oxcm.eq.'c'.and. 6
O
     #           otype.eq.'cc12').or.oxcm.eq.'m') then
O
            pa0= 0.d0 
O
            pb0= 1.d0*
O
            np0= 9)
O
            call g05faf(pa0,pb0,np0,prx0)+

O
*p(
O
*----- NC19/NC24/NC48/NC64, XS is SIGAN 

O
*
O
5
O
            if(otype.eq.'nc24'.or.otype.eq.'nc19'.or.e7
O
     #         otype.eq.'nc48'.or.otype.eq.'nc64') then 
O
            do i0=1,2f
O
            do j0=1,2f
O
             ix0= i0
O
             it0= j0
O
             ip= 2*(i0-1)+j0
O
             do i=1,npos
O
                 xshmx(i)= 0.d0 
O
             enddo
O
             npts= jnpts
O
             itrans= 0
O
             nrand= inrand
O
             ifail= 0 %
O
             if(otype.eq.'nc64') then*
O
              do j=1,npos
O

O
                  ifl(j)= 1 
O
              enddo
O

O
              call g05cbf(0)
O
              iint= 0 @
O
              call d01gcf(ndim,wtoxsn64,wtoregion,npts,vk,nrand,0
O
     #                    itrans,sig,esig,ifail)
O
             else
O

O
              call g05cbf(0)>
O
              call d01gcf(ndim,wtoxsn,wtoregion,npts,vk,nrand,0
O
     #                    itrans,sig,esig,ifail)
O
             endif
O
             sigan(ip)= sig 
O
*     "
O
             if(ostop.eq.'s') then
O
              do j=1,npos 
O
                  do i=1,ifl(j)l)
O
                      axstr(i)= stry(j,i)2
O
                  enddo2
O
                  ifail= 0
O
                  mm1= 1
O
                  mm2= ifl(j) 
O
                  ord= 'a'9
O
                  call m01daf(axstr,mm1,mm2,ord,ir,ifail)r
O
                  do i=1,ifl(j)v-
O
                      stry(j,ir(i))= axstr(i)r
O
                  enddo7&
O
                  if(ifl(j).eq.2) then.
O
                     xshmx0(j)= stry(j,ifl(j))+
O
                  else if(ifl(j).eq.1) then1$
O
                     xshmx0(j)= 0.d0
O
                  else!
O
                     ifmx= ifl(j)r"
O
                     ifm= ifl(j)-1&
O
 1133                if(ifm.eq.1) then/
O
                        xshmx0(j)= stry(j,ifmx))
O
                     elsee!
O
                        avr= 0.d0
O
"
O
                        avr2= 0.d0"
O
                        do i=1,ifm.
O
                            avr= avr+stry(j,i):
O
                            avr2= avr2+stry(j,i)*stry(j,i)
O
                        enddo*$
O
                        avr= avr/ifm&
O
                        avr2= avr2/ifm;
O
                        std2= ifm/(ifm-1.d0)*(avr2-avr*avr)q'
O
                        std= sqrt(std2) 3
O
                        dev= (stry(j,ifmx)-avr)/std ,
O
                        if(dev.lt.6.d0) then#
O
                           isort= 1 
O
                        else#
O
                           isort= 02
O
                        endif +
O
                        if(isort.eq.0) then %
O
                           ifm= ifm-1 '
O
                           ifmx= ifmx-10%
O
                           go to 11334
O
                        else1
O
                           xshmx0(j)= stry(j,ifm)d
O
                        endif 
O
                     endif
O
                  endif*
O
              enddo(
O
             else
O

O
              do i=1,npos
O
%
O
                  xshmx0(i)= xshmx(i)c
O
              enddo
O

O
             endif

O
*1
O
             ibound= 0
O
             do i=1,npos'
O
              call wtomxregion(i,bl,bu)*
O
              istop= 0#
O
              if(ostop.eq.'i') then(
O
               do j=1,9
O
.
O
                   if((xmxh(i,j).lt.bl(j)).or.0
O
     #                (xmxh(i,j).gt.bu(j))) then$
O
                      istop= istop+1
O
                   endif
O
               enddo
O
              endif
O
1
O
              if(istop.gt.0.or.ostop.eq.'s') thenr%
O
               xsmxn(i,ip)= xshmx0(i)p
O
              else
O
               irst= 0)
O
               if(xshmx0(i).eq.0.d0) then !
O
                  xshmx0(i)= 1.d0 
O
                  do j=1,92
O
                      xmx0(j)= 0.5d0*(bl(j)+bu(j))
O
                  enddo 
O
               elsep
O
                  do j=1,9(
O
                      xmx0(j)= xmxh(i,j)
O
                  enddo=
O
               endif
O
 1930          ifail= 1f
O
               ireg= i6
O
               call e04jaf(ndim,ibound,bl,bu,xmx0,fmx,/
O
     #                     iw,liw,wmx,lw,ifail)6"
O
               if(ifail.eq.2) then$
O
                  if(irst.eq.0) then
O
                     irst= 1
O
                     go to 1930-
O
                  else$
O
                     tmx1= xshmx0(i))
O
                     tmx2= -fmx*xshmx0(i) 2
O
                     xsmxn(i,ip)= dmax1(tmx1,tmx2)
O
                  endifw'
O
               else if(ifail.gt.5) thena(
O
                  xsmxn(i,ip)= xshmx0(i)
O
               else -
O
                  xsmxn(i,ip)= -fmx*xshmx0(i)d
O
               endif
O
              endifa
O
             enddo

O
*h
O
             if(ip.eq.1) then '
O
                print 1999,ip,sigan(ip)*3
O
                open(nin,file='dummy',status='new')n
O
                do i=1,9,
O
                    write(nin,fmt=*) prx0(i)
O
                enddox
O
                do i=1,nposn7
O
                 write(nin,fmt=*) sigan(ip),xsmxn(i,ip) 
O
                enddo 
O
                close(nin)
O
             else.'
O
                print 1999,ip,sigan(ip)fC
O
                open(nin,file='dummy',status='old',access='append')(
O
                do i=1,nposg7
O
                 write(nin,fmt=*) sigan(ip),xsmxn(i,ip) 
O
                enddof
O
                close(nin)
O
             endif
O
            enddo 
O
            enddom

O
*
O

O
*---- NC68, XS is SIGAHB

O
*t)
O
            else if(otype.eq.'nc68') then*

O
**"
O
            if(omssm.eq.'n') then 
O
               k0mx= 5&
O
            else if(omssm.eq.'y') then
O
               k0mx= 9
O
            endif.
O
            do i0=1,2 
O
            do j0=1,2)
O
            do k0=1,k0mx
O
             ix0= i0
O
             it0= j0
O
             ip0= k0(
O
             ip= k0mx*(2*(i0-1)+j0-1)+k0
O
             do i=1,npos
O
                 xshmx(i)= 0.d0*
O
             enddo
O
             npts= jnpts
O
             itrans= 0
O
             nrand= inrand
O
             ifail= 0t
O
             do j=1,npos
O
                 ifl(j)= 1
O
             enddo
O
             call g05cbf(0)*"
O
             if(omssm.eq.'n') thenA
O
              call d01gcf(ndim,wtoxsnh64,wtoregion,npts,vk,nrand,(0
O
     #                    itrans,sig,esig,ifail)'
O
             else if(omssm.eq.'y') thenfB
O
              call d01gcf(ndim,wtoxsnsh64,wtoregion,npts,vk,nrand,0
O
     #                    itrans,sig,esig,ifail)
O
             endif
O
             sigahb(ip)= sig
O
*     "
O
             if(ostop.eq.'s') then
O
              do j=1,nposl
O
                  do i=1,ifl(j) )
O
                      axstr(i)= stry(j,i)u
O
                  enddon
O
                  ifail= 0
O
                  mm1= 1
O
                  mm2= ifl(j)/
O
                  ord= 'a'9
O
                  call m01daf(axstr,mm1,mm2,ord,ir,ifail)s
O
                  do i=1,ifl(j)n-
O
                      stry(j,ir(i))= axstr(i)
O

O
                  enddo &
O
                  if(ifl(j).eq.2) then.
O
                     xshmx0(j)= stry(j,ifl(j))+
O
                  else if(ifl(j).eq.1) then $
O
                     xshmx0(j)= 0.d0
O
                  else!
O
                     ifmx= ifl(j)s"
O
                     ifm= ifl(j)-1&
O
 1233                if(ifm.eq.1) then/
O
                        xshmx0(j)= stry(j,ifmx) 
O
                     else(!
O
                        avr= 0.d0 "
O
                        avr2= 0.d0"
O
                        do i=1,ifm.
O
                            avr= avr+stry(j,i):
O
                            avr2= avr2+stry(j,i)*stry(j,i)
O
                        enddoi$
O
                        avr= avr/ifm&
O
                        avr2= avr2/ifm;
O
                        std2= ifm/(ifm-1.d0)*(avr2-avr*avr)m'
O
                        std= sqrt(std2)
O
3
O
                        dev= (stry(j,ifmx)-avr)/std2,
O
                        if(dev.lt.6.d0) then#
O
                           isort= 1p
O
                        else#
O
                           isort= 00
O
                        endife+
O
                        if(isort.eq.0) then(%
O
                           ifm= ifm-1p'
O
                           ifmx= ifmx-1 %
O
                           go to 1233 
O
                        else1
O
                           xshmx0(j)= stry(j,ifm)6
O
                        endifd
O
                     endif
O
                  endif 
O
              enddo6
O
             else(
O
              do i=1,npos %
O
                  xshmx0(i)= xshmx(i) 
O
              enddo 
O
             endif

O
*+
O
             ibound= 0
O
             do i=1,npos'
O
              call wtomxregion(i,bl,bu)v
O
              istop= 0#
O
              if(ostop.eq.'i') then 
O
               do j=1,9+.
O
                   if((xmxh(i,j).lt.bl(j)).or.0
O
     #                (xmxh(i,j).gt.bu(j))) then$
O
                      istop= istop+1
O
                   endif
O
               enddo
O
              endif)1
O
              if(istop.gt.0.or.ostop.eq.'s') then6&
O
               xsmxhb(i,ip)= xshmx0(i)
O
              else
O
               irst= 0)
O
               if(xshmx0(i).eq.0.d0) thenl!
O
                  xshmx0(i)= 1.d0 
O
                  do j=1,93
O
                      xmxh0(j)= 0.5d0*(bl(j)+bu(j))4
O
                  enddo(
O
               else 
O
                  do j=1,9)
O
                      xmxh0(j)= xmxh(i,j)
O

O
                  enddo6
O
               endif
O
 6901          ifail= 1
O

O
               ireg= i8
O
               call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh,/
O
     #                     iw,liw,wmx,lw,ifail)e"
O
               if(ifail.eq.2) then$
O
                  if(irst.eq.0) then
O
                     irst= 1
O
                     go to 6901f
O
                  else$
O
                     tmx1= xshmx0(i)*
O
                     tmx2= -fmxh*xshmx0(i)3
O
                     xsmxhb(i,ip)= dmax1(tmx1,tmx2) 
O
                  endif/'
O
               else if(ifail.gt.5) then )
O
                  xsmxhb(i,ip)= xshmx0(i)l
O
               else1/
O
                  xsmxhb(i,ip)= -fmxh*xshmx0(i)(
O
               endif
O
              endif0
O
             enddo
O
             if(ip.eq.1) then (
O
                print 1999,ip,sigahb(ip)3
O
                open(nin,file='dummy',status='new') 
O
                do i=1,9,
O
                    write(nin,fmt=*) prx0(i)
O
                enddos
O
                do i=1,npos 9
O
                 write(nin,fmt=*) sigahb(ip),xsmxhb(i,ip) 
O
                enddo1
O
                close(nin)
O
             else (
O
                print 1999,ip,sigahb(ip)C
O
                open(nin,file='dummy',status='old',access='append')v
O
                do i=1,npos 9
O
                 write(nin,fmt=*) sigahb(ip),xsmxhb(i,ip)3
O
                enddo 
O
                close(nin)
O
             endif
O
            enddo 
O
            enddo 
O
            enddo4

O
*5
O
*---- NC26/NC33, XS is SIGAHQ5

O
*4<
O
            else if(otype.eq.'nc33'.or.otype.eq.'nc26') then

O
*.
O
            do i0=1,2 
O
            do j0=1,2
O

O
            do k0=1,35
O
             ix0= i0
O
             it0= j0
O
             ip0= k0%
O
             ip= 3*(2*(i0-1)+j0-1)+k0 
O
             do i=1,npos
O
                 xshmx(i)= 0.d0 
O
             enddo
O
             npts= jnpts
O
             itrans= 0
O
             nrand= inrand
O
             ifail= 0 
O
             call g05cbf(0)v%
O
             if(otype.eq.'nc33') then)A
O
              call d01gcf(ndim,wtoxsnh32,wtoregion,npts,vk,nrand, 0
O
     #                    itrans,sig,esig,ifail)*
O
             else if(otype.eq.'nc26') thenA
O
              call d01gcf(ndim,wtoxsnh26,wtoregion,npts,vk,nrand,
O
0
O
     #                    itrans,sig,esig,ifail)
O
             endif
O
             sigahq(ip)= sig
O
             do i=1,npos$
O
                 xshmx0(i)= xshmx(i)
O
             enddo
O
             ibound= 0
O
             do i=1,npos'
O
              call wtomxregion(i,bl,bu)v
O
              istop= 0
O
              do j=1,9-
O
                  if((xmxh(i,j).lt.bl(j)).or. /
O
     #               (xmxh(i,j).gt.bu(j))) then #
O
                     istop= istop+1.
O
                  endif 
O
              enddo
O
1
O
              if(istop.gt.0.or.ostop.eq.'s') then &
O
               xsmxhq(i,ip)= xshmx0(i)
O
              else
O
               irst= 0)
O
               if(xshmx0(i).eq.0.d0) then !
O
                  xshmx0(i)= 1.d0 
O
                  do j=1,93
O
                      xmxh0(j)= 0.5d0*(bl(j)+bu(j)) 
O
                  enddo 
O
               else.
O
                  do j=1
)
O
                      xmxh0(j)= xmxh(i,j) 
O
                  enddo
O

O
               endif
O
 1901          ifail= 1v
O
               ireg= i8
O
               call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh,/
O
     #                     iw,liw,wmx,lw,ifail) "
O
               if(ifail.eq.2) then$
O
                  if(irst.eq.0) then
O
                     irst= 1
O
                     go to 1901e
O
                  else$
O
                     tmx1= xshmx0(i)*
O
                     tmx2= -fmxh*xshmx0(i)3
O
                     xsmxhq(i,ip)= dmax1(tmx1,tmx2) 
O
                  endif0'
O
               else if(ifail.gt.5) then0)
O
                  xsmxhq(i,ip)= xshmx0(i) 
O
               else./
O
                  xsmxhq(i,ip)= -fmxh*xshmx0(i)
O

O
               endif
O
              endifn
O
             enddo
O
             if(ip.eq.1) then (
O
                print 1999,ip,sigahq(ip)3
O
                open(nin,file='dummy',status='new')-
O
                do i=1,9,
O
                    write(nin,fmt=*) prx0(i)
O
                enddof
O
                do i=1,nposs9
O
                 write(nin,fmt=*) sigahq(ip),xsmxhq(i,ip)
O

O
                enddo.
O
                close(nin)
O
             elses(
O
                print 1999,ip,sigahq(ip)C
O
                open(nin,file='dummy',status='old',access='append')a
O
                do i=1,nposd9
O
                 write(nin,fmt=*) sigahq(ip),xsmxhq(i,ip),
O
                enddoi
O
                close(nin)
O
             endif
O
            enddo 
O
            enddo
O

O
            enddo 

O
*
O

O
*---- NC50, XS is SIGAH 

O
* )
O
            else if(otype.eq.'nc50') then(

O
*
O

O
            do i0=1,2(
O
            do j0=1,2p
O
            do k0=1,2 
O
             ix0= i0
O
             it0= j0
O
             ip0= k0%
O
             ip= 2*(2*(i0-1)+j0-1)+k0 
O
             do i=1,npos
O
                 xshmx(i)= 0.d0v
O
             enddo
O
             npts= jnpts
O
             itrans= 0
O
             nrand= inrand
O
             ifail= 0 
O
             call g05cbf(0) @
O
             call d01gcf(ndim,wtoxsnh49,wtoregion,npts,vk,nrand,/
O
     #                   itrans,sig,esig,ifail) 
O
             sigah(ip)= sig 
O
             do i=1,npos$
O
                 xshmx0(i)= xshmx(i)
O
             enddo
O
             ibound= 0
O
             do i=1,npos'
O
              call wtomxregion(i,bl,bu) 
O
              istop= 0
O
              do j=1,9-
O
                  if((xmxh(i,j).lt.bl(j)).or.,/
O
     #               (xmxh(i,j).gt.bu(j))) then(#
O
                     istop= istop+1 
O
                  endif 
O
              enddo 1
O
              if(istop.gt.0.or.ostop.eq.'s') thend%
O
               xsmxh(i,ip)= xshmx0(i) 
O
              else
O
               irst= 0)
O
               if(xshmx0(i).eq.0.d0) thenh!
O
                  xshmx0(i)= 1.d0w
O
                  do j=1,93
O
                      xmxh0(j)= 0.5d0*(bl(j)+bu(j))t
O
                  enddo 
O
               else 
O
                  do j=1,9)
O
                      xmxh0(j)= xmxh(i,j)s
O
                  enddo 
O
               endif
O
 1201          ifail= 1s
O
               ireg= i8
O
               call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh,/
O
     #                     iw,liw,wmx,lw,ifail)
O
"
O
               if(ifail.eq.2) then$
O
                  if(irst.eq.0) then
O
                     irst= 1
O
                     go to 1201 
O
                  else$
O
                     tmx1= xshmx0(i)*
O
                     tmx2= -fmxh*xshmx0(i)2
O
                     xsmxh(i,ip)= dmax1(tmx1,tmx2)
O
                  endifi'
O
               else if(ifail.gt.5) thend(
O
                  xsmxh(i,ip)= xshmx0(i)
O
               else).
O
                  xsmxh(i,ip)= -fmxh*xshmx0(i)
O
               endif
O
              endif1
O
             enddo
O
             if(ip.eq.1) thenm'
O
                print 1999,ip,sigah(ip) 3
O
                open(nin,file='dummy',status='new')s
O
                do i=1,9,
O
                    write(nin,fmt=*) prx0(i)
O
                enddo 
O
                do i=1,npos
O
7
O
                 write(nin,fmt=*) sigah(ip),xsmxh(i,ip)#
O
                enddo1
O
                close(nin)
O
             else '
O
                print 1999,ip,sigah(ip)=C
O
                open(nin,file='dummy',status='old',access='append')N
O
                do i=1,npos 7
O
                 write(nin,fmt=*) sigah(ip),xsmxh(i,ip) 
O
                enddo4
O
                close(nin)
O
             endif
O
            enddoo
O
            enddo 
O
            enddo 

O
* &
O
*-----CC12/NC25/NC32/MX43, XS is SIGAH

O
* :
O
            else if(otype.eq.'nc25'.or.otype.eq.'nc32'.or.<
O
     #              otype.eq.'cc12'.or.otype.eq.'mx43') then

O
*
O
$
O
            if(otype.eq.'mx43') then#
O
               if(ofs.eq.'ll') thenq
O
                  k0mx= 1 9
O
               else if(ofs.eq.'qq'.and.opglu.eq.'n') then 
O
                  k0mx= 1g9
O
               else if(ofs.eq.'qq'.and.opglu.eq.'y') then(
O
                  k0mx= 2s
O
               endif)
O
            else if(otype.eq.'nc32') then $
O
               if(opglu.eq.'n') then
O
                  k0mx= 1d)
O
               else if(opglu.eq.'y') then
O

O
                  k0mx= 2t
O
               endif
O
            else
O
               k0mx= 2
O
            endif 

O
* 
O
            do i0=1,2h
O
            do j0=1,21
O
            do k0=1,k0mx
O
             ix0= i0
O
             it0= j0
O
             ip0= k0(
O
             ip= k0mx*(2*(i0-1)+j0-1)+k0
O
             do i=1,npos
O
                 xshmx(i)= 0.d0j
O
             enddo
O
             npts= jnpts
O
             itrans= 0
O
             nrand= inrand
O
             ifail= 0j%
O
             if(otype.eq.'nc25') then=
O
              call g05cbf(0)A
O
              call d01gcf(ndim,wtoxsnh24,wtoregion,npts,vk,nrand, 0
O
     #                    itrans,sig,esig,ifail)*
O
             else if(otype.eq.'nc32') then
O
              call g05cbf(0)@
O
              call d01gcf(ndim,wtoxsn32,wtoregion,npts,vk,nrand,0
O
     #                    itrans,sig,esig,ifail)*
O
             else if(otype.eq.'cc12') then
O
              call g05cbf(0)?
O
              call d01gcf(ndim,wtoxssc,wtoregion,npts,vk,nrand, 0
O
     #                    itrans,sig,esig,ifail)*
O
             else if(otype.eq.'mx43') then
O
               iint= 0
O
               call g05cbf(0))A
O
               call d01gcf(ndim,wtoxsm43,wtoregion,npts,vk,nrand, 1
O
     #                     itrans,sig,esig,ifail)m
O
             endif
O
             sigah(ip)= sigv
O
             do i=1,npos$
O
                 xshmx0(i)= xshmx(i)
O
             enddo
O
             ibound= 0

O
* 
O
             do i=1,npos'
O
              call wtomxregion(i,bl,bu) 
O
              istop= 0
O
              do j=1,9-
O
                  if((xmxh(i,j).lt.bl(j)).or. /
O
     #               (xmxh(i,j).gt.bu(j))) thent#
O
                     istop= istop+1 
O
                  endif 
O
              enddo=1
O
              if(istop.gt.0.or.ostop.eq.'s') then %
O
               xsmxh(i,ip)= xshmx0(i) 
O
              else
O
               irst= 0)
O
               if(xshmx0(i).eq.0.d0) then !
O
                  xshmx0(i)= 1.d0*
O
                  do j=1,93
O
                      xmxh0(j)= 0.5d0*(bl(j)+bu(j)) 
O
                  enddoi
O
               elsed
O
                  do j=1,9)
O
                      xmxh0(j)= xmxh(i,j),
O
                  enddow
O
               endif
O
 1202          ifail= 1 
O
               ireg= i8
O
               call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh,/
O
     #                     iw,liw,wmx,lw,ifail) "
O
               if(ifail.eq.2) then$
O
                  if(irst.eq.0) then
O
                     irst= 1
O
                     go to 1202 
O
                  else$
O
                     tmx1= xshmx0(i)*
O
                     tmx2= -fmxh*xshmx0(i)2
O
                     xsmxh(i,ip)= dmax1(tmx1,tmx2)
O
                  endif '
O
               else if(ifail.gt.5) then (
O
                  xsmxh(i,ip)= xshmx0(i)
O
               elsej.
O
                  xsmxh(i,ip)= -fmxh*xshmx0(i)
O
               endif
O
              endif 
O
             enddo
O
             if(ip.eq.1) thend'
O
                print 1999,ip,sigah(ip)i3
O
                open(nin,file='dummy',status='new') 
O
                do i=1,9,
O
                    write(nin,fmt=*) prx0(i)
O
                enddo 
O
                do i=1,npos 7
O
                 write(nin,fmt=*) sigah(ip),xsmxh(i,ip)s
O
                enddo 
O
                close(nin)
O
             else '
O
                print 1999,ip,sigah(ip) C
O
                open(nin,file='dummy',status='old',access='append')x
O
                do i=1,npos 7
O
                 write(nin,fmt=*) sigah(ip),xsmxh(i,ip) 
O
                enddo=
O
                close(nin)
O
             endif
O
            enddom
O
            enddo 
O
            enddo 

O
* 
O
*---- NC21, XS is SIGAH 

O
*d)
O
            else if(otype.eq.'nc21') then 

O
* 
O
            do i0=1,2i
O
            do j0=1,2 
O
            do k0=1,2,
O
             ix0= i0
O
             it0= j0
O
             ip0= k0%
O
             ip= 2*(2*(i0-1)+j0-1)+k0n
O
             do i=1,npos
O
                 xshmx(i)= 0.d0n
O
             enddo
O
             npts= jnpts
O
             itrans= 0
O
             nrand= inrand
O
             ifail= 0 
O
             call g05cbf(0) @
O
             call d01gcf(ndim,wtoxsnh19,wtoregion,npts,vk,nrand,/
O
     #                   itrans,sig,esig,ifail)i
O
             sigah(ip)= sig(
O
             do i=1,npos$
O
                 xshmx0(i)= xshmx(i)
O
             enddo
O
             ibound= 0
O
             do i=1,npos'
O
              call wtomxregion(i,bl,bu)t
O
              istop= 0
O
              do j=1,9-
O
                  if((xmxh(i,j).lt.bl(j)).or. /
O
     #               (xmxh(i,j).gt.bu(j))) then=#
O
                     istop= istop+1o
O
                  endifj
O
              enddoo1
O
              if(istop.gt.0.or.ostop.eq.'s') then %
O
               xsmxh(i,ip)= xshmx0(i)i
O
              else
O
               irst= 0)
O
               if(xshmx0(i).eq.0.d0) then
O
!
O
                  xshmx0(i)= 1.d0t
O
                  do j=1,93
O
                      xmxh0(j)= 0.5d0*(bl(j)+bu(j)) 
O
                  enddop
O
               elsei
O
                  do j=1,9)
O
                      xmxh0(j)= xmxh(i,j) 
O
                  enddo 
O
               endif
O
 1203          ifail= 1n
O
               ireg= i8
O
               call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh,/
O
     #                     iw,liw,wmx,lw,ifail)1"
O
               if(ifail.eq.2) then$
O
                  if(irst.eq.0) then
O
                     irst= 1
O
                     go to 1203p
O
                  else$
O
                     tmx1= xshmx0(i)*
O
                     tmx2= -fmxh*xshmx0(i)2
O
                     xsmxh(i,ip)= dmax1(tmx1,tmx2)
O
                  endif
O
'
O
               else if(ifail.gt.5) then (
O
                  xsmxh(i,ip)= xshmx0(i)
O
               elser.
O
                  xsmxh(i,ip)= -fmxh*xshmx0(i)
O
               endif
O
              endif1
O
             enddo
O
             if(ip.eq.1) thenr'
O
                print 1999,ip,sigah(ip) 3
O
                open(nin,file='dummy',status='new')h
O
                do i=1,9,
O
                    write(nin,fmt=*) prx0(i)
O
                enddox
O
                do i=1,npos 7
O
                 write(nin,fmt=*) sigah(ip),xsmxh(i,ip) 
O
                enddo3
O
                close(nin)
O
             else '
O
                print 1999,ip,sigah(ip) C
O
                open(nin,file='dummy',status='old',access='append') 
O
                do i=1,npos 7
O
                 write(nin,fmt=*) sigah(ip),xsmxh(i,ip) 
O
                enddo 
O
                close(nin)
O
             endif
O
            enddo 
O
            enddoi
O
            enddo 
O
            endif
O

O
         endif

O
* 
O
       else if(osm.eq.'g') then !
O
          if(itmx.gt.itmxmx) then2E
O
             print*,' Number of events exceedes the maximum allowed ' 
O
             stopi
O
          endif 
O
         if(oxcm.eq.'c') then./
O
            open(nin,file='dummy',status='old')
O

O
            do i=1,9(
O
                read(nin,fmt=*) prx0n(i)
O
            enddo 
O
            do ip=1,4 
O
             do i=1,npos1
O
              read(nin,fmt=*) siga(ip),xsmx(i,ip)
O

O
             enddo
O
            enddo 
O
            close(nin)
O
            nint= 4 A
O
            call wtovolume(npos,nint,xsmx,oxsmx,gbl,gbu,vol,volt)y1
O
         else if(oxcm.eq.'n'.or.oxcm.eq.'m') then 5
O
            if(otype.eq.'nc24'.or.otype.eq.'nc19'.or. 7
O
     #         otype.eq.'nc48'.or.otype.eq.'nc64') then
O
0
O
             open(nin,file='dummy',status='old')
O
             do i=1,9 )
O
                 read(nin,fmt=*) prx0n(i) 
O
             enddo
O
             do ip=1,4
O
              do i=1,npos 4
O
               read(nin,fmt=*) sigan(ip),xsmxn(i,ip)
O
              enddo
O

O
             enddo
O
             close(nin).
O
             nint= 4D
O
             call wtovolume(npos,nint,xsmxn,oxsmxn,gbl,gbu,vol,volt)
O
            else%
O
             if(otype.eq.'nc68') then 1
O
              open(nin,file='dummy',status='old')s
O
              do i=1,9*
O
                  read(nin,fmt=*) prx0n(i)
O
              enddo=#
O
              if(omssm.eq.'n') thend
O
                 ipmx= 20 (
O
              else if(omssm.eq.'y') then
O
                 ipmx= 36m
O
              endifb
O
              do ip=1,ipmx
O
               do i=1,npos7
O
                read(nin,fmt=*) sigahb(ip),xsmxhb(i,ip))
O
               enddo
O
              enddo 
O
              close(nin)=
O
             else if(otype.eq.'nc33'.or.otype.eq.'nc26') then01
O
              open(nin,file='dummy',status='old') 
O
              do i=1,9*
O
                  read(nin,fmt=*) prx0n(i)
O
              enddof
O
              do ip=1,12
O
               do i=1,npos7
O
                read(nin,fmt=*) sigahq(ip),xsmxhq(i,ip) 
O
               enddo
O
              enddo 
O
              close(nin)
O
             else 1
O
              open(nin,file='dummy',status='old') 
O
              do i=1,9*
O
                  read(nin,fmt=*) prx0n(i)
O
              enddos&
O
              if(otype.eq.'mx43') then%
O
                 if(ofs.eq.'ll') then(
O
                    ipmx= 4 ;
O
                 else if(ofs.eq.'qq'.and.opglu.eq.'n') thent
O
                    ipmx= 49;
O
                 else if(ofs.eq.'qq'.and.opglu.eq.'y') thens
O
                    ipmx= 8i
O
                 endif+
O
              else if(otype.eq.'nc32') thenn&
O
                 if(opglu.eq.'n') then
O
                    ipmx= 4 +
O
                 else if(opglu.eq.'y') thend
O
                    ipmx= 8)
O
                 endif
O
              else
O
                 ipmx= 8
O
              endifl

O
*d
O
              do ip=1,ipmx
O
               do i=1,npos5
O
                read(nin,fmt=*) sigah(ip),xsmxh(i,ip))
O
               enddo
O
              enddo 
O

            close(nin)
O
             endif%
O
             if(otype.eq.'nc68') then 
O
              nint= ipmxB
O
              call wtovolume(npos,nint,xsmxhb,oxsmxhb,gblhb,gbuhb,'
O
     #                    volhb,volhtb) =
O
             else if(otype.eq.'nc33'.or.otype.eq.'nc26') then 
O
              nint= 12B
O
              call wtovolume(npos,nint,xsmxhq,oxsmxhq,gblhq,gbuhq,'
O
     #                    volhq,volhtq) 
O
             else=
O
              nint= ipmxC
O
              call wtovolume(npos,nint,xsmxh,oxsmxh,gblh,gbuh,volh,d 
O
     #                    volht)
O
             endif
O
            endif 
O
         endif

O
*3
O
*-----generation for CC starts

O
*(
O
         nlim= npos-1n
O
         do i= 1,itmx#
O
             amps1(i)= 0.d0g
O
             amps2(i)= 0.d0s
O
         enddo
O
         if(oxcm.eq.'c') thend

O
*c
O
         it= 1.
O
         sigt= siga(1)+siga(2)+siga(3)+siga(4)'
O
         sigt2= siga(2)+siga(3)+siga(4) 
O
         sigt3= siga(3)+siga(4)s
O
         tsig= sigt=
O
 1504    wxs= g05caf(dmy)h
O
         wxs0= siga(1)/sigt 
O
         if(wxs.lt.wxs0) then
O

O
          ix0= 1
O
          it0= 1
O
 1500     jp= 1i
O
          volp= volt(1) 
O
 3000     test0= vol(jp,1)/volp9
O
          pru0= g05caf(dmy)i 
O
          if(pru0.lt.test0) then
O
           do l=1,9b 
O
            aprx(l)= g05caf(dmy)6
O
            prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+ 
O
     #               gbl(jp,l,1)
O
           enddo
O
          else
O
           if(jp.lt.nlim) then 
O
            volp= volp-vol(jp,1)
O
            jp= jp+1
O
            go to 3000     e
O
           else 
O
            jp= jp+1
O
            do l=1,9!
O
             aprx(l)= g05caf(dmy))7
O
             prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+ !
O
     #                gbl(jp,l,1)j
O
            enddo 
O
           endif
O
          endif 
O
          pru= g05caf(pru)
O
          xt0= wtoxsc(ndim,prx)l
O
          xt= xt0/siga(1)i"
O
          prh= oxsmx(jp,1)/siga(1)
O
          prtst= xt-pru*prh#/
O
          if(prtst.ge.0.d0.and.xt.ne.0.d0) then 
O
           do i=1,ninv
O
            prxa(i,it)= xap(i)
O
           enddo
O
           amps(it)= xt0
O
           it= it+1 
O
          else  
O
            go to 1500
O
          endif=
         else 
O
          wxs2= g05caf(dmy)x
O
          wxs0= siga(2)/sigt2x
O
          if(wxs2.lt.wxs0) then 
O
             ix0= 10
O
             it0= 2 
O
 1501        jp= 1
O
             volp= volt(2)"
O
 3001        test0= vol(jp,2)/volp
O
             pru0= g05caf(dmy)#
O
             if(pru0.lt.test0) thene
O
              do l=1,9#
O
               aprx(l)= g05caf(dmy)i9
O
               prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+
O
#
O
     #                  gbl(jp,l,2),
O
              enddo 
O
             else
O
!
O
              if(jp.lt.nlim) then)#
O
               volp= volp-vol(jp,2) 
O
               jp= jp+1
O

O
               go to 3001     
O
              else
O
               jp= jp+1n
O
               do l=1,9s$
O
                aprx(l)= g05caf(dmy):
O
                prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+$
O
     #                   gbl(jp,l,2)
O
               enddo
O
              endif 
O
             endif
O
             pru= g05caf(pru) "
O
             xt0= wtoxsc(ndim,prx)
O
             xt= xt0/siga(2)%
O
             prh= oxsmx(jp,2)/siga(2)d
O
             prtst= xt-pru*prh2
O
             if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O
              do i=1,ninv !
O
               prxa(i,it)= xap(i)=
O
              enddok
O
              amps(it)= xt0
O

O
              it= it+1
O
             else 
O
              go to 1501
O
             endif
O
          else
O
             wxs3= g05caf(dmy) 
O
             wxs0= siga(3)/sigt3"
O
             if(wxs3.lt.wxs0) then
O
              ix0= 2
O
              it0= 1
O
 1502         jp= 1 
O
              volp= volt(3) #
O
 3002         test0= vol(jp,3)/volpg
O
              pru0= g05caf(dmy) $
O
              if(pru0.lt.test0) then
O
               do l=1,9 $
O
                aprx(l)= g05caf(dmy):
O
                prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+$
O
     #                   gbl(jp,l,3)
O
               enddo
O
              else"
O
               if(jp.lt.nlim) then$
O
                volp= volp-vol(jp,3)
O
                jp= jp+1
O
                go to 3002     (
O
               else 
O
                jp= jp+1
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy)f;
O
                 prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+)%
O
     #                    gbl(jp,l,3) 
O
                enddo 
O
               endif
O
              endif 
O
              pru= g05caf(pru)#
O
              xt0= wtoxsc(ndim,prx) 
O
              xt= xt0/siga(3))&
O
              prh= oxsmx(jp,3)/siga(3)
O
              prtst= xt-pru*prh 3
O
              if(prtst.ge.0.d0.and.xt.ne.0.d0) thens
O
               do i=1,ninv"
O
                prxa(i,it)= xap(i)
O
               enddo
O
               amps(it)= xt0
O
               it= it+1,
O
              else #
O
                go to 1502
O
              endif 
O
           else.
O
              ix0= 2
O
              it0= 2
O
 1503         jp= 1 
O
              volp= volt(4) #
O
 3003         test0= vol(jp,4)/volpl
O
              pru0= g05caf(dmy)h$
O
              if(pru0.lt.test0) then
O
               do l=1,9 $
O
                aprx(l)= g05caf(dmy):
O
                prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+$
O
     #                   gbl(jp,l,4)
O
               enddo
O
              else"
O
               if(jp.lt.nlim) then$
O
                volp= volp-vol(jp,4)
O
                jp= jp+1
O
                go to 3003     e
O
               elser
O
                jp= jp+1
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy) ;
O
                 prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+ %
O
     #                    gbl(jp,l,4)i
O
                enddo 
O
               endif
O
              endif
O

O
              pru= g05caf(pru)#
O
              xt0= wtoxsc(ndim,prx) 
O
              xt= xt0/siga(4),&
O
              prh= oxsmx(jp,4)/siga(4)
O
              prtst= xt-pru*prhp3
O
              if(prtst.ge.0.d0.and.xt.ne.0.d0) then 
O
               do i=1,ninv"
O
                prxa(i,it)= xap(i)
O
               enddo
O
               amps(it)= xt0
O
               it= it+1 
O
              else  
O
                go to 1503
O
              endifG
O
           endif
O
          endifq
O
         endif
O
         if(it.le.itmx) then
O
            go to 1504
O
         endif

O
*h(
O
*-----generation for Higgs/NC/Mix starts

O
*
O
1
O
         else if(oxcm.eq.'n'.or.oxcm.eq.'m') then 

O
* 2
O
         if(otype.eq.'nc24'.or.otype.eq.'nc19'.or.4
O
     #      otype.eq.'nc48'.or.otype.eq.'nc64') thenA
O
         if(otype.eq.'nc64'.and.ofs.eq.'qq'.and.opglu.eq.'n'.and. 
O
     #      oint.eq.'y') thenf
O
            iselsp= 1t
         elset
O
            iselsp= 0(
O
         endif
O
          it= 1 3
O
          sigt= sigan(1)+sigan(2)+sigan(3)+sigan(4)e+
O
          sigt2= sigan(2)+sigan(3)+sigan(4) "
O
          sigt3= sigan(3)+sigan(4)
O
          tsig= sigt
O
 1564     wxs= g05caf(dmy)
O
          wxs0= sigan(1)/sigt,
O
          if(wxs.lt.wxs0) then
O
           ix0= 1
O

O
           it0= 10
O
 1560      jp= 1
O
           volp= volt(1) 
O
 3060      test0= vol(jp,1)/volp
O
           pru0= g05caf(dmy)!
O
           if(pru0.lt.test0) theno
O
            do l=1,9!
O
             aprx(l)= g05caf(dmy) B
O
             prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+gbl(jp,l,1)
O
            enddo 
O
           elseg
O
            if(jp.lt.nlim) thenc!
O
             volp= volp-vol(jp,1)n
O
             jp= jp+1 
O
             go to 3060     
O
            else
O
             jp= jp+1 
O
             do l=1,9l"
O
              aprx(l)= g05caf(dmy)C
O
              prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+gbl(jp,l,1)t
O
             enddo
O
            endif 
O
           endif
O
           pru= g05caf(pru)5#
O
           if(otype.eq.'nc64') thenm
O
            iint= 0t#
O
            xt0= wtoxsn64(ndim,prx)  
O
            if(iselsp.eq.1) then
O
             iint= 1$
O
             xt1= wtoxsn64(ndim,prx)
O
             iint= 2$
O
             xt2= wtoxsn64(ndim,prx)
O
            endifi
O
            xt= xt0/sigan(1)
O
           elsei!
O
            xt0= wtoxsn(ndim,prx) 
O
            xt= xt0/sigan(1)
O
           endif%
O
           prh= oxsmxn(jp,1)/sigan(1)(
O
           prtst= xt-pru*prh0
O
           if(prtst.ge.0
0.and.xt.ne.0.d0) then
O
            do i=1,ninvl 
O
             prxa(i,it)= xaph(i)
O
            enddo 
O
            amps(it)= xt0  
O
            if(iselsp.eq.1) then
O
               amps1(it)= xt1i
O
               amps2(it)= xt2 
O
            endifs
O
            it= it+1
O
           else 
O
             go to 1560 
O
           endif
O
          else
O
           wxs2= g05caf(dmy)
O
           wxs0= sigan(2)/sigt2l 
O
           if(wxs2.lt.wxs0) then
O
              ix0= 1
O
              it0= 2
O
 1561         jp= 1=
O
              volp= volt(2)j#
O
 3061         test0= vol(jp,2)/volp 
O
              pru0= g05caf(dmy) $
O
              if(pru0.lt.test0) then
O
               do l=1,9
O
$
O
                aprx(l)= g05caf(dmy):
O
                prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+$
O
     #                   gbl(jp,l,2)
O
               enddo
O
              else"
O
               if(jp.lt.nlim) then$
O
                volp= volp-vol(jp,2)
O
                jp= jp+1
O
                go to 3061     .
O
               else 
O
                jp= jp+1
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy) ;
O
                 prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+x%
O
     #                    gbl(jp,l,2) 
O
                enddo 
O
               endif
O
              endifi
O
              pru= g05caf(pru)&
O
              if(otype.eq.'nc64') then
O
               iint= 0&
O
               xt0= wtoxsn64(ndim,prx)#
O
               if(iselsp.eq.1) then
O

O
                iint= 1 '
O
                xt1= wtoxsn64(ndim,prx) 
O
                iint= 2
O
'
O
                xt2= wtoxsn64(ndim,prx)t
O
               endif
O
               xt= xt0/sigan(2) 
O
              else$
O
               xt0= wtoxsn(ndim,prx)
O
               xt= xt0/sigan(2) 
O
              endif*(
O
              prh= oxsmxn(jp,2)/sigan(2)
O
              prtst= xt-pru*prh)3
O
              if(prtst.ge.0.d0.and.xt.ne.0.d0) theni
O
               do i=1,ninv#
O
                prxa(i,it)= xaph(i)s
O
               enddo
O
               amps(it)= xt0#
O
               if(iselsp.eq.1) thens 
O
                  amps1(it)= xt1 
O
                  amps2(it)= xt2
O
               endif
O
               it= it+1 
O
              else n
O
               go to 1561 
O
              endif 
O
           elseq
O
              wxs3= g05caf(dmy)o"
O
              wxs0= sigan(3)/sigt3#
O
              if(wxs3.lt.wxs0) then0
O
               ix0= 2 
O
               it0= 10
O
 1562          jp= 1
O
               volp= volt(3)$
O
 3062          test0= vol(jp,3)/volp 
O
               pru0= g05caf(dmy)%
O
               if(pru0.lt.test0) thenr
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy) ;
O
                 prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+1%
O
     #                    gbl(jp,l,3) 
O
                enddo,
O
               else #
O
                if(jp.lt.nlim) then %
O
                 volp= volp-vol(jp,3)s
O
                 jp= jp+1
O
 
O
                 go to 3062     
O
                else
O
                 jp= jp+1b
O
                 do l=1,9=&
O
                  aprx(l)= g05caf(dmy)<
O
                  prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+&
O
     #                     gbl(jp,l,3)
O
                 enddo
O
                endif
O

O
               endif
O
               pru= g05caf(pru)t'
O
               if(otype.eq.'nc64') theni
O
                iint= 0 '
O
                xt0= wtoxsn64(ndim,prx) $
O
                if(iselsp.eq.1) then
O
                 iint= 1(
O
                 xt1= wtoxsn64(ndim,prx)
O
                 iint= 2(
O
                 xt2= wtoxsn64(ndim,prx)
O
                endifs 
O
                xt= xt0/sigan(3)
O
               elsem%
O
                xt0= wtoxsn(ndim,prx)n 
O
                xt= xt0/sigan(3)
O
               endif)
O
               prh= oxsmxn(jp,3)/sigan(3)0 
O
               prtst= xt-pru*prh3
O
              if(prtst.ge.0.d0.and.xt.ne.0.d0) then 
O
                do i=1,ninvh$
O
                 prxa(i,it)= xaph(i)
O
                enddo 
O
                amps(it)= xt0g$
O
                if(iselsp.eq.1) then!
O
                   amps1(it)= xt1
O
!
O
                   amps2(it)= xt2h
O
                endif 
O
                it= it+1
O
               else 
O
                 go to 1562e
O
               endif
O
            else
O
               ix0= 2(
O
               it0= 2r
O
 1563          jp= 1
O
               volp= volt(4)$
O
 3063          test0= vol(jp,4)/volp 
O
               pru0= g05caf(dmy)%
O
               if(pru0.lt.test0) thenr
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy)t;
O
                 prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+t%
O
     #                    gbl(jp,l,4)n
O
                enddoo
O
               else #
O
                if(jp.lt.nlim) thenx%
O
                 volp= volp-vol(jp,4) 
O
                 jp= jp+1  
O
                 go to 3063     
O
                else
O
                 jp= jp+1l
O
                 do l=1,9s&
O
                  aprx(l)= g05caf(dmy)<
O
                  prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+&
O
     #                     gbl(jp,l,4)
O
                 enddo
O
                endifd
O
               endif
O
               pru= g05caf(pru) '
O
               if(otype.eq.'nc64') thenf
O
                iint= 0 '
O
                xt0= wtoxsn64(ndim,prx) $
O
                if(iselsp.eq.1) then
O
                 iint= 1(
O
                 xt1= wtoxsn64(ndim,prx)
O
                 iint= 2(
O
                 xt2= wtoxsn64(ndim,prx)
O
                endif1 
O
                xt= xt0/sigan(4)
O
               else %
O
                xt0= wtoxsn(ndim,prx)  
O
                xt= xt0/sigan(4)
O
               endif)
O
               prh= oxsmxn(jp,4)/sigan(4)  
O
               prtst= xt-pru*prh4
O
               if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O
                do i=1,ninv,$
O
                 prxa(i,it)= xaph(i)
O
                enddo 
O
                amps(it)= xt0.$
O
                if(iselsp.eq.1) then!
O
                   amps1(it)= xt1t!
O
                   amps2(it)= xt2u
O
                endif 
O
                it= it+1
O
               else 
O
                 go to 1563e
O
               endif
O
            endif 
O
           endif
O
          endif(
O
          if(it.le.itmx) then)
O
             go to 1564 
O
          endif
O


O
* &
O
         else if(otype.eq.'nc68') then
O
          it= 1c
O
          sigt= 0.d0
O
          do i=1,ipmxv!
O
             sigt= sigt+sigahb(i) 
O
          enddo'
O
          tsig= sigt
O
          sigp= sigt
O
 8517     do i0=1,2
O

O
           do j0=1,2
O
            do k0=1,k0mx
O
             ix0= i0
O
             it0= j0
O
             ip0= k0(
O
             ip= k0mx*(2*(i0-1)+j0-1)+k0
O
             wxs= g05caf(dmy)'"
O
             wxs0= sigahb(ip)/sigp!
O
             if(wxs.lt.wxs0) then 
O
 8505         jp= 1 
O
              volhp= volhtb(ip) '
O
 8005         test0= volhb(jp,ip)/volhpp
O
              pru0= g05caf(dmy) $
O
              if(pru0.lt.test0) then
O
               do l=1,9y$
O
                aprx(l)= g05caf(dmy)8
O
                prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))*.
O
     #                  aprx(l)+gblhb(jp,l,ip)
O
               enddo
O
              else"
O
               if(jp.lt.nlim) then)
O
                volhp= volhp-volhb(jp,ip)*
O
                jp= jp+1
O
                go to 8005      
O
               else 
O
                jp= jp+1
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy) 9
O
                 prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))*(/
O
     #                   aprx(l)+gblhb(jp,l,ip)e
O
                enddo 
O
               endif
O
              endif 
O
              pru= g05caf(dmy)#
O
              if(omssm.eq.'n') thenn)
O
                 xt0= wtoxsnh64(ndim,prx) #
O
                 xt= xt0/sigahb(ip)o(
O
              else if(omssm.eq.'y') then*
O
                 xt0= wtoxsnsh64(ndim,prx)#
O
                 xt= xt0/sigahb(ip) 
O
              endif.,
O
              prh= oxsmxhb(jp,ip)/sigahb(ip)
O
              prtst= xt-pru*prhy3
O
              if(prtst.ge.0.d0.and.xt.ne.0.d0) then 
O
                 do i=1,ninv(
O
                     prxa(i,it)= xaph(i)
O
                 enddo
O
                 amps(it)= xt0
O
                 it= it+1e
O
              else p
O
                 go to 8505 
O
              endif 
O
             endif$
O
             if(ip.le.(ipmx-1)) then#
O
              sigp= sigp-sigahb(ip) 
O
             elsep
O
 8506         jp= 1l
O
              volhp= volhtb(ip)h'
O
 8006         test0= volhb(jp,ip)/volhp 
O
              pru0= g05caf(dmy) $
O
              if(pru0.lt.test0) then
O
               do l=1,9t$
O
                aprx(l)= g05caf(dmy)8
O
                prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))*.
O
     #                  aprx(l)+gblhb(jp,l,ip)
O
               enddo
O
              else"
O
               if(jp.lt.nlim) then)
O
                volhp= volhp-volhb(jp,ip)
O

O
                jp= jp+1
O
                go to 8006     
O

O
               elseC
O
                jp= jp+1
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy)
O
9
O
                 prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))*x/
O
     #                   aprx(l)+gblhb(jp,l,ip)g
O
                enddo(
O
               endif
O
              endif+
O
              pru= g05caf(dmy)#
O
              if(omssm.eq.'n') then0'
O
               xt0= wtoxsnh64(ndim,prx)g!
O
               xt= xt0/sigahb(ip))(
O
              else if(omssm.eq.'y') then(
O
               xt0= wtoxsnsh64(ndim,prx)!
O
               xt= xt0/sigahb(ip)9
O
              endifc,
O
              prh= oxsmxhb(jp,ip)/sigahb(ip)
O
              prtst= xt-pru*prhp3
O
              if(prtst.ge.0.d0.and.xt.ne.0.d0) theng
O
                 do i=1,ninv(
O
                     prxa(i,it)= xaph(i)
O
                 enddo
O
                 amps(it)= xt0
O
                 it= it+1
O

O
              else 
O

O
                 go to 8506
O

O
              endif 
O
             endif
O
            enddo
O

O
           enddo
O
          enddo
O

O
          if(it.le.itmx) then1
O
             go to 8517
O

O
          endif 

O
* 9
O
         else if(otype.eq.'nc33'.or.otype.eq.'nc26') then 
O
          it= 1 
O
          sigt= 0.d0
O
          do i=1,12d!
O
             sigt= sigt+sigahq(i)i
O
          enddoo
O
          tsig= sigt
O
          sigt2= sigt-sigahq(1)  
O
          sigt3= sigt2-sigahq(2) 
O
          sigt4= sigt3-sigahq(3) 
O
          sigt5= sigt4-sigahq(4) 
O
          sigt6= sigt5-sigahq(5) 
O
          sigt7= sigt6-sigahq(6) 
O
          sigt8= sigt7-sigahq(7) 
O
          sigt9= sigt8-sigahq(8)!
O
          sigt10= sigt9-sigahq(9) #
O
          sigt11= sigt10-sigahq(10) 
O
 2517     wxs= g05caf(dmy)
O
          wxs0= sigahq(1)/sigt
O
          if(wxs.lt.wxs0) then
O
           ix0= 1l
O
           it0= 1t
O
           ip0= 1p
O
 2505      jp= 1
O
           volhp= volhtq(1)i#
O
 4005      test0= volhq(jp,1)/volhp 
O
           pru0= g05caf(dmy)!
O
           if(pru0.lt.test0) thenx
O
            do l=1,9!
O
             aprx(l)= g05caf(dmy) 3
O
             prx(l)= (gbuhq(jp,l,1)-gblhq(jp,l,1))* *
O
     #               aprx(l)+gblhq(jp,l,1)
O
            enddol
O
           else 
O
            if(jp.lt.nlim) then %
O
             volhp= volhp-volhq(jp,1) 
O
             jp= jp+1=
O
             go to 4005     
O
            else
O
             jp= jp+1
O

O
             do l=1,9 "
O
              aprx(l)= g05caf(dmy)4
O
              prx(l)= (gbuhq(jp,l,1)-gblhq(jp,l,1))*,
O
     #                 aprx(l)+gblhq(jp,l,1)
O
             enddo
O
            endif 
O
           endif
O
           pru= g05caf(dmy)x#
O
           if(otype.eq.'nc33') thenp$
O
            xt0= wtoxsnh32(ndim,prx)
O
            xt= xt0/sigahq(1)e(
O
           else if(otype.eq.'nc26') then$
O
            xt0= wtoxsnh26(ndim,prx)
O
            xt= xt0/sigahq(1)k
O
           endif'
O
           prh= oxsmxhq(jp,1)/sigahq(1) 
O
           prtst= xt-pru*prh0
O
           if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O
              do i=1,ninv)%
O
                  prxa(i,it)= xaph(i) 
O
              enddox
O
              amps(it)= xt02
O
              it= it+1
O
           else 
O
              go to 2505
O
           endif
O
          else
O
           wxs2= g05caf(dmy) 
O
           wxs0= sigahq(2)/sigt2 
O
           if(wxs2.lt.wxs0) then
O
            ix0= 1
O
            it0= 1
O
            ip0= 2
O
 2506       jp= 1j
O
            volhp= volhtq(2)$
O
 4006       test0= volhq(jp,2)/volhp
O
            pru0= g05caf(dmy) "
O
            if(pru0.lt.test0) then
O
             do l=1,9 "
O
              aprx(l)= g05caf(dmy)4
O
              prx(l)= (gbuhq(jp,l,2)-gblhq(jp,l,2))*+
O
     #                aprx(l)+gblhq(jp,l,2) 
O
             enddo
O
            else 
O
             if(jp.lt.nlim) then&
O
              volhp= volhp-volhq(jp,2)
O
              jp= jp+1
O
              go to 4006      
O
             else
O

O
              jp= jp+1
O
              do l=1,9#
O
               aprx(l)= g05caf(dmy) 5
O
               prx(l)= (gbuhq(jp,l,2)-gblhq(jp,l,2))*3-
O
     #                  aprx(l)+gblhq(jp,l,2) 
O
              enddor
O
             endif
O
            endifa
O
            pru= g05caf(dmy)$
O
            if(otype.eq.'nc33') then%
O
             xt0= wtoxsnh32(ndim,prx)
O

O
             xt= xt0/sigahq(2))
O
            else if(otype.eq.'nc26') then#%
O
             xt0= wtoxsnh26(ndim,prx) 
O
             xt= xt0/sigahq(2)
O
            endif
O
(
O
            prh= oxsmxhq(jp,2)/sigahq(2)
O
            prtst= xt-pru*prh
O
1
O
            if(prtst.ge.0.d0.and.xt.ne.0.d0) then 
O
               do i=1,ninv&
O
                   prxa(i,it)= xaph(i)
O
               enddo
O
               amps(it)= xt0
O
               it= it+1 
O
            else ,
O
               go to 2506 
O
            endif,
O
           else 
O
            wxs3= g05caf(dmy) !
O
            wxs0= sigahq(3)/sigt3 !
O
            if(wxs3.lt.wxs0) then
O

O
             ix0= 1p
O
             it0= 1o
O
             ip0= 3 
O
 2507        jp= 1
O
             volhp= volhtq(3) %
O
 4007        test0= volhq(jp,3)/volhp)
O
             pru0= g05caf(dmy)#
O
             if(pru0.lt.test0) thenr
O
              do l=1,9#
O
               aprx(l)= g05caf(dmy) 5
O
               prx(l)= (gbuhq(jp,l,3)-gblhq(jp,l,3))* ,
O
     #                 aprx(l)+gblhq(jp,l,3)
O
              enddo
O

O
             elset!
O
              if(jp.lt.nlim) thenx'
O
               volhp= volhp-volhq(jp,3)r
O
               jp= jp+1p
O
               go to 4007     
O
              else
O
               jp= jp+1 
O
               do l=1,9 $
O
                aprx(l)= g05caf(dmy)6
O
                prx(l)= (gbuhq(jp,l,3)-gblhq(jp,l,3))*.
O
     #                   aprx(l)+gblhq(jp,l,3)
O
               enddo
O
              endif 
O
             endif
O
             pru= g05caf(dmy) %
O
             if(otype.eq.'nc33') then-&
O
              xt0= wtoxsnh32(ndim,prx)
O
              xt= xt0/sigahq(3)x*
O
             else if(otype.eq.'nc26') then&
O
              xt0= wtoxsnh26(ndim,prx)
O
              xt= xt0/sigahq(3)6
O
             endif)
O
             prh= oxsmxhq(jp,3)/sigahq(3).
O
             prtst= xt-pru*prh2
O
             if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O
                do i=1,ninv '
O
                    prxa(i,it)= xaph(i)i
O
                enddoi
O
                amps(it)= xt0t
O
                it= it+1
O
             else 
O
                go to 2507
O
             endif
O
            else
O
             wxs4= g05caf(dmy)"
O
             wxs0= sigahq(4)/sigt4"
O
             if(wxs4.lt.wxs0) then
O
              ix0= 1
O
              it0= 2
O
              ip0= 1
O
 2508         jp= 1 
O
              volhp= volhtq(4)&
O
 4008         test0= volhq(jp,4)/volhp
O
              pru0= g05caf(dmy) $
O
              if(pru0.lt.test0) then
O
               do l=1,9r$
O
                aprx(l)= g05caf(dmy)6
O
                prx(l)= (gbuhq(jp,l,4)-gblhq(jp,l,4))*-
O
     #                  aprx(l)+gblhq(jp,l,4)=
O
               enddo
O
              else"
O
               if(jp.lt.nlim) then(
O
                volhp= volhp-volhq(jp,4)
O
                jp= jp+1
O
                go to 4008     j
O
               elsep
O
                jp= jp+1
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy)u7
O
                 prx(l)= (gbuhq(jp,l,4)-gblhq(jp,l,4))*t/
O
     #                    aprx(l)+gblhq(jp,l,4) 
O
                enddoh
O
               endif
O
              endif 
O
              pru= g05caf(dmy)&
O
              if(otype.eq.'nc33') then'
O
               xt0= wtoxsnh32(ndim,prx)t 
O
               xt= xt0/sigahq(4)+
O
              else if(otype.eq.'nc26') then '
O
               xt0= wtoxsnh26(ndim,prx)  
O
               xt= xt0/sigahq(4)
O
              endift*
O
              prh= oxsmxhq(jp,4)/sigahq(4)
O
              prtst= xt-pru*prh13
O
              if(prtst.ge.0.d0.and.xt.ne.0.d0) then 
O
                 do i=1,ninv(
O
                     prxa(i,it)= xaph(i)
O
                 enddo
O
                 amps(it)= xt0
O
                 it= it+1 
O
              else  
O
                 go to 2508g
O
              endif 
O
             elsel
O
              wxs5= g05caf(dmy)
O
#
O
              wxs0= sigahq(5)/sigt5 #
O
              if(wxs5.lt.wxs0) then 
O
               ix0= 1 
O
               it0= 2 
O
               ip0= 2 
O
 2509          jp= 1
O
               volhp= volhtq(5) '
O
 4009          test0= volhq(jp,5)/volhp  
O
               pru0= g05caf(dmy)%
O
               if(pru0.lt.test0) thenx
O
                do l=1,9%
O
                 aprx(l)= g05caf(dmy))7
O
                 prx(l)= (gbuhq(jp,l,5)-gblhq(jp,l,5))* .
O
     #                   aprx(l)+gblhq(jp,l,5)
O
                enddo 
O
               elseo#
O
                if(jp.lt.nlim) then1)
O
                 volhp= volhp-volhq(jp,5) 
O
                 jp= jp+1  
O
                 go to 4009     
O
                else
O
                 jp= jp+1 
O
                 do l=1,9l&
O
                  aprx(l)= g05caf(dmy)8
O
                  prx(l)= (gbuhq(jp,l,5)-gblhq(jp,l,5))*0
O
     #                     aprx(l)+gblhq(jp,l,5)
O
                 enddo
O
                endifn
O
               endif
O
               pru= g05caf(dmy) '
O
               if(otype.eq.'nc33') thenl(
O
                xt0= wtoxsnh32(ndim,prx)!
O
                xt= xt0/sigahq(5)p,
O
               else if(otype.eq.'nc26') then(
O
                xt0= wtoxsnh26(ndim,prx)!
O
                xt= xt0/sigahq(5)n
O
               endif+
O
               prh= oxsmxhq(jp,5)/sigahq(5)  
O
               prtst= xt-pru*prh4
O
               if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O
                  do i=1,ninvr)
O
                      prxa(i,it)= xaph(i)x
O
                  enddo 
O
                  amps(it)= xt0 
O
                  it= it+1
O
               else 
O
                  go to 2509
O
               endif
O
              else 
O
               wxs6= g05caf(dmy)$
O
               wxs0= sigahq(6)/sigt6$
O
               if(wxs6.lt.wxs0) then
O
                ix0= 1
O
                it0= 2
O
                ip0= 3
O
 2510           jp= 1  
O
                volhp= volhtq(6)(
O
 4010           test0= volhq(jp,6)/volhp!
O
                pru0= g05caf(dmy) &
O
                if(pru0.lt.test0) then
O
                 do l=1,9 &
O
                  aprx(l)= g05caf(dmy)8
O
                  prx(l)= (gbuhq(jp,l,6)-gblhq(jp,l,6))*/
O
     #                    aprx(l)+gblhq(jp,l,6) 
O
                 enddo
O
                else$
O
                 if(jp.lt.nlim) then*
O
                  volhp= volhp-volhq(jp,6)
O
                  jp= jp+1!
O
                  go to 4010     e
O
                 elsef
O
                  jp= jp+1
O
                  do l=1,9'
O
                   aprx(l)= g05caf(dmy) 9
O
                   prx(l)= (gbuhq(jp,l,6)-gblhq(jp,l,6))* 1
O
     #                      aprx(l)+gblhq(jp,l,6) 
O
                  enddom
O
                 endif
O
                endifj 
O
                pru= g05caf(dmy)(
O
                if(otype.eq.'nc33') then)
O
                 xt0= wtoxsnh32(ndim,prx) "
O
                 xt= xt0/sigahq(6)-
O
                else if(otype.eq.'nc26') then
O
)
O
                 xt0= wtoxsnh26(ndim,prx)t"
O
                 xt= xt0/sigahq(6)
O
                endifn,
O
                prh= oxsmxhq(jp,6)/sigahq(6)!
O
                prtst= xt-pru*prh 5
O
                if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O

O
                   do i=1,ninv*
O
                       prxa(i,it)= xaph(i)
O
                   enddo 
O
                   amps(it)= xt0
O
                   it= it+1 
O
                else  
O
                   go to 2510 
O
                endifp
O
               else !
O
                wxs7= g05caf(dmy)e%
O
                wxs0= sigahq(7)/sigt7 %
O
                if(wxs7.lt.wxs0) then 
O
                 ix0= 2 
O
                 it0= 1 
O
                 ip0= 11
O
 2511            jp= 1!
O
                 volhp= volhtq(7) )
O
 4011            test0= volhq(jp,7)/volhp
O
"
O
                 pru0= g05caf(dmy)'
O
                 if(pru0.lt.test0) then
O

O
                  do l=1,9'
O
                   aprx(l)= g05caf(dmy) 9
O
                   prx(l)= (gbuhq(jp,l,7)-gblhq(jp,l,7))*=0
O
     #                     aprx(l)+gblhq(jp,l,7)
O
                  enddo(
O
                 else(%
O
                  if(jp.lt.nlim) theno+
O
                   volhp= volhp-volhq(jp,7)t
O
                   jp= jp+1b"
O
                   go to 4011     
O
                  else
O
                   jp= jp+1 
O
                   do l=1,9
O
(
O
                    aprx(l)= g05caf(dmy):
O
                    prx(l)= (gbuhq(jp,l,7)-gblhq(jp,l,7))*2
O
     #                       aprx(l)+gblhq(jp,l,7)
O
                   enddo
O
                  endif 
O
                 endif!
O
                 pru= g05caf(dmy) )
O
                 if(otype.eq.'nc33') then**
O
                  xt0= wtoxsnh32(ndim,prx)#
O
                  xt= xt0/sigahq(7) .
O
                 else if(otype.eq.'nc26') then*
O
                  xt0= wtoxsnh26(ndim,prx)#
O
                  xt= xt0/sigahq(7) 
O
                 endif-
O
                 prh= oxsmxhq(jp,7)/sigahq(7)i"
O
                 prtst= xt-pru*prh6
O
                 if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O
                    do i=1,ninv +
O
                        prxa(i,it)= xaph(i)n
O
                    enddot!
O
                    amps(it)= xt0 
O
                    it= it+1
O
                 else 
O
                    go to 2511
O
                 endif
O
                else"
O
                 wxs8= g05caf(dmy)&
O
                 wxs0= sigahq(8)/sigt8&
O
                 if(wxs8.lt.wxs0) then
O
                  ix0= 2
O
                  it0= 1
O
                  ip0= 2
O
 2512             jp= 1 "
O
                  volhp= volhtq(8)*
O
 4012             test0= volhq(jp,8)/volhp#
O
                  pru0= g05caf(dmy) (
O
                  if(pru0.lt.test0) then
O
                   do l=1,9 (
O
                    aprx(l)= g05caf(dmy):
O
                    prx(l)= (gbuhq(jp,l,8)-gblhq(jp,l,8))*1
O
     #                      aprx(l)+gblhq(jp,l,8) 
O
                   enddo
O
                  else&
O
                   if(jp.lt.nlim) then,
O
                    volhp= volhp-volhq(jp,8)
O
                    jp= jp+1#
O
                    go to 4012     
O

O
                   else1
O
                    jp= jp+1
O
                    do l=1,9)
O
                     aprx(l)= g05caf(dmy)0;
O
                     prx(l)= (gbuhq(jp,l,8)-gblhq(jp,l,8))*03
O
     #                        aprx(l)+gblhq(jp,l,8)w
O
                    enddo 
O
                   endif
O
                  endif
O
"
O
                  pru= g05caf(dmy)*
O
                  if(otype.eq.'nc33') then+
O
                   xt0= wtoxsnh32(ndim,prx)f$
O
                   xt= xt0/sigahq(8)/
O
                  else if(otype.eq.'nc26') then +
O
                   xt0= wtoxsnh26(ndim,prx)p$
O
                   xt= xt0/sigahq(8)
O
                  endif .
O
                  prh= oxsmxhq(jp,8)/sigahq(8)#
O
                  prtst= xt-pru*prh 7
O
                  if(prtst.ge.0.d0.and.xt.ne.0.d0) then  
O
                     do i=1,ninv,
O
                         prxa(i,it)= xaph(i)
O
                     enddo"
O
                     amps(it)= xt0
O
                     it= it+1h
O
                  else )
O
                     go to 2512(
O
                  endif 
O
                 else #
O
                  wxs9= g05caf(dmy)
O
'
O
                  wxs0= sigahq(9)/sigt9 '
O
                  if(wxs9.lt.wxs0) then 
O
                   ix0= 2 
O
                   it0= 1s
O
                   ip0= 3l
O
 2513              jp= 1#
O
                   volhp= volhtq(9),+
O
 4013              test0= volhq(jp,9)/volhp $
O
                   pru0= g05caf(dmy))
O
                   if(pru0.lt.test0) thenr
O
                    do l=1,9)
O
                     aprx(l)= g05caf(dmy) ;
O
                     prx(l)= (gbuhq(jp,l,9)-gblhq(jp,l,9))*x2
O
     #                       aprx(l)+gblhq(jp,l,9)
O
                    enddo 
O
                   else '
O
                    if(jp.lt.nlim) then -
O
                     volhp= volhp-volhq(jp,9) 
O
                     jp= jp+1 $
O
                     go to 4013     
O
                    else
O
                     jp= jp+1=
O
                     do l=1,9t*
O
                      aprx(l)= g05caf(dmy)<
O
                      prx(l)= (gbuhq(jp,l,9)-gblhq(jp,l,9))*4
O
     #                         aprx(l)+gblhq(jp,l,9)
O
                     enddo
O
                    endif*
O
                   endif#
O
                   pru= g05caf(dmy) +
O
                   if(otype.eq.'nc33') then(,
O
                    xt0= wtoxsnh32(ndim,prx)%
O
                    xt= xt0/sigahq(9)10
O
                   else if(otype.eq.'nc26') then,
O
                    xt0= wtoxsnh26(ndim,prx)%
O
                    xt= xt0/sigahq(9) 
O
                   endif/
O
                   prh= oxsmxhq(jp,9)/sigahq(9) $
O
                   prtst= xt-pru*prh8
O
                   if(prtst.ge.0.d0.and.xt.ne.0.d0) then!
O
                      do i=1,ninvr-
O
                          prxa(i,it)= xaph(i)t
O
                      enddox#
O
                      amps(it)= xt0x
O
                      it= it+1
O
                
 else  
O
                      go to 2513
O
                   endif
O
                  else%
O
                   wxs10= g05caf(dmy)x*
O
                   wxs0= sigahq(10)/sigt10)
O
                   if(wxs10.lt.wxs0) thena
O
                    ix0= 2
O
                    it0= 2
O
                    ip0= 1
O
 2514               jp= 1e%
O
                    volhp= volhtq(10) -
O
 4014               test0= volhq(jp,10)/volhp
O
%
O
                    pru0= g05caf(dmy) *
O
                    if(pru0.lt.test0) then
O
                     do l=1,9 *
O
                      aprx(l)= g05caf(dmy)>
O
                      prx(l)= (gbuhq(jp,l,10)-gblhq(jp,l,10))*4
O
     #                        aprx(l)+gblhq(jp,l,10)
O
                     enddo
O
                    else(
O
                     if(jp.lt.nlim) then/
O
                      volhp= volhp-volhq(jp,10)s
O
                      jp= jp+1%
O
                      go to 4014     3
O
                     else=
O
                      jp= jp+1
O
                      do l=1,9+
O
                       aprx(l)= g05caf(dmy)s?
O
                       prx(l)= (gbuhq(jp,l,10)-gblhq(jp,l,10))*t5
O
     #                         aprx(l)+gblhq(jp,l,10)1
O
                      enddo 
O
                     endif
O
                    endifn$
O
                    pru= g05caf(dmy),
O
                    if(otype.eq.'nc33') then-
O
                     xt0= wtoxsnh32(ndim,prx)h'
O
                     xt= xt0/sigahq(10)m1
O
                    else if(otype.eq.'nc26') then -
O
                     xt0= wtoxsnh26(ndim,prx) '
O
                     xt= xt0/sigahq(10)p
O
                    endif 2
O
                    prh= oxsmxhq(jp,10)/sigahq(10)%
O
                    prtst= xt-pru*prht9
O
                    if(prtst.ge.0.d0.and.xt.ne.0.d0) thenj"
O
                       do i=1,ninv.
O
                           prxa(i,it)= xaph(i)
O
                       enddo$
O
                       amps(it)= xt0
O
                       it= it+1j
O
                    else  !
O
                       go to 2514 
O
                    endifi
O
                   else &
O
                    wxs11= g05caf(dmy)+
O
                    wxs0= sigahq(11)/sigt11s*
O
                    if(wxs11.lt.wxs0) then
O
                     ix0= 2e
O
                     it0= 2t
O
                     ip0= 2 
O
 2515                jp= 1&
O
                     volhp= volhtq(11).
O
 4015                test0= volhq(jp,11)/volhp&
O
                     pru0= g05caf(dmy)+
O
                     if(pru0.lt.test0) then 
O
                      do l=1,9+
O
                       aprx(l)= g05caf(dmy) ?
O
                       prx(l)= (gbuhq(jp,l,11)-gblhq(jp,l,11))*t5
O
     #                         aprx(l)+gblhq(jp,l,11)s
O
                      enddos
O
                     else )
O
                      if(jp.lt.nlim) then10
O
                       volhp= volhp-volhq(jp,11)
O
                       jp= jp+1l&
O
                       go to 4015     
O
                      else
O
                       jp= jp+1t
O
                       do l=1,9 ,
O
                        aprx(l)= g05caf(dmy)@
O
                        prx(l)= (gbuhq(jp,l,11)-gblhq(jp,l,11))*6
O
     #                          aprx(l)+gblhq(jp,l,11)
O
                       enddo
O
                      endifh
O
                     endif%
O
                     pru= g05caf(dmy)0-
O
                     if(otype.eq.'nc33') then .
O
                      xt0= wtoxsnh32(ndim,prx)(
O
                      xt= xt0/sigahq(11)2
O
                     else if(otype.eq.'nc26') then.
O
                      xt0= wtoxsnh26(ndim,prx)(
O
                      xt= xt0/sigahq(11)
O
                     endif3
O
                     prh= oxsmxhq(jp,11)/sigahq(11) &
O
                     prtst= xt-pru*prh:
O
                     if(prtst.ge.0.d0.and.xt.ne.0.d0) then#
O
                        do i=1,ninvp/
O
                            prxa(i,it)= xaph(i) 
O
                        enddoq%
O
                        amps(it)= xt0- 
O
                        it= it+1
O
                     else "
O
                        go to 2515
O
                     endif
O
                    else
O
                     ix0= 2
O

O
                     it0= 2 
O
                     ip0= 3 
O
 2516                jp= 1&
O
                     volhp= volhtq(12).
O
 4016                test0= volhq(jp,12)/volhp&
O
                     pru0= g05caf(dmy)+
O
                     if(pru0.lt.test0) then 
O
                      do l=1,9+
O
                       aprx(l)= g05caf(dmy) ?
O
                       prx(l)= (gbuhq(jp,l,12)-gblhq(jp,l,12))*
O
5
O
     #                         aprx(l)+gblhq(jp,l,12)l
O
                      enddo 
O
                     else )
O
                      if(jp.lt.nlim) then 0
O
                       volhp= volhp-volhq(jp,12)
O
                       jp= jp+1t&
O
                       go to 4016     
O
                      else
O
                       jp= jp+1+
O
                       do l=1,9 ,
O
                        aprx(l)= g05caf(dmy)@
O
                        prx(l)= (gbuhq(jp,l,12)-gblhq(jp,l,12))*7
O
     #                           aprx(l)+gblhq(jp,l,12) 
O
                       enddo
O
                      endif 
O
                     endif%
O
                     pru= g05caf(dmy) -
O
                     if(otype.eq.'nc33') thent.
O
                      xt0= wtoxsnh32(ndim,prx)(
O
                      xt= xt0/sigahq(12)2
O
                     else if(otype.eq.'nc26') then.
O
                      xt0= wtoxsnh26(ndim,prx)(
O
                      xt= xt0/sigahq(12)
O
                     endif3
O
                     prh= oxsmxhq(jp,12)/sigahq(12)s&
O
                     prtst= xt-pru*prh:
O
                     if(prtst.ge.0.d0.and.xt.ne.0.d0) then!
O
                      do i=1,ninv *
O
                       prxa(i,it)= xaph(i)
O
                      enddo #
O
                      amps(it)= xt0e
O
                      it= it+1
O
                     else  
O
                      go to 2516
O
                     endif
O
                    endif
O

O
                   endif
O
                  endif 
O
                 endif
O
                endif4
O
               endif
O
              endif 
O
             endif
O
            endif 
O
           endif
O
          endif 
O
          if(it.le.itmx) then 
O
             go to 2517 
O
          endif)

O
*g
         elsel

O
*j1
O
         if((otype.eq.'mx43'.and.ofs.eq.'ll').or.lB
O
     #      (otype.eq.'mx43'.and.ofs.eq.'qq'.and.opglu.eq.'n').or.4
O
     #      (otype.eq.'nc32'.and.opglu.eq.'n')) then

O
*
O
A
O
         if(otype.eq.'mx43'.and.ofs.eq.'qq'.and.opglu.eq.'n'.and. 
O
     #      oint.eq.'y') thenp
O
            iselsp= 1 
         else 
O
            iselsp= 0c
O
         endif

O
* 
O
         it= 1
O
         sigt= 0.d0p
O
         do i=1,4  
O
             sigt= sigt+sigah(i)
O
         enddo
O
         tsig= sigt 
O
         sigt2= sigt-sigah(1)i
O
         sigt3= sigt2-sigah(2)
O
 9513    wxs= g05caf(dmy)e
O
         wxs0= sigah(1)/sigt
O
         if(wxs.lt.wxs0) then 
O
          ix0= 1
O
          it0= 1
O
          ip0= 1
O
 9505     jp= 1e
O
          volhp= volht(1)x!
O
 9005     test0= volh(jp,1)/volhp=
O
          pru0= g05caf(dmy)  
O
          if(pru0.lt.test0) then
O
           do l=1,9  
O
            aprx(l)= g05caf(dmy)0
O
            prx(l)= (gbuh(jp,l,1)-gblh(jp,l,1))*(
O
     #              aprx(l)+gblh(jp,l,1)
O
           enddo
O
          else
O
           if(jp.lt.nlim) then#
O
            volhp= volhp-volh(jp,1)i
O
            jp= jp+1
O
            go to 9005      
O
           else 
O
            jp= jp+1
O
            do l=1,9!
O
             aprx(l)= g05caf(dmy) 1
O
             prx(l)= (gbuh(jp,l,1)-gblh(jp,l,1))*w*
O
     #                aprx(l)+gblh(jp,l,1)
O
            enddo 
O
           endif
O
          endifj
O
          pru= g05caf(dmy)"
O
          if(otype.eq.'nc32') then$
O
             xt0= wtoxsn32(ndim,prx)
O
             xt= xt0/sigah(1)('
O
          else if(otype.eq.'mx43') thenl
O
             iint= 0$
O
             xt0= wtoxsm43(ndim,prx)
O
             xt= xt0/sigah(1)(
O
          endif#$
O
          prh= oxsmxh(jp,1)/sigah(1)
O
          prtst= xt-pru*prh /
O
          if(prtst.ge.0.d0.and.xt.ne.0.d0) then)
O
             do i=1,ninv$
O
                 prxa(i,it)= xaph(i)
O
             enddo
O
             amps(it)= xt0!
O
             if(iselsp.eq.1) then 
O
    
          iint= 1 '
O
                xt1= wtoxsm43(ndim,prx)l
O
                iint= 2 '
O
                xt2= wtoxsm43(ndim,prx),
O
                amps1(it)= xt1
O
                amps2(it)= xt2
O
             endif
O
             it= it+1 
O
          else 
O

O
             go to 9505f
O
          endif 
         elseq
O
          wxs2= g05caf(dmy) 
O
          wxs0= sigah(2)/sigt2
O
          if(wxs2.lt.wxs0) then 
O
           ix0= 1f
O
           it0= 2h
O
           ip0= 1 
O
 9506      jp= 1
O
           volhp= volht(2)"
O
 9006      test0= volh(jp,2)/volhp
O
           pru0= g05caf(dmy)!
O
           if(pru0.lt.test0) thent
O
            do l=1,9!
O
             aprx(l)= g05caf(dmy).1
O
             prx(l)= (gbuh(jp,l,2)-gblh(jp,l,2))* )
O
     #               aprx(l)+gblh(jp,l,2) 
O
            enddo 
O
           elsem
O
            if(jp.lt.nlim) thent$
O
             volhp= volhp-volh(jp,2)
O
             jp= jp+1 
O
             go to 9006     
O
            else
O
             jp= jp+1)
O
             do l=1,9s"
O
              aprx(l)= g05caf(dmy)2
O
              prx(l)= (gbuh(jp,l,2)-gblh(jp,l,2))**
O
     #                aprx(l)+gblh(jp,l,2)
O
             enddo
O
            endif=
O
           endif
O
           pru= g05caf(dmy)v#
O
           if(otype.eq.'nc32') theny%
O
              xt0= wtoxsn32(ndim,prx)h
O
              xt= xt0/sigah(2)(
O
           else if(otype.eq.'mx43') then%
O
              xt0= wtoxsm43(ndim,prx)g
O
              xt= xt0/sigah(2)
O
           endif%
O
           prh= oxsmxh(jp,2)/sigah(2)
O

O
           prtst= xt-pru*prh0
O
           if(prtst.ge.0.d0.and.xt.ne.0.d0) then
O
              do i=1,ninv %
O
                  prxa(i,it)= xaph(i) 
O
              enddo 
O
              amps(it)= xt0 "
O
              if(iselsp.eq.1) then
O
                 iint= 1(
O
                 xt1= wtoxsm43(ndim,prx)
O
                 iint= 2(
O
                 xt2= wtoxsm43(ndim,prx)
O
                 amps1(it)= xt1 
O
                 amps2(it)= xt2 
O
              endif 
O
              it= it+1
O
           else 
O
              go to 9506
O
           endif
O
          else
O
           wxs3= g05caf(dmy)
O
           wxs0= sigah(3)/sigt3( 
O
           if(wxs3.lt.wxs0) then
O
            ix0= 2
O
            it0= 1
O
            ip0= 1
O
 9507       jp= 1x
O
            volhp= volht(3) #
O
 9007       test0= volh(jp,3)/volhpx
O
            pru0= g05caf(dmy) "
O
            if(pru0.lt.test0) then
O
             do l=1,9d"
O
              aprx(l)= g05caf(dmy)2
O
              prx(l)= (gbuh(jp,l,3)-gblh(jp,l,3))**
O
     #                aprx(l)+gblh(jp,l,3)
O
             enddo
O
            else 
O
             if(jp.lt.nlim) then%
O
              volhp= volhp-volh(jp,3) 
O
              jp= jp+1
O
              go to 9007      
O
             elsey
O
              jp= jp+1
O
              do l=1,9#
O
               aprx(l)= g05caf(dmy) 3
O
               prx(l)= (gbuh(jp,l,3)-gblh(jp,l,3))* ,
O
     #                  aprx(l)+gblh(jp,l,3)
O
              enddop
O
             endif
O
            endifo
O
            pru= g05caf(dmy)$
O
            if(otype.eq.'nc32') then&
O
               xt0= wtoxsn32(ndim,prx)
O
               xt= xt0/sigah(3) )
O
            else if(otype.eq.'mx43') thenp
O
               iint= 0&
O
               xt0= wtoxsm43(ndim,prx)
O
               xt= xt0/sigah(3) 
O
            endif &
O
            prh= oxsmxh(jp,3)/sigah(3)
O
            prtst= xt-pru*prh 1
O
            if(prtst.ge.0.d0.and.xt.ne.0.d0) then 
O
               do i=1,ninv&
O
                   prxa(i,it)= xaph(i)
O
               enddo
O
               amps(it)= xt0#
O
               if(iselsp.eq.1) thenx
O
                  iint= 1 )
O
                  xt1= wtoxsm43(ndim,prx),
O
                  iint= 2 )
O
                  xt2= wtoxsm43(ndim,prx)  
O
                  amps1(it)= xt1 
O
                  amps2(it)= xt2
O
               endif
O
               it= it+1i
O
            else t
O
               go to 9507=
O
            endif)
O
           else 
O
            ix0= 2
O
            it0= 2
O
            ip0= 1
O
 9512       jp= 1 
O
            volhp= volht(4)r#
O
 9012       test0= volh(jp,4)/volhp7
O
            pru0= g05caf(dmy) "
O
            if(pru0.lt.test0) then
O
               do l=1,9p'
O
    
O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O














































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































                prx(l)= (gbuh(jp,l,2)-gblh(jp,l,2))*
     #                   aprx(l)+gblh(jp,l,2)
               enddo
              endif
             endif
             pru= g05caf(dmy)
             if(otype.eq.'nc50') then
                xt0= wtoxsnh49(ndim,prx)
             else if(otype.eq.'nc25') then
                xt0= wtoxsnh24(ndim,prx)
             else if(otype.eq.'nc21') then
                xt0= wtoxsnh19(ndim,prx)
             else if(otype.eq.'nc32') then
                xt0= wtoxsn32(ndim,prx)
             else if(otype.eq.'cc12') then
                xt0= wtoxssc(ndim,prx)
             else if(otype.eq.'mx43') then
                iint= 0
                xt0= wtoxsm43(ndim,prx)
             endif
             xt= xt0/sigah(2)
             prh= oxsmxh(jp,2)/sigah(2)
             prtst= xt-pru*prh
             if(prtst.ge.0.d0.and.xt.ne.0.d0) then
                do i=1,ninv
                    prxa(i,it)= xaph(i)
                enddo
                amps(it)= xt0
                it= it+1
             else 
                go to 1508
             endif
            else
             wxs5= g05caf(dmy)
             wxs0= sigah(7)/sigt5
             if(wxs5.lt.wxs0) then
              ix0= 2
              it0= 2
              ip0= 1
 1509         jp= 1
              volhp= volht(7)
 3009         test0= volh(jp,7)/volhp
              pru0= g05caf(dmy)
              if(pru0.lt.test0) then
               do l=1,9
                aprx(l)= g05caf(dmy)
                prx(l)= (gbuh(jp,l,7)-gblh(jp,l,7))*
     #                   aprx(l)+gblh(jp,l,7)
               enddo
              else
               if(jp.lt.nlim) then
                volhp= volhp-volh(jp,7)
                jp= jp+1
                go to 3009     
               else
                jp= jp+1
                do l=1,9
                 aprx(l)= g05caf(dmy)
                 prx(l)= (gbuh(jp,l,7)-gblh(jp,l,7))*
     #                    aprx(l)+gblh(jp,l,7)
                enddo
               endif
              endif
              pru= g05caf(dmy)
              if(otype.eq.'nc50') then
                 xt0= wtoxsnh49(ndim,prx)
              else if(otype.eq.'nc25') then
                 xt0= wtoxsnh24(ndim,prx)
              else if(otype.eq.'nc21') then
                 xt0= wtoxsnh19(ndim,prx)
              else if(otype.eq.'nc32') then
                 xt0= wtoxsn32(ndim,prx)
              else if(otype.eq.'cc12') then
                 xt0= wtoxssc(ndim,prx)
              else if(otype.eq.'mx43') then
                 iint= 0
                 xt0= wtoxsm43(ndim,prx)
              endif
              xt= xt0/sigah(7)
              prh= oxsmxh(jp,7)/sigah(7)
              prtst= xt-pru*prh
              if(prtst.ge.0.d0.and.xt.ne.0.d0) then
                 do i=1,ninv
                     prxa(i,it)= xaph(i)
                 enddo
                 amps(it)= xt0
                 it= it+1
              else 
                 go to 1509
              endif
             else
              wxs6= g05caf(dmy)
              wxs0= sigah(6)/sigt6
              if(wxs6.lt.wxs0) then
               ix0= 2
               it0= 1
               ip0= 2
 1510          jp= 1
               volhp= volht(6)
 3010          test0= volh(jp,6)/volhp
               pru0= g05caf(dmy)
               if(pru0.lt.test0) then
                do l=1,9
                 aprx(l)= g05caf(dmy)
                 prx(l)= (gbuh(jp,l,6)-gblh(jp,l,6))*
     #                    aprx(l)+gblh(jp,l,6)
                enddo
               else
                if(jp.lt.nlim) then
                 volhp= volhp-volh(jp,6)
                 jp= jp+1
                 go to 3010     
                else
                 jp= jp+1
                 do l=1,9
                  aprx(l)= g05caf(dmy)
                  prx(l)= (gbuh(jp,l,6)-gblh(jp,l,6))*
     #                     aprx(l)+gblh(jp,l,6)
                 enddo
                endif
               endif
               pru= g05caf(dmy)
               if(otype.eq.'nc50') then
                  xt0= wtoxsnh49(ndim,prx)
               else if(otype.eq.'nc25') then
                  xt0= wtoxsnh24(ndim,prx)
               else if(otype.eq.'nc21') then
                  xt0= wtoxsnh19(ndim,prx)
               else if(otype.eq.'nc32') then
                  xt0= wtoxsn32(ndim,prx)
               else if(otype.eq.'cc12') then
                  xt0= wtoxssc(ndim,prx)
               else if(otype.eq.'mx43') then
                  iint= 0
                  xt0= wtoxsm43(ndim,prx)
               endif
               xt= xt0/sigah(6)
               prh= oxsmxh(jp,6)/sigah(6)
               prtst= xt-pru*prh
               if(prtst.ge.0.d0.and.xt.ne.0.d0) then
                 do i=1,ninv
                     prxa(i,it)= xaph(i)
                 enddo
                 amps(it)= xt0
                 it= it+1
               else 
                  go to 1510
               endif
              else
               wxs7= g05caf(dmy)
               wxs0= sigah(4)/sigt7
               if(wxs7.lt.wxs0) then
                ix0= 1
                it0= 2
                ip0= 2
 1511           jp= 1
                volhp= volht(4)
 3011           test0= volh(jp,4)/volhp
                pru0= g05caf(dmy)
                if(pru0.lt.test0) then
                 do l=1,9
                  aprx(l)= g05caf(dmy)
                  prx(l)= (gbuh(jp,l,4)-gblh(jp,l,4))*
     #                     aprx(l)+gblh(jp,l,4)
                 enddo
                else
                 if(jp.lt.nlim) then
                  volhp= volhp-volh(jp,4)
                  jp= jp+1
                  go to 3011     
                 else
                  jp= jp+1
                  do l=1,9
                   aprx(l)= g05caf(dmy)
                   prx(l)= (gbuh(jp,l,4)-gblh(jp,l,4))*
     #                      aprx(l)+gblh(jp,l,4)
                  enddo
                 endif
                endif
                pru= g05caf(dmy)
                if(otype.eq.'nc50') then
                   xt0= wtoxsnh49(ndim,prx)
                else if(otype.eq.'nc25') then
                   xt0= wtoxsnh24(ndim,prx)
                else if(otype.eq.'nc21') then
                   xt0= wtoxsnh19(ndim,prx)
                else if(otype.eq.'nc32') then
                   xt0= wtoxsn32(ndim,prx)
                else if(otype.eq.'cc12') then
                   xt0= wtoxssc(ndim,prx)
                else if(otype.eq.'mx43') then
                   iint= 0
                   xt0= wtoxsm43(ndim,prx)
                endif
                xt= xt0/sigah(4)
                prh= oxsmxh(jp,4)/sigah(4)
                prtst= xt-pru*prh
                if(prtst.ge.0.d0.and.xt.ne.0.d0) then
                   do i=1,ninv
                       prxa(i,it)= xaph(i)
                   enddo
                   amps(it)= xt0
                   it= it+1
                else 
                   go to 1511
                endif
               else
                ix0= 2
                it0= 2
                ip0= 2
 1512           jp= 1
                volhp= volht(8)
 3012           test0= volh(jp,8)/volhp
                pru0= g05caf(dmy)
                if(pru0.lt.test0) then
                 do l=1,9
                  aprx(l)= g05caf(dmy)
                  prx(l)= (gbuh(jp,l,8)-gblh(jp,l,8))*
     #                     aprx(l)+gblh(jp,l,8)
                 enddo
                else
                 if(jp.lt.nlim) then
                  volhp= volhp-volh(jp,8)
                  jp= jp+1
                  go to 3012     
                 else
                  jp= jp+1
                  do l=1,9
                   aprx(l)= g05caf(dmy)
                   prx(l)= (gbuh(jp,l,8)-gblh(jp,l,8))*
     #                      aprx(l)+gblh(jp,l,8)
                  enddo
                 endif
                endif
                pru= g05caf(dmy)
                if(otype.eq.'nc50') then
                   xt0= wtoxsnh49(ndim,prx)
                else if(otype.eq.'nc25') then
                   xt0= wtoxsnh24(ndim,prx)
                else if(otype.eq.'nc21') then
                   xt0= wtoxsnh19(ndim,prx)
                else if(otype.eq.'nc32') then
                   xt0= wtoxsn32(ndim,prx)
                else if(otype.eq.'cc12') then
                   xt0= wtoxssc(ndim,prx)
                else if(otype.eq.'mx43') then
                   iint= 0
                   xt0= wtoxsm43(ndim,prx)
                endif
                xt= xt0/sigah(8)
                prh= oxsmxh(jp,8)/sigah(8)
                prtst= xt-pru*prh
                if(prtst.ge.0.d0.and.xt.ne.0.d0) then
                   do i=1,ninv
                       prxa(i,it)= xaph(i)
                   enddo
                   amps(it)= xt0
                   it= it+1
                else 
                   go to 1512
                endif
               endif
              endif
             endif
            endif
           endif
          endif
         endif
         if(it.le.itmx) then
            go to 1513
         endif
*
*-----CC or Higgs/NC generation ends
*
           endif
          endif
         endif
* 
*-----CC or Higgs storage starts
*
         if(ostore.eq.'d') then
          if(oxcm.eq.'c') then
             xint= 10.d0
             xintl= 25.d0
             xintwp= 2.d0
             xinte= 25.d0
             xint4j= 40.d0
             do i=1,im
                 exi= 75.d0+(i-1.d0)*xint/im
                 exf= 75.d0+i*xint/im
                 exli= -15.d0+(i-1.d0)*xintl/im
                 exlf= -15.d0+i*xintl/im
                 exwpi= -1.d0+(i-1.d0)*xintwp/im
                 exwpf= -1.d0+i*xintwp/im
                 exei= (i-1.d0)*xinte/im
                 exef= i*xinte/im
                 ex4ji= 90.d0+(i-1.d0)*xint4j/im
                 ex4jf= 75.d0+i*xint4j/im
                 exc(i)= 0.d0
                 excl(i)= 0.d0
                 excwp(i)= 0.d0
                 exce(i)= 0.d0
                 exc4j(i)= 0.d0
                 do j=1,itmx
                     cvv= prxa(1,j)*prxa(2,j)
                     extst(j)= cvv*prxa(4,j)
                     if(extst(j).gt.0.d0) then
                      extst(j)= sqrt(extst(j))*rs
                      if(exf.gt.extst(j).and.extst(j).gt.exi) then
                         exc(i)= exc(i)+1.d0
                      endif
                      he1= prxa(3,j)+prxa(6,j)+prxa(7,j)
                      ht1= prxa(9,j)
                      if((prxa(2,j)*(he1-ht1)).gt.0.d0) then
                       qtrf(j)= log(prxa(2,j)*(he1-ht1))
                       if(exlf.gt.qtrf(j).and.qtrf(j).gt.exli) then
                        excl(i)= excl(i)+1.d0
                       endif
                      endif
                     endif
                     hewp= 0.5d0*(prxa(2,j)*(1.d0+prxa(4,j)-
     #                     prxa(3,j))-(prxa(2,j)-prxa(1,j))*
     #                     (1.d0-prxa(8,j)))
                     hbwps= 1.d0-cvv*prxa(4,j)/hewp/hewp
                     if(hbwps.gt.0.d0) then
                        hbwp= sqrt(hbwps)
                        cwpst(j)= (1.d0-prxa(1,j)*(1.d0-prxa(8,j))/
     #                             hewp)/hbwp
                        if(exwpf.gt.cwpst(j).and.cwpst(j).gt.
     #                     exwpi) then
                           excwp(i)= excwp(i)+1.d0
                        endif
                     endif
                     exetst(j)= prxa(2,j)*(prxa(3,j)+prxa(6,j)+
     #                          prxa(7,j))-(prxa(2,j)-prxa(1,j))*
     #                          prxa(9,j)*rs
                     if(exef.gt.exetst(j).and.exetst(j).gt.exei) then
                        exce(i)= exce(i)+1.d0
                     endif
                     xiv1= sqrt(cvv*prxa(3,j))*rs
                     xiv2= sqrt(cvv*prxa(4,j))*rs
                     xiv3= sqrt(cvv*prxa(5,j))*rs
                     xiv4= sqrt(cvv*prxa(6,j))*rs
                     xiv5= sqrt(cvv*prxa(7,j))*rs
                     xiv6= sqrt(cvv)*sqrt(1.d0-prxa(3,j)-
     #                     prxa(4,j)-prxa(5,j)-prxa(6,j)-
     #                     prxa(7,j))*rs
                     div1= abs(xiv1-xiv2)
                     div2= abs(xiv3-xiv4)
                     div3= abs(xiv5-xiv6)
                     divall= dmin1(div1,div2,div3)
                     if(divall.eq.div1) then
                        sum4j= xiv1+xiv2
                     else if(divall.eq.div2) then
                        sum4j= xiv3+xiv4
                     else if(divall.eq.div3) then
                        sum4j= xiv5+xiv6
                     endif
                     if(ex4ji.lt.sum4j.and.sum4j.lt.ex4jf) then
                        exc4j(i)= exc4j(i)+1
                     endif
                 enddo
             enddo
*
             ern= xint/im
             open(nout,file='mco',status='new')
             do i=1,im
                 exb= 75.d0+(i-1.d0)*xint/im
                 ecx= 75.d0+(i-0.5d0)*xint/im
                 write(nout,fmt=2500) exb,ecx,exc(i)/itmx/ern*tsig
             enddo
             close(nout)      
             ern4j= xint4j/im
             open(nout,file='mco4j',status='new')
             do i=1,im
                 exb4j= 90.d0+(i-1.d0)*xint4j/im
                 ecx4j= 90.d0+(i-0.5d0)*xint4j/im
                 write(nout,fmt=2500) exb4j,ecx4j,exc4j(i)/itmx/
     #                                ern4j*tsig
             enddo
             close(nout)      
             ernl= xintl/im
             open(nout,file='mcol',status='new')
             do i=1,im
                 exb= -15.d0+(i-1.d0)*xintl/im
                 ecx= -15.d0+(i-0.5d0)*xintl/im
                 write(nout,fmt=2500) exb,ecx,excl(i)/itmx/ernl*tsig
             enddo
             close(nout)      
             ernwp= xintwp/im
             open(nout,file='mcowp',status='new')
             do i=1,im
                 exb= -1.d0+(i-1.d0)*xintwp/im
                 ecx= -1.d0+(i-0.5d0)*xintwp/im
                 write(nout,fmt=2500) exb,ecx,
     #                                excwp(i)/itmx/ernwp*tsig
             enddo
             close(nout)      
             erne= xinte/im
             open(nout,file='mcoe',status='new')
             do i=1,im
                 exeb= (i-1.d0)*xinte/im
                 ecex= (i-0.5d0)*xinte/im
                 write(nout,fmt=2500) exeb,ecex,exce(i)/itmx/erne*tsig
             enddo
             close(nout)      
          else if(oxcm.eq.'n') then
             xint= 60.d0
             xintv= 100.d0
             xintqt= 120.d0
             do i=1,im
                 exi= 60.d0+(i-1.d0)*xint/im
                 exf= 60.d0+i*xint/im
                 exiv= 30.d0+(i-1.d0)*xintv/im
                 exfv= 30.d0+i*xintv/im
                 if(otype.eq.'nc68') then
                  excm(i)= 0.d0
                  do j=1,itmx
                   bthr= 50.d0
                   hvv= prxa(1,j)*prxa(2,j)
                   extst1(j)= sqrt(hvv*prxa(3,j))*rs
                   extst2(j)= sqrt(hvv*prxa(4,j))*rs
                   extst3(j)= sqrt(hvv*prxa(5,j))*rs
                   extst4(j)= sqrt(hvv*prxa(6,j))*rs
                   extst5(j)= sqrt(hvv*prxa(7,j))*rs
                   if(oev.eq.'a1') then
                    secr= 1.d0-prxa(3,j)-prxa(4,j)-prxa(5,j)-
     #                    prxa(6,j)-prxa(7,j)
                    if(secr.gt.0.d0) then
                     extst6(j)= sqrt(hvv*secr)*rs
                     if(exf.gt.extst1(j).and.extst1(j).gt.exi.
     #                  and.extst2(j).gt.bthr) then
                        excm(i)= excm(i)+1.d0
                     else if(exf.gt.extst2(j).and.
     #                       extst2(j).gt.exi.and.
     #                       extst1(j).gt.bthr) then
                        excm(i)= excm(i)+1.d0
                     else if(exf.gt.extst3(j).and.
     #                       extst3(j).gt.exi.and.
     #                       extst4(j).gt.bthr) then
                        excm(i)= excm(i)+1.d0
                     else if(exf.gt.extst4(j).and.
     #                       extst4(j).gt.exi.and.
     #                       extst3(j).gt.bthr) then
                        excm(i)= excm(i)+1.d0
                     else if(exf.gt.extst5(j).and.
     #                       extst5(j).gt.exi.and.
     #                       extst6(j).gt.bthr) then
                        excm(i)= excm(i)+1.d0
                     else if(exf.gt.extst6(j).and.
     #                       extst6(j).gt.exi.and.
     #                       extst5(j).gt.bthr) then
                        excm(i)= excm(i)+1.d0
                     endif
                    endif
                   else if(oev.eq.'a2') then
                    rim1= extst1(j)
                    rim2= extst2(j)
                    rim3= extst3(j)
                    rim4= extst4(j)
                    rim5= extst5(j)
                    secr= 1.d0-prxa(3,j)-prxa(4,j)-prxa(5,j)-
     #                    prxa(6,j)-prxa(7,j)
                    if(secr.gt.0.d0) then
                     extst6(j)= sqrt(hvv*secr)*rs
                     rim6= extst6(j)
                     iev= ievs(exi,exf,rim1,rim2,rim3,rim4,rim5,rim6)
                     if(iev.eq.1) then
                        excm(i)= excm(i)+1.d0
                     endif
                    endif
                   endif
                  enddo
                 endif
                 exc1(i)= 0.d0
                 exc2(i)= 0.d0
                 do j=1,itmx
                  hvv= prxa(1,j)*prxa(2,j)
                  extst1(j)= sqrt(hvv*prxa(4,j))*rs
                  extst2(j)= sqrt(hvv*prxa(3,j))*rs
                  if(exf.gt.extst1(j).and.extst1(j).gt.exi) then
                     exc1(i)= exc1(i)+1.d0
                  endif
                  if(exf.gt.extst2(j).and.extst2(j).gt.exi) then
                     exc2(i)= exc2(i)+1.d0
                  endif
                 enddo
             enddo
             ern= xint/im
             ern= tsig/itmx/ern
             if(otype.eq.'nc68') then
              open(nout,file='mcombb',status='new')
              do i=1,im
                  exb= 60.d0+(i-1.d0)*xint/im
                  ecx= 60.d0+(i-0.5d0)*xint/im
                  write(nout,fmt=2500) exb,ecx,excm(i)*ern
              enddo
              close(nout)      
             else
              open(nout,file='mcombb',status='new')
              do i=1,im
                  exb= 60.d0+(i-1.d0)*xint/im
                  ecx= 60.d0+(i-0.5d0)*xint/im
                  write(nout,fmt=2500) exb,ecx,exc1(i)*ern,
     #                                exc2(i)*ern
              enddo
              close(nout)      
             endif
*
             if(ipr.eq.31.or.ipr.eq.34.or.ipr.eq.35.or.
     #          ipr.eq.10.or.ipr.eq.13.or.ipr.eq.14) then
                icl= 1
             else if(ipr.eq.30.or.ipr.eq.32.or.ipr.eq.33.or.
     #               ipr.eq.7.or.ipr.eq.19.or.ipr.eq.25) then
                icl= 2
             endif
             do i=1,im
                 exi= -1.d0+(i-1.d0)*2.d0/im
                 exf= -1.d0+i*2.d0/im
                 exim= 60.d0+(i-1.d0)*xint/im
                 exfm= 60.d0+i*xint/im
                 exiv= 50.d0+(i-1.d0)*xintv/im
                 exfv= 50.d0+i*xintv/im
                 exiqt= (i-1.d0)*xintqt/im
                 exfqt= i*xintqt/im
                 exc1(i)= 0.d0
                 exc2(i)= 0.d0
                 exc3(i)= 0.d0
                 exc4(i)= 0.d0
                 exc5(i)= 0.d0
                 exc6(i)= 0.d0
                 exc7(i)= 0.d0
                 do j=1,itmx
                  hvv= prxa(1,j)*prxa(2,j)
                  he1= prxa(3,j)+prxa(6,j)+prxa(7,j)
                  he2= 1.d0-prxa(4,j)-prxa(6,j)-prxa(7,j)
                  he3= prxa(4,j)+prxa(5,j)+prxa(7,j)
                  he4= 1.d0-prxa(3,j)-prxa(5,j)-prxa(7,j)
                  ht1= prxa(9,j)
                  ht2= prxa(8,j)-prxa(9,j)
                  ht3= prxa(10,j)
                  ht4= 1.d0-prxa(8,j)-prxa(10,j)
                  hxp= prxa(2,j)
                  hxm= prxa(1,j) 
                  hct1= 1.d0-2.d0*hxm*ht1/(hxp*he1-(hxp-hxm)*ht1)
                  hct2= 1.d0-2.d0*hxm*ht2/(hxp*he2-(hxp-hxm)*ht2)
                  hct3= 1.d0-2.d0*hxm*ht3/(hxp*he3-(hxp-hxm)*ht3)
                  hct4= 1.d0-2.d0*hxm*ht4/(hxp*he4-(hxp-hxm)*ht4)
                  hct12= 1.d0-0.5d0*hvv*prxa(3,j)/he1/he2
                  hct34= 1.d0-0.5d0*hvv*prxa(4,j)/he3/he4
                  he12= 1.d0+prxa(3,j)-prxa(4,j)
                  he34= 1.d0+prxa(4,j)-prxa(3,j)
                  hem= 0.5d0*(hxp*he12-(hxp-hxm)*prxa(8,j))
                  hbem= 1.d0-hvv*prxa(3,j)/hem/hem
                  if(hbem.le.0.d0) then
                     isrm= 0
                  else
                     isrm= 1
                     hbem= sqrt(hbem)
                     hctsumm= (1.d0-hxm*prxa(8,j)/hem)/hbem
                  endif
                  hep= 0.5d0*(hxp*he34-(hxp-hxm)*(1.d0-prxa(8,j)))
                  hbep= 1.d0-hvv*prxa(4,j)/hep/hep
                  if(hbep.le.0.d0) then
                     isrp= 0
                  else
                     isrp= 1
                     hbep= sqrt(hbep)
                     hctsump= (1.d0-hxm*(1.d0-prxa(8,j))/hep)/hbep
                  endif
                  rmmm= s*(1.d0-2.d0*hem+hvv*prxa(3,j))
                  rmmp= s*(1.d0-2.d0*hep+hvv*prxa(4,j))
                  rmmm= sqrt(rmmm)
                  rmmp= sqrt(rmmp)
                  if(icl.eq.2) then
                   tq2= -hvv*prxa(3,j)+hem*hem
                   beqs= 1.d0-hvv*prxa(3,j)/hem/hem
                   if(beqs.lt.0.d0.or.tq2.lt.0.d0) then
                      iqt= 0
                   endif
                   beq= sqrt(beqs)
                   tq= sqrt(tq2)
                   ctq= (1.d0-hxm*prxa(8,j)/hem)/beq
                   ctqf= 1.d0-ctq*ctq
                   if(ctqf.lt.0.d0) then
                      iqt= 0
                   endif
                   tqt= tq*sqrt(ctqf)*rs
                   iqt= 1
                  endif
                  if(icl.eq.1) then
                   if((exf.gt.hct3.and.hct3.gt.exi).or.
     #                (exf.gt.hct4.and.hct4.gt.exi)) then
                      exc1(i)= exc1(i)+1.d0
                   endif
                   if(exf.gt.hct34.and.hct34.gt.exi) then
                      exc2(i)= exc2(i)+1.d0
                   endif
                   if((isrp.eq.1).and.(exf.gt.hctsump.and.
     #                 hctsump.gt.exi)) then
                       exc3(i)= exc3(i)+1.d0
                   endif
                   if(exfm.gt.rmmm.and.rmmm.gt.exim) then
                      exc4(i)= exc4(i)+1.d0
                   endif
                   if((isrm.eq.1).and.(exfm.gt.rmmm.and.
     #                 rmmm.gt.exim).and.(0.8d0.gt.hctsumm.and.
     #                 hctsumm.gt.-0.8d0)) then
                       exc5(i)= exc5(i)+1.d0
                   endif
                   if(exfv.gt.(hep*rs).and.(hep*rs).gt.exiv) then
                      exc6(i)= exc6(i)+1.d0
                   endif
                  else if(icl.eq.2) then
                   if((exf.gt.hct1.and.hct1.gt.exi).or.
     #                (exf.gt.hct2.and.hct2.gt.exi)) then
                      exc1(i)= exc1(i)+1.d0
                   endif
                   if(exf.gt.hct12.and.hct12.gt.exi) then
                      exc2(i)= exc2(i)+1.d0
                   endif
                   if((isrm.eq.1).and.(exf.gt.hctsumm.and.
     #                 hctsumm.gt.exi)) then
                       exc3(i)= exc3(i)+1.d0
                   endif
                   if(exfm.gt.rmmp.and.rmmp.gt.exim) then
                      exc4(i)= exc4(i)+1.d0
                   endif
                   if((isrp.eq.1).and.(exfm.gt.rmmp.and.
     #                 rmmp.gt.exim).and.(0.8d0.gt.hctsump.and.
     #                 hctsump.gt.-0.8d0)) then
                       exc5(i)= exc5(i)+1.d0
                   endif
                   if(exfv.gt.(hem*rs).and.(hem*rs).gt.exiv) then
                      exc6(i)= exc6(i)+1.d0
                   endif
                   if((iqt.eq.1).and.(exfv.gt.tqt.and.
     #                 tqt.gt.exiv)) then
                       exc7(i)= exc7(i)+1.d0
                   endif
                  endif
                 enddo
             enddo
*
             if(otype.ne.'nc68') then
              ern= 2.d0/im
              ern= tsig/itmx/ern
              open(nout,file='mcocbb',status='new')
              do i=1,im
                  exb= -1.d0+(i-1.d0)*2.d0/im
                  ecx= -1.d0+(i-0.5d0)*2.d0/im
                  write(nout,fmt=2501) exb,ecx,exc1(i)*ern,
     #                                 exc2(i)*ern,exc3(i)*ern
              enddo
              close(nout)      
              open(nout,file='mcomm',status='new')
              ern= xint/im
              ern= tsig/itmx/ern
              do i=1,im
                  exb= 60.d0+(i-1.d0)*xint/im
                  ecx= 60.d0+(i-0.5d0)*xint/im
                  write(nout,fmt=2500) exb,ecx,exc4(i)*ern,
     #                                 exc5(i)*ern
              enddo
              close(nout)      
              ern= xintv/im
              ern= tsig/itmx/ern
              open(nout,file='mcoev',status='new')
              do i=1,im
                  exb= 30.d0+(i-1.d0)*xintv/im
                  ecx= 30.d0+(i-0.5d0)*xintv/im
                  write(nout,fmt=2500) exb,ecx,exc6(i)*ern
              enddo
              close(nout)      
              if(icl.eq.2) then
               ern= xintqt/im
               ern= tsig/itmx/ern
               open(nout,file='mcoqt',status='new')
               do i=1,im
                   exb= (i-1.d0)*xintqt/im
                   ecx= (i-0.5d0)*xintqt/im
                   write(nout,fmt=2500) exb,ecx,exc7(i)*ern
               enddo
               close(nout)      
              endif
             endif
*
          endif
          print 2503,itmx
         else if(ostore.eq.'m') then
          open(nout,file='mom',status='new')
          if(iosf.eq.0) then 
             nhep= 4
          else
             nhep= 6
          endif
*
              if(iosf.gt.0) then
                 idhep(5)= 22
                 idhep(6)= 22
              endif
              if(ipr.eq.1) then
                 idhep(1)= 16
                 idhep(2)= -15
                 idhep(3)= 13
                 idhep(4)= -14
              else if(ipr.eq.2) then
                 if(ockm.eq.'n') then
                  idhep(1)= 2
                  idhep(2)= -1
                  idhep(3)= 13
                  idhep(4)= -14
                 else
                  if(ickm.eq.1) then
                   idhep(1)= 2
                   idhep(2)= -1
                   idhep(3)= 13
                   idhep(4)= -14
                  else if(ickm.eq.2) then
                   idhep(1)= 2
                   idhep(2)= -3
                   idhep(3)= 13
                   idhep(4)= -14
                  else if(ickm.eq.3) then
                   idhep(1)= 2
                   idhep(2)= -5
                   idhep(3)= 13
                   idhep(4)= -14
                  else if(ickm.eq.4) then
                   idhep(1)= 4
                   idhep(2)= -1
                   idhep(3)= 13
                   idhep(4)= -14
                  else if(ickm.eq.5) then
                   idhep(1)= 4
                   idhep(2)= -3
                   idhep(3)= 13
                   idhep(4)= -14
                  else if(ickm.eq.6) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 13
                   idhep(4)= -14
                  endif
                 endif
              else if(ipr.eq.3) then
                 if(ockm.eq.'n') then
                  idhep(1)= 4
                  idhep(2)= -3
                  idhep(3)= 1
                  idhep(4)= -2
                 else
                  if(ickm.eq.1) then
                   idhep(1)= 4
                   idhep(2)= -3
                   idhep(3)= 1
                   idhep(4)= -2
                  else if(ickm.eq.2) then
                   idhep(1)= 4
                   idhep(2)= -1
                   idhep(3)= 1
                   idhep(4)= -2
                  else if(ickm.eq.3) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 1
                   idhep(4)= -2
                  else if(ickm.eq.4) then
                   idhep(1)= 4
                   idhep(2)= -3
                   idhep(3)= 3
                   idhep(4)= -2
                  else if(ickm.eq.5) then
                   idhep(1)= 4
                   idhep(2)= -1
                   idhep(3)= 3
                   idhep(4)= -2
                  else if(ickm.eq.6) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 3
                   idhep(4)= -2
                  else if(ickm.eq.7) then
                   idhep(1)= 4
                   idhep(2)= -3
                   idhep(3)= 5
                   idhep(4)= -2
                  else if(ickm.eq.8) then
                   idhep(1)= 4
                   idhep(2)= -1
                   idhep(3)= 5
                   idhep(4)= -2
                  else if(ickm.eq.9) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 5
                   idhep(4)= -2
                  else if(ickm.eq.10) then
                   idhep(1)= 2
                   idhep(2)= -3
                   idhep(3)= 1
                   idhep(4)= -2
                  else if(ickm.eq.11) then
                   idhep(1)= 2
                   idhep(2)= -5
                   idhep(3)= 1
                   idhep(4)= -2
                  else if(ickm.eq.12) then
                   idhep(1)= 2
                   idhep(2)= -5
                   idhep(3)= 3
                   idhep(4)= -2
                  else if(ickm.eq.13) then
                   idhep(1)= 4
                   idhep(2)= -3
                   idhep(3)= 1
                   idhep(4)= -4
                  else if(ickm.eq.14) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 1
                   idhep(4)= -4
                  else if(ickm.eq.15) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 3
                   idhep(4)= -4
                  endif
                 endif
              else if(ipr.eq.4) then
                 idhep(1)= 14
                 idhep(2)= -13
                 idhep(3)= 11
                 idhep(4)= -12
              else if(ipr.eq.5) then
                 if(ockm.eq.'n') then
                  idhep(1)= 2
                  idhep(2)= -1
                  idhep(3)= 11
                  idhep(4)= -12
                 else
                  if(ickm.eq.1) then
                   idhep(1)= 2
                   idhep(2)= -1
                   idhep(3)= 11
                   idhep(4)= -12
                  else if(ickm.eq.2) then
                   idhep(1)= 2
                   idhep(2)= -3
                   idhep(3)= 11
                   idhep(4)= -12
                  else if(ickm.eq.3) then
                   idhep(1)= 2
                   idhep(2)= -5
                   idhep(3)= 11
                   idhep(4)= -12
                  else if(ickm.eq.4) then
                   idhep(1)= 4
                   idhep(2)= -1
                   idhep(3)= 11
                   idhep(4)= -12
                  else if(ickm.eq.5) then
                   idhep(1)= 4
                   idhep(2)= -3
                   idhep(3)= 11
                   idhep(4)= -12
                  else if(ickm.eq.6) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 11
                   idhep(4)= -12
                  endif
                 endif
              else if(ipr.eq.6) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 16
                 idhep(4)= -16
              else if(ipr.eq.7) then
                 idhep(1)= 1
                 idhep(2)= -1
                 idhep(3)= 14
                 idhep(4)= -14
              else if(ipr.eq.8) then
                 idhep(1)= 2
                 idhep(2)= -2
                 idhep(3)= 14
                 idhep(4)= -14
              else if(ipr.eq.9) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 15
                 idhep(4)= -15
              else if(ipr.eq.10) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 1
                 idhep(4)= -1
              else if(ipr.eq.11) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 2
                 idhep(4)= -2
              else if(ipr.eq.12) then
                 idhep(1)= 14
                 idhep(2)= -14
                 idhep(3)= 16 
                 idhep(4)= -16
              else if(ipr.eq.13) then
                 idhep(1)= 3
                 idhep(2)= -3
                 idhep(3)= 2
                 idhep(4)= -2
              else if(ipr.eq.14) then
                 idhep(1)= 1
                 idhep(2)= -1
                 idhep(3)= 3
                 idhep(4)= -3
              else if(ipr.eq.15) then
                 idhep(1)= 2
                 idhep(2)= -2
                 idhep(3)= 4
                 idhep(4)= -4
              else if(ipr.eq.16) then
                 idhep(1)= 14
                 idhep(2)= -14
                 idhep(3)= 12
                 idhep(4)= -12
              else if(ipr.eq.17) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 12
                 idhep(4)= -12
              else if(ipr.eq.18) then
                 idhep(1)= 2
                 idhep(2)= -2
                 idhep(3)= 12
                 idhep(4)= -12
              else if(ipr.eq.19) then
                 idhep(1)= 1
                 idhep(2)= -1
                 idhep(3)= 12
                 idhep(4)= -12
              else if(ipr.eq.20) then
                 idhep(1)= 14
                 idhep(2)= -13
                 idhep(3)= 13
                 idhep(4)= -14
              else if(ipr.eq.21) then
                 if(ockm.eq.'n') then
                  idhep(1)= 2
                  idhep(2)= -1
                  idhep(3)= 1
                  idhep(4)= -2
                 else
                  if(ickm.eq.1) then
                   idhep(1)= 2
                   idhep(2)= -1
                   idhep(3)= 1
                   idhep(4)= -2
                  else if(ickm.eq.2) then
                   idhep(1)= 2
                   idhep(2)= -3
                   idhep(3)= 3
                   idhep(4)= -2
                  else if(ickm.eq.3) then
                   idhep(1)= 2
                   idhep(2)= -5
                   idhep(3)= 5
                   idhep(4)= -2
                  else if(ickm.eq.4) then
                   idhep(1)= 4
                   idhep(2)= -1
                   idhep(3)= 1
                   idhep(4)= -4
                  else if(ickm.eq.5) then
                   idhep(1)= 4
                   idhep(2)= -3
                   idhep(3)= 3
                   idhep(4)= -4
                  else if(ickm.eq.6) then
                   idhep(1)= 4
                   idhep(2)= -5
                   idhep(3)= 5
                   idhep(4)= -4
                  endif
                 endif
              else if(ipr.eq.22) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 11
                 idhep(4)= -11
              else if(ipr.eq.23) then
                 idhep(1)= 14
                 idhep(2)= -14
                 idhep(3)= 11
                 idhep(4)= -11
              else if(ipr.eq.24) then
                 idhep(1)= 2
                 idhep(2)= -2
                 idhep(3)= 11
                 idhep(4)= -11
              else if(ipr.eq.25) then
                 idhep(1)= 1
                 idhep(2)= -1
                 idhep(3)= 11
                 idhep(4)= -11
              else if(ipr.eq.26) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 13
                 idhep(4)= -13
              else if(ipr.eq.27) then
                 idhep(1)= 14
                 idhep(2)= -14
                 idhep(3)= 14
                 idhep(4)= -14
              else if(ipr.eq.28) then
                 idhep(1)= 2
                 idhep(2)= -2
                 idhep(3)= 2
                 idhep(4)= -2
              else if(ipr.eq.29) then
                 idhep(1)= 1
                 idhep(2)= -1
                 idhep(3)= 1
                 idhep(4)= -1
              else if(ipr.eq.30) then
                 idhep(1)= 5
                 idhep(2)= -5
                 idhep(3)= 14
                 idhep(4)= -14
              else if(ipr.eq.31) then
                 idhep(1)= 13
                 idhep(2)= -13
                 idhep(3)= 5
                 idhep(4)= -5
              else if(ipr.eq.32) then
                 idhep(1)= 5
                 idhep(2)= -5
                 idhep(3)= 12
                 idhep(4)= -12
              else if(ipr.eq.33) then
                 idhep(1)= 5
                 idhep(2)= -5
                 idhep(3)= 11
                 idhep(4)= -11
              else if(ipr.eq.34) then
                 idhep(1)= 2
                 idhep(2)= -2
                 idhep(3)= 5
                 idhep(4)= -5
              else if(ipr.eq.35) then
                 idhep(1)= 1
                 idhep(2)= -1
                 idhep(3)= 5
                 idhep(4)= -5
              else if(ipr.eq.36) then
                 idhep(1)= 5
                 idhep(2)= -5
                 idhep(3)= 5
                 idhep(4)= -5
              else if(ipr.eq.37) then
                 idhep(1)= 15
                 idhep(2)= -15
                 idhep(3)= 5
                 idhep(4)= -5
              else if(ipr.eq.38) then
                 idhep(1)= 15
                 idhep(2)= -16
                 idhep(3)= 4
                 idhep(4)= -3
              endif
*
          do i=1,itmx
              gxp= prxa(2,i)
              gxm= prxa(1,i)
              ge1= prxa(3,i)+prxa(6,i)+prxa(7,i)
              ge2= 1.d0+prxa(3,i)-prxa(4,i)-ge1
              ge3= prxa(4,i)+prxa(5,i)+prxa(7,i)
              ge4= 1.d0+prxa(4,i)-prxa(3,i)-ge3
              gt1= prxa(9,i)
              gt2= prxa(8,i)-prxa(9,i)
              gt3= prxa(10,i)
              gt4= 1.d0-prxa(8,i)-gt3
*
              phep(4,1)= gxp*ge1-(gxp-gxm)*gt1
              phep(4,2)= gxp*ge2-(gxp-gxm)*gt2
              phep(4,3)= gxp*ge3-(gxp-gxm)*gt3
              phep(4,4)= gxp*ge4-(gxp-gxm)*gt4
              phep(3,1)= gxp*ge1-(gxp+gxm)*gt1
              phep(3,2)= gxp*ge2-(gxp+gxm)*gt2
              phep(3,3)= gxp*ge3-(gxp+gxm)*gt3
              phep(3,4)= gxp*ge4-(gxp+gxm)*gt4
              phep(2,1)= 0.d0
              phep(1,1)= phep(4,1)*phep(4,1)-phep(3,1)*phep(3,1)
              phep(1,1)= sqrt(phep(1,1))
              phep(1,2)= (-phep(3,1)*phep(3,2)+phep(4,1)*phep(4,2)-
     #                 2.d0*gxp*gxm*prxa(3,i))/phep(1,1)
          
O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O