Homec4science
Diffusion COSOlver (master) Browse

cosolver/src/gkcoulomb_mod.f90
5b66a09e0f20master

Last Change
Raw File

gkcoulomb_mod.f90
View Options

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 gkcoulomb_coeffs
use coulomb
use T4T5
IMPLICIT NONE
PUBLIC
! pointer
PROCEDURE( coulombGK ) , POINTER :: coulombselfgk => null()
PROCEDURE( coulombGK ) , POINTER :: coulombTGK => null()
PROCEDURE( coulombGK ) , POINTER :: coulombFGK => null()
! pointer
PROCEDURE( tabljpmf_fct), POINTER :: tabljpmf
ABSTRACT INTERFACE
subroutine coulombgk( l,k,j,n,np,rCaa)
use fmzm
IMPLICIT NONE
INTEGER,INTENT(in) :: l,k
INTEGER,INTENT(in) :: j,n,np
TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
end subroutine coulombgk
END INTERFACE
ABSTRACT INTERFACE
function tabljpmf_fct(l,j,p,m,f) result(XOUT)
use fmzm
IMPLICIT NONE
INTEGER,intent(in) :: l,j,p,m,f
TYPE(FM):: XOUT
end function tabljpmf_fct
END INTERFACE
CONTAINS
!
subroutine set_gkcoulomb
! set pointer to routine
implicit none
IF( me .eq. 0) print*, 'Set gk coulomb ....'
IF( .not. GKC_USE_FULL) THEN
if( me .eq. 0) print*, 'Use prec-computed gk-coulomb ....'
IF( GKC_USE_OPT) THEN
! optimized version
coulombselfgk => coulombgk_opt_prec
coulombTGK => coulombgk_opt_prec
coulombFGK => coulombgk_opt_prec
ELSE ! not optimized and not working with the new version of the code (see n summation that is now contained in coulombgk_opt_prec
coulombselfgk => coulombselfgk_prec ! pre-computed
coulombTGK => coulombTGK_prec
coulombFGK => coulombFGK_prec
ENDIF
ELSE
if( me .eq. 0) print*, 'Use full gk-coulomb ....'
coulombselfgk => coulombgk_full
coulombTGK => coulombgk_full
coulombFGK => coulombgk_full
call init_T4T5(.false. , 0)
ENDIF
!
end subroutine set_gkcoulomb
!
SUBROUTINE coulombSelfGK_prec(l,k,j,n,np,rCaa)
!
! construct the like-species gyrokinetic Coulomb Operator.
!
! Note: - use pre-computed Tljpmf coefficients and CurlyNpjmn computed using gkcoulomb_suite
!
!
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, mmax , n_, np_
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_, prodkern_
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=FMlen) :: 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)
prodkern_ = kern_ * kernp_
n_ = n
np_ = 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)
mmax = p
IF( LWL .and. p .gt. 2) THEN
mmax = 2
ENDIF
mloop: DO m=0, mmax ! p ! ... collisional harmonics
IF( m .eq. 0) THEN
bapowm = TO_FM(1)
ELSE
bapowm = 2._dp*TO_FM(bafm_/TO_FM(2))**(2*m)
ENDIF
floop: DO f=0,k+n_
bardmnkf_ = ALL2L(n_,k,f,m)
IF( bardmnkf_ .ne. ZEROFM) THEN
IF( p .le. l+2*f+m) THEN
Tljpmfstr_ = Tljpmf_STR( f, m, p, j, l)
Tljpmf_ = TO_FM(Tljpmfstr_)
IF(Tljpmf_ .ne. ZEROFM) THEN
DO icol =1,numba(Pmaxa,Jmaxa)
! column major
RowNpjstr_ = curlyNpjmn_STR(icol, np_, m, j, p)
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* prodkern_
ENDDO
ENDIF
ENDIF ! ... if p <= l + 2*f + m
ENDIF
ENDDO floop
ENDDO mloop
ENDDO ploop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombSelfGK_prec
!
!___________________________________________________________________________
!
SUBROUTINE coulombgk_opt_prec(l, k, j, n, np, rCab)
!
! construct the Gyrokinetic Coulomb Operator between a and b
!
! Note: - use pre-computed Tljpmf coefficients and CurlyNpjmn computed using gkcoulomb_suite
! - see L. Lechot TPIV reports (not yet completed)
!
IMPLICIT NONE
!
! input vars.
INTEGER, INTENT(in) :: l,k
INTEGER, INTENT(in) :: j, np, n
TYPE(FM), DIMENSION(:), INTENT(inout):: rCab
! local variables
INTEGER :: p,f,m,r, icol, mmax , n_, np_
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=FMlen) :: Tljpmfstr_, RowNpjstr_
! to store temp. vars
TYPE(FM), save :: tmp_
CHARACTER(len= FMlen) :: ALL2Lstr_
!
CALL FM_ENTER_USER_ROUTINE
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
kernp_ = kerneln(bpfm_,np)
np_ = 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)
mmax = p
mloop_: DO m=0, mmax
IF( m .eq. 0) THEN
bapowm = TO_FM(1)
ELSE
bapowm = 2._dp* (TO_FM(bfm_/2)**TO_FM(m)) * (TO_FM(bpfm_/2)**TO_FM(m))
ENDIF!
tmp_ = TO_FM(0)
do n_ = 0, naFLR ! this summation was orignially contained outside
kern_ = kerneln(bfm_,n_)
floop_: DO f=0,k+n_
IF( p .le. l+2*f+m) THEN ! use dot product
ALL2Lstr_ = array_ALL2L(f,m,k,n_)
bardmnkf_ = TO_FM(ALL2Lstr_)
IF( bardmnkf_ .ne. ZEROFM) THEN
Tljpmfstr_ = Tljpmf_STR( f, m, p, j, l)
Tljpmf_ = TO_FM(Tljpmfstr_)
tmp_ = tmp_ + bardmnkf_*Tljpmf_ * &
kern_* FACTORIAL(TO_FM(n_))/FACTORIAL(TO_FM(n_ +m))
endif
endif
enddo floop_
enddo
DO icol =1,numba(Pmaxa,Jmaxa)
RowNpjstr_ = curlyNpjmn_STR(icol, np_, m, j, p)
RowNpj_ = TO_FM(RowNpjstr_)
rCab(icol) = rCab(icol) + NX0*RowNpj_/FACTORIAL(TO_FM(np_ +m))/SPJ*NX1*FACTORIAL(TO_FM(np_))*tmp_ *bapowm * kernp_
ENDDO
ENDDO mloop_
ENDDO ploop_
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombgk_opt_prec
!
!___________________________________________________________________________
!
SUBROUTINE coulombGK_full(l,k,j,n,np,rCaa)
!
! construct the like-species GK Coulomb Operator.
!
!
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, mmax
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_, T5_,X1
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,Lptd_,Tljpmf_
!
CALL FM_ENTER_USER_ROUTINE
!
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)
mmax = p
IF( LWL .and. p .gt. 2) THEN
mmax = 2
ENDIF
mloop: DO m=0, mmax ! ... 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( bardmnkf_ .ne. ZEROFM) THEN
IF( p .le. l+2*f+m) THEN
Tljpmf_ = Tabljpmf(l,j,p,m,f)
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,np + m + h)
! Change order : n + m + h <=Jmax
dmnhs1_ = ALLX2L(np,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_
ENDDO
ENDIF
ENDDO
ENDDO
ENDIF
ENDIF
ENDIF
ENDDO floop
ENDDO mloop
ENDDO ploop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombGK_full
!
!___________________________________________________________________________
!
!
SUBROUTINE coulombTGK_prec(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 , mmax
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_, prodkern_
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm, Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=FMlen) :: 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)
prodkern_ = kern_ * kernp_
!
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)
mmax = p
IF( LWL .and. p .gt. 2) THEN
mmax = 2
ENDIF
mloop: DO m=0, mmax ! ... 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
Tljpmfstr_ = Tljpmf_STR( f, m, p, j, l)
Tljpmf_ = TO_FM( Tljpmfstr_)
IF(Tljpmf_ .ne. ZEROFM) THEN
DO icol =1,numba(Pmaxa,Jmaxa)
RowNpjstr_ = curlyNpjmn_STR(icol, np, m, j, p)
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*prodkern_
ENDDO
ENDIF
ENDIF ! ... if p < l + 2*f + m
ENDDO floop
ENDDO mloop
ENDDO ploop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombTGK_prec
!
!___________________________________________________________________________
!
SUBROUTINE coulombFGK_prec(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, mmax
TYPE(FM),SAVE :: NX0,NX1,SPJ,kern_,kernp_, prodkern_
TYPE(FM), SAVE :: bardmnkf_,dmnhs1_,bapowm,bbpowm, Lptd_,Tljpmf_,RowNpj_
CHARACTER(len=FMlen) :: 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)
prodkern_ = kern_ * kernp_
!
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)
mmax = p
IF( LWL .and. p .gt. 2) THEN
mmax = 2
ENDIF
mloop: DO m=0, mmax
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
Tljpmfstr_ = Tljpmf_STR( f, m, p, j, l)
Tljpmf_ = TO_FM(Tljpmfstr_)
IF(Tljpmf_ .ne. ZEROFM) THEN
DO icol =1,numbb(Pmaxb,Jmaxb)
RowNpjstr_ = curlyNpjmn_STR(icol, np, m, j, p)
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*prodkern_
ENDDO
ENDIF
ENDIF ! ... if p < l + 2*f + m
ENDDO floop
ENDDO mloop
ENDDO ploop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombFGK_prec
!
!___________________________________________________________________________
!
! routines needed to compute gk coulomb on the fly
!
SUBROUTINE set_tabljpmf(coll_type)
! set Tabljpmf
implicit none
character(1), intent(in ) :: coll_type
select case( coll_type)
case('s') ! self
tabljpmf => Tabljpmf_self
case('t') ! test
tabljpmf => Tabljpmf_test
case('f') ! field
tabljpmf => Tabljpmf_field
end select
!
END SUBROUTINE
FUNCTION Tabljpmf_self(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)
!T5inv = TO_FM(2)
IF(T5inv .ne. ZEROFM) THEN
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)
XOUT = XOUT + Lptd*( nuT_ + nuF_)*T5inv
enddo
ENDIF
enddo
!
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION Tabljpmf_self
FUNCTION Tabljpmf_test(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, 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))
nuT_ = nusTabpjgd(p,j,p,d)
XOUT = XOUT + Lptd* nuT_ *T5inv
enddo
ENDIF
enddo
!
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION Tabljpmf_test
FUNCTION Tabljpmf_field(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_
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))
nuF_ = nusFabpjgd(p,j,p,d)
XOUT = XOUT + Lptd* nuF_ *T5inv
enddo
ENDIF
enddo
!
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION Tabljpmf_field
! 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=FMlen) :: 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=FMlen) :: 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=FMlen) :: 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
!
!___________________________________________________________________
!
!
END MODULE gkcoulomb
c4science · Help