c October 2006 c PRECOM.for - Main program file for PRECO-2006 c program preco c c **************************************************************** c * PRECO-2006 * c * by Constance Kalbach Walker * c * Replaces PRECO-2000 released in 2001 * c * * c * o Description of reactions with complex particle channels * c * significantly revised. * c * o Default assumptions with regard to isospin conservation * c * added. c **************************************************************** c c Closed form preequilibrium reaction calculations c in the two-component exciton model formalism c c Shell effects included thru shell-shifted equispacing model c Collective preequilbrium pairing option c Isospin conservation option with new default options c Direct excitation of collective states a la Kalka '89 c Revised surface effects defaults, June 1999 c c Subroutines for phenomenological descriptions of c 1. Direct knockout and inelastic scattering with complex c particle clusters c 2. Direct nucleon transfer c 3. Angular distributions using exponentials in cos(theta) c c Direct reaction formalism has been verified for nucleon, c deuteron, alpha and some He-3 induced reactions. c Projectile breakup is, however, not yet included, so breakup c channels will be underestimated. c (no direct reaction calculations for heavier projectiles) c (no default inverse or reaction cross sections for emitted c particles with A>4) c c calls to: ANGEL, COMPOUND, CROSS, EMMY, INEL, KNOCK, NUTRA, PAULI, c PRODUCT, REPENT, SETISO, SIGPAR, SPEC, SPECOUT c * * * * * * * * * * c * PRAISE THE LORD * c * * * * * * * * * * character*64 title character*1 newpg common kprint(6), mlo(2), mhi(2), p1(14), ref(2), 1 tau2(13), tau(13,13) common /alph/ title common /angels/ elastf(51), elcoll(51), kang, tomsc(51), 1 tomsd(51), xs(6) common /direct/ betalam(10), coul(7), elam(10), elast(51), fracel, 1 lam(10), lamax, scoll(51), sigbar(7), sigcni, siggeom(51), 2 snock(6,51), snutra(6,51), wcoll(10) 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 /hauser/ csres0, csres1(9,52), csres2(3,52,52), csr(9), 1 eres0, eres1(9,52), eres2(3,52,52) common /pair/ ceq(3,7,2), cpre(3,7,2), gres(3,7,2), icol, kin, 1 xncrit(7,2) common /par/ xl0(6), xl1(6), xm0(6), xm1(6) ,xn0(6), xn1(6), 1 xn2(6), xp0(6), xp1(6), xp2(6) common /repeat/ bfis(3), ecut, ifiss, irecyc, kpmax, nemin, 1 neprev, nepsin, einit, rgamma, rzcn(3), rncn(3) 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 common /wmain/ w(3,13,14,51), strong(9,10,51), seq(6,51), 1 seq2(2,51), steq(6,51), steq2(2,51), xiso(9,10,51), xisoeq(2,51), 2 spre(2,51) c c * * * * * * * * * * c * INITIALIZATION * c * * * * * * * * * * c irea=5 iwri=6 open(irea,file='precom.dat') open(iwri,file='precom.out',status='new') newpg=char(12) length=60 c ** length is the number of lines per page of output. c ** newpg forces a new page of output. Currently this is done c ** with an ASCII FF (form feed). mxdimp=13 mxdim2=9 mxdime=51 c ** mxdim's are the dimensions of the particle number (primary and c ** secondary emission) and energy indices c mxdimp.le.14 or fac needs to be changed to real*8 so it c doesn't overflow. c mxdimp.le.24 or the dimension of fac needs to be increased. c If mxdim's are changed, the common blocks need to be changed c along with the dimension statements in ANGEL and EMMY. c --------------------------------------------------------- c 1. Read in global parameters used for all problems in the c input file and set their default values. c --------------------------------------------------------- read (irea,52) ihfp, vfull, ehole1, ehole2 ihf = iabs(ihfp) c ** ihfp=0 specifies that PRECO is being used alone. This c ** bypasses the generation of arrays of cross section and c ** energy for the various residual nuclei populated. c ** If PRECO is to be incorporated into a Hauser-Feshbach model c ** code then set ihf=iabs(ihfp)=1 c ** to generate arrays of residual nuclei at various excitation c ** energies populated by preequilibrium emission. These are c ** used as input for evaporation calc'ns. c ** If ihf=2, then the arrays calculated include c ** equilibrium emission and give production cross sections at c ** different final excitation energies. c ** All these arrays are stored in common block /hauser/. c ** A negative value of ihfp means that compact arrays will be c ** printed for the secondary residual nuclei whenever an equally c ** spaced emeission energy grid is used and the minimum energy c ** is equal to the energy spacing. This combines all cross c ** section going to a given nucleus at a given excitation energy c ** regardless of the energies of the individual emitted particles c ** c ** vfull is the full potential well depth and defaults to 38 c ** MeV. c ** ehole1 and ehole2 are effective well depths after 1st c ** interaction, involving the creation of a proton(1) or c ** neutron(2) particle-hole pair. c ** If ehole1=0, systematic values with surface effects are used, c ** and vfull should be read in as 38 MeV or 0 (defaulting to 38). c ** Calculations currently start with h=0, p=jin. if (vfull.le.0.) vfull=38. if (ehole1.gt.vfull) ehole1=vfull if (ehole2.gt.vfull) ehole2=vfull read (irea,40) fpp,fnp,fnn, gnor c ** fpp, fnp, fnn are M-square normalizations. c ** gnor gives default s.p. state densities: gpi=Z/gnor; gnu=N/gnor c ** It is also used in calculating the washout of shell structure c ** with N or Z even when gpi and gnu are read in. c ** Finally it is used in the normalization of the mean square M.E. read (irea,40) frange, fspan, rgamma c ** frange multiplies D/2d to determine the range of N or Z for c ** the washout of shell structure effects. c ** fspan*t is the span for averaging over S2-ESM s.p. states to c ** get the equilibrium level density parameter. The averagying c ** interval is centered in the middle of the shell gap. c ** rgamma is the requirement on sigp/sign for secondary emission if (fpp.le.0.) fpp=2200 if (fnp.le.0.) fnp=900 if (fnn.le.0.) fnn=900 if (gnor.le.0) gnor=15. if (frange.le.0) frange = 2.0 if (fspan.le.0) fspan = 2.3 if (rgamma.le.0.) rgamma = 0.005 c ----------------------------------------------------------- c 2. Set numbers to be checked against for overflow/underflow c protection. c ----------------------------------------------------------- flow = 1.e-18 spill = 1.e+18 c ------------------------------------------------------------ c 3. Set up data for emitted particle types c (kp=a+z of ejectile for kp=1 to 6. exit channel label) c (kp=7 is for composite nucleus / entrance channel) c ------------------------------------------------------------ jpout(1) = 0 jpout(2) = 1 jpout(3) = 1 jpout(4) = 1 jpout(5) = 2 jpout(6) = 2 jpout(7) = 0 do 2 kp=1,6 jnout(kp) = kp - 2*jpout(kp) 2 nsd(kp) = 2 jnout(7) = 0 nsd(3) = 3 nsd(6) = 1 xs(1) = 0. xs(2) = 0. xs(3) = 2.22 xs(4) = 8.48 xs(5) = 7.72 xs(6) = 28.30 c ** These are the Z, N, (2s+1) and internal binding energy of c ** the particles in the kp channels. c ** ipout(kc) and inout(kc), calculated below, give the Z and N c ** for the previously emitted particle (if any) leading to c ** the nucleus from which preequilibrium emission is being c ** considered. kc=1,2,3 for initial CN, n resid., and p resid. ipout(1) = 0 ipout(2) = 0 ipout(3) = 1 inout(1) = 0 inout(2) = 1 inout(3) = 0 c ---------------------------------------------------------- c 4. Store parameters for approximate inverse cross sections c ---------------------------------------------------------- call sigpar c ------------------------------- c 5. Generate table of factorials c ------------------------------- fac(1)=1. ne=2*mxdimp+3 do 10 n=1,ne x=n 10 fac(n+1)=fac(n)*x c ------------------------------------------------------- c 6. Zero arrays not fully initialized with finite values c ------------------------------------------------------- 20 do 23 ne=1,mxdime eps(ne) = 0. sigin(7,ne) = 0. tomsd(ne) = 0. tomsc(ne) = 0. elast(ne) = 0. elastf(ne) = 0. elcoll(ne) = 0. scoll(ne) = 0. do 21 kp=1,2 preeq2(kp,ne) = 0. spre(kp,ne) = 0. seq2(kp,ne) = 0. 21 steq2(kp,ne) = 0. do 22 kp=1,6 preeq1(kp,ne) = 0. ttot(kp,ne) = 0. seq(kp,ne) = 0. steq(kp,ne) = 0. weiss(kp,ne) = 0. weist(kp,ne) = 0. snutra(kp,ne) = 0. snock(kp,ne) = 0. csres1(kp,ne) = 0. eres1(kp,ne) = 0. 22 sigin(kp,ne) = 0. do 23 kpp=1,3 csres1(kpp+6,ne) = 0. eres1(kpp+6,ne) = 0. do 23 ne2 = 1,mxdime csres2(kpp,ne,ne2) = 0. 23 eres2(kpp,ne,ne2) = 0. csres0 = 0. eps(mxdime+1) = 0. do 24 np=1,mxdimp tau2(np) = 0. p1(np) = 0. do 24 npp=1,mxdimp tau(np,npp) = 0. gamhp(np,npp) = 0. gamhn(np,npp) = 0. gampn(np,npp) = 0. gamnp(np,npp) = 0. rho(np,npp) = 0. 24 clos(np,npp) = 0. fwd(1) = 0. fwd(2) = 0. do 32 kp=1,7 rose(kp) = 0. twot(kp)=0. tz(kp) = 0. do 28 kc=1,3 ct(kc,kp,1) = 1. ct(kc,kp,2) = 0. 28 esym(kc,kp,3) = 0. do 30 ia=1,2 do 29 kc=1,3 esym(kc,kp,ia) = 0. 29 gres(kc,kp,ia) = 0. dodd(kp,ia) = 0. mdnear(kp,ia) = 0. 30 mdfar(kp,ia) = 0. do 32 np=1,mxdimp thresh(kp,np,1) = 0. thresh(kp,np,2) = 0. epauli(kp,np,1) = 0. 32 epauli(kp,np,2) = 0. do 36 ne=1,mxdime do 34 npp=1,mxdimp+1 do 34 np=1,mxdimp do 34 kp=1,3 34 w(kp,np,npp,ne) = 0. c ** W should be w(6,mxdimp,mxdimp,mxdime) but nearly half of grid c ** would be unused. Thus the array is compacted to c ** w(3,mxdimp,mxdimp+1,mxdime). c ** For kp.gt.3 dummy variables are used: kp goes to mdp = kp-3 c np1 goes to mp1 = mxdimp-np1+1 c nppi1 goes to mppi1 = mxdimp-nppi1+2 c ** A similar scheme is used for the arrays strong and xiso below xisoeq(1,ne) = 0. xisoeq(2,ne) = 0. do 36 npp=1,mxdim2+1 do 36 np=1,mxdim2 strong(np,npp,ne) = 0. 36 xiso(np,npp,ne) = 0. ntmin(1) = 0 ntmin(2) = 0 ntmax(1) = 0 ntmax(2) = 0 do 38 kc = 1,2 fiss(kc) = 0. 38 bfis(kc) = 0. c c * * * * * * * * * * * * * * * * * * * * c * INPUT and Preliminary Calculations * c * * * * * * * * * * * * * * * * * * * * c c ------------------------- c 1. Reaction specification c ------------------------- read (irea,41) title 41 format(1a64) write (*,43) title 43 format (1x,1a64) read (irea,40) e, ben(1,7), c, fteq einit = e if (ihf.gt.0) eres0 = e iso = c + 0.001 c ** iso=0 means t is mixed. c ** iso=1 means t is conserved, default symmetry energies. c ** iso>1 means use default cond'ns for t conservation and fteq. c ** fteq is the fraction of isospin conservation at equilibrium. read (irea,40) rzz,rnn, rzmag,rnmag, db(1),db(2) 40 format(7f10.2) c ** Z and N of target; Z and N of closed shells; p and n shell gaps c ** If rzmag < 0, then default shell parameters are used. c ** These are generated after projectile information is read in. c ** If rzmag and/or rnmag = 0, no shell effects are used for c ** that type of nucleon. atar = rzz + rnn jrz = rzz + 0.001 jrn = rnn + 0.001 ishell = 0 c ** ishell=0 for no shell effects; ishell=1 if shells considered. if (rzmag + rnmag.gt.1) then ishell = 1 if (rzmag.le.0.) db(1)=0. if (rnmag.le.0.) db(2)=0. read (irea,42) mlo(1), mhi(1), mlo(2), mhi(2) 42 format(bn,7i5) c ** mlo and mhi are degeneracies of piled up levels in S2-ESM. c ** These are later stored in arrays mdnear and mdfar. c ** Default is S2-ESM values which are generated in c ** subroutine pauli. Shell model values can be input. end if read (irea,42) jpin,jnin, multi c ** Z and N of projectile and control variable for 2ndary emission. c ** If multi=1, then secondary nucleon emission is calculated. jzcn = jrz + jpin jncn = jrn + jnin c ** Now set default values of rzmag, rnmag, and shell gaps. if (rzmag.lt.-0.001) then mlo(1) = 0 mlo(2) = 0 mhi(1) = 0 mhi(2) = 0 db(1) = 0. db(2) = 0. rzmag = 0. rnmag = 0. if (jzcn.le.12) then if (jzcn.ge.6) then rzmag = 8. db(1) = 2.5 end if else if (jzcn.le.24) then if (jzcn.ge.17) then rzmag = 20. db(1) = 1.6 end if else if (jzcn.le.33) then rzmag = 28. db(1) = 1.5 else if (jzcn.le.59) then if (jzcn.ge.43) then rzmag = 50. db(1) = 2.0 end if else if (jzcn.le.92) then if (jzcn.ge.74) then rzmag = 82. db(1) = 1.4 end if end if if (jncn.le.12) then if (jncn.ge.6) then rnmag = 8. db(2) = 3.1 end if else if (jncn.le.25) then if (jncn.ge.16) then rnmag = 20. db(2) = 2.4 end if else if (jncn.le.33) then rnmag = 28. db(2) = 1.3 else if (jncn.le.60) then if (jncn.ge.42) then rnmag = 50. db(2) = 2.3 end if else if (jncn.le.95) then if (jncn.ge.71) then rnmag = 82. db(2) = 1.9 end if else if (jncn.le.142) then if (jncn.ge.112) then rnmag = 126. db(2) = 1.6 end if end if if (rzmag+rnmag.gt.0.001) ishell=1 end if if (ishell.gt.0) then range(1) = frange * db(1) * rzmag / (gnor+gnor) range(2) = frange * db(2) * rnmag / (gnor+gnor) c ** Range(ia) is the range of Z or N over which shell effects c ** disappear. end if jin = jpin+jnin kin = jin+jpin if (kin.le.0) kin = 7 c ** This would correspond to an indident gamma ray, but the model c ** has never been studied for this case, and there are no default c ** inverse cross sections. xin = jin xpin = jpin xnin = jnin np0 = max0(jin,1) nh0 = np0 - jin acom(1) = atar + xin acom(2) = atar + xin -1. acom(3) = acom(2) elab = (e - ben(1,7)) * atar / (atar+xin) kc = 1 kpmax = 6 read (irea,40) (ben(1,kp),kp=1,6), bfis(1) c ** binding energies of emitted particles in original composite c ** nucleus and fission barrier ifiss = 0 if (bfis(1).gt.0.) ifiss = 1 if(multi.gt.0) then read (irea,40) ben(2,1), ben(2,2), ben(3,1), ben(3,2), 1 bfis(2), bfis(3) c ** neutron and proton binding energies and fission barriers c ** for secondary emission if (ifiss.gt.0) then if (bfis(2).le.0.) bfis(2) = bfis(1) if (bfis(3).le.0.) bfis(3) = bfis(1) end if end if read (irea,42) kang, (kprint(kp),kp=1,6) c ** If kang>0 the double differential cross sections will c ** be printed. Otherwise, only the energy differential cross c ** sections are included in the output. c ** If kprint(kp)=1 results for kp emission are printed c ** The default (if all kprints = 0) is to print all. ktot = 0 do 44 kp=1,6 44 ktot = ktot + kprint(kp) if (ktot.le.0) then do 46 kp=1,6 46 kprint(kp) = 1 end if c ---------------------------------------- c 2. Collective Resonance State Parameters c ---------------------------------------- do 47 i=1,10 lamax= i read (irea,52) lam(i), elam(i), betalam(i), wcoll(i) c ** Multipolarity, excitation energy, deformation parameter c ** and gaussian width parameter or the ith collective state. c ** If the multipolarity of a state is negative, then a Lorentzian c ** line shape is used instead of a Gaussian. c ** A negative width parameter signals that this is the last c ** collective state to be read in. if (wcoll(i).lt.0.) then wcoll(i) = -wcoll(i) go to 48 else if (wcoll(i).eq.0.) then if (i.eq.1) then wcoll(i) = 1. else wcoll(i) = wcoll(i-1) end if end if 47 continue c ** NOTE: wcoll is the width parameter of a gaussian c ** (or lorentzian) and should be no smaller than c ** about half the spacing of the emission energy grid. c ------------------------------------------------------------ c 3. Composite Nucleus Formation Cross Section; c CN and geometric cross sections for collective resonances c ------------------------------------------------------------ 48 read (irea,40) sigcni c ** zero invokes defaults and is recommended unless jin=0 or jin>4 if (jin.gt.0) then if (jin.lt.5) then neps1(1,kin) = 2 eps(2) = e - ben(1,7) call cross(1,kin) if (sigcni.eq.0.) sigcni = sigin(kin,2) siggeom(1) = siggeom(2) end if end if c ** Subroutine cross is called even if the default inverse cross c ** section isn't used, because cross also calculates the c ** corresponding geometric cross section, which is needed in the c ** collective model calculations. c ** Geometric cross sections at resonances are stored inside cross. c ----------------------------------------------- c 4. Emission Energies and Inverse Cross Sections c ----------------------------------------------- read (irea,52) neps,eps(2),deleps 52 format(bn,1i5,4f10.2) c ** Number of emission energies to be considered; lowest energy c ** and energy increment (if default inverse cross sections are c ** used). c ** If eps(2)<0 then energies and inverse cross sections are read. c ** The energy eps(1) is used as the lower limit in energy c ** "integrations" while eps(neps+2) is the upper limit. if (neps.ge.mxdime) then write (iwri,53) mxdime-1 53 format(' too many energies; maximum is',1i4) go to 1000 end if if (neps.le.0) neps=2 nepsin = neps ihfcomp = 0 if (ihfp.lt.0.) ihfcomp=1 c ** ihfcomp=1 specifies compact printout of the H-F or c ** production cross section arrays, if this is possible. c ** Below the code checks for the necessary conditions: c ** an equally spaced emission energy grid with the minimum c ** energy equal to the energy spacing. if (eps(2).le.0.) then nsigin = 1 c ** Parameter identifying the source of sigin: c ** nsigin=1 implies values are to be read in c ** nsigin=0 implies default values read (irea,40) eps(2),(sigin(kp,2),kp=1,6) do 54 ne=3,neps+1 read (irea,40) eps(ne),(sigin(kp,ne),kp=1,6) if (eps(ne)-eps(ne-1).ne.eps(2)) ihfcomp=0 54 continue do 57 kp=1,6 emu = e - ben(1,kp) if (kp.eq.kin) emu = emu + 4.*wcoll(1) do 56 ne=2,neps+1 ner = ne if (eps(ne).gt.emu) go to 57 56 continue ner = ner+1 57 neps1(1,kp) = ner-1 c ** note: ner is introduced to avoid problems with compilers c ** that reset loop indices to some strange value after normal c ** loop termination. else nsigin = 0 if (deleps.le.0.) deleps=1. if (eps(2).ne.deleps) ihfcomp=0 do 58 ne=3,neps+1 58 eps(ne) = eps(ne-1) + deleps do 59 kp=1,6 xne = (e-ben(1,kp)-eps(2)) / deleps + 1. if (kp.eq.kin) xne = xne + 4.*wcoll(1)/deleps ne = xne neps1(1,kp) = min0(ne,neps) + 1 59 call cross(1,kp) end if eps(1) = 0. eps(neps+2) = eps(neps+1)*2. - eps(neps) if (multi.eq.0) ihfcomp=0 c -------------------------------------------- c 5. Single Particle State Densities; gpi, gnu c -------------------------------------------- read (irea,40) gpi,gnu,gpfis,gnfis c ** These are for the composite nucleus in its g.s. and at top c ** of fission barrier. c ** Zero values invoke default of Z/gnor and N/gnor c ** Now set Z and N of composite nucleus. rzcn(1) = rzz + xpin rncn(1) = rnn + xnin if (gpi.gt.0.) then c ** Option to read in s.p. state densities for residual nuclei c ** Blank lines causes them to scale with Z and N. read (irea,40) (gres(1,kp,1),kp=1,6) read (irea,40) (gres(1,kp,2),kp=1,6) else gpi = rzcn(1)/gnor gnu = rncn(1)/gnor end if if (gres(1,1,1).le.0.) then do 64 kp=1,6 xpout = jpout(kp) xnout = jnout(kp) gres(1,kp,1) = gpi * (rzcn(1)-xpout) / rzcn(1) 64 gres(1,kp,2) = gnu * (rncn(1)-xnout) / rncn(1) end if gres(1,7,1) = gpi gres(1,7,2) = gnu if (gpfis.le.0.) then gpfis = gpi gnfis = gnu end if c ------------------------------------------------ c 6. Preliminary Quantities for Secondary Emission c ------------------------------------------------ if (multi.gt.0) then rzcn(2) = rzcn(1) rncn(2) = rncn(1)-1. rzcn(3) = rzcn(1)-1. rncn(3) = rncn(1) gres(2,7,1) = gres(1,1,1) gres(2,7,2) = gres(1,1,2) gres(3,7,1) = gres(1,2,1) gres(3,7,2) = gres(1,2,2) gres(2,1,1) = gres(1,1,1) gres(2,1,2) = gres(1,4,2) gres(3,1,1) = gres(1,3,1) gres(3,1,2) = gres(1,3,2) gres(2,2,1) = gres(1,3,1) gres(2,2,2) = gres(1,3,2) gres(3,2,1) = gres(1,5,1) gres(3,2,2) = gres(1,2,2) end if c -------------------------------- c 7. Pairing Condensation Energies c -------------------------------- read (irea,52) ipair,dppi1,dpnu1,dppi2,dpnu2 read (irea,52) iepair,depi1,denu1,depi2,denu2 c ** ipair = -3 simple pairing with default gaps c ** = -2 simple pairing with input gaps (including zero) c ** = -1 simple pairing with eqb. gaps c ** = 0 simple pairing with input gaps (including zero) c ** = 1 collective pairing with eqb. gaps c ** = 2 collective pairing w input gaps (including 0) c ** = 3 collective pairing w default gaps c ** (if default not used for equilibrium) c ** iepair = 0 input equilibrium gaps (including zero) c ** = 1 default equilibrium gaps c ** dppi, dpnu, depi, and denu are traditional pairing energies c ** (as, e.g., from Gilbert and Cameron) for even configurations c ** suffix 1 is for cn or jpout(jnout)=0,1 c ** suffix 2 is for jpout(jnout)=2 if (iepair.gt.0) then evenz = 2. * int(rzcn(1)*0.5+0.01) evenn = 2. * int(rncn(1)*0.5+0.01) depi1 = 9.69 * evenz**(-0.567) denu1 = 9.83 * evenn**(-0.552) depi2 = 9.69 * (evenz-2.)**(-0.567) denu2 = 9.83 * (evenn-2.)**(-0.552) end if ip = ipair if (ip.lt.0) ip=-ip if (ip.eq.1) then dppi1 = depi1 dpnu1 = denu1 dppi2 = depi2 dpnu2 = denu2 else if (ip.ge.3) then evenz = 2. * int(rzcn(1)*0.5+0.01) evenn = 2. * int(rncn(1)*0.5+0.01) dppi1 = 9.69 * evenz**(-0.567) dpnu1 = 9.83 * evenn**(-0.552) dppi2 = 9.69 * (evenz-2.)**(-0.567) dpnu2 = 9.83 * (evenn-2.)**(-0.552) end if c ** Proton Pairing Energies for primary emission. m = rzz+xpin + 0.1 me = (rzz+xpin)*0.5 + 0.1 enod = m - me-me enod1 = 1 - m + me+me c ** enod=0 for even z; enod=1 for odd z in the composite nucleus. cpre(1,7,1) = enod1 * dppi1 cpre(1,1,1) = cpre(1,7,1) cpre(1,2,1) = enod * dppi1 cpre(1,3,1) = cpre(1,2,1) cpre(1,4,1) = cpre(1,2,1) cpre(1,5,1) = enod1 * dppi2 cpre(1,6,1) = cpre(1,5,1) ceq(1,7,1) = enod1 * depi1 ceq(1,1,1) = ceq(1,7,1) ceq(1,2,1) = enod * depi1 ceq(1,3,1) = ceq(1,2,1) ceq(1,4,1) = ceq(1,2,1) ceq(1,5,1) = enod1 * depi2 ceq(1,6,1) = ceq(1,5,1) c ** Neutron Pairing Energies for primary emission. m = rnn+xnin + 0.1 me = (rnn+xnin)*0.5 + 0.1 enod = m - me-me enod1 = 1 - m + me+me cpre(1,7,2) = enod1 * dpnu1 cpre(1,1,2) = enod * dpnu1 cpre(1,2,2) = cpre(1,7,2) cpre(1,3,2) = cpre(1,1,2) cpre(1,4,2) = enod1 * dpnu2 cpre(1,5,2) = cpre(1,1,2) cpre(1,6,2) = cpre(1,4,2) ceq(1,7,2) = enod1 * denu1 ceq(1,1,2) = enod * denu1 ceq(1,2,2) = ceq(1,7,2) ceq(1,3,2) = ceq(1,1,2) ceq(1,4,2) = enod1 * denu2 ceq(1,5,2) = ceq(1,1,2) ceq(1,6,2) = ceq(1,4,2) c ** Pairing Energies for Secondary Emission if (multi.gt.0) then cpre(2,7,1) = cpre(1,1,1) cpre(2,7,2) = cpre(1,1,2) cpre(3,7,1) = cpre(1,2,1) cpre(3,7,2) = cpre(1,2,2) cpre(2,1,1) = cpre(1,1,1) cpre(2,1,2) = cpre(1,4,2) cpre(3,1,1) = cpre(1,3,1) cpre(3,1,2) = cpre(1,3,2) cpre(2,2,1) = cpre(1,3,1) cpre(2,2,2) = cpre(1,3,2) cpre(3,2,1) = cpre(1,5,1) cpre(3,2,2) = cpre(1,2,2) ceq(2,7,1) = ceq(1,1,1) ceq(2,7,2) = ceq(1,1,2) ceq(3,7,1) = ceq(1,2,1) ceq(3,7,2) = ceq(1,2,2) ceq(2,1,1) = ceq(1,1,1) ceq(2,1,2) = ceq(1,4,2) ceq(3,1,1) = ceq(1,3,1) ceq(3,1,2) = ceq(1,3,2) ceq(2,2,1) = ceq(1,3,1) ceq(2,2,2) = ceq(1,3,2) ceq(3,2,1) = ceq(1,5,1) ceq(3,2,2) = ceq(1,2,2) end if icol = 1 x = cpre(1,kin,1) + cpre(1,kin,2) if (x.le.0.) icol = 0 if (ipair.lt.1) icol = 0 c icol greater than 0 invokes collective pairing in c.n. c (using target fermi level to allow for passives) c ---------------------------------------------- c 8. Isospin Parameters for Primary Calculations c (Main work done in subroutine SETISO) c ---------------------------------------------- nt(1) = 0 nt(2) = 0 nt(3) = 0 c ** nt(kc) is T-T(g.s.) in the composite nucleus designated by kc itest = jrn-jrz test = itest tzeff = abs(test) * 0.5 80 if (iso.gt.0) then call setiso(1,rzcn(1),rncn(1)) if (iso.gt.1) then ttest = e / esym(1,7,2) if (ttest.ge.4.) then iso = 0 do 84 kp = 1,6 ct(1,kp,1) = 1. 84 ct(1,kp,2) = 0. go to 90 else iso = 1 fteq = 6. * (ttest-3.5) fteq = 0.4 / (1.+exp(fteq)) end if end if if (e.le.esym(1,7,1)) then write(iwri,88) 88 format(' excit. energy too low') go to 470 end if ct(1,7,1) = ct(1,kin,1) ct(1,7,2) = 0. twot7 = twot(7) twotkin = twot(kin) c ** ct(kc,kp,it) is the square of the isospin coupling c ** Clebsch-Gordan coefficient for the lowest T considered c ** (it=1) and the next higher value in the emission of a c ** particle of type kp from nucleus kc. c ** twot(kc) is 2Tz for the indicated nucleus. end if c --------------------------------------------------------------- c 9. Coulomb Barriers and Channel Averaged Inverse Cross Sections c (used in knockout and inelastic subroutines) c --------------------------------------------------------------- 90 do 94 kp=1,6 sigbar(kp) = 0. athrd = kp - jpout(kp) athrd = (acom(1)-athrd)**0.3333 xpout = jpout(kp) coul(kp) = 0.75 * (rzz+xpin-xpout) * xpout / athrd emu = e - ben(1,kp) c ** If T conservation is added to KNOCK & INEL, emu needs an Esym c ** correction, and more than one emu per channel may be needed. c ** Currently, only T mixed calc'ns are done in KNOCK & INEL. c ** These are scaled by the entrance channel Ct value. esum = 0. do 92 ne=2,neps1(1,kp) if (eps(ne).lt.coul(kp)) go to 92 if (eps(ne).gt.emu) go to 94 eint = (eps(ne+1)-eps(ne-1)) *0.5 esum = esum+eint sigbar(kp) = sigbar(kp) + sigin(kp,ne)*eint 92 continue 94 if (esum.gt.0.) sigbar(kp) = sigbar(kp)/esum coul(7) = 0. if (kin.le.0) go to 100 if (kin.gt.6) go to 100 coul(7) = coul(kin) sigbar(7) = sigbar(kin) c ----------------------------------------------------- c 10. Minimum and Maximum Values of the Particle Number c ----------------------------------------------------- 100 gpi = gres(kc,7,1) gnu = gres(kc,7,2) size = (gpi+gnu) * (e-esym(kc,7,1)) if (size.gt.0.) nup = sqrt(size)*0.75 + 0.5 c ** Above condition always satisfied except for some 2-ary emiss. if (kc.eq.1) then ndwn = max0(jin,1) c ** Currently np0 is also max(jin,1), but could be changed. nup = min0(mxdimp-1,nup) else ndwn = max0(jin-1,1) nup = min0(mxdim2-1,nup) end if ndwn1 = ndwn+1 ndwn2 = ndwn+2 nup1 = nup+1 nup2 = nup+2 c --------------------------------------------------- c 11. Preliminary Calculations for Collective Pairing c --------------------------------------------------- if (icol.le.0) go to 120 do 112 ia=1,2 x = cpre(1,kin,ia) / gres(1,kin,ia) del0 = 2.*sqrt(x) xncrit(7,ia) = 0.792 * del0 * gres(kc,7,ia) do 112 kp=1,kpmax 112 xncrit(kp,ia) = 0.792 * del0 * gres(kc,kp,ia) c ** xncrit gives the critical exciton numbers for protons and c ** neutrons in the nucleus of the kp channel. c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * PAULI CORRECTION FUNCTIONS in the S2-ESM with Simple Pairing * c * (corrections for collective pairing done in subroutines * c * COMPOUND and OMEGA) * c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c 120 iz = rzcn(kc) - rzmag in = rncn(kc) - rnmag do 126 kp=1,7 if (kc.gt.1) then if (kp.gt.2) then if (kp.lt.7) go to 126 c ** This bypasses calculations for d,t,h,alpha in c ** secondary emission where they are not considered. end if end if ds(kp,1) = 1./gres(kc,kp,1) ds(kp,2) = 1./gres(kc,kp,2) idif(kp,1) = iz - jpout(kp) idif(kp,2) = in - jnout(kp) nphd(kc,kp,1) = jpin - jpout(kp) - ipout(kc) nphd(kc,kp,2) = jnin - jnout(kp) - inout(kc) if (ishell.gt.0) then do 124 ia=1,2 if (idif(kp,ia).gt.0) then mdnear(kp,ia) = mhi(ia) mdfar(kp,ia) = mlo(ia) else mdnear(kp,ia) = mlo(ia) mdfar(kp,ia) = mhi(ia) end if 124 continue end if c ** Above loop stores the input shell degeneracies in the c ** appropriate arrays for later use in pauli. c ** If default degeneracies are used, c ** they are stored in mdnear & mdfar inside pauli. call pauli(kc,kp) 126 continue c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * FINITE WELL DEPTH and Surface Effect corrections * c * (surface effect parameters based on ESM assumption) * c * (with isospin conserved, e-esym used in these calcs) * c * This section also bypasses preeq calcns for * c * secondary emission if EBfis. ef = e - bfis(kc) nfish = ef + 1. afish = 1.645 * (gpfis+gnfis) asp = 1.645 * (gnu+gpi) umatch = 2.5 + 150./acom(kc) temp = sqrt(afish*umatch) - 1.25/umatch temp = 1. / temp const = 2. * sqrt(afish*umatch) x = (afish*umatch)**0.25 * umatch const = exp(const) / x x = exp((umatch+ceq(kc,7,1)+ceq(kc,7,2))/temp) const = const / x ref(1) = 2. * sqrt(asp*e) imax = 0 if (iso.gt.0) then if (ee.gt.0.) then imax = 1 ref(2) = 2. * sqrt(asp*ee) else fwt = 0. end if end if do 284 k=1,nfish xk1 = k-1 eff = (ef-xk1) efg = e do 280 i=0,imax c ** i=0 for the T-mixed fission rates; i=1 is for T-conserved if (i.eq.1) then eff = eff - esym(kc,7,1) efg = e - esym(kc,7,1) end if if (eff.le.0.) go to 284 eff = eff - ceq(kc,7,1) - ceq(kc,7,2) if (eff.ge.umatch) then quo = (efg/eff)**1.25 * (asp/afish)**0.25 ex = 2. * sqrt(afish*eff) else eff = eff + ceq(kc,7,1) + ceq(kc,7,2) quo = const * efg**1.25 * asp**0.25 ex = eff / temp end if ex = ex - ref(i+1) ex = quo * exp(ex) if (ex.lt.flow) ex=0. if (ex.gt.spill) ex=0. if (i.eq.0) then fw = ex else fwt = ex end if 280 continue if(iso.gt.0) fishwt = fishwt + fwt 284 fishw = fishw + fw fishw = fishw * 0.242 fishwt = fishwt * 0.242 if (kc.gt.1) then if (ee.le.0.) then c ** Only T-mixed equilibrium emission is possible. c ** y is generated just above stmt 142 where secondary c ** preequilibrium calculations are bypassed when E1 is optional end if clos(np1,npilo1) = p1(npilo1) + p1(npilo2)*gampn(np1,npilo2) clos(np1,npihi1) = p1(npihi1) + p1(npihi)*gamnp(np1,npihi) if (npihi.lt.npilo2) go to 308 do 306 nppi1=npilo2,npihi nppi = nppi1-1 nppi2 = nppi1+1 x = p1(nppi2)*gampn(np1,nppi2) + p1(nppi)*gamnp(np1,nppi) 306 clos(np1,nppi1) = x + p1(nppi1) 308 continue 310 nup1 = np nup = nup1 - 1 nup2 = nup1 + 1 c ** Previous stmts reduce the maximum particle number when strength c ** is exhausted by preeq emission before the old max is reached. c ** Ptot is the strength left & available for equilibrium emission. c ** For secondary emission, strength from primary preeq emission c ** from states with p>mxdim2 must be added into ptot. c ** This strength is in the array spre(kp,ne) c ** and is included next, along with any unused c ** strength in strong(np1,nppi1,neprev) when secondary preeq c ** emission is cut off before mxdim2 is reached. if (kc.gt.1) then ptot = ptot + spre(kc-1,neprev) if (nup2.lt.mxdim2) then do 312 np1=nup2+1,mxdim2 npihi1 = np1 - max0(jnin-inout(kc),0) do 312 nppi1=npilo1,npihi1 312 ptot = ptot + strong(np1,nppi1,neprev) c ** Note: STRONG from np1=nup2 is added in the last pass c ** thru the loop for the closed form rxn calc'ns. end if end if c ** Next calculate the total equilibrium emission rates (all c ** particle types and energies for the isospin mixed (fweis) c ** and isospin conserved (fweist) parts of the cross section. fweis = 0. fweist = 0. do 314 ne=2,neps+1 eint = (eps(ne+1)-eps(ne-1)) * 0.5 if (iso.gt.0) then fw = weist(1,ne) + weist(2,ne) if (kc.eq.1) then fw = fw + weist(3,ne) + weist(4,ne) fw = fw + weist(5,ne) + weist(6,ne) end if fweist = fweist + fw*eint end if fw = weiss(1,ne) + weiss(2,ne) if (kc.eq.1) then fw = fw + weiss(3,ne) + weiss(4,ne) fw = fw + weiss(5,ne) + weiss(6,ne) end if 314 fweis = fweis + fw*eint fweis = fweis + fishw c ** Now convert fweis and fweist into normalization factors c ** that will multiply the specific equilibrium emission rates c ** in order to get the corresponding cross sections. Thus they c ** become the available strength / total eqb. emission rate. c ** The entrance channel cross section and isospin coupling c ** coefficent are applied separately. if (kc.gt.1) eint = (eps(neprev+1)-eps(neprev-1)) * 0.5 if (fweis.gt.0.) then x = (1.-fteq) * ptot if (kc.gt.1) then c ** Next stmts add primary direct reaction strength that was c ** not directed to secondary preeq. emission (as snutra was). if (nt(kc).eq.ntmin(kc-1)) then x = x + seq(kc-1,neprev)*eint sdir = snock(kc-1,neprev) if (kc-1.eq.kin) sdir = sdir + scoll(neprev) x = x + sdir*eint/(sigcni*ct(1,7,1)) end if end if fweis = x / fweis end if if (iso.gt.0) then fweist = fweist + fishwt if (fweist.gt.0.) then x = fteq * ptot if (kc.gt.1) then x = x + steq(kc-1,neprev)*xisoeq(kc-1,neprev)*eint end if fweist = x / fweist end if end if c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * DIRECT REACTION CALCULATIONS and Strength Depletion * c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c if (kc.gt.1) go to 400 avail = 1. used = 0. do 334 kp=1,6 if (jin.gt.4) go to 350 jout = jpout(kp) + jnout(kp) jweak = 1 c ** jweak=1 specifies a weakly bound ejectile (d, t, or h) if (jout.gt.4) go to 350 if (jout.eq.4) jweak=0 if (jout.eq.1) jweak=0 if (jout.eq.jin) go to 332 call nutra(kp,multi,mxdim2) if (jweak.eq.0) call knock(kp) c ** NOTE: As currently programmed knock is called for nucleon or c ** alpha emission, even for d, t or He-3 projectiles which are c ** loosely bound. Projectile cluster survival is assumed and c ** this needs to be checked. Use with caution. go to 334 332 if (jpout(kp).eq.jpin) call inel(kp) c ** NOTE: inel is called even for d, t, and He-3 projectiles which c ** are loosely bound, yet inel assumes projectile cluster c ** survival. Use with caution for these reactions until data c ** can be studied. call nutra(kp,multi,mxdim2) 334 continue do 340 ne=2,neps+1 x = 0. do 338 kp=1,6 338 x = x + snutra(kp,ne) + snock(kp,ne) x = x + scoll(ne) 340 used = used + x*(eps(ne+1)-eps(ne-1))*0.5 sigpre = sigcni*ct(1,7,1) - used c ** sigpre is the cross section entering the equilibration process. c ** avail is the ratio of this to the total rxn cross section. avail = sigpre / (sigcni*ct(1,7,1)) c c * * * * * * * * * * * * * * * * c * PRINT CLOSED FORM STRENGTHS * c * * * * * * * * * * * * * * * * c 350 write (iwri,351) newpg,title 351 format (a1,' ',a64) write (iwri,352) 352 format (' Reaction Strength Passing Thru Configurations') jzmag = rzmag jnmag = rnmag write (iwri,354) jrz,jrn, jpin,jnin, jzmag,jnmag 354 format (' Target Z=',i3,' N=',i3,t27,'Projectile z=' 1,1i2,' n=',1i2,t52,'Shells Z=',i3,' N=',i3) ptot=1.-ptot write (iwri,356) e, sigcni, db(1),db(2) 356 format(' Excit Energy=',f7.2,t27,'Rxn Cross Sect.=',f6.1 1,t52,'Shell gaps=',2f6.2) 358 format (1h ) c write (iwri,360) gpi,gnu, fpp*1.e-6,fnn*1.e-6, fnp*1.e-6 c360 format (' gpi=',f6.3,' gnu=',f6.3,t27,'Kpp=',f5.2,' Knn=',f5.2 c 1,t52,'Knp=',f5.2) write (iwri,360) gpi, gnu, fpp, fnn, fnp 360 format (' gpi=',f6.3,' gnu=',f6.3,t27,'Kpp=',f5.0,' Knn=',f5.0 1,t52,'Knp=',f5.0) if (dppi1+dpnu1+depi1+denu1.gt.0.) then write (iwri,361) ipair,iepair 361 format (' ipair=',i2,t27,'iepair=',i2) else write (iwri,362) 362 format (' no pairing') end if write( iwri,363) vfull,ehole(1),ehole(2) 363 format (' V(central)=',f6.2,t27,'Veff(1)=',f6.2 1,t52,'Veff(2)=',f6.2) if (iso.gt.0) write(iwri,364) nt(1), fteq, esym(1,7,1), ct(1,7,1) 364 format (' T=Tz+',i1,' fteq=',f4.2,t27,'Esym(cn)=',f5.2 1,t52,'Ct2=',f5.3) write (iwri,366) used, sigpre*ptot, ptot 366 format (' Direct sigma=',f7.2,t27,'Preeq sigma=',f7.2 1,t52,'Frac Preeq=',f6.3) c ** Frac Preeq is the fraction of the cross section surviving the c ** direct rxns that is emitted in the exciton model. c ** The overall fraction of preequilibrium emission c ** (direct+exciton) is avail*FracPreeq + 1-avail. write (iwri,358) write (iwri,368)(i,i=npilo1-1,npilo1+5) 368 format (8h p,ppi=,1i2,6i9,9h sum) write (iwri,358) do 374 np1=ndwn1,nup1 closum = 0. npihi1 = np1-jnin np = np1-1 do 370 nppi1=npilo1,npihi1 370 closum=closum+clos(np1,nppi1) write (iwri,372)np,(clos(np1,i),i=npilo1,npilo1+6),closum 372 format (1i3,8f9.5) c ** Renormalize all primary emission strengths for depletion c ** from Direct Reactions (This is carried over for 2ndary emiss.) do 374 nppi1=npilo1,npihi1 374 clos(np1,nppi1) = clos(np1,nppi1)*avail fweis = fweis * avail fweist = fweist * avail c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * Calculate PARTICLE ENERGY SPECTRA * c * (Still will need to be normalized by sigin and Ct^2) * c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c 400 do 402 np1=ndwn1,nup1 npihi1 = np1 - max0(jnin-inout(kc),0) do 402 nppi1=npilo1,npihi1 402 clos(np1,nppi1) = clos(np1,nppi1) * tau(np1,nppi1) 404 if (kc.eq.1) then csres0 = csres0 + sigcni * ct(1,7,1) end if do 408 kp=1,kpmax if (ihf.gt.0) then c ** To generate Hauser-Feshbach (or production cross section) arrays, c ** spec must always be called. Exciton model c ** (and eqb) contributions to these array are generated there. c ** Direct reaction contributions are included here. c ** Arrays for secondary emission residuals are finalized later. call spec(kc,kp,multi,ihf,neprev,sigcni) if (kc.eq.1) then do 406 ne=2,neps1(1,kp) eint = (eps(ne+1)-eps(ne-1)) *0.5 csdir = snock(kp,ne) + snutra(kp,ne) if (kp.eq.kin) csdir = csdir + scoll(ne) csres0 = csres0 - csdir*eint csres1(kp,ne) = csres1(kp,ne) + csdir eres = e - ben(1,kp) - eps(ne) 406 if (eres.ge.0) eres1(kp,ne) = eres end if else if (multi.gt.0) then if (kc.eq.1) then if (kp.lt.3) then call spec(1,kp,multi,0,neprev,sigcni) else if (kprint(kp).gt.0) call spec(1,kp,multi,0,neprev,sigcni) end if else if (kprint(kp).gt.0) call spec(kc,kp,multi,0,neprev,sigcni) end if else if (kprint(kp).gt.0) call spec(kc,kp,multi,0,neprev,sigcni) end if 408 continue if (ifiss.gt.0) fiss(kc) = fiss(kc) + fishw*fweis + fishwt*fweist c c * * * * * * * * * * c * RECYCLE Options * c * * * * * * * * * * c ----------------------------------------- c 1. Recycle for Multiple Particle Emission c ----------------------------------------- 420 if (multi.gt.0) then if (kc.eq.1) then neprev = 2 irecyc = 0 nemin = 2 end if call repent(kc,neps) c ** Subroutine repent does initialization of variables for c ** a new pass thru the main program with a new CN or E if (irecyc.eq.0) then c ** This is first U and T in current primary resid. nucl. go to 100 else if (irecyc.eq.1) then c ** This is a subsequent U and/or T in this nucleus. go to 140 end if end if c ------------------------- c 2. Print Particle Spectra c ------------------------- do 452 kp=1,6 if (kprint(kp).gt.0) then call specout(kp,newpg,multi) if (kang.gt.0) call angel(kp,length,newpg) end if 452 continue c ----------------------------------------- c 3. Recycle to New Composite Isospin Value c ----------------------------------------- c ** If isospin is mixed, then this problem is done. if (iso.le.0) go to 470 c ** If isospin is conserved and if T is already T>, then done. if (nt(1).gt.0) go to 470 c ** If it is a T< composite state then a T> state is possible c ** only for T(targ,g.s) > T(cn,g.s.). For N>Z this means p and c ** He-3 projectiles. For Z>N it means n and t projectiles. if (twotkin.le.twot7) go to 470 nt(1)=1 do 466 ne=1,mxdime preeq2(1,ne) = 0. preeq2(2,ne) = 0. spre(1,ne) = 0. spre(2,ne) = 0. seq2(1,ne) = 0. seq2(2,ne) = 0. steq2(1,ne) = 0. steq2(2,ne) = 0. elast(ne) = 0. elastf(ne) = 0. elcoll(ne) = 0. scoll(ne) = 0. do 462 npp=1,mxdimp+1 do 462 np=1,mxdimp do 462 kp=1,3 462 w(kp,np,npp,ne) = 0. do 464 kp=1,6 preeq1(kp,ne) = 0. seq(kp,ne) = 0. steq(kp,ne) = 0. weiss(kp,ne) = 0. weist(kp,ne) = 0. snock(kp,ne)=0. 464 snutra(kp,ne)=0. do 466 npp=1,mxdim2+1 do 466 np=1,mxdim2 strong(np,npp,ne) = 0. 466 xiso(np,npp,ne) = 0. fiss(1) = 0. fiss(2) = 0. fiss(3) = 0. go to 80 c ------------------------------------------------------------------ c 4. Finalize and print residual nucleus arrays for Hauser-Feshbach c (ihf=1) or production cross section (ihf=2) calculations c ------------------------------------------------------------------ 470 if (ihf.gt.0) call product(ihf,ihfcomp,multi,neps) c ---------------------------------------------- c 5. Recyle to new problem (0 or 1) or exit (-1) c ---------------------------------------------- read (irea,42) newp if (newp.ge.0.) go to 20 1000 continue c c * * * * * * * * c * Hallelujah * c * * * * * * * * c end