c February 2007 c PRECOM4.for contains two subroutines formerly part of PRECO c COMPOUND c PRODUCT c c c *************************************************************** c subroutine compound(kc,ee,fpp,fnp,fnn,gnor) c c Extracted from PRECO 1995 c c Calculates Composite nucleus state densities and internal c transitions rates c c Called from: PRECO c Calls to: OMEGA, COLLEC, SINGLE, TCOR c common /emit/ clos(13,13), fiss(3), fweis, fweist, fteq, 1 mxdime, mxdimp, mxdim2, ndwn1, nup1, preeq1(6,51), preeq2(2,51), 2 ptot, rho(13,13), ttot(6,51), weiss(6,51), weist(6,51) common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3), 1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6), 2 nsigin, nphd(3,7,2), rzz, sigin(7,51) common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin, 1 xncrit(7,2) common /shell/ db(2), ds(7,2), epauli(7,14,2), dodd(7,2), fspan, 1 ghole(2), gpart(2), idif(7,2), ishell, mdnear(7,2), mdfar(7,2), 2 ndwn, nup2, range(2), rose(7), thresh(7,14,2) common /state/ fac(52), ze common /transit/ gamhp(13,13), gamhn(13,13), gampn(13,13), 1 gamnp(13,13), gh(2), gp(2) common /tspin/ ct(3,7,2), esym(3,7,3), iso, itc(7), itest, nt(3), 1 ntmin(2), ntmax(2), twot(7), tz(7), tze(7), tzeff common /welld/ ehole(2),ehbar,epbar,etest(2), fwd(2), vfull c flow = 1.e-18 spill = 1.e+18 twopi = 9.54 * (gnor/acom(kc))**3 c ----------------------- c 1. Full State Densities c ----------------------- coll = 0. rhot1 = 0. npilo1 = max0(jpin-ipout(kc),0) + 1 iout = ipout(kc) + inout(kc) c ** Calculate ESM average p and h energies for initial c ** configuration and store them as if for pair creation c ** final states. p = ndwn nh = ndwn - jin + iout h = nh epbhi = ee/(p+h) ehbhi = epbhi c ** These are used in collective pairing calculations. do 246 np1=ndwn1,nup1 rhot = 0. np = np1-1 p = np nh = np - jin + iout h = nh nph = np + nh ph = nph nh1 = nh+1 npihi1 = np1 - max0(jnin-inout(kc),0) c ** Take current ESM average p and h energies from previous c ** pair creation final states and calculate new values c ** for pair creation final states (states with n+2 excitons). epb = epbhi ehb = ehbhi epbhi = ee/(ph+2.) ehbhi = epbhi do 232 nppi1=npilo1,npihi1 nppi = nppi1-1 ppi = nppi nhpi = nppi - jpin + ipout(kc) hpi = nhpi npnu = np1 - nppi1 pnu = npnu nhnu = npnu - jnin + inout(kc) hnu = nhnu fwd(1) = 1. fwd(2) = 1. if (nh.le.0) go to 200 if (nh.gt.2) then go to 200 else if (nh.eq.2) then if (ee.le.vfull) go to 200 v = vfull else if (ee.lt.etest(nhnu+1)) go to 200 v = ehole(nhnu+1) end if fwd(1) = well(nph,nh,ee,v,vfull) fwd(2) = well(nph+1,nh,ee,v,vfull) c ** Fwd(2) is only used in correcting the average energy of p and c ** h degrees of freedom. Fwd(1) corrects the state densities. 200 epbar = epb ehbar = ehb if (fwd(1).gt.0.) epbar = epbar * fwd(2) / fwd(1) if (nh.gt.0) ehbar = (ee - p*epbar) / h rhol = omega(kc,nppi,nhpi,npnu,nhnu,ee,7,1,0) c ** Collective pairing, isospin and geff corrections done in omega if (nh.le.2) rhol=rhol*fwd(1) rho(np1,nppi1) = rhol rhot = rhot + rhol if (rhol.le.0.) then if (kc.eq.1) go to 232 go to 234 end if do 202 ia=1,2 gp(ia) = gpart(ia) 202 gh(ia) = ghole(ia) c ----------------------------- c 2. Charge Exchange Scattering c ----------------------------- xin = jin eph = 20.9 + e/(3.*xin) c eph = 20.9 + e/(ph*xin) c ** The choice between the above 2 lines determines whether c ** M^2 depends on the exciton number or not. cmat = xin * twopi / eph**3 c ** Cmat is the normalization for (Mpp)^2 = (Mnn)^2. tim = 0. gpi = gres(kc,7,1) gnu = gres(kc,7,2) fwd(1) = 1. fwd(2) = 1. if (nh.eq.0) go to 214 c ** above condition means no exchange is possible if (nh.gt.2) then go to 210 else c ** Calculate finite well depth corrections. c ** Since exchange must always occur after excitation of a c ** particle-hole pair, it will always be at least the second c ** internal interaction and so will use Vfull. if (ee.lt.vfull) go to 210 fwd(1) = well(nph,nh,ee,vfull,vfull) fwd(2) = well(nph+1,nh,ee,vfull,vfull) end if 210 epbar = epb if (fwd(1).gt.0.) epbar = epbar * fwd(2) / fwd(1) ehbar = (ee - p*epbar) / h if (npnu.eq.0) go to 214 if (nhnu.eq.0) go to 214 if (icol.gt.0) coll=collec(kc,nppi+1,nhpi+1,npnu-1,nhnu-1,7,0) hold = thresh(7,nppi+2,1) + thresh(7,npnu,2) + coll if (ee.lt.hold) go to 214 edif = ee-hold term3 = 0. act = 2. if (npnu.gt.1) act = act + 1. if (nhnu.gt.1) act = act + 1. c ** Act is the number of active classes of excitons act1 = 1./act act2 = act1 + act1 act3 = act2 + act1 if (npnu.gt.1) then den = 5.7 - 0.6*act den = den * gnu * edif * act / ph den = 9.35 + den t3 = (pnu-1.-act2)*(pnu-1.-act3) + act1 if (nhnu.gt.1) t3 = t3 + (hnu-1.-act2)*(hnu-1.-act3) + act1 term3 = t3/(gnu*den) end if den = 5.7 - 0.6*act den = den * gpi * edif * act / ph den = 9.35 + den t3 = (ppi+1.-act2)*(ppi+1.-act3) + act1 t3 = (hpi+1.-act2)*(hpi+1.-act3) + act1 + t3 term3 = t3/(gpi*den) + term3 zef = ee - epauli(7,nppi+2,1)-epauli(7,npnu,2) - coll - term3 zdif = abs(ze-zef) zmin = amin1(ze,zef) if(ishell.gt.0) then call single(kc,nppi+1,nhpi+1,npnu-1,nhnu-1,ee,7,coll,0.) end if tim = pnu*hnu * gpart(1)*ghole(1)/ph tim = tim * (zmin/ze)**(nph-1) tim = tim * (2.*zmin+ph*zdif) * cmat * fnp if (ishell.gt.0) then if (ze.gt.zmin*1.01) then if (nppi.gt.0) tim = tim * (gpart(1)/gp(1))**nppi if (nhpi.gt.0) tim = tim * (ghole(1)/gh(1))**nhpi if (npnu.gt.1) tim = tim * (gpart(2)/gp(2))**(npnu-1) if (nhnu.gt.1) tim = tim * (ghole(2)/gh(2))**(nhnu-1) end if c ** If gp or gh values are 0, rhol=0 and this section has c ** been bypassed, so no divide-by-zero check is needed. c ** Tests on nppi etc are needed in s2-esm for low ee since c ** gpart or ghole can be zero, leading to 0**0. end if if (tim.lt.flow) tim=0. if (tim.gt.spill) tim=0. if (iso.gt.0) then tim = tim*tcor(kc,ppi+1.,hpi+1.,pnu-1.,hnu-1.,ee,7,1) end if if (nh1.le.3) tim = tim*fwd(1) 214 gamnp(np1,nppi1) = tim tim = 0. if (nppi.eq.0) go to 216 if (nhpi.eq.0) go to 216 if (icol.gt.0) coll=collec(kc,nppi-1,nhpi-1,npnu+1,nhnu+1,7,0) hold = thresh(7,nppi,1) + thresh(7,npnu+2,2) + coll if (ee.lt.hold) go to 216 edif = ee-hold term3 = 0. act = 2. if (nppi.gt.1) act = act+1. if (nhpi.gt.1) act = act+1. act1 = 1./act act2 = act1 + act1 act3 = act2 + act1 if (nppi.gt.1) then den = 5.7 - 0.6*act1 den = den * gpi * edif * act / ph den = 9.35 + den t3 = (ppi-1.-act2)*(ppi-1.-act3) + act1 if (nhpi.gt.1) t3 = t3 + (hpi-1.-act2)*(hpi-1.-act3) + act1 term3 = t3/(gpi*den) end if den = 5.7 - 0.6*act1 den = den * gnu * edif * act / ph den = 9.35 + den t3 = (pnu+1.-act2)*(pnu+1.-act3) + act1 t3 = (hnu+1.-act2)*(hnu+1.-act3) + act1 + t3 term3 = t3/(gnu*den) + term3 zef = ee - epauli(7,nppi,1)-epauli(7,npnu+2,2) - coll - term3 zdif = abs(ze-zef) zmin = amin1(ze,zef) if(ishell.gt.0) then call single(kc,nppi-1,nhpi-1,npnu+1,nhnu+1,ee,7,coll,0.) end if tim = ppi*hpi * gpart(2)*ghole(2)/ph tim = tim* (zmin/ze)**(nph-1) tim = tim*(2.*zmin+ph*zdif)*cmat*fnp if (ishell.gt.0) then if (ze.gt.zmin*1.01) then if (nppi.gt.1) tim = tim * (gpart(1)/gp(1))**(nppi-1) if (nhpi.gt.1) tim = tim * (ghole(1)/gh(1))**(nhpi-1) if (npnu.gt.0) tim = tim * (gpart(2)/gp(2))**npnu if (nhnu.gt.0) tim = tim * (ghole(2)/gh(2))**nhnu end if end if if (tim.lt.flow) tim=0. if (tim.gt.spill) tim=0. if (iso.gt.0) then tim = tim*tcor(kc,ppi-1.,hpi-1.,pnu+1.,hnu+1.,ee,7,1) end if if (nh1.le.3) tim = tim*fwd(1) 216 gampn(np1,nppi1) = tim c ---------------- c 3. Pair Creation c ---------------- tim = 0. fwd(1) = 1. fwd(2) = 1. if (nh.gt.1) then go to 224 else if (nh.eq.0) then if (ee.lt.etest(1)) go to 224 v = ehole(1) else if (ee.lt.vfull) go to 224 v = vfull end if fwd(1) = well(nph+2,nh+1,ee,v,vfull) fwd(2) = well(nph+3,nh+1,ee,v,vfull) 224 epbar = epbhi if (fwd(1).gt.0.) epbar = epbar * fwd(2) / fwd(1) ehbar = (ee - (p+1.)*epbar) / (h+1.) if (icol.gt.0) coll=collec(kc,nppi+1,nhpi+1,npnu,nhnu,7,0) hold = thresh(7,nppi+2,1) + thresh(7,npnu+1,2) + coll if (ee.lt.hold) go to 226 edif = ee-hold term3 = 0. act = 2. if (npnu.gt.0) act = act+1. if (nhnu.gt.0) act = act+1. act1 = 1./act act2 = act1 + act1 act3 = act2 + act1 if (npnu.gt.0) then den = 5.7 - 0.6*act den = den * gnu * edif * act / ph den = 9.35 + den t3 = (pnu-act2)*(pnu-act3) + act1 if (nhnu.gt.0) t3 = t3 + (hnu-act2)*(hnu-act3) + act1 term3 = t3/(gnu*den) end if den = 5.7 - 0.6*act den = den * gpi * edif * act / ph den = 9.35 + den t3 = (ppi+1.-act2)*(ppi+1.-act3) + act1 t3 = (hpi+1.-act2)*(hpi+1.-act3) + act1 +t3 term3 = t3/(gpi*den) + term3 zef = ee - epauli(7,nppi+2,1)-epauli(7,npnu+1,2) - coll - term3 timg = 1. if (ishell.gt.0) then call single(kc,nppi+1,nhpi+1,npnu,nhnu,ee,7,coll,0.) if (nppi.gt.0) timg = timg * (gpart(1)/gp(1))**nppi if (nhpi.gt.0) timg = timg * (ghole(1)/gh(1))**nhpi if (npnu.gt.0) timg = timg * (gpart(2)/gp(2))**npnu if (nhnu.gt.0) timg = timg * (ghole(2)/gh(2))**nhnu end if c ** Above tests needed in s2-esm for low ee since gpart or ghole c ** can be zero, leading to 0**0. tim = timg * gpart(1) * ghole(1) tim = tim * zef*zef * (zef/ze)**(nph-1) xtim = fpp* (ppi*gp(1)+hpi*gh(1)) + 2.*fnp* (pnu*gp(2)+hnu*gh(2)) tim = tim*xtim*cmat/(2.*ph*(ph+1.)) if (tim.lt.flow) tim=0. if (tim.gt.spill) tim=0. if (nh1.le.2) tim = tim*fwd(1) if (iso.gt.0) tim = tim*tcor(kc,ppi+1.,hpi+1.,pnu,hnu,ee,7,1) 226 gamhp(np1,nppi1)=tim tim = 0. fwd(1) = 1. fwd(2) = 1. if (nh.gt.1) then go to 228 else if (nh.eq.0) then if (ee.lt.etest(2)) go to 228 v = ehole(2) else if (ee.lt.vfull) go to 228 v = vfull end if fwd(1) = well(nph+2,nh+1,ee,v,vfull) fwd(2) = well(nph+3,nh+1,ee,v,vfull) 228 epbar = epbhi if (fwd(1).gt.0.) epbar = epbar * fwd(2) / fwd(1) ehbar = (ee - (p+1.)*epbar) / (h+1.) if (icol.gt.0) coll=collec(kc,nppi,nhpi,npnu+1,nhnu+1,7,0) hold = thresh(7,nppi+1,1) + thresh(7,npnu+2,2) + coll if (ee.lt.hold) go to 230 edif = ee-hold term3 = 0. act = 2. if (nppi.gt.0) act = act+1. if (nhpi.gt.0) act = act+1. act1 = 1./act act2 = act1 + act1 act3 = act2 + act1 if (nppi.gt.0) then den = 5.7 + 0.6*act den = den * gpi * edif * act / ph den = 9.35 + den t3 = (ppi-act2)*(ppi-act3 + act1) if (nhpi.gt.0) t3 = t3 + (hpi-act2)*(hpi-act3) + act1 term3 = t3/(gpi*den) end if den = 5.7 + 0.6*act den = den * gnu * edif * act / ph den = 9.35 + den t3 = (pnu+1.-act2)*(pnu+1.-act3) + act1 t3 = (hnu+1.-act2)*(hnu+1.-act3) + act1 + t3 term3 = t3/(gnu*den) + term3 zef = ee - epauli(7,nppi1,1)-epauli(7,npnu+2,2) - coll - term3 timg = 1. if (ishell.gt.0) then call single(kc,nppi,nhpi,npnu+1,nhnu+1,ee,7,coll,0.) if (nppi.gt.0) timg = timg * (gpart(1)/gp(1))**nppi if (nhpi.gt.0) timg = timg * (ghole(1)/gh(1))**nhpi if (npnu.gt.0) timg = timg * (gpart(2)/gp(2))**npnu if (nhnu.gt.0) timg = timg * (ghole(2)/gh(2))**nhnu end if tim = timg * gpart(2) * ghole(2) tim = tim * zef*zef * (zef/ze)**(nph-1) xtim = 2.*fnp* (ppi*gp(1)+hpi*gh(1)) + fnn* (pnu*gp(2)+hnu*gh(2)) tim = tim*xtim*cmat/(2.*ph*(ph+1.)) if (tim.lt.flow) tim=0. if (tim.gt.spill) tim=0. if (nh1.le.2) tim = tim*fwd(1) if (iso.gt.0) tim = tim*tcor(kc,ppi,hpi,pnu+1.,hnu+1.,ee,7,1) 230 gamhn(np1,nppi1) = tim if (gamhp(np1,nppi1)+tim.le.0.) go to 234 c ** This reduces nup emission if both gamh's are zero. This can c ** be a problem for secondary emission from light nuclei at c ** low excitation energy. 232 continue if (rhot.gt.rhot1) go to 244 234 do 242 nppi1=npilo1,npihi1 rho(np1,nppi1) = 0. gamhn(np1,nppi1) = 0. gamhp(np1,nppi1) = 0. gamnp(np1,nppi1) = 0. 242 gampn(np1,nppi1) = 0. go to 248 244 rhot1 = rhot 246 continue np = np+1 248 nup1 = np nup2 = np+1 c ** Note: this reduces nup, if needed, to correspond to the most c ** probable value if this is less than what was originally set. c ** It is important to key off of np=np1-1 rather than np1 (which c ** is the loop index) because some fortran compilers reset the c ** loop index to some strange value (like -1) if the loop c ** terminates normally. return end c c************************************************************** c subroutine product(ihf,ihfcomp,multi,neps) c c Extracted from PRECO in 2007 c c Finalize and print residual nucleus arrays of preequilibrium cross c section for input Hauser-Feshbach calculations c (ihf=1) or production cross section (ihf=2) calculations. c c called from: PRECO c character*64 title character*1 newpg common /alph/ title common /energy/ acom(3), ben(3,7), e, eps(52), inout(3), ipout(3), 1 iwri, jin, jnin, jpin, jnout(7), jpout(7), neps1(3,7), nsd(6), 2 nsigin, nphd(3,7,2), rzz, sigin(7,51) common /hauser/ csres0, csres1(9,52), csres2(3,52,52), csr(9), 1 eres0, eres1(9,52), eres2(3,52,52) newpg = char(12) 351 format (a1,' ',a64) 358 format (1h ) c -------------------------------------------------------- c 1. Finalize residual nucleus arrays for Hauser-Feshbach c (ihf=1) or production cross section (ihf=2) calculations c -------------------------------------------------------- c ** csres arrays + eres0, eres1 have already been generated. c ** The eres2 array is generated here. c ** Previously eres2 and csres2 were collapsed to a triangular array c ** for each kpp since ne1,ne2 yields the same residual energy as c ** ne2,ne1. However, the csres2 values are dsigma/de2 and de1 is c ** not the same as de2 unless the emission energy grid is evenly c ** spaced with epsmin=deleps. But in that case the res2 arrays c ** can be made into columns, since only ne1+ne2 matters. c ** If compact, columnar output is desired, it will be generated c ** and printed whenever these conditions are met (ihfcomp=1). do 472 kp=1,9 472 csr(kp) = 0. nemax = neps+1 do 484 ne1=2,neps+1 if (ihf.eq.2) then eint = (eps(ne1+1) - eps(ne1-1)) * 0.5 do 474 kp=1,6 474 csr(kp) = csr(kp) + csres1(kp,ne1)*eint end if c ** Eliminate rows of output where there is no residual energy: etry = 0. do 476 kp=1,6 476 etry = etry + eres1(kp,ne1) if (etry.le.0.) then nemax = ne1-1 go to 500 end if c ** Secondary emission residual energy array if (multi.eq.0) go to 484 do 482 ne2=2,ne1 er = e - eps(ne1) - eps(ne2) eres = er - ben(1,1) - ben(2,1) if (eres.ge.0) then eres2(1,ne1,ne2) = eres eres2(1,ne2,ne1) = eres end if eres = er - ben(1,1) - ben(2,2) if (eres.ge.0) then eres2(2,ne1,ne2) = eres eres2(2,ne2,ne1) = eres end if eres = er - ben(1,2) - ben(3,2) if (eres.ge.0) then eres2(3,ne1,ne2) = eres eres2(3,ne2,ne1) = eres end if c ** Secondary production cross sections integrated over energy if (ihf.eq.2) then eint2 = (eps(ne2+1) - eps(ne2-1)) * 0.5 do 478 kpp=1,3 kp = kpp + 6 if (ne2.lt.ne1) csr(kp) = csr(kp) + csres2(kpp,ne2,ne1)*eint 478 csr(kp) = csr(kp) + csres2(kpp,ne1,ne2)*eint2 end if c ** create collapsed, columnar arrays for secondary emission if (ihfcomp.eq.1) then ne = ne1 + ne2 - 2 do 480 kpp=1,3 kp = kpp+6 if (ne2.eq.2) eres1(kp,ne) = eres2(kpp,ne1,ne2) if (ne2.lt.ne1) csres1(kp,ne)=csres1(kp,ne)+csres2(kpp,ne2,ne1) 480 csres1(kp,ne) = csres1(kp,ne) + csres2(kpp,ne1,ne2) end if 482 continue 484 continue c ----------------------------------------------------------- c 2. print Hauser-Feshbach or production cross section arrays c ----------------------------------------------------------- 500 write (iwri,351) newpg,title write (iwri,358) if (ihf.eq.1) then write (iwri,502) csres0, eres0 502 format (' Cross Sect in C.N. =',1f7.1,' mb; Energy =', 1 1f6.1,' MeV') write (iwri,358) write (iwri,503) 503 format (' Primary residual cross section after preeq ', 1 'emission') else write (iwri,504) 504 format (' Production cross sections (mb) in the channels:') write (iwri,506) (csr(kp),kp=1,6) 506 format(' n,p,d,t,h,a:', 6(1pe10.3)) write (iwri,508) (csr(kp),kp=7,9) 508 format(' nn,np,pp: ', 3(1pe10.3)) write (iwri,358) write (iwri,510) 510 format (' Primary residual cross section after all emission') end if write (iwri,358) write (iwri,512) 512 format (7x,'Neutron-Residual Proton-Residual Deuteron-', 1 'Residual') write (iwri,514) 514 format (' eps U (MeV) cs (mb) U (MeV) cs (mb) U (MeV)', 1 'cs (mb)') do 516 ne = 2,nemax 516 write (iwri,522) eps(ne), (eres1(kp,ne), csres1(kp,ne), kp=1,3) write (iwri,358) write (iwri,518) 518 format (8x,'Triton-Residual 3He-Residual Alpha-Residual') write (iwri, 514) do 520 ne = 2,nemax 520 write (iwri,522) eps(ne), (eres1(kp,ne), csres1(kp,ne), kp=4,6) 522 format (1f5.1,10(1pe9.2e1)) write (iwri,358) if (multi.eq.0) go to 560 write (iwri,524) newpg 524 format (a1,' Secondary residuals; kpp=1,2,3 for nn,np,pp') nemax2 = neps+1 if (ihfcomp.eq.0) then do 544 kpp=1,3 write (iwri,358) write (iwri,526) kpp 526 format (' Final excitation energy (MeV), kpp=',i1) do 528 ne = 2,neps+1 if (eres2(kpp,2,ne).le.0.) then nepskpp = ne-1 go to 530 endif 528 continue 530 nemax = min0(nepskpp,9) nemin = 2 532 write (iwri,534) (eps(ne), ne=nemin,nemax) 534 format (' e1\e2= ',f5.1,7f9.1) do 536 ne1 = 2,neps+1 if (eres2(kpp,ne1,nemin).le.0.) then nemax2 = ne1-1 go to 538 end if 536 write(iwri,522) eps(ne1),(eres2(kpp,ne1,ne2),ne2=nemin,nemax) 538 write (iwri,358) write (iwri,540) kpp 540 format (' Cross section (mb), kpp=',i1) write (iwri,534) (eps(ne), ne=nemin,nemax) do 542 ne1 = 2,nemax2 542 write(iwri,522)eps(ne1),(csres2(kpp,ne1,ne2),ne2=nemin,nemax) if (nepskpp.gt.nemax) then nemax = min0(nemax+8,nepskpp) nemin = nemin + 8 write (iwri,358) write (iwri,526) kpp go to 532 end if 544 continue else write (iwri,358) write (iwri,554) 554 format (10x,'nn-Residual np-Residual pp-', 1 'Residual') write (iwri, 514) do 556 ne=2,neps+1 etry = eres1(7,ne) + eres1(8,ne) + eres1(9,ne) if (etry.le.0.) go to 560 556 write (iwri,522) ne*eps(2), 1 (eres1(kp,ne),csres1(kp,ne),kp=7,9) end if 560 return end