Skip to content

Fixes #383: Implemented required updates to enhance code robustness, … #384

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 1 commit into
base: main
Choose a base branch
from
Open

Fixes #383: Implemented required updates to enhance code robustness, … #384

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
41 changes: 28 additions & 13 deletions src/acecm.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ module acecm
private

! main ace container array
integer,parameter,public::nxss=50000000
integer,parameter,public::nxss=100000000
real(kr),dimension(:),allocatable,public::xss

!--Public routines
Expand All @@ -30,13 +30,13 @@ subroutine mtname(mt,name,izai)
integer::i
character(10),dimension(500),parameter::hndf=(/&
'total ','elastic ','nonelastic','inelastic ','(n,x) ',&
'(n,1/2*1) ','(n,1/2*2) ','(n,1/2*3) ','(n,1/2*4) ','(n,x) ',&
'(n,1/2*1) ','(n,1/2*2) ','(n,1/2*3) ','(n,1/2*4) ','(n,cont) ',&
'(n,2nd) ','(n,x) ','(n,x) ','(n,x) ','(n,x) ',&
'(n,2n) ','(n,3n) ','fission ','(n,f) ','(n,n*f) ',&
'(n,2nf) ','(n,n*)a ','(n,n*)3a ','(n,2n)a ','(n,3n)a ',&
'(n,2n)iso ','(n,abs) ','(n,n*)p ','(n,n*)2a ','(n,2n)2a ',&
'(n,x) ','(n,n*)d ','(n,n*)t ','(n,n*)he3 ','(n,n*)d2a ',&
'(n,n*)t2a ','(n,4n) ','(n,3nf) ','(n,x) ','(n,x) ',&
'(n,2nf) ','(n,na) ','(n,n3a) ','(n,2na) ','(n,3na) ',&
'(n,2n)iso ','(n,abs) ','(n,np) ','(n,n2a) ','(n,2n2a) ',&
'(n,x) ','(n,nd) ','(n,nt) ','(n,nhe3) ','(n,nd2a) ',&
'(n,nt2a) ','(n,4n) ','(n,3nf) ','(n,x) ','(n,x) ',&
'(n,2np) ','(n,3np) ','(n,x) ','(n,n2p) ','(n,npa) ',&
'(n,2/2*1) ','(n,2/2*2) ','(n,2/2*3) ','(n,2/2*4) ','(n,n*0) ',& !50
'(n,n*1) ','(n,n*2) ','(n,n*3) ','(n,n*4) ','(n,n*5) ',&
Expand All @@ -49,7 +49,7 @@ subroutine mtname(mt,name,izai)
'(n,n*36) ','(n,n*37) ','(n,n*38) ','(n,n*39) ','(n,n*40) ',&
'(n,n*c) ','(n,x) ','(n,x) ','(n,x) ','(n,x) ',&
'(n,x) ','(n,x) ','(n,x) ','(n,n*)gma ','(n,x) ',& !100
'(n,parab) ','(n,gma) ','(n,p) ','(n,d) ','(n,t) ',&
'(n,disap) ','(n,gma) ','(n,p) ','(n,d) ','(n,t) ',&
'(n,he3) ','(n,a) ','(n,2a) ','(n,3a) ','(n,x) ',&
'(n,2p) ','(n,pa) ','(n,t2a) ','(n,d2a) ','(n,pd) ',&
'(n,pt) ','(n,da) ','(n,x) ','(n,x) ','(n,dest) ',&
Expand Down Expand Up @@ -207,10 +207,25 @@ subroutine mtname(mt,name,izai)

!--alternate name when processing incident charged particle files
if (izai.gt.1) then
if (mt.gt.3) then
if (izai.eq.1001) then
name(1:2)='(p'
elseif (izai.eq.1002) then
name(1:2)='(d'
elseif (izai.eq.1003) then
name(1:2)='(t'
elseif (izai.eq.2003) then
name(1:2)='(s'
elseif (izai.eq.2004) then
name(1:2)='(a'
else
name(1:2)='(z'
endif
endif
if (mt.eq.4) then
name='(z,n) '
name(3:10)=',n) '
elseif (mt.eq.5) then
name='(z,x) '
name(3:10)=',x) '
elseif ((izai.eq.1001.and.mt.eq.103).or.&
(izai.eq.1002.and.mt.eq.104).or.&
(izai.eq.1003.and.mt.eq.105).or.&
Expand Down Expand Up @@ -697,16 +712,16 @@ subroutine advance_to_locator(nout,l,locator)
! internals
character(66)::text
if (l.lt.locator) then
write(text,'(''expected xss index ('',i6,'') greater than '',&
&''current index ('',i6,'')'')') locator, l
write(text,'(''expected xss index ('',i9,'') > '',&
&''current index ('',i9,'')'')') locator, l
call mess('advance',text,'xss array was padded accordingly')
do while (l.lt.locator)
call typen(l,nout,1)
l=l+1
enddo
else if (l.gt.locator) then
write(text,'(''expected xss index ('',i6,'') less than '',&
&''current index ('',i6,'')'')') locator, l
write(text,'(''expected xss index ('',i9,'') < '',&
&''current index ('',i9,'')'')') locator, l
call error('advance',text,'this may be a serious problem')
endif
return
Expand Down
3 changes: 2 additions & 1 deletion src/endf.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -525,8 +525,9 @@ subroutine tpidio(nin,nout,nscr,a,nb,nw)
integer::inin,inout,inscr,i

!--input.
inin=iabs(nin)
if (inin.ne.0) rewind inin
if (nin.lt.0) then
inin=iabs(nin)
read(inin) math,mfh,mth,nb,nw,(a(i),i=1,17)
else if (nin.gt.0) then
read(nin,'(16a4,a2,i4,i2,i3,i5)') (hb(i),i=1,17),math,mfh,mth
Expand Down
7 changes: 4 additions & 3 deletions src/leapr.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -423,7 +423,7 @@ subroutine leapr
isym=0
if (ncold.ne.0) isym=1
if (isabt.eq.1) isym=isym+2

! Based on endout, to write the actual TSL data, the max number of entries
! needed in scr is either 8+2*nalpha, or 8+2*nedge. However, we have no way
! of knowing how many comment lines were added to the leaper input. The
Expand Down Expand Up @@ -522,7 +522,7 @@ subroutine contin(temp,itemp,np,maxn)
if (maxn.gt.maxnphon) then
call timer(time)
write(nsyse,'(/'' performing phonon expansion sum'',&
&37x,f8.1,''s'')'),time
&37x,f8.1,''s'')')time
endif
do n=2,maxn
npn=np+npl-1
Expand Down Expand Up @@ -559,7 +559,7 @@ subroutine contin(temp,itemp,np,maxn)
if (maxn.gt.maxnphon) then
call timer(time)
write(nsyse,'(/'' done with phonon expansion sum'',&
&38x,f8.1,''s'')'),time
&38x,f8.1,''s'')')time
endif

!--print out start of sct range for each beta
Expand Down Expand Up @@ -1139,6 +1139,7 @@ real(kr) function terps(sd,nsd,delta,be)
terps=0
if (be.gt.delta*nsd) return
i=int(be/delta)
if (i.lt.0) return
if (i.lt.nsd-1) then
bt=i*delta
btp=bt+delta
Expand Down
2 changes: 1 addition & 1 deletion src/main.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ program njoy
!-----------------------------------------------------------------------
!
! NJOY Nuclear Data Processing System
! Version 2016.20
! Version 2016 (check vers.f90)
!
!-----------------------------------------------------------------------
!
Expand Down
4 changes: 2 additions & 2 deletions src/plotr.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -383,9 +383,9 @@ subroutine plotr
real(kr)::st1,st2,st3,st4,stmax,stmin,ststp,temp,dleth,dener
real(kr)::sigig,flag,factx1,facty1,temper2,reset,enxt,zz,z2,itypx
real(kr)::z(15)
integer,parameter::mmax=20000 !same in plotr and viewr
integer,parameter::mmax=50000 !same in plotr and viewr
integer,parameter::nwamax=45000
integer,parameter::maxaa=200000
integer,parameter::maxaa=500000
real(kr),dimension(nwamax)::a
real(kr),dimension(maxaa)::aa
real(kr),dimension(mmax)::x,y,b,dxm,dxp,dym,dyp
Expand Down
3 changes: 3 additions & 0 deletions src/util.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,6 +193,8 @@ subroutine openz(lun,new)
if (nun.ge.10.and.nun.le.19) then
! scratch units
age='scratch'
inquire(nun,opened=there)
if (there) close (nun)
open(nun,form=for,status=age)
else
! regular units
Expand Down Expand Up @@ -377,6 +379,7 @@ real(kr) function sigfig(x,ndig,idig)
if (x.ne.zero) then
aa=log10(abs(x))
ipwr=int(aa)
if ((ndig.eq.9).and.(aa.gt.zero).and.(aa-ipwr.lt.1.0e-6_kr)) ipwr=ipwr-1
if (aa.lt.zero) ipwr=ipwr-1
ipwr=ndig-1-ipwr
ii=nint(x*ten**ipwr+ten**(ndig-11))
Expand Down
11 changes: 6 additions & 5 deletions src/viewr.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -339,8 +339,8 @@ subroutine viewr
real(kr)::factx,facty
real(kr)::xx,yy,zz
character(80)::text
integer,parameter::mmax=20000 !same in plotr and viewr
integer,parameter::maxaa=500000
integer,parameter::mmax=50000 !same in plotr and viewr
integer,parameter::maxaa=5000000
real(kr),dimension(15)::z
real(kr),dimension(maxaa)::aa
real(kr),dimension(mmax)::x,y,b,dxm,dxp,dym,dyp
Expand Down Expand Up @@ -689,7 +689,7 @@ subroutine viewr
idone=1
else
if (iskip.eq.0) then
if (l+5000.ge.maxaa) then
if (l+mmax.ge.maxaa) then
call mess('viewr','too much 3d data',&
'data truncated')
iskip=1
Expand Down Expand Up @@ -1270,16 +1270,17 @@ subroutine set3d(iplot,xyz,nxyz)
real(kr)::xyz(nxyz)
! internals
integer,parameter::length=2000
integer,parameter::lenxyz=10000
integer::i,j,nn,k,ncurv,major,minor,itop,ibot,l,n
real(kr)::wt,xlo,xhi,ylo,yhi,zlo,zhi,xn,yn,zn,yy
real(kr)::top,bot,xop,yop,zop
integer::lll(length)
real(kr)::x(2000),y(2000),z(2000)
real(kr)::x(lenxyz),y(lenxyz),z(lenxyz)
real(kr),parameter::big=1.e10_kr
real(kr),parameter::d0=.001e0_kr
real(kr),parameter::d3=.301e0_kr
real(kr),parameter::d7=.699e0_kr
integer::kmax=1999
integer::kmax=lenxyz-1 ! x(:),y(:),z(:) arrays dimension - 1
real(kr),parameter::zero=0
real(kr),parameter::one=1
real(kr),parameter::ten=10
Expand Down
25 changes: 20 additions & 5 deletions src/wimsr.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ module wimsm
integer,parameter::nymax=100
real(kr)::yield(nymax)
integer::ifisp(nymax)
integer::ifiss,nfiss
integer::ifiss,ifissr,nfiss
integer::isg
integer::ixs
real(kr),dimension(:),allocatable::snu
Expand Down Expand Up @@ -434,7 +434,7 @@ subroutine resint
integer::i,jg,is,jtemp,jfiss,nl,nz,ntw,l,ig,kg,lim,iadd
integer::loca,jz,loc,locn,jtem,iterm,ioff,it,index
integer::indexl,iz
real(kr)::xid,siglam,sigb,siga,sig
real(kr)::xid,siglam,sigb,siga,sigs,sig
real(kr),dimension(:),allocatable::sabs
real(kr),dimension(:),allocatable::snux
real(kr),dimension(:),allocatable::snsf
Expand All @@ -447,6 +447,7 @@ subroutine resint
real(kr),parameter::zero=0

!--allocate storage.
ifissr=0
if (ires.eq.0) go to 510
write(nsyso,'(/'' ***resonance integrals***'')')
ntsr=ires*nsigz*nrg
Expand Down Expand Up @@ -519,6 +520,7 @@ subroutine resint
if (mfh.ne.3) go to 300
if (mth.eq.1) go to 162
if (mth.eq.2) go to 162
if ((mth.gt.50).and.(mth.le.91)) go to 162
if (mth.lt.18.or.mth.gt.150) go to 300
if (mth.gt.21.and.mth.lt.102.and.mth.ne.38) go to 300
162 continue
Expand All @@ -537,7 +539,8 @@ subroutine resint
jg=kg-nfg
lim=nsigz
if (nsigz.gt.nz) lim=nz
if (mth.ne.1.and.mth.ne.2.and.mth.ne.18.and.mth.ne.102) go to 300
if (mth.ne.1.and.mth.ne.2.and.mth.ne.18.and.mth.ne.102 &
.and.((mth.le.50).or.(mth.gt.91))) go to 300

!--absorption
if (mth.eq.102) then
Expand Down Expand Up @@ -578,8 +581,8 @@ subroutine resint
sabs(iadd-jz+1)=sabs(iadd-jz+1)+scr(nl*jz+loca)
enddo

!--elastic
else if (mth.eq.2) then
!--elastic or inelastic scattering
else if ((mth.eq.2).or.((mth.gt.50).and.(mth.le.91))) then
iadd=nsigz+nsigz*(jtemp-1+ires*(jg-1))
loca=l+lz+nl*(nz-1)
do jz=1,lim
Expand Down Expand Up @@ -661,13 +664,16 @@ subroutine resint
sigb=sigz(iz)+siglam
siga=sabs(loc+iz-1)
sig=snux(loc+iz-1)
sigs=elas(loc+iz-1)
sabs(loc+iz-1)=sigb*siga/(sigb+siga)
snux(loc+iz-1)=sigb*sig/(sigb+siga)
if (iverw.ne.4) elas(loc+iz-1)=sigb*sigs/(sigb+siga)
enddo
enddo
enddo

!--write out results
ifissr=jfiss
call rsiout(xid,jfiss,fa,sabs,snux,flxr,sigz,elas)

!--resint is finished.
Expand Down Expand Up @@ -929,6 +935,7 @@ subroutine xsecs
call tpidio(ngendf,0,0,scr,nb,nw)
isg=0
if (sgref.lt.dilinf) isg=1
if (iverw.eq.5) isg=1
jtemp=0
i318=0
jfisd=0
Expand All @@ -944,11 +951,15 @@ subroutine xsecs
xxi=1+log(alf)*alf/(1-alf)

!--loop over temperatures
do i=1,ngnd
snu(i)=0
enddo
140 continue
do i=1,ngnd
csp1(i)=0
spot(i)=sigp
xi(i)=xxi
abs1(i)=0
abs2(i)=0
sf0(i)=0
sfi(i)=0
Expand All @@ -957,6 +968,9 @@ subroutine xsecs
chi(i)=0
scat(i)=0
sn2n(i)=0
ab0(i)=0
xtr(i)=0
sdp(i)=0
l1e(i)=ngnd
l2e(i)=1
l1(i)=ngnd
Expand Down Expand Up @@ -2040,6 +2054,7 @@ subroutine wimout(jcc)

!--write material identifier data
awt=awr*amassn
if((ifis.eq.3).and.(ifissr.ne.3)) ifis=2
write(nout,'(i6,1p,e15.8,5i6)')&
ident,awt,iznum,ifis,ntemp,nrestb,isof

Expand Down
Loading