Homec4science
Diffusion COSOlver (improved_sugama) Browse

cosolver/src/coll_sugama.h
cd1286e2d410improved_sugama

Last Change
Raw File

coll_sugama.h
View Options

! This file contains routines to compute the Hermite-Laguerre drift and gyrokinetic
! Sugama collision operator
!
! Ref: Sugama et. al., Phsyics of Plasmas 16, 11 (2009)
!
!_________________________________________________________________
!
!
! DRIFT-KINETIC FIELD/BACK-REACTION SUGAMA
!
!
!_________________________________________________________________
SUBROUTINE SugamaTDK(l,k,rCabT)
!
! Construct the drift-kinetic field Sugama collision operator
!
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabT
! CabS[f]^lk
CALL SugamaTCabDK(l,k,rCabT)
! CPa[f]^lk
CALL SugamaTCPabDK(l,k,rCabT)
! PaC[f]^lk
CALL SugamaTPCabDK(l,k,rCabT)
! PaCPa[f]^lk
CALL SugamaTPCPabDK(l,k,rCabT)
!
END SUBROUTINE SugamaTDK
!__________________________________________________
SUBROUTINE SugamaTCabDK(l,k,rCabT)
!
! Construct CabS[s]^lk of the test part of the Sugama collision operator
!
! Ref: Sugama et al., Physics of Plasmas 16, 112503 (2009)
! COMPLETED
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
! Declare local variables
INTEGER :: p,h,j,d,t,g
TYPE(FM),SAVE :: NX0,NX1,NX2,SPJ,Lphd,T4,T4inv,nusS
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_, T4invstr_
!
CALL FM_ENTER_USER_ROUTINE
!
irow = numba(l,k)
!
NX0 = 1/SQRT((TO_FM(2)**l)*FACTORIAL(TO_FM(l)))
!
jloop: DO j=0,jamaxx
ploop: DO p=0,l + 2*k
SPJ = sigmapj(p,j)
NX1 = (TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM(2*p +1)
tloop: DO t=0, k + floor(real(l)/2)
T4invstr_ = GET_T4(p,t,l,k) ! ... (T^-1)^pt_lk = ... T^lk_pt
T4inv =TO_FM(T4invstr_)*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)
IF( T4inv .ne. TO_FM('0.0')) THEN
gloop: DO g=0, min(p + 2*j,Pmaxa)
NX2 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0, min(j + floor(real(p)/2),Jmaxa)
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
icol = numba(g,h)
IF( T4 .ne. TO_FM('0.0')) THEN
dloop: DO d=0, t
Lphd=Lpjl(real(p,dp),real(t,dp),real(d,dp))
nusS = nusSabpjd(p,j,d)
rCabT(icol) = rCabT(icol) + NX0*NX1*T4inv*NX2*T4*Lphd*nusS/SPJ
ENDDO dloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
CALL FM_EXIT_USER_ROUTINE
!
!
END SUBROUTINE SugamaTCabDK
!_________________________________________________________________________________
!
SUBROUTINE SugamaTCPabDK(l,k,rCabT)
!
! Construct the CabSPa[f]^lk
!
!
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: l,k
TYPE(FM), DIMENSION(:), INTENT(INOUT) :: rCabT
!
! LOCAL VARS.
!
INTEGER :: icol
TYPE(FM) :: NX0,NX1,I1lkpara,Ilk,coeff
!
CALL FM_ENTER_USER_ROUTINE
!
NX0 = 1/ SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
NX1 = - (1 + chifm_**2)
!
coeff = thetaabfm_ -TO_FM(1)
!
IF(coeff .ne. TO_FM(0)) THEN
!
I1lkpara = IabSparalk(l,k)
Ilk = IabSlk(l,k)
!
! Perturbed parallel velocity u_a[f]^lk\cdot e^{10}
!
rCabT(numba(1,0)) = rCabT(numba(1,0)) + coeff*SQRT(TO_FM(2))*NX0*NX1*I1lkpara
!
!
! Perturbed temperature deltaT_a[f]^lk
!
rCabT(numba(2,0)) = rCabT(numba(2,0)) + coeff*Ilk*NX0*SQRT(TO_FM(2))/TO_FM(3)
rCabT(numba(0,1)) = rCabT(numba(0,1)) - coeff*Ilk*2*NX0/TO_FM(3)
!
!
ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
!
END SUBROUTINE SugamaTCPabDK
!_________________________________________________________________________________
!
SUBROUTINE SugamaTPCPabDK(l,k,rCabT)
!
! Construct the PaCabSPa[f]^lk
!
! COMPLETED
IMPLICIT NONE
!
INTEGER, INTENT(in) :: l,k
TYPE(FM), DIMENSION(:), INTENT(INOUT) :: rCabT
!
! LOCAL VARS.
!
INTEGER :: irow,icol
TYPE(FM) :: NX0,NX1,alphapara,alphaE,coeff
!
irow = numba(l,k)
!
CALL FM_ENTER_USER_ROUTINE
!
coeff = (thetaabfm_ -TO_FM(1))**2
!
IF (coeff .ne. TO_FM(0)) THEN
!
alphapara = - SQRT(TO_FM(2))*chifm_/SQRT((TO_FM(1) + chifm_**2))*4/3/SQRT(PIFM)
alphaE = - 2*chifm_/( TO_FM(1) + chifm_**2)**(TO_FM(3)/2)*4/3/SQRT(PIFM)
!
IF(l .eq. 1 .and. k .eq. 0) THEN
!
! Perturbed parallel velocity u_a[f]^lk\cdot e^{10}
!
rCabT(numba(1,0)) = rCabT(numba(1,0)) + coeff*alphapara/SQRT(TO_FM(2))
!
ENDIF
!
IF( l .eq. 2 .and. k .eq. 0) THEN
!
! Perturbed temperature deltaT_a[f]^lk
!
rCabT(numba(2,0)) = rCabT(numba(2,0)) + coeff*alphaE/TO_FM(3)
rCabT(numba(0,1)) = rCabT(numba(0,1)) - coeff*2*alphaE/TO_FM(3)/SQRT(TO_FM(2))
!
ENDIF
!
IF( l .eq. 0 .and. k .eq. 1) THEN
!
! Perturbed temperature deltaT_a[f]^lk
!
rCabT(numba(2,0)) = rCabT(numba(2,0)) - coeff*SQRT(TO_FM(2))*alphaE/TO_FM(3)
rCabT(numba(0,1)) = rCabT(numba(0,1)) + coeff*2*alphaE/TO_FM(3)
!
ENDIF
!
ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
!
END SUBROUTINE SugamaTPCPabDK
!_________________________________________________________________________________
!
SUBROUTINE SugamaTPCabDK_V0(l,k,rCabT)
!
! Construct PC_{ab}[f]^lk of the test part of the Sugama collision operator
!
! Note: Version 0.0
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:), INTENT(inout) :: rCabT
! Declare local variables
INTEGER :: p,h,j,g
TYPE(FM),SAVE :: NX0,NX1,NX2,SPJ,L100,T4,nuS1j0,nuS0j0,nuS0j1,L010_,L011_,nus_,coeff
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_, T4invstr_
!
CALL FM_ENTER_USER_ROUTINE
!
irow = numba(l,k)
!
coeff = (thetaabfm_ -TO_FM(1))
!
IF(coeff .ne. TO_FM(0)) THEN
!
IF(l .eq. 1 .and. k .eq. 0) THEN
!
jloop: DO j=0,jamaxx
! ___________ SPEED FUNCTIONS _____________
nuS1j0 = nusSabpjd(1,j,0)
L100 = Lpjl(real(1,dp),real(0,dp),real(0,dp))
!__________________________________________
NX0 = 2*FACTORIAL(TO_FM(1)/2)*FACTORIAL(TO_FM(j))/FACTORIAL(TO_FM(2*j + 3)/2)*SQRT(TO_FM(2)/SQRT(PIFM))
gloop: DO g=0, min(1 + 2*j,Pmaxa)
NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
hloop: DO h=0, min(j,Jmaxa)
T4str_ = GET_T4(1,j,g,h) !T^gh_1j
T4 = TO_FM(T4str_)
icol = numba(g,h)
rCabT(icol) = rCabT(icol) + coeff*T4*NX0*NX1/2*L100*nuS1j0
ENDDO hloop
ENDDO gloop
ENDDO jloop
!
ENDIF
! !
IF(l .eq. 2 .and. k .eq. 0) THEN
DO j=0,jamaxx
! ___________ SPEED FUNCTIONS _____________
nuS0j0 = nusSabpjd(0,j,0)
L010_ = Lpjl(real(0,dp),real(1,dp),real(0,dp))
L011_ = Lpjl(real(0,dp),real(1,dp),real(1,dp))
nuS0j1 = nusSabpjd(0,j,1)
!________________________
NX0 = - FACTORIAL(TO_FM(1)/2)*FACTORIAL(TO_FM(j))/FACTORIAL(TO_FM(2*j + 1)/2)*SQRT(TO_FM(2))/3
DO g=0, min( 2*j,Pmaxa)
NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
DO h=0, min(j,Jmaxa)
T4str_ = GET_T4(0,j,g,h) !T^gh_0j
T4 = TO_FM(T4str_)
icol = numba(g,h)
rCabT(icol) = rCabT(icol) + coeff*T4*NX1*NX0*(L010_*nuS0j0 + L011_*nuS0j1)
ENDDO
ENDDO
ENDDO
ENDIF
! !
! !
IF(l .eq. 0 .and. k .eq. 0) THEN
!
DO j=0,jamaxx
! ___________ SPEED FUNCTIONS _____________
nuS0j0 = nusSabpjd(0,j,0)
nuS0j1 = nusSabpjd(0,j,1)
L010_ = Lpjl(real(0,dp),real(1,dp),real(0,dp))
L011_ = Lpjl(real(0,dp),real(1,dp),real(1,dp))
!________________________
NX0 = FACTORIAL(TO_FM(1)/2)*FACTORIAL(TO_FM(j))/FACTORIAL(TO_FM(2*j + 1)/2)*TO_FM(2)/3
DO g=0, min( 2*j,Pmaxa)
NX1 = SQRT(TO_FM(2)**g*FACTORIAL(TO_FM(g)))
DO h=0, min(j,Jmaxa)
T4str_ = GET_T4(0,j,g,h) !T^gh_0j
T4 = TO_FM(T4str_)
icol = numba(g,h)
rCabT(icol) = rCabT(icol) + coeff*T4*NX1*NX0*(L010_*nuS0j0 + L011_*nuS0j1)
ENDDO
ENDDO
ENDDO
ENDIF
!
!
ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaTPCabDK_V0
!
!_________________________________________________________________
!
SUBROUTINE SugamaTPCabDK(l,k,rCabT)
!
! Construct PC_{ab}[f]^lk of the test part of the Sugama collision operator
!
! Note: New verision 1.0
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:), INTENT(inout) :: rCabT
!
! Local var
INTEGER :: icol,p,j
TYPE(FM),SAVE :: I1abSparapj_,NX0,IabESpj_,coeff
!
CALL FM_ENTER_USER_ROUTINE
!
coeff = (thetaabfm_ -TO_FM(1))
!
IF(coeff .ne. TO_FM(0)) THEN
!
IF (l .eq. 1 .and. k .eq. 0 ) THEN
DO p=0,Pmaxa
NX0 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
DO j=0,Jmaxa
I1abSparapj_ = IabSparalk(p,j)
icol = numba(p,j)
rCabT(icol) = rCabT(icol) - coeff*(1 + chifm_**2)*SQRT(TO_FM(2))*NX0*I1abSparapj_
ENDDO
ENDDO
ENDIF
!
IF( l .eq. 2 .and. k .eq. 0) THEN
!
DO p=0,Pmaxa
NX0 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
DO j=0,Jmaxa
IabESpj_ = IabSlk(p,j)
icol = numba(p,j)
rCabT(icol) = rCabT(icol) + coeff*NX0*SQRT(TO_FM(2))/3*IabESpj_
ENDDO
ENDDO
!
ENDIF
!
IF(l .eq. 0 .and. k .eq. 1 ) THEN
DO p=0,Pmaxa
NX0 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
DO j=0,Jmaxa
IabESpj_ = IabSlk(p,j)
icol = numba(p,j)
rCabT(icol) = rCabT(icol) - coeff*NX0*IabESpj_*2/3
ENDDO
ENDDO
!
ENDIF
!
!
ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaTPCabDK
!_________________________________________________________________
!
!
! DRIFT-KINETIC IMPROVED SUGAMA
!
!
!_________________________________________________________________
SUBROUTINE SugamaImprovedTDK(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)
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_
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedTDK
!_________________________________________________________________
!
SUBROUTINE SugamaImprovedFDK(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)
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_
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedFDK
!_________________________________________________________________
!
!
SUBROUTINE SugamaImprovedSelfDK(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_{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
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_)
ENDDO hloop
ENDDO gloop
ENDDO iloop
ENDIF
ENDDO jloop
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedSelfDK
!_________________________________________________________________
!
!
! DRIFT-KINETIC FIELD/BACK-REACTION SUGAMA
!
!
!_________________________________________________________________
SUBROUTINE SugamaFDK_V1(l,k,rCabF)
!
! Construct the drift-kinetic back-reaction Sugama collision operator
!
!
! Note: - not used anymore
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabF
!
! LOCAL VAR.
INTEGER :: irow,icol
TYPE(FM), SAVE :: gammaab_,etaab_,T41010,X0,curlyVab_,curlyWab_,T40110,T40101,T40111,T40120,T40121
TYPE(FM), DIMENSION(:),ALLOCATABLE,SAVE :: rCbaT_
TYPE(FM),SAVE :: buffchi,buffthetaab
!
!
CALL FM_ENTER_USER_ROUTINE
!
call allocate_array(rCbaT_,numbb(0,0),numbb(Pmaxb,Jmaxb))
!
buffchi = chifm_
buffthetaab = thetaabfm_
gammaab_ = - chifm_*(taufm_ + chifm_**2)/((1 + chifm_**2)**(TO_FM(3)/2))*sigmafm_
etaab_ = - 3*chifm_*(taufm_ + chifm_**2)/((1 + chifm_**2)**(TO_FM(5)/2))*taufm_
!
!
T41010 = TO_FM(GET_T4(1,0,1,0)) !
T40110 = TO_FM(GET_T4(0,1,1,0))
T40101 = TO_FM(GET_T4(0,1,0,1))
T40111 = TO_FM(GET_T4(0,1,1,1))
T40120 = TO_FM(GET_T4(0,1,2,0))
T40121 = TO_FM(GET_T4(0,1,2,1))
X0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
curlyVab_ = curlyVab(l,k)
curlyWab_ = curlyWab(l,k)
!
!
! Invert model to compute CTba[fb] collisional moments
! Invert: a -> b, b-> a
call invert_ab
!
rCbaT_(:) = TO_FM(0)
call SugamaTDK(1,0,rCbaT_)
rCabF(:) = rCabF(:) + rCbaT_(:)*T41010*2*SQRT(TO_FM(2))*X0*curlyVab_/gammaab_/buffchi*(1 + buffchi**2)*buffthetaab &
+ rCbaT_(:)*SQRT(TO_FM(2))*curlyWab_*X0*T40110/etaab_
!
rCbaT_(:) = TO_FM(0)
call SugamaTDK(2,0,rCbaT_)
rCabF(:) = rCabF(:) + rCbaT_(:)*2*SQRT(TO_FM(2))*curlyWab_*X0*T40120/etaab_
!
rCbaT_(:) = TO_FM(0)
call SugamaTDK(0,1,rCbaT_)
rCabF(:) = rCabF(:) + rCbaT_(:)*curlyWab_*X0*T40101/etaab_
!
rCbaT_(:) = TO_FM(0)
call SugamaTDK(1,1,rCbaT_)
rCabF(:) = rCabF(:) + rCbaT_(:)*curlyWab_*X0*SQRT(TO_FM(2))*T40111/etaab_
!
rCbaT_(:) = TO_FM(0)
call SugamaTDK(2,1,rCbaT_)
rCabF(:) = rCabF(:) + rCbaT_(:)*curlyWab_*X0*2*SQRT(TO_FM(2))*T40121/etaab_
! Invert: b -> a, a-> b
call invert_ab
!
! Free memory
!
deallocate(rCbaT_)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaFDK_V1
!
!___________________________________________________________________
!
SUBROUTINE SugamaFDK(l,k,rCabF)
!
! Construct the drift-kinetic back-reaction Sugama collision operator
!
!
! Note: - Computes the field part using explicity the speed integrals E_{ba}^k and e_{ba}^l in
! \mathcal{V}_{ba}^{pj } and \mathcal{W}_{ba}^{pj }
! - there is an additional factor of 2 added to match the full DK Coulomb, need to be investigated further
!
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabF
!
! local vars.
INTEGER ::p,j,icol
TYPE(FM), SAVE :: curlyVablk_,curlyVbapj_,curlyWablk_,curlyWbapj_
TYPE(FM), SAVE :: NX0,NX1,NX2,NX3,NX4
CALL FM_ENTER_USER_ROUTINE
!
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
NX2 = (sigmafm_**2)*(chifm_**3)*( (TO_FM(1) + chifm_**2)**(1.5_dp)/sigmafm_/chifm_/(taufm_ + chifm_**2))*4/3/SQRT(PIFM) ! ... m_b/(tau^S_{ab} \tau_{ba}^S)/gamma_{ab} = ...
!
NX3 = 4*thetaabfm_*thetabafm_ ! ... need to check this additionnal factor 2 missing to match coulomb in the limit ma/mb<<1 !!!
!
NX4= (1 + chifm_**2)**(2.5_dp)/3/chifm_/(taufm_ + chifm_**2)/taufm_*3*SQRT(PIFM)/4 ! ... T_b/eta_{ab} = ...
!
!
curlyVablk_ = curlyVab(l,k)
curlyWablk_ = curlyWab(l,k)
!
ploop: DO p=0,Pmaxb
NX1 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
jloop: DO j=0,Jmaxb
icol = numbb(p,j)
curlyVbapj_ = curlyVba(p,j)
curlyWbapj_ = curlyWba(p,j)
rCabF(icol) = rCabF(icol) + NX0*NX1*NX3*NX2*( TO_FM(1) + TO_FM(1)/(chifm_**2))*( TO_FM(1) + chifm_**2)*curlyVablk_*curlyVbapj_/chifm_ &
+ 2*NX4*NX0*NX1*curlyWbapj_*curlyWablk_ ! ... factor 2 important for energy conservation
ENDDO jloop
ENDDO ploop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaFDK
!_________________________________________________________________
!
!
! DRIFT KINETIC SELF SUGAMA
!
!
!_________________________________________________________________
SUBROUTINE SugamaSelfDK(l,k,rCaa)
!
! computes the Drift-kinetic self-collision Sugama
!
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCaa
CALL FM_ENTER_USER_ROUTINE
!
!
call SugamaTDK(l,k,rCaa)
call SugamaFDK(l,k,rCaa)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaSelfDK
!_________________________________________________________________
!
!
! GYROKINETIC SUGAMA TEST PART
!
!
!_________________________________________________________________
SUBROUTINE SugamaTGK(l,k,rCabT)
!
! Computes the test-part of the gyrokinetic Sugama collision operator
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
CALL FM_ENTER_USER_ROUTINE
!
! Drift kinetic part
CALL SugamaTCabDK(l,k,rCabT)
!
! FLR corrections
CALL Rablk(l,k,rCabT) ! ... ok !
! Add GK (theta_ab - 1)*(Cab[Pa[f]] + Pa[Cab[f]])
!
IF(thetaabfm_ .ne. TO_FM(1)) THEN ! ... if theta_{ab} =1, do nothing ...
!
CALL UparalkcurlyRabpara1(l,k,rCabT) ! ... ok !
CALL UperplkcurlyRabpara2(l,k,rCabT) ! ... ok !
CALL TperplkcurlyRabpara3(l,k,rCabT) ! ... ok !
CALL Fab0paralkupara(l,k,rCabT) ! ... ok !
CALL Fab1paralkuperp(l,k,rCabT) ! ... ok !
CALL FabElkTa(l,k,rCabT) ! ... ok !
! add Gk (theta - 1)**2 PaCaPa: Matches the drift-kinetic part
CALL Uparalkupara(l,k,rCabT) ! .... ok !
CALL Uperplkuperp(l,k,rCabT) ! .... ok !
CALL curlyTperpTa(l,k,rCabT) ! .... GK checked !
!
ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaTGK
!_________________________________________________________________
!
SUBROUTINE Rablk(l,k,rCabT)
!
! Computes the FLR corrections to GK Sugama test-part
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! Local vars.
!
INTEGER :: p, j,g,h,r,s,s1,h1
TYPE(FM),SAVE :: Lrss1,Lghh1,T4inv1,T4inv2,NX0,NX1
TYPE(FM), SAVE :: alpha0_,NX
CHARACTER(len=T4len) :: T4str_
INTEGER :: icol
!
!
CALL FM_ENTER_USER_ROUTINE
!
!
NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
NX = - TO_FM(4)/SQRT(PIFM)*bafm_**2/4
!
alpha0_ = TO_FM(0)
!
ploop: DO p=0, Pmaxa
NX1 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
jloop: DO j=0,Jmaxa
icol = numba(p,j)
gloop: DO g=0, l + 2*k
hloop: DO h=0, k + floor(real(l)/2)
!______________________________________________
! INVERSE BASIS TRANSFORMATION
T4str_ = GET_T4(g,h,l,k) ! ... T^{lk}_{gh}
T4inv1 = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*TO_FM(2*g + 1)/2/FACTORIAL(TO_FM(2*g+2*h +1)/2)*TO_FM(T4str_)
IF(T4inv1 .ne. TO_FM(0)) THEN
rloop: DO r=0, p + 2*j
sloop: DO s=0, j + floor(real(p)/2)
T4str_ = GET_T4(r,s,p,j) ! ... T^{pj}_{rs}
T4inv2 = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(s))*TO_FM(2*r + 1)/2/FACTORIAL(TO_FM(2*r+2*s +1)/2)*TO_FM(T4str_)
IF(T4inv2 .ne. TO_FM(0)) THEN
s1loop: DO s1=0, s
Lrss1 = Lpjl(real(r,dp),real(s,dp),real(s1,dp))
h1loop: DO h1=0,h
Lghh1 = Lpjl(real(g,dp),real(h,dp),real(h1,dp))
! compute speed functions
alpha0_ = TO_FM(0)
IF( r .eq. g) THEN
alpha0_ = (1/TO_FM(2*g+1) )*( nuabparak(g + 2 + s1 +h1) + nuabDk( g+2 +s1 +h1 ))
ENDIF
alpha0_ = alpha0_ + X2prodLegendre(r,g)/2*(nuabDk( (r +g)/2 +2 +s1 +h1 ) - nuabparak((g+r)/2 + 2 + s1 +h1))
rCabT(icol) = rCabT(icol) + NX*NX0*NX1*T4inv1*T4inv2*Lrss1*Lghh1*alpha0_
ENDDO h1loop
ENDDO s1loop
ENDIF
ENDDO sloop
ENDDO rloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDDO jloop
ENDDO ploop
!
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Rablk
!_________________________________________________________________
!
SUBROUTINE UparalkcurlyRabpara1(l,k,rCabT)
!
! Computes \mathcal{U}_{\parallel}^{lk}\mathcal{R}_{ab}^{\parallel 1}
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! Local vars.
integer ::p,j,n,s,h,h1,icol
TYPE(FM), SAVE :: NX0,L1hh1,d0njs,T4inv,curlyUpara_,kern_,nupara_,thetaabm1_,NX
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
thetaabm1_ = (thetaabfm_ -TO_FM(1))
!
NX = - 8*(1 + chifm_**2)/3/SQRT(PIFM)*thetaabm1_
!
IF( l .eq. 1) THEN
!
! FLR Effects
curlyUpara_ = curlyUparalk(l,k)
ploop:DO p=0,Pmaxa
NX0 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
jloop: DO j=0,Jmaxa
icol = numba(p,j)
nloop: DO n=0,naFLR
kern_ = kerneln(bafm_,n)
sloop: DO s=0,n + j
hloop: DO h=0, s + floor(real(p)/2)
T4str_ = GET_T4(1,h,p,s) ! ...
T4inv = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(h))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*h +1)/2)*TO_FM(T4str_)
IF(T4inv .ne. TO_FM(0)) THEN
d0njs = ALLX2L(n,j,s,0)
h1loop: DO h1=0,h
L1hh1 = Lpjl(real(1,dp),real(h,dp),real(h1,dp))
nupara_ = nuabparak(3+h1)
rCabT(icol) = rCabT(icol) + NX*L1hh1*T4inv*d0njs*NX0*nupara_*curlyUpara_*kern_
ENDDO h1loop
ENDIF
ENDDO hloop
ENDDO sloop
ENDDO nloop
ENDDO jloop
ENDDO ploop
!
!
ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE UparalkcurlyRabpara1
!_________________________________________________________________
!
SUBROUTINE UperplkcurlyRabpara2(l,k,rCabT)
!
! Computes \mathcal{U}_{\perp}^{lk}\mathcal{R}_{ab}^{\parallel 2}
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! Local vars.
integer ::p,j,n,s,h,h1,icol
TYPE(FM), SAVE :: NX0,L0hh1,d1njs,T4inv,kern_,nupara_,thetaabm1_,curlyUperp,NX
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
!
thetaabm1_ = (thetaabfm_ -TO_FM(1))
!
NX = - 8*(1 + chifm_*chifm_)/SQRT(PIFM)*thetaabm1_
!
!
!
IF ( l .eq. 0 ) THEN
!
curlyUperp = curlyUperplk(l,k)
!
ploop:DO p=0,Pmaxa
NX0 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
jloop: DO j=0,Jmaxa
icol = numba(p,j)
nloop: DO n=0,naFLR
kern_ = kerneln(bafm_,n)
sloop: DO s=0,n + j +1
hloop: DO h=0, s + floor(real(p)/2)
T4str_ = GET_T4(0,h,p,s) ! ...
T4inv = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(h))*TO_FM( 1)/2/FACTORIAL(TO_FM(2*h +1)/2)*TO_FM(T4str_)
IF(T4inv .ne. TO_FM(0)) THEN
d1njs = ALLX2L(n,j,s,1)
h1loop: DO h1=0,h
L0hh1 = Lpjl(real(0,dp),real(h,dp),real(h1,dp))
nupara_ = nuabparak(2+h1)
rCabT(icol) = rCabT(icol) + NX*bafm_/(n+1)/2*L0hh1*T4inv*d1njs*NX0*nupara_*kern_*curlyUperp
!
ENDDO h1loop
ENDIF
ENDDO hloop
ENDDO sloop
ENDDO nloop
ENDDO jloop
ENDDO ploop
!
ENDIF
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE UperplkcurlyRabpara2
!_________________________________________________________________
!
SUBROUTINE TperplkcurlyRabpara3(l,k,rCabT)
!
! Computes \mathcal{T}_{\perp}^{lk} \mathcal{R}_{ab}^{\parallel 3}
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! Local vars.
!
integer ::p,j,n,s,h,h1,icol
TYPE(FM), SAVE :: curlyTlk, NX0,L0hh1,d0njs,T4inv,kern_,thetaabm1_,NX
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
thetaabm1_ = thetaabfm_ -TO_FM(1)
!
IF( l .eq. 2 .or. l .eq. 0 ) THEN
curlyTlk = curlyTperplk(l,k)
!
ploop:DO p=0,Pmaxa
NX0 = 1/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
jloop: DO j=0,Jmaxa
icol = numba(p,j)
nloop: DO n=0,naFLR
kern_ = kerneln(bafm_,n)
sloop: DO s=0,n + j
hloop: DO h=0, s + floor(real(p)/2)
T4str_ = GET_T4(0,h,p,s) ! ...
T4inv = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(h))*TO_FM( 1)/2/FACTORIAL(TO_FM(2*h +1)/2)*TO_FM(T4str_)
IF(T4inv .ne. TO_FM(0)) THEN
d0njs = ALLX2L(n,j,s,0)
h1loop: DO h1=0,h
L0hh1 = Lpjl(real(0,dp),real(h,dp),real(h1,dp))
rCabT(icol) = rCabT(icol) -thetaabm1_*8*L0hh1*d0njs*NX0*T4inv*kern_*2*(Erfkab_array(h1) - chifm_*(1 + chifm_**2)*ekab_array(h1+1))/chifm_**2/3/SQRT(PIFM)*curlyTlk
!
ENDDO h1loop
ENDIF
ENDDO hloop
ENDDO sloop
ENDDO nloop
ENDDO jloop
ENDDO ploop
!
ENDIF
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE TperplkcurlyRabpara3
!_________________________________________________________________
!
SUBROUTINE Fab0paralkupara(l,k,rCabT)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! Local vars.
TYPE(FM),SAVE :: kern_,Fasb0paralk_,NX0
INTEGER :: n,icol
!
CALL FM_ENTER_USER_ROUTINE
!
Fasb0paralk_ = curlyF0paralk(l,k)
!
NX0 = thetaabfm_ - 1
!
nloop: DO n=0,naFLR
kern_ = kerneln(bafm_,n)
rCabT(numba(1,n)) = rCabT(numba(1,n)) + kern_/SQRT(TO_FM(2))*NX0*Fasb0paralk_
ENDDO nloop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Fab0paralkupara
!_________________________________________________________________
!
SUBROUTINE Fab1paralkuperp(l,k,rCabT)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! local vars.
INTEGER :: n
TYPE(FM),SAVE :: thetaabm1_,Fab1paralk
!
CALL FM_ENTER_USER_ROUTINE
!
thetaabm1_ = thetaabfm_ -TO_FM(1)
!
Fab1paralk = curlyF1paralk(l,k)
!
DO n=0,naFLR
rCabT(numba(0,n)) = rCabT(numba(0,n)) + bafm_/2*kerneln(bafm_,n)*thetaabm1_*Fab1paralk
IF( n+1 .le. naFLR) THEN
rCabT(numba(0,n+1)) = rCabT(numba(0,n+1)) - bafm_/2*kerneln(bafm_,n)*thetaabm1_*Fab1paralk
ENDIF
ENDDO
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Fab1paralkuperp
!_________________________________________________________________
!
SUBROUTINE FabElkTa(l,k,rCabT)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabT
!
! local vars.
TYPE(FM),SAVE :: curlyFabElk_,thetaabm1_
INTEGER :: n
!
CALL FM_ENTER_USER_ROUTINE
!
curlyFabElk_ = curlyFElk(l,k)
!
thetaabm1_ = thetaabfm_ - TO_FM(1)
!
DO n=0,naFLR
rCabT(numba(0,n)) = rCabT(numba(0,n)) + TO_FM(4)/3*kerneln(bafm_,n)*TO_FM(n)*thetaabm1_*curlyFabElk_
rCabT(numba(2,n)) = rCabT(numba(2,n)) + SQRT(TO_FM(2))/3*kerneln(bafm_,n)*thetaabm1_*curlyFabElk_
IF(n-1 .ge. 0) THEN
rCabT(numba(0,n-1)) = rCabT(numba(0,n-1)) - TO_FM(2)/3*kerneln(bafm_,n)*TO_FM(n)*thetaabm1_*curlyFabElk_
ENDIF
IF( n+1 .le. naFLR) THEN
rCabT(numba(0,n+1)) = rCabT(numba(0,n+1)) - TO_FM(2)/3*kerneln(bafm_,n)*(n+1)*thetaabm1_*curlyFabElk_
ENDIF
ENDDO
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE FabElkTa
!_________________________________________________________________
!
SUBROUTINE Uparalkupara(l,k,rCabT)
!
! computes \mathcal{U}_{\parallel}^{lk} u_\parallel a
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:), INTENT(INOUT) :: rCabT
!
! local vars
!
INTEGER :: n
TYPE(FM), SAVE :: curlyUparalk_,thetaabm1_,NX0
!
CALL FM_ENTER_USER_ROUTINE
!
thetaabm1_ = (thetaabfm_ -TO_FM(1))
!
curlyUparalk_ = curlyUparalk(l,k)
!
NX0 = - 2*thetaabm1_**2*chifm_/SQRT(TO_FM(1 + chifm_**2))*4/SQRT(PIFM)/3
!
DO n=0,naFLR
rCabT(numba(1,n)) = rCabT(numba(1,n)) + NX0/SQRT(TO_FM(2))*curlyUparalk_*kerneln(bafm_,n)
ENDDO
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Uparalkupara
!_________________________________________________________________
!
SUBROUTINE Uperplkuperp(l,k,rCabT)
!
! computes \mathcal{U}_{\perp}^{lk} u_{\perp a}
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:), INTENT(INOUT) :: rCabT
!
! local vars
!
INTEGER :: n
TYPE(FM), SAVE :: curlyUperplk_,thetaabm1_,NX0
!
CALL FM_ENTER_USER_ROUTINE
!
thetaabm1_ = (thetaabfm_ -TO_FM(1))
!
NX0 = - 2*thetaabm1_**2*chifm_/SQRT(TO_FM(1 + chifm_**2))*4/3/SQRT(PIFM)
!
curlyUperplk_ = curlyUperplk(l,k)
!
nloop: DO n=0,naFLR
rCabT(numba(0,n)) = rCabT(numba(0,n)) + NX0*bafm_*curlyUperplk_/2*kerneln(bafm_,n)
IF(n+1 .le. naFLR) THEN
rCabT(numba(0,n+1)) = rCabT(numba(0,n+1)) - NX0*bafm_*curlyUperplk_/2*kerneln(bafm_,n)
ENDIF
ENDDO nloop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Uperplkuperp
!_________________________________________________________________
!
SUBROUTINE curlyTperpTa(l,k,rCabT)
!
! computes \mathcal{T}_{\perp}^{lk} T_{ a}
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:), INTENT(INOUT) :: rCabT
!
! local vars
!
INTEGER :: n
TYPE(FM), SAVE :: curlyTperplk_,thetaabm1_,NX0,NX1
!
CALL FM_ENTER_USER_ROUTINE
!
thetaabm1_ = (thetaabfm_ -TO_FM(1))
!
curlyTperplk_ = curlyTperplk(l,k)
!
NX1 = TO_FM(1)/(1 + chifm_**2)
NX0 = - 2*thetaabm1_**2*chifm_/SQRT(TO_FM(1 + chifm_**2))*4/SQRT(PIFM)/3
!
nloop: DO n=0,naFLR
rCabT(numba(2,n)) = rCabT(numba(2,n)) + NX1*NX0*SQRT(TO_FM(2))/3*kerneln(bafm_,n)*curlyTperplk_
rCabT(numba(0,n)) = rCabT(numba(0,n)) + TO_FM(n)*NX1*NX0*4*kerneln(bafm_,n)*curlyTperplk_/3
IF(n .ge. 1) THEN
rCabT(numba(0,n-1)) = rCabT(numba(0,n-1)) - NX1*NX0*n*2/3*kerneln(bafm_,n)*curlyTperplk_
ENDIF
IF(n + 1 .le. naFLR) THEN
rCabT(numba(0,n+1)) = rCabT(numba(0,n+1)) - NX1*NX0*(n+1)*kerneln(bafm_,n)*2/3*curlyTperplk_
ENDIF
ENDDO nloop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE curlyTperpTa
!_________________________________________________________________
!
!
! GYROKINETIC SUGAMAM FIELD (BACK-REACTION) PART
!
!
!_________________________________________________________________
SUBROUTINE SugamaFGK(l,k,rCabF)
!
! Computes the field (back-reaction) part of the gyrokinetic Sugama collision operator
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCabF
!
!
CALL overlinecurlyVparalkVparaab(l,k,rCabF) ! ... ok !
CALL overlinecurlyVperplkVperpaab(l,k,rCabF) ! .... ok !
CALL overlinecurlyWlkWab(l,k,rCabF) ! ...
!
END SUBROUTINE SugamaFGK
!_________________________________________________________________
!
SUBROUTINE overlinecurlyVparalkVparaab(l,k,rCabF)
!
! computes - \frac{1}{\sqrt{2^l l !}} \overline{\mathcal{V}}_{\parallel ab}^{lk} \overline{V}_{\parallel ab} (*)
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:), INTENT(INOUT) :: rCabF
!
! Local vars.
TYPE(FM), SAVE :: NX0,overlinecurlyVparalk_
!
CALL FM_ENTER_USER_ROUTINE
!
!
NX0 = - 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l))) ! ... includes (-1) factor from (*) [see above]
!
overlinecurlyVparalk_ = overlinecurlyVparalk(l,k)
!
rCabF(:) = rCabF(:) + NX0*overlinecurlyVparalk_*roverlineVparaab(:)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE overlinecurlyVparalkVparaab
!_________________________________________________________________
!
!
SUBROUTINE overlinecurlyVperplkVperpaab(l,k,rCabF)
!
! computes - \frac{1}{\sqrt{2^l l !}} \overline{\mathcal{V}}_{\perp ab}^{lk}\overline{V}_{\perp ab} \overline{W}_{ab} \frac{\overline{\mathcal{W}}_{ab}^{lk}}{\sqrt{2^l l!}},
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:), INTENT(INOUT) :: rCabF
!
! Local vars.
INTEGER:: p,j,n,s,icol
TYPE(FM),SAVE :: buffchi_,gammaab_,overlinecurlyVperplk_,Ibaps_
TYPE(FM),SAVE :: NX0,NX1,d1njs_,kern_,mbovergammaab_,NX
CALL FM_ENTER_USER_ROUTINE
!
!
overlinecurlyVperplk_ = overlinecurlyVperplk(l,k)
!
NX0 = - 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
rCabF(:) =rCabF(:) + NX0*overlinecurlyVperplk_*roverlineVperpab(:)
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE overlinecurlyVperplkVperpaab
!_________________________________________________________________
!
SUBROUTINE overlinecurlyWlkWab(l,k,rCabF)
!
! computes - \overline{W}_{ab} \frac{\overline{\mathcal{W}}_{ab}^{lk}}{\sqrt{2^l l!}},
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:), INTENT(INOUT) :: rCabF
!
! local vars.
!
INTEGER :: p,j,n,s,icol
TYPE(FM),SAVE :: overlinecurlyWab_, NX0,etaab_,NX1,d0njs_,curlyWbaps_,Tboveretaab_,NX
!
CALL FM_ENTER_USER_ROUTINE
!
overlinecurlyWab_ = overlinecurlyWlk(l,k)
!
NX0 = - 1/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
! Note : add additional 2 factor to match DK Sugama and to conserve energy in the limit k_\perp =0
!
rCabF(:) = rCabF(:) + 2*NX0*overlinecurlyWab_*roverlineWab(:)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE overlinecurlyWlkWab
!_________________________________________________________________
!
!
! GYROKINETIC SUGAMA SELF (TEST + BACK-REACTION)
!
!
!_________________________________________________________________
SUBROUTINE SugamaSelfGK(l,k,rCaa)
!
! computes the Gyrokinetic self-collision Sugama
!
IMPLICIT NONE
!
INTEGER, intent(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCaa
CALL FM_ENTER_USER_ROUTINE
!
!
CALL SugamaTGK(l,k,rCaa)
!!CALL SugamaFDK(l,k,rCaa)
CALL SugamaFGK(l,k,rCaa)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaSelfGK
!_________________________________________________________________
!
!
! GYROKINETIC IMPROVED FIELD SUGAMA
!
!
!_________________________________________________________________
SUBROUTINE SugamaImprovedTGK(l,k,rCabT)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),dimension(:),INTENT(INOUT) :: rCabT
!
CALL FM_ENTER_USER_ROUTINE
!
! CODE HERE ...... !!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedTGK
!
!__________________________________________
!
SUBROUTINE SugamaImprovedFGK(l,k,rCabF)
!
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),dimension(:),INTENT(INOUT) :: rCabF
!
CALL FM_ENTER_USER_ROUTINE
!
! CODE HERE ...... !!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaImprovedFGK
c4science · Help