*comdeck zc
double-precision declaration and named constants.
*if def,CHEAP
floating-point entities on 32-bit computers.
implicit double precision (a-h,o-z)
code name and version number.
character kod*8,ver*5
parameter (kod='mcnpx',ver='2.1.5')
processor-dependent named constants.
mdas is the initial size of dynamically allocated storage.
on systems where memory adjustment is not available, set mdas
large enough for your biggest problem.
*if def,POINTERS
parameter (mdas=20808)
*if -def,POINTERS
parameter (mdas=4000000)
ndp2 = numeric storage units per floating-point value.
huge = a very large number.
*if -def,CHEAP
parameter (ndp2=1,huge=1e123)
*if def,CHEAP
parameter (ndp2=2,huge=1d37)
array dimensions.
i/o unit numbers.
general constants.
include MIPT and MLPT parameter definitions, etc.
parameter(maxi=34,maxv=20,maxw=3,mcoh=55,mcpu=200,
+ minc=21,mink=200,m_lca=9,m_lcb=8,m_lea=8,m_leb=4,
*if def,MESHTAL
+ mipt=35,maxg=20,maxgv=10,
*if -def,MESHTAL
+ mipt=34,
+ mlpt=22,mjsf=9,mkft=9,mktc=22,mlgc=1000,mpb=5,mpng=21,mseb=301,
+ mspare=3,mtop=49,mwng=25,mxdt=20,mxdx=5,mxlv=10,nbmx=100,
+ ncolor=10,
*if def,MESHTAL
+ ndef=16,
*if -def,MESHTAL
+ ndef=15,
+ novr=5,nptr=13,nsp=602,nsp12=nsp+12,npkey=6,
+ maxsec=5,iemaxn=31,iemaxp=26,msigg=6,mxso=2,
+ ntp=201,nspt=nsp+ntp+7)
parameter(iui=31,iuo=32,iur=33,iux=34,iud=35,iup=37,ius=38,
+ iu1=39,iu2=40,iusw=41,iusr=42,iusc=43,iuc=44,iut=45,iuz=46,
+ iuk=47,iu3=48,iu4=49,iupw=50,iupc=51,iupx=52,iuw=53,
*if def,MESHTAL
For HISTP, the unit number nhstp=55 is defined elsewhere.
parameter(
7 aneut=1.,avogad=6.244d23,
8 avgdn=1.d-24*avogad/aneut,dfdmp=-15.,dftint=100.,
8 emax_def=1.0d+02,emax_top=1.0d+03,hundred=1.0d+02,
9 euler=.,p_1m12=1.0d-12,
9 half=0.5d0,two=2.d0,three=3.d0,ten=10.d0,p_1m6=1.0d-06,
9 fscon=137.0393d0,hsll=1.d-30,one=1.d0,pie=3.8d0,
9 planck=4.,slite=299.7925d0,third=one/three,zero=0.)
c ----------------------------------------------------------------------
*comdeck vv
tables and character commons.
common /tables/ gpt(mipt),charge(mipt),decay(mipt),
1 ecut_min(mipt),ecut_def(mipt),
1 qfiss(23),rkt(mtop),talb(8,2),vco(mcoh),vic(minc),wco(mcoh),
2 jrwb(16,5),jsf(mjsf),mfiss(22),nvs(maxv),itty,jtty
character common -- character variables and arrays.
character aid*80,aid1*80,aids*80,chcd*10,exms*80,hbln(maxv,2)*3,
1 hblw(maxw)*3,hcs(2)*7,hdpath*80,hdpth*80,hft(mkft)*3,hfu(2)*11,
2 hip*(mipt),hmes*69,hnp(mipt)*11,hovr*8,hsd(2)*10,hsub*31,ibin*8,
3 idtm*19,idtms*19,ilbl(8)*8,iname*8,klin*80,kods*8,ksf(29)*3,
4 loddat*28,lods*8,msub(ndef)*8,probid*19,probs*19,rfq(11)*57,
5 ufil(3,6)*11,vers*5
common /charcm/ aid,aid1,aids,chcd,exms,hbln,hblw,hcs,hdpath,
1 hdpth,hft,hfu,hip,hmes,hnp,hovr,hsd,hsub,ibin,idtm,idtms,ilbl,
2 iname,klin,kods,ksf,loddat,lods,msub,probid,probs,rfq,ufil,vers
names of files
character*8 inp,outp,runtpe,mctal,wssa,ptrac,comout,srctp,plotm,
1 rssa,xsdir,com,wwinp,
*if def,MESHTAL
+ dumn1,dumn2,isub(ndef)
common /charcm/ inp,outp,runtpe,mctal,wssa,ptrac,comout,srctp,
1 plotm,rssa,xsdir,com,wwinp,
*if def,MESHTAL
+ dumn1,dumn2
equivalence (isub,inp)
c ----------------------------------------------------------------------
*comdeck cm
common blocks for all sections.
*if def,MULTT.and.VMS,1
cpar$ private /pblcom/,/tskcom/
whole-block declarations.
V2.1.4: nfixcm decreased by 1 & increased by mipt for efac(mipt).
V2.1.4: lfixcm increased by mipt for in_erg(mipt).
parameter (nfixcm=m_lcb+m_leb+maxi+3*maxv+mtop+mipt*(25+7*mxdx)+
+ 2*mxdt+nsp+105)
parameter (lfixcm=m_lca+m_lea+
*if def,HISTP
*if def,MESHTAL
+ 4*maxg+maxg*maxgv+6+
+ 3*mxdt+mink+13*mipt+4*maxv+301)
dimension gfixcm(nfixcm),jfixcm(lfixcm)
equivalence (bbrem,gfixcm),(ibad,jfixcm)
parameter (nvarcm=107*mipt+140)
parameter (lvarcm=
*if def,RADIOG
+ mipt*(1+8*mxdx)+mcpu+329)
dimension gvarcm(nvarcm),jvarcm(lvarcm)
equivalence (cpk,gvarcm),(idum,jvarcm)
parameter (nephcm=29,lephcm=nptr+novr+ncolor+
*if def,MESHTAL
dimension gephcm(nephcm),jephcm(lephcm)
equivalence (cp0,gephcm),(ichan,jephcm)
logical lockl
********************
statically allocated common
*********************
fixed common -- constant after the problem is initiated.
common /fixcom/ bbrem(mtop),bcw(2,3),bnum,calph(maxi),coincd,
1 cut_n,ddg(2,mxdt),ddx(mipt,2,mxdx),dxw(mipt,3),dxx(mipt,5,mxdx),
1 ecf(mipt),efac(mipt),
1 emcf(mipt),emx(mipt),enum,espl(mipt,10),fnw,rim,hsb(nsp),
+ optlcb(m_lcb),optleb(m_leb),
+ rssp,rnok,rnfb,rnfs,rngb,rngs,rnmult,snit,srv(3,maxv),
2 tco(mipt),thgf(0:50),wc1(mipt),wc2(mipt),wwg(7),wwp(mipt,5),
+ wwm(26),zfixcm
common /fixcom/
3 ibad,icw,idefv(maxv),ides,idrc(mxdt),ifft,igm,igww,ikz,img,imt,
4 indt,ink(mink),iphot,iplt,ipty(mipt),isb,issw,istrg,
*if def,HISTP
4 ivdd(maxv),ivdis(maxv),ivord(maxv),jgm(mipt),jtlx,junf,kf8,kfl,
5 knods,knrm,kpt(mipt),ktls,kufil(2,6),lfcdg,lfcdj,
+ lcaopt(m_lca),leaopt(m_lea),
5 locdt(2,mxdt),lvcdg,lvcdj,lxs,mbnk,mcal,mct,mgegbt(mipt),
6 mgm(mipt+1),mgww(mipt+1),mix,mjss,mlaj,mlja,mrkp,mrl,msd,msrk,
7 mtasks,nlat,nsrc,in_erg(mipt),
6 mww(mipt+1),mxa,mxafs,mxe,mxe1,mxf,mxj,mxt,mxtr,mxxs,ndet(mipt),
7 ndnd,ndtt,ndx(mipt),nee_max,nee(mipt),isanti,
7 ngww(mipt),nhb,nilr,nilw,nips,niss,njsr,
8 njss,nkxs,nlev,nlja,nmat,nmxf,nnpos,nocoh,nord,np1,npikmt,npn,
8 ipert,mnnm,mxfp,nmaz,nmip,nmzu,npert,nrcd,nrss,nsph,nsr,nsrck,
+ nstrid,ntal,nvec,nwang,nwgeom,nwgm,nww(mipt),nwwm,nxnx
*if def,MESHTAL
common /fixcom/ nmesh,mtlflg,mspflg,medflg,mpdflg,mdxflg
offsets for virtual arrays in dynamically allocated storage.
common /fixcom/ lara,lasp,lava,lavz,
1 lawc,lawn,lcmg,lden,ldpt,ldrs,ldxp,
1 leaa,lear,leba,lebd,lebl,lebt,lech,ledg,leee,leek,legg,lelp,
2 lesa,lewg,lfim,lflc,lfme,lfmg,
2 laa9,larg,lexn,lexp,lesp,lfmf,lhsg,lpcp,lsgg,lsgx,lzz9,lxso,
2 lfor,lfrc,lfst,lftt,lgmg,lgvl,
3 lgwt,lpbr,lpbt,lpkn,lpmg,lpru,lpxr,lqav,lqax,lqcn,lrho,lrpt,
4 lscf,lsmg,lspf,lsqq,lsso,ltbt,ltds,ltmp,ltrf,ltth,lvcl,lvec,
5 lvol,lwgm,lwwe,lwwf,lwwk,lxnm,lx85,lipa,lipb,lpptb,lipt,liax,
5 lmel,lnel,liss,litd,lixl,liza,ljar,ljmd,ljpt,ljsc,ljss,lixs,
7 ljtf,ljun,ljvc,ljxs,lkcp,lksc,lksd,lkst,lksu,lktp,lkxs,llaf,
8 llat,llca,llfc,llft,llja,llme,llmt,llct,llph,llst,llsc,lmat,
9 lmbd,lmbi,lmfl,lmzp,lmzu,lncl,lngm,lnht,lnmt,lnpt,lnpq,lnsb,
1 lnsf,lnty,lnxs,lddm,lddn,ldec,ldxc,ldxd,lfeb,lflx,lfso,lgww,
2 lpac,lpan,lpcc,lpwb,lrkp,lshs,lstt,ltfc,lwns,
3 lise,ljfq,llaj,llcj,llse,lmaz,lndp,lndr,lnhs,lnpw,lnsl,lntb,
4 lscr,ldrc,lemi,lfdd,lgnr,lpik,lptb,lptr,lpts,lrng,lrtc,ltgp,
5 lifl,lpc2,lixc,ljfl,ljft,ljmt,lkmt,lktc,lkxd,llbb,
*if def,MESHTAL
+ ltgd(maxg),lgd1(maxg),lgd2(maxg),lgd3(maxg),lemsh(maxg*maxgv),
+ llgtsk,lljtsk,ltal,lgbn,lbnk,lxss,lexs,mfixcm
variable common -- variable but required for a continue run.
arrays that are backed up when a track is lost.
common /varcom/ cpk,cts,dbcn(30),dmp,eacc(mipt,3),racc(mipt,2),
1 osum(3),coll(mipt),
1 osum2(3,3),pax(6,16,mipt),prn,rani,ranj,rdum(50),rijk,rkk,
2 rlt(4,2),rnr,rsum(3),rsum2(3,3),smul(3),sumk(3),tmav(mipt,3),
3 twac,twss,wcs1(mipt),wcs2(mipt),wgts(2),wt0,wssi(7),
5 idum(50),inif,ist,ist0,ixak,ixak0,jrad,kc8,kcsf,kct,kcy,knod,
6 ksdef,kcz,lost(2),nbal(mcpu),nbhwm,nbov,nbt(mipt),
7 ndmp,nerr,netb(2),nfer,npc(20),npd,npnm,npp,nppm,nps,npsout,npsr,
8 nqss,nqsw,nrnh(3),nrrs,nrsw,nsa,nsa0,nskk,nss,nss0,nssi(8),ntc,
9 ntc1,ntss,nwer,nwsb,nwse,nwsg(3),nwst,nwws(2,99),
+ nziy(8,mxdx,mipt),
*if def,RADIOG
+ notal,notrn,ntprt(100),monod,npsmg,
ephemeral common -- not needed after the current run.
common /ephcom/ cp0,cp1,cp2,cp3,cpa,ctme,fpi,freq,
1 ssb(10),tdc,tlc,trm,wnvp(4),xhom,xnum,yhom,
3 ichan,ics,idmp,ifile,iln,iln1,inform,iovr,inpd,iptr,iptra(nptr),
4 irup,itask,iterm,itfxs,itotnu,iuou,jchar,jfcn,jgf,jgxa(2),
5 jgxo(2),jovr(novr),jtasks,jtfc,jvp,kbp,kcolor(ncolor+6),
6 kdbnps,kmplot,komout,konrun,kprod,krflg,krtm,ksr,ktfile,
7 lfatl,lfll,locki,lspeed,mazf(3),mcolor,mdc,mmkdb,mnk,
8 mpc,mrm,nde,nkrp,nomore,nrc,nst,ntasks,
*if def,MESHTAL
+ nugd(maxg),
9 mephcm,lockl,mynum
pbl common -- particle description
*if def,unicos.and.multt,1
task common /pblcom/ xxx,yyy,zzz,uuu,vvv,www,erg,wgt,tme,vel,dls,
*if -def,unicos.or.multt,1
common /pblcom/ xxx,yyy,zzz,uuu,vvv,www,erg,wgt,tme,vel,dls,
1 dxl,dtc,elc(mipt),fiml(mipt),fismg,wtfasv,rnk,spare(mspare),
2 totmp,zpblcm,
3 xxx9(mpb),yyy9(mpb),zzz9(mpb),uuu9(mpb),vvv9(mpb),www9(mpb),
4 erg9(mpb),wgt9(mpb),tme9(mpb),vel9(mpb),dls9(mpb),dxl9(mpb),
5 dtc9(mpb),elc9(mpb,mipt),fiml9(mpb,mipt),fismg9(mpb),
6 wtfas9(mpb),rnk9(mpb),spare9(mpb,mspare),totmp9(mpb),
7 zpb9cm(mpb),
1 npa,icl,jsu,ipt,iex,node,idx,ncp,jgp,jan,lev,iii,jjj,kkk,iap,
1 iexp,mtp,mpblcm,
3 npa9(mpb),icl9(mpb),jsu9(mpb),ipt9(mpb),iex9(mpb),node9(mpb),
4 idx9(mpb),ncp9(mpb),jgp9(mpb),jan9(mpb),
4 lev9(mpb),iii9(mpb),jjj9(mpb),
5 kkk9(mpb),iap9(mpb),iexp9(mpb),mtp9(mpb),
6 mpb9cm(mpb)
parameter (npblcm=mspare+2*mipt+17,lpblcm=17)
dimension gpblcm(npblcm+1),jpblcm(lpblcm+1),
1 gpb9cm(mpb,npblcm+1),jpb9cm(mpb,lpblcm+1)
equivalence (xxx,gpblcm),(npa,jpblcm),(xxx9,gpb9cm),(npa9,jpb9cm)
general task common.
*if def,unicos.and.multt,1
task common /tskcom/ amfp,ang(3),cbwf,cmult,colout(2,11),cpv,ddet,
*if -def,unicos.or.multt,1
common /tskcom/ amfp,ang(3),cbwf,cmult,colout(2,11),cpv,ddet,
1 colltc(mipt),deb,dti(mlgc),eacctc(mipt,3),racctc(mipt,2),
1 eg0,ergace,paxtc(6,16,mipt),
*if def,RADIOG
+ deltas,deltat,
1 pfp,ple,pmf,psc,qpl,ranb,ranitc,ranjtc,rans,rijktc,rlttc(4,2),
2 rnrtc,rnrtc0,sff(3,maxv),siga,smultc(3),ssr,stp,sumktc(3),
2 tmavtc(mipt,3),totgp1,totm,tpd(7),tpp(20),ttn,udt(10,0:mxlv),
3 udtr(10*mxlv+10),udts(10*mxlv+10),udtt(10*mxlv+10),uold(3),
3 vtr(3),wcs1tc(mipt),wcs2tc(mipt),wgtstc(2),ycn,
4 iax,ibc,ibe,ibs,ibt,ibu,iclp(5,0:mxlv),idet,iet,imd,ipsc,irt,
5 isic(maxv),ital,iti(mlgc),ixcos,ixre,jap,jbd,jbnk,jev,jtls,kdb,
5 jlock,kqss,ktask,levp,lgc(mlgc+1),lsb,mbb,mkc,mpan,nbhwtc,nbnk,
6 nbttc(mipt),nch(mipt),ngp,nlaj,nlse,nlt,npb,npsrtc,
6 npstc,nrnhtc(3),nter,ntii,ntx,ntyn,nziytc(8,mxdx,mipt),
*if def,unicos.and.multt,1
task common /itskpt/ kddm,kddn,kdec,kdxc,kdxd,kfeb,kflx,kgww,
*if -def,unicos.or.multt,1
common /itskpt/ kddm,kddn,kdec,kdxc,kdxd,kfeb,kflx,kgww,
1 kpac,kpan,kpcc,kpwb,kwns,kise,klaj,klcj,klse,kndp,kmaz,
2 kndr,kdrc,kfdd,kgnr,kpik,kptb,krtc,ktgp,kifl,kpc2,kjfl,
3 kjft,kktc,knhs,kshs,kstt,ktal,kgbn,kbnk,mtskpt
parameter (ntskcm=
*if def,RADIOG
+ 107*mipt+40*(mxlv+1)+3*maxv+mlgc+102)
parameter (ltskcm=mipt*(2+8*mxdx)+5*(mxlv+1)+2*mlgc+maxv+47)
parameter (ltskpt=39)
dimension gtskcm(ntskcm),jtskcm(ltskcm),udt1(10*mxlv+10),
1 ktskpt(ltskpt)
equivalence (amfp,gtskcm),(iax,jtskcm),(udt,udt1),(kddm,ktskpt)
********************
dynamically allocated common
********************
*if -def,POINTERS
dac contains the arrays themselves.
common /dac/ das(mdas/ndp2)
*if def,POINTERS
*if def,LP64,1
integer*8 kdy
dac contains only the pointer.
common /dac/ kdy
fixed dynamically allocated common.
dimension aaafd(2),ara(1),asp(1),ava(1),avz(1),
1 awc(1),awn(1),cmg(1),den(1),dptb(3,1),drs(mipt,1),
1 dxcp(0:mxdx,mipt,1),eaa(1),ear(1),eba(mtop,1),
2 ebd(mtop,1),ebl(1),ebt(mtop,1),ech(mpng,mwng,1),edg(1),
+ eee(mipt,1),eek(1),egg(maxi,1),elp(mipt,1),
*if def,MESHTAL
+ emsh(1),
+ esa(1),ewwg(1),fim(mipt,1),flc(1),fme(1),fmg(1),
9 aa9(1),arg(1),esxln(iemaxn,1),esxlp(iemaxp,1),esxmp(1),fmf(1),
9 hsigg(mlpt,1),pcap(1),sigg(msigg,1),sigmx(mlpt,1),zz9(1),
9 xsiso(mxso,1),
4 for(mipt,1),frc(1),
5 fst(1),ftt(1),gmg(1),gvl(1),gwt(1),pbr(1),pbt(1),pkn(1),
6 pmg(1),pru(1),pxr(1),qav(mipt,1),qax(mipt,1),qcn(1),rho(1),
7 rptb(1),scf(1),smg(1),spf(4,2),sqq(12,1),sso(1),tbt(1),
8 tds(1),tmp(1),trf(17,0:1),tth(1),vcl(3,7,1),vec(3,1),
9 vol(1),wgm(1),wwe(1),wwf(1),wwk(1),xnm(1),xse85(10,1)
dimension iiifd(1),ipan(1),iptal(8,5,1),iptb(2+2*npkey,1),
1 kpptb(mipt,1),iazdex(1),mel(1),nel(1),
2 iss(1),itds(1),ixl(3,1),iza(1),jasr(1),jmd(1),jptal(8,1),jscn(1),
3 jss(1),jtf(8,1),jun(1),jvc(1),jxs(32,1),kcp(1),ksc(1),
4 ksd(21,1),ixs(10*maxsec,1),
4 kst(1),ksu(1),ktp(mipt,1),kxs(1),laf(3,3),lat(2,1),
5 lca(1),lfcl(1),lft(mkft,1),lja(1),lme(mipt,1),lmt(1),
6 locct(mipt,1),locph(1),locst(mipt,1),lsc(1),mat(1),mazp(3,1),
7 mazu(1),mbd(1),mbi(1),mfl(3,1),ncl(1),ngmfl(1),nhtfl(1),
8 nmt(1),npq(1),nptb(1),nsb(1),nsf(1),nty(1),nxs(16,1)
*if -def,POINTERS
equivalence (das,aaafd,ara,asp,ava,avz,
1 awc,awn,cmg,den,dptb,drs,dxcp,
+ eaa,ear,eba,ebd,ebl,ebt,ech,edg,eee,eek,egg,elp,
*if def,MESHTAL,1
+ esa,ewwg,fim,flc,fme,fmg,
2 aa9,arg,esxln,esxlp,esxmp,fmf,hsigg,pcap,sigg,sigmx,zz9,xsiso,
2 for,frc,fst,ftt,gmg,gvl,gwt,pbr,pbt,pkn,pmg,pru,
3 pxr,qav,qax,qcn,rho,rptb,scf,smg,spf,sqq,sso,tbt,tds,tmp,trf,
4 tth,vcl,vec,vol,wgm,wwe,wwf,wwk,xnm,xse85,
1 iiifd,ipan,iptal,iptb,kpptb,iazdex,mel,nel,
1 iss,itds,ixl,iza,jasr,jmd,jptal,jscn,
2 jss,jtf,jun,jvc,jxs,kcp,ksc,ksd,ixs,kst,ksu,ktp,kxs,laf,lat,lca,
3 lfcl,lft,lja,lme,lmt,locct,locph,locst,lsc,mat,mazp,mazu,mbd,
4 mbi,mfl,ncl,ngmfl,nhtfl,nmt,npq,nptb,nsb,nsf,nty,nxs)
*if def,POINTERS
pointer (kdy,aaafd),(kdy,ara),(kdy,asp),(kdy,ava),(kdy,avz),
1 (kdy,awc),(kdy,awn),
1 (kdy,cmg),(kdy,den),(kdy,dptb),(kdy,drs),(kdy,dxcp),(kdy,eaa),
2 (kdy,ear),(kdy,eba),(kdy,ebd),(kdy,ebl),(kdy,ebt),(kdy,ech),
+ (kdy,edg),(kdy,eee),(kdy,eek),(kdy,egg),(kdy,elp),
*if def,MESHTAL,1
+ (kdy,emsh),
+ (kdy,esa),(kdy,ewwg),(kdy,fim),(kdy,flc),(kdy,fme),(kdy,fmg),
4 (kdy,aa9),(kdy,arg),(kdy,esxln),(kdy,esxlp),(kdy,esxmp),
4 (kdy,fmf),(kdy,hsigg),(kdy,pcap),(kdy,sigg),(kdy,sigmx),
4 (kdy,zz9),(kdy,xsiso),(kdy,for),
5 (kdy,frc),(kdy,fst),(kdy,ftt),(kdy,gmg),(kdy,gvl),(kdy,gwt),
6 (kdy,pbr),(kdy,pbt),(kdy,pkn),(kdy,pmg),(kdy,pru),(kdy,pxr),
7 (kdy,qav),(kdy,qax),(kdy,qcn),(kdy,rho),(kdy,rptb),(kdy,scf),
8 (kdy,smg),(kdy,spf),(kdy,sqq),(kdy,sso),(kdy,tbt),(kdy,tds),
9 (kdy,tmp),(kdy,trf),(kdy,tth),(kdy,vcl),(kdy,vec),(kdy,vol),
1 (kdy,wgm),(kdy,wwe),(kdy,wwf),(kdy,wwk),(kdy,xnm),(kdy,xse85)
pointer (kdy,iiifd),(kdy,ipan),(kdy,iptal),(kdy,iptb),
1 (kdy,kpptb),(kdy,iazdex),(kdy,mel),(kdy,nel),(kdy,iss),
1 (kdy,itds),(kdy,ixl),(kdy,iza),(kdy,jasr),(kdy,jmd),(kdy,jptal),
2 (kdy,jscn),(kdy,jss),(kdy,jtf),(kdy,jun),(kdy,jvc),(kdy,jxs),
3 (kdy,kcp),(kdy,ksc),(kdy,ksd),
3 (kdy,ixs),(kdy,kst),(kdy,ksu),(kdy,ktp),
4 (kdy,kxs),(kdy,laf),(kdy,lat),(kdy,lca),(kdy,lfcl),(kdy,lft),
5 (kdy,lja),(kdy,lme),(kdy,lmt),(kdy,locct),(kdy,locph),
6 (kdy,locst),(kdy,lsc),(kdy,mat),(kdy,mazp),(kdy,mazu),
7 (kdy,mbd),(kdy,mbi),(kdy,mfl),(kdy,ncl),(kdy,ngmfl),
8 (kdy,nhtfl),(kdy,nmt),(kdy,npq),(kdy,nptb),(kdy,nsb),
9 (kdy,nsf),(kdy,nty),(kdy,nxs)
variable dynamically allocated common.
dimension aaavd(1),ddm(2,1),ddn(24,1),dec(3,1),dxc(3,1),
1 dxd(mipt,24,mxdx),febl(2,1),flx(1),fso(1),gww(2,9,1),
2 pac(mipt,10,1),pan(2,6,1),pcc(3,1),pwb(mipt,19,1),rkpl(19,1),
3 shsd(nspt,1),stt(ntp,1),tfc(6,20,1),wns(2,1),
4 iiivd(1),isef(2,1),jfq(8,0:1),laj(1),lcaj(1),lse(1),maze(1),
5 ndpf(6,1),ndr(1),nhsd(nsp12,1),npsw(1),nsl(2+4*mipt,1),ntbb(4,1)
*if -def,POINTERS
equivalence (das,aaavd,ddm,ddn,dec,dxc,dxd,febl,flx,fso,gww,
1 pac,pan,pcc,pwb,rkpl,tfc,wns,shsd,stt,iiivd,isef,jfq,laj,
2 lcaj,lse,maze,ndpf,ndr,nhsd,npsw,nsl,ntbb)
*if def,POINTERS
pointer (kdy,aaavd),(kdy,ddm),(kdy,ddn),(kdy,dec),(kdy,dxc),
1 (kdy,dxd),(kdy,febl),(kdy,flx),(kdy,fso),(kdy,gww),(kdy,pac),
2 (kdy,pan),(kdy,pcc),(kdy,pwb),(kdy,rkpl),(kdy,shsd),(kdy,stt),
+ (kdy,tfc),(kdy,wns),
4 (kdy,iiivd),(kdy,isef),(kdy,jfq),(kdy,laj),(kdy,lcaj),(kdy,lse),
5 (kdy,maze),(kdy,ndpf),(kdy,ndr),(kdy,nhsd),(kdy,npsw),(kdy,nsl),
6 (kdy,ntbb)
ephemeral dynamically allocated common.
dimension scr(1),
*if def,RADIOG,1
+ drc(18,1),
*if -def,RADIOG,1
+ drc(16,1),
*if def,MESHTAL
+ gdata(1),gd1(1),gd2(1),gd3(1),
+ emi(1),fdd(2,1),genr(1),pik(1),ptr(1),ptb(4,1),pts(1),
+ rng(mipt,1),rtc(10,1),tgp(1),ifl(1),ipac2(1),ixc(61,1),jfl(1),
2 jft(1),jmt(1),kmt(3,1),ktc(2,1),kxd(1),lbb(1)
*if -def,POINTERS
equivalence (das,scr,drc,
*if def,MESHTAL,1
+ gdata,gd1,gd2,gd3,
+ emi,fdd,genr,pik,ptb,ptr,pts,rng,rtc,tgp,
1 ifl,ipac2,ixc,jfl,jft,jmt,kmt,ktc,kxd,lbb)
*if def,POINTERS
pointer (kdy,scr),(kdy,drc),
*if def,MESHTAL,1
+ (kdy,gdata),(kdy,gd1),(kdy,gd2),(kdy,gd3),
+ (kdy,emi),(kdy,fdd),(kdy,genr),
1 (kdy,pik),(kdy,ptb),(kdy,ptr),(kdy,pts),(kdy,rng),(kdy,rtc),
2 (kdy,tgp),(kdy,ifl),(kdy,ipac2),(kdy,ixc),(kdy,jfl),(kdy,jft),
3 (kdy,jmt),(kdy,kmt),(kdy,ktc),(kdy,kxd),(kdy,lbb)
tallies, bank, and cross-sections in dac.
dimension tal(1),gbnk(1),ibnk(1),exs(1)
*if def,XS64,1
dimension xss(1)
*if -def,XS64,1
real xss(1)
*if -def,POINTERS,1
equivalence (das,tal,gbnk,ibnk,xss,exs)
*if def,POINTERS,1
pointer (kdy,tal),(kdy,gbnk),(kdy,ibnk),(kdy,xss),(kdy,exs)
dynamically allocated common for imcn.
dimension jtr(1),awt(1),bbv(1),prb(1),rtp(1),sfb(1),
1 ipnt(2,mktc,0:1),jasw(1),kaw(1),kdr(1),kdup(mipt+2,1),
1 kmm(1),ktr(1),
2 lxd(mipt,1),mfm(1),nlv(1),nslr(2+4*mipt,1),
3 aras(2,1),atsa(2,1),rscrn(2,1),rsint(2,1),scfq(5,1),vols(2,1),
4 iint(1),icrn(3,1),ljav(1),ljsv(1),lsat(1)
*if -def,POINTERS
equivalence (das,jtr,awt,bbv,prb,rtp,sfb,ipnt,jasw,kaw,kdr,kdup,
1 kmm,ktr,lxd,mfm,nlv,nslr,
2 aras,atsa,rscrn,rsint,scfq,vols,iint,icrn,ljav,ljsv,lsat)
*if def,POINTERS
pointer (kdy,jtr),(kdy,awt),(kdy,bbv),(kdy,prb),(kdy,rtp),
1 (kdy,sfb),(kdy,ipnt),(kdy,jasw),(kdy,kaw),(kdy,kdr),(kdy,kdup),
2 (kdy,kmm),(kdy,ktr),(kdy,lxd),(kdy,mfm),(kdy,nlv),(kdy,nslr),
3 (kdy,aras),(kdy,atsa),(kdy,rscrn),(kdy,rsint),(kdy,scfq),
4 (kdy,vols),(kdy,iint),(kdy,icrn),(kdy,ljav),(kdy,ljsv),(kdy,lsat)
dynamically allocated common for plot.
dimension amx(4,4,1),coe(6,2,1),crs(1),jst(2,1),kcl(102,1),kfm(1),
1 lcl(1),lsg(1),ncs(1),plb(1),qmx(3,3,2,1),zst(1)
*if -def,POINTERS,1
equivalence (das,amx,coe,crs,jst,kcl,kfm,lcl,lsg,ncs,plb,qmx,zst)
*if def,POINTERS
pointer (kdy,amx),(kdy,coe),(kdy,crs),(kdy,jst),(kdy,kcl),
1 (kdy,kfm),(kdy,lcl),(kdy,lsg),(kdy,ncs),(kdy,plb),(kdy,qmx),
2 (kdy,zst)
dynamically allocated common for mcplot.
dimension ab1(1),ab2(1),erb(1),kxspie(1),kxspen(1),kxspxs(1),
1 kxspnx(1),mcc(1),ord(1),xcc(1),ycc(1)
real xrr(1),yrr(1)
*if -def,POINTERS
equivalence (das,ab1,ab2,erb,kxspie,kxspen,kxspxs,kxspnx,
1 mcc,ord,xcc,xrr,ycc,yrr)
*if def,POINTERS
pointer (kdy,ab1),(kdy,ab2),(kdy,erb),(kdy,kxspie),(kdy,kxspen),
1 (kdy,kxspxs),(kdy,kxspnx),(kdy,mcc),(kdy,ord),
1 (kdy,xcc),(kdy,xrr),(kdy,ycc),(kdy,yrr)
*comdeck gs
common block for the gks simulation subroutines.
common /gkssim/ als,chite(5),chup(2),psize(4),rgb(100),wsf,xlf,
1 xrt,ybt,ytp,icol,jta(2),kls,lvd,ltype,munit,nctext,nltext,
2 npages,nxnorm,nynorm
c ----------------------------------------------------------------------
*comdeck jc
common blocks for imcn
V2.1.4: Eliminates the za, zb, and zc cards, reducing nkcd by 3.
HISTP Implementation: nkcd includes 1 for histp card.
RADIOG Implementation: nkcd includes 3 for the following cards:
fi, notrn, and talnp.
MESHTAL Implementation: nkcd includes 7 for the following cards:
rmesh, cmesh, cora, corb, corc, ergsh, and mshmf.
parameter (nkcd=94,ntalmx=100,mopts=5)
common /imccom/ ajsh,bbb(4,4),fes(33),ritm,swtm,swtx,
1 ica,icn,icx,ifip(mipt),iitm,ioid,ipl,irc,irs,jsd(4,33),jui,kitm,
2 krq(7,nkcd),ktl(ntalmx,2),licc,likef,lit,lrt,ltd,m1c,m2c,m3c,m4c,
2 m5c,m6c,m7c,m8c,m9c,m10c,
*if def,RADIOG
3 mipts,mkcp,mlaf,mssc,mxit,mxmel1,
3 naw,ncomp,ncpar(mipt,nkcd),ncparf,ncrn,
4 ndp(ntalmx),ndup(3),nii,nitm,niwr,njsw,njsx,nmat1,nmfm,
5 novol,nqp(mipt),nqw,nsc,nsjv,nsv,ntl(0:ntalmx),
6 ktt(ntalmx,mipt),nwc,nxsc
offsets for virtual arrays in dynamically allocated storage.
common /imccom/ ljtr,lawt,lbbv,lprb,lrtp,lsfb,lipn,ljaw,lkaw,lkdr,
1 llxd,lmfm,ldup,lkmm,lktr,lnlv,lnsr,lars,lats,lrsc,lrsn,lscq,
2 lvls,liin,licr,llav,llsv,llsa
character common
character cnm(nkcd)*5,hitm*80,hlin*80,hmopt(mopts)*5,hpbl(24)*7,
1 hptb(npkey)*7,hptr(nptr)*7,ich*5
common /jmccom/ cnm,hitm,hlin,hmopt,hpbl,hptb,hptr,ich
c ----------------------------------------------------------------------
*comdeck lkoff
turn off multitasking lock.
*if -def,MULTT,1
*if def,MULTT
*if def,UNICOS,1
if(jlock.eq.1)call lockoff(locki)
*if def,VMS,3
if(jlock.eq.1)then
cpar$ lockoff lockl
if(jlock.eq.1)jlock=-1
*comdeck lkon
turn on multitasking lock.
*if -def,MULTT,1
*if def,MULTT
*if def,UNICOS,1
if(jlock.eq.-1)call lockon(locki)
*if def,VMS,3
if(jlock.eq.-1)then
cpar$ lockon lockl
if(jlock.eq.-1)jlock=1
*comdeck lm
common blocks for electron landau treatment
parameter (mlam=5001)
common /lancom/ eqlm(mlam)
c ----------------------------------------------------------------------
*comdeck lt
common blocks for electron data
parameter (mlanc=1591)
common /lancut/ avlm(mlanc),flam(mlanc)
c ----------------------------------------------------------------------
*comdeck mb
common block for message passing subroutines
common /msgcom/ idbuf,itid(mcpu),nsub
c ----------------------------------------------------------------------
*comdeck mp
common blocks for mcplot.
parameter (isbm=100,jsbm=100,mclevs=20,nkeys=57)
save /mplcom/,/zchar/
common /mplcom/ abhi(2),ablo(2),clev(mclevs),contur(3),cthick,
1 dissf(3),pht(2),plim(4),scalf(2,3),tensn,xlg,xst,xyzmn(3),
2 xyzmx(3),ylg,yst,yval,ibl(8,2),icut(2),ifree(2),ipct,iprpts,
3 itfc,itik(2),ititle(7),j3d,jlbl(2,8),jlim(2),kbin(8,2),klbl(43),
4 koplot,kurv,kxsplt,kxspma,kxspmt,kxspar,kxsmat,kxspkm,kxsptp,
5 kxspu(43),lax,legalc(nkeys),legalm(nkeys),legalx(nkeys),legend,
6 lnstyl,lpert,lput,mscal,mtal,nclev,noerbr,nfree,nonorm,npt(2)
offsets for virtual arrays in dynamically allocated storage.
common /mplcom/ lab1,lab2,lerb,lmcc,lord,lxcc,lxie,lxen,
1 lxnx,lxrr,lxxs,lycc,lyrr
character common.
character hlbl(43)*40,hxspu(15)*40,keys(nkeys)*8,titles(7)*40,
1 xspttl*10
common /zchar/ hlbl,hxspu,keys,titles,xspttl
c ----------------------------------------------------------------------
*comdeck pc
common blocks for plot
parameter (mclb=19,mplm=100,nkeyp=28,res=1./1500.)
save /pltcom/,/qltcom/
common /pltcom/ basis(9),exsav(2),extent(2),origin(3),orsav(3),
1 plmx(4,4),pxx(4,4),sch,sclabl(4),
2 icolor(mplm),icurs,icurs1,iqc,jloc,jscal,lcolor,ncel,ncrs,nlb,
3 nplb,nxp
offsets for virtual arrays in dynamically allocated storage.
common /pltcom/ lamx,lcoe,lcrs,ljst,lkcl,lkfm,llcl,llsg,lncs,lplb,
1 lqmx,lzst
character common.
character keyp(nkeyp)*8,hcolor(ncolor)*12
common /qltcom/ keyp,hcolor
c ----------------------------------------------------------------------
c ----------------------------------------------------------------------
c___________________________________________________________________
c The PRPR preprocessable deck is prepared from the *.F and *.h files
c obtained from LANL, by P. Vertes, KFKI Atomic Energy Research
c Institute, Budapest, Hungary, September, 2000.
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
*deck av abvals
*if def,MCPLOT
subroutine abvals(la,iv)
extract npt(iv) abscissa values of free variable ifree(iv) from
tds and put them in ab1(la).
la may be lab1 or lab2.
branch for an xs plot.
if(kxsplt.ne.0)go to 120
use bin numbers if the variable is not c, e, or t.
if(ifree(iv).ge.6)go to 20
*if def,RADIOG,1
if(ifree(iv).eq.4.and.jptal(ljpt+8,mtal).gt.2)go to 20
do 10 i=1,npt(iv)
10 ab1(la+i)=i
icut(iv)=1
abhi(iv)=npt(iv)+1
ablo(iv)=0.
use upper bin bounds for c, e, and t.
20 do 30 i=1,npt(iv)
30 ab1(la+i)=tds(iptal(lipt+ifree(iv),1,mtal)+i)
abhi(iv)=ab1(la+npt(iv))
define left insertion point.
ablo(iv)=0.
if(ifree(iv).eq.6)ablo(iv)=-1.
if(ifree(iv).ne.7.or.ktfile.ne.1)go to 50
ablo(iv)=huge
do 40 i=1,mipt
40 if(ktp(lktp+i,mtal).ne.0)ablo(iv)=min(ablo(iv),ecf(i))
50 if(ifree(iv).eq.8)ablo(iv)=min(zero,ab1(la+1))
*if def,RADIOG,2
if((ifree(iv).eq.4.or.ifree(iv).eq.6).and.jptal(ljpt+8,mtal).gt.2)
+ ablo(iv)=tds(iptal(lipt+ifree(iv),1,mtal))
if(nfree.eq.1.and.lax.ge.3.and.ablo(iv).le.0.)
1 ablo(iv)=max(zero+1e-14,.00001*ab1(la+1))
determine position of cutoff.
do 60 i=1,npt(iv)
icut(iv)=i
60 if(ablo(iv).lt.ab1(la+i))go to 80
write(jtty,70)
70 format(38h cutoff is larger than all bin values.)
normalize 2d ordinate values by bin width.
80 if(ktfile.ne.1.and.iptal(lipt+ifree(iv),2,mtal).ne.0)return
if(nfree.eq.2)go to 100
*if def,RADIOG,2
if(jptal(ljpt+2,mtal).eq.5.and.(ifree(iv).eq.4.or.ifree(iv).eq.6))
+ go to 92
if(nonorm.ne.0)go to 100
do 90 i=icut(iv),npt(iv)
if(i.eq.icut(iv))d=ab1(la+i)-ablo(iv)
if(i.ne.icut(iv))d=ab1(la+i)-ab1(la+i-1)
if(ifree(iv).eq.6)d=2.*pie*d
90 if(d.ne.0.)ord(lord+i)=ord(lord+i)/d
*if def,RADIOG,16
if nonorm.ne.0 multiply ordinate by bin area of grid.
92 if(nonorm.eq.0)go to 100
if(ifree(iv).eq.6)ib=4
if(kbin(ib,1).eq.iptal(lipt+ib,4,mtal))go to 94
d1=tds(iptal(lipt+ib,1,mtal)+kbin(ib,1))-
1 tds(iptal(lipt+ib,1,mtal)+kbin(ib,1)-1)
94 d1=tds(iptal(lipt+ib,1,mtal)+iptal(lipt+ib,3,mtal))-
1 tds(iptal(lipt+ib,1,mtal))
96 do 98 i=icut(iv),npt(iv)
if(i.eq.icut(iv))d=d1*(ab1(la+i)-ablo(iv))
if(i.ne.icut(iv))d=d1*(ab1(la+i)-ab1(la+i-1))
98 ord(lord+i)=d*ord(lord+i)
bin center the abscissa values.
100 do 110 iz=icut(iv)+1,npt(iv)
i=npt(iv)+icut(iv)+1-iz
110 ab1(la+i)=.5*(ab1(la+i)+ab1(la+i-1))
ab1(la+icut(iv))=.5*(ab1(la+icut(iv))+ablo(iv))
multigroup plots.
npt(iv)=2*ne.
120 if(kxsptp.ne.6.and.kxsptp.ne.7)go to 140
ne=nxs(lnxs+5,kxspie(lxie+1))
do 130 i=1,ne
kp=kxspen(lxen+1)+i
ab1(la+npt(iv)-2*i+1)=xss(kp)-.5*xss(kp+ne)
130 ab1(la+npt(iv)-2*i+2)=xss(kp)+.5*xss(kp+ne)
create energy grid for non-multigroup xs plots.
if(kxspar.lt.3.or.kxsptp.ne.8)go to 160
set up grid for electrons.
in=nee(3)-1
do 150 j=1,in
150 ab1(la+j)=eee(leee+3,k)
get energy grid for single nuclide case.
if(kxsptp.eq.4.and.nxs(lnxs+5,kxspie(lxie+1)).eq.4.and.
1 kxspmt.le.2)kf=2
if(kxsplt.gt.1)go to 200
if(kxsptp.eq.2.or.kf.ne.0)go to 180
npt(iv)=kxspnx(lxnx+1)
do 170 i=1,npt(iv)
ab1(la+i)=xss(kxspen(lxen+1)+i)
170 if(kxspar.eq.2)ab1(la+i)=exp(xss(kxspen(lxen+1)+i))
drxs or s(a,b) elastic / energy.
180 npt(iv)=2*kxspnx(lxnx+1)
do 190 i=1,npt(iv)
190 ab1(la+i)=xss(kxspen(lxen+1)+(i+1)/2)
multiple nuclides on different grids.
loop through nuclides then through energies.
set up dummy pointer in upper half of kxspnx array.
lower half of kxspnx array has number of points in xss array.
200 do 210 j=1,kxsplt
210 kxspnx(lxnx+mxe1+j)=1
find next minimum energy of given isotopes.
eo=-huge for photons.
220 en=huge
do 230 j=1,kxsplt
if(kxspen(lxen+j).le.0.or.kxspxs(lxxs+j).le.0)go to 230
if(kxspnx(lxnx+mxe1+j).gt.kxspnx(lxnx+j))go to 230
e=xss(kxspen(lxen+j)+kxspnx(lxnx+mxe1+j))
if(e.gt.en)go to 230
if(e.eq.en)kxspnx(lxnx+mxe1+k)=kxspnx(lxnx+mxe1+k)+1
230 continue
if(en.eq.huge)go to 240
kxspnx(lxnx+mxe1+k)=kxspnx(lxnx+mxe1+k)+1
if(en.le.eo)go to 220
ab1(la+in)=en
if(kxspar.eq.2)ab1(la+in)=exp(en)
kxsptp=nty(lnty+kxspie(lxie+k))
drxs or s(a,b) elastic / energy.
if(kxsptp.ne.2.and.kf.ne.k)go to 220
ab1(la+in)=en-1.e-5*en
ab1(la+in)=en
240 npt(iv)=in
determine left insertion point and set upper energy.
250 abhi(iv)=ab1(la+npt(iv))
do 260 i=1,npt(iv)
260 if(lax.lt.3.or.ab1(la+i).gt.0.)go to 270
270 icut(iv)=i
ablo(iv)=ab1(la+i)
subroutine acecas(ls,ip,q,ia,ka,id,kd)
sample the emission parameters from the appropriate law data.
the subroutine takes in all the information necessary to
sample (or debug inability to sample) the emission
energy and scattering angle in the laboratory coordinate
Preconditions:
explicitly passed in variables:
ls - the current particle index in the colout array
ip - the ipt particle type for the incident particle
q - the q value for the reaction being sampled
ia - the first word of the relevant AND block in the xss array
ka - the offset to the first word of the table in AND block
id - the first word of the relevant DLW block in the xss array
kd - the offset to the first word of the table in DLW block
implicitly uses variables:
awn(lawn+iex) - the awr of the current target isotope
colout(1,ls) - the sampled emission energy
colout(2,ls) - the sampled emission scattering angle
erg - the incident particle energy in the lab system
gpt - array of particle awr's
iex - the table index of the current target isotope
ipsc - index to indicate what kind of law was sampled
ipt - the ipt particle type of the emitted particle
ixl - the ZAID storage array
ixre - the current reaction index being sampled
kdb - fatal error flag
mtp - the reaction mt number being sampled
ntyn - the coordinate system of the sampled parameters
tpd(1,2) - storage for correlated energy/angle parameters
wgt - the weight of the emitted particle
xss - data array containing sampling distributions etc.
read only input variables:
ls, ip, q, ia, ka, id, kd, awn, erg, gpt, iex, ipt, ixl,
ixre, mtp, ntyn and all xss(i)
Postconditions:
returns the sample emission parameters in the lab system:
colout(1,ls) - the emission energy
colout(2,ls) - the emission scattering angle
Law 4/44/61 makes use of a biased distribution which affects
the outgoing particle weight
kdb is set on fatal error (inability to sample reasonable
emission parameters) during call to expirx
modified variables: colout(1,ls), colout(2,ls), ipsc, kdb,
tpd and wgt
parameter (ep=0.000001)
character ht*10
c***********************************************************************
select the law.
10 t1=rang()
20 n=id+nint(xss(n-1))
25 if(nint(xss(n-1)).eq.0)go to 30
t1=t1-acefcn(n+2,erg,ln)
if(t1.ge.0.)go to 20
c***********************************************************************
use the selected law to sample the energy [and possibly angle].
if law samples without error, go to sample angle or
coordinate transform as appropriate
30 colout(1,ls)=-huge
lw=nint(xss(n))
iw=id-1+nint(xss(n+1))
if(lw.eq.1) go to 40
if(lw.eq.2) go to 60
if(lw.eq.3) go to 70
if(lw.eq.4) go to 80
if(lw.eq.5) go to 160
if(lw.eq.7) go to 170
if(lw.eq.9) go to 190
if(lw.eq.11)go to 210
if(lw.eq.22)go to 230
if(lw.eq.24)go to 242
if(lw.eq.33)go to 70
if(lw.eq.44)go to 80
if(lw.eq.61)go to 80
if(lw.eq.66)go to 245
if(lw.eq.67)go to 255
law 1 (From ENDF Law 1) -- tabular equi-probable energy bins.
40 call acetbl(iw,ic,r,ln)
nt=nint(xss(iw+ln))
iw=iw+ln+nt*(ic-1)
k=int(rang()*(nt-1)+1)
if(r.ne.0.)go to 50
sample from single table.
colout(1,ls)=xss(iw+k)+fr*(xss(iw+k+1)-xss(iw+k))
sample by scaled interpolation between tables.
50 t1=xss(iw+1)+r*(xss(iw+1+nt)-xss(iw+1))
if(rang().le.r)i=i+nt
colout(1,ls)=t1+(xss(i+k)+fr*(xss(i+k+1)-xss(i+k))-xss(i+1))*
1 (xss(iw+nt)+r*(xss(iw+2*nt)-xss(iw+nt))-t1)/(xss(i+nt)-xss(i+1))
law 2 -- discrete photon lines.
this law applies to photon production only and should
not be used by any other emission particle type
60 if(ipt.ne.2)go to 300
colout(1,ls)=xss(iw+1)
if(nint(xss(iw)).eq.2)colout(1,ls)=colout(1,ls)+
1 erg*awn(lawn+iex)/(awn(lawn+iex)+1.)
law 3 & 33 (From ENDF Law 3) -- level scattering.
70 colout(1,ls)=xss(iw+1)*(erg-xss(iw))
law 4 (From ENDF Law 1) -- continuous erg tabular distribution
law 44 (From ENDF Law 1) -- Kalbach-87 correlated formalism
law 61 (From ENDF Law 1) -- correlated tab energy-angle dist
80 call acetbl(iw,ic,r,ln)
nr=nint(xss(iw))
lb=iw-1+ln+ic
lc=id+nint(xss(lb))
xss(lc,ld,lf) is an overloaded variable.
it contains the
number of points in the emission distribution and if part
of the continuum has been expunged, it contains 0.5 times
the cumulative probability of the portion expunged.
np=int(xss(lc)+ep)
jj=nint(xss(lc-1))
if(jj.lt.9999)nd=jj/10
if(r.eq.0.)go to 90
ld=id+nint(xss(lb+1))
mp=int(xss(ld)+ep)
t1=xss(lc+nd+1)+r*(xss(ld+nd+1)-xss(lc+nd+1))
t2=xss(lc+np)+r*(xss(ld+mp)-xss(lc+np))
if(ra.ge.r)go to 90
jj=nint(xss(ld-1))
if(jj.lt.9999.and.jj/10.ne.nd)call expirx(1,'acecas',
1 'wrong number of discrete lines for law 4/44')
if(jj.gt.10000)nd=0
90 if(jj.eq.9999)return
jj=jj-nd*10
check for a hit in the discrete part.
use histogram or corresponding-point interpolation.
if(nd.eq.0)go to 97
do 95 ih=1,nd
cc=xss(lc+2*np+ih)+r*(xss(ld+2*mp+ih)-xss(lc+2*np+ih))
95 if(cc.ge.r1)go to 96
if(np.eq.nd.and.mp.eq.nd)go to 152
96 t=xss(lc+ih)+r*(xss(ld+ih)-xss(lc+ih))
if(t.lt.0.)t=erg*awn(lawn+iex)/(awn(lawn+iex)+1.)-t
if(lw.ne.44)go to 152
tpd(1)=xss(lc+ih+3*np)+r*(xss(ld+ih+3*mp)-xss(lc+ih+3*np))
tpd(2)=xss(lc+ih+4*np)+r*(xss(ld+ih+4*mp)-xss(lc+ih+4*np))
handle a hit in the continuous part.
use histogram or unit-base interpolation inside table.
use scaled interpolation between tables.
97 np=int(xss(lf)+ep)
if(nd.ne.0)r1=1.-(1.-r1)*(1.-xss(lf+2*np+nd))/(1.-cc)
adjust for the energy cutoff if necessary.
see subroutine expung for law 4/44 distribution. 0 < wc >>>>
law 5 (From ENDF Law 5) -- general evaporation spectrum.
160 t1=acefcn(iw,erg,ln)
i=iw+ln+1+int(rang()*(nint(xss(iw+ln))-1))
colout(1,ls)=t1*(xss(i)+rang()*(xss(i+1)-xss(i)))
law 7 (From ENDF Law 7) -- simple Maxwell fission spectrum.
170 t1=acefcn(iw,erg,ln)
t3=erg-xss(iw+ln)
if(t3.le.0.)go to 260
180 t4=rang()**2
t2=t4+rang()**2
if(t2.gt.1.)go to 180
t2=log(rang())*t4/t2
colout(1,ls)=-t1*(t2+log(rang()))
reject if outside range 0 ... e-u
if(colout(1,ls).gt.t3)go to 180
law 9 (From ENDF Law 9) -- evaporation spectrum.
190 t1=acefcn(iw,erg,ln)
t2=erg-xss(iw+ln)
if(t2.le.0.)go to 260
200 fr=rang()
colout(1,ls)=-t1*log(fr*rang())
reject if outside range 0 ... e-u
if(colout(1,ls).gt.t2)go to 200
law 11 (From ENDF Law 11) -- energy dependent Watt spectrum.
210 t1=acefcn(iw,erg,ln)
t2=acefcn(iw+ln,erg,lb)
if(erg.le.xss(iw+ln+lb))go to 260
t5=sqrt((1.+.125*t1*t2)**2-1.)+1.+.125*t1*t2
220 t=-log(rang())
colout(1,ls)=t1*t5*t
if(((1.-t5)*(1.+t)-log(rang()))**2.gt.t2*colout(1,ls))go to 220
law 22 (From UK Law 2) -- tabular linear functions.
230 call acetbl(iw,ic,r,ln)
ie=id-1+nint(xss(iw+ln+ic-1))
nf=nint(xss(ie))
240 iw=iw+1
fr=fr-xss(iw)
if(fr.ge.0.)go to 240
if(iw.gt.ie+nf)go to 235
colout(1,ls)=xss(iw+2*nf)*(erg-xss(iw+nf))
law 24 (From UK Law 6) -- tabular energy multipliers.
242 call acetbl(iw,ic,r,ln)
i=iw+ln+1+nint(xss(iw+ln))*(ic-1)+int(rang()*(nint(xss(iw+ln))-1))
colout(1,ls)=erg*(xss(i)+rang()*(xss(i+1)-xss(i)))
law 66 (From ENDF Law 6) -- n-body phase space distribution.
245 nb=nint(xss(iw))
ap=xss(iw+1)
if(ipt.gt.2)ap=ap*gpt(1)/gpt(ipt)
246 r=rang()**2
s=r+rang()**2
if(s.gt.1.)go to 246
x=-r*log(s)/s-log(rang())
248 r=rang()**2
s=r+rang()**2
if(s.gt.1.)go to 248
go to(252,251,250)nb-2
250 p=p*rang()*rang()
251 p=p*rang()
252 y=-r*log(s)/s-log(p)
aw=awn(lawn+iex)
colout(1,ls)=t*((ap-1.)/ap)*(erg*aw/(aw+1.)+q)
colout(2,ls)=2.*rang()-1.
if(ntyn.ge.0)go to 295
law 67 (endf/b-vi law 7) -- correleted energy-angle scatter.
255 call acetbl(iw,ic,r,ln)
cs=acecos(ixcos,ia,ka)
colout(2,ls)=cs
colout(1,ls)=acecs6(0,id,iw,ic,r,cs)
if(ntyn.le.0)go to 295
if(colout(1,ls).lt.0.)go to 300
c***********************************************************************
if not correlated energy-angle, calculate the cosine.
260 colout(2,ls)=acecos(ixcos,ia,ka)
c***********************************************************************
adjust if energy and cosine are given in center-of-mass system.
280 if(colout(1,ls).lt.0.)go to 300
ergace=colout(1,ls)
if(ntyn.ge.0)go to 290
Formulas below are from p. 2 of X-6:RES-93-68.
These formulas assume two-body kinematics.
Seamon's formulas specify atomic weight ratios (to neutron)
For incident particle (a), AWR=GPT(IPT_INCIDENT)/GPT(1)
For exiting particle (b), AWR=GPT(IPT)/GPT(1)
For target (A), AWR=AWN(LAWN+IEX)
a1=gpt(ip)/gpt(1)
a2=gpt(ipt)/gpt(1)
a3=awn(lawn+iex)
t1=colout(1,ls)
t2=erg*a1*a2/(a3+a1)**2
t3=2.*sqrt(a2*a1*erg*colout(1,ls))*colout(2,ls)/(a3+a1)
t4=t1+t2+t3
s1=colout(2,ls)*sqrt(colout(1,ls)/t4)
s2=sqrt(a1*a2*erg/t4)/(a3+a1)
colout(1,ls)=t4
colout(2,ls)=s1+s2
c***********************************************************************
resample energy up to 100 times if > emx(1).
290 if(colout(1,ls).le.emx(ipt))return
call errprn(1,netb(2),4,erg,zero,'erg',' ',
1 'energy of particle from inelastic collision > emx')
if(nx.lt.100)go to 10
call expirx(1,'acecas','erg>emx 100 times in one collision.')
c***********************************************************************
print debug information for cross-section table errors.
295 colout(1,ls)=huge
300 call zaid(2,ht,ixl(lixl+1,iex))
*call lkon
write(iuo,310)ht,erg,ixre,mtp,ntyn,lw,colout(1,ls)
310 format(/30h error in cross-section table ,a10/12h energy in =,
1 1pe12.4,5x,16hreaction index =,i3,5x,4hmt =,i4,5x,4hty =,i4,
2 5x,5hlaw =,i3,5x,12henergy out =,1pe12.4)
if(colout(1,ls).eq.-huge)call expirx(1,'acecas',
1 'an inappropriate or non-existent law was selected.')
if(colout(1,ls).lt.0.)call expirx(1,'acecas',
1 'emission energy was negative.')
if(colout(1,ls).eq.huge)call expirx(1,'acecas',
1 'faulty cross-section data.')
if(colout(1,ls).gt.erg)call expirx(1,'acecas',
1 'emission energy exceeds incident energy.')
*deck ac acecol
subroutine acecol(jq)
sample an elastic or inelastic neutron collision.
jq=0 normally.
jq=2 for kcode sampling of fission energy only.
character ht*10
10 if(jq.eq.2)go to 80
sample elastic vs. inelastic collision.
if(nxs(lnxs+5,iex).eq.0)go to 60
ic=ktc(kktc+2,iex)+jxs(ljxs+1,iex)-1
n=nxs(lnxs+3,iex)
el=xss(ic+3*n)+rtc(krtc+2,iex)*(xss(ic+3*n+1)-xss(ic+3*n))
if(erg.ne.eg0)go to 20
sr=rtc(krtc+4,iex)-rtc(krtc+3,iex)-rtc(krtc+8,iex)-el
20 sr=xss(ic+n)-xss(ic+2*n)+rtc(krtc+2,iex)*
1 (xss(ic+n+1)-xss(ic+n)-xss(ic+2*n+1)+xss(ic+2*n))-el
if(lfcl(llfc+icl).eq.0)go to 30
l=jxs(ljxs+21,iex)
if(l.eq.0)go to 30
j=l+2+ktc(kktc+2,iex)-nint(xss(l))
sr=sr-xss(j)-rtc(krtc+2,iex)*(xss(j+1)-xss(j))
30 if(awn(lawn+iex)*erg.gt.500.*tbt(ltbt+iex))go to 50
a2=awn(lawn+iex)*erg/tbt(ltbt+iex)
if(a2.lt.4.)go to 40
el=el*a2/(a2+.5)
40 a=sqrt(a2)
el=el*a/(thgf(i)+(b-i)*(thgf(i+1)-thgf(i)))
50 r=rang()*(sr+el)-el
if(r.ge.0.)go to 100
****************************
elastic case
****************************
sample the neutron output direction and calculate its energy.
c=acecos(ixcos,jxs(ljxs+9,iex),nint(xss(jxs(ljxs+8,iex))))
if(awn(lawn+iex).lt.1.)go to 70
t1=1.+awn(lawn+iex)*(awn(lawn+iex)+2.*c)
colout(1,1)=erg*t1/(1.+awn(lawn+iex))**2
colout(2,1)=(1.+awn(lawn+iex)*c)/sqrt(t1)
special calculation for hydrogen.
70 colout(1,1)=.5*erg*(1.+c)
colout(2,1)=sqrt(.5+.5*c)
**************************
inelastic case.
***************************
sample the reaction, using cumulative partial cross sections.
80 if(erg.ne.eg0)go to 90
r=rang()*rtc(krtc+8,iex)
90 l=jxs(ljxs+21,iex)
j=l+2+ktc(kktc+2,iex)-nint(xss(l))
r=rang()*(xss(j)+rtc(krtc+2,iex)*(xss(j+1)-xss(j)))
do 120 ixre=1,nxs(lnxs+5,iex)
if(lfcl(llfc+icl).eq.0)go to 110
if(jq.ne.2.eqv.nint(xss(jxs(ljxs+5,iex)+ixre-1)).eq.19)go to 120
calculate the cross section for reaction whose index is ixre.
ie,ne,(xs(i),i=1,ne)
xs(1) corresponds to es(ie).
110 is=jxs(ljxs+7,iex)+nint(xss(jxs(ljxs+6,iex)+ixre-1))
ic=ktc(kktc+2,iex)+1-nint(xss(is-1))
if(ic.lt.1.or.ic.gt.nint(xss(is)).or.ic.eq.nint(xss(is)).and.
1 rtc(krtc+2,iex).ne.0.)go to 120
t=xss(ic+is)
if(rtc(krtc+2,iex).ne.0.)t=t+rtc(krtc+2,iex)*
1 (xss(ic+is+1)-xss(ic+is))
if(t.eq.0.)go to 120
if(t.gt.tt)kx=ixre
tt=max(t,tt)
if(r.lt.0.)go to 150
120 continue
handle failure of the reaction cross sections to add up.
if(kx.eq.0)go to 60
call zaid(2,ht,ixl(lixl+1,iex))
call errprn(1,npnm,4,erg,zero,'erg',' ',
1 'no reaction mt found. collision resampled. zaid = '//ht)
get the reaction type number.
150 ntyn=nint(xss(jxs(ljxs+5,iex)+ixre-1))
mtp=nint(xss(jxs(ljxs+3,iex)+ixre-1))
q=xss(jxs(ljxs+4,iex)+ixre-1)
ia=jxs(ljxs+9,iex)
ka=nint(xss(jxs(ljxs+8,iex)+ixre))
id=jxs(ljxs+11,iex)
kd=nint(xss(jxs(ljxs+10,iex)+ixre-1))
check that the energy is within the band for this reaction.
if(nty(lnty+iex).ne.2)go to 160
l=jxs(ljxs+11,iex)+2+nint(xss(jxs(ljxs+10,iex)+ixre-1))
ie=2*nint(xss(l))+l
if(erg.le.xss(ie+2).or.erg.ge.xss(ie+1+nint(xss(ie+1))))go to 180
sample the energy and scattering angle of emerging neutrons.
if(abs(ntyn).ge.100)go to 190
cmult=abs(ntyn)
if(ntyn.eq.19)cmult=aint(acenu()+rang())
if(jq.eq.0)ns=cmult
do 170 i=1,ns
call acecas(i,1,q,ia,ka,id,kd)
170 if(kdb.ne.0)return
try again, up to 100 times, if the energy is out of the
reaction band.
180 if(erg.le.xss(jxs(ljxs+1,iex)).or.erg.gt.xss(jxs(ljxs+1,iex)+
1 nxs(lnxs+3,iex)-1))go to 160
ndr(kndr+iex)=ndr(kndr+iex)+1
if(nt.le.100)go to 10
call expirx(1,'acecol',
1 'energy was not within the band for the reaction 100 times.')
high energy reaction other than fission with energy-dependent
multiplicity.
kalbach-87 (law 44) endf/b-vi.
the location of the multiplicity table relative to dlw is
abs(ntyn)-100.
acefcn gets non-fission reaction multiplicity.
190 l=jxs(ljxs+11,iex)+abs(ntyn)-101
cmult=acefcn(l,erg,ln)
if(cmult.lt.1.)go to 210
cmult=aint(cmult+rang())
do 200 i=1,int(cmult)
call acecas(i,1,q,ia,ka,id,kd)
200 if(kdb.ne.0)return
multiplicity
0 indicates Angular Law 1 binned data at xss(ia-1+lmu)
0 indicates a binned distribution sampled
Isotropic Angular Distribution
if(lm.ne.0)go to 20
10 call angiso(t)
Law 1 Equiprobable Binned Angular Distribution
20 if(lm.lt.0)go to 30
it=ia-1+lm
call anglw1(it,t)
Law 2 Tabular Probability Angular Distribution
30 it=ia-1-lm
call anglw2(it,t)
40 acecos=t
*deck acecp acecp
subroutine acecp
generate and bank charged particles from a neutron collision
(very similar to ACEGAM for neutron-induced photons).
character ht*10
dimension yd(2)
data f/1.d0/
uold(1)=uuu
uold(2)=vvv
uold(3)=www
do loop over number of other secondary particles
do 10 i=1,nxs(lnxs+7,iex)
ip=nint(xss(jxs(ljxs+30,iex)+i-1))
if(fiml(ip).eq.0.)go to 10
determine total production cross section for this particle (t)
is=ixs(lixs+10*(i-1)+1,iex)+1
ic=ktc(kktc+1,iex)+1-nint(xss(is-1))
if(ic.lt.1.or.ic.gt.nint(xss(is)).or.ic.eq.nint(xss(is)).and.
1 rtc(krtc+1,iex).ne.0.) go to 10
t=xss(ic+is)
if(rtc(krtc+1,iex).ne.0.) t=t+rtc(krtc+1,iex)*
1 (xss(ic+is+1)-xss(ic+is))
if(t.eq.0.) go to 10
r is the ratio of total particle production cross section
to total neutron cross section and represents the average
number of particles produced per collision.
r=t/rtc(krtc+5,iex)
produce 'n' of this type of secondary particles.
from nearest two integers to r.
weight of each particle is
simply the neutron weight.
n=int(r+rang())
if(n.eq.0) go to 10
if(npb.gt.mpb)then
call expirx(1,'acecp','pbl stack overflow.')
do 20 j=1,npblcm
20 gpb9cm(npb,j)=gpblcm(j)
do 30 j=1,lpblcm
30 jpb9cm(npb,j)=jpblcm(j)
do 40 j=1,n
sample neutron reaction that produced the particle
if(nint(xss(jxs(ljxs+31,iex)+i-1)).eq.1) go to 70
b=rang()*t
do 80 ixre=1,nint(xss(jxs(ljxs+31,iex)+i-1))
l=ixs(lixs+10*(i-1)+5,iex)+
1 nint(xss(ixs(lixs+10*(i-1)+4,iex)+ixre-1))+1
if(nint(xss(l-2)).ne.12) call expirx(1,'acecp',
1 'mftype=12 expected but not found')
mf = 12; partial cross section = yield * neutron cross section
do 90 jy=1,2
This call to ACEFCN incorporates the same f "safety factor" as
is used in ACEGAM. This is probably there to avoid problems at
discontinuities in yield tables. It is unknown whether it is
strictly needed or not.
90 yd(jy)=acefcn(l,xss(jxs(ljxs+1,iex)+ktc(kktc+1,iex)+jy-2)*f,ln)
if(yd(1).eq.0..and.yd(2).eq.0.) go to 80
ix=nint(xss(l-1))
is=jxs(ljxs+7,iex)+nint(xss(jxs(ljxs+6,iex)+ix-1))
RCL writes:
For the logic of the next 4 lines of the code, I originally
copied 4 lines from ACEGAM (ag4xq.9 through ag.133).
looking at the logic contained therein, I don't agree with
the treatment if we are looking up a cross section at an
energy beyond the last tabulated value for the particular
ACEGAM, in such cases, always uses the last value
in the table.
I think it is better to use the logic found
in ACECOL (ac.76 through ac.78).
That is what I have tried
to implement here.
ic=ktc(kktc+1,iex)-nint(xss(is-1))+1
if(ic.gt.0.and.ic.le.nint(xss(is))) x1=xss(is+ic)
if(ic.gt.0.and.ic.le.nint(xss(is))) x2=xss(is+ic)
s=x1*yd(1)+rtc(krtc+1,iex)*(x2*yd(2)-x1*yd(1))
if(s.eq.0.) go to 80
if(s.gt.tt) kx=ixre
tt=max(s,tt)
if(a.gt.b) go to 70
80 continue
handle failure of reaction cross sections to sum
to the total.
treatment is same as in acegam,
whereby a warning is issued (1st time) and the reaction
having the largest cross section is used.
call zaid(2,ht,ixl(lixl+1,iex))
call errprn(1,npxm,3,zero+ipt,erg,'ipt','erg',
1 'no particle-production mt found in acecp. zaid = '//ht)
70 continue
call acecas to sample secondary energy
ntyn=nint(xss(ixs(lixs+10*(i-1)+3,iex)+ixre-1))
mtp=nint(xss(ixs(lixs+10*(i-1)+2,iex)+ixre-1))
ia=ixs(lixs+10*(i-1)+7,iex)
ka=nint(xss(ixs(lixs+10*(i-1)+6,iex)+ixre-1))
id=ixs(lixs+10*(i-1)+9,iex)
kd=nint(xss(ixs(lixs+10*(i-1)+8,iex)+ixre-1))
call acecas(1,1,zero,ia,ka,id,kd)
call rotas(colout(2,1),uold,uuu,lev,irt)
erg=colout(1,1)
vel=slite*sqrt(erg*(erg+2.*gpt(ipt)))/(erg+gpt(ipt))
load particle in bank.
paxtc(1,11,ipt)=paxtc(1,11,ipt)+1.
paxtc(2,11,ipt)=paxtc(2,11,ipt)+wgt
paxtc(3,11,ipt)=paxtc(3,11,ipt)+wgt*erg
pwb(kpwb+ipt,18,icl)=pwb(kpwb+ipt,18,icl)+wgt
call bankit(30)
40 continue
do 50 j=1,lpblcm
50 jpblcm(j)=jpb9cm(npb,j)
do 60 j=1,npblcm
60 gpblcm(j)=gpb9cm(npb,j)
10 continue
c deck ae4a
function acecs6(ii,id,iw,jc,r,cs)
sample energy acecs6 given angle cs, law 67 (endf/b-vi law 7).
ii=0 acecs6 called from acecas (transport): save iw, jc, r
and random numbers generate
ii=1 acecs6 called from calcps (next-event estimator):
use iw, jc, r, and random numbers saved from transport.
dimension el(2),eh(2),rnb(5),lx(2)
save id0,iw0,ic0,rr0,rnb
if(ii.ne.0)go to 10
go to appropriate table for incident energy
10 jw=iw0+2*nint(xss(iw0))+1
ne=nint(xss(jw))
lx(1)=id0+nint(xss(jw+ne+ic0))-1
if(rr0.eq.0.)go to 20
lx(2)=id0+nint(xss(jw+ne+ic0+1))-1
if(ii.eq.0)rnb(1)=rang()
if(rnb(1).le.rr0)ix=2
20 do 110 m=1,2
if(lx(m).eq.0)go to 110
mu=nint(xss(le))
nm=nint(xss(le+1))
find appropriate table for sampled cosine
do 30 iq=1,nm-1
30 if(xss(le+iq+2).ge.cs)go to 40
40 if(mu.eq.1)go to 50
if(ii.eq.0)rnb(ir)=rang()
if(rnb(ir).le.(cs-xss(le+iq+1))/(xss(le+iq+2)-
1 xss(le+iq+1)))iq=iq+1
sample from tabulated energy distribution
50 lb=id0+nint(xss(le+nm+iq+1))
jj=nint(xss(lb-1))
np=nint(xss(lb))
el(m)=xss(lb+1)
eh(m)=xss(lb+np)
if(m.ne.ix)go to 110
if(ii.eq.0)rnb(ir)=rang()
ic=lb+2*np+1
ib=lb+3*np
call bnsrch(rnb(ir),ic,ib,ig)
80 l=ic-2*np
fa=xss(l+np)
if(jj.eq.1)go to 90
bb=(xss(l+np+1)-fa)/(xss(l+1)-ea)
if(bb.eq.0.)go to 90
t=ea+(sqrt(max(zero,fa**2+2.*bb*(rnb(ir)-xss(ic))))-fa)/bb
90 t=ea+(rnb(ir)-xss(ic))/fa
100 acecs6=t
110 continue
use scaled interpolation between energies
if(rr0.eq.0.)return
t1=el(1)+rr0*(el(2)-el(1))
t2=eh(1)+rr0*(eh(2)-eh(1))
acecs6=t1+(acecs6-el(ix))*(t2-t1)/(eh(ix)-el(ix))
function acefcn(it,eg,ln)
interpolate value in table at xss(it) for energy eg
Preconditions:
xss(it) is the first word of an appropriate table
table data are in the format:
nr - number of regions
nbt(i=1..nr) - number of points in int. region i
int(i=1..nr) - interpolation scheme for region i
(nbt and int don't exist if nr is zero and a
linear-linear int. scheme is used across all points)
ne - number of energy points
e(j=1..ne) - energy values
f(j=1..ne) - function values on which to interpolate
the interpolation values are picked using energy value eg
it, eg and all xss(i) are read only variables
Postconditions:
returns value f(eg) appropriately interpolated
return value ln is the length of the table,
i.e. xss(it+ln-1) = f(ne)
if eg is not within the bounds of the table e(i)
i.e. bnsrch returns with a non-zero value of ig
the extreme edge boundary value is returned
acefcn = f(1)
ln and acefcn are modified return values
get key parameter values for table
nr=nint(xss(it))
ie=it+2*nr+1
nf=nint(xss(ie))
ln=2*(nr+nf+1)
binary search for the location of eg in the table.
call bnsrch(eg,il,ih,ig)
set up the parameters for interpolation
ea=xss(il)
eb=xss(ih)
fa=xss(il+nf)
fb=xss(ih+nf)
if outside bounds of table, use extreme edge value
if(ig.ne.0)go to 30
find out which kind of interpolation should be used.
if(nr.eq.0)go to 40
do 10 n=1,nr
10 if(ih-ie.le.nint(xss(it+n)))go to 20
interpolate between table entries.
20 go to(30,40,50,60,70)nint(xss(it+nr+n))
histogram interpolation
30 acefcn=fa
linear-linear interpolation
40 acefcn=fa+(fb-fa)*(eg-ea)/(eb-ea)
log-linear interpolation
50 acefcn=fa+(fb-fa)*log(eg/ea)/log(eb/ea)
linear-log interpolation
60 acefcn=fa*(fb/fa)**((eg-ea)/(eb-ea))
log-log interpolation
70 acefcn=fa*(fb/fa)**(log(eg/ea)/log(eb/ea))
*deck ap acefpt
subroutine acefpt
sum the type 19 fission cross sections, put the totals after
the rest of the cross sections, and recalculate the locators.
this table is required if the problem has any type 7
has any cells where fission i is a
mode n p problem with some expanded photon-
has any fm cards with pseudo-reactions -6, -7, or -8; or is a
kcode problem.
find min and max energy indices over all fission reactions.
do 10 j=1,nxs(lnxs+5,iex)
if(nint(xss(jxs(ljxs+5,iex)+j-1)).ne.19)go to 10
is=jxs(ljxs+7,iex)+nint(xss(jxs(ljxs+6,iex)+j-1))
mi=min(mi,nint(xss(is-1)))
ma=max(ma,nint(xss(is-1))+nint(xss(is))-1)
10 continue
if(ma.eq.0)return
form fsig list beginning at location fis=end+1
kf=jxs(ljxs+22,iex)+2
xss(kf-1)=mi
xss(kf)=ma-mi+1
mc=kf+ma-mi+5
*if def,XS64
mc=mc*ndp2
if(lfll.lt.mc)call chgmem(mdas,lfll,mc,'acefpt a')
zero the fission cross section storage for accumulation.
do 20 i=1,ma-mi+1
20 xss(kf+i)=0.
accumulate fission cross sections for all fission reactions in
corresponding locations of fsig list.
do 40 j=1,nxs(lnxs+5,iex)
if(nint(xss(jxs(ljxs+5,iex)+j-1)).ne.19)go to 40
is=jxs(ljxs+7,iex)+nint(xss(jxs(ljxs+6,iex)+j-1))
ir=kf+nint(xss(is-1))-mi
do 30 i=1,nint(xss(is))
30 xss(ir+i)=xss(ir+i)+xss(is+i)
40 continue
assign new values to fis, end, and table length.
jxs(ljxs+21,iex)=jxs(ljxs+22,iex)+1
jxs(ljxs+22,iex)=jxs(ljxs+22,iex)+ma-mi+3
nxs(lnxs+1,iex)=nxs(lnxs+1,iex)+ma-mi+3
*deck ag acegam
subroutine acegam
generate and bank photons from a neutron collision.
character ht*10
dimension en(30),yd(2)
data f/1.d0/
data en / 1.39e-10, 1.52e-7, 4.14e-7, 1.13e-6, 3.06e-6,
1 8.32e-6, 2.26e-5, 6.14e-5, 1.67e-4, 4.54e-4, 1.235e-3, 3.35e-3,
2 9.12e-3, 2.48e-2, 6.76e-2, 0.184, 0.303, 0.500, 0.823, 1.353,
3 1.738, 2.232, 2.865, 3.68, 6.07, 7.79, 10.0, 12.0, 13.5, 15.0/
return if kcode problem is not settled.
initialize.
if(kc8.gt.0)return
with pikmt, select the nuclide and generate exactly one photon.
if(npikmt.eq.0)go to 30
t=rang()*totgp1
mk=mat(lmat+icl)
do 10 m=jmd(ljmd+mk),jmd(ljmd+mk+1)-1
t=t-tgp(ktgp+lme(llme+1,m))*fme(lfme+m)
10 if(t.lt.0.)go to 20
call expirx(1,'acegam','no isotope found')
20 iex=lme(llme+1,m)
gp=wgt*totgp1/totm
kp=ipan(lipa+icl)+m-jmd(ljmd+mk)
nt=nxs(lnxs+15,iex)
calculate the probability of production of a photon.
30 l=jxs(ljxs+12,iex)+ktc(kktc+1,iex)
gp=(xss(l-1)+rtc(krtc+1,iex)*(xss(l)-xss(l-1)))/rtc(krtc+5,iex)
if(gp.eq.0.)return
find the neutron energy group in the en table.
if(jxs(ljxs+13,iex).ne.0)go to 90
40 if(ib-ic.eq.1)go to 60
ih=(ic+ib)/2
if(erg.lt.en(ih))go to 50
60 lg=jxs(ljxs+12,iex)+nxs(lnxs+3,iex)+(ic-1)*20-1
do 70 it=1,20
70 if(xss(lg+it).gt.elc(2))go to 80
80 gp=gp*(1.-(it-1)*.05)
decide how many photons to make.
90 gp=wgt*gp
if(idx.eq.0)go to 100
t1=max(dxw(1,3),-dxw(1,3)*wgt9(1))
100 if(nww(2).ne.0)go to 110
t1=max(gwt(lgwt+icl),-gwt(lgwt+icl)*wgt9(1))*fiml9(1,1)/fiml(1)
110 t1=huge
do 120 i=1,nww(2)
t=wwval(2,icl,i,1)
if(t.le.0.)t1=min(t1,gp)
120 if(t.gt.0.)t1=min(t,t1)
if(gp.lt.t1)go to 140
if(t1.ne.0.)ni=min(10,int(gp/(5.*t1)+1.))
140 if(gp.le.t1*rang())return
generate photons and write them to the bank.
150 npb=npb+1
if(npb.gt.mpb)then
call expirx(1,'acegam','pbl stack overflow.')
do 160 i=1,npblcm
160 gpb9cm(npb,i)=gpblcm(i)
do 170 i=1,lpblcm
170 jpb9cm(npb,i)=jpblcm(i)
do 180 i=1,ndx(2)
180 if((xxx-dxx(2,1,i))**2+(yyy-dxx(2,2,i))**2+(zzz-dxx(2,3,i))**2.lt.
1 dxx(2,5,i))idx=i
do 450 k=1,ni
sample from full photon-production data, if any.
if(jxs(ljxs+13,iex).eq.0)go to 380
sample the reaction that produced the photon.
if(nxs(lnxs+6,iex).eq.1)go to 350
if(nt.ne.0)go to 190
l=jxs(ljxs+12,iex)+ktc(kktc+1,iex)
r=rang()*(xss(l-1)+rtc(krtc+1,iex)*(xss(l)-xss(l-1)))
do 290 ixre=1,nxs(lnxs+6,iex)
if(nt.eq.0)go to 220
do 200 ik=1,nt
200 if(ixre.eq.nint(xss(jxs(ljxs+29,iex)+2*ik-2)))go to 210
210 pik(kpik+ik)=0.
220 l=jxs(ljxs+15,iex)+nint(xss(jxs(ljxs+14,iex)+ixre-1))+1
if(nint(xss(l-2)).eq.13)go to 250
mf=12 -- partial cross section = yield * neutron cross section.
do 230 jy=1,2
230 yd(jy)=acefcn(l,xss(jxs(ljxs+1,iex)+ktc(kktc+1,iex)+jy-2)*f,ln)
if(yd(1).eq.0..and.yd(2).eq.0.)go to 290
ix=nint(xss(l-1))
if(ix.eq.md)go to 240
is=jxs(ljxs+21,iex)+1
if(ix.gt.0)is=jxs(ljxs+7,iex)+nint(xss(jxs(ljxs+6,iex)+ix-1))
ic=min(ktc(kktc+1,iex)-nint(xss(is-1))+1,nint(xss(is)))
if(ic.gt.0)x1=xss(is+ic)
ic=min(ic+1,nint(xss(is)))
if(ic.gt.0)x2=xss(is+ic)
240 t=x1*yd(1)+(x2*yd(2)-x1*yd(1))*rtc(krtc+1,iex)
mf=13 -- partial cross section is direct from table.
250 ic=ktc(kktc+1,iex)-nint(xss(l-1))+1
if(ic.lt.1)go to 290
if(rtc(krtc+1,iex).ne.0.)go to 260
if(ic.gt.nint(xss(l)))go to 290
t=xss(ic+l)
260 if(ic.ge.nint(xss(l)))go to 290
t=xss(ic+l)+rtc(krtc+1,iex)*(xss(ic+l+1)-xss(ic+l))
270 if(t.eq.0.)go to 290
if(nt.eq.0)go to 280
pik(kpik+ik)=t*xss(jxs(ljxs+29,iex)+ik*2-1)
r=r+pik(kpik+ik)
280 if(t.gt.tt)kx=ixre
tt=max(t,tt)
if(s.gt.r)go to 350
290 continue
if(nt.eq.0)go to 320
if(r.le.0.)go to 450
if(nt.eq.1)go to 310
t=rang()*r
do 300 ik=1,nt
t=t-pik(kpik+ik)
300 if(t.lt.0.)go to 310
call expirx(1,'acegam','no photon mt selected.')
310 ixre=nint(xss(jxs(ljxs+29,iex)+ik*2-2))
wgt=wgt*r/(tgp(ktgp+iex)*xss(jxs(ljxs+29,iex)+ik*2-1))
handle failure of the reaction cross sections to add up.
320 call zaid(2,ht,ixl(lixl+1,iex))
call errprn(1,nppm,4,erg,zero,'erg',' ',
1 'no photon-production mt found in acegam. zaid = '//ht)
sample the energy and direction of the photon.
350 erg=es
mtp=nint(xss(jxs(ljxs+13,iex)+ixre-1))
ia=jxs(ljxs+17,iex)
if(jxs(ljxs+16,iex).ne.0)then
ka=nint(xss(jxs(ljxs+16,iex)+ixre-1))
id=jxs(ljxs+19,iex)
kd=nint(xss(jxs(ljxs+18,iex)+ixre-1))
call acecas(1,1,zero,ia,ka,id,kd)
if(kdb.ne.0)return
if(ixcos.ne.0)ipsc=8
if photon energy below cell cutoff, ignore it.
if(colout(1,1).lt.elc(2))go to 450
call rotas(colout(2,1),uold,uuu,lev,irt)
erg=colout(1,1)
if(nint(xss(jxs(ljxs+13,iex)+ixre-1)).lt.18000)go to 390
if(nint(xss(jxs(ljxs+13,iex)+ixre-1)).gt.19999)go to 390
if(lfcl(llfc+icl).eq.2)go to 450
if(nsr.ne.71)go to 390
if(kcy.le.ikz)go to 390
do 360 j=1+mjss,mjss+nilw
360 if(jss(ljss+j).eq.ncl(lncl+icl))go to 370
370 call sufwrt(j,zero)
otherwise sample from equi-probable energies table.
380 if(it.eq.20)erg=xss(lg+20)
if(it.ne.20)erg=xss(lg+it+int(rang()*(21-it)))
call isos(uuu,lev)
russian roulette photons below weight window.
390 if(nww(2).le.1.or.idx.ne.0)go to 420
do 400 i=1,nww(2)
400 if(erg.lt.wwe(lwwe+i+mww(2)))go to 410
410 t=wwval(2,icl,i,1)
if(wgt.ge.t)go to 420
if(wgt.lt.rang()*t)go to 450
store photon-production tally information.
420 paxtc(1,11,2)=paxtc(1,11,2)+1.
paxtc(2,11,2)=paxtc(2,11,2)+wgt
paxtc(3,11,2)=paxtc(3,11,2)+wgt*erg
if(nmaz.ne.0)call rslmaz(1,1)
pwb(kpwb+2,13,icl)=pwb(kpwb+2,13,icl)+wgt
pan(kpan+2,4,kp)=pan(kpan+2,4,kp)+1.
pan(kpan+2,5,kp)=pan(kpan+2,5,kp)+wgt
pan(kpan+2,6,kp)=pan(kpan+2,6,kp)+wgt*erg
pcc(kpcc+1,icl)=pcc(kpcc+1,icl)+1.
pcc(kpcc+2,icl)=pcc(kpcc+2,icl)+wgt
pcc(kpcc+3,icl)=pcc(kpcc+3,icl)+wgt*erg
do 430 i=1,15
430 if(erg.gt.ebl(lebl+i))go to 440
440 febl(kfeb+1,i)=febl(kfeb+1,i)+1.
febl(kfeb+2,i)=febl(kfeb+2,i)+wgt
tally detectors, create dxtran particles, and bank the photon.
if(locph(llph+icl).ne.0)call talph(icl,erg)
if(kdb.ne.0)return
if(ndet(2).ne.0)call tallyd
if(kdb.ne.0)return
if(ndx(2).gt.1.or.ndx(2).eq.1.and.idx.eq.0)call dxtran
if(kdb.ne.0)return
call bankit(8)
450 if(kdb.ne.0)return
do 460 i=1,lpblcm
460 jpblcm(i)=jpb9cm(npb,i)
do 470 i=1,npblcm
470 gpblcm(i)=gpb9cm(npb,i)
if(npikmt.ne.0)iex=io
*deck an acenu
function acenu()
calculate the average number of neutrons from fission.
the data for all real nuclides are such that acenu will always
be greater than 1 and almost always greater than 2.
data (for polynomial function):
lnu,nc,(c(i),i=1,nc)
find out which kind of calculation to do.
l=jxs(ljxs+2,iex)
if(nint(xss(l)).ne.1)go to 20
calculate acenu by evaluating a polynomial in energy.
n=nint(xss(l+1))-1
acenu=xss(l+n+2)
do 10 i=1,n
10 acenu=acenu*erg+xss(l+n+2-i)
calculate acenu by interpolation in a table.
20 acenu=acefcn(l+1,erg,i)
subroutine acetbl(it,il,r,ln)
returns the interpolation parameters and table length
Preconditions:
xss(it) is the first word of an appropriate interpolated
energy region with data in the format:
nr - number of regions
nbt(i=1..nr) - number of points in int. region i
int(i=1..nr) - interpolation scheme for region i
(nbt and int don't exist if nr is zero and a
linear-linear int. scheme is used across all points)
ne - number of energy points
e(j=1..ne) - energy values
the interpolation values are picked using erg
it, erg and all xss(i) are read only variables
Postconditions:
return value il is the lower indice on which to interpolate
return value r is the interpolation factor such that
e(il) + r * ( e(ih) - e(il) )
acetbl ignores interpolation schemes other than
linear-linear or histogram
return value ln is the length of the table,
i.e. xss(it+ln-1) = e(ne)
if erg is not within the bounds of the table e(i)
i.e. bnsrch returns with a non-zero value of ig
the extreme edge boundary is returned in il
and the interpolation factor is returned as zero
ii = 1 and r = 0.0
il, r and ln are modified return values
nr=nint(xss(it))
ie=it+2*nr+1
nf=nint(xss(ie))
ln=2*(nr+1)+nf
binary search for the location of the energy in the table.
call bnsrch(erg,il,ih,ig)
if(ig.ne.0)go to 40
calculate interpolation fraction r
use histogram interpolation if int(i) = 1
use linear-linear interpolation for all else
if(nr.eq.0)go to 30
do 10 n=1,nr
10 if(ih-ie.le.nint(xss(it+n)))go to 20
20 if(nint(xss(it+nr+n)).eq.1)go to 40
30 if(erg-xss(il).lt.1.e-6*(xss(ih)-xss(il)))go to 40
r=(erg-xss(il))/(xss(ih)-xss(il))
40 il=il-ie
*deck at acetot
subroutine acetot(mm,mk)
calculate the neutron cross sections in material mk.
mm=1 for call from wtmult, otherwise zero.
get the temperature ttn of cell icl.
if(mcal.ne.0)go to 310
if(mm.ne.0)go to 50
l=ltmp+(icl-1)*mxt
if(tme.ge.tth(ltth+mxt))go to 40
if(tme.le.tth(ltth+1))go to 30
do 10 n=2,mxt
10 if(tme.lt.tth(ltth+n))go to 20
20 ttn=tmp(l+n-1)+(tmp(l+n)-tmp(l+n-1))*
1 (tme-tth(ltth+n-1))/(tth(ltth+n)-tth(ltth+n-1))
30 ttn=tmp(l+1)
40 ttn=tmp(l+mxt)
calculate cross sections for each nuclide in the material.
50 totm=0.
do 270 m=jmd(ljmd+mk),jmd(ljmd+mk+1)-1
iex=lme(llme+1,m)
find the energy bin in the cross-section table.
if(rtc(krtc+6,iex).eq.erg)go to 130
n=nxs(lnxs+3,iex)
if(nxs(lnxs+16,iex).eq.lc)go to 100
ic=jxs(ljxs+1,iex)
if(erg.lt.rtc(krtc+6,iex))ib=ic+ktc(kktc+1,iex)
if(erg.gt.rtc(krtc+6,iex))ic=ic+ktc(kktc+1,iex)-1
60 if(ib-ic.eq.1)go to 80
ih=(ic+ib)/2
if(erg.lt.xss(ih))go to 70
get the bin index, the interpolation fraction, and
the absorption and total cross sections.
80 ktc(kktc+1,iex)=ic-jxs(ljxs+1,iex)+1
if(nty(lnty+iex).eq.2)go to 90
rtc(krtc+1,iex)=max(zero,min(one,(erg-xss(ic))/(xss(ib)-xss(ic))))
rtc(krtc+3,iex)=xss(ic+2*n)+rtc(krtc+1,iex)*
1 (xss(ic+2*n+1)-xss(ic+2*n))
rtc(krtc+4,iex)=xss(ic+n)+rtc(krtc+1,iex)*(xss(ib+n)-xss(ic+n))
90 lc=nxs(lnxs+16,iex)
id=ktc(kktc+1,iex)
100 ktc(kktc+1,iex)=id
ic=id+jxs(ljxs+1,iex)-1
110 rtc(krtc+3,iex)=xss(ic+2*n)
rtc(krtc+4,iex)=xss(ic+n)
120 rtc(krtc+6,iex)=erg
rtc(krtc+7,iex)=-1.
130 rtc(krtc+8,iex)=0.
adjust the cross sections for thermal effects if necessary.
if(erg.lt.esa(lesa+lmt(llmt+m)))go to 160
if(rtc(krtc+7,iex).eq.ttn)go to 260
rtc(krtc+7,iex)=ttn
rtc(krtc+5,iex)=rtc(krtc+4,iex)
if(tbt(ltbt+iex).ne.0.)go to 260
a2=awn(lawn+iex)*erg
if(a2.gt.500.*ttn)go to 260
if(a2.lt.4.*ttn)go to 140
f=.5*ttn/a2
140 a=sqrt(a2/ttn)
f=(thgf(i)+(b-i)*(thgf(i+1)-thgf(i)))/a-1.
150 i=ktc(kktc+1,iex)+jxs(ljxs+1,iex)-1+3*nxs(lnxs+3,iex)
rtc(krtc+5,iex)=rtc(krtc+5,iex)+f*(xss(i)+rtc(krtc+1,iex)*
1 (xss(i+1)-xss(i)))
calculate the s(a,b) inelastic scattering cross section.
160 it=lmt(llmt+m)
if(rtc(krtc+6,it).eq.erg)go to 250
n=nint(xss(jxs(ljxs+1,it)))
ic=jxs(ljxs+1,it)+1
ib=jxs(ljxs+1,it)+n
if(erg.lt.rtc(krtc+6,it))ib=ic+ktc(kktc+1,it)
if(erg.gt.rtc(krtc+6,it))ic=ic+ktc(kktc+1,it)-1
rtc(krtc+6,it)=erg
170 if(ib-ic.eq.1)go to 190
ih=(ib+ic)/2
if(erg.lt.xss(ih))go to 180
190 ktc(kktc+1,it)=ic-jxs(ljxs+1,it)
rtc(krtc+1,it)=(erg-xss(ic))/(xss(ib)-xss(ic))
rtc(krtc+8,it)=xss(ic+n)+rtc(krtc+1,it)*(xss(ic+n+1)-xss(ic+n))
rtc(krtc+7,it)=rtc(krtc+8,it)
calculate the s(a,b) elastic scattering cross section.
if(jxs(ljxs+4,it).eq.0)go to 250
n=nint(xss(jxs(ljxs+4,it)))
ic=jxs(ljxs+4,it)+1
ib=jxs(ljxs+4,it)+n
if(nxs(lnxs+5,it).ne.4)go to 200
if(erg.le.xss(ic))go to 250
if(erg.ge.xss(ib))go to 220
200 if(ib-ic.eq.1)go to 230
ih=(ib+ic)/2
if(erg.lt.xss(ih))go to 210
230 ktc(kktc+2,it)=ic-jxs(ljxs+4,it)
if(nxs(lnxs+5,it).ne.4)go to 240
rtc(krtc+7,it)=rtc(krtc+7,it)+xss(ic+n)/erg
240 rtc(krtc+4,it)=(erg-xss(ic))/(xss(ib)-xss(ic))
rtc(krtc+7,it)=rtc(krtc+7,it)+xss(ic+n)+rtc(krtc+4,it)*
1 (xss(ib+n)-xss(ic+n))
250 rtc(krtc+5,iex)=rtc(krtc+7,it)+rtc(krtc+3,iex)
accumulate the total cross section, totm, and the absorption
cross section, siga.
260 totm=totm+rtc(krtc+5,iex)*fme(lfme+m)
270 siga=siga+rtc(krtc+3,iex)*fme(lfme+m)
if(mm.ne.0)return
set up for photon production biasing.
if(npikmt.eq.0)go to 290
if(gwt(lgwt+icl).eq.-1e6)go to 290
do 280 m=jmd(ljmd+mk),jmd(ljmd+mk+1)-1
iex=lme(llme+1,m)
tgp(ktgp+iex)=0.
if(nxs(lnxs+15,iex).lt.0)go to 280
j=jxs(ljxs+12,iex)+ktc(kktc+1,iex)
qu=xss(j-1)
if(nty(lnty+iex).eq.1)qu=qu+rtc(krtc+1,iex)*(xss(j)-xss(j-1))
tgp(ktgp+iex)=qu
totgp1=totgp1+tgp(ktgp+iex)*fme(lfme+m)
280 continue
calculate the fission cross section and nubar.
290 if(lfcl(llfc+icl).eq.0)return
do 300 m=jmd(ljmd+mk),jmd(ljmd+mk+1)-1
iex=lme(llme+1,m)
l=jxs(ljxs+21,iex)
if(l.eq.0)go to 300
if(ktc(kktc+1,iex).lt.nint(xss(l)))go to 300
j=l+2+ktc(kktc+1,iex)-nint(xss(l))
rtc(krtc+8,iex)=xss(j)+rtc(krtc+1,iex)*(xss(j+1)-xss(j))
rtc(krtc+10,iex)=acenu()
300 continue
compute multigroup cross sections.
310 totm=0.
do 330 m=jmd(ljmd+mk),jmd(ljmd+mk+1)-1
iex=lme(llme+1,m)
rtc(krtc+3,iex)=xss(jxs(ljxs+6,iex)-1+jgp)
rtc(krtc+5,iex)=xss(jxs(ljxs+2,iex)-1+jgp)
rtc(krtc+8,iex)=0.
if(nxs(lnxs+11,iex).eq.0)go to 320
if(jxs(ljxs+7,iex).ne.0)stp=stp+fme(lfme+m)*
1 xss(jxs(ljxs+7,iex)-1+jgp)
if(jxs(ljxs+8,iex).ne.0)pfp=pfp+fme(lfme+m)*
1 xss(jxs(ljxs+8,iex)-1+jgp)
320 if(lfcl(llfc+icl).eq.0)go to 330
if(jxs(ljxs+3,iex).eq.0)go to 330
rtc(krtc+8,iex)=xss(jxs(ljxs+3,iex)-1+jgp)
rtc(krtc+10,iex)=xss(jxs(ljxs+4,iex)-1+jgp)
330 totm=totm+rtc(krtc+5,iex)*fme(lfme+m)
*deck az action
subroutine action
print cell activity, weight balance, and nuclide activity.
TEMPORARY RESTRICTION to n, p, e PWB tables.
dimension kl(4),l3(3)
character hc(2)*10,ht*6,hz*10
data kl/1,6,13,20/,l3/17,19,18/
data hc/'collisions',' substeps'/
print the table of activity in each cell of the problem.
do 70 ip=1,mipt
if(kpt(ip).eq.0)go to 70
if(charge(ip).eq.zero)write(iuo,10)hnp(ip),hc(1),hc(1)
if(charge(ip).ne.zero)write(iuo,10)hnp(ip),hc(2),hc(2)
10 format(1h1,a11,22h activity in each cell,70x,15hprint table 126//
1 17x,6htracks,5x,10hpopulation,3x,a10,3x,
2 a10,5x,6hnumber,8x,4hflux,8x,7haverage,6x,7haverage/
3 8x,4hcell,4x,8hentering,31x,8h* weight,2(5x,8hweighted),
4 3x,5htrack,19h weight
track mfp/52x,13h(per history),4x,
5 6henergy,7x,6henergy,5x,10h(relative),6x,4h(cm)/)
do 20 j=1,8
20 tpp(j)=0.
do 50 i=1,mxa
if(fim(lfim+ip,i).eq.0.)go to 50
do 30 j=1,4
tpp(j)=tpp(j)+pac(lpac+ip,j,i)
30 tpp(j+4)=0.
if(pac(lpac+ip,9,i).ne.0.)tpp(5)=pac(lpac+ip,5,i)/pac(lpac+ip,9,i)
if(pac(lpac+ip,10,i).ne.0.)tpp(6)=pac(lpac+ip,6,i)/
1 pac(lpac+ip,10,i)
if(pac(lpac+ip,7,i).ne.0.)tpp(7)=pac(lpac+ip,10,i)*fim(lfim+ip,i)
1 /pac(lpac+ip,7,i)
if(pac(lpac+ip,10,i).ne.0.)tpp(8)=pac(lpac+ip,8,i)/
1 pac(lpac+ip,10,i)
write(iuo,40)i,ncl(lncl+i),(pac(lpac+ip,j,i),j=1,3),
1 pac(lpac+ip,4,i)*fpi,(tpp(j),j=5,8)
40 format(2i6,f13.0,tl1,f14.0,tl1,f14.0,tl1,1h ,1p5e13.4)
50 continue
write(iuo,60)tpp(1),tpp(2),tpp(3),tpp(4)*fpi
60 format(/5x,5htotal,f15.0,tl1,f14.0,tl1,f14.0,tl1,1pe14.4)
if(ink(128).ne.0)call mapmaz(ip)
70 continue
print the weight balance tables.
if(ink(130).eq.0)go to 200
TEMPORARY RESTRICTION to n, p, e PWB tables.
do 190 ip=1,3
if(kpt(ip).eq.0)go to 190
do 80 i=1,kl(4)-1
80 tpp(i)=0.
do 180 jw=1,3
if(jw.eq.1)write(iuo,90)hnp(ip)
90 format(1h1,a11,40h weight balance in each cell -- external,
1 7h events,45x,15hprint table 130//8x,4hcell,4x,8hentering,6x,
2 6hsource,7x,6henergy,8x,4htime,7x,7hexiting,7x,5htotal/
3 43x,6hcutoff,7x,6hcutoff/)
if(jw.eq.2)write(iuo,100)hnp(ip)
100 format(1h1,a11,40h weight balance in each cell -- variance,
1 17h reduction events,35x,15hprint table 130//8x,4hcell,5x,
2 6hweight,8x,4hcell,8x,6hweight,7x,6henergy,7x,6hdxtran,7x,
3 6hforced,4x,11hexponential,5x,5htotal/17x,6hwindow,5x,
4 10himportance,5x,6hcutoff,5x,10himportance,16x,9hcollision,4x,
5 9htransform/)
if(jw.eq.3.and.ip.eq.1)write(iuo,110)
110 format(1h1,47hneutron weight balance in each cell -- physical,
1 7h events,49x,15hprint table 130//8x,4hcell,5x,6h(n,xn),6x,
2 7hfission,6x,7hcapture,6x,7hloss to,6x,7hloss to,7x,5htotal/
3 56x,6h(n,xn),6x,7hfission/)
if(jw.eq.3.and.ip.eq.2)write(iuo,120)
120 format(1h1,46hphoton weight balance in each cell -- physical,
1 7h events,50x,15hprint table 130//8x,4hcell,6x,4hfrom,8x,
2 6hbrems-,5x,9hp-annihi-,5x,8helectron,2x,12hfluorescence,4x,
3 7hcapture,8x,4hpair,8x,5htotal/16x,8hneutrons,4x,9hstrahlung,6x,
4 6hlation,7x,6hx-rays,31x,10hproduction/)
if(jw.eq.3.and.ip.eq.3)write(iuo,130)
130 format(1h1,48helectron weight balance in each cell -- physical,
1 7h events,48x,15hprint table 130//8x,4hcell,6x,4hpair,7x,
2 7hcompton,7x,6hphoto-,7x,6hphoton,6x,8helectron,5x,8hknock-on,
3 6x,5htotal/15x,10hproduction,5x,6hrecoil,6x,8helectron,6x,
4 5hauger,8x,5hauger/)
j2=kl(jw+1)-1
if(jw.eq.3)j2=l3(ip)
do 160 i=1,mxa
if(fim(lfim+ip,i).eq.0.)go to 160
do 140 j=j1,j2
tpp(j)=tpp(j)+fpi*pwb(lpwb+ip,j,i)
t3=max(t3,fpi*pwb(lpwb+ip,j,i))
extra parenthesis required for convex compiler optimization.
140 t1=(t1+fpi*pwb(lpwb+ip,j,i))
if(abs(t1).lt.1.e-10*t3)t1=0.
t4=max(t4,t1)
write(iuo,150)i,ncl(lncl+i),(fpi*pwb(lpwb+ip,j,i),j=j1,j2),t1
150 format(2i6,1p8e13.}