MODULE sierk_stuff ! original subroutine of Arnie Sierk !srk This subroutine returns the barrier height bfis, the ground-state !srk energy E_rot, in MeV, and the angular momentum at which the fission !srk barrier disappears, Jmax, in units of h-bar, !srk when called with integer arguments z, the atomic number, !srk a, the atomic mass number, and J, the angular momentum in units !srk of h-bar, (Planck's constant divided by 2*pi). ! !srk The fission barrier for J = 0 is calculated from a 7th order !srk fit in two variables to 638 calculated fission barriers for z value !srk from 20 to 110. These 638 barriers are fit with an rms deviation !srk 0.10 MeV by this 49-parameter function. !srk If barfit is called with (z,a) values outside the range of the !srk the barrier height is set to 999 MeV, and a message is !srk printed on the default output file. ! !srk For J values not equal to zero, the values of !srk J at which the barrier is 80% and 20% of the J=0 value are !srk respectively fit to 20-parameter functions of Z and A, over a more !srk restricted range of A values, than is the case for J = 0. !srk The value of J where the barrier disappears, Jmax, for 61 nuclei !srk is fit to a 35-parameter function of Z and A, with the same range !srk Z and A values as J-80 and J-20. ! Once again, if an (iz,ia) pair is outside of the range of ! validity of the fit, the barrier value is set to 0.0 and a message ! is printed. These three values (Bfis(J=0),J-80, and J-20) and the ! constraints of Bfis = 0 and d(Bfis)/dJ = 0 at J = Jmax and J = 0 ! lead to a fifth-order fit to Bfis(J) for J> L-20. The first three ! constraints lead to a third-order fit for the region J < J-20. ! ! ! The ground-state energies are calculated from a 175-parameter ! fit in Z, A, and L to 329 ground-state energies for 36 different ! Z and A values. !srk (The range of Z and A is the same as for J-80, J-20, and J-max) ! !srk The calculated barriers from which the fits were !srk made were calculated in 1983-1984 by A. J. Sierk of Los Alamos !srk National Laboratory Group T-9, using Yukawa-plus-exponential double !srk folded nuclear energy, exact Coulomb diffuseness corrections, !srk and diffuse-matter moments of inertia. The parameters of the model !srk are those derived by Moller and Nix in 1979: !srk r-0 = 1.16 fm, as = 21.13 MeV, kappa-s = 2.3 a = 0.68 fm. !srk The diffuseness of the matter and charge distributions used !srk corresponds to a surface diffuseness parameter (defined by Myers) !srk of 0.99 fm. The calculated barriers for J = 0 are !srk accurate to a little less than 0.1 MeV; the output from this !srk subroutine is a little less accurate. Worst errors may be as large !srk as 0.5 MeV; characteristic uncertainty is in the range of 0.1-0.2 !srk MeV. !srk Jmax may be in error by up to 1.0 for light nuclei, !srk but is probably good to about 0.2-0.3 for heavy nuclei which do !srk not have prolate ground states (Z >75 or so). ! The values of egs are generally approximated to within ! about 0.1-0.2 MeV; the largest deviation is about 0.5 MeV, ! near L-I for light nuclei.c ! The rms deviation of Lmax from the 61 input values is 0.31 ! h-bar. The approximate value is nearly always within ! 0.5 h-bar of the calculated one. ! !srk written by A. J. Sierk, LANL T-9 !srk Version 1.0 February, 1984 ! Version 1.1 January, 1985 Improved coefficients in egs and Lmax ! Version 1.2 September, 1985 Improved Lmax, egs coefficients ! Version 1.21 June, 1986 minor changes made USE gem_kind IMPLICIT NONE PRIVATE PUBLIC :: jmax_sierk, bfis_sierk, yrast_sierk, m_inertia_sierk, bs_sierk PUBLIC :: fisrot PRIVATE :: lpoly REAL (KIND=r4), PRIVATE :: amin,amax REAL (KIND=r8), PRIVATE :: aa,zz,J20,J80,bfis,JJ,JMAX REAL (KIND=r8), PRIVATE :: zzz,x,y,q,qa,qb,aj,ak,a1,a2,E_rot REAL (KIND=r8), DIMENSION(7), PRIVATE :: pa, pz REAL (KIND=r8), DIMENSION(9), PRIVATE :: pl REAL (KIND=r8), DIMENSION(5), PRIVATE :: paa REAL (KIND=r8), DIMENSION(8), PRIVATE :: pzz REAL (KIND=r8), PRIVATE :: aizro,ai70,aimax,ai95,aimin,bizro,bi70,bimax,bi95 REAL (KIND=r8), PRIVATE :: aimax2,ai952,ff1,ff2,fg1,fg2,aimaxh,aimidh REAL (KIND=r8), PRIVATE :: saizro,sai70,sai95,sai952,simax2,sbimax,sbi70,sbi95 REAL (KIND=r8), PRIVATE :: sbizro,q1,q2,q3,q4,q5,q6,gam,saimax,bb,gam2,aa2 REAL (KIND=r8), PRIVATE :: bb2,aa3,saimid REAL (KIND=r8), PRIVATE :: bb3,gam3,alpha,beta,silt,sjlt,silt2,saimin,saimx REAL (KIND=r8), PRIVATE :: sigt,sjgt,sigt2,f1,f3,f1m,f2,f4,f2m,aa4,bb4 REAL (KIND=r4), PRIVATE :: Z,A REAL (KIND=r8), DIMENSION(5,4), PARAMETER, PRIVATE :: emncof = RESHAPE((/& -9.01100e2_r8,-1.40818e3_r8, 2.77000e3_r8,-7.06695e2_r8, 8.89867e2_r8, & 1.35355e4_r8,-2.03847e4_r8, 1.09384e4_r8,-4.86297e3_r8,-6.18603e2_r8, & -3.26367e3_r8, 1.62447e3_r8, 1.36856e3_r8, 1.31731e3_r8, 1.53372e2_r8, & 7.48863e3_r8,-1.21581e4_r8, 5.50281e3_r8,-1.33630e3_r8, 5.05367e-2_r8/),& (/5,4/)) REAL (KIND=r8), DIMENSION(5,4), PARAMETER, PRIVATE :: elmcof = RESHAPE((/& 1.84542e3_r8,-5.64002e3_r8, 5.66730e3_r8,-3.15150e3_r8, 9.54160e2_r8, & -2.24577e3_r8, 8.56133e3_r8,-9.67348e3_r8, 5.81744e3_r8,-1.86997e3_r8, & 2.79772e3_r8,-8.73073e3_r8, 9.19706e3_r8,-4.91900e3_r8, 1.37283e3_r8, & -3.01866e1_r8, 1.41161e3_r8,-2.85919e3_r8, 2.13016e3_r8,-6.49072e2_r8/),& (/5,4/)) REAL (KIND=r8), DIMENSION(7,5), PARAMETER, PRIVATE :: emxcof = RESHAPE((/& -4.10652732e6_r8, 1.00064947e7_r8,-1.09533751e7_r8, 7.84797252e6_r8, & -3.78574926e6_r8, 1.12237945e6_r8,-1.77561170e5_r8, 1.08763330e7_r8, & -2.63758245e7_r8, 2.85472400e7_r8,-2.01107467e7_r8, 9.48373641e6_r8, & -2.73438528e6_r8, 4.13247256e5_r8,-8.76530903e6_r8, 2.14250513e7_r8, & -2.35799595e7_r8, 1.70161347e7_r8,-8.23738190e6_r8, 2.42447957e6_r8, & -3.65427239e5_r8, 6.30258954e6_r8,-1.52999004e7_r8, 1.65640200e7_r8, & -1.16695776e7_r8, 5.47369153e6_r8,-1.54986342e6_r8, 2.15409246e5_r8, & -1.45539891e6_r8, 3.64961835e6_r8,-4.21267423e6_r8, 3.24312555e6_r8, & -1.67927904e6_r8, 5.23795062e5_r8,-7.66576599e4_r8/),(/7,5/)) REAL (KIND=r8), DIMENSION(7,7), PARAMETER, PRIVATE :: elzcof = RESHAPE((/ & 5.11819909e5_r8,-1.30303186e6_r8, 1.90119870e6_r8,-1.20628242e6_r8, & 5.68208488e5_r8, 5.48346483e4_r8,-2.45883052e4_r8,-1.13269453e6_r8, & 2.97764590e6_r8,-4.54326326e6_r8, 3.00464870e6_r8,-1.44989274e6_r8, & -1.02026610e5_r8, 6.27959815e4_r8, 1.37543304e6_r8,-3.65808988e6_r8, & 5.47798999e6_r8,-3.78109283e6_r8, 1.84131765e6_r8, 1.53669695e4_r8, & -6.96817834e4_r8,-8.56559835e5_r8, 2.48872266e6_r8,-4.07349128e6_r8, & 3.12835899e6_r8,-1.62394090e6_r8, 1.19797378e5_r8, 4.25737058e4_r8, & 3.28723311e5_r8,-1.09892175e6_r8, 2.03997269e6_r8,-1.77185718e6_r8, & 9.96051545e5_r8,-1.53305699e5_r8,-1.12982954e4_r8, 4.15850238e4_r8, & 7.29653408e4_r8,-4.93776346e5_r8, 6.01254680e5_r8,-4.01308292e5_r8, & 9.65968391e4_r8,-3.49596027e3_r8,-1.82751044e5_r8, 3.91386300e5_r8, & -3.03639248e5_r8, 1.15782417e5_r8,-4.24399280e3_r8,-6.11477247e3_r8, & 3.66982647e2_r8/),(/7,7/)) REAL (KIND=r8), DIMENSION(88), PARAMETER, PRIVATE :: egscof_a = (/ & -1.781665232e6_r8,-2.849020290e6_r8, 9.546305856e5_r8, 2.453904278e5_r8, & 3.656148926e5_r8, 4.358113622e6_r8, 6.960182192e6_r8,-2.381941132e6_r8, & -6.262569370e5_r8,-9.026606463e5_r8,-4.804291019e6_r8,-7.666333374e6_r8, & 2.699742775e6_r8, 7.415602390e5_r8, 1.006008724e6_r8, 3.505397297e6_r8, & 5.586825123e6_r8,-2.024820713e6_r8,-5.818008462e5_r8,-7.353683218e5_r8, & -1.740990985e6_r8,-2.759325148e6_r8, 1.036253535e6_r8, 3.035749715e5_r8, & 3.606919356e5_r8, 5.492532874e5_r8, 8.598827288e5_r8,-3.399809581e5_r8, & -9.852362945e4_r8,-1.108872347e5_r8,-9.229576432e4_r8,-1.431344258e5_r8, & 5.896521547e4_r8, 1.772385043e4_r8, 1.845424227e4_r8, 4.679351387e6_r8, & 7.707630513e6_r8,-2.718115276e6_r8,-9.845252314e5_r8,-1.107173456e6_r8, & -1.137635233e7_r8,-1.870617878e7_r8, 6.669154225e6_r8, 2.413451470e6_r8, & 2.691480439e6_r8, 1.237627138e7_r8, 2.030222826e7_r8,-7.334289876e6_r8, & -2.656357635e6_r8,-2.912593917e6_r8,-8.854155353e6_r8,-1.446966194e7_r8, & 5.295832834e6_r8, 1.909275233e6_r8, 2.048899787e6_r8, 4.290642787e6_r8, & 6.951223648e6_r8,-2.601557110e6_r8,-9.129731614e5_r8,-9.627344865e5_r8, & -1.314924218e6_r8,-2.095971932e6_r8, 8.193066795e5_r8, 2.716279969e5_r8, & 2.823297853e5_r8, 2.131536582e5_r8, 3.342907992e5_r8,-1.365390745e5_r8, & -4.417841315e4_r8,-4.427025540e4_r8,-3.600471364e6_r8,-5.805932202e6_r8, & 1.773029253e6_r8, 4.064280430e5_r8, 7.419581557e5_r8, 8.829126250e6_r8, & 1.422377198e7_r8,-4.473342834e6_r8,-1.073350611e6_r8,-1.845960521e6_r8, & -9.781712604e6_r8,-1.575666314e7_r8, 5.161226883e6_r8, 1.341287330e6_r8, & 2.083994843e6_r8, 7.182555931e6_r8, 1.156915972e7_r8,-3.941330542e6_r8/) REAL (KIND=r8), DIMENSION(87), PARAMETER, PRIVATE :: egscof_b = (/ & -1.108259560e6_r8,-1.543982755e6_r8,-3.579820035e6_r8,-5.740079339e6_r8, & 2.041827680e6_r8, 5.981648181e5_r8, 7.629263278e5_r8, 1.122573403e6_r8, & 1.777161418e6_r8,-6.714631146e5_r8,-1.952833263e5_r8,-2.328129775e5_r8, & -1.839672155e5_r8,-2.871137706e5_r8, 1.153532734e5_r8, 3.423868607e4_r8, & 3.738902942e4_r8, 2.421750735e6_r8, 4.107929841e6_r8,-1.302310290e6_r8, & -5.267906237e5_r8,-6.197966854e5_r8,-5.883394376e6_r8,-9.964568970e6_r8, & 3.198405768e6_r8, 1.293156541e6_r8, 1.506909314e6_r8, 6.387411818e6_r8, & 1.079547152e7_r8,-3.517981421e6_r8,-1.424705631e6_r8,-1.629099740e6_r8, & -4.550695232e6_r8,-7.665548805e6_r8, 2.530844204e6_r8, 1.021187317e6_r8, & 1.141553709e6_r8, 2.182540324e6_r8, 3.646532772e6_r8,-1.228378318e6_r8, & -4.813626449e5_r8,-5.299974544e5_r8,-6.518758807e5_r8,-1.070414288e6_r8, & 3.772592079e5_r8, 1.372024952e5_r8, 1.505359294e5_r8, 9.952777968e4_r8, & 1.594230613e5_r8,-6.029082719e4_r8,-2.023689807e4_r8,-2.176008230e4_r8, & -4.902668827e5_r8,-8.089034293e5_r8, 1.282510910e5_r8,-1.704435174e4_r8, & 8.876109934e4_r8, 1.231673941e6_r8, 2.035989814e6_r8,-3.727491110e5_r8, & 4.071377327e3_r8,-2.375344759e5_r8,-1.429330809e6_r8,-2.376692769e6_r8, & 5.216954243e5_r8, 7.268703575e4_r8, 3.008350125e5_r8, 1.114306796e6_r8, & 1.868800148e6_r8,-4.718718351e5_r8,-1.215904582e5_r8,-2.510379590e5_r8, & -5.873353309e5_r8,-9.903614817e5_r8, 2.742543392e5_r8, 9.055579135e4_r8, & 1.364869036e5_r8, 1.895325584e5_r8, 3.184776808e5_r8,-9.500485442e4_r8, & -3.406036086e4_r8,-4.380685984e4_r8,-2.969272274e4_r8,-4.916872669e4_r8, & 1.596305804e4_r8, 5.741228836e3_r8, 6.669912421e3_r8/) REAL (KIND=r8), DIMENSION(5,7,5), PARAMETER, PRIVATE :: egscof = & RESHAPE((/egscof_a,egscof_b/),(/5,7,5/)) REAL (KIND=r8), DIMENSION(6,5), PARAMETER, PRIVATE :: aizroc = RESHAPE((/& 2.34441624e4_r8,-5.88023986e4_r8, 6.37939552e4_r8,-4.79085272e4_r8, & 2.27517867e4_r8,-5.35372280e3_r8,-4.19782127e4_r8, 1.09187735e5_r8, & -1.24597673e5_r8, 9.93997182e4_r8,-4.95141312e4_r8, 1.19847414e4_r8, & 4.18237803e4_r8,-1.05557152e5_r8, 1.16142947e5_r8,-9.00443421e4_r8, & 4.48976290e4_r8,-1.10161792e4_r8,-8.27172333e3_r8, 2.49194412e4_r8, & -3.39090117e4_r8, 3.33727886e4_r8,-1.98040399e4_r8, 5.37766241e3_r8, & 5.79695749e2_r8,-1.61762346e3_r8, 2.14044262e3_r8,-3.55379785e3_r8, & 3.25502799e3_r8,-1.15583400e3_r8/),(/6,5/)) REAL (KIND=r8), DIMENSION(6,5), PARAMETER, PRIVATE :: ai70c = RESHAPE((/& 3.11420101e4_r8,-7.54335155e4_r8, 7.74456473e4_r8,-4.79993065e4_r8, & 2.23439118e4_r8,-4.81961155e3_r8,-7.24025043e4_r8, 1.72276697e5_r8, & -1.72027101e5_r8, 1.03891065e5_r8,-4.83180786e4_r8, 1.08040504e4_r8, & 7.14932917e4_r8,-1.72792523e5_r8, 1.75814382e5_r8,-1.07245918e5_r8, & 4.86163223e4_r8,-1.10623761e4_r8,-2.87206866e4_r8, 6.76667976e4_r8, & -6.50167483e4_r8, 3.67161268e4_r8,-1.74755753e4_r8, 4.67495427e3_r8, & 1.67914908e4_r8,-3.97304542e4_r8, 3.81446552e4_r8,-2.04628156e4_r8, & 7.20091899e3_r8,-1.49978283e3_r8/),(/6,5/)) REAL (KIND=r8), DIMENSION(6,5), PARAMETER, PRIVATE :: ai95c = RESHAPE((/& -6.17201449e5_r8, 1.45561724e6_r8,-1.47514522e6_r8, 9.37798508e5_r8, & -3.74435017e5_r8, 7.81254880e4_r8, 1.24304280e6_r8,-2.94179116e6_r8, & 3.00170753e6_r8,-1.92737183e6_r8, 7.79238772e5_r8,-1.64803784e5_r8, & -1.49648799e6_r8, 3.52658199e6_r8,-3.56784327e6_r8, 2.26413602e6_r8, & -9.02243251e5_r8, 1.88619658e5_r8, 7.27293223e5_r8,-1.72140677e6_r8, & 1.75634889e6_r8,-1.12885888e6_r8, 4.57150814e5_r8,-9.74833991e4_r8, & -3.75965723e5_r8, 8.83032946e5_r8,-8.87134867e5_r8, 5.58350462e5_r8, & -2.20433857e5_r8, 4.62178756e4_r8/),(/6,5/)) REAL (KIND=r8), DIMENSION(6,5), PARAMETER, PRIVATE :: aimaxc = RESHAPE((/& -1.07989556e6_r8, 2.54617598e6_r8,-2.56762409e6_r8, 1.62814115e6_r8, & -6.39575059e5_r8, 1.34017942e5_r8, 2.17095357e6_r8,-5.13081589e6_r8, & 5.19610055e6_r8,-3.31651644e6_r8, 1.31229476e6_r8,-2.77511450e5_r8, & -2.66020302e6_r8, 6.26593165e6_r8,-6.31060776e6_r8, 3.99082969e6_r8, & -1.56447660e6_r8, 3.25613262e5_r8, 1.29464191e6_r8,-3.05746938e6_r8, & 3.09487138e6_r8,-1.97160118e6_r8, 7.79696064e5_r8,-1.63704652e5_r8, & -7.13073644e5_r8, 1.67482279e6_r8,-1.67984330e6_r8, 1.05446783e6_r8, & -4.10928559e5_r8, 8.43774143e4_r8/),(/6,5/)) REAL (KIND=r8), DIMENSION(6,5), PARAMETER, PRIVATE :: ai952c = RESHAPE((/& -7.37473153e5_r8, 1.73682827e6_r8,-1.75850175e6_r8, 1.11320647e6_r8, & -4.41842735e5_r8, 9.02463457e4_r8, 1.49541980e6_r8,-3.53222507e6_r8, & 3.59762757e6_r8,-2.29652257e6_r8, 9.21077757e5_r8,-1.90079527e5_r8, & -1.80243593e6_r8, 4.24319661e6_r8,-4.29072662e6_r8, 2.71416936e6_r8, & -1.07624953e6_r8, 2.20863711e5_r8, 8.86920591e5_r8,-2.09589683e6_r8, & 2.13507675e6_r8,-1.36546686e6_r8, 5.48868536e5_r8,-1.14532906e5_r8, & -4.62131503e5_r8, 1.08555722e6_r8,-1.09187524e6_r8, 6.87308217e5_r8, & -2.70986162e5_r8, 5.61637883e4_r8/),(/6,5/)) REAL (KIND=r8), DIMENSION(6,5), PARAMETER, PRIVATE :: aimax2c = RESHAPE((/& -1.16343311e6_r8, 2.74470544e6_r8,-2.77664273e6_r8, 1.76933559e6_r8, & -7.02900226e5_r8, 1.49345081e5_r8, 2.36929777e6_r8,-5.60655122e6_r8, & 5.70413177e6_r8,-3.66528765e6_r8, 1.47006527e6_r8,-3.15794626e5_r8, & -2.82646077e6_r8, 6.66086824e6_r8,-6.72677653e6_r8, 4.27484625e6_r8, & -1.69427298e6_r8, 3.58429081e5_r8, 1.39112772e6_r8,-3.29007553e6_r8, & 3.34544584e6_r8,-2.14723142e6_r8, 8.61118401e5_r8,-1.84500129e5_r8, & -7.21329917e5_r8, 1.69371794e6_r8,-1.69979786e6_r8, 1.07037781e6_r8, & -4.20662028e5_r8, 8.80728361e4_r8/),(/6,5/)) REAL (KIND=r8), DIMENSION(4,4), PARAMETER, PRIVATE :: aimax3c = RESHAPE((/& -2.88270282e3_r8, 5.30111305e3_r8,-3.07626751e3_r8, 6.56709396e2_r8, & 5.84303930e3_r8,-1.07450449e4_r8, 6.24110631e3_r8,-1.33480875e3_r8, & -4.20629939e3_r8, 7.74058373e3_r8,-4.50256063e3_r8, 9.65788439e2_r8, & 1.23820134e3_r8,-2.28228958e3_r8, 1.33181316e3_r8,-2.87363568e2_r8/), & (/4,4/)) REAL (KIND=r8), DIMENSION(4,4), PARAMETER, PRIVATE :: aimax4c = RESHAPE((/& -3.34060345e3_r8, 6.26384099e3_r8,-3.77635848e3_r8, 8.57180868e2_r8, & 6.76377873e3_r8,-1.26776571e4_r8, 7.64206952e3_r8,-1.73406840e3_r8, & -4.74821371e3_r8, 8.89857519e3_r8,-5.36266252e3_r8, 1.21614216e3_r8, & 1.46369384e3_r8,-2.74251101e3_r8, 1.65205435e3_r8,-3.74262365e2_r8/), & (/4,4/)) REAL (KIND=r8), DIMENSION(6,4), PARAMETER, PRIVATE :: bizroc = RESHAPE((/& 5.88982505e2_r8,-1.35630904e3_r8, 1.32932125e3_r8,-7.78518395e2_r8, & 2.73122883e2_r8,-3.49600841e1_r8,-9.67701343e2_r8, 2.24594418e3_r8, & -2.24303790e3_r8, 1.35440047e3_r8,-4.96538939e2_r8, 6.66791793e1_r8, & 1.17090267e3_r8,-2.71181535e3_r8, 2.67008958e3_r8,-1.58801770e3_r8, & 5.66896359e2_r8,-8.21530057e1_r8,-3.83031864e2_r8, 9.05191483e2_r8, & -9.30560410e2_r8, 5.96618532e2_r8,-2.34403480e2_r8, 3.97909172e1_r8/), & (/6,4/)) REAL (KIND=r8), DIMENSION(6,4), PARAMETER, PRIVATE :: bi70c = RESHAPE((/& 2.32414810e3_r8,-5.42381778e3_r8, 5.40202710e3_r8,-3.26923144e3_r8, & 1.18318943e3_r8,-1.93186467e2_r8,-4.38084778e3_r8, 1.03523570e4_r8, & -1.05573803e4_r8, 6.59901160e3_r8,-2.47601209e3_r8, 4.19497260e2_r8, & 4.35377377e3_r8,-1.01728647e4_r8, 1.01311246e4_r8,-6.14038462e3_r8, & 2.21957562e3_r8,-3.62854365e2_r8,-1.84533539e3_r8, 4.41613298e3_r8, & -4.59403284e3_r8, 2.95951225e3_r8,-1.14630148e3_r8, 2.02702459e2_r8/), & (/6,4/)) REAL (KIND=r8), DIMENSION(6,4), PARAMETER, PRIVATE :: bi95c = RESHAPE((/& 1.55359266e3_r8,-3.58209715e3_r8, 3.50693744e3_r8,-2.03992913e3_r8, & 7.05498010e2_r8,-1.49075519e2_r8,-2.86876240e3_r8, 6.77107086e3_r8, & -6.90300614e3_r8, 4.20246063e3_r8,-1.50290693e3_r8, 3.13662258e2_r8, & 2.60138185e3_r8,-5.95414919e3_r8, 5.70261588e3_r8,-3.17188958e3_r8, & 9.89207911e2_r8,-1.76320647e2_r8,-1.75198402e3_r8, 4.16635208e3_r8, & -4.25212424e3_r8, 2.59953301e3_r8,-9.09813362e2_r8, 1.51070448e2_r8/), & (/6,4/)) REAL (KIND=r8), DIMENSION(6,4), PARAMETER, PRIVATE :: bimaxc = RESHAPE((/& 4.17708254e3_r8,-8.59358778e3_r8, 6.46392215e3_r8,-8.84972189e2_r8, & -1.59735594e3_r8, 1.39662071e3_r8,-1.56318394e4_r8, 3.54574417e4_r8, & -3.35945173e4_r8, 1.65495998e4_r8,-3.32021998e3_r8,-1.46150905e3_r8, & 1.41292811e4_r8,-3.11818487e4_r8, 2.77454429e4_r8,-1.19628827e4_r8, & 1.28008968e3_r8, 1.66111636e3_r8,-1.92878152e4_r8, 4.56505796e4_r8, & -4.66413277e4_r8, 2.89229633e4_r8,-1.07284346e4_r8, 1.50513815e3_r8/), & (/6,4/)) REAL (KIND=r8), DIMENSION(100), PARAMETER, PRIVATE :: b_a = (/& -17605.486, 7149.518, 23883.053, 21651.272, -348.402, & 67617.841, -43739.724, -72255.975, -79130.794, -44677.571, & -38284.438, 1629.405, 49928.158, 41804.076, -4744.661, & 32124.997, -30258.577, -45416.684, -40347.682, -20718.220, & -7009.141, -12274.997, 5703.797, 4375.711, -2519.378, & 44633.464, -19664.988, -59871.761, -55426.914, -834.542, & -168222.448, 111291.480, 179939.092, 198092.569, 111347.893, & 97361.913, -9404.674, -125757.922, -108003.071, 7218.634, & -79483.915, 77325.943, 113433.935, 100793.473, 51134.715, & 17974.287, 28084.202, -15187.024, -12244.758, 4662.638, & -52775.944, 26372.246, 69061.418, 66412.802, 4456.150, & 191743.968, -131672.398, -204896.966, -227867.388, -126609.907, & -115704.145, 22053.324, 146315.778, 131364.763, 1176.787, & 89552.374, -92216.043, -129901.381, -115369.145, -56853.433, & -21565.916, -26595.626, 19490.177, 16758.541, -2102.079, & 43089.492, -24205.277, -54145.651, -54916.885, -6982.880, & -148893.417, 105958.013, 157835.761, 178691.634, 97246.541, & 94988.919, -28100.912, -115999.016, -110379.031, -10165.655, & -68277.650, 75226.302, 100957.904, 89635.423, 41887.718, & 17720.511, 14837.634, -17476.245, -15762.519, -1297.194/) REAL (KIND=r8), DIMENSION(100), PARAMETER, PRIVATE :: b_b = (/& -25673.002, 15584.030, 30291.575, 33008.500, 6161.347, & 83671.721, -60540.161, -86493.251, -101211.404, -53525.729, & -56872.290, 21837.845, 65632.326, 67209.922, 11367.245, & 37417.011, -43832.496, -56023.110, -49995.410, -21521.767, & -10406.752, -4483.654, 11314.537, 10435.311, 2206.541, & 11115.697, -6672.659, -11696.751, -14183.580, -3586.280, & -33914.223, 23453.203, 32630.486, 40847.174, 21223.389, & 24704.794, -10226.568, -25497.313, -29082.833, -7190.865, & -14734.982, 17457.735, 21457.946, 19706.266, 7745.219, & 4283.242, 214.103, -5058.390, -4715.987, -1303.557, & -3196.271, 1579.967, 2718.777, 3997.275, 1361.541, & 9152.527, -5182.345, -7412.465, -10822.583, -5804.143, & -7177.709, 2317.015, 5874.736, 8208.479, 2858.079, & 3919.329, -4068.562, -4982.296, -5065.564, -1948.636, & -1147.488, 265.063, 1372.992, 1333.320, 416.369, & 442.331, -219.051, -337.663, -575.255, -232.634, & -1217.351, 590.169, 854.438, 1468.791, 811.481, & 1031.157, -227.306, -651.183, -1174.558, -528.960, & -536.534, 413.429, 533.977, 663.064, 278.203, & 154.548, -54.623, -169.817, -183.006, -63.999/) REAL (KIND=r8), DIMENSION(5,5,8), PARAMETER, PRIVATE :: b =& RESHAPE((/b_a,b_b/),(/5,5,8/)) CONTAINS SUBROUTINE jmax_sierk (z_in,a_in,Jmax_s) ! ! ! DUMMY ARGUMENTS: ! ! input: Z (REAL (KIND=4) ::) proton number of nucleus ! A (REAL (KIND=4) ::) nucleon number of nucleus ! J (REAL (KIND=4) ::) spin of nucleus (hbar) ! output: ! bfis (REAL (KIND=4) ::) - symmetric fission barrier of nucleus ! Jmax (REAL (KIND=4) ::) - angular momentum at which bfis=0 MeV ! ! ! C H A N G E L O G ! ! Date | Name | Description ! ----------------+-------+----------------------------------------------------- ! 21-11-89 | RJC | I have beautified this routine of Sierks and made all ! the double precision variables single precision. ! ----------------+-------+----------------------------------------------------- ! 30-8-90 | RJC | corrected two errors of no probable consequence ! ! error 1/ there were two variables Z in routine ! ! 2/ there was a el instead of J in if statement ! ----------------+-------+----------------------------------------------------- ! REAL (KIND=r4), INTENT(IN) :: Z_in,A_in REAL (KIND=r4), INTENT(OUT), OPTIONAL :: Jmax_s INTEGER :: i,k z = z_in A = a_in IF (z < 19.0 .OR. z > 111.0) THEN WRITE (UNIT=*, FMT=*) "z=",Z, "out of limits" IF (PRESENT(Jmax_s)) jmax_s = -1.0 RETURN END IF amin = 1.2*z + 0.01*z**2.0 amax = 5.8*z - 0.024*z**2.0 IF (a < amin .OR. a > amax) THEN WRITE (UNIT=*, FMT=*) "A out of limits" IF (PRESENT(Jmax_s)) jmax_s = -1.0 RETURN END IF aa = 2.5e-3_r8 *REAL(a,KIND=r8) zz = 1.0e-2_r8 *REAL(z,KIND=r8) jmax= 0.0 CALL lpoly (zz,7,pz) CALL lpoly (aa,7,pa) DO i=1,5 DO k=1,7 Jmax = Jmax + emxcof(k,i)*pz(k)*pa(i) END DO END DO IF (PRESENT(jmax_s)) jmax_s = jmax RETURN END SUBROUTINE jmax_sierk SUBROUTINE bfis_sierk(j,bfis_s) REAL (KIND=r4), INTENT(IN) :: J REAL (KIND=r4), INTENT(OUT) :: bfis_s INTEGER :: i,k bfis_s = 999.0 IF (z > 102.0 .AND. J > 0.0) THEN WRITE (UNIT=*, FMT=*) "z and J out of limits" RETURN END IF !bfis = 0 !DO i=1,7 ! DO k=1,7 ! bfis = bfis + elzcof(k,i)*pz(k)*pa(i) ! END DO !END DO bfis = DOT_PRODUCT(pz,MATMUL(elzcof,pa)) IF (J < 1) THEN bfis_s = bfis RETURN END IF J80 = 0.0 J20 = 0.0 DO i=1,4 DO k=1,5 J80 = J80 + elmcof(k,i)*pz(k)*pa(i) J20 = J20 + emncof(k,i)*pz(k)*pa(i) END DO END DO IF (J <= J20) THEN q = 0.2/(J20**2*J80**2*(J20-J80)) qa = q*(4.0*J80**3-J20**3) qb = -q*(4.0*J80**2-J20**2) bfis = bfis*(1.0+qa*REAL(J,KIND=r8)**2.0+qb*REAL(J,KIND=r8)**3.0) ELSE x = J20/Jmax y = J80/Jmax aj = (-20.0*x**5+25.0*x**4-4.0)*(y-1.0)**2*y*y ak = (-20.0*y**5+25.0*y**4-1.0)*(x-1.0)**2*x*x q = 0.2/((y-x)*((1.0-x)*(1.0-y)*x*y)**2) qa = q*(aj*y-ak*x) qb = -q*(aj*(2.0*y+1.0)-ak*(2.0*x+1.0)) zzz = REAL(J,KIND=r8)/Jmax a1 = 4.0*zzz**5-5.0*zzz**4 + 1.0 a2 = qa*(2.0*zzz+1.0) bfis = bfis*(a1+(zzz-1.0)*(a2+qb*zzz)*zzz*zzz*(zzz-1.0)) END IF IF (bfis <= 0.0) bfis = 0.0 IF (J > Jmax) bfis = 0.0 bfis_s = bfis RETURN END SUBROUTINE bfis_sierk SUBROUTINE yrast_sierk(J,E_rot_s) REAL (KIND=r4), INTENT(IN) :: J REAL (KIND=r4), INTENT(OUT) :: E_rot_s INTEGER :: k,l,m E_rot_s = 999.0 IF (z > 102.0 .AND. J > 0.0) THEN WRITE (UNIT=*, FMT=*) "z and J out of limits" RETURN END IF IF (A < amin .OR. A > amax) THEN WRITE (UNIT=*, FMT=*) "A is out of range" RETURN END IF E_rot = 0.0 JJ = REAL(J,KIND=r8)/Jmax CALL lpoly (jj,9,pl) !srk Now calculate rotating ground-state energy ! IF (J > Jmax) RETURN DO k=1,5 DO l=1,7 DO m=1,5 E_rot = E_rot & + egscof(m,l,k)*pz(l)*pa(k)*pl(2*m-1) END DO END DO END DO IF (E_rot < 0.0) E_rot = 0.0 E_rot_S = E_rot RETURN END SUBROUTINE yrast_sierk SUBROUTINE M_inertia_sierk(J,M_inertia_1,M_inertia_2,M_inertia_3) REAL (KIND=r4), INTENT(IN) :: J REAL (KIND=r4), INTENT(OUT), OPTIONAL :: m_inertia_1,m_inertia_2,m_inertia_3 INTEGER :: l,k JJ = REAL(J,KIND=r8)/Jmax aizro = 0.0_r8 ai70 = 0.0_r8 aimax = 0.0_r8 ai95 = 0.0_r8 aimin = 0.0_r8 bizro = 0.0_r8 bi70 = 0.0_r8 bimax = 0.0_r8 bi95 = 0.0_r8 aimax2 = 0.0_r8 ai952 = 0.0_r8 DO l=1,6 DO k = 1,5 aizro = aizro + aizroc(l,k) *pz(l)*pa(k) ai70 = ai70 + ai70c(l,k) *pz(l)*pa(k) ai95 = ai95 + ai95c(l,k) *pz(l)*pa(k) aimax = aimax + aimaxc(l,k) *pz(l)*pa(k) ai952 = ai952 + ai952c(l,k) *pz(l)*pa(k) aimax2 = aimax2 + aimax2c(l,k)*pz(l)*pa(k) END DO DO k = 1,4 bizro = bizro + bizroc(l,k)*pz(l)*pa(k) bi70 = bi70 + bi70c(l,k) *pz(l)*pa(k) bi95 = bi95 + bi95c(l,k) *pz(l)*pa(k) bimax = bimax + bimaxc(l,k)*pz(l)*pa(k) END DO END DO ff1 = 1.0 ff2 = 0.0 fg1 = 1.0 fg2 = 0.0 IF (Z > 70) THEN aimaxh = 0.0 aimidh = 0.0 DO l=1,4 DO k=1,4 aimaxh = aimaxh + aimax3c(l,k)*pz(l)*pa(k) aimidh = aimidh + aimax4c(l,k)*pz(l)*pa(k) END DO END DO IF (z > 80) ff1 = 0.0 IF (z >= 80) fg1 = 0.0 IF (bimax > 0.95) fg1 = 0.0 IF (aimaxh > aimax) ff1 = 0.0 ff2 = 1.0 - ff1 fg2 = 1.0 - fg1 aimax = aimax*ff1 + ff2*aimaxh aimax2 = aimax2*ff1 + ff2*aimidh END IF bimax = bimax*fg1 + aimidh*fg2 saizro = MAX(aizro,0.0_r8) sai70 = MAX(ai70,0.0_r8) sai95 = MAX(ai95,0.0_r8) saimax = MAX(aimax,0.0_r8) sai952 = MAX(ai952,0.0_r8) simax2 = MAX(aimax2,0.0_r8) sbimax = MAX(bimax,0.0_r8) sbi70 = MAX(bi70,0.0_r8) sbi95 = MAX(bi95,0.0_r8) sbizro = MAX(bizro,0.0_r8) q1 = -3.148849569 q2 = 4.465058752 q3 = -1.316209183 q4 = 2.26129233 q5 = -4.94743352 q6 = 2.68614119 gam = - 20.0*LOG(ABS(saizro-sai95)/ABS(saizro-saimax)) aa = q1*saizro + q2*sai70 + q3*sai95 bb = q4*saizro + q5*sai70 + q6*sai95 gam2 = - 20.0*LOG(ABS(saizro-sai952)/ABS(saizro-simax2)) aa2 = q1*saizro + q2*sai70 + q3*sai952 bb2 = q4*saizro + q5*sai70 + q6*sai952 aa3 = q1*sbizro + q2*sbi70 + q3*sbi95 bb3 = q4*sbizro + q5*sbi70 + q6*sbi95 gam3 = 60.0 alpha = pi*(jj - 0.7) beta = 5.0*pi*(jj - 0.9) silt = saizro + aa*jj**2.0 + bb*jj**4.0 sjlt = sbizro + aa3*jj**2.0 + bb3*jj**4.0 silt2 = saizro + aa2*jj**2.0 + bb2*jj**4.0 IF (jj <= 0.70) THEN saimin = sjlt saimx = silt saimid = silt2 ELSE sigt = saizro + (saimax-saizro)*EXP(gam*(jj-1.0)) sjgt = sbi95 + (sbimax-sbi95)*EXP(gam3*(jj-1.0)) sigt2 = saizro + (simax2-saizro)*EXP(gam2*(jj-1.0)) IF (jj <= 0.95) THEN f1 = silt*(COS(alpha))**2 + sigt*(SIN(alpha))**2 saimx = f1 f3 = sjlt*(COS(alpha))**2 + sjgt*(SIN(alpha))**2 saimin = f3 f1m = silt2*(COS(alpha))**2 + sigt2*(SIN(alpha))**2 saimid = f1m ELSE f2 = silt*(COS(beta))**2 + sigt*(SIN(beta))**2 saimx = f2 f4 = sjlt*(COS(beta))**2 + sjgt*(SIN(beta))**2 saimin = f4 f2m = silt2*(COS(beta))**2 + sigt2*(SIN(beta))**2 saimid = f2m END IF END IF IF (ff2 > 0.01 .AND. fg2 > 0.01) THEN q1 = 4.001600640 q2 = 0.960784314 q3 = 2.040816327 aa3 = q1*sai70 - q2*saimax - (1.0+q3)*saizro bb3 = -q1*sai70 + (1.0+q2)*saimax + q3*saizro aa4 = q1*sai70 - q2*simax2 - (1.0+q3)*saizro bb4 = -q1*sai70 + (1.0+q2)*simax2 + q3*saizro saimx = saizro + aa3*jj**2.0 + bb3*jj**4.0 saimid = saizro + aa4*jj**2.0 + bb4*jj**4.0 END IF IF (saimid > saimx) saimid = saimx saimin = MAX(saimin,0.0_r8) IF (PRESENT(M_inertia_1)) m_inertia_1 = saimx IF (PRESENT(M_inertia_2)) M_inertia_2 = saimid IF (PRESENT(M_inertia_3)) M_inertia_3 = saimin RETURN END SUBROUTINE M_inertia_sierk SUBROUTINE bs_sierk (j,bs_s) REAL (kind=r4), INTENT (IN) :: j REAL (kind=r4), INTENT (OUT) :: bs_s INTEGER :: iz,ia,ill REAL (KIND=r8) :: xa,bs xa = a/320.0_r8 JJ = REAL(J,KIND=r8)/Jmax CALL lpoly (jj,9,pl) CALL lpoly (xa,5,paa) CALL lpoly (zz,8,pzz) bs = 0.0_r8 DO iz = 1,8 DO ia = 1,5 DO ill = 1,5 bs = bs + b(ill,ia,iz)*pzz(iz)*paa(ia)*pl(2*ill-1) END DO END DO END DO bs_s = REAL(bs) END SUBROUTINE bs_sierk SUBROUTINE lpoly (x,n,pl) ! ! FUNCTIONAL DESCRIPTION: ! !srk this subroutine calculates the ordinary Legendre Polynomials of !srk order 0 to n-1 of argument x and stores them in the vector !srk pl. They are calculated by recursion relation from the first two !srk polynomials. ! !srk written by A. J. Sierk LANL T-9 February,1984 ! !srk NOTE: pl and x must be REAL (KIND=8) :: on 32-bit computers! ! ! DUMMY ARGUMENTS: ! ! input: x (REAL (KIND=4) ::) - argument of Legendre Polynomial ! n (INTEGER (kind=4) ::) - n-1 is highest order of polynomial ! output: pl(1:20) (REAL (KIND=4) ::) - value of Legendre Polynomials ! INTEGER, INTENT(IN) :: n REAL (KIND=r8), INTENT(IN) :: x REAL (KIND=r8), DIMENSION(:), INTENT(OUT) :: pl INTEGER :: i pl(1) = 1.0 pl(2) = x DO i=3,n pl(i) = (REAL(2*i-3,KIND=r8)*x*pl(i-1) & - REAL(i-2,KIND=r8)*pl(i-2))/REAL(i-1,KIND=r8) END DO RETURN END SUBROUTINE lpoly SUBROUTINE fisrot (Z,A,al,bf) ! Cohen-Plasil-Swiatecki rotating liquid drop fission barrier ! Provided by courtesy of Dr. F. Plasil, Oak Ridge National Lab. ! REAL (KIND=r4), INTENT(IN) :: a,z,al REAL (KIND=r4), INTENT(OUT) :: bf REAL (KIND=r4), PARAMETER, DIMENSION(6,11) :: x1b = RESHAPE ( & (/0.28,.243,.221,.208,.195,.18,.211,.186,.17,.1506,.136,.12, & 0.152,.131,.1155,.096,.0795,.0625,.09725,.0795,.065,.0506,.0375,& 0.0253,.05771,.0455,.03414,.0235,.014,.0065,.03325,.0235,.0153,.0081, & 0.001,.0,.01625,.009,.0032,.0,.0,.0,.0071,.0,.0,.0,.0,.0,.0,.0, & 0.0,.0,.0,.0,.0,.0,.0,.0,.0,.0,0.,0.,0.,0.,0.,0./),(/6,11/)) REAL (KIND=r4), PARAMETER, DIMENSION(6,11) :: x2b = RESHAPE ( & (/.18,.1695,.1515,.133,.1155,.0949,.1495,.1363,.1165,.099, & .0815,.0594,.12,.1032,.0864,.0678,.0469,.028,.09,.0725,.0556,.037, & .019,.0057,.0625,.045,.0304,.016,.005,0.,.0406,.0264,.0151,.0052, & 0.,0.,.0253,.0144,.0027,0.,0.,0.,.0141,.006,0.,0.,0.,0.,.0065,.0008, & 0.,0.,0.,0.,.002,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./),(/6,11/)) REAL (KIND=r4), PARAMETER, DIMENSION(10,20) :: x3b = RESHAPE ( & (/.0949,.0755,.0564,.0382,.0223,.0121,.00588,.00242,.00069, & .0001,.0873,.0684,.049,.0306,.0162,.0074,.00267,.00055,0.,0.,.0801, & .061,.0418,.0235,.0108,.00373,.00071,0.,0.,0.,.073,.054,.035,& .0178,.0062,.00125,0.,0.,0.,0.,.0661,.047,.0284,.012,.0025,0.,0.,0., & 0.,0.,.0594,.0404,.022,.0065,0.,0.,0.,0.,0.,0.,.0528,.034,.0159, & .002,0.,0.,0.,0.,0.,0.,.0465,.0277,.01,0.,0.,0.,0.,0.,0.,0.,.0401, & .0217,.0044,0.,0.,0.,0.,0.,0.,0.,.0339,.0158,.00024,0.,0.,0.,0., & 0.,0.,0.,.028,.0106,0.,0.,0.,0.,0.,0.,0.,0.,.0219,.0064,0.,0.,0., & 0.,0.,0.,0.,0.,.0164,.0025,0.,0.,0.,0.,0.,0.,0.,0.,.0122,0.,0.,0., & 0.,0.,0.,0.,0.,0.,.0085,0.,0.,0.,0.,0.,0.,0.,0.,0.,.0057,0.,0., & 0.,0.,0.,0.,0.,0.,0.,.0035,0.,0.,0.,0.,0.,0.,0.,0.,0.,.0016,0.,0., & 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,0.,0.,0./),(/10,20/)) REAL (KIND=r4) :: an,paren,eso,cx,bx,dx,by,cy,dy,b2,b1 INTEGER :: ix,iy an=a-z paren=1.-1.7826*((an-z)/a)**2 eso=17.9439*paren*a**.66667 x=0.019655*z*(z/a)/paren y=1.9254*al**2/(paren*a**2.3333) ix=20.*x+.999 cx=ix bx=20.*x+.999 dx=bx-cx IF (x <= 0.25) THEN by=10.*y+.999 IF (by.GT.9.) by=9. IF (by.LT.1.) by=1. iy=by cy=iy dy=by-cy b2=(x1b(ix+1,iy+1)-x1b(ix,iy+1))*dx+x1b(ix,iy+1) b1=(x1b(ix+1,iy)-x1b(ix,iy))*dx+x1b(ix,iy) bf=(b2-b1)*dy+b1 ELSE IF (x <= .5) THEN by=20.*y+.999 IF (by.GT.11.) by=11. IF (by.LT.1.) by=1. ix=ix-5 iy=by cy=iy dy=by-cy b1=(x2b(ix+1,iy)-x2b(ix,iy))*dx+x2b(ix,iy) b2=(x2b(ix+1,iy+1)-x2b(ix,iy+1))*dx+x2b(ix,iy+1) bf=(b2-b1)*dy+b1 ELSE x=MIN(x,0.95) ix=20.*x+.999 ix=ix-10 by=100.*y+.999 IF (by.GT.19.) by=19. IF (by.LT.1.) by=1. iy=by cy=iy dy=by-cy b1=(x3b(ix+1,iy)-x3b(ix,iy))*dx+x3b(ix,iy) b2=(x3b(ix+1,iy+1)-x3b(ix,iy+1))*dx+x3b(ix,iy+1) bf=(b2-b1)*dy+b1 END IF bf=bf*eso RETURN END SUBROUTINE fisrot END MODULE sierk_stuff