Homec4science
Diffusion COSOlver (improved_sugama) Browse

cosolver/src/coulomb_pol_mod.f90
cd1286e2d410improved_sugama

Last Change
Raw File

coulomb_pol_mod.f90
View Options

module coulomb_pol
!
! Gyrokinetic polarization terms for the gyrokinetic Coulomb collision operator
!
USE prec_const
USE basic
USE basic_mpi
USE model
USE FMZM
USE coeff
USE basis_transformation
USE speed_functions
use load_prec_coeffs
IMPLICIT NONE
PUBLIC
contains
!===============================================================================
SUBROUTINE coulombSelfGK_polarization(l,k,rCaa)
! add polarization terms to the self GK Coulomb collision operator
!
! Note: - adiabatic electron are assumed.
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
! Energy and FLR indices
TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
!
! loc. vars
INTEGER :: p,m,f,j,n,icol
TYPE(FM) , SAVE :: kern_,NX0,NX1,bapowm,SPJ,pipjm_,Tljpmf_,bardmnkf_,XOUT
CHARACTER(len=T5len) :: Tljpmfstr_
!
CALL FM_ENTER_USER_ROUTINE
!
!
XOUT = TO_FM('0')
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
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))
jloop: DO j=0, JEmaxx
SPJ = sigmapj(p,j)
mloop: DO m=0,p
IF( m .eq. 0) THEN
bapowm = TO_FM(1)
ELSE
bapowm = 2*TO_FM(bafm_/2)**(m)
ENDIF
!! pipjm_ = Pipjm(p,j,m)
pipjm_ = TO_FM(PipjmSTR(p,j,m))
IF( pipjm_ .ne. TO_FM(0)) THEN
nloop: DO n=0,naFLR
kern_ = kerneln(bafm_,n)
floop: DO f=0,n+k
bardmnkf_ = ALL2L(n,k,f,m)
IF( p .le. l + m + 2*f) THEN
Tljpmfstr_ = GET_Tljpmf(l,j,p,m,f)
Tljpmf_ = TO_FM(Tljpmfstr_)
IF(Tljpmf_ .ne. TO_FM('0')) THEN
XOUT = XOUT + NX0*NX1*Tljpmf_*pipjm_*bardmnkf_*bapowm*kern_/SPJ*FACTORIAL(TO_FM(n))/FACTORIAL(TO_FM(n+m))
ENDIF
ENDIF
ENDDO floop
ENDDO nloop
ENDIF
ENDDO mloop
ENDDO jloop
ENDDO ploop
!
! Coupling due to electrostatic potential
!
npoissonloop: DO n=0,Jmaxa
icol = numba(0,n)
kern_ = kerneln(bafm_,n)
rCaa(icol) = rCaa(icol) + XOUT*kern_/TO_FM(tauie)/poisson_coeff
ENDDO npoissonloop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE coulombSelfGK_polarization
!===============================================================================
FUNCTION Pipjm(p,j,m) RESULT(XOUT)
! computes the polarization FLR terms contained in the spherical harmonics particle moments, \Pi^{pjm}(b_a)
IMPLICIT NONE
!
INTEGER, INTENT(in) :: p,j,m
TYPE(FM) :: XOUT
!
! loc. vars.
INTEGER :: h,s,sp
TYPE(FM), SAVE :: T5_,kerns_,kernsp_,bapow_,dshspm_
CHARACTER(len=T5len) :: T5str_
!
CALL FM_ENTER_USER_FUNCTION(XOUT)
!
XOUT = TO_FM('0.0')
!
bapow_ = TO_FM(bafm_/2)**m
!
hloop: DO h=0,j + floor(real(p+m)/2)
T5STR_ = GET_T5(p,j,m,0,h)
T5_ = TO_FM(T5STR_)
IF(T5_ .ne. TO_FM('0')) THEN
sloop: DO s=0,naFLR
kerns_ = kerneln(bafm_,s)
sploop: DO sp=0, s + m + h
dshspm_ = ALLX2L(s,h,sp,m)
kernsp_ = kerneln(bafm_,sp)
XOUT = XOUT + T5_/FACTORIAL(TO_FM( s+ m))*FACTORIAL(TO_FM(s))*bapow_*kerns_*kernsp_*dshspm_
ENDDO sploop
ENDDO sloop
ENDIF
ENDDO hloop
!
CALL FM_EXIT_USER_FUNCTION(XOUT)
!
END FUNCTION Pipjm
!
!===============================================================================
!
SUBROUTINE evaluate_poisson
! evaluates the poisson coefficient, ion polarization and electron adiabatique response
!
implicit none
!
! loc. vars.
integer :: n
TYPE(FM), save :: kern_,XOUT
!
CALL FM_ENTER_USER_ROUTINE
!
poisson_coeff = TO_FM(1) ! electron adiabatic response
!
XOUT = TO_FM('0')
!
! ion polarization
DO n=0,2*naFLR
kern_ = kerneln(bafm_,n)
XOUT = XOUT + kern_*kern_
ENDDO
!
poisson_coeff = poisson_coeff + (TO_FM(1) - XOUT)/TO_FM(tauei)
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE evaluate_poisson
!
!===============================================================================
!
end module coulomb_pol
c4science · Help