! FDMNES subroutines

! Post treatment of the spectra calculated by FDMNES
! 1) Sommation with energy shift of varius spectra
! 2) Convolution by a lorentzienne for XANES
!              L(x) = (1/(pi*b)) * 1 / ( 1 + ( (x-a)/b )**2 )
! Integration over energy for Dafs

subroutine Convolution(Analyzer,bav_open,Bormann,Circular,Conv_done,Convolution_out, &
        Delta_edge,E_cut_imp,E_cut_man,Ecent,Elarg,Epsii_ref,Epsii_ref_man,Estart,Fit_cal,Gamma_hole,Gamma_hole_man, &
        Gamma_max,Gamma_max_man,Header,ical,icheck,indice_par,iscratchconv,itape1,kw_conv,Length_line,n_col_max, &
        n_shift_core,ngamh,ngroup_par,nkw_conv,nomfich,nomfichbav,npar,nparm,param,Scan_a,Shift_core,Typepar,ncal)

  use declarations
  implicit none

  integer:: Dafs_exp_type, eof, i, i_conv, i_surface, i_bulk, i_z, ical, icheck, ie, ie1, ie2, ifich, igr, &
    ii, init, ip, ipar, ipl, ipr, ipr1, ipr2, is, iscr, iscratchconv, istop, istat, itape1, j, &
    jfich, jpl, js, jseuil, k, Length, Length_line, mfich, n, n_bir, n_col, n_col_max, n_couple_weight_co, &
    n_energ_tr, n_mat_pol, n_selec_core, n_shift_core, n_Stokes, n_bulk_site, n_z_SRXRD, n_z_SRXRD_abs, &
    n_z_SRXRD_abs_bulk, n_z_SRXRD_bulk, ncal, ne2, nef, nelor, nen2, nenerg, nenerge, nes, nes_in, nfich, &
    ngamh, ngroup_par, ninit0, ninit1, ninitm, nj, nkw_conv, nnombre, np_stokes, nparm, nphim, npldafs, npldafs_b, &
    npldafs_NR, npldafs_t, npldafs_th, NR, nseuil, nxan, nw

! njp : Points beyond the energy range to make the border effect less strong
  integer, parameter:: njp = 500

  integer, dimension(3):: hkl_S
  integer, dimension(10):: num_core
  integer, dimension(ngroup_par):: npar
  integer, dimension(ngroup_par,nparm):: indice_par
  integer, dimension(:), allocatable:: i_done, indf, ne, ninit, nphi, Numat
  integer, dimension(:,:), allocatable:: Couple_co, hkl_dafs, nsup, ne_initl
  integer, dimension(:), allocatable:: n_div_fpp, nL_header

  character(len=1):: Diese
  character(len=9):: keyword
  character(len=132):: identmot, mot, mots
  character(len=Length_word):: nomab   
  character(len=Length_name):: chemin, Convolution_out, Convolution_out_L, File_scan_out, nomfich, nomfichbav, File_name, Word
  character(len=9), dimension(nkw_conv):: kw_conv
  character(len=9), dimension(ngroup_par,nparm) :: typepar
  character(len=Length_word), dimension(:), allocatable:: nom_col
  character(len=Length_name), dimension(:), allocatable:: Convolution_out_all, File_in, File_name_temp, File_scan_in, &
                                                          File_scan_out_all
  character(len=Length_word), dimension(:), allocatable:: Stokes_name
  
  complex(kind=db):: cf
  complex(kind=db), dimension(:), allocatable:: Bragg_abs_T, Bragg_abs_TT, dampl, dph_t, f0_th, Sum_Bragg_NR_Bulk, Sum_Bragg_NR_T
  complex(kind=db), dimension(:,:), allocatable:: Sum_Bragg_abs_f0, Sum_Bragg_abs_f0_bulk, Bragg_abs_T_scan, Sum_Bragg_NR_T_scan, &
                                                  Sum_Bragg_nonabs_f, Sum_Bragg_nonabs_f_bulk
  complex(kind=db), dimension(:,:,:), allocatable:: Ad, Adafs, As, Mu_m, Mu_mat_comp, Mu_t
  complex(kind=db), dimension(:,:,:,:), allocatable:: As_SRXRD, As_SRXRD_bulk, mu, mus

  logical:: Abs_before, Ampl_phase, Analyzer, Another_one, Arc, bav_open, Bormann, Check_birefringence, Dafs, &
    Dafs_bio, Check_conv, chem, Circular, Conv_done, Cor_abs, decferm, Deuxieme, Double_cor, &
    E_cut_man, E_cut_param, Energphot, Epsii_ref_man, Extrap, Fermip, Fit_cal, &
    Forbidden, fprim, fprime_atom, Full_self_abs, Gamma, Gamma_hole_man, Gamma_hole_param, Gamma_max_man, Gamma_var, &
    Gaussian_default, Green_int, Header, Just_total, Magn, Multi_file, no_extrap, Non_resonant_bulk, nxan_lib, Scan_a, &
    Scan_true, Seah, Self_abs, Signal_Sph, SRXRD, Step_conv, Stokes, Stokes_Dafs, Stokes_xan, Sup_sufix, Tenseur, Tenseur_car, &
    Thomson, Transpose_file, U_iso_man, XES
    
  logical, dimension(:), allocatable:: Bulk_site, run_done, Skip_run

  real(kind=db):: a, a1, a2, a3, a4, Abs_cap, Abs_U_iso_inp, alambda, Asea, b, b1, b2, b3, b4, bba, bbb, c, c_micro, &
    conv_mbarn_nelec, ct_epsilon, ct_nelec, d, d_dead, de_obj, de1, de2, Delta_edge, Delta_gauss, &
    E, E_cut_imp, E_obj, E_step_conv, e1m, Ecent, E_cut, E_cut_orig, Eintmax, Elarg, Eph, Epsii_ref, Esmin, &
    Estart, f0_forward, fac, fpp0, Gamm, Gamma_h, Gamma_max, &
    natomsym, p1, p2, Low_Energy_Step, Pdt, Pdt_bulk, S0_2, Sample_thickness, Shift_U_iso,  &
    Surface_ref, Tab_width, Volume_maille, Volume_maille_bulk, Weight_co

  real(kind=db), dimension(3):: angxyz, axyz
  real(kind=db), dimension(10):: Gamma_hole
  real(kind=db), dimension(ngroup_par,nparm) :: param
  real(kind=db), dimension(n_shift_core):: Shift_core
  real(kind=db), dimension(:), allocatable:: Abs_NR_SRXRD, Abs_NR_SRXRD_abs, Abs_NR_SRXRD_bulk, Abs_NR_SRXRD_abs_bulk, Abs_U_iso, &
    Angle, bb, betalor, e1, e2, Efermip, Elor, En_fermi, Energ, &
    Energe, Energ_tr, Eph1, Ephm, Ephoton, Es, Es_temp, Eseuil, fpp_avantseuil, fi, fr, L_dafs, Length_abs, lori, & 
    lorix, lorr, lorrx, natomsym_f, Pds, p1f, p2f, Pos_z_SRXRD, Pos_z_SRXRD_abs, Pos_z_SRXRD_bulk, Pos_z_SRXRD_abs_bulk, &
    Shift_site, V0muf, Yr, Yi
  real(kind=db), dimension(:,:), allocatable:: Epsii, Mu_tt, mua_r, mua_i, Shift_init, Stokes_param, Xa, Xanes, Xs, Xs_SRXRD, &
                                               Xs_SRXRD_bulk
  real(kind=db), dimension(:,:,:), allocatable:: Icirc, Icirccor, Icor, Icircdcor, Idcor, Mu_mat
  real(kind=db), dimension(:,:,:,:), allocatable:: Mus_mat

! Put to .true. when TDDFT
  Another_one = .false.

  do  ! end of loop at the end of the routine

    call Initialization_conv(Abs_before,Abs_U_iso_inp,Ampl_phase,Arc,Asea,bav_open,Check_conv,Chem,Convolution_out, &
             d_dead,Dafs,Dafs_exp_type,decferm,Delta_gauss,Deuxieme,Double_cor,E_cut,E_cut_imp,E_cut_param,Eintmax,f0_forward, &
             Fermip,File_scan_out,Forbidden,fprim,fprime_atom,Full_self_abs,Gamma_hole_param,Gamma_var,Gaussian_default,Green_int, &
             hkl_S,icheck,jseuil,Just_total,Magn,Multi_file,n_couple_weight_co,n_Stokes,nelor,nfich,npldafs,npldafs_t,npldafs_th, &
             nseuil,num_core,nxan_lib,S0_2,Sample_thickness,scan_true,seah,self_abs,Shift_U_iso,Signal_Sph,Stokes,Stokes_dafs, &
             Stokes_xan,Surface_ref,Tenseur,Tenseur_car,Thomson,Transpose_file,U_iso_man,Volume_maille,Volume_maille_bulk,XES, &
             Weight_co)
  
! -- Lecture --------------------------------------------

    allocate( File_name_temp(10000) ) 
     
    call Dimension_lect_conv(Circular,File_name_temp,itape1,kw_conv,Length_line,Multi_file,n_couple_weight_co,n_energ_tr, &
                   n_Stokes,nelor,nfich,nkw_conv,npldafs_th)

    allocate( run_done(nfich) )
    allocate( skip_run(nfich) )
    allocate( i_done(nfich) )
    allocate( Stokes_param(5,n_Stokes) )
    allocate( Stokes_name(n_Stokes) )
    allocate( V0muf(nfich) )

    rewind(itape1)

    Stokes_param(1:3,:) = 0._db
! Analyzer: when the 2 angles are zero, no analyzer
! Rotation angle de of the analyzer
    Stokes_param(4,:) = 0._db
! Bragg angle of the analyzer
    Stokes_param(5,:) = 0._db
    Stokes_name(:) = 'no_name'

    run_done(:) = .false.
    Skip_run(:) = .false. 
    do ifich = 1,nfich
      i_done(ifich) = ifich
    end do
    mfich = nfich

    boucle_jj: do ii = 1,1000000
      n = nnombre(itape1,132)
      read(itape1,'(A)',iostat=eof) mot
      if( eof /= 0 ) exit boucle_jj
      keyword = identmot(mot,9)
      k = 0
      if( keyword == 'run_done' ) then
        do i = 1,mfich
          read(itape1,*,iostat=eof) j, is
          if( eof /= 0 ) exit boucle_jj
          if( j == 0 ) then
            run_done(i) = .true.
          else
            k = k + 1
            if( is == 0 ) Skip_run(k) = .true.
          endif
          if( run_done(i) ) then
            nfich = nfich - 1
            do j = i+1,mfich
              i_done(j) = i_done(j) - 1
            end do
          endif
        end do
        exit
      endif
    end do boucle_jj

    if( nfich == 0 ) then
      nfich = 1
      mfich = 1
      allocate( File_in(nfich) )
      Word = ' '
      Word = nomfich
      Length = len_trim(Word)
      Word(Length+1:Length+4) = '.txt'
      File_in(1) = Word
    else
      allocate( File_in(nfich) )
    endif

    allocate( Efermip(nfich) )
    allocate( En_fermi(nfich) )
    allocate( Eseuil(nfich) )
    allocate( File_scan_in(nfich) )
    allocate( fpp_avantseuil(nfich) )
    allocate( n_div_fpp(nfich) )
    allocate( ne(nfich) )
    allocate( ninit(nfich) )
    allocate( nL_header(nfich) )
    allocate( Numat(nfich) )
    allocate( Pds(nfich) )
    allocate( Shift_site(nfich) )

    Pds(1) = 1._db
    Efermip(1) = 0._db
    ninit(:) = 1
    ninitm = 1
    n_selec_core = 0
    Shift_site(:) = 0._db

    if( ncal > 1 .and. ical == 1 ) then
      do jfich = 1,mfich
        iscr = 100 + jfich
        open( iscr, status = 'scratch' )
      end do
    endif

    allocate( Elor(nelor) )
    allocate( betalor(nelor) )
    allocate( f0_th(npldafs_th) )
    allocate( Energ_tr(n_energ_tr) )
    allocate(  Couple_co(n_couple_weight_co,2) )

    call Lecture_conv( Abs_before, Abs_U_iso_inp, Ampl_phase, Another_one, Arc,asea, Betalor, Check_birefringence ,Chem, &
         Convolution_out, Couple_co, d_dead, Dafs_exp_type, Decferm, Delta_gauss, Deuxieme, Double_cor, E_cut, E_cut_man, &
         E_step_conv, Ecent, Efermip, Eintmax, Elarg, Elor, Energ_tr, f0_th, Fermip, File_in, File_name, File_name_temp, &
         File_scan_in, File_scan_out, Forbidden, fprim, Gaussian_default, Gamma_hole, Gamma_hole_man, Gamma_max, Gamma_max_man, &
         Gamma_var, hkl_S, itape1, Just_total, kw_conv, mfich, Multi_file, n_couple_weight_co, n_energ_tr, n_selec_core,&
         n_stokes, nelor, nfich, nkw_conv, no_extrap, npldafs_th, num_core, nxan_lib, Pds, Run_done, S0_2, Sample_thickness, &
         Scan_a, Scan_true, Seah, Shift_site, Shift_U_iso, Step_conv, Stokes, Stokes_name, Stokes_param, Thomson, Transpose_file, &
         U_iso_man, Weight_co, XES )

    if( Fit_cal ) Just_total = .true.
 
    if( nelor == 0 ) deallocate( betalor, Elor )
    deallocate( File_name_temp )

  allocate( File_scan_out_all(nfich) ) 
  allocate( Convolution_out_all(nfich) )
  Convolution_out_all(:) = ' ' 
  File_scan_out_all(:) = ' '
  
  call Conv_out_name(Abs_U_iso_inp,bav_open,check_conv,Chem,Chemin,convolution_out,Convolution_out_all,Dafs_bio,Delta_gauss, &
                        Deuxieme,E_cut,E_cut_man,Epsii_ref,Epsii_ref_man,File_in,File_scan_in,File_scan_out,File_scan_out_all, &
                        Gamma_hole(1),Gamma_hole_man,Gamma_max,Gamma_max_man,Length_line,nfich,nomfichbav,Scan_a,Scan_true, &
                        U_iso_man,XES)

  if( .not. ( Seah .or. Arc ) .and. nelor == 0 ) then
    nelor = 1
    allocate( Elor(nelor) )
    allocate( betalor(nelor) )
    Elor(1) = 0._db
    betalor(1) = 0._db
  endif

! -- Table dimensions -------------------------------------

  ninitm = 1
  nL_header(:) = 0
  do ifich = 1,nfich
    open(2, file = File_in(ifich), status='old', iostat=istat)
    n = nnombre(2,Length_line)
    if( n == 0 ) then
      do
        read(2,'(a1)') Diese
        if( Diese /= '#' ) then
          backspace(2)
          n = nnombre(2,Length_line)
          exit
        endif
        nL_header(ifich) = nL_header(ifich) + 1 
      end do
    endif
    if( n > 8 ) read(2,*) Eseuil(ifich), Numat(ifich), nseuil, jseuil, fpp_avantseuil(ifich), V0muf(ifich), En_fermi(ifich), ninit0
    if( n /= ninit0 + 11 .and. n /= ninit0 + 12 .and. n /= ninit0 + 13 .and. n /= ninit0 + 14 ) then
      call write_error
      do ipr = 6,9,3
        write(ipr,'(///,A/)') ' Old format in the first line of not convoluted file !'
      end do
      stop
    endif
    ninitm = max( ninitm, ninit0 )
    Close(2)
  end do

  if( n_selec_core /= 0 ) then
    do i = 1,n_selec_core
      if( num_core(i) > ninitm ) then
        call write_error
        do ipr = 6,9,3
          write(ipr,130) num_core(1:n_selec_core)
          write(ipr,135) ninitm
        end do
        stop
      endif
    end do
  endif
  
  call n_div_fpp_avanseuil(Eseuil,n_div_fpp,nfich)

  allocate( Abs_U_iso(nfich) )
  allocate( Bulk_site(nfich) )
  allocate( Eph1(nfich) )
  allocate( Ephm(nfich) )
  allocate( Epsii(ninitm,nfich) )
  allocate( natomsym_f(nfich) )
  allocate( ne_initl(ninitm,nfich) )
  allocate( nsup(ninitm,nfich) )
  allocate( Shift_init(ninitm,nfich) )

  Abs_U_iso(:) = 0._db
  Bulk_site(:) = .false.
  Epsii(:,:) = 0._db
  Shift_init(:,:) = 0._db

  call Dimension_file(Abs_U_iso,Bulk_site,Cor_abs,Eintmax,En_Fermi,Eph1,Ephm,Epsii,f0_forward,File_in, &
          File_scan_in,fprim,Full_self_abs,Green_int,jseuil,Length_line,Magn,n_bulk_site,n_col_max,n_mat_pol, &
          n_z_SRXRD,n_z_SRXRD_abs,n_z_SRXRD_abs_bulk,n_z_SRXRD_bulk,natomsym_f,ne,nfich,ninit1,ninit,ninitm,nL_header, &
          Non_resonant_bulk,nphim,npldafs,npldafs_NR,nseuil,nxan,nxan_lib,Sample_thickness,Scan_true,Self_abs,Signal_Sph,SRXRD, &
          Surface_ref,Tenseur,Tenseur_car,V0muf,Volume_maille,Volume_maille_bulk)

  if( .not. Cor_abs ) Double_cor = .false.
  
  if( jseuil >= 4 ) no_extrap = .true.

! Modification when fitting
  if( Fit_cal ) then

    do igr = 2,ngroup_par
      istop = 0
      do ipar = 1,npar(igr)
        if( typepar(igr,ipar) /= 'shift' .and. typepar(igr,ipar) /= 'weight' ) cycle
        if( indice_par(igr,ipar) > nfich ) then
          call write_error
          do ipr = 6,9,3
            write(ipr,145) typepar(igr,ipar), indice_par(igr,ipar), nfich
          end do
          istop = 1
         endif
      end do
    end do
    if( istop == 1 ) stop

    do igr = 2,ngroup_par
      do ipar = 1,npar(igr)
        select case( typepar(igr,ipar) )
          case('abs_u_iso')
            Abs_U_iso_inp = param(igr,ipar)
            U_iso_man = .true.
          case('aseah')
            asea = param(igr,ipar)
          case('ecent')
            Ecent = param(igr,ipar)
            Ecent = Ecent / rydb
          case('gaussian')
            Delta_gauss = param(igr,ipar)
          case('elarg')
            Elarg = param(igr,ipar)
            Elarg = Elarg / rydb
          case('gamma_max')
            Gamma_max = param(igr,ipar)
            Gamma_max = Gamma_max / rydb
          case('gamma_hol')
            Gamma_hole(1) = param(igr,ipar)
            Gamma_hole(1) = Gamma_hole(1) / rydb
            Gamma_hole_param = .true.
          case('e_cut','efermi')
            E_cut = param(igr,ipar)
            E_cut_param = .true.
          case('weight_co')
            Weight_co = param(igr,ipar)
          case('shift')
            if( .not. run_done( indice_par(igr,ipar) ) ) then
              ifich = i_done( indice_par(igr,ipar) )
              Shift_site( ifich ) = param(igr,ipar)
            endif
          case('weight')
            if( .not. run_done( indice_par(igr,ipar) ) ) then
              ifich = i_done( indice_par(igr,ipar) )
              Pds( ifich ) = param(igr,ipar)
            endif
        end select
      end do
    end do
  endif

  deallocate( i_done )

  if( Cor_abs .or. Dafs_bio ) then
    npldafs_b = npldafs
  else
    npldafs_b = 0
  endif
  if( Analyzer ) then
    npldafs_t = npldafs
  else
    npldafs_t = npldafs / 2
  endif

  allocate( hkl_dafs(3,npldafs_b) )

  Dafs = npldafs > 0
   
  allocate( Bragg_abs_T(npldafs) )
  allocate( Bragg_abs_TT(npldafs) )
  Bragg_abs_TT(:) = (0._db,0._db)
  allocate( fi(npldafs) ); allocate( fr(npldafs) )
  allocate( Sum_Bragg_NR_T(npldafs) ) 
  Sum_Bragg_NR_T(:) = (0._db, 0._db )
  allocate( Sum_Bragg_NR_Bulk(npldafs_NR) ) 
  Sum_Bragg_NR_Bulk(:) = (0._db, 0._db )
  allocate( nphi(npldafs) )
  if( SRXRD ) then
    allocate( Abs_NR_SRXRD(n_z_SRXRD) )
    allocate( Abs_NR_SRXRD_abs(n_z_SRXRD_abs) )
    allocate( Abs_NR_SRXRD_bulk(n_z_SRXRD_bulk) )
    allocate( Abs_NR_SRXRD_abs_bulk(n_z_SRXRD_abs_bulk) )
    allocate( Pos_z_SRXRD(n_z_SRXRD) )
    allocate( Pos_z_SRXRD_abs(n_z_SRXRD_abs) )
    allocate( Pos_z_SRXRD_bulk(n_z_SRXRD_bulk) )
    allocate( Pos_z_SRXRD_abs_bulk(n_z_SRXRD_abs_bulk) )
    allocate( L_dafs(npldafs) )
    allocate( Sum_Bragg_abs_f0(npldafs,n_z_SRXRD_abs) )
    allocate( Sum_Bragg_abs_f0_bulk(npldafs,n_z_SRXRD_abs_bulk) )
    allocate( Sum_Bragg_nonabs_f(npldafs,n_z_SRXRD) )
    allocate( Sum_Bragg_nonabs_f_bulk(npldafs,n_z_SRXRD_bulk) )
    allocate( Length_abs(npldafs) )
  endif

! Reading of Thomson part
  if( npldafs > 0 ) then

    if( Thomson ) then
      n = min( npldafs, npldafs_th )
      Sum_Bragg_NR_T(1:n) = f0_th(1:n)
    else
      Pdt = 0._db; Pdt_bulk = 0._db
      i_bulk = 0
      do ifich = 1,nfich
        if( Bulk_site(ifich) ) i_bulk = i_bulk + 1
        if( Skip_run(ifich) ) cycle
        open(2, file = File_in(ifich), status='old', iostat=istat)
        if( istat /= 0) call write_open_error(File_in(ifich),istat,1)
        do i = 1,nL_header(ifich)+1
          read(2,*)
        end do
        n = nnombre(2,Length_line)
        if( n > 0 ) then
          read(2,*) ( fr(ipl), fi(ipl), ipl = 1,npldafs )
          Sum_Bragg_NR_T(:) = Sum_Bragg_NR_T(:) + Pds(ifich) * cmplx( fr(:), fi(:), db )
          Pdt = Pdt + Pds(ifich)
        endif
        Close(2)
      end do
      if( abs(Pdt) > 1e-10_db ) Sum_Bragg_NR_T(:) = Sum_Bragg_NR_T(:) / Pdt
    endif

  endif

  deallocate( f0_th )

  if( E_cut_man .or. E_cut_param ) En_Fermi(:) = E_cut

  if( Gaussian_default ) Delta_gauss = max( Eseuil(1) / 10000, 0.0001_db )

  Delta_gauss = Delta_gauss / sqrt( 8 * log(2._db) )

  Shift_site(:) = Shift_site(:) / rydb
  E_step_conv = E_step_conv / Rydb
  En_fermi(:) = En_fermi(:) / rydb
  Eph1(:) = Eph1(:) / rydb
  Ephm(:) = Ephm(:) / rydb
  Epsii(:,:) = Epsii(:,:) / rydb
  Eseuil(:) = Eseuil(:) / rydb
  V0muf(:) = V0muf(:) / rydb
  Delta_gauss = Delta_gauss / rydb
  Shift_U_iso = Shift_U_iso / rydb
  if( Fermip ) Efermip(1:nfich) = Efermip(1:nfich) / rydb
  if( .not. ( Seah .or. Arc ) ) Elor(:) = Elor(:) - En_fermi(1)

! Elaboration of energy grid

  call Shift_eval(delta_edge,Energphot,Eph1,Ephm,Epsii,Epsii_ref,Epsii_ref_man,Eseuil,Esmin, &
                        n_selec_core,n_shift_core,nfich,ninit1,ninit,ninitm,num_core,Shift_core,Shift_init,Shift_site)
 
  Low_Energy_Step = 0.5_db / rydb
  nes_in = 100000
  allocate( Es_temp(nes_in) )
   
  call Output_Energy_Grid( E_step_conv, Energphot, Eph1, Es_temp, Eseuil, Esmin, Estart, File_in, Length_line, &
                       Low_Energy_Step, n_selec_core, ne, ne_initl, nes, nes_in, nfich, ninit, ninitm, nL_header, nsup, &
                       num_core, Shift_init, Step_conv )  
  allocate( Es(nes) )
  Es(1:nes) = Es_temp(1:nes)
  deallocate( Es_temp )

  np_stokes = n_stokes * npldafs / 4
  allocate( angle(nphim) )
  allocate( As(nes,nphim,npldafs) )
  allocate( Icor(nes,nphim,npldafs) )
  allocate( Icirccor(nes,nphim,np_stokes) )
  allocate( Icircdcor(nes,nphim,np_stokes) )
  allocate( Icirc(nes,nphim,np_stokes) )
  allocate( Idcor(nes,nphim,npldafs) )
  allocate( indf(nes) )
  allocate( Sum_Bragg_NR_T_scan(nphim,npldafs) )
  allocate( p1f(nes) )
  allocate( p2f(nes) )
  allocate( Bragg_abs_T_scan(nphim,npldafs) )
  allocate( Xs(nes,nxan) )
  if( SRXRD ) then
    allocate( Xs_SRXRD(nes,n_z_SRXRD_abs) )
    allocate( Xs_SRXRD_bulk(nes,n_z_SRXRD_abs_bulk) )
    allocate( As_SRXRD(nes,nphim,npldafs,n_z_SRXRD_abs) )
    allocate( As_SRXRD_bulk(nes,nphim,npldafs,n_z_SRXRD_abs_bulk) )
  endif
  if( Check_birefringence ) then
    n_bir = 4
  else
    n_bir = 1
  endif
  allocate( Mus_mat(nes,n_stokes,n_bir,n_mat_pol) )
  allocate( Mu_mat_comp(nes,4,n_mat_pol) )
  allocate( mus(nes,nphim,npldafs_b,2) )

  Stokes = n_stokes > 0
  Stokes_xan = Stokes .and. n_mat_pol > 0
  Stokes_Dafs = Stokes .and. .not. Stokes_xan

  if( Circular ) then
! The analysor is supposed perfect
    i = n_Stokes - 3
! Sigma - Sigma
    Stokes_param(3,i) = 1._db; Stokes_param(4,i) = 0._db
    Stokes_param(5,i) = pi / 4
    i = i + 1
! Sigma - Pi
    Stokes_param(3,i) = 1._db; Stokes_param(4,i) = pi / 2
    Stokes_param(5,i) = pi / 4
    i = i + 1
    Stokes_param(5,i) = pi / 4
    Stokes_param(3,i) = -1._db; Stokes_param(4,i) = 0._db
    i = i + 1
    Stokes_param(3,i) = -1._db; Stokes_param(4,i) = pi / 2
    Stokes_param(5,i) = pi / 4
  endif
  
  if( Tenseur ) then
    n_col = 2*npldafs
  elseif( Bormann ) then
    n_col = nxan + n_mat_pol * n_stokes + 2*npldafs
  else
    if( .not. Analyzer ) then
      n_col = nxan + n_mat_pol * ( n_stokes * n_bir + 2 ) + npldafs / 2
    else
      n_col = nxan + n_mat_pol * ( n_stokes * n_bir + 2 ) + npldafs
    endif
    if( Dafs_bio ) then
      if( icheck > 1 ) then
         n_col = n_col + 4*npldafs
         nw = 5
      else
         nw = 1
      endif
    else
      if( fprim ) n_col = n_col + 2*npldafs
      if( Ampl_phase ) n_col = n_col + 2*npldafs
      if( .not. Analyzer ) then
        if( Self_abs ) n_col = n_col + 2*npldafs + npldafs / 2
        if( Full_self_abs ) n_col = n_col + 2*npldafs + npldafs / 2    ! 4 mu + Ic
        if( Double_cor ) n_col = n_col + npldafs / 2
        if( Stokes_Dafs ) n_col = n_col + np_stokes / 2
        if( Stokes .and. Full_self_abs ) n_col = n_col + np_stokes / 2
        if( Stokes .and. Double_cor ) n_col = n_col + np_stokes / 2
      else
        if( Self_abs ) n_col = n_col + 2*npldafs + npldafs
        if( Full_self_abs ) n_col = n_col + 4*npldafs + npldafs    ! 4 mu + Ic
        if( Double_cor ) n_col = n_col + npldafs
        if( Stokes_Dafs ) n_col = n_col + np_stokes
        if( Stokes .and. Full_self_abs ) n_col = n_col + np_stokes
        if( Stokes .and. Double_cor ) n_col = n_col + np_stokes
      endif
    endif
  endif

  allocate( nom_col(n_col) )

! Preparation is done --------------------------------------

  Sup_sufix = ninit(1) > 1
   
  call Col_name(Ampl_phase,Analyzer,Bormann,Cor_abs,Dafs_bio,Double_cor,File_in(nfich),File_scan_in(nfich),fprim,Full_self_abs, &
      hkl_dafs,Length_line,n_bir,n_col,n_mat_pol,n_stokes,nL_header(nfich),nom_col,npldafs,npldafs_b,nxan,Self_abs, &
      Signal_sph,Stokes,Stokes_name,Stokes_param,Sup_sufix,Tenseur,XES)

! This loop is for the case of output for all indata files  
  boucle_conv_file: do i_conv = 0,nfich
  
    Convolution_out_L = Convolution_out
    
    if( i_conv > 0 ) then
      if ( Just_total .or. nfich == 1 ) exit
      Convolution_out_L = Convolution_out_all(i_conv)
      File_scan_out = File_scan_out_all(i_conv)
    endif

  if( SRXRD ) then
    Xs_SRXRD(:,:) = 0._db
    Xs_SRXRD_bulk(:,:) = 0._db
    As_SRXRD(:,:,:,:) = (0._db,0._db)
    As_SRXRD_bulk(:,:,:,:) = (0._db,0._db)
  endif

  Xs(:,:) = 0._db
  Mus_mat(:,:,:,:) = 0._db
  Mu_mat_comp(:,:,:) = (0._db,0._db)
  if( Dafs ) As(:,:,:) = (0._db,0._db)
  if( Cor_abs ) mus(:,:,:,:) = (0._db,0._db)
  Sum_Bragg_NR_T_scan(:,:) = ( 0._db, 0._db )
  angle(:) = 0._db
  Bragg_abs_T_scan(:,:) = ( 0._db, 0._db ) 
 
  ifich = 0
  i_bulk = 0
  i_surface = 0
  do jfich = 1,mfich
    if( Bulk_site(jfich) ) then
      i_bulk = i_bulk + 1
    else
      i_surface = i_surface + 1
    endif
    if( run_done(jfich) ) cycle
    ifich = ifich + 1

    if( .not. Just_total .and. i_conv > 0 .and. i_conv /= ifich ) cycle
     
    if( Fermip ) then
      E_cut_orig = Efermip(ifich)
    else
      E_cut_orig = En_fermi(ifich)
    endif

    do init = 1,ninit(ifich)

      if( n_selec_core /= 0 ) then
        do i = 1,n_selec_core
          if( init == num_core(i) ) exit
        end do
        if( i > n_selec_core ) cycle
      endif

      nenerg = ne_initl(init,ifich)

      allocate( Adafs(nenerg,nphim,npldafs) )
      Adafs(:,:,:) = (0._db, 0._db)
      if( Cor_abs ) then
        allocate( mu(nenerg,nphim,npldafs,2) )
        mu(:,:,:,:) = (0._db, 0._db)
      endif
      allocate( Xanes(nenerg,nxan) )
      allocate( Mu_mat(nenerg,6,n_mat_pol) )
      Xanes(:,:) = 0._db
      Mu_mat(:,:,:) = 0._db

! -- Lecture -----------------------------------------------------------

      open(2, file = File_in(ifich), status='old', iostat=istat)

      do i = 1,nL_header(ifich) + 1
        read(2,*)
      end do

      natomsym = natomsym_f(ifich)
      
      n = nnombre(2,Length_line)
      if( n > 0 ) then
        read(2,*)      ! This line has already been red
        read(2,*) ( fr(ipl), fi(ipl), ipl = 1,npldafs )
        Bragg_abs_T(:) = cmplx( fr(:), fi(:),db )
        Bragg_abs_TT(:) = Bragg_abs_TT(:) + Pds(ifich) * Bragg_abs_T(:)
        if( Cor_abs ) then
          read(2,*) axyz(:), angxyz(:), ( hkl_dafs(:,ipl), ipl = 1,npldafs)
          angxyz(:) = angxyz(:) * radian
          axyz(:) = axyz(:) / bohr
        elseif( SRXRD ) then
          if( Bulk_site(ifich) ) then
            read(2,*) Pos_z_SRXRD_abs_bulk(i_bulk), Abs_NR_SRXRD_abs_bulk(i_bulk), Abs_cap, n
            read(2,*) ( fr(ipl), fi(ipl), ipl = 1,npldafs )
            Sum_Bragg_abs_f0_bulk(1:npldafs,i_bulk) = cmplx( fr(1:npldafs), fi(1:npldafs), db )
          else
            read(2,*) Pos_z_SRXRD_abs(i_surface), Abs_NR_SRXRD_abs(i_surface), Abs_cap, n, NR
            read(2,*) ( fr(ipl), fi(ipl), ipl = 1,npldafs )
            Sum_Bragg_abs_f0(1:npldafs,i_surface) = cmplx( fr(1:npldafs), fi(1:npldafs), db )
          endif 
          do i_z = 1,n
            if( Bulk_site(ifich) ) then
              read(2,*) Pos_z_SRXRD_bulk(i_z), Abs_NR_SRXRD_bulk(i_z)
              read(2,*) ( fr(ipl), fi(ipl), ipl = 1,npldafs )
              Sum_Bragg_nonabs_f_bulk(1:npldafs,i_z) = cmplx( fr(1:npldafs), fi(1:npldafs), db )
            else
              read(2,*) Pos_z_SRXRD(i_z), Abs_NR_SRXRD(i_z)
              read(2,*) ( fr(ipl), fi(ipl), ipl = 1,npldafs )
              Sum_Bragg_nonabs_f(1:npldafs,i_z) = cmplx( fr(1:npldafs), fi(1:npldafs), db )
            endif 
          end do
          if( NR == -1 ) then
            read(2,*) ( fr(ipl), fi(ipl), ipl = 1,npldafs )
            Sum_Bragg_NR_bulk(1:npldafs) = cmplx( fr(1:npldafs), fi(1:npldafs), db )
          endif
          read(2,*) Length_abs(:)
          if( Bulk_site(ifich) ) read(2,*) L_dafs(:)
          n = nnombre(2,Length_line)
        endif
      elseif( Tenseur ) then
        Sum_Bragg_NR_T(:) = ( 0._db, 0._db )
        Bragg_abs_T(:) = ( 0._db, 0._db )
        Bragg_abs_T(1) = ( 1._db, 0._db )
        if( Tenseur_car ) then
          Bragg_abs_T(4) = ( 1._db, 0._db )
          Bragg_abs_T(6) = ( 1._db, 0._db )
        endif
        Bragg_abs_TT(:) = Bragg_abs_TT(:) + Pds(ifich) * Bragg_abs_T(:)
      endif

      read(2,*)

      allocate( Energ(nenerg) )

      if( Dafs ) then

        nphi(:) = 1

        if( Cor_abs ) allocate( mua_r(npldafs,2) )
        if( Full_self_abs ) allocate( mua_i(npldafs,2) )

        do ie = nsup(init,ifich)+1,nenerg
          if( Tenseur .and. .not. Magn ) then
            Read(2,*) Energ(ie), ((fr(ipl), ipl = 1,npldafs), j = 1,init)
            if( .not. scan_true ) Adafs(ie,1,:) = cmplx(fr(:), 0._db,db)
          elseif( nxan > 0 .or. n_mat_pol > 0 ) then
            if( Full_self_abs ) then
              Read(2,*) Energ(ie),( ( Xanes(ie,ipl), ipl = 1,nxan ), ( Mu_mat(ie,:,ipl), ipl = 1,n_mat_pol ), &
                                   ( fr(ipl), fi(ipl), (mua_r(ipl,i), mua_i(ipl,i), i = 1,2), ipl = 1,npldafs ), j = 1,init)
            elseif( self_abs ) then
              Read(2,*) Energ(ie),(( Xanes(ie,ipl), ipl = 1,nxan ), ( Mu_mat(ie,:,ipl), ipl = 1,n_mat_pol ), &
                                  ( fr(ipl), fi(ipl), (mua_r(ipl,i), i = 1,2), ipl = 1,npldafs ), j = 1,init)
            else
              Read(2,*) Energ(ie), (( Xanes(ie,ipl), ipl = 1,nxan ), ( Mu_mat(ie,:,ipl), ipl = 1,n_mat_pol ), &
                                  ( fr(ipl), fi(ipl), ipl = 1,npldafs ), j = 1,init )
            endif
          else
            if( Full_self_abs ) then
              Read(2,*) Energ(ie), ( ( fr(ipl), fi(ipl), ( mua_r(ipl,i), mua_i(ipl,i), i = 1,2), ipl = 1,npldafs ), j = 1,init)
            elseif( self_abs ) then
              Read(2,*) Energ(ie), ( ( fr(ipl), fi(ipl), ( mua_r(ipl,i), i = 1,2 ), ipl = 1,npldafs ), j = 1,init)
            else
              Read(2,*) Energ(ie), ( ( fr(ipl), fi(ipl), ipl = 1,npldafs), j = 1,init )
            endif
          endif
          if( .not. scan_true ) then
            Adafs(ie,1,:) = cmplx( fr(:), fi(:), db)
            if( Full_self_abs ) mu(ie,1,:,:) = cmplx( mua_r(:,:), mua_i(:,:), db)
            if( Self_abs ) mu(ie,1,:,:) = cmplx( mua_r(:,:), 0._db, db)
          endif
        end do

        if( Cor_abs ) deallocate( mua_r)
        if( Full_self_abs ) deallocate( mua_i )

      else

        do ie = nsup(init,ifich)+1,nenerg
          Read(2,*) Energ(ie), (( Xanes(ie,ipl), ipl = 1,nxan ), ( Mu_mat(ie,:,ipl), ipl = 1,n_mat_pol ), j = 1,init)
        end do

      endif

      close(2)

      if( Scan_true ) then
        open(2, file = File_scan_in(ifich), status='old',iostat=istat)

        do ipl = 1,npldafs
          read(2,*) nphi(ipl)
        end do
        
        do ie = nsup(init,ifich)+1,nenerg
          read(2,*) mots
          do ipl = 1,npldafs
            if( Dafs_bio ) then
              read(2,*) hkl_dafs(:,ipl)
            else
              read(2,*) mots
            endif
            do i = 1,nphi(ipl)
              if( Self_abs ) then
                read(2,*) angle(i), ( a, b, a3, a4, j = 1,init ), ( c, d, c, d, j = init+1,ninit(ifich) ), &
                          a1, b1, a2, b2
                mu(ie,i,ipl,1) = cmplx( a3, 0._db, db )
                mu(ie,i,ipl,2) = cmplx( a4, 0._db, db )
              elseif( Full_self_abs ) then
                read(2,*) angle(i), ( a, b, a3, b3, a4, b4, j = 1,init ), ( c, d, c, d, c, d, j = init+1,ninit(ifich) ), &
                          a1, b1, a2, b2
                mu(ie,i,ipl,1) = cmplx( a3, b3, db )
                mu(ie,i,ipl,2) = cmplx( a4, b4, db )
              else
                if( Dafs_bio ) then
                  read(2,*) ( a, b, j = 1,init ), ( c, d, j = init+1,ninit(ifich) ), a1, b1, a2, b2
                else
                  read(2,*) angle(i), ( a, b, j = 1,init ), ( c, d, j = init+1,ninit(ifich) ), a1, b1, a2, b2
                endif
              endif
              Adafs(ie,i,ipl) = cmplx( a, b, db)
              Sum_Bragg_NR_T_scan(i,ipl) = cmplx( a1, b1, db)
              Bragg_abs_T_scan(i,ipl) = cmplx( a2, b2, db)
            end do
          end do
        end do
        close(2)

      endif

      Energ(:) = Energ(:) / rydb

      if( abs( S0_2 - 1._db ) > eps10 ) then
        Xanes(:,:) = S0_2 * Xanes(:,:)
        Mu_mat(:,:,:) = S0_2 * Mu_mat(:,:,:)
        if( Dafs ) then
          Adafs(:,:,:) = S0_2 * Adafs(:,:,:)
          if( Cor_abs ) mua_r(:,:) = S0_2 * mua_r(:,:)
          if( Full_self_abs ) mua_i(:,:) = S0_2 * mua_i(:,:)
        endif
      endif

      do ie = nsup(init,ifich),1,-1
        Energ(ie) = Energ(ie+1) - Low_Energy_Step
      end do

      if( Decferm ) then
      ! Fermi does not follow the shift of the energy
        if( Energphot ) then
          E_cut = E_cut_orig - Shift_init(init,ifich) + Eseuil(ifich) 
        else
          E_cut = E_cut_orig - Shift_init(init,ifich) + Eseuil(ifich) - Esmin 
        endif
      elseif( Energphot ) then
        E_cut = E_cut_orig + Eseuil(ifich) 
      else
        E_cut = E_cut_orig 
      endif

      allocate( Mu_m(nenerg,4,n_mat_pol) )

      if( .not. ( Conv_done .or. Skip_run(ifich) ) ) then

        nenerge = nenerg + njp

        if( Cor_abs ) then
          allocate( Ad(nenerg,nphim,3*npldafs) )
        else
          allocate( Ad(nenerg,nphim,npldafs) )
        endif
        allocate( Xa(nenerge,nxan) )

        if( Seah .or. Arc ) then
          nelor = nenerg
          allocate( Elor(nelor) )
          allocate( betalor(nelor) )
          Elor(:) = Energ(:)
        endif

        if( XES .or. Green_int ) then
          betalor(:) = 0._db
        elseif( seah ) then
          call seahdench(asea,E_cut,Gamma_max,nelor,Elor,betalor)
        elseif( Arc ) then
          call gammarc(Ecent,Elarg,Gamma_max,E_cut,nelor,Elor,betalor)
        endif

        if( .not. ( Gamma_hole_man .or. Gamma_hole_param ) ) then
          if( ninit(ifich) == 1 .or. ninit1 == ninit(ifich) .or. init <= ninit1 ) then
            js = jseuil
          else
            js = jseuil + 1
          endif
          Gamma_h = Tab_width(Eseuil(ifich),js,nseuil,Numat(ifich))
        elseif( ngamh == 1 .or. ninit(ifich) == 1 ) then
          Gamma_h = Gamma_hole(1)
        elseif( ngamh == ninit(ifich) ) then
          Gamma_h = Gamma_hole(init)
        elseif( init <= ninit1 ) then
          Gamma_h = Gamma_hole(1)
        else
          Gamma_h = Gamma_hole(2)
        endif
        betalor(:) = betalor(:) + Gamma_h

        if( i_conv == 0 ) then
          if( Check_conv ) then
            ipr1 = 3
          else
            ipr1 = 6
          endif
          ipr2 = 6
          do ipr = ipr1,ipr2,3
            if( ifich == 1 .and. ( init == 1 .or. init == num_core(1) ) )  then
              if( seah ) then
                Write(ipr,'(/A)') ' Seah-Dench model'
                Write(ipr,150) Gamma_max*rydb, asea
              elseif( Arc ) then
                Write(ipr,'(/A)') ' Arctangent model'
                Write(ipr,160) Gamma_max*rydb, Ecent*rydb, Elarg*rydb
              endif
            endif
            if( nfich == 1 .and. ninit(ifich) == 1 ) then
              Write(ipr,170) Gamma_h*rydb, E_cut_orig*rydb, Shift_init(init,ifich)*rydb
            elseif( nfich == 1 ) then
              Write(ipr,172) Gamma_h*rydb, E_cut_orig*rydb, Shift_init(init,ifich)*rydb, init
            elseif( ninit(ifich) == 1 ) then
              Write(ipr,174) Gamma_h*rydb, E_cut_orig*rydb, Shift_init(init,ifich)*rydb, ifich
            else
              Write(ipr,176) Gamma_h*rydb, E_cut_orig*rydb, Shift_init(init,ifich)*rydb, ifich, init
            endif
            if( U_iso_man .and. Abs_U_iso_inp > eps10 ) then
              Write(ipr,'(a11,f10.7)') '    U_iso =', Abs_U_iso_inp 
            elseif( Abs_U_iso(ifich) > eps10 .and. .not. U_iso_man ) then
              Write(ipr,'(a11,f10.7,a22)') '    U_iso =', Abs_U_iso(ifich),'  coming from cif file' 
            endif
          end do
        endif

        if( i_conv == 0 .and. ( ifich == 1 .or. Esmin /= Eseuil(ifich) ) &
                .and. ( init == 1 .or. init == num_core(1) ) .and. Gamma_max > eps10 ) then
          de_obj = ( Elor(nelor) - Elor(1) ) / 15
          E_obj = Elor(1) - 0.00001_db
          do ie = 1,nelor
            Gamm = betalor(ie) - Gamma_h
            if( Gamm > 0._db ) exit
          end do
          n = ie - 1
          do ie = 1,nelor
            if( ie > 1 .and. ie < n ) cycle
            if( Energphot ) then
              E = Elor(ie) - Eseuil(ifich) - V0muf(ifich)
            else
              E = Elor(ie) - V0muf(ifich)
            endif
            Gamm = betalor(ie) - Gamma_h
            if( Gamm > 0._db ) then
              alambda = sqrt( 2 / ( sqrt( E**2 + Gamm**2 ) - E ) )
            else
              alambda = 0._db
            endif
            E = E + V0muf(ifich)
            do ipr = ipr1,ipr2,3
              if( ie == 1 ) then
                if( Energphot ) then
                  Write(ipr,'(A)') '    E - E_edge   Width_(eV) lambda_(A)'
                else
                  Write(ipr,'(A)') '      E_(eV)     Width_(eV) lambda_(A)'
                endif
              endif
              if( ipr == 6 ) then
                if( Elor(ie) < E_obj .and. ie /= nelor ) cycle
                E_obj = E_obj + de_obj
              endif
              if( alambda * bohr < 10000000 ) then
                Write(ipr,180) E * Rydb, betalor(ie) * Rydb, alambda * bohr
              else
                Write(ipr,185) E * Rydb, betalor(ie) * Rydb, alambda * bohr
              endif
            end do
          end do
        endif
  
        allocate( dampl(nenerg) )
        allocate( bb(nenerge) )
        allocate( Energe(nenerge) )
        allocate( Ephoton(nenerg) )
        allocate( e1(nenerge) )
        allocate( e2(nenerge) )
        Ephoton(:) = 0._db
        allocate( lori(nenerge) )
        allocate( lorix(nenerge) )
        allocate( lorrx(nenerge) )
        if( Dafs .or. Stokes_xan ) allocate( lorr(nenerge) )

! Calculation of the coefficients of the lorentzian
        call cflor(bb,betalor,E_cut,Energe,Elor,Energ,ie1,ie2,nef,nelor,nenerg,nenerge,XES)
        if( Energphot ) then
          Ephoton(:) = Energ(:) + Shift_init(init,ifich)
        else
          Ephoton(:) = Energ(:) + Shift_init(init,ifich) + Esmin
        endif
  
  ! Temperature effect in the Debye model
        if( ( U_iso_man .and. Abs_U_iso_inp > eps10 ) .or. ( Abs_U_iso(ifich) > eps10 .and. .not. U_iso_man ) ) then

        select case( ninit(ifich) )
          case(1)
            n = ninit(ifich)
          case(2)
            if( ninit1 == 2 ) then
              n = ninit(ifich)
            else
              n = 1
            endif
          case(4,6,10)
            if( ninit1 /= ninit(ifich) ) then
              if( init <= ninit1 ) then 
                n = ninit1
              else
                n = ninit(ifich) - ninit1 
              endif 
            else
              n = ninit(ifich)
            endif
          end select

 ! natomsym is real
          fac = natomsym / n

          call Debye_effect(Abs_U_iso,Abs_U_iso_inp,Energ,Energphot,Ephoton,Eseuil,fac,ifich,jseuil,n_col,nef,nenerg,nfich, &
                        nom_col,Numat,nxan,Shift_U_iso,U_iso_man,V0muf,Xanes)
          fac = 1._db / n

          if( Dafs ) call Debye_effect_a(Abs_U_iso,Abs_U_iso_inp,Adafs,Bragg_abs_T,Energ,Energphot,Ephoton,Eseuil,fac,ifich, &
                        jseuil,nenerg,nfich,.true.,nphim,npldafs,Numat,Shift_U_iso,U_iso_man,V0muf) 

          fac = natomsym * 100 / ( n * Volume_maille )
  
          if( Cor_abs ) then
            allocate( Mu_t(nenerg,nphim,npldafs) )
            allocate( dph_t(npldafs) )
            dph_t(:) = ( 1._db, 0._db )
            do i = 1,2
              Mu_t(:,:,:) = Mu(:,:,:,i) 
              call Debye_effect_a(Abs_U_iso,Abs_U_iso_inp,Mu_t,dph_t,Energ,Energphot,Ephoton,Eseuil,fac,ifich,jseuil, &
                        nenerg,nfich,.false.,nphim,npldafs,Numat,Shift_U_iso,U_iso_man,V0muf)
              Mu(:,:,:,i) = Mu_t(:,:,:) 
            end do
            deallocate( dph_t, Mu_t )
          endif 

          if( Stokes_xan ) then
            allocate( Mu_tt(nenerg,npldafs) )
            do i = 1,6,5
              Mu_tt(:,:) = Mu_mat(:,i,:) 
              call Debye_effect(Abs_U_iso,Abs_U_iso_inp,Energ,Energphot,Ephoton,Eseuil,fac,ifich,jseuil,n_col,nef,nenerg,nfich, &
                        nom_col,Numat,npldafs,Shift_U_iso,U_iso_man,V0muf,Mu_tt)
              Mu_mat(:,i,:) = Mu_tt(:,:)  
            end do  
            deallocate( Mu_tt )
          endif

        endif

        if( no_extrap ) then
          Extrap = .false.
        else
          Extrap = Stokes_xan .or. Cor_abs
          if( Dafs .and. .not. Extrap ) then
            do ipl = 1,npldafs
              if( abs( Bragg_abs_T(ipl) ) < 1.e-10_db ) cycle
              Extrap = .true.
              exit
            end do
          endif
        endif

        fprime_atom = icheck > 2 .and. ifich == 1 .and. i_conv == 0 .and. ( init == 1 .or. init == num_core(1) )
        if( Extrap .or. fprime_atom .or. Cor_abs .or. Stokes_xan ) then
          call Extrapat(bb(nenerg),dampl,Ephoton,Eseuil(ifich),Extrap,fpp0,fprime_atom,icheck,nenerg,Numat(ifich))
           dampl(:) = dampl(:) / ninit(ifich)
           fpp0 = fpp0 / ninit(ifich)
        endif

! Convolution by lorentzian
        if( Dafs ) Ad(1:nenerg,:,1:npldafs) = Adafs(1:nenerg,:,1:npldafs)
        if( Cor_abs ) then
          do i = 1,2
            do ipl = 1,npldafs
              Ad(1:nenerg,:,i*npldafs+ipl) = mu(1:nenerg,:,ipl,i)
            end do
          end do
        endif

        do ipl = 1,nxan
          Xa(1:nenerg,ipl) = Xanes(1:nenerg,ipl)
! Extrapolation
          do ie = nenerg+1,nenerge
            Xa(ie,ipl) = Xa(nenerg,ipl)
          end do
        end do

        gamma = .false.
        do ie = nef,nenerge
          if( abs( bb(ie) ) > 1.e-10_db ) then
            gamma = .true.
            exit
          endif
        end do

        ne2 = ie2

        do ie = ie1,ie2
          if( ie == ie1 ) then
            if( XES ) then
              e1(ie) = 1.5 * Energe(1) - 0.5 * Energe(2)
            else
              e1(ie) = E_cut
              if( ie == 1 ) then
                e1m = 1.5 * Energe(1) - 0.5 * Energe(2)
                e1(ie) = max( e1(ie), e1m )
              endif
            endif
          else
            e1(ie) = 0.5 * ( Energe(ie) + Energe(ie-1) )
          endif
          if( ie == ne2 ) then
            e2(ie) = 1.5 * Energe(ne2) - 0.5 * Energe(ne2-1)
          else
            e2(ie) = 0.5 * ( Energe(ie) + Energe(ie+1) )
          endif
        end do

        if( XES ) then
          nen2 = ie2
        else
          nen2 = nenerg
        endif
   
        if( Gamma .and. .not. Green_int ) then
   
          do ie = 1,nenerg
   
            do j = ie1,ne2
   
              if( .not. Gamma_var ) then
                bba = bb(ie)
                if( bba >= 0._db ) then
                  bbb = max( bba, 1.E-08_db )
                else
                  bbb = bba
                endif
                de2 = ( e2(j) - Energe(ie) ) / bbb
                de1 = ( e1(j) - Energe(ie) ) / bbb
                lorix(j) = atan( de1 ) - atan( de2 )
                if( ( Dafs .or. Stokes_xan ) .and. j <= nenerg ) lorrx(j) = 0.5 * log( (1 + de1**2) / (1 + de2**2) )
              endif
              if( Dafs .or. Gamma_var ) then
                bba = bb(j)
                if( bba >= 0._db ) then
                  bbb = max( bba, 1.E-08_db )
                else
                  bbb = bba
                endif
                de2 = ( e2(j) - Energe(ie) ) / bbb
                de1 = ( e1(j) - Energe(ie) ) / bbb
                lori(j) = atan( de1 ) - atan( de2 )
                if( Gamma_var ) then
                  lorix(j) = lori(j)
                  lorrx(j) = lorr(j)
                endif
                if( ( Dafs .or. Stokes_xan ) .and. j <= nenerg ) lorr(j) = 0.5 * log( (1 + de1**2) / (1 + de2**2) )
              endif
            end do

            do ipl = 1,nxan
              Xanes(ie,ipl) = - sum( lorix(ie1:ne2) * Xa(ie1:ne2,ipl) ) / pi
            end do

            if( Stokes_xan ) then
              do ipl = 1,n_mat_pol
                Mu_m(ie,1,ipl) = sum( cmplx( lorrx(ie1:nen2), lorix(ie1:nen2), db) * Mu_mat(ie1:nen2,1,ipl) ) / pi
                Mu_m(ie,2,ipl) = sum( cmplx( lorrx(ie1:nen2), lorix(ie1:nen2), db) &
                                    * cmplx( Mu_mat(ie1:nen2,2,ipl), Mu_mat(ie1:nen2,3,ipl), db ) ) / pi
                Mu_m(ie,3,ipl) = sum( cmplx( lorrx(ie1:nen2), lorix(ie1:nen2), db) &
                                    * cmplx( Mu_mat(ie1:nen2,4,ipl), Mu_mat(ie1:nen2,5,ipl), db ) ) / pi
                Mu_m(ie,4,ipl) = sum( cmplx( lorrx(ie1:nen2), lorix(ie1:nen2), db) * Mu_mat(ie1:nen2,6,ipl) ) / pi
              end do
            endif
     
            do ipl = 1,npldafs
              do ip = 1,nphi(ipl)

                if( Tenseur ) then

 ! Here Ad is in number of electron, the convolution needs division by pi
                  Adafs(ie,ip,ipl) = sum( cmplx( lorr(ie1:nen2), lori(ie1:nen2),db ) * real( Ad(ie1:nen2,ip,ipl),db ) ) / pi

                else

                  Adafs(ie,ip,ipl) = sum( cmplx(lorr(ie1:nen2),lori(ie1:nen2),db) * Ad(ie1:nen2,ip,ipl) )

! Here the imaginary part becomes the Kramers Koenig of the real part
                  if( Cor_abs ) then
                    do i = 1,2
                      jpl = ipl + i*npldafs
                      mu(ie,ip,ipl,i) = sum( cmplx( lorrx(ie1:nen2), lorix(ie1:nen2), db) * Ad(ie1:nen2,ip,jpl) ) / pi
                    end do
                  endif

                endif

              end do
            end do

          end do

        else

          if( Dafs .or. Stokes_xan ) then

            do ie = 1,nenerg

              do j = ie1,nen2
                de2 = e2(j) - Energe(ie)
                de1 = e1(j) - Energe(ie)
                if( j == ie ) then
                  if( abs( de2 + de1 ) < 1.e-10_db ) then
                    lorr(j) = 0._db
                  elseif( abs( de1 ) < 1.e-10_db ) then
                    lorr(j) = - log( 4._db )
                  else
                    lorr(j) = log( - de1 / de2 )
                  endif
                else
                  lorr(j) = log( de1 / de2 )
                endif
              end do

              do ipl = 1,npldafs
                do ip = 1,nphi(ipl)

                  if( Tenseur .and. .not. Green_int ) then
 ! Here Ad is in number of electron, the convolution needs division by pi
                    Adafs(ie,ip,ipl) = sum( lorr(ie1:nen2) * real( Ad(ie1:nen2,ip,ipl),db ) ) / pi
                    if( ie >= ie1 ) Adafs(ie,ip,ipl) = Adafs(ie,ip,ipl) - img * real( Ad(ie,ip,ipl),db )
                  else
                    if( .not. Green_int ) then
                      Adafs(ie,ip,ipl) = sum( lorr(ie1:nen2) * Ad(ie1:nen2,ip,ipl) )
                      if( ie >= ie1 ) Adafs(ie,ip,ipl) = Adafs(ie,ip,ipl) + img * pi * Ad(ie,ip,ipl)
                    endif
! For the correction by absorption, the second xanes becomes the Kramers Koenig of the frst one
                    if( Cor_abs ) then
                      do i = 1,2
                        jpl = ipl + i*npldafs
                        mu(ie,ip,ipl,i) = sum( lorr(ie1:nen2) * Ad(ie1:nen2,ip,jpl)) / pi
                        if( ie >= ie1 ) mu(ie,ip,ipl,i) = mu(ie,ip,ipl,i) - img * Ad(ie,ip,jpl)
                      end do
                    endif
                  endif
                end do
              end do

              if( Stokes_xan ) then
                do ipl = 1,n_mat_pol
                  Mu_m(ie,1,ipl) = sum( lorr(ie1:nen2) * Mu_mat(ie1:nen2,1,ipl) ) / pi
                  if( ie >= ie1 ) Mu_m(ie,1,ipl) = Mu_m(ie,1,ipl) - img * Mu_mat(ie,1,ipl)
                  Mu_m(ie,2,ipl) = sum( lorr(ie1:nen2) &
                             * cmplx( Mu_mat(ie1:nen2,2,ipl), Mu_mat(ie1:nen2,3,ipl), db ) ) / pi
                  if( ie >= ie1 ) Mu_m(ie,2,ipl) = Mu_m(ie,2,ipl) - img * cmplx( Mu_mat(ie,2,ipl), Mu_mat(ie,3,ipl), db )
                  Mu_m(ie,3,ipl) = sum( lorr(ie1:nen2) &
                             * cmplx( Mu_mat(ie1:nen2,4,ipl), Mu_mat(ie1:nen2,5,ipl), db ) ) / pi
                  if( ie >= ie1 ) Mu_m(ie,3,ipl) = Mu_m(ie,3,ipl) - img * cmplx( Mu_mat(ie,4,ipl), Mu_mat(ie,5,ipl), db )
                  Mu_m(ie,4,ipl) = sum( lorr(ie1:nen2) * Mu_mat(ie1:nen2,6,ipl) ) / pi
                  if( ie >= ie1 ) Mu_m(ie,4,ipl) = Mu_m(ie,4,ipl) - img * Mu_mat(ie,6,ipl)
                end do
              endif
          
            end do

          endif

          if( XES ) then
            do ipl = 1,nxan
              Xanes(nef+1:nenerg,ipl) = 0._db
            end do
          elseif( .not. Green_int ) then
            do ipl = 1,nxan
              Xanes(1:nef-1,ipl) = 0._db
            end do
          endif
          if( Stokes_xan ) then
            do ipl = 1,n_mat_pol
              Mu_m(1:nef-1,1:4,ipl) = 0._db
            end do
          endif
     
        endif

! Transformation to crystallo convention with f" > 0
        if( Dafs ) then
          Adafs(:,:,:) = conjg( Adafs(:,:,:) )
          if( Cor_abs ) mu(:,:,:,:) = Conjg( mu(:,:,:,:) )
        endif
        if( Stokes_xan ) Mu_m(:,:,:) = Conjg( Mu_m(:,:,:) )

! One wants f' and f" in number of electrons
        if( Tenseur_car ) then
          do ie = 1,nenerg
            Eph = Ephoton(ie)
            ct_nelec = conv_mbarn_nelec(Eph) / pi
            adafs(ie,:,:) = ct_nelec * adafs(ie,:,:)
          end do
        endif

        if( Extrap ) then
          do ipl = 1,npldafs
            if( Signal_Sph .and. ipl > 2 ) exit
            do ip = 1,nphi(ipl)
              if( ip == 1 ) then
                cf = Bragg_abs_T(ipl)
              else
                cf = Bragg_abs_T_scan(ip,ipl)
              endif
              if( Green_int ) then
                do ie = 1,nenerg
                  Adafs(ie,ip,ipl) = Adafs(ie,ip,ipl) + img * fpp0 * cf
                end do
              else
                do ie = 1,nenerg
                  Adafs(ie,ip,ipl) = Adafs(ie,ip,ipl) + ( dampl(ie) + img * fpp0 ) * cf
                end do
              endif
            end do
          end do
        endif

! One transforms to crystallo convention with f" > 0
        if( Bormann ) Adafs(:,:,:) = conjg( Adafs(:,:,:) )

        if( Cor_abs ) then
          do ipl = 1,npldafs
            if( Full_self_abs .and. ( mod(ipl,4) == 2 .or. mod(ipl,4) == 3 ) ) cycle
            do ie = 1,nenerg ! Unit cell volume is in A^3
              fac = natomsym * 100 / ( Volume_maille * Conv_mbarn_nelec(Ephoton(ie)) )
              mu(ie,:,ipl,:) = mu(ie,:,ipl,:) + fac * dampl(ie)
            end do
          end do
        endif
      
        if( Stokes_xan ) then
          do ie = 1,nenerg
            fac = natomsym * 100 / ( Volume_maille * Conv_mbarn_nelec(Ephoton(ie)) )
            Mu_m(ie,1,:) = Mu_m(ie,1,:) + fac * dampl(ie)
            Mu_m(ie,4,:) = Mu_m(ie,4,:) + fac * dampl(ie)
          end do
! Real part is absorption
        endif

! Real part is absorption
        if( Cor_abs ) mu(:,:,:,:) = img * Conjg( mu(:,:,:,:) )
        if( Stokes_xan ) Mu_m(:,:,:) = img * Conjg( Mu_m(:,:,:) )
        
! Addition of the constant term corresponding to the edges of lower energy (but on dichroic terms)
        if( i_conv == 0 .and. Abs_before .and. init == 1 ) then
          do i = 1,nxan
            nomab = nom_col(i)
            nomab = adjustl(nomab)
            if( nomab(1:3) == 'dic' .or. nomab(1:3) == 'Dic' ) cycle
            Xanes(:,i) = Xanes(:,i) + fpp_avantseuil(ifich) / n_div_fpp(ifich)
          end do
        endif

        if( ncal > 1 ) then
          iscr = 100 + jfich
          rewind( iscr )
          do ie = 1,nenerg
            write(iscr,*) Xanes(ie,1:nxan)
          end do
          if( Dafs ) then
            do ie = 1,nenerg
              write(iscr,*) ( ( Adafs(ie,ip,ipl), ip = 1,nphi(ipl) ), ipl = 1,npldafs )
            end do
          endif
        endif

        deallocate( Ad )
        deallocate( bb )
        deallocate( dampl )
        deallocate( Energe )
        deallocate( Ephoton )
        deallocate( e1 )
        deallocate( e2 )
        deallocate( lori )
        deallocate( lorix, lorrx )
        if( Dafs .or. Stokes_xan ) deallocate( lorr )
        deallocate( Xa )
        if( Seah .or. Arc ) then
          deallocate( Elor )
          deallocate( betalor )
        endif
 
      else

        iscr = 100 + jfich
        rewind( iscr )
        do ie = 1,nenerg
          read(iscr,*) Xanes(ie,1:nxan)
        end do
        if( Dafs ) then
          do ie = 1,nenerg
            read(iscr,*) ( ( Adafs(ie,ip,ipl), ip = 1,nphi(ipl) ), ipl = 1,npldafs )
          end do
        endif

      endif
      
      Energ(1:nenerg) = Energ(1:nenerg) + Shift_init(init,ifich)
      
      fac = 1._db ! (Rydb)
      
      do ie = 1,nes
        E = Es(ie)
        E = max(E ,Energ(1) )
        E = min(E ,Energ(nenerg) )
        do i = 2,nenerg
          if( Energ(i) > E - eps10 ) exit
        end do
        p1 = ( E - Energ(i-1) ) / ( Energ(i) - Energ(i-1) )
        p2 = 1 - p1
        if( Es(ie) < Energ(1) - fac .or. Es(ie) > Energ(nenerg) + fac &
                                    .or. ( Es(ie) > Energ(nenerg) + eps10 .and. Es(nes) - Es(ie) > fac ) ) then
          p1f(ie) = 0._db
          p2f(ie) = 0._db
        else  
          p1f(ie) = Pds(ifich) * p1
          p2f(ie) = Pds(ifich) * p2
        endif
        indf(ie) = i
      end do

      do ie = 1,nes
        i = indf(ie)
        if( SRXRD ) then
          if( Bulk_site(ifich) ) then
            As_SRXRD_bulk(ie,:,:,i_bulk) = As_SRXRD_bulk(ie,:,:,i_bulk) + p2f(ie) * Adafs(i-1,:,:) + p1f(ie) * Adafs(i,:,:)
            Xs_SRXRD_bulk(ie,i_bulk) = Xs_SRXRD_bulk(ie,i_bulk) + p2f(ie) * Xanes(i-1,nxan) + p1f(ie) * Xanes(i,nxan)
          else
            As_SRXRD(ie,:,:,i_surface) = As_SRXRD(ie,:,:,i_surface) + p2f(ie) * Adafs(i-1,:,:) + p1f(ie) * Adafs(i,:,:)
            Xs_SRXRD(ie,i_surface) = Xs_SRXRD(ie,i_surface) + p2f(ie) * Xanes(i-1,nxan) + p1f(ie) * Xanes(i,nxan)
          endif
        else
          if( Dafs ) As(ie,:,:) = As(ie,:,:) + p2f(ie) * Adafs(i-1,:,:) + p1f(ie) * Adafs(i,:,:)
          if( Cor_abs ) mus(ie,:,:,:) = mus(ie,:,:,:) + p2f(ie) * mu(i-1,:,:,:) + p1f(ie) * mu(i,:,:,:)

          if( Stokes_xan ) Mu_mat_comp(ie,:,:) = Mu_mat_comp(ie,:,:) + p2f(ie) * Mu_m(i-1,:,:) + p1f(ie) * Mu_m(i,:,:)

        endif
 
!        do ipl = 1,nxan
!          if( i == 2 .or. i == nenerg ) then
!            Xs(ie,ipl) = Xs(ie,ipl) + p2f(ie) * Xanes(i-1,ipl) + p1f(ie) * Xanes(i,ipl)
!          else
!            Xs(ie,ipl) = Xs(ie,ipl) + f_interp3(Es(ie),Energ(i-2),Energ(i-1),Energ(i),Energ(i+1), &
!                                                          Xanes(i-2,ipl),Xanes(i-1,ipl),Xanes(i,ipl),Xanes(i+1,ipl))       
!          endif
!        end do        
        Xs(ie,:) = Xs(ie,:) + p2f(ie) * Xanes(i-1,:) + p1f(ie) * Xanes(i,:)
      end do

      deallocate( Adafs )
      deallocate( Energ )
      deallocate( Mu_mat )
      deallocate( Mu_m )
      deallocate( Xanes )
      if( Cor_abs ) deallocate( mu )
      
    end do ! end of loop over initial states (init)

  end do ! end of loop over files

  if( .not. ( Seah .or. Arc ) ) then
    deallocate( Elor )
    deallocate( betalor )
  endif

! For SRXRD, one applies self-absorption and truncation for the Bulk
  if( SRXRD ) call Total_amplitude_SRXRD(Abs_cap,Abs_NR_SRXRD,Abs_NR_SRXRD_abs,Abs_NR_SRXRD_abs_bulk,Abs_NR_SRXRD_bulk,As, &
      As_SRXRD,As_SRXRD_bulk,Es,icheck,L_dafs,Length_abs,n_z_SRXRD,n_z_SRXRD_abs,n_z_SRXRD_abs_bulk,n_z_SRXRD_bulk,nes, &
      Non_resonant_bulk,nphim,npldafs,npldafs_NR,Pos_z_SRXRD,Pos_z_SRXRD_abs,Pos_z_SRXRD_bulk,Pos_z_SRXRD_abs_bulk, &
      Sum_Bragg_abs_f0,Sum_Bragg_abs_f0_bulk,Sum_Bragg_nonabs_f,Sum_Bragg_nonabs_f_bulk,Sum_Bragg_NR_Bulk,Sum_Bragg_NR_T, &
      Xs_SRXRD,Xs_SRXRD_bulk)
  
  if( DAFS .and. .not. ( Forbidden .or. Tenseur .or. SRXRD ) ) then
    do ipl = 1,npldafs
      do ip = 1,nphi(ipl)
        if( nphi(ipl) == 1 ) then
          if( Bormann ) then
            As(:,ip,ipl) = As(:,ip,ipl) + conjg( Sum_Bragg_NR_T(ipl) )
          else
            As(:,ip,ipl) = As(:,ip,ipl) + Sum_Bragg_NR_T(ipl)
          endif
        else
          if( Bormann ) then
            As(:,ip,ipl) = As(:,ip,ipl) + conjg( Sum_Bragg_NR_T_scan(ip,ipl) )
          else
            As(:,ip,ipl) = As(:,ip,ipl) + Sum_Bragg_NR_T_scan(ip,ipl)
          endif
        endif
      end do
    end do
  endif
  
! In fact one calculates epsilon - 1 = -4*pi*Atom_density*r0/k2
! We are in atomic units r0 --> r0 / a0 = alfa^2
! k = 1/2 alfa E
! Volume_maille = unit cell volume in A^3
! ct = -4*pi * alfa^2  * 4  * bohr**3 / ( alfa^2 * E^2 * Volume_maille )
!    = -16*pi  * bohr**3 / ( E^2 * Volume_maille )
  if( Bormann ) then
    do ie = 1,nes
      if( Energphot ) then
        Eph = Es(ie) + Esmin
      else
        Eph = Es(ie)
      endif
      ct_epsilon = - 16 * pi * bohr**3 / ( Volume_maille * Eph**2 )
      As(ie,:,:) = ct_epsilon * As(ie,:,:)
    end do
  endif

  if( Cor_abs ) then
! Conversion from Megabarn to micrometer^-1
    c_micro = 100 / Volume_maille
! mus is already in micrometre -1
    do ipl = 1,npldafs
      if( Full_self_abs .and. ( mod(ipl,4) == 2 .or. mod(ipl,4) == 3 ) ) cycle
      mus(:,:,ipl,:) = mus(:,:,ipl,:) + c_micro * cmplx( fpp_avantseuil(ifich), f0_forward, db )
    end do
  endif

  if( Stokes_xan ) then
! Conversion from Megabarn to micrometer^-1
    c_micro = 100 / Volume_maille
! Mu_mat_comp is already in micrometer -1
    do i = 1,4,3
      Mu_mat_comp(:,i,:) = Mu_mat_comp(:,i,:) + c_micro * cmplx( fpp_avantseuil(ifich), f0_forward, db)
    end do
  endif

! Convolution by a Gaussian
  if( abs(Delta_gauss) > 1.e-10_db ) then

    allocate( Yr(nes) )
    allocate( Yi(nes) )
! number of points added beyond the range in both directions    
    nj = 10

    do ipl = 1,nxan
      Yr(:) = Xs(:,ipl)
      call gaussi(Delta_gauss,Es,nes,nj,Yr)
      Xs(:,ipl) = Yr(:)
    end do

    do ipl = 1,npldafs
      do ip = 1,nphi(ipl)
        Yr(:) = real( As(:,ip,ipl), db )
        call gaussi(Delta_gauss,Es,nes,nj,Yr)
        Yi(:) = aimag( As(:,ip,ipl) )
        call gaussi(Delta_gauss,Es,nes,nj,Yi)
        As(:,ip,ipl) = cmplx( Yr(:), Yi(:),db )

        if( Cor_abs ) then
          do i = 1,2
            Yr(:) = real( mus(:,ip,ipl,i), db )
            call gaussi(Delta_gauss,Es,nes,nj,Yr)
            Yi(:) = aimag( mus(:,ip,ipl,i) )
            call gaussi(Delta_gauss,Es,nes,nj,Yi)
            mus(:,ip,ipl,i) = cmplx( Yr(:), Yi(:),db )
          end do
        endif

      end do
    end do

    deallocate( Yr )
    deallocate( Yi )

  endif

  if( Stokes_dafs ) call Cal_stokes(As,Icirc,nes,n_stokes,np_stokes,nphi,nphim,npldafs,Stokes_param)

  if( Stokes_xan ) then
    if( icheck > 1 ) then
      do ipl = 1,n_mat_pol
        write(3,189) ipl
        write(3,190) ( ipl, i = 1,8 )  
        do ie = 1,nes
          write(3,'(f10.3,1p,121e13.5)') Es(ie)*rydb, Mu_mat_comp(ie,:,ipl) 
        end do
      end do
    endif
    call Cal_MuTrans(icheck,Es,Mu_mat_comp,Mus_mat,n_bir,n_stokes,n_mat_pol,nes,Sample_thickness,Stokes_param)
  endif

  if( Cor_abs ) &
    call Corr_abs(angxyz,As,axyz,c_micro,d_dead,Double_cor,Eseuil(1),fpp_avantseuil(1),hkl_dafs,hkl_S,Icor,Icirccor, &
          Icircdcor,Idcor,mus,n_stokes,nes,np_stokes,nphi,nphim,npldafs,Self_abs,Stokes_param)

! 2D diffraction case
  if( Surface_ref > eps10 )  As(:,:,:) = As(:,:,:) / Surface_ref

!---- Writing -------------------------------------------------

    call Writing_output_conv(Ampl_Phase,Analyzer,Angle,As,Bormann,Convolution_out_L,Cor_abs,Couple_co,Dafs_bio, &
      Dafs_exp_type,Double_cor,E_cut,Es,Energ_tr,Eseuil,File_scan_out,Fit_cal,fprim,Full_self_abs, &
      Gamma_h,Gamma_max,Header,icheck,Icirc,Icirccor,Icircdcor,Icor,Idcor,iscratchconv,Mu_mat_comp, &
      Mus,Mus_mat,n_bir,n_col,nes,n_couple_weight_co,n_energ_tr,n_mat_pol,n_stokes,nfich,ninitm, &
      nom_col,np_stokes,nphi,nphim,npldafs,npldafs_b,npldafs_t,nxan,Scan_true,Shift_init,Stokes,Tenseur,Transpose_file, &
      Weight_co,Xs)

  end do boucle_conv_file
 
    deallocate( Abs_U_iso, angle, As, Bragg_abs_T, Bragg_abs_T_scan, Bragg_abs_TT, Bulk_site, Convolution_out_all, Couple_co )
    deallocate( Efermip, En_fermi, Energ_tr, Eph1, Ephm, Epsii, Es, Eseuil )
    deallocate( fi, File_in, File_scan_in, File_scan_out_all, fpp_avantseuil, fr, hkl_dafs )
    deallocate( Icirc, Icirccor, Icircdcor, Icor, Idcor, indf )
    deallocate( Mu_mat_comp,  mus, Mus_mat, n_div_fpp, natomsym_f, ne, ne_initl, ninit, nL_header, nom_col, nphi, nsup, Numat )
    deallocate( p1f, p2f, Pds, Run_done, Shift_init, Shift_site, Stokes_name, Stokes_param, Skip_run, Sum_Bragg_NR_Bulk, &
                Sum_Bragg_NR_T, Sum_Bragg_NR_T_scan, V0muf, Xs )
    if( SRXRD ) then
      deallocate( Abs_NR_SRXRD, Abs_NR_SRXRD_abs, Abs_NR_SRXRD_bulk, Abs_NR_SRXRD_abs_bulk, L_dafs, Length_abs, &
                          Pos_z_SRXRD, Pos_z_SRXRD_abs, Pos_z_SRXRD_bulk, Pos_z_SRXRD_abs_bulk )
      deallocate( Sum_Bragg_abs_f0, Sum_Bragg_abs_f0_bulk, Sum_Bragg_nonabs_f, Sum_Bragg_nonabs_f_bulk )
    endif

    if( ncal > 1 .and. ical == ncal+1 ) then
      do jfich = 1,mfich
        iscr = 100 + jfich
        Close( iscr )
      end do
    endif

    if( .not. Another_one ) exit

  end do

  return

  110 format(///' Line not understood in the indata file :'//1x,A//, &
          ' If it is supposed to be a keyword check the spelling.',/ &
          ' If it is a line containing numbers, check:',/ &
          '   - How many numbers are supposed to be there.',/ &
          '   - Numbers must be separated by spaces.',/ &
          '   - Are there extra characters (points, tabulations...)'//)
  112 format(//' Error opening the file:',//3x,A,//  &
           5x,'It does not exist but the following file is found existing: ',//3x,A,// &
           5x,'Modify the input file and eventualy add other indata files with the convenient indexes !',//)
  130 format(//' Selec_core keyword is used with the numbers =',10i3)
  135 format(' At list one of them is larger than the number of core states (or edges when non magnetic) =',i2,/ &
             ' what is forbidden !'/)
  145 format(///' Error under the keyword Par_',a6,'in the indata file:',// &
                ' The wanted file is the number',i3,' !',/ &
                ' There are only',i3,' files in the job !'//)
  150 format('    Gamma_max  =',f7.2,',  Aseah =',f7.2)
  160 format('    Gamma_max  =',f7.2,',  Ecent =',f7.2,', Elarg =',f7.2)
  170 format('    Gamma_hole =',f7.2,', E_cut =',f8.3,', Shift =',f8.3,' eV')
  172 format('    Gamma_hole =',f7.2,', E_cut =',f8.3,', Shift =',f8.3,' eV,  init =',i3)
  174 format('    Gamma_hole =',f7.2,', E_cut =',f8.3,', Shift =',f8.3,' eV,  site', i3)
  176 format('    Gamma_hole =',f7.2,', E_cut =',f8.3,', Shift =',f8.3,' eV,  site', i3,', init =',i3)
  178 format(/'      E_(eV)      Width_(eV) lambda_(A)')
  179 format(/' E_(eV) - E_edge  Width_(eV) lambda_(A)')
  180 format(3f12.3)
  185 format(2f12.3,1p,e12.4)
  189 format(/' Complex linear absorption coefficients (micrometer^-1) for the orientation matrix number =',i3)
  190 format(/'   Energy   ',' Mu_ss_r_',i1,'    Mu_ss_i_',i1,'    Mu_sp_r_',i1,'    Mu_sp_i_',i1,'    Mu_ps_r_',i1, &
                          '    Mu_ps_i_',i1,'    Mu_pp_r_',i1,'    Mu_pp_i_',i1)
  200 format(f10.3,1p,240e13.5)
  210 format(f7.1,1p,3e13.5)
end

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

subroutine Initialization_conv(Abs_before,Abs_U_iso_inp,Ampl_phase,Arc,Asea,bav_open,Check_conv,Chem,Convolution_out, &
             d_dead,Dafs,Dafs_exp_type,decferm,Delta_gauss,Deuxieme,Double_cor,E_cut,E_cut_imp,E_cut_param,Eintmax,f0_forward, &
             Fermip,File_scan_out,Forbidden,fprim,fprime_atom,Full_self_abs,Gamma_hole_param,Gamma_var,Gaussian_default,Green_int, &
             hkl_S,icheck,jseuil,Just_total,Magn,Multi_file,n_couple_weight_co,n_Stokes,nelor,nfich,npldafs,npldafs_t,npldafs_th, &
             nseuil,num_core,nxan_lib,S0_2,Sample_thickness,scan_true,seah,self_abs,Shift_U_iso,Signal_Sph,Stokes,Stokes_dafs, &
             Stokes_xan,Surface_ref,Tenseur,Tenseur_car,Thomson,Transpose_file,U_iso_man,Volume_maille,Volume_maille_bulk,XES, &
             Weight_co)

  use declarations
  implicit none
  
  integer:: Dafs_exp_type, i, icheck, jseuil, n_couple_weight_co, n_Stokes, nelor, nfich, npldafs, npldafs_t, npldafs_th, nseuil

  integer, dimension(3):: hkl_S
  integer, dimension(10):: num_core
  
  character(len=Length_name):: Convolution_out, File_scan_out
  
  logical:: Abs_before, Ampl_phase, Arc, bav_open, Check_birefringence, Check_conv, Chem, Dafs, decferm, &
    Deuxieme, Double_cor, E_cut_param, Fermip, Forbidden, fprim, fprime_atom, Full_self_abs, Gamma_hole_param, Gamma_var, &
    Gaussian_default, Green_int, Just_total, Magn, Multi_file, no_extrap, nxan_lib, scan_true, seah, self_abs, Signal_Sph, Stokes, &
    Stokes_dafs, Stokes_xan, Tenseur, Tenseur_car, Thomson, Transpose_file, U_iso_man, XES 
  
  real(kind=db):: Abs_U_iso_inp, Asea, d_dead, Delta_gauss, E_cut, E_cut_imp, Eintmax, f0_forward, S0_2, Sample_thickness, &
    Shift_U_iso, Surface_ref, Volume_maille, Volume_maille_bulk, Weight_co

  Abs_before = .false.
  Abs_U_iso_inp = 0._db
  Ampl_phase = .false.
  Arc = .true.
  Asea = 0.2_db  ! Slope of Gamma at the origin in Seah Dench model
  Check_birefringence = .false.
  chem = .false.
  Check_conv = .false.
  Convolution_out = ' '
  d_dead = 0._db
  Dafs = .false.
  Dafs_exp_type = 4
  decferm = .false.
  Delta_gauss = 0._db
  Deuxieme = .false.
  Double_cor = .false.
  E_cut = E_cut_imp * rydb
  E_cut_param = .false.
  Eintmax = 1000000._db
  f0_forward = 0._db
  Fermip = .false.
  File_scan_out = ' '
  Forbidden = .false.
  fprim = .false.
  fprime_atom = .false.
  Full_self_abs = .false.
  Gamma_hole_param = .false.
  Gamma_var = .false.
  Gaussian_default = .false.
  Green_int = .false.
  hkl_S(:) = 0
  jseuil = 1
  Just_total = .true.
  Magn = .false.
  Multi_file = .false.
  n_couple_weight_co = 0
  n_Stokes = 0
  nelor = 0
  nfich = 0
  npldafs = 0
  npldafs_t = 0
  npldafs_th = 0
  no_extrap = .false. 
  nseuil = -1
  do i = 1,10
    num_core(i) = i
  end do
  nxan_lib = .false.
  S0_2 = 1._db
  Sample_thickness = 0._db
  scan_true = .false.
  seah = .false.
  self_abs = .false.
  Signal_Sph = .false.
  Shift_U_iso = 0._db
  Stokes = .false.
  Stokes_dafs = .false.
  Stokes_xan = .false.
  Surface_ref = 0._db
  Tenseur = .false.
  Tenseur_car = .false.
  Thomson = .false.
  Transpose_file = .false.
  U_iso_man = .false.
  Volume_maille = 0._db
  Volume_maille_bulk = 0._db
  XES = .false.
  Weight_co = 0.5_db

  if( bav_open .or. icheck > 1 ) Check_conv = .true.
  
  return
end

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

! Dimension of the tables for lecture

subroutine Dimension_lect_conv(Circular,File_name_temp,itape1,kw_conv,Length_line,Multi_file,n_couple_weight_co,n_energ_tr, &
                   n_Stokes,nelor,nfich,nkw_conv,npldafs_th)

  use declarations
  implicit none
  
  integer:: eof, i, ie, ifich, ii, istat, itape1, j, Length, Length_line, n, nelor, n_couple_weight_co, n_energ_tr, n_Stokes, &
    nfich, nkw_conv, nnombre, npldafs_th

  character(len=9):: keyword, mot9
  character(len=132):: identmot, mot
  character(len=Length_name):: File_name, Word, Word_b
  character(len=9), dimension(nkw_conv):: kw_conv
  character(len=Length_name), dimension(10000):: File_name_temp

  logical:: Circular, Multi_file
  
  n_energ_tr = 0
      
  Rewind(itape1)

  boucle_ii: do ii = 1,1000000

    n = nnombre(itape1,132)
    read(itape1,'(A)',iostat=eof) mot
    if( eof /= 0 ) exit boucle_ii

    keyword = identmot(mot,9)

    select case(keyword)

      case('calculati')

        nfich = 0
        n = nnombre(itape1,132)
        boucle_i: do i = 1,1000000
          n = nnombre(itape1,132)
          read(itape1,'(A)',iostat=eof) File_name
          if( eof /= 0 ) exit boucle_i
          mot9 = identmot(File_name,9)
          do j = 1,nkw_conv
            if( mot9 /= kw_conv(j) ) cycle
            backspace(itape1)
            exit boucle_i
          end do
         if( mot9 == 'run_done' ) exit boucle_i
          if( n == 0 ) then
            nfich = nfich + 1
            Length = len_trim( File_name )
            if( Length > 4 ) then
              if( File_name(Length-3:Length) /= '.txt' ) File_name(Length+1:Length+4) = '.txt'
              File_name_temp(nfich) = File_name
            endif
          endif
       end do boucle_i

! looking for several indata files with same main name      
        if( nfich == 1 ) then
          Rewind( itape1 )
          do
            read(itape1,'(A)') mot
            keyword = identmot(mot,9)
            if( keyword /= 'calculati' ) cycle
            n = nnombre(itape1,132)
            read(itape1,'(A)',iostat=eof) File_name
            File_name = adjustl(File_name)
            Word = File_name
            Length = len_trim( File_name )
            if( Length > 4 ) then
              if( File_name(Length-3:Length) /= '.txt' ) File_name(Length+1:Length+4) = '.txt'
            endif
            open(2, file = File_name, status='old',iostat=istat)
            if( istat == 0 ) exit
            Word(Length+1:Length+1) = '_'
            ifich = 0
            do i = 1,10000
              File_name = Word
              call ad_number(i,File_name,Length_name)
              Length = len_trim( File_name )
              File_name(Length+1:Length+4) = '.txt'
              open(2, file = File_name, status='old',iostat=istat)
              if( istat /= 0 ) then
  ! for bulk case
                Word_b = File_name
                Word_b(Length+1:Length+4) = '_   '  
                do j = 1,10
                  File_name = Word_b 
                  call ad_number(j,File_name,Length_name)
                  Length = len_trim( File_name )
                  File_name(Length+1:Length+4) = '.txt'
                  open(2, file = File_name, status='old',iostat=istat)
                  if( istat /= 0 ) then
                    cycle
                  else
                    Multi_file = .true.
                    ifich = ifich + 1 
                    File_name_temp(ifich) = File_name
                  endif
                end do
                cycle
              endif
              Multi_file = .true.
              ifich = ifich + 1 
              File_name_temp(ifich) = File_name
            end do
            nfich = ifich
            exit
          end do
          if( nfich == 0 ) then
            Length = len_trim(Word)
            Word(Length:Length) = ' '
            call write_open_error(Word,istat,1)
          endif
        endif

      case('table')

        n = nnombre(itape1,132)
        if( n == 1 ) read(itape1,*)
        nelor = 0
        do ie = 1,10000
          n = nnombre(itape1,132)
          if( n == 0 ) exit
          read(itape1,*)
          nelor = nelor + 1
        end do

      case('stokes')
      
        do i = 1,1000000
          n = nnombre(itape1,132)
          if( n == 0 ) exit
          n_Stokes = n_Stokes + 1
          read(itape1,*)
        end do

      case('thomson')

        npldafs_th = nnombre(itape1,Length_line) / 2

      case('transpose')
        n_energ_tr = nnombre(itape1,132)
         
      case('weight_co')
        n = nnombre(itape1,132)
        if( n == 0 ) then
          call write_err_form(itape1,keyword)
        elseif( n == 1 ) then
          read(itape1,*,iostat=eof)
          n = nnombre(itape1,132)
          if( n == 0 .or. mod(n,2) == 1 ) call write_err_form(itape1,keyword)
        else
          n = nnombre(itape1,132)
        endif
        n_couple_weight_co = n / 2

    end select

  end do boucle_ii

  if( Circular ) n_Stokes = n_Stokes + 4

  return
end

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

! Lecture of the parameters for convolution

subroutine Lecture_conv( Abs_before, Abs_U_iso_inp, Ampl_phase, Another_one, Arc,asea, Betalor, Check_birefringence, Chem, &
         Convolution_out, Couple_co, d_dead, Dafs_exp_type, Decferm, Delta_gauss, Deuxieme, Double_cor, E_cut, E_cut_man, &
         E_step_conv, Ecent, Efermip, Eintmax, Elarg, Elor, Energ_tr, f0_th, Fermip, File_in, File_name, File_name_temp, &
         File_scan_in, File_scan_out, Forbidden, fprim, Gaussian_default, Gamma_hole, Gamma_hole_man, Gamma_max, Gamma_max_man, &
         Gamma_var, hkl_S, itape1, Just_total, kw_conv, mfich, Multi_file, n_couple_weight_co, n_energ_tr, n_selec_core, &
         n_stokes, nelor, nfich, nkw_conv, no_extrap, npldafs_th, num_core, nxan_lib, Pds, Run_done, S0_2, Sample_thickness, &
         Scan_a, Scan_true, Seah, Shift_site, Shift_U_iso, Step_conv, Stokes, Stokes_name, Stokes_param, Thomson, Transpose_file, &
         U_iso_man, Weight_co, XES )

  use declarations
  implicit none
  
  integer:: Dafs_exp_type, eof, i, ii, ie, ifich, ipl, ipr, istat, itape1, j, L, Length, mfich, n, n_couple_weight_co, &
    n_energ_tr, n_selec_core, n_stokes, nelor, nfich, nkw_conv, nnombre, npldafs_th

  integer, dimension(3):: hkl_S
  integer, dimension(10):: num_core
  integer, dimension(n_couple_weight_co,2):: Couple_co

  character(len=9):: keyword, mot9, Traduction
  character(len=132):: identmot, mot, mots
  character(len=Length_name):: Chemin, Convolution_out, Word, File_name, File_scan_out
  character(len=9), dimension(nkw_conv):: kw_conv
  character(len=Length_name), dimension(10000):: File_name_temp
  character(len=Length_name), dimension(nfich):: File_in, File_scan_in
  character(len=Length_word), dimension(n_stokes):: Stokes_name

  complex(kind=db), dimension(npldafs_th):: f0_th

  logical:: Abs_before, Ampl_phase, Another_one, Arc, Check_birefringence, Chem, Decferm, Deuxieme, Double_cor, E_cut_man, &
    Fermip, Forbidden, fprim, Gamma_hole_man, Gamma_max_man, Gamma_var, Gaussian_default, Just_total, Multi_file, &
    No_extrap, nxan_lib, Scan_a, Scan_true, Seah, Stokes, Step_conv, Thomson, Transpose_file, U_iso_man, XES
  
  logical, dimension(nfich):: Run_done
  
  real(kind=db):: Abs_U_iso_inp, asea, d_dead, Delta_gauss, E_cut, E_step_conv, Ecent, Eintmax, Elarg, Gamma_max, S0_2, &
    Sample_thickness, Shift_U_iso, Weight_co 

  real(kind=db), dimension(10):: Gamma_hole
  real(kind=db), dimension(nfich):: Efermip, Pds, Shift_site
  real(kind=db), dimension(nelor):: Betalor, Elor
  real(kind=db), dimension(npldafs_th):: fi, fr
  real(kind=db), dimension(n_energ_tr):: Energ_tr
  real(kind=db), dimension(5,n_Stokes):: Stokes_param
  
  E_step_conv = -1._db
  
  no_extrap = .false.
  Step_conv = .false.
  
  rewind(itape1)

  boucle_lect: do ii = 1,1000000

    n = nnombre(itape1,132)
    read(itape1,'(A)',iostat=eof) mot
    if( eof /= 0 ) exit boucle_lect

    keyword = identmot(mot,9)
    keyword = Traduction(keyword)

    select case(keyword)

      case('run_done')
        exit

      case('abs_befor')
        Abs_before = .true.

      case('all_conv')
        Just_total = .false.

      case('amplphase')
        Ampl_phase = .true.

      case('calculati')

        if( Multi_file ) then

          n = nnombre(itape1,132)
          read(itape1,*)

          Pds(:) = 1._db
          Shift_site(:) = 0._db

          do ifich = 1,nfich
            File_in(ifich) = File_name_temp(ifich)
          end do
        
        else

          ifich = 0
          do i = 1,mfich
            n = nnombre(itape1,132)
           
            if( run_done(i) ) then
              read(itape1,*)
              n = nnombre(itape1,132)
              if( n /= 0 ) read(itape1,*)
              cycle
            endif
            ifich = ifich + 1
            read(itape1,'(A)') File_in(ifich)
            File_in(ifich) = adjustl( File_in(ifich) )
            n = nnombre(itape1,132)
            if( n == 0 ) then
              Pds(ifich) = 1._db
              Shift_site(ifich) = 0._db
            elseif( n == 1 ) then
              read(itape1,*) Pds(ifich)
              Shift_site(ifich) = 0._db
            elseif( n == 2 ) then
              read(itape1,*) Pds(ifich), Shift_site(ifich)
            else
              Fermip = .true.
              read(itape1,*) Pds(ifich), Shift_site(ifich), Efermip(ifich)
            endif
          end do
          
        endif

      case('cal_tddft')

        ifich = 0
        do i = 1,mfich
          n = nnombre(itape1,132)
          if( run_done(i) ) then
            read(itape1,*)
            n = nnombre(itape1,132)
            if( n /= 0 ) read(itape1,*)
            cycle
          endif
          ifich = ifich + 1
          if( Another_one ) then
            read(itape1,'(A)') File_in(ifich)
            File_in(ifich) = adjustl( File_in(ifich) )
            Deuxieme = .true.
          else
            read(itape1,'(A)') mot
          endif
          Pds(ifich) = 1._db
          Shift_site(ifich) = 0._db
        end do
        Another_one = .not. Another_one

      case('conv_out')

        n = nnombre(itape1,132)
        read(itape1,'(A)') convolution_out
        convolution_out = adjustl( convolution_out )

      case('dafs_exp_')

        n = nnombre(itape1,132)
        read(itape1,*) Dafs_exp_type

      case('e_step_co')

        n = nnombre(itape1,132)
        Step_conv = .true.
        if( n > 0 ) read(itape1,*) E_step_conv

      case('gaussian')

        n = nnombre(itape1,132)
        if( n == 1 .or. n == 2 ) then
          read(itape1,*) Delta_gauss
        else
          Gaussian_default = .true.
        endif

      case('fprime')

        fprim = .true.

      case('forbidden')

        Forbidden = .true.

      case('scan_file')

        scan_true = .true.
        do ifich = 1,nfich
          n = nnombre(itape1,132)
          read(itape1,'(A)') File_scan_in(ifich)
          File_scan_in(ifich) = adjustl( File_scan_in(ifich) )

          open(2, file = File_scan_in(ifich), status='old',iostat=istat)
          if( istat /= 0 ) then
            Word = File_scan_in(ifich)
            Length = len_trim(Word)
            if( Word(Length-3:Length) /= '.txt' ) then
              Word(Length+1:Length+4) = '.txt'
              Close(2) 
              open(2, file = Word, status='old', iostat=istat)
              if( istat /= 0 ) then
                Length = len_trim(Word)
                Word(Length-3:Length+2) = '_1.txt'
                open(2, file = Word, status='old', iostat=istat)
                if( istat == 0 ) then
                  call write_error
                  do ipr = 6,9,3
                    write(ipr,112) File_in(ifich), Word
                  end do
                  stop
                else
                  call write_open_error(File_in(ifich),istat,1)
                endif
              endif
            endif
            File_scan_in(ifich) = Word
          endif
        end do

      case('scan')

        Scan_a = .true.

      case('scan_conv')

        n = nnombre(itape1,132)
        read(itape1,'(A)') File_scan_out
        File_scan_out = adjustl( File_scan_out )

      case('directory')

        Chem = .true.
        n = nnombre(itape1,132)
        read(itape1,'(A)') chemin

      case('abs_b_iso')
        U_iso_man = .true.
        n = nnombre(itape1,132)
        if( n == 1 ) then
          read(itape1,*) Abs_U_iso_inp
        else
          read(itape1,*) Abs_U_iso_inp, Shift_U_iso          
        endif
        Abs_U_iso_inp = Abs_U_iso_inp / ( 8 * pi**2 ) 

      case('abs_u_iso')
        U_iso_man = .true.
        n = nnombre(itape1,132)
        if( n == 1 ) then
          read(itape1,*) Abs_U_iso_inp
        else
          read(itape1,*) Abs_U_iso_inp, Shift_U_iso          
        endif

      case('check_bir')
        Check_birefringence = .true.
        
      case('seah')

        Arc = .false.
        Seah = .true.
        n = nnombre(itape1,132)
        select case(n)
          case(0)
            continue
          case(1)
            read(itape1,*) asea
          case(2)
            read(itape1,*) asea, Gamma_max
            Gamma_max = Gamma_max / rydb
            Gamma_max_man = .true.
          case(3)
            read(itape1,*) asea, Gamma_max, Gamma_hole(1)
            Gamma_hole(1) = Gamma_hole(1) / rydb
            Gamma_max = Gamma_max / rydb
            Gamma_hole_man = .true.
            Gamma_max_man = .true.
          case default
            read(itape1,*) asea, Gamma_max, Gamma_hole(1), E_cut
            Gamma_hole(1) = Gamma_hole(1) / rydb
            Gamma_max = Gamma_max / rydb
            Gamma_hole_man = .true.
            Gamma_max_man = .true.
            E_cut_man = .true.
        end select

      case('convoluti')

        Arc = .true.
        n = nnombre(itape1,132)
        select case(n)
          case(0)
            continue
          case(1)
            read(itape1,*) Ecent
            Ecent = Ecent / rydb
            Elarg = Elarg / rydb
          case(2)
            read(itape1,*) Ecent, Elarg
            Ecent = Ecent / rydb
            Elarg = Elarg / rydb
          case(3)
            read(itape1,*) Ecent, Elarg, Gamma_max
            Ecent = Ecent / rydb
            Elarg = Elarg / rydb
            Gamma_max = Gamma_max / rydb
            Gamma_max_man = .true.
          case(4)
            read(itape1,*) Ecent, Elarg, Gamma_max, Gamma_hole(1)
            Ecent = Ecent / rydb
            Elarg = Elarg / rydb
            Gamma_hole(1) = Gamma_hole(1) / rydb
            Gamma_max = Gamma_max / rydb
            Gamma_hole_man = .true.
            Gamma_max_man = .true.
         case default
            read(itape1,*) Ecent, Elarg, Gamma_max, Gamma_hole(1), E_cut
            Ecent = Ecent / rydb
            Elarg = Elarg / rydb
            Gamma_hole(1) = Gamma_hole(1) / rydb
            Gamma_max = Gamma_max / rydb
            Gamma_hole_man = .true.
            Gamma_max_man = .true.
            E_cut_man = .true.
        end select

      case('table')

        Arc = .false.
        n = nnombre(itape1,132)
        if( n == 1 ) then
          read(itape1,*) E_cut
          E_cut_man = .true.
        endif
        do ie = 1,nelor
          n = nnombre(itape1,132)
          read(itape1,*,iostat=eof) Elor(ie), betalor(ie)
          if( eof > 0 ) call write_err_form(itape1,keyword)
        end do
        Elor(:) = Elor(:) / rydb
        betalor(:) = betalor(:) / rydb

      case('eintmax')

        n = nnombre(itape1,132)
        read(itape1,*,iostat=eof) Eintmax
        if( eof > 0 ) call write_err_form(itape1,keyword)

      case('gamma_fix')
        Gamma_var = .false.

      case('gamma_var')
        gamma_var = .true.

! To make the shift before the convolution.
      case('dec')

        Decferm = .true.

      case('no_extrap')

        no_extrap = .true.

      case('nxan_lib')

        nxan_lib = .true.

      case('thomson')

        Thomson = .true.
        read(itape1,*,iostat=eof) (fr(ipl), fi(ipl),ipl = 1,npldafs_th)
        if( eof > 0 ) call write_err_form(itape1,keyword)
        f0_th(:) = cmplx( fr(:), fi(:),db)

      case('transpose')
        Transpose_file = .true.
        read(itape1,*) Energ_tr(:)
        Energ_tr(:) = Energ_tr(:) / rydb
         
      case('s0_2')
        read(itape1,*,iostat=eof) S0_2

      case('XES')
        XES = .true.

      case('selec_cor')
        n_selec_core = nnombre(itape1,132)
        read(itape1,*,iostat=eof) num_core(1:n_selec_core)
        if( eof > 0 ) call write_err_form(itape1,keyword)

      case('surface_p')
        n = nnombre(itape1,132)
        read(itape1,*,iostat=eof) hkl_S(:)
        if( eof > 0 ) call write_err_form(itape1,keyword)

      case('sample_th')
        n = nnombre(itape1,132)
        read(itape1,*,iostat=eof) Sample_thickness

      case('stokes')
        Stokes = .true.
        do i = 1,n_stokes
          n = nnombre(itape1,132)
          if( n == 0 ) exit
          read(itape1,*,iostat=eof) Stokes_param(1:min(5,n),i)
          if( eof > 0 ) call write_err_form(itape1,keyword)
!Si l'angle de rotation de l'analyseur est defini mais pas son angle de Bragg,
!ce dernier est supposé = 45 degres correspondant a un analyseur "parfait"
          if( n == 4 ) Stokes_param(5,i) = 45._db
          Stokes_param(4:5,i) = Stokes_param(4:5,i) * pi / 180
        end do

      case('stokes_na')
        boucle_n: do i = 1,n_stokes
          read(itape1,'(A)',iostat=eof) mots
          if( eof /= 0 ) exit boucle_lect
          mot9 = identmot(mots,9)
          do j = 1,nkw_conv
            if( mot9 /= kw_conv(j) ) cycle
            backspace(itape1)
            exit boucle_n
          end do
          mots = adjustl( mots )
          Stokes_name(i) = mots(1:Length_word)
          if( Stokes_name(i) == 'noname') Stokes_name(i) = 'no_name'
        end do boucle_n

      case('dead_laye')
        n = nnombre(itape1,132)
        read(itape1,*,iostat=eof) d_dead
        if( eof > 0 ) call write_err_form(itape1,keyword)
        d_dead = d_dead / 10000

      case('double_co')
        Double_cor = .true.
        
      case('weight_co')
        n = nnombre(itape1,132)
        if( n == 0 ) then
          call write_err_form(itape1,keyword)
        elseif( n == 1 ) then
          read(itape1,*,iostat=eof) Weight_co
          n = nnombre(itape1,132)
          if( n == 0 .or. mod(n,2) == 1 ) call write_err_form(itape1,keyword)
        else
          n = nnombre(itape1,132)
        endif
        n_couple_weight_co = n / 2
        if( mod(n,2) == 1 ) then
          read(itape1,*,iostat=eof) Weight_co, ( Couple_co(i,:), i = 1, n_couple_weight_co )
        else
          read(itape1,*,iostat=eof) ( Couple_co(i,:), i = 1, n_couple_weight_co )
        endif
                                
      case default

        call write_error
        do ipr = 6,9,3
          write(ipr,110) mot
        end do
        stop

    end select

  end do boucle_lect

  if( Chem ) then
    L = len_trim(Chemin)
    Word = File_in(1)
    if( Chemin(L:L) /= '/' .and. File_in(1)(1:1) /= '/' ) then
      L = L + 1
      Chemin(L:L) = '/'
    endif
    do ifich = 1,nfich
      Word = File_in(ifich)
      Length = len_trim(Word)
      File_in(ifich) = chemin(1:L) // Word(1:Length)
    end do
  endif
  
! Test on input file names
  do ifich = 1,nfich
    open(2, file = File_in(ifich), status='old', iostat=istat)
    if( istat /= 0 ) then
      File_name = File_in(ifich)

      Length = len_trim(File_name)
      if( Length <= 4 ) call write_open_error(File_in(ifich),istat,1)

      if( File_name(Length-3:Length) /= '.txt' ) then
        File_name(Length+1:Length+4) = '.txt'
        Close(2) 
        open(2, file = File_name, status='old', iostat=istat)
        if( istat /= 0 ) then
          Length = len_trim(File_name)
          File_name(Length-3:Length+2) = '_1.txt'
          open(2, file = File_name, status='old', iostat=istat)
          if( istat == 0 ) then
           call write_error
           do ipr = 6,9,3
             write(ipr,112) File_in(ifich), File_name
           end do
           stop
         else
           call write_open_error(File_in(ifich),istat,1)
         endif
        endif
        File_in(ifich) = File_name
      else
        Length = len_trim(File_name)
        File_name(Length-3:Length+2) = '_1.txt'
        open(2, file = File_name, status='old', iostat=istat)
        if( istat == 0 ) then
          call write_error
          do ipr = 6,9,3
            write(ipr,112) File_in(ifich), File_name
          end do
          stop
        else
          call write_open_error(File_in(ifich),istat,1)
        endif
      endif
      Close(2)
    endif
  end do

  return
  110 format(///' Line not understood in the indata file :'//1x,A//, &
          ' If it is supposed to be a keyword check the spelling.',/ &
          ' If it is a line containing numbers, check:',/ &
          '   - How many numbers are supposed to be there.',/ &
          '   - Numbers must be separated by spaces.',/ &
          '   - Are there extra characters (points, tabulations...)'//)
  112 format(//' Error opening the file:',//3x,A,//  &
           5x,'It does not exist but the following file is found existing: ',//3x,A,// &
           5x,'Modify the input file and eventualy add other indata files with the convenient indexes !',//)
end

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

subroutine n_div_fpp_avanseuil(Eseuil,n_div_fpp,nfich)

  use declarations
  implicit none
  
  integer:: i, ifich, jfich, nfich

  integer, dimension(nfich):: n_div_fpp
  
  logical, dimension(nfich):: ok
  
  real(kind=db), dimension(nfich):: Eseuil
   
  ok(:) = .false.
  do ifich = 1,nfich
    if( ok(ifich) ) cycle
    i = 1
    do jfich = ifich+1,nfich
      if( ok(jfich) ) cycle
      if( abs( Eseuil(ifich) - Eseuil(jfich) ) > 1._db ) cycle
      i = i + 1
    end do 
    do jfich = ifich,nfich
      if( ok(jfich) ) cycle
      if( abs( Eseuil(ifich) - Eseuil(jfich) ) > 1._db ) cycle
      n_div_fpp(jfich) =  i
      ok(jfich) = .true.
    end do 
    
  end do
  return
end

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

! Evaluation of relative shifts between input data

  subroutine Shift_eval(delta_edge,Energphot,Eph1,Ephm,Epsii,Epsii_ref,Epsii_ref_man,Eseuil,Esmin, &
                        n_selec_core,n_shift_core,nfich,ninit1,ninit,ninitm,num_core,Shift_core,Shift_init,Shift_site)
  
  use declarations
  implicit none

  integer:: i, ifich, jfich, n, n_selec_core, n_shift_core, nfich, ninit1, ninitm
  
  integer, dimension(nfich):: ninit
  integer, dimension(10):: num_core
  
  logical:: Energphot, Epsii_ref_man
  logical, dimension(:), allocatable:: ok
  
  real(kind=db):: decal_min, delta_edge, Epsii_moy, Epsii_ref, Esmin
  real(kind=db), dimension(nfich):: decal_in, Eph1, Ephm, Eseuil, Shift_site 
  real(kind=db), dimension(n_shift_core):: Shift_core
  real(kind=db), dimension(ninitm,nfich):: Epsii, Shift_init 

! At this stage contains only the imposed value.  
  decal_in(:) = Shift_site(:)
  Shift_site(:) = 0._db
  
  Esmin = Eseuil(1)
  do ifich = 2,nfich
    Esmin = min( Esmin, Eseuil(ifich) ) 
  end do

  Energphot = ( Ephm(1) > Eseuil(1) ) .and. ( Eseuil(1) > Ephm(1) - Eph1(1) )

  if( .not. Energphot ) then
    do ifich = 1,nfich
      Shift_site(ifich) = Shift_site(ifich) + Eseuil(ifich) - Esmin
    end do
  endif

  allocate( ok(nfich) )
  ok(:) = .false.
 
  do ifich = 1,nfich
! Shift core is an option, not in the manual, to impose shift at specific sites  
    if( ifich <= n_shift_core ) then
      Shift_site(ifich) = Shift_site(ifich) + shift_core(ifich)
      cycle
    endif

    if( ok(ifich) ) cycle
    ok(ifich) = .true.

    Epsii_moy = sum( Epsii(1:ninit(ifich),ifich) )
    n = ninit(ifich) 

    do jfich = ifich+1,nfich
      if( abs( Eseuil(jfich) - Eseuil(ifich) ) > 1._db ) cycle
      n = n + ninit(jfich) 
      ok(jfich) = .true.
      Epsii_moy = Epsii_moy + sum( Epsii(1:ninit(jfich),jfich) )
    end do
    Epsii_moy = Epsii_moy / n

    if( Epsii_ref_man ) Epsii_moy = Epsii_ref

    decal_min = 0._db
    do jfich = ifich,nfich
      if( abs( Eseuil(jfich) - Eseuil(ifich) ) > 0.1_db ) cycle
      Shift_site(jfich) = Shift_site(jfich) - Epsii_moy + sum( Epsii(1:ninit(jfich),jfich) ) / ninit(jfich)
      decal_min = min( Shift_site(jfich), decal_min ) 
    end do
    if( .not. Epsii_ref_man ) then 
      do jfich = ifich,nfich
        if( abs( Eseuil(jfich) - Eseuil(ifich) ) > 1._db ) cycle
        Shift_site(jfich) = Shift_site(jfich) - decal_min
      end do
    endif
    
  end do
 
  deallocate( ok )
  
  Shift_site(:) = Shift_site(:) + decal_in(:)
  
  Shift_init(:,:) = 0._db
  do ifich = 1,nfich
    select case( ninit(ifich) )
      case(1)
        Epsii_moy = Epsii(1,ifich)
      case(2)
         if( ninit1 == 2 ) then
           Epsii_moy = sum( Epsii(1:ninit(ifich),ifich) ) / ninit(ifich)
         elseif( n_selec_core == 0 .or. n_selec_core == ninit(ifich) ) then 
           Epsii_moy = Epsii(2,ifich)
         elseif( num_core(1) == 2 ) then
           Epsii_moy = Epsii(2,ifich)
         else
           Epsii_moy = Epsii(1,ifich)
         endif
      case(4,6,10)
         if( ninit1 /= ninit(ifich) ) then
           if( n_selec_core == 0 ) then 
             Epsii_moy = sum( Epsii(ninit1+1:ninit(ifich),ifich) ) / ( ninit(ifich) - ninit1 )
           else
             n = 0
             do i = 1, n_selec_core
               if( num_core(i) > ninit1 ) n = n + 1
             end do
             if( n /= 0 ) then
               Epsii_moy = sum( Epsii(ninit1+1:ninit(ifich),ifich) ) / ( ninit(ifich) - ninit1 )
             else
               Epsii_moy = sum( Epsii(1:ninit1,ifich) ) / ninit1
             endif 
           endif
         else
           Epsii_moy = sum( Epsii(1:ninit(ifich),ifich) ) / ninit(ifich)
         endif
    end select
    Shift_init(1:ninit(ifich),ifich) = Epsii(1:ninit(ifich),ifich) - Epsii_moy + Shift_site(ifich)
    if( abs(Delta_edge) > eps10 .and. ( ninit(ifich) /= ninit1 ) ) Shift_init(1:ninit1,ifich) &
               = Shift_init(1:ninit1,ifich) + Delta_edge * Rydb
  end do

  return
end

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

! Realization of the output energy grid

subroutine Output_Energy_Grid( E_step_conv, Energphot, Eph1, Es_temp, Eseuil, Esmin, Estart, File_in, Length_line, &
                       Low_Energy_Step, n_selec_core, ne, ne_initl ,nes, nes_in, nfich, ninit, ninitm, nL_header, nsup, &
                       num_core, Shift_init, Step_conv )

  use declarations
  implicit none
  
  integer:: i, ie, ifich, ifichref, init, initlref, ipr, j, je, istat, Length_line, n, n_selec_core, ne_ref, &
            nemax, nes, nes_in, nfich, ninitm, nnombre 

  integer, dimension(10):: num_core
  integer, dimension(nfich):: ne, ninit, nL_header
  integer, dimension(ninitm,nfich):: ne_initl, nsup
  
  character(len=Length_name), dimension(nfich):: File_in
  
  logical:: Energphot, File_change, New_range, Step_conv
  logical, dimension(:,:), allocatable:: Fichdone 
  
  real(kind=db):: d, dmin, E, E_step, E_step_conv, Emax, Emin, Esmin, Estart, Estart_L, Low_Energy_Step
  real(kind=db), dimension(nes_in):: Es_temp 
  real(kind=db), dimension(nfich):: Eph1, Eseuil 
  real(kind=db), dimension(ninitm,nfich):: Shift_init 
  real(kind=db), dimension(:), allocatable:: Energ, Energ_ref 
  real(kind=db), dimension(:,:,:), allocatable:: Ef 

  ne_initl(:,:) = 0
  New_range = .true.
  nsup(:,:) = 0
  
  Estart_L = Estart
  do ifich = 1,nfich
    do init = 1,ninit(ifich)
    
      if( n_selec_core /= 0 ) then
        do i = 1,n_selec_core
          if( init == num_core(i) ) exit
        end do
        if( i > n_selec_core ) cycle
      endif
    
      if( Energphot ) then
        E = Eph1(ifich) - Eseuil(ifich) + Shift_init(init,ifich)
      else
        E = Eph1(ifich) + Shift_init(init,ifich) - Eseuil(ifich) + Esmin
      endif
      Estart_L = Min( Estart_L, E )
    end do
  end do

  do ifich = 1,nfich
    do init = 1,ninit(ifich)
    
      if( n_selec_core /= 0 ) then
        do i = 1,n_selec_core
          if( init == num_core(i) ) exit
        end do
        if( i > n_selec_core ) cycle
      endif
        
      if( Energphot ) then
        E = Eph1(ifich) - Eseuil(ifich) + Shift_init(init,ifich)
      else
        E = Eph1(ifich) + Shift_init(init,ifich) - Eseuil(ifich) + Esmin
      endif
      if( E > Estart_L - eps10 ) then
        nsup(init,ifich) = nint( ( E - Estart_L ) / Low_Energy_Step )
        ne_initl(init,ifich) = ne(ifich) + nsup(init,ifich)
      else
        ne_initl(init,ifich) = ne(ifich)
        nsup(init,ifich) = 0
      endif
    end do
  end do

  nemax = maxval( ne_initl )

  allocate( Ef(nemax,ninitm,nfich) )

  do ifich = 1,nfich

    allocate( Energ( ne(ifich) ) )
    
    open(2, file = File_in(ifich), status='old', iostat=istat)

    do i = 1,nL_header(ifich) + 1
      read(2,*)
    end do
    do
      n = nnombre(2,Length_line)
      if( n == 0 ) then
        read(2,*)
        exit
      endif
      read(2,*)
    end do

    do ie = 1,ne(ifich)
      Read(2,*) Energ(ie)
    end do

    Close(2)
    
    Energ(:) = Energ(:) / Rydb

    if( ifich == 1 ) then
      if( ne(ifich) > 1 ) then
        allocate( Energ_ref( ne(1) ) )
        Energ_ref(:) = Energ(:)
        ne_ref = ne(1)
          
        E_step = Energ(2) - Energ(1)
      else
        E_step = 0._db
      endif 
    endif
    
    do init = 1,ninit(ifich)

      if( ne_initl(init,ifich) == 0 ) cycle
      
      do ie = 1,ne(ifich)
        je = ie + nsup(init,ifich)
        Ef(je,init,ifich) = Energ(ie)
      end do

      do ie = nsup(init,ifich),1,-1
        Ef(ie,init,ifich) = Ef(ie+1,init,ifich) - Low_Energy_Step
      end do

      n = ne_initl(init,ifich)
      Ef(1:n,init,ifich) = Ef(1:n,init,ifich) + Shift_init(init,ifich)
    end do

    deallocate( Energ )
    
  end do

 ! When there are absorbing atom with different Z, 
 ! the maximum energy of the convoluted spectra corresponds to the one of maximum edge energy  
  E = Eseuil(1)
  j = 0
  do ifich = 2,nfich
    if( abs( Eseuil(ifich) - E ) < 1._db ) cycle
    if( ifich == 2 ) j = 1
    if( Eseuil(ifich) < E ) cycle
    E = Eseuil(ifich)
    j = ifich
  end do

  if( j == 0 ) then
    do init = 1,ninit(1)
      if( ne_initl(init,1) /= 0 ) exit
    end do
    Emax = Ef(ne_initl(init,1),init,1)
  else
    do init = 1,ninit(j)
      if( ne_initl(init,j) /= 0 ) exit
    end do
    Emax = Ef(ne_initl(init,j),init,j)
  endif

  Emin = Estart_L
  
  do ifich = 1,nfich
    do init = 1,ninit(ifich)
    
      if( ne_initl(init,ifich) == 0 ) cycle
         
      if( j == 0 .or. abs( Eseuil(ifich) - Eseuil(max(j,1)) ) < 1._db ) then
        if( Ef(ne_initl(init,ifich),init,ifich) - Emax > Emax - Emin ) then
          Emax = Ef(ne_initl(init,ifich),init,ifich)
        elseif( .not. Emax - Ef(ne_initl(init,ifich),init,ifich) > Ef(ne_initl(init,ifich),init,ifich) - Emin ) then  
          Emax = min( Ef(ne_initl(init,ifich),init,ifich), Emax )
        endif
      endif
    end do
  end do

! Grid elaboration

  if( Step_conv .and. E_step > eps10 ) then
  
    if( E_step_conv > eps10 ) E_step = E_step_conv

    nes = nint( ( Emax - Emin ) / E_step ) + 1
    
    do ie = 1,nes
      Es_temp(i) = Emin + ( ie - 1 ) * E_step
    end do
    
  elseif( New_range ) then

    E_step = Low_Energy_Step
        
    n = nint( ( Energ_ref(1) - Emin ) / E_step )
    ie = 0

    if( n > 0 ) then
      E = Energ_ref(1) - n * E_step
      if( E > Emin - eps10 ) then
        do je = -n, -1
          ie = ie + 1
          Es_temp(ie) = Energ_ref(1) + je * E_step
        end do  
      else
        ie = ie + 1
        Es_temp(ie) = Emin
        do je = - n + 1, 1
          ie = ie + 1
          Es_temp(ie) = Energ_ref(1) + je * E_step
        end do
      endif
    endif

    ie = ie + 1
    Es_temp(ie) = Energ_ref(1)
      
    if( ne_ref > 1 ) then
      
      do    
 
        E_step = 1000000._db

        do ifich = 1,nfich
          do init = 1,ninit(ifich)
            do je = 1,ne_initl(init,ifich) - 1
              if( Ef(je,init,ifich) > Es_temp(ie) - eps10 ) exit
            end do
            if( je == ne_initl(init,ifich) ) cycle           
            if( Ef(je+1,init,ifich) - Ef(je,init,ifich) > E_step + eps10 ) cycle
            E_step = Ef(je+1,init,ifich) - Ef(je,init,ifich) 
          end do
        end do
        
        ie = ie + 1
        Es_temp(ie) = Es_temp(ie-1) + E_step
        if( Es_temp(ie) > Emax - eps10 ) exit
        
      end do

      if( Es_temp(ie) > Emax + eps10 ) then
        nes = ie - 1
      else
        nes = ie
      endif
           
    endif
      
  else

    allocate( Fichdone(ninitm,nfich) )
    Fichdone(:,:) = .false.

    dmin = 1000000000._db
    do ifich = 1,nfich
      do init = 1,ninit(ifich)
        if( ne_initl(init,1) == 0 ) cycle
        if( Shift_init(init,ifich) > dmin - eps10 ) cycle
        ifichref = ifich
        initlref = init
        dmin = Shift_init(initlref,ifichref) 
      end do
    end do
    Fichdone(initlref,ifichref) = .true.

    je = 0
    ie = 0
    do
      je = je + 1
      if( je <= ne_initl(initlref,ifichref) ) then 
        E = Ef(je,initlref,ifichref)
        ie = ie + 1

        if( ie > nes_in ) then
          call write_error
          do ipr = 6,9,3
            write(ipr,110)
          end do
          stop
        endif

        Es_temp(ie) = E
        if( E > Emax - eps10 ) exit
      endif

      if( ie == 1 ) cycle 

      d = Shift_init(initlref,ifichref)
      dmin = 1000000000._db
      File_change = .false.
      do ifich = 1,nfich
        do init = 1,ninit(ifich)
          if( Fichdone(init,ifich) ) cycle
          if( ne_initl(init,ifich) == 0 ) cycle
          if( Shift_init(init,ifich) < d + eps10 ) then
            Fichdone(init,ifich) = .true.
            cycle
          endif  
          if( Shift_init(init,ifich) > d + eps10 .and. Shift_init(init,ifich) < dmin &
             .and. ( E > Ef(nsup(init,ifich)+1,init,ifich) - eps10 .or. je > ne_initl(initlref,ifichref) ) ) then
            File_change = .true.
            initlref = init
            ifichref = ifich
            dmin = Shift_init(init,ifich)
          endif
        end do
      end do
      if( File_change ) then
        Fichdone(initlref,ifichref) = .true.
        do je = 1,ne_initl(initlref,ifichref)
          if( Ef(je,initlref,ifichref) > Es_temp(ie-1) + eps10 ) exit
        end do
        if( je > ne_initl(initlref,ifichref) ) exit
        Es_temp(ie) = Ef(je,initlref,ifichref)  
      endif

    end do
    
    if( E > Emax + eps10 ) then
      nes = ie - 1
    else
      nes = ie
    endif

    deallocate( Fichdone )

  endif
  
  deallocate( Ef, Energ_ref )

  if( nes == 0 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,120)
    end do
    stop
  endif

  return
  110 format(///' Problem in the energy grid elaboration in convolution',/ &
                ' Too much points !' //)
  120 format(///' Taking into account the energy shifts, there is no overlap between',/ &
                ' the energy ranges of the different absorbing atoms !',/ &
                ' No summation and convolution are possible !' //)
end

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

! Name of the output file

  subroutine Conv_out_name(Abs_U_iso_inp,bav_open,check_conv,Chem,Chemin,convolution_out,Convolution_out_all,Dafs_bio,Delta_gauss, &
                        Deuxieme,E_cut,E_cut_man,Epsii_ref,Epsii_ref_man,File_in,File_scan_in,File_scan_out,File_scan_out_all, &
                        Gamma_hole,Gamma_hole_man,Gamma_max,Gamma_max_man,Length_line,nfich,nomfichbav,Scan_a,Scan_true, &
                        U_iso_man,XES)

  use declarations
  implicit none

  integer:: i, ifich, ipr, istat, jfich, k, Length, l_max, L, Length_line, long, longf, n, nfich, nnombre, ns
  
  character(len=8):: dat
  character(len=10):: mot10, tim
  character(len=50):: com_date, com_time
  character(len=Length_name):: chemin, convolution_out, File_scan_out, File_name, nomfichbav, Word
  character(len=Length_name), dimension(nfich):: convolution_out_all, File_in, File_scan_in, File_scan_out_all

  logical:: bav_open, check_conv, Chem, Dafs_bio, Deuxieme, E_cut_man, Epsii_ref_man, Gamma_hole_man, Gamma_max_man, Scan_a, &
            Scan_true, U_iso_man, XES
  
  real(kind=db):: Abs_U_iso_inp, Delta_gauss, E_cut, Epsii_ref, Gamma_hole, Gamma_max, x, y
  
  if( convolution_out == ' ' ) then
 
    do ifich = 0,nfich
      if( nfich == 1 .and. ifich > 0 ) exit

      if( ifich == 0 ) then
        l_max = 10000
        do jfich = 1,nfich
          long = len_trim( File_in(jfich) )
          if( long < l_max ) then
            l_max = long
            File_name = File_in(jfich)
          endif
        end do
      else 
        File_name = File_in(ifich)
      endif
      Length = len_trim( File_name )
      if( Length > 4 ) then
        if( File_name(Length-3:Length) == '.txt' )  File_name(Length-3:Length) = '    '     
      endif
      
      Length = len_trim( File_name )
      if( nfich > 1 .and. ifich == 0 .and. Length > 2 ) then
        if( File_name(Length-1:Length-1) == '_' ) then
          File_name(Length-1:Length) = '  '
        elseif( File_name(Length-2:Length-2) == '_' ) then
          File_name(Length-2:Length) = '   '
        endif
      endif

      do k = 1,6
      
        select case(k)
          case(1)
           if( .not. E_cut_man ) cycle
           y = E_cut
          case(2)
           if( .not. Gamma_max_man ) cycle
           y = Gamma_max * Rydb
          case(3)
           if( .not. Gamma_hole_man ) cycle
           y = Gamma_hole * Rydb
          case(4)
           if( .not. Epsii_ref_man ) cycle
           y = Epsii_ref * Rydb          
          case(5)
           if( abs( Delta_gauss ) < eps10 ) cycle
           y = Delta_gauss          
          case(6)
           if( .not. U_iso_man ) cycle
           y = Abs_U_iso_inp          
        end select
        
        x = y * 0.1_db
        do i = 0,4
          x = 10._db * x
          if( abs( x - nint(x) ) < eps10 ) exit
        end do

        select case(i)
          case(0)
            write(mot10,'(i10)') nint( y )
          case(1)
            write(mot10,'(f10.1)') y
          case(2)
            write(mot10,'(f10.2)') y
          case(3)
            write(mot10,'(f10.3)') y
          case default
            write(mot10,'(f10.4)') y
        end select
        
        mot10 = adjustl( mot10 )   
        L = len_trim( mot10 )
        if( i > 0 ) mot10(L-i:L-i) = 'p'
     
        Length = len_trim( File_name ) + 1

        select case(k)
          case(1)          
            File_name(Length:Length+2) = '_EF'
          case(2)
            File_name(Length:Length+2) = '_Gm'
          case(3)
            File_name(Length:Length+2) = '_GH'
          case(4)
            File_name(Length:Length+2) = '_Ep'
          case(5)
            File_name(Length:Length+2) = '_Ga'
          case(6)
            File_name(Length:Length+2) = '_Ui'
        end select
        Length = Length + 2
        
        File_name(Length+1:Length+L) = mot10(1:L)
        
      end do
     
      Length = len_trim( File_name ) + 1
      if( XES ) then
        File_name(Length:Length+12) = '_XES_conv.txt'
      else
        File_name(Length:Length+8) = '_conv.txt'
      endif
      
      if( ifich == 0 ) then
        convolution_out = File_name
      else
        convolution_out_all(ifich) = File_name
      endif

    end do      
  
  else
  
    if( Chem ) then
      long = len_trim(chemin)
      longf = len_trim(convolution_out)
      convolution_out = chemin(1:long) // convolution_out(1:longf)
    endif

    Length = len_trim( convolution_out )

    if( Length > 5 ) then
      if( convolution_out(Length-3:Length) /= '.txt' ) convolution_out(Length+1:Length+4) = '.txt'  
    endif

    do ifich = 1,nfich
      if( nfich == 1 ) exit
      File_name = convolution_out
      Length = len_trim( convolution_out )
      if( convolution_out(Length-8:Length) == '_conv.txt' ) then
        File_name(Length-7:Length) = '        '
        call ad_number(ifich,File_name,Length_name)
        Length = len_trim( File_name )
        if( XES ) then
          File_name(Length+1:Length+16) = '_xes_conv.txt'
        else
          File_name(Length+1:Length+9) = '_conv.txt'
        endif 
      else
        if( XES ) then
          File_name(Length-3:Length+1) = '_xes_'
        else
          File_name(Length-3:Length) = '_   '
        endif
        call ad_number(ifich,File_name,Length_name)
        Length = len_trim( File_name )
        File_name(Length+1:Length+4) = '.txt'
      endif 
      convolution_out_all(ifich) = File_name
    end do

    do ifich = 1,nfich
      if( Convolution_out /= File_in(ifich) ) cycle
      call write_error
      do ipr = 6,9,3
        write(ipr,120) File_in(ifich), File_in(ifich)
      end do
      stop
    end do
    
  endif

  if( .not. bav_open .and. Check_conv ) then
    Length = len_trim(convolution_out)
    nomfichbav = ' '
    nomfichbav(1:Length-4) = convolution_out(1:Length-4)
    nomfichbav(Length-3:Length+4) = '_bav.txt'
    open(3, file = nomfichbav, status='unknown',iostat=istat)
    if( istat /= 0 ) call write_open_error(nomfichbav,istat,1)
    bav_open = .true.
    call date_and_time( date = dat, time = tim )
    com_date = '   Date = ' // dat(7:8) // ' ' // dat(5:6) // ' ' // dat(1:4)
    com_time = '   Time = ' // tim(1:2)  // ' h ' // tim(3:4) // ' mn ' // tim(5:6) // ' s'
    write(3,'(1x,A/A/A)') Revision, com_date, com_time
    
    write(3,'(/A)') ' Convolution of the file(s):' 
    do ifich = 1,nfich
      write(3,'(2x,A)') File_in(ifich)
    end do

  endif

  if( Check_conv .and. .not. Deuxieme ) write(3,110)

  if( Scan_a .and. .not. Scan_true ) then
    do ifich = 1,nfich
      Word = ' '
      File_name = File_in(ifich)
      Length = len_trim(File_name) - 4
      if( nfich == 1 ) then
        ns = 0
      else
        if( File_name(Length-2:Length-2) == '_' ) then
          ns = 3
        elseif( File_name(Length-1:Length-1) == '_' ) then
          ns = 2
        else
          ns = 0
        endif
      endif
      Word(1:Length-ns) = File_name(1:Length-ns)
      Word(Length-ns+1:Length-ns+5) = '_scan'
      if( ns > 0 ) Word(Length-ns+6:Length+5) = File_name(Length-ns+1:Length)
      Word(Length+6:Length+9) = '.txt'
      File_scan_in(ifich) = Word
    end do
    Scan_true = .true.
  endif

  Dafs_bio = .false.
  if( Scan_true ) then
    open(2, file = File_scan_in(1), status='old',iostat=istat)
    n = nnombre(2,Length_line)
    read(2,*) n
    if( n == 4 ) Dafs_bio = .true.
    close(2)
  endif

  if( Scan_true .and. .not. Dafs_bio ) then
    if( File_scan_out /= ' ' ) then
      Length = len_trim(File_scan_out)
      if( File_scan_out(Length-3:Length-3) /= '.' ) File_scan_out(Length+1:Length+4) = '.txt' 
    else
      File_name = convolution_out
      Length = len_trim(File_name)
      File_name(Length-7:Length+5) = 'scan_conv.txt'
      File_scan_out = File_name
      do ifich = 1,nfich
        if( nfich == 1 ) exit
        File_name = convolution_out_all(ifich)
        Length = len_trim(File_name)
        File_name(Length-7:Length+5) = 'scan_conv.txt'
        File_scan_out_all(ifich) = File_name
      end do
    endif
  endif

  return
  110 format(/' ---- Convolution ------',100('-'))
  112 format(//' Error opening the file:',//3x,A,//  &
           5x,'It does not exist but the following file is found existing: ',//3x,A,// &
           5x,'Modify the input file and eventualy add other indata files with the convenient indexes !',//)
  120 format(///' The output file has the same name than one of the input files:',/ &
                3x,A,// &
                ' Choose another one with keyword "Conv_out" !',//)
end

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

! Reading of dimensions

subroutine Dimension_file(Abs_U_iso,Bulk_site,Cor_abs,Eintmax,En_Fermi,Eph1,Ephm, Epsii,f0_forward,File_in, &
          File_scan_in,fprim,Full_self_abs,Green_int,jseuil,Length_line,Magn,n_bulk_site,n_col_max,n_mat_pol, &
          n_z_SRXRD,n_z_SRXRD_abs,n_z_SRXRD_abs_bulk,n_z_SRXRD_bulk,natomsym_f,ne,nfich,ninit1,ninit,ninitm,nL_header, &
          Non_resonant_bulk,nphim,npldafs,npldafs_NR,nseuil,nxan,nxan_lib,Sample_thickness,Scan_true,Self_abs,Signal_Sph,SRXRD, &
          Surface_ref,Tenseur,Tenseur_car,V0muf,Volume_maille,Volume_maille_bulk)

  use declarations
  implicit none

  integer:: eof, i, ie, ifich, ipl, ipr, istat, jseuil, Length, Length_line, n, n_bulk_site, n_col_max, n_mat_pol, n_mat_pol1, &
    n_z_SRXRD, n_z_SRXRD_abs, n_z_SRXRD_abs_bulk, n_z_SRXRD_bulk, nfich, ninit1, ninitm, nnombre, nphi, nphim, nphim1, &
    npldafs, npldafs_NR, npldafs1, nseuil, Numat, nxan, nxan1
  integer, dimension(nfich):: ne, ninit, nL_header

  character(len=15):: nomab
  character(len=132):: mot
  character(len=Length_name), dimension(nfich):: File_in, File_scan_in
  character(len=15), dimension(:), allocatable:: nomcol

  logical:: Cor_abs, fprim, Full_self_abs, Green_int, Magn, Non_resonant_bulk, nxan_lib, Scan_true, Self_abs, Signal_Sph, SRXRD, &
            Tenseur, Tenseur_car
            
  logical, dimension(nfich):: Bulk_site

  real(kind=db):: Eintmax, Eph, Eseuil, f0_forward, fpp_avantseuil, Sample_thickness, Surface_ref, Volume, &
                  Volume_maille, Volume_maille_bulk, x  
  real(kind=db), dimension(ninitm,nfich):: Epsii
  real(kind=db), dimension(nfich):: Abs_U_iso, En_Fermi, Eph1, Ephm, natomsym_f, V0muf

  SRXRD = .false.
  n_z_SRXRD = 0
  n_z_SRXRD_bulk = 0
  n_bulk_site = 0
  natomsym_f(:) = 0._db
  Non_resonant_bulk = .false.
  npldafs_NR = 0
 
  do ifich = 1,nfich

    n_mat_pol = 0

    open(2, file = File_in(ifich), status='old', iostat=istat)

    do i = 1,nL_header(ifich)
      read(2,*)
    end do
    n = nnombre(2,Length_line)
    if( n == 11 + ninit(ifich) ) then
      read(2,*) Eseuil, Numat, nseuil, jseuil, fpp_avantseuil, V0muf(ifich), En_fermi(ifich), ninit(ifich), &
        ninit1, Epsii(1:ninit(ifich),ifich), Volume, f0_forward
    elseif( n == 12 + ninit(ifich) ) then
      read(2,*) Eseuil, Numat, nseuil, jseuil, fpp_avantseuil, V0muf(ifich), En_fermi(ifich), ninit(ifich), &
        ninit1, Epsii(1:ninit(ifich),ifich), Volume, f0_forward, natomsym_f(ifich)
    elseif( n == 13 + ninit(ifich) ) then
      read(2,*) Eseuil, Numat, nseuil, jseuil, fpp_avantseuil, V0muf(ifich), En_fermi(ifich), ninit(ifich), &
        ninit1, Epsii(1:ninit(ifich),ifich), Volume, Surface_ref, f0_forward, natomsym_f(ifich)
    else
      read(2,*) Eseuil, Numat, nseuil, jseuil, fpp_avantseuil, V0muf(ifich), En_fermi(ifich), ninit(ifich), &
        ninit1, Epsii(1:ninit(ifich),ifich), Volume, Surface_ref, f0_forward, natomsym_f(ifich), Abs_U_iso(ifich)
    endif

    if( ninit1 < 0 ) then
      Green_int = .true.
      ninit1 = abs( ninit1 )
    else
      Green_int = .false.
    endif
    
    n = nnombre(2,Length_line)
   
    if( n >= n_col_max ) then
      call write_error
      do ipr = 6,9,3
        write(ipr,110) n_mat_pol
      end do
      stop
    endif
    
    npldafs = n / 2
    
    if( n /= 0 ) then    
      read(2,*)
      read(2,*)
    endif

    n = nnombre(2,Length_line)
    if( n == 5 ) then

      SRXRD = .true.
      read(2,*) x, x, x, n, i
      if( i == 1 ) then 
        Bulk_site(ifich) = .true.
        n_bulk_site = n_bulk_site + 1
        n_z_SRXRD_bulk = max( n, n_z_SRXRD_bulk ) 
      else
        if( i == -1 ) then
          Non_resonant_bulk = .true.
          npldafs_NR = npldafs
        endif
        n_z_SRXRD = max( n, n_z_SRXRD )
      endif
      do i = 1, 1 + 2 * n + 1
        read(2,*) 
      end do
      if( Non_resonant_bulk ) read(2,*) 
      if( Bulk_site(ifich) ) read(2,*)
      
    elseif( n /= 0 ) then
      read(2,*)
      n = nnombre(2,Length_line)
      read(2,'(A)') mot
      do i = 1,131
        if( mot(i:i+1) == 'mu' ) then
          Full_self_abs = .true.
          exit
        endif
      end do
      Self_abs = i > 131
      backspace(2)
    endif            

    Cor_abs = Self_abs .or. Full_self_abs

    if( Bulk_site(ifich) ) then
      Volume_maille_bulk = Volume
    else
      Volume_maille = Volume
    endif

    read(2,'(10x,a15)') nomab
    nomab = adjustl( nomab )

    n = ( nnombre(2,Length_line) - 1 ) / ninit(ifich)
    if( npldafs > 0 ) then
      if( Full_self_abs ) then
        nxan = n - 6 * npldafs
      elseif( Self_abs ) then
        nxan = n - 4 * npldafs
      else
        nxan = n - 2 * npldafs
      endif
      if( nomab(1:8) == 'Sum_dd_r' .or. nomab(1:8) == 'Sum_tot_' ) Signal_Sph = .true.
    elseif( fprim .and. nomab(1:1) == 'D' ) then
      Tenseur = .true.
      nxan = 0
      if( nomab(1:4) == 'D_xx' ) Tenseur_car = .true.
      Magn = nomab(1:6) == 'D_xx_r'
      if( Magn ) then
        npldafs = n / 2
      else
        npldafs = n
      endif
    else
      nxan = n
    endif

! number of polarization matrix
    backspace(2)
    allocate( nomcol(nxan) )
    read(2,'(10x,10000a15)') nomcol(:)
    do ipl = 1,nxan
      nomab = nomcol(ipl)
      Length = len_trim( nomab )
      if( Length <= 3 ) cycle
      do i = 1,Length-2
        if( nomab(i:i+2) /= '_ss' ) cycle
        n_mat_pol = n_mat_pol + 1
        exit
      end do
    end do
    deallocate( nomcol )
    
    if( n_mat_pol /= 0 .and. Sample_thickness < eps10 ) then
      call write_error
      do ipr = 6,9,3
        if( n_mat_pol == 1 ) then
          write(ipr,120) n_mat_pol
        else
          write(ipr,122) n_mat_pol
        endif 
      end do
      stop
    endif

    nxan = nxan - 6 * n_mat_pol
     
    if( ifich == 1 ) then
      n_mat_pol1 = n_mat_pol
      nxan1 = nxan
      npldafs1 = npldafs
    else
      if( nxan_lib ) nxan = min(nxan1,nxan)
      if( nxan1 /= nxan .or. npldafs1 /= npldafs ) then
        call write_error
        do ipr = 6,9,3
          write(ipr,130) nxan1, npldafs1, n_mat_pol1, ifich, nxan, npldafs, n_mat_pol
        end do
        stop
      endif
    endif

    do ie = 1,10000
      Read(2,*,iostat=eof) Eph
      if( eof /= 0 ) exit
      if( ie == 1 ) eph1(ifich) = eph
      Ephm(ifich) = eph
      if( eph > Eintmax ) exit
    end do
    ne(ifich) = ie - 1

    if( ne(ifich) < 2 ) then
      call write_error
      do ipr = 6,9,3
        write(ipr,140) File_in(ifich)
      end do
      stop
    endif

    Close(2)
  end do

  n_z_SRXRD_abs_bulk = n_bulk_site
  if( SRXRD ) then
    n_z_SRXRD_abs = nfich - n_bulk_site
  else
    n_z_SRXRD_abs = 0
  endif 

  if( npldafs > 0 ) then 
    nphim = 1
  else
    nphim = 0
  endif
  if( Scan_true ) then
    do ifich = 1,nfich
      open(2, file = File_scan_in(ifich), status='old',iostat=istat)
      if( istat /= 0 ) call write_open_error(File_scan_in(ifich),istat,1)
      n = nnombre(2,Length_line)
      do ipl = 1,npldafs
        read(2,*) nphi
        nphim = max( nphim, nphi )
      end do
      Close(2)
      
      if( ifich == 1 ) then
        nphim1 = nphim
      elseif( nphim /= nphim1 ) then
        call write_error
        do ipr = 6,9,3
          write(ipr,150) nphim1, nphim
        end do
        stop
      endif
 
    end do
  endif

  return
  110 format(///' The maximum of number of columns =',i6,', is reached in the files to convolute.', & 
               ' Use the keyword Length_line to increase this number.', &
               ' Some compilers such as ftn95 does not allow higher number.',// ) 
  120 format(///'     In convolution, there is', i2,' polarization matrix to calculate',/&
                '     A sample tickness must be given in micrometer with keyword "Sample_tickness" !',// )
  122 format(///'     In convolution, there are', i3,' polarization matrix to calculate',/&
                '     A sample tickness must be given in micrometer with keyword "Sample_tickness" !',// )
  130 format(///'     In convolution, indata files are different,',/ &
                '  file   1 : number of xanes =',i4,', number of reflexions =',i4,', number of polarization matrix =',i4,/ &
                '  file',i4,' : number of xanes =',i4,', number of reflexions =',i4,', number of polarization matrix =',i4)
  140 format(//' For the convolution, the number of energy must be greater than one !'/, &
               ' It is not the case in the file:'//A//)
  150 format(///'     Input files are different,',/ &
                '  file 1 : nphim =',i4,'  file ',i2,' nphim =',i4)
end

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

subroutine Col_name(Ampl_phase,Analyzer,Bormann,Cor_abs,Dafs_bio,Double_cor,File_in,File_scan_in,fprim,Full_self_abs, &
      hkl_dafs,Length_line,n_bir,n_col,n_mat_pol,n_stokes,nL_header,nom_col,npldafs,npldafs_b,nxan,Self_abs,Signal_sph, &
      Stokes,Stokes_name,Stokes_param,Sup_sufix,Tenseur,XES)

  use declarations
  implicit none

  integer:: i, i_bir, i_hk, ii, index_in, index_out, ipl, istat, j, k, kk, Length, Length_line, m, n, n_bir, n_col, n_col_o, &
    n_col_in, n_mat_pol, n_stokes, nc, nL_header, nnombre, npldafs, npldafs_b, nxan
  integer, dimension(npldafs):: nphi

  character(len=Length_name):: File_in, File_scan_in
  character(len=Length_word):: nomab, nomac
  character(len=length_line):: motl
  character(len=Length_word), dimension(n_stokes):: Stokes_name
  character(len=Length_word), dimension(n_col):: nom_col
  character(len=Length_word), dimension(2*n_col):: nom_col_o
  character(len=Length_word), dimension(:), allocatable:: nom_col_in

  integer, dimension(3,npldafs_b):: hkl_dafs
  integer, dimension(npldafs):: n_index_hk

  logical:: Ampl_phase, Analyzer, Bormann, Cor_abs, Dafs_bio, Double_cor, fprim, Full_self_abs, Self_abs, Signal_sph, &
            Stokes, Sup_sufix, Tenseur, XES

  real(kind=db), dimension(5,n_stokes):: Stokes_param

  open(2, file = File_in, status='old', iostat=istat)
  if( istat /= 0 ) call write_open_error(File_in,istat,1)

  do i = 1,nL_header + 1
    read(2,*)
  end do
  do
    n = nnombre(2,Length_line)
    if( n == 0 ) exit
    read(2,*)
  end do

  if( Dafs_bio ) then
    n_col_in = 1 + nxan
  else
    n_col_in = 1 + nxan + 6 * n_mat_pol + 2 * npldafs
    if( Self_abs ) n_col_in = n_col_in + 2*npldafs
    if( Full_self_abs ) n_col_in = n_col_in + 4*npldafs
  endif

  allocate( nom_col_in(n_col_in) )

  read(2,'(A)') motl

  call Extract_word(Length_line,motl,nom_col_in,n_col_in,n_col_in)
  if( Sup_sufix .or. nxan /= 0 .or. n_mat_pol /= 0 ) then
    if( nxan == 0 .and. n_mat_pol == 0 ) then
      nc = n_col_in - 1
    else
      nc = nxan + 6 * n_mat_pol 
    endif 

    index_out = 0
    index_in = 1   ! Energy
    do ipl = 1,nc
      index_out = index_out + 1
      index_in = index_in + 1
      nomab = nom_col_in(index_in)
      Length = len_trim( nomab )
      if( nomab(Length-1:Length-1) == '_' ) nomab(Length-1:Length) = '  '
      if( nomab(Length-2:Length-1) == '_L' ) nomab(Length-2:Length) = '   '
      if( nomab(Length-2:Length-1) == '_M' ) nomab(Length-2:Length) = '   '
      if( nomab(Length-2:Length-1) == '_N' ) nomab(Length-2:Length) = '   '
      if( nomab(Length-2:Length-1) == '_O' ) nomab(Length-2:Length) = '   '

      if( index_in > nxan + 1 .and. index_in <= nxan + 6 * n_mat_pol + 1 ) then

        if( mod(index_in-nxan-1,6) /= 1 ) then
          index_out = index_out - 1
          cycle
        endif
        Length = len_trim( nomab )
        do k = 1,Length
          if( nomab(k:k) == '_' ) exit
        end do
        do m = k+1,Length
          nomab(m:m) = ' '
        end do
        nomac = nomab
        nomac(k:k) = ' '
        Length = len_trim( nomab )
        nomab(Length+1:Length+2) = 'ss'
        nom_col_o(index_out) = nomab
        nomab(Length+1:Length+2) = 'pp'
        index_out = index_out + 1
        nom_col_o(index_out) = nomab
        do n = 1,n_stokes
          if( Stokes_name(n) /= 'no_name' ) then
            nomab = Stokes_name(n)
          else
           nomab = nomac
           call ad_number(n,nomab,Length_word)
          endif
          index_out = index_out + 1
          nom_col_o(index_out) = nomab
          Length = min( Length_word-3, len_trim( nomab ) )
          do i_bir = 2,n_bir
            if( i_bir == 2 ) then
              nomab(Length+1:Length+3) = 'nB'
            elseif( i_bir == 3 ) then
              nomab(Length+1:Length+3) = 'nD'
            else
              nomab(Length+1:Length+3) = 'nO'
            endif
            index_out = index_out + 1
            nom_col_o(index_out) = nomab
          end do
        end do
      else 
        if( Length > 4 ) then
          if( nomab(Length-3:Length-2) == '_r' .or. nomab(Length-3:Length-2) == '_i' ) nomab(Length-3:Length-2) = '  '
        endif
        if( nxan == 0 ) then
          nom_col_in(index_in) = nomab
        else
          nom_col_o(index_out) = nomab
        endif
      endif
      nomab = adjustl( nom_col_o(index_out) )
      if( XES .and. nomab(1:7) == '<xanes>' ) then
        nomab(1:7) = '  <xes>'
        call center_word(nomab,Length_word)
        nom_col_o(index_out) = nomab
      endif
       
    end do
  endif

  i = nxan + n_stokes * n_mat_pol
  j = nxan + 6 * n_mat_pol + 1

  if( Dafs_bio ) then

! When using Dafs_bio, there is necessarily a scan file with 4 polarizations value
! We use this file to get the hkl_dafs indexes
    open(2, file = File_scan_in, status='old',iostat=istat)

    do ipl = 1,npldafs
      read(2,*) nphi(ipl)
    end do
        
    read(2,*)
    do ipl = 1,npldafs
      read(2,*) hkl_dafs(:,ipl)
      do k = 1,nphi(ipl)
        read(2,*) 
      end do
    end do
 
    Close(2)
 
    if( n_col == nxan + 6 * n_mat_pol + npldafs ) then
      n = 1
    else
      n = 5
    endif

    do ipl = 1,npldafs
      do j = 1,n
        i = i + 1
        nomab = '('
        do k = 1,3
          if( j == 4 .or. j == 5 ) then
            call ad_number(-hkl_dafs(k,ipl),nomab,Length_word)
          else
            call ad_number(hkl_dafs(k,ipl),nomab,Length_word)
          endif
          Length = len_trim(nomab) + 1
          if( Length > Length_word ) exit
          if( k == 3 ) then
            nomab(Length:Length) = ')'
          else
            nomab(Length:Length) = ','
          endif
        end do
        if( Length+1 <= Length_word ) then 
          if( j == 2 .or. j == 4 ) then
            nomab(Length+1:Length+1) = 's'
          elseif( j == 3 .or. j == 5 ) then
            nomab(Length+1:Length+1) = 'p'
          endif
        endif
        call center_word(nomab,Length_word)
        nom_col_o(i) = nomab
      end do
    end do

  else

    do ipl = 1,npldafs

      if( Tenseur ) then

        i = i + 1
        j = j + 1
        nomab = nom_col_in(j)
        Length = len_trim( nomab )
        nomab(Length+1:Length+1) = 'p'
        nom_col_o(i) = nomab

        i = i + 1
        nomab(Length+1:Length+1) = 's'
        nom_col_o(i) = nomab

        cycle

      elseif( Bormann ) then

        i = i + 1
        j = j + 2
        nomab = nom_col_in(j)
        nomab(1:1) = 'e'
        nom_col_o(i) = nomab

        i = i + 1
        Length = min( len_trim( nomab ) + 1, Length_word )
        nomab(Length:Length) = 'i'
        nom_col_o(i) = nomab

        cycle

      endif

! step 2, for nom_col_in, because real part and imaginary part
      j = j + 2
      nomab = nom_col_in(j)
      if( Signal_sph ) then
        Length = len_trim( nomab )
        nomab(2:Length+1) = nomab(1:Length)
      endif
      nomab(1:1) = 'I'
      i = i + 1
      nom_col_o(i) = nomab

      if( Ampl_phase ) then
        i = i + 1
        nomab = nom_col_in(j)
        Length = len_trim( nomab )
        do k = 1,Length-1
          nomab(k:k) = nomab(k+1:k+1)
        end do
        nomab(Length:Length) = 'A'
        nom_col_o(i) = nomab

        i = i + 1
        nomab(Length:Length) = 'P'
        nom_col_o(i) = nomab
      endif

      if( fprim ) then
        i = i + 1
        nomab = nom_col_in(j)
        Length = len_trim( nomab )
        do k = 1,Length-1
          nomab(k:k) = nomab(k+1:k+1)
        end do
        nomab(Length:Length) = 'p'
        nom_col_o(i) = nomab

        i = i + 1
        nomab(Length:Length) = 's'
        nom_col_o(i) = nomab
      endif

      if( Cor_abs ) then
        i = i + 1
        nomab = nom_col_in(j)
        Length = len_trim( nomab )
        do k = min(Length+1,Length_word),2,-1
          nomab(k:k) = nomab(k-1:k-1)
        end do
        nomab(1:2) = 'Ic'
        nom_col_o(i) = nomab
      endif

      if( Double_cor ) then
        i = i + 1
        nomab = nom_col_o(i-1)
        nomab(2:2) = 'd'
        nom_col_o(i) = nomab
      endif

      if( Cor_abs ) then
        do ii = 1,4     ! the mu
          if( .not. Full_self_abs .and. ii > 2 ) exit
          i = i + 1
          j = j + 1
          nom_col_o(i) = nom_col_in(j)
        end do
      endif

      if( Stokes .and. mod(ipl,4) == 0 ) then
        do ii = 1,n_stokes

          i = i + 1

          nomab = ' '

          if( Stokes_name(ii) /= 'no_name' ) then

            nomab = Stokes_name(ii)
            Length = len_trim( nomab )
            nomab(2:Length+1) = nomab(1:Length)
            nomab(1:1) = 'I'
            Length = Length + 1

          else

            if( Full_self_abs .and. ii == 1 ) j = j - 4
            nomab = nom_col_in(j)
            nomab(1:1) = 'I'
            Length = len_trim( nomab )
            do k = 2,Length-1
              if( nomab(k:k+1) == 'pp' ) exit
            end do
            if( abs( Stokes_param(1,ii) - 1 ) < 0.0001 ) then
              nomab(k:k) = 's'
            elseif( abs( Stokes_param(1,ii) + 1 ) < 0.0001 ) then
              nomab(k:k) = 'p'
            elseif( abs( Stokes_param(3,ii) - 1 ) < 0.0001 ) then
              nomab(k:k) = 'r'
            elseif( abs( Stokes_param(3,ii) + 1 ) < 0.0001 ) then
              nomab(k:k) = 'l'
            else
              nomab(k:k) = 'a'
            endif

            k = k + 1

            if( abs( Stokes_param(5,ii) ) < 0.0001 ) then
              nomab(k:Length-1) = nomab(k+1:Length)
              nomab(Length:Length) = ' '
            elseif( abs( Stokes_param(4,ii) ) < 0.0001 ) then
              nomab(k:k) = 's'
            elseif( abs( Stokes_param(4,ii) - pi/2 ) < 0.0001 ) then
              nomab(k:k) = 'p'
            else
              nomab(k:k) = 'a'
            endif

          endif

          nom_col_o(i) = nomab

          if( Cor_abs ) then
            i = i + 1
            do k = min(Length+1,Length_word),2,-1
              nomab(k:k) = nomab(k-1:k-1)
            end do
            nomab(1:2) = 'Ic'
            nom_col_o(i) = nomab
          endif

          if( Double_cor ) then
            i = i + 1
            nomab(2:2) = 'd'
            nom_col_o(i) = nomab
          endif

        end do
        if( Full_self_abs ) j = j + 4
      endif

    end do

  endif

  Close(2)
  
  deallocate( nom_col_in )

  if( Analyzer ) then
   
    nom_col(1:n_col) = nom_col_o(1:n_col)

  else
    n_col_o = i
  
    i_hk = 0
    k = 0
    kk = 0
    do i = 1,n_col_o
    
      if( i <= k ) cycle
      
      nomab = adjustl( nom_col_o(i) )
      if( nomab(1:2) /= 'I(' .and. nomab(1:2) /= 'I_' ) then
        k = k + 1
        kk = kk + 1
        nom_col(kk) = nom_col_o(i)  
        cycle
      endif
      Length = len_trim( nomab ) - 1
      
      do j = i+1,n_col_o
        nomac = adjustl( nom_col_o(j) )
  
        if( nomab(1:Length) /= nomac(1:Length) ) cycle
        i_hk = i_hk + 1
        n_index_hk(i_hk) = j - i
        k = k + 2 * n_index_hk(i_hk) 
        exit
      end do
      if( j > n_col_o ) then
        kk = kk + 1 
        nom_col(kk) = nom_col_o(i)  
      else
      
        do j = i, i + n_index_hk(i_hk) - 1
          kk = kk + 1  
          nomab = nom_col_o(j)
         if( nomab(1:2) /= 'mu' ) then 
            Length = len_trim( nomab )    ! one does not change the name of mu
            nomab(Length:Length) = ' '
          endif
          nom_col(kk) = nomab
        end do
      
      endif  

    end do

  endif
   
  return
end

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

! Tables giving the core level width 
! Edges K, L1, L2 and L3: comes from M. O. Krause et J. H. Oliver, J. Phys. Chem. Ref. Data 8, 329 (1979)
! but for Z = 1, up to 9 at K edge: from Elam 2002

function Tab_width(Eseuil,jseuil,nseuil,Numat)

  use declarations
  implicit none
  
  integer, parameter:: nam = 110
  
  integer:: jseuil, nseuil, Numat
  
  real(kind=db):: Es, Eseuil, Gamma_hole, Tab_width
  real(kind=db), dimension(nam):: GH_K1
  real(kind=db), dimension(10:nam):: GH_L1, GH_L2, GH_L3

  data GH_K1/ &
     0.02_db,  0.027_db, 0.036_db, 0.048_db, 0.065_db, 0.087_db, 0.116_db, 0.156_db, 0.209_db, 0.24_db, &
     0.30_db,  0.36_db,  0.42_db,  0.48_db,  0.53_db,  0.59_db,  0.64_db,  0.68_db,  0.74_db,  0.81_db, &
     0.86_db,  0.94_db,  1.01_db,  1.08_db,  1.16_db,  1.25_db,  1.33_db,  1.44_db,  1.55_db,  1.67_db, &
     1.82_db,  1.96_db,  2.14_db,  2.33_db,  2.52_db,  2.75_db,  2.99_db,  3.25_db,  3.52_db,  3.84_db, &
     4.14_db,  4.52_db,  4.91_db,  5.33_db,  5.77_db,  6.24_db,  6.75_db,  7.28_db,  7.91_db,  8.49_db, &
     9.16_db,  9.89_db, 10.6_db,  11.4_db,  12.3_db,  13.2_db,  14.1_db,  15.1_db,  16.2_db,  17.3_db,  &
    18.5_db,  19.7_db,  21.0_db,  22.3_db,  23.8_db,  25.2_db,  26.8_db,  28.4_db,  30.1_db,  31.9_db,  &
    33.7_db,  35.7_db,  37.7_db,  39.9_db,  42.1_db,  44.4_db,  46.8_db,  49.3_db,  52.0_db,  54.6_db,  &
    57.4_db,  60.4_db,  63.4_db,  66.6_db,  69.8_db,  73.3_db,  76.8_db,  80.4_db,  84.1_db,  88.0_db,  &
    91.9_db,  96.1_db, 100._db,  105._db,  109._db,  114._db,  119._db,  124._db,  129._db,  135._db,   &
   140._db,  145._db,  150._db,  156._db,  162._db,  168._db,  174._db,  181._db,  187._db,  193._db    /

  data GH_L1/                                                                                  0.1_db,  &
     0.2_db,   0.41_db,  0.73_db,  1.03_db,  1.26_db,  1.49_db,  1.58_db,  1.63_db,  1.92_db,  2.07_db, &
     2.21_db,  2.34_db,  2.41_db,  2.54_db,  2.62_db,  2.76_db,  2.79_db,  2.89_db,  3.06_db,  3.28_db, &
     3.38_db,  3.53_db,  3.79_db,  3.94_db,  4.11_db,  4.28_db,  4.44_db,  4.67_db,  4.71_db,  4.78_db, &
     3.94_db,  4.25_db,  4.36_db,  4.58_db,  4.73_db,  4.93_db,  4.88_db,  4.87_db,  5.00_db,  2.97_db, &
     3.13_db,  3.32_db,  3.46_db,  3.64_db,  3.78_db,  3.92_db,  4.06_db,  4.21_db,  4.34_db,  4.52_db, &
     4.67_db,  4.80_db,  4.91_db,  5.05_db,  5.19_db,  5.25_db,  5.33_db,  5.43_db,  5.47_db,  5.53_db, &
     5.54_db,  5.63_db,  5.58_db,  5.61_db,  6.18_db,  7.25_db,  8.30_db,  9.39_db, 10.5_db,  11.3_db,  &
    12.0_db,  12.2_db,  12.4_db,  12.6_db,  12.8_db,  13.1_db,  13.3_db,  13.4_db,  13.6_db,  13.7_db,  &
    14.3_db,  14.0_db,  14.0_db,  13.5_db,  13.3_db,  13.6_db,  13.8_db,  14.0_db,  14.3_db,  14.4_db,  &
    14.8_db,  15.1_db,  15.9_db,  16.2_db,  16.5_db,  16.8_db,  17.0_db,  17.2_db,  17.6_db,  18.1_db   /

  data GH_L2/                                                                                  0.0_db,  &
     0.0_db,  0.001_db, 0.004_db, 0.015_db, 0.032_db, 0.054_db, 0.083_db, 0.126_db, 0.152_db,  0.17_db, &
     0.19_db,  0.24_db,  0.26_db,  0.29_db,  0.34_db,  0.37_db,  0.43_db,  0.52_db,  0.62_db,  0.72_db, &
     0.83_db,  0.95_db,  1.03_db,  1.13_db,  1.21_db,  1.31_db,  1.43_db,  1.54_db,  1.65_db,  1.78_db, &
     1.87_db,  1.97_db,  2.08_db,  2.23_db,  2.35_db,  2.43_db,  2.57_db,  2.62_db,  2.72_db,  2.84_db, &
     3.00_db,  3.12_db,  3.25_db,  3.40_db,  3.51_db,  3.57_db,  3.68_db,  3.80_db,  3.89_db,  3.97_db, &
     4.06_db,  4.15_db,  4.23_db,  4.32_db,  4.43_db,  4.55_db,  4.66_db,  4.73_db,  4.79_db,  4.82_db, &
     4.92_db,  5.02_db,  5.15_db,  5.33_db,  5.48_db,  5.59_db,  5.69_db,  5.86_db,  6.00_db,  6.17_db, &
     6.32_db,  6.48_db,  6.67_db,  6.83_db,  7.01_db,  7.20_db,  7.47_db,  7.68_db,  7.95_db,  8.18_db, &
     8.75_db,  9.32_db,  9.91_db, 10.5_db,  10.9_db,  11.4_db,  11.8_db,  12.2_db,  12.7_db,  13.1_db,  &
    13.6_db,  14.0_db,  14.4_db,  14.9_db,  15.5_db,  15.9_db,  16.4_db,  17.0_db,  17.5_db,  18.1_db   /

  data GH_L3/                                                                                  0.0_db,  &
     0.0_db,  0.001_db, 0.004_db, 0.014_db, 0.033_db, 0.054_db, 0.087_db, 0.128_db, 0.156_db,  0.17_db, &
     0.19_db,  0.22_db,  0.24_db,  0.27_db,  0.32_db,  0.36_db,  0.43_db,  0.48_db,  0.56_db,  0.65_db, &
     0.76_db,  0.82_db,  0.94_db,  1.00_db,  1.08_db,  1.17_db,  1.27_db,  1.39_db,  1.50_db,  1.57_db, &
     1.66_db,  1.78_db,  1.91_db,  2.00_db,  2.13_db,  2.25_db,  2.40_db,  2.50_db,  2.65_db,  2.75_db, &
     2.87_db,  2.95_db,  3.08_db,  3.13_db,  3.25_db,  3.32_db,  3.41_db,  3.48_db,  3.60_db,  3.65_db, &
     3.75_db,  3.86_db,  3.91_db,  4.01_db,  4.12_db,  4.17_db,  4.26_db,  4.35_db,  4.48_db,  4.60_db, &
     4.68_db,  4.80_db,  4.88_db,  4.98_db,  5.04_db,  5.16_db,  5.25_db,  5.31_db,  5.41_db,  5.50_db, &
     5.65_db,  5.81_db,  5.98_db,  6.13_db,  6.29_db,  6.41_db,  6.65_db,  6.82_db,  6.98_db,  7.13_db, &
     7.33_db,  7.43_db,  7.59_db,  7.82_db,  8.04_db,  8.26_db,  8.55_db,  8.75_db,  9.04_db,  9.33_db, &
     9.61_db,  9.90_db, 10.1_db,  10.5_db,  10.8_db,  11.2_db,  11.6_db,  12.0_db,  12.3_db,  12.7_db    /

  select case(nseuil)

    case(1)   ! Seuil K
    
      if( Numat < 10 .or. Numat > nam ) then
   ! Formule ajustee sur les data de Krause et Oliver (correspond donc a une extrapolation)
        Es = Eseuil * Rydb
        if( Numat == 1 ) then
          Es = Es + 13.6_db
        elseif( Numat == 2 ) then
          Es = Es + 24.6_db
        endif
        Gamma_hole = 0.13007_db + 0.00011544_db * Es + 9.2262E-09_db * Es**2 - 3.2344E-14_db * Es**3
      else
        Gamma_hole = GH_K1( Numat )
      endif
      
    case(2)

      select case(jseuil)

        case(1)   ! Seuil L1
          if( Numat < 10 ) then
            Gamma_hole = 0._db
          elseif( Numat >= 10 .and. Numat <= nam ) then
            Gamma_hole = GH_L1( Numat )
          else
            Gamma_hole = 20._db
          endif

        case(2)   ! Seuil L2
          if( Numat < 10 ) then
            Gamma_hole = 0._db
          elseif( Numat >= 10 .and. Numat <= nam ) then
            Gamma_hole = GH_L2( Numat )
          else
            Gamma_hole = 19._db
          endif
! Anciennes valeurs avant 13/2/2015
!    Z < 26     Gamma_hole = 1.699*p  + 0.5001, p = ( Numat - 1 ) / 25.    
!     Fe        Gamma_hole = 1.5     
!     Co        Gamma_hole = 1.3    
!     Ni        Gamma_hole = 1.1
!  Z = 29,30,31 Gamma_hole = -0.9*p + 2.2,  p = ( Numat - 26 ) / 5.     
!  31 < Z < 61  Gamma_hole = 2.7*p + 1.3,   p = ( Numat - 31 ) / 29.     
!  60 < Z < 81  Gamma_hole = 1.5*p + 4.0,   p = ( Numat - 60 ) / 20.
!  Z > 80       Gamma_hole = 5*p + 5.5,     p = ( Numat - 80 ) / 15.

        case(3)   ! Seuil L3
          if( Numat < 10 ) then
            Gamma_hole = 0._db
          elseif( Numat >= 10 .and. Numat <= nam ) then
            Gamma_hole = GH_L3( Numat )
          else
            Gamma_hole = 13._db
          endif
! Anciennes valeurs avant 13/2/2015, comme L2 - 0.5
      end select

! For the edges M, N, O, values are not known
    case default

      select case(jseuil)

        case(0) ! Optic

          Gamma_hole = 0._db

        case(1,2,3)

          Gamma_hole = 0.2_db

        case(4,5,6,7)

          Gamma_hole = 0.1_db

    end select

  end select

  Tab_width = Gamma_hole / Rydb

  return
end

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

subroutine Seahdench(A,E_cut,Gamma_max,nelor,Elor,betalor)

  use declarations
  implicit none
  
  integer:: ie, nelor

  real(kind=db):: A, E, E_cut, Ep, lambda, Gamma_max
  real(kind=db), dimension(nelor):: Elor, betalor

  do ie = 1,nelor
    E = Elor(ie) - E_cut
    Ep = E
    if( E > 0._db .and. Ep > 0._db .and. ( A > 0._db .or. Gamma_max > 0._db ) ) then
      lambda = 0._db
      if( A > 0._db ) lambda = lambda + 1 / ( A * sqrt(Ep) )
      if( Gamma_max > 0._db ) lambda = lambda + sqrt(Ep) / Gamma_max
      betalor(ie) = sqrt(E) / lambda
    else
      betalor(ie) = 0._db
    endif
  end do

  return
end

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

subroutine gammarc(Ecent,Elarg,Gamma_max,E_cut,nelor, Elor,betalor)

  use declarations
  implicit none

  integer ie, nelor

  real(kind=db):: E, Ec, Ecent, E_cut, El, Elarg, Gamma_max, p
  real(kind=db), dimension(nelor):: Elor, betalor

  Ec = max( Ecent, 1.E-10_db )
  El = max( Elarg, 1.E-10_db )
  p = ( pi / 3 ) * Gamma_max / El

  do ie = 1,nelor
    E = Elor(ie) - E_cut
    if ( E <= 0._db ) then
      betalor(ie) = 0._db
    else
      betalor(ie) = Gamma_max * ( 0.5 + atan( p*(E/Ec - (Ec/E)**2)) / pi)
    endif
  end do

  return
end

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

! Convolution by a gaussian
! nj : Points au dela de la gamme pour diminuer les effets de bord de la convolution
  
subroutine Gaussi(Delta_gauss,Energ,nenerg,nj,Y)

  use declarations
  implicit none

  integer:: i, ie, je, n, ne1, ne2, nenerg, nj
  
  real(kind=db):: def, Delta_gauss, E, Efac, fac, p, pas, pdt, Yint
  real(kind=db), dimension(nenerg):: Energ, Y
  real(kind=db), dimension(1-nj:nenerg+nj):: dE, Ef, Gaus(1-nj:nenerg+nj), Xa(1-nj:nenerg+nj)

  Ef(1:nenerg) = Energ(1:nenerg)
  Xa(1:nenerg) = Y(1:nenerg)

! Creation des points au dela de la gamme
  ne1 = 1 - nj
!      ne2 = nenerg + nj
  def = Ef(2) - Ef(1)
  do ie = 0,ne1,-1
    Ef(ie) = Ef(ie+1) - def
    Xa(ie) = Xa(1)
  end do
  ne2 = nenerg

  de(ne1) = Ef(ne1+1) - Ef(ne1)
  do ie = ne1+1,ne2-1
    de(ie) = 0.5_db * ( Ef(ie+1) - Ef(ie-1) )
  end do
  de(ne2) = ( Ef(ne2) - Ef(ne2-1) )

  if( abs(Delta_gauss) < 1.e-10_db ) return

  do ie = 1,nenerg

    gaus(:) = 0._db
    Pdt = 0._db
    do je = ne1,ne2
      n = max( int( 10 * de(je) / Delta_gauss ), 1 )
      pas = de(je) / ( n + 1 )
      if ( je == ne1 ) then
        E = Ef(je) - 0.5_db * ( Ef(je+1) -  Ef(je) )
      else
        E = Ef(je) - 0.5_db * ( Ef(je) -  Ef(je-1) )
      endif
      do i = 1,n
        E = E + pas
        if( ( E < Ef(je) .and. je /= ne1 ) .or. je == ne2 ) then
          p = ( E - Ef(je-1) ) / ( Ef(je) - Ef(je-1) )
          Yint = (1 - p ) * Xa(je-1) + p * Xa(je)
        else
          p = ( E - Ef(je) ) / ( Ef(je+1) - Ef(je) )
          Yint = (1 - p ) * Xa(je) + p * Xa(je+1)
        endif
        fac = -0.5_db * ( ( E - Ef(ie) ) / Delta_gauss )**2
        if( fac > -600._db ) then
          efac = exp( fac )
          gaus(je) = gaus(je) + efac * Yint
          Pdt = Pdt + efac * de(je) / n
        endif
      end do
      gaus(je) = ( gaus(je) / n ) * de(je)
    end do
    Y(ie) = sum( gaus(ne1:ne2) ) / Pdt

  end do

  Where( Y < 1.e-99_db ) Y = 0._db
  return
end

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

! Calculation of the coefficients of the lorentzian

subroutine cflor(bb,betalor,E_cut,Eph,Elor,Energ,ie1,ie2,nef,nelor,nenerg,nenerge,XES)

  use declarations
  implicit none

  integer:: ie, ie1, ie2, j, nef, nelor, nenerg, nenerge
  
  logical:: XES

  real(kind=db):: def, E_cut, p
  real(kind=db), dimension(nenerge):: bb, Eph
  real(kind=db), dimension(nenerg):: Energ
  real(kind=db), dimension(nelor):: betalor, Elor

  Eph(1:nenerg) = Energ(1:nenerg)

! Creation des points au dela de la gamme
  def = Eph(2) - Eph(1)
  def = Eph(nenerg) - Eph(nenerg-1)
  do ie = nenerg+1,nenerge
    Eph(ie) = Eph(ie-1) + def
  end do

! Les etats en dessous de Fermi sont occupes
  do ie = 1,nenerg - 1
    if( Eph(ie) > E_cut - 1.e-10_db ) exit
  end do
  nef = ie

  if( XES ) nef = max(1, nef-1)

  if( XES ) then
    ie1 = 1
    ie2 = nef
  else
    ie1 = nef
    ie2 = nenerge
  endif

  do ie = 1,nenerge
    if( Eph(ie) <= Elor(1) ) then
      bb(ie) = betalor(1)
    elseif( Eph(ie) >= Elor(nelor) ) then
      bb(ie) = betalor(nelor)
    else
      do j = 2,nelor
        if( Elor(j) >= Eph(ie) ) exit
      end do
      p = ( Eph(ie) - Elor(j-1) ) / ( Elor(j) - Elor(j-1) )
      bb(ie) = p * betalor(j) + ( 1 - p ) * betalor(j-1)
    endif
    bb(ie) = 0.5_db * bb(ie) ! bb est en fait Gamma / 2
  end do

  return
end

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

! Calculation of the extrapolation after the grid for the atomic scattering

subroutine Extrapat(bb_nenerg,dampl,Ephoton,Eseuil,Extrap,fpp0,fprime_atom,icheck,nenerg,Numat)

  use declarations
  implicit none

  integer, parameter:: nim = 10000

  integer:: i, icheck, ie, nenerg, ni, Numat
  
  character(len=2):: Chemical_symbol

  complex(kind=db):: dampl(nenerg)

  logical:: Extrap, fprime_atom

  real(kind=db):: alfa, bb_nenerg, bb2, dde, E1, Eph, Ephi, Ephm, Eseuil, f, fp, fpp, fpp0
  real(kind=db):: ei(nim), ei1(nim), ei2(nim), Ephoton(nenerg), fppn(nim)

  if( fprime_atom ) then
    Eph = max( Eseuil - 200 / rydb, 10 / rydb )
    dde = 0.5 /rydb
    Ephm = Eseuil + 5000. / rydb
    Ephi = Eseuil + 500. / rydb
    write(3,110) Chemical_symbol(Numat)
    write(3,'(/A)') '   Energy_eV       f_prim      f_second'
    do i = 1,nim
      if( Eph < Ephi ) then
        Eph = Eph + dde
      else
        Eph = Eph + 10 * dde
      endif
      if( Eph > Ephm ) exit
      call fprime(Numat,Eph,fpp,fp)
      write(3,'(3f13.6)') Eph*rydb, fp, fpp
    end do
    if( .not. Extrap ) return
  endif

  Eph = Eseuil - 2 / rydb
  call fprime(Numat,Eph,fpp0,fp)

  Eph = Ephoton(nenerg)
  Ephm = Eseuil + 10000._db / rydb
  if( nenerg > 1 ) then
    dde = Ephoton(nenerg) - Ephoton(nenerg-1)
  else
    dde = 2 / rydb
  endif

  alfa = 1.02_db
  f = 0.5_db * alfa
  do i = 1,nim
    dde = alfa * dde
    Eph = Eph + dde
    if( Eph > ephm ) exit
    call fprime(Numat,Eph,fpp,fp)
    Ei(i) = Eph
    Ei2(i) = Ei(i) + f * dde
    Ei1(i) = Ei(i) - 0.5_db * dde
    fppn(i) = - ( fpp - fpp0 ) / pi
  end do
  ni = min(i-1,nim)

  bb2 = bb_nenerg**2

  do ie = 1,nenerg
    E1 = Ephoton(ie)
    if( abs(bb_nenerg) < 1.e-10_db ) then
      dampl(ie) = sum( fppn(1:ni) * log( ( Ei2(1:ni) - E1 ) / (Ei1(1:ni) - E1 ) ) )
    else
      dampl(ie) = sum( fppn(1:ni) * ( 0.5*log( ( ( Ei2(1:ni) - E1 )**2 + bb2 ) / ( ( Ei1(1:ni) - E1 )**2 + bb2 ) ) &
              - img * ( atan( ( Ei2(1:ni) - E1 ) / bb_nenerg ) - atan( ( Ei1(1:ni) - E1 ) / bb_nenerg ) ) ) )
    endif
  end do

  if( icheck > 1 ) then
    write(3,'(/A)') '  Energy      dampl_r     dampl_i'
    do ie = 1,nenerg
      write(3,'(f10.3,1p,2e13.5)') Ephoton(ie)*rydb, dampl(ie)
    end do
  endif

  return
  110 format(/1x,a2,' atomic scattering scattering amplitude')
end

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

subroutine extract_word(Length_line,motl,word,n_word,n_dim)

  use declarations
  
  implicit none
  
  integer:: i, i_word, ii, j, Length_line, n_dim, n_word

  character(len=Length_line):: motl
  character(len=length_word):: mot
  character(len=length_word), dimension(n_dim):: word

  i_word = 0
  i = 0
  do ii = 1,Length_line
    i = i + 1
    if( i > Length_line ) exit
    if( motl(i:i) == ' ' ) cycle
    i_word = i_word + 1
    mot = ' '
    do j = i,Length_line
      if( motl(j:j) == ' ' ) exit
      mot(j-i+1:j-i+1) = motl(j:j)
    end do
    i = j
    word(i_word) = mot
    if( i_word == n_word ) exit
  end do

  return
end

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

! Calculation of transmission

subroutine Cal_MuTrans(icheck,Es,Mu_mat_comp,Mus_mat,n_bir,n_stokes,n_mat_pol,nes,Sample_thickness,Stokes_param)

  use declarations
  implicit none

  integer:: i, i_bir, icheck, indp, ie, ipl, n_bir, n_mat_pol, n_stokes, nes

  complex(kind=db):: Ch, diag, exp_TS, exp_TS_i, Mu_iso, mu_pp, mu_ps, mu_sp, mu_ss, Sh, Tau, Trace

  complex(kind=db), dimension(2,2):: AM, Mu, Mat, Mat_pol
  complex(kind=db), dimension(nes,4,n_mat_pol):: Mu_mat_comp

  real(kind=db):: Cos_eta, Cos_2Bra, P1, P2, P3, Sample_thickness, Sin_eta

  real(kind=db), dimension(nes):: Es, Intensity
  real(kind=db), dimension(2,2):: Analyseur
  real(kind=db), dimension(nes,3):: Stokes_P
  real(kind=db), dimension(5,n_stokes):: Stokes_param
  real(kind=db), dimension(nes,n_stokes,n_bir,n_mat_pol):: Mus_mat

  if( icheck > 1 ) write(3,100)
  
  do ipl = 1,n_mat_pol

    do indp = 1,n_stokes   ! Boucle sur la polarisation

! Stokes parameter
      P1 = Stokes_param(1,indp)
      P2 = Stokes_param(2,indp)
      P3 = Stokes_param(3,indp)

! Polarization matrix
      Mat_pol(1,1) = 1._db + P1;  Mat_pol(1,2) = P2 - img*P3
      Mat_pol(2,1) = P2 + img*P3; Mat_pol(2,2) = 1._db - P1
      Mat_pol = Mat_pol / 2

! All is with conjugate complex (Bragg, F, ..)
      Mat_pol = Conjg( Mat_pol)

! Rotation angle of analyzor
      Cos_eta = Cos( Stokes_param(4,indp) )
      Sin_eta = Sin( Stokes_param(4,indp) )
! Bragg angle of analyzor
      Cos_2Bra = Cos( 2 * Stokes_param(5,indp) )

! Analyzor matrix
      Analyseur(1,1) = Cos_eta
      Analyseur(1,2) = - Sin_eta
      Analyseur(2,1) = Cos_2Bra * Sin_eta
      Analyseur(2,2) = Cos_2Bra * Cos_eta

      if( icheck > 1 ) then
        write(3,110) ipl, indp, P1, P2, P3, Stokes_param(4:5,indp) * 180 / pi
        do i = 1,2
          write(3,120) Mat_pol(i,:), Analyseur(i,:)
        end do 
      endif
      
      do ie = 1,nes

        do i_bir = 1,n_bir
        
          if( i_bir == 1 ) then 
            mu_ss = Mu_mat_comp(ie,1,ipl) 
            mu_sp = Mu_mat_comp(ie,2,ipl) 
            mu_ps = Mu_mat_comp(ie,3,ipl) 
            mu_pp = Mu_mat_comp(ie,4,ipl) 
            Mu_iso = exp( -0.25_db * ( mu_pp + mu_ss ) * Sample_thickness )
          elseif( i_bir == 2 ) then   ! Without birefringence
            mu_ss = cmplx( real( Mu_mat_comp(ie,1,ipl), db ), 0._db, db )  
            mu_sp = cmplx( real( Mu_mat_comp(ie,2,ipl), db ), 0._db, db ) 
            mu_ps = cmplx( real( Mu_mat_comp(ie,3,ipl), db ), 0._db, db ) 
            mu_pp = cmplx( real( Mu_mat_comp(ie,4,ipl), db ), 0._db, db )
          elseif( i_bir == 3 ) then   ! Without dichroism
            mu_ss = cmplx(0._db, aimag( Mu_mat_comp(ie,1,ipl) ), db )  
            mu_sp = cmplx(0._db, aimag( Mu_mat_comp(ie,2,ipl) ), db ) 
            mu_ps = cmplx(0._db, aimag( Mu_mat_comp(ie,3,ipl) ), db ) 
            mu_pp = cmplx(0._db, aimag( Mu_mat_comp(ie,4,ipl) ), db )
          else                        ! No off-diagonal components
            mu_ss = Mu_mat_comp(ie,1,ipl)  
            mu_sp = ( 0._db, 0._db ) 
            mu_ps = ( 0._db, 0._db ) 
            mu_pp = Mu_mat_comp(ie,4,ipl)
          endif 

          if( abs( mu_pp - mu_ss ) < 1.e-20_db .and. abs( mu_sp * mu_ps ) < 1.e-20_db ) then
            Mu(1,1) = 1._db
            Mu(2,1) = 0._db
            Mu(1,2) = 0._db
            Mu(2,2) = 1._db
          else
            Tau = 0.25_db * sqrt( (mu_pp - mu_ss)**2 + 4 * mu_sp * mu_ps )
            exp_TS = exp( Tau * Sample_thickness )
            exp_TS_i = 1 / exp_TS
            Ch  = 0.5_db * ( Exp_TS + Exp_Ts_i )  
            Sh = 0.5_db * ( Exp_TS - Exp_Ts_i ) 
            diag = 0.25_db * Sh * ( mu_pp - mu_ss ) / Tau
            Mu(1,1) = ch + diag
            Mu(2,1) = - 0.5_db * Sh * mu_ps / Tau
            Mu(1,2) = - 0.5_db * Sh * mu_sp / Tau
            Mu(2,2) = ch - diag
          endif
          
          Mu(:,:) =  Mu_iso * Mu(:,:)
     
          if( icheck > 1 .and. ie == 1 ) then
            if( i_bir == 1 ) then
              write(3,'(/A)') '   Mu Matrix'
            elseif( i_bir == 2 ) then
              write(3,'(/A)') '   Mu Matrix, without birefringence'
            elseif( i_bir == 3 ) then
              write(3,'(/A)') '   Mu Matrix, without dichroism'
            else
              write(3,'(/A)') '   Mu Matrix, without off-diagonal components'
            endif
            do i = 1,2
              write(3,130) Mu(i,:)
            end do 
          endif
       
          AM = Matmul( Analyseur, Mu )
     
          if( icheck > 1 .and. ie == 1 ) then
            write(3,'(/A)') '   A.Mu Matrix'
            do i = 1,2
              write(3,130) AM(i,:)
            end do 
          endif
     
          Mat = Transpose( Conjg( AM ) )
          Mat = Matmul( AM, Matmul( Mat_pol, Mat ) )
          Trace = real( Mat(1,1) + Mat(2,2), db )
     
          if( icheck > 1 .and. ie == 1 ) then
            write(3,'(/A)') '   A.Mu.P.Tr(Conj(A.Mu)) Matrix'
            do i = 1,2
              write(3,130) Mat(i,:)
            end do 
          endif
     
          Mus_mat(ie,indp,i_bir,ipl) = Trace

          if( i_bir == 1 .and. icheck > 0 ) then
 ! Writing of polarization matrix for outgoing beam

            Mat = Matmul( Mu, Matmul( Mat_pol, Transpose( Conjg( Mu ) ) ) )
 
 ! Because before we take the complex conjugate, so we go back to the conventional way for the Jones matrix          
            Mat(:,:) = conjg( Mat(:,:) )
 
            Trace = real( Mat(1,1) + Mat(2,2), db )
            Mat(:,:) = Mat(:,:) / Trace
            Intensity(ie) = Trace
            Stokes_P(ie,1) = real( Mat(1,1) - Mat(2,2), db )
            Stokes_P(ie,2) = real( Mat(2,1) + Mat(1,2), db )
            Stokes_P(ie,3) = aimag( Mat(2,1) - Mat(1,2) )
          endif
     
        end do
      end do ! end of loop on ie (energy)
      
      if( icheck > 0 ) then
        write(3,140)
        write(3,'(10x,1p,3e13.5)') Stokes_P(ie,1:3)
        write(3,150)
        do ie = 1,nes
          write(3,'(f10.3,1p,3e13.5,1x,2e13.5)') Es(ie) * rydb, Stokes_p(ie,:), Intensity(ie), &
                                                          - Log( Intensity(ie) ) / Sample_thickness
        end do
      endif
      
    end do
  end do

  return
 100 format(/' ---- Transmission calculation and polarization rotation --------',60('-')) 
 110 format(/' Wave vector direction',i2,', Stoke parameters',i3, / &
             ' P1, P2, P3 =',3f9.5,', Analyzer angles =',2f11.5, // & 
             '                  Polarization matrix                    Analyzer Matrix')
 120 format(5x,2(1x,2f10.7),5x,2f11.7)
 130 format(1p,5x,2(1x,2e13.5))
 140 format(/' Polarization parameters of the incoming beam (S1=1(-1) = sigma(pi), S2 = 1(-1) = +(-)(s+p)/2, S3=1(-1) Circ. R(L)', &
            /16x,'S1',11x,'S2',11x,'S3')
 150 format(/' Polarization matrix of the outgoing beam',/'    Energy',6x,'S1',11x,'S2',11x,'S3',9x,'Intensity',7x,'mu')
end

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

! Calculation of intensity versus Stokes parameters 

subroutine Cal_stokes(As,Icirc,nes,n_stokes,np_stokes, nphi,nphim,npldafs,Stokes_param)

  use declarations
  implicit none

  integer:: ie, indp, ip, ipl, is, nes, n_stokes, np_stokes, nphim, npldafs
  integer, dimension(npldafs):: nphi

  complex(kind = db):: fps, fpp, fsp, fss, Trace

  complex(kind=db), dimension(2,2):: AF, F, Mat, Mat_pol
  complex(kind=db), dimension(nes,nphim,npldafs):: As

  real(kind=db):: Cos_eta, Cos_2Bra, P1, P2, P3, Sin_eta

  real(kind=db), dimension(2,2):: Analyseur
  real(kind=db), dimension(5,n_stokes):: Stokes_param
  real(kind=db), dimension(nes,nphim,np_stokes):: Icirc

  is = 0

  do ipl = 1,npldafs,4

    do indp = 1,n_stokes   ! Loop on polarization

      is = is + 1

! Stokes parameters
      P1 = Stokes_param(1,indp)
      P2 = Stokes_param(2,indp)
      P3 = Stokes_param(3,indp)

! Polarisation matrix
      Mat_pol(1,1) = 1._db + P1;  Mat_pol(1,2) = P2 - img*P3
      Mat_pol(2,1) = P2 + img*P3; Mat_pol(2,2) = 1._db - P1
      Mat_pol = Mat_pol / 2

! AA is with conjugate complex (Bragg, F, ..)
      Mat_pol = Conjg( Mat_pol)

! Rotation angle of analyzor
      Cos_eta = Cos( Stokes_param(4,indp) )
      Sin_eta = Sin( Stokes_param(4,indp) )
! Bragg angle of analyzor
      Cos_2Bra = Cos( 2 * Stokes_param(5,indp) )

! Analyzor matrix
      Analyseur(1,1) = Cos_eta
      Analyseur(1,2) = - Sin_eta
      Analyseur(2,1) = Cos_2Bra * Sin_eta
      Analyseur(2,2) = Cos_2Bra * Cos_eta

      do ip = 1,nphi(ipl)

        do ie = 1,nes

          fss = As(ie,ip,ipl)
          fsp = As(ie,ip,ipl+1)
          fps = As(ie,ip,ipl+2)
          fpp = As(ie,ip,ipl+3)

          F(1,1) = fss; F(1,2) = fps
          F(2,1) = fsp; F(2,2) = fpp

          AF = Matmul( Analyseur, F )

          Mat = Transpose( Conjg( AF ) )
          Mat = Matmul( AF, Matmul( Mat_pol, Mat ) )
          Trace = Mat(1,1) + Mat(2,2)

          Icirc(ie,ip,is) = Real ( Trace , db )

        end do

      end do
    end do
  end do

  return
end

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

! Calculation of intensity corrected by self-absorption taking into account birefringence

subroutine Corr_abs(angxyz,As,axyz,c_micro,d_dead,Double_cor,Eseuil,fpp_avantseuil,hkl_dafs,hkl_S,Icor,Icirccor, &
          Icircdcor,Idcor,mus,n_stokes,nes,np_stokes,nphi,nphim, npldafs,Self_abs,Stokes_param)

  use declarations
  implicit none

  integer:: i, ie, indp, ip, ipl, ipl_pas, is, j, k, n_stokes, nes, np_stokes, nphim, npldafs
  integer, dimension(3):: hkl, hkl_S
  integer, dimension(npldafs):: nphi
  integer, dimension(3,npldafs):: hkl_dafs

  logical:: Check, Double_cor, Self_abs

  complex(kind = db):: dmu_i, dmu_s, Expo, fps, fpp, fsp, fss, mu_i, mu_pp_i, mu_pp_s, mu_ps_i, mu_ps_s, mu_s, &
       mu_sp_i, mu_sp_s, mu_ss_i, mu_ss_s, mups_i, mups_s, musp_i, musp_s, t_i, t_s, tau_i, tau_s, Trace

  complex(kind=db), dimension(2):: Ex_i, Ex_s
  complex(kind=db), dimension(4):: Ex
  complex(kind=db), dimension(2,2):: F, M1, M2, M3, Mat_pol
  complex(kind=db), dimension(2,2,2):: Tr_i, Tr_s
  complex(kind=db), dimension(2,2,4):: ATFT, TFT
  complex(kind=db), dimension(nes,nphim,npldafs):: As
  complex(kind=db), dimension(nes,nphim,npldafs,2):: mus

  real(kind=db):: alfa_S, Attenuation, Brag_A, c_micro, corr_ref, Correction, Cos_eta, Cos_2Bra, csi, css, d_dead, Eseuil, eta, &
    fpp_avantseuil, mu_i_m, mu_s_m, mum, mum_0, P1, P1A, P2, P3, Resultat, Sin_eta, theta_B

  real(kind=db), dimension(3):: angxyz, axyz
  real(kind=db), dimension(2,2):: Analyseur
  real(kind=db), dimension(5,n_stokes):: Stokes_param
  real(kind=db), dimension(nes,nphim,npldafs):: Icor, Idcor
  real(kind=db), dimension(nes,nphim,np_stokes):: Icirccor, Icircdcor

  Check = .false.

  if( Self_abs ) then
    ipl_pas = 1   ! not Full_self_abs
  else
    ipl_pas = 4
  endif

  is = 0
  
  do ipl = 1,npldafs,ipl_pas

    hkl(:) = hkl_dafs(:,ipl)
    call Angle_Bragg(alfa_S,angxyz,axyz,Eseuil,hkl,hkl_S,theta_B)

    csi = 1 / sin( theta_B + alfa_S )
    css = 1 / sin( theta_B - alfa_S )
    mum_0 = c_micro * fpp_avantseuil * ( csi + css )
    corr_ref = mum_0 / exp( - mum_0 * d_dead )

    do ip = 1,nphi(ipl)

      do ie = 1,nes

        if( Self_abs .or. Double_cor ) mum = real( mus(ie,ip,ipl,1), db )*csi + real( mus(ie,ip,ipl,2), db )*css
        if( Double_cor ) Correction = mum_0 / mum

        if( self_abs ) then
          Attenuation = exp( - mum * d_dead ) * corr_ref
          Icor(ie,ip,ipl) = ( Attenuation / mum ) * abs( As(ie,ip,ipl) )**2
          if( Double_cor ) Idcor(ie,ip,ipl) = Icor(ie,ip,ipl) * Correction
          cycle
        endif

        mu_ss_i = mus(ie,ip,ipl,1)
        mu_ss_s = mus(ie,ip,ipl,2)
        mu_sp_i = mus(ie,ip,ipl+1,1)
        mu_sp_s = mus(ie,ip,ipl+1,2)
        mu_ps_i = mus(ie,ip,ipl+2,1)
        mu_ps_s = mus(ie,ip,ipl+2,2)
        mu_pp_i = mus(ie,ip,ipl+3,1)
        mu_pp_s = mus(ie,ip,ipl+3,2)

        fss = As(ie,ip,ipl)
        fsp = As(ie,ip,ipl+1)
        fps = As(ie,ip,ipl+2)
        fpp = As(ie,ip,ipl+3)

        mu_i = ( mu_ss_i + mu_pp_i ) * csi / 4
        mu_s = ( mu_ss_s + mu_pp_s ) * css / 4

        mum = 0.5_db*real( mu_ss_i + mu_pp_i, db )*csi + 0.5_db*real( mu_ss_s + mu_pp_s, db )*css

        Attenuation = exp( - mum * d_dead ) * corr_ref

        tau_i = sqrt( ( mu_pp_i - mu_ss_i )**2 + 4*mu_sp_i*mu_ps_i )
        tau_s = sqrt( ( mu_pp_s - mu_ss_s )**2 + 4*mu_sp_s*mu_ps_s )

        t_i = csi * tau_i / 4
        t_s = css * tau_s / 4

        if( abs(tau_i) > eps10 ) then
          mups_i = 2*mu_ps_i / tau_i
          musp_i = 2*mu_sp_i / tau_i
          dmu_i = ( mu_pp_i - mu_ss_i ) / tau_i
        else
          mups_i = (0._db, 0._db)
          musp_i = (0._db, 0._db)
          dmu_i = (0._db, 0._db)
        endif

        if( abs(tau_s) > eps10 ) then
          mups_s = 2*mu_ps_s / tau_s
          musp_s = 2*mu_sp_s / tau_s
          dmu_s = ( mu_pp_s - mu_ss_s ) / tau_s
        else
          mups_s = (0._db, 0._db)
          musp_s = (0._db, 0._db)
          dmu_s = (0._db, 0._db)
        endif

! Transmittance matrices
! They have cut in 2 parts. The total transmitance matrix is given by
!   T = 0.5*(exp(Tau*l)*Tr(1) + exp(-Tau*l)*Tr(2) )
        Tr_i(1,1,1) = 1 + dmu_i; Tr_i(1,2,1) = - musp_i
        Tr_i(2,1,1) = - mups_i;  Tr_i(2,2,1) = 1 - dmu_i

        Tr_i(1,1,2) = 1 - dmu_i; Tr_i(1,2,2) = musp_i
        Tr_i(2,1,2) = mups_i;    Tr_i(2,2,2) = 1 + dmu_i

        Tr_i(:,:,:) = Tr_i(:,:,:) / 2

        Tr_s(1,1,1) = 1 + dmu_s; Tr_s(1,2,1) = - musp_s
        Tr_s(2,1,1) = - mups_s;  Tr_s(2,2,1) = 1 - dmu_s

        Tr_s(1,1,2) = 1 - dmu_s; Tr_s(1,2,2) = musp_s
        Tr_s(2,1,2) = mups_s;    Tr_s(2,2,2) = 1 + dmu_s

        Tr_s(:,:,:) = Tr_s(:,:,:) / 2

! Exponent of the exponentials

        Ex_i(1) = - mu_i + t_i;  Ex_i(2) = - mu_i - t_i

        Ex_s(1) = - mu_s + t_s;  Ex_s(2) = - mu_s - t_s

! Scattering amplitude

        F(1,1) = fss; F(1,2) = fps
        F(2,1) = fsp; F(2,2) = fpp

! Product
        k = 0
        do i = 1,2
          do j = 1,2
            k = k + 1
            M1(:,:) = Tr_s(:,:,i)
            M2(:,:) = Tr_i(:,:,j)
            M3 = Matmul( M1, Matmul( F, M2 ) )
            TFT(:,:,k) = M3
            Ex(k) = Ex_s(i) + Ex_i(j)
          end do
        end do

        if( check ) then
          Write(6,'(A)') ' tau_i tau_s'
          Write(6,'(1p,8e13.5)') tau_i, tau_s, t_i, t_s
          Write(6,'(A)') ' Tr_i'
          do i = 1,2
            Write(6,'(1p,8e13.5)') Tr_i(i,:,1), Tr_i(i,:,2)
          end do
          Write(6,'(A)') ' Tr_s'
          do i = 1,2
            Write(6,'(1p,8e13.5)') Tr_s(i,:,1), Tr_s(i,:,2)
          end do
          Write(6,'(A)') ' Ex_i'
            Write(6,'(1p,8e13.5)') Ex_i(:)
          Write(6,'(A)') ' Ex_s'
            Write(6,'(1p,8e13.5)') Ex_s(:)
          Write(6,'(A)') ' F'
          do i = 1,2
            Write(6,'(1p,8e13.5)') F(i,:)
          end do
          do k = 1,4
            Write(6,'(a9,i2)') ' TFT, k =', k
            Write(6,'(1p,8e13.5)') TFT(1,:,k), Ex(k)
            Write(6,'(1p,8e13.5)') TFT(2,:,k)
          end do
        endif

        do indp = 1,4+n_stokes  ! Loop over polarization and analyzor 

          select case(indp)

            case(1)
              P1 = 1._db; P2 = 0._db; P3 = 0._db
              eta = 0._db; Brag_A = pi / 4

            case(2)
              P1 = 1._db; P2 = 0._db; P3 = 0._db
              eta = pi / 2; Brag_A = pi / 4

            case(3)
              P1 = -1._db; P2 = 0._db; P3 = 0._db
              eta = 0._db; Brag_A = pi / 4

            case(4)
              P1 = -1._db; P2 = 0._db; P3 = 0._db
              eta = pi / 2; Brag_A = pi / 4

            case default

!              is = ipl + indp - 5

              eta = Stokes_param(4,indp-4)
              Brag_A = Stokes_param(5,indp-4)

! Stokes parameters
              P1 = Stokes_param(1,indp-4)
              P2 = Stokes_param(2,indp-4)
              P3 = Stokes_param(3,indp-4)

! Parametre de Stokes P1 equivalent pour l'analyseur
              P1A = Cos( 2 * eta ) * Sin( 2 * Brag_A )

          end select

! Polarization matrix
          Mat_pol(1,1) = 1._db + P1;  Mat_pol(1,2) = P2 - img*P3
          Mat_pol(2,1) = P2 + img*P3; Mat_pol(2,2) = 1._db - P1
          Mat_pol = Mat_pol / 2

! Tout est avec le complexe conjugue (Bragg, F, ..)
          Mat_pol = Conjg( Mat_pol)

! Angle de rotation de l'analyseur
          Cos_eta = Cos( eta )
          Sin_eta = Sin( eta )
! Angle de Bragg de l'analyseur
          Cos_2Bra = Cos( 2 * Brag_A )

! Matrice Analyseur
          Analyseur(1,1) = Cos_eta
          Analyseur(1,2) = - Sin_eta
          Analyseur(2,1) = Cos_2Bra * Sin_eta
          Analyseur(2,2) = Cos_2Bra * Cos_eta

          if( Check ) then
            Write(6,'(/A)') ' Analyseur'
            Write(6,'(2f8.5)') Analyseur(1,:)
            Write(6,'(2f8.5)') Analyseur(2,:)
          endif

          do i = 1,4
            M1(:,:) = TFT(:,:,i)
            M2 = Matmul( Analyseur, M1 )
            ATFT(:,:,i) = M2
          end do

          Resultat = 0._db
          do i = 1,4
            do j = 1,4
              M1(:,:) = ATFT(:,:,i)
              M2(:,:) = ATFT(:,:,j)
              M2 = Transpose( Conjg( M2 ) )
              M3 = Matmul( M1, Matmul( Mat_pol, M2 ) )
              Trace = M3(1,1) + M3(2,2)
              Expo = Ex(i) + Conjg( Ex(j) )
! Signe moins car integrale de 0 a l'infini = - Primitive(0)
              Resultat = Resultat - Real ( Trace / Expo, db )
            end do
          end do
          Resultat = Resultat * Attenuation

          if( check ) then
            Write(6,'(/A)') ' Polar'
            Write(6,'(4f8.3)') Mat_pol(1,:)
            Write(6,'(4f8.3)') Mat_pol(2,:)
            Write(6,'(/A)') ' Attenuation, Resultat'
            Write(6,'(1p,5e13.5)') Attenuation, Resultat
          endif

          select case(indp)

            case(1)  ! sigma - sigma
              Icor(ie,ip,ipl) = Resultat
              if( Double_cor ) Idcor(ie,ip,ipl) = Resultat * Corr_ref / Real( mu_ss_i*csi + mu_ss_s*css, db)
            case(2)  ! pi - sigma
              Icor(ie,ip,ipl+2) = Resultat
              if( Double_cor ) Idcor(ie,ip,ipl+2) = Resultat * Corr_ref / Real( mu_pp_i*csi + mu_ss_s*css, db)
            case(3)  ! sigma - pi
              Icor(ie,ip,ipl+1) = Resultat
              if( Double_cor ) Idcor(ie,ip,ipl+1) = Resultat * Corr_ref / Real( mu_ss_i*csi + mu_pp_s*css, db)
            case(4)  ! pi - pi
              Icor(ie,ip,ipl+3) = Resultat
                if( Double_cor ) Idcor(ie,ip,ipl+3) = Resultat * Corr_ref / Real( mu_pp_i*csi + mu_pp_s*css, db)

            case default

              is = ( ( ipl + 3 ) / 4 - 1 ) * n_stokes + indp - 4
              
              Icirccor(ie,ip,is) = Resultat

              if( Double_cor ) then
                mu_i_m = ( 1 + P1 ) * Real( mu_ss_i, db ) + ( 1 - P1 ) * Real( mu_pp_i, db )

                mu_s_m = ( 1 + P1A ) * Real( mu_ss_s, db ) + ( 1 - P1A ) * Real( mu_pp_s, db )

                Icircdcor(ie,ip,is) = Resultat * Corr_ref * 2 / ( mu_i_m * csi + mu_s_m * css )
              endif

          end select

        end do ! Boucle indp

      end do
    end do
  end do

  return
end

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

! Effect of temperature using the Debye model

subroutine Debye_effect(Abs_U_iso,Abs_U_iso_inp,Energ,Energphot,Ephoton,Eseuil,fac,ifich,jseuil,n_col,nef,nenerg,nfich, &
                        nom_col,Numat,nxan,Shift_U_iso,U_iso_man,V0muf,Xanes)

  use declarations
  implicit none

  character(len=Length_word):: nomab
  character(len=Length_word), dimension(n_col):: nom_col

  integer:: i, ie, ie0, ifich, jseuil, Length, n_col, nef, nenerg, nfich, nxan, Z
  integer, dimension(nfich):: Numat
  
  logical:: Energphot, U_iso_man

  real(kind=db):: Abs_U_iso_inp, c, Conv_Mbarn_nelec, Delta_E, E, E_medium, Eph, f, fac, fp, fpp, fpp0, &
                  p, Shift_BE, Shift_U_iso, Xanes_atom_av
  real(kind=db), dimension(nxan):: Ratio, Xanes_av 
  real(kind=db), dimension(nenerg):: Energ, Ephoton, Xanes_atom 
  real(kind=db), dimension(nenerg,nxan):: Xanes 
  real(kind=db), dimension(nfich):: Abs_U_iso, Eseuil, V0muf 

  Z = Numat(ifich) 
  
!  Abs_U_iso = <u^2> is in A^2
  c = 0.04 * Rydb * ( m_electron * e_electron / hbar**2 )

  do i = 1,nxan
    nomab = nom_col(i)
    nomab = adjustl(nomab)
    if( nomab /= 'XANES_atom' ) cycle
    Xanes_atom(:) = Xanes(:,i)
    exit   
  end do

 ! Calculation of atomic spectra when not existing
  if( i > nxan ) then

    call Shift_from_fprime_binding_energy(Eseuil(ifich),fpp0,jseuil,Shift_BE,Z) 
    do ie = 1,nenerg
      Eph = Ephoton(ie) - Shift_U_iso + Shift_BE
      call fprime(Z,Eph,fpp,fp)
      f = fac / conv_mbarn_nelec(Eph)
      Xanes_atom(ie) = f * ( fpp - fpp0 )
    end do
  
  endif
  
! This part is also used for X-ray Raman. In that case, the Xanes_atom must be renormalized to
! a Raman_atom. Because here one does not know if it is Raman or XANES one checks the normalization.
  Ratio(:) = 1._db
  nomab = nom_col(1)
  nomab = adjustl(nomab)
  if( nomab(1:2) == 'q=' ) then  ! Test to check if it Raman X
    E_medium = Energ(nenerg) / 2
    do ie = 1,nenerg
      if( Energ(ie) > E_medium ) exit
    end do
    ie0 = ie
    if(  ie0 > 2 .and. ie0 < nenerg - 2 ) then
      Xanes_atom_av = 0._db
      Xanes_av(:) = 0._db
      do ie = ie0,nenerg
        if( ie == ie0 ) then
          Delta_E = ( Energ(ie+1) - Energ(ie) ) / 2 + Energ(ie) - E_medium
        elseif( ie == nenerg ) then
          Delta_E = ( Energ(ie) - Energ(ie-1) ) / 2
        else
          Delta_E = ( Energ(ie+1) - Energ(ie-1) ) / 2
        endif
        Xanes_atom_av = Xanes_atom_av + Delta_E * Xanes_atom(ie)
        if( ie < nef ) cycle  ! it is before convolution, thus before E_fermi, it is zero
        do i = 1,nxan
          nomab = nom_col(i)
          nomab = adjustl(nomab)
          Length = len_trim(nomab)
          if( nomab(Length-1:Length) == 'l0' .or. nomab(Length-1:Length) == 'l2' .or. nomab(Length-1:Length) == 'l3' ) cycle
          Xanes_av(i) = Xanes_av(i) + Delta_E * Xanes(ie,i)   
        end do
      end do
      Xanes_atom_av = Xanes_atom_av / E_medium
      do i = 1,nxan
        nomab = nom_col(i)
        nomab = adjustl(nomab)
        Length = len_trim(nomab)
        if( nomab(Length-1:Length) == 'l0' .or. nomab(Length-1:Length) == 'l2' .or. nomab(Length-1:Length) == 'l3' ) then
          Ratio(i) = 0._db
        else
          Xanes_av(i) = Xanes_av(i) / E_medium
          Ratio(i) = Xanes_av(i) / Xanes_atom_av
        endif
      end do
    endif
  endif

  do ie = 1,nenerg
    if( Energphot ) then
      E = Energ(ie) - Eseuil(ifich) - V0muf(ifich)
    else
      E = Energ(ie) - V0muf(ifich)
    endif
    if( E < eps10 ) cycle

    if( U_iso_man ) then 
      p = exp( - c * Abs_U_iso_inp * E )
    else
      p = exp( - c * Abs_U_iso(ifich) * E )
    endif
    do i = 1,nxan
      nomab = nom_col(i)
      nomab = adjustl(nomab)
      Length = len_trim(nomab)
      if( nomab(Length-1:Length) == 'l0' .or. nomab(Length-1:Length) == 'l2' .or. nomab(Length-1:Length) == 'l3' ) then
        cycle
      elseif( nomab(1:3) == 'dic' .or. nomab(1:3) == 'Dic' ) then
        Xanes(ie,i) = p * Xanes(ie,i)
      else
        Xanes(ie,i) = p * Xanes(ie,i) + (1 - p) * Ratio(i) * Xanes_atom(ie)
      endif
    end do
  end do
           
  return
end

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

! fprime.f90 has tabulated edge energy different from fdmnes (see fprime_data).
! This routine calculate the shift between them

subroutine Shift_from_fprime_binding_energy(Eseuil,fpp0,jseuil,Shift_BE,Z) 

  use declarations
  implicit none

  integer:: i, is, j, jseuil, Z
  
  logical:: Special
  
  real(kind=db):: Binding_energy, Delta, Eph, Eseuil, fp, fpp0, fpp1, Shift_BE 

! Special case for L1 Oxygen case
  if( Z == 8 .and. abs( Eseuil*rydb - 29.79_db ) < 1._db ) then
    Binding_energy = 23.5_db / Rydb 
    Shift_BE = Binding_energy - Eseuil
    return
  endif
    
! This edge is special because between M5 and M4, fpp is increasing
  Special = ( Z == 91 .or. Z == 92 ) .and. jseuil == 4 
  Special = .false.
  if( Special ) then
    is = - 1
  else
    is = 1
  endif
  
  Eph = Eseuil - is
  Eph = max( Eph,  0.5_db )
  call fprime(Z,Eph,fpp0,fp)

  Delta = 0.8_db / Rydb
  do j = 1,3
    Loop_i: do i = 1,100
      Eph = Eph + is * Delta
      call fprime(Z,Eph,fpp1,fp)
      if( ( fpp1 > fpp0 .and. .not. Special ) .or. ( fpp1 < fpp0 .and. Special ) ) then
        Eph = Eph - is * Delta
        Delta = Delta / 2
        exit Loop_i
      else
        Binding_energy = Eph
        fpp0 = fpp1
      endif
    end do Loop_i
  end do
    
  Shift_BE = Binding_energy - Eseuil
   
  return
end

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

subroutine Debye_effect_a(Abs_U_iso,Abs_U_iso_inp,Adafs,Bragg_abs_T,Energ,Energphot,Ephoton,Eseuil,fac,ifich,jseuil, &
                        nenerg,nfich,Dafs_cal,nphim,npldafs,Numat,Shift_U_iso,U_iso_man,V0muf) 
  use declarations
  implicit none

  integer:: i, icheck, ie, ifich, jseuil, nenerg, nfich, nphim, npldafs, Z
  integer, dimension(nfich):: Numat

  complex(kind=db), dimension(nenerg,nphim,npldafs):: Adafs 
  complex(kind=db), dimension(npldafs):: Bragg_abs_T 
  complex(kind=db), dimension(nenerg):: Dafs_atom 
  
  logical:: Dafs_cal, Energphot, U_iso_man

  real(kind=db):: Abs_U_iso_inp, c, Conv_Mbarn_nelec, E, Eph, f, fac, fp, fpp, fpp0, p, Shift_BE, Shift_U_iso
  real(kind=db), dimension(nenerg):: Energ, Ephoton 
  real(kind=db), dimension(nfich):: Abs_U_iso, Eseuil, V0muf 

  Z = Numat(ifich) 
  
!  Abs_U_iso = <u^2> is in A^2
  c = 0.04 * Rydb * ( m_electron * e_electron / hbar**2 )

 ! Calculation of atomic spectra when not existing

  call Shift_from_fprime_binding_energy(Eseuil(ifich),fpp0,jseuil,Shift_BE,Z) 

  f = fac / pi
  
  do ie = 1,nenerg
    Eph = Ephoton(ie) - Shift_U_iso + Shift_BE
    call fprime(Z,Eph,fpp,fp)
    if( .not. Dafs_cal ) f = fac / ( pi * Conv_Mbarn_nelec(Eph) ) 
    Dafs_atom(ie) = f * cmplx( fpp - fpp0, 0._db, db )
  end do

  icheck = 1
  if( icheck > 1 ) then
    open(99, file = 'Debye_out.txt')
    write(99,'(A)') '   Energy     Dafs_r       Dafs_i      Dafs_a_r    Dafs_a_i'
    do ie = 1,nenerg
      write(99,'(f10.3,1p,4e13.5)') Energ(ie)*Rydb, Adafs(ie,1,5), Dafs_atom(ie)*conjg(Bragg_abs_T(5))
    end do
    stop
  endif

  do ie = 1,nenerg
    if( Energphot ) then
      E = Energ(ie) - Eseuil(ifich) - V0muf(ifich)
    else
      E = Energ(ie) - V0muf(ifich)
    endif
    if( E < eps10 ) cycle
    
    if( U_iso_man ) then 
      p = exp( - c * Abs_U_iso_inp * E )
    else
      p = exp( - c * Abs_U_iso(ifich) * E )
    endif
    do i = 1,npldafs
      Adafs(ie,:,i) = p * Adafs(ie,:,i) + (1 - p) * conjg( Bragg_abs_T(i) ) * Dafs_atom(ie)
    end do
  end do
   
  return
end

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

subroutine Angle_Bragg(alfa_S,angxyz,axyz,Eseuil,hkl,hkl_S, theta_B)

  use declarations
  implicit none

  integer:: i, j
  integer, dimension(3):: hkl, hkl_S

  real(kind=db):: alfa_S, detmat, dhkl, Eseuil, fac, konde, qkn, qkn_S, theta_B, vol
  real(kind=db), dimension(3):: angxyz, axyz, cosdir, hklred, qk, qk_S, vx, vy, vz, wx, wy, wz
  real(kind=db), dimension(0:3):: det
  real(kind=db), dimension(3,3):: Cubmat, mat

  character(len=5):: Struct

  call cal_cubmat(angxyz,Cubmat,struct)

  konde = 0.5_db * alfa_sf * Eseuil

  cosdir(:) = cos( angxyz(:) )

! Base du reseau direct, exprimee dans une base orthonormee
  vx(:) = Cubmat(:,1) * axyz(:)
  vy(:) = Cubmat(:,2) * axyz(:)
  vz(:) = Cubmat(:,3) * axyz(:)

! wx, wy, wz : base du reseau reciproque exprimee dans la meme base
  call Prodvec(Wx,vy,vz)

  vol = sum( wx(:) * vx(:) )
  Wx(:) = Wx(:) / vol
  call Prodvec(Wy,vz,vx)
  Wy(:) = Wy(:) / vol
  call Prodvec(Wz,vx,vy)
  Wz(:) = Wz(:) / vol

  qk_S(:) = hkl_S(1) * wx(:) + hkl_S(2) * wy(:) + hkl_S(3) * wz(:)
  qkn_S = sqrt( sum( qk_S(:)**2 ) )

  if( abs(qkn_S) > eps10 ) qk_S(:) = qk_S(:) / qkn_S

  if( hkl(1) == 0 .and. hkl(2) == 0 .and. hkl(3) == 0 ) then

    dhkl = 0._db
    theta_B = 0._db

  else

    hklred(:) = hkl(:) / axyz(:)
    do i = 0,3
      do j = 1,3
        mat(j,j) = 1._db
      end do
      mat(1,2) = cosdir(3)
      mat(1,3) = cosdir(2)
      mat(2,3) = cosdir(1)
      mat(2,1) = mat(1,2)
      mat(3,1) = mat(1,3)
      mat(3,2) = mat(2,3)
      if( i > 0 ) mat(:,i) = hklred(:)
      det(i) = detmat(mat)
    end do
! Distance interplan
    dhkl = sqrt( det(0) / sum( hklred(1:3) * det(1:3) ) )
    fac = abs( pi / ( konde * dhkl ) )

    if( fac > 1._db ) then
      call write_error
      write(6,120) hkl(:)
      stop
    endif

    theta_B = asin( pi / ( konde * dhkl ) )

  endif

! Calcul angle avec plan de surface

  if( qkn_S > eps10 ) then

    qk(:) = hkl(1) * wx(:) + hkl(2) * wy(:) + hkl(3) * wz(:)
    qkn = sqrt( sum( qk(:)**2 ) )
    if( abs(qkn) > eps10 ) then
      qk(:) = qk(:) / qkn
      alfa_S = acos( sum( qk(:)*qk_S(:) ) )
    else
      alfa_S = 0._db
    endif
  else
    alfa_S = 0._db
  endif

  return
  120 format(//' The reflection number',i3,' : (h,k,l) = (',3i3,') does not exist at this energy !'//)
end

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

! Calculation of scttaring amplitudes foror SRXRD.
! One applies self-absorption and truncation for the Bulk
! The possible occupancy < 1 are already included in the absorption in Write_coabs and Abs_NR_layer routines

subroutine Total_amplitude_SRXRD(Abs_cap,Abs_NR_SRXRD,Abs_NR_SRXRD_abs,Abs_NR_SRXRD_abs_bulk,Abs_NR_SRXRD_bulk,As, &
      As_SRXRD,As_SRXRD_bulk,Es,icheck,L_dafs,Length_abs,n_z_SRXRD,n_z_SRXRD_abs,n_z_SRXRD_abs_bulk,n_z_SRXRD_bulk,nes, &
      Non_resonant_bulk,nphim,npldafs,npldafs_NR,Pos_z_SRXRD,Pos_z_SRXRD_abs,Pos_z_SRXRD_bulk,Pos_z_SRXRD_abs_bulk, &
      Sum_Bragg_abs_f0,Sum_Bragg_abs_f0_bulk,Sum_Bragg_nonabs_f,Sum_Bragg_nonabs_f_bulk,Sum_Bragg_NR_Bulk,Sum_Bragg_NR_T, &
      Xs_SRXRD,Xs_SRXRD_bulk)

  use declarations
  implicit none
  
  integer:: i, icheck, ie, ipl, j, n_z_SRXRD, n_z_SRXRD_abs, n_z_SRXRD_abs_bulk, n_z_SRXRD_bulk, nes, nphim, npldafs, npldafs_NR
  
  logical:: Non_resonant_bulk

  complex(kind=db), dimension(npldafs):: Cap_contrib, Sum_Bragg_NR_T
  complex(kind=db), dimension(npldafs_NR):: Sum_Bragg_NR_Bulk
  complex(kind=db), dimension(nes):: Trs
  complex(kind=db), dimension(nes,nphim,npldafs):: As
  complex(kind=db), dimension(nes,nphim,npldafs):: As_bulk
  complex(kind=db), dimension(nes,nphim,npldafs,n_z_SRXRD_abs):: As_SRXRD
  complex(kind=db), dimension(nes,nphim,npldafs,n_z_SRXRD_abs_bulk):: As_SRXRD_bulk
  complex(kind=db), dimension(npldafs,n_z_SRXRD_abs):: Sum_Bragg_abs_f0
  complex(kind=db), dimension(npldafs,n_z_SRXRD_abs_bulk):: Sum_Bragg_abs_f0_bulk
  complex(kind=db), dimension(npldafs,n_z_SRXRD):: Sum_Bragg_nonabs_f
  complex(kind=db), dimension(npldafs,n_z_SRXRD_bulk):: Sum_Bragg_nonabs_f_bulk

  real(kind=db):: Abs_cap, Abs_z, Pos_z
  
  real(kind=db), dimension(n_z_SRXRD):: Abs_NR_SRXRD, Pos_z_SRXRD
  real(kind=db), dimension(n_z_SRXRD_abs):: Abs_NR_SRXRD_abs, Pos_z_SRXRD_abs
  real(kind=db), dimension(n_z_SRXRD_abs_bulk):: Abs_NR_SRXRD_abs_bulk, Pos_z_SRXRD_abs_bulk
  real(kind=db), dimension(n_z_SRXRD_bulk):: Abs_NR_SRXRD_bulk, Pos_z_SRXRD_bulk
  real(kind=db), dimension(npldafs):: L_dafs, Length_abs
  real(kind=db), dimension(nes):: Abs_z_a, Es
  real(kind=db), dimension(nes,n_z_SRXRD_abs):: Xs_SRXRD
  real(kind=db), dimension(nes,n_z_SRXRD_abs_bulk):: Xs_SRXRD_bulk

! Surface

  if( Abs_cap > 0._db ) then
    Cap_contrib(:) = Sum_Bragg_NR_T(:) * exp( - 0.5_db * Abs_cap * Length_abs(:) ) 
    do ie = 1,nes 
      As(ie,1,:) =  Cap_contrib(:)  
    end do
  endif
 
  do i = 1,n_z_SRXRD + n_z_SRXRD_abs

    Abs_z_a(:) = 0._db
    Abs_z = Abs_cap

    if( i > n_z_SRXRD ) then
      Pos_z = Pos_z_SRXRD_abs(i-n_z_SRXRD)
    else
      Pos_z = Pos_z_SRXRD(i)
    endif 
    do j = 1,n_z_SRXRD
      if( Pos_z_SRXRD(j) < Pos_z - eps10 ) then
        cycle
      elseif( abs( Pos_z_SRXRD(j) - Pos_z ) < eps10 ) then
        Abs_z = Abs_z + 0.5_db * Abs_NR_SRXRD(j)
      else
        Abs_z = Abs_z + Abs_NR_SRXRD(j)
      endif
    end do

    do j = 1,n_z_SRXRD_abs
      if( Pos_z_SRXRD_abs(j) < Pos_z - eps10 ) then
        cycle
      elseif( abs( Pos_z_SRXRD_abs(j) - Pos_z ) < eps10 ) then
        Abs_z = Abs_z + 0.5_db * Abs_NR_SRXRD_abs(j)
        Abs_z_a(:) = Abs_z_a(:) + 0.5_db * Xs_SRXRD(:,j)  
      else
        Abs_z = Abs_z + Abs_NR_SRXRD_abs(j)
        Abs_z_a(:) = Abs_z_a(:) + Xs_SRXRD(:,j)  
      endif
    end do

    if( icheck > 2 ) then
      write(3,110) i, Abs_z
      write(3,'(/A)' ) '   Energy        Abs'
      do ie = 1,nes
        write(3,'(f10.5,1pe13.5)') Es(ie)*rydb, Abs_z_a(ie) 
      end do
    endif

    Abs_z_a(:) = Abs_z + Abs_z_a(:)
      
    if( i > n_z_SRXRD ) then

      j = i - n_z_SRXRD

      do ipl = 1,npldafs
        As(:,1,ipl) = As(:,1,ipl) + ( Sum_Bragg_abs_f0(ipl,j) + As_SRXRD(:,1,ipl,j) ) * exp( - Abs_z_a(:) * Length_abs(ipl) )  
      end do
     
    else
  
      do ipl = 1,npldafs
        As(:,1,ipl) = As(:,1,ipl) + Sum_Bragg_nonabs_f(ipl,i) * exp( - Abs_z_a(:) * Length_abs(ipl) )  
      end do
     
    endif

    if( icheck > 2 ) then
      write(3,*)
      if( i > n_z_SRXRD ) then
        write(3,120) 'Sum_Bragg_abs_f0 =', Sum_Bragg_abs_f0(1:min(npldafs,5),j) 
      else
        write(3,120) 'Sum_Bragg_nabs_f =', Sum_Bragg_nonabs_f(1:min(npldafs,5),i) 
      endif
      write(3,120)   'Length_abs       =', Length_abs(1:min(npldafs,5))

      write(3,'(/A)')   '   Energy       As'
      do ie = 1,nes
        write(3,'(f10.5,1p,50(1x,2e13.5))') Es(ie)*rydb, As(ie,1,1:min(npldafs,5)) 
      end do
    endif
        
  end do 

! Bulk
  As_bulk(:,:,:) = ( 0._db, 0._db )
  
  if( .not. Non_resonant_bulk ) then
  
    do i = 1,n_z_SRXRD_bulk + n_z_SRXRD_abs_bulk

      Abs_z_a(:) = 0._db
      Abs_z = 0._db
      
      if( i > n_z_SRXRD_bulk ) then
        Pos_z = Pos_z_SRXRD_abs_bulk(i-n_z_SRXRD_bulk)
      else
        Pos_z = Pos_z_SRXRD_bulk(i)
      endif 
      do j = 1,n_z_SRXRD_bulk
        if( Pos_z_SRXRD_bulk(j) < Pos_z - eps10 ) then
          cycle
        elseif( abs( Pos_z_SRXRD_bulk(j) - Pos_z ) < eps10 ) then
          Abs_z = Abs_z + 0.5_db * Abs_NR_SRXRD_bulk(j)
        else
          Abs_z = Abs_z + Abs_NR_SRXRD_bulk(j)
        endif
      end do

      do j = 1,n_z_SRXRD_abs_bulk
        if( Pos_z_SRXRD_abs_bulk(j) < Pos_z - eps10 ) then
          cycle
        elseif( abs( Pos_z_SRXRD_abs_bulk(j) - Pos_z ) < eps10 ) then
          Abs_z = Abs_z + 0.5_db * Abs_NR_SRXRD_abs_bulk(j)
          Abs_z_a(:) = Abs_z_a(:) + 0.5_db * Xs_SRXRD_bulk(:,j)  
        else
          Abs_z = Abs_z + Abs_NR_SRXRD_abs_bulk(j)
          Abs_z_a(:) = Abs_z_a(:) + Xs_SRXRD_bulk(:,j)  
        endif
      end do

      if( icheck > 2 ) then
        write(3,110) i, Abs_z
        write(3,'(/A)' ) '   Energy        Abs    Bulk'
        do ie = 1,nes
          write(3,'(f10.5,1pe13.5)') Es(ie)*rydb, Abs_z_a(ie) 
        end do
      endif

      Abs_z_a(:) = Abs_z + Abs_z_a(:)

      if( i > n_z_SRXRD_bulk ) then

        j = i - n_z_SRXRD_bulk

        do ipl = 1,npldafs
          As_bulk(:,1,ipl) = As_bulk(:,1,ipl) + ( Sum_Bragg_abs_f0_bulk(ipl,j) + As_SRXRD_bulk(:,1,ipl,j) ) &
                                              * exp( - Abs_z_a(:) * Length_abs(ipl) )  
        end do
       
      else
    
        do ipl = 1,npldafs
          As_bulk(:,1,ipl) = As_bulk(:,1,ipl) + Sum_Bragg_nonabs_f_bulk(ipl,i) * exp( - Abs_z_a(:) * Length_abs(ipl) )  
        end do
       
      endif

      if( icheck > 2 ) then
        write(3,*)
        if( i > n_z_SRXRD ) then
          write(3,125) 'Sum_Bragg_abs_f0_bulk =', Sum_Bragg_abs_f0_bulk(1:min(npldafs,5),j) 
        else
          write(3,125) 'Sum_Bragg_nabs_f_bulk =', Sum_Bragg_nonabs_f_bulk(1:min(npldafs,5),i) 
        endif
        write(3,125)   'Length_abs            =', Length_abs(1:min(npldafs,5))

        write(3,'(/A)')   '   Energy       As_bulk'
        do ie = 1,nes
          write(3,'(f10.5,1p,50(1x,2e13.5))') Es(ie)*rydb, As_bulk(ie,1,1:min(npldafs,5)) 
        end do
      endif
                  
    end do
    
! Total absorption from bulk

    Abs_z = sum( Abs_NR_SRXRD_bulk(:) ) + sum( Abs_NR_SRXRD_abs_bulk(:) )
    Abs_z_a(:) = Abs_z
    do i = 1, n_z_SRXRD_abs_bulk
      Abs_z_a(:) = Abs_z_a(:) + Xs_SRXRD_bulk(:,i)
    end do

    if( icheck > 2 ) then
      write(3,140) Abs_z
      write(3,'(/A)' ) '   Energy        total Abs    Bulk'
      do ie = 1,nes
        write(3,'(f10.5,1pe13.5)') Es(ie)*rydb, Abs_z_a(ie) 
      end do
    endif

    do ipl = 1,npldafs 
      Trs(:) = 1 / ( 1 - exp( - Abs_z_a(:) * Length_abs(ipl) - img * 2 * pi * L_dafs(ipl) ) )
      As_bulk(:,1,ipl) = Trs(:) * As_bulk(:,1,ipl)

      if( ipl > 5 .or. icheck < 3 ) cycle 
      write(3,'(/a69,i3)' ) '   Energy             Trs                      As_bulk          ipl =', ipl
      do ie = 1,nes
        write(3,'(f10.5,1p,2(1x,2e13.5))') Es(ie)*rydb, Trs(ie),  As_bulk(ie,1,ipl) 
      end do

    end do

  endif
  
! Total absorption from surface

  Abs_z = Abs_cap + sum( Abs_NR_SRXRD(:) ) + sum( Abs_NR_SRXRD_abs(:) )
  Abs_z_a(:) = Abs_z
  do i = 1, n_z_SRXRD_abs
    Abs_z_a(:) = Abs_z_a(:) + Xs_SRXRD(:,i)
  end do

  if( icheck > 2 ) then
    write(3,150) Abs_z
    write(3,'(/A)' ) '   Energy        total Abs    surface'
    do ie = 1,nes
      write(3,'(f10.5,1pe13.5)') Es(ie)*rydb, Abs_z_a(ie) 
    end do
  endif

  if( Non_resonant_bulk ) then
    do ipl = 1,npldafs 
      As_bulk(:,1,ipl) = Sum_Bragg_NR_Bulk(ipl) * exp( - Abs_z_a(:) * Length_abs(ipl) )
    end do
  else
    do ipl = 1,npldafs 
      As_bulk(:,1,ipl) = As_bulk(:,1,ipl) * exp( - Abs_z_a(:) * Length_abs(ipl) )
    end do
  endif

  As(:,:,:) = As(:,:,:) + As_bulk(:,:,:)
    
  if( icheck > 1 ) then
    write(3,'(/A)')   '   Energy       As_bulk   with truncature and self-absorption from the surface'
    do ie = 1,nes
      write(3,'(f10.5,1p,50(1x,2e13.5))') Es(ie)*rydb, As_bulk(ie,1,1:min(npldafs,5)) 
    end do

    write(3,'(/A)')   '   Energy       As        with truncature'
    do ie = 1,nes
      write(3,'(f10.5,1p,50(1x,2e13.5))') Es(ie)*rydb, As(ie,1,1:min(npldafs,5)) 
    end do
  endif
    
  return
 110 format(/' Atom =',i3,', Abs_z =',1p,e13.5,' Mbarn')
 120 format(2x,a18,1p,50e13.5) 
 125 format(2x,a23,1p,50e13.5) 
 140 format(/' Total absorption from bulk, Abs_z =',1p,e13.5,' Mbarn')
 150 format(/' Total absorption from surface, Abs_z =',1p,e13.5,' Mbarn')
end

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

! Convolution by a gaussian

subroutine Main_gaussian(File_in,itape,File_out)

  use declarations
  implicit none

  integer:: eof, i, ier, ipr, istat, itape, Length, ligne, n, n_text, ncol, nnombre, nx

  character(len=9):: keyword
  character(len=132):: identmot, mot, Traduction
  character(len=1320):: Title
  character(len=Length_name):: File_in, File_out

  real(kind=db):: a, sigma
  real(kind=db), dimension(:), allocatable:: x
  real(kind=db), dimension(:,:), allocatable:: y, z
 
  Rewind(itape)
  
  i = 0
  do ligne = 1,1000000

    n = nnombre(itape,132)
    read(itape,'(A)',iostat=eof) mot
    if( eof /= 0 ) exit

    keyword = identmot(mot,9)
    keyword = Traduction(keyword)
    
    if( keyword == 'end' ) exit

    if( keyword == 'conv_gaus' ) then
      n = nnombre(itape,132)
      read(itape,*,iostat=ier) sigma
      if( ier > 0 ) call write_err_form(itape,keyword)
      i = i + 1
    endif
    
  end do
  
  if( i /= 1 .or. File_in == ' ' ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,'(A)') ' Both "File_in" and "Conv_gaus" keywords must be in the indata file !'
    end do
    stop
  endif

  if( sigma < -eps10 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,110) sigma
    end do
    stop
  endif
  
! Dimensions

  Open(20, file = File_in, status='old', iostat=istat )
      
  do i = 1,10000
    n = nnombre(20,10000)
    if( n == 0 ) then
      read(20,'(A)') mot
    else
      exit
    endif
  end do
  n_text = i - 1
  
  ncol = n - 1

  do i = 1,1000000
    read(20,*,iostat=istat) a
    if( istat /= 0 ) exit
  end do
  nx = i - 1

  allocate( x(nx) ); allocate( y(nx,ncol) ); allocate( z(nx,ncol) )

! Reading

  Rewind(20)

  do i = 1, n_text - 1
    read(20,*)
  end do
  if( n_text > 0 ) read(20,'(A)') Title

  do i = 1,nx
    read(20,*,iostat=istat) x(i), y(i,:)
    if( istat /= 0 ) then
      call write_error
      do ipr = 6,9,3
        write(ipr,120) File_in
      end do
      stop
    endif
  end do

  Close(20) 
  
  call Conv_gaussian(ncol,nx,sigma,x,y,z)

! Writing
  
  Length = len_trim( File_out )
  if( Length > 4 ) then
    if( File_out(Length-3:Length-3) /= '.' ) File_out(Length+1:Length+4) = '.txt'
  endif

  Open(20, file = File_out )

  if( n_text > 0 ) then
    Length = len_trim(Title)
    write(20,'(4x,A)') Title(1:Length)
  endif
  
  do i = 1,nx
    write(20,130) x(i), z(i,:)
  end do
  
  Close(20)
    
  deallocate( x, y, z )
    
  return
  110 format(//' In the gaussian broadening, the width =',1p,e11.3,' is negative, what is not possible !'//)
  120 format(//' In the file:',/A,/'  the columns must have the same number of elements',//)
  130 format(f15.7,1p,1000e13.5) 
end

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

subroutine Conv_gaussian(ncol,nx,sigma,x,y,z)

  use declarations
  implicit none

  integer:: i, j, ncol, nx
  
  real(kind=db):: alfa, f, g, g_sum, mu, sigma
  real(kind=db), dimension(nx):: dx, x
  real(kind=db), dimension(nx,ncol):: y, z

  if( abs( sigma ) < eps10 ) then
    z(:,:) = y(:,:)
    return
  endif
   
  alfa = 0.5_db / sigma**2
  f = 1 / ( sigma * sqrt( 2._db * pi ) )
  
  z(:,:) = 0._db

  dx(1) = ( x(2) - x(1) ) / 2
  do i = 2,nx-1
    dx(i) = ( x(i+1) - x(i-1) ) / 2
  end do
  dx(nx) = ( x(nx) - x(nx-1) ) / 2
    
  do i = 1,nx
  
    mu = x(i)
    
    g_sum = 0
    
    do j = 1,nx
    
      g = f * exp( - alfa * ( x(j) - mu )**2 ) * dx(j)
      g_sum = g_sum + g
      
      z(i,:) = z(i,:) + g * y(j,:)
      
    end do
    
    z(i,:) = z(i,:) / g_sum
    
  end do

  return
end

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

! Writing

Subroutine Writing_output_conv(Ampl_Phase,Analyzer,Angle,As,Bormann,Convolution_out_L,Cor_abs,Couple_co,Dafs_bio, &
      Dafs_exp_type,Double_cor,E_cut,Es,Energ_tr,Eseuil,File_scan_out,Fit_cal,fprim,Full_self_abs, &
      Gamma_h,Gamma_max,Header,icheck,Icirc,Icirccor,Icircdcor,Icor,Idcor,iscratchconv,Mu_mat_comp, &
      Mus,Mus_mat,n_bir,n_col,nes,n_couple_weight_co,n_energ_tr,n_mat_pol,n_stokes,nfich,ninitm, &
      nom_col,np_stokes,nphi,nphim,npldafs,npldafs_b,npldafs_t,nxan,Scan_true,Shift_init,Stokes,Tenseur,Transpose_file, &
      Weight_co,Xs)

  use declarations
  implicit none
  
  integer:: Dafs_exp_type, i, i_hk, i1, icheck, ie, ipas, ip, ipl, iscratchconv, jpl, j, kpl, Length, n, n_bir, n_col, &
    n_couple_weight_co, n_energ_tr, n_energ_tr_t, n_header_info, n_mat_pol, n_signal, n_stokes, nes, nfich, ninitm, np_stokes, &
    nphim, npldafs, npldafs_b, npldafs_t, nw, nxan

  integer, dimension(npldafs):: nphi
  integer, dimension(n_couple_weight_co,2):: Couple_co

  character(len=Length_word):: nomab   
  character(len=Length_name):: Convolution_out_L, File_scan_out
  character(len=Length_word), dimension(n_col):: nom_col

  character(len=Length_name), dimension(:), allocatable:: Header_info
  
  complex(kind=db), dimension(nes,nphim,npldafs):: As
  complex(kind=db), dimension(nes,4,n_mat_pol):: Mu_mat_comp
  complex(kind=db), dimension(nes,nphim,npldafs_b,2):: Mus
  
  logical:: Ampl_Phase, Analyzer, Bormann, Cor_abs, Dafs_bio, Double_cor, First_E, Fit_cal, fprim, Full_self_abs, Header, &
            New_file, Scan_true, Stokes, Tenseur, Transpose_file

  real(kind=db):: a, E_cut, Gamma_h, Gamma_max, Im_pi, Im_sig, Ip_pi, Ip_sig, Weight_co
  
  real(kind=db), dimension(nphim):: Angle
  real(kind=db), dimension(nfich):: Eseuil
  real(kind=db), dimension(n_col):: Tens
  real(kind=db), dimension(nes):: Es
  real(kind=db), dimension(n_energ_tr):: Energ_tr
  real(kind=db), dimension(nes,nxan):: Xs
  real(kind=db), dimension(ninitm,nfich):: Shift_init
  real(kind=db), dimension(nes,n_stokes,n_bir,n_mat_pol):: Mus_mat
  real(kind=db), dimension(nes,nphim,npldafs):: Icor, Idcor
  real(kind=db), dimension(nes,nphim,np_stokes):: Icirc, Icirccor, Icircdcor

  real(kind=db), dimension(:), allocatable:: Energ_tr_t   
  real(kind=db), dimension(:,:), allocatable:: Signal 

  if( Double_cor ) then
    n_signal = 3 * nes
    i1 = n_col - npldafs_t - 4*npldafs + 1
    ipas = 5
  elseif( Cor_abs ) then
    n_signal = 2*nes
    i1 = n_col - npldafs_t - 3*npldafs + 1
    ipas = 4
  else
    n_signal = nes
    i1 = n_col - npldafs_t + 1
    ipas = 1
  endif
  
  if( Transpose_file ) allocate( Signal(npldafs_t,n_signal) )
  
  First_E = .true.
  
  do ie = 1,nes

    do ipl = 1,nxan
      Tens(ipl) = Xs(ie,ipl)
    end do

    jpl = nxan

    do ipl = 1,n_mat_pol
      jpl = jpl + 1
      Tens(jpl) = real( Mu_mat_comp(ie,1,ipl), db )   ! sigma-sigma
      jpl = jpl + 1
      Tens(jpl) = real( Mu_mat_comp(ie,4,ipl), db )   ! pi-pi
      do i = 1,n_stokes
        do j = 1,n_bir
          jpl = jpl + 1
          Tens(jpl) = Mus_mat(ie,i,j,ipl)
        end do
      end do
    end do

    if( Dafs_bio ) then
      if( icheck > 1 ) then
        nw = 5
      else
        nw = 1
      endif

! There are 4 polarization states corresponding to Q and -Q, both with sigma out and pi out.
! The intensity is the difference between I(Q) = I(Q)_sigma + I(Q)_pi, and I(-Q)
      do ipl = 1,npldafs
        Ip_sig = real( As(ie,1,ipl), db )**2 + aimag( As(ie,1,ipl) )**2
        Ip_pi =  real( As(ie,2,ipl), db )**2 + aimag( As(ie,2,ipl) )**2
        Im_sig = real( As(ie,3,ipl), db )**2 + aimag( As(ie,3,ipl) )**2
        Im_pi =  real( As(ie,4,ipl), db )**2 + aimag( As(ie,4,ipl) )**2
        do i = 1,nw
          jpl = jpl + 1
          select case(i)
            case(1)
              select case( Dafs_exp_type )
                case(1)
                  Tens(jpl) = Ip_sig + Ip_pi
                case(2)
                  Tens(jpl) = sqrt( Ip_sig + Ip_pi )
                case(3)
                  Tens(jpl) = Ip_sig + Ip_pi - ( Im_sig + Im_pi )
                case default
                  Tens(jpl) = sqrt( Ip_sig + Ip_pi ) - sqrt( Im_sig + Im_pi )
              end select
            case(2)
              Tens(jpl) = Ip_sig
            case(3)
              Tens(jpl) = Ip_pi
            case(4)
              Tens(jpl) = Im_sig
            case(5)
              Tens(jpl) = Im_pi
          end select
        end do
      end do

    else

      i_hk = 1
      kpl = 0
 
      do ipl = 1,npldafs
        if( .not. ( Tenseur .or. Bormann ) ) then
          if( Analyzer .or. mod(ipl,2) == 1 ) jpl = jpl + 1
          if( Analyzer ) then
            Tens(jpl) = abs( As(ie,1,ipl) )**2
            if( Transpose_file ) Signal(ipl,ie) = Tens(jpl)
          elseif( mod(ipl,2) == 1 ) then
            Tens(jpl) = abs( As(ie,1,ipl) )**2 + abs( As(ie,1,ipl+1) )**2
            kpl = kpl + 1
            if( Transpose_file ) Signal(kpl,ie) = Tens(jpl)
          endif
        endif
        if( Ampl_Phase ) then
          jpl = jpl + 1
          Tens(jpl) = abs( As(ie,1,ipl) )
          jpl = jpl + 1
          Tens(jpl) = atan2( aimag( As(ie,1,ipl) ), real( As(ie,1,ipl) ) )
        endif
        if( fprim .or. Tenseur .or. Bormann ) then
          jpl = jpl + 1
          Tens(jpl) = real( As(ie,1,ipl), db )
          jpl = jpl + 1
          Tens(jpl) = aimag( As(ie,1,ipl) )
          if( Tenseur .or. Bormann ) cycle
        endif
        if( Cor_abs ) then
          if( Analyzer .or. mod(ipl,2) == 1 ) jpl = jpl + 1
          if( Analyzer ) then
            Tens(jpl) = Icor(ie,1,ipl)
          elseif( mod(ipl,2) == 1 ) then
            Tens(jpl) = Icor(ie,1,ipl) + Icor(ie,1,ipl+1)
          endif
        endif
        if( Double_cor ) then
          if( Analyzer .or. mod(ipl,2) == 1 ) jpl = jpl + 1
          if( Analyzer ) then
            Tens(jpl) = Idcor(ie,1,ipl)
            if( Transpose_file ) Signal(ipl,2*nes+ie) = Tens(jpl)
          elseif( mod(ipl,2) == 1 ) then
            Tens(jpl) = Idcor(ie,1,ipl) + Idcor(ie,1,ipl+1)
          endif
        endif
        if( Cor_abs ) then
          if( Analyzer .or. mod(ipl,2) == 1 ) then
            do i = 1,2
              jpl = jpl + 1
              Tens(jpl) = real( Mus(ie,1,ipl,i), db )
              if( Full_self_abs ) then
                jpl = jpl + 1
                Tens(jpl) = aimag( Mus(ie,1,ipl,i) )
              endif
            end do
          endif
        endif
        if( Stokes .and. mod(ipl,4) == 0 ) then
          j = ( ipl/4 - 1 )*n_stokes
          do i = 1,n_stokes
            if( Analyzer .or. mod(i,2) == 1 ) jpl = jpl + 1
            if( Analyzer ) then
              Tens(jpl) = Icirc(ie,1,j+i)
            elseif( mod(i,2) == 1 ) then
              Tens(jpl) = Icirc(ie,1,j+i) + Icirc(ie,1,j+i+1)
            endif
            if( Cor_abs ) then
              if( Analyzer .or. mod(i,2) == 1 ) jpl = jpl + 1
              if( Analyzer ) then
                Tens(jpl) = Icirccor(ie,1,j+i)
              elseif( mod(j+i,2) == 1 ) then
                Tens(jpl) = Icirccor(ie,1,j+i) + Icirccor(ie,1,j+i+1)
               endif
            endif
            if( Double_cor ) then
              if( Analyzer .or. mod(i,2) == 1 ) jpl = jpl + 1
              if( Analyzer ) then
                Tens(jpl) = Icircdcor(ie,1,j+i)
              elseif( mod(j+i,2) == 1 ) then
                Tens(jpl) = Icircdcor(ie,1,j+i) + Icircdcor(ie,1,j+i+1)
              endif
            endif
          end do
        endif
      end do

    endif

! Weighted average between columns
    do i = 1,n_couple_weight_co
      ipl = Couple_co(i,1)
      jpl = Couple_co(i,2)
      a = Tens(ipl)
      Tens(ipl) = Weight_co * Tens(ipl) + ( 1 - Weight_co ) * Tens(jpl)  
      Tens(jpl) = ( 1 - Weight_co ) * a + Weight_co * Tens(jpl)  
    end do 
    
    if( Fit_cal ) then
! Command the writing in a temporary file
      n = -iscratchconv
    else
      n = 1
    endif
    New_file = First_E
    if( Header .and. First_E .and. .not. Fit_cal ) then
      if( nfich == 1 .and. ninitm == 1 ) then
        n_header_info = 1
      else
        n_header_info = 2
      endif
      allocate( Header_info(n_header_info) )
      call Construct_header(E_cut,Eseuil(1),Gamma_h,Gamma_max,Header_info,n_header_info,ninitm,nfich,Shift_init)

      call Write_Header(Convolution_out_L,Header_info,n_header_info)
      New_file = .false.
      deallocate( Header_info )
    endif
 
    call Write_out_stand(Es(ie),Convolution_out_L,First_E,n_col,jpl,New_file,nom_col,n,Tens)
     
    First_E = .false.
  end do
  
  if( Transpose_file ) then
    if( n_energ_tr == 0 ) then
      n_energ_tr_t = nes
    else
      n_energ_tr_t = n_energ_tr
    endif
    allocate( Energ_tr_t(n_energ_tr_t) )
    if( n_energ_tr == 0 ) then
      Energ_tr_t(:) = Es(:)
    else
      Energ_tr_t(:) = Energ_tr(:)
    endif
    call Write_transpose(Convolution_out_L,Energ_tr_t,Es,n_col,n_energ_tr_t,n_signal,nes,nom_col,npldafs_t,Signal)
    deallocate( Signal )
  endif

  if( Scan_true .and. .not. Dafs_bio ) then

    Open(7, file = File_scan_out)

! One skips the '_0'
    do i = 1,n_col
      nomab = adjustl( nom_col(i) )
      if( nomab(1:1) /= 'I' ) cycle
      Length = len_trim(nomab)
      if( nomab(Length-1:Length) == '_0' ) nomab(Length-1:Length) = '  '
      nom_col(i) = nomab
    end do

    do ie = 1,nes
      Write(7,*)
      Write(7,200) Es(ie) * rydb
      jpl = nxan

      do ipl = 1,npldafs

        if( Analyzer .or. mod(ipl,2) == 1 ) then
          jpl = jpl + 1
          Write(7,*) nom_col(jpl)
        endif
        do ip = 1,nphi(ipl)
          if( Analyzer ) then
            Write(7,210) angle(ip), real( As(ie,ip,ipl), db )**2 + aimag( As(ie,ip,ipl) )**2
          elseif( mod(ipl,2) == 1 ) then
            Write(7,210) angle(ip), real( As(ie,ip,ipl), db )**2 + aimag( As(ie,ip,ipl) )**2 &
                                  + real( As(ie,ip,ipl+1), db )**2 + aimag( As(ie,ip,ipl+1) )**2
          endif
        end do

        if( Ampl_phase ) jpl = jpl + 2
        if( fprim ) jpl = jpl + 2

        if( Cor_abs ) then
          if( Analyzer .or. mod(ipl,2) == 1 ) then
            jpl = jpl + 1
            Write(7,*) nom_col(jpl)
          endif
          do ip = 1,nphi(ipl)
            if( Analyzer ) then
              Write(7,210) angle(ip), Icor(ie,ip,ipl)
            elseif( mod(ipl,2) == 1 ) then
              Write(7,210) angle(ip), Icor(ie,ip,ipl) + Icor(ie,ip,ipl+1)
            endif  
          end do
          if( Double_cor  ) then
            if( Analyzer .or. mod(ipl,2) == 1 ) then
              jpl = jpl + 1
              Write(7,*) nom_col(jpl)
            endif
            do ip = 1,nphi(ipl)
              if( Analyzer ) then
                Write(7,210) angle(ip), Idcor(ie,ip,ipl)
              elseif( mod(ipl,2) == 1 ) then
                Write(7,210) angle(ip), Idcor(ie,ip,ipl) + Idcor(ie,ip,ipl+1)
              endif  
            end do
          endif
          if( Analyzer .or. mod(ipl,2) == 1 ) then
            jpl = jpl + 2  ! To skip mu
            if( Full_self_abs ) jpl = jpl + 2  ! To skip mu
          endif
        endif

        if( Stokes .and. mod(ipl,4) == 0 ) then
          j = ( ipl/4 - 1 ) * n_stokes
          do i = 1,n_stokes
            if( Analyzer .or. mod(i,2) == 1 ) then
              jpl = jpl + 1
              Write(7,*) nom_col(jpl)
            endif
            do ip = 1,nphi(ipl)
              if( Analyzer ) then
                Write(7,210) angle(ip), Icirc(ie,ip,j+i)
              elseif( mod(i,2) == 1 ) then
                Write(7,210) angle(ip), Icirc(ie,ip,j+i) + Icirc(ie,ip,j+i+1)
              endif  
            end do
            if( Analyzer .or. mod(i,2) == 1 ) then
              jpl = jpl + 1
              Write(7,*) nom_col(jpl)
            endif
            if( Analyzer ) then
              do ip = 1,nphi(ipl)
                Write(7,210) angle(ip), Icirccor(ie,ip,j+i)
              end do
            elseif( mod(i,2) == 1 ) then
              do ip = 1,nphi(ipl)
                Write(7,210) angle(ip), Icirccor(ie,ip,j+i) + Icirccor(ie,ip,j+i+1)
              end do
            endif
            if( Double_cor ) then
              if( Analyzer .or. mod(i,2) == 1 ) then
                jpl = jpl + 1
                Write(7,*) nom_col(jpl)
              endif
              do ip = 1,nphi(ipl)
                if( Analyzer ) then
                  Write(7,210) angle(ip), Icircdcor(ie,ip,j+i)
                elseif( mod(i,2) == 1 ) then
                  Write(7,210) angle(ip), Icircdcor(ie,ip,j+i) + Icircdcor(ie,ip,j+i+1)
                endif  
              end do
            endif
          end do
        endif

      end do
    end do
    
    Close(7)

  endif

  return
  200 format(f10.3,1p,240e13.5)
  210 format(f7.1,1p,3e13.5)
end

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

subroutine Construct_header(E_cut,Eseuil,Gamma_h,Gamma_max,Header_info,n_header_info,ninitm,nfich,Shift_init)

  use declarations
  implicit none
  
  integer:: i, i_header, ifich, init, L, n_header_info, nfich, ni, ninitm
  
  character(len=13):: Word13
  character(len=Length_name):: Word
  character(len=Length_name), dimension(n_header_info):: Header_info

  real(kind=db):: E, E_cut, Gamma_h, Gamma_max
  real(kind=db), dimension(nfich):: Eseuil
  real(kind=db), dimension(ninitm,nfich):: Shift_init
  
  Header_info(:) = ' '
  
  do i_header = 1,n_header_info

    Word = ' '
    Word(1:1) = '#'
    L = 3 - 14
    if( i_header == 1 ) then
      ni = 4
    else
      ni = ninitm * nfich
    endif
    
    do i = 1,ni
      L = L + 14
      if( i_header == 1 ) then
        select case(i)
          case(1)
            E = Eseuil(1) * Rydb
          case(2)
            E = Gamma_h * Rydb
          case(3)
            E = Gamma_max * Rydb
          case(4)
            E = E_cut * Rydb
        end select
      else
        ifich = 1 + ( i - 1 ) / ninitm
        init = i - ( ifich - 1 ) * ninitm
        E = Shift_init(init,ifich) * Rydb
      endif
      if( abs( nint( 10*E ) - 10*E ) < eps10 ) then
        write(Word13,'(f13.1)') E
      elseif( abs( nint( 100*E ) - 100*E ) < eps10 ) then
        write(Word13,'(f13.2)') E
      elseif( abs( nint( 1000*E ) - 1000*E ) < eps10 .or. E > 999.99_db ) then
        write(Word13,'(f13.3)') E
      elseif( abs( nint( 1000*E ) - 1000*E ) < eps10 ) then
        write(Word13,'(f13.4)') E
      else
        write(Word13,'(f13.5)') E
      endif
      Word13 = adjustr(Word13)
      Word(L:L+12) = Word13
    
    end do

    L = L+15
    if( i_header == 1 ) then
      Word(L:L+37) = '= E_edge, Gamma_hole, Gamma_max, E_cut'
    else
      Word(L:L+7) = '= Shifts'
    endif  
    Header_info(i_header) = Word
   
  end do

  return
end

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

subroutine Write_transpose(Convolution_out,Energ_tr,Es,n_col,n_energ_tr,n_signal,nes,nom_col,npldafs,Signal)

  use declarations
  implicit none
  
  integer:: i, i_col, i_cor, i_hk, i1, i2, ie, index, index_hk, index_max, ipas, ipr, j, Length, l2, l3, le, length_hk, LL, &
            n_col, n_cor, n_hk_tr, n_energ_tr, n_signal, nes, npldafs
  integer, dimension(:), allocatable:: hk_length, n_index_hk
  integer, dimension(n_energ_tr):: index_ie

  logical:: Double_cor, Python
  
  character(len=Length_word-9):: Word_m9, Word_m9_b
  character(len=Length_word-9), dimension(:), allocatable:: hk
  character(len=Length_word):: Word, Word_b, Word_c, Word_d
  character(len=Length_word), dimension(n_col):: nom_col
  character(len=Length_word), dimension(:,:,:), allocatable:: E_string, l_value
  character(len=Length_word), dimension(:,:,:,:), allocatable:: Signal_out
  character(len=Length_name):: Convolution_out, Convolution_tr

  real(kind=db):: E
  real(kind=db), dimension(nes):: Es
  real(kind=db), dimension(n_energ_tr):: Energ_tr, p1
  real(kind=db), dimension(npldafs,n_signal):: Signal
 
  Convolution_tr = Convolution_out
  Length = len_trim(Convolution_tr)
  Convolution_tr(Length-3:Length+3) = '_tr.txt'
  
  ipr = 4
  
  open( ipr, file = Convolution_tr )

  if( n_signal == 3*nes ) then
    Double_cor = .true.
    i1 = n_col - 5*npldafs + 1
    ipas = 5
    n_cor = 3
  elseif( n_signal == 2*nes ) then
    Double_cor = .false.
    i1 = n_col - 4*npldafs + 1
    ipas = 4
    n_cor = 2
  else
    Double_cor = .false.
    i1 = n_col - npldafs + 1
    ipas = 1
    n_cor = 1
  endif

  do ie = 1,n_energ_tr
    do i = 2,nes
      if( Es(i) > Energ_tr(ie) - eps10 ) exit
    end do
    i = min( i, nes )
    index_ie(ie) = i 
    p1(ie) = ( Energ_tr(ie) - Es(i-1) ) / ( Es(i) - Es(i-1) )
  end do

  do j = 1,2
  
    i_hk = 0
    index_hk = 0
    Word_m9_b = ' '
    do i_col = i1,n_col,ipas
      Word = ' '
      Word = adjustl( nom_col(i_col) )
      do i = 1,Length_word
        if( Word(i:i) /= '(' .and. Word(i:i) /= '_' ) cycle

        Python = Word(i:i) == '_' 

        Word_m9 = ' '
        if( Word(i+1:i+1) == '-' .and. Word(i+3:i+3) == '-' ) then
          Word_m9(1:4) = Word(i+1:i+4)
        elseif( Word(i+1:i+1) == '-' .or. Word(i+2:i+2) == '-' ) then
          Word_m9(1:3) = Word(i+1:i+3)
        else
          Word_m9(1:2) = Word(i+1:i+2)  
        endif
        length_hk = len_trim( Word_m9 )
        exit
      end do

      do i = 3,Length_word
        if( Word(i:i) == ')' .or. ( Python .and. Word(i:i) == '_' ) ) exit
      end do
      
      Le = len_trim( Word ) - i 
      LL = len_trim( Word_m9 )
      Le = min( Length_word - 9 - LL, Le ) 
      Word_m9(LL+1:LL+Le) = Word(i+1:i+Le)
 
      if( Word_m9 == Word_m9_b ) then    
        index_hk = index_hk + 1
        if( j == 2 ) n_index_hk(i_hk) = index_hk 
      else  
        i_hk = i_hk + 1
        index_hk = 1
        if( j == 2 ) then
          hk(i_hk) = Word_m9
          hk_length(i_hk) = length_hk
        endif
        Word_m9_b = Word_m9 
      endif
    end do
    
    if( j == 1 ) then
      n_hk_tr = i_hk
      allocate( n_index_hk(n_hk_tr) )
      n_index_hk(:) = 1
      allocate( hk(n_hk_tr) )
      allocate( hk_length(n_hk_tr) )
    endif
    
  end do
  
  index_max = 0  
  do i_hk = 1,n_hk_tr
    index_max = max( index_max, n_index_hk(i_hk) )
  end do
  
  allocate( l_value(index_max,n_cor,n_hk_tr) )
  allocate( Signal_out(index_max,n_energ_tr,n_cor,n_hk_tr) )
  allocate( E_string(n_energ_tr,n_cor,n_hk_tr) )
  Signal_out(:,:,:,:) = ' '
  l_value(:,:,:) = ' '

  index = 0
  index_hk = 0
  i_hk = 1

  do i_col = i1,n_col,ipas
    index = index + 1
    index_hk = index_hk + 1
    Word = ' '
    Word = nom_col(i_col)
    i1 = 0
    i2 = Length_word + 1
    do i = 1,Length_word
      if( i1 == 0 .and. ( Word(i:i) == '(' .or. Word(i:i) == '_' ) ) then
        i1 = i + hk_length(i_hk) + 1
      elseif( Word(i:i) == ')' .or. ( Python .and. Word(i:i) == '_' )  ) then
        i2 = i-1
        exit
      endif
    end do
    Word_b = ' '
    Word_b(1:i2-i1+1) = Word(i1:i2)
    l_value(index_hk,:,i_hk) = adjustr( Word_b )

    do i_cor = 1,n_cor
      do ie = 1,n_energ_tr
        j = ( i_cor - 1 ) * nes + index_ie(ie)
        Word_b = ' '
        Select case( Length_word )
          case(15)
            write(Word_b,'(1p,e15.7)') ( 1 - p1(ie) ) * Signal(index,j-1) + p1(ie) * Signal(index,j)
          case(16)
            write(Word_b,'(1p,e16.8)') ( 1 - p1(ie) ) * Signal(index,j-1) + p1(ie) * Signal(index,j)
          case(17)
            write(Word_b,'(1p,e17.9)') ( 1 - p1(ie) ) * Signal(index,j-1) + p1(ie) * Signal(index,j)
        end select
        Signal_out(index_hk,ie,i_cor,i_hk) = Word_b
      end do
    end do 
    
    if( index_hk == n_index_hk(i_hk) ) then
      index_hk = 0
      i_hk = i_hk + 1
      if( i_hk > n_hk_tr ) exit
    endif   
  end do 

  Length = 1
  Word_b = ''
  Word_b(1:1) = 'I'
  do i_cor = 1,n_cor
    if( i_cor == 2 ) then
      Word_b(2:2) = 'c'
      Length = 2
    elseif( i_cor == 3 ) then
      Word_b(2:2) = 'd'
      Length = 2
    endif
    do i_hk = 1,n_hk_tr
      Word_d = Word_b
      Word_m9 = hk(i_hk)
      le = len_trim( Word_m9 ) 
      Word_d(Length+1:Length+1) = '_'
      Word_d(Length+2:Length+1+le) = Word_m9(1:le) 
      l2 = len_trim(Word_d) + 1
      if( Python ) then
        Word_d(l2:l2) = '_'
      else
        Word_d(l2:l2) = '('
      endif
      do i = 1,n_energ_tr
        E = Energ_tr(i)*rydb
        Word_c = Word_d
        Word = ' '
        if( abs( nint( E ) - E ) < eps10 ) then
          write(Word,'(i6)') nint( E )
        elseif( abs( nint( 10*E ) - 10*E ) < eps10 ) then
          if( Length_word == 15 ) then
            write(Word,'(f15.1)') Energ_tr(i)*rydb
          elseif( Length_word == 17 ) then
            write(Word,'(f17.1)') Energ_tr(i)*rydb
          endif
        elseif( abs( nint( 100*E ) - 100*E ) < eps10 ) then
          if( Length_word == 15 ) then
            write(Word,'(f15.2)') Energ_tr(i)*rydb
          elseif( Length_word == 17 ) then
            write(Word,'(f17.2)') Energ_tr(i)*rydb
          endif
        elseif( abs( nint( 1000*E ) - 1000*E ) < eps10 ) then
           if( Length_word == 15 ) then
            write(Word,'(f15.3)') Energ_tr(i)*rydb
          elseif( Length_word == 17 ) then
            write(Word,'(f17.3)') Energ_tr(i)*rydb
          endif
        elseif( abs( nint( 1000*E ) - 1000*E ) < eps10 ) then
          if( Length_word == 15 ) then
            write(Word,'(f15.4)') Energ_tr(i)*rydb
          elseif( Length_word == 17 ) then
            write(Word,'(f17.4)') Energ_tr(i)*rydb
          endif
        endif
        Word = adjustl(Word)
        l3 = len_trim(Word)
        Word_c(l2+1:l2+l3) = Word(1:l3)
        if( l2+l3 < 14 .and. .not. Python ) Word_c(l2+l3+1:l2+l3+1) = ')'
        Word = adjustr( Word_c )
        E_string(i,i_cor,i_hk) = Word
      end do
    end do

  end do
  
  select case(Length_word)
    case(15)  
      write(ipr,110) (( '       L_'//hk(i_hk), E_string(:,i_cor,i_hk), i_cor = 1,n_cor ), i_hk = 1,n_hk_tr )
    case(17)  
      write(ipr,120) (( '       L_'//hk(i_hk), E_string(:,i_cor,i_hk), i_cor = 1,n_cor ), i_hk = 1,n_hk_tr )
  end select

  do index = 1,index_max
    select case(Length_word)
      case(15)  
        write(ipr,110) ( ( l_value(index,i_cor,i_hk), Signal_out(index,:,i_cor,i_hk), i_cor = 1,n_cor ), i_hk = 1,n_hk_tr )
      case(17)  
        write(ipr,120) ( ( l_value(index,i_cor,i_hk), Signal_out(index,:,i_cor,i_hk), i_cor = 1,n_cor ), i_hk = 1,n_hk_tr )
    end select
  end do 
  
  Close(ipr)
  
  deallocate( E_string, hk, hk_length, l_value, n_index_hk, Signal_out )

  return
  110 format(10000a15)
  120 format(10000a17)
end

