! FDMNES subroutine

! Set of subroutines for procedure optic

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

subroutine main_optic(Abs_U_iso,angxyz,Allsite,axyz,Bragg_abs,Cartesian_tensor, &
          Core_resolved,Dafs,Dafs_bio,dv0bdcF,E_cut,E_cut_imp,E_cut_man,E_Fermi,Eclie,Eneg,Energ_s,Ephot_min, &
          Eseuil,Full_potential,Full_self_abs,Header,hkl_dafs,Hubb_a,Hubb_d,icheck, &
          iabsorig,ip_max,ip0,isigpi,isymeq, &
          jseuil,ldip,Length_abs,lmax_pot,lmax_probe,lmaxabs_t,lmaxat0,lmaxfree,lmoins1,loct,lplus1,lqua, &
          ltypcal,m_hubb,Matper,Moyenne,mpinodes,mpirank,mpirank0,msymdd,msymddi,msymdq, &
          msymdqi,msymdo,msymdoi,msymoo,msymooi,msymqq,msymqqi,multi_0,Multipole,n_bulk_sup,n_multi_run, &
          n_oo,n_rout,n_tens_max,natomsym,nbseuil,ncolm,ncolr,ncolt,nenerg_s,ninit1,ninitr,nlm_pot, &
          nlm_probe,nlmamax,nomabs,nomfich,nomfich_cal_conv,nomfich_s,nphi_dafs,nphim,npldafs,nplr,nplrm, &
          nr,nseuil,nspin,nspino,nspinp,numat_abs,nxanout,pdp,phdt,pol,poldafse,poldafss, &
          r,Relativiste,Renorm,Rmtg,Rmtsd,rot_atom_abs,Rot_cluster,Self_abs,Spherical_signal, &
          Spherical_tensor,Spinorbite,Sum_Bragg_abs_f0,Sum_Bragg_nonabs_f,Sum_Bragg_NR_T,Taull_tdd,Taux_eq, &
          Time_rout,V_intmax,V_hubb,V0muf,Vcato,vec,vecdafse,vecdafss,Vhbdc,Volume_maille,VxcbdcF,Vxcato,Ylm_comp)

  use declarations
  implicit none

  integer, parameter:: n_z_SRXRD = 0
  integer, parameter:: n_z_SRXRD_abs = 0
  integer, parameter:: n_z_SRXRD_max_abs = 0
  integer, parameter:: n_max = 1

  integer:: iabsorig, icheck_s, ie, ie_computer, ie_q, ie_s, ie_t, ip_max, ip0, &
    isp, iss1, iss2, je, jseuil, l0_nrixs, lm1, lm2, lmax, lmax_pot, &
    lmax_probe, lmaxabs_t, m_hubb, mpinodes, mpirank, mpirank0, &
    lmax_nrixs, lmaxat0, multi_0, n_abs_rgh, n_bulk_sup, n_Ec, n_multi_run, n_oo, n_rout, &
    n_tens_max, n_V, natomsym, nbseuil, ncolm, ncolr, ncolt, &
    nenerg, nenerg_s, nge, ninit1, ninitr, nlm_pot, nlm_probe, nlm_p_fp, &
    nlmamax, nphim, npldafs, nplr, nplrm, nq_nrixs, nr, ns_dipmag, nseuil, nspin, &
    nspino, nspinp, numat_abs, nxanout

  integer, dimension(30):: icheck
  integer, dimension(natomsym):: isymeq
  integer, dimension(npldafs):: nphi_dafs
  integer, dimension(2,npldafs):: isigpi
  integer, dimension(3):: ldip
  integer, dimension(3,3):: lqua, msymdd, msymddi
  integer, dimension(3,3,3):: loct, msymdq, msymdqi
  integer, dimension(3,3,3,3)::  msymdo, msymdoi, msymqq, msymqqi
  integer, dimension(3,n_oo,3,n_oo):: msymoo, msymooi
  integer, dimension(n_z_SRXRD_abs):: n_zc_SRXRD_abs
  integer, dimension(n_z_SRXRD_max_abs,n_z_SRXRD_abs):: index_z_SRXRD_abs

  character(len=Length_name):: nomfich, nomfich_s
  character(len=5), dimension(nplrm):: ltypcal
  character(len=length_word), dimension(ncolm):: nomabs
  character(len=Length_name), dimension(n_multi_run+n_bulk_sup):: nomfich_cal_conv

  complex(kind=db):: f_avantseuil
  complex(kind=db), dimension(3,nplrm):: pol
  complex(kind=db), dimension(natomsym,npldafs):: Bragg_abs
  complex(kind=db), dimension(npldafs,nphim):: Sum_Bragg_NR_Bulk, Sum_Bragg_NR_T
  complex(kind=db), dimension(npldafs,nphim,n_max):: phdt
  complex(kind=db), dimension(npldafs,n_max):: Sum_Bragg_abs_f0
  complex(kind=db), dimension(npldafs,n_z_SRXRD):: Sum_Bragg_nonabs_f
  complex(kind=db), dimension(3,npldafs,nphim):: poldafse, poldafss
  complex(kind=db), dimension(3,3,ninitr,0:mpinodes-1):: T_E1E1, T_E1E1_m, T_E1E1_t, T_E1E1_m_t, T_E1SP, T_E1SP_m, T_E1SP_t, &
                     T_E1SP_m_t, T_E1M1, T_E1M1_m, T_M1M1, T_M1M1_m, T_E1M1_t, T_E1M1_m_t, T_M1M1_t, T_M1M1_m_t, T_SPSP, &
                     T_SPSP_m, T_SPSP_t, T_SPSP_m_t
  complex(kind=db), dimension(3,3,3,ninitr,0:mpinodes-1):: T_E1E2, T_E1E2_m, T_E1E2_t, T_E1E2_m_t, T_E1M2, T_E1M2_m, T_E1M2_t, &
                                                           T_E1M2_m_t 
  complex(kind=db), dimension(3,3,3,3,ninitr,0:mpinodes-1):: T_E1E3, T_E1E3_m, T_E1E3_t, T_E1E3_m_t, &
                     T_E2E2, T_E2E2_m, T_E2E2_t, T_E2E2_m_t, T_E2M2, T_E2M2_m, T_E2M2_t, T_E2M2_m_t
  complex(kind=db), dimension(3,n_oo,3,n_oo,ninitr,0:mpinodes-1):: T_E3E3, T_E3E3_m, T_E3E3_t, T_E3E3_m_t
  complex(kind=db), dimension(-m_hubb:m_hubb,-m_hubb:m_hubb,nspinp,nspinp):: V_hubb
  complex(kind=db), dimension(nenerg_s,nlmamax,nspinp,nlmamax,nspinp):: Taull_tdd
  complex(kind=db), dimension(:,:,:,:,:), allocatable:: S_nrixs, Taull_abs
  complex(kind=db), dimension(:,:,:,:,:,:,:), allocatable:: Chi

  logical:: Abs_in_Bulk_roughness, Allsite, Cartesian_tensor, Core_resolved, Dafs, Dafs_bio, E_cut_man, &
    E1E1, E1E2, E1E3, E1M1, E1M2, E1SP, E2E2, E2M2, E3E3, Eneg, Energphot, FDM_comp, Final_optic, Final_tddft, First_E, &
    Full_potential, Full_self_abs, Green_int, Header, Hubb_a, Hubb_d, lmaxfree, lmoins1, lplus1, M1M1, Matper, &
    Moyenne, Non_resonant_Bulk, Relativiste, Renorm, Self_abs, Solsing, &
    Spherical_signal, Spherical_tensor, Spinorbite, SPSP, Xan_atom, Ylm_comp

  logical, dimension(n_Multipole_max):: Multipole

  real(kind=sg):: time

  real(kind=db):: Abs_cap, Abs_U_iso, delta_E, E_cut, E_cut_imp, E_cut_optic, E_Fermi, Eclie, Ecmax, Ephot_min, &
    p, Rmtg, Rmtsd, Surface_ref, V_intmax, V0muf, Vhbdc, Volume_maille

  real(kind=db), dimension(3):: angxyz, axyz
  real(kind=db), dimension(3):: tp
  real(kind=db), dimension(n_rout):: Time_rout
  real(kind=db), dimension(nspin):: dv0bdcF, VxcbdcF
  real(kind=db), dimension(nenerg_s):: Energ_s
  real(kind=db), dimension(nbseuil):: Eseuil
  real(kind=db), dimension(ninitr):: Epsii, sec_atom
  real(kind=db), dimension(nr):: r
  real(kind=db), dimension(natomsym) :: Taux_eq
  real(kind=db), dimension(3,3):: rot_atom_abs, Rot_cluster
  real(kind=db), dimension(nplrm,2):: pdp
  real(kind=db), dimension(3,nplrm):: vec
  real(kind=db), dimension(3,npldafs):: hkl_dafs
  real(kind=db), dimension(npldafs):: Length_abs
  real(kind=db), dimension(n_tens_max*ninitr,0:natomsym):: Int_tens
  real(kind=db), dimension(nr,nlm_pot):: Vcato
  real(kind=db), dimension(3,npldafs,nphim):: vecdafse, vecdafss
  real(kind=db), dimension(nr,nlm_pot,nspin) :: Vxcato
  real(kind=db), dimension(n_z_SRXRD):: Abs_SRXRD, Pos_z_SRXRD
  real(kind=db), dimension(n_z_SRXRD_abs):: Abs_SRXRD_abs, Pos_z_SRXRD_abs

  real(kind=db), dimension(:), allocatable:: Eimag, En, Energ, Enervide
  real(kind=db), dimension(:,:), allocatable:: Ecinetic, V0bdc
  real(kind=db), dimension(:,:,:,:), allocatable:: Vrato

  if( icheck(1) > 0 ) write(3,100)

  Abs_cap = 0._db
  Abs_in_Bulk_roughness = .false.
  n_abs_rgh = 0
  Epsii(:) = 0._db
  Green_int = .false.
  FDM_comp = .false.
  n_V = 1
  n_Ec = 2
  ns_dipmag = 2  ! corresponds here at the 2, occupied and non-occupied state, energies, not the spin-flip in the E1M1 transition 
  Non_resonant_Bulk = .false.
  Surface_ref = 0._db
  
  if( Hubb_a .or. Full_potential ) then
    nlm_p_fp = nlm_probe
  else
    nlm_p_fp = 1
  endif

  allocate( Taull_abs(nlm_probe,nspinp,nlm_probe,nspinp,2) )
  Taull_abs(:,:,:,:,:) = (0._db, 0._db)
  allocate( En(n_Ec) )
  allocate( V0bdc(nspin,n_V) )
  
  if( E_cut_man ) then
    E_cut_optic = E_cut_imp
  else
    E_cut_optic = E_cut
  endif

  E1E1 = Multipole(1); E1E2 = Multipole(2); E1E3 = Multipole(3)
  E1M1 = Multipole(4); E1M2 = Multipole(5)
  E2E2 = Multipole(6);
  E3E3 = Multipole(7); M1M1 = Multipole(8)
  E1SP = Multipole(11); SPSP = Multipole(12)
  E2M2 = Multipole(13)

  f_avantseuil = (0._db, 0._db)
  Energphot = .false.
  Final_optic = .true.
  Final_tddft = .false.
  sec_atom(:) = 0._db

  allocate( Ecinetic(nspin,n_Ec) )
  allocate( Enervide(n_Ec) )

! Elaboration of energy grid 
  call dim_grille_optic(E_cut_optic,Energ_s,icheck(29),mpirank0,nenerg,nenerg_s)

  allocate( Energ(nenerg) )
  allocate( Eimag(nenerg) )

  call grille_optic(E_cut_optic,Energ,Energ_s,icheck(29),nenerg,nenerg_s)

  Eimag(:) = 0._db
  Solsing = .false.

  Xan_atom = .false.

  allocate( Vrato(nr,nlm_pot,nspin,n_V) )

  do isp = 1,nspin
    Vrato(1:nr,:,isp,1) = Vcato(1:nr,:) + Vxcato(1:nr,:,isp)
    V0bdc(isp,1) = Vhbdc + VxcbdcF(isp) + dV0bdcF(isp)   
  end do

  nge = ( nenerg - 1 ) / mpinodes + 1

  First_E = .true.
  
! Loop over photon energy
  boucle_energ: do je = 1,nge

    ie = ( je - 1 ) * mpinodes + mpirank + 1
    
    if( ie <= nenerg ) then

      if( Energ(ie) < Ephot_min - eps10 ) cycle

      T_E1E1(:,:,:,mpirank) = (0._db,0._db)
      T_E1SP(:,:,:,mpirank) = (0._db,0._db)
      T_M1M1(:,:,:,mpirank) = (0._db,0._db)
      T_E1M1(:,:,:,mpirank) = (0._db,0._db)
      T_SPSP(:,:,:,mpirank) = (0._db,0._db)
      T_E1E2(:,:,:,:,mpirank) = (0._db,0._db)
      T_E1M2(:,:,:,:,mpirank) = (0._db,0._db)
      T_E2E2(:,:,:,:,:,mpirank) = (0._db,0._db)
      T_E1E3(:,:,:,:,:,mpirank) = (0._db,0._db)
      T_E1E1_m(:,:,:,mpirank) = (0._db,0._db)
      T_E1SP_m(:,:,:,mpirank) = (0._db,0._db)
      T_SPSP_m(:,:,:,mpirank) = (0._db,0._db)
      T_M1M1_m(:,:,:,mpirank) = (0._db,0._db)
      T_E1M1_m(:,:,:,mpirank) = (0._db,0._db)
      T_E1E2_m(:,:,:,:,mpirank) = (0._db,0._db)
      T_E2E2_m(:,:,:,:,:,mpirank) = (0._db,0._db)
      T_E1E3_m(:,:,:,:,:,mpirank) = (0._db,0._db)
      if( E3E3 ) T_E3E3(:,:,:,:,:,mpirank) = (0._db,0._db)
      if( E3E3 ) T_E3E3_m(:,:,:,:,:,mpirank) = (0._db,0._db)

! Loop over valence electron energy (below Fermi level)
      do ie_s = 1,nenerg_s

        if( mpirank0 == 0 ) then
          call CPU_TIME(time)
          tp(1) = real(time,db)
        endif

        En(1) = Energ_s(ie_s)
        En(2) = Energ_s(ie_s) + Energ(ie)
        if( En(1) > E_cut_optic - eps10 ) exit
        if( En(2) < E_cut_optic - eps10 .or. En(2) > Energ_s(nenerg_s) + eps10 ) cycle

! n_Ec = 2
! ie_t = 1: index of valence level energy
! ie_t = 2: index of conduction level energy
        do ie_t = 1,n_Ec
          Enervide(ie_t) = En(ie_t) + E_Fermi
          Ecinetic(:,ie_t) = Enervide(ie_t) - V0bdc(:,1)
          if( .not. Eneg ) Ecinetic(:,ie_t) = max( Ecinetic(:,ie_t), Eclie )
        end do

        Ecmax = maxval( Ecinetic )
        call clmax(Ecmax,Rmtg,lmaxat0,lmax,numat_abs,lmaxfree)
        lmax = min(lmax,lmaxabs_t)

! Calculation of cartesian tensor

        icheck_s = max( icheck(29), icheck(20) )

        do ie_t = 1,n_Ec
          do ie_q = 2,nenerg_s
            if( Energ_s(ie_q) > En(ie_t) - eps10 ) exit
          end do
          p = ( En(ie_t) - Energ_s(ie_q - 1) ) / ( Energ_s(ie_q) - Energ_s(ie_q - 1) )
          do lm1 = 1,nlm_probe
            do iss1 = 1,nspinp
              do lm2 = 1,nlm_probe
                do iss2 = 1,nspinp
                  if( .not. Spinorbite .and. iss1 /= iss2 ) cycle
! In practice interpolation is only on the conduction electron (ie_t = 2 ) and when steps are different
                  Taull_abs(lm1,iss1,lm2,iss2,ie_t) = p * Taull_tdd(ie_q,lm1,iss1,lm2,iss2) &
                                                              + (1 - p ) * Taull_tdd(ie_q-1,lm1,iss1,lm2,iss2)
                end do
              end do
            end do
          end do
        end do

        if( ie_s == nenerg_s ) then
          delta_E = Energ_s(ie_s) - Energ_s(ie_s-1)
        elseif( ie_s == 1 ) then
          delta_E = Energ_s(ie_s+1) - Energ_s(ie_s)
        elseif(  Energ_s(ie_s+1) > E_cut_optic ) then
          delta_E = E_cut_optic - 0.5_db * ( Energ_s(ie_s) + Energ_s(ie_s-1) )
        else
          delta_E = 0.5_db * ( Energ_s(ie_s+1) - Energ_s(ie_s-1) )
        endif
! The square-root of the step energy is set to the square in subroutine tenseur_car_optic
        Taull_abs(:,:,:,:,:) = sqrt( Delta_E ) * Taull_abs(:,:,:,:,:) 
 
        ns_dipmag = 0
        allocate( Chi(nlm_probe*nspino,nlm_probe*nspino,nlm_probe*nspino,nlm_probe*nspino,2,2,ns_dipmag) )

        call Tenseur_car_optic( Chi, Core_resolved, Ecinetic, Eimag(ie), Energ(ie), Enervide, Final_tddft, &
                  Full_potential, Green_int, Hubb_a, Hubb_d, icheck_s, ip_max, ip0, Lmax_probe, Lmax_pot, &
                  Ldip, Lmoins1, Loct, Lplus1, Lqua, m_hubb, mpinodes, mpirank, mpirank0, &
                  msymdd, msymddi, msymdq, msymdqi, msymdo, msymdoi, msymoo, msymooi, msymqq, msymqqi, Multipole, &
                  n_Ec, n_oo, n_V, ns_dipmag, ninitr, nlm_pot, nlm_probe, &
                  nlm_p_fp, nr, nspin, nspino, nspinp, numat_abs, r, Relativiste, Renorm, Rmtg, Rmtsd, Rot_atom_abs, Rot_cluster, &  
                  T_E1E1_t, T_E1E1_m_t, T_E1E3_t, T_E1E3_m_t, T_E1E2_t, T_E1E2_m_t, T_E1SP_t, T_E1SP_m_t, &
                  T_E1M1_t, T_E1M1_m_t, T_M1M1_t, T_M1M1_m_t, T_E3E3_t, T_E3E3_m_t, T_E2E2_t, T_E2E2_m_t, T_SPSP_t, T_SPSP_m_t, &   
                  Solsing, Spinorbite, Taull_abs, V_hubb, V_intmax, V0bdc, Vrato, Ylm_comp )

        deallocate( Chi )

        if( E1E1 ) T_E1E1(:,:,:,mpirank) = T_E1E1(:,:,:,mpirank) + T_E1E1_t(:,:,:,mpirank)
        if( E1E1 ) T_E1E1_m(:,:,:,mpirank) = T_E1E1_m(:,:,:,mpirank) + T_E1E1_m_t(:,:,:,mpirank)
        if( E1SP ) T_E1SP(:,:,:,mpirank) = T_E1SP(:,:,:,mpirank) + T_E1SP_t(:,:,:,mpirank)
        if( E1SP ) T_E1SP_m(:,:,:,mpirank) = T_E1SP_m(:,:,:,mpirank) + T_E1SP_m_t(:,:,:,mpirank)
        if( SPSP ) T_SPSP(:,:,:,mpirank) = T_SPSP(:,:,:,mpirank) + T_SPSP_t(:,:,:,mpirank)
        if( SPSP ) T_SPSP_m(:,:,:,mpirank) = T_SPSP_m(:,:,:,mpirank) + T_SPSP_m_t(:,:,:,mpirank)
        if( M1M1 ) T_M1M1(:,:,:,mpirank) = T_M1M1(:,:,:,mpirank) + T_M1M1_t(:,:,:,mpirank)
        if( M1M1 ) T_M1M1_m(:,:,:,mpirank) = T_M1M1_m(:,:,:,mpirank) + T_M1M1_m_t(:,:,:,mpirank)
        if( E1M1 ) T_E1M1(:,:,:,mpirank) = T_E1M1(:,:,:,mpirank) + T_E1M1_t(:,:,:,mpirank)
        if( E1M1 ) T_E1M1_m(:,:,:,mpirank) = T_E1M1_m(:,:,:,mpirank) + T_E1M1_m_t(:,:,:,mpirank)
        if( E1E2 ) T_E1E2(:,:,:,:,mpirank) = T_E1E2(:,:,:,:,mpirank) + T_E1E2_t(:,:,:,:,mpirank)
        if( E1M2 ) T_E1M2(:,:,:,:,mpirank) = T_E1M2(:,:,:,:,mpirank) + T_E1M2_t(:,:,:,:,mpirank)
        if( E1E2 ) T_E1E2_m(:,:,:,:,mpirank) = T_E1E2_m(:,:,:,:,mpirank) + T_E1E2_m_t(:,:,:,:,mpirank)
        if( E1M2 ) T_E1M2_m(:,:,:,:,mpirank) = T_E1M2_m(:,:,:,:,mpirank) + T_E1M2_m_t(:,:,:,:,mpirank)
        if( E2E2 ) T_E2E2(:,:,:,:,:,mpirank) = T_E2E2(:,:,:,:,:,mpirank) + T_E2E2_t(:,:,:,:,:,mpirank)
        if( E2E2 ) T_E2E2_m(:,:,:,:,:,mpirank) = T_E2E2_m(:,:,:,:,:,mpirank) + T_E2E2_m_t(:,:,:,:,:,mpirank)
        if( E2M2 ) T_E2M2(:,:,:,:,:,mpirank) = T_E2M2(:,:,:,:,:,mpirank) + T_E2M2_t(:,:,:,:,:,mpirank)
        if( E2M2 ) T_E2M2_m(:,:,:,:,:,mpirank) = T_E2M2_m(:,:,:,:,:,mpirank) + T_E2M2_m_t(:,:,:,:,:,mpirank)
        if( E1E3 ) T_E1E3(:,:,:,:,:,mpirank) = T_E1E3(:,:,:,:,:,mpirank) + T_E1E3_t(:,:,:,:,:,mpirank)
        if( E1E3 ) T_E1E3_m(:,:,:,:,:,mpirank) = T_E1E3_m(:,:,:,:,:,mpirank) + T_E1E3_m_t(:,:,:,:,:,mpirank)
        if( E3E3 ) T_E3E3(:,:,:,:,:,mpirank) = T_E3E3(:,:,:,:,:,mpirank) + T_E3E3_t(:,:,:,:,:,mpirank)
        if( E3E3 ) T_E3E3_m(:,:,:,:,:,mpirank) = T_E3E3_m(:,:,:,:,:,mpirank) + T_E3E3_m_t(:,:,:,:,:,mpirank)

        if( mpirank0 == 0 ) then
          call CPU_TIME(time)
          tp(2) = real(time,db)
          Time_rout(11) = Time_rout(11) + tp(2) - tp(1)
        endif

      end do ! End of loop over electron energies

    endif  ! arrival when ie > nenerg

    if( mpinodes > 1 ) then
      l0_nrixs = 0; lmax_nrixs = 0; nq_nrixs = 0
      allocate( S_nrixs(nq_nrixs,(lmax_nrixs+1)**2,(lmax_nrixs+1)**2,ninitr,0:mpinodes-1) )
      call MPI_RECV_all( L0_nrixs, lmax_nrixs, mpinodes, mpirank, mpirank0, Multipole, n_oo, ninitr, &
                        nq_nrixs, S_nrixs, T_E1E1, T_E1E2, T_E1E3, T_E1M1, T_E1M2, T_E1SP, T_E2E2, T_E2M2, &
                        T_E3E3, T_M1M1, T_SPSP )
      deallocate( S_nrixs )
    endif

    if( mpirank0 /= 0 ) cycle
    call CPU_TIME(time)
    tp(2) = real(time,db)

    icheck_s = max( icheck(29), icheck(21) )

    do ie_computer = 0,mpinodes-1

      ie = ( je - 1 ) * mpinodes + ie_computer + 1

      if( ie > nenerg ) exit

      call Write_coabs( Abs_cap, Abs_in_Bulk_roughness, Abs_SRXRD, Abs_SRXRD_abs, Abs_U_iso, Allsite, angxyz, axyz, &
                Bragg_abs, .false., Cartesian_tensor, &
                Core_resolved, Dafs, Dafs_bio, E_cut_optic, Energ, Energphot, Epsii, Eseuil, Final_tddft, First_E, &
                f_avantseuil, Full_self_abs, Green_int, Header, hkl_dafs, iabsorig, icheck_s, ie, ie_computer, &
                index_z_SRXRD_abs, Int_tens, isigpi, isymeq, jseuil, Length_abs, ltypcal, Matper, Moyenne, mpinodes, &
                multi_0, Multipole, n_abs_rgh, n_bulk_sup, n_multi_run, n_z_SRXRD, n_z_SRXRD_abs, &
                n_z_SRXRD_max_abs, n_zc_SRXRD_abs, n_max, n_oo, n_tens_max, natomsym, nbseuil, ncolm, ncolr, ncolt, &
                nenerg, ninit1, ninitr, nomabs, nomfich, nomfich_cal_conv, nomfich_s, Non_resonant_Bulk, nphi_dafs, &
                npldafs, nphim, nplr, nplrm, nseuil, nspinp, numat_abs, nxanout, pdp, phdt, pol ,poldafse, poldafss, &
                Pos_z_SRXRD, Pos_z_SRXRD_abs, sec_atom, T_E1E1, T_E1E1_m, T_E1E2, T_E1E2_m, T_E1E3, T_E1E3_m, T_E1M1, T_E1M1_m, &
                T_E1M2, T_E1M2_m, T_E1SP, T_E1SP_m, T_E2E2, T_E2E2_m, T_E2M2, T_E2M2_m, T_E3E3, T_E3E3_m, T_M1M1, T_M1M1_m, &
                T_SPSP, T_SPSP_m, Self_abs, &
                Spherical_signal, Spherical_tensor, Spinorbite, Sum_Bragg_abs_f0, Sum_Bragg_nonabs_f, Sum_Bragg_NR_Bulk, &
                Sum_Bragg_NR_T, Surface_ref, Taux_eq, V0muf, Vecdafse, Vecdafss, Vec, Volume_maille, Xan_atom ) 

      First_E = .false.
      
    end do

    call CPU_TIME(time)
    tp(3) = real(time,db)
    Time_rout(13) = Time_rout(13) + tp(3) - tp(2)

  end do boucle_energ

  deallocate( En, Taull_abs, V0bdc, Vrato )
  deallocate( Ecinetic, Eimag, Energ, Enervide )

  return
  100 format(/1x,120('-')//,' Cycle Optic')
end

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

! Subroutine giving the dimension of the photon energy table in optics

subroutine dim_grille_optic(E_cut,Energ_s,icheck,mpirank0,nenerg,nenerg_s)

  use declarations
  implicit none

  integer:: icheck, mpirank0, nenerg, nenerg_s

  real(kind=db), intent(in):: E_cut
  real(kind=db), dimension(nenerg_s), intent(in):: Energ_s

  integer:: ie, ipr, ne

  if( E_cut < Energ_s(1) .and. mpirank0 == 0 ) then
    call write_error
    do ipr = 6,9,3
      if( ipr == 3 .and. icheck == 0 ) cycle
      write(ipr,110) E_cut*rydb, Energ_s(1)*rydb
    end do
    stop
  endif

  if( E_cut > Energ_s(nenerg_s) .and. mpirank0 == 0 ) then
    call write_error
    do ipr = 6,9,3
      if( ipr == 3 .and. icheck == 0 ) cycle
      write(ipr,120) E_cut*rydb, Energ_s(nenerg_s)*rydb
    end do
    stop
  endif

  do ie = 1,nenerg_s
    if( Energ_s(ie) > E_cut + eps10 ) exit
  end do
  ne = ie - 1

  nenerg = nenerg_s - ne

  return
  110 format(///' E_Fermi =',f7.3,' eV < First energy =',f7.3,' eV',// ' Increase the energy range of calculation !'//)
  120 format(///' E_Fermi =',f7.3,' eV > Last energy =',f7.3,' eV',// ' Increase the energy range of calculation !'//)
end

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

! Subroutine giving the photon energy table for optics.

subroutine grille_optic(E_cut,Energ,Energ_s,icheck,nenerg, nenerg_s)

  use declarations
  implicit none

  integer,intent(in):: icheck, nenerg, nenerg_s

  real(kind=db), intent(in):: E_cut
  real(kind=db), dimension(nenerg_s), intent(in):: Energ_s
  real(kind=db), dimension(nenerg), intent(out):: Energ

  integer:: ie, ne

  do ie = 1,nenerg_s
    if( Energ_s(ie) > E_cut + eps10 ) exit
  end do
  ne = ie - 1

  do ie = 1,nenerg
    Energ(ie) = Energ_s(ie+ne) - Energ_s(ne)
  end do

  if( icheck > 2 ) then
    write(3,100)
    write(3,110)
    write(3,120) Energ(:)*rydb
  end if

  return
  100 format(/' ---- Grille_optic -------',100('-'))
  110 format(/'The energy grid for the optic calculation:')
  120 format(5f13.7)
end

