Homec4science
Diffusion COSOlver (new_coulomb) Browse

cosolver/src/improved_sugama_mod.f90
83b8f38f6c51new_coulomb

Last Change
Raw File

improved_sugama_mod.f90
View Options

MODULE improved_sugama
!
! 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 coulomb
use sugama
IMPLICIT NONE
PUBLIC
CONTAINS
SUBROUTINE evaluate_CurlIpjk
! computes CurlIparaapjk, CurlIperpapjk, CurlIparabpjk, CurlIperpbpjk
! for p=0,p=MAX(Pmaxe,Pmaxi)
! for j=0,j=MAX(Jmaxe,Jmaxi)
! for k=0,k=max(ImpSugamaJmax, ImpSugamaKmax)
integer :: p,j,k
character(len=CurlI_IS_len) :: val_str
call FM_ENTER_USER_ROUTINE
do j=0,MAX(Jmaxe,Jmaxi)
do p=0,MAX(Pmaxe,Pmaxi)
do k=0,MAX(ImpSugamaJmax, ImpSugamaKmax)
CALL FMFORM('F60.40',CurlIparapjk(p,j,k,bafm_),val_str)
! CurlIparaapjk_array_str(p,j,k) = val_str
CurlIparaapjk_array_str(k, p,j) = val_str
!! print '(" - CurlI perp a p="I2"/"I2" j="I2"/"I2" k="I2"/"I2)', p, Pmaxa, j, Jmaxa, k, MAX(ImpSugamaJmax, ImpSugamaKmax)
CALL FMFORM('F60.40',CurlIperppjk(p,j,k,bafm_),val_str)
! CurlIperpapjk_array_str(p,j,k) = val_str
CurlIperpapjk_array_str(k, p,j) = val_str
!print '(" - CurlI para b p="I2"/"I2" j="I2"/"I2" k="I2"/"I2)', p, Pmaxa, j, Jmaxa, k, MAX(ImpSugamaJmax, ImpSugamaKmax)
CALL FMFORM('F60.40',CurlIparapjk(p,j,k,bbfm_),val_str)
! CurlIparabpjk_array_str(p,j, k ) = val_str
CurlIparabpjk_array_str(k, p,j ) = val_str
!print '(" - CurlI perp b p="I2"/"I2" j="I2"/"I2" k="I2"/"I2)', p, Pmaxa, j, Jmaxa, k, MAX(ImpSugamaJmax, ImpSugamaKmax)
CALL FMFORM('F60.40',CurlIperppjk(p,j,k,bbfm_),val_str)
! CurlIperpbpjk_array_str(p,j,k) = val_str
CurlIperpbpjk_array_str(k, p,j) = val_str
end do
end do
end do
call FM_EXIT_USER_ROUTINE
!
END SUBROUTINE evaluate_CurlIpjk
! Friction Coeffs
! SUBROUTINE compute_friction_coefficients()
! IMPLICIT NONE
! ! Local vars
! INTEGER, PARAMETER :: imax = 5, jmax = 5
! INTEGER :: i, j
! TYPE(FM) :: NX0
! ! Contains 1/tauS
! TYPE(FM), DIMENSION(0:imax, 0:jmax) :: Mee, Mii, Mei, Mie, MSee, MSii, MSei, MSie
! TYPE(FM), DIMENSION(0:imax, 0:jmax) :: Nee, Nii, Nei, Nie, NSee, NSii, NSei, NSie
! CALL FM_ENTER_USER_ROUTINE
! NX0 = TO_FM(4)/TO_FM(3)/SQRT(PIFM)
! WRITE(*,*) '============= Friction coefficients ============='
! PRINT *, '===== Precomputing ====='
! print *, ' - Electron-Electron'
! CALL set_model('ee')
! CALL set_basic('ee')
! CALL evaluate_auxval_field
! DO i=0, imax
! DO j=0, jmax
! Mee(i,j) = MabLlk(i,j)*NX0
! Nee(i,j) = NabLlk(i,j)*NX0
! MSee(i,j) = MabSlk(i,j)*NX0
! NSee(i,j) = NabSlk(i,j)*NX0
! END DO
! END DO
! print *, ' - Ion-Ion'
! CALL set_model('ii')
! CALL set_basic('ii')
! CALL evaluate_auxval_field
! DO i=0, imax
! DO j=0, jmax
! Mii(i,j) = MabLlk(i,j)*NX0
! Nii(i,j) = NabLlk(i,j)*NX0
! MSii(i,j) = MabSlk(i,j)*NX0
! NSii(i,j) = NabSlk(i,j)*NX0
! END DO
! END DO
! print *, ' - Electron-Ion'
! CALL set_model('ei')
! CALL set_basic('ei')
! CALL evaluate_auxval_field
! DO i=0, imax
! DO j=0, jmax
! Mei(i,j) = MabLlk(i,j)*NX0
! Nei(i,j) = NabLlk(i,j)*NX0
! MSei(i,j) = MabSlk(i,j)*NX0
! NSei(i,j) = NabSlk(i,j)*NX0
! END DO
! END DO
! print *, ' - Ion-Electron'
! CALL set_model('ie')
! CALL set_basic('ie')
! CALL evaluate_auxval_field
! DO i=0, imax
! DO j=0, jmax
! Mie(i,j) = MabLlk(i,j)*NX0
! Nie(i,j) = NabLlk(i,j)*NX0
! MSie(i,j) = MabSlk(i,j)*NX0
! NSie(i,j) = NabSlk(i,j)*NX0
! END DO
! END DO
! PRINT *, '===== Computing coeffs ====='
! print *, ' - Electron-Electron'
! CALL set_model('ee')
! CALL set_basic('ee')
! CALL evaluate_auxval_field
! DO i = 1, imax
! DO j= 1, jmax
! PRINT '("l_ee_"I1"_"I1" = "E17.10)', i, j, TO_DP(Mee(i-1,j-1) + Mei(i-1,j-1) + Nee(i-1,j-1))
! PRINT '("lS_ee_"I1"_"I1" = "E17.10)', i, j, TO_DP(MSee(i-1,j-1) + MSei(i-1,j-1) + NSee(i-1,j-1))
! END DO
! END DO
! print *, ' - Ion-Ion'
! CALL set_model('ii')
! CALL set_basic('ii')
! CALL evaluate_auxval_field
! DO i = 1, imax
! DO j= 1, jmax
! PRINT '("l_ii_"I1"_"I1" = "E17.10)', i, j, TO_DP(Mii(i-1,j-1) + Mie(i-1,j-1) + Nii(i-1,j-1))
! PRINT '("lS_ii_"I1"_"I1" = "E17.10)', i, j, TO_DP(MSii(i-1,j-1) + MSie(i-1,j-1) + NSii(i-1,j-1))
! END DO
! END DO
! print *, ' - Electron-Ion'
! CALL set_model('ei')
! CALL set_basic('ei')
! CALL evaluate_auxval_field
! DO i = 1, imax
! DO j= 1, jmax
! PRINT '("l_ei_"I1"_"I1" = "E17.10)', i, j, TO_DP(Nei(i-1,j-1))
! PRINT '("lS_ei_"I1"_"I1" = "E17.10)', i, j, TO_DP(NSei(i-1,j-1))
! END DO
! END DO
! print *, ' - Ion-Electron'
! CALL set_model('ie')
! CALL set_basic('ie')
! CALL evaluate_auxval_field
! DO i = 1, imax
! DO j= 1, jmax
! PRINT '("l_ie_"I1"_"I1" = "E17.10)', i, j, TO_DP(Nie(i-1,j-1))
! PRINT '("lS_ie_"I1"_"I1" = "E17.10)', i, j, TO_DP(NSie(i-1,j-1))
! END DO
! END DO
! WRITE(*,*) '============= End Friction coefficients ============='
! CALL FM_EXIT_USER_ROUTINE
! END SUBROUTINE compute_friction_coefficients
! 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 = 20, kmax = 20
! 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_ab
! 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_ab
! 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: "I2"/"I2", "I2"/"I2)', l, lmax, k, kmax
! PRINT '("Mab_"I2"_"I2" - Mab_"I2"_"I2" = "E20.10)', l, k, k, l, TO_DP(Mab(l,k) - Mab(k,l))
! PRINT '(">Nab_"I2"_"I2" = "E20.10)', l, k, TO_DP(Nab(l,k))
! PRINT '(">Nba_"I2"_"I2" = "E20.10)', k, l, TO_DP(Nba(k,l))
! PRINT '("Nab_"I2"_"I2" - tau*chi*Nba_"I2"_"I2" = "E20.10)', 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: "I2"/"I2", "I2"/"I2)', l, lmax, k, kmax
! PRINT '(">NSab_"I2"_"I2" = "E20.10)', l, k, TO_DP(NabS(l,k))
! PRINT '("MSab_"I2"_"I2" - Mab_"I2"_"I2" = "E20.10)', l, k, l, k, TO_DP(MabS(l,k) - Mab(k,l))
! PRINT '("NSab_"I2"_"I2" - Nab_"I2"_0*Nab_0_"I2"/Nab_0_0 = "E20.10)', 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: "I2"/"I2", "I2"/"I2)', l, lmax, 0, kmax
! PRINT '("NSab_"I2"_0 - Nab_"I2"_0 = "E20.10)', l, l, TO_DP(NabS(l,0) - Nab(l,0))
! END DO lloop2_1
! lloop2_2: DO k=0, kmax
! PRINT '("Loop: "I2"/"I2", "I2"/"I2)', 0, lmax, k, kmax
! PRINT '("NSab_0_"I2" - Nab_0_"I2" = "E20.10)', k, k, TO_DP(NabS(0,k) - Nab(0,k))
! END DO lloop2_2
! END IF
! WRITE(*,*) "# Temperature Ta != Tb"
! PRINT '("Loop: "I2"/"I2", "I2"/"I2)', 0, lmax, 0, kmax
! PRINT '("Mab_0_0 + Nab_0_0 ="E20.10)', TO_DP(Mab(0,0) + Nab(0,0))
! PRINT '("Mab_0_0 - MSab_0_0 ="E20.10)', TO_DP(Mab(0,0) - MabS(0,0))
! PRINT '("Mab_0_0 + NSab_0_0 ="E20.10)', TO_DP(Mab(0,0) + NabS(0,0))
! PRINT '("Nab_0_0 + MSab_0_0 ="E20.10)', TO_DP(Nab(0,0) + MabS(0,0))
! PRINT '("Nab_0_0 - NSab_0_0 ="E20.10)', TO_DP(Nab(0,0) - NabS(0,0))
! PRINT '("MSab_0_0 + NSab_0_0 ="E20.10)', TO_DP(MabS(0,0) + NabS(0,0))
! kloop2_3: DO k=0, kmax
! PRINT '("Loop: "I2"/"I2", "I2"/"I2)', 0, lmax, k, kmax
! PRINT '("MSab_0_"I2" + tau*chi*NSba_0_"I2" = "E20.10)', k, k, TO_DP(MabS(0,k) + taufm_*chifm_*NbaS(0,k))
! PRINT '("Mab_0_"I2" + tau*chi*Nba_0_"I2" = "E20.10)', 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: "I2"/"I2", "I2"/"I2)', l, lmax, k, kmax
! PRINT '("MSab_"I2"_"I2" - MSab_"I2"_"I2" = "E20.10)', l, k, k, l, TO_DP(MabS(l,k) - MabS(k,l))
! PRINT '("NSab_"I2"_"I2" - tau*chi*NSba_"I2"_"I2" = "E20.10)', l, k, k, l, TO_DP(NabS(l,k) - taufm_*chifm_*NbaS(k,l))
! PRINT '("NSab_"I2"_"I2" - MSab_0_"I2"*NSab_0_"I2"/Mab_0_0 = "E20.10)', l,k, l, k, TO_DP(NabS(l,k) - MabS(0,l)*NabS(0,k)/Mab(0,0))
! PRINT '("NSab_"I2"_"I2" - NSab_"I2"_0*NSab_0_"I2"/Nab_0_0 = "E20.10)', 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: "I2"/"I2", "I2"/"I2)', l, lmax, k, kmax
! PRINT '("DeltaMab_"I2"_"I2" = "E20.10)', l, k, TO_DP(DeltaMab(l,k))
! PRINT '("DeltaNab_"I2"_"I2" - (Nab_0_0*Nab_"I2"_"I2"-Nab_"I2"_0*Nab_0_"I2")/Nab_0_0 = "E20.10)', 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: "I2"/"I2", "I2"/"I2)', 0, lmax, 0, kmax
! PRINT '("DeltaNab_0_0 ="E20.10)', TO_DP(DeltaNab(0,0))
! lloop3_1: DO l=0, lmax
! PRINT '("Loop: "I2"/"I2", "I2"/"I2)', l, lmax, 0, kmax
! PRINT '("DeltaNab_"I2"_0 ="E20.10)', l, TO_DP(DeltaNab(l,0))
! END DO lloop3_1
! lloop3_2: DO k=0, kmax
! PRINT '("Loop: "I2"/"I2", "I2"/"I2)', 0, lmax, k, kmax
! PRINT '("DeltaNab_0_"I2" ="E20.10)', 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
IF( ADD_IS_TO_OS) CALL SugamaTDK(l,k, rCabT)
! add correction matrix elements \Delta C_{ab}^{FLS}
IF( l .ge. 1 .or. k .ge. 1) THEN
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
jloop: DO j=0, MIN(k + floor(real(l)/2), ImpSugamaJmax)
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, ImpSugamaKmax
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
IF( ADD_IS_TO_OS ) CALL SugamaFDK(l, k, rCabF)
!
! add correction matrix elements \Delta C_{ab}^{FLS}
IF( l .ge. 1 .or. k .ge. 1) THEN
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
!
jloop: DO j=0, MIN(k + floor(real(l)/2), ImpSugamaJmax)
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,ImpSugamaKmax
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) / chifm_ ! added 1 / chifm_ term
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
IF( addtest) CALL ImprovedSugamaTDK(l,k,rCaa)
IF( addback) CALL ImprovedSugamaFDK(l,k,rCaa)
!
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) :: DeltaMabnt_, IPara_a_lkn_, IPerp_a_lkn_
TYPE(FM) :: IPara_a_pjt_, IPerp_a_pjt_
!
CALL FM_ENTER_USER_ROUTINE
! Construct GK Test Sugama
IF( ADD_IS_TO_OS ) 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, ImpSugamaJmax
! IPara_a_lkn_ = TO_FM(CurlIparaapjk_array_str(l,k,n))
! IPerp_a_lkn_ = TO_FM(CurlIperpapjk_array_str(l,k,n))
IPara_a_lkn_ = TO_FM(CurlIparaapjk_array_str(n, l,k))
IPerp_a_lkn_ = TO_FM(CurlIperpapjk_array_str(n, l,k))
cn_ = CJ(n)
jloop: DO j=0, Jmaxa
ploop: DO p=0, Pmaxa
tloop: DO t=0, ImpSugamaKmax
DeltaMabnt_ = DeltaMlk(n,t)
ct_ = CJ(t)
! IPara_a_pjt_ = TO_FM(CurlIparaapjk_array_str(p,j,t))
! IPerp_a_pjt_ = TO_FM(CurlIperpapjk_array_str(p,j,t))
IPara_a_pjt_ = TO_FM(CurlIparaapjk_array_str(t,p,j))
IPerp_a_pjt_ = TO_FM(CurlIperpapjk_array_str(t,p ,j))
icol = numba(p,j)
rCabT(icol) = rCabT(icol) + NX0 *cn_*ct_ * DeltaMabnt_ *(IPara_a_pjt_*IPara_a_lkn_ + IPerp_a_pjt_*IPerp_a_lkn_)
ENDDO tloop
ENDDO ploop
ENDDO jloop
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
IF( ADD_IS_TO_OS) CALL SugamaFGK(l,k,rCabF)
! Add correction matrix \Delta C_{ab}^{FLS}
NX0 = TO_FM(8)/TO_FM(3)/SQRT(PIFM)
nloop: DO n=0, ImpSugamaJmax
! IPara_a_lkn_ = TO_FM(CurlIparaapjk_array_str(l,k,n))
! IPerp_a_lkn_ = TO_FM(CurlIperpapjk_array_str(l,k,n))
IPara_a_lkn_ = TO_FM(CurlIparaapjk_array_str(n, l,k))
IPerp_a_lkn_ = TO_FM(CurlIperpapjk_array_str(n, l,k))
cn_ = CJ(n)
jloop: DO j=0, Jmaxa
ploop: DO p=0, Pmaxa
tloop: DO t=0, ImpSugamaKmax
DeltaNabnt_ = DeltaNlk(n,t)
ct_ = CJ(t)
! IPara_b_pjt_ = TO_FM(CurlIparabpjk_array_str(p,j,t))
! IPerp_b_pjt_ = TO_FM(CurlIperpbpjk_array_str(p,j,t))
IPara_b_pjt_ = TO_FM(CurlIparabpjk_array_str(t, p,j))
IPerp_b_pjt_ = TO_FM(CurlIperpbpjk_array_str(t, p,j))
icol = numbb(p,j)
rCabF(icol) = rCabF(icol) + NX0*cn_*ct_*DeltaNabnt_*(IPara_b_pjt_*IPara_a_lkn_ + IPerp_b_pjt_*IPerp_a_lkn_)/chifm_
ENDDO tloop
ENDDO ploop
ENDDO jloop
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
CALL ImprovedSugamaTGK(l,k,rCaa)
CALL ImprovedSugamaFGK(l,k,rCaa)
CALL FM_EXIT_USER_ROUTINE
END SUBROUTINE ImprovedSugamaSelfGK
END MODULE improved_sugama
c4science · Help