*cpl module *lst module *ident up1 */ heatr -- 21 mar 91 -- allow for more mf6 reactions. */ also, isotropy message should not be */ printed if mf6 is present. *d heatr.336 common/hmfsix/i6,mt6(320),i6g,mt6g(320) *d heatr.400,401 if (mf4.eq.4.and.mt4.eq.mtd) go to 140 if (mf4.eq.6.and.mt4.eq.mtd) go to 140 *d heatr.403 write(strng, 1 '(''mf4 and 6 missing, isotropy assumed for mt '',i3)') mtd *d heatr.417,418 if (mf4.eq.4.and.mt4.eq.mtd) go to 140 if (mf4.eq.6.and.mt4.eq.mtd) go to 140 *d heatr.420 write(strng, 1 '(''mf4 and 6 missing, isotropy assumed for mt '',i3)') mtd *i heatr.429 if (i6.gt.320) call error('hinit','too many mf6 reactions',' ') *d heatr.570 common/hmfsix/i6,mt6(320),i6g,mt6g(320) *d heatr.3131 common/hmfsix/i6,mt6(320),i6g,mt6g(320) */ heatr -- 21 mar 91 -- fix mgam flag. */ 10 means mf6 gammas, but no mf12 or 13 gammas. *d heatr.462 if (zap.eq.0) mgam=10+mod(mgam,10) *ident up2 */ reconr -- 26 mar 91 -- simplify mlbw and make sure that */ duplicate j values are treated correctly. *d reconr.494 *d reconr.705,706 jnow=jnow+6 *d reconr.741,742 if (mode.le.3) go to 370 *d reconr.758,795 *i reconr.2174 den=4.*spi+2. *d reconr.2189,2190 *d reconr.2200,2202 c reconstruct statistical factors sum=0. fl=float(ll) ajmin=abs(abs(spi-fl)-0.5) ajmax=spi+fl+0.5 nj=ajmax-ajmin+1.001 aj=ajmin do 52 i=1,nj gj(i)=(2.*aj+1.)/den aj=aj+1. 52 sum=sum+gj(i) diff=2*fl+1.-sum *d reconr.2211 j=a(inow+1)-ajmin+1.001 *d reconr.2245 sigp(2)=sigp(2)+2.*diff*cos2p inow=in */ reconr -- 26 mar 91 -- install dunford's patch to treat */ duplicate j values for reich-moore. *i reconr.2338 if (fl.gt.spi-0.5) then jjl=0 else jjl=1 endif *d reconr.2431 u11i=p2*(uno-two*ri(1,1))+two*p1*si(1,1) *i reconr.2433 if (jj.gt.jjl.and.jj.lt.numj) then termn=termn+two*gj*(1-p1) termt=termt+two*gj*(1-p1) endif *ident up3 */ plotr -- 4 apr 91 -- use more points when plotting */ 3-d energy distributions. *d plotr.931 2525 ne2=200 */ plotr -- 4 apr 91 -- last point in 3d ang plot missing *d plotr.824,825 if (en.lt..999*xmin) go to 2460 if (en.gt.1.001*xmax) go to 2460 */ plotr -- 4 apr 91 -- add 3d ang dist truncation message *i plotr.870 if (locn+3*nmu.gt.maxaa) call mess('plotr', 1 'too much 3d angular distribution data', 2 'energy range truncated') */ plotr -- 4 apr 91 -- protect storage for 3d energy plots *i plotr.306 data nw7max/6000/ *d plotr.932,933 locn=0 if (law.eq.7) locn=locn+nw7max *d plotr.938 if (mfd.eq.6) 1 call gety6(1,0.,en,idis,f2,nin,a,2,law,aa,nw7max) *d plotr.940,941 if (ei.lt..999*xmin) go to 2545 if (ei.gt.1.001*xmax) go to 2545 if (locn+3*ne2.gt.maxaa) go to 2545 *d plotr.956 if (mfd.eq.6) 1 call gety6(1,e2,en,idis,f2,nin,a,2,law,aa,nw7max) *d plotr.972 if (mfd.eq.6) 1 call gety6(1,1.e8,en,idis,f2,nin,a,2,law,aa,nw7max) *i plotr.973 if (locn+3*ne2.gt.maxaa) call mess('plotr', 1 'too much energy distribution data', 2 'energy range truncated') *d plotr.1883 subroutine gety6(i6,x,xnext,idis,y6,itape,a,lep,law,aa,maa) *d plotr.1901 if (law.eq.7) call fixl7(itape,aa,maa,a,nb,nw) *d plotr.1967 subroutine fixl7(itape,aa,maa,a,nb,nw) *i plotr.1991 if (loc.gt.maa) call error('fixl7', 1 'not enough storage to convert file 7',' ') */ plotr -- 4 apr 91 -- echo input to output */ to help in finding input errors *i plotr.312 call timer(time) write(nsyso,10) time if (ntty.gt.0) write(ntty,10) time *i plotr.328 write(nsyso,15) xpage,ypage,istyle,size *i plotr.350 write(nsyso,20) iplot,factx,facty,xll,yll,xur,yur *i plotr.369 write(nsyso,25) t1,t2 *i plotr.387 write(nsyso,30) itype,jtype,igrid,ileg if (ileg.eq.1) write(nsyso,31) xtag,ytag *i plotr.423 write(nsyso,35) xleft,xright,xstep,xl *i plotr.458 write(nsyso,36) ybot,ytop,ystep,yl *i plotr.485 write(nsyso,37) rbot,rtop,rstep,rl *i plotr.512 write(nsyso,40) iverf,nin,matd,mfd,mtd,temper,nth,ntp,nkh *i plotr.527 write(nsyso,45) icon,isym,idash,ithick *i plotr.540 if (ileg.eq.1) write(nsyso,50) aleg *i plotr.551 write(nsyso,51) aleg,xtag,ytag,xpoint *i plotr.570 write(nsyso,55) xv,yv,zv,x3,y3,z3 *i plotr.1111 write(nsyso,60) nform *i plotr.1161 10 format(/' plotr...plot endf, pendf, gendf, or input data', 1 22x,f8.1,1hs) 15 format(/40h xpage ................................ ,f10.3/ 1 40h ypage ................................ ,f10.3/ 2 40h istyle ............................... ,i10/ 3 40h size ................................. ,f10.3) 20 format(/50h -------------------------------------------------// 1 40h iplot ................................ ,i10/ 2 40h factx ................................ ,1p,e10.2/ 2 40h facty ................................ ,1p,e10.2/ 3 40h xll .................................. ,f10.3/ 4 40h yll .................................. ,f10.3/ 5 40h xur .................................. ,f10.3/ 6 40h yur .................................. ,f10.3) 25 format(/12x,a/12x,a) 30 format(/40h itype ................................ ,i10/ 1 40h jtype ................................ ,i10/ 2 40h igrid ................................ ,i10/ 3 40h ileg ................................. ,i10) 31 format(/40h xtag ................................. ,1p,e10.2/ 1 40h ytag ................................. ,1p,e10.2) 35 format(/40h xlo .................................. ,1p,e10.2/ 1 40h xhi .................................. ,1p,e10.2/ 2 40h xstep ................................ ,1p,e10.2/ 3 12h xlabel ... ,a) 36 format(/40h ylo .................................. ,1p,e10.2/ 1 40h yhi .................................. ,1p,e10.2/ 2 40h ystep ................................ ,1p,e10.2/ 3 12h ylabel ... ,a) 37 format(/40h rlo .................................. ,1p,e10.2/ 1 40h rhi .................................. ,1p,e10.2/ 2 40h rstep ................................ ,1p,e10.2/ 3 12h rlabel ... ,a) 40 format(/40h iverf ................................ ,i10/ 1 40h nin .................................. ,i10/ 2 40h matd ................................. ,i10/ 3 40h mfd .................................. ,i10/ 4 40h mtd .................................. ,i10/ 5 40h temp ................................. ,1p,e10.2/ 6 40h nth .................................. ,i10/ 7 40h ntp .................................. ,i10/ 8 40h nkh .................................. ,i10) 45 format(/40h icon ................................. ,i10/ 1 40h isym ................................. ,i10/ 2 40h idash ................................ ,i10/ 3 40h ithick ............................... ,i10) 50 format(/12h legend ... ,a) 51 format(/12h tag ... ,a/ 1 40h xtag ................................. ,1p,e10.2/ 2 40h ytag ................................. ,1p,e10.2/ 3 40h xpoint ............................... ,1p,e10.2) 55 format(/40h 3d viewpoint ......................... ,f10.3/ 1 40h ,f10.3/ 2 40h ,f10.3/ 3 40h 3d workbox ........................... ,f10.3/ 4 40h ,f10.3/ 5 40h ,f10.3) 60 format(/40h data nform ........................... ,i10) *ident up4 */ heatr -- 8 apr 91 -- fix interpolation range problem in sixbar *d heatr.1982 305 if (e.lt.ehi) go to 400 *d heatr.1984 if (nne.eq.ne.and.e.eq.ehi) go to 400 *d heatr.2035 go to 305 *ident up5 */ mixr -- 16 apr 91 -- fix bad error message *d mixr.133 if (math.lt.0) call error('mixr','mat and temp not found',' ') */ mixr -- 16 apr 91 -- fix interpolation problem *d mixr.404 if (ip.lt.ip2) go to 125 *d mixr.412 ip1=ip-1 *ident up6 */ moder -- 19 apr 91 -- allow file 40 to be processed *d moder.177 260 if (mf.gt.40) go to 490 *ident up7 */ reconr -- 19 apr 91 -- add a grid point for low-energy */ resonance evaluations (protects */ reconstruction stack from overflow) *i reconr.904 if (el.lt..0253) then nodes=nodes+1 a(ienode+nodes-1)=.0253 endif */ reconr -- 18 apr 91 -- fix some problems with */ discontinuities and thresholds *i reconr.1277 if (mth.eq.2) go to 210 if (mth.eq.18.or.mth.eq.19) go to 210 if (mth.eq.102) go to 210 *d reconr.1291 if (er.gt.sigfig(abs(eg),ndigit,-1)) go to 222 *i reconr.1293 er=sigfig(er,ndigit,1) *i reconr.1295 go to 225 222 er=abs(eg) *ident up8 */ matxsr -- 24 apr 91 -- fix length of file id block *d matxsr.745 nwds=6 *i matxsr.751 ia(icont+5)=0 */ matxsr -- 24 apr 91 -- increase container size *d matxsr.397 common/mstore/a(50000) *d matxsr.408 isiza=50000 *d matxsr.538 common/mstore/a(50000) *d matxsr.718 common/mstore/a(50000) *d matxsr.884 common/mstore/a(50000) *d matxsr.1449 common/mstore/a(50000) *d matxsr.1668 common/mstore/a(50000) *d matxsr.1776 common/mstore/a(50000) *d matxsr.1917 common/mstore/a(50000) *d matxsr.2010 common/mstore/a(50000) *d matxsr.2177 common/mstore/a(50000) *d matxsr.2519 common/mstore/a(10000),b(25000) */ matxsr -- 24 apr 91 -- fix indexing in indexm */ change to 80-column format *d matxsr.2539,2541 hn=ha(l1h) hu=ha(l1h+1) hs=ha(l1h+2) *d matxsr.2561 loc=l1h-1 *d matxsr.2703,2704 45 format(//21h --------------------) *d matxsr.2706,2707 50 format(6x,'vectors: ',8a8/(16x,8a8)) 55 format(6x,'matrices: ',8a8/(16x,8a8)) */ matxsr -- 24 apr 91 -- allow for more reactions *d matxsr.1456 common/hollr/hvps(300),hmtx(300) *d matxsr.1673 common/hollr/hvps(300),hmtx(300) *d matxsr.1783 common/hollr/hvps(300),hmtx(300) *d matxsr.2180 dimension ngrp(10),htype(300) *d matxsr.2182 dimension nfg(300),nlg(300),nmg(300) */ matxsr -- 29 apr 91 -- fix paging errors. */ keep last block smaller than maxw. *d matxsr.646,647 if (nwds+nw.lt.maxw) go to 210 *i matxsr.653 if (nwds.eq.0) go to 225 irec3=irec3+1 call reed(nscrt3,irec3,a(ivdat),nwds,0) irec=irec+1 call rite(nmatx,irec,a(ivdat),nwds,0) *d matxsr.683,684 if (nwds+nw*lord.lt.maxw) go to 270 *i matxsr.691 if (nwds.eq.0) go to 275 irec3=irec3+1 call reed(nscrt3,irec3,a(imdat),nwds,0) irec=irec+1 call rite(nmatx,irec,a(imdat),nwds,0) *i matxsr.1722 ivnow=ivdat *d matxsr.1726,1727 if (nwds+nw.lt.maxw) go to 295 *d matxsr.1729 call rite(nscrt3,irec3,a(ivnow),nwds,0) ivnow=ivnow+nwds *i matxsr.1732 if (nwds.eq.0) go to 310 irec3=irec3+1 call rite(nscrt3,irec3,a(ivnow),nwds,0) ivnow=ivnow+nwds nritev=nritev+1 *d matxsr.1735 310 irec2=irec2+1 *i matxsr.2018 c c ***retrieve pointers call findex('rinp',irinp,a) irlim=irinp+lirinp call findex('cdat',icdat,a) call findex('mdat',imdat,a) *d matxsr.2034,2036 if (num+n.lt.maxw) go to 101 imax=i-1 *d matxsr.2046,2049 *d matxsr.2147 imax1=imax+1 do 300 i=i0,imax1 if (i.gt.imax) go to 290 *d matxsr.2151,2153 if (nwds+np.lt.maxw) go to 300 290 if (nwds.eq.0) go to 300 *d matxsr.2305 n1d1=n1d+1 do 155 ik=1,n1d1 if (ik.gt.n1d) go to 128 *d matxsr.2310,2316 if (nwds+nw.lt.maxw) go to 125 128 if (nwds.eq.0) go to 155 *d matxsr.2331 nwblk=nblk-mblk+1 *d matxsr.2333 do 140 jn=1,nwblk *i matxsr.2349 irct=0 nmg(1)=1 *i matxsr.2350 irct=irct+1 nfg(irct)=nfgi nlg(irct)=nlgi nmg(irct+1)=nwds+1 htype(irct)=ha(l3h-1+ik) *d matxsr.2383 noutg1=noutg+1 do 170 i=1,noutg1 if (i.gt.noutg) go to 168 *d matxsr.2386,2388 if (nwds+nw.lt.maxw) go to 170 168 if (nwds.eq.0) go to 170 *d matxsr.2637,2638 if (nwds+nw.lt.maxw) go to 125 *i matxsr.2641 if (nwds.gt.0) irec=irec+1 *d matxsr.2663,2664 if (nwds+nw.lt.maxw) go to 160 *i matxsr.2667 if (nwds.gt.0) irec=irec+1 *ident up9 */ resxsr -- 18 jun 91 -- fix material control and blocking *d resxsr.396,397 ia(iscr+mult+n+1)=jx/n ia(iscr+mult+n+2)=je *d resxsr.405,410 if (nwds+nn.lt.nblok) go to 415 irec=irec+1 *i resxsr.416 if (nwds.eq.0) go to 400 irec=irec+1 call rite(nscr,irec,a(iscr),nwds,0) *d resxsr.468,469 nreac=ia(iscr+mult+ntemp(i)+1) nener=ia(iscr+mult+ntemp(i)+2) *ident up10 */ gaminr -- 11 jul 91 -- nint ngg and fix eq in gpanel *d gaminr.547 ngg=nint(egg(1)) *d gaminr.799 if (eq.gt.ehigh) eq=ehigh *i gaminr.836 if (enext.le.ehi) enext=rndoff*ehi */ gaminr -- 11 jul 91 -- add vitamin-j group structure *i gaminr.63 c * 10 vitamin-j 42-group structure * *i gaminr.502 c 10 vitamin-j 42-group structure *d gaminr.510 dimension eg8(39),eg10(43) *i gaminr.532 data eg10/1.0e3,1.0e4,2.0e4,3.0e4,4.5e4,6.0e4,7.0e4,7.5e4,1.0e5, 1 1.50e5,2.00e5,3.00e5,4.00e5,4.50e5,5.10e5,5.12e5,6.00e5, 2 7.00e5,8.00e5,1.00e6,1.33e6,1.34e6,1.50e6,1.66e6,2.00e6, 3 2.50e6,3.00e6,3.50e6,4.00e6,4.50e6,5.00e6,5.50e6,6.00e6, 4 6.50e6,7.00e6,7.50e6,8.00e6,1.00e7,1.20e7,1.40e7,2.00e7, 5 3.00e7,5.00e7/ *d gaminr.539 go to (100,200,300,400,500,600,700,800,800,850),igg *i gaminr.608 c c ***vitamin-j 42-group structure 850 ngg=42 do 860 ig=1,43 860 egg(ig)=eg10(ig) go to 900 *i gaminr.620 if (igg.eq.10) write(nsyso,95) *i gaminr.634 95 format(/46h gamma group structure......vitamin-j 42-group) *ident up11 */ groupr -- 11 jul 91 -- add vitamin-j photon group structure *i groupr.142 c * 10 vitamin-j 42-group structure * *i groupr.1521 c 10 vitamin-j 42-group structure *d groupr.1529 dimension eg8(39),eg10(43) *i groupr.1551 data eg10/1.0e3,1.0e4,2.0e4,3.0e4,4.5e4,6.0e4,7.0e4,7.5e4,1.0e5, 1 1.50e5,2.00e5,3.00e5,4.00e5,4.50e5,5.10e5,5.12e5,6.00e5, 2 7.00e5,8.00e5,1.00e6,1.33e6,1.34e6,1.50e6,1.66e6,2.00e6, 3 2.50e6,3.00e6,3.50e6,4.00e6,4.50e6,5.00e6,5.50e6,6.00e6, 4 6.50e6,7.00e6,7.50e6,8.00e6,1.00e7,1.20e7,1.40e7,2.00e7, 5 3.00e7,5.00e7/ *d groupr.1558 go to (100,200,300,400,500,600,700,800,800,850),igg *i groupr.1625 c c ***vitamin-j 42-group structure 850 ngg=42 do 860 ig=1,43 860 egg(ig)=eg10(ig) go to 900 *i groupr.1645 if (igg.eq.10) write(nsyso,95) *i groupr.1659 95 format(/46h gamma group structure......vitamin-j 42-group) *ident up12 */ njoy -- 17 jul 91 -- fix jump back from resxsr *i njoy.277 go to 100 *ident up13 */ groupr -- 18 jul 91 -- add more t-dep weight functions *d groupr.152 c * 7 same with t-dep thermal part * *i groupr.156 c * 12 vit-e with t-dep thermal part * *d groupr.1677 c 7 same with t-dep thermal part *i groupr.1681 c 12 vit-e with t-dep thermal part *d groupr.1695 dimension w8(66) *d groupr.1725,1733 *d groupr.1796 go to (110,120,130,140,150,160,160,180,190,190,200,200),iwtt *i groupr.1850 if (iwtt.gt.6) write(nsyso,16) *d groupr.1852,1859 *i groupr.1880 if (iwtt.gt.11) write(nsyso,16) *d groupr.1945 go to (100,200,300,400,100,600,600,100,100,70,700,700),iwtt *d groupr.1997 c ***with optional t dependence 600 tt=.054 if (iwtt.gt.6) tt=temp(jtemp)*8.61735e-5 bb=2*tt cc=1. if (iwtt.gt.6) cc=1.578551e-3*exp(2.)/bb**2 if (e.gt.bb) go to 610 wtf=cc*e*exp(-e/tt) *d groupr.2014 c ***with optional t dependence 700 enext=1.01*e *d groupr.2016 tt=therm if (iwtt.gt.11) tt=temp(jtemp)*8.61735e-5 cc=con1 if (iwtt.gt.11) cc=con2*exp(en1/tt)/en1**2 wtf=cc*e*exp(-e/tt) *ident up14 */ unresr -- 29 sep 91 -- increase to 10 temps and 10 sigzeros *d unresr.25,26 c * ntemp no. of temperatures (10 max) * c * nsigz no. of sigma zeroes (10 max) * *d unresr.45 dimension temp(10),sigz(10),bkgz(4),sigu(5,10) *d unresr.49 data ntmax/10/, nzmax/10/, ipr/1/ *d unresr.727,729 dimension sigf(2,10,5) dimension abns(30),tj(30,10,3),tk(30,10,3),tl(90,10,3), 1 del(3),sigm(10),t(4,10,3),yy(3),gg(5),temp(3) */ update card count in directory (raepsaet, cea-dmt) *d unresr.239 if (mfd.le.2.and.mtd.ne.152) go to 175 *ident up15 */ reconr -- 3 dec 91 -- fix channel radius (rowlands and eaton). */ fix enode name for reserv, etc. (panini). *d reconr.150 130 call reserv('enod',nodmax,ienode,a) *d reconr.220 call releas('enod',1,a) *d reconr.233 call releas('enod',-1,a) *d reconr.275 call findex('enod',ienode,a) *d reconr.498 call findex('enod',ienode,a) *i reconr.687 ascatl=ascat if (lrf.eq.3) ascatl=a(jnow+1) if (ascatl.eq.0.) ascatl=ascat *i reconr.690 ral=ra if (naps.eq.1) ral=ascatl *d reconr.727 210 rho=2.196771e-3*arat*sqrt(abs(a(jnow)))*ral *d reconr.1173 call findex('enod',ienode,a) *ident up16 */ wimsr -- 3 dec 91 -- fix name for reserv, etc. (panini) *d wimsr.668 call reserv('elst',nwelas,ielas,a) *d wimsr.695 call findex('elst',ielas,a) *d wimsr.1016 call findex('elst',ielas,a) *ident up17 */ groupr -- 3 dec 91 -- fix type of id in call to reserv *d groupr.6750 *d groupr.6756 call reserv('sc',nc,ic,c) */ groupr -- 16 dec 91 -- fix allocation for sed (laughton, aecl) *d groupr.3435 if (e.eq.0.) call reserv('sed',nwds,ised,a) */ groupr -- 16 dec 91 -- fix print of mfd=5/mtd=18 (laughton, aecl) *d groupr.2743 210 if (mfd.eq.5) go to 400 *i groupr.2847 c c ***prompt fission spectrum (mfd=5/mtd=18) 400 if (mtd.eq.455) go to 410 if (iprint.ne.1) return do 405 j=1,ng2,10 lim=j+9 if (lim.gt.ng2) lim=ng2 405 write(nsyso,65) j,(ans(1,1,i),i=j,lim) igzero=1 return */ groupr -- 16 dec 91 -- extend lower limit for tabulated secondary */ distributions to 1e-5 and upper limit to */ 50e6 as is done for analytic sections. */ noticed by laughton, aecl. *d groupr.6863 e1=eg(ig) if (ig.eq.1) e1=1.e-5 e2=eg(ig+1) if (ig.eq.ng) e2=50.e6 call intega(c(l),e1,e2,c(iraw+ic),ip,ir) *d groupr.6936 e1=eg(ig) if (ig.eq.1) e1=1.e-5 e2=eg(ig+1) if (ig.eq.ng) e2=50.e6 280 call intega(c(l),e1,e2,c(iraw+ic),ip,ir) *ident up18 */ njoy -- 3 dec 91 -- pad string parameter with spaces for */ better comparisons on vax (gauld, whiteshell) *d njoy.2749 subroutine findex(ident,index,ia) *i njoy.2760 character*(*) ident *i njoy.2764 id=ident(1:len(ident))//' ' *d njoy.2569 subroutine reserv(ident,nwords,index,ia) *i njoy.2584 character*(*) ident *i njoy.2591 id=ident(1:len(ident))//' ' *d njoy.2676 subroutine releas(ident,nwords,ia) *i njoy.2689 character*(*) ident *i njoy.2693 id=ident(1:len(ident))//' ' *ident up19 */ matxsr -- 4 dec 91 -- fix hollerith word size and pointers *i matxsr.409 ncw=8 *d matxsr.412 *d matxsr.415 *d matxsr.568 nwds=3*mult+1 *d matxsr.2225,2226 l2=1+nwds if (2*(l2/2).eq.l2) l2=l2+1 l2h=(l2-1)/mult+1 *d matxsr.2275 write(nout,380) ha(l2h),a(nex1) *i matxsr.2276 if (2*(l3/2).eq.l3) l3=l3+1 *ident up20 */ moder -- 4 dec 91 -- fix gendf processing (dean, winfrith). *d moder.251 730 if ((mf.ne.1).and.(ig.ne.ng.and.ig.ne.0)) go to 710 */ moder -- 4 dec 91 -- fix file 34 processing (muir, iaea). *i moder.1114 call contio(nin,nout,nscr,a,nb,nw) *ident up21 */ heatr -- 4 dec 91 -- prevent div by zero in sed (decher, ce) *d heatr.1626,1629 *i heatr.1688 c function sed(e,ep,theta,u) c ****************************************************************** c secondary energy distribution function for anadam. c a taylor expansion is used for small values to avoid overflow. c ****************************************************************** xeu=(e-u)/theta aeu=abs(xeu) if (aeu.lt.1.0e-2) go to 200 sed=ep*exp(-ep/theta)/(theta*theta*(1.-exp(-xeu)*(1.+xeu))) return 200 sed=ep*exp(-ep/theta)/(theta*theta*0.5*xeu*xeu*(1.-xeu)) return end *ident up22 */ thermr -- 4 dec 91 -- set sabflg for sw machines (dean, aguilar) *d thermr.1559 *if sw data sabflg/-89./ *else data sabflg/-225./ *endif */ thermr -- 4 dec 91 -- fix tolerance on test for x(i)=0 (decher, ce). *d thermr.1798 if (xn.gt.xl.and.xn.lt.(x(i)+xtol*(x(i)-xl))) go to 180 */ fix value of boltzmann's constant *d thermr.1186 data bk/8.61735e-5/ *ident up23 */ broadr -- 4 dec 91 -- don't thin lowest threshold (dean, winfrith). *i broadr.290 thnmax=fact*thnmax *ident up24 */ groupr -- 4 dec 91 -- impr. to flux calculator (winfrith, lanl, cea). *i groupr.1778 z(5)=1. *d groupr.1784 call openz(ninwt,1) *d groupr.1786 jsigz=nint(z(5)) *d groupr.1904 46 format(7x,13halpha2, sam =,1p,e12.4,0p,f8.3/ *d groupr.2120 nfv=1+nl*nsigz+1 */ check out the following, if ok, change to *d and remove the c. *i groupr.2166 c 120 b(iz+3+li)=(sigz(iz)-sam)*wtf*(1.-beta) *d groupr.2198,2199 f3=f1*alog(ejp/ej)-1. f4=1.-f2*alog(ejp/ej) *d groupr.2212,2215 g3=f1*alog(elim/ej)-(elim-ej)/(ejp-ej) g4=(elim-ej)/(ejp-ej)-f2*alog(elim/ej) b(lj+3+iz)=b(lj+3+iz)+g4*s2*b(lj+3+iz+nx)/(1.-alpha(k)) denom(iz)=denom(iz)-g3*s1/(1.-alpha(k)) *d groupr.2249,2252 g3=f1*alog(elim/ej)-(elim-ej)/(ejp-ej) g4=(elim-ej)/(ejp-ej)-f2*alog(elim/ej) b(li+3+iz)=b(li+3+iz)+(g3*s1*b(lj+3+iz) 1 +g4*s2*b(lj+3+iz+nx))/(1.-alpha(k)) *d groupr.2274 l=(iz-1)*nl+1 *d groupr.2288 indx=2+nl*(jsigz-1) *d groupr.2299 l=1+il+nl*(iz-1) *d groupr.2339 l=1+il+nl*(iz-1) *ident up25 */ acer -- 6 dec 91 -- misc acer fixes */ fix messages in etopl routines *d acer.2102 write(nsyso,1) a(iamu+j-1),a(iamu+j),a(iprob+j-1),e,mat,mf,mt *d acer.2127 c described above. if it differs from 1.0 by more than 1.0e-6 *d acer.2130 if (abs(1.-sum).gt.1.e-6) write(nsyso,3) avar,mat,mf,mt *d acer.2189 c now compute area for mu-pmu histogram. if *d acer.2217,2218 1 format(/46h message from ptleg--negative area between mu=, 1 1p,e12.5,5h and ,e12.5,3h, ,e12.5,5h, e=,e12.5/ *d acer.2220 2 format(/41h message from ptleg--mt.ne.zero--card is / *d acer.2222,2223 3 format(/53h message from ptleg--integrated area of legendre dist, 1 /,20x,28h using 1000 subintervals is ,1p,e12.6/ *d acer.2582 1 20x,27husing 1000 subintervals is ,1p,e12.6/ */ fix date in output file *d acer.3683,3684 call dater(hd) */ fix up case of no photon production *d acer.4358 600 if (mf1x(1).eq.0.and.mf1x(2).eq.0) go to 635 *d acer.4765 *i acer.4774 len2=lendp */ fix up type 3 output and input *d acer.5695 dimension i0(202),f0(202) equivalence (i0(1),f0(1)) *d acer.5758 i1(6)=41 *d acer.5794 f0(2)=zaid *d acer.5797 i0(5)=i1(20)+202+3 *d acer.5850 call wdisk(nout,i2(l+1),n+1,ll) *d acer.8027,8029 *if wordio call assign(nin,0,0) call rdisk(nin,i1(1),40,202) if (unit(nin)) 161,250,250 *endif *d acer.8074,8075 *if wordio call rdisk(nin,f2(1),n,243) if (unit(nin)) 163,250,250 *endif *d acer.8081,8082 *if wordio call rdisk(nin,f2(l+1),n,ll) if (unit(nin)) 164,250,250 *endif *i acer.8140 *if wordio 250 call error('acefix','rdisk/wdisk error',' ') *endif c *d acer.8146,8148 */ fix up dosimetry output (vontobel, psi) *d acer.7269 common /jxst/ lone,dum2,mtr,dum4,dum5,lsig,sigd,jxsd(14),end, 1 jxsd2(10) *d acer.7316 jscr=iscr call tab1io(nin,0,0,a(jscr),nb,nw) *d acer.7319 intr=nint(a(iscr+7)) 112 jscr=jscr+nw if (nb.eq.0) go to 115 call moreio(nin,0,0,a(jscr),nb,nw) go to 112 115 if (nr.eq.1.and.intr.eq.2) go to 130 *d acer.7335,7336 xss(l+ne)=a(k+1) l=l+1 *i acer.7338 l=l+ne *d acer.7346 nn=l1+mtr-2 *d acer.7350 n=l-sigd *d acer.7361 write(hz,'(f9.2,''y'')') zaid *d acer.7368 call dosout(itype,nace,ndir) *d acer.7381 common /jxst/ lone,dum2,mtr,dum4,dum5,lsig,sigd,jxsd(14),end, 1 jxsd2(10) *d acer.7395 l=nint(xss(lsig-1+j)+sigd-1) *d acer.7410,7411 s=xss(l+ne) l=l+1 *d acer.7457 common /jxst/ lone,dum2,mtr,dum4,dum5,lsig,sigd,jxsd(14),end, 1 jxsd2(10) *d acer.7459 equivalence (nxs(1),len2),(jxs(1),lone) *i acer.7503 call typen(0,nout,3) call closz(nout) */ fix one of the acer prompts (mattes) *d acer.216 1 '('' enter iopt, iprint, ntype, suff, nxtra'')') */ fix temperature input (mattes) *d acer.299 tempd=z(2) */ tighten up initial tolerance for integral thinning *d acer.1096 data fac/2./, rfact/.0025/, miter/7/ *d acer.1324 if (iter.eq.1.and.j.lt..75*npts) rmin=rmin/(fac**6) */ file 6 patches to topfil *d acer.1539 data nwmaxn/14000/ *d acer.1548 nwmax=nwmaxn *d acer.1570,1572 if (mf.eq.6.and.mt.le.91) go to 121 *d acer.1612,1617 128 ik=ik+1 *d acer.1624,1625 140 if (mf.eq.5.and.lf.eq.5) go to 185 if (mf.eq.5.and.lf.ne.1) go to 200 if (mf.eq.6.and.lf.eq.6) go to 260 *i acer.1731 if (mf.eq.6) call asend(0,nout) *i acer.1773 c c ***copy phase-space law from file 6. 260 call contio(nin,nout,0,a(iscr),nb,nw) if (ik.lt.nk) go to 128 go to 100 */ prevent writing of cont in conver if mf.gt.15 is present *i acer.2998 if (mfh.gt.15) go to 515 */ generalize transfer of mf6 photons to mf16 *d acer.3279 530 if (law.eq.0.or.law.eq.3.or.law.eq.4) go to 558 if (law.eq.6) go to 551 if (law.eq.7) go to 552 call tab2io(nin,0,0,a(iscr+6),nb,nw) *d acer.3296 go to 558 551 call contio(nin,0,0,a(iscr+6),nb,nw) go to 558 552 call tab2io(nin,0,0,a(iscr+6),nb,nw) ne=n2h do 556 ie=1,ne nmu=n2h do 554 imu=1,nmu call tab1io(nin,0,0,a(iscr+6),nb,nw) 553 if (nb.eq.0) go to 554 call moreio(nin,0,0,a(iscr+6),nb,nw) go to 553 554 continue 556 continue 558 ik=ik+1 */ define nin in gamout (vontobel) *i acer.3370 nin=ngend */ allow for larger nu tables (mattes) *d acer.3664 nwscr=10000 *i acer.3732 inow=iscr 121 inow=inow+nw if (nb.eq.0) go to 123 call moreio(nin,0,0,a(inow),nb,nw) go to 121 123 continue */ supress extra 1e-5 at start of table (mattes) *d acer.3780 */ move inactive message *d acer.3774,3776 165 esz=0 *i acer.3802 if ((esz+5*nes).gt.max3) 1 call error('acelod','insufficient storage for esz block.',' ') */ patches for file 6 in acelod *d acer.3915 c ***store angular distributions *d acer.3921 if (math.eq.-1) go to 335 if (mfh.gt.6) go to 335 if (mfh.eq.0) go to 290 if (mfh.eq.5) go to 330 *i acer.3922 if (mfh.eq.6) go to 292 *d acer.3925 go to 296 292 lct=l2h call tab1io(nin,0,0,a(iscr),nb,nw) law=l2h 296 call tab2io(nin,0,0,a(iscr),nb,nw) *d acer.3937 300 if (mfh.eq.4) go to 304 if (mfh.eq.6.and.law.eq.2) go to 302 xss(land+ir)=-1 go to 325 302 do 303 k=1,nxc if (mts(k).eq.mth.and.mfs(k).eq.6) mfs(k)=64 303 continue 304 last=next *d acer.4237,4238 *d acer.4240 *d acer.4255 c c ***store yield 530 xss(ldlw+i-1)=next-dlw+1 xss(land+i)=-1 last=next xss(next)=0 xss(next+1)=44 next=next+3 if (m.eq.1.and.jnt.eq.2) go to 540 *i acer.4264 ntyr=nint(a(iscr+5+2*m+2)) xss(tyr+i-1)=(3-2*lct)*ntyr *d acer.4267,4269 550 xss(next+n+k)=1.0 *d acer.4272,4285 c ***store energy-angle distribution *d acer.4296 xss(last+2)=next-last+1 */ vontobel fixes for mf16 photons *d acer.4670 xss(nex+2)=10 *d acer.4683,4684 xss(nex+5)=a(iscr+6+2*m)*1.e-6 xss(nex+6)=a(iscr+4+2*m+2+n)*1.e-6 *d acer.4688 *d acer.4690 lep=l2h */ allow for coupled energy-angle distributions (mattes) *d acer.6002 if (nn.le.0) go to 205 *i acer.6033 if (nout.ne.1) ly=nint(xss(tyr+i-1)) if (nout.eq.1) ly=iss(tyr+i-1) ly=iabs(ly) if (ly.gt.100) then l=ly-100+dlw-1 if (nout.ne.1) nr=nint(xss(l)) if (nout.eq.1) nr=iss(l) call typen(l,nout,1) l=l+1 if (nr.eq.0) go to 212 n=2*nr do 211 j=1,n call typen(l,nout,1) 211 l=l+1 212 if (nout.ne.1) ne=nint(xss(l)) if (nout.eq.1) ne=iss(l) call typen(l,nout,1) l=l+1 n=2*ne do 213 j=1,n call typen(l,nout,2) 213 l=l+1 endif */ missing increments in law 1 section (mattes) *i acer.6073 l=l+1 *i acer.6081 l=l+1 *i acer.6087 l=l+1 */ add branch for mf16 data (mattes) *i acer.6293 if (mftype.eq.16) go to 465 */ fix typen to prevent two output lines for i=4 (vontobel) *d acer.6501 save hl,i *d acer.6503 if (iflag.eq.3.and.nout.gt.1.and.i.lt.4) then */ reintroduce lost lines (mattes) *i acer.6654 if (emax.gt.em) emax=em c c ***process coherent reaction if requested if (mte.gt.0) go to 120 write(nout) mte,mte go to 200 120 if (ielas.gt.0) go to 150 call findf(matd,3,mte,nin) call contio(nin,0,0,a(iscr),nb,nw) */ part of thermal inelastic missing (lazo) *d acer.7202 n=nieb*(nil+2) */ thermal zaid should be a name, not a number *i acer.123 c * tname thermal zaid name (6 char max, def=za) *i acer.163 character*6 tname,tscr common/ace3a/tname *d acer.294,295 1 '(/'' enter matd, tempd, tname for thermal data'')') nz=3 *if sw nz=4 *endif *d acer.297 iz(3)=iblank iz(4)=iblank call free(nsysi,z,nz,6) *d acer.300 *if sw write(tscr,'(a4,a2)') iz(3),iz(4) *else write(tscr,'(a6)') iz(3) *endif nch=0 do 132 i=1,6 if (tscr(i:i).eq.' ') go to 133 nch=nch+1 132 continue 133 tname=' ' if (nch.gt.0) tname(7-nch:6)=tscr(1:nch) write(nsyso,27) matd,tempd,tname *d acer.419 1 40h temperature .......................... ,1p,e10.2/ 2 40h thermal name ......................... ,4x,a6) *i acer.6882 character*6 tname common/ace3a/tname *d acer.6903,6904 write(hz,'(a6,f3.2,''t'')') tname,suff *d acer.7123 10 format(1h1///////38x,4hzaid,1x,a10/39x,3hawr,f10.3/ */ increase the size of the container arrays *d acer.168 common/astore/a(20000) *d acer.176 data namax/20000/, maxnpb/50/, nidmax/27/ *d acer.188 *d acer.190,191 max1=0 max2=0 max3=400000 *d acer.432 common/astore/a(20000) *d acer.3651 common/astore/a(20000) *d acer.4800 common/astore/a(20000) *d acer.6580 common/astore/a(20000) *d acer.6885 common/astore/a(20000) *d acer.7270 common/astore/a(20000) *d acer.7544 common/astore/a(20000) *d acer.3641 common /xsst/ n3,xss(400000) *d acer.4792 common /xsst/ n3,xss(400000) *d acer.5693 common /xsst/ n3,xss(400000) *d acer.5893 common /xsst/ n3,xss(400000) *d acer.6497 common /xsst/ n3,xss(400000) *d acer.6877 common /xsst/ n3,xss(400000) *d acer.7004 common /xsst/ n3,xss(400000) *d acer.7173 common /xsst/ n3,xss(400000) *d acer.7274 common /xsst/ n3,xss(400000) *d acer.7382 common /xsst/ n3,xss(400000) *d acer.7460 common /xsst/ n3,xss(400000) *d acer.7551 common /xsst/ n3,xss(400000) *d acer.7807 common /xsst/ n3,xss(400000) *d acer.7903 common /xsst/ n3,xss(400000) *d acer.7978 common /xsst/ n3,xss(400000) */ make sure mt600,etc are available for mf6 photons *d acer.494 common/phot/ntrp,nf12s,mf12s(10),nf16s,mf16s(250) *d acer.1090 common/phot/ntrp,nf12s,mf12s(10),nf16s,mf16s(250) *d acer.1443 252 if (nf12s.eq.0) go to 255 *i acer.1446 255 if (nf16s.eq.0) go to 260 do 256 i=1,nf16s if (mt.eq.mf16s(i)) go to 215 256 continue *d acer.2664 common/phot/ntrp,nf12s,mf12s(10),nf16s,mf16s(250) *i acer.2924 common/phot/ntrp,nf12s,mf12s(10),nf16s,mf16s(250) *i acer.2962 nf16s=0 *i acer.3267 if (mth.ge.600) nf16s=nf16s+1 if (mth.ge.600) mf16s(nf16s)=mth *d acer.3649 common /phot/ ntrpp,nf12s,mf12s(10),nf16s,mf16s(250) *i acer.3692 if16s=0 *d acer.3705 ntr=ntr+nf12s+nf16s *d acer.3877 if (if12s.eq.0.and.if16s.eq.0) xss(it+j)=xss(it+j)+s *d acer.3889 if (if12s.gt.0.or.if16s.gt.0) go to 255 *d acer.3894 if (if12s.gt.nf12s) go to 256 *d acer.3898 256 if16s=if16s+1 if (if16s.gt.nf16s) go to 260 call findf(matd,3,mf16s(if16s),nin) call contio(nin,0,0,a(iscr),nb,nw) go to 205 260 if (nf12s.eq.0.and.nf16s.eq.0) go to 265 */ fix up reaction naming *i acer.5631 common/util/npage,iverf *d acer.5633 character*10 hndf(400) character*10 hndf1(50),hndf2(48),hndf3(51) character*10 hndf4(50),hndf5(50),hndf6(50) character*10 hndf7(50),hndf8(50) character*10 h719,h739,h759,h779,h799 *i acer.5635 equivalence (hndf4(1),hndf(150)) equivalence (hndf5(1),hndf(200)) equivalence (hndf6(1),hndf(250)) equivalence (hndf7(1),hndf(300)) equivalence (hndf8(1),hndf(350)) *d acer.5646 a '(n,x) ', '(n,2np) ', '(n,3np) ', '(n,x) ', *d acer.5667 5 '(n,pd) ', '(n,pt) ', '(n,x) ', '(n,x) ', *d acer.5675 d '(n,x) ', '(n,x) ', '(n,x) '/ data hndf4/'(n,p*0) ', 1 '(n,p*1) ', '(n,p*2) ', '(n,p*3) ', '(n,p*4) ', 2 '(n,p*5) ', '(n,p*6) ', '(n,p*7) ', '(n,p*8) ', 3 '(n,p*9) ', '(n,p*10) ', '(n,p*11) ', '(n,p*12) ', 4 '(n,p*13) ', '(n,p*14) ', '(n,p*15) ', '(n,p*16) ', 5 '(n,p*17) ', '(n,p*18) ', '(n,p*19) ', '(n,p*20) ', 6 '(n,p*21) ', '(n,p*22) ', '(n,p*23) ', '(n,p*24) ', 7 '(n,p*25) ', '(n,p*26) ', '(n,p*27) ', '(n,p*28) ', 8 '(n,p*29) ', '(n,p*30) ', '(n,p*31) ', '(n,p*32) ', 9 '(n,p*33) ', '(n,p*34) ', '(n,p*35) ', '(n,p*36) ', a '(n,p*37) ', '(n,p*38) ', '(n,p*39) ', '(n,p*40) ', b '(n,p*41) ', '(n,p*42) ', '(n,p*43) ', '(n,p*44) ', c '(n,p*45) ', '(n,p*46) ', '(n,p*47) ', '(n,p*48) ', d '(n,p*c) '/ data hndf5/'(n,d*0) ', 1 '(n,d*1) ', '(n,d*2) ', '(n,d*3) ', '(n,d*4) ', 2 '(n,d*5) ', '(n,d*6) ', '(n,d*7) ', '(n,d*8) ', 3 '(n,d*9) ', '(n,d*10) ', '(n,d*11) ', '(n,d*12) ', 4 '(n,d*13) ', '(n,d*14) ', '(n,d*15) ', '(n,d*16) ', 5 '(n,d*17) ', '(n,d*18) ', '(n,d*19) ', '(n,d*20) ', 6 '(n,d*21) ', '(n,d*22) ', '(n,d*23) ', '(n,d*24) ', 7 '(n,d*25) ', '(n,d*26) ', '(n,d*27) ', '(n,d*28) ', 8 '(n,d*29) ', '(n,d*30) ', '(n,d*31) ', '(n,d*32) ', 9 '(n,d*33) ', '(n,d*34) ', '(n,d*35) ', '(n,d*36) ', a '(n,d*37) ', '(n,d*38) ', '(n,d*39) ', '(n,d*40) ', b '(n,d*41) ', '(n,d*42) ', '(n,d*43) ', '(n,d*44) ', c '(n,d*45) ', '(n,d*46) ', '(n,d*47) ', '(n,d*48) ', d '(n,d*c) '/ data hndf6/'(n,t*0) ', 1 '(n,t*1) ', '(n,t*2) ', '(n,t*3) ', '(n,t*4) ', 2 '(n,t*5) ', '(n,t*6) ', '(n,t*7) ', '(n,t*8) ', 3 '(n,t*9) ', '(n,t*10) ', '(n,t*11) ', '(n,t*12) ', 4 '(n,t*13) ', '(n,t*14) ', '(n,t*15) ', '(n,t*16) ', 5 '(n,t*17) ', '(n,t*18) ', '(n,t*19) ', '(n,t*20) ', 6 '(n,t*21) ', '(n,t*22) ', '(n,t*23) ', '(n,t*24) ', 7 '(n,t*25) ', '(n,t*26) ', '(n,t*27) ', '(n,t*28) ', 8 '(n,t*29) ', '(n,t*30) ', '(n,t*31) ', '(n,t*32) ', 9 '(n,t*33) ', '(n,t*34) ', '(n,t*35) ', '(n,t*36) ', a '(n,t*37) ', '(n,t*38) ', '(n,t*39) ', '(n,t*40) ', b '(n,t*41) ', '(n,t*42) ', '(n,t*43) ', '(n,t*44) ', c '(n,t*45) ', '(n,t*46) ', '(n,t*47) ', '(n,t*48) ', t '(n,t*c) '/ data hndf7/'(n,he3*0) ', 1 '(n,he3*1) ', '(n,he3*2) ', '(n,he3*3) ', '(n,he3*4) ', 2 '(n,he3*5) ', '(n,he3*6) ', '(n,he3*7) ', '(n,he3*8) ', 3 '(n,he3*9) ', '(n,he3*10)', '(n,he3*11)', '(n,he3*12)', 4 '(n,he3*13)', '(n,he3*14)', '(n,he3*15)', '(n,he3*16)', 5 '(n,he3*17)', '(n,he3*18)', '(n,he3*19)', '(n,he3*20)', 6 '(n,he3*21)', '(n,he3*22)', '(n,he3*23)', '(n,he3*24)', 7 '(n,he3*25)', '(n,he3*26)', '(n,he3*27)', '(n,he3*28)', 8 '(n,he3*29)', '(n,he3*30)', '(n,he3*31)', '(n,he3*32)', 9 '(n,he3*33)', '(n,he3*34)', '(n,he3*35)', '(n,he3*36)', a '(n,he3*37)', '(n,he3*38)', '(n,he3*39)', '(n,he3*40)', b '(n,he3*41)', '(n,he3*42)', '(n,he3*43)', '(n,he3*44)', c '(n,he3*45)', '(n,he3*46)', '(n,he3*47)', '(n,he3*48)', h '(n,he3*c) '/ data hndf8/'(n,a*0) ', 1 '(n,a*1) ', '(n,a*2) ', '(n,a*3) ', '(n,a*4) ', 2 '(n,a*5) ', '(n,a*6) ', '(n,a*7) ', '(n,a*8) ', 3 '(n,a*9) ', '(n,a*10) ', '(n,a*11) ', '(n,a*12) ', 4 '(n,a*13) ', '(n,a*14) ', '(n,a*15) ', '(n,a*16) ', 5 '(n,a*17) ', '(n,a*18) ', '(n,a*19) ', '(n,a*20) ', 6 '(n,a*21) ', '(n,a*22) ', '(n,a*23) ', '(n,a*24) ', 7 '(n,a*25) ', '(n,a*26) ', '(n,a*27) ', '(n,a*28) ', 8 '(n,a*29) ', '(n,a*30) ', '(n,a*31) ', '(n,a*32) ', 9 '(n,a*33) ', '(n,a*34) ', '(n,a*35) ', '(n,a*36) ', a '(n,a*37) ', '(n,a*38) ', '(n,a*39) ', '(n,a*40) ', b '(n,a*41) ', '(n,a*42) ', '(n,a*43) ', '(n,a*44) ', c '(n,a*45) ', '(n,a*46) ', '(n,a*47) ', '(n,a*48) ', d '(n,a*c) '/ data h719/'(n,p*c)x '/ data h739/'(n,d*c)x '/ data h759/'(n,t*c)x '/ data h779/'(n,he3*c)x'/ data h799/'(n,a*c)x '/ *d acer.5677 if (iverf.ge.6) then i=mt if (i.ge.600) i=i-450 name=hndf(i) else if (mt.lt.150) then name=hndf(mt) else if (mt.ge.700.and.mt.lt.718) then name=hndf(mt-550) else if (mt.eq.718) then name=hndf(199) else if (mt.eq.719) then name=h719 else if (mt.ge.720.and.mt.lt.738) then name=hndf(mt-520) else if (mt.eq.738) then name=hndf(249) else if (mt.eq.739) then name=h739 else if (mt.ge.740.and.mt.lt.758) then name=hndf(mt-490) else if (mt.eq.758) then name=hndf(299) else if (mt.eq.759) then name=h759 else if (mt.ge.760.and.mt.lt.779) then name=hndf(mt-460) else if (mt.eq.778) then name=hndf(349) else if (mt.eq.779) then name=h779 else if (mt.ge.780.and.mt.lt.798) then name=hndf(mt-430) else if (mt.eq.798) then name=hndf(399) else if (mt.eq.799) then name=h799 endif endif *ident up26 */ dtfr -- 9 dec 91 -- fix pointer for photon prod xsec *d dtfr.443 660 loca=lz+ip+nl*((jz-1)+nz) */ fix reaction list for pointwise plotting *i dtfr.874 common/util/npage,iverf *i dtfr.879 common/dtf4/dum(4),iphph *d dtfr.926,929 135 if (iphph.eq.1) then if (iverf.ge.6) then if (mt.lt.501.or.mt.gt.522) go to 100 else if (mt.lt.501.or.mt.gt.602) go to 100 endif else if (iverf.ge.6) then if (mt.gt.150.and.mt.lt.301) go to 100 if (mt.gt.450.and.mt.lt.600) go to 100 if (mt.gt.849) go to 100 else if (mt.gt.150.and.mt.lt.301) go to 100 if (mt.gt.450.and.mt.lt.700) go to 100 if (mt.gt.799) go to 100 endif endif *ident up27 */ matxsr -- 21 dec 91 -- allow bcd output, modify format. */ remove print and index in favor of bbc. c * the bbc code can be used to convert between binary * c * and bcd modes, list, index, or modify a matxs file. * c * * *d matxsr.19 *d matxsr.27 *d matxsr.34 *i matxsr.77 c c (bcd format changed 12/21/91) *d matxsr.138 cb format(4h 0v ,a8,1h*,2a8,1h*,i6) *d matxsr.157 cb format(6h 1d ,6i6) *d matxsr.178 cb format(4h 2d /(9a8)) *d matxsr.195 cb format(4h 3d ,4x,8a8/(9a8)) hprt,htype,hmatn *d matxsr.241 cb format(4h 4d ,8x,1p,5e12.5/(6e12.5)) *d matxsr.257 cb format(4h 5d ,a8,1p,2e12.5/(2e12.5,5i6)) *d matxsr.280 cb format(4h 6d ,4x,8a8/(9a8)) hvps *d matxsr.307 cb format(4h 7d ,8x,1p,5e12.5/(6e12.5)) *d matxsr.325 cb format(4h 8d ,4x,a8/(12i6)) hmtx,lord,jconst, *d matxsr.345 cb format(4h 9d ,8x,1p,5e12.5/(6e12.5)) *d matxsr.374 cb format(4h10d ,8x,1p,5e12.5/(6e12.5)) *d matxsr.400,402 dimension hz(1),iz(1) equivalence (z(1),hz(1)),(z(1),iz(1)) data blank/8h / *d matxsr.434 *d matxsr.456 call openz(nmatx,1) *d matxsr.462,473 nwds=2*mult+1 if (ntty.gt.0) write(ntty, 1 '(/'' enter matxs version number and user id.'')') z(1)=0. z(2)=blank if (mult.gt.1) z(3)=blank z(2+mult)=blank if (mult.gt.1) z(3+mult)=blank call free(nsysi,z(1),nwds,ncw) ivers=nint(z(1)) nw=nwds-1 do 110 i=1,nw 110 iz(i)=iz(i+1) huse(1)=hz(1) huse(2)=hz(2) write(nsyso,50) ivers,huse *d matxsr.478,484 *d matxsr.486,489 call closz(nmatx) *d matxsr.521,523 50 format(40h file version number .................. ,i10/ 1 40h user id .............................. ,2a8) *d matxsr.574 if (nmatx.lt.0) call rite(-nmatx,irec,a(iscrt),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 0v '',a8,1h*,2a8,1h*,i6)') 1 hfile,huse(1),huse(2),ivers *d matxsr.582 if (nmatx.lt.0) call rite(-nmatx,irec,a(icont),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 1d '',6i6)') 1 (ia(icont+i-1),i=1,nwds) *d matxsr.588 ihholl=(iholl-1)/mult if (nmatx.lt.0) call rite(-nmatx,irec,a(iholl),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 2d ''/(9a8))') 1 (ha(ihholl+i),i=1,nholl) *d matxsr.594 if (nmatx.lt.0) call rite(-nmatx,irec,a(ifild),nwds,0) if (nmatx.gt.0) then ndr=npart+ntype+nmat nir=npart+2*ntype+2*nmat ihfild=(ifild-1)/mult write(nmatx,'('' 3d '',4x,8a8/(9a8))') (ha(ihfild+i),i=1,ndr) write(nmatx,'(12i6)') (ia(ifild+ndr*mult+i-1),i=1,nir) endif *d matxsr.603 if (nmatx.lt.0) call rite(-nmatx,irec,a(lout),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 4d '',1p,8x,5e12.5/(6e12.5))') 1 (a(lout+j-1),j=1,nwds) *d matxsr.618 if (nmatx.lt.0) call rite(-nmatx,irec,a(imatc),nwds,0) if (nmatx.gt.0) then ihmatc=(imatc-1)/mult+1 write(nmatx,'('' 5d '',a8,1p,e12.5)') 1 ha(ihmatc),a(imatc+mult) do 175 i=1,nsubj ll=imatc+mult+1+6*(i-1) write(nmatx,'(1p,2e12.5,4i6)') a(ll),a(ll+1), 1 ia(ll+2),ia(ll+3),ia(ll+4),ia(ll+5) 175 continue endif *d matxsr.636 if (nmatx.lt.0) call rite(-nmatx,irec,a(ivcon),nwds,0) if (nmatx.gt.0) then nw=2*n1d ihvcon=(ivcon-1)/mult write(nmatx,'('' 6d '',4x,8a8/(9a8))') 1 (ha(ihvcon+i),i=1,n1d) write(nmatx,'(12i6)') (ia(ivcon+n1d*mult+i-1),i=1,nw) endif *d matxsr.651 if (nmatx.lt.0) call rite(-nmatx,irec,a(ivdat),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 7d '',8x,1p,5e12.5/(6e12.5))') 1 (a(ivdat+i-1),i=1,nwds) *d up8.65 if (nmatx.lt.0) call rite(-nmatx,irec,a(ivdat),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 7d '',8x,1p,5e12.5/(6e12.5))') 1 (a(ivdat+i-1),i=1,nwds) *d matxsr.669 if (nmatx.lt.0) call rite(-nmatx,irec,a(imcon),nwmcr,0) if (nmatx.gt.0) then ihmcon=(imcon-1)/mult+1 nw=nwmcr-mult write(nmatx,'('' 8d '',4x,a8)') ha(ihmcon) write(nmatx,'(12i6)') (ia(imcon+mult+i-1),i=1,nw) endif *d matxsr.689 if (nmatx.lt.0) call rite(-nmatx,irec,a(imdat),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 9d '',8x,1p,5e12.5/(6e12.5))') 1 (a(imdat+i-1),i=1,nwds) *d up8.73 if (nmatx.lt.0) call rite(-nmatx,irec,a(imdat),nwds,0) if (nmatx.gt.0) write(nmatx,'('' 9d '',8x,1p,5e12.5/(6e12.5))') 1 (a(imdat+i-1),i=1,nwds) *d matxsr.696 if (nmatx.lt.0) call rite(-nmatx,irec,a(icdat),nwdc,0) if (nmatx.gt.0) write(nmatx,'(''10d '',8x,1p,5e12.5/(6e12.5))') 1 (a(icdat+i-1),i=1,nwdc) *d matxsr.2164,2710 *i matxsr.2850 c *ident up28 */ reconr -- 6 mar 91 -- avoid the reserved name "save" *d reconr.965 tsave=x(j) *d reconr.967 x(i)=tsave */ fix the value of boltzmann's constant *d reconr.2004 data rpi/1.7724585/, bk/8.61735e-5/ *d reconr.3125 data bk/8.61735e-5/, rpi/1.7724585/, ki/1/ *ident up29 */ heatr -- 5 mar 92 -- avoid the reserved name "save" *d heatr.2978 tsave=a(ieg+i-1) *d heatr.2980,2981 a(ieg+ii-1)=tsave tsave=a(iy+i-1) *d heatr.2983,2984 a(iy+ii-1)=tsave tsave=a(ies+i-1) *d heatr.2986 a(ies+ii-1)=tsave *ident up30 */ groupr -- 5 mar 92 -- fix format for message *d groupr.446 381 write(strng,'(''illegal mfd ='',i3)') mfd */ groupr -- 5 mar 92 -- avoid the variable name "save" (cea-dmt) *d groupr.6383 tsave=a(ieg+i-1) *d groupr.6385,6386 a(ieg+ii-1)=tsave tsave=a(ies+i-1) *d groupr.6388,6389 a(iy+ii-1)=tsave tsave=a(ies+i-1) *d groupr.6391 a(ies+ii-1)=tsave */ groupr -- 5 mar 92 -- fix tabulated mf6 option (konieczny) *d groupr.4802 if (llo.eq.ilo.or.llo.gt.mlo) go to 320 *d groupr.4827 jj=lhi-nchi+2 *d groupr.4829 if (chi(jj).gt.qp(i)) go to 345 *d groupr.4833 jj=lhi+2 *d groupr.4835 if (chi(jj).gt.qp(i)) go to 355 *i groupr.4844 do 380 l=1,nl 380 call terp1(elo,term1(l),ehi,term2(l),e,term(l),intt) */ groupr -- 5 mar 92 -- fix dimensions (raepsaet, cea-dmt) *d groupr.4264 dimension cnow(10),term(1),p(7) *d groupr.4409 dimension cnow(10),p(8) *d groupr.4671 dimension term(nl),clo(10),chi(10) *d groupr.7430 dimension c(10) */ groupr -- 5 mar 92 -- supress excess upscatter messages *i groupr.6878 nupm=0 *d groupr.6957 340 if (sigup.gt..1.and.nupm.lt.10) then nupm=nupm+1 *i groupr.6959 if (nupm.eq.10) 1 call mess('getsed','additional messages supressed',' ') *ident up31 */ gaminr -- 5 mar 92 -- include endf-6 gheat *i gaminr.1022 if (mtd.eq.522) ff(1,2)=e *ident up32 */ acer -- 10 mar 92 -- convert legendre law1 to kalbach form *i acer.4242 awp=a(iscr+1) *d acer.4288,4289 if (lang.lt.1.or.(lang.gt.2.and.lang.lt.11).or.lang.gt.13) 1 call error('acelod', 2 'only lang=1, 2, or 11-13 allowed for endf-6 file 6 neutrons', 3 ' ') *i acer.4309 aq4=awr+1.-awp beta=sqrt(awr*aq4/awp) dodmup=1. dodmum=1. if (lct.eq.1) then dodmup=(beta/(beta+1.))**2 dodmum=(beta/(beta-1.))**2 endif dompl=dodmup+dodmum domin=dodmup-dodmum c c ***loop over incident energies *i acer.4317 ee=c2h*1.e-6 *i acer.4328 c c ***loop over secondary energies *i acer.4331 bzro=xss(nexd+ki+n)*1.e-6 *d acer.4341 *d acer.4343 if (lang.ne.2) go to 576 c c ***distribution given in kalbach format xss(ki+3*n+nexd)=a(iscr+8+ncyc*(ki-1)) xss(ki+4*n+nexd)=0.040*(ep+esep) go to 580 576 if (lang.ne.1) go to 577 c c ***convert legendre distribution to kalbach form iso=1 do 578 ik3=1,na if (a(iscr+7+ik3+ncyc*(ki-1)).ne.0.) iso=0 578 continue if (bzro.eq.0.) iso=1 if (iso.eq.0) go to 579 c c isotropic case xss(ki+3*n+nexd)=0. xss(ki+4*n+nexd)=1.e-5 go to 580 c c nonisotropic case 579 sfe=0.5 sfo=0.0 do 581 ii1=1,na bbi=(2.*ii1+1.)*0.5*a(iscr+7+ii1+ncyc*(ki-1))/bzro if (ii1.eq.2*(ii1/2)) then sfe=sfe+bbi else sfo=sfo+2*bbi endif 581 continue if (lct.eq.1) then fbarl=sfe delfl=sfo fbarcm=0.5*(fbarl*dompl+0.5*delfl*domin) delfcm=fbarl*domin+0.5*delfl*dompl else fbarcm=sfe delfcm=sfo endif c don't accept unreasonable distributions if (delfcm.lt.0.) delfcm=0. call fndar1(akal,rkal,fbarcm,delfcm,ee,ep) xss(ki+3*n+nexd)=rkal xss(ki+4*n+nexd)=akal go to 580 c ***convert tabulated distribution to kalbach form 577 emu1=a(iscr+8+ncyc*(ki-1)) emu2=a(iscr+6+na+ncyc*(ki-1)) fbl=a(iscr+9+ncyc*(ki-1)) ffl=a(iscr+7+na+ncyc*(ki-1)) if (abs(emu1+1.).lt.1.e-3.and.abs(emu2-1.).lt.1.e-3) go to 582 fbcm=fbl ffcm=ffl if (lct.eq.1) then emu1=((emu1**2-1.)+emu1*sqrt(emu1**2-1.+beta**2))/beta emu2=((emu2**2-1.)+emu2*sqrt(emu2**2-1.+beta**2))/beta dodmup=beta*beta*(beta+emu2)/sqrt((beta*beta+ 1 2.*beta*emu2+1.)**3) dodmum=beta*beta*(beta*emu1)/sqrt((beta*beta+ 1 2.*beta*emu1+1.)**3) fbcm=fbl*dodmum ffcm=ffl*dodmup endif call fndar2(akak,rkal,fbcm,ffcm,emu1,emu2,ee,ep) call mess('acelod','tabulated angular distribution', 1 'does not extend over entire cosine range.') go to 580 582 if (lct.eq.1) then fbarl=0.5*(ffl+fbl) delfl=ffl-fbl fbarcm=0.5*(fbarl*dompl+0.5*delfl*domin) delfcm=fbarl*domin+0.5*delfl*dompl else fbarcm=0.5*(ffl+fbl) delfcm=ffl-fbl endif call fndar1(akal,rkal,fbarcm,delfcm,ee,ep) xss(ki+3*n+nexd)=rkal xss(ki+4*n+nexd)=akal c c ***continue secondary-energy loop *i acer.4350 c c ***continue incident-energy loop *i acer.4779 c subroutine fndar1(akal,rkal,fbarc,delfc,ee,ep) c ****************************************************************** c this subroutine finds the kalbach angular distribution parameters c akal and rkal from the average of the forward and backward c amplitudes fbarc, and the difference of the forward and backward c amplitudes delfc, both defined in the cm frame. c writen by a. j. sierk, lanl, 1 march 1990. c ****************************************************************** iq1=0 akal=2.*fbarc*(1.-4.*fbarc**2/3.) c c ***use newton's method to converge to a solution of fofa=0, c ***where fofa=akal/tanh(akal)-2*fbarc 110 fofa=akal/tanh(akal)-2.*fbarc if (abs(fofa).lt.1.e-10) go to 120 s2a=sinh(2.*akal) c2a=cosh(2.*akal) dela=(2.*fbarc*(c2a-1.)-akal*s2a)/(s2a-2.*akal) akal=akal+dela if (akal.lt.0.) akal=-akal iq1=iq1+1 if (iq1.lt.30) go to 110 call mess('fndar1','loop to find kalbach a not converged',' ') return 120 if (akal.ge.0.) go to 130 call mess('fndar1','kalbach a is negative',' ') return 130 rkal=0. if (akal.ne.0.) rkal=delfc/akal if (rkal.lt.0..or.rkal.gt.1.) 1 call mess('fndar1','kalbach r is unreasonable',' ') return end c subroutine fndar2(akal,rkal,fbcm,ffcm,emu1,emu2,ee,ep) c ****************************************************************** c this subroutine finds the kalbach angular distribution parameters c akal and rkal from tabulated angular distributions using the lang c =11--13 options. this subroutine is only used when the table c does not cover the entire range (-1,1) in cos(theta-cm). c ****************************************************************** iq1=0 fb=fbcm+ffcm dmu=emu2-emu1 akal=fb*(1.-fb**2/3.) 110 cam=cosh(akal*emu1) cap=cosh(akal*emu2) sam=sinh(akal*emu1) sap=sinh(akal*emu2) sa=sinh(akal) ca=cosh(akal) cad=cosh(akal*dmu) sad=sinh(akal*dmu) fofa=ffcm*sam-fbcm*sap+0.5*akal*sad/sa if (abs(fofa).le.1.e-10) go to 120 denom=sad*(sa-akal*ca)+akal*dmu*cad*sa+ 1 2.*sa*sa*(ffcm*emu1*cam-fbcm*emu2*cap) dela=(2.*(fbcm*sap-ffcm*sam)*sa-akal*sad)*sa/denom akal=akal+dela if (iq1.lt.30) go to 110 call mess('fndar2','loop to find kalbach a not converged',' ') return 120 if (akal.ge.0.) go to 130 call mess('fndar2','kalbach a is negative',' ') return 130 rkal=0. if (akal.ne.0.) rkal=2.*sa*(ffcm*cam-fbcm*cap)/(akal*s2a) if (rkal.lt.0..or.rkal.gt.1.) 1 call mess('fndar2','kalbach r has an unreasonable value',' ') return end *ident up33 */ unresr -- 11 mar 91 -- fix multi-isotope case *d unresr.398 *d unresr.406 *d unresr.431 if (ier.lt.ner) go to 115 go to 110 *ident up34 */ wimsr -- 11 mar 91 -- fix bugs (mann, hanford) *d wimsr.914,915 if (a(nzmax*jtem+iterm-jz+nsigz+locsf).gt.0.) then a(nzmax*jtem+iterm+locnus)=a(nzmax*jtem+iterm+locsf)*a(locn) 1 /a(nzmax*jtem+iterm-jz+nsigz+locsf) else a(nzmax*jtem+iterm+locnus)=0. endif *d wimsr.968 a(loc+nsigz-iz)=a(loc+nsigz-iz)*a(isigz+nsigz-iz) *d wimsr.1164 ident=nint(a(idnt)) */ wimsr -- 15 may 92 -- fix problems with fission matrix *d wimsr.1687 a(locn)=a(locn)+a(icspc+k-1)*a(loca) *d wimsr.1689 locc=locchi+ngnd-ig2c *d wimsr.1691 cnorm=cnorm+a(icspc+k-1)*a(loca-1)*a(loca) *d wimsr.1697 345 a(icspc+ig2lo+i-2)=a(l+lz+nl*nz*(i-1)) *ident up35 */ gaminr -- 16 apr 92 -- increase container array to allow p9 *d gaminr.81 common/storeg/a(50000) *d gaminr.95 iamax=50000 *i gaminr.1171 if (unow.gt.1.0) unow=1.0 *ident up36 */ matxsr -- 21 apr 92 -- make sure of word boundaries *i matxsr.564 iscrt=(iscrt/2)*2+1 *d matxsr.619 next=imatc+nwds+mult-1 *ident up37 */ groupr -- 22 apr 92 -- fix problem with vit-e temperature */ dependent weighting (sartori) *d groupr.1932 if (iwtt.ne.2.and.iwtt.ne.3.and.iwtt.ne.6.and.iwtt.ne.7 1 .and.iwtt.ne.10.and.iwtt.ne.11.and.iwtt.ne.12) *ident up38 */ reconr -- 15 may 92 -- fix problem with u.r. fission widths */ for jef (nordborg) *d reconr.631 if (ener.lt.el.or.ener.ge.eh) go to 150 *ident up39 */ broadr -- 15 jun 92 -- fix problem at resolved range boundary *d broadr.264 if (fact*enext.lt.thnmax) thnmax=fact*enext *d up23.4 *ident up40 */ acer -- 15 jun 92 -- be careful with ang. dist. loop *d acer.3982 go to 290 335 call skiprz(nin,-1) *i acer.463 call atend(mscr,0) */ acer -- 16 jun 92 -- be careful with partial fission reactions *d acer.582 if (mfd.ge.4.and.mfd.le.6.and.mtd.eq.19) mt19=1 *d acer.1080 *d acer.1374 if (mt19.eq.1.and.mt.eq.18) go to 260 if (mt19.eq.0.and. 1 (mt.eq.19.or.mt.eq.20.or.mt.eq.21.or.mt.eq.38)) go to 260 *d acer.1565 if (mt19.eq.1.and.mt.eq.18) go to 105 if (mt19.eq.0.and. 1 (mt.eq.19.or.mt.eq.20.or.mt.eq.21.or.mt.eq.38)) go to 105 */ acer -- 17 jun 92 -- fix conversion of real to integer *d acer.287 if (thin(4).eq.2.) npts=nint(thin(2)) *d acer.1120,1122 ithopt=nint(thin(4)) if (ithopt.eq.2) iwtt=nint(thin(1)) if (ithopt.eq.2) npts=nint(thin(2)) *ident up41 */ acer -- 13 jul 92 -- fix problem with law44 for photons *d acer.4669,up25.234 xss(nex+1)=44 xss(nex+2)=nex-dlwp+10 *d acer.5436 if (law.ne.4.and.law.ne.44) go to 515 *i acer.6397 if (law.eq.44) go to 580 *d acer.6434 c ***law 4 and law 44 */ acer -- 17 jul 92 -- fix mf6/mt51-90 problems *d acer.1618 call tab1io(nin,0,0,a(iscr),nb,nw) *i acer.1620 if (lf.eq.3) iso=1 if (lf.eq.3) a(iscr+3)=2 call tab1io(0,nout,0,a(iscr),nb,nw) *i acer.1729 go to 183 *d up25.133 183 if (mf.eq.6) call asend(0,nout) */ acer -- 22 jul 92 -- fix problem with sorting of thresholds *d acer.675 340 if (jethr.gt.1) call aordr(jethr,100,a(iethr)) */ acer -- 23 jul 92 -- fix problem with the use of moreio *d acer.3277 call moreio(0,nout,0,a(iscr+6),nbt,nw) *d acer.3282 if (zap.eq.0.) call tab2io(0,nout,0,a(iscr+6),nbt,nw) *d acer.3287 if (zap.eq.0.) call listio(0,nout,0,a(iscr+6),nbt,nw) *d acer.3293 call moreio(0,nout,0,a(iscr+6),nbt,nw) *d acer.3299 if (igm.eq.0) go to 516 mfh=16 call asend(nout,0) */ acer -- 18 aug 92 -- be more careful with conv check *d acer.1281,1282 185 if (rsum.ge.rmin.or.rlin.ge.rmin) go to 187 *i acer.1284 if (i.ge.ne.or.en.ge.1.05*es) go to 186 */ acer -- 18 aug 92 -- add positive threshold test *i acer.1092 character*52 messs *i acer.1384 thrx=(awr+1.)*(-c2h)/awr if (thrx.lt.0.) thrx=0. test=0. if (thrx.ne.0.) test=thresh/thrx write(messs,'(i6,1p3e15.7)') mt,thresh,thrx,test if (test.lt.1..and.test.ne.0.) 1 call mess('unionx','threshold error',messs) */ acer -- 18 aug 92 -- fix error in up25 *d up25.245,246 if (nout.ne.0) ly=nint(xss(tyr+i-1)) if (nout.eq.0) ly=iss(tyr+i-1) */ acer -- 18 aug 92 -- date should be right justified in field of 10 */ dater now expects a character argument *i acer.3635 character hdt*8 *d up25.24 call dater(hdt) hd=' '//hdt *i acer.6872 character hdt*8 *d acer.6905,6906 call dater(hdt) hd=' '//hdt *i acer.7265 character hdt*8 *d acer.7363,7364 call dater(hdt) hd=' '//hdt *i acer.7539 character hdt*8 *d acer.7688,7689 call dater(hdt) hd=' '//hdt *ident up42 */ groupr -- 1 aug 92 -- extend kalbach option */ to higher energies and other particles *i groupr.4410 common/proj/izap,awrp common/mfsix/zap,awp,lip,law,lct,aprime *d groupr.4499 280 iza2=nint(zap) aa=bach(izap,iza2,izat,enow,ep) *d groupr.4521 function bach(iza1,iza2,izat,e,ep) *d groupr.4523,4524 c compute the kalbach a parameter *d groupr.4543,4544 aa=mod(iza,1000) if (aa.eq.0.) then *d groupr.4546 call error('bach',strng,' ') *d groupr.4548,4550 za=iza/1000 ac=aa+mod(iza1,1000) zc=za+iza1/1000 ab=ac-mod(iza2,1000) zb=zc-iza2/1000 na=nint(aa-za) *d groupr.4553 sa=15.68*(ac-aa) 1 -28.07*((nc-zc)**2/ac-(na-za)**2/aa) 2 -18.56*(ac**0.66667-aa**0.66667) 3 +33.22*((nc-zc)**2/ac**1.33333-(na-za)**2/aa**1.33333) 4 -0.717*(zc**2/ac**.33333-za**2/aa**.33333) 5 +1.211*(zc**2/ac-za**2/aa) if (iza1.eq.1002) sa=sa-2.22 if (iza1.eq.1003) sa=sa-8.48 if (iza1.eq.2003) sa=sa-7.72 if (iza1.eq.2004) sa=sa-28.3 sb=15.68*(ac-ab) *d groupr.4559 if (iza2.eq.1002) sb=sb-2.22 if (iza2.eq.1003) sb=sb-8.48 if (iza2.eq.2003) sb=sb-7.72 if (iza2.eq.2004) sb=sb-28.3 ecm=aa*e/ac ea=ecm*1.e-6+sa eb=ep*1.e-6+sb x1=eb if (ea.gt.130.) x1=130.*eb/ea x3=eb if (ea.gt.41.) x3=41.*eb/ea fa=1. if (iza1.eq.2004) fa=0. fb=1. if (iza2.eq.1) fb=.5 if (iza2.eq.2004) fb=2. bach=0.04*x1+1.8e-6*x1**3+6.7e-7*fa*fb*x3**4 *ident up43 */ heatr -- 1 aug 92 -- extend kalbach option */ to higher energies and other particles *d heatr.2506 280 iza2=nint(zap) aa=bacha(1,iza2,izat,enow,ep) *d heatr.2528 function bacha(iza1,iza2,izat,e,ep) *d heatr.2530,2531 c compute the kalbach a parameter *d heatr.2568,2572 aa=mod(iza,1000) if (aa.eq.0.) then write(strng,'(''dominant isotope not known for '',i8)') iza call error('bacha',strng,' ') endif za=iza/1000 ac=aa+mod(iza1,1000) zc=za+iza1/1000 ab=ac-mod(iza2,1000) zb=zc-iza2/1000 na=nint(aa-za) *d heatr.2575 sa=15.68*(ac-aa) 1 -28.07*((nc-zc)**2/ac-(na-za)**2/aa) 2 -18.56*(ac**0.66667-aa**0.66667) 3 +33.22*((nc-zc)**2/ac**1.33333-(na-za)**2/aa**1.33333) 4 -0.717*(zc**2/ac**.33333-za**2/aa**.33333) 5 +1.211*(zc**2/ac-za**2/aa) if (iza1.eq.1002) sa=sa-2.22 if (iza1.eq.1003) sa=sa-8.48 if (iza1.eq.2003) sa=sa-7.72 if (iza1.eq.2004) sa=sa-28.3 sb=15.68*(ac-ab) *d heatr.2581 if (iza2.eq.1002) sb=sb-2.22 if (iza2.eq.1003) sb=sb-8.48 if (iza2.eq.2003) sb=sb-7.72 if (iza2.eq.2004) sb=sb-28.3 ecm=aa*e/ac ea=ecm*1.e-6+sa eb=ep*1.e-6+sb x1=eb if (ea.gt.130.) x1=130.*eb/ea x3=eb if (ea.gt.41.) x3=41.*eb/ea fa=1. if (iza1.eq.2004) fa=0. fb=1. if (iza2.eq.1) fb=.5 if (iza2.eq.2004) fb=2. bacha=0.04*x1+1.8e-6*x1**3+6.7e-7*fa*fb*x3**4 *ident up44 */ acer -- 1 aug 92 -- extend kalbach option */ to higher energies *d acer.3669 *d up32.36 xss(ki+4*n+nexd)=bachaa(1,1,iza,ee,ep) *d acer.6522 function bachaa(iza1,iza2,izat,e,ep) *d acer.6524,6525 c compute the kalbach a parameter *d acer.6553,6554 aa=mod(iza,1000) if (aa.eq.0.) then *d acer.6556 call error('bachaa',strng,' ') *d acer.6558,6560 za=iza/1000 ac=aa+mod(iza1,1000) zc=za+iza1/1000 ab=ac-mod(iza2,1000) zb=zc-iza2/1000 na=nint(aa-za) *d acer.6563 sa=15.68*(ac-aa) 1 -28.07*((nc-zc)**2/ac-(na-za)**2/aa) 2 -18.56*(ac**0.66667-aa**0.66667) 3 +33.22*((nc-zc)**2/ac**1.33333-(na-za)**2/aa**1.33333) 4 -0.717*(zc**2/ac**.33333-za**2/aa**.33333) 5 +1.211*(zc**2/ac-za**2/aa) if (iza1.eq.1002) sa=sa-2.22 if (iza1.eq.1003) sa=sa-8.48 if (iza1.eq.2003) sa=sa-7.72 if (iza1.eq.2004) sa=sa-28.3 sb=15.68*(ac-ab) *d acer.6569 if (iza2.eq.1002) sb=sb-2.22 if (iza2.eq.1003) sb=sb-8.48 if (iza2.eq.2003) sb=sb-7.72 if (iza2.eq.2004) sb=sb-28.3 ecm=aa*e/ac ea=ecm+sa eb=ep+sb x1=eb if (ea.gt.130.) x1=130.*eb/ea x3=eb if (ea.gt.41.) x3=41.*eb/ea fa=1. if (iza1.eq.2004) fa=0. fb=1. if (iza2.eq.1) fb=.5 if (iza2.eq.2004) fb=2. bachaa=0.04*x1+1.8e-6*x1**3+6.7e-7*fa*fb*x3**4 *ident up45 */ reconr -- 18 aug 92 -- fix duplicate e-sigma check *d reconr.1357 260 ir=ir+1 *d reconr.1377 if (srnext.eq.sr) go to 260 *ident up46 */ groupr -- 19 aug 92 -- fix some c.p.-related labels *d groupr.874 340 if (mfd.ne.32) go to 350 *i groupr.921 if (mfd.eq.8) mfd=6 if (mfd.eq.17.or.mfd.eq.18) mfd=16 if (mfd.eq.31) mfd=21 if (mfd.eq.32) mfd=22 if (mfd.eq.33) mfd=23 if (mfd.eq.34) mfd=24 if (mfd.eq.35) mfd=25 if (mfd.eq.36) mfd=26 *d groupr.964 data nreac/60/, npart/11/, nproj/7/ *d groupr.968 do 120 i=1,nproj *d groupr.1064 do 140 i=1,npart *ident up47 */acer -- 21 aug 92 -- fix type3 with no detailed g.p. *d acer.5786 i1(20)=nxs(1)+41 if (i1(17).eq.100) i1(20)=i1(20)+20 *d acer.5810 137 if (i1(17).ne.100) go to 141 l=nxs(1)+1 *d acer.5838 141 n=2*nxs(3)+1 *d acer.5848 n=nxs(1)-2*nxs(3) if (i1(17).eq.100) n=n+20 *d acer.8079 n=nxs(1)-2*nxs(3) if (iabs(mgpt).eq.100) n=n+20 *d acer.8085 164 if (iabs(mgpt).ne.100) go to 190 *d acer.8100 190 jxs(22)=nxs(1) *d acer.8113 210 if (nohk.eq.1) go to 220 */ acer -- 21 aug 92 -- fix type3 comment field *d acer.8065 write (hk,'(7a8)') (i1(23+i),i=1,7) */ acer -- 24 sep 92 -- fix typo affecting law=4 distributions *d acer.6103 270 if (nout.ne.1) nr=nint(xss(l)) */ acer -- 7 oct 92 -- fix problem with discrete photon law *d acer.4515,4516 ef=xss(esz+ie) el=xss(esz+ie+n-1) */ acer -- 7 oct 92 -- discontinuities were not being handled */ correctly for photon sections with nk=1 *i acer.3011 c copy total yield or xs to a scratch file *i acer.3027 c supplement grid with discontinuities *i acer.3043 c compute corrected total on the new grid *d acer.3082,3084 c handle the case nk=1 for mf12 or mf13 180 jtot=itot call tab1io(nin,nout,nscr,a(jtot),nb,nw) np=n2h nr=n1h nwt=6+2*nr+2*np if (nwt.gt.nwtot) 1 call error('convr','storage exceeded for photon data',' ') *d acer.3086,3087 jtot=jtot+nw call moreio(nin,nout,nscr,a(jtot),nb,nw) *d acer.3090 if (mf.ne.13) go to 185 a(iethr+nethr)=a(itot+6+2*nr) *d acer.3092 c check for discontinuities 185 l=itot+5+2*nr do 186 i=1,np if (i.eq.1.or.a(l+2*i-1).ne.a(l+2*i-3)) go to 186 nedis=nedis+1 if (nedis.gt.nned) 1 call error('convr','storage exceeded for edis',' ') a(idisc+nedis-1)=a(l+2*i-3) 186 continue call tosend(nin,nout,nscr,a(iscr)) */ acer -- 7 oct 92 -- move all thinning to unionx. */ don't remove discontinuities or thresholds. */ make sure there is a sharp break at */ discontinuities. *i acer.1123 iskp=0 if (ithopt.eq.1) iskp=nint(thin(3)) if (ithopt.eq.1) e1=thin(1) if (ithopt.eq.1) e2=thin(2) *d acer.1144,1146 c c ***watch for discontinuities and thresholds iedis=1 egl=1.e10 egh=1.e10 if (iedis.le.nedis) then egd=a(idisc-1+iedis) egl=sigfig(egd,ndigit,-1) egh=sigfig(egd,ndigit,+1) endif jethr=1 eet=0. if (jethr.le.nethr) eet=a(iethr-1+jethr) edl=0. edh=0. iskip=0 *d acer.1165,1171 enxt=0. i=0 120 i=i+1 e=enext *d acer.1173,1187 if (e.ge.eet) then i=iskp jethr=jethr+1 eet=0. if (jethr.le.nethr) eet=a(iethr-1+jethr) endif if (e.eq.egl) then enext=egh i=iskp go to 150 endif if (e.eq.egh) then i=iskp egl=1.e10 egh=1.e10 iedis=iedis+1 if (iedis.le.nedis) then egd=a(idisc-1+iedis) egl=sigfig(egd,ndigit,-1) egh=sigfig(egd,ndigit,+1) endif go to 140 endif if (e.eq.edl) then enext=edh i=iskp go to 150 endif if (e.eq.edh) then i=iskp edl=0. edh=0. go to 140 endif if (enext.ge.egl) then enext=egl go to 150 endif 140 if (idisc.eq.0) go to 150 edl=sigfig(enext,ndigit,-1) edh=sigfig(enext,ndigit,+1) enext=edl 150 c(1)=e *d acer.1190,1192 *i acer.1197 if (iskp.gt.0.and.e.ge.e1.and.e.le.e2.and.i.lt.iskp) go to 120 i=0 j=j+1 jt=j if (enext.ge.1.e10) jt=-jt *d acer.1199,1201 *i acer.1232 jethr=1 ee=0. if (jethr.le.nethr) ee=a(iethr-1+jethr) iedis=1 egl=1.e10 egh=1.e10 if (iedis.le.nedis) then eg=a(idis-1+iedis) egl=sigfig(eg,ndigit,-1) egh=sigfig(eg,ndigit,+1) endif *d acer.1244,1247 *d acer.1265 if (i.ne.nsave*ne/10) go to 182 *i acer.1273 c check for gamma discontinuities 182 if (ee.eq.egl) go to 186 if (ee.eq.egh) then iedis=jethr+1 egl=1.e10 egh=1.e10 if (iedis.le.nedis) then eg=a(idisc-1+iedis) egl=sigfig(eg,ndigit,-1) egh=sigfig(eg,ndigit,+1) endif go to 186 endif *i acer.1337 e=c(1) *d acer.1497,1498 *d acer.2671,2672 *d acer.2686,2707 c ***use the energy grid from unionx isave=iold iold=inew inew=isave j=0 100 j=j+1 call finda(j,c,2,inew,a(ibufn),nbuf) c(2)=0. jt=j if (c(1).ge..99999*elast) jt=-j call loada(jt,c,2,iold,a(ibufo),nbuf) if (jt.gt.0) go to 100 *i acer.2712 mfd=3 *d acer.3657,3663 *d acer.3781,3785 *d acer.3789,3798 k=k+1 *d acer.3801 */ acer -- 9 oct 92 -- provide messages for kalbach conversion. *i acer.3635 common/mainio/nsysi,nsyso,nsyse,ntty *i up32.41 write(nsyso,'(/'' converting legendre to kalbach:'', 1 '' mt ='',i4)') mth *d up32.75 if (delfcm.lt.0.) then write(strng,'('' for mt='',i3,'' e='',1p,e12.4, 1 '' eprime='',e12.4)') ee,ep call mess('acelod', 1 'converted cm distribution unreasonable',strng) delfcm=0. endif *i up32.172 s2a=sinh(2.*akal) *ident up48 */ njoy -- 24 nov 92 -- revise driver to avoid hollerith *d njoy.86,87 c * mopt six character module name, e.g. reconr. the * *d njoy.110,132 dimension z(2) character*6 module data blank/4h / *d njoy.165,194 *if sw nz=2 z(2)=blank call free(nsysi,z,nz,6) write(module,'(a4,a2)') z(1),z(2) *else nz=1 call free(nsysi,z,nz,6) write(module,'(a6)') z(1) *endif if (module.eq.'stop') then go to 1000 *d njoy.198,200 else if (module.eq.'reconr') then call reconr *d njoy.204,206 else if (module.eq.'broadr') then call broadr *d njoy.210,212 else if (module.eq.'unresr') then call unresr *d njoy.216,218 else if (module.eq.'heatr') then call heatr *d njoy.222,224 else if (module.eq.'thermr') then call thermr *d njoy.228,230 else if (module.eq.'groupr') then call groupr *d njoy.234,236 else if (module.eq.'gaminr') then call gaminr *d njoy.240,242 else if (module.eq.'errorr') then call errorr *d njoy.246,248 else if (module.eq.'covr') then call covr *d njoy.252,254 else if (module.eq.'moder') then call moder *d njoy.258,260 else if (module.eq.'dtfr') then call dtfr *d njoy.264,266 else if (module.eq.'ccccr') then call ccccr *d njoy.270,272 else if (module.eq.'matxsr') then call matxsr *d njoy.276,up12.4 else if (module.eq.'resxsr') then call resxsr *d njoy.281,283 else if (module.eq.'acer') then call acer *d njoy.287,289 else if (module.eq.'powr') then call powr *d njoy.293,295 else if (module.eq.'wimsr') then call wimsr *d njoy.299,301 else if (module.eq.'plotr') then call plotr *d njoy.305,307 else if (module.eq.'mixr') then call mixr *d njoy.311,312 else if (module.eq.'purr') then call purr c c ***illegal module name else write(nsyso,'(/'' ***error in njoy***illegal module name---'', 1 a)') module if (iopt.eq.0) go to 1000 write(ntty,'(/'' illegal module name ---'',a)') module write(ntty,'('' try again.'')') endif */ njoy -- 24 nov 92 -- add save statements *i njoy.2369 save nwtot,nr,np,lt,ip1,ip2,ip,ir,nb,nbt,xlast,ylast *i njoy.2467 save nwtot,nr,np,lt,ip1,ip2,ip,ir,nb,nbt,xlast,ylast *ident up49 */ broadr -- 24 nov 92 -- remove some cft77 caution messages *i broadr.291 tempin=temp1 *d broadr.303 npage=int((namax-50-nw-ib)/(4.*(nreac+1))) *ident up50 */ unresr -- 24 nov 92 -- remove some cft77 caution messages */ and some things found by cflint *d unresr.68 if (nin.lt.0.and.nout.gt.0) call error('unresr', *d unresr.111 write(nsyso,25) (temp(i),i=1,ntemp) *d unresr.233 ncds=2+(nint(a(ib+10))-1)/6 *d unresr.235 *d unresr.365 data wide/3.0/ *d unresr.597 nwds=inow-iarry+1 *d unresr.1010,1011 *d unresr.1026,1027 *ident up51 */ heatr -- 24 nov 92 -- remove some cft77 caution messages *i heatr.1039 save rtm,zx,ax,denom,z,aw1fac *d heatr.1047,1048 if (mtd.eq.102) then tm=939.512e6*(awr+1.) rtm=1./tm endif *i heatr.1111 save imiss,thresh,awp,afact,arat save en,cn,el,cl,z,daml,damn,enx *d heatr.1125 *d heatr.1130,1131 *i heatr.1132 mfd=4 nld=20 lcd=2 *i heatr.1182 mfd=4 *i heatr.1183 lcd=2 *i heatr.1229 save rel,fl *i heatr.1257 save nktot,awr,awfac,aw1fac,z,za,mtt save iraw,nbt,nnt,nne,ne,int save elo,ehi,flo,fhi,dlo,dhi *i heatr.1489 save lnu,ip,ir *i heatr.1814 save iraw,nne,ne,law,lang,lep,int save elo,ehi,flo,fhi,dlo,dhi *i heatr.2281 save small,ndnow,ncnow,xc,elmax *i heatr.2416 save efirst,nl,inow,mnow,lnow,ncnow,izat *i heatr.2594 save cn,ex,eimax *i heatr.2639 save iso,iraw,ir,nne,ne,int save elo,ehi,nlo,nhi *i heatr.3422 save jstart,jend,istart,ilo,ihi,nnt,ne,nne,int save elo,ehi,flo,fhi,glo,ghi,hlo,hhi *i heatr.3605 save rtm */ heatr -- 11 mar 93 -- fix up some problems found by cflint *d heatr.85 *d heatr.380 do 140 i=1,nw,6 *d heatr.581 nw=npage+50 *d heatr.1241 df=dam *d heatr.1292 *d heatr.1384,1385 *d heatr.1488 *d heatr.1490 *d heatr.1508,1509 *d heatr.1798 fl=f *d heatr.1812 *d heatr.2056 *d heatr.2296 *d heatr.2435 *d heatr.2468 *d heatr.2708 if (nne.eq.ne) *d heatr.3216 call findf(matd,mfd,mtd,nin) *d heatr.3423 *d heatr.3508 call tabbar(flo,a(istart),lf) *d heatr.3573 *d heatr.3628 dimension c(30),ncds(30) *d heatr.3659 *ident up52 */ thermr -- 24 nov 92 -- remove some cft77 caution messages *i thermr.697 save k,recon,scon,rndoff,rndup */ thermr -- 11 mar 93 -- fix up problems found by cflint *d thermr.844 *d thermr.1235 *d thermr.1353 *d thermr.1846 *d thermr.1850 *ident up53 */ groupr -- 24 nov 92 -- remove some cft77 caution messages *i groupr.1924 save ip,ir,ipl,step *i groupr.2918 save ip *i groupr.2976 save lnu,ip,ir *i groupr.3063 save lfs,mt *i groupr.3361 save ifirst *i groupr.3785 save nne,ne,int save jlo,nlo,ndlo,nclo,elo save jhi,nhi,ndhi,nchi,ehi save pspmax,langn,lepn *i groupr.4268 save small,xc,ndnow,ncnow,elmax *d groupr.4282 *i groupr.4410 save ifirst,nl,inow,lnow,mnow,ncnow *i groupr.4673 save xend save ilo,llo,mlo,nclo,xlo save ihi,lhi,mhi,nchi,xhi *d groupr.4786,4787 *i groupr.4877 save iecl,iech,ecl,ech,ecn,nqp0 *d groupr.4946 npo=ld+nl+int(alog(300/awr)) *i groupr.5143 save iso,lidp,ltt,iraw,ir,nne,ne,int save elo,ehi,nlo,nhi *i groupr.5286 save nwt,ir,ncyc,ltt,nu,ne,ie,nw1,nw2 save llo,nlo,elo save lhi,nhi,ehi save ibrag,jbrag *i groupr.5673 save li *i groupr.5996 save awr *i groupr.6749 save econ,nktot,nbt,nnt,iraw,nupm save elo,ehi,istart,jstart,nne,ne,int *i groupr.7016 save new,nr,np,ip2,ir2,loct,theta,xc,locb,ca,rc,bot *i groupr.7432 save cn,eimax,ex *ident up54 */ errorr -- 24 nov 92 -- remove some cft77 caution messages *i errorr.424 c *d errorr.940,950 *i errorr.1460 save matd,mfrd,matlst c *i errorr.3731 save ir,ip *ident up55 */ covr -- 24 nov 92 -- remove some cft77 caution messages *i covr.698 save ixpo,igprnt *d covr.1574 if (ia.eq.0) go to 110 *d covr.1761 nlmax=int(10./dylin) *ident up56 */ matxsr -- 24 nov 92 -- fix missing path *d matxsr.1511 if (mt.ne.1.and.mt.ne.501) go to 134 */ matxsr -- 11 mar 93 -- increase number of materials allowed *d matxsr.542 1jinp(10),joutp(10),hmatn(250),matno(250),matgg(250) *d matxsr.722 1jinp(10),joutp(10),hmatn(250),matno(250),matgg(250) *d matxsr.892 1jinp(10),joutp(10),hmatn(250),matno(250),matgg(250) */ matxsr -- 11 mar 93 -- fix some problems found by cflint *d matxsr.505 *d matxsr.535 *d matxsr.662 *d matxsr.723 *d matxsr.1086 ng2=l1h *d matxsr.907,908 data hgsct/6hgscat /, hnthr/6hntherm/ *d matxsr.1237 if (math.eq.0.and.mfh.eq.0.and.mth.eq.0) go to 110 *d matxsr.1451,1453 *d matxsr.1674 *d matxsr.1691 nfg=ning+1 *d matxsr.1722 nwds=0 *d matxsr.1802 *d matxsr.1828 *d matxsr.1919 *d matxsr.1980 *d matxsr.2018 *d matxsr.2107 *d matxsr.2719,2720 *d matxsr.2728 if (nin.lt.0) go to 140 *d matxsr.2843 *ident up57 */ moder -- 28 nov 92 -- don't compare */ real variables containing hollerith *d moder.46 dimension ia(1) dimension iz(20) equivalence (a(1),ia(1)) equivalence (z(1),iz(1)) data iblank/4h / *d moder.77 120 iz(1)=iblank *d moder.118 call tpidio(nin,no,nscr,ia,nb,nw) *d moder.121 if (ia(1).ne.iblank) go to 160 */ moder -- 28 nov 92 -- remove previous change that */ damages gendf fission matrices *d up20.4 730 if (mf.ne.1.and.ig.ne.ng) go to 710 *ident up58 */ plotr -- 28 nov 92 -- fix error and cft77 cautions *d up3.128 51 format(/12h tag ...... ,a/ *i plotr.331 i3d=0 *d plotr.502,505 *i plotr.512 if (iverf.eq.0) go to 260 if (nin.ne.ninl) call closz(ninl) if (nin.ne.ninl) call openz(nin,0) ninl=nin *i plotr.796 ltt=0 *i plotr.808 emin=0. emax=0. lang=0 *i plotr.1066 nl=0 nz=0 *i plotr.1574 save y *i plotr.1621 j=0 *d plotr.1853 n1=int(zmin/dz) *d plotr.1861 n2=int(zmax/dz) *i plotr.1896 save nwtot,nep,ncyc,lt,ip1,ip2,ip *d plotr.1924 ln=ncyc*(ip-ip1)+lt *ident up59 */ njoy -- 25 nov 92 -- make free recognize uppercase e */ and ignore leading zero in exponents *i njoy.448 integer ee *i njoy.451 data ee/1hE/ *d njoy.607 340 if (char.ne.e.and.char.ne.ee.and. 1 char.ne.plus.and.char.ne.minus) go to 350 *d njoy.610 if (char.eq.e.or.char.eq.ee) go to 120 *d njoy.631,632 410 if (expo) go to 415 z(j)=z(j)+(n-1)*fact *i njoy.635 415 if (power.eq.0..and.n.eq.1) go to 120 power=power+(n-1)*fact fact=fact/10 ipow=ipow+1 go to 120 */ njoy -- 11 mar 93 -- fix up problems found by cflint *d up48.8 *d njoy.800 *d njoy.810 inin=iabs(nin) *d njoy.896 inin=iabs(nin) *d njoy.986 inin=iabs(nin) *d njoy.1085 *d njoy.1092 inin=iabs(nin) *d njoy.1169 inin=iabs(nin) *d njoy.1242 *d njoy.1248 inin=iabs(nin) *d njoy.1292 inin=iabs(nin) *d njoy.1362 inin=iabs(nin) */ njoy -- 1 apr 93 -- fix pointer problem in scana *d njoy.2107 *ident up60 */ groupr -- 21 dec 92 -- fix problem with hnab that leads */ to errors for e-dep watt spectrum. *i groupr.7317 if (dabs(aa-bb).ge.1.0) go to 110 *i groupr.7320 if (abs(aa-bb).ge.1.0) go to 110 *i groupr.7327 if (dabs(aa-bb).ge.1.0) go to 320 *i groupr.7330 if (abs(aa-bb).ge.1.0) go to 320 *d groupr.7337 *d groupr.7339 a=dabs(bb) b=dabs(aa) *d groupr.7341 a=abs(bb) b=abs(aa) *d groupr.7345 205 continue *d groupr.7354 h=(b-a)*pow2(1) *i groupr.7356 asq=a*a *d groupr.7358 if (asq.gt.120.) asq=120. con=dexp(-asq)*resqpi/pow2(n+1) *d groupr.7360 con=exp(-asq)*resqpi/pow2(n+1) */ groupr 21 dec 92 -- remove some nuisance messages *d groupr.6192 314 if (za2.lt.0..and. 1 mth.ne.18.and.mth.ne.19.and.mth.ne.20.and. 2 mth.ne.21.and.mth.ne.38) then *ident up61 */ broadr -- 21 dec 92 -- fix problem in hnabb for sw mach. *d broadr.1235,1236 a=dabs(bb) b=dabs(aa) */ broadr -- 11 mar 93 -- fix up problems found by cflint *d broadr.480 *d broadr.594 *d broadr.597 *d broadr.911 *d broadr.968 x=0.0 */ broadr -- 15 mar 93 -- fix problem with restart *d broadr.201,204 if (iverf.ge.5) call contio(nin,no,nscr1,a(iscr),nb,nw) if (iverf.ge.6) call contio(nin,no,nscr1,a(iscr),nb,nw) call hdatio(nin,no,nscr1,a(iscr),nb,nw) temp=a(iscr) write(nsyso,65) temp *d broadr.214,216 145 call tomend(nin,no,nscr1,a(iscr)) *ident up62 */ groupr -- 22 dec 92 -- fix groupr to handle multiple */ emissions in file 6 (f-19). */ we'll need more storage. *d groupr.233 common/gstore/a(100000) *d groupr.250 iamax=100000 *d groupr.353 if (nsigz.gt.1) then nw=nflmax*(3+nsigz) if (nw.eq.0) nw=1 call reserv('flxc',nw,iflxc,a) call genflx(a,a(iflxc)) call releas('flxc',0,a) endif *d groupr.2041 subroutine genflx(a,b) *d groupr.2073 *d groupr.2090 dimension b(1) *d groupr.2128 *i groupr.3783 dimension iyss(3),izss(3),jjss(3),jloss(500) *i up53.17 save idis,iyss,izss,jjss,jloss data ncmax/350/ *i groupr.3796 idis=0 *i groupr.3817 nss=1 jjss(1)=1 iyss(1)=iy elo=1.e10 *d groupr.3885,3886 140 if (law.ge.2.and.law.le.5) go to 194 *d groupr.3889,3972 c ***read in data for law 1 if (law.ne.1) go to 160 izss(nss)=l iz=l call tab2io(nin,0,0,c(l),nb,nw) lang=nint(c(l+2)) lep=nint(c(l+3)) ne=nint(c(l+5)) nbt=nint(c(l+6)) int=nint(c(l+7)) nn=0 l=l+nw nnn=0 do 155 ie=1,ne ilo=l jlo=ilo ii=jjss(nss)+ie-1 if (ii.gt.500) call error('getmf6', 1 'too many subsection energy points',' ') jloss(ii)=jlo call listio(nin,0,0,c(l),nb,nw) if (c(l+1).lt.elo) elo=c(l+1) nn=nint(c(l+5)) if (ie.gt.1.and.nn.ne.nnn.and.int.ge.11.and.int.le.20) then int=int+10 c(iz+7)=int call mess('getmf6', 1 'bad grids for corresponding-point interpolation-', 2 'changing to unit-base interpolation') endif nnn=nn 145 l=l+nw if (l.gt.iy+nc) call error('getmf6','storage exceeded',' ') if (nb.eq.0) go to 147 call moreio(nin,0,0,c(l),nb,nw) go to 145 147 if (lct.ne.2) go to 155 call cm2lab(ilo,jlo,l,c,nl,lang,lep,max) do 150 i=jlo,l 150 c(ilo+i-jlo)=c(i) l=ilo+l-jlo 155 continue go to 190 c c ***read in data for law 6. convert them to law 1. 160 if (law.ne.6) go to 170 call contio(nin,0,0,c(l),nb,nw) apsx=c(l) npsx=c(l+5) c(l)=0. c(l+1)=0. c(l+2)=0 c(l+3)=0 lang=0 lep=2 lct=2 irr=nint(c(iy+4)) npp=nint(c(iy+5)) enow=c(iy+6+2*irr)/1.2 pspmax=c(iy+4+2*irr+2*npp) izss(nss)=l l=l+8 ne=0 162 ne=ne+1 ii=jjss(nss)+ne-1 if (ii.gt.500) call error('getmf6', 1 'too many subsection energy points',' ') ilo=l jloss(ii)=ilo c(ilo)=apsx enow=enow*1.2 if (enow.gt.pspmax) enow=pspmax c(ilo+1)=enow if (enow.lt.elo) elo=enow c(ilo+2)=0 c(ilo+3)=0 c(ilo+4)=0 c(ilo+5)=npsx l=l+6 call cm2lab(ilo,jlo,l,c,nl,lang,lep,max) do 163 i=jlo,l 163 c(ilo+i-jlo)=c(i) l=ilo+l-jlo if (enow.lt.pspmax) go to 162 l=izss(nss) c(l+2)=lang c(l+3)=lep c(l+4)=1 c(l+5)=ne c(l+6)=ne c(l+7)=22 law=1 c(iy+3)=law go to 190 c c ***read in data for law 7 170 izss(nss)=l call tab2io(nin,0,0,c(l),nb,nw) lang=nint(c(l+2)) lep=nint(c(l+3)) ne=nint(c(l+5)) l=l+nw do 185 ie=1,ne ilo=l jlo=ilo ii=jjss(nss)+ie-1 if (ii.gt.500) call error('getmf6', 1 'too many subsection energy points',' ') jloss(ii)=ilo call tab2io(nin,0,0,c(l),nb,nw) nmu=nint(c(l+5)) l=l+nw l1=l do 177 i=1,nmu call tab1io(nin,0,0,c(l),nb,nw) if (c(l+1).lt.elo) elo=c(l+1) 172 l=l+nw if (l.gt.iy+nc) call error('getmf6','storage exceeded',' ') if (nb.eq.0) go to 175 call moreio(nin,0,0,c(l),nb,nw) go to 172 175 continue 177 continue call ll2lab(ilo,jlo,l,c,nl,max) do 180 i=jlo,l 180 c(ilo+i-jlo)=c(i) l=ilo+l-jlo 185 continue c c ***check for multiple particle subsections 190 if (ik.eq.nk) go to 195 ik=ik+1 call tab1io(nin,0,0,c(l),nb,nw) l1=l 192 l=l+nw if (l.gt.iy+nc) call error('getmf6','storage exceeded',' ') if (nb.eq.0) go to 193 call moreio(nin,0,0,c(l),nb,nw) go to 192 193 if (c(l1).ne.zap) go to 195 law=c(l1+3) nss=nss+1 if (nss.gt.3) call error('getmf6', 1 'too many subsections for one particle',' ') iyss(nss)=l1 jjss(nss)=jjss(nss-1)+ne go to 140 c c ***discrete two-body scattering 194 izss(nss)=l call getdis(ed,elo,idisc,yldd,ans,nl,ng2,iglo,nq,nin,c(l)) l=l+ncmax c c ***initialization complete 195 continue if (nss.gt.1) call mess('getmf6', 1 'there are multiple subsections in mf6', 2 'for this emitted particle') nc=l-ic call releas('sc',nc,c) enext=elo idisc=0 return *i groupr.3975 idisc=0 *d groupr.3979,3980 200 do 205 j=1,ng *i groupr.3984 yld=0. iss=1 207 iy=iyss(iss) iz=izss(iss) lip=nint(c(iy+2)) law=nint(c(iy+3)) if (law.ge.2.and.law.le.5) go to 500 *d groupr.3994 yld=yld+pe lang=nint(c(iz+2)) lep=nint(c(iz+3)) lepn=lep if (law.eq.1.and.lct.eq.2) lepn=2 if (law.eq.7) lepn=2 langn=lang if (law.eq.1.and.lct.eq.2) langn=1 if (law.eq.7) langn=1 *d groupr.4006,4092 c ***find desired energy panel nne=nint(c(iz+5)) nbt=nint(c(iz+6)) int=nint(c(iz+7)) jj=jjss(iss) do 230 ie=1,nne jlo=jloss(jj+ie-1) jhi=jloss(jj+ie) if (ed.ge.c(jlo+1).and.ed.le.c(jhi+1)) go to 300 230 continue go to 450 c c ***set up loop over secondary energy *d groupr.4107 ans(l,jg)=ans(l,jg)+pe*0.5*(en-ep)*(term(l)+terml(l)) *i groupr.4114 ndlo=nint(c(jlo+2)) nplo=nint(c(jlo+5)) nclo=nint(c(jlo+4))/nplo nphi=nint(c(jhi+5)) nchi=nint(c(jhi+4))/nphi *d groupr.4131 1 ans(1,j)=ans(1,j)+pe*g0 *i groupr.4140 c c ***loop back to process additional particles, if any. iss=iss+1 if (iss.le.nss) go to 207 *d groupr.4144 500 ilo=iz *d groupr.4151 if (yld.eq.0.) return *d groupr.4161 if (abs(test-yld).gt..02 *i groupr.6628 if (mth.eq.iabs(mf18(imf18-1))) go to 790 *i groupr.6638 if (mth.eq.iabs(mf6(imf6-1))) go to 790 *i groupr.6648 if (mth.eq.iabs(mf6p(1,imf21-1))) go to 790 *i groupr.6658 if (mth.eq.iabs(mf6p(2,imf22-1))) go to 790 *i groupr.6668 if (mth.eq.iabs(mf6p(3,imf23-1))) go to 790 *i groupr.6678 if (mth.eq.iabs(mf6p(4,imf24-1))) go to 790 *i groupr.6688 if (mth.eq.iabs(mf6p(5,imf25-1))) go to 790 *i groupr.6698 if (mth.eq.iabs(mf6p(6,imf26-1))) go to 790 *d groupr.7444 if (eimax.le.0.) eimax=100. *d groupr.7455 *ident up63 */ plotr -- 26 feb 93 -- add more 2-d and 3-d options for gendf data */ also fix up things found by cflint *i plotr.178 c c * special meaning for nth for gendf mf=6 data * c * nth=1 plot 2-d spectrum for group 1 * c * nth=2 plot 2-d spectrum for group 2 * c * etc. * c * no special flags are needed for mf=6 3d plots * *d plotr.287 common/setup8/ex3(400),ey3(400) *d plotr.290 dimension a(1200),z(15),aa(40000) *d plotr.304,306 data max/2500/, nwamax/1200/ data maxaa/40000/ data maxx3/400/, maxy3/400/ data eps/1.e-7/, delta/1.26/ *d plotr.310 call gplot(1hu,13ht2 njoy plotr,13) *d plotr.315 *d plotr.598 *d plotr.669 if (mmf.eq.6) call gety6(ntp,0.,enext,idis,yf,nin,a,lep,lf,b,1) *d plotr.678 if (mmf.eq.6) call gety6(ntp,enow,enext,idis,yf,nin,a,lep,lf,b,1) *d plotr.693 if (mmf.eq.6) call gety6(ntp,enow,enext,idis,yf,nin,a,lep,lf,b,1) *d plotr.714 if (mmf.eq.6) call gety6(ntp,enow,enext,idis,yf,nin,a,lep,lf,b,1) *i plotr.770 go to 2307 *d plotr.791,792 *i plotr.804 nmum=nmu+2 *d plotr.826 if (locn+3*nmum.gt.maxaa) go to 2460 if (i+3.gt.maxx3) go to 2460 *d plotr.831 do 2418 imu=1,nmum *d plotr.833 locn=locn+nmum *i plotr.835 ey3(1)=1. aa(locn+1)=alog10(base) locn=locn+1 *i plotr.837 ey3(1+imu)=amu *d plotr.855,857 2449 if (sum.lt.base) sum=base *i plotr.863 ey3(nmum)=-1. aa(locn+1)=alog10(base) locn=locn+1 *d plotr.866 do 2452 imu=1,nmum *d plotr.868 locn=locn+nmum *d up3.12 if (locn+3*nmum.gt.maxaa.or.i+3.gt.maxx3) 1 call mess('plotr', *d plotr.899 call set3d(iplot,aa,nmum,nmax,a) *d up3.5 2525 ne2=198 ne2m=ne2+2 *d up3.27 if (locn+3*ne2m.gt.maxaa) go to 2545 if (il+3.gt.maxx3) go to 2545 *d plotr.946 do 2528 i2=1,ne2m *d plotr.948 locn=locn+ne2m *i plotr.950 if (y3.gt.0.) ey3(1)=ymin if (y3.lt.0.) ey3(1)=ymax aa(locn+1)=-12. locn=locn+1 *i plotr.952 if (y3.gt.0.) ey3(1+i2)=e2 if (y3.lt.0.) ey3(1+ne2-i2+1)=e2 *i plotr.964 if (y3.gt.0.) ey3(ne2m)=ymax if (y3.lt.0.) ey3(ne2m)=ymin aa(locn+1)=-12. locn=locn+1 *d plotr.967 do 2542 i2=1,ne2m *d plotr.969 locn=locn+ne2m *d up3.35 if (locn+3*ne2m.gt.maxaa.or.il+3.gt.maxx3) 1 call mess('plotr', *d plotr.982,983 do 2560 i2=1,ne2m l=i2+ne2m*(i1-1) *d plotr.1008 if (ystp.ne.0.) go to 2580 *d plotr.1023 call set3d(iplot,aa,ne2m,nmax,a) *d plotr.1051,1052 450 continue *d plotr.1057 *d plotr.1065,1066 jbase=jbase+1 if (i3d.eq.1) go to 3400 *i plotr.1068 ig1=0 ig=ngn *d plotr.1090,1091 if (mth.ne.mtd) go to 455 if (mfd.eq.6.and.mfh.eq.3) go to 495 if (mfh.ne.mfd) go to 455 if (mfh.eq.6) go to 490 if (mfh.ne.3) go to 455 *i plotr.1104 490 if (ig.ne.nth) go to 455 do 492 j=2,ng2 jj=ig2lo+j-2 jnow=2*(j-2) x(jnow+1)=a(locngn+jj-1)*factx x(jnow+2)=a(locngn+jj)*factx s=a(jbase+6+nl*nz*(j-1)) dener=a(locngn+jj)-a(locngn+jj-1) s=s/dener s=s/sigig y(jnow+1)=s*facty y(jnow+2)=s*facty if (s.ne.0.) n=jnow+2 492 continue go to 455 495 if (ig.ne.nth) go to 455 sigig=a(jbase+6+nl*nz) go to 455 c c ***gendf 3d plots 3400 call contio(nin,0,0,a(jbase),nb,nw) if (math.le.0) then call mess('plotr', 1 'desired gendf matrix not found',' ') go to 110 end if if (mfh.eq.mfd.and.mth.eq.mtd) go to 3410 call tosend(nin,0,0,a(jbase)) go to 3400 3410 ng1=0 nl=nint(a(jbase+2)) nz=nint(a(jbase+3)) ngl=nint(a(jbase+5)) xmin=xleft xmax=xright if (xmax.eq.0.) xmax=a(locngn+ngl) if (xmin.eq.0.) xmin=a(locngn) xstp=xstep emin=0. emax=1.e10 ymin=ybot ymax=ytop if (mfd.eq.16) then if (ymax.eq.0.) ymax=a(locngg+ngg) if (ymin.eq.0.) ymin=a(locngg) else if (ymax.eq.0.) ymax=a(locngn+ngn) if (ymin.eq.0.) ymin=a(locngn) endif ystp=ystep if (itype.ne.2.and.itype.ne.4) go to 3414 3413 call algplt(ymin,ymax,abs(y3),yor,ycy) if (ycy.gt.0.42) go to 3414 ymin=ymin*10. go to 3413 3414 locn=0 zmax=-12. 3420 ng1=ng1+1 call listio(nin,0,0,a(jbase),nb,nw) l=jbase 3419 l=l+nw if (nb.eq.0) go to 3418 call moreio(nin,0,0,a(l),nb,nw) go to 3419 3418 ng2=nint(a(jbase+2)) ig2lo=nint(a(jbase+3)) ig=nint(a(jbase+5)) if (ng1+2.gt.maxx3) call error('plotr', 1 'too many incident groups for 3d gendf plot',' ') ex3(ng1)=a(locngn+ig-1) if (ng1.eq.1) emin=ex3(ng1) ngm=ngn+1 if (mfd.eq.16) ngm=ngg+1 if (locn+2*ngm+4.gt.maxaa) call error('plotr', 1 'too many data for 3d gendf plot',' ') ey3(1)=ymin ey3(2)=ey3(1) j2=2 do 3421 i2=1,ngm e2=a(locngn+i2-1) if (mfd.eq.16) e2=a(locngg+i2-1) if (e2.le.ymin) go to 3421 if (e2.ge.ymax) go to 3421 if (j2+4.gt.maxy3) call error('plotr', 1 'too many secondary groups for 3d gendf plot',' ') ey3(j2+1)=e2 ey3(j2+2)=e2 j2=j2+2 3421 continue ey3(j2+1)=ymax ey3(j2+2)=ymax j2=j2+2 do 3422 i2=1,j2 aa(locn+i2)=-12. 3422 aa(locn+j2+i2)=-12. locn=locn+j2 ng1=ng1+1 ex3(ng1)=a(locngn+ig-1) k2=1 do 3450 j=2,ng2 ig2=ig2lo+j-2 eglo=a(locngn+ig2-1) if (mfd.eq.16) eglo=a(locngg+ig2-1) eghi=a(locngn+ig2) if (mfd.eq.16) eghi=a(locngg+ig2) if (eghi.le.ymin) go to 3450 if (eglo.gt.ymax) go to 3450 f2=a(jbase+6+nl*nz*(j-1))/(eghi-eglo) if (f2.lt.1.e-12) f2=1.e-12 f2=alog10(f2) if (f2.gt.zmax) zmax=f2 k2=k2+1 aa(locn+k2)=f2 k2=k2+1 aa(locn+k2)=f2 3450 continue k2=k2+1 locn=locn+j2 if (ig.lt.ngl) go to 3420 ng1=ng1+1 ex3(ng1)=a(locngn+ig) ng1=ng1+1 ex3(ng1)=1.001*a(locngn+ig) do 3431 i2=1,j2 aa(locn+i2)=aa(locn+i2-2*j2) 3431 aa(locn+i2+j2)=aa(locn+i2-j2) locn=locn+2*j2 nmax=ng1 emax=ex3(ng1) c fix up z axis i=zmax if (i.lt.zmax) i=i+1 zmax=i zmin=zmax-3.3 i=zmin if (i.gt.zmin) i=i-1 zmin=i do 3440 i1=1,nmax do 3440 i2=1,j2 l=i2+j2*(i1-1) if (aa(l).lt.zmin) aa(l)=zmin 3440 continue rbot=zmin rtop=zmax rstep=1. c fix up x axis if (itype.ne.3.and.itype.ne.4) go to 3468 if (xmin.eq.0.) xmin=1.e-5 xmin=alog10(emin) xmax=alog10(emax) xstp=1. do 3466 i1=1,nmax 3466 ex3(i1)=alog10(ex3(i1)) 3468 if (xstp.ne.0.) go to 3470 xmin=emin xmax=emax call ascale(4,xmin,xmax,major,minor) xstp=(xmax-xmin)/major 3470 continue xleft=xmin xright=xmax xstep=xstp c fix up y axis if (itype.ne.2.and.itype.ne.4) go to 3478 ymin=alog10(ymin) ymax=alog10(ymax) ystp=1. do 3477 i2=1,j2 3477 ey3(i2)=alog10(ey3(i2)) 3478 if (ystp.ne.0.) go to 3480 call ascale(2,ymin,ymax,major,minor) ystp=(ymax-ymin)/major 3480 continue ybot=ymin ytop=ymax ystep=ystp c set up axis labels rl=probp nr=nprobp jtype=2 yl=sece ny=nsece if (nplot.ne.0.and.iplot.eq.1) call endpl(0) if (iplot.eq.-1) call endgr(0) call set3d(iplot,aa,j2,nmax,a) nplot=1 go to 110 *d plotr.1161 *d plotr.1285 *d plotr.1289 *d plotr.1669,1670 dimension zmat(n3y,n3x),work(1) *d plotr.1706 call surtrn('both') *d plotr.1720 call ylog(0.,1.,zor,zcy) *i plotr.1760 c function x3dmat(i) c ****************************************************************** c function for rescaling the x axis for 3d plots c ****************************************************************** common/setup8/ex3(400),ey3(400) x3dmat=ey3(i) return end *d plotr.1763 c ****************************************************************** *d plotr.1765,1767 c ****************************************************************** common/setup8/ex3(400),ey3(400) y3dmat=ex3(i) */ plotr -- 15 mar 93 -- move frame in window slightly *d plotr.1308,1309 xpos=xll+4.7*hlab ypos=yll+4.7*hlab *ident up64 */ groupr -- 1 mar 93 -- fix error in up30 */ that causes bad mt50+n gamma yields *d up30.10 tsave=a(iy+i-1) */ groupr -- 1 mar 93 -- change e interpolation in getaed */ this one tracks peaks better (zrh) *d groupr.5287 *d up53.37,39 save l2,l3,llo,nlo,elo,lhi,nhi,ehi save ibrag,jbrag,nbrag *d groupr.5397,5401 235 if (loc-iaes.gt.nwmax) *d groupr.5403 nw2=loc-l3 *d groupr.5322,5324 *d groupr.5327 lhi=llo nw1=loc-l2 *d groupr.5425,5432 egp1=0. egp2=0. *d groupr.5439,5448 egp1=egp-e+elo egp2=egp+ehi-e if (egp1.le.eg1) then egp1=egp egp2=egp end if *d groupr.5604,5662 */ groupr -- 1 mar 93 -- remove unused variables *d groupr.2460 data naed/14/ *d groupr.2977,2978 *d groupr.2990 *d groupr.2992 *d groupr.3200 *d groupr.3202 *d groupr.3275 *d groupr.3786,3793 *d groupr.3797 *d up53.22 save efirst,nl,inow,lnow,mnow,ncnow *d groupr.4428 *d groupr.4461 *d groupr.5305 *d groupr.5408 *d groupr.5414,5417 *d groupr.5434 *d groupr.5682 *d groupr.5917 *d groupr.5919 *d groupr.5922 *d groupr.6559,6565 *d groupr.6802,6803 *d groupr.6857 *d groupr.6866 *d groupr.7076 */ groupr -- 3 mar 93 -- fix getsed to work for all combinations */ of analytic and tabulated subsections *d groupr.6742,6743 *d up53.45,up53.46 save econ,nktot,nupm,loc *d groupr.6751 data nkmax/20/ *d groupr.6763,6764 *d groupr.6768 enext=1.e10 *i groupr.6769 if (ik.gt.nk) go to 190 *d groupr.6777,6778 lf=nint(c(il+ic+3)) nr=nint(c(il+ic+4)) *d groupr.6780 if (elo.lt.enext) enext=elo *d groupr.6782,6798 if (lf.eq.1) go to 162 *d groupr.6806,6807 nr=nint(c(l+4)) np=nint(c(l+5)) *d groupr.6822,6823 nr=nint(c(l+4)) np=nint(c(l+5)) *d groupr.6837,6845 160 il=l-ic go to 100 c c ***tabulated subsection 162 call tab2io(nin,0,0,c(l),nb,nw) ne=nint(c(l+5)) nnt=l+6-ic nbt=nint(c(nnt+ic)) int=nint(c(nnt+ic+1)) l=l+nw nne=0 c c ***read and average spectrum for each incident energy 165 il=l-ic l=l+ng+1 iraw=l call tab1io(nin,0,0,c(l),nb,nw) *i groupr.6846 if (l.gt.nc) call error('getsed','storage exceeded',' ') *d groupr.6850,6854 180 nne=nne+1 c(il+ic)=c(iraw+1) if (nne.eq.1) elo=c(iraw+1) if (nne.gt.1) ehi=c(iraw+1) *d groupr.6860,6862 l=il+ig+ic *d up17.29 if (ig.eq.1) e1=0. *d up17.31 if (ig.eq.ng) e2=1.e8 *d up17.32 call intega(c(l),e1,e2,c(iraw),ip,ir) *i groupr.6864 ll=l-ng-1 *i groupr.6870 l=iraw *d groupr.6874,6875 if (nne.lt.ne) go to 165 go to 100 c c ***initialization complete 190 nc=l-ic *d groupr.6877 *d groupr.6889 do 290 ik=1,nk *d groupr.6902,6943 if (pe.eq.0.) go to 290 nr=nint(c(lnow+4)) np=nint(c(lnow+5)) mnow=lnow+6+2*nr+2*np if (lf.eq.1) go to 230 c c ***analytic subsection. compute sed at e. ip=2 ir=1 do 220 ig=1,ng e1=eg(ig) if (ig.eq.1) e1=0. e2=eg(ig+1) if (ig.eq.ng) e2=1.e8 call anased(s,ed,ethi,idisc,e1,e2,c(lnow),ip,ir) s=s*pe sed(ikt,ig)=sed(ikt,ig)+s 220 continue if (ethi.eq.enext.and.idisc.gt.idis) idis=idisc if (ethi.lt.enext) idis=idisc if (ethi.lt.enext) enext=ethi go to 290 c c ***tabulated subsection. interpolate for sed at e. 230 nne=0 nr=nint(c(mnow+4)) ne=nint(c(mnow+5)) nnow=mnow+6+2*nr nnt=6 nbt=nint(c(mnow+nnt)) int=nint(c(mnow+nnt+1)) 240 nne=nne+1 if (nne.gt.ne) go to 300 if (nne.le.nbt) go to 250 nnt=nnt+2 nbt=nint(c(mnow+nnt)) int=nint(c(mnow+nnt+1)) 250 llo=nnow+(ng+1)*(nne-1) elo=c(llo) lhi=nnow+(ng+1)*nne ehi=c(lhi) if (ed.gt.ehi) go to 240 do 260 ig=1,ng call terp1(elo,c(llo+ig),ehi,c(lhi+ig),ed,s,int) sed(ikt,ig)=sed(ikt,ig)+s*pe 260 continue *d groupr.6949,6951 if (mfd.eq.15.or.mtd.eq.38) go to 280 if (mtd.gt.17.and.mtd.lt.22) go to 280 do 270 kg=1,ng *d groupr.6953 if (ed.gt.eg(ig)*1.001) go to 280 *d groupr.6954,6956 sed(ikt,ig-1)=sed(ikt,ig-1)+sed(ikt,ig) sigup=sed(ikt,ig) 270 sed(ikt,ig)=0. *d up30.43 280 if (sigup.gt..1.and.nupm.lt.10) then *d groupr.6963,6979 c c ***continue loop over subsections 290 continue if (mtd.lt.18.or.(mtd.gt.21.and.mtd.ne.38)) return econ=1.1*ed if (ed.lt.1e4) econ=1.e4 *d groupr.6986,6988 300 if (int.eq.1) go to 320 do 310 ig=1,ng 310 sed(1,ig)=0. *d groupr.6994,6996 320 do 330 ig=1,ng 330 sed(1,ig)=ehi *d groupr.7044 if (epl.gt.de.or.de.le.0.) return */ groupr -- 6 mar 93 -- fix more cft warnings and cautions */ also fix up things found by cflint *d groupr.238,239 data blank/4h /, iblank/4h / *d groupr.291 if (ngout1.eq.0.or.ngout2.eq.0) go to 170 *d groupr.472 *d groupr.640 if (matb.eq.0) go to 650 *i groupr.809 save ir *d groupr.1164,1165 *d groupr.1525,1526 *d groupr.1646,1647 do 910 ig=1,ngg 910 write(nsyso,5) ig,egg(ig),egg(ig+1) *d groupr.1910 *i groupr.2110 net2=0 *i groupr.2283 sigzj=1.e10 indx=1 *i groupr.2336 fac=0. *d groupr.2402 if (abs(tempn-temp(jtemp)).le.eps) go to 150 *d groupr.2405 if (flag.ge.0.) go to 130 *d groupr.2420 if (abs(e-en).le.eps) go to 240 *i groupr.2551 save matl,mfl,mtl,ig1w,iglow,eqw save elast,flst,slst *d groupr.2698,2699 *i groupr.2703 save irat,itr *d groupr.2894 *d groupr.3057 *d groupr.3351 *d groupr.3357 *d groupr.3359 *d groupr.3664 *d groupr.3784 *d up53.15,16 save jlo,elo,jhi,ehi *d groupr.4266 *i groupr.4308 *d groupr.4469 epnn=0. *d groupr.4511,4517 350 if (lep.gt.1) go to 370 if (epnext.eq.1.e10) go to 370 if (ep.ge.dn*dn*epnext) go to 360 epnext=dn*epnext go to 370 360 epnext=up*epnext 370 f6ddx=t *i groupr.5036 sigc=0. *d groupr.5144 *d groupr.5228 if (nne.eq.ne) *d groupr.5541 do 520 i=1,nw2 *i groupr.5881 ipm=0 *i groupr.6160 mt0=0 *d groupr.6260 if (mfh.ne.12) go to 113 *d groupr.6801 if (lf.ne.3) go to 125 *i groupr.6813 jp=0 */ groupr -- 15 mar 93 -- fix t-dependent claw *d up37.5,6 if (iwtt.eq.1.or.iwtt.eq.4.or.iwtt.eq.5.or.iwtt.eq.8. 1 or.iwtt.eq.9.or.iwtt.eq.10) *ident up65 */ acer -- 11 mar 93 -- more changes for type 3 at lanl *d acer.5747 open(nout) *d acer.5750 open(nout) *d acer.5802 if (udisk(nout)) 136,910,910 *d acer.5807 if (udisk(nout)) 137,910,910 *d acer.5842 if (udisk(nout)) 142,910,910 *d acer.5851 if (udisk(nout)) 143,910,910 *d up25.45 open(nin) *d up25.47 if (udisk(nin)) 161,250,250 *d up25.52 if (udisk(nin)) 163,250,250 *d up25.57 if (udisk(nin)) 164,250,250 */ acer -- 11 mar 93 -- fix problems found by cflint *d up25.341 data namax/20000/, nidmax/27/ *d acer.385,386 *d acer.412,413 *d acer.642 *d acer.650 *d acer.814 nen=-ne *d acer.1158 *d acer.1248 rsum=0. *d acer.1275 if (ee.eq.0..or.en.lt.ee) go to 185 *d up47.92 *d acer.1534 *d acer.1557 *d acer.1595 *d up25.16 *d acer.2158 *d acer.2221 *d acer.2356 k=k+1 *d acer.2396 k=k+1 *d acer.2619 *d acer.2802 if (nk.eq.1) go to 230 *d acer.2966 *d acer.3001 call tofend(nin,nout,nscr,a(iscr)) *d acer.3361 *d acer.3412 *d acer.3562 mf=15 *d acer.3647,3648 *d acer.3650 *d acer.3672 *d acer.3677,3680 *d acer.3800 if (enext.lt.9.e9) go to 170 *d acer.4374 negn=0 *d acer.4406,4408 *d acer.4411 *d acer.4777,4778 *d acer.4803 *d acer.4809,4811 *d acer.4814 *d acer.4817 *i acer.5507 return *d acer.5687 *d acer.5855 *d acer.6581 *d acer.6777 fract=wt(1) *d acer.6878 *d acer.6892,6895 *d acer.6900 *d acer.6907 *d acer.6937 *d acer.7012 data dashes/6h------/ *d acer.7547 *d acer.7715,7716 *d acer.7792 *ident up66 */ gaminr -- 11 mar 93 -- fix problems found by cflint *d gaminr.88 data blank/4h /, iblank/4h /, nz/1/, ipr/1/ *d gaminr.128 if (ngam1.eq.0.or.ngam2.eq.0) go to 150 *d gaminr.238 *d gaminr.381 if (matb.eq.0) go to 500 *d gaminr.506,507 *d gaminr.621,622 do 910 ig=1,ngg 910 write(nsyso,5) ig,egg(ig),egg(ig+1) *d gaminr.661 write(nsyso,10) *d gaminr.995 *d gaminr.1204 ff(1,ng+1)=siginc *ident up67 */ reconr -- 11 mar 92 -- fix some problems found by cflint *d reconr.377 *d reconr.487 *d reconr.492 *d reconr.701 *d reconr.1163,1164 *d reconr.2017 *d reconr.2083 ex=2.*(e-erp)/gtt *d reconr.2165 *d reconr.2285 *ident up68 */ acer -- 16 mar 93 -- allow multiple distributions in mf6 (f19) */ add mf6,lf6 (phase-space law) */ add mf6,lf7 (angle-energy law) *d up25.125 if (mf.eq.6.and.mt.le.91) go to 111 *i acer.1577 c c ***check for multiple neutron emissions in mf6 111 nk=n1h ik=0 112 ik=ik+1 if (ik.gt.nk) go to 119 call tab1io(nin,0,0,a(iscr),nb,nw) zap=c1h lf=l2h if (abs(zap-1.).gt.0.001) go to 119 113 if (nb.eq.0) go to 114 call moreio(nin,0,0,a(iscr),nb,nw) go to 113 114 if (lf.eq.0.or.lf.eq.3.or.lf.eq.4) go to 112 if (lf.ne.6) go to 115 call contio(nin,0,0,a(iscr),nb,nw) go to 112 115 if (lf.eq.7) go to 116 call tab2io(nin,0,0,a(iscr),nb,nw) ne=n2h do 1115 ie=1,ne call listio(nin,0,0,a(iscr),nb,nw) 1114 if (nb.eq.0) go to 1115 call moreio(nin,0,0,a(iscr),nb,nw) go to 1114 1115 continue go to 112 116 call tab2io(nin,0,0,a(iscr),nb,nw) nmu=n2h do 1118 imu=1,nmu call tab2io(nin,0,0,a(iscr),nb,nw) ne=n2h do 1117 ie=1,ne call tab1io(nin,0,0,a(iscr),nb,nw) 1116 if (nb.eq.0) go to 1117 call moreio(nin,0,0,a(iscr),nb,nw) go to 1116 1117 continue 1118 continue go to 112 119 ik=ik-1 mtd=mt call findf(matd,6,mtd,nin) call contio(nin,0,0,a(iscr),nb,nw) n1h=ik go to 121 *d acer.1582 *d acer.1586 *i acer.1627 if (mf.eq.6.and.lf.eq.7) ltt=3 *d acer.1666 153 if (ltt.eq.3) go to 156 if (ltt.eq.1) call listio(nin,0,0,a(l),nb,nw) *i acer.1672 156 call tab2io(nin,0,0,a(l),nb,nw) l=l+nw l2=l nmu=n2h l=l+2*nmu sum=0. do 1160 imu=1,nmu l1=l call tab1io(nin,0,0,a(l),nb,nw) 1157 l=l+nw if (nb.eq.0) go to 1158 call moreio(nin,0,0,a(l),nb,nw) go to 1157 1158 e1=0. e2=1.e9 ir=1 ip=2 call intega(f,e1,e2,a(l1),ip,ir) a(l2+2*imu-2)=a(l1+1) a(l2+2*imu-1)=f nr=a(l1+4) nep=a(l1+5) ll=l1+5+2*nr do 1159 iep=1,nep 1159 a(ll+2*iep)=a(ll+2*iep)/f if (imu.eq.1) go to 1160 ll=l2+2*imu-2 sum=sum+0.5*(a(ll)-a(ll-2))*(a(ll+1)+a(ll-1)) 1160 continue nw=l-loc(j) *i acer.1694 if (ltt.eq.3) go to 164 *i acer.1701 164 l=iscr call tab1io(0,0,nscr,a(l),nb,nw) l=l+nw do 1166 i=1,nmu call tab1io(0,0,nscr,a(l),nb,nw) 1165 l=l+nw if (nb.eq.0) go to 1166 call moreio(0,0,nscr,a(l),nb,nw) go to 1165 1166 continue *d acer.1730 *d up41.22 *i acer.1726 if (ltt.eq.3) go to 190 *d acer.1732,1733 183 if (ik.lt.nk) go to 128 call asend(nout,0) go to 105 *d acer.1749 190 call pttab(ltt,a(iscr),nscr,nout) *d acer.1751 call asend(nout,0) go to 105 *d acer.1773 call asend(nout,0) go to 105 *d up25.139 call asend(nout,0) go to 105 *d acer.2275 subroutine pttab(ltt,a,nin,nout) *i acer.2297 if (ltt.eq.3) nmu=np if (ltt.eq.3) intmu=nint(a(8)) *d acer.2490 go to 590 *i acer.2470 if (ltt.eq.3) a(3)=intmu if (ltt.eq.3) a(4)=nmu *i acer.2564 c c copy rest of the mf6,lf7 subsubsection 590 if (ltt.ne.3) go to 100 do 600 imu=1,nmu call tab1io(nin,nout,0,a,nb,nw) 595 if (nb.eq.0) go to 600 call moreio(nin,nout,0,a,nb,nw) go to 595 600 continue c c continue eprime loop *i up25.155 call tab2io(nin,0,0,a(iscr+6),nb,nw) *i acer.3652 dimension loc(5) *i up25.204 if (mfh.eq.6.and.law.eq.7) go to 304 *i acer.3957 if (law.eq.7) nmu=nint(a(iscr+3)) *i acer.3970 if (law.ne.7) go to 320 do 319 i=1,nmu call tab1io(nin,0,0,a(iscr),nb,nw) 318 if (nb.eq.0) go to 319 call moreio(nin,0,0,a(iscr),nb,nw) go to 318 319 continue *d acer.4234 c ***generalized yield from file 6 *i acer.4239 nk=nint(a(iscr+4)) *i acer.4241 ivar=0. ik=0 521 ik=ik+1 if (ik.gt.nk) go to 524 loc(ik)=jscr *i acer.4242 zap=a(jscr) if (abs(zap-1.).gt.0.001) go to 524 *d up32.4 awp=a(jscr+1) *d acer.4243 law=nint(a(jscr+3)) *d acer.4251,4252 522 jscr=jscr+nw if (nb.eq.0) go to 523 *d acer.4254 go to 522 523 do 1514 j=2,n ii=loc(ik)+5+2*m+2*j if (abs(a(ii)-a(ii-2)).gt.0.001*a(ii)+1.e-6) ivar=1 1514 continue ii=loc(ik)+5+2*m+2 if (ivar.eq.0.and.abs(a(ii)-nint(a(ii))).gt.0.001) ivar=2 call tab2io(nin,0,0,a(jscr),nb,nw) m=nint(a(jscr+4)) n=nint(a(jscr+5)) do 1520 j=1,n call listio(nin,0,0,a(jscr),nb,nw) nmu=l2h 1515 if (nb.eq.0) go to 1516 call moreio(nin,0,0,a(jscr),nb,nw) go to 1515 1516 if (law.ne.7) go to 1520 do 1518 k=1,nmu call tab1io(nin,0,0,a(jscr),nb,nw) 1517 if (nb.eq.0) go to 1518 call moreio(nin,0,0,a(jscr),nb,nw) go to 1517 1518 continue 1520 continue go to 521 524 ik=ik-1 if (ik.gt.1) write(nsyso, 1 '(/'' multiple mf6 subsections found for mt='',i3)') mth if (ivar.eq.1) write(nsyso, 1 '(/'' energy-dependent mf6 yields found for mt='',i3)') mth if (ivar.eq.2) write(nsyso, 1 '(/'' non-integer mf6 yields found for mt='',i3)') mth if (ivar.gt.0) go to 1525 c c ***constant integer yield or yields yield=0. do 1524 j=1,ik ii=loc(j)+7+2*m yield=yield+a(ii) 1524 continue ntyr=(3-2*lct)*nint(yield) xss(tyr+i-1)=(3-2*lct)*ntyr go to 1532 c c ***generalized yield 1525 ntyr=100+next-dlw+1 xss(tyr+i-1)=(3-2*lct)*ntyr igyl=next xnext=0. k=-1 1526 k=k+1 xx=xnext yy=0. xnext=1.e10 ishift=500 do 1527 j=1,ik ii=loc(j) ip=2 ir=1 call terpa(y,xx,xn,idis,a(ii),ip,ir) if (k.eq.0.and.xn.lt.xnext) xnext=1.0001*xn if (k.gt.0.and.xn.lt.xnext.and.xn.gt.1.001*xx) xnext=xn yy=yy+y 1527 continue if (k.le.0) go to 1526 xss(next+2*k+ishift)=xx xss(next+2*k+1+ishift)=yy if (2*k+2.ge.ishift) call error('acelod', 1 'storage exceeded for generalized yield',' ') if (xnext.lt.1.e10) go to 1526 xss(next)=0 xss(next+1)=k do 1530 j=1,k eyl=xss(next+2*j+ishift) if (j.eq.1) eyl=eyl/1.00001 if (j.eq.k) eyl=1.00001*eyl xss(next+1+j)=eyl*1.e-6 xss(next+1+k+j)=xss(next+2*j+1+ishift) 1530 continue next=next+2*k+2 c c ***return to the beginning of the section 1532 call findf(matd,6,mt,nin) call contio(nin,0,0,a(iscr),nb,nw) xss(ldlw+i-1)=next-dlw+1 if (law.ne.7) xss(land+i)=-1 ik=0 529 ik=ik+1 if (ik.gt.nk) go to 592 jscr=iscr call tab1io(nin,0,0,a(jscr),nb,nw) zap=a(iscr) if (abs(zap-1.).gt.0.001) go to 592 awp=a(iscr+1) law=nint(a(iscr+3)) m=nint(a(iscr+4)) n=nint(a(iscr+5)) jnt=nint(a(iscr+7)) 531 jscr=jscr+nw if (nb.eq.0) go to 530 call moreio(nin,0,0,a(jscr),nb,nw) go to 531 *d up25.216,217 530 if (ik.gt.1) xss(last)=next-dlw+1 *i up25.221 if (ivar.gt.0) go to 1550 *d up25.224,225 *d up25.227 550 xss(next+n+k)=a(iscr+5+2*m+2*k)/yield *i acer.4270 go to 554 1550 ngyl=nint(xss(igyl+1)) xss(next)=0 next=next+1 xss(next)=ngyl do 1551 j=1,ngyl eyl=xss(igyl+1+j)*1.e6 if (j.eq.1) eyl=1.00001*eyl if (j.eq.ngyl) eyl=eyl/1.00001 gyl=xss(igyl+1+ngyl+j) ir=1 ip=2 call terpa(y,eyl,en,idis,a(iscr),ip,ir) if (j.eq.1) eyl=eyl/1.00001 if (j.eq.ngyl) eyl=1.00001*eyl xss(next+j)=eyl*1.e-6 xss(next+ngyl+j)=y/gyl 1551 continue next=next+1+2*ngyl 554 continue *d up32.6,7 if (law.eq.1.and.(lang.lt.1.or.(lang.gt.2.and.lang.lt.11) 1 .or.lang.gt.13)) call error('acelod', *i acer.4310 if (law.eq.7) go to 587 *i acer.4350 go to 590 c c ***law 7, angle-energy format 587 jscr=iscr call tab1io(nin,0,0,a(jscr),nb,nw) intmu=l1h nmu=l2h ee=c2h*1.e-6 xss(next+j)=ee xss(nextn+j)=nexd-dlw+1 xss(nexd)=intmu xss(nexd+1)=nmu lein=nexd nexd=nexd+2 mus=nexd nexd=nexd+2*nmu do 1589 imu=1,nmu jscr=iscr call tab1io(nin,0,0,a(jscr),nb,nw) npep=n2h intep=nint(a(iscr+7)) 1584 jscr=jscr+nw if (nb.eq.0) go to 1585 call moreio(nin,0,0,a(jscr),nb,nw) go to 1584 1585 xss(mus+imu-1)=c2h xss(mus+nmu+imu-1)=nexd-dlw+1 xss(nexd)=intep xss(nexd+1)=npep nexd=nexd+1 xss(nexd+1+2*npep)=0. do 1588 ki=1,npep xss(nexd+ki)=a(iscr+8+2*(ki-1))*1.e-6 xss(nexd+npep+ki)=a(iscr+9+2*(ki-1))*1.e6 if (ki.eq.1) go to 1588 if (intep.eq.1) xss(nexd+2*npep+ki)=xss(nexd+2*npep+ki-1) 1 +a(iscr+9+2*(ki-2)) 2 *(a(iscr+8+2*(ki-1))-a(iscr+8+2*(ki-2))) if (intep.eq.2) xss(nexd+2*npep+ki)=xss(nexd+2*npep+ki-1) 1 +0.5*(a(iscr+9+2*(ki-2))+a(iscr+9+2*(ki-1))) 2 *(a(iscr+8+2*(ki-1))-a(iscr+8+2*(ki-2))) 1588 continue nexd=nexd+3*npep+1 1589 continue *d acer.4353 go to 529 592 call tosend(nin,0,0,a(iscr)) *i up32.9 if (lang.ne.2) then write(nsyso,'(/'' converting to kalbach:'', 1 '' mt ='',i4)') mth xss(tyr+i-1)=-xss(tyr+i-1) endif *d up47.207,208 *d up47.212 1 '' eprime='',e12.4)') mth,ee,ep *i up32.111 if (delfcm.lt.0.) then write(strng,'('' for mt='',i3,'' e='',1p,e12.4, 1 '' eprime='',e12.4)') mth,ee,ep call mess('acelod', 1 'converted cm distribution unreasonable',strng) delfcm=0. endif *d acer.4244,4245 if (law.ne.1.and.law.ne.6.and.law.ne.7) call error('acelod', 1 'illegal law for endf6 file6 neutrons',' ') *d up25.220 if (law.eq.1) then xss(next+1)=44 else if (law.eq.7) then xss(next+1)=67 else xss(next+1)=66 endif *i acer.4270 if (law.eq.6) go to 1590 *b acer.4354 go to 595 c c ***phase-space distribution 1590 call contio(nin,0,0,a(iscr),nb,nw) apsx=a(iscr) npsx=nint(a(iscr+5)) xss(next)=npsx xss(next+1)=apsx xss(next+2)=2 next=next+3 lcnt=next step1=10.**0.2 step2=0.02 xx=1.e-5 n=1 1591 n=n+1 if (xx.lt.0.09999) then xx=xx*step1 else xx=xx+step2 endif if (xx.lt.1.00001) go to 1591 nn=n xss(next)=nn xl=0. pl=0. yn=0. n=1 xss(next+n)=xl xss(next+nn+n)=pl xss(next+2*nn+n)=yn xx=1.e-5 1592 n=n+1 if (xx.gt..99999) xx=1.0 pn=sqrt(xx)*(1.0-xx)**(1.5*npsx-4) yn=yn+0.5*(xx-xl)*(pn+pl) xss(next+n)=xx xss(next+nn+n)=pn xss(next+2*nn+n)=yn xl=xx pl=pn if (xx.lt.0.09999) then xx=xx*step1 else xx=xx+step2 endif if (xx.lt.1.00001) go to 1592 sum=yn do 1593 k=1,nn xss(next+nn+k)=xss(next+nn+k)/sum 1593 xss(next+2*nn+k)=xss(next+2*nn+k)/sum next=next+1+3*nn *i acer.5114 if (law.eq.66) go to 381 if (law.eq.67) go to 382 *i acer.5260 c c ***law 66 (phase-space) 381 continue npsx=nint(xss(l)) apsx=xss(l+1) intt=nint(xss(l+2)) nn=nint(xss(l+3)) loci=l+3 write(nsyso,896) npsx,apsx,intt write(nsyso,897) (xss(j+loci),xss(j+nn+loci), 1 xss(j+2*nn+loci),j=1,nn) go to 385 c c ***law 67 (angle-energy) 382 continue j=nint(xss(l)) write(nsyso,685) j if (j.eq.0) go to 383 write(nsyso,690) (nint(xss(l+i)),i=1,j) l=l+j write(nsyso,695) (nint(xss(l+i)),i=1,j) l=l+j 383 l=l+1 ne=nint(xss(l)) write(nsyso,700) ne do 1385 ie=1,ne ein=xss(l+ie) lein=nint(xss(ie+ne+l)+dlw) intmu=nint(xss(lein)) write(nsyso,898) ein,intmu nmu=nint(xss(lein+1)) locmu=lein+1 do 1382 k=1,nmu amu=xss(locmu+k) ll=nint(xss(locmu+nmu+k))+dlw intep=nint(xss(ll)) ll=ll+1 npep=nint(xss(ll)) write(nsyso,900) amu,intep write(nsyso,901) (xss(ll+j),xss(ll+npep+j), 1 xss(ll+2*npep+j),j=1,npep) 1382 continue 1385 continue go to 385 *i acer.5625 896 format(12x,'npsx =',i2/ 1 12x,'apsx =',f10.4/ 2 12x,'intt =',i2) 897 format(1x,3(8x,6henergy,11x,3hpdf,11x,3hcdf)/ 1 1x,3(8x,6h------,11x,3h---,11x,3h---)/(1x,1p,9e14.6)) 898 format(/4x,'ein =',f10.5,4x,'intmu =',i2) 900 format(/5x,'mu =',f10.5,4x,'intep =',i2) 901 format(9x,6henergy,11x,3hpdf,11x,3hcdf/ 1 9x,6h------,11x,3h---,11x,3h---/(1x,1p,3e14.6)) *i acer.6068 if (law.eq.66) go to 416 if (law.eq.67) go to 418 *i acer.6255 c c ***law 66 416 call typen(l,nout,1) l=l+1 call typen(l,nout,2) l=l+1 call typen(l,nout,1) l=l+1 if (nout.ne.1) nn=nint(xss(l)) if (nout.eq.1) nn=iss(l) n=3*nn l=l+1 do 417 k=1,n call typen(l,nout,2) 417 l=l+1 go to 420 c c ***law 67 418 if (nout.ne.1) nr=nint(xss(l)) if (nout.eq.1) nr=iss(l) call typen(l,nout,1) l=l+1 if (nr.eq.0) go to 1412 n=2*nr do 1411 j=1,n call typen(l,nout,1) 1411 l=l+1 1412 if (nout.ne.1) ne=nint(xss(l)) if (nout.eq.1) ne=iss(l) call typen(l,nout,1) l=l+1 do 1414 j=1,ne call typen(l,nout,2) 1414 l=l+1 do 1415 j=1,ne call typen(l,nout,1) 1415 l=l+1 do 1425 j=1,ne call typen(l,nout,1) l=l+1 if (nout.ne.1) nmu=nint(xss(l)) if (nout.eq.1) nmu=iss(l) call typen(l,nout,1) l=l+1 do 1418 k=1,nmu call typen(l,nout,2) 1418 l=l+1 do 1419 k=1,nmu call typen(l,nout,1) 1419 l=l+1 do 1420 k=1,nmu call typen(l,nout,1) l=l+1 if (nout.ne.1) nep=nint(xss(l)) if (nout.eq.1) nep=iss(l) call typen(l,nout,1) l=l+1 nn=3*nep do 1421 n=1,nn call typen(l,nout,2) 1421 l=l+1 1420 continue 1425 continue go to 420 */ acer -- 30 mar 93 -- fix photon law *d acer.4666 c law=4 *d up41.4 xss(nex+1)=4 *d up25.237 xss(nex+6)=a(iscr+4+2*m+2*n)*1.e-6 */ acer --1 apr 93 --fix discontinuity handling *d up47.137 140 if (idis.eq.0) go to 150 *d up47.148 if (enext.gt..99e10) jt=-jt *d acer.1202 if (jt.gt.0) go to 120 *d up47.158 eg=a(idisc-1+iedis) *ident up69 */ reconr -- 30 mar 93 -- fix reich-moore problem *d up2.34 if (ll.ne.0.and.(fl.gt.spi-.5.and.fl.le.spi)) then */ reconr -- 29 apr 93 -- make sure res boundaries are clean *i reconr.1158 common/reslim/eresl,eresr,eresh,eresu,eresm *i reconr.1338 if (kr+ir.ge.n2h) go to 255 *d reconr.1452 if (ig.eq.1.or.ig.eq.ngo) go to 415 egt=sigfig(eg,ndigit,0) if (egt.eq.eresl) go to 430 if (egt.eq.eresr) go to 430 if (egt.eq.eresu) go to 430 if (egt.eq.eresm) go to 430 if (egt.eq.eresh) go to 430 415 egl=eg *ident up70 */ wimsr -- 2 feb 92 -- major rework of wimsr */ to include suggestions from u of i */ and aguilar of inin. *i wimsr.13 c * * *i wimsr.21 c * igroup group option * c * 0=69 groups (default) * c * 9=user's choice * c * * c * card 2a (igroup.eq.9 only) * c * ngnd number of groups * c * nfg number of fast groups * c * nrg number of resonance groups * c * igref reference group (default is last fast group) * c * * *d wimsr.23,42 c * mat endf mat number of the material to be processed * c * nfid identification of material for the wims library * c * rdfid identification number for the resonance data * c * iburn burnup data option * c * 0=no burnup data is provided (default) * c * 1=burnup data is provided in cards 5 and 6 * c * * c * card 4 * *d wimsr.44 c * (0=all found on input tape) * *d wimsr.46 c * (0=all found on input tape) * *d wimsr.65,87 c * * c * following two cards are for iburn gt 0 only * c * card 5 * c * ntis no. of time-dependent isotopes * c * efiss energy released per fission * c * * c * card 6 (repeat this card ntis times) * c * identa ident of fission product isotope * c * yield fission yield of identa from burnup of mat * c * * c * card 7 * c * lambda resonance-group goldstein lambdas (13 for * c * default 69-group structure, nrg otherwise). * *d wimsr.92 common/wunits/ngendf,nout,nscr0,nscr1,nscr2,nscr3,nscr4 *d wimsr.94 *d wimsr.96,97 common/wim1/ngnd,nfg,nrg,igref,iprint, 1 ntemp,nsigz,ntmax,tempr(10),nwscr *d wimsr.99,103 common/wim3/nfid,rdfid common/wim4/iu,nfiss,isof,jfiss common/wim6/mat,iznum,awr,sgrf,sigp,mti,mtc,ip1opt,inorf,ires common/burnup/yield(50),ifisp(50) dimension z(15) *d wimsr.108,112 10 format(/' wimsr...create data for wims',40x,f8.1,1hs) ntmax=10 nwscr=7500 namax=100000 *i wimsr.115 nscr0=10 nscr1=11 nscr2=12 nscr3=13 nscr4=14 *d wimsr.118 *d wimsr.126,129 if (ntty.gt.0) write(ntty, 1 '(/'' enter iprint[0], iverw[4], igroup[0]'')') nz=3 *i wimsr.131 z(3)=0. *d wimsr.135,137 igroup=nint(z(3)) write(nsyso,'(/ 1 '' input gendf unit ..................... '',i10/ 2 '' output unit .......................... '',i10/ 3 '' wims version ......................... '',i10/ 4 '' group structure ...................... '',i10)') 5 ngendf,nout,iprint,iverw,igroup if (igroup.eq.0) then *d wimsr.140,144 *d wimsr.145,154 nrg=13 igref=nfg else if (igroup.eq.9) then if (ntty.gt.0) write(ntty, 1 '(/'' enter ngnd, nfg, nrg, igref[nfg]'')') nz=4 z(4)=0 call free(nsysi,z,nz,4) ngnd=nint(z(1)) nfg=nint(z(2)) nrg=nint(z(3)) igref=nint(z(4)) if (igref.eq.0) igref=nfg *d wimsr.156,188 write(nsyso,'(/ 1 '' ngnd ................................. '',i10/ 2 '' nfg .................................. '',i10/ 3 '' nrg .................................. '',i10/ 4 '' igref ................................ '',i10)') 5 ngnd,nfg,nrg,igref ixs=1 call openz(nout,1) call repoz(nout) call reserv('scr',nwscr,iscr,a) *d wimsr.191,195 if (ntty.gt.0) write(ntty, 1 '(/''enter mat, nfid, rdfid, iburn[0]'')') nz=4 z(4)=0 call free(nsysi,z,nz,4) mat=nint(z(1)) nfid=nint(z(2)) rdfid=z(3) iburn=nint(z(4)) write(nsyso,'(/ 1 '' endf mat number ...................... '',i10/ 2 '' wims material identifier ............. '',i10/ 3 '' wims resonance identifier ............ '',f12.1/ 4 '' burn data (0=no, 1=yes, def=0) ....... '',i10)') 5 mat,nfid,rdfid,iburn *d wimsr.196,199 if (ntty.ne.0) write(ntty, *d wimsr.201,322 2 '' mti, mtc, ip1opt[1], norf[0], isof[0], ifprod[0].'')') nz=11 z(8)=1 z(9)=0 z(10)=0 z(11)=0 call free(nsysi,z,nz,4) ntemp=nint(z(1)) nsigz=nint(z(2)) sgref=z(3) ires=nint(z(4)) sigp=z(5) mti=nint(z(6)) mtc=nint(z(7)) ip1opt=nint(z(8)) inorf=nint(z(9)) isof=nint(z(10)) ifprod=nint(z(11)) write(nsyso,'(/ 1 '' no. temperatures ..................... '',i10/ 2 '' no. sigma zeroes ..................... '',i10/ 3 '' reference sigma zero ................. '',1p,e10.2/ 4 '' resonance absorber flag (0=no, 1=yes) '',i10/ 5 '' pot. scatt. cross section ............ '',f10.2/ 6 '' thermal inelastic mt ................. '',i10/ 7 '' thermal elastic mt ................... '',i10/ 8 '' p1 matrix option (0=yes, 1=no) ....... '',i10/ 9 '' resonance fission option (0=yes, 1=no) '',i10/ a '' fission spectrum option (0=no, 1=yes) '',i10/ b '' fission product indicator ............ '',i10)') c ntemp,nsigz,sgref,ires,sigp,mti,mtc, d ip1opt,inorf,isof,ifprod sgrf=sgref c c ***input burn data yield(1)=0.0 ifisp(1)=nfid if (iburn.eq.0) go to 166 if (ntty.gt.0) write(ntty,'(/''enter no. of time-dependent'', 1 '' nuclides and energy released per fission'')') nz=2 call free(nsysi,z,nz,4) ntis=nint(z(1)) efiss=z(2) if (ntty.gt.0) write(ntty,'(/'' enter fission-product ident and'', 1 '' yield'',i2,''times'')') ntis jcc=ntis*2+4 yield(4)=efiss ifisp(4)=0 ij=1 do 150 i=1,ntis nz=2 call free(nsysi,z,nz,4) if (i.eq.1.or.i.eq.3) ij=ij+1 ifisp(ij)=nint(z(1)) yield(ij)=z(2) ij=ij+1 150 continue 166 continue *d wimsr.324,344 *d wimsr.346,356 call reserv('glam',nrg,iglam,a) if (ntty.gt.0) write(ntty,'(/'' enter '',i2,'' lambdas)'')') nrg jscr=iscr *d wimsr.364 c ***read in the cross section data *d wimsr.359 *d wimsr.361,362 190 a(iglam+j-1)=a(jscr+j-1) *d wimsr.365,379 *d wimsr.380,391 call wminit(a) *d wimsr.395,397 *d wimsr.402,404 *d wimsr.409 if (iburn.eq.0) jcc=2 if (nout.ne.0) call wimout(jcc,nfid,rdfid) *d wimsr.413,414 call usag(a) *i wimsr.416 15 format(69x,f8.1,1hs/1x,7(10h**********),7h*******) *i wimsr.417 16 format(/1x,10h**********,f8.1,1hs,11h **********) *d wimsr.420,456 *d wimsr.459,468 subroutine wminit(a) *d wimsr.473,476 common/wunits/ngendf,nout,nscr0,nscr1,nscr2,nscr3,nscr4 common/wim1/ngnd,nfg,nrg,igref,iprint, 1 ntemp,nsigz,ntmax,tempr(10),nwscr *d wimsr.477,478 common/wim4/iu,nfiss,isof,ifiss common/wim6/mat,iznum,awr,sgrf,sigp,mti,mtc,ip1opt,inorf,ires *i wimsr.482 c ***read gendf tape for groups and flags *d wimsr.486,489 *d wimsr.490 *d wimsr.493,494 *d wimsr.497 *d wimsr.498,510 *d wimsr.512,522 if (math.ne.mat) call error('wminit', 1 'desired material is not on gendf tape',' ') *d wimsr.524,525 iznum=iza/1000 awr=c2h*1.0086652 *d wimsr.527,528 if (iu.eq.1) then nw=nrg+nfg call reserv('uff',nw,iuff,a) endif *i wimsr.531 tempr(1)=c1h *d wimsr.534 if (ngn.ne.ngnd) call error('wminit', 1 'incorrect group structure',' ') *d wimsr.537,547 *d wimsr.556,557 *d wimsr.584,618 150 continue c c ***check the number of temperatures and sigma zeroes jtemp=1 160 call contio(ngendf,0,0,a(iscr),nb,nw) if (math.ne.mat) go to 170 jtemp=jtemp+1 call listio(ngendf,0,0,a(iscr),nb,nw) tempr(jtemp)=c1h call tomend(ngendf,0,0,a(iscr)) go to 160 170 if (ntemp.eq.0) ntemp=jtemp if (ntemp.gt.jtemp) ntemp=jtemp if (nsigz.eq.0) nsigz=nz if (nsigz.gt.nz) nsigz=nz write(nsyso,'(/'' processing'',i3,'' temperatures''/ 1 '' processing'',i3,'' sigma zeroes'')') ntemp,nsigz nwflx=ntemp*nsigz*nrg call reserv('flux',nwflx,iflux,a) c c ***print out group structure if (iprint.gt.0) write(nsyso,'(/ 1 '' group structure''/'' ---------------''/(1p,6e12.4))') 2 (a(iegb-1+i),i=1,ngnd1) *d wimsr.620,630 *d wimsr.637 common/wunits/ngendf,nout,nscr0,nscr1,nscr2,nscr3,nscr4 *d wimsr.640,643 common/wim1/ngnd,nfg,nrg,igref,iprint, 1 ntemp,nsigz,ntmax,tempr(10),nwscr common/wim6/mat,iznum,awr,sgrf,sigp,mti,mtc,ip1opt,inorf,ires *d wimsr.645 *d wimsr.651 if (ires.eq.0) go to 510 *d wimsr.654,661 20 format(/' ***resonance integrals***') ntsr=ntemp*nsigz*nrg call reserv('sabs',ntsr,locabs,a) call reserv('snus',ntsr,locnus,a) call reserv('nsf',ntsr,locsf,a) nwflxr=ntemp*nsigz *d wimsr.663 call reserv('sigz',nsigz,isigz,a) *d wimsr.665,667 call reserv('abs',ntemp,iabs,a) nwelas=ntemp*nsigz*nrg *d wimsr.670,671 nwfa=ntemp*nsigz *d wimsr.673,683 call findex('glam',iglam,a) call findex('scr',iscr,a) *d wimsr.685,686 *d wimsr.687,689 *d wimsr.690,703 *d wimsr.707,713 c ***read gendf data *d wimsr.717,725 *d wimsr.726,729 *d wimsr.730,732 *d wimsr.735 nwflx=ntemp*nsigz*nrg *i wimsr.736 a(i-1+locabs)=0. a(i-1+locnus)=0. a(i-1+locsf)=0.