Skip to content

Feature reconr dla #374

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dlopezaldama wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Feature reconr dla #374

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c7c97b8
Fixes #373: Linearization of cross sections for incident charged part…
dlopezaldama Aug 14, 2025
173036b
Update reference tapes and output files for tests 03,04,06,14,48,50,5…
dlopezaldama Aug 15, 2025
5c3cbb9
Update reference tapes for tests 84 & 85
dlopezaldama Aug 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
164 changes: 100 additions & 64 deletions src/reconr.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,10 @@ module reconm
real(kr)::zai,el,eh,err,errmax,errint,za,awr,tempr,q18
real(kr)::tempi
real(kr)::elis,sta,efmax
integer::lis,lis0,nfor,lrel,nver
integer::lis,lis0,nfor,lrel,nver,nsub
integer::lfw,mata,itype,lrp,lfi,lssf,lrx
integer,parameter::nmtmax=10
integer::mtr4,mtr18,mtr(nmtmax),mtrt(nmtmax),nmtr
integer::mtr3,mtr4,mtr18,mtr(nmtmax),mtrt(nmtmax),nmtr
integer::mt103,mt104,mt105,mt106,mt107
integer::mpmin,mpmax,mdmin,mdmax,mtmin,mtmax,m3min,m3max,m4min,m4max
integer::nxc,ngo,mtr522,ncards
Expand Down Expand Up @@ -131,7 +131,7 @@ subroutine reconr
use samm ! provides s2sammy,desammy
use mainio ! provides nsysi,nsyso,nsyse
! internals
integer::i,nb,nw,nrtot,nwscr,n6,nx,nsub,iold,inew,ngrid
integer::i,nb,nw,nrtot,nwscr,n6,nx,iold,inew,ngrid
integer::nendf,npend,nscr1,nscr2,nscr3,nscr4,intunr
real(kr)::time
real(kr)::rlabel(17),z(17)
Expand Down Expand Up @@ -509,8 +509,10 @@ subroutine anlyzd(nin,nx)

!--select desired redundant reactions.
igam=0
maxres=0
nxn=0
nmtr=0
mtr3=0
mtr4=0
mtr18=0
mt103=0
Expand Down Expand Up @@ -546,67 +548,81 @@ subroutine anlyzd(nin,nx)
if (nmtr.ge.nmtmax)&
call error('anlyzd','too many redundant reactions.',' ')
mfi=nint(dict(i+2))
if (mfi.eq.2) then
mti=nint(dict(i+3))
if (mti.eq.151) maxres=12*nint(dict(i+4))
mti=nint(dict(i+3))
if (mfi.eq.2.and. mti.eq.151) then
maxres=12*nint(dict(i+4))
else if (mfi.eq.3) then
nxn=nxn+1
if (nmtr.eq.0) mtr(1)=1
if (nmtr.eq.0) nmtr=1
mti=nint(dict(i+3))
if (mti.eq.19) then
nmtr=nmtr+1
mtr(nmtr)=18
mtr18=1
endif
if (mti.ge.51.and.mti.le.91.and.mtr4.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=4
mtr4=1
endif
if (mti.ge.mpmin.and.mti.le.mpmax.and.mt103.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=103
mt103=1
endif
if (mti.ge.mdmin.and.mti.le.mdmax.and.mt104.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=104
mt104=1
endif
if (mti.ge.mtmin.and.mti.le.mtmax.and.mt105.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=105
mt105=1
endif
if (mti.ge.m3min.and.mti.le.m3max.and.mt106.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=106
mt106=1
endif
if (mti.ge.m4min.and.mti.le.m4max.and.mt107.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=107
mt107=1
endif
if (nsub.eq.10) then
if (nmtr.eq.0) then
mtr(1)=1
nmtr=1
endif
if (mti.eq.19) then
nmtr=nmtr+1
mtr(nmtr)=18
mtr18=1
endif
if (mti.eq.3.and.mtr3.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=3
mtr3=1
endif
if (mti.ge.51.and.mti.le.91.and.mtr4.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=4
mtr4=1
endif
if (mti.ge.mpmin.and.mti.le.mpmax.and.mt103.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=103
mt103=1
endif
if (mti.ge.mdmin.and.mti.le.mdmax.and.mt104.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=104
mt104=1
endif
if (mti.ge.mtmin.and.mti.le.mtmax.and.mt105.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=105
mt105=1
endif
if (mti.ge.m3min.and.mti.le.m3max.and.mt106.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=106
mt106=1
endif
if (mti.ge.m4min.and.mti.le.m4max.and.mt107.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=107
mt107=1
endif
endif
else if (mfi.eq.10) then
nxn=nxn+1
else if (mfi.eq.12) then
nxn=nxn+1
mti=nint(dict(i+3))
if (mti.eq.3) nmtr=nmtr+1
if (mti.eq.3) mtr(nmtr)=3
if (mti.eq.3.and.mtr3.eq.0) then
nmtr=nmtr+1
mtr(nmtr)=3
mtr3=1
endif
else if (mfi.eq.13) then
nxn=nxn+1
else if (mfi.eq.23) then
nxn=nxn+1
igam=1
if (iverf.lt.6) then
nsub=3
zain=0
awin=0
endif
if (nmtr.eq.0) mtr(1)=501
if (nmtr.eq.0) nmtr=1
if (nmtr.eq.0) then
mtr(1)=501
nmtr=1
endif
if (mtr522.eq.0) then
mti=nint(dict(i+3))
if (mti.ge.534.and.mti.le.572) then
Expand Down Expand Up @@ -650,7 +666,7 @@ subroutine anlyzd(nin,nx)
mts(1)=451
ncs(1)=ncards+2
nxc=1
if (igam.eq.1) return
if (igam.eq.1.or.maxres.eq.0) return
mfs(2)=2
mts(2)=151
ncs(2)=4
Expand Down Expand Up @@ -1815,10 +1831,16 @@ subroutine lunion(nin,nout,ngrid,ngr)
allocate(bufn(nbuf))
allocate(x(ndim))
allocate(y(ndim))
if (allocated(scr)) deallocate(scr)
allocate(scr(npage+50))

! this value fits 1/v to within err
stpmax=1+sqrt(ovfact*err)
if (nsub.eq.10.or.nsub.eq.3) then
! this value fits 1/v to within err for neutrons and photoatomic data
stpmax=1+sqrt(ovfact*err)
else
! value used for the union grid of incident charged-particles in acefc
stpmax=1.2
endif

!--copy nodes from the global area to *old* tape
iold=14
Expand Down Expand Up @@ -1879,7 +1901,8 @@ subroutine lunion(nin,nout,ngrid,ngr)
if (mfh.eq.12) scr(5)=1
if (mfh.eq.13) scr(5)=1
if (mfh.ne.3) go to 180
if (mth.eq.1.or.mth.eq.3) go to 150
if (mth.eq.1.and.nsub.ne.0) go to 150
if (mth.eq.3.and.mtr3.gt.0) go to 150
if (mth.eq.4.and.mtr4.gt.0) go to 150
if (mth.eq.103.and.mt103.gt.0) go to 150
if (mth.eq.104.and.mt104.gt.0) go to 150
Expand Down Expand Up @@ -1916,19 +1939,20 @@ subroutine lunion(nin,nout,ngrid,ngr)
if (mth.eq.19.and.mtr18.gt.0) q18=qx
if (awin.ne.0) then
thrx=-qx*(awrx+1)/awrx
thrxx=sigfig(thrx,7,+1)
else
thrx=-qx
thrxx=thrx
endif
thr6=thrx
if (thr6.lt.zero) thr6=0
if (thrx.le.zero) go to 190
thrxx=sigfig(thrx,7,+1)
l=7+2*nint(scr(5))
if (scr(l).ge.thrxx) go to 190
if (scr(l+1).ne.zero) then
write(text,'(''xsec nonzero at threshold for mt='',i3)') mth
call mess('lunion',text,'adjusted using jump in xsec')
endif
if (scr(l).ge.thrxx) go to 190
thrx=thrxx
write(nsyso,'(/&
&'' changed threshold from'',1p,e13.6,'' to'',e13.6,&
Expand Down Expand Up @@ -1973,7 +1997,7 @@ subroutine lunion(nin,nout,ngrid,ngr)
207 continue
ernext=scr(ibase+ir*2+1)
srnext=scr(ibase+ir*2+2)
if (sr.lt.ssmall.and.srnext.lt.ssmall.and.ir.lt.npr-1) go to 205
if (sr.lt.ssmall.and.srnext.lt.ssmall.and.(lr+ir).lt.npr-1) go to 205
! check for initial discontinuity
if (abs(er-ernext).gt.small*er) go to 210
er=sigfig(er,7,0)
Expand Down Expand Up @@ -2262,7 +2286,7 @@ subroutine resxs(ngrid,nout,nscr,nrtot)
real(kr),dimension(:),allocatable::bufr,bufg,bufl
real(kr),dimension(:),allocatable::x,y
real(kr),dimension(:,:),allocatable::sigs
integer,parameter::ndim=30
integer,parameter::ndim=50
real(kr),parameter::half=0.5e0_kr
real(kr),parameter::estp=4.1e0_kr
real(kr),parameter::small=1.e-10_kr
Expand Down Expand Up @@ -4657,7 +4681,7 @@ subroutine emerge(nin,ngrid,ngp,nres,nrtot,iold,inew,nscr)
! internals
integer::nneg,ntot,i,in,ig,inn,nss,iss,nb,nw,idis,it
integer::imtr,itt,k,istart,iend,j,ib,isave,ir,ith
real(kr)::er,eg,en,e,thresh,sn,sg,enext,awrx,qx
real(kr)::er,eg,en,e,thresh,sn,sg,enext,awrx,qx,onethr
real(kr)::res(nsig+1),tot(10)
real(kr)::aa(1)
real(kr),dimension(:),allocatable::bufo,bufn,bufg,bufr
Expand All @@ -4673,12 +4697,14 @@ subroutine emerge(nin,ngrid,ngp,nres,nrtot,iold,inew,nscr)
allocate(bufn(nbuf))
allocate(bufg(nbufg))
allocate(bufr(nbufr))
if (allocated(scr)) deallocate(scr)
allocate(scr(npage+50))
nneg=0
ntot=nmtr+1
do i=1,nmtr
tot(i+1)=0
enddo
onethr=sigfig(one,7,0)

!--assign scratch units.
iold=14
Expand Down Expand Up @@ -4730,6 +4756,7 @@ subroutine emerge(nin,ngrid,ngp,nres,nrtot,iold,inew,nscr)
nsc=0
210 continue
call contio(nin,0,nscr,scr,nb,nw)
if (mfh.ne.23.and.awin.ne.zero) awrx=c2h/awin
nss=1
if (mfh.eq.10) nss=n1h
iss=nss
Expand All @@ -4739,7 +4766,6 @@ subroutine emerge(nin,ngrid,ngp,nres,nrtot,iold,inew,nscr)
220 continue
if (mth.eq.0) go to 210
if (iss.eq.nss) then ! needed for MF10
if (mfh.ne.23.and.awin.ne.zero) awrx=c2h/awin
qx=c2h
if (awin.ne.0) then
thr6=-qx*(awrx+1)/awrx
Expand Down Expand Up @@ -4791,7 +4817,9 @@ subroutine emerge(nin,ngrid,ngp,nres,nrtot,iold,inew,nscr)
sn=0
if (thresh-eg.gt.test*thresh.and.itype.eq.0) go to 370
call gety1(eg,enext,idis,sn,nin,scr)
if (thresh.gt.one.and.abs(thresh-eg).lt.test*thresh) sn=0
! set zero cross section at threshold for neutron and photoatomic data
if ((nsub.eq.10.or.nsub.eq.3).and.&
& thresh.gt.onethr.and.abs(thresh-eg).lt.test*thresh) sn=0
! backgrounds in a range of unresolved-smooth overlap
! are arbitrarily assigned to the unresolved component
! this only applies to total, elastic, fission and capture.
Expand Down Expand Up @@ -4831,7 +4859,8 @@ subroutine emerge(nin,ngrid,ngp,nres,nrtot,iold,inew,nscr)
endif
sn=sigfig(sn,7,0)
tot(2)=sn
if (ith.eq.0.and.sn.gt.zero) ith=in
! negative values are allowed for charged particle scattering if LTP>2
if (ith.eq.0.and.(sn.gt.zero.or.(nsub.ge.10010.and.mth.eq.2))) ith=in
inn=in
if (ig.eq.ngo) inn=-in
call loada(inn,tot,2,ngrid,bufg,nbufg)
Expand Down Expand Up @@ -5006,9 +5035,11 @@ subroutine recout(iold,nscr,nrtot)
real(kr),parameter::zero=0

!--initialize.
if (allocated(scr)) deallocate(scr)
allocate(scr(npage+50))
allocate(bufo(nbuf))
allocate(bufn(nbuf))
if (allocated(bufl)) deallocate(bufl)
allocate(bufl(nbufl))
i152=0
if (nunr.gt.0) i152=1
Expand All @@ -5018,7 +5049,7 @@ subroutine recout(iold,nscr,nrtot)
! none found increment nxc since we'll insert a dummy mf3, mt1
! section into the output tape later.
no3=0
if (nint(zain).eq.1) then
if (nsub.eq.10) then
do i=1,nxc
if (mfs(i).eq.3)no3=1
enddo
Expand Down Expand Up @@ -5055,7 +5086,7 @@ subroutine recout(iold,nscr,nrtot)
scr(2)=efmax
scr(3)=lrel
scr(4)=0
scr(5)=int(10*zain)
scr(5)=nint(10*zain)
scr(6)=nver
call contio(0,nout,0,scr,nb,nw)
endif
Expand Down Expand Up @@ -5089,7 +5120,7 @@ subroutine recout(iold,nscr,nrtot)
nc=nc+17
enddo
nxcc=nxc
if (nint(zain).eq.1.and.no3.eq.0) nxcc=nxcc-1
if (nsub.eq.10.and.no3.eq.0) nxcc=nxcc-1
scr(6)=nxcc+nmtr+i152
call hdatio(0,nout,0,scr,nb,nw)
no2=1
Expand Down Expand Up @@ -5119,7 +5150,7 @@ subroutine recout(iold,nscr,nrtot)
imtr=imtr+1
enddo
endif
if (no3.eq.0.and.nint(zain).eq.1.and.mfs(i).gt.3)then
if (no3.eq.0.and.nsub.eq.10.and.mfs(i).gt.3)then
no3=-1
dict(j+1)=0
dict(j+2)=0
Expand Down Expand Up @@ -5257,6 +5288,11 @@ subroutine recout(iold,nscr,nrtot)
if(mfl.eq.0) go to 272
goto 270
endif

!--no redundant reaction calculation for incident
! charged particles or photonuclear data
if (nsub.ge.10010.or.nsub.eq.0) goto 270

mth=1
if (awin.eq.0) then
mth=3
Expand Down
Loading