module collision
! contains the Hermite-Laguerre collision operators. Solved using COSOlver.
IMPLICIT NONE

PUBLIC :: compute_TColl
PUBLIC :: LenardBernstein_e, LenardBernstein_i!, LenardBernstein GK
PUBLIC :: DoughertyGK_e, DoughertyGK_i!, Dougherty GK
PUBLIC :: load_COSOlver_mat
PUBLIC :: apply_COSOlver_mat_e, apply_COSOlver_mat_i

CONTAINS

  !******************************************************************************!
  !! Lenard Bernstein collision operator for electrons
  !******************************************************************************!
  SUBROUTINE LenardBernstein_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
    USE fields, ONLY: moments_e
    USE grid,   ONLY: parray_e, jarray_e, kxarray, kyarray
    USE basic
    USE model,  ONLY: sigmae2_taue_o2, nu_ee
    USE time_integration, ONLY : updatetlevel
    IMPLICIT NONE
    INTEGER,     INTENT(IN)    :: ip_,ij_,ikx_,iky_,iz_
    COMPLEX(dp), INTENT(OUT)   :: TColl_

    REAL(dp)    :: j_dp, p_dp, be_2

    !** Auxiliary variables **
    p_dp      = REAL(parray_e(ip_),dp)
    j_dp      = REAL(jarray_e(ij_),dp)
    ! be_2      = (kxarray(ikx_)**2 + kyarray(iky_)**2) * sigmae2_taue_o2 ! this is (be/2)^2

    !** Assembling collison operator **
    ! -nuee (p + 2j) Nepj
    TColl_ = -nu_ee * (p_dp + 2._dp*j_dp)*moments_e(ip_,ij_,ikx_,iky_,iz_,updatetlevel)

  END SUBROUTINE LenardBernstein_e

    !******************************************************************************!
    !! Lenard Bernstein collision operator for electrons
    !******************************************************************************!
    SUBROUTINE LenardBernstein_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
      USE fields, ONLY: moments_i
      USE grid,   ONLY: parray_i, jarray_i, kxarray, kyarray
      USE basic
      USE model,  ONLY: sigmai2_taui_o2, nu_i
      USE time_integration, ONLY : updatetlevel
      IMPLICIT NONE
      INTEGER,     INTENT(IN)    :: ip_,ij_,ikx_,iky_,iz_
      COMPLEX(dp), INTENT(OUT)   :: TColl_

      REAL(dp)    :: j_dp, p_dp, bi_2

      !** Auxiliary variables **
      p_dp      = REAL(parray_i(ip_),dp)
      j_dp      = REAL(jarray_i(ij_),dp)
      ! bi_2      = (kxarray(ikx_)**2 + kyarray(iky_)**2) * sigmai2_taui_o2 ! this is (bi/2)^2

      !** Assembling collison operator **
      ! -nuii (p + 2j) Nipj
      TColl_ = -nu_i * (p_dp + 2._dp*j_dp)*moments_i(ip_,ij_,ikx_,iky_,iz_,updatetlevel)

    END SUBROUTINE LenardBernstein_i

  !******************************************************************************!
  !! Doughtery gyrokinetic collision operator for electrons
  !******************************************************************************!
  SUBROUTINE DoughertyGK_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
    USE fields, ONLY: moments_e, phi
    USE grid,   ONLY: parray_e, jarray_e, kxarray, kyarray, Jmaxe, ip0_e, ip1_e, ip2_e
    USE array,  ONLY: kernel_e
    USE basic
    USE model,  ONLY: sigmae2_taue_o2, qe_taue, nu_ee
    USE time_integration, ONLY : updatetlevel
    IMPLICIT NONE
    INTEGER,     INTENT(IN)    :: ip_,ij_,ikx_,iky_,iz_
    COMPLEX(dp), INTENT(OUT)   :: TColl_

    COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_
    COMPLEX(dp) :: nadiab_moment_0j
    REAL(dp)    :: Knp0, Knp1, Knm1
    INTEGER     :: in_
    REAL(dp)    :: n_dp, j_dp, p_dp, be_, be_2

    !** Auxiliary variables **
    p_dp      = REAL(parray_e(ip_),dp)
    j_dp      = REAL(jarray_e(ij_),dp)
    be_2      = (kxarray(ikx_)**2 + kyarray(iky_)**2) * sigmae2_taue_o2 ! this is (be/2)^2
    be_       = 2_dp*SQRT(be_2) ! this is be

    !** Assembling collison operator **
    ! Velocity-space diffusion (similar to Lenard Bernstein)
    ! -nuee (p + 2j + b^2/2) Nepj
    TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*be_2)*moments_e(ip_,ij_,ikx_,iky_,iz_,updatetlevel)

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
      ! Get adiabatic moment
      TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*be_2) * qe_taue * Kernel_e(ij_,ikx_,iky_,iz_)*phi(ikx_,iky_,iz_)
        !** build required fluid moments **
        n_     = 0._dp
        upar_  = 0._dp; uperp_ = 0._dp
        Tpar_  = 0._dp; Tperp_ = 0._dp
        DO in_ = 1,jmaxe+1
          n_dp = REAL(in_-1,dp)
          ! Store the kernels for sparing readings
          Knp0 =  Kernel_e(in_  ,ikx_,iky_,iz_)
          Knp1 =  Kernel_e(in_+1,ikx_,iky_,iz_)
          Knm1 =  Kernel_e(in_-1,ikx_,iky_,iz_)
          ! Nonadiabatic moments (only different from moments when p=0)
          nadiab_moment_0j   = moments_e(ip0_e,in_  ,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
          ! Density
          n_     = n_     + Knp0 * nadiab_moment_0j
          ! Perpendicular velocity
          uperp_ = uperp_ + be_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
          ! Parallel temperature
          Tpar_  = Tpar_  + Knp0 * (SQRT2*moments_e(ip2_e,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
          ! Perpendicular temperature
          Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
        ENDDO
      T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
      ! Add energy restoring term
      TColl_ = TColl_ + T_* 4._dp *  j_dp          * Kernel_e(ij_  ,ikx_,iky_,iz_)
      TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_,iz_)
      TColl_ = TColl_ - T_* 2._dp *  j_dp          * Kernel_e(ij_-1,ikx_,iky_,iz_)
      TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_,iz_) - Kernel_e(ij_-1,ikx_,iky_,iz_))
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ELSEIF( p_dp .eq. 1 ) THEN ! kronecker p1
      !** build required fluid moments **
      upar_  = 0._dp
      DO in_ = 1,jmaxe+1
        ! Parallel velocity
         upar_  = upar_  + Kernel_e(in_,ikx_,iky_,iz_) * moments_e(ip1_e,in_,ikx_,iky_,iz_,updatetlevel)
      ENDDO
      TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,iz_)
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ELSEIF( p_dp .eq. 2 ) THEN ! kronecker p2
      !** build required fluid moments **
      n_     = 0._dp
      upar_  = 0._dp; uperp_ = 0._dp
      Tpar_  = 0._dp; Tperp_ = 0._dp
      DO in_ = 1,jmaxe+1
        n_dp = REAL(in_-1,dp)
        ! Store the kernels for sparing readings
        Knp0 =  Kernel_e(in_  ,ikx_,iky_,iz_)
        Knp1 =  Kernel_e(in_+1,ikx_,iky_,iz_)
        Knm1 =  Kernel_e(in_-1,ikx_,iky_,iz_)
        ! Nonadiabatic moments (only different from moments when p=0)
        nadiab_moment_0j = moments_e(ip0_e,in_,ikx_,iky_,iz_,updatetlevel) + qe_taue*Knp0*phi(ikx_,iky_,iz_)
        ! Density
        n_     = n_     + Knp0 * nadiab_moment_0j
        ! Parallel temperature
        Tpar_  = Tpar_  + Knp0 * (SQRT2*moments_e(ip2_e,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
        ! Perpendicular temperature
        Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
      ENDDO
      T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
      TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_,iz_)
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ENDIF
    ! Multiply by electron-electron collision coefficient
    TColl_ = nu_ee * TColl_

  END SUBROUTINE DoughertyGK_e
  !******************************************************************************!
  !! Doughtery gyrokinetic collision operator for ions
  !******************************************************************************!
  SUBROUTINE DoughertyGK_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
    USE fields, ONLY: moments_i, phi
    USE grid,   ONLY: parray_i, jarray_i, kxarray, kyarray, Jmaxi, ip0_i, ip1_i, ip2_i
    USE array,  ONLY: kernel_i
    USE basic
    USE model,  ONLY: sigmai2_taui_o2, qi_taui, nu_i
    USE time_integration, ONLY : updatetlevel
    IMPLICIT NONE
    INTEGER,     INTENT(IN)    :: ip_,ij_,ikx_,iky_,iz_
    COMPLEX(dp), INTENT(OUT)   :: TColl_

    COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_
    COMPLEX(dp) :: bi_, bi_2
    COMPLEX(dp) :: nadiab_moment_0j
    REAL(dp)    :: Knp0, Knp1, Knm1
    INTEGER     :: in_
    REAL(dp)    :: n_dp, j_dp, p_dp

    !** Auxiliary variables **
    p_dp      = REAL(parray_i(ip_),dp)
    j_dp      = REAL(jarray_i(ij_),dp)
    bi_2      = (kxarray(ikx_)**2 + kyarray(iky_)**2) * sigmai2_taui_o2 ! this is (bi/2)^2
    bi_       = 2_dp*SQRT(bi_2) ! this is be

    !** Assembling collison operator **
    ! Velocity-space diffusion (similar to Lenard Bernstein)
    ! -nui (p + 2j + b^2/2) Nipj
    TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*bi_2)*moments_i(ip_,ij_,ikx_,iky_,iz_,updatetlevel)

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
      ! Get adiabatic moment
      TColl_ = TColl_ - (p_dp + 2._dp*j_dp + 2._dp*bi_2) * qi_taui * Kernel_i(ij_,ikx_,iky_,iz_)*phi(ikx_,iky_,iz_)
        !** build required fluid moments **
        n_     = 0._dp
        upar_  = 0._dp; uperp_ = 0._dp
        Tpar_  = 0._dp; Tperp_ = 0._dp
        DO in_ = 1,jmaxi+1
          n_dp = REAL(in_-1,dp)
          ! Store the kernels for sparing readings
          Knp0 =  Kernel_i(in_  ,ikx_,iky_,iz_)
          Knp1 =  Kernel_i(in_+1,ikx_,iky_,iz_)
          Knm1 =  Kernel_i(in_-1,ikx_,iky_,iz_)
          ! Nonadiabatic moments (only different from moments when p=0)
          nadiab_moment_0j   = moments_i(ip0_i,in_  ,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
          ! Density
          n_     = n_     + Knp0 * nadiab_moment_0j
          ! Perpendicular velocity
          uperp_ = uperp_ + bi_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
          ! Parallel temperature
          Tpar_  = Tpar_  + Knp0 * (SQRT2*moments_i(ip2_i,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
          ! Perpendicular temperature
          Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
        ENDDO
        T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
      ! Add energy restoring term
      TColl_ = TColl_ + T_* 4._dp *  j_dp          * Kernel_i(ij_  ,ikx_,iky_,iz_)
      TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_,iz_)
      TColl_ = TColl_ - T_* 2._dp *  j_dp          * Kernel_i(ij_-1,ikx_,iky_,iz_)
      TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_,iz_) - Kernel_i(ij_-1,ikx_,iky_,iz_))
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ELSEIF( p_dp .eq. 1 ) THEN ! kxonecker p1
      !** build required fluid moments **
      upar_  = 0._dp
      DO in_ = 1,jmaxi+1
        ! Parallel velocity
         upar_  = upar_  + Kernel_i(in_,ikx_,iky_,iz_) * moments_i(ip1_i,in_,ikx_,iky_,iz_,updatetlevel)
      ENDDO
      TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,iz_)
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ELSEIF( p_dp .eq. 2 ) THEN ! kxonecker p2
      !** build required fluid moments **
      n_     = 0._dp
      upar_  = 0._dp; uperp_ = 0._dp
      Tpar_  = 0._dp; Tperp_ = 0._dp
      DO in_ = 1,jmaxi+1
        n_dp = REAL(in_-1,dp)
        ! Store the kernels for sparing readings
        Knp0 =  Kernel_i(in_  ,ikx_,iky_,iz_)
        Knp1 =  Kernel_i(in_+1,ikx_,iky_,iz_)
        Knm1 =  Kernel_i(in_-1,ikx_,iky_,iz_)
        ! Nonadiabatic moments (only different from moments when p=0)
        nadiab_moment_0j = moments_i(ip0_i,in_,ikx_,iky_,iz_,updatetlevel) + qi_taui*Knp0*phi(ikx_,iky_,iz_)
        ! Density
        n_     = n_     + Knp0 * nadiab_moment_0j
        ! Parallel temperature
        Tpar_  = Tpar_  + Knp0 * (SQRT2*moments_i(ip2_i,in_,ikx_,iky_,iz_,updatetlevel) + nadiab_moment_0j)
        ! Perpendicular temperature
        Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
      ENDDO
      T_  = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
      TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_,iz_)
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ENDIF
    ! Multiply by ion-ion collision coefficient
    TColl_ = nu_i * TColl_

  END SUBROUTINE DoughertyGK_i

  !******************************************************************************!
  !! compute the collision terms in a (Np x Nj x Nkx x Nky) matrix all at once
  !******************************************************************************!
  SUBROUTINE compute_TColl
    USE fields
    USE grid
    USE array
    USE basic
    USE prec_const
    USE time_integration
    USE model
    USE utility
    IMPLICIT NONE
    COMPLEX(dp), DIMENSION(1:total_np_e)   :: local_sum_e, buffer_e, total_sum_e
    COMPLEX(dp), DIMENSION(ips_e:ipe_e) :: TColl_distr_e
    COMPLEX(dp), DIMENSION(1:total_np_i)   :: local_sum_i, buffer_i, total_sum_i
    COMPLEX(dp), DIMENSION(ips_i:ipe_i) :: TColl_distr_i
    COMPLEX(dp) :: TColl
    INTEGER :: ikxs_C, ikxe_C, ikys_C, ikye_C

    ! Execution time start
    CALL cpu_time(t0_coll)

    IF (ABS(CO) .GE. 2) THEN !compute only if COSOlver matrices are used
      DO ikx = ikxs,ikxe
        DO iky = ikys,ikye
          DO iz = izs,ize
            ! Electrons
            DO ij = 1,Jmaxe+1
              ! Loop over all p to compute sub collision term
              DO ip = 1,total_np_e
                CALL apply_COSOlver_mat_e(ip,ij,ikx,iky,iz,TColl)
                local_sum_e(ip) = TColl
              ENDDO
              IF (num_procs_p .GT. 1) THEN
                ! Sum up all the sub collision terms on root 0
                CALL MPI_REDUCE(local_sum_e, buffer_e, total_np_e, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
                ! distribute the sum over the process among p
                CALL MPI_SCATTERV(buffer_e, counts_np_e, displs_np_e, MPI_DOUBLE_COMPLEX,&
                                  TColl_distr_e, local_np_e, MPI_DOUBLE_COMPLEX,&
                                  0, comm_p, ierr)
              ELSE
                TColl_distr_e = local_sum_e
              ENDIF
              ! Write in output variable
              DO ip = ips_e,ipe_e
                TColl_e(ip,ij,ikx,iky,iz) = TColl_distr_e(ip)
              ENDDO
            ENDDO
            ! Ions
            DO ij = 1,Jmaxi+1
              DO ip = 1,total_np_i
                CALL apply_COSOlver_mat_i(ip,ij,ikx,iky,iz,TColl)
                local_sum_i(ip) = TColl
              ENDDO
              IF (num_procs_p .GT. 1) THEN
                ! Reduce the local_sums to root = 0
                CALL MPI_REDUCE(local_sum_i, buffer_i, total_np_i, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
                ! buffer contains the entire collision term along p, we scatter it between
                ! the other processes (use of scatterv since Pmax/Np is not an integer)
                CALL MPI_SCATTERV(buffer_i, counts_np_i, displs_np_i, MPI_DOUBLE_COMPLEX,&
                                  TColl_distr_i, local_np_i, MPI_DOUBLE_COMPLEX, &
                                  0, comm_p, ierr)
              ELSE
                TColl_distr_i = local_sum_i
              ENDIF
              ! Write in output variable
              DO ip = ips_i,ipe_i
                TColl_i(ip,ij,ikx,iky,iz) = TColl_distr_i(ip)
              ENDDO
            ENDDO
          ENDDO
        ENDDO
      ENDDO
    ENDIF

    ! Execution time end
    CALL cpu_time(t1_coll)
    tc_coll = tc_coll + (t1_coll - t0_coll)
  END SUBROUTINE compute_TColl

  !******************************************************************************!
  !!!!!!! Compute ion collision term
  !******************************************************************************!
  SUBROUTINE apply_COSOlver_mat_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
    USE fields, ONLY: moments_e, moments_i
    USE grid
    USE array
    USE basic
    USE time_integration, ONLY: updatetlevel
    USE utility
    USE model, ONLY: CO, nu_e, nu_ee, CLOS
    IMPLICIT NONE

    INTEGER,     INTENT(IN)  :: ip_, ij_ ,ikx_, iky_, iz_
    COMPLEX(dp), INTENT(OUT) :: TColl_

    INTEGER     :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C
    p_int = parray_e_full(ip_); j_int = jarray_e_full(ij_);

    IF (CO .GT. 0) THEN ! GK operator (k-dependant)
      ikx_C = ikx_; iky_C = iky_
    ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for every k)
      ikx_C = 1;   iky_C = 1
    ENDIF

    TColl_ = 0._dp ! Initialization of the local sum

    ! sum the electron-self and electron-ion test terms
    ploopee: DO ip2 = ips_e,ipe_e
      p2_int = parray_e(ip2)
      jloopee: DO ij2 = ijs_e,ije_e
        j2_int = jarray_e(ij2)
        IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))&
        TColl_ = TColl_ + moments_e(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
           *( nu_e  * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C) &
             +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C))
      ENDDO jloopee
    ENDDO ploopee

    ! sum the electron-ion field terms
    ploopei: DO ip2 = ips_i,ipe_i
      p2_int = parray_i(ip2)
      jloopei: DO ij2 = ijs_i,ije_i
        j2_int = jarray_i(ij2)
        IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))&
        TColl_ = TColl_ + moments_i(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
          *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int),ikx_C, iky_C))
      END DO jloopei
    ENDDO ploopei

  END SUBROUTINE apply_COSOlver_mat_e

  !******************************************************************************!
  !!!!!!! Compute ion collision term
  !******************************************************************************!
  SUBROUTINE apply_COSOlver_mat_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
    USE fields, ONLY : moments_e, moments_i
    USE grid
    USE array
    USE basic
    USE time_integration, ONLY : updatetlevel
    USE utility
    USE model, ONLY: CO, nu_i, nu_ie, CLOS
    IMPLICIT NONE
    INTEGER,     INTENT(IN)    :: ip_, ij_ ,ikx_, iky_, iz_
    COMPLEX(dp), INTENT(OUT)   :: TColl_

    INTEGER     :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C
    p_int = parray_i_full(ip_); j_int = jarray_i_full(ij_);

    IF (CO .GT. 0) THEN ! GK operator (k-dependant)
      ikx_C = ikx_; iky_C = iky_
    ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for every k)
      ikx_C = 1;   iky_C = 1
    ENDIF

    TColl_ = 0._dp ! Initialization
    ! sum the ion-self and ion-electron test terms
    ploopii: DO ip2 = ips_i,ipe_i
      p2_int = parray_i(ip2)
      jloopii: DO ij2 = ijs_i,ije_i
        j2_int = jarray_i(ij2)
        IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))&
        TColl_ = TColl_ + moments_i(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
            *( nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C) &
              +nu_i  * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C))
      ENDDO jloopii
    ENDDO ploopii

    ploopie: DO ip2 = ips_e,ipe_e ! sum the ion-electron field terms
      p2_int = parray_e(ip2)
      jloopie: DO ij2 = ijs_e,ije_e
        j2_int = jarray_e(ij2)
        IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))&
        TColl_ = TColl_ + moments_e(ip2,ij2,ikx_,iky_,iz_,updatetlevel) &
          *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int), ikx_C, iky_C))
      ENDDO jloopie
    ENDDO ploopie

  END SUBROUTINE apply_COSOlver_mat_i

    !******************************************************************************!
    !!!!!!! Load the collision matrix coefficient table from COSOlver results
    !******************************************************************************!
    SUBROUTINE load_COSOlver_mat ! Load a sub matrix from iCa files (works for pmaxa,jmaxa<=P_full,J_full)
      use futils
      use initial_par
      USE grid
      USE array, ONLY: Ceepj, Ciipj, CeipjF, CeipjT, CiepjF, CiepjT
      USE basic
      USE time_integration, ONLY : updatetlevel
      USE utility
      USE model, ONLY: CO, NON_LIN, sigmae2_taue_o2, sigmai2_taui_o2
      IMPLICIT NONE
      ! Indices for row and columns of the COSOlver matrix (4D compressed 2D matrices)
      INTEGER :: irow_sub, irow_full, icol_sub, icol_full
      INTEGER :: fid ! file indexation

      INTEGER :: ip_e, ij_e, il_e, ik_e, ikps_C, ikpe_C                  ! indices for electrons loops
      REAL(dp), DIMENSION(2) :: dims_e
      INTEGER :: pdime, jdime                                            ! dimensions of the COSOlver matrices
      REAL(dp), DIMENSION(:,:),   ALLOCATABLE :: Ceepj_full, CeipjT_full ! To load the entire matrix
      REAL(dp), DIMENSION(:,:),   ALLOCATABLE :: CeipjF_full             ! ''
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ceepj__kp, CeipjT_kp    ! To store the coeff that will be used along kperp
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: CeipjF_kp               ! ''
      INTEGER :: ip_i, ij_i, il_i, ik_i                                  ! same for ions
      INTEGER,  DIMENSION(2) :: dims_i
      INTEGER :: pdimi, jdimi                                            ! dimensions of the COSOlver matrices
      REAL(dp), DIMENSION(:,:),   ALLOCATABLE :: Ciipj_full, CiepjT_full ! .
      REAL(dp), DIMENSION(:,:),   ALLOCATABLE :: CiepjF_full             ! .
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ciipj__kp, CiepjT_kp    ! .
      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: CiepjF_kp               ! .
      INTEGER  :: NFLR

      REAL(dp), DIMENSION(:),     ALLOCATABLE :: kp_grid_mat             ! kperp grid of the matrices
      INTEGER  :: ikp_next, ikp_prev, nkp_mat, ikp_mat
      REAL(dp) ::  kp_next,  kp_prev, kperp_sim, kperp_mat, zerotoone, be_2, bi_2

      CHARACTER(len=128) :: var_name, kperp_string, ikp_string

      LOGICAL     :: CO_AA_ONLY = .false. ! Flag to remove ei ie collision

      !! Some terminal info
      IF (CO .EQ. 2) THEN
        IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Sugama matrix ==='
      ELSEIF(CO .EQ. 3) THEN
        IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK pitch angle matrix ==='
      ELSEIF(CO .EQ. 4) THEN
        IF (my_id .EQ. 0) WRITE(*,*) '=== Load GK Coulomb matrix ==='
      ELSEIF(CO .EQ. -2) THEN
        IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Sugama matrix ==='
      ELSEIF(CO .EQ. -3) THEN
        IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK pitch angle matrix ==='
      ELSEIF(CO .EQ. -4) THEN
        IF (my_id .EQ. 0) WRITE(*,*) '=== Load DK Coulomb matrix ==='
      ENDIF

      ! Opening the compiled cosolver matrices results
      if(my_id.EQ.0)write(*,*) mat_file
      CALL openf(mat_file,fid, 'r', 'D', mpicomm=comm_p);

      ! Get matrices dimensions (polynomials degrees and kperp grid)
      CALL getarr(fid, '/dims_e', dims_e) ! Get the electron polynomial degrees
      pdime = dims_e(1); jdime = dims_e(2);
      CALL getarr(fid, '/dims_i', dims_i) ! Get the ion      polynomial degrees
      pdimi = dims_i(1); jdimi = dims_i(2);
      IF ( ((pdime .LT. pmaxe) .OR. (jdime .LT. jmaxe)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Pe,Je Matrix too small !!'
      IF ( ((pdimi .LT. pmaxi) .OR. (jdimi .LT. jmaxi)) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Pi,Ji Matrix too small !!'

      CALL getsize(fid, '/coordkperp', nkp_mat) ! Get the dimension kperp grid of the matrices
      CALL allocate_array(kp_grid_mat, 1,nkp_mat)
      CALL getarr(fid, '/coordkperp', kp_grid_mat)
      IF (NON_LIN) THEN ! check that we have enough kperps mat
        IF ( (kp_grid_mat(nkp_mat) .LT. two_third_kpmax) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!'
      ELSE
        IF ( (kp_grid_mat(nkp_mat) .LT. kp_max) .AND. (my_id .EQ. 0)) WRITE(*,*) '!! Matrix kperp grid too small !!'
      ENDIF

      IF (CO .GT. 0) THEN ! GK operator (k-dependant)
        ikps_C = 1; ikpe_C = nkp_mat
      ELSEIF (CO .LT. 0) THEN ! DK operator (only one mat for all k)
        ikps_C = 1; ikpe_C = 1
      ENDIF

      CALL allocate_array(  Ceepj__kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C)
      CALL allocate_array(  CeipjT_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C)
      CALL allocate_array(  CeipjF_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C)
      CALL allocate_array(  Ciipj__kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C)
      CALL allocate_array(  CiepjT_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C)
      CALL allocate_array(  CiepjF_kp, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikps_C,ikpe_C)

      DO ikp = ikps_C,ikpe_C ! Loop over everz kperp values
        ! we put zeros if kp>2/3kpmax because thoses frequencies are filtered through AA
        IF(kp_grid_mat(ikp) .GT. two_third_kpmax .AND. NON_LIN) THEN
          CiepjT_kp(:,:,ikp) = 0._dp
          CiepjF_kp(:,:,ikp) = 0._dp
          CeipjT_kp(:,:,ikp) = 0._dp
          CeipjF_kp(:,:,ikp) = 0._dp
          Ceepj__kp(:,:,ikp) = 0._dp
          Ciipj__kp(:,:,ikp) = 0._dp
        ELSE
          ! Kperp value in string format
          IF (CO .GT. 0) THEN
            write(ikp_string,'(i5.5)') ikp-1
          ELSE
            write(ikp_string,'(i5.5)') 0
          ENDIF
          !!!!!!!!!!!! E-E matrices !!!!!!!!!!!!
          ! get the self electron colision matrix
          ! Allocate space for storing full collision matrix
          CALL allocate_array(  Ceepj_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
          ! Naming of the array to load (kperp dependant)
          WRITE(var_name,'(a,a)') TRIM(ADJUSTL(ikp_string)),'/Caapj/Ceepj'
          CALL getarr(fid, var_name, Ceepj_full) ! get array (moli format)
          ! Fill sub array with the usefull polynmial degrees only
          DO ip_e = 0,pmaxe ! Loop over rows
          DO ij_e = 0,jmaxe
                irow_sub  = (jmaxe +1)*ip_e + ij_e +1
                irow_full = (jdime +1)*ip_e + ij_e +1
                DO il_e = 0,pmaxe ! Loop over columns
                DO ik_e = 0,jmaxe
                      icol_sub  = (jmaxe +1)*il_e + ik_e +1
                      icol_full = (jdime +1)*il_e + ik_e +1
                      Ceepj__kp (irow_sub,icol_sub,ikp) = Ceepj_full (irow_full,icol_full)
                ENDDO
                ENDDO
          ENDDO
          ENDDO
          DEALLOCATE(Ceepj_full)

          !!!!!!!!!!!!!!! I-I matrices !!!!!!!!!!!!!!
          ! get the self electron colision matrix
          CALL allocate_array(  Ciipj_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
          WRITE(var_name,'(a,a,a)') TRIM(ADJUSTL(ikp_string)),'/Caapj/Ciipj'
          CALL getarr(fid, var_name, Ciipj_full) ! get array (moli format)
          ! Fill sub array with only usefull polynmials degree
          DO ip_i = 0,Pmaxi ! Loop over rows
          DO ij_i = 0,Jmaxi
                irow_sub  = (Jmaxi +1)*ip_i + ij_i +1
                irow_full = (jdimi +1)*ip_i + ij_i +1
                DO il_i = 0,Pmaxi ! Loop over columns
                DO ik_i = 0,Jmaxi
                      icol_sub  = (Jmaxi +1)*il_i + ik_i +1
                      icol_full = (jdimi +1)*il_i + ik_i +1
                      Ciipj__kp (irow_sub,icol_sub,ikp) = Ciipj_full (irow_full,icol_full)
                ENDDO
                ENDDO
          ENDDO
          ENDDO
          DEALLOCATE(Ciipj_full)

          IF(abs(CO) .NE. 3) THEN ! Pitch angle is only applied on like-species
            !!!!!!!!!!!!!!! E-I matrices !!!!!!!!!!!!!!
            ! Get test and field e-i collision matrices
            CALL allocate_array( CeipjT_full, 1,(pdime+1)*(jdime+1), 1,(pdime+1)*(jdime+1))
            CALL allocate_array( CeipjF_full, 1,(pdime+1)*(jdime+1), 1,(pdimi+1)*(jdimi+1))
            WRITE(var_name,'(a,a)') TRIM(ADJUSTL(ikp_string)),'/Ceipj/CeipjT'
            CALL getarr(fid, var_name, CeipjT_full)
            WRITE(var_name,'(a,a)') TRIM(ADJUSTL(ikp_string)),'/Ceipj/CeipjF'
            CALL getarr(fid, var_name, CeipjF_full)
            ! Fill sub array with only usefull polynmials degree
            DO ip_e = 0,pmaxe ! Loop over rows
            DO ij_e = 0,jmaxe
                  irow_sub  = (jmaxe +1)*ip_e + ij_e +1
                  irow_full = (jdime +1)*ip_e + ij_e +1
                  DO il_e = 0,pmaxe ! Loop over columns
                  DO ik_e = 0,jmaxe
                        icol_sub  = (jmaxe +1)*il_e + ik_e +1
                        icol_full = (jdime +1)*il_e + ik_e +1
                        CeipjT_kp(irow_sub,icol_sub,ikp) = CeipjT_full(irow_full,icol_full)
                  ENDDO
                  ENDDO
                  DO il_i = 0,pmaxi ! Loop over columns
                  DO ik_i = 0,jmaxi
                        icol_sub  = (Jmaxi +1)*il_i + ik_i +1
                        icol_full = (jdimi +1)*il_i + ik_i +1
                        CeipjF_kp(irow_sub,icol_sub,ikp) = CeipjF_full(irow_full,icol_full)
                  ENDDO
                  ENDDO
            ENDDO
            ENDDO
            DEALLOCATE(CeipjF_full)
            DEALLOCATE(CeipjT_full)

            !!!!!!!!!!!!!!! I-E matrices !!!!!!!!!!!!!!
            ! get the Test and Back field electron ion collision matrix
            CALL allocate_array( CiepjT_full, 1,(pdimi+1)*(jdimi+1), 1,(pdimi+1)*(jdimi+1))
            CALL allocate_array( CiepjF_full, 1,(pdimi+1)*(jdimi+1), 1,(pdime+1)*(jdime+1))
            WRITE(var_name,'(a,a,a)') TRIM(ADJUSTL(ikp_string)),'/Ciepj/CiepjT'
            CALL getarr(fid, var_name, CiepjT_full)
            WRITE(var_name,'(a,a,a)') TRIM(ADJUSTL(ikp_string)),'/Ciepj/CiepjF'
            CALL getarr(fid, var_name, CiepjF_full)
            ! Fill sub array with only usefull polynmials degree
            DO ip_i = 0,Pmaxi ! Loop over rows
            DO ij_i = 0,Jmaxi
                  irow_sub  = (Jmaxi +1)*ip_i + ij_i +1
                  irow_full = (jdimi +1)*ip_i + ij_i +1
                  DO il_i = 0,Pmaxi ! Loop over columns
                  DO ik_i = 0,Jmaxi
                        icol_sub  = (Jmaxi +1)*il_i + ik_i +1
                        icol_full = (jdimi +1)*il_i + ik_i +1
                        CiepjT_kp(irow_sub,icol_sub,ikp) = CiepjT_full(irow_full,icol_full)
                  ENDDO
                  ENDDO
                  DO il_e = 0,pmaxe ! Loop over columns
                  DO ik_e = 0,jmaxe
                        icol_sub  = (jmaxe +1)*il_e + ik_e +1
                        icol_full = (jdime +1)*il_e + ik_e +1
                        CiepjF_kp(irow_sub,icol_sub,ikp) = CiepjF_full(irow_full,icol_full)
                  ENDDO
                  ENDDO
            ENDDO
            ENDDO
            DEALLOCATE(CiepjF_full)
            DEALLOCATE(CiepjT_full)
          ELSE
            CeipjT_kp = 0._dp; CeipjF_kp = 0._dp; CiepjT_kp = 0._dp; CiepjF_kp = 0._dp;
          ENDIF
        ENDIF
      ENDDO
      CALL closef(fid)

      IF (CO .GT. 0) THEN ! Interpolation of the kperp matrix values on kx ky grid
        IF (my_id .EQ. 0 ) WRITE(*,*) '...Interpolation from matrices kperp to simulation kx,ky...'
        DO ikx = ikxs,ikxe
          DO iky = ikys,ikye
            kperp_sim = SQRT(kxarray(ikx)**2+kyarray(iky)**2) ! current simulation kperp

            ! Find the interval in kp grid mat where kperp_sim is contained
            ! Loop over the whole kp mat grid to find the smallest kperp that is
            ! larger than the current kperp_sim (brute force...)
            DO ikp=1,nkp_mat
              ikp_mat   = ikp ! the first indice of the interval (k0)
              kperp_mat = kp_grid_mat(ikp)
              IF(kperp_mat .GT. kperp_sim) EXIT ! a matrix with kperp2 > current kx2 + ky2 has been found
            ENDDO
            ! Interpolation
            ! interval boundaries
            ikp_next  = ikp_mat     !index of k1 (larger than kperp_sim thanks to previous loop)
            ikp_prev  = ikp_mat - 1 !index of k0 (smaller neighbour to interpolate inbetween)
            ! write(*,*) kp_grid_mat(ikp_prev), '<', kperp_sim, '<', kp_grid_mat(ikp_next)

            ! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0)
            zerotoone = (kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev))

            ! Linear interpolation between previous and next kperp matrix values
            Ceepj (:,:,ikx,iky) = (Ceepj__kp(:,:,ikp_next) - Ceepj__kp(:,:,ikp_prev))*zerotoone + Ceepj__kp(:,:,ikp_prev)
            CeipjT(:,:,ikx,iky) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev)
            CeipjF(:,:,ikx,iky) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev)
            Ciipj (:,:,ikx,iky) = (Ciipj__kp(:,:,ikp_next) - Ciipj__kp(:,:,ikp_prev))*zerotoone + Ciipj__kp(:,:,ikp_prev)
            CiepjT(:,:,ikx,iky) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev)
            CiepjF(:,:,ikx,iky) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev)
          ENDDO
        ENDDO
      ELSE ! DK -> No kperp dep, copy simply to final collision matrices
        Ceepj (:,:,1,1) = Ceepj__kp (:,:,1)
        CeipjT(:,:,1,1) = CeipjT_kp(:,:,1)
        CeipjF(:,:,1,1) = CeipjF_kp(:,:,1)
        Ciipj (:,:,1,1) = Ciipj__kp (:,:,1)
        CiepjT(:,:,1,1) = CiepjT_kp(:,:,1)
        CiepjF(:,:,1,1) = CiepjF_kp(:,:,1)
      ENDIF
      ! Deallocate auxiliary variables
      DEALLOCATE (Ceepj__kp ); DEALLOCATE (CeipjT_kp); DEALLOCATE (CeipjF_kp)
      DEALLOCATE (Ciipj__kp ); DEALLOCATE (CiepjT_kp); DEALLOCATE (CiepjF_kp)

      IF( CO_AA_ONLY ) THEN
        CeipjF = 0._dp;
        CeipjT = 0._dp;
        CiepjF = 0._dp;
        CiepjT = 0._dp;
      ENDIF

      IF (my_id .EQ. 0) WRITE(*,*) '============DONE==========='

    END SUBROUTINE load_COSOlver_mat
    !******************************************************************************!

end module collision
