diff --git a/src/gkcoulomb_mod.f90 b/src/gkcoulomb_mod.f90
index 36597d0..57a0fcd 100644
--- a/src/gkcoulomb_mod.f90
+++ b/src/gkcoulomb_mod.f90
@@ -1,718 +1,722 @@
MODULE gkcoulomb
!
! Linearized GK Coulomb Collision operator between species a and b
!
USE prec_const
USE basic
USE basic_mpi
USE model
USE FMZM
USE coeff
USE basis_transformation
USE speed_functions
use load_prec_coeffs
IMPLICIT NONE
PUBLIC
CONTAINS
!
!
SUBROUTINE coulombSelfGK(l,k,j,n,np,rCaa)
!
! construct the like-species GYROKINETIC Coulomb Operator.
!
! Note: - use pre-computed Tljpmf coefficients.
! - routine tested and compare against coulombSelfGK_V2 using the pre-computed coefficients Tljpmf (DK Coulomb retrived in the kperp =0 limit)
!
IMPLICIT NONE
!
! input vars.
INTEGER,INTENT(in) :: l,k
INTEGER,INTENT(in) :: j,n,np
TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
! local variables
INTEGER :: p,f,m,r, icol
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=T5len) :: Tljpmfstr_,RowNpjstr_
!
CALL FM_ENTER_USER_ROUTINE
!
! PRECOMPUTE NORMALIZATION COEFFICIENTS
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
kern_ = kerneln(bafm_,n)
kernp_ =kerneln(bafm_,np)
!
ploop: DO p=0,PMMAXX
NX1 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
SPJ = sigmapj(p,j)
mloop: DO m=0, p ! p ! ... collisional harmonics
IF( m .eq. 0) THEN
bapowm = TO_FM(1)
ELSE
bapowm = 2._dp*TO_FM(bafm_/2)**(2*m)
ENDIF
floop: DO f=0,k+n
- bardmnkf_ = ALL2L(n,k,f,m)
+ bardmnkf_ = ALL2L(n,k,f,m)
+ IF( bardmnkf_ .ne. ZEROFM) THEN
IF( p .le. l+2*f+m) THEN
! get Tljpmf coeffs; pre computed
Tljpmfstr_ = GET_Tljpmf(l,j,p,m,f)
Tljpmf_ = TO_FM(Tljpmfstr_)
IF(Tljpmf_ .ne. ZEROFM) THEN
DO icol =1,numba(Pmaxa,Jmaxa)
RowNpjstr_ = curlyNpjmn_STR(p,j,m,np,icol)
RowNpj_ = TO_FM(RowNpjstr_)
rCaa(icol) = rCaa(icol)+ NX0*RowNpj_/FACTORIAL(TO_FM(np +m))/FACTORIAL(TO_FM(n +m))/SPJ*NX1*FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))*Tljpmf_*bardmnkf_*bapowm*kern_*kernp_
ENDDO
ENDIF
ENDIF ! ... if p <= l + 2*f + m
+ ENDIF
ENDDO floop
ENDDO mloop
ENDDO ploop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombSelfGK
!
!___________________________________________________________________________
!
SUBROUTINE coulombSelfGK_fulleval(l,k,j,n,np,rCaa)
!
! construct the like-species GYROKINETIC Coulomb Operator. Evaluate the full expression on the fly
!
!
IMPLICIT NONE
!
! input vars.
INTEGER,INTENT(in) :: l,k
INTEGER,INTENT(in) :: j,n,np
TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
! local variables
INTEGER :: p,f,m,r, icol, g ,h, s1
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_, T5_,X1
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=T5len) :: Tljpmfstr_,RowNpjstr_
!
CALL FM_ENTER_USER_ROUTINE
!
! PRECOMPUTE NORMALIZATION COEFFICIENTS
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
kern_ = kerneln(bafm_,n)
kernp_ =kerneln(bafm_,np)
!
ploop: DO p=0,PMMAXX
NX1 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
SPJ = sigmapj(p,j)
mloop: DO m=0, p ! ... collisional harmonics
IF( m .eq. 0) THEN
bapowm = TO_FM(1)
ELSE
bapowm = 2._dp*TO_FM(bafm_/2)**(2*m)
ENDIF
floop: DO f=0,k+n
bardmnkf_ = ALL2L(n,k,f,m)
- IF( p .le. l+2*f+m) THEN
- ! get Tljpmf coeffs; pre computed
- Tljpmf_ = Tabljpmf(l,j,p,m,f)
-!! print*, l,j,p,m,f,TO_DP(Tljpmf_)
- IF(Tljpmf_ .ne. ZEROFM) THEN
- DO g=0,min(Pmaxa,p + 2*j + m)
- X1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
- DO h =0, j + floor(real(p + m )/2)
- T5_ = TO_FM('0.0')
- T5_ = get_T5( p,j,m,g,h)
- IF(T5_ .ne. ZEROFM) THEN
- DO s1 =0, min(Jmaxa,n + m + h)
- ! Change order : n + m + h <=Jmax
- dmnhs1_ = ALLX2L(n,h,s1,m)
- icol = numba(g,s1)
- rCaa(icol) = rCaa(icol)+ NX0/FACTORIAL(TO_FM(np +m))/FACTORIAL(TO_FM(n +m))/SPJ*NX1*FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))*Tljpmf_*bardmnkf_*bapowm*kern_*kernp_* X1*T5_*dmnhs1_
-!! print*, TO_DP(rCaa(icol))
- ENDDO
-
- ENDIF
+ IF( bardmnkf_ .ne. ZEROFM) THEN
+ IF( p .le. l+2*f+m) THEN
+ ! get Tljpmf coeffs; pre computed
+ Tljpmf_ = Tabljpmf(l,j,p,m,f)
+ !! print*, l,j,p,m,f,TO_DP(Tljpmf_)
+ IF(Tljpmf_ .ne. ZEROFM) THEN
+ DO g=0,min(Pmaxa,p + 2*j + m)
+ X1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
+ DO h =0, j + floor(real(p + m )/2)
+ T5_ = TO_FM('0.0')
+ T5_ = get_T5( p,j,m,g,h)
+ IF(T5_ .ne. ZEROFM) THEN
+ DO s1 =0, min(Jmaxa,n + m + h)
+ ! Change order : n + m + h <=Jmax
+ dmnhs1_ = ALLX2L(n,h,s1,m)
+ icol = numba(g,s1)
+ rCaa(icol) = rCaa(icol)+ NX0/FACTORIAL(TO_FM(np +m))/FACTORIAL(TO_FM(n +m))/SPJ*NX1*FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))*Tljpmf_*bardmnkf_*bapowm*kern_*kernp_* X1*T5_*dmnhs1_
+ !! print*, TO_DP(rCaa(icol))
+ ENDDO
+
+ ENDIF
+ ENDDO
ENDDO
- ENDDO
- ENDIF
- ENDIF ! ... if p < l + 2*f + m
+ ENDIF
+ ENDIF ! ... if p < l + 2*f + m
+ ENDIF
ENDDO floop
ENDDO mloop
ENDDO ploop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombSelfGK_fulleval
!
!___________________________________________________________________________
!
SUBROUTINE coulombTGK(l,k,j,n,np,rCabT)
! evaluate the test part of the GK Coulomb operator
IMPLICIT NONE
!
! Input vars.
INTEGER,INTENT(in) :: l,k
INTEGER,INTENT(in) :: j,n,np
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
! local vars
INTEGER :: p,f,m,r, icol
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm, Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=T5len) :: Tljpmfstr_,RowNpjstr_
CALL FM_ENTER_USER_ROUTINE
! PRECOMPUTE NORMALIZATION COEFFICIENTS
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
kern_ = kerneln(bafm_,n)
kernp_ =kerneln(bafm_,np)
!
ploop: DO p=0,PMMAXX
NX1 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
SPJ = sigmapj(p,j)
mloop: DO m=0, p ! ... collisional harmonics
IF( m .eq. 0) THEN
bapowm = TO_FM(1)
ELSE
bapowm = 2._dp*TO_FM(bafm_/2)**(2*m)
ENDIF
floop: DO f=0,k+n
bardmnkf_ = ALL2L(n,k,f,m)
IF( p .le. l+2*f+m) THEN
! get Tljpmf coeffs; pre computed
Tljpmfstr_ = GET_Tljpmf(l,j,p,m,f)
Tljpmf_ = TO_FM(Tljpmfstr_)
IF(Tljpmf_ .ne. ZEROFM) THEN
DO icol =1,numba(Pmaxa,Jmaxa)
RowNpjstr_ = curlyNpjmn_STR(p,j,m,np,icol)
RowNpj_ = TO_FM(RowNpjstr_)
rCabT(icol) = rCabT(icol)+ NX0*RowNpj_/FACTORIAL(TO_FM(np +m))/FACTORIAL(TO_FM(n +m))/SPJ*NX1*FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))*Tljpmf_*bardmnkf_*bapowm*kern_*kernp_
ENDDO
ENDIF
ENDIF ! ... if p < l + 2*f + m
ENDDO floop
ENDDO mloop
ENDDO ploop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombTGK
!
!___________________________________________________________________________
!
SUBROUTINE coulombFGK(l,k,j,n,np,rCabF)
! evaluate the field part of the GK Coulomb operator
IMPLICIT NONE
!
! Input vars.
INTEGER,INTENT(in) :: l,k
INTEGER,INTENT(in) :: j,n,np
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabF
! local vars
INTEGER :: p,f,m,r, icol
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,bbpowm, Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=T5len) :: Tljpmfstr_,RowNpjstr_
CALL FM_ENTER_USER_ROUTINE
! PRECOMPUTE NORMALIZATION COEFFICIENTS
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
kern_ = kerneln(bafm_,n)
kernp_ =kerneln(bbfm_,np)
!
ploop: DO p=0,PMMAXX
NX1 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
SPJ = sigmapj(p,j)
mloop: DO m=0, p ! ... collisional harmonics
IF( m .eq. 0) THEN
bapowm = TO_FM(1)
bbpowm = TO_FM(1)
ELSE
bapowm = 2._dp*TO_FM(bafm_/2)**(m)
bbpowm = TO_FM(bbfm_/2)**(m)
ENDIF
floop: DO f=0,k+n
bardmnkf_ = ALL2L(n,k,f,m)
IF( p .le. l+2*f+m) THEN
! get Tljpmf coeffs; pre computed
Tljpmfstr_ = GET_Tljpmf(l,j,p,m,f)
Tljpmf_ = TO_FM(Tljpmfstr_)
IF(Tljpmf_ .ne. ZEROFM) THEN
DO icol =1,numbb(Pmaxb,Jmaxb)
RowNpjstr_ = curlyNpjmn_STR(p,j,m,np,icol)
RowNpj_ = TO_FM(RowNpjstr_)
rCabF(icol) = rCabF(icol)+ NX0*RowNpj_/FACTORIAL(TO_FM(np +m))/FACTORIAL(TO_FM(n +m))/SPJ*NX1*FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))*Tljpmf_*bardmnkf_*bapowm*bbpowm*kern_*kernp_
ENDDO
ENDIF
ENDIF ! ... if p < l + 2*f + m
ENDDO floop
ENDDO mloop
ENDDO ploop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombFGK
!
!___________________________________________________________________________
!
! routines needed to compute gk coulomb on the fly
!
FUNCTION Tabljpmf(l,j,p,m,f) result(XOUT)
! compute Tabljpmf coefficient
implicit none
!
INTEGER,intent(in) :: l,j,p,m,f
TYPE(FM):: XOUT
TYPE(FM),SAVE :: T5inv,Lptd,nuF_,nuT_
INTEGER :: t,d
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM(0)
!
do t= 0,f + floor(real(l+m)/2)
T5inv = GET_T5inv(l,f,p,t,m)
IF(T5inv .ne. ZEROFM) THEN
do d =0,t
Lptd = Lpjl(real(p,dp),real(t,dp),real(d,dp))
!! if( coll_type .eq. 'self') THEN
nuT_ = nusTabpjgd(p,j,p,d)
nuF_ = nusFabpjgd(p,j,p,d)
XOUT = XOUT + Lptd*( nuT_ + nuF_)*T5inv
!! ELSEIF( coll_type .eq. 'test' ) THEN
! nuT_ = nusTabpjgd(p,j,p,d)
! XOUT = XOUT + Lptd*nuT_*T5inv
! ! ELSEIF ( coll_type .eq. 'field' ) THEN
! nuF_ = nusFabpjgd(p,j,p,d)
! XOUT = XOUT + Lptd*nuF_*T5inv
! ENDIF
enddo
ENDIF
enddo
!
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION Tabljpmf
!
!___________________________________________________________________________
!
! TEST ROUTINES - NOT USED ANYMORE
! SUBROUTINE coulombTGK(l,k,j,rCabT)
! !
! ! construct the test-part of the GYROKINETIC Coulomb Operator with the full collisional harmonics by energy components
! !
! !
! IMPLICIT NONE
! !
! INTEGER,INTENT(in) :: l,k,j
! TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
! ! Declare local variables
! INTEGER :: p,n,g,h,s1,np,t,f,d,m,mmaxx
! TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T5,T5inv,SPJ,kern_,kernp_
! TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,nuT_,bapowm,Lptd_,Tabmpfjl_
! INTEGER :: icol
! CHARACTER(len=T5len) :: T5str_
! !
! CALL FM_ENTER_USER_ROUTINE
! !
! ! PRECOMPUTE NORMALIZATION COEFFICIENTS
! NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
! !
! ! mmaxx = 0
! !
! nloop: DO n=0,naFLR
! kern_ = kerneln(bafm_,n)
! floop: DO f=0,k+n
! ploop: DO p=0,l+2*f
! NX2 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
! SPJ = sigmapj(p,j)
! IF (p .eq. 0) mmaxx =0
! mloop: DO m=0, p ! ... collisional harmonics
! IF( m .eq. 0) THEN
! bapowm = TO_FM(1)
! ELSE
! bapowm = 2*(bafm_/2)**(2*m)
! ENDIF
! bardmnkf_ = ALL2L(n,k,f,m)
! Tabmpfjl_ = Tabmpfjl(m,p,f,j,l)
! IF(Tabmpfjl_ .ge. TO_FM('1E-40') .or. bardmnkf_ .ne. TO_FM('1E-40')) THEN
! gloop: DO g=0, min(Pmaxa,p+ 2*j)
! NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
! hloop: DO h=0, j + floor(real(p)/2)
! T5STR_ = GET_T5(p,j,m,g,h)
! T5 = TO_FM(T5STR_)
! IF(T5 .ge. TO_FM('1E-40')) THEN
! nploop: DO np=0,naFLR
! kernp_ = kerneln(bafm_,np)
! NX3 = FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))/FACTORIAL(TO_FM(n + m ))/FACTORIAL(TO_FM(np + m))
! s1loop: DO s1 =0, min(Jmaxa,np + m + h)
! dmnhs1_ = ALLX2L(np,h,s1,m)
! icol = numba(g,s1)
! rCabT(icol) = rCabT(icol ) + NX0*NX1*NX2*NX3*kern_*kernp_*dmnhs1_*bardmnkf_*T5*Tabmpfjl_/SPJ*bapowm/FACTORIAL(TO_FM(n + m ))/FACTORIAL(TO_FM(np + m))
! ENDDO s1loop
! ENDDO nploop
! ENDIF ! ... end if T5 .ne. 0
! ENDDO hloop
! ENDDO gloop
! ENDIF
! ENDDO mloop
! ENDDO ploop
! ENDDO floop
! ENDDO nloop
! !
! !
! CALL FM_EXIT_USER_ROUTINE
! !
! END SUBROUTINE coulombTGK
!
!___________________________________________________________________________
!
! SUBROUTINE coulombFGK(l,k,rCabF)
! !
! ! construct the field-part of the GYROKINETIC Coulomb Operator with the full collisional harmonics
! !
! !
! IMPLICIT NONE
! !
! INTEGER,INTENT(in) :: l,k
! TYPE(FM),DIMENSION(:),INTENT(inout):: rCabF
! ! Declare local variables
! INTEGER :: p,j,n,g,h,s1,np,t,f,d,m
! TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T5,T5inv,SPJ,kern_,kernp_
! TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,nuF_,bapowm,Lptd_
! INTEGER :: icol
! CHARACTER(len=T5len) :: T5str_
! !
! CALL FM_ENTER_USER_ROUTINE
! !
! ! PRECOMPUTE NORMALIZATION COEFFICIENTS
! NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
! !
! !
! jloop: DO j=0,JEMaxx
! nloop: DO n=0,nbFLR
! kern_ = kerneln(bafm_,n)
! floop: DO f=0,k+n
! ploop: DO p=0,l+2*f
! NX2 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
! SPJ = sigmapj(p,j)
! mloop: DO m=0, p ! ... collisional harmonics
! IF( m .eq. 0) THEN
! bapowm = TO_FM(1)
! ELSE
! bapowm = 2*(bafm_/2)**(2*m)
! ENDIF
! bardmnkf_ = ALL2L(n,k,f,m)
! tloop: DO t=0, f + floor(real(l)/2)
! T5STR_ = GET_T5(p,t,m,l,f)
! T5inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t)) &
! *TO_FM(2*p + 1)/2/FACTORIAL(TO_FM( 2*p + 2*t +1 )/2)/FACTORIAL(TO_FM(p + m))*FACTORIAL(TO_FM(p-m))*TO_FM(T5STR_)
! IF(T5inv .ne. TO_FM('0')) THEN
! gloop: DO g=0, min(Pmaxb,p+ 2*j)
! NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
! hloop: DO h=0, j + floor(real(p)/2)
! T5STR_ = GET_T5(p,j,m,g,h)
! T5 = TO_FM(T5STR_)
! IF(T5 .ne. TO_FM('0')) THEN
! dloop: DO d=0,t
! Lptd_ = Lpjl(real(p,dp),real(t,dp),real(d,dp))
! nuF_ = nusFabpjgd(p,j,p,d)
! nploop: DO np=0,nbFLR
! kernp_ = kerneln(bafm_,np)
! NX3 = FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))/FACTORIAL(TO_FM(n + m ))/FACTORIAL(TO_FM(np + m))
! s1loop: DO s1 =0, min(Jmaxb,np + m + h)
! dmnhs1_ = ALLX2L(np,h,s1,m)
! icol = numbb(g,s1)
! rCabF(icol) = rCabF(icol ) + NX0*NX1*NX2*NX3*kern_*kernp_*dmnhs1_*bardmnkf_*T5*T5inv*Lptd_*nuF_/SPJ*bapowm
! ENDDO s1loop
! ENDDO nploop
! ENDDO dloop
! ENDIF ! ... end if T5 .ne. 0
! ENDDO hloop
! ENDDO gloop
! ENDIF ! ... end if T5inv .ne.0
! ENDDO tloop
! ENDDO mloop
! ENDDO ploop
! ENDDO floop
! ENDDO nloop
! ENDDO jloop
! !
! !
! CALL FM_EXIT_USER_ROUTINE
! !
! END SUBROUTINE coulombFGK
!
!___________________________________________________________________
!
! SUBROUTINE coulombSelfGK(l,k,rCaa)
! !
! ! construct the like-species GYROKINETIC Coulomb Operator with the full collisional harmonics
! !
! !
! IMPLICIT NONE
! !
! INTEGER,INTENT(in) :: l,k
! TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
! ! Declare local variables
! INTEGER :: p,j,n,g,h,s1,np,t,f,d,m
! TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T5,T5inv,SPJ,kern_,kernp_
! TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,nuT_,bapowm,Lptd_,nuF_
! INTEGER :: icol
! CHARACTER(len=T5len) :: T5str_
! !
! CALL FM_ENTER_USER_ROUTINE
! !
! ! PRECOMPUTE NORMALIZATION COEFFICIENTS
! NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
! !
! !
! jloop: DO j=0,JEmaxx
! nloop: DO n=0,naFLR
! kern_ = kerneln(bafm_,n)
! floop: DO f=0,k+n
! ploop: DO p=0,l+2*f
! NX2 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
! SPJ = sigmapj(p,j)
! mloop: DO m=0, p ! ... collisional harmonics
! IF( m .eq. 0) THEN
! bapowm = TO_FM(1)
! ELSE
! bapowm = 2*(bafm_/2)**(2*m)
! ENDIF
! bardmnkf_ = ALL2L(n,k,f,m)
! tloop: DO t=0, f + floor(real(l)/2)
! T5STR_ = GET_T5(p,t,m,l,f)
! T5inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t)) &
! *TO_FM(2*p + 1)/2/FACTORIAL(TO_FM( 2*p + 2*t +1 )/2)/FACTORIAL(TO_FM(p + m))*FACTORIAL(TO_FM(p-m))*TO_FM(T5STR_)
! IF(T5inv .ne. TO_FM('0')) THEN
! gloop: DO g=0, min(Pmaxa,p+ 2*j)
! NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
! hloop: DO h=0, j + floor(real(p)/2)
! T5STR_ = GET_T5(p,j,m,g,h)
! T5 = TO_FM(T5STR_)
! IF(T5 .ne. TO_FM('0')) THEN
! dloop: DO d=0,t
! Lptd_ = Lpjl(real(p,dp),real(t,dp),real(d,dp))
! nuT_ = nusTabpjgd(p,j,p,d)
! nuF_ = nusFabpjgd(p,j,p,d)
! nploop: DO np=0,naFLR
! kernp_ = kerneln(bafm_,np)
! NX3 = FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))/FACTORIAL(TO_FM(n + m ))/FACTORIAL(TO_FM(np + m))
! s1loop: DO s1 =0, min(Jmaxa,np + m + h)
! dmnhs1_ = ALLX2L(np,h,s1,m)
! icol = numba(g,s1)
! rCaa(icol) = rCaa(icol ) + NX0*NX1*NX2*NX3*kern_*kernp_*dmnhs1_*bardmnkf_*T5*T5inv*Lptd_*( nuT_ + nuF_)/SPJ*bapowm
! ENDDO s1loop
! ENDDO nploop
! ENDDO dloop
! ENDIF ! ... end if T5 .ne. 0
! ENDDO hloop
! ENDDO gloop
! ENDIF ! ... end if T5inv .ne.0
! ENDDO tloop
! ENDDO mloop
! ENDDO ploop
! ENDDO floop
! ENDDO nloop
! ENDDO jloop
! !
! !
! CALL FM_EXIT_USER_ROUTINE
! !
! END SUBROUTINE coulombSelfGK
!
!___________________________________________________________________
!
! SUBROUTINE coulombSelfGK_V2(l,k,j,n,np,rCaa)
! !
! ! construct the like-species GYROKINETIC Coulomb Operator with the full collisional harmonics for a given energy, j, and FLR indices, n and np.
! !
! ! Note: aiming to be
! IMPLICIT NONE
! !
! INTEGER,INTENT(in) :: l,k
! ! Energy and FLR indices
! INTEGER,INTENT(in) :: j,n,np
! TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
! ! Declare local variables
! INTEGER :: p,g,h,s1,t,f,d,m,r
! TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T5,T5inv,SPJ,kern_,kernp_,XX
! TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,nuT_,bapowm,Lptd_,nuF_
! INTEGER :: icol
! CHARACTER(len=T5len) :: T5str_
! !
! CALL FM_ENTER_USER_ROUTINE
! !
! ! PRECOMPUTE NORMALIZATION COEFFICIENTS
! NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
! !
! !
! kern_ = kerneln(bafm_,n)
! floop: DO f=0,k+n
! ploop: DO p=0,PMMAXX
! mloop: DO m=0, p ! ... collisional harmonics
! IF( m .eq. 0) THEN
! bapowm = TO_FM(1)
! ELSE
! bapowm = 2*(bafm_/2)**(2*m)
! ENDIF
! bardmnkf_ = ALL2L(n,k,f,m)
! rloop: DO r=0,l+m +2*f
! IF( r .eq. p) THEN
! !
! !! ========== Compute Tljpmf ========
! XX = TO_FM(0)
! tloop: DO t=0, f + floor(real(l + m )/2)
! T5inv = GET_T5inv(l,f,p,t,m)
! IF(T5inv .ne. TO_FM('0')) THEN
! dloop: DO d=0,t
! Lptd_ = Lpjl(real(p,dp),real(t,dp),real(d,dp))
! nuT_ = nusTabpjgd(p,j,p,d)
! nuF_ = nusFabpjgd(p,j,p,d)
! XX = XX + Lptd_*( nuT_ + nuF_ )*T5inv
! ENDDO dloop
! ENDIF
! ENDDO tloop
! !==================================
! gloop: DO g=0, min(Pmaxa,p+ m + 2*j)
! NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
! hloop: DO h=0, j + floor(real(p+m)/2)
! T5STR_ = GET_T5(p,j,m,g,h)
! T5 = TO_FM(T5STR_)
! IF(T5 .ne. TO_FM('0')) THEN
! kernp_ = kerneln(bafm_,np)
! s1loop: DO s1 =0, min(Jmaxa,np + m + h)
! dmnhs1_ = ALLX2L(np,h,s1,m)
! icol = numba(g,s1)
! rCaa(icol) = rCaa(icol ) + NX1*kern_*kernp_*dmnhs1_*T5*XX*bapowm
! ENDDO s1loop
! ENDIF ! ... end if T5 .ne. 0
! ENDDO hloop
! ENDDO gloop
! ENDIF
! ENDDO rloop
! ENDDO mloop
! ENDDO ploop
! ENDDO floop
! !
! !
! CALL FM_EXIT_USER_ROUTINE
! !
! END SUBROUTINE coulombSelfGK_V2
!
!___________________________________________________________________
!
! SUBROUTINE coulombSelfGK_V3(l,k,j,n,np,rCaa)
! !
! ! construct the like-species GYROKINETIC Coulomb Operator with the full collisional harmonics for a given energy, j, and FLR indices, n and np.
! !
! ! Note: - use pre-computed Tljpmf coefficients.
! ! - routine tested and compare against coulombSelfGK_V2 using the pre-computed coefficients Tljpmf (DK Coulomb retrived in the kperp =0 limit)
! !
! IMPLICIT NONE
! !
! INTEGER,INTENT(in) :: l,k
! ! Energy and FLR indices
! INTEGER,INTENT(in) :: j,n,np
! TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
! ! Declare local variables
! INTEGER :: p,g,h,s1,f,m
! TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T5,T5inv,SPJ,kern_,kernp_
! TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,Lptd_,Tljpmf_
! INTEGER :: icol
! CHARACTER(len=T5len) :: T5str_,Tljpmfstr_
! !
! CALL FM_ENTER_USER_ROUTINE
! !
! ! PRECOMPUTE NORMALIZATION COEFFICIENTS
! NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
! !
! !
! kern_ = kerneln(bafm_,n)
! floop: DO f=0,k+n
! ploop: DO p=0,l+2*f
! NX2 = TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p +1))
! SPJ = sigmapj(p,j)
! mloop: DO m=0, p ! ... collisional harmonics
! IF( m .eq. 0) THEN
! bapowm = TO_FM(1)
! ELSE
! bapowm = 2*(bafm_/2)**(2*m)
! ENDIF
! bardmnkf_ = ALL2L(n,k,f,m)
! ! get Tljpmf coeffs; pre computed
! Tljpmfstr_ = GET_Tljpmf(l,j,p,m,f)
! Tljpmf_ = TO_FM(Tljpmfstr_)
! IF(Tljpmf_ .ne. TO_FM('0')) THEN
! gloop: DO g=0, min(Pmaxa,p+ 2*j)
! NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
! hloop: DO h=0, j + floor(real(p)/2)
! T5STR_ = GET_T5(p,j,m,g,h)
! T5 = TO_FM(T5STR_)
! IF(T5 .ne. TO_FM('0')) THEN
! kernp_ = kerneln(bafm_,np)
! NX3 = FACTORIAL(TO_FM(n))*FACTORIAL(TO_FM(np))/FACTORIAL(TO_FM(n + m ))/FACTORIAL(TO_FM(np + m))
! s1loop: DO s1 =0, min(Jmaxa,np + m + h)
! dmnhs1_ = ALLX2L(np,h,s1,m)
! icol = numba(g,s1)
! rCaa(icol) = rCaa(icol ) + NX0*NX1*NX2*NX3*kern_*kernp_*dmnhs1_*bardmnkf_*T5*Tljpmf_/SPJ*bapowm
! ENDDO s1loop
! ENDIF
! ENDDO hloop
! ENDDO gloop
! ENDIF
! ENDDO mloop
! ENDDO ploop
! ENDDO floop
! !
! !
! CALL FM_EXIT_USER_ROUTINE
! !
! END SUBROUTINE coulombSelfGK_V3
!
!___________________________________________________________________
!
END MODULE gkcoulomb