Homec4science
Diffusion COSOlver (improved_sugama) Browse

cosolver/src/sugama_mod.f90
cd1286e2d410improved_sugama

Last Change
Raw File

sugama_mod.f90
View Options

MODULE Sugama
!
! Linearized DK/GK Sugamam Collision operator between species a and b
!
! 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 16, 11 (2009)
USE prec_const
USE basic
USE model
USE FMZM
USE coeff
USE basis_transformation
USE speed_functions
USE auxval
use coulomb_parts
IMPLICIT NONE
PUBLIC
CONTAINS
!_________________________________________________________________
!
!
! 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
! CPa[f]^lk
CALL SugamaTCPabDK(l,k,rCabT)
! PaC[f]^lk
CALL SugamaTPCabDK(l,k,rCabT)
! PaCPa[f]^lk
CALL SugamaTPCPabDK(l,k,rCabT)
!
rCabT(:) = TO_FM(2)*rCabT(:)
!
! CabS[f]^lk: isothermal part of test Coulomb
CALL Coulomb_isothermalTDK(l,k,rCabT)
!
END SUBROUTINE SugamaTDK
!__________________________________________________
!
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,JEMaxx
! ___________ 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,JEmaxx
! ___________ 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,JEmaxx
! ___________ 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 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 }
! - Added factor of 2 added to match the full DK Coulomb
!
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 = 2*thetaabfm_*thetabafm_
!
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)
! Momentum
rCabF(icol) = rCabF(icol) + TO_FM(2)*NX0*NX1*NX3*NX2*( TO_FM(1) + TO_FM(1)/(chifm_**2))*( TO_FM(1) + chifm_**2)*curlyVablk_*curlyVbapj_/chifm_
! Energy
rCabF(icol) = rCabF(icol) + TO_FM(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
!
!
! 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
!
!
rCabT(:) = rCabT(:)*TO_FM(2)
!
! Drift kinetic part: isothermal part of linearized Coulomb
CALL Coulomb_isothermalTDK(l,k,rCabT)
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaTGK
!
!_________________________________________________________________
!
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._dp
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} (*)
!
! Note : add additional 2 factor to match DK Sugama
!
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) ! ... checked
!
rCabF(:) = rCabF(:) + 2*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
!7
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:), INTENT(INOUT) :: rCabF
!
! Local vars.
INTEGER :: icol
TYPE(FM),SAVE :: overlinecurlyVperplk_
TYPE(FM),SAVE :: NX0
CALL FM_ENTER_USER_ROUTINE
!
!
overlinecurlyVperplk_ = overlinecurlyVperplk(l,k)
!
NX0 = -TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
DO icol=1,numbb(Pmaxb,Jmaxb)
rCabF(icol) =rCabF(icol) + 2*NX0*overlinecurlyVperplk_*roverlineVperpab(icol)
ENDDO
!
!
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!}},
!
! Note : add additional 2 factor to match DK Sugama
!
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)))
!
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 SugamaFGK(l,k,rCaa)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE SugamaSelfGK
!
END MODULE Sugama
c4science · Help