moments_eq_rhs.F90

SUBROUTINE moments_eq_rhs

  USE basic
  USE time_integration
  USE array
  USE fields
  USE fourier_grid
  USE model

  use prec_const
  IMPLICIT NONE

  INTEGER     :: ip,ij, ikr,ikz ! loops indices
  REAL(dp)    :: ip_dp, ij_dp
  REAL(dp) :: kr, kz, kperp2
  REAL(dp) :: taue_qe_etaB, taui_qi_etaB
  REAL(dp) :: kernelj, kerneljp1, kerneljm1, b_e2, b_i2 ! Kernel functions and variable
  REAL(dp) :: factj, xb_e2, xb_i2 ! Auxiliary variables
  REAL(dp) :: xNapj, xNapp2j, xNapm2j, xNapjp1, xNapjm1 ! factors depending on the pj loop
  REAL(dp) :: xphij, xphijp1, xphijm1, xColl
  COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
  COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs

  !Precompute species dependant factors
  taue_qe_etaB  = tau_e/q_e * eta_B
  xb_e2         = sigma_e**2 * tau_e!/2.0 ! species dependant factor of the Kernel, squared
  taui_qi_etaB  = tau_i/q_i * eta_B
  xb_i2         = sigma_i**2 * tau_i!/2.0 ! species dependant factor of the Kernel, squared

  !!!!!!!!! Electrons moments RHS !!!!!!!!!
  Tphi          = 0 ! electrostatic potential term
  ploope : DO ip = ips_e, ipe_e ! This loop is from 1 to pmaxe+1
    ip_dp = REAL(ip-1.,dp) ! REAL index is one minus the loop index (0 to pmax)

    ! N_e^{p+2,j} multiplicator
    IF (ip+2 .LE. ipe_e) THEN
      xNapp2j = -taue_qe_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
    ELSE
      xNapp2j = 0.
    ENDIF

    ! N_e^{p-2,j} multiplicator
    IF (ip-2 .GE. ips_e) THEN
      xNapm2j = -taue_qe_etaB * SQRT(ip_dp*(ip_dp - 1.))
    ELSE
      xNapm2j = 0.
    ENDIF

    jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxe+1
      ij_dp = REAL(ij-1.,dp) ! REAL index is one minus the loop index (0 to jmax)

      ! N_e^{p,j+1} multiplicator
      IF (ij+1 .LE. ije_e) THEN
        xNapjp1 = taue_qe_etaB * (ij_dp + 1.)
      ELSE
        xNapjp1 = 0.
      ENDIF

      ! N_e^{p,j-1} multiplicator
      IF (ij-1 .GE. ijs_e) THEN
        xNapjm1 = taue_qe_etaB * ij_dp
      ELSE
        xNapjm1 = 0.
      ENDIF

      ! N_e^{pj} multiplicator
      xNapj   = -taue_qe_etaB * 2.*(ip_dp + ij_dp + 1.)

      ! Collision operator (DK Lenard-Bernstein basis)
      xColl = ip_dp + 2.*ij_dp
      ! ... adding Dougherty terms
      IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) ) THEN ! kronecker p0 * j1
        TColl01 = 2.0/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
        - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
        TColl20 = 0.0; TColl10 = 0.0;
      ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) ) THEN ! kronecker p2 * j0
        TColl20 = -SQRT2/3.0*(SQRT2*moments_e(3,1,ikr,ikz,updatetlevel)&
        - 2.0*moments_e(1,2,ikr,ikz,updatetlevel))
        TColl10 = 0.0; TColl01 = 0.0;
      ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
        TColl10 = moments_e(2,1,ikr,ikz,updatetlevel)
        TColl20 = 0.0; TColl01 = 0.0;
      ELSE
        TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
      ENDIF

      ! phi multiplicator for different kernel numbers
      IF (ip .EQ. 1) THEN !(kronecker delta_p^0)
        xphij   =  (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
        xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
        xphijm1 = -(eta_T - eta_B)* ij_dp
        factj   = REAL(Factorial(ij-1),dp)
      ELSE IF (ip .EQ. 3) THEN !(kronecker delta_p^2)
        xphij   =  (eta_T/SQRT2 - SQRT2*eta_B)
        factj   = REAL(Factorial(ij-1),dp)
      ELSE
        xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
        factj = 1
      ENDIF

      !write(*,*) '(ip,ij) = (', ip,',', ij,')'

      ! Loop on kspace
      krloope : DO ikr = ikrs,ikre
        kzloope : DO ikz = ikzs,ikze
          kr     = krarray(ikr)   ! Poloidal wavevector
          kz     = kzarray(ikz)   ! Toroidal wavevector
          kperp2 = kr**2 + kz**2  ! perpendicular wavevector
          b_e2   = kperp2 * xb_e2 ! Bessel argument

          !! Compute moments and mixing terms
          ! term propto N_e^{p,j}
          TNapj = moments_e(ip,ij,ikr,ikz,updatetlevel) * xNapj
          ! term propto N_e^{p+2,j}
          IF (ip+2 .LE. ipe_e) THEN
            TNapp2j = moments_e(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
          ELSE
            TNapp2j = 0.
          ENDIF
          ! term propto N_e^{p-2,j}
          IF (ip-2 .GE. ips_e) THEN
            TNapm2j = moments_e(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
          ELSE
            TNapm2j = 0.
          ENDIF
          ! xterm propto N_e^{p,j+1}
          IF (ij+1 .LE. ije_e) THEN
            TNapjp1 = moments_e(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
          ELSE
            TNapjp1 = 0.
          ENDIF
          ! term propto N_e^{p,j-1}
          IF (ij-1 .GE. ijs_e) THEN
            TNapjm1 = moments_e(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
          ELSE
            TNapjm1 = 0.
          ENDIF

          ! Collision term completed (DK Dougherty)
          TColl = -nu * (xColl * moments_e(ip,ij,ikr,ikz,updatetlevel) &
                          + TColl01 + TColl10 + TColl20)

          !! Electrical potential term
          Tphi = 0
          IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0)
            kernelj    = b_e2**(ij-1) * exp(-b_e2)/factj
            kerneljp1  = kernelj * b_e2  /(ij_dp + 1.)
            kerneljm1  = kernelj * ij_dp / b_e2
            Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
          ENDIF

          ! Sum of all precomputed terms
          moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
              imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl

        END DO kzloope
      END DO krloope

    END DO jloope
  END DO ploope

  !!!!!!!!! Ions moments RHS !!!!!!!!!
  Tphi = 0 ! electrostatic potential term

  ploopi : DO ip = ips_i, ipe_i
    ip_dp = REAL(ip-1.,dp)

    ! x N_i^{p+2,j}
    IF (ip+2 .LE. ipe_i) THEN
      xNapp2j = -taui_qi_etaB * SQRT((ip_dp + 1.) * (ip_dp + 2.))
    ELSE
      xNapp2j = 0.
    ENDIF

    ! x N_i^{p-2,j}
    IF (ip-2 .GE. ips_i) THEN
      xNapm2j = -taui_qi_etaB * SQRT(ip_dp * (ip_dp - 1.))
    ELSE
      xNapm2j = 0.
    ENDIF

    jloopi : DO ij = ijs_i, ije_i
      ij_dp = REAL(ij-1.,dp)

      ! x N_i^{p,j+1}
      IF (ij+1 .LE. ije_i) THEN
        xNapjp1 = taui_qi_etaB * (ij_dp + 1.)
      ELSE
        xNapjp1 = 0.
      ENDIF

      ! x N_i^{p,j-1}
      IF (ij-1 .GE. ijs_i) THEN
        xNapjm1 = taui_qi_etaB * ij_dp
      ELSE
        xNapjm1 = 0.
      ENDIF

      ! x N_i^{pj}
      xNapj   = -taui_qi_etaB * 2.*(ip_dp + ij_dp + 1.)

      ! Collision term completed (DK Dougherty)
      TColl = -nu * (xColl * moments_i(ip,ij,ikr,ikz,updatetlevel) &
                      + TColl01 + TColl10 + TColl20)

      ! ... adding Dougherty terms
      IF ( (ip .EQ. 1) .AND. (ij .EQ. 2) ) THEN ! kronecker p0 * j1
        TColl01 = 2.0/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
        - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
        TColl20 = 0.0; TColl10 = 0.0;
      ELSEIF ( (ip .EQ. 3) .AND. (ij .EQ. 1) ) THEN ! kronecker p2 * j0
        TColl20 = -SQRT2/3.0*(SQRT2*moments_i(3,1,ikr,ikz,updatetlevel)&
        - 2.0*moments_i(1,2,ikr,ikz,updatetlevel))
        TColl10 = 0.0; TColl01 = 0.0;
      ELSEIF ( (ip .EQ. 2) .AND. (ij .EQ. 1) ) THEN ! kronecker p1 * j0
        TColl10 = moments_i(2,1,ikr,ikz,updatetlevel)
        TColl20 = 0.0; TColl01 = 0.0;
      ELSE
        TColl10 = 0.0; TColl20 = 0.0; TColl01 = 0.0;
      ENDIF

      ! x phi
      IF (ip .EQ. 1) THEN !(krokecker delta_p^0)
        xphij   =  (eta_n + 2.*ij_dp*eta_T - 2.*eta_B*(ij_dp+1.) )
        xphijp1 = -(eta_T - eta_B)*(ij_dp+1.)
        xphijm1 = -(eta_T - eta_B)* ij_dp
        factj   = REAL(Factorial(ij-1),dp)
      ELSE IF (ip .EQ. 3) THEN !(krokecker delta_p^2)
        xphij   =  (eta_T/SQRT2 - SQRT2*eta_B)
        factj   = REAL(Factorial(ij-1),dp)
      ELSE
        xphij = 0.; xphijp1 = 0.; xphijm1 = 0.
        factj = 1.
      ENDIF

      ! Loop on kspace
      krloopi : DO ikr = ikrs,ikre
        kzloopi : DO ikz = ikzs,ikze
          kr     = krarray(ikr)   ! Poloidal wavevector
          kz     = kzarray(ikz)   ! Toroidal wavevector
          kperp2 = kr**2 + kz**2  ! perpendicular wavevector
          b_i2   = kperp2 * xb_i2 ! Bessel argument

          !! Compute moments and mixing terms
          ! term propto N_i^{p,j}
          TNapj = moments_i(ip,ij,ikr,ikz,updatetlevel) * xNapj
          ! term propto N_i^{p+2,j}
          IF (ip+2 .LE. ipe_i) THEN
            TNapp2j = moments_i(ip+2,ij,ikr,ikz,updatetlevel) * xNapp2j
          ELSE
            TNapp2j = 0.
          ENDIF
          ! term propto N_i^{p-2,j}
          IF (ip-2 .GE. ips_i) THEN
            TNapm2j = moments_i(ip-2,ij,ikr,ikz,updatetlevel) * xNapm2j
          ELSE
            TNapm2j = 0.
          ENDIF
          ! xterm propto N_i^{p,j+1}
          IF (ij+1 .LE. ije_i) THEN
            TNapjp1 = moments_i(ip,ij+1,ikr,ikz,updatetlevel) * xNapjp1
          ELSE
            TNapjp1 = 0.
          ENDIF
          ! term propto N_i^{p,j-1}
          IF (ij-1 .GE. ijs_i) THEN
            TNapjm1 = moments_i(ip,ij-1,ikr,ikz,updatetlevel) * xNapjm1
          ELSE
            TNapjm1 = 0.
          ENDIF

          ! Collision term completed (Dougherty)
          TColl = -nu* (xColl * moments_i(ip,ij,ikr,ikz,updatetlevel))

          !! Electrical potential term
          Tphi = 0
          IF ( (ip .eq. 1) .or. (ip .eq. 3) ) THEN ! 0 otherwise (krokecker delta_p^0, delta_p^2)
            kernelj    = b_i2**(ij-1) * exp(-b_i2)/factj
            kerneljp1  = kernelj * b_i2  /(ij_dp + 1.)
            kerneljm1  = kernelj * ij_dp / b_i2
            Tphi = (xphij*Kernelj + xphijp1*Kerneljp1 + xphijm1*Kerneljm1) * phi(ikr,ikz)
          ENDIF

          ! Sum of all precomputed terms
          moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
              imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 + Tphi) + TColl

        END DO kzloopi
      END DO krloopi

    END DO jloopi
  END DO ploopi

END SUBROUTINE moments_eq_rhs