*ident up_groupr */-------------------------------------------------------------- */ correction in getsed to include the cases */ of noncartesian interpolation laws */ 18.04.06/V.Sinitsa/ENEA-Bologna */-------------------------------------------------------------- *d groupr.8647 external intega,releas,findex,terpa,anased,terp1,terp2,mess *d groupr.8754,8824 c ***read spectrum for each incident energy do while (nne.lt.ne) call tab1io(nin,0,0,c(l),nb,nw) l=l+nw do while (nb.ne.0) if (l.gt.nc) call error('getsed', & 'storage exceeded',' ') call moreio(nin,0,0,c(l),nb,nw) l=l+nw enddo nne=nne+1 enddo c c ***determine econst lk=loc(ik)+ic nr=nint(c(lk+4)) np=nint(c(lk+5)) lk=lk+6+2*(nr+np) nr=nint(c(lk+4)) ne=nint(c(lk+5)) ir=1 nbt=nint(c(lk+4+2*ir)) int=nint(c(lk+5+2*ir)) lhi=lk+6+2*nr do i=2,ne if (i.gt.nbt) then ir=ir+1 nbt=nint(c(lk+4+2*ir)) int=nint(c(lk+5+2*ir)) endif llo=lhi elo=c(llo+1) nrl=nint(c(llo+4)) npl=nint(c(llo+5)) lhi=llo+6+2*(nrl+npl) ehi=c(lhi+1) ipl=2 irl=1 iph=2 irh=1 do ig=1,ng e1=eg(ig) if (ig.eq.1) e1=0. e2=eg(ig+1) if (ig.eq.ng) e2=ebig call intega(clo,e1,e2,c(llo),ipl,irl) call intega(chi,e1,e2,c(lhi),iph,irh) delta=abs(chi-clo) if (delta.gt.eps*clo) then if (int.eq.1.and.ehi.lt.econst) econst=ehi if (int.ne.1) then ee=elo+eps*clo*(ehi-elo)/delta if (ee.lt.econst) econst=ee endif endif enddo enddo endif enddo c c ***initialization complete nc=l-ic econ=0 call releas('sc',nc,c) if (lf.eq.1) call reserv('tab2',5000,lint,c) nupm=0 return endif *i groupr.8829 if (lf.eq.1) call findex('tab2',lint,c) *d groupr.8878,8907 nne=0 nr=nint(c(mnow+4)) ne=nint(c(mnow+5)) nnt=6 nbt=nint(c(mnow+nnt)) int=nint(c(mnow+nnt+1)) nnow=mnow+6+2*nr lhi=nnow ehi=0 do while (nne.lt.ne.and.ed.gt.ehi*(1+small).and. & iout.eq.0) nne=nne+1 if (nne.gt.ne) then iout=1 else if (nne.gt.nbt) then nnt=nnt+2 nbt=nint(c(mnow+nnt)) int=nint(c(mnow+nnt+1)) endif llo=lhi elo=c(llo+1) nrl=nint(c(llo+4)) npl=nint(c(llo+5)) lhi=llo+6+2*(nrl+npl) ehi=c(lhi+1) endif enddo if (iout.eq.0) then ipl=2 iph=2 irl=1 irh=1 if (int.gt.10) & call terp2(c(llo),c(lhi),ed,c(lint),int) do ig=1,ng e1=eg(ig) if (ig.eq.1) e1=0. e2=eg(ig+1) if (ig.eq.ng) e2=ebig if (int.lt.10) then call intega(clo,e1,e2,c(llo),ipl,irl) call intega(chi,e1,e2,c(lhi),iph,irh) call terp1(elo,clo,ehi,chi,ed,s,int) else call intega(s,e1,e2,c(lint),ipl,irl) end if sed(ikt,ig)=sed(ikt,ig)+s*pe enddo *i groupr.8966 c subroutine terp2(tabl,tabh,e,tab,int) c ****************************************************************** c terp2: two-dimensional interpolation c input: c tabl,tabh - two consecutive tab1 records, c surrounding e-point c e - incident neutron energy point c int - interpolation law c output: c tab - tab1 record interpolated to e-point c ***************************************************************** *if sw implicit real*8 (a-h,o-z) *endif common/mainio/nsysi,nsyso,nsyse,ntty dimension tabl(2,*),tabh(2,*),tab(2,*) external terp1 data c1/0.999999/ isch=int/10 ints=int-isch*10 inte=ints el=tabl(2,1) eh=tabh(2,1) tab(1,1)=0. tab(2,1)=e tab(1,2)=0. tab(2,2)=0. tab(2,3)=0. nrl=nint(tabl(1,3)) nel=nint(tabl(2,3)) nrh=nint(tabh(1,3)) neh=nint(tabh(2,3)) c ---------------------------------------------------- c preliminary definition of the interpolation laws ne=max(nel,neh) if (nrh.ge.nrl) then nr=nrh do ir=1,nr tab(1,3+ir)=tabh(1,3+ir) tab(2,3+ir)=tabh(2,3+ir) end do else nr=nrl do ir=1,nr tab(1,3+ir)=tabl(1,3+ir) tab(2,3+ir)=tabl(2,3+ir) end do end if me=3+nr mel=3+nrl meh=3+nrh if (inte.eq.3) then inte=5 else if (inte.eq.4) then inte=2 end if c if (isch.eq.1) then c ------------------------ method of corresponding points do i=1,ne il=mel+min(i,nel) ih=meh+min(i,neh) ii=me+i call terp1(el,tabl(1,il),eh,tabh(1,ih),e,tab(1,ii),inte) call terp1(el,tabl(2,il),eh,tabh(2,ih),e,tab(2,ii),ints) end do c --- tab(1,3)=nr tab(2,3)=ne tab(1,3+nr)=ne return end if c --------- cartesian and unit base interpolation c ------ define scale factors rl,rh,r if (isch.eq.0) then c ---------------------- cartesian interplolation rl=1. rh=1. r=1. else c ----------------------- unit base interpolation rl=tabl(1,mel+nel) rh=tabh(1,meh+neh) call terp1(el,rl,eh,rh,e,r,inte) end if ne=nint(tab(2,3)) if (ne.eq.0) then iel=1 ieh=1 ei=0. do while (iel.lt.nel.or.ieh.lt.neh) eli=tabl(1,mel+iel)/rl ehi=tabh(1,meh+ieh)/rh if (c1*eli.lt.ehi) then if (iel.lt.nel) then ei=eli if (eli.ge.c1*ehi) ieh=min(ieh+1,neh) iel=min(iel+1,nel) else ei=ehi ieh=min(ieh+1,neh) end if else if (ieh.lt.neh) then ei=ehi if (ehi.ge.c1*eli) iel=min(iel+1,nel) ieh=min(ieh+1,neh) else ei=eli iel=min(iel+1,nel) end if end if ne=ne+1 tab(1,me+ne)=ei tab(2,me+ne)=0. end do ne=ne+1 tab(1,me+ne)=1. tab(2,me+ne)=0. end if c ---------------------------- elj=tabl(1,4+nrl)/rl ehj=tabh(1,4+nrh)/rh jl=1 jh=1 npl=3+nrl nph=3+nrh np=3+nr irl=1 krl=nint(tabl(1,4)) intl=nint(tabl(2,4)) irh=1 krh=nint(tabh(1,4)) inth=nint(tabh(2,4)) do ie=1,ne ei=tab(1,np+ie) do while (jl.eq.1.or.(ei.gt.c1*elj.and.jl.lt.nel)) il=jl jl=jl+1 eli=elj elj=tabl(1,npl+jl)/rl if (jl.gt.krl) then if (nrl.ge.nrh) tab(1,3+irl)=ie irl=irl+1 krl=nint(tabl(1,3+irl)) intl=nint(tabl(2,3+irl)) end if end do pli=0. if (ei.gt.c1*eli.and.c1*ei.lt.elj) & call terp1(eli,rl*tabl(2,npl+il),elj,rl*tabl(2,npl+jl), & ei,pli,intl) do while (jh.eq.1.or.(ei.gt.c1*ehj.and.jh.lt.neh)) ih=jh jh=jh+1 ehi=ehj ehj=tabh(1,nph+jh)/rh if (jh.gt.krh) then if (nrh.ge.nrl) tab(1,3+irh)=ie irh=irh+1 krh=nint(tabh(1,3+irh)) inth=nint(tabh(2,3+irh)) end if end do phi=0. if (ei.gt.c1*ehi.and.c1*ei.lt.ehj) & call terp1(ehi,rh*tabh(2,nph+ih),ehj,rh*tabh(2,nph+jh), & ei,phi,inth) c ------ if (ie.eq.ne) then pli=tabl(2,npl+jl) phi=tabh(2,nph+jh) end if call terp1(el,pli,eh,phi,e,pei,ints) tab(1,np+ie)=ei*r tab(2,np+ie)=pei/r end do tab(1,3)=nr tab(2,3)=ne tab(1,3+nr)=ne c write(nsyso,*) ' tab: e=',e c write(nsyso,'(1p6e12.5)') (tab(1,i),tab(2,i),i=1,3+nr+ne) return end */ --------------------------------------------------end up_groupr