Project CMSSW displayed by LXR CMSSW_14_1_X_2024-07-25-2300/​GeneratorInterface/​ReggeGribovPartonMCInterface/​src/​nuclib_21.f	Source navigation Diff markup Identifier search General search
	Version: CMSSW_14_1_X_2024-07-25-2300 CMSSW_14_1_X_2024-07-24-2300 CMSSW_14_1_X_2024-07-23-2300 CMSSW_14_1_X_2024-07-22-2300 CMSSW_14_1_X_2024-07-21-2300 CMSSW_14_1_X_2024-07-19-2300 Revert   Change

File indexing completed on 2024-04-06 12:14:07

 C==========================================================================
 C.  Library of programs for the generation of nucleus-nucleus interactions
 C.  and the study of nucleus-induced cosmic ray showers
 C.
 C.  September 2001  changes in FPNI, and SIGMA_INI,
 C.                  new SIGMA_PP, SIGMA_PPI, SIGMA_KP  (R. Engel)
 C.
 C.  may  1996       small bug  corrected by Dieter Heck in NUC_CONF
 C.
 C.  march 1996      small modification to the superposition code
 C.
 C.  February 1996   change to FPNI to give an interaction length
 C.                   also  at very low energy
 C.
 C.  Version 1.01  september 1995
 C.       (small corrections P.L.)
 C.       the random number generator is called as S_RNDM(0)
 C.  ------------------------------------------------------
 C.  Version 1.00  April 1992
 C.
 C.  Authors:
 C.
 C.     J. Engel
 C.     T.K Gaisser
 C.     P.Lipari
 C.     T. Stanev
 C.
 C.  This set of routines  when used in  the simulation of cosmic ray
 C.  showers have only three  "contact points" with the "external world"
 C.
 C.    (i) SUBROUTINE NUC_NUC_INI
 C.        (no  calling arguments)
 C.         to be called once during general initialization
 C.    (ii) SUBROUTINE HEAVY (IA, E0)
 C.         where IA (integer) is the mass number of the projectile
 C.         nucleus  and E0 (TeV) is the energy per nucleon
 C.         The output (positions of first interaction for the IA
 C.         nucleons of the projectile) is  contained in  the common block:
 C.           COMMON /C1STNC/ XX0(60),XX(60),YY(60),AX(60),AY(60)
 C.         In detail:
 C.             XX0(j)   (g cm-2) =  position of interaction
 C.             XX(j) (mm)    x-distance from shower axis
 C.             YY(j) (mm)    y-distance from shower axis
 C.             AX(j) (radiants)  Theta_x with respect to original direction
 C.             AY(j) (radiants)  Theta_y with respect to original direction
 C.
 C.    (iii)  FUNCTION FPNI (E,L)
 C.           Interaction length in air.
 C.           E (TeV) is the energy of the particle, L is the particle
 C.           code (NOTE: "Sibyll" codes are used : L =1-18)
 C.           WANRNING : The nucleus-nucleus cross section
 C.           tabulated in the program are "matched" to the p-Air
 C.           cross section calculated  with this FUNCTION, in other words
 C.           they are both calculated with the same input pp cross section
 C==========================================================================
 
       SUBROUTINE NUC_NUC_INI
 C...Initialization for the generation of nucleus-nucleus interactions
 C.  INPUT : E0 (TeV) Energy per nucleon of the beam nucleus
 C........................................................................
       SAVE
 
       CALL NUC_GEOM_INI                       ! nucleus profiles
       CALL SIGMA_INI                          ! initialize pp cross sections
 
       RETURN
       END
 
 
       SUBROUTINE FRAGM1 (IA,NW, NF, IAF)
 C...Nuclear Fragmentation
 C.  total dissolution of nucleus
 C..........................................
       SAVE
 
       DIMENSION IAF(60)
       NF = IA-NW
       DO J=1,NF
          IAF(J) = 1
       ENDDO
       RETURN
       END
 C->
       SUBROUTINE FRAGM2 (IA,NW, NF, IAF)
 C...Nuclear Fragmentation
 C.  Spectator in one single fragment
 C..........................................
       SAVE
 
       DIMENSION IAF(60)
       IF (IA-NW .GT. 0)  THEN
          NF = 1
          IAF(1) = IA-NW
       ELSE
          NF = 0
       ENDIF
       RETURN
       END
 
 C====================================================================
 C...Code of fragmentation  of spectator nucleons
 C.  based on Jon Engel  abrasion-ablation algorithms
 C====================================================================
 
       BLOCK DATA FRAG_DATA
       SAVE
 
 C...Data for the fragmentation of  nucleus  projectiles
       COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
       DATA (NFLAGG(I, 1),I=1,10)  /
      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
       DATA (NFLAGG(I, 2),I=1,10)  /
      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
       DATA (NFLAGG(I, 3),I=1,10)  /
      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
       DATA (NFLAGG(I, 4),I=1,10)  /
      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
       DATA (NFLAGG(I, 5),I=1,10)  /
      +    0,  0,  0,  0,  0,  0,  0,  0,  0,  0 /
       DATA (NFLAGG(I, 6),I=1,10)  /
      +    0,  0,  0,  0,  0,  0,  0,  1,  1,  1 /
       DATA (NFLAGG(I, 7),I=1,10)  /
      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
       DATA (NFLAGG(I, 8),I=1,10)  /
      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
       DATA (NFLAGG(I, 9),I=1,10)  /
      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
       DATA (NFLAGG(I,10),I=1,10)  /
      +    1,  1,  1,  1,  1,  1,  1,  1,  1,  1 /
       DATA (A(I, 1, 1),I=1,10)  /
      +  .438E-01,.172    ,.283    ,.511    ,.715    ,.920    ,1.19    ,
      +  1.37    ,1.65    ,2.14     /
       DATA (A(I, 1, 2),I=1,10)  /
      +  .147E-01,.249E-01,.439E-01,.592E-01,.776E-01,.886E-01,.108    ,
      +  .117    ,.126    ,.128     /
       DATA (A(I, 1, 3),I=1,10)  /
      +  .216E-02,.627E-02,.834E-02,.108E-01,.144E-01,.152E-01,.196E-01,
      +  .200E-01,.210E-01,.224E-01 /
       DATA (A(I, 1, 4),I=1,10)  /
      +  .593E-01,.653E-01,.116    ,.145    ,.184    ,.204    ,.234    ,
      +  .257    ,.271    ,.248     /
       DATA (A(I, 1, 5),I=1,10)  /
      +  .000E+00,.918E-02,.362E-02,.805E-02,.436E-02,.728E-02,.466E-02,
      +  .707E-02,.932E-02,.130E-01 /
       DATA (A(I, 1, 6),I=1,10)  /
      +  .000E+00,.180E-02,.247E-02,.208E-02,.224E-02,.214E-02,.226E-02,
      +  .233E-02,.230E-02,.194E-02 /
       DATA (A(I, 1, 7),I=1,10)  /
      +  .000E+00,.106E-02,.703E-03,.687E-03,.739E-03,.674E-03,.819E-03,
      +  .768E-03,.756E-03,.720E-03 /
       DATA (A(I, 1, 8),I=1,10)  /
      +  .000E+00,.000E+00,.188E-02,.130E-02,.138E-02,.117E-02,.124E-02,
      +  .119E-02,.111E-02,.829E-03 /
       DATA (A(I, 1, 9),I=1,10)  /
      +  .000E+00,.000E+00,.302E-03,.258E-03,.249E-03,.208E-03,.248E-03,
      +  .222E-03,.210E-03,.187E-03 /
       DATA (A(I, 1,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.235E-03,.222E-03,.172E-03,.181E-03,
      +  .166E-03,.152E-03,.124E-03 /
       DATA (A(I, 1,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.238E-03,.179E-03,.145E-03,.156E-03,
      +  .138E-03,.129E-03,.111E-03 /
       DATA (A(I, 1,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.368E-03,.400E-03,.255E-03,.262E-03,
      +  .221E-03,.182E-03,.112E-03 /
       DATA (A(I, 1,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.753E-04,.712E-04,.527E-04,
      +  .537E-04,.538E-04,.487E-04 /
       DATA (A(I, 1,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.103E-03,.589E-04,.578E-04,
      +  .468E-04,.385E-04,.269E-04 /
       DATA (A(I, 1,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.444E-04,.372E-04,
      +  .318E-04,.284E-04,.218E-04 /
       DATA (A(I, 1,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.487E-04,.473E-04,
      +  .338E-04,.243E-04,.122E-04 /
       DATA (A(I, 1,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.121E-04,.117E-04,
      +  .932E-05,.792E-05,.583E-05 /
       DATA (A(I, 1,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.147E-04,
      +  .101E-04,.756E-05,.496E-05 /
       DATA (A(I, 1,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.755E-05,
      +  .612E-05,.505E-05,.341E-05 /
       DATA (A(I, 1,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .630E-05,.444E-05,.282E-05 /
       DATA (A(I, 2, 1),I=1,10)  /
      +  .269    ,.510    ,.738    ,1.12    ,1.46    ,1.83    ,2.22    ,
      +  2.57    ,3.00    ,3.67     /
       DATA (A(I, 2, 2),I=1,10)  /
      +  .121    ,.133    ,.190    ,.234    ,.293    ,.332    ,.395    ,
      +  .431    ,.468    ,.502     /
       DATA (A(I, 2, 3),I=1,10)  /
      +  .227E-01,.374E-01,.474E-01,.578E-01,.722E-01,.794E-01,.960E-01,
      +  .102    ,.110    ,.120     /
       DATA (A(I, 2, 4),I=1,10)  /
      +  .287    ,.196    ,.270    ,.314    ,.373    ,.408    ,.462    ,
      +  .498    ,.529    ,.523     /
       DATA (A(I, 2, 5),I=1,10)  /
      +  .000E+00,.433E-01,.218E-01,.384E-01,.263E-01,.385E-01,.298E-01,
      +  .405E-01,.504E-01,.671E-01 /
       DATA (A(I, 2, 6),I=1,10)  /
      +  .000E+00,.151E-01,.177E-01,.159E-01,.173E-01,.173E-01,.187E-01,
      +  .196E-01,.201E-01,.191E-01 /
       DATA (A(I, 2, 7),I=1,10)  /
      +  .000E+00,.457E-02,.607E-02,.610E-02,.677E-02,.670E-02,.784E-02,
      +  .787E-02,.806E-02,.803E-02 /
       DATA (A(I, 2, 8),I=1,10)  /
      +  .000E+00,.000E+00,.702E-02,.536E-02,.558E-02,.510E-02,.554E-02,
      +  .546E-02,.538E-02,.489E-02 /
       DATA (A(I, 2, 9),I=1,10)  /
      +  .000E+00,.000E+00,.190E-02,.199E-02,.205E-02,.191E-02,.221E-02,
      +  .214E-02,.213E-02,.204E-02 /
       DATA (A(I, 2,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.226E-02,.219E-02,.195E-02,.208E-02,
      +  .204E-02,.203E-02,.194E-02 /
       DATA (A(I, 2,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.213E-02,.195E-02,.175E-02,.191E-02,
      +  .183E-02,.179E-02,.166E-02 /
       DATA (A(I, 2,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.588E-03,.186E-02,.137E-02,.141E-02,
      +  .128E-02,.117E-02,.947E-03 /
       DATA (A(I, 2,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.554E-03,.562E-03,.454E-03,
      +  .485E-03,.505E-03,.509E-03 /
       DATA (A(I, 2,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.490E-03,.533E-03,.531E-03,
      +  .476E-03,.437E-03,.369E-03 /
       DATA (A(I, 2,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.427E-03,.382E-03,
      +  .358E-03,.340E-03,.294E-03 /
       DATA (A(I, 2,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.239E-03,.298E-03,
      +  .238E-03,.196E-03,.134E-03 /
       DATA (A(I, 2,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.299E-04,.893E-04,
      +  .796E-04,.744E-04,.683E-04 /
       DATA (A(I, 2,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.127E-03,
      +  .107E-03,.916E-04,.720E-04 /
       DATA (A(I, 2,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.397E-04,
      +  .630E-04,.565E-04,.461E-04 /
       DATA (A(I, 2,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .511E-04,.459E-04,.402E-04 /
       DATA (A(I, 3, 1),I=1,10)  /
      +  .708    ,1.02    ,1.41    ,1.91    ,2.42    ,3.00    ,3.53    ,
      +  4.09    ,4.71    ,5.57     /
       DATA (A(I, 3, 2),I=1,10)  /
      +  .397    ,.410    ,.539    ,.648    ,.795    ,.910    ,1.06    ,
      +  1.17    ,1.29    ,1.42     /
       DATA (A(I, 3, 3),I=1,10)  /
      +  .845E-01,.122    ,.157    ,.190    ,.232    ,.262    ,.307    ,
      +  .335    ,.366    ,.402     /
       DATA (A(I, 3, 4),I=1,10)  /
      +  .210    ,.379    ,.450    ,.490    ,.574    ,.636    ,.709    ,
      +  .769    ,.820    ,.849     /
       DATA (A(I, 3, 5),I=1,10)  /
      +  .000E+00,.102    ,.675E-01,.104    ,.858E-01,.115    ,.102    ,
      +  .129    ,.154    ,.194     /
       DATA (A(I, 3, 6),I=1,10)  /
      +  .000E+00,.392E-01,.615E-01,.593E-01,.649E-01,.674E-01,.735E-01,
      +  .779E-01,.817E-01,.828E-01 /
       DATA (A(I, 3, 7),I=1,10)  /
      +  .000E+00,.539E-02,.222E-01,.238E-01,.269E-01,.280E-01,.320E-01,
      +  .334E-01,.350E-01,.361E-01 /
       DATA (A(I, 3, 8),I=1,10)  /
      +  .000E+00,.000E+00,.838E-02,.130E-01,.133E-01,.131E-01,.141E-01,
      +  .144E-01,.149E-01,.152E-01 /
       DATA (A(I, 3, 9),I=1,10)  /
      +  .000E+00,.000E+00,.228E-02,.647E-02,.688E-02,.687E-02,.772E-02,
      +  .786E-02,.811E-02,.824E-02 /
       DATA (A(I, 3,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.664E-02,.828E-02,.802E-02,.845E-02,
      +  .869E-02,.902E-02,.930E-02 /
       DATA (A(I, 3,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.338E-02,.735E-02,.710E-02,.767E-02,
      +  .767E-02,.776E-02,.756E-02 /
       DATA (A(I, 3,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.280E-03,.262E-02,.349E-02,.342E-02,
      +  .322E-02,.312E-02,.291E-02 /
       DATA (A(I, 3,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.618E-03,.161E-02,.138E-02,
      +  .148E-02,.155E-02,.166E-02 /
       DATA (A(I, 3,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.313E-03,.128E-02,.161E-02,
      +  .150E-02,.144E-02,.134E-02 /
       DATA (A(I, 3,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.645E-03,.118E-02,
      +  .115E-02,.111E-02,.103E-02 /
       DATA (A(I, 3,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.117E-03,.497E-03,
      +  .581E-03,.501E-03,.401E-03 /
       DATA (A(I, 3,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.115E-04,.997E-04,
      +  .202E-03,.203E-03,.206E-03 /
       DATA (A(I, 3,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.877E-04,
      +  .242E-03,.263E-03,.226E-03 /
       DATA (A(I, 3,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.158E-04,
      +  .881E-04,.152E-03,.136E-03 /
       DATA (A(I, 3,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .358E-04,.997E-04,.117E-03 /
       DATA (A(I, 4, 1),I=1,10)  /
      +  .945    ,1.29    ,1.40    ,1.98    ,2.73    ,3.17    ,3.77    ,
      +  4.29    ,4.78    ,5.54     /
       DATA (A(I, 4, 2),I=1,10)  /
      +  .581    ,.599    ,.645    ,.839    ,1.10    ,1.25    ,1.47    ,
      +  1.64    ,1.78    ,1.99     /
       DATA (A(I, 4, 3),I=1,10)  /
      +  .127    ,.182    ,.202    ,.264    ,.344    ,.387    ,.455    ,
      +  .504    ,.549    ,.611     /
       DATA (A(I, 4, 4),I=1,10)  /
      +  .183    ,.464    ,.351    ,.444    ,.642    ,.659    ,.772    ,
      +  .830    ,.882    ,.930     /
       DATA (A(I, 4, 5),I=1,10)  /
      +  .000E+00,.122    ,.803E-01,.136    ,.134    ,.173    ,.164    ,
      +  .203    ,.239    ,.300     /
       DATA (A(I, 4, 6),I=1,10)  /
      +  .000E+00,.393E-01,.766E-01,.872E-01,.108    ,.111    ,.123    ,
      +  .132    ,.139    ,.145     /
       DATA (A(I, 4, 7),I=1,10)  /
      +  .000E+00,.416E-02,.289E-01,.360E-01,.454E-01,.477E-01,.549E-01,
      +  .583E-01,.618E-01,.654E-01 /
       DATA (A(I, 4, 8),I=1,10)  /
      +  .000E+00,.000E+00,.761E-02,.157E-01,.214E-01,.205E-01,.233E-01,
      +  .241E-01,.255E-01,.271E-01 /
       DATA (A(I, 4, 9),I=1,10)  /
      +  .000E+00,.000E+00,.238E-02,.803E-02,.123E-01,.123E-01,.140E-01,
      +  .145E-01,.153E-01,.160E-01 /
       DATA (A(I, 4,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.695E-02,.150E-01,.154E-01,.166E-01,
      +  .172E-01,.181E-01,.192E-01 /
       DATA (A(I, 4,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.355E-02,.104E-01,.143E-01,.156E-01,
      +  .158E-01,.164E-01,.165E-01 /
       DATA (A(I, 4,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.112E-03,.276E-02,.568E-02,.736E-02,
      +  .684E-02,.691E-02,.661E-02 /
       DATA (A(I, 4,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.740E-03,.222E-02,.339E-02,
      +  .352E-02,.382E-02,.409E-02 /
       DATA (A(I, 4,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.369E-03,.160E-02,.322E-02,
      +  .375E-02,.375E-02,.355E-02 /
       DATA (A(I, 4,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.750E-03,.190E-02,
      +  .298E-02,.319E-02,.299E-02 /
       DATA (A(I, 4,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.260E-03,.673E-03,
      +  .117E-02,.156E-02,.126E-02 /
       DATA (A(I, 4,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.283E-05,.131E-03,
      +  .363E-03,.618E-03,.690E-03 /
       DATA (A(I, 4,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.205E-03,
      +  .378E-03,.709E-03,.844E-03 /
       DATA (A(I, 4,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.654E-05,
      +  .150E-03,.341E-03,.527E-03 /
       DATA (A(I, 4,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .957E-04,.197E-03,.406E-03 /
       DATA (A(I, 5, 1),I=1,10)  /
      +  1.16    ,1.70    ,2.19    ,2.79    ,3.33    ,3.90    ,4.49    ,
      +  5.07    ,5.66    ,6.38     /
       DATA (A(I, 5, 2),I=1,10)  /
      +  .779    ,.899    ,1.09    ,1.28    ,1.51    ,1.71    ,1.96    ,
      +  2.18    ,2.39    ,2.62     /
       DATA (A(I, 5, 3),I=1,10)  /
      +  .167    ,.263    ,.334    ,.408    ,.482    ,.548    ,.632    ,
      +  .700    ,.767    ,.840     /
       DATA (A(I, 5, 4),I=1,10)  /
      +  .203    ,.565    ,.845    ,.867    ,.906    ,.961    ,1.08    ,
      +  1.13    ,1.21    ,1.25     /
       DATA (A(I, 5, 5),I=1,10)  /
      +  .000E+00,.129    ,.152    ,.237    ,.208    ,.268    ,.258    ,
      +  .312    ,.368    ,.450     /
       DATA (A(I, 5, 6),I=1,10)  /
      +  .000E+00,.460E-01,.126    ,.174    ,.182    ,.188    ,.208    ,
      +  .219    ,.233    ,.239     /
       DATA (A(I, 5, 7),I=1,10)  /
      +  .000E+00,.289E-02,.380E-01,.611E-01,.788E-01,.845E-01,.974E-01,
      +  .103    ,.111    ,.117     /
       DATA (A(I, 5, 8),I=1,10)  /
      +  .000E+00,.000E+00,.137E-01,.223E-01,.374E-01,.436E-01,.488E-01,
      +  .488E-01,.524E-01,.547E-01 /
       DATA (A(I, 5, 9),I=1,10)  /
      +  .000E+00,.000E+00,.162E-02,.114E-01,.198E-01,.263E-01,.315E-01,
      +  .323E-01,.348E-01,.364E-01 /
       DATA (A(I, 5,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.149E-01,.240E-01,.320E-01,.428E-01,
      +  .436E-01,.469E-01,.493E-01 /
       DATA (A(I, 5,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.562E-02,.194E-01,.290E-01,.408E-01,
      +  .460E-01,.492E-01,.500E-01 /
       DATA (A(I, 5,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.476E-04,.106E-01,.134E-01,.191E-01,
      +  .227E-01,.264E-01,.253E-01 /
       DATA (A(I, 5,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.281E-02,.679E-02,.879E-02,
      +  .123E-01,.165E-01,.190E-01 /
       DATA (A(I, 5,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.542E-04,.847E-02,.125E-01,
      +  .144E-01,.173E-01,.192E-01 /
       DATA (A(I, 5,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.652E-02,.982E-02,
      +  .129E-01,.159E-01,.192E-01 /
       DATA (A(I, 5,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.109E-03,.688E-02,
      +  .751E-02,.845E-02,.905E-02 /
       DATA (A(I, 5,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.823E-06,.237E-02,
      +  .318E-02,.446E-02,.569E-02 /
       DATA (A(I, 5,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.604E-03,
      +  .610E-02,.673E-02,.827E-02 /
       DATA (A(I, 5,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.716E-06,
      +  .412E-02,.519E-02,.617E-02 /
       DATA (A(I, 5,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .710E-03,.543E-02,.674E-02 /
       DATA (A(I, 6, 1),I=1,10)  /
      +  1.36    ,2.08    ,2.67    ,3.30    ,3.94    ,4.62    ,5.18    ,
      +  3.60    ,3.64    ,3.95     /
       DATA (A(I, 6, 2),I=1,10)  /
      +  1.07    ,1.33    ,1.58    ,1.82    ,2.10    ,2.44    ,2.74    ,
      +  1.78    ,1.73    ,1.80     /
       DATA (A(I, 6, 3),I=1,10)  /
      +  .158    ,.276    ,.402    ,.506    ,.609    ,.700    ,.802    ,
      +  .638    ,.629    ,.658     /
       DATA (A(I, 6, 4),I=1,10)  /
      +  .308    ,.739    ,1.02    ,1.12    ,1.26    ,1.35    ,1.57    ,
      +  1.94    ,1.71    ,1.55     /
       DATA (A(I, 6, 5),I=1,10)  /
      +  .000E+00,.217    ,.183    ,.324    ,.276    ,.395    ,.393    ,
      +  .558    ,.602    ,.681     /
       DATA (A(I, 6, 6),I=1,10)  /
      +  .000E+00,.658E-01,.251    ,.267    ,.299    ,.326    ,.386    ,
      +  .452    ,.475    ,.409     /
       DATA (A(I, 6, 7),I=1,10)  /
      +  .000E+00,.198E-02,.774E-01,.136    ,.149    ,.164    ,.187    ,
      +  .210    ,.238    ,.256     /
       DATA (A(I, 6, 8),I=1,10)  /
      +  .000E+00,.000E+00,.290E-01,.122    ,.139    ,.128    ,.129    ,
      +  .137    ,.147    ,.167     /
       DATA (A(I, 6, 9),I=1,10)  /
      +  .000E+00,.000E+00,.699E-03,.617E-01,.750E-01,.801E-01,.905E-01,
      +  .974E-01,.105    ,.122     /
       DATA (A(I, 6,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.310E-01,.112    ,.127    ,.140    ,
      +  .143    ,.155    ,.176     /
       DATA (A(I, 6,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.277E-02,.889E-01,.143    ,.150    ,
      +  .175    ,.184    ,.208     /
       DATA (A(I, 6,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.202E-04,.343E-01,.959E-01,.109    ,
      +  .115    ,.112    ,.116     /
       DATA (A(I, 6,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.186E-02,.435E-01,.512E-01,
      +  .744E-01,.856E-01,.103     /
       DATA (A(I, 6,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.144E-04,.427E-01,.786E-01,
      +  .911E-01,.993E-01,.108     /
       DATA (A(I, 6,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.466E-02,.518E-01,
      +  .848E-01,.109    ,.119     /
       DATA (A(I, 6,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.655E-05,.330E-01,
      +  .586E-01,.617E-01,.594E-01 /
       DATA (A(I, 6,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.228E-06,.328E-02,
      +  .190E-01,.301E-01,.454E-01 /
       DATA (A(I, 6,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.218E-04,
      +  .272E-01,.501E-01,.707E-01 /
       DATA (A(I, 6,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.146E-06,
      +  .441E-02,.378E-01,.556E-01 /
       DATA (A(I, 6,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .160E-03,.204E-01,.679E-01 /
       DATA (A(I, 7, 1),I=1,10)  /
      +  .522    ,.862    ,1.14    ,1.40    ,1.70    ,1.94    ,2.26    ,
      +  2.48    ,2.72    ,3.95     /
       DATA (A(I, 7, 2),I=1,10)  /
      +  .314    ,.450    ,.588    ,.692    ,.834    ,.936    ,1.09    ,
      +  1.18    ,1.28    ,1.80     /
       DATA (A(I, 7, 3),I=1,10)  /
      +  .814E-01,.147    ,.189    ,.226    ,.272    ,.302    ,.351    ,
      +  .378    ,.406    ,.658     /
       DATA (A(I, 7, 4),I=1,10)  /
      +  .252    ,.864    ,1.01    ,.851    ,.837    ,.774    ,.763    ,
      +  .757    ,.748    ,1.55     /
       DATA (A(I, 7, 5),I=1,10)  /
      +  .000E+00,.225    ,.180    ,.276    ,.193    ,.240    ,.190    ,
      +  .228    ,.259    ,.681     /
       DATA (A(I, 7, 6),I=1,10)  /
      +  .000E+00,.485E-01,.272    ,.273    ,.253    ,.216    ,.206    ,
      +  .197    ,.191    ,.409     /
       DATA (A(I, 7, 7),I=1,10)  /
      +  .000E+00,.137E-02,.752E-01,.137    ,.152    ,.134    ,.125    ,
      +  .119    ,.116    ,.256     /
       DATA (A(I, 7, 8),I=1,10)  /
      +  .000E+00,.000E+00,.220E-01,.155    ,.175    ,.155    ,.116    ,
      +  .977E-01,.858E-01,.167     /
       DATA (A(I, 7, 9),I=1,10)  /
      +  .000E+00,.000E+00,.326E-03,.695E-01,.881E-01,.106    ,.897E-01,
      +  .782E-01,.706E-01,.122     /
       DATA (A(I, 7,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.261E-01,.124    ,.131    ,.156    ,
      +  .141    ,.121    ,.176     /
       DATA (A(I, 7,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.785E-03,.864E-01,.130    ,.170    ,
      +  .182    ,.172    ,.208     /
       DATA (A(I, 7,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.896E-05,.225E-01,.105    ,.126    ,
      +  .126    ,.135    ,.116     /
       DATA (A(I, 7,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.542E-03,.427E-01,.553E-01,
      +  .744E-01,.980E-01,.103     /
       DATA (A(I, 7,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.515E-05,.377E-01,.831E-01,
      +  .985E-01,.104    ,.108     /
       DATA (A(I, 7,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.285E-02,.495E-01,
      +  .871E-01,.106    ,.119     /
       DATA (A(I, 7,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.110E-05,.284E-01,
      +  .588E-01,.657E-01,.594E-01 /
       DATA (A(I, 7,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.722E-07,.176E-02,
      +  .170E-01,.305E-01,.454E-01 /
       DATA (A(I, 7,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.148E-05,
      +  .213E-01,.492E-01,.707E-01 /
       DATA (A(I, 7,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.323E-07,
      +  .722E-02,.359E-01,.556E-01 /
       DATA (A(I, 7,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .461E-05,.155E-01,.679E-01 /
       DATA (A(I, 8, 1),I=1,10)  /
      +  .630    ,.974    ,1.29    ,1.58    ,1.89    ,2.16    ,2.49    ,
      +  2.75    ,3.02    ,3.95     /
       DATA (A(I, 8, 2),I=1,10)  /
      +  .328    ,.459    ,.613    ,.735    ,.879    ,.994    ,1.15    ,
      +  1.27    ,1.38    ,1.80     /
       DATA (A(I, 8, 3),I=1,10)  /
      +  .748E-01,.121    ,.164    ,.197    ,.235    ,.265    ,.310    ,
      +  .339    ,.370    ,.658     /
       DATA (A(I, 8, 4),I=1,10)  /
      +  .194    ,.211    ,.337    ,.344    ,.339    ,.351    ,.390    ,
      +  .419    ,.442    ,1.55     /
       DATA (A(I, 8, 5),I=1,10)  /
      +  .000E+00,.869E-01,.725E-01,.113    ,.810E-01,.106    ,.951E-01,
      +  .120    ,.143    ,.681     /
       DATA (A(I, 8, 6),I=1,10)  /
      +  .000E+00,.288E-01,.102    ,.922E-01,.857E-01,.845E-01,.932E-01,
      +  .983E-01,.102    ,.409     /
       DATA (A(I, 8, 7),I=1,10)  /
      +  .000E+00,.668E-03,.533E-01,.575E-01,.493E-01,.482E-01,.539E-01,
      +  .558E-01,.582E-01,.256     /
       DATA (A(I, 8, 8),I=1,10)  /
      +  .000E+00,.000E+00,.205E-01,.808E-01,.510E-01,.409E-01,.406E-01,
      +  .394E-01,.389E-01,.167     /
       DATA (A(I, 8, 9),I=1,10)  /
      +  .000E+00,.000E+00,.999E-04,.647E-01,.385E-01,.325E-01,.325E-01,
      +  .316E-01,.314E-01,.122     /
       DATA (A(I, 8,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.169E-01,.834E-01,.611E-01,.565E-01,
      +  .533E-01,.519E-01,.176     /
       DATA (A(I, 8,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.107E-03,.769E-01,.922E-01,.805E-01,
      +  .745E-01,.711E-01,.208     /
       DATA (A(I, 8,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.180E-05,.143E-01,.983E-01,.775E-01,
      +  .627E-01,.541E-01,.116     /
       DATA (A(I, 8,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.157E-04,.346E-01,.507E-01,
      +  .479E-01,.455E-01,.103     /
       DATA (A(I, 8,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.752E-06,.248E-01,.721E-01,
      +  .728E-01,.611E-01,.108     /
       DATA (A(I, 8,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.686E-04,.356E-01,
      +  .731E-01,.791E-01,.119     /
       DATA (A(I, 8,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.838E-07,.151E-01,
      +  .470E-01,.567E-01,.594E-01 /
       DATA (A(I, 8,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.759E-08,.400E-04,
      +  .193E-01,.313E-01,.454E-01 /
       DATA (A(I, 8,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.385E-07,
      +  .921E-02,.353E-01,.707E-01 /
       DATA (A(I, 8,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.219E-08,
      +  .348E-03,.226E-01,.556E-01 /
       DATA (A(I, 8,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .212E-07,.149E-01,.679E-01 /
       DATA (A(I, 9, 1),I=1,10)  /
      +  .736    ,1.13    ,1.49    ,1.82    ,2.20    ,2.49    ,2.86    ,
      +  3.17    ,3.49    ,3.95     /
       DATA (A(I, 9, 2),I=1,10)  /
      +  .339    ,.492    ,.658    ,.789    ,.958    ,1.08    ,1.25    ,
      +  1.37    ,1.50    ,1.80     /
       DATA (A(I, 9, 3),I=1,10)  /
      +  .680E-01,.110    ,.150    ,.180    ,.222    ,.247    ,.289    ,
      +  .318    ,.349    ,.658     /
       DATA (A(I, 9, 4),I=1,10)  /
      +  .110    ,.104    ,.157    ,.156    ,.210    ,.205    ,.246    ,
      +  .274    ,.300    ,1.55     /
       DATA (A(I, 9, 5),I=1,10)  /
      +  .000E+00,.379E-01,.347E-01,.477E-01,.486E-01,.576E-01,.569E-01,
      +  .732E-01,.893E-01,.681     /
       DATA (A(I, 9, 6),I=1,10)  /
      +  .000E+00,.223E-01,.354E-01,.312E-01,.436E-01,.400E-01,.489E-01,
      +  .548E-01,.600E-01,.409     /
       DATA (A(I, 9, 7),I=1,10)  /
      +  .000E+00,.338E-03,.149E-01,.142E-01,.215E-01,.188E-01,.248E-01,
      +  .278E-01,.307E-01,.256     /
       DATA (A(I, 9, 8),I=1,10)  /
      +  .000E+00,.000E+00,.553E-02,.862E-02,.150E-01,.106E-01,.145E-01,
      +  .165E-01,.181E-01,.167     /
       DATA (A(I, 9, 9),I=1,10)  /
      +  .000E+00,.000E+00,.375E-04,.641E-02,.111E-01,.792E-02,.112E-01,
      +  .127E-01,.140E-01,.122     /
       DATA (A(I, 9,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.112E-01,.200E-01,.127E-01,.176E-01,
      +  .200E-01,.220E-01,.176     /
       DATA (A(I, 9,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.244E-04,.261E-01,.162E-01,.232E-01,
      +  .263E-01,.287E-01,.208     /
       DATA (A(I, 9,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.455E-06,.635E-02,.121E-01,.186E-01,
      +  .201E-01,.207E-01,.116     /
       DATA (A(I, 9,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.146E-05,.922E-02,.116E-01,
      +  .145E-01,.165E-01,.103     /
       DATA (A(I, 9,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.135E-06,.128E-01,.202E-01,
      +  .215E-01,.220E-01,.108     /
       DATA (A(I, 9,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.237E-05,.229E-01,
      +  .259E-01,.271E-01,.119     /
       DATA (A(I, 9,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.100E-07,.534E-02,
      +  .210E-01,.193E-01,.594E-01 /
       DATA (A(I, 9,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.915E-09,.847E-06,
      +  .119E-01,.125E-01,.454E-01 /
       DATA (A(I, 9,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.298E-08,
      +  .101E-01,.242E-01,.707E-01 /
       DATA (A(I, 9,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.196E-09,
      +  .243E-05,.234E-01,.556E-01 /
       DATA (A(I, 9,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .575E-09,.364E-02,.679E-01 /
       DATA (A(I,10, 1),I=1,10)  /
      +  .959    ,1.46    ,1.92    ,2.34    ,2.80    ,3.24    ,3.64    ,
      +  4.05    ,4.48    ,3.95     /
       DATA (A(I,10, 2),I=1,10)  /
      +  .343    ,.516    ,.692    ,.836    ,1.01    ,1.16    ,1.31    ,
      +  1.46    ,1.61    ,1.80     /
       DATA (A(I,10, 3),I=1,10)  /
      +  .512E-01,.837E-01,.115    ,.138    ,.169    ,.195    ,.220    ,
      +  .245    ,.270    ,.658     /
       DATA (A(I,10, 4),I=1,10)  /
      +  .274E-01,.361E-01,.510E-01,.562E-01,.703E-01,.828E-01,.877E-01,
      +  .996E-01,.111    ,1.55     /
       DATA (A(I,10, 5),I=1,10)  /
      +  .000E+00,.850E-02,.875E-02,.118E-01,.124E-01,.170E-01,.154E-01,
      +  .194E-01,.237E-01,.681     /
       DATA (A(I,10, 6),I=1,10)  /
      +  .000E+00,.345E-02,.519E-02,.533E-02,.691E-02,.842E-02,.844E-02,
      +  .987E-02,.113E-01,.409     /
       DATA (A(I,10, 7),I=1,10)  /
      +  .000E+00,.722E-04,.130E-02,.135E-02,.189E-02,.240E-02,.235E-02,
      +  .281E-02,.331E-02,.256     /
       DATA (A(I,10, 8),I=1,10)  /
      +  .000E+00,.000E+00,.283E-03,.272E-03,.394E-03,.557E-03,.480E-03,
      +  .616E-03,.775E-03,.167     /
       DATA (A(I,10, 9),I=1,10)  /
      +  .000E+00,.000E+00,.457E-05,.122E-03,.192E-03,.275E-03,.225E-03,
      +  .292E-03,.373E-03,.122     /
       DATA (A(I,10,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.119E-03,.185E-03,.278E-03,.201E-03,
      +  .274E-03,.364E-03,.176     /
       DATA (A(I,10,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.140E-05,.129E-03,.200E-03,.137E-03,
      +  .188E-03,.252E-03,.208     /
       DATA (A(I,10,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.207E-07,.307E-04,.518E-04,.278E-04,
      +  .421E-04,.608E-04,.116     /
       DATA (A(I,10,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.306E-07,.252E-04,.111E-04,
      +  .188E-04,.295E-04,.103     /
       DATA (A(I,10,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.321E-08,.220E-04,.104E-04,
      +  .162E-04,.243E-04,.108     /
       DATA (A(I,10,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.770E-08,.632E-05,
      +  .105E-04,.162E-04,.119     /
       DATA (A(I,10,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.117E-09,.199E-05,
      +  .321E-05,.492E-05,.594E-01 /
       DATA (A(I,10,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.888E-11,.323E-09,
      +  .106E-05,.192E-05,.454E-01 /
       DATA (A(I,10,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.174E-10,
      +  .131E-05,.218E-05,.707E-01 /
       DATA (A(I,10,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.994E-12,
      +  .233E-09,.104E-05,.556E-01 /
       DATA (A(I,10,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  .144E-11,.724E-06,.679E-01 /
       DATA (AE(I, 1, 1),I=1,10)  /
      +  7.27    ,6.29    ,7.76    ,6.70    ,8.17    ,7.34    ,8.70    ,
      +  8.02    ,7.37    ,6.18     /
       DATA (AE(I, 1, 2),I=1,10)  /
      +  7.41    ,7.52    ,8.14    ,8.20    ,8.96    ,9.05    ,9.96    ,
      +  10.0    ,10.1    ,9.86     /
       DATA (AE(I, 1, 3),I=1,10)  /
      +  7.72    ,7.69    ,9.17    ,8.99    ,10.6    ,10.5    ,12.1    ,
      +  12.1    ,12.0    ,11.5     /
       DATA (AE(I, 1, 4),I=1,10)  /
      +  7.90    ,8.48    ,9.50    ,9.94    ,10.8    ,11.4    ,12.2    ,
      +  12.8    ,13.3    ,13.8     /
       DATA (AE(I, 1, 5),I=1,10)  /
      +  .000E+00,8.52    ,9.59    ,10.1    ,11.1    ,11.8    ,12.7    ,
      +  13.3    ,13.8    ,14.4     /
       DATA (AE(I, 1, 6),I=1,10)  /
      +  .000E+00,9.00    ,10.7    ,11.7    ,13.2    ,14.2    ,15.6    ,
      +  16.5    ,17.3    ,18.0     /
       DATA (AE(I, 1, 7),I=1,10)  /
      +  .000E+00,9.01    ,11.1    ,11.9    ,14.3    ,15.0    ,17.4    ,
      +  18.0    ,18.6    ,18.8     /
       DATA (AE(I, 1, 8),I=1,10)  /
      +  .000E+00,.000E+00,11.2    ,12.4    ,14.5    ,15.7    ,17.6    ,
      +  18.8    ,19.9    ,20.9     /
       DATA (AE(I, 1, 9),I=1,10)  /
      +  .000E+00,.000E+00,11.4    ,12.7    ,15.5    ,16.6    ,19.3    ,
      +  20.2    ,21.1    ,21.7     /
       DATA (AE(I, 1,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,13.2    ,15.8    ,17.3    ,19.9    ,
      +  21.2    ,22.4    ,23.2     /
       DATA (AE(I, 1,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,13.2    ,16.3    ,17.8    ,20.8    ,
      +  22.1    ,23.3    ,24.2     /
       DATA (AE(I, 1,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,13.4    ,16.2    ,18.2    ,21.0    ,
      +  22.8    ,24.4    ,25.9     /
       DATA (AE(I, 1,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,16.5    ,18.4    ,21.6    ,
      +  23.2    ,24.8    ,26.2     /
       DATA (AE(I, 1,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,16.7    ,19.0    ,22.3    ,
      +  24.3    ,26.1    ,27.4     /
       DATA (AE(I, 1,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,19.1    ,22.8    ,
      +  24.7    ,26.6    ,28.2     /
       DATA (AE(I, 1,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,19.2    ,23.0    ,
      +  25.3    ,27.5    ,29.5     /
       DATA (AE(I, 1,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,19.6    ,23.3    ,
      +  25.6    ,27.8    ,29.6     /
       DATA (AE(I, 1,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,23.6    ,
      +  26.2    ,28.5    ,30.4     /
       DATA (AE(I, 1,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,23.7    ,
      +  26.3    ,28.8    ,31.0     /
       DATA (AE(I, 1,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  26.5    ,29.2    ,31.5     /
       DATA (AE(I, 2, 1),I=1,10)  /
      +  8.74    ,8.16    ,9.25    ,8.45    ,9.46    ,8.90    ,9.83    ,
      +  9.38    ,8.96    ,8.15     /
       DATA (AE(I, 2, 2),I=1,10)  /
      +  8.96    ,9.30    ,9.95    ,10.0    ,10.8    ,10.9    ,11.7    ,
      +  11.8    ,11.9    ,11.8     /
       DATA (AE(I, 2, 3),I=1,10)  /
      +  9.44    ,9.66    ,11.0    ,11.0    ,12.3    ,12.5    ,13.7    ,
      +  13.9    ,14.0    ,13.8     /
       DATA (AE(I, 2, 4),I=1,10)  /
      +  8.86    ,9.81    ,10.8    ,11.2    ,12.0    ,12.6    ,13.4    ,
      +  14.0    ,14.5    ,15.1     /
       DATA (AE(I, 2, 5),I=1,10)  /
      +  .000E+00,10.2    ,11.4    ,12.0    ,12.9    ,13.6    ,14.5    ,
      +  15.1    ,15.7    ,16.3     /
       DATA (AE(I, 2, 6),I=1,10)  /
      +  .000E+00,10.7    ,12.5    ,13.5    ,15.1    ,16.0    ,17.5    ,
      +  18.3    ,19.2    ,19.9     /
       DATA (AE(I, 2, 7),I=1,10)  /
      +  .000E+00,11.5    ,12.9    ,13.9    ,16.1    ,17.0    ,19.1    ,
      +  19.8    ,20.6    ,21.0     /
       DATA (AE(I, 2, 8),I=1,10)  /
      +  .000E+00,.000E+00,12.4    ,13.8    ,15.9    ,17.2    ,19.1    ,
      +  20.3    ,21.4    ,22.3     /
       DATA (AE(I, 2, 9),I=1,10)  /
      +  .000E+00,.000E+00,13.4    ,14.5    ,17.1    ,18.3    ,20.9    ,
      +  21.9    ,23.0    ,23.7     /
       DATA (AE(I, 2,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,14.9    ,17.5    ,19.1    ,21.6    ,
      +  22.9    ,24.1    ,25.0     /
       DATA (AE(I, 2,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,15.0    ,18.0    ,19.6    ,22.4    ,
      +  23.8    ,25.2    ,26.2     /
       DATA (AE(I, 2,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,16.2    ,17.3    ,19.4    ,22.2    ,
      +  24.0    ,25.7    ,27.2     /
       DATA (AE(I, 2,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,17.8    ,19.8    ,22.9    ,
      +  24.6    ,26.2    ,27.7     /
       DATA (AE(I, 2,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,19.1    ,20.4    ,23.7    ,
      +  25.7    ,27.6    ,29.1     /
       DATA (AE(I, 2,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,20.5    ,24.1    ,
      +  26.1    ,28.1    ,29.9     /
       DATA (AE(I, 2,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,20.9    ,23.9    ,
      +  26.4    ,28.7    ,30.7     /
       DATA (AE(I, 2,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,22.4    ,24.2    ,
      +  26.7    ,29.0    ,30.9     /
       DATA (AE(I, 2,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,24.8    ,
      +  27.3    ,29.7    ,31.8     /
       DATA (AE(I, 2,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,26.1    ,
      +  27.3    ,29.9    ,32.3     /
       DATA (AE(I, 2,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  27.4    ,30.1    ,32.6     /
       DATA (AE(I, 3, 1),I=1,10)  /
      +  11.0    ,11.0    ,11.7    ,11.3    ,11.9    ,11.4    ,12.1    ,
      +  11.7    ,11.5    ,11.0     /
       DATA (AE(I, 3, 2),I=1,10)  /
      +  11.2    ,12.0    ,12.7    ,12.9    ,13.6    ,13.7    ,14.4    ,
      +  14.6    ,14.7    ,14.6     /
       DATA (AE(I, 3, 3),I=1,10)  /
      +  12.1    ,12.6    ,13.7    ,13.9    ,15.0    ,15.2    ,16.3    ,
      +  16.5    ,16.7    ,16.7     /
       DATA (AE(I, 3, 4),I=1,10)  /
      +  12.6    ,11.3    ,12.4    ,13.0    ,13.8    ,14.2    ,15.0    ,
      +  15.6    ,16.1    ,16.6     /
       DATA (AE(I, 3, 5),I=1,10)  /
      +  .000E+00,12.6    ,13.7    ,14.4    ,15.3    ,16.0    ,16.8    ,
      +  17.5    ,18.1    ,18.6     /
       DATA (AE(I, 3, 6),I=1,10)  /
      +  .000E+00,14.0    ,14.6    ,15.8    ,17.4    ,18.4    ,19.8    ,
      +  20.6    ,21.5    ,22.2     /
       DATA (AE(I, 3, 7),I=1,10)  /
      +  .000E+00,16.0    ,15.2    ,16.3    ,18.3    ,19.3    ,21.1    ,
      +  22.0    ,22.8    ,23.5     /
       DATA (AE(I, 3, 8),I=1,10)  /
      +  .000E+00,.000E+00,15.6    ,15.1    ,17.2    ,18.6    ,20.6    ,
      +  21.8    ,22.9    ,23.8     /
       DATA (AE(I, 3, 9),I=1,10)  /
      +  .000E+00,.000E+00,17.8    ,16.3    ,18.8    ,20.1    ,22.5    ,
      +  23.6    ,24.7    ,25.6     /
       DATA (AE(I, 3,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,17.5    ,19.0    ,20.7    ,23.1    ,
      +  24.5    ,25.8    ,26.8     /
       DATA (AE(I, 3,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,19.2    ,19.4    ,21.1    ,23.8    ,
      +  25.4    ,26.8    ,28.0     /
       DATA (AE(I, 3,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,20.7    ,19.6    ,19.7    ,22.4    ,
      +  24.4    ,26.2    ,27.9     /
       DATA (AE(I, 3,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,21.6    ,20.4    ,23.2    ,
      +  25.1    ,26.9    ,28.5     /
       DATA (AE(I, 3,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,23.5    ,22.0    ,23.8    ,
      +  26.1    ,28.1    ,29.9     /
       DATA (AE(I, 3,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,23.7    ,24.2    ,
      +  26.3    ,28.5    ,30.4     /
       DATA (AE(I, 3,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,25.4    ,24.8    ,
      +  25.6    ,28.1    ,30.5     /
       DATA (AE(I, 3,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,26.9    ,26.8    ,
      +  26.1    ,28.4    ,30.8     /
       DATA (AE(I, 3,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,28.8    ,
      +  27.6    ,29.0    ,31.5     /
       DATA (AE(I, 3,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,30.5    ,
      +  29.2    ,28.9    ,31.5     /
       DATA (AE(I, 3,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  31.0    ,30.0    ,31.7     /
       DATA (AE(I, 4, 1),I=1,10)  /
      +  13.0    ,13.2    ,14.8    ,14.2    ,14.2    ,14.1    ,14.5    ,
      +  14.4    ,14.3    ,14.0     /
       DATA (AE(I, 4, 2),I=1,10)  /
      +  13.5    ,14.5    ,16.1    ,15.9    ,16.0    ,16.3    ,16.8    ,
      +  17.0    ,17.1    ,17.2     /
       DATA (AE(I, 4, 3),I=1,10)  /
      +  14.9    ,15.3    ,17.2    ,17.1    ,17.5    ,17.8    ,18.6    ,
      +  18.9    ,19.1    ,19.3     /
       DATA (AE(I, 4, 4),I=1,10)  /
      +  15.1    ,13.5    ,16.4    ,16.7    ,16.4    ,17.3    ,17.8    ,
      +  18.5    ,19.0    ,19.6     /
       DATA (AE(I, 4, 5),I=1,10)  /
      +  .000E+00,15.6    ,17.5    ,17.7    ,17.8    ,18.6    ,19.2    ,
      +  19.9    ,20.3    ,21.1     /
       DATA (AE(I, 4, 6),I=1,10)  /
      +  .000E+00,18.0    ,18.4    ,19.2    ,19.8    ,20.9    ,22.0    ,
      +  23.1    ,23.6    ,24.7     /
       DATA (AE(I, 4, 7),I=1,10)  /
      +  .000E+00,27.4    ,19.1    ,19.8    ,20.7    ,21.8    ,23.2    ,
      +  24.4    ,24.9    ,25.9     /
       DATA (AE(I, 4, 8),I=1,10)  /
      +  .000E+00,.000E+00,18.9    ,18.9    ,19.3    ,21.1    ,22.5    ,
      +  24.0    ,24.7    ,26.0     /
       DATA (AE(I, 4, 9),I=1,10)  /
      +  .000E+00,.000E+00,21.1    ,19.7    ,20.7    ,22.3    ,24.0    ,
      +  25.6    ,26.3    ,27.7     /
       DATA (AE(I, 4,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,21.0    ,21.1    ,22.9    ,24.6    ,
      +  26.5    ,27.3    ,29.0     /
       DATA (AE(I, 4,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,21.3    ,22.4    ,23.1    ,25.0    ,
      +  27.1    ,27.9    ,29.8     /
       DATA (AE(I, 4,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,36.6    ,21.5    ,22.2    ,23.1    ,
      +  25.6    ,26.8    ,29.1     /
       DATA (AE(I, 4,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,22.9    ,23.1    ,23.7    ,
      +  26.2    ,27.3    ,29.6     /
       DATA (AE(I, 4,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,30.5    ,23.6    ,25.0    ,
      +  26.9    ,28.2    ,30.7     /
       DATA (AE(I, 4,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,25.4    ,26.2    ,
      +  27.2    ,28.3    ,31.0     /
       DATA (AE(I, 4,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,24.5    ,25.9    ,
      +  27.4    ,27.6    ,30.7     /
       DATA (AE(I, 4,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,43.3    ,28.4    ,
      +  27.5    ,27.9    ,30.9     /
       DATA (AE(I, 4,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,27.2    ,
      +  29.1    ,29.0    ,31.4     /
       DATA (AE(I, 4,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,51.3    ,
      +  30.6    ,29.5    ,31.4     /
       DATA (AE(I, 4,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  28.8    ,30.6    ,32.4     /
       DATA (AE(I, 5, 1),I=1,10)  /
      +  15.0    ,14.9    ,15.5    ,15.4    ,15.9    ,15.8    ,16.2    ,
      +  16.2    ,16.1    ,15.9     /
       DATA (AE(I, 5, 2),I=1,10)  /
      +  15.4    ,16.1    ,17.0    ,17.4    ,18.0    ,18.2    ,18.7    ,
      +  18.9    ,19.0    ,19.1     /
       DATA (AE(I, 5, 3),I=1,10)  /
      +  17.1    ,17.2    ,18.3    ,18.7    ,19.3    ,19.6    ,20.3    ,
      +  20.6    ,20.8    ,20.9     /
       DATA (AE(I, 5, 4),I=1,10)  /
      +  14.7    ,14.8    ,15.0    ,16.0    ,17.0    ,17.7    ,18.1    ,
      +  19.0    ,19.4    ,20.0     /
       DATA (AE(I, 5, 5),I=1,10)  /
      +  .000E+00,16.7    ,17.6    ,18.1    ,18.6    ,19.2    ,19.7    ,
      +  20.4    ,20.8    ,21.2     /
       DATA (AE(I, 5, 6),I=1,10)  /
      +  .000E+00,17.8    ,18.2    ,19.2    ,20.0    ,21.0    ,21.9    ,
      +  23.0    ,23.6    ,24.3     /
       DATA (AE(I, 5, 7),I=1,10)  /
      +  .000E+00,35.2    ,18.9    ,20.3    ,20.6    ,21.5    ,22.6    ,
      +  23.7    ,24.2    ,24.7     /
       DATA (AE(I, 5, 8),I=1,10)  /
      +  .000E+00,.000E+00,16.4    ,18.9    ,18.8    ,19.6    ,20.7    ,
      +  22.3    ,23.1    ,23.9     /
       DATA (AE(I, 5, 9),I=1,10)  /
      +  .000E+00,.000E+00,33.9    ,19.8    ,20.3    ,20.7    ,21.9    ,
      +  23.4    ,24.1    ,24.8     /
       DATA (AE(I, 5,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,18.0    ,20.0    ,21.4    ,22.0    ,
      +  23.8    ,24.6    ,25.4     /
       DATA (AE(I, 5,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,26.4    ,20.4    ,21.2    ,22.3    ,
      +  23.8    ,24.7    ,25.5     /
       DATA (AE(I, 5,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,41.7    ,18.2    ,19.8    ,21.1    ,
      +  22.6    ,23.4    ,24.6     /
       DATA (AE(I, 5,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,22.5    ,20.0    ,21.7    ,
      +  22.8    ,23.7    ,24.7     /
       DATA (AE(I, 5,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,54.1    ,19.9    ,21.9    ,
      +  23.2    ,24.3    ,25.3     /
       DATA (AE(I, 5,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,21.2    ,22.2    ,
      +  23.6    ,24.9    ,25.5     /
       DATA (AE(I, 5,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,44.9    ,21.9    ,
      +  23.8    ,25.2    ,25.6     /
       DATA (AE(I, 5,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,47.8    ,22.7    ,
      +  23.8    ,24.9    ,26.3     /
       DATA (AE(I, 5,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,35.5    ,
      +  23.9    ,25.9    ,26.6     /
       DATA (AE(I, 5,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,64.3    ,
      +  24.1    ,25.7    ,27.1     /
       DATA (AE(I, 5,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  34.0    ,25.7    ,27.7     /
       DATA (AE(I, 6, 1),I=1,10)  /
      +  16.6    ,16.5    ,16.8    ,16.7    ,17.0    ,16.5    ,16.7    ,
      +  18.3    ,18.9    ,19.0     /
       DATA (AE(I, 6, 2),I=1,10)  /
      +  16.2    ,16.6    ,17.2    ,17.4    ,17.9    ,17.4    ,17.7    ,
      +  20.7    ,22.0    ,22.6     /
       DATA (AE(I, 6, 3),I=1,10)  /
      +  18.9    ,18.7    ,18.8    ,18.6    ,18.9    ,18.6    ,18.9    ,
      +  21.0    ,22.3    ,22.9     /
       DATA (AE(I, 6, 4),I=1,10)  /
      +  18.3    ,12.7    ,14.2    ,15.0    ,15.7    ,16.1    ,16.3    ,
      +  16.5    ,17.9    ,19.0     /
       DATA (AE(I, 6, 5),I=1,10)  /
      +  .000E+00,15.7    ,15.1    ,15.3    ,16.5    ,16.4    ,16.4    ,
      +  17.0    ,18.3    ,19.4     /
       DATA (AE(I, 6, 6),I=1,10)  /
      +  .000E+00,22.9    ,14.9    ,15.2    ,16.2    ,16.9    ,17.4    ,
      +  18.2    ,19.5    ,21.1     /
       DATA (AE(I, 6, 7),I=1,10)  /
      +  .000E+00,40.7    ,18.4    ,15.9    ,17.1    ,17.7    ,18.9    ,
      +  19.5    ,20.3    ,21.1     /
       DATA (AE(I, 6, 8),I=1,10)  /
      +  .000E+00,.000E+00,23.3    ,16.2    ,16.3    ,17.3    ,18.7    ,
      +  19.5    ,20.3    ,21.1     /
       DATA (AE(I, 6, 9),I=1,10)  /
      +  .000E+00,.000E+00,49.2    ,19.0    ,19.1    ,19.4    ,20.2    ,
      +  20.8    ,21.6    ,22.0     /
       DATA (AE(I, 6,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,27.2    ,21.2    ,20.8    ,21.4    ,
      +  22.3    ,22.8    ,23.3     /
       DATA (AE(I, 6,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,45.6    ,25.0    ,22.8    ,23.9    ,
      +  23.6    ,24.3    ,24.4     /
       DATA (AE(I, 6,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,45.8    ,29.7    ,25.1    ,25.3    ,
      +  25.3    ,26.0    ,26.3     /
       DATA (AE(I, 6,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,42.7    ,29.0    ,28.0    ,
      +  27.0    ,27.2    ,27.6     /
       DATA (AE(I, 6,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,62.0    ,32.0    ,30.0    ,
      +  29.8    ,29.5    ,29.6     /
       DATA (AE(I, 6,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,44.5    ,34.4    ,
      +  32.7    ,31.5    ,31.8     /
       DATA (AE(I, 6,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,75.6    ,37.1    ,
      +  34.6    ,34.4    ,34.4     /
       DATA (AE(I, 6,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,51.2    ,45.2    ,
      +  39.0    ,37.5    ,36.4     /
       DATA (AE(I, 6,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,74.9    ,
      +  42.3    ,39.9    ,38.3     /
       DATA (AE(I, 6,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,69.5    ,
      +  50.7    ,42.3    ,41.4     /
       DATA (AE(I, 6,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  66.3    ,48.0    ,43.4     /
       DATA (AE(I, 7, 1),I=1,10)  /
      +  27.0    ,25.8    ,26.3    ,26.2    ,26.7    ,26.7    ,27.1    ,
      +  27.1    ,27.2    ,19.0     /
       DATA (AE(I, 7, 2),I=1,10)  /
      +  29.1    ,28.9    ,29.7    ,30.3    ,31.0    ,31.4    ,32.0    ,
      +  32.3    ,32.7    ,22.6     /
       DATA (AE(I, 7, 3),I=1,10)  /
      +  31.6    ,29.7    ,30.9    ,31.4    ,32.5    ,33.1    ,34.0    ,
      +  34.6    ,35.1    ,22.9     /
       DATA (AE(I, 7, 4),I=1,10)  /
      +  27.4    ,19.9    ,20.8    ,22.8    ,24.6    ,26.4    ,28.2    ,
      +  29.6    ,30.8    ,19.0     /
       DATA (AE(I, 7, 5),I=1,10)  /
      +  .000E+00,24.6    ,24.1    ,25.0    ,27.2    ,28.7    ,30.7    ,
      +  31.8    ,32.9    ,19.4     /
       DATA (AE(I, 7, 6),I=1,10)  /
      +  .000E+00,35.6    ,25.2    ,25.6    ,27.9    ,30.4    ,32.7    ,
      +  34.6    ,36.3    ,21.1     /
       DATA (AE(I, 7, 7),I=1,10)  /
      +  .000E+00,45.4    ,30.9    ,28.2    ,29.0    ,31.2    ,34.0    ,
      +  35.8    ,37.4    ,21.1     /
       DATA (AE(I, 7, 8),I=1,10)  /
      +  .000E+00,.000E+00,38.2    ,29.6    ,29.4    ,30.3    ,33.2    ,
      +  35.5    ,37.6    ,21.1     /
       DATA (AE(I, 7, 9),I=1,10)  /
      +  .000E+00,.000E+00,59.3    ,34.5    ,33.7    ,32.9    ,35.4    ,
      +  37.6    ,39.6    ,22.0     /
       DATA (AE(I, 7,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,44.5    ,37.8    ,37.5    ,37.2    ,
      +  39.0    ,41.4    ,23.3     /
       DATA (AE(I, 7,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,67.0    ,43.6    ,42.0    ,40.8    ,
      +  41.4    ,43.0    ,24.4     /
       DATA (AE(I, 7,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,49.9    ,50.9    ,44.6    ,43.9    ,
      +  44.2    ,44.2    ,26.3     /
       DATA (AE(I, 7,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,67.2    ,50.5    ,48.7    ,
      +  48.1    ,47.2    ,27.6     /
       DATA (AE(I, 7,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,68.1    ,55.2    ,52.3    ,
      +  51.5    ,51.6    ,29.6     /
       DATA (AE(I, 7,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,68.7    ,58.6    ,
      +  56.5    ,55.7    ,31.8     /
       DATA (AE(I, 7,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,89.3    ,62.9    ,
      +  60.0    ,59.1    ,34.4     /
       DATA (AE(I, 7,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,56.0    ,72.9    ,
      +  66.3    ,64.2    ,36.4     /
       DATA (AE(I, 7,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,105.    ,
      +  71.3    ,68.3    ,38.3     /
       DATA (AE(I, 7,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,73.4    ,
      +  76.8    ,72.4    ,41.4     /
       DATA (AE(I, 7,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  107.    ,79.9    ,43.4     /
       DATA (AE(I, 8, 1),I=1,10)  /
      +  35.5    ,35.3    ,35.7    ,35.7    ,36.3    ,36.3    ,36.7    ,
      +  36.7    ,36.7    ,19.0     /
       DATA (AE(I, 8, 2),I=1,10)  /
      +  40.6    ,41.4    ,41.9    ,42.3    ,43.2    ,43.5    ,44.0    ,
      +  44.3    ,44.5    ,22.6     /
       DATA (AE(I, 8, 3),I=1,10)  /
      +  45.4    ,45.7    ,46.4    ,47.0    ,48.1    ,48.7    ,49.4    ,
      +  49.8    ,50.2    ,22.9     /
       DATA (AE(I, 8, 4),I=1,10)  /
      +  43.9    ,44.3    ,43.4    ,45.1    ,47.3    ,48.7    ,49.6    ,
      +  50.5    ,51.3    ,19.0     /
       DATA (AE(I, 8, 5),I=1,10)  /
      +  .000E+00,49.3    ,49.6    ,50.5    ,53.2    ,54.2    ,55.4    ,
      +  56.1    ,56.8    ,19.4     /
       DATA (AE(I, 8, 6),I=1,10)  /
      +  .000E+00,59.1    ,53.0    ,55.4    ,58.0    ,60.0    ,61.2    ,
      +  62.5    ,63.6    ,21.1     /
       DATA (AE(I, 8, 7),I=1,10)  /
      +  .000E+00,54.5    ,57.1    ,59.2    ,62.3    ,64.4    ,66.0    ,
      +  67.3    ,68.5    ,21.1     /
       DATA (AE(I, 8, 8),I=1,10)  /
      +  .000E+00,.000E+00,65.9    ,62.1    ,65.1    ,67.6    ,69.4    ,
      +  71.1    ,72.6    ,21.1     /
       DATA (AE(I, 8, 9),I=1,10)  /
      +  .000E+00,.000E+00,72.2    ,67.1    ,70.5    ,73.1    ,75.1    ,
      +  76.8    ,78.4    ,22.0     /
       DATA (AE(I, 8,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,80.1    ,75.0    ,78.0    ,80.0    ,
      +  82.1    ,83.9    ,23.3     /
       DATA (AE(I, 8,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,94.5    ,82.2    ,82.8    ,85.1    ,
      +  87.3    ,89.2    ,24.4     /
       DATA (AE(I, 8,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,56.8    ,92.5    ,87.2    ,89.4    ,
      +  91.9    ,94.1    ,26.3     /
       DATA (AE(I, 8,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,116.    ,96.2    ,94.4    ,
      +  97.0    ,99.2    ,27.6     /
       DATA (AE(I, 8,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,78.1    ,104.    ,102.    ,
      +  102.    ,105.    ,29.6     /
       DATA (AE(I, 8,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,128.    ,111.    ,
      +  109.    ,110.    ,31.8     /
       DATA (AE(I, 8,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,104.    ,118.    ,
      +  117.    ,115.    ,34.4     /
       DATA (AE(I, 8,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,64.4    ,138.    ,
      +  124.    ,122.    ,36.4     /
       DATA (AE(I, 8,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,133.    ,
      +  133.    ,132.    ,38.3     /
       DATA (AE(I, 8,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,83.6    ,
      +  146.    ,139.    ,41.4     /
       DATA (AE(I, 8,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  166.    ,147.    ,43.4     /
       DATA (AE(I, 9, 1),I=1,10)  /
      +  43.3    ,43.2    ,43.6    ,43.8    ,44.1    ,44.3    ,44.7    ,
      +  44.8    ,44.8    ,19.0     /
       DATA (AE(I, 9, 2),I=1,10)  /
      +  50.9    ,51.4    ,52.0    ,52.6    ,53.1    ,53.6    ,54.2    ,
      +  54.5    ,54.7    ,22.6     /
       DATA (AE(I, 9, 3),I=1,10)  /
      +  58.0    ,58.4    ,59.3    ,60.1    ,60.7    ,61.5    ,62.3    ,
      +  62.7    ,63.1    ,22.9     /
       DATA (AE(I, 9, 4),I=1,10)  /
      +  62.0    ,63.9    ,63.7    ,65.7    ,65.5    ,67.5    ,68.2    ,
      +  68.9    ,69.7    ,19.0     /
       DATA (AE(I, 9, 5),I=1,10)  /
      +  .000E+00,72.2    ,72.5    ,74.2    ,74.2    ,76.1    ,77.0    ,
      +  77.8    ,78.6    ,19.4     /
       DATA (AE(I, 9, 6),I=1,10)  /
      +  .000E+00,80.4    ,80.5    ,83.1    ,83.0    ,85.5    ,86.8    ,
      +  88.1    ,89.2    ,21.1     /
       DATA (AE(I, 9, 7),I=1,10)  /
      +  .000E+00,63.4    ,88.5    ,91.3    ,91.1    ,94.0    ,95.8    ,
      +  97.3    ,98.6    ,21.1     /
       DATA (AE(I, 9, 8),I=1,10)  /
      +  .000E+00,.000E+00,98.8    ,98.6    ,97.8    ,102.    ,104.    ,
      +  106.    ,108.    ,21.1     /
       DATA (AE(I, 9, 9),I=1,10)  /
      +  .000E+00,.000E+00,84.1    ,107.    ,107.    ,111.    ,113.    ,
      +  116.    ,117.    ,22.0     /
       DATA (AE(I, 9,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,116.    ,115.    ,119.    ,122.    ,
      +  125.    ,127.    ,23.3     /
       DATA (AE(I, 9,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,111.    ,123.    ,127.    ,131.    ,
      +  134.    ,137.    ,24.4     /
       DATA (AE(I, 9,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,65.6    ,136.    ,135.    ,140.    ,
      +  143.    ,146.    ,26.3     /
       DATA (AE(I, 9,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,146.    ,144.    ,149.    ,
      +  152.    ,155.    ,27.6     /
       DATA (AE(I, 9,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,88.7    ,152.    ,158.    ,
      +  162.    ,165.    ,29.6     /
       DATA (AE(I, 9,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,181.    ,167.    ,
      +  171.    ,174.    ,31.8     /
       DATA (AE(I, 9,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,117.    ,174.    ,
      +  180.    ,183.    ,34.4     /
       DATA (AE(I, 9,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,72.0    ,201.    ,
      +  189.    ,192.    ,36.4     /
       DATA (AE(I, 9,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,151.    ,
      +  198.    ,201.    ,38.3     /
       DATA (AE(I, 9,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,95.2    ,
      +  220.    ,210.    ,41.4     /
       DATA (AE(I, 9,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  192.    ,217.    ,43.4     /
       DATA (AE(I,10, 1),I=1,10)  /
      +  62.1    ,62.1    ,62.6    ,62.9    ,63.3    ,63.3    ,64.0    ,
      +  64.0    ,64.0    ,19.0     /
       DATA (AE(I,10, 2),I=1,10)  /
      +  75.1    ,75.4    ,76.3    ,76.8    ,77.6    ,77.9    ,78.8    ,
      +  79.0    ,79.3    ,22.6     /
       DATA (AE(I,10, 3),I=1,10)  /
      +  87.5    ,88.3    ,89.4    ,90.2    ,91.3    ,91.9    ,93.0    ,
      +  93.5    ,93.9    ,22.9     /
       DATA (AE(I,10, 4),I=1,10)  /
      +  104.    ,104.    ,105.    ,106.    ,107.    ,108.    ,109.    ,
      +  110.    ,110.    ,19.0     /
       DATA (AE(I,10, 5),I=1,10)  /
      +  .000E+00,122.    ,122.    ,123.    ,124.    ,125.    ,126.    ,
      +  127.    ,128.    ,19.4     /
       DATA (AE(I,10, 6),I=1,10)  /
      +  .000E+00,138.    ,139.    ,140.    ,142.    ,143.    ,144.    ,
      +  146.    ,147.    ,21.1     /
       DATA (AE(I,10, 7),I=1,10)  /
      +  .000E+00,85.3    ,158.    ,159.    ,161.    ,162.    ,164.    ,
      +  166.    ,167.    ,21.1     /
       DATA (AE(I,10, 8),I=1,10)  /
      +  .000E+00,.000E+00,176.    ,177.    ,179.    ,181.    ,183.    ,
      +  184.    ,186.    ,21.1     /
       DATA (AE(I,10, 9),I=1,10)  /
      +  .000E+00,.000E+00,114.    ,199.    ,201.    ,202.    ,205.    ,
      +  206.    ,207.    ,22.0     /
       DATA (AE(I,10,10),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,218.    ,219.    ,220.    ,224.    ,
      +  225.    ,226.    ,23.3     /
       DATA (AE(I,10,11),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,150.    ,238.    ,238.    ,243.    ,
      +  244.    ,245.    ,24.4     /
       DATA (AE(I,10,12),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,85.8    ,255.    ,255.    ,261.    ,
      +  262.    ,263.    ,26.3     /
       DATA (AE(I,10,13),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,195.    ,272.    ,279.    ,
      +  279.    ,280.    ,27.6     /
       DATA (AE(I,10,14),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,115.    ,290.    ,296.    ,
      +  297.    ,298.    ,29.6     /
       DATA (AE(I,10,15),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,263.    ,313.    ,
      +  314.    ,315.    ,31.8     /
       DATA (AE(I,10,16),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,150.    ,330.    ,
      +  331.    ,332.    ,34.4     /
       DATA (AE(I,10,17),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,90.0    ,319.    ,
      +  349.    ,349.    ,36.4     /
       DATA (AE(I,10,18),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,196.    ,
      +  366.    ,367.    ,38.3     /
       DATA (AE(I,10,19),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,122.    ,
      +  387.    ,384.    ,41.4     /
       DATA (AE(I,10,20),I=1,10)  /
      +  .000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,.000E+00,
      +  247.    ,401.    ,43.4     /
       DATA (ERES(I, 1),I=1,10)  / 10*0./
       DATA (ERES(I, 2),I=1,10)  / 10*0./
       DATA (ERES(I, 3),I=1,10)  / 10*0./
       DATA (ERES(I, 4),I=1,10)  / 10*0./
       DATA (ERES(I, 5),I=1,10)  / 10*0./
       DATA (ERES(I, 6),I=1,10)  /
      +     0.000,   0.000,   0.000,   0.000,   0.000,   0.000,   0.000,
      +     2.780,   2.880,   2.890 /
       DATA (ERES(I, 7),I=1,10)  /
      +     1.500,   2.460,   2.510,   2.610,   2.700,   2.920,   3.070,
      +     3.200,   3.330,   2.890 /
       DATA (ERES(I, 8),I=1,10)  /
      +     4.470,   4.350,   4.390,   4.550,   4.660,   4.890,   4.980,
      +     5.100,   5.220,   2.890 /
       DATA (ERES(I, 9),I=1,10)  /
      +     7.480,   7.380,   7.370,   7.480,   7.510,   7.630,   7.660,
      +     7.750,   7.820,   2.890 /
       DATA (ERES(I,10),I=1,10)  /
      +    15.270,  15.190,  15.200,  15.370,  15.380,  15.430,  15.540,
      +    15.590,  15.630,   2.890 /
       END
 C->
       SUBROUTINE FRAGM (IAT,IAP, NW,B, NF, IAF)
 C...Nuclear Fragmentation, Abrasion-ablation model,
 C...Based on Jon Engel's routines ABRABL
 C...This most recent version adds for all prefragment
 C...masses > 10 the model calculation for the fragment
 C...mass distribution and the energy carried by the fragment
 C...of W. Friedmann
 C...The average values are used to implement the model
 C...in the montecarlo fashion / TSS, Dec '91
 C.
 C.  INPUT: IAP = mass of incident nucleus
 C.         IAT = mass of target   nucleus
 C.         NW = number of wounded nucleons in the beam nucleus
 C.         B  = impact parameter in the interaction
 C.
 C.  OUTPUT : NF = number of fragments  of the spectator nucleus
 C.           IAF(1:NF) = mass number of each fragment
 C.           PF(3,60) in common block /FRAGMENTS/ contains
 C.           the three momentum components (MeV/c) of each
 C.           fragment in the projectile frame
 C..............................................................
       SAVE
 
       COMMON /FRAGMENTS/ PPP(3,60)
       COMMON /FRAGMOD/A(10,10,20),AE(10,10,20),ERES(10,10),NFLAGG(10,10)
       DIMENSION IAF(60)
       DIMENSION AA(10), EAA(10)
       DATA AA/10.,15.,20.,25.,30.,35.,40.,45.,50.,56./
       DATA EAA/1.,2.,4.,6.,8.,10.,12.,16.,20.,30/
       AP=IAP
       AT=IAT
       NPF = IAP - NW
       IF (NPF .EQ. 0) THEN
          NF = 0
          RETURN
       ENDIF
 
       EB = ESTAR(AP,AT, B)
       EBP = ESTARP (NPF, NW)
 C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
       EB = EB + EBP
 C TOTAL E* IS THE SUM OF THE TWO COMPONENTS
 
 C.....Prefragment transverse momentum (MeV/nucleon)...
             FK = FERMK(AP)
 C FERMI MOMENTUM OF THE PROJECTILE NUCLEUS
             IF (NW .LT. IAP) THEN
             SIG = FK*SQRT(NW*NPF/(AP-1.))/3.162
 C GAUSSIAN SIGMA IN ALL THREE DIRECTION
             ELSE
             SIG = FK/3.162
 C THIS IS NOT CORRECT, TOO LARGE !!!!!!!!!!!!!!
             ENDIF
 ctp060203             PPFX = SIG*GASDEV(0)/NPF
 ctp060203             PPFY = SIG*GASDEV(0)/NPF
              PPFX = SIG*GASDEV(NW)/NPF
              PPFY = SIG*GASDEV(IAP)/NPF
 C THREE MOMENTUM COMPONENTS PER NUCLEON FOR THE PREFRAGMENT
 
 C.............Crude model for small prefragment mass .......
             IF (NPF .LT. 10) THEN
                  CALL EVAP(NPF, EB, EPS, NNUC, NALP)
 C                  EPS IS THE KINETIC ENERGY CARRIED BY THE EVAPORATED NUCLEONS
                ETOT = 938. + EPS
                  PP = SQRT((ETOT*ETOT - 8.79844E5)/3.)
 C                  AVERAGE MOMENTUM OF EVAPORATED NUCLEONS IN EACH DIRECTION
                  NUC = NPF - NNUC - 4*NALP
                  NF = 0
                  IF (NUC .GT. 0) THEN
                     NF = NF + 1
                     IAF(NF) = NUC
                     PPP(1,NF) = NUC*PPFX
                     PPP(2,NF) = NUC*PPFY
                  ENDIF
                  IF (NALP .NE. 0) THEN
                  DO I=1,NALP
                    NF = NF + 1
                     IAF(NF) = 4
                    CALL SINCO(S1,C1)
                    CALL SINCO(S2,C2)
                    PXE = 4.*PP*S1*S2
                    PYE = 4.*PP*S1*C2
                    PPP(1,NF) = 4.*PPFX + PXE
                    PPP(2,NF) = 4.*PPFY + PYE
                    PPP(1,1) = PPP(1,1) - PXE
                    PPP(2,1) = PPP(2,1) - PYE
                  ENDDO
                  ENDIF
                  IF (NNUC .NE. 0) THEN
                  DO I=1,NNUC
                     NF = NF + 1
                     IAF(NF) = 1
                     CALL SINCO(S1,C1)
                     CALL SINCO(S2,C2)
                     PXE = PP*S1*S2
                     PYE = PP*S1*C2
                     PPP(1,NF) = 4.*PPFX + PXE
                     PPP(2,NF) = 4.*PPFY + PYE
                     PPP(1,1) = PPP(1,1) - PXE
                     PPP(2,1) = PPP(2,1) - PYE
                  ENDDO
                  ENDIF
                  RETURN
             ENDIF
 
 C.........More refined model calculation .............
       JA = NPF/5 -1
       IF (JA .LT. 10) THEN
       IF ((NPF - AA(JA)) .GT. (AA(JA+1)-NPF)) JA = JA + 1
       ENDIF
       ARAT = FLOAT(NPF)/AA(JA)
       DO J=1,10
       IF (EB .LT. EAA(J)) GO TO 29
       ENDDO
       JE = 10
       GO TO 39
    29      JE = J
    39      IF (JE .GT. 1 .AND. JE .NE. 10) THEN
       IF ((EB - EAA(J-1)) .LT. (EAA(J)-EB)) JE = J - 1
       ENDIF
       ERAT = EB/EAA(JE)
         IF (EB .LT. 1.) THEN
         ERAT = EB
         ENDIF
 C INTERPOLATE BETWEEN EB=0. (NOTHING HAPPENS) AND EB = 1. MeV
 
          IF (JA .EQ. 10 .AND. JE .GT. 6) THEN
             WRITE(*,*)' JA=',JA,',   JE=',JE
          ENDIF
    43    ESUM = 0.
       NSUM = 0
       JF = 0
       DO J=20,1,-1
       FR =  A(JA, JE, J)*ARAT*ERAT
       N1 = 1 + FR
       FR1 = FR/FLOAT(N1)
       DO K=1, N1
       IF (S_RNDM(k) .LT. FR1) THEN
       JF = JF + 1
       IAF(JF) = J
       NSUM = NSUM + J
       EKIN = ERAT*AE(JA,JE, J)
          IF (EKIN .GT. 0.) THEN
          ESUM = ESUM + EKIN
          ETOT = 938.*IAF(JF) + EKIN
 ctp         PP = SQRT(2.*(ETOT*ETOT - IAF(JF)**2*8.79844E5)/3.)  !numerical precision problem
          PP = SQRT(2.*EKIN*(ETOT + IAF(JF)*938.)/3.)
          CALL SINCO(S1,C1)
          CALL SINCO(S2,C2)
          PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
          PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
          ENDIF
          IF (NSUM .GT. NPF) THEN
 C         WRITE(*,*)' WARNING, NSUM=', NSUM,',  NPF=',NPF
 C         WRITE(*,*)'  ARAT =', ARAT
          GO TO 43
         ELSE
         IF (NSUM .EQ. NPF) THEN
         GO TO 44
         ENDIF
         ENDIF
       ENDIF
       ENDDO
       ENDDO
       IF (NFLAGG(JA,JE) .EQ. 0) THEN
 C 'THE RESIDUE' IS A NUCLEAR FRAGMENT
       JF = JF + 1
       IAF(JF) = NPF - NSUM
       F1 = NPF*EB - ESUM
       IF (F1 .LT. 0.) F1 = 0.
 C GIVE THE REST OF EB TO THE FRAGMENT
       EKIN = F1
          IF (EKIN .GT. 0.) THEN
          ETOT = 938.*IAF(JF) + EKIN
 ctp         PP = SQRT(2.*(ETOT*ETOT - IAF(JF)**2*8.79844E5)/3.)  !numerical precision problem
          PP = SQRT(2.*EKIN*(ETOT + IAF(JF)*938.)/3.)
          CALL SINCO(S1,C1)
          CALL SINCO(S2,C2)
          PPP(1,JF) = PP*S1*S2 + IAF(JF)*PPFX
          PPP(2,JF) = PP*S1*C2 + IAF(JF)*PPFY
          ENDIF
       ELSE
 C 'THE RESIDUE' CONSISTS OF SPECTATOR NUCLEONS
       N1 = NPF - NSUM
       DO K=1,N1
       JF = JF + 1
       IAF(JF) = 1
       EKIN = ERAT*ERES(JA,JE)
          IF (EKIN .GT. 0.) THEN
          ETOT = 938.*IAF(JF) + EKIN
 ctp         PP = SQRT(2.*(ETOT*ETOT - IAF(JF)**2*8.79844E5)/3.)  !numerical precision problem
          PP = SQRT(2.*EKIN*(ETOT + IAF(JF)*938.)/3.)
          CALL SINCO(S1,C1)
          CALL SINCO(S2,C2)
          PPP(1,JF) = PP*S1*S2 + PPFX
          PPP(2,JF) = PP*S1*C2 + PPFY
          ENDIF
       ENDDO
       ENDIF
   44  NF = JF
       RETURN
       END
 C->
       FUNCTION ESTARP (NPF, NW)
 C CONTRIBUTION TO E* FROM ENERGY DEPOSITED BY SECONDARIES
 C VERY NAIVE VERSION INCORPORATING HUEFFNER'S IDEAS
       SAVE
 
       APF = NPF
       F1 = 15.3/APF**0.666666666
 C AVERAGE KINETIC ENERGY/NUCLEON IN PREFRAGMENT (MeV)
 C PER PATHLENGTH EQUAL TO THE PREFRAGMENT RADIUS
       ESTARP = 0.
       DO I=1,NW
       IF (S_RNDM(i) .GT. 0.5) THEN
 ctp060203      F2 = F1*RDIS(0)
       F2 = F1*RDIS()
       ESTARP = ESTARP + F2
       ENDIF
       ENDDO
 C SAMPLE RANDOMLY PER WOUNDED NUCLEON, x NW
       RETURN
       END
 
 ctp060203      function rdis(Idum)
       function rdis()
       SAVE
 
       dimension probr(20)
       data probr/
      *      0.10000, 0.15748, 0.21778, 0.28605, 0.36060,
      *      0.43815, 0.51892, 0.60631, 0.70002, 0.79325,
      *      0.88863, 0.98686, 1.10129, 1.21202, 1.32932,
      *      1.44890, 1.57048, 1.70139, 1.83417, 2.00000/
       nr = 20.*S_RNDM(0) + 1
       if (nr .eq. 1) then
       f1 = 0.
       else
       f1 = probr(nr-1)
       endif
       dr = probr(nr) - f1
       rdis = f1 + dr*S_RNDM(0)
       return
       end
 
 
       function estar(ap,at,b)
       implicit real*8(a-h,o-z)
       SAVE
 
       real*4 ap,at,b,estar
       sigma=4.5  !total n-n cross section in fm**2
       rt=.82*at**.3333 !target radius
       rp=.82*ap**.3333 !projectile radius
       alpha=rt**2/rp**2
       beta=b**2/rt**2
       f=at*sigma/(3.14159*rt**2)
       alf = log(f)
       alalf = log(alpha)
       gfac=0
       gfac1=0
       s1=0.
       s2=0.
       s3=0.
       ii=1
       do n=0,10 ! This limit may not need to be so high.
          if(n.ge.2) then
             gfac1=gfac
             gfac=gfac+log(float(n))
          endif
          g0=n*alf -n*beta*alpha/(n+alpha)+alalf
          g1=g0-log(alpha+n)-gfac
          g2=(n+2)*log(f)-(n+2)*beta*alpha/(n+2+alpha)
      >      +log(n+2+alpha+beta*alpha**2)-3*log(n+2+alpha)-gfac
          g3=g0-2*log(n+alpha)-gfac1
          ii=-ii
          s1=s1+ii*exp(g1)
          s2=s2+ii*exp(g2)
          if(n.ge.1) s3=s3+ii*exp(g3)
       enddo
 
       pb=s1
       e1b=197.**2/(2*938.*rp**2*pb) *s2
 c      a=b*(s3/pb-1)
 c      a=-b*s3/pb
 c      e2b=-.5* 938. * (41./(ap**.333))**2 * a**2 /(197.**2)
 c      estar=e1b+e2b
       estar = e1b
       return
       end
 
       subroutine evap(npf,eb,eps,nnuc,nalp)
       SAVE
 
       eps=7.5+sqrt(8*eb)
       n=min(npf*int(eb/eps),npf)
       nalp=n/5
       nnuc=n-4*nalp
       return
       end
 C->
       FUNCTION FERMK(A)
       SAVE
 
       DIMENSION AA(6), FK(6)
       DATA AA/4., 6., 12., 24., 40., 57./
       DATA FK/130.,169.,221.,235.,251.,260./
       DO I=2,4
       IF (A .LT. AA(I)) GO TO 25
       ENDDO
       I = 5
    25      F11 = AA(I-1)
       F12 = AA(I)
       F13 = AA(I+1)
       F21 = FK(I-1)
       F22 = FK(I)
       F23 = FK(I+1)
       FERMK = QUAD_INT(A,F11,F12,F13, F21,F22,F23)
       RETURN
       END
 
 C========================================================================
 C. Multiple interaction structure
 C========================================================================
 
       SUBROUTINE INT_NUC (IA, IB, SIG0, SIGEL)
 C...Compute with a montecarlo code  the  "multiple interaction structure"
 C.  of a nucleus-nucleus interaction
 C.
 C.  INPUT : IA            = mass of target nucleus
 C.          IB            = mass of projectile nucleus
 C.          SIG0 (mbarn)  = inelastic pp cross section
 C.          SIGEL(mbarn)  = elastic pp cross section
 C.
 C.  OUTPUT : in common block /CNUCMS/
 C.           B = impact parameter (fm)
 C.           BMAX = maximum impact parameter for generation
 C.           NTRY = number of "trials" before one interaction
 C.           NA = number of wounded nucleons in A
 C.           NB =    "        "        "     in B
 C.           NI = number of nucleon-nucleon inelastic interactions
 C.           NAEL = number of elastically scattered nucleons in  A
 C.           NBEL =    "         "           "          "    in  B
 C.           JJA(J)  [J=1:IA]   = number of inelastic interactions
 C.                                of J-th nucleon of nucleus A
 C.           JJB(J)  [J=1:IB]   = number of inelastic interactions
 C.                                of J-th nucleon of nucleus B
 C.           JJAEL(J)  [J=1:IA]   = number of elastic interactions
 C.                                of J-th nucleon of nucleus A
 C.           JJBEL(J)  [J=1:IB]   = number of elastic interactions
 C.                                of J-th nucleon of nucleus B
 C.           JJINT(J,K)  [J=1:NB, K=1:NA]  (0 = no interaction)
 C.                                         (1 = interaction )
 C.                                         between nucleon J of A and K of B
 C-----------------------------------------------------------------------------
       SAVE
 
       PARAMETER (IAMAX=56)
       COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
      +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
      +         ,JJAEL(IAMAX), JJBEL(IAMAX)
       DIMENSION XA(IAMAX), YA(IAMAX), XB(IAMAX), YB(IAMAX)
       DATA PI /3.1415926/
       SIGT = SIG0 + SIGEL
       R2  = 0.1 * SIG0/PI
       R2T = 0.1 * SIGT/PI
       BMAX = 15.                             ! fm
       NTRY = 0
       CALL NUC_CONF (IA, XA, YA)
       CALL NUC_CONF (IB, XB, YB)
       NI = 0
       NIEL = 0
       DO JA=1,IA
          JJA(JA) = 0
          JJAEL(JA) = 0
       ENDDO
       DO JB=1,IB
          JJB(JB) = 0
          JJBEL(JB) = 0
          DO JA=1,IA
             JJINT(JB,JA) = 0
          ENDDO
       ENDDO
 1000  B = BMAX*SQRT(S_RNDM(0))
       PHI = 2.*PI*S_RNDM(0)
       BX = B*COS(PHI)
       BY = B*SIN(PHI)
       NTRY = NTRY+1
       DO JA=1,IA
          DO JB=1,IB
             S = (XA(JA)-XB(JB)-BX)**2 + (YA(JA)-YB(JB)-BY)**2
             IF (S .LT. R2)  THEN
                NI = NI + 1
                JJA(JA) = JJA(JA)+1
                JJB(JB) = JJB(JB)+1
                JJINT(JB,JA) = 1
             ELSE IF (S .LT. R2T)  THEN
                NIEL = NIEL + 1
                JJAEL(JA) = JJAEL(JA)+1
                JJBEL(JB) = JJBEL(JB)+1
             ENDIF
          ENDDO
       ENDDO
       IF (NI + NIEL .EQ. 0)  GOTO 1000
       NA = 0
       NB = 0
       NAEL = 0
       NBEL = 0
       DO JA=1,IA
          IF (JJA(JA) .GT. 0)  THEN
             NA = NA + 1
          ELSE
             IF (JJAEL(JA) .GT. 0)  NAEL = NAEL+1
          ENDIF
       ENDDO
       DO JB=1,IB
          IF (JJB(JB) .GT. 0)  THEN
             NB = NB + 1
          ELSE
             IF (JJBEL(JB) .GT. 0)  NBEL = NBEL+1
          ENDIF
       ENDDO
       RETURN
       END
 
        SUBROUTINE NUC_CONF (IA, XX, YY)
 C...This routine generates the configuration  of a nucleus
 C.  need an initialization call to NUC_GEOM_INI
 C.
 C.  INPUT  : IA = mass number of the nucleus
 C.  OUTPUT : XX(1:IA), YY(1:IA) (fm) = position in impact parameter
 C.                                     space of the IA nucleons
 C...................................................................
       SAVE
 
       PARAMETER (IAMAX=56)
       DIMENSION XX(IAMAX), YY(IAMAX)
       PARAMETER (NB=401)
       COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
       DATA PI /3.1415926/
       DO J=1,IA
          Z = S_RNDM(j)
          JZ = INT((Z-ZMIN)/DZ)+1
 CDH
          JZ = MIN(JZ,400)
          T = (Z-ZMIN)/DZ - FLOAT(JZ-1)
          B = BBZ(JZ,IA)*(1.-T) + BBZ(JZ+1,IA)*T
          PHI = 2.*PI*S_RNDM(jz)
          XX(J) = B*COS(PHI)
          YY(J) = B*SIN(PHI)
       ENDDO
       RETURN
       END
 
       SUBROUTINE NUC_GEOM_INI
 C...Initialize all nucleus profiles
       SAVE
 
       PARAMETER (NB=401)
       PARAMETER (IAMAX=56)
       COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
       COMMON /CPROFA/ ZMIN, DZ, BBZ(NB,IAMAX)
       DIMENSION FFB(NB), GGB(NB)
       DATA PI /3.1415926/
       CALL SHELL_INI
       CALL WOOD_SAXON_INI
       DO IA= 2,IAMAX
            JA = IA
          CALL NUC_PROFIL(JA)
          DO K=1,NB
            FFB(K) = BB(K)*TB(K) * (2.*PI)
          ENDDO
          GGB(1) = 0.
          GGB(NB) = 1.
          DO K=2,NB-1
            GGB(K) = GGB(K-1) + FFB(K-1)*DB
          ENDDO
          CALL INVERT_ARRAY(GGB,0.,DB,NB, BBZ(1,IA), ZMIN, DZ)
       ENDDO
       RETURN
       END
 
       SUBROUTINE NUC_PROFIL (JA)
 C...Compute the profile function T(b)
 C.  normalised as INT[d2b T(b) = 1]
 C.  INPUT : JA = integer mass number of nucleus
 C...............................................
       SAVE
 
       PARAMETER (NB=401)
       EXTERNAL DENSA
       REAL DENSA
       COMMON /CC01/  B
       COMMON /CCDA/ JJA
       COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
       BMAX = 7.5
       DB = BMAX/FLOAT(NB-1)
       JJA = JA
       A = JA
       DO JB=1,NB
         B = DB*FLOAT(JB-1)
         BB(JB) = B
         IF (JA .LE. 18)  THEN
             TB(JB) = PROFNUC (B, JA)
          ELSE
             TB(JB) = 2.*GAUSS (DENSA,0.,BMAX)
          ENDIF
       ENDDO
       RETURN
       END
 
       SUBROUTINE NUC1_PROFIL (AA)
 C...Compute the profile function T(b)
 C.  normalised as INT[d2b T(b) = 1]
 C.  INPUT : AA = mass number of nucleus
 C...............................................
       SAVE
 
       PARAMETER (NB=401)
       EXTERNAL DENSA
       REAL DENSA
       COMMON /CC01/  B
       COMMON /CPROF/ DB, BMAX, BB(NB), TB(NB), A
       A = AA
       IA1 = INT(AA)
       IA2 = IA1 + 1
       U = AA - FLOAT(IA1)
       BMAX = 7.5
       DB = BMAX/FLOAT(NB-1)
       DO JB=1,NB
          B = DB*FLOAT(JB-1)
          BB(JB) = B
          IF (A .LE. 18.)  THEN
              T1 = PROFNUC (B, IA1)
              T2 = PROFNUC (B, IA2)
           ELSE
              JJA = IA1
              T1 = 2.*GAUSS (DENSA,0.,BMAX)
              JJA = IA2
              T2 = 2.*GAUSS (DENSA,0.,BMAX)
           ENDIF
           TB(JB) = (1.-U)*T1  + U*T2
       ENDDO
       RETURN
       END
 
 C===========================================================================
 C.   Code about nuclear densities
 C===========================================================================
 
       FUNCTION DENS_NUC (R, JA)
 C....Nuclear density (normalised to 1)
 C.   for a nucleus of mass number JA
 C.   INPUT R = radial coordinate  (fm)
 C.         JA = integer mass number
 C.  OUTPUT (fm**-3)
 C--------------------------------------------------------
       SAVE
 
       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
       IF (JA .GT. 18)  THEN
          DENS_NUC = WOOD_SAXON(R,JA)
       ELSE IF (JA .NE. 4)  THEN
          DENS_NUC = HELIUM(R)
       ELSE
          DENS_NUC = SHELL(R,JA)
       ENDIF
       RETURN
       END
 
       FUNCTION WOOD_SAXON (R, JA)
 C....Wood-Saxon nuclear density (normalised to 1)
 C.   for a nucleus of mass number A.
 C.   INPUT R =  (fm)
 C.         JA = mass number
 C.   OUTPUT (fm**-3)
 C------------------------------------------------------
       SAVE
 
       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
       WOOD_SAXON = CC0(JA)/(1.+EXP((R-RR0(JA))/AA0(JA)))
       RETURN
       END
 
       FUNCTION HELIUM (R)
 C... Helium density from Barrett and Jackson
 C.   INPUT R = r coordinate (fm)
 C.   OUTPUT (fm**-3)
 C........................................................
       SAVE
 
       DATA R0 /0.964/, CA /0.322/   ! fm
       DATA W /0.517/, CC /5.993224E-02/
       HELIUM = CC*(1.+W*(R/R0)**2)/(1. + EXP((R-R0)/CA))
       RETURN
       END
 
       FUNCTION SHELL (R,JA)
 C...Density in the shell model
       COMMON /CSHELL/ RR0(18), RR02(18)
       SAVE
 
       DATA PI /3.1415926/
       R0 = RR0(JA)
       C1 = MIN(1.,4./FLOAT(JA))
       CS = 1./(R0**3*PI**(1.5))
       CP = 2.*CS/3.
       FS = EXP(-(R/R0)**2)
       FP = (R/R0)**2 * FS
       SHELL = C1*CS*FS + (1.-C1)*CP*FP
       RETURN
       END
 
       FUNCTION PROFNUC (B, JA)
 C...This function return
 C.  the profile T(b) for a nucleus of mass number A
 C.  INPUT B = impact parameter (GeV**-1)
 C.        JA = integer mass number
 C.  OUTPUT  (fm**-2)
 C.
 C.  The  density of the nucleus is the `shell model density'
 C.  the parameter r0 must beinitialized in the common block
 C.............................................................
       SAVE
 
       COMMON /CSHELL/ RR0(18), RR02(18)
       DATA PI /3.1415926/
       B2 = B*B
       ARG = B2/RR02(JA)
       TS = EXP(-ARG)
       TP = TS*(2.*B2+RR02(JA))/(3.*RR02(JA))
       CS = MIN(1.,4./FLOAT(JA))
       PROFNUC = (CS*TS + (1.-CS)*TP)/(PI*RR02(JA))
       RETURN
       END
 
       SUBROUTINE SHELL_INI
 C...Initialize the parameter  of the shell model
 C.  for the nuclei with    6 < A < 18
 C..............................................
       SAVE
 
       COMMON /CSHELL/ RR0(18), RR02(18)
       DIMENSION RR(18)
 C...Data on Sqrt[<r**2>]  in fermi
       DATA RR /0.81,2.095,1.88,1.674, -1.,2.56,2.41,-1.,2.519,2.45
      +          ,2.37, 2.460, 2.440, 2.54, 2.58, 2.718, 2.662,2.789 /
       DO JA=1,18
          A = FLOAT(JA)
          RMED = RR(JA)
          IF (RMED .LE. 0.)   RMED = 0.5*(RR(JA-1) + RR(JA+1))
          C = MAX(1.5,(5./2. - 4./A) )
          R0 = RMED/SQRT(C)
          RR0 (JA) = R0
          RR02(JA) = R0*R0
       ENDDO
       RETURN
       END
 C->
       SUBROUTINE WOOD_SAXON_INI
       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
       SAVE
 
       DATA PI /3.1415926/
 C...Wood-Saxon parameters from  table 6.2   of Barrett and Jackson
       RR0 (19) = 2.59
       AA0 (19) = 0.564
       RR0 (20) = 2.74
       AA0 (20) = 0.569
       RR0 (22) = 2.782
       AA0 (22) = 0.549
       RR0 (24) = 2.99
       AA0 (24) = 0.548
       RR0 (27) = 2.84
       AA0 (27) = 0.569
       RR0 (28) = 3.14
       AA0 (28) = 0.537
       RR0 (29) = 3.77
       AA0 (29) = 0.52
       RR0 (48) = 3.912
       AA0 (48) = 0.5234
       RR0 (56) = 3.98
       AA0 (56) = 0.569
       DO J=19, 56
          IF (RR0(J) .LE. 0.)  THEN
             RR0(J) = 1.05*FLOAT(J)**0.333333
             AA0(J) = 0.545
          ENDIF
          CC0(J)=3./(4.*PI*RR0(J)**3)/(1.+((AA0(J)*PI)/RR0(J))**2)
       ENDDO
       RETURN
       END
 
       FUNCTION DENSA (Z)
 C....Woods Saxon nuclear density (normalised to 1)
 C.   for a nucleus of mass number A.
 C.   INPUT z = z coordinate (fm)
 C.         JA = integer mass number
 C.         B (in common /CC01/)  impact parameter  (fm)
 C.  OUTPUT (fm**-3)
 C--------------------------------------------------------
       SAVE
 
       COMMON /CC01/  B
       COMMON /CCDA/ JA
       COMMON /CWOOD/ RR0(19:56), AA0(19:56), CC0(19:56)
       R = SQRT (Z*Z + B*B)
       DENSA = CC0(JA)/(1.+EXP((R-RR0(JA))/AA0(JA)))
       RETURN
       END
 
 C==========================================================================
 C. Cross sections
 C==========================================================================
 
       SUBROUTINE SIGMA_AIR (IB,SIG0,SIGEL,NINT,
      +                            SIGMA,DSIGMA,SIGQE,DSIGQE)
 C...Compute with a montecarlo method the "production"
 C.  and "quasi-elastic" cross section for
 C.  a nucleus-air  interaction
 C.
 C.  INPUT : IB            = mass of projectile nucleus
 C.          SIG0 (mbarn)  = inelastic pp cross section
 C.          NINT            = number  of interactions to generate
 C.  OUTPUT : SIGMA (mbarn) = "production" cross section
 C.           DSIGMA   "    = error
 C.           SIGQE    "    = "quasi-elastic" cross section
 C.           DSIGQE   "    = error
 C.           additional output is in the common block  /CPROBAB/
 C..........................................................................
       SAVE
 
       PARAMETER (IAMAX=56)
       PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
       COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX),
      +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
      +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
      +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
       COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
      +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
      +         ,JJAEL(IAMAX), JJBEL(IAMAX)
       DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
       DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
       DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
       DATA WOX /0.346/
       DATA PI /3.1415926/
       R2 = 0.1 * SIG0/PI
       BMAX = 15.                             ! fm
       SIGMA0 = PI*BMAX*BMAX*10.              ! mbarn
       IA = 16
       DO J=1,IA
          MMA(J) = 0
          M1AEL(J) = 0
          M2AEL(J) = 0
       ENDDO
       DO J=1,IB
          MMB(J) = 0
          M1BEL(J) = 0
          M2BEL(J) = 0
       ENDDO
       DO J=1,IA*IB
          MMI(J) = 0
       ENDDO
       NN = 0
       M = 0
       DO KK=1,NINT
          IA = 14 + 2*INT((1.+WOX)*S_RNDM(kk))
          CALL INT_NUC (IA, IB, SIG0, SIGEL)
          NN = NN + NTRY
          MMI(NI) = MMI(NI) + 1
          MMA(NA) = MMA(NA)+1
          MMB(NB) = MMB(NB)+1
          IF (NI .GT. 0)  THEN
             M = M+1
             M1AEL(NAEL) = M1AEL(NAEL)+1
             M1BEL(NBEL) = M1BEL(NBEL)+1
          ELSE
             M2AEL(NAEL) = M2AEL(NAEL)+1
             M2BEL(NBEL) = M2BEL(NBEL)+1
          ENDIF
       ENDDO
       MQE = NINT - M
       SIGMA  = SIGMA0 * FLOAT(M)/FLOAT(NN)
       DSIGMA = SIGMA0 * SQRT(FLOAT(M))/FLOAT(NN)
       SIGQE  = SIGMA0 * FLOAT(MQE)/FLOAT(NN)
       DSIGQE = SIGMA0 * SQRT(FLOAT(MQE))/FLOAT(NN)
       DO J=1,IA
          PROBA(J) = FLOAT(MMA(J))/FLOAT(M)
          DPROBA(J) = SQRT(FLOAT(MMA(J)))/FLOAT(M)
       ENDDO
       DO J=1,IB
          PROBB(J) = FLOAT(MMB(J))/FLOAT(M)
          DPROBB(J) = SQRT(FLOAT(MMB(J)))/FLOAT(M)
       ENDDO
       DO J=1,IA*IB
          PROBI(J) = FLOAT(MMI(J))/FLOAT(M)
          DPROBI(J) = SQRT(FLOAT(MMI(J)))/FLOAT(M)
       ENDDO
       DO J=0,IA
          P1AEL(J) = FLOAT(M1AEL(J))/FLOAT(M)
          DP1AEL(J) = SQRT(FLOAT(M1AEL(J)))/FLOAT(M)
          P2AEL(J) = FLOAT(M2AEL(J))/FLOAT(MQE)
          DP2AEL(J) = SQRT(FLOAT(M2AEL(J)))/FLOAT(MQE)
       ENDDO
       DO J=0,IB
          P1BEL(J) = FLOAT(M1BEL(J))/FLOAT(M)
          DP1BEL(J) = SQRT(FLOAT(M1BEL(J)))/FLOAT(M)
          P2BEL(J) = FLOAT(M2BEL(J))/FLOAT(MQE)
          DP2BEL(J) = SQRT(FLOAT(M2BEL(J)))/FLOAT(MQE)
       ENDDO
       RETURN
       END
 C->
       SUBROUTINE SIGMA_MC (IA,IB,SIG0,SIGEL,NINT,
      +                            SIGMA,DSIGMA,SIGQE,DSIGQE)
 C...Compute with a montecarlo method the "production"
 C.  and "quasi-elastic" cross section for
 C.  a nucleus-nucleus interaction
 C.
 C.  INPUT : IA            = mass of target nucleus
 C.          IB            = mass of projectile nucleus
 C.          SIG0 (mbarn)  = inelastic pp cross section
 C.          NINT            = number  of interactions to generate
 C.  OUTPUT : SIGMA (mbarn) = "production" cross section
 C.           DSIGMA   "    = error
 C.           SIGQE    "    = "quasi-elastic" cross section
 C.           DSIGQE   "    = error
 C.           additional output is in the common block  /CPROBAB/
 C.           Prob(n_A), Prob(n_B), Prob(n_int)
 C..........................................................................
       SAVE
 
       PARAMETER (IAMAX=56)
       PARAMETER (IAMAX2=3136)          ! IAMAX*IAMAX
       COMMON  /CPROBAB/ PROBA(IAMAX), DPROBA(IAMAX),
      +   PROBB(IAMAX), DPROBB(IAMAX), PROBI(IAMAX2), DPROBI(IAMAX2),
      +   P1AEL(0:IAMAX),DP1AEL(0:IAMAX),P1BEL(0:IAMAX), DP1BEL(0:IAMAX),
      +   P2AEL(0:IAMAX),DP2AEL(0:IAMAX),P2BEL(0:IAMAX), DP2BEL(0:IAMAX)
       COMMON /CNUCMS/ B, BMAX, NTRY, NA, NB, NI, NAEL, NBEL
      +         ,JJA(IAMAX), JJB(IAMAX), JJINT(IAMAX,IAMAX)
      +         ,JJAEL(IAMAX), JJBEL(IAMAX)
       DIMENSION  MMA(0:IAMAX), MMB(0:IAMAX), MMI(0:IAMAX2)
       DIMENSION  M1AEL(0:IAMAX), M1BEL(0:IAMAX)
       DIMENSION  M2AEL(0:IAMAX), M2BEL(0:IAMAX)
       DATA PI /3.1415926/
       R2 = 0.1 * SIG0/PI
       BMAX = 15.                             ! fm
       SIGMA0 = PI*BMAX*BMAX*10.              ! mbarn
       DO J=1,IA
          MMA(J) = 0
          M1AEL(J) = 0
          M2AEL(J) = 0
       ENDDO
       DO J=1,IB
          MMB(J) = 0
          M1BEL(J) = 0
          M2BEL(J) = 0
       ENDDO
       DO J=1,IA*IB
          MMI(J) = 0
       ENDDO
       NN = 0
       M = 0
       DO KK=1,NINT
          CALL INT_NUC (IA, IB, SIG0, SIGEL)
          NN = NN + NTRY
          MMI(NI) = MMI(NI) + 1
          MMA(NA) = MMA(NA)+1
          MMB(NB) = MMB(NB)+1
          IF (NI .GT. 0)  THEN
             M = M+1
             M1AEL(NAEL) = M1AEL(NAEL)+1
             M1BEL(NBEL) = M1BEL(NBEL)+1
          ELSE
             M2AEL(NAEL) = M2AEL(NAEL)+1
             M2BEL(NBEL) = M2BEL(NBEL)+1
          ENDIF
       ENDDO
       MQE = NINT - M
       SIGMA  = SIGMA0 * FLOAT(M)/FLOAT(NN)
       DSIGMA = SIGMA0 * SQRT(FLOAT(M))/FLOAT(NN)
       SIGQE  = SIGMA0 * FLOAT(MQE)/FLOAT(NN)
       DSIGQE = SIGMA0 * SQRT(FLOAT(MQE))/FLOAT(NN)
       DO J=1,IA
          PROBA(J) = FLOAT(MMA(J))/FLOAT(M)
          DPROBA(J) = SQRT(FLOAT(MMA(J)))/FLOAT(M)
       ENDDO
       DO J=1,IB
          PROBB(J) = FLOAT(MMB(J))/FLOAT(M)
          DPROBB(J) = SQRT(FLOAT(MMB(J)))/FLOAT(M)
       ENDDO
       DO J=1,IA*IB
          PROBI(J) = FLOAT(MMI(J))/FLOAT(M)
          DPROBI(J) = SQRT(FLOAT(MMI(J)))/FLOAT(M)
       ENDDO
       DO J=0,IA
          P1AEL(J) = FLOAT(M1AEL(J))/FLOAT(M)
          DP1AEL(J) = SQRT(FLOAT(M1AEL(J)))/FLOAT(M)
          P2AEL(J) = FLOAT(M2AEL(J))/FLOAT(MQE)
          DP2AEL(J) = SQRT(FLOAT(M2AEL(J)))/FLOAT(MQE)
       ENDDO
       DO J=0,IB
          P1BEL(J) = FLOAT(M1BEL(J))/FLOAT(M)
          DP1BEL(J) = SQRT(FLOAT(M1BEL(J)))/FLOAT(M)
          P2BEL(J) = FLOAT(M2BEL(J))/FLOAT(MQE)
          DP2BEL(J) = SQRT(FLOAT(M2BEL(J)))/FLOAT(MQE)
       ENDDO
       RETURN
       END
 
 C=============================================================
 C.  Cross sections
 C=============================================================
 
 
       SUBROUTINE SIG_H_AIR (SSIG, SLOPE, ALPHA,  SIGT, SIGEL, SIGQE)
 C...Subroutine to compute hadron-air cross sections
 C.  according to:
 C.  R.J. Glauber and G.Matthiae  Nucl.Phys. B21, 135, (1970)
 C.
 C.  Air is a linear combination of Nitrogen and oxygen
 C.
 C.  INPUT :  SSIG  (mbarn) total pp cross section
 C.           SLOPE (GeV**-2)  elastic scattering slope for pp
 C.           ALPHA    real/imaginary part of the forward pp elastic
 C.                                               scattering amplitude
 C.  OUTPUT : SIGT  = Total cross section
 C.           SIGEL = Elastic cross section
 C.           SIGQEL  = Elastic + Quasi elastic cross section
 C......................................................................
       SAVE
 
       DATA  FOX /0.257/
       CALL GLAUBER(14,SSIG,SLOPE,ALPHA,SIG1,SIGEL1,SIGQE1)
       CALL GLAUBER(16,SSIG,SLOPE,ALPHA,SIG2,SIGEL2,SIGQE2)
       SIGT  = (1.-FOX)*SIG1   + FOX*SIG2
       SIGEL = (1.-FOX)*SIGEL1 + FOX*SIGEL2
       SIGQE = (1.-FOX)*SIGQE1 + FOX*SIGQE2
       RETURN
       END
 
       SUBROUTINE GLAUBER(JA,SSIG,SLOPE,ALPHA,SIGT,SIGEL,SIGQEL)
 C...Subroutine to compute hadron-Nucleus cross sections
 C.  according to:
 C.  R.J. Glauber and G.Matthiae  Nucl.Phys. B21, 135, (1970)
 C.
 C.  This formulas assume that the target nucleus  density is
 C.  modeled by a shell-model form.  A reasonable range of models
 C   is  4 < JA < 18
 C.
 C.  INPUT :  A = mass number of the nucleus
 C.           SSIG  (mbarn) total pp cross section
 C.           SLOPE (GeV**-2)  elastic scattering slope for pp
 C.           ALPHA    real/imaginary part of the forward pp elastic
 C.                                               scattering amplitude
 C.  OUTPUT : SIGT  = Total cross section
 C.           SIGEL = Elastic cross section
 C.           SIGQEL  = Elastic + Quasi elastic cross section
 C.
 C. Internally  everything is computed in GeV (length = GeV**-1)
 C......................................................................
       SAVE
 
       COMMON /CA0SH/ R0, R02
       COMPLEX  ZZ, ZS, ZP, ZC
       DIMENSION RR(18)
       DATA CMBARN /0.389385/
       DATA PI /3.1415926/
       DATA BMAX /50./            ! GeV**-1
       DATA NB /100/
 C...data on Sqrt[<r**2>] (fm). (A=5,8 are not correct).
 C   from Barett and Jackson
       DATA RR /0.81,2.095,1.88,1.674, 2.56,2.56,2.41,2.5,2.519,2.45
      +          ,2.37, 2.460, 2.440, 2.54, 2.58, 2.718, 2.662,2.789 /
       A = FLOAT(JA)
 C...Parameter of shell model density
       R0 = RR(JA)/0.197/SQRT(5./2. - 4./A)         ! GeV**-1
       R02 = R0*R0
       SIG = SSIG/CMBARN                           ! GeV**-2
       DB = BMAX/FLOAT(NB)
       SUM = 0.
       SUM1 = 0.
       SUM2 = 0.
       DO JB=1,NB
          B = DB*(FLOAT(JB)-0.5)
          GS = GLAUBGS (B,SLOPE, SIG)
          GP = GLAUBGP (B,SLOPE, SIG)
          XS = (1.- GS)
          YS =  GS*ALPHA
          ZS = CMPLX(XS,YS)
          XP = (1.- GP)
          YP =  GP*ALPHA
          ZP = CMPLX(XP,YP)
          ZZ = ZS**4. * ZP**(A-4.)
          X = REAL (ZZ)
          Y = AIMAG(ZZ)
          ZC = CMPLX(X,-Y)
          SUM = SUM + (1.-X)*B
          SUM1 = SUM1 + ((1.-X)**2 + Y**2)*B
          OMS = OMEGAS(B,SIG,SLOPE,ALPHA)
          OMP = OMEGAP(B,SIG,SLOPE,ALPHA)
          OM = (1.- 2.*GS + OMS)**4. * (1. -2.*GP + OMP)**(A-4.)
          SUM2 = SUM2 + (1.-2.*X + OM)*B
       ENDDO
       SIGT =   SUM  * DB * 4.*PI * CMBARN
       SIGEL =  SUM1 * DB * 2.*PI * CMBARN
       SIGQEL = SUM2 * DB * 2.*PI * CMBARN
       RETURN
       END
 C->
       FUNCTION GLAUBGS (B,SLOPE, SIG)
       SAVE
 
       COMMON /CA0SH/ A0, A02
       DATA PI /3.1415926/
       GAMMA2 = A02/4. + 0.5*SLOPE
       ARG = B**2/(4.*GAMMA2)
       GLAUBGS = SIG/(8.*PI*GAMMA2) * EXP(-ARG)
       RETURN
       END
 C->
       FUNCTION GLAUBGP (B,SLOPE, SIG)
       SAVE
 
       COMMON /CA0SH/ A0, A02
       DATA PI /3.1415926/
       GAMMA2 = A02/4. + 0.5*SLOPE
       ARG = B**2/(4.*GAMMA2)
       C1 = 1.- A02/(6.*GAMMA2)*(1.-ARG)
       GLAUBGP = SIG/(8.*PI*GAMMA2) *  C1 * EXP(-ARG)
       RETURN
       END
 C->
       FUNCTION OMEGAS (B, SIG, SLOPE, RHO)
       SAVE
 
       COMMON /CA0SH/ A0, A02
       DATA PI /3.1415926/
       ETA2 = 0.25*(A02 + SLOPE)
       F02 = SIG*SIG*(1.+RHO*RHO)/(16.*PI**2)
       ARG = -B*B/(4.*ETA2)
       OMEGAS = F02/(4.*ETA2*SLOPE) *EXP(ARG)
       RETURN
       END
 C->
       FUNCTION OMEGAP (B, SIG, SLOPE, RHO)
       SAVE
 
       COMMON /CA0SH/ A0, A02
       DATA PI /3.1415926/
       ETA2 = 0.25*(A02 + SLOPE)
       F02 = SIG*SIG*(1.+RHO*RHO)/(16.*PI**2)
       ARG = -B*B/(4.*ETA2)
       OMEGAP=F02/(4.*ETA2*SLOPE)*(1.-A02/(6.*ETA2)*(1.+ARG))*EXP(ARG)
       RETURN
       END
 
 C------------------------------------------------------------------------
 C.  Fit of Block and Cahn to pp and pbar-p cross sections
 C------------------------------------------------------------------------
 
       SUBROUTINE BLOCK(SQS,SIG1,SIG2,SLOP1,SLOP2,
      +                 RHO1,RHO2,SIGEL1,SIGEL2)
 C...p-p and pbar-p cross sections
 C.  Parametrization of  Block and Cahn
 C
 C.  INPUT  : SQS   (GeV)  = c.m. energy
 C.
 C.  OUPUT : SIG1 (mbarn)    = pp  total  cross section
 C.          SLOP1 (GeV**2)  = slope of elastic scattering
 C.          RHO1            = Real/Imaginary part of the amplitude
 C.                            for forward elastic  scattering (pp)
 C.          SIGEL1 (mbarn)  = pp  elastic scattering  cross section
 C.          [1 -> 2   : pp -> pbar p]
 C-----------------------------------------------------------------------
       SAVE
 
       DATA PI /3.1415926/
       DATA CMBARN /0.389385/
       S = SQS*SQS
       CALL FPLUS  (S, FR, FI)
       CALL FMINUS (S, GR, GI)
       SIG1 = FI-GI
       SIG2 = FI+GI
       RHO1 = (FR-GR)/(FI-GI)
       RHO2 = (FR+GR)/(FI+GI)
       CALL SSLOPE (S, BP, BM)
       SLOP1 = BP - GI/FI*(BM-BP)
       SLOP2 = BP + GI/FI*(BM-BP)
       SIGEL1 = SIG1**2*(1.+RHO1**2)/(16.*PI*SLOP1)/CMBARN
       SIGEL2 = SIG2**2*(1.+RHO2**2)/(16.*PI*SLOP2)/CMBARN
       RETURN
       END
 
       SUBROUTINE FPLUS (S, FR, FI)
       SAVE
 
       COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
       COMPLEX Z1, Z2, Z3
       DATA PI /3.1415926/
       F1 = LOG(S/S0)
       Z1 = CMPLX(F1,-PI/2.)
       Z1 = Z1*Z1
       Z2 = 1. + A0*Z1
       Z3 = Z1/Z2
       F2 = CC*S**(AMU-1.)
       F3 = 0.5*PI*(1.-AMU)
       FI = AA + F2*COS(F3) + BETA*REAL(Z3)
       FR = -BETA*AIMAG(Z3)+F2*SIN(F3)
       RETURN
       END
 
       SUBROUTINE FMINUS (S, FR, FI)
       SAVE
 
       COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
       DATA PI /3.1415926/
       F1 = S**(ALPHA-1.)
       F2 = 0.5*PI*(1.-ALPHA)
       FR = -DD*F1*COS(F2)
       FI = -DD*F1*SIN(F2)
       RETURN
       END
 
       SUBROUTINE SSLOPE (S, BP, BM)
       SAVE
 
       COMMON /BLOCKD/ CP, DP, EP, CM, DM
       AL = LOG(S)
       BP = CP + DP*AL + EP*AL*AL
       BM = CM + DM*AL
       RETURN
       END
 
       SUBROUTINE BLOCK_INI
 C...Parameters of fit IFIT=1 of Block and Cahn
       SAVE
 
       COMMON /BLOCKC/ AA, BETA, S0, CC, AMU, DD, ALPHA, A0
       COMMON /BLOCKD/ CP, DP, EP, CM, DM
       AA = 41.74
       BETA = 0.66
       S0 = 338.5
       CC = 0.
       AMU = 0.
       DD = -39.37
       ALPHA = 0.48
       A0 = 0.
       CP = 10.90
       DP = -0.08
       EP = 0.043
       CM = 23.27
       DM = 0.93
       RETURN
       END
 
 C=============================================================
 C.  Nucleus-nucleus cross sections
 C=============================================================
 
       SUBROUTINE SIGNUC_INI (IA,E0)
 C...This subroutine receives in INPUT E0 (TeV)
 C.  energy per nucleon and computes the cross sections
 C.  and interactions lengths for  all nuclei
 C.  with A  between 2 and IA
 C.  The output is contained in common block /CLENNN/
 C.
 C.  Attention: the tabulated cross sections are obtained with
 C.  new p-p cross sections as used in SIBYLL 2x,
 C.  in addition field dimensions changed (RE 04/2000)
 C.
 C........................................................
       COMMON /CLENNN/ SSIGNUC(60), ALNUC(60)
       DIMENSION SIGMA(6,56), SIGQE(6,56)
       DIMENSION AA(6)
       DATA NE /6/, AMIN /1./, DA /1./
       DATA AA /1.,2.,3.,4.,5.,6./
       DATA AVOG /6.0221367E-04/
       DATA ATARGET /14.514/               ! effective masss of air
 C...Data on `inelastic-production' nucleus-air cross section
       DATA (SIGMA(J, 2),J=1,6) /  392., 434., 507., 617., 739., 866./
       DATA (SIGMA(J, 3),J=1,6) /  459., 501., 579., 691., 819., 943./
       DATA (SIGMA(J, 4),J=1,6) /  493., 532., 608., 721., 853., 975./
       DATA (SIGMA(J, 5),J=1,6) /  589., 632., 716., 839., 986.,1122./
       DATA (SIGMA(J, 6),J=1,6) /  695., 746., 838., 974.,1125.,1275./
       DATA (SIGMA(J, 7),J=1,6) /  710., 759., 857., 986.,1141.,1292./
       DATA (SIGMA(J, 8),J=1,6) /  751., 808., 902.,1039.,1192.,1344./
       DATA (SIGMA(J, 9),J=1,6) /  794., 847., 944.,1082.,1241.,1395./
       DATA (SIGMA(J,10),J=1,6) /  808., 857., 958.,1095.,1248.,1411./
       DATA (SIGMA(J,11),J=1,6) /  808., 858., 958.,1094.,1255.,1409./
       DATA (SIGMA(J,12),J=1,6) /  856., 909.,1001.,1144.,1306.,1462./
       DATA (SIGMA(J,13),J=1,6) /  868., 916.,1016.,1166.,1324.,1486./
       DATA (SIGMA(J,14),J=1,6) /  913., 961.,1063.,1214.,1386.,1543./
       DATA (SIGMA(J,15),J=1,6) /  937., 995.,1101.,1240.,1411.,1573./
       DATA (SIGMA(J,16),J=1,6) /  998.,1037.,1157.,1304.,1474.,1644./
       DATA (SIGMA(J,17),J=1,6) /  991.,1044.,1154.,1310.,1472.,1647./
       DATA (SIGMA(J,18),J=1,6) / 1041.,1101.,1208.,1359.,1544.,1718./
       DATA (SIGMA(J,19),J=1,6) / 1094.,1151.,1255.,1431.,1614.,1785./
       DATA (SIGMA(J,20),J=1,6) / 1141.,1193.,1314.,1485.,1667.,1842./
       DATA (SIGMA(J,21),J=1,6) / 1165.,1222.,1341.,1498.,1694.,1880./
       DATA (SIGMA(J,22),J=1,6) / 1159.,1219.,1333.,1506.,1687.,1863./
       DATA (SIGMA(J,23),J=1,6) / 1211.,1265.,1381.,1553.,1738.,1918./
       DATA (SIGMA(J,24),J=1,6) / 1221.,1285.,1402.,1571.,1757.,1947./
       DATA (SIGMA(J,25),J=1,6) / 1245.,1314.,1429.,1595.,1792.,1993./
       DATA (SIGMA(J,26),J=1,6) / 1270.,1334.,1447.,1624.,1811.,2005./
       DATA (SIGMA(J,27),J=1,6) / 1251.,1314.,1431.,1603.,1790.,1984./
       DATA (SIGMA(J,28),J=1,6) / 1289.,1351.,1465.,1654.,1839.,2029./
       DATA (SIGMA(J,29),J=1,6) / 1412.,1482.,1617.,1796.,2003.,2204./
       DATA (SIGMA(J,30),J=1,6) / 1351.,1409.,1526.,1711.,1895.,2105./
       DATA (SIGMA(J,31),J=1,6) / 1361.,1422.,1548.,1734.,1932.,2132./
       DATA (SIGMA(J,32),J=1,6) / 1381.,1435.,1574.,1751.,1952.,2140./
       DATA (SIGMA(J,33),J=1,6) / 1400.,1466.,1593.,1767.,1971.,2173./
       DATA (SIGMA(J,34),J=1,6) / 1414.,1479.,1606.,1796.,1985.,2193./
       DATA (SIGMA(J,35),J=1,6) / 1425.,1503.,1621.,1800.,2008.,2203./
       DATA (SIGMA(J,36),J=1,6) / 1444.,1516.,1643.,1822.,2031.,2225./
       DATA (SIGMA(J,37),J=1,6) / 1460.,1529.,1655.,1844.,2042.,2246./
       DATA (SIGMA(J,38),J=1,6) / 1484.,1547.,1676.,1861.,2071.,2269./
       DATA (SIGMA(J,39),J=1,6) / 1494.,1557.,1700.,1886.,2085.,2291./
       DATA (SIGMA(J,40),J=1,6) / 1510.,1585.,1710.,1897.,2116.,2292./
       DATA (SIGMA(J,41),J=1,6) / 1532.,1594.,1724.,1922.,2116.,2317./
       DATA (SIGMA(J,42),J=1,6) / 1540.,1617.,1735.,1930.,2144.,2334./
       DATA (SIGMA(J,43),J=1,6) / 1559.,1627.,1764.,1944.,2158.,2363./
       DATA (SIGMA(J,44),J=1,6) / 1570.,1641.,1774.,1949.,2172.,2367./
       DATA (SIGMA(J,45),J=1,6) / 1587.,1651.,1784.,1983.,2187.,2396./
       DATA (SIGMA(J,46),J=1,6) / 1609.,1675.,1810.,1997.,2206.,2412./
       DATA (SIGMA(J,47),J=1,6) / 1611.,1688.,1814.,2011.,2236.,2436./
       DATA (SIGMA(J,48),J=1,6) / 1630.,1700.,1833.,2024.,2240.,2457./
       DATA (SIGMA(J,49),J=1,6) / 1645.,1717.,1843.,2047.,2251.,2474./
       DATA (SIGMA(J,50),J=1,6) / 1655.,1729.,1858.,2051.,2272.,2488./
       DATA (SIGMA(J,51),J=1,6) / 1668.,1742.,1875.,2076.,2288.,2504./
       DATA (SIGMA(J,52),J=1,6) / 1684.,1757.,1891.,2091.,2308.,2515./
       DATA (SIGMA(J,53),J=1,6) / 1697.,1773.,1908.,2103.,2328.,2544./
       DATA (SIGMA(J,54),J=1,6) / 1713.,1785.,1917.,2125.,2335.,2533./
       DATA (SIGMA(J,55),J=1,6) / 1726.,1803.,1939.,2131.,2355.,2570./
       DATA (SIGMA(J,56),J=1,6) / 1757.,1831.,1970.,2168.,2388.,2620./
 C...Data on `quasi-elastic' nucleus-air cross section
       DATA (SIGQE(J, 2),J=1,6) /   41.,  39.,  54.,  99., 159., 229./
       DATA (SIGQE(J, 3),J=1,6) /   43.,  40.,  55., 101., 167., 238./
       DATA (SIGQE(J, 4),J=1,6) /   40.,  40.,  55., 102., 168., 243./
       DATA (SIGQE(J, 5),J=1,6) /   45.,  44.,  61., 113., 179., 262./
       DATA (SIGQE(J, 6),J=1,6) /   52.,  49.,  67., 121., 195., 278./
       DATA (SIGQE(J, 7),J=1,6) /   50.,  48.,  69., 121., 194., 276./
       DATA (SIGQE(J, 8),J=1,6) /   52.,  50.,  70., 125., 199., 282./
       DATA (SIGQE(J, 9),J=1,6) /   54.,  50.,  69., 126., 206., 289./
       DATA (SIGQE(J,10),J=1,6) /   53.,  50.,  70., 127., 204., 288./
       DATA (SIGQE(J,11),J=1,6) /   51.,  50.,  69., 125., 203., 285./
       DATA (SIGQE(J,12),J=1,6) /   52.,  51.,  70., 127., 210., 293./
       DATA (SIGQE(J,13),J=1,6) /   53.,  50.,  71., 129., 210., 293./
       DATA (SIGQE(J,14),J=1,6) /   52.,  53.,  72., 134., 212., 301./
       DATA (SIGQE(J,15),J=1,6) /   54.,  52.,  75., 135., 219., 306./
       DATA (SIGQE(J,16),J=1,6) /   57.,  56.,  76., 136., 222., 313./
       DATA (SIGQE(J,17),J=1,6) /   55.,  54.,  76., 137., 221., 311./
       DATA (SIGQE(J,18),J=1,6) /   59.,  55.,  79., 142., 223., 320./
       DATA (SIGQE(J,19),J=1,6) /   63.,  58.,  81., 143., 231., 323./
       DATA (SIGQE(J,20),J=1,6) /   63.,  60.,  82., 149., 237., 332./
       DATA (SIGQE(J,21),J=1,6) /   64.,  61.,  84., 150., 236., 333./
       DATA (SIGQE(J,22),J=1,6) /   62.,  61.,  84., 145., 237., 329./
       DATA (SIGQE(J,23),J=1,6) /   65.,  62.,  83., 149., 240., 341./
       DATA (SIGQE(J,24),J=1,6) /   64.,  61.,  84., 152., 244., 338./
       DATA (SIGQE(J,25),J=1,6) /   64.,  61.,  86., 153., 249., 343./
       DATA (SIGQE(J,26),J=1,6) /   66.,  64.,  88., 156., 251., 344./
       DATA (SIGQE(J,27),J=1,6) /   65.,  61.,  84., 155., 244., 345./
       DATA (SIGQE(J,28),J=1,6) /   65.,  63.,  87., 153., 249., 347./
       DATA (SIGQE(J,29),J=1,6) /   72.,  64.,  91., 161., 260., 362./
       DATA (SIGQE(J,30),J=1,6) /   67.,  63.,  90., 157., 252., 355./
       DATA (SIGQE(J,31),J=1,6) /   67.,  66.,  90., 161., 256., 353./
       DATA (SIGQE(J,32),J=1,6) /   68.,  65.,  89., 161., 256., 355./
       DATA (SIGQE(J,33),J=1,6) /   69.,  66.,  89., 161., 257., 360./
       DATA (SIGQE(J,34),J=1,6) /   69.,  63.,  89., 160., 262., 359./
       DATA (SIGQE(J,35),J=1,6) /   68.,  65.,  90., 166., 259., 360./
       DATA (SIGQE(J,36),J=1,6) /   68.,  66.,  92., 167., 257., 362./
       DATA (SIGQE(J,37),J=1,6) /   71.,  67.,  92., 166., 261., 365./
       DATA (SIGQE(J,38),J=1,6) /   70.,  67.,  92., 167., 265., 367./
       DATA (SIGQE(J,39),J=1,6) /   69.,  68.,  96., 163., 261., 370./
       DATA (SIGQE(J,40),J=1,6) /   71.,  67.,  92., 167., 266., 371./
       DATA (SIGQE(J,41),J=1,6) /   72.,  68.,  93., 167., 267., 373./
       DATA (SIGQE(J,42),J=1,6) /   71.,  69.,  96., 167., 271., 370./
       DATA (SIGQE(J,43),J=1,6) /   72.,  70.,  95., 168., 265., 369./
       DATA (SIGQE(J,44),J=1,6) /   70.,  68.,  98., 170., 270., 372./
       DATA (SIGQE(J,45),J=1,6) /   72.,  69.,  97., 167., 272., 379./
       DATA (SIGQE(J,46),J=1,6) /   75.,  69.,  98., 173., 272., 370./
       DATA (SIGQE(J,47),J=1,6) /   73.,  70.,  93., 170., 275., 377./
       DATA (SIGQE(J,48),J=1,6) /   73.,  71.,  96., 170., 275., 379./
       DATA (SIGQE(J,49),J=1,6) /   75.,  72.,  96., 171., 275., 379./
       DATA (SIGQE(J,50),J=1,6) /   74.,  73.,  97., 174., 275., 380./
       DATA (SIGQE(J,51),J=1,6) /   73.,  68.,  97., 173., 281., 388./
       DATA (SIGQE(J,52),J=1,6) /   75.,  71.,  99., 175., 275., 390./
       DATA (SIGQE(J,53),J=1,6) /   75.,  71.,  98., 180., 279., 383./
       DATA (SIGQE(J,54),J=1,6) /   74.,  71.,  99., 177., 280., 386./
       DATA (SIGQE(J,55),J=1,6) /   75.,  71., 100., 179., 278., 394./
       DATA (SIGQE(J,56),J=1,6) /   76.,  70., 101., 175., 284., 390./
 
       ASQS = 0.5*LOG10(1.876E+03*E0)
       JE = MIN(INT((ASQS-AMIN)/DA)+1,NE-2)
       DO JA=2,IA
          ABEAM = FLOAT(JA)
          S1 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
      +                   SIGMA(JE,JA),SIGMA(JE+1,JA),SIGMA(JE+2,JA))
          S2 = QUAD_INT(ASQS, AA(JE),AA(JE+1),AA(JE+2),
      +                   SIGQE(JE,JA),SIGQE(JE+1,JA),SIGQE(JE+2,JA))
          SSIGNUC(JA) = S1 + S2
          ALNUC(JA) = ATARGET/(AVOG*SSIGNUC(JA))
       ENDDO
       ALNUC(1) = FPNI(E0, 13)
       SSIGNUC(1) = ATARGET/(AVOG*ALNUC(1))
 
       RETURN
       END
 
 
 C=======================================================================
 C.  General utilities
 C.======================================================================
 
       FUNCTION QUAD_INT (R,X0,X1,X2,V0,V1,V2)
 C...Quadratic interpolation
       SAVE
 
       R0=R-X0
       R1=R-X1
       R2=R-X2
       S0=X0-X1
       S1=X0-X2
       S2=X1-X2
       QUAD_INT = V0*R1*R2/(S0*S1)-V1*R0*R2/(S0*S2)+V2*R0*R1/(S1*S2)
       RETURN
       END
 
       FUNCTION GAUSS (FUN, A,B)
 C...Returns the  8 points Gauss-Legendre integral
 C.  of function FUN from A to B
 C...........................................................
       SAVE
 
       DIMENSION X(8), W(8)
       DATA X / .0950125098, .2816035507, .4580167776, .6178762444
      1          ,.7554044083, .8656312023, .9445750230, .9894009349/
       DATA W / .1894506104, .1826034150, .1691565193, .1495959888
      1          ,.1246289712, .0951585116, .0622535239, .0271524594/
       XM = 0.5*(B+A)
       XR = 0.5*(B-A)
       SS = 0.
       DO J=1,8
         DX = XR*X(J)
         SS = SS + W(J) * (FUN(XM+DX) + FUN(XM-DX))
       ENDDO
       GAUSS = XR*SS
       RETURN
       END
 
 ctp060203      FUNCTION GASDEV(Idum)
       FUNCTION GASDEV(Idum)
 C...Gaussian deviation
       SAVE
 
 ctp060123      SAVE GSET
       DATA ISET/0/
       IF (ISET.EQ.0) THEN
 1       V1=2.*S_RNDM(ISET)-1.
         V2=2.*S_RNDM(Idum)-1.
         R=V1**2+V2**2
         IF(R.GE.1.)GO TO 1
         FAC=SQRT(-2.*LOG(R)/R)
         GSET=V1*FAC
         GASDEV=V2*FAC
         ISET=1
       ELSE
         GASDEV=GSET
         ISET=0
       ENDIF
       RETURN
       END
 
       subroutine invert_array (yy, xmin, dx, n, xnew, ymin, dy)
 C..    This subroutine receives one   array
 C      of n y values in input yy(1:n)
 C      that correspond to  equispaced values of x_j = xmin + dx*(j-1)
 C
 C      and "reverse" the array returning an array of  x values
 C      xnew (1:n) that  corresponds to equispaced values of y
 C      The relation is assumed monotonous but can be
 C      increasing or decreasing
 C..............................................................
       SAVE
 
       dimension  yy(n), xnew (n)
       ymin = yy(1)
       ymax = yy(n)
       dy = (ymax - ymin)/float(n-1)
       xnew (1) = xmin
       xnew (n) = xmin + dx*float(n-1)
       k0 = 1
       do j=2,n-1
          y = ymin + float(j-1)*dy
          do k=k0,n
             if((yy(k) .gt. y) .eqv. (yy(n) .gt. yy(1))) goto 100
          enddo
 100      y2 = yy(k)
          y1 = yy(k-1)
          k0 = k-1
          x1 = xmin + dx*float(k-2)
          x2 = x1+dx
          xnew (j)  = x1 + dx* (y-y1)/(y2-y1)
       enddo
       return
       end
 C->
       SUBROUTINE SINCO(S,C)
       SAVE
 
       DATA PI /3.1415926/
       F = 2.*PI*S_RNDM(0)
       C = COS (F)
       S = SIN (F)
       RETURN
       END
 
 
 
 C=============================================================
 C.  Cross sections for cascade calculations (FPNI)
 C=============================================================
 
 
 
       SUBROUTINE SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
 C-----------------------------------------------------------------------
 C...p-p cross sections
 C.
 C.  this routine serves the purpose to calculate cascades with different
 C.  cross sections
 C.
 C. INPUT: E0 = Laboratory Energy  (TeV)
 C.
 C. OUTPUT: SIGT = total cross section
 C.         SIGEL = elastic cross section
 C.         SIGINEL = inelastic cross section
 C.         SLOPE = slope of elastic scattering (GeV**-2)
 C.         RHO = Imaginary/Real part of forward elastic amplitude
 C.
 C.  (old cross section tables end at 10^6 GeV)
 C-----------------------------------------------------------------------
       SAVE
 
       DIMENSION SSIG0(51)
       DIMENSION SIGDIF(3)
       DATA ICSPA /0/
       DATA PI /3.1415926/
       DATA CMBARN /0.389385/
 
 C...p-p inelastic cross sections (mbarn)
       DATA (SSIG0(J),J=1,51) /
      +      32.05,    32.06,    32.08,    32.13,    32.22,    32.36,
      +      32.56,    32.85,    33.24,    33.75,    34.37,    35.14,
      +      36.05,    37.12,    38.37,    39.78,    41.36,    43.13,
      +      45.07,    47.18,    49.47,    51.91,    54.54,    57.28,
      +      60.15,    63.15,    66.28,    69.48,    72.80,    76.22,
      +      79.71,    83.27,    86.87,    90.55,    94.26,    98.05,
      +     101.89,   105.75,   109.71,   113.65,   117.60,   121.55,
      +     125.53,   129.56,   133.60,   137.70,   141.77,   145.84,
      +     149.92,   154.02,   158.15/
 
 
       SQS = SQRT(2000.*0.938*E0)
 
 *  old standard NUCLIB/SIBYLL model
 
       IF(ICSPA.EQ.-1) THEN
 
         AL = LOG10(SQS)
         if(AL.le.1.) then
           SIGINEL = SSIG0(1)
         else
           J1 = (AL - 1.)*10. + 1
           J1 = min(J1,50)
           T = (AL-1.)*10. - FLOAT(J1-1)
           SIGINEL = SSIG0(J1)*(1.-T) + SSIG0(J1+1)*T
         endif
         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
         R = SIGEL1/SIGT1
         RHO = RHO1
         SIGT  = SIGINEL/(1.-R)
         SIGEL = SIGINEL*R/(1.-R)
         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO1**2) /CMBARN
 
 *  cross section as calculated in SIBYLL
 
       ELSE IF(ICSPA.EQ.0) THEN
 
         CALL SIB_SIGMA_HP(1,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
 
 *  Donnachie-Landshoff  (sig-tot)
 
       ELSE IF(ICSPA.EQ.1) THEN
 
         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
      +             SIGEL1,SIGEL2)
         R = SIGEL1/SIGT1
         RHO = RHO1
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 21.7*S**DELDL+56.08*S**EPSDL
         SIGEL = R*SIGT
         SIGINEL = SIGT-SIGEL
         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO**2) /CMBARN
 
 *  Donnachie-Landshoff (sig-tot and sig-el)
 
       ELSE IF(ICSPA.EQ.2) THEN
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 21.7*S**DELDL+56.08*S**EPSDL
         IMODEL = 1
         IF(IMODEL.EQ.1) THEN
           ALPHAP = 0.25D0
           SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
         ELSE IF(IMODEL.EQ.2) THEN
           ALPHAP = 0.3D0
           SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
         ENDIF
         SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
         SIGINEL = SIGT-SIGEL
         RHO = 0.
 
 *  geometrical scaling with Donnachie-Landshoff sig-tot
 
       ELSE IF(ICSPA.EQ.3) THEN
 
         R = 0.17D0
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 21.7*S**DELDL+56.08*S**EPSDL
 
         SIGEL = R*SIGT
         SIGINEL = SIGT-SIGEL
         SLOPE = SIGT**2/(16*PI*SIGEL)/CMBARN
         RHO = 0.
 
       ENDIF
 
       RETURN
       END
 
 
       SUBROUTINE SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
 C-----------------------------------------------------------------------
 C...pi-p cross sections
 C.
 C.  this routine serves the purpose to calculate cascades with different
 C.  cross sections
 C.
 C. INPUT: E0 = Laboratory Energy  (TeV)
 C.
 C. OUTPUT: SIGT = total cross section
 C.         SIGEL = elastic cross section
 C.         SIGINEL = inelastic cross section
 C.         SLOPE = slope of elastic scattering (GeV**-2)
 C.         RHO = Imaginary/Real part of forward elastic amplitude
 C.
 C.  (old cross section tables end at 10^6 GeV)
 C-----------------------------------------------------------------------
       SAVE
 
       DIMENSION SSIG0(51)
       DIMENSION SIGDIF(3)
       DATA ICSPA /0/
       DATA PI /3.1415926/
       DATA CMBARN /0.389385/
 
 C...pi-p inelastic cross sections (mbarn)
       DATA (SSIG0(J),J=1,51) /
      +      20.76,    20.78,    20.81,    20.88,    20.98,    21.13,
      +      21.33,    21.61,    21.96,    22.39,    22.92,    23.56,
      +      24.31,    25.18,    26.18,    27.32,    28.60,    30.04,
      +      31.64,    33.40,    35.34,    37.43,    39.72,    42.16,
      +      44.77,    47.56,    50.53,    53.66,    56.99,    60.50,
      +      64.17,    68.03,    72.05,    76.27,    80.67,    85.27,
      +      90.08,    95.04,   100.27,   105.65,   111.21,   116.94,
      +     122.87,   129.03,   135.37,   141.93,   148.62,   155.49,
      +     162.48,   169.60,   176.94/
 
       SQS = SQRT(2000.*0.938*E0)
 
 *  old standard NUCLIB/SIBYLL model
 
       IF(ICSPA.EQ.-1) THEN
 
         AL = LOG10(SQS)
         if(AL.le.1.) then
           SIGINEL = SSIG0(1)
         else
           J1 = (AL - 1.)*10. + 1
           J1 = min(J1,50)
           T = (AL-1.)*10. - FLOAT(J1-1)
           SIGINEL = SSIG0(J1)*(1.-T) + SSIG0(J1+1)*T
         endif
         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
         R = SIGEL1/SIGT1
         RHO = RHO1
         SIGT  = SIGINEL/(1.-R)
         SIGEL = SIGINEL*R/(1.-R)
         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO1**2) /CMBARN
 
 *  cross section as calculated in SIBYLL
 
       ELSE IF(ICSPA.EQ.0) THEN
 
         CALL SIB_SIGMA_HP(2,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
 
 *  Donnachie-Landshoff  (sig-tot)
 
       ELSE IF(ICSPA.EQ.1) THEN
 
         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
      +             SIGEL1,SIGEL2)
         R = SIGEL1/SIGT1
         RHO = RHO1
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 13.63*S**DELDL+(36.02+27.56)/2.*S**EPSDL
         SIGEL = R*SIGT
         SIGINEL = SIGT-SIGEL
         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO**2) /CMBARN
 
 *  Donnachie-Landshoff (sig-tot and sig-el)
 
       ELSE IF(ICSPA.EQ.2) THEN
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 13.63*S**DELDL+(36.02+27.56)/2.*S**EPSDL
         IMODEL = 1
         IF(IMODEL.EQ.1) THEN
           ALPHAP = 0.25D0
           SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
         ELSE IF(IMODEL.EQ.2) THEN
           ALPHAP = 0.3D0
           SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
         ENDIF
         SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
         SIGINEL = SIGT-SIGEL
         RHO = 0.
 
 *  geometrical scaling with Donnachie-Landshoff sig-tot
 
       ELSE IF(ICSPA.EQ.3) THEN
 
         R = 0.17D0
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 13.63*S**DELDL+(36.02+27.56)/2.*S**EPSDL
 
         SIGEL = R*SIGT
         SIGINEL = SIGT-SIGEL
         SLOPE = SIGT**2/(16*PI*SIGEL)/CMBARN
         RHO = 0.
 
       ENDIF
 
       RETURN
       END
 
 
       SUBROUTINE SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
 C-----------------------------------------------------------------------
 C...K-p cross sections
 C.
 C.  this routine serves the purpose to calculate cascades with different
 C.  cross sections
 C.
 C.  if old cross sections are selected then sigma_pi = sigma_K
 C.
 C. INPUT: E0 = Laboratory Energy  (TeV)
 C.
 C. OUTPUT: SIGT = total cross section
 C.         SIGEL = elastic cross section
 C.         SIGINEL = inelastic cross section
 C.         SLOPE = slope of elastic scattering (GeV**-2)
 C.         RHO = Imaginary/Real part of forward elastic amplitude
 C.
 C.  (old cross section tables end at 10^6 GeV)
 C-----------------------------------------------------------------------
       SAVE
 
       DIMENSION SSIG0(51)
       DIMENSION SIGDIF(3)
       DATA ICSPA /0/
       DATA PI /3.1415926/
       DATA CMBARN /0.389385/
 
 C...pi-p inelastic cross sections (mbarn)
       DATA (SSIG0(J),J=1,51) /
      +      20.76,    20.78,    20.81,    20.88,    20.98,    21.13,
      +      21.33,    21.61,    21.96,    22.39,    22.92,    23.56,
      +      24.31,    25.18,    26.18,    27.32,    28.60,    30.04,
      +      31.64,    33.40,    35.34,    37.43,    39.72,    42.16,
      +      44.77,    47.56,    50.53,    53.66,    56.99,    60.50,
      +      64.17,    68.03,    72.05,    76.27,    80.67,    85.27,
      +      90.08,    95.04,   100.27,   105.65,   111.21,   116.94,
      +     122.87,   129.03,   135.37,   141.93,   148.62,   155.49,
      +     162.48,   169.60,   176.94/
 
       SQS = SQRT(2000.*0.938*E0)
 
 *  old standard NUCLIB/SIBYLL model
 
       IF(ICSPA.EQ.-1) THEN
 
         AL = LOG10(SQS)
         if(AL.le.1.) then
           SIGINEL = SSIG0(1)
         else
           J1 = (AL - 1.)*10. + 1
           J1 = min(J1,50)
           T = (AL-1.)*10. - FLOAT(J1-1)
           SIGINEL = SSIG0(J1)*(1.-T) + SSIG0(J1+1)*T
         endif
         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,SIGEL1,SIGEL2)
         R = SIGEL1/SIGT1
         RHO = RHO1
         SIGT  = SIGINEL/(1.-R)
         SIGEL = SIGINEL*R/(1.-R)
         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO1**2) /CMBARN
 
 *  cross section as calculated in SIBYLL
 
       ELSE IF(ICSPA.EQ.0) THEN
 
         CALL SIB_SIGMA_HP(3,SQS,SIGT,SIGEL,SIGINEL,SIGDIF,SLOPE,RHO)
 
 *  Donnachie-Landshoff  (sig-tot)
 
       ELSE IF(ICSPA.EQ.1) THEN
 
         CALL BLOCK(SQS,SIGT1,SIGT2,SLOP1,SLOP2,RHO1,RHO2,
      +             SIGEL1,SIGEL2)
         R = SIGEL1/SIGT1
         RHO = RHO1
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 11.82*S**DELDL+(26.36+ 8.15)/2.*S**EPSDL
         SIGEL = R*SIGT
         SIGINEL = SIGT-SIGEL
         SLOPE = SIGT**2/(SIGEL * 16.*PI) * (1.+RHO**2) /CMBARN
 
 *  Donnachie-Landshoff (sig-tot and sig-el)
 
       ELSE IF(ICSPA.EQ.2) THEN
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 11.82*S**DELDL+(26.36+ 8.15)/2.*S**EPSDL
         IMODEL = 1
         IF(IMODEL.EQ.1) THEN
           ALPHAP = 0.25D0
           SLOPE = 8.5D0+2.D0*ALPHAP*LOG(S)
         ELSE IF(IMODEL.EQ.2) THEN
           ALPHAP = 0.3D0
           SLOPE = 8.D0+2.D0*ALPHAP*LOG(S)
         ENDIF
         SIGEL = SIGT**2/(16.D0*PI*SLOPE*CMBARN)
         SIGINEL = SIGT-SIGEL
         RHO = 0.
 
 *  geometrical scaling with Donnachie-Landshoff sig-tot
 
       ELSE IF(ICSPA.EQ.3) THEN
 
         R = 0.17D0
 
         DELDL = 0.0808
         EPSDL = -0.4525
         S = SQS*SQS
         SIGT = 11.82*S**DELDL+(26.36+ 8.15)/2.*S**EPSDL
 
         SIGEL = R*SIGT
         SIGINEL = SIGT-SIGEL
         SLOPE = SIGT**2/(16*PI*SIGEL)/CMBARN
         RHO = 0.
 
       ENDIF
 
       RETURN
       END
 
 
 
       SUBROUTINE SIGMA_INI
 C-----------------------------------------------------------------------
 C.  Initialize the cross section and interaction lengths on air
 C-----------------------------------------------------------------------
       SAVE
 
       COMMON /CSAIR/ NSQS, ASQSMIN, ASQSMAX, DASQS,
      &               SSIG0(61,3),SSIGA(61,3),ALINT(61,3)
 
       DATA AVOG /6.0221367E-04/
       ATARGET = 14.514
 
       CALL BLOCK_INI
 
 C...Loop on c.m. energy
       NSQS = 61
       SQSMIN = 10.
       SQSMAX = 1.E+07
       ASQSMIN = LOG10(SQSMIN)
       ASQSMAX = LOG10(SQSMAX)
       DASQS = (ASQSMAX-ASQSMIN)/FLOAT(NSQS-1)
       DO J=1,NSQS
          ASQS = ASQSMIN + DASQS*FLOAT(J-1)
          SQS = 10.**ASQS
          E0 = SQS*SQS/(2.*0.938) * 1.E-03       ! TeV
 C...p-air
          CALL SIGMA_PP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
          CALL SIG_H_AIR (SIGT, SLOPE, RHO, SSIGT, SSIGEL, SSIGQE)
          SSIGA(J,1) = SSIGT-SSIGQE
          SSIG0(J,1) = SIGINEL
          ALINT(J,1) = 1./(AVOG*SSIGA(J,1)/ATARGET)
 C...pi-air
          CALL SIGMA_PIP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
          CALL  SIG_H_AIR (SIGT, SLOPE, RHO, SSIGT, SSIGEL, SSIGQE)
          SSIGA(J,2) = SSIGT-SSIGQE
          SSIG0(J,2) = SIGINEL
          ALINT(J,2) = 1./(AVOG*SSIGA(J,2)/ATARGET)
 C...K-air
          CALL SIGMA_KP (E0, SIGT, SIGEL, SIGINEL, SLOPE, RHO)
          CALL  SIG_H_AIR (SIGT, SLOPE, RHO, SSIGT, SSIGEL, SSIGQE)
          SSIGA(J,3) = SSIGT-SSIGQE
          SSIG0(J,3) = SIGINEL
          ALINT(J,3) = 1./(AVOG*SSIGA(J,3)/ATARGET)
       ENDDO
 
 *     WRITE(6,'(1X,A)')
 *    &  'SIGMA_INI: NUCLIB interaction lengths (p-air, pi-air, K-air)'
 *     DO J=1,NSQS
 *        SQS = 10.**(ASQSMIN + DASQS*FLOAT(J-1))
 *        WRITE(6,'(1X,1P,4E12.3)') SQS,ALINT(J,1),ALINT(J,2),ALINT(J,3)
 *     ENDDO
 
       RETURN
       END
 
 
       FUNCTION FPNI (E,Linp)
 C-----------------------------------------------------------------------
 C...This function  returns the interaction length
 C.  of an hadronic particle travelling in air
 C.
 C.  INPUT:   E (TeV)   particle energy
 C.           Linp      particle code
 C.  OUTPUT:  FPNI      (g cm-2)
 C-----------------------------------------------------------------------
       SAVE
 
       COMMON /CSAIR/ NSQS, ASQSMIN, ASQSMAX, DASQS,
      &               SSIG0(61,3),SSIGA(61,3),ALINT(61,3)
 
       DIMENSION KK(6:14)
       DATA KK /3*2, 4*3, 2*1/
 
       SQS = SQRT(2000.*E*0.937)                        ! GeV
       AL = LOG10 (SQS)
       L = abs(Linp)
       IF (AL .LE. ASQSMIN)  THEN
          FPNI = ALINT(1,KK(L))
       ELSE
          T = (AL-ASQSMIN)/DASQS
          J = INT(T)
          J = MIN(J,NSQS-2)
          T = T-FLOAT(J)
          FPNI = ((1.-T)*ALINT(J+1,KK(L)) + T*ALINT(J+2,KK(L)))
       ENDIF
 
       RETURN
       END
 
 
 
       SUBROUTINE INT_LEN_INI
 C-----------------------------------------------------------------------
 C...Initialize the interaction lengths from NUCLIB
 C-----------------------------------------------------------------------
       SAVE
 
       CALL NUC_GEOM_INI                 ! nucleus profiles
       CALL SIGMA_INI                    ! initialize cross sections
 
       RETURN
       END