diff --git a/src/model_mod.f90 b/src/model_mod.f90 index aed7a22..0509413 100644 --- a/src/model_mod.f90 +++ b/src/model_mod.f90 @@ -1,344 +1,345 @@ module model ! phyical parameters use prec_const use FMZM use basic IMPLICIT NONE PRIVATE ! physical parameters real(dp), public, protected :: nu = 1._dp ! electron-ion collision freq. real(dp), public, protected :: mu =1._dp ! sqr(m_e/m_i) real(dp), public, protected :: tau = 1._dp ! T_i0/T_e0 real(dp), public, protected :: kperp = 1._dp ! perpendicular wavenumber normalized to cs ! electron-ion mass/temperature ratio real(dp), public, protected :: sigmaei, sigmaaa,sigmaie real(dp), public, protected :: tauei, tauie, tauaa real(dp), public, protected :: sigma_ = 1._dp ! real(dp), public, protected :: tau_ = 1._dp ! real(dp), public, protected :: chi_ = 1._dp ! real(dp), public, protected :: be_ = 0._dp ! perpendicular wavenumber normalized to electron thermal Larmor radius real(dp), public, protected :: bi_ = 0._dp ! perpendicular wavenumber normalized to electron thermal Larmor radius real(dp), public, protected :: ba_ = 0._dp real(dp), public, protected :: bb_ = 0._dp ! Multiple-Precision FM type TYPE(FM), public, protected :: sigmafm_ ! mass ratio TYPE(FM), public, protected :: taufm_ ! temperature ratio TYPE(FM), public, protected :: chifm_ ! SQRT(tau/sigma) = vTa/vTb TYPE(FM), public, protected :: thetaabfm_ ! Sugama theta_{ab} coefficient [see Eq. (29) in Sugama et al. (2009)] TYPE(FM), public, protected :: thetabafm_ ! Sugma theta_{ba} coefficient TYPE(FM), public, protected :: bafm_ TYPE(FM), public, protected :: bbfm_ ! ! Poisson TYPE(FM), public :: poisson_coeff ! ion polarization and adiabatic electron coeff. ! label *.h5 file character(len=15), public :: strnu character(len=15), public:: strmu character(len=15), public :: strkperp character(len=15), public :: strtau character(len=1), PUBLIC :: strETEST, strITEST,strEBACK,strIBACK,strESELF,strISELF,strGKE,strGKI ! Test/Back-Reaction Model Collision Operator: 0 => Coulomb INTEGER,PUBLIC :: ETEST=0 ! Electron Test Model INTEGER,PUBLIC :: ITEST=0 ! Ion Test Model INTEGER,PUBLIC :: EBACK=0 ! Electron Back-Reaction Model INTEGER,PUBLIC :: IBACK=0 ! Ion Back-Reaction Model INTEGER, PUBLIC :: ESELF=0 ! Electron self Model INTEGER, PUBLIC :: ISELF=0 ! Ion Self Model INTEGER, PUBLIC :: GKE =0 ! Gyrokinetic Electrons INTEGER, PUBLIC :: GKI =0 ! Gyrokinetic Ions ! Colliding species CHARACTER(len=2), PUBLIC, PROTECTED :: coll_ab ! Sugama Flag: used in AUXVAL for error functions necessary for back-reaction parts LOGICAL, PUBLIC, PROTECTED :: IFSUGAMA = .FALSE. ! Improved Sugama Flag: used in AUXVAL and MEMORY for precomputation of CurlI(para/perp)(a/b)pjk in GK mode LOGICAL, PUBLIC, PROTECTED :: IFIMPROVEDSUGAMA = .FALSE. ! Coulomb flag LOGICAL, PUBLIC, PROTECTED :: IFCOULOMB = .FALSE. ! Polarization flag LOGICAL, PUBLIC, PROTECTED :: IFPOL = .FALSE. ! Zeroth order Hirshman Sigmar Clarke flag logical, PUBLIC, protected :: IF0HSC = .false. ! GK Abel flag: used in AUXVAL for error functions necessary for back-reaction parts logical, PUBLIC, protected :: IFABEL = .false. ! DK/GK test or/and field: used only if GKI and/or GKE ! Note: implented only for Abel collision operator logical, PUBLIC, protected :: DKBACK = .false. ! use DK back reaction logical, PUBLIC, protected :: DKTEST = .false. ! use DK test part ! Flag turn on/off Test/Field parts logical, public, protected :: ADDBACK =.true. ! IF 1 add field part, 0 otherwise logical, public, protected :: ADDTEST=.true. ! IF 1 add test part, 0 otherwise ! Flag turn on/off spatial diff in test part (implemetend of GK Abel) logical, public, protected :: only_gk_part = .false. ! Gyrokinetic flag LOGICAL, PUBLIC, PROTECTED :: IFGK = .FALSE. public :: model_readinputs, set_model,set_model_fm,operatormodel_readinput,invert_model contains !================================================================================ subroutine model_readinputs ! Read inputs physical input parameters use basic, only: lu_in,ierr use prec_const implicit none character(len=8) :: fmt = '(f6.4)' NAMELIST /MODEL_PAR/ nu,mu,tau,kperp IF( me .ne. 0) GOTO 111 read(lu_in,MODEL_PAR) write(*,MODEL_PAR) 111 CONTINUE ! We compute the normalized collisional matrix nu = 1._dp ! ! Broadcast from PE 0 to ! CALL MPI_bcast(nu,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(mu,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(tau,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(kperp,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) ! convert to string write(strnu,fmt) nu write(strmu,fmt) mu write(strtau,fmt) tau write(strkperp,fmt) kperp ! set electron-ion mass/temperature ratio sigmaei = mu ! defined as me/mi !!!!!! sigmaie = 1/sigmaei sigmaaa = 1._dp tauei = 1/tau tauie = tau tauaa = 1._dp ! Set electron/ion perpendicular wavenumber be_ = kperp*SQRT(2._dp)*SQRT(sigmaei) bi_ = kperp*SQRT(2._dp)*SQRT(tauei) end subroutine model_readinputs !================================================================================ subroutine model_outputinputs ! output physical input parameters to *.h5 file use basic, only: lu_in,me use prec_const implicit none ! code here end subroutine model_outputinputs !================================================================================ subroutine set_model(ab) ! set the species mass/temperature ratios implicit none character(len=2), intent(in) :: ab ! Ion-Ion collisions if(ab .eq. 'ii') then chi_ = 1._dp tau_ = 1._dp sigma_ = 1._dp ba_ = bi_ bb_ = bi_ coll_ab = 'ii' end if ! Electron- Electron collisions if(ab .eq. 'ee') then chi_ = 1._dp tau_ = 1._dp sigma_ = 1._dp ba_ = be_ bb_ = be_ coll_ab = 'ee' end if ! Electron-ion collisions if(ab .eq. 'ei') then chi_ = SQRT(tauei/sigmaei) tau_ = tauei sigma_ = sigmaei ba_ = be_ bb_ = bi_ coll_ab = 'ei' end if ! Ion-electron collisions if(ab .eq. 'ie') then chi_ = SQRT(tauie/sigmaie) tau_ = tauie sigma_ = sigmaie ba_ = bi_ bb_ = be_ coll_ab = 'ie' end if ! set model parameter to FM type CALL set_model_fm ! end subroutine set_model !================================================================================ subroutine set_model_fm ! set mass/temperature ratio to FM type implicit none ! CALL FM_ENTER_USER_ROUTINE ! Set temperature and mass ratios to FM type sigmafm_ = TO_FM(sigma_) taufm_ = TO_FM(tau_) chifm_ = SQRT(taufm_/sigmafm_) ! Set perpendicular waves numbers to FM type bafm_ = TO_FM(ba_) bbfm_ = TO_FM(bb_) ! ! ! Set Sugama parameter thetaabfm_ = SQRT(TO_FM( taufm_ + chifm_**2)/(TO_FM(1) + chifm_**2)) thetabafm_ = SQRT(TO_FM( TO_FM(1)/taufm_ + TO_FM(1)/chifm_**2)/(TO_FM(1) + TO_FM(1)/chifm_**2)) ! CALL FM_EXIT_USER_ROUTINE end subroutine set_model_fm !================================================================================ subroutine invert_model ! inverse ab --> ba ! Used in the Sugama collision operator IMPLICIT NONE ! CALL FM_ENTER_USER_ROUTINE ! sigmafm_ = 1/sigmafm_ taufm_ = 1/taufm_ chifm_ = 1/chifm_ ! ! Set Sugama parameter thetaabfm_ = SQRT( (TO_FM(1) + sigmafm_)/(TO_FM(1) + sigmafm_/taufm_)) + thetabafm_ = SQRT( (TO_FM(1) + TO_FM(1)/sigmafm_)/(TO_FM(1) + taufm_/sigmafm_)) ! CALL FM_EXIT_USER_ROUTINE ! END SUBROUTINE invert_model !================================================================================ SUBROUTINE operatormodel_readinput ! Read Model Operator from input file USE BASIC IMPLICIT NONE ! Namelist NAMELIST /OPERATOR_MODEL/ ETEST,ITEST,EBACK,IBACK,ESELF,ISELF,GKE,GKI,DKTEST,DKBACK,ADDBACK,ADDTEST,only_gk_part ! IF(me .ne. 0) GOTO 111 READ(lu_in,OPERATOR_MODEL) WRITE(*,OPERATOR_MODEL) ! 111 CONTINUE ! ! MPI COMM CALL MPI_bcast(ETEST,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(EBACK,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(ITEST,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(IBACK,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(ESELF,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(ISELF,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(GKI,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(GKE,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(DKTEST,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(DKBACK,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(ADDBACK,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast(ADDTEST,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr) CALL MPI_bcast( only_gk_part,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr) ! ! write(strETEST,'(I1)') ETEST write(strITEST,'(I1)') ITEST write(strEBACK,'(I1)') EBACK write(strIBACK,'(I1)') IBACK write(strESELF,'(I1)') ESELF write(strISELF,'(I1)') ISELF write(strGKE,'(I1)') GKE write(strGKI,'(I1)') GKI ! Set different flags ! If Sugama COOP IF(ETEST .eq. 3 .or. EBACK .eq. 3 .or. ITEST .eq. 3 .or. IBACK .eq. 3 .or. ESELF .eq. 3 .or. ISELF .eq. 3) THEN IFSUGAMA = .true. ENDIF ! If Improved Sugama COOP IF(ETEST .eq. 5 .or. EBACK .eq. 5 .or. ITEST .eq. 5 .or. IBACK .eq. 5 .or. ESELF .eq. 5 .or. ISELF .eq. 5) THEN IFSUGAMA = .true. IFIMPROVEDSUGAMA = .true. imaxx = JEmaxx ENDIF ! IF zeroth order Hirshman-Sigmar-Clarke (0HSC) IF(ETEST .eq. 6 .or. EBACK .eq. 6 .or. ITEST .eq. 6 .or. IBACK .eq. 6 .or. ESELF .eq. 6 .or. ISELF .eq. 6) THEN IF0HSC = .true. imaxx = JEmaxx ENDIF ! IF GK abel IF(ISELF .eq. 7 ) THEN IFABEL = .true. ENDIF ! IF Coulomb operator IF(ETEST .eq. 1 .or. EBACK .eq. 1 .or. ITEST .eq. 1 .or. IBACK .eq. 1 .or. ESELF .eq. 1 .or. ISELF .eq. 1) THEN IFCOULOMB = .true. ENDIF ! IF gyrokinetic IF(GKE .eq. 1 .or. GKI .eq. 1) THEN IFGK = .true. ENDIF ! ! IF GK Polarization IF(ISELF .eq. 9) THEN IFPOL = .true. ENDIF ! END SUBROUTINE operatormodel_readinput END MODULE model diff --git a/src/speed_functions_mod.f90 b/src/speed_functions_mod.f90 index 216f1e7..8515693 100644 --- a/src/speed_functions_mod.f90 +++ b/src/speed_functions_mod.f90 @@ -1,1781 +1,1781 @@ MODULE speed_functions ! ! contains the closed formulas used to compute the speed functions of the test and back-reaction parts of COOPs ! ! ! Author : B. J. Frei, SPC/Lausanne, 2019 ! USE prec_const USE model USE FMZM USE coeff USE basic use array use basis_transformation IMPLICIT NONE CONTAINS !To promote cleanliness, use .h include files include 'coulomb_speed_functions.h' !----------------------------------------------------------- !> @brief !> Computes the analytical closed form of even power velocity integral of error function derivative ! !> !> \f[ e_{ab}^{k}= \int_0^\infty d s_a s_a^{2k} erf'(s_b) e^{-s_a^2} =\frac{(k-1/2)!}{(-1/2)!(1+\chi^2)^{k+1/2}}. \f] !> !> Note: this integral converges only if k > -1 !> !> @param[in] k even velocity power !> !> @param[out] XOUT !----------------------------------------------------------- FUNCTION ekab(k) RESULT(XOUT) ! ! compute int_0^\infty d x Erf'(s_a) s_b^{2k} e^{-x^2) [see Eq.(60) in Ji & Held (2006)] ! ! Returns 0 if k < 0 ! IMPLICIT NONE INTEGER, INTENT(in) :: k TYPE(FM) :: XOUT CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! IF( k .ge. 0) THEN XOUT = FACTORIAL(TO_FM(2*k - 1 )/2)/FACTORIAL(TO_FM(-1)/2)/(TO_FM(1._dp + chifm_**2))**(real(k,dp) + 0.5_dp) ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION ekab !======================================================================! !----------------------------------------------------------- !> @brief !> Computes the analytical closed form of odd power velocity power integral of error function ! !> !> \f[ E_{ab}^{k}= \int_0^\infty d s_a s_a^{2k+1} erf(s_b)e^{-s_a^2}= \frac{\chi}{2} \sum_{j=0}^k \frac{k!}{j!}e_{ab}^j \f] ! ! Note: Only for k > -2, otherwise integral does not converge ! !> @param[in] k even velocity power !> !> @param[out] XOUT !----------------------------------------------------------- FUNCTION Erfkab(k) RESULT(XOUT) ! ! computes int_0^\infty d x Erf(s_a) s_b^{2k+1} e^{-x^2) [see Eq.(61) in Ji & Held (2006)] ! ! ! Note: Returns 0 if k < -1. ! IMPLICIT NONE ! INTEGER, INTENT(in) :: k TYPE(FM) :: XOUT TYPE(FM), SAVE:: ekab_ INTEGER :: i CALL FM_ENTER_USER_FUNCTION(XOUT) XOUT = TO_FM('0.0') ekab_ = TO_FM('0.0') ! ! ===== IF k > -2 ===== IF( k .ge. 0) THEN DO i=0, k ekab_ = ekab(i) XOUT = XOUT + FACTORIAL(TO_FM(k))/FACTORIAL(TO_FM(i))*ekab_ END DO ! XOUT = XOUT*chifm_/2 ! ENDIF ! IF( k .eq. -1 ) THEN XOUT = ASINH(chifm_) ENDIF ! ! ======== NOT DEFINED RETURN 0 ===== ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! end FUNCTION Erfkab !======================================================================! !======================================================================! ! ! Sugama Collision Operator ! !======================================================================! !======================================================================! FUNCTION nusSabparapjd(p,j,d) RESULT(XOUT) ! ! Compute the velocity integrals over the energy diffusion frequency, ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: p,j,d TYPE(FM) :: XOUT ! ! LOCAL variables INTEGER :: j1 TYPE(FM), SAVE :: Lpjj1_,alphaErfpj10,alphaErfpj11,alphaepj10,alphaepj11 TYPE(FM), SAVE :: Erfkab0_,Erfkab1_,ekab0_,ekab1_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') DO j1 =0, j Lpjj1_= Lpjl(real(p,dp),real(j,dp),real(j1,dp)) ! Error function integrals ! n =0 Erfkab0_ = Erfkab_array(p + d + j1 - 1) ekab0_ = ekab_array(p + d + j1 ) ! n=1 Erfkab1_ = Erfkab_array(p + d + j1 - 2) ekab1_ = ekab_array(p + d + j1 - 1) ! Compute coeff. ! n=0 alphaErfpj10 = - sigmafm_*TO_FM(2*(p + 2*j1))/taufm_ alphaepj10 = 2*TO_FM(2*j1 +p)*(TO_FM(1) + chifm_**2 )/chifm_ !n=1 alphaErfpj11 = - sigmafm_*TO_FM(p + 2*j1)*TO_FM(2 - p - 2*j1 )/taufm_ alphaepj11 = - TO_FM(p + 2*j1)*TO_FM( 2*j1 + p -2 )/chifm_ XOUT = XOUT + Lpjj1_ * ( alphaErfpj10*Erfkab0_ + ekab0_*alphaepj10 & + alphaErfpj11*Erfkab1_ + ekab1_*alphaepj11) ENDDO ! XOUT = 4*XOUT/SQRT(PIFM) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nusSabparapjd ! !======================================================================! FUNCTION nusabDpjd(p,j,d) RESULT(XOUT) ! ! Compute the velocity integrals over the deflection frequency ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: p,j,d TYPE(FM) :: XOUT ! ! local variable INTEGER :: j1 TYPE(FM), SAVE :: X0,X1,X2,Lpjj1_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! DO j1 =0,j Lpjj1_= Lpjl(real(p,dp),real(j,dp),real(j1,dp)) X0 = Erfkab_array( p + d + j1 -1) X1 = Erfkab_array(p + d + j1 -2 ) X2 = ekab_array( p + d + j1 -1) XOUT = XOUT + Lpjj1_*( X0 - X1/2/chifm_**2 + X2/2/chifm_ ) ENDDO ! XOUT = XOUT/SQRT(PIFM)*4 ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nusabDpjd ! !======================================================================! FUNCTION nusSabpjd(p,j,d) result(XOUT) ! Compute the velocity integrals over the deflection + energy diffusion frequency ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: p,j,d TYPE(FM) :: XOUT ! ! LOCAL VAR. TYPE(FM),SAVE :: X1,X2 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! X1 = nusabDpjd(p,j,d) X2 = nusSabparapjd(p,j,d) XOUT = X2- TO_FM(p)*TO_FM((p+1))*X1 ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nusSabpjd ! !======================================================================! FUNCTION IabSlk(l,k) RESULT(XOUT) ! ! Computes I_{ab}^{Slk} appearing in CP_{ab}[f]^lk ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! LOCAL VARS. INTEGER :: h,h1 TYPE(FM), SAVE :: L0hh1_, T4inv,Erfkab_,ekab_ CHARACTER(len=T4len) :: T4invstr_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! XOUT = TO_FM('0') ! DO h=0,k + floor(real(l/2)) T4invstr_ = GET_T4(0,h,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(h))*TO_FM( 1)/2/FACTORIAL(TO_FM(2*h +1)/2) DO h1=0,h L0hh1_= Lpjl(real(0,dp),real(h,dp),real(h1,dp)) Erfkab_ = Erfkab_array(h1) ekab_ = ekab_array(h1+1) XOUT = XOUT + L0hh1_*T4inv*( Erfkab_ - chifm_*(1 + chifm_**2)*ekab_) ENDDO ENDDO ! XOUT = - XOUT*8/SQRT(PIFM)/chifm_**2 ! ... check the coefficient in the front ! XOUT = - XOUT*6/chifm_**2 ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION IabSlk ! !======================================================================! FUNCTION IabSparalk(l,k) RESULT(XOUT) ! ! Computes I_{ab}^{S||lk} appearing in CP_{ab}[f]^lk and in PC_{ab}[f]^lk ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! LOCAL VARS. INTEGER :: h,h1 TYPE(FM), SAVE :: L1hh1_, T4inv,Erfkab_,ekab_ CHARACTER(len=T4len) :: T4invstr_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! XOUT = TO_FM('0') ! IF ( l .ge. 1 .or. k .ge. 1) THEN DO h=0,k + floor(real(l/2)) T4invstr_ = GET_T4(1,h,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(h))*TO_FM( 1 + 2)/2/FACTORIAL(TO_FM(2 +2*h +1)/2) DO h1=0,h L1hh1_= Lpjl(real(1,dp),real(h,dp),real(h1,dp)) Erfkab_ = Erfkab_array(h1) ekab_ = ekab_array(h1+1) XOUT = XOUT + L1hh1_*T4inv*( Erfkab_ - chifm_*ekab_)/chifm_**2 ENDDO ENDDO ! XOUT = XOUT*4/SQRT(PIFM)/3 ! ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION IabSparalk ! !======================================================================! ! FUNCTION curlyVab(l,k) RESULT(XOUT) ! ! Computes \mathcal{V}_{ab}^lk appearing in the back-reaction velocity integral of Sugama operator ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! ! LOCAL VARS. INTEGER ::h, h1 TYPE(FM), SAVE :: T4inv,Erfkab_,ekab_,L1hh1_ CHARACTER(len=T4len) :: T4invstr_ CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0') ! ! IF (l .ge. 1 .or. k .ge. 1)THEN ! hloop: DO h=0, k +floor(real(l)/2._dp) T4invstr_ = GET_T4(1,h,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(h))*TO_FM( 1 + 2)/2/FACTORIAL(TO_FM(2 +2*h +1)/2) IF(T4inv .ne. TO_FM('0')) THEN h1loop: DO h1=0,h L1hh1_ = Lpjl(real(1,dp),real(h,dp),real(h1,dp)) Erfkab_ = Erfkab_array(h1) ekab_ = ekab_array(h1+1) XOUT = XOUT + L1hh1_*T4inv*TO_FM(4)/TO_FM(3)/SQRT(PIFM)*( TO_FM(3)*SQRT(PIFM)/TO_FM(4)/chifm_**2*( Erfkab_ - chifm_*ekab_ ) & + chifm_*(thetaabfm_ -TO_FM(1))*FACTORIAL(TO_FM(3 + 2*h1)/2)/TO_FM(2)/(TO_FM(1) + chifm_**2)**(1.5_dp)) ENDDO h1loop ENDIF ENDDO hloop ! ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyVab ! !=======================================================================! ! FUNCTION curlyVba(l,k) RESULT(XOUT) ! ! Computes \mathcal{V}_{ba}^lk appearing in the back-reaction velocity integral of Sugama operator ! ! Note: - uses E_{ba}^k and e_{ba}^k arrays instead precomputed if IFSUGAMA ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! ! LOCAL VARS. INTEGER ::h, h1 TYPE(FM), SAVE :: T4inv,Erfkba_,ekba_,L1hh1_ CHARACTER(len=T4len) :: T4invstr_ CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! ! IF (l .ge. 1 .or. k .ge. 1)THEN ! hloop: DO h=0, k + floor(real(l)/2._dp) T4invstr_ = GET_T4(1,h,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(h))*TO_FM( 1 + 2)/2/FACTORIAL(TO_FM(2 +2*h +1)/2) IF( T4inv .ne. ZEROFM) THEN h1loop: DO h1=0,h L1hh1_ = Lpjl(real(1,dp),real(h,dp),real(h1,dp)) Erfkba_ = Erfkba_array(h1) ekba_ = ekba_array(h1+1) XOUT = XOUT + L1hh1_*T4inv*TO_FM(4)/TO_FM(3)/SQRT(PIFM)*( TO_FM(3)*SQRT(PIFM)/TO_FM(2)/TO_FM(2)*chifm_**(2)*( Erfkba_ - ekba_/chifm_ ) & + (thetabafm_ -TO_FM(1))*FACTORIAL(TO_FM(3 + 2*h1)/2)/chifm_/TO_FM(2)/(TO_FM(1) + TO_FM(1)/chifm_**(2))**(1.5_dp)) ENDDO h1loop ENDIF ENDDO hloop ! ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyVba ! !======================================================================! ! FUNCTION curlyWab(l,k) RESULT(XOUT) ! ! Computes \mathcal{W}_{ab}^lk appearing in the back-reaction velocity integral of Sugama operator ! ! IMPLICIT NONE ! ! INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! ! LOCAL VARS. INTEGER ::h, h1 TYPE(FM), SAVE :: T4inv,L0hh1_ CHARACTER(len=T4len) :: T4invstr_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! hloop: DO h=0, k + floor(real(l)/2._dp) T4invstr_ = GET_T4(0,h,l,k) ! ... (T^-1)^0h_lk = ... T^lk_0h T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(h))*TO_FM( 1)/2/FACTORIAL(TO_FM(2*h +1)/2) IF( T4inv .ne. ZEROFM) THEN h1loop: DO h1=0,h L0hh1_ = Lpjl(real(0,dp),real(h,dp),real(h1,dp)) XOUT = XOUT - thetaabfm_*L0hh1_*T4inv*TO_FM(6)/(chifm_**2)*( Erfkab_array(h1) -chifm_*(TO_FM(1) + chifm_**2)*ekab_array(h1 +1) ) ENDDO h1loop ENDIF ENDDO hloop ! IF(l .eq. 2 .and. k .eq. 0) THEN XOUT = XOUT - thetaabfm_*4*chifm_*(thetaabfm_ - TO_FM(1))/(TO_FM(1) + chifm_**2)**(1.5_dp) ENDIF ! IF(l .eq. 0 .and. k .eq. 1) THEN XOUT = XOUT + thetaabfm_*TO_FM(2)*chifm_*(thetaabfm_ - TO_FM(1))/(TO_FM(1) + chifm_**2)**(1.5_dp) ENDIF ! ! XOUT = TO_FM(4)*XOUT/TO_FM(3)/SQRT(PIFM) ! .... rescale sugama collisional time ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyWab ! !======================================================================! ! FUNCTION curlyWba(l,k) RESULT(XOUT) ! ! Computes \mathcal{W}_{ba}^lk appearing in the back-reaction velocity integral of Sugama operator ! ! Note: - uses E_{ba}^k and e_{ba}^k arrays instead precomputed if IFSUGAMA ! - need to multiply by a factor tau_{ab}/tau_{ba} ! IMPLICIT NONE ! ! INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! ! LOCAL VARS. INTEGER ::h, h1 TYPE(FM), SAVE :: T4inv,L0hh1_,NX0 CHARACTER(len=T4len) :: T4invstr_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! ! tau_{ab}/tau_{ba} ! NX0 = sigmafm_**2*chifm_**3 ! hloop: DO h=0, k + floor(real(l)/2._dp) T4invstr_ = GET_T4(0,h,l,k) ! ... (T^-1)^0h_lk = ... T^lk_0h T4inv =TO_FM(T4invstr_)*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) IF(T4inv .ne. ZEROFM) THEN h1loop: DO h1=0,h L0hh1_ = Lpjl(real(0,dp),real(h,dp),real(h1,dp)) XOUT = XOUT - thetabafm_*L0hh1_*T4inv*6*(chifm_**2)*( Erfkba_array(h1) - (TO_FM(1) + TO_FM(1)/(chifm_**2))*ekba_array(h1 +1)/chifm_ ) ENDDO h1loop ENDIF ENDDO hloop ! ! ! T4invstr_ = GET_T4(0,0,l,k) ! ... (T^-1)^0h_lk = ... T^lk_0h ! T4inv =TO_FM(T4invstr_)*SQRT(PIFM)*TO_FM(2)**l*FACTORIAL(TO_FM(l))**TO_FM( 1)/2/FACTORIAL(TO_FM(1)/2) ! ! ! XOUT = XOUT + (thetabafm_-1)*thetabafm_*3/PIFM/chifm_/(( 1 + 1/chifm_**(2))**(3/2))*T4inv ! ! IF(l .eq. 2 .and. k .eq. 0) THEN XOUT = XOUT - thetabafm_*4*(thetabafm_ - TO_FM(1))/(TO_FM(1) + 1/(chifm_**2))**(1.5_dp)/chifm_ ENDIF ! IF(l .eq. 0 .and. k .eq. 1) THEN XOUT = XOUT + thetabafm_*2*(thetabafm_ - TO_FM(1))/(TO_FM(1) + 1/(chifm_**2))**(1.5_dp)/chifm_ ENDIF ! XOUT = XOUT*NX0 ! ... rescale by tau_{ab}/tau_{ba} ! ! XOUT = 4*XOUT/3/SQRT(PIFM) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyWba ! !=======================================================================! FUNCTION nuabparak(k) RESULT(XOUT) ! ! computes \overline{nu}_{ab}^{\parallel k} = \int d s s^{2k} \nu_{ab}^{\parallel}(v) e^{-s^2}. ! ! Note : - the integral converges only for k >= 1 ! - k =1 value is only for GK and like-species collisions ! IMPLICIT NONE ! INTEGER, INTENT(in) :: k TYPE(FM) :: XOUT ! ! Local vars TYPE(FM),SAVE :: Erfkab_,ekab_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! XOUT = TO_FM(0) ! IF(k .eq. 1) THEN ! only for like species !! XOUT =TO_FM('0.532839975353552023569079399229905769541511547115312662423384') ENDIF ! IF(k .ge. 2) THEN Erfkab_ = Erfkab_array(k-3) ekab_ = ekab_array(k-2) ! XOUT = (Erfkab_ - chifm_*ekab_)/chifm_**2 ENDIF ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nuabparak ! !=======================================================================! ! FUNCTION nuabDk(k) RESULT(XOUT) ! ! computes \overline{nu}_{ab}^{D k} = \int d s s^{2k} \nu_{ab}^{D}(v) e^{-s^2} ! ! Note : maybe need to add the factor 4/sqrt(pi) ! IMPLICIT NONE ! INTEGER, INTENT(in) :: k TYPE(FM) :: XOUT ! ! Local vars TYPE(FM),SAVE :: Erfkabm2_,Erfkabm3_,ekabm2_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! XOUT = TO_FM(0) ! IF(k .eq. 1 ) THEN ! ... k=1 is not closed; used Mathematica instead. ! XOUT = TO_FM('0.614953599342767013448069625364839424257404554703979079541604') ! value for inter-species only !!!!!! ! ELSE ! ... for k > 1 ! Erfkabm2_ = Erfkab_array(k-2) Erfkabm3_ = Erfkab_array(k-3) ekabm2_ = ekab_array(k-2) ! XOUT = Erfkabm2_ - Erfkabm3_/TO_FM(2)/chifm_**2 + ekabm2_/TO_FM(2)/chifm_ ! ENDIF ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nuabDk ! !=======================================================================! ! ! FUNCTION HlLkCabTS(l,k) RESULT(XOUT) ! ! computes \int d v H_l Lk C_{ab}^{TS}[F_{aM } m_a v_|| / T_a] appearing in the ! ! IMPLICIT NONE ! INTEGER, INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! Local vars INTEGER :: h,h1 TYPE(FM),SAVE :: L0hh1_,T4inv CHARACTER(len=T4len) :: T4str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! XOUT = TO_FM(0) ! hloop: DO h=0,k + floor(real(l/2)) T4str_ = GET_T4(0,h,l,k) ! ... T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*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 h1loop: DO h1 =0,h L0hh1_ = Lpjl(real(0,dp),real(h,dp),real(h1,dp)) XOUT = XOUT + L0hh1_*3*SQRT(PIFM)/4/chifm_**2*( Erfkab_arraY(h1-1) - chifm_*ekab_array(h1) )*T4inv ENDDO h1loop ! XOUT = XOUT - T4inv*chifm_*(thetaabfm_ -1)/(1 + chifm_**2)**(TO_FM(3)/2)*FACTORIAL(TO_FM(2*h - 3)/2)/2/SQRT(PIFM)/FACTORIAL(TO_FM(h)) ENDIF ! ENDDO hloop ! XOUT = - 2*thetaabfm_*(1 + chifm_**2)*4/SQRT(PIFM)*XOUT ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION HlLkCabTS !___________________________________________________ ! FUNCTION curlyIabls(l,s) RESULT(XOUT) ! ! computes \mathcal{I}_{ab}^ls ! ! IMPLICIT NONE ! INTEGER, INTENT(in) :: l,s TYPE(FM) :: XOUT ! ! Local vars INTEGER :: h,h1 TYPE(FM),SAVE :: L0hh1_,T4inv CHARACTER(len=T4len) :: T4str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! XOUT = TO_FM('0') ! 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)/2/FACTORIAL(TO_FM(2*h +1)/2)*TO_FM(T4str_) IF(T4inv .ne. TO_FM('0')) THEN h1loop: DO h1 =0,h L0hh1_ = Lpjl(real(0,dp),real(h,dp),real(h1,dp)) XOUT = XOUT + T4inv*L0hh1_/chifm_/chifm_*(Erfkab_array(h1-1) - chifm_*ekab_array(h1) ) ENDDO h1loop ! ENDIF ! ENDDO hloop ! XOUT = - TO_FM(4)*THETAABFM_*(TO_FM(1) + CHIFM_*CHIFM_)/SQRT(PIFM)*XOUT ! IF(L .EQ. 0 .AND. S .EQ. 0) THEN XOUT = XOUT - TO_FM(4)*CHIFM_*THETAABFM_*(THETAABFM_-TO_FM(1))/SQRT((TO_FM(1) + CHIFM_*CHIFM_))/TO_FM(3)/SQRT(PIFM) ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyIabls !___________________________________________________ ! !___________________________________________________ ! FUNCTION curlyIbals(l,s) RESULT(XOUT) ! ! computes \mathcal{I}_{ba}^ls ! ! IMPLICIT NONE ! INTEGER, INTENT(in) :: l,s TYPE(FM) :: XOUT ! ! Local vars INTEGER :: h,h1 TYPE(FM),SAVE :: L0hh1_,T4inv CHARACTER(len=T4len) :: T4str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! XOUT = TO_FM('0.0') ! 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)/2/FACTORIAL(TO_FM(2*h +1)/2)*TO_FM(T4str_) IF(T4inv .ne. TO_FM(0)) THEN h1loop: DO h1 =0,h L0hh1_ = Lpjl(real(0,dp),real(h,dp),real(h1,dp)) XOUT = XOUT + T4inv*L0hh1_*chifm_*chifm_*(Erfkba_array(h1-1) - ekba_array(h1)/chifm_ ) ENDDO h1loop ! ENDIF ! ENDDO hloop ! XOUT = -TO_FM(4)*thetabafm_*(TO_FM(1) + TO_FM(1)/chifm_/chifm_)/SQRT(PIFM)*XOUT*sigmafm_*sigmafm_*chifm_*chifm_*chifm_ ! IF(l .eq. 0 .and. s .eq. 0) THEN XOUT = XOUT - TO_FM(4)/chifm_*thetabafm_*(thetabafm_- TO_FM(1))/SQRT((TO_FM(1) + TO_FM(1)/chifm_/chifm_))/TO_FM(3)/SQRT(PIFM)*sigmafm_**2*chifm_**3 ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyIbals !___________________________________________________ ! FUNCTION curlyUparalk(l,k) RESULT(XOUT) ! ! Computes \mathcal{U}_{\parallel}^{lk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! IF( l .eq. 1 ) THEN XOUT = kerneln(bafm_,k)/SQRT(TO_FM(2)) ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyUparalk !_________________________________________________________________ ! FUNCTION curlyUperplk(l,k) RESULT(XOUT) ! ! Computes \mathcal{U}_{\perp}^{lk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! IF(l .eq. 0) THEN XOUT = bafm_/2*(kerneln(bafm_,k) - kerneln(bafm_,k-1)) ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyUperplk !_________________________________________________________________ ! FUNCTION curlyTperplk(l,k) RESULT(XOUT) ! ! Computes \mathcal{T}_{\perp}^{lk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! IF ( l .eq. 2 ) THEN XOUT = kerneln(bafm_,k)/SQRT(TO_FM(2)) ELSEIF( l .eq. 0) THEN XOUT = (2*k)*kerneln(bafm_,k) - k*kerneln(bafm_,k-1 ) - (k+1)*kerneln(bafm_,k+1) ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyTperplk !_________________________________________________________________ ! FUNCTION curlyF0paralk(l,k) RESULT(XOUT) ! ! Computes \mathcal{F}_{ab 0 \parallel}^{Slk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! ! local vars ! INTEGER :: n,s,r,t,t1 TYPE(FM),SAVE :: L1tt1_,d0nks_,T4inv,nupara_,kern_,NX0 CHARACTER(len=T4len) :: T4str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! NX0 = 1/SQRT(TO_FM(2)**l*FACTORIAL(l)) ! nloop: DO n=0,naFLR kern_ = kerneln(bafm_, n) sloop: DO s=0, n+k d0nks_ = ALLX2L(n,k,s,0) IF( s .ge. 1 .or. l .ge. 1) THEN tloop: DO t=0, s +floor(real(l/2)) T4str_ = GET_T4(1,t,l,s) ! ... T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*t +1)/2)*TO_FM(T4str_) IF(T4inv .ne. TO_FM(0)) THEN t1loop: DO t1=0,t L1tt1_ = Lpjl(real(1,dp),real(t,dp),real(t1,dp)) nupara_ = nuabparak(t1+3) XOUT = XOUT + NX0*L1tt1_*d0nks_*T4inv*nupara_*kern_ ENDDO t1loop ENDIF ENDDO tloop ENDIF ENDDO sloop ENDDO nloop ! XOUT =- 8*XOUT*( 1 + chifm_**2)/3/SQRT(PIFM) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyF0paralk !_________________________________________________________________ ! FUNCTION curlyF1paralk(l,k) RESULT(XOUT) ! ! Computes \mathcal{F}_{ab 1 \parallel}^{Slk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! local vars ! INTEGER :: n,s,r,t,t1 TYPE(FM),SAVE :: L0tt1_,d1nks_,T4inv,nupara_,kern_,NX0 CHARACTER(len=T4len) :: T4str_ ! ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! NX0 = 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+1 d1nks_ = ALLX2L(n,k,s,1) tloop: DO t=0, s +floor(real(l)/2) T4str_ = GET_T4(0,t,l,s) ! ... T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t))*TO_FM( 1)/2/FACTORIAL(TO_FM(2*t +1)/2)*TO_FM(T4str_) IF(T4inv .ne. TO_FM(0)) THEN t1loop: DO t1=0,t L0tt1_ = Lpjl(real(0,dp),real(t,dp),real(t1,dp)) nupara_ = nuabparak(t1+2) XOUT = XOUT + NX0*L0tt1_*d1nks_*T4inv*nupara_*kern_*bafm_/2/(n+1) ENDDO t1loop ENDIF ENDDO tloop ENDDO sloop ENDDO nloop ! XOUT =- 8*XOUT*( 1 + chifm_**2)/SQRT(PIFM) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyF1paralk !_________________________________________________________________ ! FUNCTION curlyFElk(l,k) RESULT(XOUT) ! ! Computes \mathcal{F}_{ab E}^{Slk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! local vars ! INTEGER :: n,s,r,t,t1 TYPE(FM),SAVE :: L0tt1_,d0nks_,T4inv,nupara_,kern_,NX0 CHARACTER(len=T4len) :: T4str_ ! ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! NX0 = 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) tloop: DO t=0, s +floor(real(l)/2) T4str_ = GET_T4(0,t,l,s) ! ... T4inv = SQRT(PIFM)*(TO_FM(2)**l)*FACTORIAL(TO_FM(l))*FACTORIAL(TO_FM(t))*TO_FM( 1)/2/FACTORIAL(TO_FM(2*t +1)/2)*TO_FM(T4str_) IF(T4inv .ne. TO_FM(0)) THEN t1loop: DO t1=0, t L0tt1_ = Lpjl(real(0,dp),real(t,dp),real(t1,dp)) XOUT = XOUT + NX0*T4inv*kern_*L0tt1_*d0nks_*(Erfkab_array(t1) - chifm_*(1+chifm_**2)*ekab_array(t1+1))/chifm_**2 ENDDO t1loop ENDIF ENDDO tloop ENDDO sloop ENDDO nloop ! XOUT = - 8*XOUT/SQRT(PIFM) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyFElk !_________________________________________________________________ ! FUNCTION overlinecurlyVparalk(l,k) RESULT(XOUT) ! ! Computes \overline{\mathcal{V}}_{\parallel ab}^{lk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! local vars ! INTEGER :: n,s TYPE(FM),SAVE ::NX1,kern_,d0nks_,curlyVls_ ! ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! NX1 = - TO_FM(4)*(TO_FM(1) + chifm_**2)*thetaabfm_/TO_FM(3)/SQRT(PIFM) ! nloop: DO n=0,naFLR +1 kern_ = kerneln(bafm_,n) sloop: DO s=0,n+k d0nks_ = ALLX2L(n,k,s,0) curlyVls_ = curlyVab(l,s) XOUT = XOUT + d0nks_*kern_*curlyVls_*NX1 ENDDO sloop ENDDO nloop ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION overlinecurlyVparalk !_________________________________________________________________ ! FUNCTION overlinecurlyVperplk(l,k) RESULT(XOUT) ! ! Computes \overline{\mathcal{V}}_{\perp ab}^{lk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! local vars ! INTEGER :: n,s,f TYPE(FM),SAVE ::kern_,d1nks_,curlyIabls_,ba_,X1 ! ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! nloop: DO n=0,naFLR kern_ = kerneln(bafm_,n) floop1: DO f =0, n + k curlyIabls_ = curlyIabls(l,f) X1 = ALL2L(n,k,f,0) XOUT = XOUT + bafm_/TO_FM(2)*kern_*X1*curlyIabls_ ENDDO floop1 floop2: DO f =0, n + k +1 curlyIabls_ = curlyIabls(l,f) X1 = ALL2L(n+1,k,f,0) XOUT = XOUT - bafm_/TO_FM(2)*kern_*X1*curlyIabls_ ENDDO floop2 ENDDO nloop ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION overlinecurlyVperplk !_________________________________________________________________ ! ! FUNCTION overlinecurlyWlk(l,k) RESULT(XOUT) ! ! Computes \overline{\mathcal{W}}_{ab}^{lk} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l,k TYPE(FM) :: XOUT ! ! local vars ! INTEGER :: n,s TYPE(FM),SAVE ::kern_,d0nks_,curlyWabls_ ! ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0.0') ! ! nloop: DO n=0,naFLR kern_ = kerneln(bafm_,n) !! kern_ = kerneln(TO_FM('0.0'),n) sloop: DO s=0,n+k d0nks_ = ALLX2L(n,k,s,0) curlyWabls_= curlyWab(l,s) XOUT = XOUT + d0nks_*kern_*curlyWabls_ ENDDO sloop ENDDO nloop ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION overlinecurlyWlk ! !_______________________________________________ ! ! IMPROVED DRIFT-KINETIC SUGAMA OPERATOR !_______________________________________________ ! ! FUNCTION DeltaMlk(l,k) RESULT(XOUT) ! Compute the test correction matrix element \Delta M_{ab}^{lk} = M_{ab}^{Llk} - M_{ab}^{Slk} ! where M_{ab}^{Alk} is the Braginskii matrix element for the collision operator model $A$. ! Used in the improved Sugamam collision operator model. ! ! Note: The Sugama matrix element (l,k) are symmetric, a consquence of the self-adjointess relation of ! IMPLICIT NONE INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! Local vars. ! TYPE(FM), SAVE :: MabLlk_, MabSlk_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! _______________ COULOMB BRAGINKSII MATRIX ______________ MabLlk_ = MabLlk(l,k) ! ! _______________ SUGAMA BRAGINKSII MATRIX ______________ MabSlk_ = MabSlk(l,k) ! ! Test Correction Matrix element !______________________________________________ ! XOUT = MabLlk_ - MabSlk_ ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION DeltaMlk ! !_________________________________________________________________ ! FUNCTION MabLlk(l,k) RESULT(XOUT) ! Computes the test correction matrix element M_{ab}^{Llk} if the linearized test Landau coll. op. ! ! IMPLICIT NONE INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! Local vars. TYPE(FM), SAVE :: Lll1,nuT1k1l1 INTEGER :: l1 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) DO l1=0,l Lll1=Lpjl(real(1,dp),real(l,dp),real(l1,dp)) nuT1k1l1 = nusTabpjgd(1,k,1,l1) XOUT = XOUT + 2*Lll1*nuT1k1l1/3 ENDDO XOUT = TO_FM(3)*SQRT(PIFM)*XOUT/TO_FM(4) ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION MabLlk ! !_________________________________________________________________ ! FUNCTION MabSlk(l,k) RESULT(XOUT) ! Compute the test correction matrix element M_{ab}^{Slk} ! IMPLICIT NONE INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! Local vars. TYPE(FM), SAVE :: MabSlk0,MabSlk1,MabSlk2,nusS,Lll1 INTEGER :: l1 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! MabSlk0 = TO_FM(0) MabSlk1 = TO_FM(0) MabSlk2 = TO_FM(0) ! DO l1=0,l Lll1=Lpjl(real(1,dp),real(l,dp),real(l1,dp)) nusS = nusSabpjd(1,k,l1) MabSlk0 = MabSlk0 + 2*Lll1*nusS/3 ENDDO ! ! MabSlk1 = - 8*(thetaabfm_ - 1)/3/SQRT(PIFM)*(1 + chifm_**2)*(uparak(l)*overlinenuSabk3(k)+ uparak(k)*overlinenuSabk3(l)) ! ... this has been corrected with an additional 8/3/SQRT(PI) factor ! ! MabSlk2 = - chifm_*(thetaabfm_ -1)**2/SQRT(TO_FM(1 + chifm_**2))*uparak(l)*uparak(k) MabSlk2 = 4*MabSlk2/SQRT(PIFM)/3 ! ... rescale the collision frequency ! - XOUT = MabSlk0 + MabSlk1 + MabSlk2 + XOUT = MabSlk0 + TO_FM(2)*MabSlk1 + TO_FM(2)*MabSlk2 XOUT = TO_FM(3)*SQRT(PIFM)*XOUT/TO_FM(4) ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION MabSlk !______________________________________________ ! FUNCTION DeltaNlk(l,k) RESULT(XOUT) ! Compute the field correction matrix element \Delta N_{ab}^{lk} = N_{ab}^{Llk} - N_{ab}^{Slk} where N_{ab }^{Alk} is the ! Braginskii matrix element for the collision operator model $A$. ! Used in the improved Sugamam collision operator model. ! IMPLICIT NONE INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! Local vars. TYPE(FM), SAVE :: NabLlk_,NabSlk_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT =TO_FM(0) ! ! Coulomb Correction Matrix NabLlk_ = NabLlk(l,k) ! ! Sugama Correction Matric NabSlk_ = NabSlk(l,k) ! XOUT = NabLlk_ - NabSlk_ ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION DeltaNlk !____________________________________________ ! FUNCTION NabLlk(l,k) RESULT(XOUT) ! Compute the field correction matrix element N_{ab}^{Llk} with the linearized test part of Coulomb ! IMPLICIT NONE INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! Local vars. TYPE(FM), SAVE :: Lll1,nuF1k1l1 INTEGER :: l1 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! _______________ COULOMB BRAGINKSII MATRIX ______________ ! XOUT = TO_FM(0) ! DO l1=0,l ! Lll1=Lpjl(real(1,dp),real(l,dp),real(l1,dp)) ! nuF1k1l1 = nusFabpjgd(1,k,1,l1) ! XOUT = XOUT + TO_FM(2)*Lll1*nuF1k1l1/TO_FM(3)/(chifm_**(2*l1)) ! ENDDO XOUT = TO_FM(0) DO l1=0,l Lll1=Lpjl(real(1,dp),real(l,dp),real(l1,dp)) nuF1k1l1 = nusFabpjgd(1,k,1,l1) XOUT = XOUT + TO_FM(2)*chifm_*Lll1*nuF1k1l1/TO_FM(3) ENDDO XOUT = TO_FM(3)*SQRT(PIFM)*XOUT/TO_FM(4) ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION NabLlk ! !______________________________________________________ ! FUNCTION NabSlk(l,k) RESULT(XOUT) ! Compute the field correction matrix element N_{ab}^{Slk} ! IMPLICIT NONE INTEGER,INTENT(IN) :: l,k TYPE(FM) :: XOUT ! Local vars. TYPE(FM), SAVE :: curlyVabl_,curlyVbak_,NX0 INTEGER :: l1 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! !____________ SUGAMA BRAGKINSKII MATRIX ELEMENT _________ ! NX0 = - 3*SQRT(PIFM)/4*( 1 + chifm_**2)**1.5_dp/chifm_/sigmafm_/( taufm_ + chifm_**2)! ... m_b / \gamma_{ab} curlyVabl_ = curlyVabl(l) curlyVbak_ = curlyVbal(k) ! !XOUT = -2*curlyVabl_*curlyVbak_*NX0 ! ... There should be an extra factor 2 ... why ? XOUT = -curlyVabl_*curlyVbak_*NX0 !XOUT = TO_FM(4)*XOUT/TO_FM(3)/SQRT(PIFM) ! ... rescale the collision freq. !XOUT = TO_FM(3)*SQRT(PIFM)/TO_FM(4) XOUT = TO_FM(3)*SQRT(PIFM)*XOUT/TO_FM(4) XOUT = TO_FM(2)*XOUT !XOUT = XOUT/chifm_ ! Why 1/chi ? ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION NabSlk !_________________________________________________________________ ! FUNCTION overlinenuSabk3(k) RESULT(XOUT) ! Term sum_{l=0}^k L_{kl}^1 \bar\nu_{ab}^{// l+3} ! IMPLICIT NONE ! INTEGER,INTENT(in) :: k TYPE(FM) :: XOUT ! INTEGER :: k1 TYPE(FM), SAVE :: Lkk1,nu_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! DO k1=0,k Lkk1=Lpjl(real(1,dp),real(k,dp),real(k1,dp)) nu_ = nuabparak(k1 + 3) XOUT = XOUT + Lkk1*nu_ ENDDO ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION overlinenuSabk3 !_________________________________________________________________ ! FUNCTION uparak(k) RESULT(XOUT) ! ! computes the equilibrium parallel flow vector ! ! IMPLICIT NONE ! INTEGER,INTENT(in) :: K TYPE(FM) :: XOUT ! ! local vars. INTEGER :: k1 TYPE(FM), SAVE :: Lkk1 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! ! DO k1=0,k ! Lkk1=Lpjl(real(1,dp),real(k,dp),real(k1,dp)) ! XOUT = XOUT + Lkk1*FACTORIAL(TO_FM(2*k1 + 1 )/2) ! ENDDO ! IF( k .eq. 0) THEN XOUT = TO_FM(1) ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION uparak ! !_________________________________________________________________ ! FUNCTION curlyVabl(l) RESULT(XOUT) ! ! computes \mathcal{V}_{ab}^l used in the correction matrix element of the field part of the Sugama operator ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l TYPE(FM) :: XOUT ! ! local vars. INTEGER :: l1 TYPE(FM), SAVE :: NX0,NX1,Lll1 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) NX1 = TO_FM(0) ! NX0 = - 8*thetaabfm_*(1 + chifm_**2)/3/SQRT(PIFM) ! DO l1 = 0,l Lll1=Lpjl(real(1,dp),real(l,dp),real(l1,dp)) XOUT = XOUT + Lll1*( TO_FM(1)/chifm_**2 & *(Erfkab_array(l1) - chifm_*ekab_array(l1+1) )) ENDDO ! IF( l .eq. 0) THEN ! uparak(l) == 1 iff l=0, else uparak(l) == 0 NX1 = chifm_*(thetaabfm_ -1)/(1 + chifm_**2)**(1.5_dp)/TO_FM(2) ENDIF XOUT = NX0*(XOUT + NX1) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyVabl ! !________________________________________________ ! FUNCTION curlyVbal(l) RESULT(XOUT) ! ! computes \mathcal{V}_{ba}^l ! IMPLICIT NONE ! INTEGER,INTENT(in) :: l TYPE(FM) :: XOUT ! ! local vars. INTEGER :: l1 TYPE(FM), SAVE :: NX0,NX1,Lll1 ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) NX1 = TO_FM(0) ! NX0 = - 8*thetabafm_*(1 + TO_FM(1)/chifm_**2)/3/SQRT(PIFM)*chifm_**3*sigmafm_**2 ! DO l1 = 0,l Lll1=Lpjl(real(1,dp),real(l,dp),real(l1,dp)) XOUT = XOUT + Lll1*( chifm_**2 & *(Erfkba_array(l1) - ekba_array(l1 +1 )/chifm_ )) ENDDO ! IF( l .eq. 0) THEN ! uparak(l) == 1 iff l=0, else uparak(l) == 0 NX1 = TO_FM(1)/chifm_*(thetabafm_ -1)/(1 + TO_FM(1)/chifm_**2)**(1.5_dp)/TO_FM(2) ENDIF XOUT = NX0*(XOUT + NX1) ! !XOUT = TO_FM(4)*XOUT/3/SQRT(PIFM)*chifm_**3*sigmafm_**2 ! ... rescale tau_{ab}^S collision freq. ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION curlyVbal ! !_______________________________________________ ! ! IMPROVED GYROKINETIC SUGAMA OPERATOR !_______________________________________________ ! FUNCTION CurlIparapjk(p,j,k, bs) RESULT(XOUT) ! ! Compute \mathcal I_{\paralell s}^{pjk} ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: p, j, k TYPE(FM), INTENT(IN) :: bs TYPE(FM) :: XOUT ! Local vars INTEGER :: n, s TYPE(FM), SAVE :: NX0, T4inv, kern_, d0njs CHARACTER(len=T4len) :: T4str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) XOUT = TO_FM(0) NX0 = TO_FM(2)/TO_FM(3)/SQRT(PIFM)*FACTORIAL(TO_FM(2*k+1)/2)/FACTORIAL(TO_FM(k))/SQRT(TO_FM(2)**p * FACTORIAL(p)) nloop: DO n=0, naFLR kern_ = kerneln(bs, n) sloop: DO s=0, n+j IF( (p .ge. 1 .or. s .ge. 1) .and. (s+p/2 .ge. k) ) THEN T4str_ = GET_T4(1,k,p,s) T4inv = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(k))*TO_FM(2 + 1)/2/FACTORIAL(TO_FM(2+2*k +1)/2)*TO_FM(T4str_) d0njs = ALLX2L(n,j,s,0) XOUT = XOUT+NX0*T4inv*kern_*d0njs END IF END DO sloop END DO nloop CALL FM_EXIT_USER_FUNCTION(XOUT) END FUNCTION CurlIparapjk FUNCTION CurlIperppjk(p,j,k,bs) RESULT(XOUT) ! ! Compute \mathcal I_{\perp s}^{pjk} ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: p, j, k TYPE(FM), INTENT(IN) :: bs TYPE(FM) :: XOUT ! Local vars INTEGER :: n, s, r, q TYPE(FM), SAVE :: NX0, NX1, T4inv, kern_, d1njs, Lkq CHARACTER(len=T4len) :: T4str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) XOUT = TO_FM(0) NX0 = TO_FM(1)/SQRT(PIFM)/SQRT(TO_FM(2)**p * FACTORIAL(p))*bs nloop: DO n=0, naFLR kern_ = kerneln(bs, n) sloop: DO s=0, n+j+1 d1njs = ALLX2L(n,j,s,1) rloop: DO r=0, s+p/2 T4str_ = GET_T4(0,r,p,s) T4inv = SQRT(PIFM)*(TO_FM(2)**p)*FACTORIAL(TO_FM(p))*FACTORIAL(TO_FM(k))*TO_FM(1)/2/FACTORIAL(TO_FM(2*k +1)/2)*TO_FM(T4str_) qloop: DO q=0, k Lkq=Lpjl(real(1,dp),real(k,dp),real(q,dp)) NX1 = FACTORIAL(TO_FM(2*1+1)/2)/(TO_FM(n+1))*GENERALIZED_BINOMIAL(TO_FM(r-q-1), TO_FM(r)) XOUT = XOUT + NX0*NX1*kern_*T4inv*d1njs*Lkq END DO qloop END DO rloop END DO sloop END DO nloop CALL FM_EXIT_USER_FUNCTION(XOUT) END FUNCTION CurlIperppjk FUNCTION GENERALIZED_BINOMIAL(x,y) RESULT(XOUT) IMPLICIT NONE ! TYPE(FM), INTENT(IN) :: x, y TYPE(FM) :: XOUT ! Local vars TYPE(FM), SAVE :: x_, y_, sign_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) x_ = x y_ = y sign_ = TO_FM(1) IF( x .lt. 0 .or. y .lt. 0) THEN IF(y .ge. 0) THEN sign_ = TO_FM(-1)**y x_ = -x + y - 1 y_ = y ELSEIF(y .le. x) THEN sign_ = TO_FM(-1)**(x-y) x_ = -y - 1 y_ = x - y ELSE sign_ = TO_FM(0) END IF END IF IF(x .lt. y .or. sign_ .eq. TO_FM(0)) THEN XOUT = TO_FM(0) ELSE XOUT = sign_*FACTORIAL(x)/FACTORIAL(y)/FACTORIAL(x-y) END IF CALL FM_EXIT_USER_FUNCTION(XOUT) END FUNCTION GENERALIZED_BINOMIAL ! FUNCTION curlyIparaalkj(l,k,j) RESULT(XOUT) ! ! ! ! ! IMPLICIT NONE ! ! ! INTEGER,INTENT(in) :: l,k,j ! TYPE(FM) :: XOUT ! ! ! INTEGER :: n,s,h ! TYPE(FM),SAVE :: kern_,d0nks_,NX0,NX1,T4inv ! CHARACTER(len=T4len) :: T4str_ ! ! ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! ! XOUT = TO_FM(0) ! ! ! IF( l .ge. 1 .or. k .ge. 1 )THEN ! ! ! NX0 = 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) ! hloop: DO h =0, s + floor(real(l)/2) ! IF( h .eq. j) THEN ! 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)) THEN ! XOUT = XOUT + T4inv*NX0*kern_*d0nks_*FACTORIAL(TO_FM(2*j + 3)/2)/FACTORIAL(TO_FM(j)) ! ENDIF ! ENDIF ! ENDDO hloop ! ENDDO sloop ! ENDDO nloop ! ! ! ENDIF ! ! ! XOUT = 2*XOUT/3/SQRT(PIFM) ! ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! ! ! END FUNCTION curlyIparaalkj ! !_______________________________________________ ! ! ! FUNCTION curlyIperpalkj(l,k,j) RESULT(XOUT) ! ! ! ! ! IMPLICIT NONE ! ! ! INTEGER,INTENT(in) :: l,k,j ! TYPE(FM) :: XOUT ! ! ! ! local vars. ! INTEGER :: n,s,h,j1 ! TYPE(FM), SAVE :: kern_,d1nks_,T4inv,L1jj1,binom_,NX0 ! CHARACTER(len=T4len) :: T4str_ ! ! ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! ! ! XOUT = TO_FM(0) ! ! ! ! ! ! ! NX0 = 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 +1 ! d1nks_ = ALLX2L(n,k,s,1) ! hloop: DO h =0, s + floor(real(l)/2) ! 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)/2/FACTORIAL(TO_FM(2*h +1)/2)*TO_FM(T4str_) ! IF(T4inv .ne. TO_FM(0)) THEN ! j1loop: DO j1=0,j ! L1jj1=Lpjl(real(1,dp),real(j,dp),real(j1,dp)) ! ! compute binomial coeff ! binom_ = BinomialFM(h - j1 -1,h) ! XOUT = XOUT + T4inv*NX0*bafm_*kern_/( n +1)/SQRT(PIFM)*L1jj1*d1nks_*FACTORIAL((2*j1 + 1 )/2)*binom_ ! ENDDO j1loop ! ENDIF ! ENDDO hloop ! ENDDO sloop ! ENDDO nloop ! ! ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! ! ! END FUNCTION curlyIperpalkj !_______________________________________________ ! ! FUNCTION Tabmpfjl(m,p,f,j,l) RESULT(XOUT) ! ! evaluate T_ab(m,p,f,j,l) for the gyrokinetic full Coulomb ! IMPLICIT NONE ! INTEGER,INTENT(in)::m,p,f,j,l TYPE(FM) :: XOUT ! ! local vars. ! integer :: t,d TYPE(FM),SAVE:: T5inv,Lptd CHARACTER(len=T5len) :: T5str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! 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-abs(m)))*TO_FM(T5STR_) IF(T5inv .ne. TO_FM('0')) THEN do d =0,t Lptd = Lpjl(real(p,dp),real(t,dp),real(d,dp)) XOUT = XOUT + Lptd*nusTabpjgd(p,j,p,d)*T5inv enddo ENDIF enddo ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION Tabmpfjl FUNCTION Fabmpfjl(m,p,f,j,l) RESULT(XOUT) ! ! evaluate F_ab(m,p,f,j,l) for the gyrokinetic full Coulomb ! IMPLICIT NONE ! INTEGER,INTENT(in)::m,p,f,j,l TYPE(FM) :: XOUT ! ! local vars. ! integer :: t,d TYPE(FM),SAVE:: T5inv,Lptd CHARACTER(len=T5len) :: T5str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! 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-abs(m)))*TO_FM(T5STR_) IF(T5inv .ne. TO_FM('0')) THEN do d =0,t Lptd = Lpjl(real(p,dp),real(t,dp),real(d,dp)) XOUT = XOUT + Lptd*nusFabpjgd(p,j,p,d)*T5inv enddo ENDIF enddo ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION Fabmpfjl ! FUNCTION Selfmpfjl(m,p,f,j,l) RESULT(XOUT) ! ! evaluate S_aa(m,p,f,j,l) for the gyrokinetic full Coulomb ! IMPLICIT NONE ! INTEGER,INTENT(in)::m,p,f,j,l TYPE(FM) :: XOUT ! ! local vars. ! integer :: t,d TYPE(FM),SAVE:: T5inv,Lptd CHARACTER(len=T5len) :: T5str_ ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM(0) ! 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-abs(m)))*TO_FM(T5STR_) IF(T5inv .ne. TO_FM('0')) THEN do d =0,t Lptd = Lpjl(real(p,dp),real(t,dp),real(d,dp)) XOUT = XOUT + Lptd*(nusFabpjgd(p,j,p,d) + nusTabpjgd(p,j,p,d))*T5inv enddo ENDIF enddo ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION Selfmpfjl ! !_______________________________________________ ! ! ZEROTH-ORDER Hirshman-Sigmar !_______________________________________________ FUNCTION nuabSDk(k) result(XOUT) ! Slowing down speed integrated collision frequency ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: k TYPE(FM) :: XOUT ! ! LOCAL VAR. ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! IF(k .ge. 1) THEN ! ... k >= 1 XOUT = taufm_/chifm_**2*(1 + TO_FM(1)/sigmafm_)*(Erfkab_array(k-2) - chifm_*ekab_array(k-1)) ENDIF ! IF(k .eq. 0) THEN ! ... k =0 XOUT = TO_FM('1.06567995070710404713815879845981153908302309423062532484677') ENDIF ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nuabSDk !_______________________________________________ ! FUNCTION nubaSDk(k) result(XOUT) ! Slowing down speed integrated collision frequency IMPLICIT NONE ! INTEGER,INTENT(IN) :: k TYPE(FM) :: XOUT ! ! LOCAL VAR. ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = chifm_**2/taufm_*(1 + sigmafm_)*(Erfkba_array(k-2) - ekba_array(k-1)/chifm_) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nubaSDk ! !_______________________________________________ ! ! Abel collision frequuency !_______________________________________________ ! ! FUNCTION nuEk(k) result(XOUT) ! like-species energy frequency ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: k TYPE(FM) :: XOUT ! ! LOCAL VAR. ! CALL FM_ENTER_USER_FUNCTION(XOUT) ! XOUT = TO_FM('0') XOUT = TO_FM(2)*(nuabSDk(k) -nuabDk(k) ) - nuabparak(k) ! CALL FM_EXIT_USER_FUNCTION(XOUT) ! END FUNCTION nuEk ! END MODULE speed_functions