MODULE sortsub
USE gem_kind
IMPLICIT NONE
PRIVATE

PUBLIC :: write_file_info, get_event, read_file_header, open_event_file
PUBLIC :: start_sort, name

INTEGER, PARAMETER, PUBLIC :: max_frag=70
TYPE (identity), DIMENSION(max_frag), PUBLIC :: frag
INTEGER, PUBLIC :: n_frag,n_binary
REAL (KIND=r4), PUBLIC :: weight

CONTAINS
SUBROUTINE write_file_info(unit_out)

![subroutine_header_comments]

INTEGER, INTENT(IN) :: unit_out


WRITE (UNIT=unit_out, FMT="(TR2,A40)") para%reaction

WRITE (UNIT=unit_out, FMT="(TR2,a19,a9)")&
"date of gemini run=",para%day

WRITE (UNIT=unit_out, FMT="(TR2,A28,f10.3,a8)")&
"number of events per l-wave=",para%R_num,"*(2*L+1)"

WRITE (UNIT=unit_out, FMT="(TR2,a22,i3,TR2,a4,i3,/TR2,a19,f7.2,TR1,a3)")&
"Compund nucleus:- Zcn=",para%Zcn, "Acn=",para%Acn,&
"Excitation energy =", para%Excn, "MeV"
IF (para%dexcn > 0) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A31,F7.3)")&
    "spread in excitation energy= +-",para%dexcn
END IF


WRITE (UNIT=unit_out, FMT="(TR2,a47,i3,TR2,i3)")&
"minimum and maximum L waves simulated by Gemini="&
,para%lmin, para%lmax

If(para%N_weight > 0)then
 WRITE (UNIT=unit_out, FMT="(TR2,a24,I3,A)")&
 "Weighting was turned on for ",para%N_weight,"binary decays"
ELSE
 WRITE (UNIT=unit_out, FMT="(TR2,A)")&
 "NO Weighting"
end if


IF (para%I_angle) THEN
 WRITE (UNIT=unit_out, FMT="(TR2,a)")&
 "Angles and Velocities of all fragments calculated"
ELSE
 WRITE (UNIT=unit_out, FMT="(TR2,a)")&
 "No angles and velocities calculated"
END IF


IF (para%IMF_option == 1 .AND. para%kramer /= 0) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A,F8.3)")&
    "Kramers factor=",para%kramer
END IF

IF (para%t_delay > 1000.0 .and. para%ratio /= 1.0&
.and. para%imf_option == 0) THEN
  WRITE (UNIT=unit_out, FMT="(TR2,A38)")&
 "evaporation from deformed system only"
ELSE IF (para%t_delay > 0.0) THEN
  WRITE (UNIT=unit_out, FMT="(TR2,A17,TR1,F9.3,TR1,A24)")&
 "Fission delay of",para%t_delay,"x10-21 seconds included"
  WRITE (UNIT=unit_out, FMT="(TR2,A30,TR1,F9.3,TR1,A24)")&
 "second moment of eta dependence=",para%sig_delay&
 ," included"
  IF (para%sharp_delay) THEN
      WRITE (UNIT=unit_out, FMT="(TR2,A)")&
      "Fission and Imf emission totally suppressed during delay time"
  ELSE
      WRITE (UNIT=unit_out, FMT="(TR2,A)")&
      "Decay width ramps up during delay time"
  END IF
ELSE
  WRITE (UNIT=unit_out, FMT="(TR2,A16)") "No fission delay"
END IF

IF (para%imf_option == 0) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "No IMF emission, only evaporation"
ELSE IF (para%imf_option == 1) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "evaporation and symmetric fission only"
ELSE IF (para%IMF_option == 2) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "All asymmetric divisions possible"
END IF

IF (para%imf_option >= 1) THEN
    if(para%b_scale <= 0.0) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Sierks fission barriers used"
    ELSE
    WRITE (UNIT=unit_out, FMT="(TR2,A,TR2,F5.2)")&
    "Sierks fission barrier scaled by ",para%b_scale
    END IF 
END IF 


SELECT CASE (para%aden_type)
CASE (0)
    WRITE (UNIT=unit_out, FMT="(TR2,A23,F5.2)")&
     "level density constant=",para%aden_0
CASE (-1)
    WRITE (UNIT=unit_out, FMT="(TR2,A40)")&
    "Toke + Swiatecki level density parameter"
CASE (-2)
    WRITE (UNIT=unit_out, FMT="(TR2,A32)")&
    "Ignatyuk level density parameter"
CASE (-3)
    WRITE (UNIT=unit_out, FMT="(TR2,A50)")&
    "Gottschalk and Ledergerber level density parameter"
CASE (-4)
    WRITE (UNIT=unit_out, FMT="(TR2,A45)")&
    "Temperature dependent level density parameter"
CASE (-6)
    WRITE (UNIT=unit_out, FMT="(TR2,A45)")&
    "Lestones Temperature dependent level density parameter"
CASE (-9)
    WRITE (UNIT=unit_out, FMT="(TR2,A,F5.2)")&
    "Fineman's Temp dependent level-density parameter, kappa=", &
    para%aden_0
CASE DEFAULT
    WRITE (UNIT=unit_out, FMT="(TR2,A,I5)")&
    "aden_type=", para%aden_type
END SELECT

IF (para%imf_option == 1)&
WRITE (UNIT=unit_out, FMT="(TR2,A,F6.3)") "af/an=",para%a_scale


IF (para%tl_IWBC) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A29)")&
    "Transmission coeff. from IWBC"
ELSE
    WRITE (UNIT=unit_out, FMT="(TR2,A35)")&
    "Sharp cut-off Transmission coeff."
END IF


IF (para%ratio == 1.0 .OR. para%ratio == 0.0) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A,F5.2)")&
    "spherical transmission coef ,ratio=",para%ratio
ELSE
    WRITE (UNIT=unit_out, FMT="(TR2,A,F5.2)")&
    "deformed transmission coef ,ratio=",para%ratio
END IF

!IF (para%tl_iwbc) THEN
!    WRITE (UNIT=unit_out, FMT="(TR2,A,E8.2)")&
!    "threshold probabilty for evaporation=",para%threshold
!END IF 


!IF (para%collective_enhancement) THEN
!    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
!    "rotationally-symmetric-deformed level densities"
!ELSE
!    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
!    "spherical level densities"
!END IF

IF (para%lestone) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Lestone modification to the transition state formulism"
ELSE
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Standard transitional state formulism"
END IF


IF (para%quantum) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Quantum treatment of spin projection"
ELSE
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Semi-classical treatment of spin projection and emission angles"
END IF

IF (para%z_imf_min == 3) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)") "evaporation Z=0,1,2"
ELSE IF (para%Z_imf_min == 4) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)") "evaporation Z=0,1,2,3"
ELSE IF (para%Z_imf_min == 5) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)") "evaporation Z=0,1,2,3,4"
END IF


IF (para%mass_option == 0) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "liquid drop masses used"
ELSE IF (para%Mass_Option == 1) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "experimental masses used"
ELSE IF (para%Mass_Option == 2) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "experimental masses with shell fade-out used"
END IF

IF (para%polarization) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Polarization of protons considered"
END IF


IF (para%exotic_index == 0) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Only boring fragments evaporated"
ELSE IF (para%exotic_index == 1) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Some mildly exciting fragments evaporated"
ELSE IF (para%exotic_index == 2) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Some real exotic fragments evaporated"
ELSE 
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "Some real exotic fragments evaporated"
END IF

if (para%IMF_option == 2) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A,TR2,F6.3)")&
    "af/an=",para%a_scale
    WRITE (UNIT=unit_out, FMT=*)&
    "barrier scaled by",para%b_scale
END IF 

IF (para%time_flag) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "emission times are written in the event file"
END IF

IF (para%time_flag) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "fragment spins written in the event file"
END IF


IF (para%k_sum) THEN
    WRITE (UNIT=unit_out, FMT="(TR2,A)")&
    "K and shape distributions considered"

END IF 

WRITE (UNIT=unit_out, FMT="(TR2,A,F7.2,A)")"E2 strength=",&
para%E2_strength, "W.u."

WRITE (UNIT=unit_out, FMT="(TR2,A,F7.2,A)")"E1 strength=",&
para%E1_strength, "W.u."


RETURN 
END SUBROUTINE write_file_info



SUBROUTINE get_event (l,I_read_error)

![subroutine_header_comments]
INTEGER, INTENT(OUT) :: l
INTEGER, INTENT(OUT) :: I_read_error
INTEGER :: i
INTEGER (KIND=2), DIMENSION(600) :: array
INTEGER (KIND=2) :: n_record
REAL (KIND=r4) :: partr4
INTEGER :: ipt

READ (UNIT=2, IOSTAT=I_read_error) n_record,(array(i),i=1,n_record-1)
IF (I_read_error /= 0) RETURN
ipt = 1
n_frag = array(ipt)
ipt = ipt + 1
n_binary = array(ipt)
ipt = ipt + 1
l = array(ipt)
weight = TRANSFER(array(ipt+1:ipt+2),partr4)

 ipt = ipt + 2

DO i=1,n_frag
    frag(i)%A = array(i+ipt)
    frag(i)%Z = IBITS(frag(i)%A,8,8)
    IF (frag(i)%z == 255) frag(i)%z = -1
    frag(i)%A = IBITS(frag(i)%A,0,8)
END DO
ipt = ipt + n_frag

FORALL (i=1:n_frag) frag(i)%id = array(ipt+i)
ipt = ipt + n_frag


if (para%i_angle) THEN
    FORALL (i=1:n_frag)
       frag(i)%vel%v = TRANSFER(array(ipt+2*i-1:ipt+2*i),partr4)
    END FORALL
    ipt = ipt + 2*n_frag

    FORALL (i=1:n_frag)
       frag(i)%vel%theta = TRANSFER(array(ipt+2*i-1:ipt+2*i),partr4)
    END FORALL
    ipt = ipt + 2*n_frag

    FORALL (i=1:n_frag)
       frag(i)%vel%phi = TRANSFER(array(ipt+2*i-1:ipt+2*i),partr4)
    END FORALL
    ipt = ipt + 2*n_frag
END IF

if (para%polarization) THEN
   forall (i=1:n_frag)  frag(i)%M = REAL(array(i+ipt) )/2.0
   ipt = ipt + n_frag
END IF

if (para%time_flag) THEN
    FORALL (i=1:n_frag)
       frag(i)%time = TRANSFER(array(ipt+2*i-1:ipt+2*i),partr4)
    END FORALL
    ipt = ipt + 2*n_frag
END IF 

if (para%j_flag .OR. para%polarization) THEN
   forall (i=1:n_frag)  frag(i)%J = REAL(array(i+ipt) )/2.0
   ipt = ipt + n_frag
END IF


if (para%ex_flag) THEN
   forall (i=1:n_frag)  frag(i)%exx = REAL(array(i+ipt) )/10.0
   ipt = ipt + n_frag
END IF

RETURN  
END SUBROUTINE get_event


SUBROUTINE read_file_header(unit_in,I_read_error)

![subroutine_header_comments]
INTEGER, INTENT(IN) :: unit_in
INTEGER, INTENT(OUT) :: i_read_error
REAL (KIND=r4), save :: r_num=0.0

READ (UNIT=unit_in, IOSTAT=I_read_error) para
para%r_num = para%r_num + r_num
r_num = para%r_num

RETURN 
END SUBROUTINE read_file_header


SUBROUTINE open_event_file(unit_in,I_read_error,unit_out)

![subroutine_header_comments]
INTEGER, INTENT(IN) :: unit_in,unit_out
INTEGER, INTENT(OUT) :: I_read_error
INTEGER :: i
CHARACTER(LEN=4), DIMENSION(3), PARAMETER :: EXT=(/".inp",".evt",".out"/)
CHARACTER(LEN=30) :: FILE_name
CHARACTER (LEN=34) :: evtfile,infile,outfile
LOGICAL, SAVE :: be=.FALSE.

WRITE (UNIT=*, FMT=*)"input file"

READ (UNIT=5, FMT="(a)", IOSTAT=I_read_error) file_name
IF (I_read_error /= 0) THEN
   WRITE (UNIT=*, FMT=*)"WHAT!!!"
   RETURN
END IF

file_name = ADJUSTL(file_name) !remove leading blanks
outfile = TRIM(file_name)//ext(3)
evtfile = TRIM(file_name)//ext(2)

WRITE (UNIT=*, FMT=*)"file=",evtfile

OPEN (                   &
UNIT=unit_in,            &
 ACCESS = "SEQUENTIAL", &
 FORM = "UNFORMATTED",  &
 STATUS = "OLD",        &
 position = "REWIND",   &
 ACTION = "READ",       &
 SHARED,                &
 FILE = evtfile,        &
 IOSTAT=I_read_error)   

IF (I_read_error /= 0) THEN
  WRITE (UNIT=*, FMT=*)"Could not open event file",evtfile
  RETURN
END IF


IF (unit_out >= 1) THEN 
   IF (.NOT. be)THEN !if an outfile is already open, doen"t open a new one
     OPEN (             &
     UNIT=unit_out,     &
     STATUS = "NEW",    &
     ACTION = "WRITE",  &
     FILE = outfile)
     be = .TRUE.
   END IF
   WRITE (UNIT=unit_out,FMT=*)"evtfile=",evtfile
END IF


RETURN 
END SUBROUTINE open_event_file


SUBROUTINE start_sort(unit_in,unit_out,I_read_error)

![subroutine_header_comments]
INTEGER, INTENT(IN) :: unit_in,unit_out
INTEGER, INTENT(OUT) :: I_read_error

CALL open_event_file(unit_in,I_read_error,unit_out)
if (i_read_error /= 0)RETURN
CALL read_file_header(unit_in,I_read_error)
IF (i_read_error == 0) THEN
    call write_file_info(6)
    IF (unit_out >= 1) THEN
        call write_file_info(unit_out)
    END IF
END IF
RETURN
END SUBROUTINE start_sort

FUNCTION name(Z) RESULT(funct_out)
INTEGER, INTENT(IN) :: z

CHARACTER(LEN=2) :: funct_out
CHARACTER(LEN=2), DIMENSION(103), PARAMETER :: element = (/&
" H","He","Li","Be"," B"," C"," N"," O"," F","Ne",  &
"Na","Mg","Al","Si"," P"," S","Cl","Ar"," K","Ca",  &
"Sc","Ti"," V","Cr","Mn","Fe","Co","Ni","Cu","Zn",  &
"Ga","Ge","As","Se","Br","Kr","Rb","Sr"," Y","Zr",  &
"Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn",  &
"Sb","Te"," I","Xe","Cs","Ba","La","Ce","Pr","Nd",  &
"Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb",  &
"Lu","Hf","Ta"," W","Re","Os","Ir","Pt","Au","Hg",  &
"Tl","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th",  &
"Pa"," U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm",  &
"Md","No","Lw"/)
IF (z > 0 .AND. Z < 104) THEN
    funct_out = element(Z)
ELSE IF (z == 0) THEN
    funct_out = " n"
ELSE
    funct_out = "  "
  END IF
RETURN
END FUNCTION name
   
END MODULE sortsub