diff --git a/src/improved_sugama_mod.f90 b/src/improved_sugama_mod.f90
index 9c59be6..9bc87f6 100644
--- a/src/improved_sugama_mod.f90
+++ b/src/improved_sugama_mod.f90
@@ -1,683 +1,690 @@
MODULE ImprovedSugama
!
! Linearized DK/GK Improved Sugamam Collision operator between species a and b
!
! !! TODO(Sam): Corriger ces annotations
! Note: as defined below, the Sugama operator differs from a 2 factor from the original paper due to the definition of the isothermal part of the linearized Coulomb
!
! Ref: Sugama et. al., Phsyics of Plasmas 26 (2019)
USE prec_const
USE basic
USE model
USE FMZM
USE coeff
USE basis_transformation
USE speed_functions
USE auxval
use coulomb_parts
use sugama
IMPLICIT NONE
PUBLIC
CONTAINS
! TESTS
SUBROUTINE test_SugamaRelations()
IMPLICIT NONE
! Local vars
INTEGER :: l, k
TYPE(FM), SAVE :: MabLlk_, MabSlk_
TYPE(FM), SAVE :: NabLlk_,NabSlk_
TYPE(FM), SAVE :: DeltaMablk_, DeltaNablk_
TYPE(FM), SAVE :: val
INTEGER, PARAMETER :: lmax = 5, kmax = 5
!INTEGER, PARAMETER :: lmax = 0, kmax = 0
TYPE(FM), DIMENSION(0:lmax, 0:kmax) :: Mab, Mba, Nab, Nba, MabS, MbaS, NabS, NbaS, DeltaMab, DeltaNab
CALL FM_ENTER_USER_ROUTINE
lloop_ab: DO l=0, lmax
kloop_ab: DO k=0, kmax
Mab(l,k) = MabLlk(l,k)
Nab(l,k) = NabLlk(l,k)
MabS(l,k) = MabSlk(l,k)
NabS(l,k) = NabSlk(l,k)
DeltaMab(l,k) = DeltaMlk(l,k)
DeltaNab(l,k) = DeltaNlk(l,k)
END DO kloop_ab
END DO lloop_ab
CALL invert_basic
CALL invert_model
CALL evaluate_auxval_erf
lloop_ba: DO l=0, lmax
kloop_ba: DO k=0, kmax
Mba(l,k) = MabLlk(l,k)
Nba(l,k) = NabLlk(l,k)
MbaS(l,k) = MabSlk(l,k)
NbaS(l,k) = NabSlk(l,k)
END DO kloop_ba
END DO lloop_ba
CALL invert_model
CALL invert_basic
CALL evaluate_auxval_erf
val = TO_FM(0)
WRITE(*,*) '====== Tests Sugama Relations ======='
WRITE(*,*) "## Landau"
IF( tau .eq. 1.0 ) THEN
WRITE(*,*) "# Temperature Ta = Tb"
lloop: DO l=0, lmax
kloop: DO k=0, kmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', l, lmax, k, kmax
PRINT '("Mab_"I1"_"I1" - Mab_"I1"_"I1" = "E11.4)', l, k, k, l, TO_DP(Mab(l,k) - Mab(k,l))
PRINT '(">Nab_"I1"_"I1" = "E17.10)', l, k, TO_DP(Nab(l,k))
PRINT '(">Nba_"I1"_"I1" = "E17.10)', k, l, TO_DP(Nba(k,l))
PRINT '("Nab_"I1"_"I1" - tau*chi*Nba_"I1"_"I1" = "E11.4)', l, k, k, l, TO_DP(Nab(l,k) - taufm_*chifm_*Nba(k,l))
END DO kloop
END DO lloop
END IF
WRITE(*,*) "## Sugama"
IF( tau .eq. 1.0 ) THEN
WRITE(*,*) "# Temperature Ta = Tb"
lloop2: DO l=0, lmax
kloop2: DO k=0, kmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', l, lmax, k, kmax
PRINT '(">NSab_"I1"_"I1" = "E17.10)', l, k, TO_DP(NabS(l,k))
PRINT '("MSab_"I1"_"I1" - Mab_"I1"_"I1" = "E11.4)', l, k, l, k, TO_DP(MabS(l,k) - Mab(k,l))
PRINT '("NSab_"I1"_"I1" - Nab_"I1"_0*Nab_0_"I1"/Nab_0_0 = "E11.4)', l, k, l, k, TO_DP(NabS(l,k) - Nab(l,0)*Nab(0,k)/Nab(0,0))
END DO kloop2
END DO lloop2
lloop2_1: DO l=0, lmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', l, lmax, 0, kmax
PRINT '("NSab_"I1"_0 - Nab_"I1"_0 = "E11.4)', l, l, TO_DP(NabS(l,0) - Nab(l,0))
END DO lloop2_1
lloop2_2: DO k=0, kmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', 0, lmax, k, kmax
PRINT '("NSab_0_"I1" - Nab_0_"I1" = "E11.4)', k, k, TO_DP(NabS(0,k) - Nab(0,k))
END DO lloop2_2
END IF
WRITE(*,*) "# Temperature Ta != Tb"
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', 0, lmax, 0, kmax
PRINT '("Mab_0_0 + Nab_0_0 ="E11.4)', TO_DP(Mab(0,0) + Nab(0,0))
PRINT '("Mab_0_0 - MSab_0_0 ="E11.4)', TO_DP(Mab(0,0) - MabS(0,0))
PRINT '("MSab_0_0 + NSab_0_0 ="E11.4)', TO_DP(MabS(0,0) + NabS(0,0))
kloop2_3: DO k=0, kmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', 0, lmax, k, kmax
PRINT '("MSab_0_"I1" + tau*chi*NSba_0_"I1" = "E11.4)', k, k, TO_DP(MabS(0,k) + taufm_*chifm_*NbaS(0,k))
PRINT '("Mab_0_"I1" + tau*chi*Nba_0_"I1" = "E11.4)', k, k, TO_DP(Mab(0,k) + taufm_*chifm_*Nba(0,k))
END DO kloop2_3
lloop4: DO l=0, lmax
kloop4: DO k=0, kmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', l, lmax, k, kmax
PRINT '("MSab_"I1"_"I1" - MSab_"I1"_"I1" = "E11.4)', l, k, k, l, TO_DP(MabS(l,k) - MabS(k,l))
PRINT '("NSab_"I1"_"I1" - tau*chi*NSba_"I1"_"I1" = "E11.4)', l, k, k, l, TO_DP(NabS(l,k) - taufm_*chifm_*NbaS(k,l))
PRINT '("NSab_"I1"_"I1" - MSab_0_"I1"*NSab_0_"I1"/Mab_0_0 = "E11.4)', l,k, l, k, TO_DP(NabS(l,k) - MabS(0,l)*NabS(0,k)/Mab(0,0))
PRINT '("NSab_"I1"_"I1" - NSab_"I1"_0*NSab_0_"I1"/Nab_0_0 = "E11.4)', l,k, l, k, TO_DP(NabS(l,k) - NabS(l,0)*NabS(0,k)/Nab(0,0))
END DO kloop4
END DO lloop4
WRITE(*,*) "## Sugama+"
IF( tau .eq. 1.0 ) THEN
WRITE(*,*) "# Temperature Ta = Tb"
lloop3: DO l=0, lmax
kloop3: DO k=0, kmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', l, lmax, k, kmax
PRINT '("DeltaMab_"I1"_"I1" = "E11.4)', l, k, TO_DP(DeltaMab(l,k))
PRINT '("DeltaNab_"I1"_"I1" - (Nab_0_0*Nab_"I1"_"I1"-Nab_"I1"_0*Nab_0_"I1")/Nab_0_0 = "E11.4)', l, k, l, k, l, k, TO_DP(DeltaNab(l,k) - (Nab(0,0)*Nab(l,k)-Nab(l,0)*Nab(0,k))/Nab(0,0))
END DO kloop3
END DO lloop3
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', 0, lmax, 0, kmax
PRINT '("DeltaNab_0_0 ="E11.4)', TO_DP(DeltaNab(0,0))
lloop3_1: DO l=0, lmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', l, lmax, 0, kmax
PRINT '("DeltaNab_"I1"_0 ="E11.4)', l, TO_DP(DeltaNab(l,0))
END DO lloop3_1
lloop3_2: DO k=0, kmax
PRINT '("Loop: "I1"/"I1", "I1"/"I1)', 0, lmax, k, kmax
PRINT '("DeltaNab_0_"I1" ="E11.4)', k, TO_DP(DeltaNab(0,k))
END DO lloop3_2
END IF
WRITE(*,*) '====== End Tests ======='
CALL FM_EXIT_USER_ROUTINE
END SUBROUTINE test_SugamaRelations
!_________________________________________________________________
!
! DRIFT-KINETIC TEST PART
!_________________________________________________________________
SUBROUTINE ImprovedSugamaTDK(l,k,rCabT)
!
! Add the test correction matrix to DK Sugama. The improved Sugama collision
! operator reproduces the same friction-flow relations as the lineared DK full
! Coulomb collision operator [see Eqs. (12) and (13) in Ref. [1]]
!
! Ref[1]: Sugama et al, Phys. Plasmas 26, 102108 (2019)
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabT
!
!
! local vars.
INTEGER :: j,i,g,h,icol
TYPE(FM), SAVE :: NX0,NX1,T4,T4inv,NX2,cj_,ci_
TYPE(FM), SAVE :: DeltaMablk_
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
! construct DK Test Sugama
CALL SugamaTDK(l,k, rCabT)
! add correction matrix elements \Delta C_{ab}^{FLS}
IF( l .ge. 1 .or. k .ge. 1) THEN
!
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
jloop: DO j=0, k + floor(real(l)/2)
NX1 = 4*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j))/3/SQRT(PIFM)
cj_ = CJ(j)
T4str_ = GET_T4(1,j,l,k) ! ... T^{lk}_{gh}
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(j))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*j +1)/2)*TO_FM(T4str_)
IF( T4inv .ne. TO_FM(0)) THEN
iloop: DO i=0,imaxx ! .... energy moment
ci_ = CJ(i)
! Complute test correction matrix
DeltaMablk_= DeltaMlk(j,i) ! Delta M_{ab}^{ji}
gloop: DO g =0,min(2*i+1,Pmaxa)
NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0,min(i,Jmaxa)
T4str_ = GET_T4(1,i,g,h) ! ... T^{lk}_{gh}
T4 = TO_FM(T4str_)
icol = numba(g,h)
rCabT(icol) = rCabT(icol) + NX0*NX1*NX2*T4inv*T4*cj_*ci_*DeltaMablk_*TO_FM(4)/TO_FM(3)/SQRT(PIFM)
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE ImprovedSugamaTDK
!
!_________________________________________________________________
!
! DRIFT-KINETIC FIELD PART
!_________________________________________________________________
!
SUBROUTINE ImprovedSugamaFDK(l,k,rCabF)
! Add field correction matrix to DK Sugama. The improved Sugama collision operator
! reproduces the same friction-flow relations as the lineared DK full Coulomb collision
! operator [see Eqs. (12) and (13) in Ref. [1]]
!
! Ref[1]: Sugama et al, Phys. Plasmas 26, 102108 (2019)
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabF
!
!
! local vars.
INTEGER :: j,i,g,h,icol
TYPE(FM), SAVE :: NX0,NX1,T4,T4inv,NX2,cj_,ci_
TYPE(FM), SAVE :: DeltaNablk_
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
! construct Sugama Field part
CALL SugamaFDK(l,k, rCabF)
!
! add correction matrix elements \Delta C_{ab}^{FLS}
IF( l .ge. 1 .or. k .ge. 1) THEN
!
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
jloop: DO j=0, k + floor(real(l)/2)
NX1 = 4*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j))/3/SQRT(PIFM)
cj_ = CJ(j)
T4str_ = GET_T4(1,j,l,k) ! ... T^{lk}_{gh}
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(j))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*j +1)/2)*TO_FM(T4str_)
IF( T4inv .ne. TO_FM(0)) THEN
iloop: DO i=0,imaxx ! ... energy moment
ci_ = CJ(i)
! Complute test correction matrix
DeltaNablk_= DeltaNlk(j,i) ! Delta N_{ab}^{ji}
gloop: DO g =0,min(2*i +1,Pmaxb)
NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0,min(i,Jmaxb)
T4str_ = GET_T4(1,i,g,h) ! ... T^{lk}_{gh}
T4 = TO_FM(T4str_)
icol = numbb(g,h)
rCabF(icol) = rCabF(icol) + NX0*NX1*NX2*T4inv*T4*cj_*ci_*DeltaNablk_*TO_FM(4)/TO_FM(3)/SQRT(PIFM)
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE ImprovedSugamaFDK
!
!_________________________________________________________________
!
! DRIFT-KINETIC SELF
!_________________________________________________________________
!
SUBROUTINE ImprovedSugamaSelfDK(l,k,rCaa)
! Add correction matrix to self DK Sugama. The improved Sugama collision operator
! reproduces the same friction-flow relations as the lineared DK full Coulomb collision
! operator [see Eqs. (12) and (13) in Ref. [1]]
!
! Ref[1]: Sugama et al, Phys. Plasmas 26, 102108 (2019)
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCaa
!
!
! local vars.
INTEGER :: j,i,g,h,icol
TYPE(FM), SAVE :: NX0,NX1,T4,T4inv,NX2,cj_,ci_
TYPE(FM), SAVE :: DeltaMaalk_,DeltaNaalk_
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
- !
- ! Construct self DK Sugama
- CALL SugamaSelfDK(l,k,rCaa)
-
- !
- ! add correction matrix elements \Delta C_{aa}^{FLS}
- IF( l .ge. 1 .or. k .ge. 1) THEN
- !
- NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
- !
- !
- jloop: DO j=0, k + floor(real(l)/2)
- NX1 = 4*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j))/3/SQRT(PIFM)
- cj_ = CJ(j)
- T4str_ = GET_T4(1,j,l,k) ! ... T^{lk}_{gh}
- T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(j))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*j +1)/2)*TO_FM(T4str_)
- IF( T4inv .ne. TO_FM(0)) THEN
-
- iloop: DO i=0,imaxx
- ci_ = CJ(i)
-
- ! Complute test and field correction matrix
- DeltaNaalk_= DeltaNlk(j,i)
- DeltaMaalk_= DeltaMlk(j,i)
- gloop: DO g =0,min(2*i +1,Pmaxa)
- NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
- hloop: DO h=0,min(i,Jmaxa)
- T4str_ = GET_T4(1,i,g,h) ! ... T^{lk}_{gh}
- T4 = TO_FM(T4str_)
- icol = numba(g,h)
- rCaa(icol) = rCaa(icol) + NX0*NX1*NX2*T4inv*T4*cj_*ci_*(DeltaNaalk_ + DeltaMaalk_)*TO_FM(4)/TO_FM(3)/SQRT(PIFM)
- ENDDO hloop
- ENDDO gloop
- ENDDO iloop
- ENDIF
- ENDDO jloop
- ENDIF
+ CALL ImprovedSugamaTDK(l,k,rCaa)
+ CALL ImprovedSugamaFDK(l,k,rCaa)
+
+ ! !
+ ! ! Construct self DK Sugama
+ ! CALL SugamaSelfDK(l,k,rCaa)
+
+ ! !
+ ! ! add correction matrix elements \Delta C_{aa}^{FLS}
+ ! IF( l .ge. 1 .or. k .ge. 1) THEN
+ ! !
+ ! NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
+ ! !
+ ! !
+ ! jloop: DO j=0, k + floor(real(l)/2)
+ ! NX1 = 4*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j))/3/SQRT(PIFM)
+ ! cj_ = CJ(j)
+ ! T4str_ = GET_T4(1,j,l,k) ! ... T^{lk}_{gh}
+ ! T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(j))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*j +1)/2)*TO_FM(T4str_)
+ ! IF( T4inv .ne. TO_FM(0)) THEN
+
+ ! iloop: DO i=0,imaxx
+ ! ci_ = CJ(i)
+
+ ! ! Complute test and field correction matrix
+ ! DeltaNaalk_= DeltaNlk(j,i)
+ ! DeltaMaalk_= DeltaMlk(j,i)
+
+ ! gloop: DO g =0,min(2*i +1,Pmaxa)
+ ! NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
+ ! hloop: DO h=0,min(i,Jmaxa)
+ ! T4str_ = GET_T4(1,i,g,h) ! ... T^{lk}_{gh}
+ ! T4 = TO_FM(T4str_)
+ ! icol = numba(g,h)
+ ! rCaa(icol) = rCaa(icol) + NX0*NX1*NX2*T4inv*T4*cj_*ci_*(DeltaNaalk_ + DeltaMaalk_)*TO_FM(4)/TO_FM(3)/SQRT(PIFM)
+ ! ENDDO hloop
+ ! ENDDO gloop
+ ! ENDDO iloop
+ ! ENDIF
+ ! ENDDO jloop
+ ! ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE ImprovedSugamaSelfDK
!
!_________________________________________________________________
!
! GYROKINETIC TEST PART
!_________________________________________________________________
!
SUBROUTINE ImprovedSugamaTGK(l,k,rCabT)
!
! Add the test correction matrix to GK Sugama. The improved Sugama collision
! operator reproduces the same friction-flow relations as the lineared GK full
! Coulomb collision operator [see Eqs. (12) and (13) in Ref. [1]]
!
! Ref[1]: Sugama et al, Phys. Plasmas 26, 102108 (2019)
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),dimension(:),INTENT(INOUT) :: rCabT
!
! Local variables
INTEGER :: n, t, p, j, icol
TYPE(FM), SAVE :: NX0, cn_, ct_
TYPE(FM), SAVE :: DeltaMabnt_, IPara_a_lkn_, IPerp_a_lkn_
TYPE(FM), SAVE :: IPara_a_pjt_, IPerp_a_pjt_
!
CALL FM_ENTER_USER_ROUTINE
! Construct GK Test Sugama
CALL SugamaTGK(l,k,rCabT)
! Add correction matrix \Delta C_{ab}^{TLS}
NX0 = TO_FM(8)/TO_FM(3)/SQRT(PIFM)
nloop: DO n=0, ImpSugamaJmaxa
IPara_a_lkn_ = CurlIparapjk(l,k,n,bafm_)
IPerp_a_lkn_ = CurlIperppjk(l,k,n,bafm_)
cn_ = CJ(n)
tloop: DO t=0, ImpSugamaKmaxa
DeltaMabnt_ = DeltaMlk(n,t)
ct_ = CJ(t)
ploop: DO p=0, Pmaxa
jloop: DO j=0, Jmaxa
IPara_a_pjt_ = CurlIparapjk(p,j,t,bafm_)
IPerp_a_pjt_ = CurlIperppjk(p,j,t,bafm_)
icol = numba(p,j)
rCabT(icol) = rCabT(icol) + NX0*cn_*ct_*DeltaMabnt_*(IPara_a_lkn_*IPara_a_pjt_ + IPerp_a_lkn_*IPerp_a_pjt_)
ENDDO jloop
ENDDO ploop
ENDDO tloop
ENDDO nloop
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE ImprovedSugamaTGK
!
!_________________________________________________________________
!
! GYROKINETIC FIELD PART
!_________________________________________________________________
!
SUBROUTINE ImprovedSugamaFGK(l,k,rCabF)
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),dimension(:),INTENT(INOUT) :: rCabF
!
! Local variables
INTEGER :: n, t, p, j, icol
TYPE(FM), SAVE :: NX0, cn_, ct_
TYPE(FM), SAVE :: DeltaNabnt_, IPara_a_lkn_, IPerp_a_lkn_
TYPE(FM), SAVE :: IPara_b_pjt_, IPerp_b_pjt_
!
CALL FM_ENTER_USER_ROUTINE
! Construct GK Test Sugama
CALL SugamaTGK(l,k,rCabF)
! Add correction matrix \Delta C_{ab}^{FLS}
NX0 = TO_FM(8)/TO_FM(3)/SQRT(PIFM)
nloop: DO n=0, ImpSugamaJmaxa
IPara_a_lkn_ = CurlIparapjk(l,k,n,bafm_)
IPerp_a_lkn_ = CurlIperppjk(l,k,n,bafm_)
cn_ = CJ(n)
tloop: DO t=0, ImpSugamaKmaxa
DeltaNabnt_ = DeltaNlk(n,t)
ct_ = CJ(t)
ploop: DO p=0, Pmaxa
jloop: DO j=0, Jmaxa
IPara_b_pjt_ = CurlIparapjk(p,j,t,bbfm_)
IPerp_b_pjt_ = CurlIperppjk(p,j,t,bbfm_)
- icol = numba(p,j)
+ icol = numbb(p,j)
rCabF(icol) = rCabF(icol) + NX0*cn_*ct_*DeltaNabnt_*(IPara_a_lkn_*IPara_b_pjt_ + IPerp_a_lkn_*IPerp_b_pjt_)/chifm_
ENDDO jloop
ENDDO ploop
ENDDO tloop
ENDDO nloop
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE ImprovedSugamaFGK
!
!_________________________________________________________________
!
! GYROKINETIC SELF
!_________________________________________________________________
!
SUBROUTINE ImprovedSugamaSelfGK(l,k,rCaa)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),dimension(:),INTENT(INOUT) :: rCaa
!
! Local variables
INTEGER :: n, t, p, j, icol
TYPE(FM), SAVE :: NX0, cn_, ct_
TYPE(FM), SAVE :: DeltaMabnt_, DeltaNabnt_
TYPE(FM), SAVE :: IPara_a_lkn_, IPerp_a_lkn_
TYPE(FM), SAVE :: IPara_a_pjt_, IPerp_a_pjt_
TYPE(FM), SAVE :: IPara_b_pjt_, IPerp_b_pjt_
!
CALL FM_ENTER_USER_ROUTINE
- ! Construct GK Self Sugama
- CALL SugamaSelfGK(l,k,rCaa)
-
- ! Add correction matrix \Delta C_{ab}^{SelfLS}
- NX0 = TO_FM(8)/TO_FM(3)/SQRT(PIFM)
- nloop: DO n=0, ImpSugamaJmaxa
- IPara_a_lkn_ = CurlIparapjk(l,k,n,bafm_)
- IPerp_a_lkn_ = CurlIperppjk(l,k,n,bafm_)
- cn_ = CJ(n)
- tloop: DO t=0, ImpSugamaKmaxa
- DeltaMabnt_ = DeltaMlk(n,t)
- DeltaNabnt_ = DeltaNlk(n,t)
- ct_ = CJ(t)
- ploop: DO p=0, Pmaxa
- jloop: DO j=0, Jmaxa
+ CALL ImprovedSugamaTGK(l,k,rCaa)
+ CALL ImprovedSugamaFGK(l,k,rCaa)
+
+ ! ! Construct GK Self Sugama
+ ! CALL SugamaSelfGK(l,k,rCaa)
+
+ ! ! Add correction matrix \Delta C_{ab}^{SelfLS}
+ ! NX0 = TO_FM(8)/TO_FM(3)/SQRT(PIFM)
+ ! nloop: DO n=0, ImpSugamaJmaxa
+ ! IPara_a_lkn_ = CurlIparapjk(l,k,n,bafm_)
+ ! IPerp_a_lkn_ = CurlIperppjk(l,k,n,bafm_)
+ ! cn_ = CJ(n)
+ ! tloop: DO t=0, ImpSugamaKmaxa
+ ! DeltaMabnt_ = DeltaMlk(n,t)
+ ! DeltaNabnt_ = DeltaNlk(n,t)
+ ! ct_ = CJ(t)
+ ! ploop: DO p=0, Pmaxa
+ ! jloop: DO j=0, Jmaxa
- PRINT *, 'n=', n, ' t=', t, ' p=', p, ' j=', j
+ ! PRINT *, 'n=', n, ' t=', t, ' p=', p, ' j=', j
- IPara_a_pjt_ = CurlIparapjk(p,j,t,bafm_)
- IPerp_a_pjt_ = CurlIperppjk(p,j,t,bafm_)
- IPara_b_pjt_ = CurlIparapjk(p,j,t,bbfm_)
- IPerp_b_pjt_ = CurlIperppjk(p,j,t,bbfm_)
+ ! IPara_a_pjt_ = CurlIparapjk(p,j,t,bafm_)
+ ! IPerp_a_pjt_ = CurlIperppjk(p,j,t,bafm_)
+ ! IPara_b_pjt_ = CurlIparapjk(p,j,t,bbfm_)
+ ! IPerp_b_pjt_ = CurlIperppjk(p,j,t,bbfm_)
- icol = numba(p,j)
- rCaa(icol) = rCaa(icol) + NX0*cn_*ct_*DeltaMabnt_*(IPara_a_lkn_*IPara_a_pjt_ + IPerp_a_lkn_*IPerp_a_pjt_)
- rCaa(icol) = rCaa(icol) + NX0*cn_*ct_*DeltaNabnt_*(IPara_a_lkn_*IPara_b_pjt_ + IPerp_a_lkn_*IPerp_b_pjt_)/chifm_
-
- ENDDO jloop
- ENDDO ploop
- ENDDO tloop
- ENDDO nloop
+ ! icol = numba(p,j)
+ ! rCaa(icol) = rCaa(icol) + NX0*cn_*ct_*DeltaMabnt_*(IPara_a_lkn_*IPara_a_pjt_ + IPerp_a_lkn_*IPerp_a_pjt_)
+ ! rCaa(icol) = rCaa(icol) + NX0*cn_*ct_*DeltaNabnt_*(IPara_a_lkn_*IPara_b_pjt_ + IPerp_a_lkn_*IPerp_b_pjt_)/chifm_
+
+ ! ENDDO jloop
+ ! ENDDO ploop
+ ! ENDDO tloop
+ ! ENDDO nloop
CALL FM_EXIT_USER_ROUTINE
END SUBROUTINE ImprovedSugamaSelfGK
!
!_______________________________________________________________
!
! Usefull functions
!_______________________________________________________________
!
FUNCTION CurlIparapjk(p,j,k, bs) RESULT(XOUT)
!
! Compute \mathcal I_{\paralell s}^{pjk}
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: p, j, k
TYPE(FM), INTENT(IN) :: bs
TYPE(FM) :: XOUT
! Local vars
INTEGER :: n, s
TYPE(FM), SAVE :: NX0, T4inv, kern_, d0njs
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM(0)
NX0 = TO_FM(2)/TO_FM(3)/SQRT(PIFM)*FACTORIAL(TO_FM(2*k+1)/2)/FACTORIAL(TO_FM(k))/SQRT(TO_FM(2)**p * FACTORIAL(p))
nloop: DO n=0, naFLR
kern_ = kerneln(bs, n)
sloop: DO s=0, n+j
IF( (p .ge. 1 .or. s .ge. 1) .and. (s+p/2 .ge. k) ) THEN
T4str_ = GET_T4(1,k,p,s)
T4inv = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(k))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*k +1)/2)*TO_FM(T4str_)
d0njs = ALLX2L(n,j,s,0)
XOUT = XOUT+NX0*T4inv*kern_*d0njs
END IF
END DO sloop
END DO nloop
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION CurlIparapjk
FUNCTION CurlIperppjk(p,j,k,bs) RESULT(XOUT)
!
! Compute \mathcal I_{\perp s}^{pjk}
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: p, j, k
TYPE(FM), INTENT(IN) :: bs
TYPE(FM) :: XOUT
! Local vars
INTEGER :: n, s, r, q
TYPE(FM), SAVE :: NX0, NX1, T4inv, kern_, d1njs, Lkq
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
XOUT = TO_FM(0)
NX0 = TO_FM(1)/SQRT(PIFM)/SQRT(TO_FM(2)**p * FACTORIAL(p))*bs
nloop: DO n=0, naFLR
kern_ = kerneln(bs, n)
sloop: DO s=0, n+j+1
d1njs = ALLX2L(n,j,s,1)
rloop: DO r=0, s+p/2
T4str_ = GET_T4(0,r,p,s)
T4inv = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(k))*TO_FM(1)/2/FACTORIAL(TO_FM(2*k +1)/2)*TO_FM(T4str_)
qloop: DO q=0, k
Lkq=Lpjl(real(1,dp),real(k,dp),real(q,dp))
NX1 = FACTORIAL(TO_FM(2*1+1)/2)/(TO_FM(n+1))*GENERALIZED_BINOMIAL(TO_FM(r-q-1), TO_FM(r))
XOUT = XOUT + NX0*NX1*kern_*T4inv*d1njs*Lkq
END DO qloop
END DO rloop
END DO sloop
END DO nloop
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION CurlIperppjk
FUNCTION GENERALIZED_BINOMIAL(x,y) RESULT(XOUT)
IMPLICIT NONE
!
TYPE(FM), INTENT(IN) :: x, y
TYPE(FM) :: XOUT
! Local vars
TYPE(FM), SAVE :: x_, y_, sign_
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
x_ = x
y_ = y
sign_ = TO_FM(1)
IF( x .lt. 0 .or. y .lt. 0) THEN
IF(y .ge. 0) THEN
sign_ = TO_FM(-1)**y
x_ = -x + y - 1
y_ = y
ELSEIF(y .le. x) THEN
sign_ = TO_FM(-1)**(x-y)
x_ = -y - 1
y_ = x - y
ELSE
sign_ = TO_FM(0)
END IF
END IF
IF(x .lt. y .or. sign_ .eq. TO_FM(0)) THEN
XOUT = TO_FM(0)
ELSE
XOUT = sign_*FACTORIAL(x)/FACTORIAL(y)/FACTORIAL(x-y)
END IF
CALL FM_EXIT_USER_FUNCTION(XOUT)
END FUNCTION GENERALIZED_BINOMIAL
END MODULE ImprovedSugama