!**********************************************
!*  Declarataion of the Communication Module  *
!**********************************************
module comm
    parameter (max_entry=10000)
    integer i_entry(0:max_entry)
    integer i_entry_max
end module comm

!*******************************************
!*     Declarataion of the Path module     *
!*******************************************
module path_module
    character*256 files_path
    character*10  re_import_test
end module path_module


!*****************************************
!*     Declarataion of the ISO Module    *
!*****************************************
module iso
   implicit none
   character*256 init_file, input_file_network, input_file_abundance, output_file_abundance
   integer max_number_isotopes, number_isotopes, current_isotope, abundance_file_format
   character*20  element_names(0:200)
   integer number_reaction_types, max_number_reactions, max_number_fissile_isotopes
   double precision max_hwz ! maximum half life - above that, isotope will be considered stable
   parameter (number_reaction_types=22)
   parameter (max_number_isotopes=10000)  ! maximum number of species
   parameter (max_number_reactions=100)  ! maximum number of reactions per species
   parameter (max_number_fissile_isotopes=200) ! maximum number of species which undergo fission
   parameter (max_hwz=1e20)  ! maximum half life - above that, isotope will be considered stable - ~3000 Gyr
   integer z(max_number_isotopes), n(max_number_isotopes)  ! Number of Protons and Neutrons of given species
   integer reactions(max_number_isotopes,-1:max_number_reactions)  ! 0: Number of reactions 1..n reaction type of given species, -1 index in fission vector
   double precision decay_constant(max_number_isotopes,0:max_number_reactions)  ! decay rates of reactions of given species in 1/s  0 ... total, 1..n partial (later normalized to 1)
   double precision decay(max_number_isotopes)  ! sum of decays per time step of given species
   double precision production(max_number_isotopes)  ! sum of production per time step of given species
   double precision level(max_number_isotopes)  ! level in keV above ground state of given species (allows isomers...)
   double precision spin(max_number_isotopes)  ! spin of given species (needed in case of fission)
   double precision level_product(max_number_isotopes,max_number_reactions)  ! level in keV above ground state of product species (allows isomers...)
   integer product_isomer(max_number_isotopes,max_number_reactions)  ! species produced by reaction on (isomer,reaction)
   double precision abundance(max_number_isotopes)  ! abundance of given species 
   double precision initial_abundance(max_number_isotopes)  ! abundance of given species 
   character*3 reaction_types(number_reaction_types)
   integer reaction_vector(2,number_reaction_types) 
!  gives change in Z,N for desired reaction - EC is (-1,1), B- is (+1,-1) etc.
   data reaction_types  /'B-',   'EC',   'N',     'P',    'A',   'BN' , 'EP',   &          ! 7
                       'BA',   'EA',                                          &          ! 9
                       '2N',   '2P',   '2A',    'B2A',                        &          ! 13
                       'B2N',  'B3N',  'B4N', 'E2P', 'BNA',   'EPA', 'IT' , '12C',   &   ! 21
                       'SF'/  ! 22
   data reaction_vector /1,-1,   -1,1,   0,-1,   -1,0,   -2,-2,  1,-2,  -2,1,   &
                       -1,-3,  -3,-1,                                         &  ! in the current file, all BA is in fact B-A
                       +0,-2,  -2,0,   -4,-4,   -3,-5,                        &
                       +1,-3,  +1,-4,  +1,-5,   -3,1,   -1,-4,   -4,1,  0,0,   -6,-6,  &
                       +0,0/  
! spontanous fission - requires special treatment
   double precision s_fission_vector(0:max_number_fissile_isotopes,max_number_isotopes)  ! fission yields per fissile isotope
   double precision time_decay, dt
   double precision kt, mb_distribution(1000)  ! Maxwell-Boltzmann Distribution, need for fission/scission
   parameter (kt=1.8d0)  ! temperature of neutron before scission in MeV
   integer index_p, index_n, index_a, index_12c, time_steps, short_lived_counter
   logical short_lived(max_number_isotopes)

end module iso


!*****************************************
!*     Declarataion of the ISO Module    *
!*****************************************
! Taken from GEFSUB.FOR
module GEFSUB_FOR

!     '    Copyright 2009,2010,2011,2012,2013,2014,2015: 
!     '       Dr. Karl-Heinz Schmidt,Rheinstrasse 4,64390 Erzhausen,Germany 
!     '       and 
!     '       Dr. Beatriz Jurado,Centre d'Etudes Nucleaires de Bordeaux-Gradignan, 
!     '       Chemin du Solarium,Le Haut Vigneau,BP 120,33175 Gradignan,Cedex, 
!     '       France 
!     ' 
!     '    This program is free software: you can redistribute it and/or modify 
!     '    it under the terms of the GNU General Public License as published by 
!     '    the Free Software Foundation,either version 3 of the License,or 
!     '    (at your option) any later version. 
!     ' 
!     '    This program is distributed in the hope that it will be useful, 
!     '    but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     '    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
!     '    GNU General Public License for more details. 
!     ' 
!     '    You should have received a copy of the GNU General Public License 
!     '    along with this program.  If not,see <http://www.gnu.org/licenses/>. 
!     ' 
!     ' 
!     /' Documentation: '/ 
!     /' (1) K.-H. Schmidt and B. Jurado,Contribution to '/ 
!     /'     ESNT Workshop "The scission process",Saclay (France),April 12-16,2010 '/ 
!     /' (2) B. Jurado and K.-H. Schmidt,Contribution to '/ 
!     /'     Seminar an fission,Gent (Belgium),May 17-20,2010 '/ 
!     /' (3) K.-H. Schmidt and B. Jurado,Contribution to '/ 
!     /'     Seminar on fission,Gent (Belgium),May 17-20,2010 '/ 
!     /' (4) B. Jurado and K.-H. Schmidt,Contribution to '/ 
!     /'     EFNUDAT Workshop,Paris (France),May 25-27,2010 '/ 
!     /' (5) K.-H. Schmidt and B. Jurado,Contribution to '/ 
!     /'     EFNUDAT Workshop,Paris (France),May 25-27,2010 '/ 
!     /' (6) K.-H. Schmidt and B. Jurado,'/ 
!     /'     Final Report to EFNUDAT,October,2010 '/ 
!     /' (7) K.-H. Schmidt and B. Jurado,Phys. Rev. Lett. 104 (2010) 21250 '/ 
!     /' (8) K.-H. Schmidt and B. Jurado,Phys. Rev. C 82 (2011) 014607 '/ 
!     /' (9) K.-H. Schmidt and B. Jurado,Phys. Rev. C 83 (2011) 061601 '/ 
!     /' (10) K.-H. Schmidt and B. Jurado,arXiv:1007.0741v1[nucl-th] (2010) '/ 
!     /' (11) K.-H. Schmidt and B. Jurado,JEF/DOC 1423,NEA of OECD,2012 '/ 
!     /' (12) K.-H. Schmidt and B. Jurado,Phys. Rev. C 86 (2012) 044322 '/ 
!     /' (13) K.-H. Schmidt,B. Jurado,Ch. Amouroux,JEFF-Report 24,NEA of OECD,2014 '/ 
!     ' 
!     ' 
!     /' Further documentation and the newest version of the GEF code are '/ 
!     /' available from                                                   '/ 
!     /' http://www.cenbg.in2p3.fr/GEF and http://www.khs-erzhausen.de/ . '/ 
!     ' 
!     ' 
!     '    The development of the GEF code has been supported by the European Union, 
!     '    EURATOM 6 in the Framework Program "European Facilities for Nuclear Data 
!     '    Measurements" (EFNUDAT),contract number FP6-036434,the Framework 
!     '    Program "European Research Infrastructure for Nuclear Data Applications 
!     '    (ERINDA),contract number FP7-269499,and by the OECD Nuclear Energy Agency. 

contains

!*****************************************
!*               U Valid                 *
!*****************************************
INTEGER*4 function U_Valid(I_Z,I_A)
      IMPLICIT NONE
      INTEGER*4 I_Z
      INTEGER*4 I_A
      INTEGER*4  Ivalid
      Ivalid = 1 
!     '   If I_A / I_Z < 210.E0/90.E0 
      IF (  I_A / I_Z .LT. 172.E0 / 80.E0 .OR. I_A / I_Z .GT. 250.E0/90.E0  ) THEN 
      Ivalid = 0 
      End If 
      IF (  I_Z .LT. 70 .OR. I_Z .GT. 120  ) THEN 
      Ivalid = 0 
      End If 
!     ' Ivalid = 1 
      U_Valid = Ivalid 

end


!*****************************************
!*              U Delta S0               *
!***************************************** 
REAL*4 function get_U_Delta_S0(I_Z,I_A)

      IMPLICIT NONE
      INTEGER*4 I_Z
      INTEGER*4 I_A
!     ' I_Z and I_A refer to the fissioning nucleus90 22 
      REAL*4  Delta
      Delta = 0.3 
      IF (  I_Z .EQ. 90 .AND. I_A .EQ. 228  )  Delta = 0.70 
!     'N 
      IF (  I_Z .EQ. 90 .AND. I_A .EQ. 230  )  Delta = 0.6 
!     'N 
      IF (  I_Z .EQ. 90 .AND. I_A .EQ. 233  )  Delta = 0.3 
!     ' 
      IF (  I_Z .EQ. 91 .AND. I_A .EQ. 228  )  Delta = 0.65 
!     ' 
      IF (  I_Z .EQ. 92 .AND. I_A .EQ. 233  )  Delta = 0.5 
!     'N 
      IF (  I_Z .EQ. 92 .AND. I_A .EQ. 234  )  Delta = 0.6 
!     'N 
      IF (  I_Z .EQ. 92 .AND. I_A .EQ. 235  )  Delta = 0.3 
      IF (  I_Z .EQ. 92 .AND. I_A .EQ. 236  )  Delta = 0.3 
!     'N 
      IF (  I_Z .EQ. 92 .AND. I_A .EQ. 237  )  Delta = 0.3 
      IF (  I_Z .EQ. 92 .AND. I_A .EQ. 238  )  Delta = 0.3 
      IF (  I_Z .EQ. 92 .AND. I_A .EQ. 239  )  Delta = 0.1 
!     ' 
      IF (  I_Z .EQ. 93 .AND. I_A .EQ. 238  )  Delta = -0.1 
!     'N 
!     ' 
      IF (  I_Z .EQ. 94 .AND. I_A .EQ. 240  )  Delta = -0.1 
!     'N 
      IF (  I_Z .EQ. 94 .AND. I_A .EQ. 241  )  Delta = -0.5 
!     'N 
      IF (  I_Z .EQ. 94 .AND. I_A .EQ. 242  )  Delta = -0.15 
!     'N 
      IF (  I_Z .EQ. 94 .AND. I_A .EQ. 243  )  Delta = -0.45 
!     'N 
      IF (  I_Z .EQ. 94  )  Delta = 0.25 
!     ' 
      IF (  I_Z .EQ. 95 .AND. I_A .EQ. 242  )  Delta = -0.35 
!     'N 
!     ' 
      IF (  I_Z .EQ. 95 .AND. I_A .EQ. 243  )  Delta = -0.1 
!     'N 
!     ' 
      IF (  I_Z .EQ. 95 .AND. I_A .EQ. 244  )  Delta = -0.1 
!     ' 
      IF (  I_Z .EQ. 96 .AND. I_A .EQ. 244  )  Delta = 0.0
!     'N 
      IF (  I_Z .EQ. 96 .AND. I_A .EQ. 246  )  Delta = -0.2 
!     'N 
      get_U_Delta_S0 = Delta 

end


!*****************************************
!*             Get Getyield              *
!***************************************** 
REAL*4 function get_Getyield(E_rel,E_ref,T_low,T_high)

      IMPLICIT NONE
      REAL*4 E_rel
      REAL*4 E_ref
      REAL*4 T_low
      REAL*4 T_high
!     /' Erel: Energy relative to the barrier '/ 
!     /' T_low: Effective temperature below barrier '/ 
!     /' T_high: Effective temperature above barrier '/ 
      REAL*4  Exp1
      REAL*4  Yield
!     ' 
      Exp1 = E_rel/T_low - E_ref/0.4 
!     ' energy far below barrier 
!     ' Subtraction of E_ref/0.4 to prevent numerical problems. 
      IF (  Exp1 .LT. -50  ) THEN 
      Yield = 0 
      Else 
      Yield = Exp(E_rel / T_high - E_ref/0.4) * 1.E0 / (1.E0 + exp(-E_rel/ (T_high*T_low/(T_high-T_low) ) ) ) 
      End If 

      get_Getyield = Max(Yield,0.0) 

end
 

!*****************************************
!*               Get F1                  *
!***************************************** 
REAL*4 function get_F1(Z_S_A)
      IMPLICIT NONE
      REAL*4 Z_S_A
!     /' Fit to the lower part of the data '/ 
      REAL*4  Result
      Result = exp(-9.05E0 + 4.58E0 * Log(Z_S_A/2.3E0)) 
      get_F1 = Result 

end


!*****************************************
!*               Get F2                  *
!*****************************************
REAL*4 function get_F2(Z_S_A)

      IMPLICIT NONE
      REAL*4 Z_S_A
!     /' Fit to the upper part of the data '/ 
      REAL*4  Result
      Result = exp(12.08E0 - 3.27E0 * Log(Z_S_A/2.3E0)) 
      get_F2 = Result 

end


!*****************************************
!*             Get Masscurv              *
!*****************************************
REAL*4 function get_Masscurv(Z,A,RL,kappa)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 RL
      REAL*4 kappa
!     /'  Fit to  Data of Fig. 7 of                                             '/ 
!     /'  "Shell effect in the symmetric-modal fission of pre-actinide nuclei"  '/ 
!     /'  S. I. Mulgin,K.-H. Schmidt,A. Grewe,S. V. Zhdanov                  '/ 
!     /'  Nucl. Phys. A 640 (1998) 375 
!     /' (From fit of the width of the mass distributions.) '/                                         '/ 
      REAL*4  RI,Result1,Result2,Result
      REAL*4  Z_square_over_A
      REAL*4  ZsqrA
      REAL*4  c_rot
      DATA c_rot/600.0/
      REAL*4 F1 
      REAL*4 F2 
!     ' 
      Z_square_over_A = Z**2/A 
      RI = (A - 2*Z)/A 
      ZsqrA = Z_square_over_A * (1.E0 - kappa * RI**2) / &
          (1.E0 - kappa * ((226.E0 - 2.E0*91.E0)/226.E0)**2) + &
           c_rot * RL**2 / A**(7.0/3.0) 
!     ' Hasse & Myers 
!     '      + 0.0017 * RL^2 
!     ' 
      Result1 = get_F1(ZsqrA) 
      Result2 = get_F2(ZsqrA) 
      Result = Min(Result1,Result2) 
      get_Masscurv = Result 
!     ' 

end


!*****************************************
!*             Get Masscurv1             *
!*****************************************
REAL*4 function get_Masscurv1(Z,A,RL,kappa)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 RL
      REAL*4 kappa
!     /'  Fit to  Data of Fig. 7 of                                             '/ 
!     /'  "Shell effect in the symmetric-modal fission of pre-actinide nuclei"  '/ 
!     /'  S. I. Mulgin,K.-H. Schmidt,A. Grewe,S. V. Zhdanov                  '/ 
!     /'  Nucl. Phys. A 640 (1998) 375 
!     /' (The left part assumed to be valid for the yields of the fission channels.) '/                                         '/ 
      REAL*4  RI,Result1,Result2,Result
!     '    Dim As Single A,A_central,Z 
      REAL*4  Z_square_over_A
      REAL*4  ZsqrA
      REAL*4  c_rot
      DATA c_rot/600.0/
      REAL*4 F1 
      REAL*4 F2 
!     ' 
!     'A_central = -28.8156 + Z * 2.86587  ' Stability line for heavy nuclei 
!     ' 
      Z_square_over_A = Z**2/A 
      RI = (A - 2*Z)/A 
      ZsqrA = Z_square_over_A * (1.E0 - kappa * RI**2) / &
          (1.E0 - kappa * ((226.E0 - 2.E0*91.E0)/226.E0)**2) + &
           c_rot * RL**2 / A**(7.0/3.0) 
!     ' Hasse & Myers 
!     '      + 0.0017 * RL^2 
!     ' 
      IF (  ZsqrA .LT. 36.0  ) THEN 
!     ' adjusted to Y(S2) in light nuclei (80<Z<92) 
      ZsqrA = ZsqrA + 0.9 * (36.0 - ZsqrA) 
      End If 
!     ' 
      Result1 = get_F1(ZsqrA) 
!     '  Result2 = get_F2(ZsqrA) 
!     '  Result = Min(Result1,Result2) 
      get_Masscurv1 = Result1 
!     ' 

end


!*****************************************
!         Get De Saddle Scission         *
!***************************************** 
REAL*4 function get_De_Saddle_Scission(Z_square_over_Athird,ESHIFTSASCI)

      IMPLICIT NONE
      REAL*4 Z_square_over_Athird
      REAL*4 ESHIFTSASCI
!     /' Energy release between saddle and scission '/ 
!     /' M. Asghar,R. W. Hasse,J. Physique C 6 (1984) 455 '/ 
      REAL*4  Result
      Result = (31.E0 - 11.E0) / (1550.E0 - 1300.E0) * &
          (Z_square_over_Athird - 1300.E0 + ESHIFTSASCI) + 11.E0 
!     ' This formula with ESHIFTSASCI = 0 is the parameterisation of the results 
!     ' of Ashgar and Hasse,JPC 6 (1984) 455,see 
!     ' F. Rejmund,A. V. Ignatyuk,A. R. Junghans,K.-H. Schmidt 
!     ' Nucl. Phys. A 678 (2000) 215 
      Result = max(Result,0.0) 
      get_De_Saddle_Scission = Result 

end


!*****************************************
!*             Get TEgidy                *
!***************************************** 
REAL*4 function get_TEgidy(A,DU,Fred)

      IMPLICIT NONE
      REAL*4 A
      REAL*4 DU
      REAL*4 Fred
!     /' Temperature parameter of the constant-temperature formula for the 
!     nuclear level density. 
!     Input parameters: A = Mass number of nucleus 
!     DU = Shell effect (corrected for pairing:P=0 for odd-A nuclei) 
!     From "Correlations between the nuclear level density parameters" 
!     Dorel Bucurescu,Till von Egidy 
!     Phys. Rev. C 72 (2005) 067304    and 
!     "Systematics of nuclear level density parameters" 
!     Dorel Bucurescu,Till von Egidy 
!     J. Phys. G: Nucl. Part. Phys. 31 (2005) S1675 and 
!     "Systematics of nuclear level density parameters" 
!     Till von Egidy,Dorel Bucurescu 
!     Phys. Rev. C 72 (2005) 044311 '/ 
      REAL*4  Temp_smooth,Temp,T_Fac
!     ' Temp_smooth = 17.45E0 / (A^0.666667E0) 
!     ' Temp = (17.45E0 - 0.51E0 * DU + 0.051 * DU^2) / (A^0.666667E0) 
      Temp_smooth = 1.0 / (0.0570 * A**0.6666667) 
      Temp = 1.0 / ( (0.0570 + 0.00193*DU) * A**0.6666667) 
!     ' from  PRC 80 (2009) 054310 
      T_Fac = Temp / Temp_smooth 
      Temp = Temp * Fred 
!     /' (For influence of deformation) '/ 
      get_TEgidy = Temp 

end


!*****************************************
!*             Get TRusanov              *
!***************************************** 
REAL*4 function get_TRusanov(E,A)

      IMPLICIT NONE
      REAL*4 E
      REAL*4 A
!     /' Fermi-gas level density,parameterisation of Rusanov et al. '/ 
      IF (  E >0  ) THEN 
      get_TRusanov = SQRT(E / (0.094E0 * A) ) 
      Else 
      get_TRusanov = 0.0 
      End If 

end


!*****************************************
!*              Get LyMass               *
!*****************************************
REAL*4 function get_LyMass(Z,A,beta)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 beta
!     ' 
!     /' liquid-drop mass,Myers & Swiatecki,Lysekil,1967  '/ 
!     /' pure liquid drop,without pairing and shell effects '/ 
!     ' 
!     /' On input:    Z     nuclear charge of nucleus        '/ 
!     /'              N     number of neutrons in nucleus    '/ 
!     /'              beta  deformation of nucleus           '/ 
!     /' On output:   binding energy of nucleus              '/ 
!     ' 
      REAL*4  pi
      PARAMETER (pi=3.14159)
      REAL*4  N
      REAL*4  alpha
      REAL*4  XCOM,XVS,XE,EL
!     ' 
      N = A - Z 
      alpha = SQRT(5.E0/(4.E0*pi)) * beta 
      XCOM = 1.E0 - 1.7826E0 * ((A - 2.E0*Z)/A)**2 
!     /' factor for asymmetry dependence of surface and volume term '/ 
      XVS = - XCOM * (15.4941E0*A - 17.9439E0*A**(2.E0/3.E0)*(1.E0+0.4E0*Alpha**2)) 
!     /' sum of volume and surface energy '/ 
      XE = Z**2 * (0.7053E0/A**(1.E0/3.E0)*(1.E0-0.2E0*Alpha**2) - 1.1529E0/A) 
      EL = XVS + XE 
!     /'   EL = EL + LyPair(Z,A); '/ 
      get_LyMass = EL 

end


!*****************************************
!*             Get LyPair                *
!*****************************************
REAL*4 function get_LyPair(Z,A)

      IMPLICIT NONE
      INTEGER*4 Z
      INTEGER*4 A
!     /' Calculates pairing energy '/ 
!     /' odd-odd nucleus:   get_LyPair = 0 '/ 
!     /' even-odd nucleus:  get_LyPair = -12/sqr(A) '/ 
!     /' even-even nucleus: get_LyPair = -2*12/sqr(A) '/ 
      REAL*4  E_PAIR
!     ' 
      E_PAIR = - 12.E0 / SQRT(REAL(A)) * ( MOD((Z+1) , 2) + MOD((A-Z+1), 2)) 
!     ' 
      get_LyPair = E_PAIR 

end

!*****************************************
!*                TFPair                 *
!***************************************** 
REAL*4 function TFPair(Z,A)

      IMPLICIT NONE
      INTEGER*4 Z
      INTEGER*4 A
!     /' Pairing energy from Thomas-Fermi model of Myers and Swiatecki '/ 
!     /' Shifted that TFPair is zero for odd-odd nuclei '/ 
      INTEGER*4  N
      REAL*4  E_Pair
      N = A - Z 
      IF (   MOD(Z,2)  .EQ. 0 .AND.  MOD(N,2)  .EQ. 0  ) THEN 
!     /' even-even '/ 
      E_Pair = - 4.8E0 / Z**0.333333E0 - 4.8E0 / N**0.333333E0 + 6.6E0 / A**0.666666E0 
      END IF 
      IF (   MOD(Z,2)  .EQ. 0 .AND.  MOD(N,2)  .EQ. 1  ) THEN 
!     /' even Z,odd N '/ 
      E_Pair = - 4.8E0 / Z**0.333333E0 + 6.6E0 / A**0.666666E0 
      END IF 
      IF (   MOD(Z,2)  .EQ. 1 .AND.  MOD(N,2)  .EQ. 0  ) THEN 
!     /' odd Z,even N '/ 
      E_Pair = - 4.8E0 / N**0.333333E0 + 6.6E0 / A**0.666666E0 
      END IF 
      IF (   MOD(Z,2)  .EQ. 1 .AND.  MOD(N,2)  .EQ. 1  ) THEN 
!     /' odd N,odd N '/ 
      E_Pair = 0.0 
      END IF 
      TFPair = E_Pair 

end


!*****************************************
!*                 Pmass                 *
!*****************************************
REAL*4 function Pmass(Z,A,beta)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 beta
!     /' Liquid-drop model of Pearson,2001 '/ 
      REAL*4  N,EA,BE
      REAL*4  avol
      DATA avol/-15.65/
      REAL*4  asf
      DATA asf/17.63/
      REAL*4  r0
      DATA r0/1.233/
      REAL*4  asym
      DATA asym/27.72/
      REAL*4  ass
      DATA ass/-25.60/
      REAL*4  alpha
      REAL*4  pi
      PARAMETER (pi=3.14159)
!     ' 
      N = A - Z 
      alpha = SQRT(5.E0/(4.E0*pi)) * beta 
      EA = avol + asf * A**(-0.333333)*(1.E0+0.4E0*Alpha**2) + &
          0.6E0 * 1.44E0 * Z**2 / (A**1.333333 * r0 )*(1.E0-0.2E0*Alpha**2) + &
          (asym + ass * A**(-0.333333)) * (N-Z)**2 / A**2 
      BE = EA * A 
      Pmass = BE 

end


!*****************************************
!*              Get FEDEFOP              *
!***************************************** 
REAL*4 function get_FEDEFOP(Z,A,beta)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 beta
!     /' According to liquid-drop model of Pearson 2001 '/ 
      REAL*4  asf
      DATA asf/17.63/
      REAL*4  r0
      DATA r0/1.233/
      REAL*4  N,Alpha
      REAL*4  pi
      PARAMETER (pi=3.14159)
!     ' 
      N = A - Z 
      alpha = SQRT(5.E0/(4.E0*pi)) * beta 
      get_FEDEFOP = asf * A**(0.666667)*(0.4E0*Alpha**2) - &
          0.6E0 * 1.44E0 * Z**2 / (A**0.333333 * r0 )*(0.2E0*Alpha**2) 

end


!*****************************************
!*           Get FEDEFOLys               *
!*****************************************
REAL*4 function get_FEDEFOLys(Z,A,beta)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 beta
      REAL*4 LYMASS 
      get_FEDEFOLys = get_LyMass(Z,A,beta) - get_LyMass(Z,A,0.0) 

end


!*****************************************
!*             Get LDMass                *
!*****************************************
REAL*4 function get_LDMass(Z,A,beta)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 beta
      REAL*4  N,BEtab
      REAL*4 LYMASS 
      REAL*4 FEDEFOLYS 
      REAL*4 BEldmTF 
      REAL*4 BEexp 
      N = A - Z 
      BEtab = get_BEldmTF(NINT(N),NINT(Z)) + 2.0 * 12.0 / SQRT(REAL(A)) - 0.00001433*Z**2.39 
!     ' The values in BEtab are the negative binding energies! 
!     ' Pairing in Thomas Fermi masses is zero for Z,N even ! 
      IF (  BEtab .EQ. 0.0  ) THEN 
      BEtab = get_LyMass(Z,A,0.0) 
!     '         Print "Warning: Binding energy of Z=";Z;",A=";A;" not in mass table,";                         " replaced by LYMASS" 
!     '         Print "I_Mode = ";I_Mode 
      End If 
      get_LDMass = BEtab + get_FEDEFOLys(Z,A,beta) 

end


!*****************************************
!*               Get BEexp               *
!*****************************************
! Taken from BEexp.FOR file
REAL*4 function get_BEexp(N,Z)

    IMPLICIT NONE
    SAVE
    INTEGER*4 N
    INTEGER*4 Z
    INTEGER*4 I
    INTEGER*4 J
    INTEGER*4 NI
    INTEGER*4 ZI
    INTEGER*4 AI
    REAL*4 R
    INTEGER*4 IFIRST
    DATA IFIRST /1/
    REAL*4, DIMENSION(0:203,0:136) :: BEexptab
    If (IFIRST.EQ.1) THEN
      print*,'get_BEexp(N,Z), was here  ONCE', BEexptab(N,Z)
      DO I = 0, 203, 1
        DO J = 0, 136, 1
          BEexptab(I,J) = -1.E11
        END DO
      END DO          
      OPEN (UNIT = 1, FILE = 'BEexp.dat', STATUS = 'OLD', ACTION = 'READ')
      DO I = 1, 3352, 1      
         READ (1,*)  ZI, AI, R
         NI = AI - ZI
         BEexptab(NI,ZI) = R
      END DO   
      CLOSE (UNIT = 1)
      IFIRST = 0
    END IF 
    get_BEexp = BEexptab(N,Z)  

end


!*****************************************
!*              Get BELdmFT              *
!*****************************************
! Taken from BELdmFT.FOR file
    REAL*4 function get_BEldmTF(N,Z)
    IMPLICIT NONE
    SAVE
    INTEGER*4 N
    INTEGER*4 Z
    INTEGER*4 I
    INTEGER*4 NI
    INTEGER*4 ZI
    REAL*4 R
    INTEGER*4 IFIRST
    DATA IFIRST /1/
    REAL*4, DIMENSION(0:203,0:136) :: BEldmTFtab
    If (IFIRST.EQ.1) THEN
      print*,'get_BEldmTF(N,Z), was here  ONCE', BEldmTFtab(N,Z)
      OPEN (UNIT = 1, FILE = 'BEldmTF.dat', STATUS = 'OLD', ACTION = 'READ')
      DO I = 1, 8293, 1      
         READ (1,*)  NI, ZI, R
         BEldmTFtab(NI,ZI) = R
      END DO   
      CLOSE (UNIT = 1)
      IFIRST = 0
    END IF 
    get_BEldmTF = BEldmTFtab(N,Z)  

end


!*****************************************
!*            Get AME2012                *
!*****************************************
REAL*4 function get_AME2012(IZ,IA)

      IMPLICIT NONE
      INTEGER*4 IZ
      INTEGER*4 IA
!     ' Masses from the 2003 mass evaluation,complemented by TF masses 
!     ' and Lysekil masses. 
      REAL*4  BEexpval
      REAL*4  Z,A,N
      INTEGER*4  INeu
      REAL*4 LYPAIR 
      REAL*4 U_SHELL 
      REAL*4 LDMASS 
      REAL*4 BEexp 
      INeu = IA - IZ 
      A = REAL(IA) 
      Z = REAL(IZ) 
      N = A - Z 
      BEexpval = get_BEexp(INeu,IZ) 
      IF (  BEexpval .GT. -1.E10  ) THEN 
      get_AME2012 = BEexpval 
      Else 
      get_AME2012 = get_LDMass(Z,A,0.0) + get_U_SHELL(IZ,IA) + get_LyPair(IZ,IA) 
      End If 

end


!*****************************************
!*             Get U Shell               *
!*****************************************
REAL*4 function get_U_SHELL(Z,A)

      IMPLICIT NONE
      INTEGER*4 Z
      INTEGER*4 A
      INTEGER*4  N
      REAL*4  Res
      REAL*4 ShellMO 
      N = A - Z 
      Res = get_ShellMO(N,Z) 
      IF (  Res .GT. 0.0  )  Res = 0.3 * Res 
!     ' KHS (12. Feb. 2012) 
!     '      ' The positive shell effects for deformed nuclei seem to be too positive 
!     ' This gives too many high-energetic prompt neutrons. 
      get_U_SHELL = Res 

end


!*****************************************
!*           Get U Shell Exp             *
!*****************************************
REAL*4 function get_U_SHELL_exp(IZ,IA)

      IMPLICIT NONE
      INTEGER*4 IZ
      INTEGER*4 IA
      REAL*4  Res
      REAL*4  Z,A
      REAL*4 LDMASS 
      REAL*4 LYPAIR 
      REAL*4 AME2012 
      Z = REAL(IZ) 
      A = REAL(IA) 

      Res = 0.5 * ( get_AME2012(IZ,IA) - get_LyPair(IZ,IA) - get_LDMass(Z,A,0.0) ) + &
          0.125 * ( get_AME2012(IZ,IA-1) - get_LyPair(IZ,IA-1) - get_LDMass(Z,A-1.0,0.0) ) + &
          0.125 * ( get_AME2012(IZ,IA+1) - get_LyPair(IZ,IA+1) - get_LDMass(Z,A+1.0,0.0) ) + &
          0.125 * ( get_AME2012(IZ+1,IA+1) - get_LyPair(IZ+1,IA+1) - &
          get_LDMass(Z+1.0,A+1.0,0.0) ) + 0.125 * ( get_AME2012(IZ-1,IA-1) - &
          get_LyPair(IZ-1,IA-1) - get_LDMass(Z-1.0,A-1.0,0.0) ) 
      get_U_SHELL_exp = Res 

end


!*****************************************
!*          Get U Shell EO Exp           *
!*****************************************
REAL*4 function get_U_SHELL_EO_exp(IZ,IA)

      IMPLICIT NONE
      INTEGER*4 IZ
      INTEGER*4 IA
!     ' Returns experimental shell and even-odd staggering 
      REAL*4  Res
      REAL*4  Z,A
      REAL*4 LDMASS 
      REAL*4 LYPAIR 
      REAL*4 AME2012 
      Z = REAL(IZ) 
      A = REAL(IA) 
      Res = get_AME2012(IZ,IA) - get_LDMass(Z,A,0.0) 
      get_U_SHELL_EO_exp = Res 

end


!*****************************************
!*             Get U MASS                *
!*****************************************
REAL*4 function get_U_MASS(Z,A)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
!     /' LD + congruence energy + shell (no pairing) '/ 
      REAL*4  BE
      REAL*4 U_SHELL 
      REAL*4 LDMASS 
!      IF (  Z .LT. 0 .OR. A .LT. 0  ) THEN 
!     '       Print "U_Mass: Z,A",Z,A 
!      End If 
      BE = get_LDMass(Z,A,0.0) + get_U_SHELL(NINT(Z),NINT(A))  
      get_U_MASS = BE 

end


!*****************************************
!*              Get ECOUL                *
!***************************************** 
REAL*4 function get_ECOUL(Z1,A1,beta1,Z2,A2,beta2,d)

      IMPLICIT NONE
      REAL*4 Z1
      REAL*4 A1
      REAL*4 beta1
      REAL*4 Z2
      REAL*4 A2
      REAL*4 beta2
      REAL*4 d
!     ' 
!     /' Coulomb potential between two nuclei                    '/ 
!     /' surfaces are in a distance of d                         '/ 
!     /' in a tip to tip configuration                           '/ 
!     ' 
!     /' approximate formulation                                 '/ 
!     /' On input: Z1      nuclear charge of first nucleus       '/ 
!     /'           A1      mass number of irst nucleus   '/ 
!     /'           beta1   deformation of first nucleus          '/ 
!     /'           Z2      nuclear charge of second nucleus      '/ 
!     /'           A2      mass number of second nucleus  '/ 
!     /'           beta2   deformation of second nucleus         '/ 
!     /'           d       distance of surfaces of the nuclei    '/ 
!     ' 
      REAL*4  N1,N2,recoul
      REAL*4  dtot
      REAL*4  r0
      DATA r0/1.16/
!     ' 
      N1 = A1 - Z1 
      N2 = A2 - Z2 
      dtot = r0 *( (Z1+N1)**0.3333333E0 * (1.E0+0.6666667E0*beta1) + &
          (Z2+N2)**0.3333333E0 * (1.E0+0.6666667E0*beta2) ) + d 
      REcoul = Z1 * Z2 * 1.44E0 / dtot 
!     ' 
      get_ECOUL = REcoul 

end


!*****************************************
!*              Beta Light               *
!*****************************************
REAL*4 function get_beta_light(Z,betaL0,betaL1)

      IMPLICIT NONE
      INTEGER*4 Z
      REAL*4 betaL0
      REAL*4 betaL1
!     /' Deformation of light fission fragment for S1 and S2 '/ 
!     /' Systematic correlation Z vs. beta for deformed shells '/ 
!     /' Z of fission fragment '/ 
      REAL*4  beta
      beta = (Z - betaL0) * betaL1/20.E0 
      get_beta_light = beta 

end


!*****************************************
!*              Beta Heavy               *
!*****************************************
REAL*4 function get_beta_heavy(Z,betaH0,betaH1)

      IMPLICIT NONE
      INTEGER*4 Z
      REAL*4 betaH0
      REAL*4 betaH1
!     /' Deformation of heavy fission fragment for S2 '/ 
!     /' Systematic correlation Z vs. beta for deformed shells '/ 
!     /' Z of fission fragment '/ 
      REAL*4  beta
      beta = (Z - betaH0) * betaH1/20.E0 
      get_beta_heavy = beta 

end


!*****************************************
!*                Z Equi                 *
!***************************************** 
REAL*4 function get_Z_equi(ZCN,A1,A2,beta1,beta2,d,Imode,POLARadd,POLARfac)

      IMPLICIT NONE
      INTEGER*4 ZCN
      INTEGER*4 A1
      INTEGER*4 A2
      REAL*4 beta1
      REAL*4 beta2
      REAL*4 d
      INTEGER*4 Imode
      REAL*4 POLARadd
      REAL*4 POLARfac

!     /' Determines the minimum potential of the scission-point configuration 
!     represented by two deformed nuclei divided by a tip distance d. 
!     A1,A2,beta1,beta2,d are fixed,Z1 is searched for and returned on output.  '/ 
!     ' 
!     /' ZCN: Z of fissioning nucleus '/ 
!     /' A1: A of first fission fragment '/ 
!     /' A2: A of second fission fragment '/ 
!     /' beta1: deformation of first fission fragment '/ 
!     /' beta2: deformation of second fission fragment '/ 
!     /' d: tip distance '/ 
!     ' 
      REAL*4  RZ_equi
      REAL*4  RA1,RA2,RZCN,RACN
      REAL*4  Z1UCD,Z2UCD
      REAL*4  re1,re2,re3,eps1,eps2,DZ_Pol
!     /' help variables '/ 
      REAL*4 ECOUL 
      REAL*4 LYMASS 
!     ' 
      RA1 = REAL(A1) 
      RA2 = REAL(A2) 
      RZCN = REAL(ZCN) 
      RACN = RA1 + RA2 
      Z1UCD = RA1 / (RA1 + RA2) * RZCN 
      Z2UCD = RZCN - Z1UCD 
      re1 = get_LyMass( Z1UCD-1.E0,RA1,beta1 ) + &
          get_LyMass( Z2UCD+1.E0,RA2,beta2 ) + &
          get_ECOUL( Z1UCD-1.E0,RA1,beta1,Z2UCD+1.E0,RA2,beta2,d ) 
      re2 = get_LyMass( Z1UCD,RA1,beta1) + get_LyMass( Z2UCD,RA2,beta2) + &
          get_ECOUL( Z1UCD,RA1,beta1,Z2UCD,RA2,beta2,d ) 
      re3 = get_LyMass( Z1UCD+1.E0,RA1,beta1 ) + &
            get_LyMass( Z2UCD-1.E0,RA2,beta2 ) + &
            get_ECOUL( Z1UCD+1.E0,RA1,beta1,Z2UCD-1.E0,RA2,beta2,d ) 
      eps2 = ( re1 - 2.E0*re2 + re3 ) / 2.E0 
      eps1 = ( re3 - re1 ) / 2.E0 
      DZ_Pol = -eps1 / ( 2.E0 * eps2 ) 
!     ' 
      IF (  DZ_Pol .GT. 2 .OR. DZ_Pol .LT. -2  )  DZ_Pol = 0 
!     ' 
      IF (  Imode .GT. 0  ) THEN 
!     /' Purely empirical enhancement of charge polarization '/ 
      DZ_POL = DZ_POL * POLARfac + POLARadd 
      End If 
!     ' 
      RZ_equi = Z1UCD + DZ_POL 
      get_Z_equi = RZ_equi 

end


!*****************************************
!*            Beta Opt Light             *
!*****************************************
subroutine Beta_opt_light(A1,A2,Z1,Z2,d,beta2_imposed,beta1_opt)

      IMPLICIT NONE
      REAL*4 A1
      REAL*4 A2
      REAL*4 Z1
      REAL*4 Z2
      REAL*4 d
      REAL*4 beta2_imposed
      REAL*4 beta1_opt
!     /' Determines the optimum deformation of the light fragment when the deformation of the 
!     heavy fragment is imposed. '/ 
!     ' 
      REAL*4  beta1,dbeta1,beta1_prev,beta1_next
      REAL*4  Uguess,Uplus,Uminus,Uprev,Unext
      INTEGER*4  I
      REAL*4 ECOUL 
      REAL*4 LYMASS 
!     ' 
!     /' List('Beta_opt_light called with '); 
!     List(A1,A2,Z1,Z2,d,beta2_imposed,beta1_opt); 
!     DCL Byes Bit(1) aligned; 
!     Call GPYES('Continue',Byes); '/ 
      beta1 = 0.5 
      dbeta1 = 0.01 
      Uguess = get_LyMass(Z1,A1,beta1) + get_LyMass(Z2,A2,beta2_imposed) + &
          get_ECOUL(Z1,A1,beta1,Z2,A2,beta2_imposed,d) 
      Uplus = get_LyMass(Z1,A1,beta1 + dbeta1) + get_LyMass(Z2,A2,beta2_imposed) + &
          get_ECOUL(Z1,A1,beta1 + dbeta1,Z2,A2,beta2_imposed,d) 
      Uminus = get_LyMass(Z1,A1,beta1 - dbeta1) + get_LyMass(Z2,A2,beta2_imposed) + &
          get_ECOUL(Z1,A1,beta1 - dbeta1,Z2,A2,beta2_imposed,d) 
      IF (  Uplus .GT. Uguess .AND. Uminus .GT. Uguess  ) THEN 
      beta1_opt = beta1 
      Else 
      IF (  Uplus .LT. Uguess  )  dbeta1 = 0.01 
      IF (  Uminus .LT. Uguess  )  dbeta1 = -0.01 
      Unext = Uguess 
      beta1_next = beta1 
      DO I = 1 , 10000 
      beta1_prev = beta1_next 
      Uprev = Unext 
      beta1_next = beta1_prev + dbeta1 
      Unext = get_LyMass(Z1,A1,beta1_next) + get_LyMass(Z2,A2,beta2_imposed) + &
          get_ECOUL(Z1,A1,beta1_next,Z2,A2,beta2_imposed,d) 
      IF (  Unext .GE. Uprev  )  Exit 
      END DO 
      beta1_opt = beta1_prev 
      END IF 

end
      

!*****************************************
!*              Beta Equi                *
!*****************************************
subroutine Beta_Equi(A1,A2,Z1,Z2,d,beta1prev,beta2prev,beta1opt,beta2opt)

      IMPLICIT NONE
      REAL*4 A1
      REAL*4 A2
      REAL*4 Z1
      REAL*4 Z2
      REAL*4 d
      REAL*4 beta1prev
      REAL*4 beta2prev
      REAL*4 beta1opt
      REAL*4 beta2opt
!     /' Determines the minimum potential of the scission-point configuration 
!     represented by two deformed nuclei,divided by a tip distance d. 
!     A1,A2,Z1,Z2,d are fixed,beta1 and beta2 are searched for and returned on output '/ 
!     ' 
      INTEGER*4  B_analytical
      DATA B_analytical/0/
!     ' Switch to use the analytical approximation 
!     ' that replaces the long numerical calculation. 
      REAL*4  x,y,xcoul
      REAL*4  xcoul236U
      DATA xcoul236U/1369.64/
!     ' 
      REAL*4  beta1,beta2
!     ' 
!     '      Dim As Double U,Uprev,Ulast,Ubest,Uopt 
      REAL*4  U,Uprev,Ulast,Ubest,Uopt
!     ' 
!     '      Dim As Double sbeta1,sbeta2 
      REAL*4  sbeta1,sbeta2
!     ' 
      INTEGER*4  N,N1,N2,Nopt
!     ' 
!     '      Dim As Double eps = 5.E-4 
      REAL*4  eps
      DATA eps/5.E-4/
!     ' 
      INTEGER*4  I
      REAL*4 LYMASS 
      REAL*4 ECOUL 
!     ' 
      IF (  B_analytical .EQ. 0  ) THEN 
!     ' Numerical algorithm 
!     ' 
      beta1 = beta1prev 
      beta2 = beta2prev 
      Uprev = get_LyMass(Z1,A1,beta1) + get_LyMass(Z2,A2,beta2) + &
          get_ECOUL(Z1,A1,beta1,Z2,A2,beta2,d) 
      Uopt = Uprev 
!     ' 
!     /' Test slope of variation of U '/ 
      beta1 = beta1prev + eps 
      U = 1.E30 
!     ' 
      beta2 = beta2prev 
!     '     For beta2 = beta2prev to 0 Step -eps 
      DO I = 1 , Int(beta2prev/eps) 
      beta2 = beta2 - eps 
      Ulast = U 
      U = get_LyMass(Z1,A1,beta1) + get_LyMass(Z2,A2,beta2) + &
          get_ECOUL(Z1,A1,beta1,Z2,A2,beta2,d) 
      IF (  U .GT. Ulast  ) THEN 
      Exit 
      Else 
      Ubest = U 
      END IF 
      END DO 
      IF (  Ubest .LT. Uopt  ) THEN 
      Uopt = Ubest 
      sbeta1 = eps 
      sbeta2 = -eps 
      END IF 
!     ' 
      U = 1.E30 
      beta2 = beta2prev 
!     '   For beta2 = beta2prev To 1 Step eps 
      DO I = 1 , Int((1 - beta2prev)/eps) 
      beta2 = beta2 + eps 
      Ulast = U 
      U = get_LyMass(Z1,A1,beta1) + get_LyMass(Z2,A2,beta2) + &
          get_ECOUL(Z1,A1,beta1,Z2,A2,beta2,d) 
      IF (  U .GT. Ulast  ) THEN 
      Exit 
      Else 
      Ubest = U 
      END IF 
      END DO 
      IF (  Ubest .LT. Uopt  ) THEN 
      Uopt = Ubest 
      sbeta1 = eps 
      sbeta2 = eps 
      End If 
!     ' 
      beta1 = beta1prev - eps 
      U = 1.E30 
      beta2 = beta2prev 
!     '   For beta2 = beta2prev To 0 Step -eps 
      DO I = 1 , Int(beta2prev/eps) 
      beta2 = beta2 - eps 
      Ulast = U 
      U = get_LyMass(Z1,A1,beta1) + get_LyMass(Z2,A2,beta2) + &
          get_ECOUL(Z1,A1,beta1,Z2,A2,beta2,d) 
      IF (  U .GT. Ulast  ) THEN 
      Exit 
      Else 
      Ubest = U 
      End If 
      END DO 
      IF (  Ubest .LT. Uopt  ) THEN 
      Uopt = Ubest 
      sbeta1 = -eps 
      sbeta2 = -eps 
      END IF 
!     ' 
      U = 1.E30 
      beta2 = beta2prev 
!     '   For beta2 = beta2prev To 1 Step eps 
      DO I = 1 , Int((1-beta2prev)/eps) 
      beta2 = beta2 + eps 
      Ulast = U 
      U = get_LyMass(Z1,A1,beta1) + get_LyMass(Z2,A2,beta2) + &
          get_ECOUL(Z1,A1,beta1,Z2,A2,beta2,d) 
      IF (  U .GT. Ulast  ) THEN 
      Exit 
      Else 
      Ubest = U 
      END IF 
      END DO 
      IF (  Ubest .LT. Uopt  ) THEN 
      Uopt = Ubest 
      sbeta1 = -eps 
      sbeta2 = eps 
      END IF 
!     ' 
!     ' 
      Ubest = get_LyMass(Z1,A1,beta1prev) + get_LyMass(Z2,A2,beta2prev) + &
          get_ECOUL(Z1,A1,beta1prev,Z2,A2,beta2prev,d) 
      U = get_LyMass(Z1,A1,beta1prev+REAL(sbeta1)) + &
          get_LyMass(Z2,A2,beta2prev+REAL(sbeta2)) + &
          get_ECOUL(Z1,A1,beta1prev+sbeta1,Z2,A2,beta2prev+REAL(sbeta2),d) 
!     ' 
!     '   L1: 
      DO N = 1 , 1000 
!     ' 
!     '   L2: 
      DO N1 = 1 , N 
      N2 = N-N1 
      beta1 = beta1prev + sbeta1*N1 
      beta2 = beta2prev + sbeta2*N2 
      U = get_LyMass(Z1,A1,beta1) + get_LyMass(Z2,A2,beta2) + &
          get_ECOUL(Z1,A1,beta1,Z2,A2,beta2,d) 
      IF (  U .LT. Ubest  ) THEN 
      Ubest = U 
      beta1opt = beta1 
      beta2opt = beta2 
      Nopt = N 
      END IF 
      END DO 
      IF (  N-Nopt .GT. 2  )  Exit 
      END DO 
!     ' 
!     ' 
      Else 
!     ' Analytical approximation 
!     ' Must be adapted if the relevant parameters of GEF are modified! 
      xcoul = (Z1 + Z2)**2 / (A1 + A2)**(1.0/3.0) 
      x = (Z1 / (Z1 + Z2))**(xcoul/xcoul236U) 
      y = 1.2512E-4 + 0.00122851*x - 0.00267707*x**2 + 0.00372901*x**3 - 0.00219903*x**4 
      beta1opt = y * xcoul 
!     ' 
      x = (Z2 / (Z1 + Z2))**(xcoul/xcoul236U) 
      y = 1.2512E-4 + 0.00122851*x - 0.00267707*x**2 + 0.00372901*x**3 - 0.00219903*x**4 
      beta2opt = y * xcoul 
!     ' 
      End If 
!     ' 

end


!*****************************************
!*             Get U Ired                *
!***************************************** 
REAL*4 function get_U_Ired(Z,A)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
!     ' Effective moment of inertia by pairing with correction for excitation energy 
      REAL*4  I_rigid_spher,IfragEff
!     ' 
      REAL*4 U_SHELL 
!     ' 
      I_rigid_spher = 1.16E0**2 * A**1.6667E0 / 103.8415E0 
!     '   IfragEff = I_rigid_spher + 0.003 * A^(4.0/3.0) * U_shell(Cint(Z),Cint(A)) 
!     '   IfragEff = I_rigid_spher + 0.005 * A^(4.0/3.0) * U_shell(Cint(Z),Cint(A)) 
!     ' reduction due to shell (Deleplanque et al. PRC 69 (2004) 044309) 
      IfragEff = 0.45 * I_rigid_spher 
!     ' Effect of superfluidity 
!     '   IfragEff = 0.65 * IfragEff   ' Average effect of superfluidity and deformation 
!     ' 
      get_U_Ired = IfragEff 

end


!*****************************************
!*                U IredFF               *
!***************************************** 
REAL*4 function U_IredFF(Z,A)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
!     ' Effective moment of inertia by pairing with correction for excitation energy 
!     ' of final fission fragments 
!     ' 
      REAL*4 U_Ired 
      REAL*4 U_I_Shell 
!     ' 
      U_IredFF = get_U_Ired(Z,A) * get_U_I_Shell(Z,A) 

end


!*****************************************
!*            Get U I Shell              *
!*****************************************
REAL*4 function get_U_I_Shell(Z,A)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      INTEGER*4  N_shells(6)
!     ' Shell effect on the effective moment of inertia 
      INTEGER*4  I
      REAL*4  dNmin,dZmin,dNsubmin
      REAL*4  Inv_add
      DATA Inv_add/0/
      REAL*4  I_inv_add_Z
      DATA I_inv_add_Z/0/
      REAL*4  I_inv_add_N
      DATA I_inv_add_N/0/
      REAL*4  I_inv_add_Nsub
      DATA I_inv_add_Nsub/0/
      N_shells(1) = 20 
      N_shells(2) = 28 
      N_shells(3) = 50 
      N_shells(4) = 82 
      N_shells(5) = 126 
      N_shells(6) = 56 
      dNmin = 100 
      dZmin = 100 
      dNsubmin = 100 
      DO I = 1 , 5 
      dZmin = Min(dZmin,Abs(N_shells(I) - Z)) 
      END DO 
!     ' 
      DO I = 1 , 5 
      dNmin = Min(dNmin,Abs(N_shells(I) - (A-Z))) 
      END DO 
!     ' 
      dNsubmin = Abs(N_shells(6) - (A-Z)) 
!     ' 
!     ' Effect of shells: 
      IF (  dZmin .LT. 10.0  ) THEN 
!     '        I_inv_add_Z = 0.33 * (6.0 * sqr(A/140.) - dZmin) * sqr(140./A) 
      I_inv_add_Z = 0.33 * (6.0 * SQRT(A/140.0) - dZmin) * (140.0/A)**1.5 
!     ' A^(-1/3) dependence: "A simple phenomenology for 2gamma+ states", 
!     ' N. V. Zamfir,D. Bucurescu,R. F. Casten,M. Ivascu, 
!     ' Phys. Lett. B 241 (1990) 463 
      I_inv_add_Z = Max(I_inv_add_Z,0.0) 
      End If 
      IF (  dNmin .LT. 10.0  ) THEN 
!     '        I_inv_add_N = 0.42 * (8.0 * sqr(A/140.) - dNmin) * sqr(140./A) 
      I_inv_add_N = 0.42 * (8.0 * SQRT(A/140.0) - dNmin) * (140.0/A)**1.5 
      I_inv_add_N = Max(I_inv_add_N,0.0) 
      End If 
      IF (  DNsubmin .LT. 6.0  ) THEN 
!     '    I_inv_add_Nsub = 1.7 * (4.0 - dNsubmin) * (1.0 - 0.32 * Abs(40.0-Z)) 
      I_inv_add_Nsub = 1.7 * (4.0 - dNsubmin) * (1.0 - 0.18 * Abs(40.0-Z)) 
!     ' N = 56 subshell only around Z = 40 
      I_inv_add_Nsub = Max(I_inv_add_Nsub,0.0) 
      End If 
      get_U_I_Shell = 1.0 / (1.0 + Max(I_inv_add_N,I_inv_add_Nsub) + I_inv_add_Z) 

end


!*****************************************
!*            Get U Alev Ld              *
!*****************************************
REAL*4 function get_U_alev_ld(Z,A)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
!     '  get_U_alev_ld = 0.073 * A + 0.095 * A^0.666667  'Ignatyuk (1970's) 
      get_U_alev_ld = 0.078 * A + 0.115 * A**0.6666667 
!     ' Ignatyuk (Bologna 2000) 
!     '  get_U_alev_ld = 0.089 * A    ' only volume term 

end


!*****************************************
!*              Get U Temp               *
!*****************************************
REAL*4 function get_U_Temp(Z,A,E,Ishell,Ipair,Tscale,Econd)

      IMPLICIT NONE
      REAL*4 Z
      REAL*4 A
      REAL*4 E
      INTEGER*4 Ishell
      INTEGER*4 Ipair
      REAL*4 Tscale
      REAL*4 Econd
!     ' Temperature (modified Gilbert-Cameron composite level density) 
!     ' KHS (10. 2. 2012) 
      REAL*4  alev
      REAL*4  Eeff0,Eeff1,Rho0,Rho1,TCT,TFG
      REAL*4  fgamma
      DATA fgamma/0.055/
      REAL*4  RShell,RPair,Res
      REAL*4 U_ALEV_LD 
      REAL*4 U_SHELL 
      REAL*4 LYPAIR 
      REAL*4 TEGIDY 
!     ' Used global parameters: Tscale 
!     '   alev = get_U_alev_ld(Z,A) * 1.1   ' Factor adjusted to high-energy prompt neutrons in U235(nth,f) 
!     '  alev = get_U_alev_ld(Z,A) * 0.8  ' " with the correction for non-constant T (FG range) 
      alev = get_U_alev_ld(Z,A) 
!     ' 
      IF (  Ishell .EQ. 1  ) THEN 
      RShell = get_U_SHELL(NINT(Z),NINT(A)) 
      Else 
      RShell = 0.0 
      End If 
      TCT = get_TEgidy(A,RShell,Tscale) 
!     ' 
      IF (  Ipair .EQ. 1  ) THEN 
      RPair = get_LyPair(NINT(Z),NINT(A)) 
      Else 
      Rpair = 0.0 
      End If 
      Eeff0 = E - Econd + RPair + Rshell*(1.0 - exp(-fgamma * E)) 
!     ' 
      IF (  Eeff0 .GT. 0.5  ) THEN 
      Eeff1 = Eeff0 + 0.1 
      Rho0 = 1.E0/Eeff0**1.25 * exp(2.E0 * SQRT(alev * Eeff0)) 
      Rho1 = 1.E0/Eeff1**1.25 * exp(2.E0 * SQRT(alev * Eeff1)) 
!     '         Rho0 = 1.E0/Eeff0 * exp(2.E0 * sqr(alev * Eeff0)) 
!     '         Rho1 = 1.E0/Eeff1 * exp(2.E0 * sqr(alev * Eeff1)) 
      TFG = 0.1E0 / (log(Rho1) - log(Rho0)) 
      Else 
      TFG = 0.0 
      End If 
      Res = TCT 
      IF (  TFG .GT. Res  )  Res = TFG 
!     ' 
!     ' If Res > 1.4 Then Res = 1.4 
!     ' 
      get_U_Temp = Res 

end


!*****************************************
!*           Get U Even Odd              *
!*****************************************
REAL*4 function get_U_Even_Odd(I_Channel,PEO)

      IMPLICIT NONE
      INTEGER*4 I_Channel
      REAL*4 PEO
!     ' Creates even-odd fluctuations 
      REAL*4  R
      IF (   MOD(I_Channel,2)  .EQ. 0  ) THEN 
      R = 1.0 + PEO 
      Else 
      R = 1.0 - PEO 
      End If 
      get_U_Even_Odd = R 

end


!*****************************************
!*              Get BFTF                 *
!***************************************** 
REAL*4 function get_BFTF(RZ,RA,I_Switch)

      IMPLICIT NONE
      REAL*4 RZ
      REAL*4 RA
      INTEGER*4 I_Switch
!     /' Fission barriers from Myers and Swiatecki,Thomas-Fermi model '/ 
!     /'  I_Switch: 0: liquid-drop; 1: with shells and pairing, 
!     2: averaged over pairing,3: with shell and pairing + pairing gap at barrier '/ 
!     ' 4: liquid-drop + g.s. shell,no Z correction 
      REAL*4  RN,RI,Rkappa,RS,RF,RX
      REAL*4  RX0
      DATA RX0/48.5428/
      REAL*4  RX1
      DATA RX1/34.15/
      REAL*4  RB
      INTEGER*4  IZ,IA
      REAL*4 U_SHELL 
      REAL*4 U_SHELL_EXP 
      REAL*4 U_SHELL_EO_EXP 
      REAL*4 LYPAIR 
!     ' 
      IZ = NINT(RZ) 
      IA = NINT(RA) 
      RN = RA - RZ 
      RI = (RN-RZ) / RA 
      Rkappa = 1.9E0 + (RZ - 80.E0) / 75.E0 
      RS = RA**0.666667E0 * (1.E0 - Rkappa * RI**2) 
      RX = RZ**2 / (RA * (1.E0 - Rkappa * RI**2)) 
      IF (  RX .LT. 30  ) THEN 
!     /' out of range '/ 
      RF = 1.E10 
      End If 
      IF (  RX .GT. RX0  ) THEN 
!     /' out of range '/ 
      RF = 0.0 
      End If 
      IF (  RX .LT. RX1 .AND. RX .GT. 30  ) THEN 
      RF = 0.595553E0 - 0.124136E0 * (RX - RX1) 
      End If 
      IF (  RX .GE. RX1 .AND. RX .LE. RX0  ) THEN 
      RF = 0.000199749 * (RX0 - RX)**3 
      End If 
      RB = RF * RS 
!     ' 
      Select CASE( I_Switch) 
      CASE( 0) 
      get_BFTF = RB 
      CASE( 1) 
!     ' including even-odd staggering due to increased pairing strength at barrier 
!     ' Tentative modification from comparison with experimental fission barriers 
!     ' (shell correction at the barrier?) 
      IF (  RZ .GT. 86.5  )  RB = RB - 0.15 * (RZ - 86.5) 
!     '    If RZ > 90 Then RB = RB + 0.3 * (RZ - 90.0) 
!     '    If RZ > 98 Then RB = RB - 0.15 * (RZ - 98.0) 
      IF (  RZ .GT. 90  )  RB = RB + 0.35 * (RZ - 90.0) 
      IF (  RZ .GT. 93  )  RB = RB + 0.15 * (RZ - 93.0) 
      IF (  RZ .GT. 95  )  RB = RB - 0.25 * (RZ - 95.0) 
!     '    get_BFTF = RB - get_U_SHELL(IZ,IA) 
!     '    get_BFTF = RB - get_U_SHELL_exp(IZ,IA) 
      get_BFTF = RB - get_U_SHELL_EO_exp(IZ,IA) + get_LyPair(IZ,IA) * 14.0/12.0 
      CASE( 2) 
!     ' averaged over even-odd staggering 
      IF (  RZ .GT. 86.5  )  RB = RB - 0.15 * (RZ - 86.5) 
      IF (  RZ .GT. 90  )  RB = RB + 0.35 * (RZ - 90.0) 
      IF (  RZ .GT. 93  )  RB = RB + 0.15 * (RZ - 93.0) 
      IF (  RZ .GT. 95  )  RB = RB - 0.25 * (RZ - 95.0) 
      get_BFTF = RB - get_U_SHELL_exp(IZ,IA) 
      CASE( 3) 
!     ' like Case 1 + pairing gap at barrier 
      IF (  RZ .GT. 86.5  )  RB = RB - 0.15 * (RZ - 86.5) 
      IF (  RZ .GT. 90  )  RB = RB + 0.35 * (RZ - 90.0) 
      IF (  RZ .GT. 93  )  RB = RB + 0.15 * (RZ - 93.0) 
      IF (  RZ .GT. 95  )  RB = RB - 0.25 * (RZ - 95.0) 
      get_BFTF = RB - get_U_SHELL_EO_exp(IZ,IA) 
      CASE( 4) 
!     ' like case 3 but without Z correction 
!     ' This is the direct description from the topographic theorem. 
      get_BFTF = RB - get_U_SHELL_exp(IZ,IA) 
      Case DEFAULT 
!     '         Print "Undefined option in BFTF" 
!     '         Sleep 
      End Select 

end


!*****************************************
!*              Get BFTFA                *
!*****************************************
REAL*4 function get_BFTFA(RZ,RA,I_Switch)

      IMPLICIT NONE
      REAL*4 RZ
      REAL*4 RA
      INTEGER*4 I_Switch
!     /' inner barrier height '/ 
      REAL*4  EA,BF0,Z4A,Z3A,DB
      REAL*4  coeff
      DATA coeff/0.5/
      REAL*4 BFTF 
      BF0 = get_BFTF(RZ,RA,I_Switch) 
!     ' Z4A = RZ^4 / RA 
!     '  EB - EA from fit to Smirenkin barriers: 
!     '  V. M. Kupriyanov,K. K. Istekov,B. I. Fursov,G. N. Smirenkin 
!     '  Sov. J. Nucl. Phys. 32 (1980) 184 
!     '  DB = -10.3517 + 1.6027E-5 * Z4A + 5.4945E-11 * Z4A^2  ' EA - EB 
!     ' 
!     '  EB - EA from fit to data from Dahlinger et al. (KHS,21. Dec. 2012) 
      Z3A = RZ**3 / RA 
      DB = -(5.40101 - 0.00666175*Z3A + 1.52531E-6*Z3A**2) 
      IF (  DB .GT. 0.0  ) THEN 
      EA = BF0 - DB 
      Else 
      EA = BF0 
      End If 
      get_BFTFA = EA 

end


!*****************************************
!*             Get BFTFB                 *
!*****************************************
REAL*4 function get_BFTFB(RZ,RA,I_Switch)

      IMPLICIT NONE
      REAL*4 RZ
      REAL*4 RA
      INTEGER*4 I_Switch
!     /' outer barrier height '/ 
      REAL*4  EB,BF0,Z4A,Z3A,DB
      REAL*4  coeff
      DATA coeff/0.5/
      REAL*4 BFTF 
      BF0 = get_BFTF(RZ,RA,I_Switch) 
!     ' Z4A = RZ^4 / RA 
!     '  EB - EA from fit to Smirenkin barriers: 
!     '  V. M. Kupriyanov,K. K. Istekov,B. I. Fursov,G. N. Smirenkin 
!     '  Sov. J. Nucl. Phys. 32 (1980) 184 
!     '   DB = -10.3517 + 1.6027E-5 * Z4A + 5.4945E-11 * Z4A^2  ' EA - EB 
!     ' 
!     '  EB - EA from fit to data from Dahlinger et al. (KHS,21. Dec. 2012) 
      Z3A = RZ**3 / RA 
      DB = -(5.40101 - 0.00666175*Z3A + 1.52531E-6*Z3A**2) 
      IF (  DB .LT. 0.0  ) THEN 
      EB = BF0 + DB 
      Else 
      EB = BF0 
      End If 
      get_BFTFB = EB 

end


!*****************************************
!*           Get Gaussintegral           *
!*****************************************
REAL*4 function get_Gaussintegral(R_x,R_sigma)

      IMPLICIT NONE
      REAL*4 R_x
      REAL*4 R_sigma
!     /' Smoothed step function. Grows from 0 to 1 around R_x 
!     with a Gauss-integral function with given sigma'/ 
      REAL*4  R_ret
!     ' Note: The variable R_sigma = standard deviation / sqr(2) ! 
      REAL*4 ERF 
      R_ret = 0.5E0 + 0.5E0 * Erf(R_x / R_sigma) 
      get_Gaussintegral = R_ret 

end


!*****************************************
!*              Get U Box                *
!*****************************************
REAL*4 function get_U_Box(x,sigma,length)

      IMPLICIT NONE
      REAL*4 x
      REAL*4 sigma
      REAL*4 length
      REAL*4  y
!     ' Note: The variable sigma = standard deviation / sqr(2) ! 
      REAL*4 GAUSSINTEGRAL 
      y = get_Gaussintegral(x+0.5*length,sigma) - get_Gaussintegral(x-0.5*length,sigma) 
      get_U_Box = y/length 

end


!*****************************************
!*              Get U Box2               *
!*****************************************
REAL*4 function get_U_Box2(x,sigma1,sigma2,length)

      IMPLICIT NONE
      REAL*4 x
      REAL*4 sigma1
      REAL*4 sigma2
      REAL*4 length
      REAL*4  y
!     ' Note: The variable sigma = standard deviation / sqr(2) ! 
      REAL*4 GAUSSINTEGRAL 
      y = get_Gaussintegral(x+0.5*length,sigma2) - get_Gaussintegral(x-0.5*length,sigma1) 
      get_U_Box2 = y/length 

end


!*****************************************
!*               U Gauss                 *
!*****************************************
REAL*4 function U_Gauss(x,sigma)

      IMPLICIT NONE
      REAL*4 x
      REAL*4 sigma
      REAL*4  y
      REAL*4  pi
      PARAMETER (pi=3.14159)
!     ' 
      y = 1.0 / (SQRT(2.0 * pi) * sigma) * exp(-x**2/ ( 2.0 * sigma**2 ) ) 
      U_Gauss = y 

end


!*****************************************
!*            Get U Gauss Mod            *
!*****************************************
REAL*4 function get_U_Gauss_mod(x,sigma)

      IMPLICIT NONE
      REAL*4 x
      REAL*4 sigma
!     ' Gaussian with Sheppard correction 
      REAL*4  y
      REAL*4  sigma_mod
      REAL*4  pi
      PARAMETER (pi=3.14159)
      sigma_mod = SQRT(sigma**2 + 1./12.) 
!     ' 
      y = 1.0 / (SQRT(2.0 * pi) * sigma_mod) * exp(-x**2/ ( 2.0 * sigma_mod**2 ) ) 
      get_U_Gauss_mod = y 

end


!*****************************************
!*            Get U LinGauss             *
!*****************************************
REAL*4 function get_U_LinGauss(x,R_Sigma)

      IMPLICIT NONE
      REAL*4 x
      REAL*4 R_Sigma
!     /' Gaussian times a linear function '/ 
!     /' Not normalized! '/ 
      REAL*4  R_Res
      IF (  R_Sigma .GT. 0.0  ) THEN 
      R_Res = x * exp(-x**2/(2.0 * R_Sigma**2)) 
      Else 
      R_Res = 0.0 
      End If 
      get_U_LinGauss = R_Res 

end


!*****************************************
!*             Get ShellMO               *
!*****************************************
! Taken from ShellMO.FOR
REAL*4 function get_ShellMO(N,Z)

      IMPLICIT NONE
      SAVE
      INTEGER*4 N
      INTEGER*4 Z
      INTEGER*4 I
      INTEGER*4 NI
      INTEGER*4 ZI
      REAL*4 R
      INTEGER*4 IFIRST
      DATA IFIRST /1/
      REAL*4, DIMENSION(0:203,0:136) :: ShellMOtab
      If (IFIRST.EQ.1) THEN
        print*,'get_ShellMO(N,Z), was here  ONCE', SHellMOtab(N,Z)
        OPEN (UNIT = 1, FILE = 'BEldmTF.dat', STATUS = 'OLD', ACTION = 'READ')
        DO I = 1, 8270, 1      
           READ (1,*)  NI, ZI, R
           ShellMOtab(NI,ZI) = R
        END DO   
        CLOSE (UNIT = 1)
        IFIRST = 0
      END IF 
      get_ShellMO = SHellMOtab(N,Z)  

end


end module GEFSUB_FOR
!===========================================


!*******************************************
!*  Declarataion of the GEFSUBdcl2 Module  *
!*******************************************
! Taken from GEFSUBdcl2.FOR
module GEFSUBdcl2

!     Output of FBtoFO from GEFSUB.BAS
      REAL*4 Getyield
      REAL*4 Masscurv
      REAL*4 Masscurv1
!     ' 
      REAL*4 De_Saddle_Scission
!     ' 
      REAL*4 TEgidy
!     ' 
      REAL*4 TRusanov
!     ' 
      REAL*4 LyMass
!     ' 
      REAL*4 LyPair
!     ' 
      REAL*4 TFPair
!     ' 
      REAL*4 Pmass
!     ' 
      REAL*4 FEDEFOLys
!     ' 
      REAL*4 FEDEFOP
!     ' 
      REAL*4 LDMass
!     ' 
      REAL*4 AME2012
!     ' 
      REAL*4 U_SHELL
!     ' 
      REAL*4 U_SHELL_exp
!     ' 
      REAL*4 U_SHELL_EO_exp
!     ' 
      REAL*4 U_MASS
!     ' 
      REAL*4 ECOUL
!     ' 
      REAL*4 beta_light
!     ' 
      REAL*4 beta_heavy
!     ' 
      REAL*4 Z_equi
!     ' 
!     ' 
!     ' 
      REAL*4 U_Ired
!     ' 
      REAL*4 U_IredFF
!     ' 
      REAL*4 U_I_Shell
!     ' 
      REAL*4 U_alev_ld
!     ' 
      REAL*4 U_Temp
!     ' 
      REAL*4 U_Even_Odd
!     ' 
      REAL*4 BFTF
      REAL*4 BFTFA
      REAL*4 BFTFB
!     ' 
      REAL*4 Gaussintegral
!     ' 
!     /' Utility functions '/ 
!     ' 
!     ' 
      REAL*4 U_Box
      REAL*4 U_Box2
      REAL*4 U_Gauss
      REAL*4 U_Gauss_mod
      REAL*4 U_LinGauss
!     ' 
      INTEGER*4 U_Valid
!     ' 
!     ' 
      REAL*4 U_Delta_S0
!     ' 
!     ' 
!     ' 
!     /' Internal variables '/ 
      REAL*4  pi
      PARAMETER (pi=3.14159)
      INTEGER*4  I_N_CN
!     /' Neutron number of fissioning nucleus '/ 
      INTEGER*4  I,J,K
      REAL*4  T_coll_Mode_1,T_coll_Mode_2,T_coll_Mode_3,T_coll_Mode_4
      REAL*4  T_asym_Mode_1,T_asym_Mode_2,T_asym_Mode_3,T_asym_Mode_4,T_asym_Mode_0
      REAL*4  Sigpol_Mode_1,Sigpol_Mode_2,Sigpol_Mode_3,Sigpol_Mode_4
      REAL*4  R_Z_Curv_S0,R_Z_Curv1_S0,R_A_Curv1_S0
      REAL*4  ZC_Mode_0,ZC_Mode_1,ZC_Mode_2,ZC_Mode_3,ZC_Mode_4,ZC_Mode_4L
      REAL*4  SigZ_Mode_0,SigZ_Mode_1,SigZ_Mode_2,SigZ_Mode_3,SigZ_Mode_4
      REAL*4  SN
      REAL*4  E_exc_S0_prov,E_exc_S1_prov,E_exc_S2_prov,E_exc_S3_prov,E_exc_S4_prov
      REAL*4  E_exc_S11_prov,E_exc_S22_prov
      REAL*4  E_exc_Barr
      REAL*4  E_LD_S1,E_LD_S2,E_LD_S3,E_LD_S4
      REAL*4  R_Shell_S1_eff,R_Shell_S2_eff,R_Shell_S3_eff,R_Shell_S4_eff
      REAL*4  Yield_Norm
      REAL*4  R_E_exc_eff
      REAL*4  R_Z_Heavy,R_Z_Light
      INTEGER*4  I_Mode
      REAL*4  T_Pol_Mode_0,T_Pol_Mode_1,T_Pol_Mode_2,T_Pol_Mode_3,T_Pol_Mode_4
      REAL*4  E_Min_Barr
      REAL*4  RI
      REAL*4  rbeta,beta1,beta2
      REAL*4  rbeta_ld,rbeta_shell
      REAL*4  ZUCD
      REAL*4  Z
      REAL*4  E_tunn
      REAL*4  beta1_opt,beta2_opt,beta1_prev,beta2_prev
      REAL*4  Z1,Z2
      INTEGER*4  IZ1,IN1,IZ2,IN2
      REAL*4  A1,A2
      INTEGER*4  IA1,IA2
      REAL*4  E_defo
      REAL*4  R_Pol_Curv_S0,R_Pol_Curv_S1,R_Pol_Curv_S2,R_Pol_Curv_S3,R_Pol_Curv_S4
      REAL*4  RA,RZ
      REAL*4  SigA_Mode_0,SigA_Mode_1,SigA_Mode_2,SigA_Mode_3,SigA_Mode_4
      REAL*4  AC_Mode_0,AC_Mode_1,AC_Mode_2,AC_Mode_3,AC_Mode_4
      REAL*4  R_A_heavy,R_A_light
      REAL*4  RZpol
      REAL*4  T_intr_Mode_0,T_intr_Mode_1_heavy,T_intr_Mode_1_light
      REAL*4  T_intr_Mode_2_heavy,T_intr_Mode_2_light
      REAL*4  T_intr_Mode_3_heavy,T_intr_Mode_3_light
      REAL*4  T_intr_Mode_4_heavy,T_intr_Mode_4_light
      REAL*4  T
      REAL*4  DU0,DU1,DU2,DU3,DU4
      REAL*4  T_low_S1_used
      REAL*4  SigA_Mode_11,SigA_Mode_22
      INTEGER*4  Ngtot
      DATA Ngtot/0/
      INTEGER*4  Nglight
      DATA Nglight/0/
      INTEGER*4  Ngheavy
      DATA Ngheavy/0/
      REAL*4  Egtot1000
      DATA Egtot1000/0/
      REAL*4  S1_enhance,S2_enhance
      REAL*4  DZ_S2_lowE
      DATA DZ_S2_lowE/0/
      INTEGER*4  I_A_CN,I_Z_CN
      REAL*4  P_I_rms_CN
!     /' rms initial angular momentum '/ 
      DATA P_I_rms_CN/0/
!     ' 
!     ' Model parameters of GEF 
!     ' 
      REAL*4  Emax_valid
!     /' Maximum allowed excitation energy '/ 
      DATA Emax_valid/100/
      REAL*4  xP_DZ_Mean_S1
      DATA xP_DZ_Mean_S1/0.05/
      REAL*4  xP_DZ_Mean_S2
      DATA xP_DZ_Mean_S2/-1.0/
      REAL*4  xP_DZ_Mean_S3
!     /' Shift of mean Z of Mode 3 '/ 
      DATA xP_DZ_Mean_S3/0/
      REAL*4  xP_DZ_Mean_S4
!     /' Shell for structure at A around 190 '/ 
      DATA xP_DZ_Mean_S4/0/
      REAL*4  xP_Z_Curv_S1
      DATA xP_Z_Curv_S1/0.30/
      REAL*4  P_Z_Curvmod_S1
!     /' Scales energy-dependent shift '/ 
      DATA P_Z_Curvmod_S1/1.75/
      REAL*4  xP_Z_CurV_S2
      DATA xP_Z_CurV_S2/0.095/
      REAL*4  S2leftmod
!     /' Asymmetry in diffuseness of S2 mass peak '/ 
      DATA S2leftmod/0.55/
      REAL*4  S2leftmod_global
      DATA S2leftmod_global/0.6/
      REAL*4  P_Z_Curvmod_S2
!     /' Scales energy-dependent shift '/ 
      DATA P_Z_Curvmod_S2/10/
      REAL*4  xP_A_Width_S2
!     /' A width of Mode 2 (box) '/ 
      DATA xP_A_Width_S2/14.5/
!     '  Dim Shared As Single _P_Z_Curv_S3 = 0.076 
      REAL*4  xP_Z_Curv_S3
      DATA xP_Z_Curv_S3/0.068/
      REAL*4  P_Z_Curvmod_S3
!     /' Scales energy-dependent shift '/ 
      DATA P_Z_Curvmod_S3/10/
      REAL*4  P_Z_Curv_SL4
      DATA P_Z_Curv_SL4/0.28/
!     '  Dim Shared As Single _P_Z_Curv_S4 = 0.025  /' Curvature in Z of Mode 4 '/ 
      REAL*4  xP_Z_Curv_S4
      DATA xP_Z_Curv_S4/0.05/
      REAL*4  P_Z_Curvmod_S4
!     /' Scales energy-dependent shift '/ 
      DATA P_Z_Curvmod_S4/10/
      REAL*4  xDelta_S0
!     /' Shell effect for SL,for individual systems '/ 
      DATA xDelta_S0/0/
      REAL*4  xP_Shell_S1
!     /' Shell effect for Mode 1 '/ 
      DATA xP_Shell_S1/-1.85/
      REAL*4  xP_Shell_S2
!     /' Shell effect for Mode 2 '/ 
      DATA xP_Shell_S2/-4.0/
      REAL*4  xP_Shell_S3
!     /' Shell effect for Mode 3 '/ 
      DATA xP_Shell_S3/-6.0/
      REAL*4  P_Shell_SL4
!     /' Shell enhancing S1 '/ 
      DATA P_Shell_SL4/-1.3/
      REAL*4  xP_Shell_S4
!     /' Shell effect for Mode 4 '/ 
      DATA xP_Shell_S4/-1.0/
      REAL*4  PZ_S3_olap_pos
!     /' Pos. of S3 shell in light fragment (in Z) '/ 
      DATA PZ_S3_olap_pos/39.7/
      REAL*4  PZ_S3_olap_curv
!     /' for width of S3 shell in light fragment '/ 
      DATA PZ_S3_olap_curv/0.008/
      REAL*4  Level_S11
!     /' Level for mode S11 '/ 
      DATA Level_S11/-1.3/
      REAL*4  Shell_fading
!     /' fading of shell effect with E* '/ 
      DATA Shell_fading/50/
      REAL*4  xT_low_S1
      DATA xT_low_S1/0.342/
      REAL*4  xT_low_S2
!     /' Slope parameter for tunneling '/ 
      DATA xT_low_S2/0.31/
      REAL*4  xT_low_S3
!     /' Slope parameter for tunneling '/ 
      DATA xT_low_S3/0.31/
      REAL*4  xT_low_S4
!     /' Slope parameter for tunneling '/ 
      DATA xT_low_S4/0.31/
      REAL*4  xT_low_SL
!     /' Slope parameter for tunneling '/ 
      DATA xT_low_SL/0.31/
      REAL*4  T_low_S11
!     /' Slope parameter for tunneling '/ 
      DATA T_low_S11/0.36/
      REAL*4  xP_att_pol
!     /' Attenuation of 132Sn shell '/ 
      DATA xP_att_pol/4.5/
      REAL*4  dE_Defo_S1
!     /' Deformation energy expense for Mode 1 '/ 
      DATA dE_Defo_S1/-2.8/
      REAL*4  dE_Defo_S2
!     /' Deformation energy expense for Mode 2 '/ 
      DATA dE_Defo_S2/0/
      REAL*4  dE_Defo_S3
!     /' Deformation energy expense for Mode 3 '/ 
      DATA dE_Defo_S3/0/
      REAL*4  dE_Defo_S4
      DATA dE_Defo_S4/0/
!     /' Deformation energy expense for Mode 4 '/ 
!     '    Dim Shared As Single betaL0 = 24.5 
!     '    Dim Shared As Single betaL1 = 0.65 
!     '    Dim Shared As Single betaL0 = 26.7 
!     '    Dim Shared As Single betaL1 = 0.71 
      REAL*4  betaL0
      DATA betaL0/25.4/
      REAL*4  betaL1
      DATA betaL1/0.69/
      REAL*4  betaH0
!     /' Offset for deformation of heavy fragment '/ 
      DATA betaH0/48.0/
      REAL*4  betaH1
      DATA betaH1/0.55/
      REAL*4  kappa
!     /' N/Z dedendence of A-asym. potential '/ 
      DATA kappa/0/
      REAL*4  TCOLLFRAC
!     /' Tcoll per energy gain from saddle to scission '/ 
      DATA TCOLLFRAC/0.04/
      REAL*4  ECOLLFRAC
      DATA ECOLLFRAC/0.055/
!     'Dim Shared As Single ECOLLFRAC = 0.1 
      REAL*4  TFCOLL
      DATA TFCOLL/0.034/
      REAL*4  TCOLLMIN
      DATA TCOLLMIN/0.12/
      REAL*4  ESHIFTSASCI_intr
!     /' Shift of saddle-scission energy '/ 
      DATA ESHIFTSASCI_intr/-58/
      REAL*4  ESHIFTSASCI_coll
!     /' Shift of saddle-scission energy '/ 
      DATA ESHIFTSASCI_coll/-90/
      REAL*4  ESHIFTSASCI_coll_global
      DATA ESHIFTSASCI_coll_global/-90/
      REAL*4  EDISSFRAC
      DATA EDISSFRAC/0.35/
      REAL*4  SIGDEFO
      DATA SIGDEFO/0.165/
      REAL*4  SIGDEFO_0
      DATA SIGDEFO_0/0.165/
      REAL*4  SIGDEFO_slope
      DATA SIGDEFO_slope/0/
      REAL*4  EexcSIGrel
      DATA EexcSIGrel/0.7/
      REAL*4  DNECK
!     /' Tip distance at scission / fm '/ 
      DATA DNECK/1/
      REAL*4  FTRUNC50
!     /' Truncation near Z=50 '/ 
      DATA FTRUNC50/1/
      REAL*4  ZTRUNC50
!     /' Z value for truncation '/ 
      DATA ZTRUNC50/50/
      REAL*4  FTRUNC28
!     /' Truncation near Z=28 '/ 
      DATA FTRUNC28/0.56/
      REAL*4  ZTRUNC28
!     /' Z value for truncation '/ 
      DATA ZTRUNC28/30.5/
      REAL*4  ZMAX_S2
!     /' Maximum Z of S2 channel in light fragment '/ 
      DATA ZMAX_S2/60/
      REAL*4  NTRANSFEREO
!     /' Steps for E sorting for even-odd effect '/ 
      DATA NTRANSFEREO/6/
      REAL*4  NTRANSFERE
!     /' Steps for E sorting for energy division '/ 
      DATA NTRANSFERE/12/
      REAL*4  Csort
!     /' Smoothing of energy sorting '/ 
      DATA Csort/0.1/
      REAL*4  PZ_EO_symm
!     /' Even-odd effect in Z at symmetry '/ 
      DATA PZ_EO_symm/2.25/
      REAL*4  PN_EO_Symm
!     /' Even-odd effect in N at symmetry '/ 
      DATA PN_EO_Symm/0.5/
      REAL*4  R_EO_THRESH
!     /' Threshold for asymmetry-driven even-odd effect'/ 
      DATA R_EO_THRESH/0.04/
      REAL*4  R_EO_SIGMA
      DATA R_EO_SIGMA/0.35/
      REAL*4  R_EO_MAX
!     /' Maximum even-odd effect '/ 
      DATA R_EO_MAX/0.40/
      REAL*4  xPOLARadd
!     /' Offset for enhanced polarization '/ 
      DATA xPOLARadd/0.32/
      REAL*4  POLARfac
!     /' Enhancement of polarization of ligu. drop '/ 
      DATA POLARfac/1/
      REAL*4  T_POL_RED
!     /' Reduction of temperature for sigma(Z) '/ 
      DATA T_POL_RED/0.01/
      REAL*4  xHOMPOL
!     /' hbar omega of polarization oscillation '/ 
      DATA xHOMPOL/2.0/
      REAL*4  ZPOL1
!     /' Extra charge polarization of S1 '/ 
      DATA ZPOL1/0/
      REAL*4  P_n_x
!     /' Enhanced inverse neutron x section '/ 
      DATA P_n_x/0/
      REAL*4  Tscale
      DATA Tscale/1/
      REAL*4  EOscale
!     /' Scaling factor for even-odd structure in yields '/ 
      DATA EOscale/1.0/
      REAL*4  Econd
      DATA Econd/2/
      INTEGER*4  Emode
!     /' 0: E over BF_B; 1: E over gs; 2: E_neutron '/ 
      DATA Emode/1/
      REAL*4  T_orbital
!     /' From orbital ang. momentum '/ 
      DATA T_orbital/0/
      REAL*4  Jscaling
!     /' General scaling of fragment angular momenta '/ 
      DATA Jscaling/1.0/
      REAL*4  Spin_odd
!     /' RMS Spin enhancement for odd Z '/ 
      DATA Spin_odd/0.4/
!     ' 
!     /' I. Properties of nuclide distributions '/ 
!     ' 
      REAL*4 , DIMENSION(0:6,2,150) :: Beta
!     ' 
      REAL*4 , DIMENSION(0:4,2,150) :: Edefo
!     ' 
      REAL*4 , DIMENSION(0:4,2,350) :: Zmean
!     ' 
      REAL*4 , DIMENSION(0:4,2,350) :: Zshift
!     ' 
      REAL*4 , DIMENSION(0:4,2,350) :: Temp
!     ' 
      REAL*4 , DIMENSION(0:4,2,350) :: TempFF
!     ' 
      REAL*4 , DIMENSION(0:4,2,350) :: Eshell
!     ' 
      REAL*4 , DIMENSION(0:6,2,350) :: PEOZ
!     ' 
      REAL*4 , DIMENSION(0:6,2,350) :: PEON
!     ' pre-neutron evaporation 
!     ' 
      REAL*4 , DIMENSION(0:6,2,350) :: EPART
!     ' 
      REAL*4 , DIMENSION(0:6,2,1:200,1:150) :: SpinRMSNZ
!     ' 
!     ' 
!     /' Masses etc. '/ 
!     ' 
      REAL*4 , DIMENSION(0:203,0:136) :: BEldmTF
!     ' 
      REAL*4 , DIMENSION(0:203,0:136) :: BEexp
!     ' 
      REAL*4 , DIMENSION(0:203,0:136) :: DEFOtab
!     ' 
      REAL*4 , DIMENSION(0:203,0:136) :: ShellMO
!     ' 
      REAL*4 , DIMENSION(0:203,0:136) :: EVOD
!     ' 
!     ' 
      REAL*4 , DIMENSION(0:200,0:150) :: NZPRE
!     ' 
      REAL*4 , DIMENSION(0:6,0:200,0:150) :: NZMPRE
!     /' Internal parameters for error analysis: '/ 
      REAL*4  P_DZ_Mean_S1
      REAL*4  P_DZ_Mean_S2
      REAL*4  P_DZ_Mean_S3
      REAL*4  P_DZ_Mean_S4
      REAL*4  P_Z_Curv_S1
      REAL*4  P_Z_Curv_S2
      REAL*4  P_A_Width_S2
      REAL*4  P_Z_Curv_S3
      REAL*4  P_Z_Curv_S4
      REAL*4  Delta_S0
      REAL*4  P_Shell_S1
      REAL*4  P_Shell_S2
      REAL*4  P_Shell_S3
      REAL*4  P_Shell_S4
      REAL*4  T_low_S1
      REAL*4  T_low_S2
      REAL*4  T_low_S3
      REAL*4  T_low_S4
      REAL*4  T_low_SL
      REAL*4  P_att_pol
      REAL*4  HOMPOL
      REAL*4  POLARadd
!     ' 
!     ' 
!     /' Control parameters: '/ 
      REAL*4  B_F
!     /' Fission barrier '/ 
      DATA B_F/0/
      REAL*4  B_F_ld
!     /' Fission barrier,liquid drop '/ 
      DATA B_F_ld/0/
      REAL*4  E_B
!     /' Outer fission barrier '/ 
      DATA E_B/0/
      REAL*4  E_B_ld
!     /' Outer fission barrier,liquid drop '/ 
      DATA E_B_ld/0/
      REAL*4  R_E_exc_Eb
!     /' Energy above outer barrier '/ 
      DATA R_E_exc_Eb/0/
      REAL*4  R_E_exc_GS
!     /' Energy above ground state '/ 
      DATA R_E_exc_GS/0/
      REAL*4  P_Z_Mean_S0
!     /' Mean Z of Mode 1 '/ 
      DATA P_Z_Mean_S0/0/
      REAL*4  P_Z_Mean_S1
!     /' Mean Z of Mode 1 '/ 
      DATA P_Z_Mean_S1/52.8/
      REAL*4  P_Z_Mean_S2
!     /' Mean Z of Mode 2 '/ 
      DATA P_Z_Mean_S2/55/
      REAL*4  P_Z_Mean_S3
!     /' Mean Z of Mode 3 '/ 
      DATA P_Z_Mean_S3/65/
      REAL*4  P_Z_Mean_S4
!     /' Mean Z of Mode 4 '/ 
      DATA P_Z_Mean_S4/42.05/
      REAL*4  NC_Mode_0
!     /' Mean N of symm. Mode '/ 
      DATA NC_Mode_0/0/
      REAL*4  NC_Mode_1
!     /' Mean N of Mode 1 '/ 
      DATA NC_Mode_1/0/
      REAL*4  NC_Mode_2
!     /' Mean N of Mode 2 '/ 
      DATA NC_Mode_2/0/
      REAL*4  NC_Mode_3
!     /' Mean N of Mode 3 '/ 
      DATA NC_Mode_3/0/
      REAL*4  NC_Mode_4
      DATA NC_Mode_4/0/
      REAL*4  B_S1
!     /' Barrier S1,relative to SL '/ 
      DATA B_S1/0/
      REAL*4  B_S2
!     /' Barrier S2,relative to SL '/ 
      DATA B_S2/0/
      REAL*4  B_S3
!     /' Barrier S3,relative to SL '/ 
      DATA B_S3/0/
      REAL*4  B_S4
      DATA B_S4/0/
      REAL*4  B_S11
!     /' Barrier S11,relative to SL '/ 
      DATA B_S11/0/
      REAL*4  B_S22
!     /' Barrier S22,relative to SL '/ 
      DATA B_S22/0/
      REAL*4  DES11ZPM
!     /' Mod. of eff. barrier due to ZPM in overlap '/ 
      DATA DES11ZPM/0/
      REAL*4  Delta_NZ_Pol
!     /' Polarization for 132Sn '/ 
      DATA Delta_NZ_Pol/0/
      REAL*4  Yield_Mode_0
!     /' Relative yield of SL '/ 
      DATA Yield_Mode_0/0/
      REAL*4  Yield_Mode_1
!     /' Relative yield of S1 '/ 
      DATA Yield_Mode_1/0/
      REAL*4  Yield_Mode_2
!     /' Relative yield of S2 '/ 
      DATA Yield_Mode_2/0/
      REAL*4  Yield_Mode_3
!     /' Relative yield of S3 '/ 
      DATA Yield_Mode_3/0/
      REAL*4  Yield_Mode_4
!     /' Relative yield of S4 '/ 
      DATA Yield_Mode_4/0/
      REAL*4  Yield_Mode_11
!     /' Relative yield of S11 '/ 
      DATA Yield_Mode_11/0/
      REAL*4  Yield_Mode_22
!     /' Relative yield of S22 '/ 
      DATA Yield_Mode_22/0/
      REAL*4  P_POL_CURV_S0
!     /' Stiffnes in N/Z '/ 
      DATA P_POL_CURV_S0/0/
      REAL*4  T_Coll_Mode_0
!     /' Effective collective temperature '/ 
      DATA T_Coll_Mode_0/0/
      REAL*4  E_Exc_S0
!     /' Energy over barrier of symmetric channel '/ 
      DATA E_Exc_S0/0/
      REAL*4  E_Exc_S1
!     /' Energy over barrier of S1 channel '/ 
      DATA E_Exc_S1/0/
      REAL*4  E_Exc_S2
!     /' Energy over barrier of S2 channel '/ 
      DATA E_Exc_S2/0/
      REAL*4  E_Exc_S3
!     /' Energy over barrier of S3 channel '/ 
      DATA E_Exc_S3/0/
      REAL*4  E_Exc_S4
!     /' Energy over barrier of S4 channel '/ 
      DATA E_Exc_S4/0/
      REAL*4  E_Exc_S11
!     /' Energy over barrier of S11 channel '/ 
      DATA E_Exc_S11/0/
      REAL*4  E_Exc_S22
!     /' Energy over barrier of S22 channel '/ 
      DATA E_Exc_S22/0/
      REAL*4  E_POT_SCISSION
!     /' Potential-energy gain saddle-scission '/ 
      DATA E_POT_SCISSION/0/
      REAL*4  EINTR_SCISSION
!     /' Intrinsic excitation energy at scission '/ 
      DATA EINTR_SCISSION/0/
      REAL*4  EeffS2
!     /' Governs S1 reduction by pairing '/ 
      DATA EeffS2/0/
      REAL*4  Sigpol_Mode_0
!     /' Width of isobaric Z distribution '/ 
      DATA Sigpol_Mode_0/0/
!     ' 
!     '  #Include Once "BEldmTF.bas" 
!     ' 
!     '  #Include Once "BEexp.bas" 
!     ' 
!     '  #Include Once "DEFO.bas" 
!     ' 
!     '  #Include Once "ShellMO.bas" 
!     ' 
!     ' 
!     ' 
!     ' 
!     ' 
      INTEGER*4  I_E_iso
!     ' numbered in sequence of increasing energy 
      REAL*4  Spin_CN
      REAL*4  Spin_pre_fission
      REAL*4  Spin_gs_light
      REAL*4  Spin_gs_heavy
!     ' 
!     /' Shell effects for the symmetric fission channel '/ 
      REAL*4  R_E_exc_used
      REAL*4  R_Z_mod
      REAL*4  T_Rusanov
      REAL*4  R_E_intr_S1,R_E_intr_S2,R_E_intr_S3
!     ' intrinsic exc. energies at barrier 
      REAL*4  R_E_intr_S4
      REAL*4 , DIMENSION(6) :: R_Att
!     ' attenuation of shell 
      REAL*4 , DIMENSION(6) :: R_Att_Sad
!     '  Dim As Single E_backshift 
!     '  E_backshift = -3 
!     ' 
!     ' 
      REAL*4  DZ_S1,DZ_S2,DZ_S3,DZ_S4
      REAL*4  EtotS2
      REAL*4  P_Z_Curv_S1_eff
      REAL*4  P_Z_Curv_S2_eff
      REAL*4  P_Z_Curv_S3_eff
      REAL*4  P_Z_Curv_S4_eff
      REAL*4  Etot,E1FG,E1ES
      REAL*4  Rincr1P,Rincr1N,Rincr2,Rincr2P,Rincr2N
      REAL*4  T1,T2,E1,E2
      REAL*4 , DIMENSION(0:6) :: E_coll_saddle
      REAL*4  Ediff
!     ' 
!     ' 
      REAL*4  DT
!     ' 
      REAL*4  AUCD
!     /' UCD fragment mass '/ 
      REAL*4  I_rigid_spher
!     /' I rigid for spherical shape '/ 
      REAL*4  I_rigid
!     /' I rigid for deformed scission shape '/ 
      REAL*4  I_eff
!     /' I with reduction due to pairing '/ 
      REAL*4  alph
!     /' deformation parameter '/ 
      REAL*4  E_exc
!     /' Excitation energy '/ 
      REAL*4  J_rms
!     /' rms angular momentum '/ 
      INTEGER*4  Ic,Jc
      REAL*4  R_Help,Zs,R_Sum
!     ' 
      REAL*4  R_Cut1,R_Cut2
      INTEGER*4  N_index,Z_index,A_index,M_index
      REAL*4  Ymin
!     ' Minimum yield to be stored 
      DATA Ymin/1.E-7/
      REAL*4  Eexc_mean,Eexc_sigma
      REAL*4  Eexc_intr,Eexc_coll
!     ' 
!      INCLUDE "RESULTS.FOR"      
!     ' 
!     ' 
      REAL*4  Rint
      REAL*4  RS
      REAL*4  RintE
      REAL*4 F1
      REAL*4 F2

end module GEFSUBdcl2


!*******************************************
!*   Declarataion of the RESULTS Module    *
!*******************************************
! Taken from RESULTS.FOR
module RESULTS

      INTEGER*4  N_cases
!     ' Number of cases in NZMkey,Etab,Jtab and Ytab 
!     ' (First dimension of NZMkey,Etab,Jtab and Ytab) 
      INTEGER*4 , DIMENSION(10000,3) :: NZMkey
!     ' Key (Mode,N,Z) for E*,spin and yield distr. of fragments 
      REAL*4 , DIMENSION(10000,1000) :: Etab
!     ' Excitation-energy distribution of fragments (0.1 MeV bins) 
      REAL*4 , DIMENSION(10000,100) :: Jtab
!     ' Spin distribution of fragments 
!     ' (0 to 100 hbar for even-A or 1/2 to 201/2 hbar for odd-A nuclei) 
      REAL*4 , DIMENSION(10000) :: Ytab

end module RESULTS


!*****************************************
!*       Init Function for Python        *
!*****************************************
!subroutine init_decay_module()

!    integer iii
!    iii = 1

!end


!*****************************************
!*              Initialize               *
!*****************************************
subroutine initialize(f_network, f_format, f_path)

    use iso
    use path_module
    implicit none

    character*256    f_network
    integer          f_format
    character*256    f_path
    character*10     re_import_keyphrase

    ! Make sure not to reimport the module if
    ! already imported and initialized.  This is to
    ! solve a problem when calling NuPyCEE multiple 
    ! times while it import the decay module each time.
    if (re_import_test .ne. 'keyphrase') then

        re_import_test = 'keyphrase'
        files_path = f_path
        input_file_network = trim(files_path) // 'decay_data/' // f_network
        abundance_file_format = f_format

        call init_variables()
        call read_elements()
        call init_network()
!        call special_fission() ! in special_fission.f95
    !    call init_n_fission() ! in n_fission.f95

        !print*, 'in init! ', files_path, f_path, re_import_test

    end if

end subroutine initialize


!*****************************************
!*              Run Decay                *
!*****************************************
subroutine run_decay(t_decay, t_steps, input_abundance)

    use iso
    use comm
    implicit none
    double precision t_decay
    integer          t_steps
    double precision input_abundance(max_number_isotopes)
    integer isot
 
    time_decay = t_decay
    time_steps = t_steps
    dt = time_decay / time_steps

    ! Re-initialize decay arrays
    do isot=0,max_entry
        i_entry(isot) = -1
    end do
    i_entry_max = 0
    do isot=1,number_isotopes
        abundance(isot) = 0.0
        initial_abundance(isot) = 0.0
    end do

    call init_abundances(input_abundance)
    call nucleosynthesis()

end


!*****************************************
!*         Get Non-Zero Indexes          *
!*****************************************
subroutine get_non_zero_indexes()

    use iso
    use comm

    i_entry_max = 0
    do isot=1,number_isotopes
      if (initial_abundance(isot) .gt. 0.0 .or. abundance(isot) .gt. 1.d-95) then
          i_entry(i_entry_max) = isot
          i_entry_max = i_entry_max + 1
      end if
    end do

end 


!*****************************************
!*            Get Iso Names              *
!*****************************************
subroutine get_iso_names(a_return)

    use iso
    implicit none
    character*6 a_return(max_number_isotopes)
    character*10 a_number
    integer temp
    integer isot
    integer ii

!f2py intent(out) a_return
    print*, 'Write function to calculate the len of non-zero isotopes'
!    do i=1,max_number_isotopes
    ii = 1
    do isot=1,number_isotopes
      if (initial_abundance(isot) .gt. 0.d0 .or. abundance(isot) .gt. 1.d-95) then
        print*, z(isot), z(isot)+n(isot), '  ',element_names(z(isot)), &
                   level(isot), abundance(isot), initial_abundance(isot)
        temp = z(isot)+n(isot)
        print*, z(isot)+n(isot)
        print*, temp
        write(a_number, '(i10)') temp
        print*, a_number
!        write(a_number, 10) temp
! 10     format (I4)
        a_return(ii) = trim(element_names(z(isot))) // '-' // trim(a_number)
        ii = ii + 1
      end if  
    end do

    do isot=1,ii
        print*, a_return(isot)
    end do

end


!*****************************************
!*            Init Variables             *
!*****************************************
subroutine init_variables()

    use iso
    implicit none
    integer i
    double precision en

    z = 0
    n = 0
!    level = 0.
    do i=1,max_number_isotopes
      level(i) = 0.0
    end do
    reactions = 0
    decay_constant = 0.
    decay = 0.
    production = 0.
    abundance = 0.
    initial_abundance = 0.
    s_fission_vector = 0.
    do i=1,1000
      en = 0.1d0 * (i - 0.5)
      mb_distribution(i) = dsqrt(en) * dexp(-en/kt)
    end do  
    mb_distribution = mb_distribution / sum(mb_distribution)

end


!*****************************************
!*            Read Elements              *
!*****************************************
subroutine read_elements()

    use iso
    use path_module
    implicit none
    integer iostats,i
    character*256 the_path

    the_path = trim(files_path) // trim('decay_data/element_symbols.txt')

    open(unit=13,file=the_path,status='unknown')
      iostats=0
      i=0
      element_names(0) = 'NN'
      do while (iostats.eq.0)
        i=i+1
        read(13,*,iostat=iostats) element_names(i)
      end do
    close(13)

end


!*****************************************
!*            Set Init Files             *
!*****************************************
subroutine set_init_files(f_network, f_format, t_decay, t_steps)

    ! Import list of parameters
    use iso

    ! Define the arguments
    implicit none
    character*256    f_network
    integer          f_format
    double precision t_decay
    integer          t_steps

    print*, 'in set_init_files'

    ! Set the value of the parameters
    input_file_network = f_network
    abundance_file_format = f_format
    time_decay = t_decay
    time_steps = t_steps

end


!*****************************************
!*            Init Abundances            *
!*****************************************
subroutine init_abundances(input_abundance)
    use iso
    implicit none
    double precision input_abundance(max_number_isotopes)
    integer iostats, ind

!    print*, 'in init_abundances'
!    input_file_abundance = 'data/iso_mass_U238.DAT'
!    input_file_abundance = 'data/iso_mass_Fe60.DAT'

!    open(unit=13,file=input_file_abundance,status='old')
!      select case (abundance_file_format)
!        case(1)  
!          call read_nugrid_abundances()
!        case default
!          print*,abundance_file_format, ' is no valid file format'  
!          stop
!      end select
!    close(13)  
!    abundance = initial_abundance

    abundance = input_abundance

end


!*****************************************
!*        Read NuGrid Abundances         *
!*****************************************
subroutine read_nugrid_abundances()

    use iso
    implicit none
    integer iostats, ind, al, zl, nl, ll, get_isotope
    double precision levell, abundancel
    character*42 c42 

    iostats=0
    do ind=1,7
      read(13,*)
    end do
    ind=0
    do while (iostats.eq.0) 
      ind = ind+1
      read(13,'(A42)',iostat=iostats) c42
      read(c42(7:11),*) zl
      read(c42(13:16),*) al
      nl = al-zl
      read(c42(19:21),*) ll
      levell = 0
      if (zl.eq.73 .and. al.eq.180 .and.  ll.eq.1) levell = 0.0771d0  ! 180Ta, isomer
      if (zl.eq.13 .and. al.eq.26  .and.  ll.eq.2) levell = 0.2283d0  ! 26Al, isomer
      if (zl.eq.36 .and. al.eq.85  .and.  ll.eq.2) levell = 0.305d0   ! 85Kr, isomer
      if (zl.eq.48 .and. al.eq.115 .and.  ll.eq.2) levell = 0.181d0    ! 115Cd, isomer
      if (zl.eq.71 .and. al.eq.176 .and.  ll.eq.2) levell = 0.1229d0   ! 176Lu, isomer
      read(c42(23:35),*) abundancel
      current_isotope = get_isotope(zl,nl,levell)
      if (current_isotope .eq. 0) then
        print*,'isotope not found: ',ind, zl, nl, levell, c42
        stop
      else
        initial_abundance(current_isotope) = abundancel
      end if  
    end do

end


!*****************************************
!*             Init Network              *
!*****************************************
subroutine init_network()

    use iso
    implicit none
    integer isomer, reac, get_isotope, zt, nt, fission_counter
    double precision norm

    !print*, 'in init_network'

    call read_network()
    norm = 0.
    fission_counter = 0
    index_p = get_isotope(1,0,norm)
    index_n = get_isotope(0,1,norm)
    index_a = get_isotope(2,2,norm)
    index_12c = get_isotope(6,6,norm)
    do isomer = 1, number_isotopes           ! loop over all species
      if (reactions(isomer,0) .gt. 0) then   ! unstable isotope
        norm = 1. / sum(decay_constant(isomer,1:reactions(isomer,0)))
        if (norm .le. 0.) then
          print *,norm
          print *,z(isomer),n(isomer)
          stop
        end if  
        decay_constant(isomer,1:reactions(isomer,0)) =  decay_constant(isomer,1:reactions(isomer,0))*norm  ! now all decay rates normalized to 1
        do reac = 1,reactions(isomer,0)      ! loop over all decay channels
          zt = z(isomer) + reaction_vector(1,reactions(isomer,reac)) ! Z of reaction product
          nt = n(isomer) + reaction_vector(2,reactions(isomer,reac)) ! N of reaction product
          product_isomer(isomer, reac) = get_isotope(zt,nt,level_product(isomer, reac))
          if (reactions(isomer,reac) .eq. 22) then
            fission_counter = fission_counter+1
            if (fission_counter .gt. max_number_fissile_isotopes) then
              print*,'number of fissile iostopes .gt. max_number_fissile_isotopes ', max_number_fissile_isotopes
              print*,'please change in the source code and recompile'
              stop
            end if  
            reactions(isomer,-1) = fission_counter                ! index for fission yields
            call init_fission(isomer)
          end if
        end do
      end if
    end do
!    print*, 'Number of fissile isotopes: ',fission_counter
!    dt = time_decay / time_steps

end


!*****************************************
!*             Read Network              *
!*****************************************
subroutine read_network()

    use iso
    implicit none
    integer iostats,io,ind, al,zl, nl, dm
    integer character_to_decay_mode, get_isotope
    double precision level_energy, hwzl, macs, target_level,dp
    double precision character_to_decay_probability,level_spin
    double precision get_spin
    character*4 decay_mode
    character*14 decay_prob
    character*159 c159

    open(unit=13,file=input_file_network,status='unknown')
      iostats=0
      ind=0
      read(13,*)
      do while (iostats.eq.0)
        ind=ind+1
        if (ind.gt.max_number_isotopes) then
          print*,'number of species greater than ', max_number_isotopes
          print*,'change in source code and recompile'
          stop
        end if  
        read(13,'(A159)',iostat=iostats) c159
        read(c159(8:11),*) zl
        read(c159(19:22),*) al
        nl = al-zl
        read(c159(37:54),*) level_energy
        level_spin = get_spin(c159(55:71))
        read(c159(72:85),*) hwzl
        if (max_hwz.lt.hwzl) hwzl = -1
        decay_mode = c159(109:112)
        decay_prob = c159(114:127)
     
        read(c159(145:159),*) target_level
        current_isotope = get_isotope(zl,nl,level_energy)
        if ((hwzl .gt. 0.) .and. (index(decay_mode,'?') .eq. 0.)) then
          ! unstable isotope, decay mode confirmed, bound nucleus
          dp         = character_to_decay_probability(decay_prob)
          dm         = character_to_decay_mode(decay_mode)
          if (current_isotope .eq. 0) then ! new unstable species
            call new_unstable(zl,nl,level_energy,level_spin,hwzl,dm,dp, target_level)
          else
            call add_unstable(dm,dp,target_level)
          end if  
        end if  
        if ((hwzl .eq. -1) ) then  ! stable isotope
          if (current_isotope .eq. 0) then ! new stable species
            call new_stable(zl,nl,level_energy)
          else 
            print*,'Warning: stable - second time??!! Entry:',ind, trim(element_names(zl)) , zl+nl
            print*,'will ignore entry'
!            stop
          end if  
        end if  
      end do 
    close(13)

end

!*****************************************
!*             Add Unstable              *
!*****************************************
subroutine add_unstable(dm,dp, target_level)

    use iso
    implicit none
    integer dm
    double precision target_level, dp

    reactions(current_isotope,0) = reactions(current_isotope,0) + 1
    reactions(current_isotope,reactions(current_isotope,0)) = dm
    decay_constant(current_isotope,reactions(current_isotope,0)) = dp  ! to be normalized later
    level_product(current_isotope,reactions(current_isotope,0)) = target_level

end


!*****************************************
!*             New Unstable              *
!*****************************************
subroutine new_unstable(zl,nl,level_energy,level_spin,hwzl,dm,dp, target_level)
    use iso
    implicit none
    integer zl,nl,dm
    double precision level_energy,hwzl,target_level, dp,level_spin

    number_isotopes                   = number_isotopes +1
    current_isotope                   = number_isotopes
    z(current_isotope)                = zl
    n(current_isotope)                = nl
    level(current_isotope)            = level_energy
    spin(current_isotope)             = level_spin
    reactions(current_isotope,0)      = 1
    reactions(current_isotope,1)      = dm
    decay_constant(current_isotope,0) = dlog(2.d0)/hwzl  ! lambda
    decay_constant(current_isotope,1) = dp  ! to be normalized later
    level_product(current_isotope,1)  = target_level

end


!*****************************************
!*        Character to Decay Mode        *
!*****************************************
integer function character_to_decay_mode(c4)

    use iso
    implicit none
    character*4 c4
    integer i

    character_to_decay_mode = 0
    do i=3,1,-1
      if (c4(i:i) .eq. ' ') then  ! get rid of leading or interspersed spaces
        c4(i:3) = c4(i+1:4)
        c4(4:4) = ' '
      end if  
    end do  
    do i=1, number_reaction_types
!      if (index(c4,trim(reaction_types(i))) .gt. 0) then
      if (trim(c4) .eq. trim(reaction_types(i)) ) then
        character_to_decay_mode = i
        exit
      end if   
    end do
    if (character_to_decay_mode .eq. 0) then
      print*,'Reaction type ', c4 ,' not found! '
      stop
    end if  

end


!*****************************************
!*    Character to Decay Probability     *
!*****************************************
double precision function character_to_decay_probability(c14)

    implicit none
    character*14 c14
    integer io

    if ((index(c14,'?') .gt. 0 ) .or. (trim(c14).eq.'') ) then
!      print*,c14
      character_to_decay_probability = 1.
    else
      c14(index(c14,'~'):index(c14,'~')) = ' ' 
      c14(index(c14,'='):index(c14,'=')) = ' ' 
      c14(index(c14,'<'):index(c14,'<')) = ' ' 
      c14(index(c14,'>'):index(c14,'>')) = ' ' 
      read(c14,*,iostat=io) character_to_decay_probability  
      if (io.ne.0) print*,'io',c14
    end if  
    if (character_to_decay_probability .eq. 0.) character_to_decay_probability = 1.d-10

end


!*****************************************
!*              Get Isotope              *
!*****************************************
integer function get_isotope(zl,nl,level_energy)

    use iso
    implicit none
    integer zl,nl,i
    double precision level_energy

    get_isotope = 0
    do i=1, number_isotopes
      if (zl.eq.z(i) .and. nl.eq.n(i) .and. level_energy.eq.level(i)) then
        get_isotope = i
        exit
      end if   
    end do

end


!*****************************************
!*               Get Spin                *
!*****************************************
double precision function get_spin(cin)

    implicit none
    integer i,j,number_symbols,char_length
    parameter(char_length=17)
    character*17 cin
    logical get_rid
    double precision n,z
    parameter(number_symbols=13)
    character*1 symbols(number_symbols)
    data symbols/ ' ', '(', ')', '[', ']', '+','-', 'G','E','L', '<','>','|'/ 

    get_spin = 0.
    do i=char_length,1,-1
      get_rid = .false.
      do j=1,number_symbols
        get_rid = get_rid .or. (cin(i:i) .eq. symbols(j))
      end do
      if (get_rid) then  ! get rid of leading or interspersed spaces and ()+-
        cin(i:char_length-1) = cin(i+1:char_length)
        cin(char_length:char_length) = ' '
      end if  
    end do  
    if (trim(cin) .ne. '') then
      i=index(cin,'/')
      if (i.eq.0) then
        read(cin,*) get_spin
      else
        read(cin(1:i-1),*) n 
        read(cin(i+1:i+1),*) z
        get_spin = n/z
      end if 
    end if

end


!*****************************************
!*             Init Fission              *
!*****************************************
subroutine init_fission(isomer)

    use iso
    use RESULTS
    use path_module
    implicit none
    character*256 fission_file
    character*3 chz, chn
    character*11 chl, chs
    integer isomer,M_index, N_index, Z_index, A_index,k,i,get_isotope,li
    integer fission_index, istatus
    double precision yield, dytab(10000)
    REAL*4    level_e, level_spin
    INTEGER*4 iso_n, iso_z, iso_a
    character*256 the_path

    the_path = trim(files_path) // trim('decay_data/fission/')

    iso_z = z(isomer)
    iso_n = n(isomer)
    level_e    = level(isomer) / 1000. ! keV -> MeV
    level_spin = spin(isomer)
    fission_index = reactions(isomer, -1)

!    iso_z = 92
!    iso_n = 143
!    iso_n = 146
!    level_e    = 0.
!    level_spin = 0.

    iso_a = iso_z+iso_n

    write(chz,'(I3.3)') iso_z
    write(chn,'(I3.3)') iso_n
    write(chl,'(ES11.5)') level_e
    write(chs,'(ES11.5)') level_spin
    fission_file = trim(the_path) // chz // '_' // chn // '_' // chl //  '_'  // chs // '.dat'
    open(unit=13,file=fission_file,status='old', iostat=istatus)
    if (istatus .ne. 0) then                      ! need to create the fission file
      print*,'Fission of: ',iso_a,element_names(iso_z), '  ', level_e,level_spin
      print*,'Creating: ',trim(fission_file)
      call GEFSUB(iso_z,iso_a,level_e,level_spin)
      dytab = Ytab
      dytab = 2./sum(dytab(1:N_cases)) * dytab   ! normalize to 2.
      call scission(fission_index,dytab)
      print*,'Neutron yield',s_fission_vector(fission_index,index_n)
      open(unit=13,file=fission_file,status='unknown')
        do li=1,number_isotopes
          if (s_fission_vector(fission_index,li) .gt. 0.0) then
            write(13,'(2I6,ES20.10)') z(li), n(li), s_fission_vector(fission_index,li)
          end if  
        end do
    else
      do 
        read(13,*,iostat=istatus) Z_index,N_index, yield
        if (istatus .eq. 0) then
          li = get_isotope(Z_index,N_index,0.0d0)
          if (li .eq. 0) then
            print*,Z_index+N_index,element_names(Z_index)
            print*,'Isotope with fission yield not in reaction network.'
            print*,'Consider deleting and recreating fission data with current network file.'
          else  
            s_fission_vector(fission_index,li) = yield
          end if  
        else
          exit
        end if    
      end do
    end if    
!    yield = sum(s_fission_vector(fission_index,1:index_n-1)) &
!          + sum(s_fission_vector(fission_index,index_n+1:number_isotopes) )
!    print*,  iso_z+iso_n,element_names(iso_z), yield
    close(13)

end


!*****************************************
!*              New Stable               *
!*****************************************
subroutine new_stable(zl,nl,level_energy)

    use iso
    implicit none
    integer zl,nl
    double precision level_energy

    number_isotopes              = number_isotopes +1
    current_isotope              = number_isotopes
    z(current_isotope)           = zl
    n(current_isotope)           = nl
    level(current_isotope)       = level_energy
    reactions(current_isotope,0) = 0

end


!*****************************************
!*               Scission                *
!*****************************************
subroutine scission(fission_index,dytab)

    use iso
    use RESULTS 
    implicit none
    integer isomer,M_index, N_index, Z_index, A_index,k,i,j,get_isotope,li,ch_sn, n_ch_n
    integer fission_index, istatus, ch_counter, li_ne, neutron_counter
    double precision neutron_yield, e_ave, e_sum, de_ave2,neutron_separation_energy
    double precision sn, e_excitation(2,1000), e_exc, en(1000),dytab(10000)
    REAL*4    level_e, level_spin
    INTEGER*4 iso_n, iso_z, iso_a

!    INCLUDE "RESULTS.FOR" 

    DO K = 1 , N_cases    
      e_excitation(1,:) = Etab(k,:)
      e_excitation(1,:) = e_excitation(1,:) / sum(e_excitation(1,:)) * dytab(K)
      e_excitation(2,:) = 0.
      M_index = NZMkey(K,1) 
      N_index = NZMkey(K,2) 
      Z_index = NZMkey(K,3) 
      A_index = N_index + Z_index 
      sn    = neutron_separation_energy(Z_index,N_index)
      li    = get_isotope(Z_index,N_index,0.0d0)
      li_ne = get_isotope(Z_index,N_index-1,0.0d0)
      if (li*li_ne.eq.0) then
        neutron_counter = 0
        do
          neutron_counter = neutron_counter + 1
          li    = get_isotope(Z_index,N_index-neutron_counter,0.0d0)
          li_ne = get_isotope(Z_index,N_index-neutron_counter-1,0.0d0)
          sn = sn + neutron_separation_energy(Z_index,N_index-neutron_counter)
          s_fission_vector(fission_index,index_n) = s_fission_vector(fission_index,index_n) + dytab(K)
          if (li*li_ne.ne.0) then
            exit
          end if
        end do   
      end if  
      ch_sn = int(sn*10. + 0.5)
      e_ave = 0.
      do i=1,1000
        if (e_excitation(1,i) .gt. 0.) then
          e_exc = 0.1d0 * (i-0.5)
          e_ave = e_ave + e_excitation(1,i) * e_exc
          if (e_exc.le.sn) then  !  below neutron separation threshold
            neutron_yield = 0.
          else  ! neutron emission
            neutron_yield  = sum(mb_distribution(ch_sn:i)) / sum(mb_distribution(1:i)) 
            n_ch_n         = i-ch_sn+1
            en(1:n_ch_n)   = mb_distribution(ch_sn:i)  / sum(mb_distribution(1:i))
            do j=1,n_ch_n
              e_excitation(2,i-j-ch_sn) = e_excitation(2,i-j-ch_sn) + e_excitation(1,i)*en(j)
            end do
          end if
          s_fission_vector(fission_index,li)      = s_fission_vector(fission_index,li)      + e_excitation(1,i)*(1-neutron_yield)
          s_fission_vector(fission_index,index_n) = s_fission_vector(fission_index,index_n) + e_excitation(1,i)*neutron_yield
          s_fission_vector(fission_index,li_ne)   = s_fission_vector(fission_index,li_ne)   + e_excitation(1,i)*neutron_yield
        end if  
      end do
   
!
!  repeat this loop until all energy is lt separation energy
!
      do 
        N_index = N_index - 1
        sn    = neutron_separation_energy(Z_index,N_index)
        ch_sn = int(sn*10. + 0.5)
        if (sum(e_excitation(2,ch_sn:1000)) .eq. 0.) exit
        e_excitation(1,:) = e_excitation(2,:) / sum(e_excitation(1,:))
        e_excitation(2,:) = 0.
        li = li_ne
        if (li.eq.0) then
          print*,'isotope after multi neutron emission not found!!!',Z_index,N_index-1
          stop
        end if
        do i=1,1000
          if (e_excitation(1,i) .gt. 0.) then
            e_exc = 0.1d0 * (i-0.5)
            if (e_exc.gt.sn) then  ! neutron emission
              neutron_yield = sum(mb_distribution(ch_sn:i)) / sum(mb_distribution(1:i))
              n_ch_n         = i-ch_sn+1
              en(1:n_ch_n)   = mb_distribution(ch_sn:i)  / sum(mb_distribution(1:i))
              do j=1,n_ch_n
                e_excitation(2,i-j-ch_sn) = e_excitation(2,i-j-ch_sn) + e_excitation(1,i)*en(j)
              end do
              s_fission_vector(fission_index,li)      = s_fission_vector(fission_index,li)      - e_excitation(1,i)*neutron_yield
              s_fission_vector(fission_index,index_n) = s_fission_vector(fission_index,index_n) + e_excitation(1,i)*neutron_yield
              s_fission_vector(fission_index,li_ne)   = s_fission_vector(fission_index,li_ne)   + e_excitation(1,i)*neutron_yield
            end if  
          end if  
        end do
      end do  
    END DO

end


!*****************************************
!*       Neutron Separation Energy       *
!*****************************************
double precision function neutron_separation_energy(z,n)
! Taken from particle_separation_energies.f95
!
!  adopted from K Vogt, T Hartmann, A Zilges, PLB 517 (2001) 255-260
!
!  return value: neutron separation energy in MeV

    implicit none
    integer z,n
    double precision a(3), a_pair, a_shell, d_pair, d_shell, nn, mass, n_over_z                         
    data a/ 6.29  , 3.43, 5.85 / ! parameters for equation as given in paper
    data a_pair/ 10.59 / ! parameters for equation as given in paper
    data a_shell/ 1.51 / ! parameters for equation as given in paper

    nn = 0.
    if (n .gt. 28)  nn = nn + 1
    if (n .gt. 50)  nn = nn + 1
    if (n .gt. 82)  nn = nn + 1
    if (n .gt. 126) nn = nn + 1
    d_shell = nn * a_shell

    mass = z+n
    if (modulo(n,2) .eq. 0 ) then    ! even N
      d_pair =  a_pair /dsqrt(mass)
    else  
      d_pair = -a_pair /dsqrt(mass)
    end if

    n_over_z = float(n)/float(z)
    neutron_separation_energy = a(1)
    neutron_separation_energy = neutron_separation_energy + a(2)*mass**(1./3.)
    neutron_separation_energy = neutron_separation_energy / n_over_z
    neutron_separation_energy = neutron_separation_energy - a(3)
    neutron_separation_energy = neutron_separation_energy + d_pair
    neutron_separation_energy = neutron_separation_energy - d_shell
    neutron_separation_energy = neutron_separation_energy 

end


!*****************************************
!*            Nucleosynthesis            *
!*****************************************
subroutine nucleosynthesis()

    use iso
    implicit none
    integer step, s_out
    double precision levell


    s_out=max(1,time_steps / 10)
    do step =1,time_steps
!      if ( modulo(step, s_out) .eq. 0) then
!        print*,step,' steps **********************************  '
!      end if  
      call make_production_rates()
      abundance = abundance + production - decay
    end do 
!    call equilibrium_test()

end 


!*****************************************
!*            Equilibrium Test           *
!*****************************************
subroutine equilibrium_test()

    use iso
    implicit none
    integer isomer, reac, get_isotope, zt, nt, fission_index
    double precision decayl

    do isomer = 1, number_isotopes  ! loop over all species
      if (abundance(isomer).gt.0.) then  ! isomer is abundant
        if (reactions(isomer,0).eq.0) then  ! isomer is stable
          print*, z(isomer)+n(isomer), element_names(z(isomer)), &
          decay_constant(isomer,0) , abundance(isomer), short_lived(isomer)
        else
          print*, z(isomer)+n(isomer), element_names(z(isomer)), &
          decay_constant(isomer,0) , abundance(isomer), short_lived(isomer)&
          , decay_constant(isomer,0) * abundance(isomer)
        end if  
      end if
    end do

end


!*****************************************
!*         Make Production Rates         *
!*****************************************
subroutine make_production_rates()

    use iso
    implicit none
    integer isomer, reac, get_isotope, zt, nt, fission_index
    double precision decayl, test_abundance
    integer short_lived_isotopes(max_number_isotopes), short1, short2, isomer_short1,isomer_short2, i
    logical hit

    decay      = 0.
    production = 0.   
    short_lived_counter = 0
    short_lived = .false.
    do isomer = 1, number_isotopes  ! loop over all species
      if (abundance(isomer).gt.0. .and. reactions(isomer,0).gt.0) then ! isomer is unstable and abundant
        if (decay_constant(isomer,0) * dt .gt. 1.d-3) then
!          decay(isomer) = (1.d0 - dexp(-decay_constant(isomer,0) * dt))*abundance(isomer)
          short_lived_counter = short_lived_counter + 1
          short_lived(isomer) = .true.
          short_lived_isotopes(short_lived_counter) = isomer
!          print*, element_names(z(isomer)), z(isomer)+n(isomer)
        else
          decay(isomer) = decay_constant(isomer,0) * dt * abundance(isomer)  ! linear approximation
          call make_production(isomer)
        end if  
      end if
    end do  
    short1 =  short_lived_counter
    do while (short_lived_counter .gt. 0)  ! loop over all short-lived species
      isomer_short1 = short_lived_isotopes(short1)
!      print*, z(isomer_short1)+n(isomer_short1),element_names(z(isomer_short1)), ' check .. ',short_lived_counter,isomer_short1
      hit = .false.
      do short2 = 1,short_lived_counter  ! loop over all short-lived species
        if (short2 .ne. short1) then
          isomer_short2 = short_lived_isotopes(short2)
          do reac = 1,reactions(isomer_short2,0)      ! loop over all decay channels
            hit = product_isomer(isomer_short2, reac) .eq. isomer_short1 
            if (hit) exit
          end do
        end if
        if (hit) exit  
      end do  
      if (.not. hit) then                                                           ! isotope is not produced by short-lived isotopes
!        print*, z(isomer_short1)+n(isomer_short1),element_names(z(isomer_short1)), ' found !! ',short_lived_counter
        decay(isomer_short1) = (1.d0 - dexp(-decay_constant(isomer_short1,0) * dt))* &
          (abundance(isomer_short1)-production(isomer_short1)/dt/&
           decay_constant(isomer_short1,0)) + production(isomer_short1)
        call make_production(isomer_short1)  
         do i=short1,short_lived_counter-1
          short_lived_isotopes(i) = short_lived_isotopes(i+1)
        end do                     
        short_lived_counter = short_lived_counter - 1
        short1              = short_lived_counter
!        stop
       else
        short1 = short1 - 1
        if (short1 .lt. 1) then
          print*,'no solution found', short1, short_lived_counter
          do i=1,short_lived_counter
            isomer_short2 = short_lived_isotopes(i)
            print*,z(isomer_short2)+n(isomer_short2),&
                   element_names(z(isomer_short2)),isomer_short2
          end do
          print*
          stop
        end if  
      end if
    end do

end


!*****************************************
!*           Make Production             *
!*****************************************
subroutine make_production(isomer)

     use iso
     implicit none
     integer isomer, reac, get_isotope, zt, nt, fission_index
     double precision decayl

     do reac = 1,reactions(isomer,0)      ! loop over all decay channels
       decayl = decay(isomer)*decay_constant(isomer, reac)
       if (reactions(isomer,reac) .ne. 22) then  ! not fission
         production(product_isomer(isomer, reac)) = &
             production(product_isomer(isomer, reac)) + decayl   
       end if
       select case (reactions(isomer,reac))
         case(3)  ! N
           production(index_n) = production(index_n) + decayl
         case(4)  ! P
           production(index_p) = production(index_p) + decayl
         case(5)  ! A
           production(index_a) = production(index_a) + decayl
         case(6)  ! BN
           production(index_n) = production(index_n) + decayl
         case(7)  ! EP
           production(index_p) = production(index_p) + decayl
         case(8)  ! BA
           production(index_a) = production(index_a) + decayl
         case(9)  ! EA
           production(index_a) = production(index_a) + decayl
         case(10)  ! 2N
           production(index_n) = production(index_n) + 2.*decayl
         case(11)  ! 2P
           production(index_p) = production(index_p) + 2.*decayl
         case(12)  ! 2A
           production(index_a) = production(index_a) + 2.*decayl
         case(13)  ! B2A
           production(index_a) = production(index_a) + 2.*decayl
         case(14)  ! B2N
           production(index_n) = production(index_n) + 2.*decayl
         case(15)  ! B3N
           production(index_n) = production(index_n) + 3.*decayl
         case(16)  ! B4N
           production(index_n) = production(index_n) + 4.*decayl
         case(17)  ! E2P
           production(index_p) = production(index_p) + 2.*decayl
         case(18)  ! BNA
           production(index_n) = production(index_n) + decayl
           production(index_a) = production(index_a) + decayl
         case(19)  ! EPA
           production(index_p) = production(index_p) + decayl
           production(index_a) = production(index_a) + decayl
         case(21)  ! 12C
           production(index_12c) = production(index_12c) + decayl
         case(22)  ! fission
           fission_index = reactions(isomer, -1)
           production = production + decayl * s_fission_vector(fission_index,:)
       end select
     end do

end


!*****************************************
!*                Output                 *
!*****************************************
subroutine output()

    use iso
    use path_module
    implicit none
    integer index1, index2, index3,get_isotope,isot
    double precision levell
    character*256 the_path

    levell = 0.
    index1 = get_isotope(36,51,levell) !87Kr
    index2 = get_isotope(37,50,levell) !87Rb 
    index3 = get_isotope(38,49,levell) !87Sr 
! 
!    print*,'Output - abu  :  ',abundance(index1),abundance(index2),abundance(index3)
    print*, 'Reaction file: ',trim(input_file_network)
    print*, 'Initial abundance file: ',trim(input_file_abundance)
    print*, 'Decay time: ',time_decay, &
         ' s           (',time_decay/365.25/24/3600,' y , ',&
         time_decay/365.25/24/3600/1.d9,' Gy)'
    print*, 'Simulation done in ', time_steps, ' steps'
    print*, '    Z     A Symbol   Level            Final-Abundance  Initial-Abundance'
    do isot=1,number_isotopes
       if (initial_abundance(isot) .gt. 0.d0 .or. abundance(isot) .gt. 1.d-95) then
         print*, z(isot), z(isot)+n(isot), '  ',&
               element_names(z(isot)), level(isot), abundance(isot), initial_abundance(isot)
       end if  
    end do  

end


!*****************************************
!*                 TEST                  *
!*****************************************
! Test function for F2Py
character*256 function get_float(init_fileee, fl_arr, NN)

    use iso

    character*256 init_fileee
    integer NN
    double precision fl_arr(NN)

!    print*, init_fileee
!    print*, init_file

    print*, fl_arr(1)

    get_float = init_fileee
    return

end
!end function get_float

!*****************************************
!*                GEFSUB                 *
!*****************************************
! Taken from GEFSUB.FOR file
subroutine GEFSUB(P_Z_CN,P_A_CN,P_E_EXC,P_J_CN)                

      use GEFSUB_FOR
      use GEFSUBdcl2
      use RESULTS
      implicit none

      INTEGER*4 P_Z_CN
      INTEGER*4 P_A_CN
      REAL*4 P_E_EXC
      REAL*4 P_J_CN

!     /' Input parameters: '/ 
!     /' Atomic number,mass number,excitation energy/MeV,spin/h_bar of CN '/ 
!     /' Results are stored in external arrays. '/ 
!     ' 
      !INCLUDE "GEFSUBdcl2.FOR" 

      ! Benoit, changed U_Delta_S0 by get_U_delta_S0
      xDelta_S0 = get_U_Delta_S0(P_Z_CN,P_A_CN) 
!     ' default values 
!     ' 
!     ' Use nominal parameter values: 
      P_DZ_Mean_S1 = xP_DZ_Mean_S1 
      P_DZ_Mean_S2 = xP_DZ_Mean_S2 
      P_DZ_Mean_S3 = xP_DZ_Mean_S3 
      P_DZ_Mean_S4 = xP_DZ_Mean_S4 
      P_Z_Curv_S1 = xP_Z_Curv_S1 
      P_Z_Curv_S2 = xP_Z_Curv_S2 
      P_A_Width_S2 = xP_A_Width_S2 
      P_Z_Curv_S3 = xP_Z_Curv_S3 
      P_Z_Curv_S4 = xP_Z_Curv_S4 
      Delta_S0 = xDelta_S0 
      P_Shell_S1 = xP_Shell_S1 
      P_Shell_S2 = xP_Shell_S2 
      P_Shell_S3 = xP_Shell_S3 
      P_Shell_S4 = xP_Shell_S4 
      T_low_S1 = xT_low_S1 
      T_low_S2 = xT_low_S2 
      T_low_S3 = xT_low_S3 
      T_low_S4 = xT_low_S4 
      T_low_SL = xT_low_SL 
      P_att_pol = xP_att_pol 
      HOMPOL = xHOMPOL 
      POLARadd = xPOLARadd 
!     ' 
      R_E_exc_used = P_E_exc 
      I_A_CN = P_A_CN 
      I_Z_CN = P_Z_CN 
!     ' 
!     /' Central Z values of fission modes '/ 
!     ' 
!     /' Fit to positions of fission channels (Boeckstiegel et al.,2008) '/ 
!     /' P_DZ_Mean_S1 and P_DZ_Mean_S2 allow for slight adjustments '/ 
!     '    Scope 
      R_Z_mod = I_Z_CN 
      ZC_Mode_0 = R_Z_mod * 0.5E0 
!     /' Central Z value of SL mode '/ 
      ZC_Mode_1 = (53.0E0 - 51.5E0) / (1.56E0 - 1.50E0) * &
          (R_Z_mod**1.3E0 / I_A_CN - 1.50E0) + 51.5E0 + P_DZ_Mean_S1 
      ZC_Mode_2 = (55.8E0 - 54.5E0) / (1.56E0 - 1.50E0) * &
          (R_Z_mod**1.3E0 / I_A_CN - 1.50E0) + 54.5E0 + P_DZ_Mean_S2 
      ZC_Mode_3 = ZC_Mode_2 + 4.5E0 + P_DZ_Mean_S3 
!     '  ZC_Mode_4 = 38.5 + P_DZ_Mean_S4  ' structure in nuclei with A around 190 for 201Tl 
!     '  ZC_Mode_4 = 35.5 + P_DZ_Mean_S4  ' for 180Hg  ( 36.2 for 208Po ) 
!     ' 
!     ' Do not delete these lines (,because this is a very good fit!): 
!     '    ZC_Mode_4 = 38.5 + (I_A_CN-I_Z_CN-110)*0.12 - (I_A_CN-I_Z_CN-110)^2 * 0.009   '                - (I_Z_CN-77)*0.34 + P_DZ_Mean_S4 
!     ' 
      ZC_Mode_4 = 38.5 + (I_A_CN-I_Z_CN-110)*0.12 - (I_A_CN-I_Z_CN-110)**2 * &
          0.009 - (I_Z_CN-77)*0.34 + P_DZ_Mean_S4 
!     ' assumption: mode position moves with Z and A (adjusted to exp. data 
!     ' of Itkis and Andreyev et al. 
!     ' 
      ZC_Mode_4L = 42.05 
!     ' enhances S1 ' 
!     '    End Scope 
!     ' 
!     ' 
!     ' 
      I_N_CN = I_A_CN - I_Z_CN 
!     /' Mean deformation at scission as a function of mass '/ 
!     ' 
!     /' Mode 0: liquid drop and mode 4: Z = 38 '/ 
      beta1_prev = 0.3 
      beta2_prev = 0.3 
      beta1_opt = beta1_prev 
      beta2_opt = beta2_prev 
      DO I = 10 , I_Z_CN - 10 
      IZ1 = I 
      Z1 = REAL(IZ1) 
      IZ2 = I_Z_CN - IZ1 
      Z2 = REAL(IZ2) 
      A1 = Z1 / REAL(I_Z_CN) * REAL(I_A_CN) 
      A2 = I_A_CN - A1 
!     ' 
      CALL Beta_Equi(A1,A2,Z1,Z2,dneck,beta1_prev,beta2_prev,beta1_opt,beta2_opt) 
!     ' 
!     'Print "Mode 0,Z1,Z2,beta1,beta2 ";Z1;" ";Z2;" ";beta1_opt,beta2_opt 
!     'Print Z1;" ";Z2;" ";beta1_opt,beta2_opt 
      Beta(0,1,IZ1) = beta1_opt 
!     /' "light" fragment '/ 
      Beta(4,1,IZ1) = beta1_opt 
      Beta(0,2,IZ2) = beta2_opt 
!     /' "heavy" fragment '/ 
      Beta(4,2,IZ2) = beta2_opt 
      beta1_prev = beta1_opt 
      beta2_prev = beta2_opt 
      E_defo = get_LyMass(Z1,A1,beta1_opt) - get_LyMass(Z1,A1,0.0) 
      Edefo(0,1,IZ1) = E_defo 
!     /' "light" fragment '/ 
      Edefo(4,1,IZ1) = E_defo 
      E_defo = get_LyMass(Z2,A2,beta2_opt) - get_LyMass(Z2,A2,0.0) 
      Edefo(0,2,IZ2) = E_defo 
!     /' "heavy" fragment '/ 
      Edefo(4,2,IZ2) = E_defo 
      END DO 
!     ' 
!     /' Mode 1: deformed shells (light) and spherical (heavy) '/ 
      DO I = 10 , I_Z_CN - 10 
      Z1 = I 
      Z2 = I_Z_CN - Z1 
      A1 = (Z1 - 0.5E0) / REAL(I_Z_CN) * REAL(I_A_CN) 
!     /' polarization roughly considered '/ 
      A2 = I_A_CN - A1 
      IF (  I_Z_CN * 0.5 .LT. ZC_Mode_1  ) THEN 
!     ' Beta_opt_light(A1,A2,Z1,Z2,dneck,0,rbeta_ld) 
!     /' nu_mean of Cf requires shells in the light fragment: '/ 
      rbeta = get_beta_light(I,betaL0,betaL1) - 0.1 
!     ' smaller than general deformation of light fragment 
!     '        (less neck influence due to spherical heavy fragment) 
      IF (  rbeta .LT. 0  )  rbeta = 0 
      Else 
      rbeta = get_beta_heavy(I,betaH0,betaH1) 
!     ' equal to S2 channel 
      IF (  rbeta .LT. 0  )  rbeta = 0 
      End If 
      Beta(1,1,I) = rbeta 
!     /' "light" fragment '/ 
      E_defo = get_LyMass(Z1,A1,rbeta) - get_LyMass(Z1,A1,0.0) 
      Edefo(1,1,I) = E_defo 
!     /' "light" fragment '/ 
      END DO 
!     ' 
      DO I = 10 , I_Z_CN - 10 
      rbeta = 0 
      Beta(1,2,I) = rbeta 
      Edefo(1,2,I) = 0 
!     /' "heavy" fragment (at S1 shell) '/ 
      END DO 
!     ' 
!     /' Mode 2: deformed shells (light and heavy) '/ 
      DO I = 10 , I_Z_CN - 10 
      Z1 = I 
      Z2 = I_Z_CN - Z1 
      A1 = (Z1 - 0.5E0) / REAL(I_Z_CN) * REAL(I_A_CN) 
!     /' polarization roughly considered '/ 
      A2 = I_A_CN - A1 
      IF (  I_Z_CN * 0.5 .LT. ZC_Mode_2  ) THEN 
!     ' Beta_opt_light(A1,A2,Z1,Z2,dneck,beta_heavy(Z2),rbeta_ld) 
      rbeta = get_beta_light(I,betaL0,betaL1) 
!     ' general deformation of light fragment 
      IF (  rbeta .LT. 0  )  rbeta = 0 
!     ' negative values replaced by 0 
      Else 
      rbeta = get_beta_heavy(I,betaH0,betaH1) 
!     ' equal to S2 channel 
      End If 
      Beta(2,1,I) = rbeta 
      E_defo = get_LyMass(Z1,A1,rbeta) - get_LyMass(Z1,A1,0.0) 
      Edefo(2,1,I) = E_defo 
      END DO 
      DO I = 10 , I_Z_CN - 10 
      rbeta = get_beta_heavy(I,betaH0,betaH1) 
!     /' "heavy" fragment (at S2 shell)'/ 
      IF (  rbeta .LT. 0  )  rbeta = 0 
!     ' negative values replaced by 0 
      Beta(2,2,I) = rbeta 
      Z1 = I 
      A1 = (Z1 + 0.5E0) / I_Z_CN * I_A_CN 
!     /' polarization roughly considered '/ 
      E_defo = get_LyMass(Z1,A1,rbeta) - get_LyMass(Z1,A1,0.0) 
      Edefo(2,2,I) = E_defo 
      END DO 
!     ' 
!     /' Mode 3 '/ 
      DO I = 10 , I_Z_CN - 10 
      Z1 = I 
      Z2 = I_Z_CN - Z1 
      A1 = (Z1 - 0.5E0) / REAL(I_Z_CN) * REAL(I_A_CN) 
!     /' polarization roughly considered '/ 
      A2 = I_A_CN - A1 
      rbeta = get_beta_light(I,betaL0,betaL1) 
      rbeta = Max(rbeta-0.10,0.0) 
!     /' for low nu-bar of lightest fragments '/ 
!     '  Beta_opt_light(A1,A2,Z1,Z2,dneck,beta_heavy(Z2,betaH0,betaH1),rbeta) 
      Beta(3,1,I) = rbeta 
      E_defo = get_LyMass(Z1,A1,rbeta) - get_LyMass(Z1,A1,0.0) 
      Edefo(3,1,I) = E_defo 
      END DO 
      DO I = 10 , I_Z_CN - 10 
      rbeta = get_beta_heavy(I,betaH0,betaH1) + 0.2 
!     /' for high nu-bar of heaviest fragments '/ 
      IF (  rbeta .LT. 0  )  rbeta = 0 
      Beta(3,2,I) = rbeta 
      Z1 = I 
      A1 = (Z1 + 0.5E0) / REAL(I_Z_CN) * REAL(I_A_CN) 
!     /' polarization roughly considered '/ 
      E_defo = get_LyMass(Z1,A1,rbeta) - get_LyMass(Z1,A1,0.0) 
      Edefo(3,2,I) = E_defo 
      END DO 
!     ' 
!     /' Mode 5: (Channel ST1 in both fragments) '/ 
      DO I = 10 , I_Z_CN - 10 
      Z1 = I 
      Z2 = I_Z_CN - Z1 
      rbeta = Beta(1,2,I) 
      IF (  rbeta .LT. 0  )  rbeta = 0 
      Beta(5,1,Int(Z1)) = rbeta 
      Beta(5,2,Int(Z1)) = rbeta 
      END DO 
!     ' 
!     /' Mode 6: (Channel ST2 in both fragments) '/ 
      DO I = 10 , I_Z_CN - 10 
      Z1 = I 
      Z2 = I_Z_CN - Z1 
      rbeta = Beta(2,2,I) 
      IF (  rbeta .LT. 0  )  rbeta = 0 
      Beta(6,1,Int(Z1)) = rbeta 
      Beta(6,2,Int(Z1)) = rbeta 
      END DO 
!     ' 
!     ' 
!     /' Mean Z as a function of mass '/ 
!     ' 
!     /' Mode 0 '/ 
      DO I = 10 , I_A_CN - 10 
      ZUCD = REAL(I) / REAL(I_A_CN) * REAL(I_Z_CN) 
      beta1 = Beta(0,1,Int(ZUCD + 0.5)) 
      beta2 = Beta(0,2,Int(I_Z_CN - ZUCD + 0.5)) 
      Z1 = get_Z_equi(I_Z_CN,I,I_A_CN - I,beta1,beta2,dneck,0,0.0,1.0) 
      Zmean(0,1,I) = Z1 
      Zshift(0,1,I) = Z1 - ZUCD 
      Zmean(0,2,I_A_CN - I) = I_Z_CN - Z1 
      Zshift(0,2,I_A_CN - I) = ZUCD - Z1 
      END DO 
!     ' 
!     /' Mode 1 '/ 
      DO I = 10 , I_A_CN - 10 
      ZUCD = REAL(I) / REAL(I_A_CN) * REAL(I_Z_CN) 
      Z = ZUCD + ZPOL1 
!     /' Charge polarisation is considered in a crude way '/ 
      beta1 = Beta(1,1,NINT(Z)) 
!     /' "light" fragment '/ 
      Z = ZUCD - ZPOL1 
      beta2 = Beta(1,2,NINT(I_Z_CN-Z)) 
!     /' "heavy" fragment  at S1 shell '/ 
      IF (  REAL(I_Z_CN) * 0.5 .LT. ZC_Mode_1  ) THEN 
      Z1 = get_Z_equi(I_Z_CN,I,I_A_CN - I,beta1,beta2,dneck,1,POLARadd,POLARfac) 
      Else 
      Z1 = get_Z_equi(I_Z_CN,I,I_A_CN - I,beta1,beta2,dneck,1,0.0,0.0) 
      End If 
      Z1 = Z1 + ZPOL1 
!     /' Charge polarization by shell '/ 
!     ' 
      IF (  I_Z_CN - Z1 .LT. 50 .AND. (I_Z_CN - Z1) .GT. Z1  ) THEN 
      Z1 = I_Z_CN - 50 
!     /' Z of mean heavy fragment not below 50 '/ 
      END IF 
!     ' 
      Zmean(1,1,I) = Z1 
      Zshift(1,1,I) = Z1 - ZUCD 
!     ' neutron-deficient 
      Zmean(1,2,I_A_CN - I) = I_Z_CN - Z1 
      Zshift(1,2,I_A_CN - I) = ZUCD - Z1 
!     ' neutron rich at shell 
      END DO 
!     ' 
!     /' Mode 2 '/ 
      DO I = 10 , I_A_CN - 10 
      ZUCD = REAL(I) / REAL(I_A_CN) * REAL(I_Z_CN) 
      Z = ZUCD 
!     /' Charge polarisation is here neglected '/ 
      beta1 = Beta(2,1,NINT(Z)) 
      beta2 = Beta(2,2,NINT(I_Z_CN-Z)) 
      IF (  REAL(I_Z_CN) * 0.5 .LT. ZC_Mode_2  ) THEN 
      Z1 = get_Z_equi(I_Z_CN,I,I_A_CN-I,beta1,beta2,dneck,2,POLARadd,POLARfac) 
      Else 
      Z1 = get_Z_equi(I_Z_CN,I,I_A_CN-I,beta1,beta2,dneck,2,0.0,0.0) 
      End If 
!     ' 
      Zmean(2,1,I) = Z1 
      Zshift(2,1,I) = Z1 - ZUCD 
!     ' neutron deficieint 
      Zmean(2,2,I_A_CN - I) = I_Z_CN - Z1 
      Zshift(2,2,I_A_CN - I) = ZUCD - Z1 
!     ' neutron rich at shell 
      END DO 
!     ' 
!     /' Mode 3 '/ 
      DO I = 10 , I_A_CN - 10 
      ZUCD = REAL(I) / REAL(I_A_CN) * REAL(I_Z_CN) 
      Z = ZUCD 
!     /' Charge polarisation is here neglected '/ 
      beta1 = Beta(3,1,NINT(Z)) 
      beta2 = Beta(3,2,NINT(I_Z_CN-Z)) 
      Z1 = get_Z_equi(I_Z_CN,I,I_A_CN - I,beta1,beta2,dneck,3,POLARadd,POLARfac) 
      Zmean(3,1,I) = Z1 
      Zshift(3,1,I) = Z1 - ZUCD 
      Zmean(3,2,I_A_CN - I) = I_Z_CN - Z1 
      Zshift(3,2,I_A_CN - I) = ZUCD - Z1 
      END DO 
!     ' 
!     /' Mode 4 (assumed to be equal to mode 0) '/ 
      DO I = 10 , I_A_CN - 10 
      Zmean(4,1,I) = Zmean(0,1,I) 
      Zshift(4,1,I) = Zshift(0,1,I) 
      Zmean(4,2,I_A_CN - I) = Zmean(0,2,I_A_CN - I) 
      Zshift(4,2,I_A_CN - I) = Zshift(0,2,I_A_CN - I) 
      END DO 
!     ' 
!     ' 
!     /' General relations between Z and A of fission channels '/ 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_0 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(0,2,NINT(RA)) 
      END DO 
      AC_Mode_0 = (ZC_Mode_0 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
!     /' mean position in mass '/ 
      NC_Mode_0 = AC_Mode_0 - ZC_Mode_0 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_1 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(1,2,NINT(RA)) 
      END DO 
      AC_Mode_1 = (ZC_Mode_1 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_1 = AC_Mode_1 - ZC_Mode_1 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_2 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(2,2,NINT(RA)) 
      END DO 
      AC_Mode_2 = (ZC_Mode_2 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_2 = AC_Mode_2 - ZC_Mode_2 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_3 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(3,2,NINT(RA)) 
      END DO 
      AC_Mode_3 = (ZC_Mode_3 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_3 = AC_Mode_3 - ZC_Mode_3 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_4 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(4,2,NINT(RA)) 
      END DO 
      AC_Mode_4 = (ZC_Mode_4 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_4 = AC_Mode_4 - ZC_Mode_4 
!     ' 
!     ' 
!     /' Potential curvatures of fission modes '/ 
!     ' 
!     ' For the width of the mass distribution (potential between saddle and scission): 
      R_Z_Curv_S0 = 8.E0 / REAL(I_Z_CN)**2 * get_Masscurv(REAL(I_Z_CN),REAL(I_A_CN),P_I_rms_CN,kappa) 
!     ' For the yields of the fission channels (potential near saddle): 
      R_Z_Curv1_S0 = 8.E0 / REAL(I_Z_CN)**2 * get_Masscurv1(REAL(I_Z_CN),REAL(I_A_CN),0.0,kappa) 
      R_A_Curv1_S0 = 8.E0 / REAL(I_A_CN)**2 * get_Masscurv1(REAL(I_Z_CN),REAL(I_A_CN),0.0,kappa) 
!     ' 
!     ' 
!     /' Energy transformation '/ 
!     ' 
      Select CASE( Emode) 
      CASE( 0) 
!     ' Energy above outer barrier given 
      R_E_exc_Eb = R_E_exc_used 
      R_E_exc_GS = R_E_exc_used + get_BFTFB(REAL(I_Z_CN),REAL(I_A_CN),1) 
      CASE( 1,3,-1) 
!     ' Energy above ground state given 
      R_E_exc_Eb = R_E_exc_used - get_BFTFB(REAL(I_Z_CN),REAL(I_A_CN),1) 
      R_E_exc_GS = R_E_exc_used 
      CASE( 2) 
!     ' kinetic energy of neutron given 
      SN = (get_U_MASS(REAL(I_Z_CN),REAL(I_A_CN-1)) + &
           get_LyPair(I_Z_CN,I_A_CN-1)) - &
           (get_U_MASS(REAL(I_Z_CN),REAL(I_A_CN)) + &
             get_LyPair(I_Z_CN,I_A_CN))
      R_E_exc_GS = R_E_exc_used + SN 
      R_E_exc_Eb = R_E_exc_GS - get_BFTFB(REAL(I_Z_CN),REAL(I_A_CN),1) 
      End Select 
!     ' 
!     ' 
!     /' Fission barriers -> global parameters '/ 
!     ' 
      B_F = get_BFTF(REAL(I_Z_CN),REAL(I_A_CN),1) 
      B_F_ld = get_BFTF(REAL(I_Z_CN),REAL(I_A_CN),0) 
      E_B = get_BFTFB(REAL(I_Z_CN),REAL(I_A_CN),1) 
      E_B_ld = get_BFTFB(REAL(I_Z_CN),REAL(I_A_CN),0) 
!     ' 
!     ' 
!     /' Barriers and excitation energies of the fission modes '/ 
!     ' 
      E_exc_S0_prov = R_E_exc_Eb 
!     ' 
!     ' 
!     /' Additional influence of N=82 assumed '/ 
      Delta_NZ_Pol = 82.E0/50.E0 - REAL(I_N_CN)/REAL(I_Z_CN) 
      R_Shell_S1_eff = P_Shell_S1 * (1.E0 - P_Att_Pol * Abs(Delta_NZ_Pol)) 
!     ' 
!     ' 
!     /' In Pu,the Z=50 shell meets Z=44 in the light fragment. '/ 
!     /' A deformed shell at Z=44 is assumed to explain the enhancement        of the S1 channel around Pu '/ 
!     /' This very same shell automatically produces the double-humped '/ 
!     /' mass distribution in 180Hg '/ 
      S1_enhance = P_Shell_SL4 + (REAL(I_Z_CN) - ZC_Mode_1 - ZC_Mode_4L)**2 * P_Z_Curv_SL4 
!     'Print "ZC_Mode_1,ZC_Mode_4",ZC_Mode_1,ZC_Mode_4 
!     'Print "Delta-Z S1-S4,S1_enhance",I_Z_CN-ZC_Mode_1 - ZC_Mode_4,S1_enhance 
      IF (  S1_enhance .GT. 0  )  S1_enhance = 0 
      R_Shell_S1_eff = R_Shell_S1_eff + S1_enhance 
!     ' 
!     /' The high TKE of S1 in 242Pu(sf) (and neighbours) is obtained by assuming '/ 
!     /' that the Z=44 shell reduces the deformation of the light fragment. '/ 
      DO I = 10 , I_Z_CN - 10 
      Z1 = I 
      A1 = (Z1 - 0.5E0) / REAL(I_Z_CN) * REAL(I_A_CN) 
!     /' polarization roughly considered '/ 
!     '      Beta(1,1,Z1) = Beta(1,1,Z1) + 0.15 * S1_enhance   /' "light" fragment '/ 
      Beta(1,1,I) = exp(S1_enhance) * Beta(1,1,I) + (1.E0-exp(S1_enhance)) * (Beta(1,1,I)-0.25) 
      Beta(1,1,I) = Max(Beta(1,1,I),0.0) 
      E_defo = get_LyMass(Z1,A1,Beta(1,1,I)) - get_LyMass(Z1,A1,0.0) 
      Edefo(1,1,I) = E_defo 
!     /' "light" fragment '/ 
      END DO 
!     ' 
!     ' Influence of S2 shell in complementary fragment 
!     ' May be called "S12 fission channel" 
      T_Asym_Mode_2 = 0.5 
      SigZ_Mode_2 = SQRT(0.5E0 * T_Asym_Mode_2/(P_Z_Curv_S2)) 
      SigA_Mode_2 = SigZ_Mode_2 * REAL(I_A_CN) / REAL(I_Z_CN) 
      S1_enhance = P_Shell_S2 * &
          get_U_Box(REAL(P_A_CN) - AC_Mode_2 - AC_Mode_1,SigA_Mode_2,P_A_Width_S2) *P_A_Width_S2 
      IF (  S1_enhance .LT. 0.01  ) THEN 
      R_Shell_S1_eff = R_Shell_S1_eff + S1_enhance 
      End If 
!     ' Modify deformation of complementary fragment in corresponding analyzer 
!     ' 
!     ' Overlap of S2 and shell in light fragment 
      R_Shell_S2_eff = P_Shell_S2 
!     '   S2_enhance = P_Shell_S4 +  '             (Csng(I_Z_CN) - ZC_Mode_2 - ZC_Mode_4)^2 * P_Z_Curv_S4 
!     '   If S2_enhance > 0 Then S2_enhance = 0 
!     '   R_Shell_S2_eff = R_Shell_S2_eff + S2_enhance 
!     ' 
!     ' Overlap of S3 and shell in light fragment 
      R_Shell_S3_eff = P_Shell_S3 * (1.E0 - PZ_S3_olap_curv * &
          (REAL(I_Z_CN) - ZC_Mode_3 - PZ_S3_olap_pos)**2) 
!     '        * (Csng(I_Z_CN) - 60.5E0 - PZ_S3_olap_pos)^2) 
      R_Shell_S3_eff = Min(R_Shell_S3_eff,0.0) 
!     ' 
!     '   R_Shell_S4_eff = 2.0 * (P_Shell_S4 + P_Z_Curv_S4*(ZC_Mode_4 - ZC_Mode_0)^2) 
      R_Shell_S4_eff = 2.0 * (P_Shell_S4 + P_Z_Curv_S4 * (ZC_Mode_4 - ZC_Mode_0)**2) 
!     ' overlap of S4 in both fragments 
      IF (  R_Shell_S4_eff .GT. P_Shell_S4  )  R_Shell_S4_eff = P_Shell_S4 
!     ' no overlap at large distance 
!     ' 
      E_ld_S1 = R_A_Curv1_S0 * (REAL(I_A_CN)/REAL(I_Z_CN)*(ZC_MODE_1 - ZC_MODE_0) )**2 
      B_S1 = E_ld_S1 + R_Shell_S1_eff 
      E_exc_S1_prov = E_Exc_S0_prov - B_S1 
!     ' 
      E_ld_S2 = R_A_Curv1_S0 * (REAL(I_A_CN)/REAL(I_Z_CN)*(ZC_MODE_2 - ZC_MODE_0) )**2 
      B_S2 = E_ld_S2 + R_Shell_S2_eff 
      E_exc_S2_prov = E_Exc_S0_prov - B_S2 
!     ' 
      E_ld_S3 = R_A_Curv1_S0 * (REAL(I_A_CN)/REAL(I_Z_CN)*(ZC_MODE_3 - ZC_MODE_0) )**2 
      B_S3 = E_ld_S3 + R_Shell_S3_eff 
      E_exc_S3_prov = E_Exc_S0_prov - B_S3 
!     ' 
      IF (  I_A_CN .LT. 220  ) THEN 
!     ' Only here S4 is close enough to symmetry to have a chance 
      E_ld_S4 = R_A_Curv1_S0 * (REAL(I_A_CN)/REAL(I_Z_CN)*(ZC_MODE_4 - ZC_MODE_0) )**2 
      B_S4 = E_ld_S4 + R_Shell_S4_eff 
      E_exc_S4_prov = E_Exc_S0_prov - B_S4 
      Else 
      B_S4 = 9999 
      E_exc_S4_prov = - 9999 
      End If 
!     ' 
!     /' Mode 11 (overlap of channel 1 in light and heavy fragment '/ 
!     /' Potential depth with respect to liquid-drop potential: B_S11 '/ 
      B_S11 = 2.E0 * (R_Shell_S1_eff + De_Defo_S1 + P_Z_Curv_S1 * &
          (ZC_Mode_1 - ZC_Mode_0)**2 ) - De_Defo_S1 
!     ' 
!     ' 
!     /' Lowering of effective barrier by lower ZPM due to larger width in 
!     partial overlap region (shells in light and heavy fragment) '/ 
      DES11ZPM = Level_S11 * Min(Abs(ZC_Mode_1 - ZC_Mode_0),4.E0*P_Z_Curv_S1) 
!     ' Print B_S11,DES11ZPM,ZC_Mode_1-ZC_Mode_0 
!     ' 
      B_S11 = B_S11 + DES11ZPM 
!     ' 
!     '  If B_S11 > R_Shell_S1_eff + 0.5E0 Then 
!     '   If B_S11 > R_Shell_S1_eff + Level_S11 Then 
!     '     B_S11 = 100   ' S1 and S11 are exclusive 
!     '   Else 
!     '     B_S11 = Min(B_S11,R_Shell_S1_eff) 
!     '   End If 
!     ' 
!     ' 
      E_exc_S11_prov = E_Exc_S0_prov - B_S11 
!     ' 
!     /' Mode 22 (overlap of channel 2 in light and heavy fragment '/ 
!     /' Potential depth with respect to liquid-drop potential: B_S22 '/ 
!     ' 
!     '   B_S22 = 2.E0 * (E_ld_S2 + P_Shell_S2)  '       + 2.E0 * P_Z_Curv_S2 * (ZC_Mode_2 - ZC_Mode_0)^2   /' Parabola '/ 
!     'Print E_ld_S2,P_Shell_S2,P_Z_Curv_S2,ZC_Mode_2,ZC_Mode_0 
      B_S22 = 2.E0 * R_Shell_S2_eff * &
          get_U_Box(REAL(P_A_CN)/2.0 - AC_Mode_2,SigA_Mode_2,P_A_Width_S2) * P_A_Width_S2 
!     ' The integral of U_Box is normalized,not the height! 
!     '    If Abs((P_A_CN/2.E0) - AC_Mode_2) > P_A_Width_S2 Then B_S22 = 9999 
      IF (  P_A_CN .LT. 226  )  B_S22 = 9999 
!     ' 
      E_exc_S22_prov = E_Exc_S0_prov - B_S22 
!     ' 
!     ' 
      E_Min_Barr = Min(0.0,B_S1) 
      E_Min_Barr = Min(E_Min_Barr,B_S2) 
      E_Min_Barr = Min(E_Min_Barr,B_S3) 
      E_Min_Barr = Min(E_Min_Barr,B_S4) 
      E_Min_Barr = Min(E_Min_Barr,B_S11) 
      E_Min_Barr = Min(E_Min_Barr,B_S22) 
!     ' 
!     /' Energy minus the height of the respective fission saddle '/ 
      E_exc_S0 = E_exc_S0_prov + E_Min_Barr - Delta_S0 
      E_exc_S1 = E_exc_S1_prov + E_Min_Barr 
      E_exc_S2 = E_exc_S2_prov + E_Min_Barr 
      E_exc_S3 = E_exc_S3_prov + E_Min_Barr 
      E_exc_S4 = E_exc_S4_prov + E_Min_Barr 
      E_exc_S11 = E_exc_S11_prov + E_Min_Barr 
      E_exc_S22 = E_exc_S22_prov + E_Min_Barr 
!     ' 
!     /' Energy above the lowest fission saddle '/ 
      E_exc_Barr = Max(E_Exc_S0,E_Exc_S1) 
      E_exc_Barr = Max(E_exc_Barr,E_Exc_S2) 
      E_exc_Barr = Max(E_exc_Barr,E_Exc_S3) 
      E_exc_Barr = Max(E_exc_Barr,E_Exc_S4) 
      E_exc_Barr = Max(E_exc_Barr,E_exc_S11) 
      E_exc_Barr = Max(E_exc_Barr,E_exc_S22) 
!     ' 
!     ' 
!     /' Collective temperature used for calculating the widths 
!     in mass asymmetry and charge polarization '/ 
!     ' 
      IF (  E_Exc_S0 .LT. 0  ) THEN 
      E_tunn = -E_Exc_S0
      Else
      E_tunn = 0
      END IF
      R_E_exc_eff = Max(0.1,E_Exc_S0) 
!     '  T_Coll_Mode_0 = TFCOLL * R_E_exc_eff + _  /' empirical,replaced by TRusanov '/ 
      T_Coll_Mode_0 = TCOLLFRAC * &
          (get_De_Saddle_Scission(REAL(I_Z_CN)**2 / REAL(I_A_CN)**0.33333E0,ESHIFTSASCI_coll) - &
           E_tunn) 
      T_Coll_Mode_0 = Max(T_Coll_Mode_0,0.0) 
!     ' 
!     ' 
!     /' Temperature description fitting to the empirical systematics of Rusanov et al. '/ 
!     /' Here from Ye. N. Gruzintsev et al.,Z. Phys. A 323 (1986) 307 '/ 
!     /' Empirical description of the nuclear temperature according to the '/ 
!     /' Fermi-gas description. Should be valid at higher excitation energies '/ 
      T_Rusanov = get_TRusanov(R_E_exc_eff,REAL(I_A_CN)) 
!     '  Print "Temperatures,(GEF,Total,Rusanov): ",T_Coll_Mode_0,TFCOLL * R_E_exc_eff,T_Rusanov 
      T_Coll_Mode_0 = Max(T_Coll_Mode_0,T_Rusanov) 
!     /' Transition vom const. temp. to Fermi gas occurs around 20 MeV by MAX function '/ 
!     '    T_Pol_Mode_0 = T_Pol_Red * T_Coll_Mode_0 
!     ' 
!     ' Application of the statistical model,intrinsic temperature at saddle 
      T_Pol_Mode_0 = get_U_Temp(0.5 * REAL(I_Z_CN),0.5 *REAL(I_A_CN),R_E_exc_eff,0,0,Tscale,Econd) 
      T_Asym_Mode_0 = SQRT(T_Coll_Mode_0**2 + (6E0*TCOLLMIN)**2) 
!     ' 
      E_pot_scission = (&
        get_De_Saddle_Scission(REAL(I_Z_CN)**2 / REAL(I_A_CN)**0.33333E0,ESHIFTSASCI_intr)-&
        E_tunn) 
!     ' 
!     ' 
      T_low_S1_used = T_low_S1 
!     ' 
      T_Coll_Mode_1 = TFCOLL * Max(E_exc_S1,0.E0) + &
          TCOLLFRAC * (get_De_Saddle_Scission(I_Z_CN**2 / &
          I_A_CN**0.33333E0,ESHIFTSASCI_coll) - E_tunn) 
      T_Coll_Mode_1 = Max(T_Coll_mode_1,0.0) 
!     '    T_Pol_Mode_1 = T_Pol_Red * T_Coll_Mode_1 
      T_Pol_Mode_1 = T_Pol_Mode_0 
      T_Asym_Mode_1 = SQRT(T_Coll_Mode_1**2 + (4.0*TCOLLMIN)**2) 
!     ' TCOLLMIN for ZPM 
!     ' 
      T_Coll_Mode_2 = TFCOLL * Max(E_exc_S2,0.E0) + &
          TCOLLFRAC * (get_De_Saddle_Scission(REAL(I_Z_CN)**2 / &
          REAL(I_A_CN)**0.33333E0,ESHIFTSASCI_coll) - E_tunn) 
      T_Coll_Mode_2 = Max(T_Coll_mode_2,0.0) 
!     '    T_Pol_Mode_2 = T_Pol_Red * T_Coll_Mode_2 
      T_Pol_Mode_2 = T_Pol_Mode_0 
      T_Asym_Mode_2 = SQRT(T_Coll_Mode_2**2 + TCOLLMIN**2) 
!     ' 
      T_Coll_Mode_3 = TFCOLL * Max(E_exc_S3,0.E0) + &
          TCOLLFRAC * (get_De_Saddle_Scission(REAL(I_Z_CN)**2 / &
          REAL(I_A_CN)**0.33333E0,ESHIFTSASCI_coll) - E_tunn) 
      T_Coll_Mode_3 = Max(T_Coll_mode_3,0.0) 
!     '    T_Pol_Mode_3 = T_Pol_Red * T_Coll_Mode_3 
      T_Pol_Mode_3 = T_Pol_Mode_0 
      T_Asym_Mode_3 = SQRT(T_Coll_Mode_3**2 + TCOLLMIN**2) 
!     ' 
      T_Coll_Mode_4 = TFCOLL * Max(E_exc_S4,0.E0) + &
          TCOLLFRAC * (get_De_Saddle_Scission(REAL(I_Z_CN)**2 / &
          REAL(I_A_CN)**0.33333E0,ESHIFTSASCI_coll) - E_tunn) 
      T_Coll_Mode_4 = Max(T_Coll_mode_4,0.0) 
!     '    T_Pol_Mode_4 = T_Pol_Red * T_Coll_Mode_4 
      T_Pol_Mode_4 = T_Pol_Mode_0 
      T_Asym_Mode_4 = SQRT(T_Coll_Mode_4**2 + 4.0*TCOLLMIN**2) 
!     ' ZPM like S1 
!     ' 
!     /' Stiffness in polarization '/ 
!     ' 
      RZ = REAL(I_Z_CN) * 0.5E0 
      RA = REAL(I_A_CN) * 0.5E0 
      beta1 = Beta(0,1,NINT(RZ)) 
      beta2 = Beta(0,2,NINT(RZ)) 
      R_Pol_Curv_S0 = ( get_LyMass( RZ - 1.E0,RA,beta1 ) + &
          get_LyMass( RZ + 1.0E0,RA,beta2 ) + get_LyMass( RZ + 1.0E0,RA,beta1 ) + &
          get_LyMass( RZ - 1.0E0,RA,beta2 ) + get_ECOUL( RZ - 1.0E0,RA,beta1,RZ + &
          1.0E0,RA,beta2,dneck) + &
          get_ECOUL( RZ + 1.0E0,RA,beta1,RZ - 1.0E0,RA,beta2,dneck) - &
          2.0E0*get_ECOUL( RZ,RA,beta1,RZ,RA,beta2,dneck) - &
          2.0E0*get_LyMass( RZ,RA,beta1 ) - &
          2.0E0*get_LyMass( RZ,RA,beta2)) * 0.5E0 
!     ' 
      P_Pol_Curv_S0 = R_Pol_Curv_S0 
!     ' 
      R_Pol_Curv_S1 = R_Pol_Curv_S0 
      R_Pol_Curv_S2 = R_Pol_Curv_S0 
      R_Pol_Curv_S3 = R_Pol_Curv_S0 
      R_Pol_Curv_S4 = R_Pol_Curv_S0 
!     ' 
!     ' 
!     ' 
!     /' Mean values and standard deviations of fission modes '/ 
!     ' 
      SIGZ_Mode_0 = SQRT(0.5E0 * T_Asym_Mode_0/R_Z_Curv_S0) 
      IF (  T_Pol_Mode_0 .GT. 1.E-2  ) THEN 
      SigPol_Mode_0 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S0 / Tanh(HOMPOL/(2.E0 * T_Pol_Mode_0))) 
      Else 
      SigPol_Mode_0 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S0) 
!     /' including influence of zero-point motion '/ 
      END IF 
!     ' 
!     ' 
      R_E_intr_S1 = Max(E_Exc_S1+get_LyPair(I_Z_CN,I_A_CN),0.0) 
      R_Att(1) = exp(-R_E_intr_S1/Shell_fading) 
      R_Att(5) = R_Att(1) 
      R_Att_Sad(1) = exp(-R_E_intr_S1/Shell_fading) 
      R_Att_Sad(5) = R_Att_Sad(1) 
      SIGZ_Mode_1 = SQRT(0.5E0 * T_Asym_Mode_1/(P_Z_Curv_S1*SQRT(R_Att(1)))) 
      IF (  T_Pol_Mode_1 .GT. 1.E-2  ) THEN 
      SigPol_Mode_1 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S1 / Tanh(HOMPOL/(2.E0 * T_Pol_Mode_1))) 
      Else 
      SigPol_Mode_1 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S1) 
      END IF 
!     ' 
      R_E_intr_S2 = Max(E_Exc_S2+get_LyPair(I_Z_CN,I_A_CN),0.0) 
      R_Att(2) = exp(-R_E_intr_S2/Shell_fading) 
      R_Att(6) = R_Att(2) 
      R_Att_Sad(2) = exp(-R_E_intr_S2/Shell_fading) 
      R_Att_Sad(6) = R_Att_Sad(2) 
      SIGZ_Mode_2 = SQRT(0.5E0 * T_Asym_Mode_2/(P_Z_Curv_S2*SQRT(R_Att(2)))) 
      IF (  T_Pol_Mode_2 .GT. 1.E-2  ) THEN 
      SigPol_Mode_2 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S2 / Tanh(HOMPOL/(2.E0 * T_Pol_Mode_2))) 
      Else 
      SigPol_Mode_2 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S2) 
      End If 
!     ' 
      R_E_intr_S3 = Max(E_exc_S3+get_LyPair(I_Z_CN,I_A_CN),0.0) 
      R_Att(3) = exp(-R_E_intr_S3/Shell_fading) 
      R_Att_Sad(3) = exp(-R_E_intr_S3/Shell_fading) 
      SIGZ_Mode_3 = SQRT(0.5E0 * T_Asym_Mode_3/(P_Z_Curv_S3*SQRT(R_Att(3)))) 
      IF (  T_Pol_Mode_3 .GT. 1.E-2  ) THEN 
      SigPol_Mode_3 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S3 / Tanh(HOMPOL/(2.E0 * T_Pol_Mode_3))) 
      Else 
      SigPol_Mode_3 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S3) 
      End if 
!     ' 
      R_E_intr_S4 = Max(E_exc_S4+get_LyPair(I_Z_CN,I_A_CN),0.0) 
      R_Att(4) = exp(-R_E_intr_S4/Shell_fading) 
      R_Att_Sad(4) = exp(-R_E_intr_S4/Shell_fading) 
      SIGZ_Mode_4 = SQRT(0.5E0 * T_Asym_Mode_4/(P_Z_Curv_S4*SQRT(R_Att(4)))) 
      IF (  T_Pol_Mode_4 .GT. 1.E-2  ) THEN 
      SigPol_Mode_4 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S4 / Tanh(HOMPOL/(2.E0 * T_Pol_Mode_4))) 
      Else 
      SigPol_Mode_4 = SQRT(0.25E0 * HOMPOL / R_Pol_Curv_S4) 
      End if 
!     ' 
!     ' 
!     ' 
!     /' Energy-dependent shift of fission channels '/ 
!     '    Scope 
      P_Z_Curv_S1_eff = P_Z_Curv_S1 * P_Z_Curvmod_S1 
      P_Z_Curv_S2_eff = P_Z_Curv_S2 * P_Z_Curvmod_S2 
      P_Z_Curv_S3_eff = P_Z_Curv_S3 * P_Z_Curvmod_S3 
      P_Z_Curv_S4_eff = P_Z_Curv_S4 * P_Z_Curvmod_S4 
!     ' 
      DZ_S1 = ZC_Mode_1 * (P_Z_Curv_S1_eff*R_Att(1) / &
          (R_Z_Curv_S0 + P_Z_Curv_S1_eff*R_Att(1)) - &
          (P_Z_Curv_S1_eff / (R_Z_Curv_S0 + P_Z_Curv_S1_eff) ) ) 
      DZ_S2 = ZC_Mode_2 * (P_Z_Curv_S2_eff*R_Att(2) / &
          (R_Z_Curv_S0 + P_Z_Curv_S2_eff*R_Att(2)) - &
          (P_Z_Curv_S2_eff / (R_Z_Curv_S0 + P_Z_Curv_S2_eff) ) ) 
      DZ_S3 = ZC_Mode_3 * (P_Z_Curv_S3_eff*R_Att(3) / &
          (R_Z_Curv_S0 + P_Z_Curv_S3_eff*R_Att(3)) - &
          (P_Z_Curv_S3_eff / (R_Z_Curv_S0 + P_Z_Curv_S3_eff) ) ) 
      DZ_S4 = SIGN(1.0,ZC_Mode_4 - ZC_Mode_0) * ZC_Mode_4 * &
          (P_Z_Curv_S4_eff*R_Att(4) / (R_Z_Curv_S0 + &
           P_Z_Curv_S4_eff*R_Att(4)) - (P_Z_Curv_S4_eff / &
          (R_Z_Curv_S0 + P_Z_Curv_S4_eff) ) ) 
!     ' 
!     ' Empirical shift of S2 channel at low excitation energy at scission 
!     ' for better reproduction of 238U(s,f) and some data for Th isotopes. 
!     ' Does not solve the problem of 229Th(nth,f). 
      EtotS2 = Max(E_Exc_S2 + EDISSFRAC * E_pot_scission,0.0) 
      IF (  EtotS2 .LT. 5.E0  ) THEN 
      DZ_S2 = DZ_S2 + (5.E0 - EtotS2) * 0.1 
      End If 
!     ' 
!     '   DZ_S1 = 0 
!     '   DZ_S2 = 0 
!     '   DZ_S3 = 0 
!     '   DZ_S4 = 0 
!     ' 
!     ' 
      P_Z_Mean_S0 = ZC_Mode_0 
      ZC_Mode_1 = ZC_Mode_1 + DZ_S1 
      P_Z_Mean_S1 = ZC_Mode_1 
!     /' Copy to global parameter '/ 
      ZC_Mode_2 = ZC_Mode_2 + DZ_S2 
      P_Z_Mean_S2 = ZC_Mode_2 
!     /'             "            '/ 
      ZC_Mode_3 = ZC_Mode_3 + DZ_S3 
      P_Z_Mean_S3 = ZC_Mode_3 
!     '   ZC_Mode_4 = ZC_Mode_4 + DZ_S4 
!     ' shift is very small,because S4 exists only close to symmetry 
      P_Z_Mean_S4 = ZC_Mode_4 
!     '    End Scope 
!     ' 
!     /' Energy dependence of charge polarization '/ 
!     /' Due to washing out of shells '/ 
!     ' 
      DO I = 10 , I_A_CN - 10 
!     ' mass number 
      DO J = 1 , 4 
!     ' fission channel 
      DO K = 1 , 2 
!     ' light - heavy group 
      Zshift(J,K,I) = Zshift(0,K,I) + (Zshift(J,K,I) - Zshift(0,K,I))*R_Att(J) 
      END DO 
      END DO 
      END DO 
!     ' 
!     ' 
!     /' Energy dependence of shell-induced deformation '/ 
!     /' Due to washing out of shells '/ 
!     /' (Under development) '/ 
!     /'For I = 10 To I_Z_CN - 10  ' mass number 
!     For J = 1 To 4           ' fission channel 
!     For K = 1 To 2         ' light - heavy group 
!     beta(J,K,I) = beta(0,K,I) + (beta(J,K,I) - beta(0,K,I))*R_Att_Sad(J) 
!     if beta(J,K,I) < 0 Then 
!     beta(J,K,I) = 0 
!     End If 
!     Z1 = I 
!     Z2 = I_Z_CN - Z1 
!     A1 = Z1 / Csng(I_Z_CN) * Csng(I_A_CN) 
!     A2 = I_A_CN - A1 
!     E_defo = Lymass(Z1,A1,beta(J,K,I)) - Lymass(Z1,A1,0.0) 
!     Edefo(J,K,I) = E_defo 
!     Next 
!     Next 
!     Next  '/ 
!     ' 
!     ' 
!     ' 
!     ' 
!     /' General relations between Z and A of fission channels '/ 
!     /' 2nd iteration '/ 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_0 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(0,2,NINT(RA)) 
      END DO 
      AC_Mode_0 = (ZC_Mode_0 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
!     /' mean position in mass '/ 
      NC_Mode_0 = AC_Mode_0 - ZC_Mode_0 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_1 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(1,2,NINT(RA)) 
      END DO 
      AC_Mode_1 = (ZC_Mode_1 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_1 = AC_Mode_1 - ZC_Mode_1 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_2 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(2,2,NINT(RA)) 
      END DO 
      AC_Mode_2 = (ZC_Mode_2 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_2 = AC_Mode_2 - ZC_Mode_2 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_3 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(3,2,NINT(RA)) 
      END DO 
      AC_Mode_3 = (ZC_Mode_3 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_3 = AC_Mode_3 - ZC_Mode_3 
!     ' 
      RZpol = 0 
      DO I = 1 , 3 
      RA = (ZC_Mode_4 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      RZpol = Zshift(4,2,NINT(RA)) 
      END DO 
      AC_Mode_4 = (ZC_Mode_4 - RZPol) * REAL(I_A_CN) / REAL(I_Z_CN) 
      NC_Mode_4 = AC_Mode_4 - ZC_Mode_4 
!     ' 
!     ' 
!     ' 
!     /' Yields of the fission modes '/ 
!     ' 
      Yield_Mode_0 = get_Getyield(E_exc_S0,E_exc_S0,T_low_SL,get_TEgidy(REAL(I_A_CN),0.E0,Tscale)) 
!     ' 
      Yield_Mode_1 = get_Getyield(E_exc_S1,E_exc_S0,T_low_S1_used,get_TEgidy(REAL(I_A_CN),R_Shell_S1_eff + dE_Defo_S1,Tscale)) 
!     ' 
      Yield_Mode_2 = get_Getyield(E_exc_S2,E_exc_S0,T_low_S2,get_TEgidy(REAL(I_A_CN),R_Shell_S2_eff + dE_Defo_S2,Tscale)) 
!     ' 
      Yield_Mode_3 = get_Getyield(E_exc_S3,E_exc_S0,T_low_S3,get_TEgidy(REAL(I_A_CN),R_Shell_S3_eff + dE_Defo_S3,Tscale)) 
!     ' 
      Yield_Mode_4 = get_Getyield(E_exc_S4,E_exc_S0,T_low_S4,get_TEgidy(REAL(I_A_CN),R_Shell_S4_eff + dE_Defo_S4,Tscale)) 
!     ' 
!     ' 
      IF (  B_S11 .GT. B_S1  ) THEN 
      Yield_Mode_11 = 0.0 
      Else 
      Yield_Mode_11 = get_Getyield(E_exc_S11,E_exc_S0,T_low_S11,get_TEgidy(REAL(I_A_CN),R_Shell_S1_eff + 2.E0 * dE_Defo_S1,Tscale)) 
      End If 
!     ' 
      IF (  B_S22 .GT. B_S2  ) THEN 
      Yield_Mode_22 = 0.0 
      Else 
      Yield_Mode_22 = get_Getyield(E_exc_S22,E_exc_S0,T_low_S2,get_TEgidy(REAL(I_A_CN),R_Shell_S2_eff,Tscale)) 
      End If 
!     ' 
!     ' 
      Yield_Norm = Yield_Mode_0 + Yield_Mode_1 + Yield_Mode_2 + Yield_Mode_3 + Yield_Mode_4 + Yield_Mode_11 + Yield_Mode_22 
      Yield_Mode_0 = Yield_Mode_0 / Yield_Norm 
      Yield_Mode_1 = Yield_Mode_1 / Yield_Norm 
      Yield_Mode_2 = Yield_Mode_2 / Yield_Norm 
      Yield_Mode_3 = Yield_Mode_3 / Yield_Norm 
      Yield_Mode_4 = Yield_Mode_4 / Yield_Norm 
      Yield_Mode_11 = Yield_Mode_11 / Yield_Norm 
      Yield_Mode_22 = Yield_Mode_22 / Yield_Norm 
!     ' 
!     ' 
!     /' Mass widhts of the fission channels '/ 
!     ' 
      SigA_Mode_0 = SigZ_Mode_0 * REAL(I_A_CN) / REAL(I_Z_CN) 
!     /' width in mass '/ 
      SigA_Mode_1 = SigZ_Mode_1 * REAL(I_A_CN) / REAL(I_Z_CN) 
      SigA_Mode_1 = Min(SigA_Mode_1,SigA_Mode_0) 
!     ' not broader than liquid-drop 
      SigA_Mode_2 = SigZ_Mode_2 * REAL(I_A_CN) / REAL(I_Z_CN) 
      SigA_Mode_2 = Min(SigA_Mode_2,SigA_Mode_0) 
!     ' not broader than liquid-drop 
      SigA_Mode_3 = SigZ_Mode_3 * REAL(I_A_CN) / REAL(I_Z_CN) 
      SigA_Mode_3 = Min(SigA_Mode_3,SigA_Mode_0) 
      SigA_Mode_4 = SigZ_mode_4 * REAL(I_A_CN) / REAL(I_Z_CN) 
      SigA_Mode_4 = Min(SigA_Mode_4,SigA_Mode_0) 
      SigA_Mode_11 = SigZ_Mode_1 * SQRT(2.E0) * REAL(I_A_CN) / REAL(I_Z_CN) 
      SigA_Mode_11 = Min(SigA_Mode_11,SigA_Mode_0) 
      SigA_Mode_22 = SigZ_Mode_2 * SQRT(2.E0) * REAL(I_A_CN) / REAL(I_Z_CN) 
      SigA_Mode_22 = Min(SigA_Mode_22,SigA_Mode_0) 
!     ' 
!     ' 
!     ' 
!     /' Shell effects of different fission channels '/ 
!     /' This is the "real" microscopic shell effect,not the effective shell-correction energy '/ 
!     /' EShell acts on the level density and determines the T parameter '/ 
!     ' 
      DO I = 1 , I_A_CN - 1 
      DO J = 0 , 4 
      EShell(J,1,I) = 0 
!     /' Shells in "light" fragment assumed to be zero '/ 
      END DO 
      DU0 = 0 
      EShell(0,2,I) = 0 
!     /' Shell = 0 in symmetric mode '/ 
      DU1 = R_Shell_S1_eff + dE_Defo_S1 
!     /' + R_A_Curv1_S1 * (AC_Mode_1 - Float(I,6))**2; '/ 
      DU1 = MIN(DU1,0.E0) 
!     /' Technical limit '/ 
      EShell(1,2,I) = DU1 
!     ' 
      DU2 = R_Shell_S2_eff + dE_Defo_S2 
!     /' + R_A_Curv1_S2 * (AC_Mode_2 - Float(I,6))**2; '/ 
      DU2 = Min(DU2,0.E0) 
!     /' Technical limit '/ 
      EShell(2,2,I) = DU2 
!     ' 
      DU3 = R_Shell_S3_eff + dE_Defo_S3 
!     /' + R_A_Curv1_S3 * (AC_Mode_3 - Float(I,6))**2; '/ 
      DU3 = Min(DU3,0.E0) 
!     /' Technical limit '/ 
      EShell(3,2,I) = DU3 
!     ' 
      DU4 = R_Shell_S4_eff + dE_Defo_S4 
!     /' + R_A_Curv1_S4 * (AC_Mode_4 - Float(I,6))**2; '/ 
      DU4 = Min(DU4,0.E0) 
!     /' Technical limit '/ 
      EShell(4,2,I) = DU4 
!     ' 
      END DO 
!     ' 
!     ' 
!     /' Intrinsic temperatures of fragments at scission '/ 
!     ' 
!     /' Mean values '/ 
      T_intr_Mode_0 = get_TEgidy(AC_Mode_0,0.0,0.8) 
      T_intr_Mode_1_heavy = get_TEgidy(AC_Mode_1,R_Shell_S1_eff + dE_Defo_S1,Tscale) 
      T_intr_Mode_1_light = get_TEgidy(REAL(I_A_CN) - AC_Mode_1,0.0,Tscale) 
      T_intr_Mode_2_heavy = get_TEgidy(AC_Mode_2,R_Shell_S2_eff + dE_Defo_S2,Tscale) 
      T_intr_Mode_2_light = get_TEgidy(REAL(I_A_CN) - AC_Mode_2,0.0,Tscale) 
      T_intr_Mode_3_heavy = get_TEgidy(AC_Mode_3,R_Shell_S3_eff + dE_Defo_S3,Tscale) 
      T_intr_Mode_3_light = get_TEgidy(REAL(I_A_CN) - AC_Mode_3,0.0,Tscale) 
      T_intr_Mode_4_heavy = get_TEgidy(AC_Mode_4,R_Shell_S4_eff + dE_Defo_S4,Tscale) 
      T_intr_Mode_4_light = get_TEgidy(REAL(I_A_CN) - AC_Mode_4,0.0,Tscale) 
!     ' 
!     ' 
!     /' Mass-dependent values of individual fragments '/ 
!     /' Mode 0 '/ 
      DO I = 1 , I_A_CN - 1 
      T = get_TEgidy(REAL(I),EShell(0,1,I),Tscale) 
      Temp(0,1,I) = T 
!     /' "light" fragment at freeze-out (somewhere before scission) '/ 
      T = get_TEgidy(REAL(I),EShell(0,2,I),Tscale) 
      Temp(0,2,I) = T 
!     /' "heavy" fragment at freeze-out (somewhere before scission) '/ 
!     ' 
      T = get_TEgidy(REAL(I),0.0,1.0) 
      TempFF(0,1,I) = T 
!     ' FF in their ground state 
      TempFF(0,2,I) = T 
!     ' FF in their ground state 
      END DO 
!     ' 
!     /' Mode 1 '/ 
      DO I = 1 , I_A_CN - 1 
      T = get_TEgidy(REAL(I),EShell(1,1,I),Tscale) 
      Temp(1,1,I) = T 
!     /' "light" fragment '/ 
      T = get_TEgidy(REAL(I),EShell(1,2,I),Tscale) 
      Temp(1,2,I) = T 
!     /' "heavy" fragment '/ 
!     ' 
      T = get_TEgidy(REAL(I),0.0,1.0) 
      TempFF(1,1,I) = T 
!     ' FF in their ground state 
      TempFF(1,2,I) = T 
!     ' FF in their ground state 
      END DO 
!     ' 
!     /' Mode 2 '/ 
      DO I = 1 , I_A_CN - 1 
      T = get_TEgidy(REAL(I),EShell(2,1,I),Tscale) 
      Temp(2,1,I) = T 
!     /' "light" fragment '/ 
      T = get_TEgidy(REAL(I),EShell(2,2,I),Tscale) 
      Temp(2,2,I) = T 
!     /' "heavy" fragment '/ 
!     ' 
!     /' The next section is introduced,because energy sorting is not strong enough, 
!     when shells are only introduced in the heavy fragment. 
!     Ad hoc assumption: For Mode 2 there are shells in both fragments of about 
!     equal size. Technically,we neglect the shells in both fragments. 
!     This has about the same effect for the energy sorting. '/ 
      T = get_TEgidy(REAL(I),0.0,Tscale) 
!     ' FF at scssion 
      Temp(2,1,I) = T 
!     /' "light" fragment '/ 
      T = get_TEgidy(REAL(I),0.0,Tscale) 
!     ' FF at scission 
      Temp(2,2,I) = T 
!     /' "heavy" fragment '/ 
!     ' 
      T = get_TEgidy(REAL(I),0.0,1.0) 
!     ' shell effect neglected 
      TempFF(2,1,I) = T 
!     ' FFs in their ground state 
      TempFF(2,2,I) = T 
!     ' FFs in their ground state 
      END DO 
!     ' 
!     /' Mode 3 '/ 
      DO I = 1 , I_A_CN -1 
      T = get_TEgidy(REAL(I),0.0,Tscale) 
      Temp(3,1,I) = T 
      T = get_TEgidy(REAL(I),0.0,Tscale) 
      Temp(3,2,I) = T 
!     ' 
      T = get_TEgidy(REAL(I),0.0,1.0) 
      TempFF(3,1,I) = T 
!     ' FF in their ground state 
      TempFF(3,2,I) = T 
!     ' FF in their ground state 
      END DO 
!     ' 
!     /' Mode 4 '/ 
      DO I = 1 , I_A_CN -1 
      T = get_TEgidy(REAL(I),0.0,Tscale) 
      Temp(4,1,I) = T 
      T = get_TEgidy(REAL(I),0.0,Tscale) 
      Temp(4,2,I) = T 
!     ' 
      T = get_TEgidy(REAL(I),0.0,1.0) 
      TempFF(4,1,I) = T 
!     ' FF in their ground state 
      TempFF(4,2,I) = T 
!     ' FF in their ground state 
      END DO 
!     ' 
!     ' 
!     /'** Intrinsic excitation energy at saddle and at scission as well as   **'/ 
!     /'** Even-odd effect in proton and neutron number for each fission mode **'/ 
      DO I_Mode = 0 , 6 
      E_coll_saddle(I_Mode) = 0 
      IF (  I_Mode .EQ. 0  )  Etot = E_exc_S0 
      IF (  I_Mode .EQ. 1  )  Etot = E_exc_S1 
      IF (  I_Mode .EQ. 2  )  Etot = E_exc_S2 
      IF (  I_Mode .EQ. 3  )  Etot = E_exc_S3 
      IF (  I_Mode .EQ. 4  )  Etot = E_exc_S4 
      IF (  I_Mode .EQ. 5  )  Etot = E_exc_S11 
      IF (  I_Mode .EQ. 6  )  Etot = E_exc_S22 
!     ' 
      IF (   MOD(I_Z_CN,2)  +  MOD(I_N_CN,2)  .EQ. 0  ) THEN 
!     /' Even-even CN '/ 
      IF (  Etot .GT. 0 .AND. Etot .LT. 2.E0 * 14.E0/SQRT(REAL(I_A_CN))  ) THEN 
      E_coll_saddle(I_Mode) = Etot 
      Etot = 0 
!     /' Excitation below the pairing gap in even-even CN goes into collective excitations '/ 
      End If 
      End If 
!     ' 
!     '    If I_Z_CN Mod 2 + I_N_CN Mod 2 = 0 Then    ' even-even 
!     '      Ediff = Min(Etot,14.0/sqr(Csng(I_A_CN))) 
!     '    End If 
!     '    If I_Z_CN Mod 2 + I_N_CN Mod 2 = 1 Then    ' even-odd or odd-even 
!     '       Ediff = Min(Etot,2.0 * 14.0/sqr(Csng(I_A_CN))) 
!     '    End If 
!     '    Ediff = Max(Ediff,0.0) 
!     '    Etot = Etot - Ediff 
!     ' 
!     ' 
      IF (  Etot .LT. 0  ) THEN 
      E_tunn = -Etot
      Else
      E_tunn = 0
      END IF
      Etot = Max(Etot,0.0) 
!     ' 
      E_pot_scission = (get_De_Saddle_Scission(REAL(I_Z_CN)**2 / REAL(I_A_CN)**0.33333E0,ESHIFTSASCI_intr) ) 
      Etot = Etot + EDISSFRAC * (E_pot_scission - E_tunn) 
!     /' All excitation energy at saddle and part of the potential-energy gain to scission 
!     go into intrinsic excitation energy at scission '/ 
!     ' 
!     ' 
!     ' 
      IF (  I_Mode .EQ. 2  ) THEN 
      EINTR_SCISSION = Etot 
!     /' (For Mode 2) Global parameter '/ 
      End If 
!     ' 
      DO IA1 = 40 , I_A_CN - 40 
!     ' 
      IA2 = I_A_CN - IA1 
      IF (  I_Mode .LE. 4  ) THEN 
      T1 = Temp(I_Mode,1,IA1) 
      T2 = Temp(I_Mode,2,IA2) 
      End If 
      IF (  I_Mode .EQ. 5  ) THEN 
      T1 = Temp(1,2,IA1) 
      T2 = Temp(1,2,IA2) 
      End If 
      IF (  I_Mode .EQ. 6  ) THEN 
      T1 = Temp(2,2,IA1) 
      T2 = Temp(2,2,IA2) 
      End If 
      DT = ABS(T2 - T1) 
!     ' 
!     /' Even-odd effect '/ 
      IF (   MOD(I_Z_CN,2)  .EQ. 0  ) THEN 
      Rincr1P = Exp(-Etot/PZ_EO_symm) 
      Else 
      Rincr1P = 0 
      End If 
      IF (   MOD(I_N_CN,2)  .EQ. 0  ) THEN 
      Rincr1N = Exp(-Etot/PN_EO_symm) 
      Else 
      Rincr1N = 0 
      End If 
      PEOZ(I_Mode,1,IA1) = Rincr1P 
      PEOZ(I_Mode,2,IA2) = Rincr1P 
      PEON(I_Mode,1,IA1) = Rincr1N 
      PEON(I_Mode,2,IA2) = Rincr1N 
!     ' 
      Rincr2 = get_Gaussintegral(DT/Etot-R_EO_Thresh,R_EO_Sigma*(DT+0.0001)) 
!     /' even-odd effect due to asymmetry '/ 
      Rincr2P = (R_EO_MAX - Rincr1P) * Rincr2 
      Rincr2N = (R_EO_MAX - Rincr1N) * Rincr2 
!     ' 
      IF (  IA1 .LT. IA2  ) THEN 
!     ' A1 is lighter 
      PEOZ(I_Mode,1,IA1) = PEOZ(I_Mode,1,IA1) + Rincr2P 
      IF (   MOD(I_Z_CN,2)  .EQ. 0  ) THEN 
      PEOZ(I_Mode,2,IA2) = PEOZ(I_Mode,2,IA2) + Rincr2P 
      Else 
      PEOZ(I_Mode,2,IA2) = PEOZ(I_Mode,2,IA2) - Rincr2P 
      End if 
      PEON(I_Mode,1,IA1) = PEON(I_Mode,1,IA1) + Rincr2N 
      IF (   MOD(I_N_CN,2)  .EQ. 0  ) THEN 
      PEON(I_Mode,2,IA2) = PEON(I_Mode,2,IA2) + Rincr2N 
      Else 
      PEON(I_Mode,2,IA2) = PEON(I_Mode,2,IA2) - Rincr2N 
      End if 
      Else 
      PEOZ(I_Mode,1,IA1) = PEOZ(I_Mode,2,IA1) 
      PEON(I_Mode,1,IA1) = PEON(I_Mode,2,IA1) 
      PEOZ(I_Mode,2,IA2) = PEOZ(I_Mode,1,IA2) 
      PEON(I_Mode,2,IA2) = PEON(I_Mode,1,IA2) 
      End If 
!     ' 
!     ' 
!     /'  Else 
!     PEOZ(I_Mode,2,IA2) =                PEOZ(I_Mode,1,IA2) + Rincr2P 
!     IF I_Z_CN Mod 2 = 0 Then 
!     PEOZ(I_Mode,1,IA1) =                 PEOZ(I_Mode,1,IA1) + Rincr2P 
!     Else 
!     PEOZ(I_Mode,1,IA1) =                 PEOZ(I_Mode,1,IA1) - Rincr2P 
!     End if 
!     PEON(I_Mode,2,IA2) =              PEON(I_Mode,2,IA2) + Rincr2N 
!     IF I_N_CN Mod 1 = 0 Then 
!     PEON(I_Mode,1,IA1) =                 PEON(I_Mode,1,IA1) + Rincr2N 
!     Else 
!     PEON(I_Mode,1,IA1) =                 PEON(I_Mode,1,IA1) - Rincr2N 
!     End if 
!     End If  '/ 
!     ' 
      PEOZ(I_Mode,1,IA1) = PEOZ(I_Mode,1,IA1) * EOscale 
      PEOZ(I_Mode,2,IA2) = PEOZ(I_Mode,2,IA2) * EOscale 
      PEON(I_Mode,1,IA1) = PEON(I_Mode,1,IA1) * EOscale 
      PEON(I_Mode,2,IA2) = PEON(I_Mode,2,IA2) * EOscale 
!     ' 
!     /' Energy sorting '/ 
!     /' E1 = Etot * get_Gaussintegral(T2-T1,0.03); '/ 
      IF (  Abs(T1-T2) .LT. 1.E-6  ) THEN 
      E1 = 0.5E0 * Etot 
      Else 
      E1ES = Csort * T1 * T2 / ( Abs(T1 - T2) ) 
      E1ES = Min(E1ES,0.5E0*Etot) 
!     /' Asymptotic value after "complete" energy sorting '/ 
      E1FG = Etot * IA1 / I_A_CN 
!     /' in Fermi-gas regime '/ 
      IF (  Etot .LT. 13  )  E1 = E1ES 
!     ' complete energy sorting 
      IF (  Etot .GE. 13 .AND. Etot .LE. 20  ) THEN 
!     ' transition region 
      E1 = E1ES + (Etot-13)/7*(E1FG-E1ES) 
      End If 
      IF (  Etot .GT. 20  )  E1 = E1FG 
!     ' Fermi-gas regime 
      End If 
      E2 = Etot - E1 
      EPART(I_Mode,1,IA1) = Max(E1,0.0) 
!     /' Mean E* in light fragment '/ 
      EPART(I_Mode,2,IA2) = Max(E2,0.0) 
!     /' Mean E* in heavy fragment '/ 
      END DO 
      END DO 
!     ' 
!     ' 
!     /'** RMS angular momentum of fission fragments **'/ 
!     /' Following Naik et al.,EPJ A 31 (2007) 195 and  '/ 
!     /' S. G. Kadmensky,Phys. At. Nucl. 71 (2008) 1193 '/ 
!     ' 
!     '   Scope 
      Spin_CN = P_J_CN 
      P_I_rms_CN = P_J_CN 
      Spin_pre_fission = SPIN_CN 
!     ' CN ground-state spin 
!     ' 
      DO IZ1 = 10 , I_Z_CN - 10 
      AUCD = Int(REAL(IZ1) * REAL(I_A_CN) / REAL(I_Z_CN)) 
      DO IA1 = Int(AUCD - 15) , Int(AUCD + 15) 
      IN1 = IA1 - IZ1 
      IF (  IA1 - IZ1 .GE. 10  ) THEN 
!     /' Rigid momentum of inertia for spherical nucleus '/ 
      I_rigid_spher = 1.16E0**2 * REAL(IA1)**1.6667E0 / 103.8415E0 
!     /' unit: hbar^2/MeV '/ 
      DO I_Mode = 0 , 6 
!     ' 
!     /' First (normally light) fission fragment: '/ 
!     ' 
      beta1 = Beta(I_Mode,1,IZ1) 
      alph = beta1 / SQRT(4.E0 * pi / 5.E0) 
      I_rigid = I_rigid_spher * (1.E0 + 0.5E0*alph + 9.E0/7.E0*alph**2) 
!     /' From Hasse & Myers,Geometrical Relationships ... '/ 
      E_exc = EPART(I_Mode,1,IA1) 
      IF (  E_exc .LT. 0  )  E_exc = 0 
      T = get_U_Temp(REAL(IZ1),REAL(IA1),E_exc,1,1,Tscale,Econd) 
!     '   T = sqr(T^2 + 0.8^2)       ' For ZPM 
!     '   T = T_orbital 
!     '   T =  sqr(T^2 + T_orbital^2) 
      IF (  T_orbital .GT. 0.1  ) THEN 
      T = T_orbital / tanh(T_orbital/T) 
!     ' T_orbital represents the ZPM 
      End If 
      I_eff = I_rigid * (1.E0 - 0.8E0 * exp(-0.693E0 * E_exc / 5.E0)) 
      J_rms = SQRT(2.E0 * I_eff * T) 
!     ' 
      J_rms = J_rms * Jscaling 
!     ' 
      IF (   MOD(IZ1,2)  .EQ. 1 .OR.  MOD(IN1,2)  .EQ. 1  )  J_rms = J_rms + Spin_odd * (REAL(IA1)/140.0)**0.66667 
!     '                * Max(0,1 - (E_exc-1)/9) /' empirical '/ 
!     /' Additional angular momentum of unpaired proton. '/ 
!     /' See also Tomar et al.,Pramana 68 (2007) 111 '/ 
!     ' 
!     ' Print Z1,I_Mode,beta1,T,E_exc,Spin_CN 
!     ' Print " ",I_rigid_spher,I_rigid,I_eff,J_rms 
!     ' 
      J_rms = SQRT(J_rms**2 + (IA1/I_A_CN * Spin_pre_fission)**2) 
!     ' 
      SpinRMSNZ(I_Mode,1,IA1-IZ1,IZ1) = J_rms 
!     ' 
!     '     Print 
!     '     Print IA1,T,E_exc,I_rigid_spher,I_rigid,I_eff,J_rms 
!     ' 
!     /' Second (normally heavy) fission fragment: '/ 
!     ' 
      beta2 = Beta(I_Mode,2,IZ1) 
      alph = beta2 / SQRT(4.E0 * pi / 5.E0) 
      I_rigid = I_rigid_spher * (1.E0 + 0.5E0*alph + 9.E0/7.E0*alph**2) 
!     /' From Hasse & Myers,Geometrical Relationships ... '/ 
      E_exc = EPART(I_Mode,2,IA1) 
      IF (  E_exc .LT. 0  )  E_exc = 0 
      T = get_U_Temp(REAL(IZ1),REAL(IA1),E_exc,1,1,Tscale,Econd) 
!     '    T = sqr(T^2 + 0.8^2)       ' For ZPM 
!     '    T = T_orbital 
!     '    T =  sqr(T^2 + T_orbital^2) 
      IF (  T_orbital .GT. 0.1  ) THEN 
      T = T_orbital / tanh(T_orbital/T) 
!     ' T_orbital represents the ZPM 
      End If 
      I_eff = I_rigid * (1.E0 - 0.8E0 * exp(-0.693E0 * E_exc / 5.E0)) 
      J_rms = SQRT(2.E0 * I_eff * T) 
!     ' 
      J_rms = J_rms * Jscaling 
!     ' 
      IF (   MOD(IZ1,2)  .EQ. 1 .OR.  MOD(IN1,2)  .EQ. 1  )  J_rms = J_rms + Spin_odd * (REAL(IA1)/140.0)**0.66667 
!     '                 * Max(0,1 - (E_exc-1)/9) /' empirical '/ 
!     /' Additional angular momentum of unpaired proton. '/ 
!     /' See also Tomar et al.,Pramana 68 (2007) 111 '/ 
!     ' 
      J_rms = SQRT(J_rms**2 + (IA1/I_A_CN * Spin_pre_fission)**2) 
!     ' 
      SpinRMSNZ(I_Mode,2,IA1-IZ1,IZ1) = J_rms 
!     ' 
!     '      Print IA1,T,E_exc,I_rigid_spher,I_rigid,I_eff,J_rms 
!     ' 
      END DO 
      ENd If 
      END DO 
      END DO 
!     '   End Scope 
!     ' 
!     ' **************************************************************** 
!     ' *** Filling arrays with results in the folding mode (GEFSUB) *** 
!     ' **************************************************************** 
!     ' 
      DO I = 10 , I_A_CN - P_Z_CN - 10 
      DO J = 10 , P_Z_CN - 10 
      DO K = 0 , 6 
      NZMPRE(K,I,J) = 0.0 
      END DO 
      END DO 
      END DO 
!     ' 
!     ' Mode 0 
      DO I = 20 , I_A_CN - 20 
      Ic = I_A_CN - I 
      R_Help = Yield_Mode_0 * (get_U_Gauss_mod(AC_Mode_0 - REAL(I),SigA_Mode_0) + &
          get_U_Gauss_mod(AC_Mode_0 - REAL(Ic),SigA_Mode_0)) 
!     ' Mass yield 
      IF (  I .LT. Ic  ) THEN 
      Zs = ZShift(0,1,I) 
      Else 
      Zs = -ZShift(0,1,Ic) 
      End If 
      DO J = 10 , P_Z_CN - 10 
      Jc = P_Z_CN - J 
      IF (  I-J .GE. 0 .AND. Ic-Jc .GE. 0 .AND. I-J .LE. 200 .AND. Ic-Jc .LE. 200  ) THEN 
      NZMPRE(0,I-J,J) = R_Help * &
          get_U_Gauss_mod(REAL(P_Z_CN)/REAL(I_A_CN)*REAL(I) + Zs - REAL(J),SigPol_Mode_0) * &
          get_U_Even_Odd(J,PEOZ(0,1,I)) * get_U_Even_Odd(I-J,PEON(0,1,I)) 
      End If 
      END DO 
      END DO 
!     ' 
!     ' Mode 1 
      DO I = 20 , I_A_CN - 20 
      Ic = I_A_CN - I 
      R_Help = Yield_Mode_1 * (get_U_Gauss_mod(AC_Mode_1 - REAL(I),SigA_Mode_1) + &
          get_U_Gauss_mod(AC_Mode_1 - REAL(Ic),SigA_Mode_1)) 
!     ' Mass yield 
      IF (  I .LT. Ic  ) THEN 
      Zs = ZShift(1,1,I) 
      Else 
      Zs = -ZShift(1,1,Ic) 
      End If 
      DO J = 10 , P_Z_CN - 10 
      Jc = P_Z_CN - J 
      IF (  I-J .GE. 0 .AND. Ic-Jc .GE. 0 .AND. I-J .LE. 200 .AND. Ic-Jc .LE. 200  ) THEN 
      NZMPRE(1,I-J,J) = R_Help * get_U_Gauss_mod(REAL(P_Z_CN)/REAL(I_A_CN)*REAL(I) + &
          Zs - REAL(J),SigPol_Mode_1)* get_U_Even_Odd(J,PEOZ(1,1,I)) * &
          get_U_Even_Odd(I-J,PEON(1,1,I)) 
      End If 
      END DO 
      END DO 
!     ' 
!     ' Mode 2 
      DO I = 20 , I_A_CN - 20 
      Ic = I_A_CN - I 
      R_Help = Yield_Mode_2 * &
          (get_U_Box2(AC_Mode_2-REAL(I),SQRT(2.0)*S2leftmod*SigA_Mode_2,SQRT(2.0)*SigA_Mode_2,P_A_Width_S2) + &
          get_U_Box2(AC_Mode_2 - REAL(Ic),SQRT(2.0)* &
          S2leftmod*SigA_Mode_2,SQRT(2.0)*SigA_Mode_2,P_A_Width_S2)) 
      IF (  I .LT. Ic  ) THEN 
      Zs = ZShift(2,1,I) 
      Else 
      Zs = -ZShift(2,1,Ic) 
      End If 
      DO J = 10 , P_Z_CN - 10 
      Jc = P_Z_CN - J 
      IF (  I-J .GE. 0 .AND. Ic-Jc .GE. 0 .AND. I-J .LE. 200 .AND. Ic-Jc .LE. 200  ) THEN 
      R_Cut1 = R_Help 
      R_Cut2 = R_Help 
      IF (  J .GT. Jc  ) THEN 
      R_Cut1 = R_Help * get_Gaussintegral(REAL(J)-ZTRUNC50,FTRUNC50*SigZ_Mode_2) 
      Else 
      R_Cut2 = R_Help * get_Gaussintegral(REAL(J)-ZTRUNC50,FTRUNC50*SigZ_Mode_2) 
      End If 
      NZMPRE(2,I-J,J) = R_Help * &
          get_U_Gauss_mod(REAL(P_Z_CN)/REAL(I_A_CN)*REAL(I)+Zs-REAL(J),SigPol_Mode_2) * &
          get_U_Even_Odd(J,PEOZ(2,1,I)) * get_U_Even_Odd(I-J,PEON(2,1,I)) 
      End If 
      END DO 
      END DO 
!     ' 
!     ' Mode 3 
      DO I = 20 , I_A_CN - 20 
      Ic = I_A_CN - I 
      R_Help = Yield_Mode_3 * (get_U_Gauss_mod(AC_Mode_3 - REAL(I),SigA_Mode_3) + &
          get_U_Gauss_mod(AC_Mode_3 - REAL(Ic),SigA_Mode_3)) 
!     ' Mass yield 
      IF (  I .LT. Ic  ) THEN 
      Zs = ZShift(3,1,I) 
      Else 
      Zs = -ZShift(3,1,Ic) 
      End If 
      DO J = 10 , P_Z_CN - 10 
      Jc = P_Z_CN - J 
      IF (  I-J .GE. 0 .AND. Ic-Jc .GE. 0 .AND. I-J .LE. 200 .AND. Ic-Jc .LE. 200  ) THEN 
      NZMPRE(3,I-J,J) = R_Help * &
          get_U_Gauss_mod(REAL(P_Z_CN)/REAL(I_A_CN)*REAL(I)+Zs-REAL(J),SigPol_Mode_3) * &
          get_U_Even_Odd(J,PEOZ(3,1,I)) * get_U_Even_Odd(I-J,PEON(3,1,I)) 
      End If 
      END DO 
      END DO 
!     ' 
!     ' Mode 4 
      DO I = 20 , I_A_CN - 20 
      Ic = I_A_CN - I 
      R_Help = Yield_Mode_4 * (get_U_Gauss_mod(AC_Mode_4 - REAL(I),SigA_Mode_4) + &
          get_U_Gauss_mod(AC_Mode_4 - REAL(Ic),SigA_Mode_4)) 
!     ' Mass yield 
      IF (  I .LT. Ic  ) THEN 
      Zs = ZShift(3,1,I) 
      Else 
      Zs = -ZShift(3,1,Ic) 
      End If 
      DO J = 10 , P_Z_CN - 10 
      Jc = P_Z_CN - J 
      IF (  I-J .GE. 0 .AND. Ic-Jc .GE. 0 .AND. I-J .LE. 200 .AND. Ic-Jc .LE. 200  ) THEN 
      NZMPRE(4,I-J,J) = R_Help * &
          get_U_Gauss_mod(REAL(P_Z_CN)/REAL(I_A_CN)*REAL(I)+Zs-REAL(J),SigPol_Mode_4) * &
          get_U_Even_Odd(J,PEOZ(4,1,I)) * get_U_Even_Odd(I-J,PEON(4,1,I)) 
      End If 
      END DO 
      END DO 
!     ' 
!     ' Mode 11 
      DO I = 20 , I_A_CN - 20 
      Ic = I_A_CN - I 
      R_Help = Yield_Mode_11 * (get_U_Gauss_mod(AC_Mode_0 - REAL(I),SigA_Mode_11) + &
          get_U_Gauss_mod(AC_Mode_0 - REAL(Ic),SigA_Mode_11)) 
!     ' Mass yield 
      DO J = 10 , P_Z_CN - 10 
      Jc = P_Z_CN - J 
      IF (  I-J .GE. 0 .AND. Ic-Jc .GE. 0 .AND. I-J .LE. 200 .AND. Ic-Jc .LE. 200  ) THEN 
      NZMPRE(5,I-J,J) = R_Help * &
          get_U_Gauss_mod(REAL(P_Z_CN)/REAL(I_A_CN)*REAL(I)-REAL(J),SigPol_Mode_0) * &
          get_U_Even_Odd(J,PEOZ(5,1,I)) * get_U_Even_Odd(I-J,PEON(5,1,I)) 
      End If 
      END DO 
      END DO 
!     ' 
!     ' Mode 22 
      DO I = 20 , I_A_CN - 20 
      Ic = I_A_CN - I 
      R_Help = Yield_Mode_22 * (get_U_Gauss_mod(AC_Mode_0 - REAL(I),SigA_Mode_22) + &
          get_U_Gauss_mod(AC_Mode_0 - REAL(Ic),SigA_Mode_22)) 
!     ' Mass yield 
      DO J = 10 , P_Z_CN - 10 
      Jc = P_Z_CN - J 
      IF (  I-J .GE. 0 .AND. Ic-Jc .GE. 0 .AND. I-J .LE. 200 .AND. Ic-Jc .LE. 200  ) THEN 
      NZMPRE(6,I-J,J) = R_Help * &
          get_U_Gauss_mod(REAL(P_Z_CN)/REAL(I_A_CN)*REAL(I)-REAL(J),SigPol_Mode_0) * &
          get_U_Even_Odd(J,PEOZ(6,1,I)) * get_U_Even_Odd(I-J,PEON(6,1,I)) 
      End If 
      END DO 
      END DO 
!     ' 
!     ' 
!     ' Normalization 
      R_Sum = 0 
      DO I = 10 , (I_A_CN - P_Z_CN) - 10 
      DO J = 10 , P_Z_CN - 10 
      NZPRE(I,J) = 0 
      DO K = 0 , 6 
      IF (  NZMPRE(K,I,J) .GT. 0  ) THEN 
      R_Sum = R_Sum + NZMPRE(K,I,J) 
      NZPRE(I,J) = NZPRE(I,J) + NZMPRE(K,I,J) 
!     ' sum of all modes 
      End If 
      END DO 
      END DO 
      END DO 
!     ' Print R_Sum 
      DO I = 10 , (I_A_CN - P_Z_CN) - 10 
      DO J = 10 , P_Z_CN - 10 
      NZPRE(I,J) = NZPRE(I,J) / R_Sum 
      DO K = 0 , 6 
      NZMPRE(K,I,J) = NZMPRE(K,I,J) / R_Sum 
      END DO 
      END DO 
      END DO 
!     ' 
!     ' Calculate and store distributions of fragment excitation energy and spin 
!     ' 
      N_cases = 0 
      DO N_index = 10 , (I_A_CN - P_Z_CN) - 10 
!     ' Neutron number 
      DO Z_index = 10 , P_Z_CN - 10 
!     ' Atomic number 
      DO M_index = 0 , 6 
      IF (  NZMPRE(M_index,N_index,Z_index) .GT. Ymin  ) THEN 
      N_cases = N_cases + 1 
      IF (  N_cases .EQ. Ubound(NZMkey,1)  ) THEN 
!     '           Print "Upper bound of NZkey reached" 
!     '           Print "Result will be incomplete" 
      End If 
      NZMkey(N_cases,1) = M_index 
!     ' Fission mode 
      NZMkey(N_cases,2) = N_index 
!     ' Neutron number of fragment 
      NZMkey(N_cases,3) = Z_index 
!     ' Atomic number of fragment 
      End If 
      END DO 
      END DO 
      END DO 
!     ' Print "N_cases  ",N_cases 
!      WRITE (*,*) "N_cases ",N_cases 
!     ' 
      DO K = 1 , N_cases 
      M_index = NZMkey(K,1) 
!     ' fission mode 
      N_index = NZMkey(K,2) 
!     ' neutron number 
      Z_index = NZMkey(K,3) 
!     ' atomic number 
      A_index = N_index + Z_index 
!     ' 
!     ' Yield 
      Ytab(K) = NZMpre(M_index,N_index,Z_index) 
!     ' 
!     ' Angular momentum: 
      DO I = 1 , 100 
      IF (  M_index .LE. 4  ) THEN 
      IF (  Z_index .LT. 0.5 * P_Z_CN  ) THEN 
      Jtab(K,I) = get_U_LinGauss(REAL(I),SpinRMSNZ(M_index,1,N_index,Z_index)/SQRT(2.0))
      Else 
      Jtab(K,I) = get_U_LinGauss(REAL(I),SpinRMSNZ(M_index,2,N_index,Z_index)/SQRT(2.0))
      End If 
      End If 
      IF (  M_index .EQ. 5  ) THEN 
      Jtab(K,I) = get_U_LinGauss(REAL(I),SpinRMSNZ(1,2,N_index,Z_index)/SQRT(2.0)) 
      End If 
      IF (  M_index .EQ. 6  ) THEN 
      Jtab(K,I) = get_U_LinGauss(REAL(I),SpinRMSNZ(2,2,N_index,Z_index)/SQRT(2.0)) 
      End If 
      END DO 
!     ' 
!     ' Normalize numerically (due to non-continuous values) 
!     '   Scope 
      Rint = 0 
      DO I = 1 , 100 
      Rint = Rint + Jtab(K,I) 
      END DO 
      IF (  Rint .GT. 0  ) THEN 
      DO I = 1 , 100 
      Jtab(K,I) = Jtab(K,I) / Rint 
      END DO 
      End If 
!     '   End Scope 
!     ' 
!     ' 
!     ' Excitation energy: 
!     ' 1. Deformation energy at scission 
      IF (  M_index .EQ. 0  ) THEN 
      IF (  Z_index .LT. 0.5 * P_Z_CN  ) THEN 
      Eexc_mean = Edefo(M_index,1,Z_index) 
      Eexc_sigma = ( &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,1,Z_index) + SIGDEFO_0) - &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,1,Z_index) )) 
      Else 
      Eexc_mean = Edefo(M_index,2,Z_index) 
      Eexc_sigma = ( &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,2,Z_index) + SIGDEFO_0) - &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,2,Z_index) )) 
      End If 
      End If 
      IF (  M_index .GT. 0 .AND. M_index .LE. 4  ) THEN 
      IF (  Z_index .LT. 0.5 * P_Z_CN  ) THEN 
      Eexc_mean = Edefo(M_index,1,Z_index) 
      RS = SIGDEFO/SQRT(R_Att_Sad(M_index)) 
      Eexc_sigma = ( &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,1,Z_index) + RS) - &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,1,Z_index) )) 
      Else 
      Eexc_mean = Edefo(M_index,2,Z_index) 
      RS = SIGDEFO/SQRT(R_Att_Sad(M_index)) 
      Eexc_sigma = ( &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,2,Z_index) + RS) - &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(M_index,2,Z_index) )) 
      End If 
      End If 
      IF (  M_index .EQ. 5  ) THEN 
      Eexc_mean = Edefo(1,2,Z_index) 
      RS = SIGDEFO/SQRT(R_Att_Sad(M_index)) 
      Eexc_sigma = ( &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(1,2,Z_index) + RS) - &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(1,2,Z_index) )) 
      End If 
      IF (  M_index .EQ. 6  ) THEN 
      Eexc_mean = Edefo(2,2,Z_index) 
      RS = SIGDEFO/SQRT(R_Att_Sad(M_index)) 
      Eexc_sigma = ( &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(2,2,Z_index) + RS) - &
          get_LyMass(REAL(Z_index),REAL(A_index),beta(2,2,Z_index) )) 
      End If 
      Eexc_mean = Max(Eexc_mean,0.0) 
!     ' 
!     ' 2. Intrinsic excitation energy at scission 
      IF (  Z_index .LT. 0.5 * REAL(P_Z_CN)  ) THEN 
      Eexc_intr = EPART(M_index,1,A_index) 
      Else 
      Eexc_intr = EPART(M_index,2,A_index) 
      End If 
      IF (  M_index .EQ. 0  ) THEN 
!     ' add shell and pairing of final fragment 
      Eexc_intr = Eexc_intr - get_AME2012(Z_index,A_index) + &
          get_LDMass(REAL(Z_index),REAL(A_index),0.) - 2.0 * 12.0 / SQRT(REAL(A_index)) 
!     ' general shift 
      End If 
      Eexc_intr = Max(Eexc_intr,0.0) 
      Eexc_mean = Eexc_mean + Eexc_intr 
      Eexc_sigma = SQRT(Eexc_sigma**2 + (EexcSIGrel * Eexc_intr)**2) 
!     ' 
!     ' 3. Pairing staggering 
      Eexc_mean = Eexc_mean - get_LyPair(Z_index,A_index) 
!     ' 
!     ' 4. Collective energy 
      Eexc_coll = 0.5 * ECOLLFRAC * &
       (get_De_Saddle_Scission(REAL(P_Z_CN)**2/REAL(I_A_CN)**0.33333E0,ESHIFTSASCI_coll)-&
        E_tunn) 
      Eexc_coll = Max(Eexc_coll,0.0) 
      Eexc_sigma = SQRT(Eexc_sigma**2 + 0.5*(EexcSIGrel*Eexc_coll)**2) 
      Eexc_mean = Eexc_mean + Eexc_coll + 0.5 * E_coll_saddle(M_index) 
!     ' 
!     ' 5. Total excitation energy distribution of fragments (all contributions summed up) 
      DO I = 0 , 1000 
!     ' 100 keV bins up to 100 MeV 
      Etab(K,I) = exp(-(0.1*REAL(I)-Eexc_mean)**2/(2.0 * Eexc_sigma)) 
      END DO 
!     ' 
!     ' Normalize excitation-energy distribution 
!     '   Scope 
      RintE = 0 
      DO I = 0 , 1000 
      RintE = RintE + Etab(K,I) 
      END DO 
      IF (  RintE .GT. 0  ) THEN 
      DO I = 0 , 1000 
      Etab(K,I) = Etab(K,I) / RintE 
      END DO 
      End If 
!     '   End Scope 
!     ' 
      END DO 

end


!*****************************************
!*            special_fission            *
!*****************************************
! Taken from special_fission.f95
subroutine special_fission()

    use iso
    use RESULTS
    use path_module
    implicit none
    integer isomer, get_isotope, n_fission_counter,fission_index, li
    integer istatus
    double precision neutron_separation_energy
    character*3 chz, chn 
    character*256 fission_file
    double precision yield, dytab(10000)
    INTEGER*4 z_t, n_t, a_t
    REAL*4    level_e, level_spin 
    character*256 the_path

    the_path = trim(files_path) // trim('decay_data/')

    n_fission_counter = 0
    fission_index     = 0
    do isomer = 1, number_isotopes           ! loop over all species
      if (z(isomer) .eq. 92 .and. n(isomer).eq.146) then            ! only 238U for now
        z_t        = z(isomer)  ! compound nucleus
        n_t        = n(isomer)  ! compound nucleus
        a_t        = n_t + z_t  ! compound nucleus
        level_e    = 0.0              ! GS
        level_spin = spin(isomer)           ! 
        n_fission_counter = n_fission_counter + 1
        write(chz,'(I3.3)') z(isomer)
        write(chn,'(I3.3)') n(isomer)
        fission_file = the_path // chz // '_' // chn // '_spontaneous_no_scission.dat'
!
!        open(unit=13,file=fission_file,status='old', iostat=istatus)
!        if (istatus .ne. 0) then                      ! need to create the fission file
          print*,'s-Fission of: ',a_t,element_names(z(isomer)), '  ', level_e,level_spin
          print*,'Creating: ',trim(fission_file)
          call GEFSUB(z_t,a_t,level_e,level_spin)
          dytab = Ytab
          dytab = 2./sum(dytab(1:N_cases)) * dytab   ! normalize to 2.
          s_fission_vector(fission_index,:) = 0.
          call special_scission(fission_index,N_cases,NZMkey,dytab,Etab)
          print*,'Neutron yield',s_fission_vector(fission_index,index_n) &
                , sum(s_fission_vector(fission_index,:))
          open(unit=13,file=fission_file,status='unknown')
          do li=1,number_isotopes
            if (s_fission_vector(fission_index,li) .gt. 0.0) then
              write(13,'(2I6,ES20.10)') z(li), n(li), s_fission_vector(fission_index,li)
            end if  
          end do
!        end if
        close(13)
      end if
    end do
    print*,'Number of isotopes with neutron-induced fission: ', n_fission_counter

end


!*****************************************
!*            Special Scission           *
!*****************************************
! Taken from special_fission.f95
subroutine special_scission(fission_index,dytab)

    use iso
    use RESULTS
    use path_module
    implicit none
    integer isomer,M_index, N_index, Z_index, A_index,k,i,j,get_isotope,li,ch_sn, n_ch_n
    integer fission_index, istatus, ch_counter, li_ne, neutron_counter
    double precision neutron_yield, e_ave, e_sum, de_ave2,neutron_separation_energy
    double precision sn, e_excitation(2,1000), e_exc, en(1000),dytab(10000)
    double precision yield_local(300)
    REAL*4    level_e, level_spin
    INTEGER*4 iso_n, iso_z, iso_a
    character*256 the_path

    the_path = trim(files_path) // trim('decay_data/debug/u238sf_A_yield_GLE_no-scission.csv')

    yield_local = 0.
    DO K = 1 , N_cases    
      e_excitation(1,:) = Etab(k,:)
      e_excitation(1,:) = e_excitation(1,:) / sum(e_excitation(1,:)) * dytab(K)
      e_excitation(2,:) = 0.
      M_index = NZMkey(K,1) 
      N_index = NZMkey(K,2) 
      Z_index = NZMkey(K,3) 
      A_index = N_index + Z_index 
      yield_local(A_index) = yield_local(A_index) + dytab(K)
      sn    = neutron_separation_energy(Z_index,N_index)
      li    = get_isotope(Z_index,N_index,0.0d0)
      li_ne = get_isotope(Z_index,N_index-1,0.0d0)
      if (li*li_ne.eq.0) then
        neutron_counter = 0
        do
          neutron_counter = neutron_counter + 1
          li    = get_isotope(Z_index,N_index-neutron_counter,0.0d0)
          li_ne = get_isotope(Z_index,N_index-neutron_counter-1,0.0d0)
          sn = sn + neutron_separation_energy(Z_index,N_index-neutron_counter)
          s_fission_vector(fission_index,index_n) = &
              s_fission_vector(fission_index,index_n) + dytab(K)
          if (li*li_ne.ne.0) then
            exit
          end if
        end do   
      end if  
      sn = 1.d12
      ch_sn = int(sn*10. + 0.5)
      e_ave = 0.
      do i=1,1000
        if (e_excitation(1,i) .gt. 0.) then
          e_exc = 0.1d0 * (i-0.5)
          e_ave = e_ave + e_excitation(1,i) * e_exc
          if (e_exc.le.sn) then  !  below neutron separation threshold
            neutron_yield = 0.
          else  ! neutron emission
            neutron_yield  = sum(mb_distribution(ch_sn:i)) / sum(mb_distribution(1:i)) 
            n_ch_n         = i-ch_sn+1
            en(1:n_ch_n)   = mb_distribution(ch_sn:i)  / sum(mb_distribution(1:i))
            do j=1,n_ch_n
              e_excitation(2,i-j-ch_sn) = e_excitation(2,i-j-ch_sn) + e_excitation(1,i)*en(j)
            end do
          end if
          s_fission_vector(fission_index,li) = &
              s_fission_vector(fission_index,li) + e_excitation(1,i)*(1-neutron_yield)
          s_fission_vector(fission_index,index_n) = &
              s_fission_vector(fission_index,index_n) + e_excitation(1,i)*neutron_yield
          s_fission_vector(fission_index,li_ne) = &
              s_fission_vector(fission_index,li_ne) + e_excitation(1,i)*neutron_yield
        end if  
      end do
   
!
!  repeat this loop until all energy is lt separation energy
!
      do 
        N_index = N_index - 1
        sn    = neutron_separation_energy(Z_index,N_index)
        ch_sn = int(sn*10. + 0.5)
        if (sum(e_excitation(2,ch_sn:1000)) .eq. 0.) exit
        e_excitation(1,:) = e_excitation(2,:) / sum(e_excitation(1,:))
        e_excitation(2,:) = 0.
        li = li_ne
        if (li.eq.0) then
          print*,'isotope after multi neutron emission not found!!!',Z_index,N_index-1
          stop
        end if
        do i=1,1000
          if (e_excitation(1,i) .gt. 0.) then
            e_exc = 0.1d0 * (i-0.5)
            if (e_exc.gt.sn) then ! neutron emission
              neutron_yield = sum(mb_distribution(ch_sn:i)) / sum(mb_distribution(1:i))
              n_ch_n         = i-ch_sn+1
              en(1:n_ch_n)   = mb_distribution(ch_sn:i)  / sum(mb_distribution(1:i))
              do j=1,n_ch_n
                e_excitation(2,i-j-ch_sn) = e_excitation(2,i-j-ch_sn) + e_excitation(1,i)*en(j)
              end do
              s_fission_vector(fission_index,li) = &
                  s_fission_vector(fission_index,li) - e_excitation(1,i)*neutron_yield
              s_fission_vector(fission_index,index_n) = &
                  s_fission_vector(fission_index,index_n) + e_excitation(1,i)*neutron_yield
              s_fission_vector(fission_index,li_ne) = &
                  s_fission_vector(fission_index,li_ne) + e_excitation(1,i)*neutron_yield
            end if  
          end if  
        end do
      end do  
    END DO
    open(unit=23,file=the_path,status='unknown') 
      do A_index=1,300
        if (yield_local(A_index).gt.0.) write(23,*) A_index,yield_local(A_index)
      end do
    close(23)

end


!*****************************************
!*             Init N Fission            *
!*****************************************
! Taken from n_fission.f95
subroutine init_n_fission()

    use iso
    use RESULTS
    use path_module
    implicit none
    integer isomer, get_isotope, n_fission_counter,fission_index, li
    integer istatus
    double precision neutron_separation_energy
    character*3 chz, chn 
    character*256 fission_file
    double precision yield, dytab(10000)
    INTEGER*4 z_t, n_t, a_t
    REAL*4    level_e, level_spin 
    character*256 the_path

    the_path = trim(files_path) // trim('decay_data/n_fission_Sn_plus_500keV/')

    n_fission_counter = 0
    fission_index     = 0
    do isomer = 1, number_isotopes ! loop over all species
      if (z(isomer) .ge. 84) then ! Z=84 is Po. Below no fission is considered to be of importance
        z_t        = z(isomer)  ! compound nucleus
        n_t        = n(isomer) + 1  ! compound nucleus
        a_t        = n_t + z_t  ! compound nucleus
        level_e    = neutron_separation_energy(z(isomer),n(isomer) + 1) + 0.5 ! neutron separation energy of target comp-nucleus + average neutron energy
        if ( level_e .lt. 0.) then
          print*,'Sn+0.1 MeV .lt. 0.)'
          level_e = 0.
        end if  
        level_spin = spin(isomer) ! in principle, it should be spin(isomer) +- spin_of_neutron (0.5, 1.5, ...)
        n_fission_counter = n_fission_counter + 1
        write(chz,'(I3.3)') z(isomer)
        write(chn,'(I3.3)') n(isomer)
        fission_file = the_path // chz // '_' // chn // '_n_induced.dat'
!
        open(unit=13,file=fission_file,status='old', iostat=istatus)
        if (istatus .ne. 0) then                      ! need to create the fission file
          print*,'n-Fission of: ',a_t,element_names(z(isomer)), '  ', level_e,level_spin
          print*,'Creating: ',trim(fission_file)
          call GEFSUB(z_t,a_t,level_e,level_spin)
          dytab = Ytab
          dytab = 2./sum(dytab(1:N_cases)) * dytab   ! normalize to 2.
          s_fission_vector(fission_index,:) = 0.
          call scission(fission_index,N_cases,NZMkey,dytab,Etab)
          print*,'Neutron yield',s_fission_vector(fission_index,index_n) &
                , sum(s_fission_vector(fission_index,:))
          open(unit=13,file=fission_file,status='unknown')
          do li=1,number_isotopes
            if (s_fission_vector(fission_index,li) .gt. 0.0) then
              write(13,'(2I6,ES20.10)') z(li), n(li), s_fission_vector(fission_index,li)
            end if  
          end do
        end if
        close(13)
      end if
    end do
    print*,'Number of isotopes with neutron-induced fission: ', n_fission_counter

end


