PROGRAM GEMINI
USE gem_kind
USE gemsub
USE functs
USE angle
USE evaporation
USE evap
USE light_stuff
USE tl_stuff
USE yrast_mod
USE mass_stuff
USE stat_mod
IMPLICIT NONE

INTEGER :: ii,time_sec
INTEGER :: l,seed(3),values(8)
REAL (KIND=r4) :: start_time,stop_time,num_residue,weight
TYPE (identity) :: cn !compound nucleus



! now the arrays frag(i)%z,frag(i)%A,frag(i)%ex,frag(i)%J contain the
! parameters of the compound system and other systems formed in the decay
! of the compound system.
! The index i=1 corresponds to the compound system, after the
! first binary division the two fragments are put into  i=2 and
! 3. sequential decays
! are put into greater i's.
! n_frag is the total number of fragments formed at any one time.



CALL CPU_TIME(start_time)

CALL DATE_AND_TIME(DATE=para%day,VALUES=values)


! generate seed for random number generator
! from system clock
time_sec = (values(5)*60+values(6))*60+values(7)
seed = 0
seed(1) = time_sec
seed(1) = IOR(seed(1),1)
!seed(1) = 10091 !for debugging purposes
CALL RANDOM_SEED(size=ii)
WRITE (UNIT=*, FMT=*) "seed has size=",ii
CALL RANDOM_SEED(PUT=seed(1:ii))
num_residue = 0.
mult0 = 0.


CALL input ()


! *****************************************************************
! read in mass
CALL prepare_mass ()

! read list of evaporated particles
IF (para%tl_iwbc) THEN
    CALL prepare_Hauser_Feshbach ()
    CALL read_tl ()
ELSE
    call prepare_evap ()
END IF 

IF (para%I_angle .OR. para%k_sum) CALL factorial ()

!if (para%ratio /= 1.0 .AND. para%ratio > 0. .AND.&
!.NOT. para%k_SUM) CALL prepare_deformed_emission (para%ratio,1.)

if (para%ratio /= 1.0 .AND.&
.NOT. para%k_SUM) CALL prepare_deformed_emission (para%ratio,1.)


IF (para%Z_imf_min > 2) CALL prepare_5li_5he ()

IF (para%imf_option == 2) call read_sad ()

! generate array of angular distributions
IF (para%I_angle) call quantum_ang_generator ()
! prepare arrays for choosing emission angle from deformed
! emitter
CALL alloc_storage

!****************************************************************

!if restarting a run set lmin to maximum l of previous
!calculation.
max_l = MAX(para%lmin,Max_L)

!loop over input channel L wave
loop_ell: DO l = max_l,para%lmax

!  determine the number of simulations to be performed.
   N_evt = NINT((2.*L+1)*para%R_num)

!  if restarting calculation, adjust n_evt for first l wave.
   IF (l == Max_L) N_evt = N_evt - Num_L
   IF (N_evt <= 0) N_evt = 0
   cascades: DO II = 1,N_EVT

       CALL fusion(l,cn)      !       input compound nucleus
                              !       parameters

       CALL event (cn,weight)    !       simulate decay of compound
                              !       nucleus

       CALL dmp(l,weight)     !       write event on event file
                              !       (*.evt)

       if (frag(n_frag)%Z > REAL(para%zcn+4)/2&
        .AND. para%tl_iwbc) THEN
          mult0 = mult0 + mult1*weight
          num_residue = num_residue + weight
       END IF 

    END DO cascades

END DO loop_ell
CALL CPU_TIME (stop_time)
WRITE (UNIT=*, FMT=*) "cpu time=",(stop_time-start_time)/60.," minutes"

CALL write_mult (num_residue)
STOP
END PROGRAM gemini