MODULE stat_mod
USE gem_kind
USE mass_stuff
IMPLICIT NONE
PRIVATE 

PUBLIC :: alloc_storage, get_level_density,collective_enhancement

INTEGER, PARAMETER, PUBLIC :: max_channels = 2000
INTEGER, PARAMETER, PUBLIC :: n_light_ch = 49
INTEGER, PARAMETER, PUBLIC :: n_light_ksum = 250000
INTEGER, PARAMETER, PUBLIC :: n_light_normal = 24000
INTEGER, PARAMETER, PUBLIC :: Z_shell = 2


TYPE, PUBLIC :: dau_info
    INTEGER :: ll,level
    REAL (KIND=r4) :: e_kinetic,l_p_s,excit
    REAL (KIND=r4) :: Ex_residue,J_residue
    REAL (KIND=r4) :: k_residue,s_particle
    REAL (KIND=r4) :: K_particle
END TYPE dau_info

TYPE (dau_info), DIMENSION(n_light_ch), PUBLIC :: dau

INTEGER, DIMENSION(max_channels), PUBLIC :: Za
INTEGER, DIMENSION(max_channels), PUBLIC :: Aa
INTEGER, PUBLIC :: I_channels
INTEGER, DIMENSION(60), PUBLIC :: evap_mode
INTEGER, DIMENSION(50), PUBLIC :: list_HF,list_HF0
INTEGER, PUBLIC :: N_light

REAL (KIND=r4), DIMENSION(max_Channels), PUBLIC :: U_therm,gamma,time
REAL (KIND=r4), PUBLIC :: Mass0,Bs
REAL (KIND=r4), PUBLIC :: pe,ee,m_inertia_0

TYPE, PUBLIC :: evap_store
    INTEGER :: spin,kk2
    REAL (KIND=r4) :: gam,lps,energy
END TYPE evap_store
TYPE (evap_store), DIMENSION(:), ALLOCATABLE, PUBLIC :: storage


CONTAINS
    SUBROUTINE alloc_storage ()
        IF (para%k_sum) THEN
            n_light = n_light_ksum
        ELSE
            n_light = n_light_normal
        END IF      
        ALLOCATE (storage(0:n_light))
    END SUBROUTINE alloc_storage

    SUBROUTINE get_level_density(U,J,aden,M_Inertia,&
                  entropy,level_density)
      
    !+
    ! 
    ! FUNCTIONAL DESCRIPTION: 
    ! 
    ! returns the daughter level density divided by the
    ! parent level density
    ! 
    ! FORMAL PARAMETERS:
    !  
    ! U - thermal excitation energy above yrast line
    ! A - mass number of nucleus
    ! entropy - entropy 
    ! M_inertia - moment of inertia
    ! aden - little "a" level density parameter
    ! J - spin of nucleus
    !  
    !-
    REAL (KIND=r4), INTENT(IN) :: U,J
    REAL (KIND=r4), INTENT(IN) :: aden,M_inertia,entropy
    REAL (KIND=r4), INTENT(OUT) :: level_density

    REAL (KIND=r4) :: sigma_2,temp
    REAL (KIND=r4) :: ee_ek,level_density_spherical

        
        ee_ek = entropy - ee
    !   choice of level density


        temp = (U/aden)**0.5 
        sigma_2 = m_inertia*temp/40.848



    !   normal fermi gas level densities for spherical nucleus
        level_density_spherical = (2.0*J+1.0)/pe*EXP(ee_ek)&
        /(1.0+U**1.25*sigma_2**1.5)/aden**0.25


        level_density = level_density_spherical


    RETURN
    END SUBROUTINE get_level_density


    SUBROUTINE collective_enhancement(IZ,IA,U,aden,inertia,level_density)
    INTEGER, INTENT(IN) :: iz,ia
    REAL (KIND=r4), INTENT(IN) :: inertia,U,aden
    REAL (KIND=r4), INTENT(INOUT) :: level_density
    REAL (KIND=r4) :: beta2,inertia2,f,e_cr,d_cr,temp,sigma2_perp

     RETURN

    CALL getmass(iz,ia,beta_two=beta2)    

!   this follows Hansen and Jenssen. 
    IF (beta2 == 0.0) RETURN
    inertia2 = inertia *(1.0+beta2/3.0)
    E_cr = 120.0*beta2**2*REAL(IA)**0.33333
    d_cr = 1400.*beta2**2/REAL(ia)**0.6666666

    f = (U-e_cr)/d_cr
    IF (f > 65.0) THEN
        f = 0.0
        RETURN
    ELSE
        f = 1.0/(1.0+EXP(f))
    END IF 

    temp = (U/aden)**0.5 
    sigma2_perp = inertia2*temp/40.848
    level_density = ((sigma2_perp-1.0)*f + 1.0)*level_density
    
    END SUBROUTINE collective_enhancement

END MODULE stat_mod