Homec4science
Diffusion COSOlver (new_coulomb) Browse

cosolver/src/pitchangle_mod.f90
83b8f38f6c51new_coulomb

Last Change
Raw File

pitchangle_mod.f90
View Options

MODULE pitchangle
!
! Linearized pitch-angle scattering (wo velocity-dependent coll. freq.) and deflection (w velocity-dependent coll. freq.)
!
USE prec_const
USE basic
USE model
USE FMZM
USE coeff
USE basis_transformation
USE speed_functions
IMPLICIT NONE
PUBLIC
CONTAINS
SUBROUTINE LorentzDK(l,k,rCabT)
! Lorentz Pitch-angle scattering Operator (with constant coll frequency)
!
! Note:
! * Remove the element p=0 from the summation since Injlkmp is not defined for p=0
USE prec_const
IMPLICIT NONE
!
!
INTEGER, INTENT(in) :: l,k
! REAL(dp), DIMENSION(:,:),INTENT(inout) :: CabT
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
!
INTEGER :: p,j,g,h,t,t1,j1
TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T4,SPJ,Lptt1,Lpjj1,T4inv
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_,T4invstr_
CALL FM_ENTER_USER_ROUTINE
!
irow = numba(l,k)
!
! Pre compute quantities
!
! Hermite normalization coeff
NX0 = 1/SQRT(FACTORIAL(TO_FM(l))*(TO_FM(2)**l))
!
jloop:DO j=0,JEmaxx
ploop:DO p=1,l+2*k
NX1 = 2*p*(p+1)*TO_FM(2)**p*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p+1))/SQRT(PIFM)
SPJ = sigmapj(p,j)
tloop: DO t=0, k + l/2
!______________________________________________
! INVERSE BASIS TRANSFORMATION
T4str_ = GET_T4(p,t,l,k) ! ...
T4inv = 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)*TO_FM(T4str_)
!______________________________________________
IF (T4INV .ne. TO_FM('0.0')) THEN
gloop:DO g=0, min(Pmaxa,p+2*j)
NX2 = SQRT(FACTORIAL(TO_FM(g))*TO_FM(2)**g)
hloop:DO h=0, min(Jmaxa,j +p/2)
!_______________________________________
! BASIS TRANSFORMATION
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
IF (T4 .ne. TO_FM('0.0')) THEN
!_______________________________________
t1loop: DO t1=0, t
Lptt1=Lpjl(real(p,dp),real(t,dp),real(t1,dp))
j1loop:DO j1=0, j
Lpjj1=Lpjl(real(p,dp),real(j,dp),real(j1,dp))
NX3 = FACTORIAL(TO_FM( p + t1 + j1 -1))
! CONSTRUCT Lorentz COLLISIONAL MATRIX
icol = numba(g,h)
rCabT(icol ) = rCabT(icol ) - NX1*NX0*NX2*NX3*T4*T4inv*Lptt1*Lpjj1/SPJ
ENDDO j1loop
ENDDO t1loop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE LorentzDK
!
!===============================================================================
!
SUBROUTINE LorentzGK(l,k,rCabT)
! Electron GYROKINETIC Lorentz Pitch-angle scattering Operator
!
! Note:
! * Remove the element p=0 from the summation since Injlkmp is not defined for p=0
!
! * Recommanded: use instead deflection operator GK
USE prec_const
IMPLICIT NONE
!
!
INTEGER, INTENT(in) :: l,k
! REAL(dp), DIMENSION(:,:),INTENT(inout) :: CabT
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
!
INTEGER :: p,j,n,g,h,s,np,f,t1,j1,t,s1
TYPE(FM),SAVE :: NX0,NX1,NX2,NX3,T4,d0nhs_,SPJ,curlyA_,bdmnkf,T4inv,Lptt1,Lpjj1,dnphs1,ker_
INTEGER :: icol,irow
CHARACTER(len=T4len) :: T4str_
CALL FM_ENTER_USER_ROUTINE
!
!
! Pre compute quantities
!
! Hermite normalization coeff
NX0 = 1/SQRT(FACTORIAL(TO_FM(l))*(TO_FM(2)**l))
!
jloop:DO j=0,JEmaxx
nloop: DO n=0,naFLR
curlyA_ = curlyAmn(0,n,bafm_)
floop: DO f=0,n +k
bdmnkf= ALL2L(n,k,f,0) ! ... overline{d}_{nkf}^0
ploop:DO p=1,l + 2*f
NX1 = 2*p*(p+1)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(p))/FACTORIAL(TO_FM(2*p))/TO_FM((2*p+1))/SQRT(PIFM)
SPJ = sigmapj(p,j)
tloop: DO t=0, f + l/2
!______________________________________________
! INVERSE BASIS TRANSFORMATION
T4str_ = GET_T4(p,t,l,f) ! ...
T4inv = 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)*TO_FM(T4str_)
!______________________________________________
IF (T4INV .ne. TO_FM('0.0')) THEN
gloop:DO g=0, min(Pmaxa,p+2*j)
NX2 = SQRT(FACTORIAL(TO_FM(g))*TO_FM(2)**g)
hloop:DO h=0, j +p/2
!_______________________________________
! BASIS TRANSFORMATION
T4str_ = GET_T4(p,j,g,h) !T^gh_pj
T4 = TO_FM(T4str_)
IF (T4 .ne. TO_FM('0.0')) THEN
!_______________________________________
nploop: DO np=0, naFLR
ker_ = kerneln(bafm_,np)
s1loop: DO s1=0,min(Jmaxa,np+h)
dnphs1 = ALLX2L(np,h,s1,0)
t1loop: DO t1=0, t
Lptt1=Lpjl(real(p,dp),real(t,dp),real(t1,dp))
j1loop:DO j1=0, j
Lpjj1=Lpjl(real(p,dp),real(j,dp),real(j1,dp))
NX3 = FACTORIAL(TO_FM( p + t1 + j1 -1))
! CONSTRUCT Lorentz COLLISIONAL MATRIX
icol = numba(g,s1)
rCabT(icol ) = rCabT(icol ) - NX1*NX0*NX2*NX3*T4*T4inv*ker_*Lptt1*curlyA_*bdmnkf*dnphs1*Lpjj1/SPJ
ENDDO j1loop
ENDDO t1loop
ENDDO s1loop
ENDDO nploop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO floop
ENDDO nloop
ENDDO jloop
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE LorentzGK
!============================================================================
SUBROUTINE Cei1FDK(l,k,rCabF)
! Drift-kinetic Leading order sqrt(m_e/m_i) momentum convservative field part of Lorentz pitch-angle operator
!
!
USE prec_const
IMPLICIT NONE
!
!
INTEGER,INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabF
!
TYPE(FM),SAVE :: NX0
INTEGER ::irow
CALL FM_ENTER_USER_ROUTINE
!
!
NX0 = curlyAlk0pareiDK(l,k)
!
rCabF(numbb(1,0)) = rCabF(numbb(1,0)) +8*NX0/3/SQRT(TO_FM(2))/chifm_/PIFM
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Cei1FDK
!=================================================================================
SUBROUTINE Cei1FGK(l,k,rCabF)
! Gyrokinetic Leading order sqrt(m_e/m_i) momentum convservative field part of Lorentz pitch-angle operator
!
!
USE prec_const
IMPLICIT NONE
!
!
INTEGER,INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabF
!
TYPE(FM),SAVE :: NX0,NX1,ker_
INTEGER ::irow,n
!
CALL FM_ENTER_USER_ROUTINE
!
irow = numba(l,k)
!
NX0 = curlyAlk0pareiGK(l,k)
NX1 = curlyAlk0perpeiGK(l,k)
nloop: DO n=0,niFLR
ker_ = kerneln(bbfm_,n)
rCabF(numbb(1,n)) = rCabF(numbb(1,n)) +ker_*8*NX0/3/SQRT(TO_FM(2))/chifm_/PIFM
rCabF(numbb(0,n)) = rCabF(numbb(0,n)) - ker_*4/SQRT(PIFM)/chifm_*bbfm_*NX1
rCabF(numbb(0,n+1)) = rCabF(numbb(0,n+1)) + ker_*4/SQRT(PIFM)/chifm_*bbfm_*NX1
!
!
ENDDO nloop
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE Cei1FGK
!================================================================================
SUBROUTINE pitchTDK(l,k,rCabT)
!
! Construct the DK pitch-angle scattering operator with velocity dependent deflection frequency \nu_{ab}^D(v)
!
!
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,nusD
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,JEmaxx
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))
nusD = nusabDpjd(p,j,d)
rCabT(icol) = rCabT(icol) - TO_FM(p)*(TO_FM(p) +1)*NX0*NX1*T4inv*NX2*T4*Lphd*nusD/SPJ
ENDDO dloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDIF
ENDDO tloop
ENDDO ploop
ENDDO jloop
!
!
! rCabT(:)= rCabT(:)/TO_FM(2)
!
CALL FM_EXIT_USER_ROUTINE
!
!
END SUBROUTINE pitchTDK
!
!================================================================================
!
SUBROUTINE pitchTGK(l,k,rCabT)
!
! Construct the GK pitch-angle scattering operator with velocity dependent deflection frequency \nu_{ab}^D(v)
!
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCabT
CALL FM_ENTER_USER_ROUTINE
! add pitch-angle (nu = nuD(v))
CALL pitchTDK(l,k,rCabT)
! FLR terms
CALL RDlk(l,k,rCabT)
!
!! rCaa(:) = rCaa(:)*TO_FM(2)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE pitchTGK
!
!================================================================================
!
SUBROUTINE pitchselfDK(l,k,rCaa)
! Add momentum conserving term to the deflection (DM) operator for like-species only.
! Refs: A.A. Schekochihin et al. (2009), see Eq. (B18)
implicit none
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCaa
!
! local vars.
INTEGER :: p,j,h,h1,icol
TYPE(FM), save:: NX0,T4inv,L1hh1,nuDh1p2,NXlk
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
! test component
!
IF(addtest) THEN
! add pitch-angle (nu = nuD(v))
CALL pitchTDK(l,k,rCaa)
!
!
ENDIF
!
! field component
!
IF(addback) call Rpara_pitch_lkDK(l,k,rCaa)
!
!
! rCaa(:) = rCaa(:)*TO_FM(2)
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE pitchselfDK
!
!================================================================================
!
SUBROUTINE pitchselfGK(l,k,rCaa)
! Add momentum conserving term to the deflection (DM) operator for like-species only.
! Refs: A. A. Schekochihin et al. (2009), see Eq. (B18)
implicit none
!
INTEGER, INTENT(IN) :: l,k
TYPE(FM), DIMENSION(:),INTENT(inout) :: rCaa
!
! local vars.
INTEGER :: p,j,h,h1,icol
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
! test component
IF(addtest) THEN
!
CALL pitchTDK(l,k,rCaa)
!
! Add FLR correction
IF( .not. DKTEST) THEN
CALL RDlk(l,k,rCaa)
ENDIF
!
ENDIF
!
! Momentum conserving term: v_\perp J1 U_\perp, v_|| J0 U_||
!
! field component
IF(addback) THEN
IF(DKBACK) THEN
call Rpara_pitch_lkDK(l,k,rCaa)
ELSE
CALL Rperp_pitch_lkGK(l,k,rCaa)
CALL Rpara_pitch_lkGK(l,k,rCaa)
ENDIF
ENDIF
!
! rCaa(:) = rCaa(:)*TO_FM(2) ! this factor comes from the spherical part of the veloctiy-space Laplacien defined withouth the 1/2 in the convention used here !
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE pitchselfGK
!
!================================================================================
!
subroutine evaluate_rUpara_pitch_DK
! evalute Upara. Only Done once
implicit none
!
! locs. vars.
INTEGER :: p,j,h,h1,icol
TYPE(FM), save:: NX0,T4inv,L1hh1,nu2,nuDh1p2
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
! Allocate array
IF(allocated(rUpara)) deallocate( rUpara)
!
CALL allocate_array(rUpara,numbb(0,0),numbb(Pmaxb,Jmaxb))
!
! contrust Uperp
nu2 = nuabDk(2)
!
ploop: DO p=0,Pmaxa
NX0 = TO_FM(1)/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
jloop: DO j=0,Jmaxa
icol = numba(p,j)
IF( p .ge. 1 .or. j .ge. 1) THEN
hloop: DO h=0,j + floor(real(p)/2._dp)
T4str_ = GET_T4(1,h,p,j) ! ...
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.0')) THEN
h1loop: DO h1= 0,h
L1hh1 = Lpjl(real(1,dp),real(h,dp),real(h1,dp))
nuDh1p2 = nuabDk(h1 +2)
rUpara(icol) = rUpara(icol) + T4inv*L1hh1*NX0*nuDh1p2
enddo h1loop
ENDIF
enddo hloop
ENDIF
enddo jloop
enddo ploop
!
rUpara(:) = rUpara(:)/TO_FM(2)/nu2
!
CALL FM_EXIT_USER_ROUTINE
!
end subroutine evaluate_rUpara_pitch_DK
!
!================================================================================
!
subroutine Rpara_pitch_lkDK(l,k,rCaa)
!
!
implicit none
INTEGER,INTENT(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCaa
!
! loc. vars.
INTEGER :: n,s,h,h1
TYPE(FM), save:: kern_, NX0,d0nks,T4inv,L1hh1,nuDh1p2,NXlk
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
NXlk = TO_FM('0')
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
IF( l .ge. 1 .or. k .ge. 1) THEN
hloop: DO h=0,k + floor(real(l)/2._dp)
T4str_ = GET_T4(1,h,l,k) ! ...
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*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.0')) THEN
h1loop: DO h1= 0,h
L1hh1 = Lpjl(real(1,dp),real(h,dp),real(h1,dp))
nuDh1p2 = nuabDk(h1 +2)
NXlk = NXlk + T4inv*L1hh1*NX0*nuDh1p2*TO_FM(8)/SQRT(PIFM)/TO_FM(3)
enddo h1loop
ENDIF
enddo hloop
ENDIF
!
rCaa(:) = rCaa(:) + rUpara(:)*NXlk* TO_FM(2)
!
CALL FM_EXIT_USER_ROUTINE
!
end subroutine Rpara_pitch_lkDK
!
!================================================================================
!
SUBROUTINE RDlk(l,k,rCaa)
!
! Computes the FLR corrections to deflection operator
!
IMPLICIT NONE
!
INTEGER,INTENT(in) :: l,k
TYPE(FM),DIMENSION(:),INTENT(inout):: rCaa
!
! 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, icolend
INTEGER :: pj(2)
!
!
CALL FM_ENTER_USER_ROUTINE
!
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
NX = - TO_FM(1)/SQRT(PIFM)*bafm_*bafm_
!
alpha0_ = TO_FM(0)
!
icolend = numba(Pmaxa,Jmaxa)
colloop: DO icol=1,icolend
pj = invnumba(icol)
p= pj(1)
j =pj(2)
NX1 = TO_FM(1)/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
gloop: DO g=0, l + 2*k
hloop: DO h=0, k + floor(real(l)/2._dp)
!______________________________________________
! 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
IF( (r .eq. g) .or. (r .eq. g+2) .or. (r .eq. g-2)) THEN
sloop: DO s=0, j + floor(real(p)/2._dp)
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_ = (TO_FM(1)/TO_FM(2*g+1) )*nuabDk( g+2 +s1 +h1 )
ENDIF
alpha0_ = alpha0_ + X2prodLegendre(r,g)/TO_FM(2)*nuabDk( (r +g)/2 +2 +s1 +h1 )
!
rCaa(icol) = rCaa(icol) + NX*NX0*NX1*T4inv1*T4inv2*Lrss1*Lghh1*alpha0_ * TO_FM(2) ! not the extra factor 2
ENDDO h1loop
ENDDO s1loop
ENDIF
ENDDO sloop
ENDIF
ENDDO rloop
ENDIF
ENDDO hloop
ENDDO gloop
ENDDO colloop
!
!
!
CALL FM_EXIT_USER_ROUTINE
!
END SUBROUTINE RDlk
!
!================================================================================
!
subroutine evaluate_rUperp_pitch_GK
!
! evaluate Uperp. Done only once
implicit none
!
! loc. vars
INTEGER :: p,j,n,s,h,h1,icol
TYPE(FM), save:: kern_, NX0,d1njs,T4inv,L0hh1,nuD2,nuDh1p1
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
! ALLOCATE ARRAY
IF( allocated( rUperp)) deallocate(rUperp)
!
call allocate_array(rUperp,numbb(0,0),numbb(Pmaxb,Jmaxb))
!
! contrust Uperp
nuD2 = nuabDk(2)
!
ploop: DO p=0,Pmaxa
NX0 = TO_FM(1)/SQRT(TO_FM(2)**p*FACTORIAL(TO_FM(p)))
jloop: DO j=0,Jmaxa
icol = numba(p,j)
nloop: DO n=0,naFLR - 1
kern_ =kerneln(bafm_,n)
sloop: DO s=0,n + j + 1
d1njs = ALLX2L(n,j,s,1)
hloop: DO h=0,s + floor(real(p)/2._dp)
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.0')) THEN
h1loop: DO h1= 0,h
L0hh1 = Lpjl(real(0,dp),real(h,dp),real(h1,dp))
nuDh1p1 = nuabDk(h1 +1)
rUperp(icol) = rUperp(icol) + T4inv*d1njs*L0hh1*bafm_&
*kern_*nuDh1p1/TO_FM(2)/TO_FM(n+1)*NX0
enddo h1loop
ENDIF
enddo hloop
enddo sloop
enddo nloop
enddo jloop
enddo ploop
!
rUperp(:) = rUperp(:)/TO_FM(2)/nuD2*TO_FM(3)
!
CALL FM_EXIT_USER_ROUTINE
!
end subroutine evaluate_rUperp_pitch_GK
!
!================================================================================
!
subroutine evaluate_rUpara_pitch_GK
! evalute Upara. Only Done once
implicit none
!
! locs. vars.
INTEGER :: p,j,n,s,h,h1,icol
TYPE(FM), save:: kern_, NX0,d0njs,T4inv,L1hh1,nuD2,nuDh1p2
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
! Allocate array
IF(allocated( rUpara)) deallocate(rUpara)
!
CALL allocate_array(rUpara,numbb(0,0),numbb(Pmaxb,Jmaxb))
!
! contrust Uperp
nuD2 = nuabDk(2)
!
ploop: DO p=0,Pmaxa
NX0 = TO_FM(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
d0njs = ALLX2L(n,j,s,0)
IF( p .ge. 1 .or. s .ge. 1) THEN
hloop: DO h=0,s + floor(real(p)/2._dp)
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.0')) THEN
h1loop: DO h1= 0,h
L1hh1 = Lpjl(real(1,dp),real(h,dp),real(h1,dp))
nuDh1p2 = nuabDk(h1 +2)
rUpara(icol) = rUpara(icol) + T4inv*d0njs*L1hh1*kern_*NX0*nuDh1p2
enddo h1loop
ENDIF
enddo hloop
ENDIF
enddo sloop
enddo nloop
enddo jloop
enddo ploop
!
rUpara(:) = rUpara(:)/TO_FM(2)/nuD2
!
CALL FM_EXIT_USER_ROUTINE
!
end subroutine evaluate_rUpara_pitch_GK
!
!================================================================================
!
subroutine Rperp_pitch_lkGK(l,k,rCabF)
!
!
implicit none
INTEGER,INTENT(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabF
! loc. vars
INTEGER :: n,s,h,h1
TYPE(FM), save:: kern_, NX0,d1nks,T4inv,L0hh1,nuDh1p1,NXlk
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
NXlk = TO_FM('0.0')
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
!
nloop: DO n=0,naFLR-1
kern_ =kerneln(bafm_,n)
sloop: DO s=0,n + k + 1
d1nks = ALLX2L(n,k,s,1)
hloop: DO h=0,s + floor(real(l)/2._dp)
T4str_ = GET_T4(0,h,l,s) !
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*TO_FM( 1)/TO_FM(2)/FACTORIAL(TO_FM(2*h +1)/2)*TO_FM(T4str_)
IF( T4inv .ne. TO_FM('0.0')) THEN
h1loop: DO h1= 0,h
L0hh1 = Lpjl(real(0,dp),real(h,dp),real(h1,dp))
nuDh1p1 = nuabDk(h1 +1)
NXlk =NXlk + T4inv*d1nks*L0hh1*bafm_*kern_*nuDh1p1/TO_FM(n+1)*NX0*TO_FM(4)/SQRT(PIFM)
enddo h1loop
ENDIF
enddo hloop
enddo sloop
enddo nloop
!
rCabF(:) = rCabF(:) + rUperp(:)*NXlk *TO_FM(2)
!
CALL FM_EXIT_USER_ROUTINE
!
end subroutine Rperp_pitch_lkGK
!
!================================================================================
!
subroutine Rpara_pitch_lkGK(l,k,rCabF)
!
!
implicit none
INTEGER,INTENT(in) :: l,k
TYPE(FM), DIMENSION(:),INTENT(INOUT):: rCabF
!
! loc. vars.
INTEGER :: n,s,h,h1
TYPE(FM), save:: kern_, NX0,d0nks,T4inv,L1hh1,nuDh1p2,NXlk
CHARACTER(len=T4len) :: T4str_
!
CALL FM_ENTER_USER_ROUTINE
!
NXlk = TO_FM('0')
!
NX0 = TO_FM(1)/SQRT(TO_FM(2)**l*FACTORIAL(TO_FM(l)))
nloop: DO n=0,naFLR
kern_ =kerneln(bafm_,n)
sloop: DO s=0,n + k
d0nks = ALLX2L(n,k,s,0)
IF( l .ge. 1 .or. s .ge. 1) THEN
hloop: DO h=0,s + floor(real(l)/2._dp)
T4str_ = GET_T4(1,h,l,s) ! ...
T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*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.0')) THEN
h1loop: DO h1= 0,h
L1hh1 = Lpjl(real(1,dp),real(h,dp),real(h1,dp))
nuDh1p2 = nuabDk(h1 +2)
NXlk = NXlk + T4inv*d0nks*L1hh1*kern_*NX0*nuDh1p2*TO_FM(8)/SQRT(PIFM)/TO_FM(3)
enddo h1loop
ENDIF
enddo hloop
ENDIF
enddo sloop
enddo nloop
!
rCabF(:) = rCabF(:) + rUpara(:)*NXlk* TO_FM(2)
!
CALL FM_EXIT_USER_ROUTINE
!
end subroutine Rpara_pitch_lkGK
!
!================================================================================
!
END MODULE pitchangle
c4science · Help