subroutine include_CoulombCollisionsLA_GC_p(spp,achunk,tt,params, &
Y_R,Y_PHI,Y_Z,Ppll,Pmu,me,flagCon,flagCol,F,P,E_PHI,ne,Te,PSIp)
TYPE(SPECIES), INTENT(INOUT) :: spp
TYPE(PROFILES), INTENT(IN) :: P
TYPE(FIELDS), INTENT(IN) :: F
TYPE(KORC_PARAMS), INTENT(INOUT) :: params
INTEGER, INTENT(IN) :: achunk
REAL(rp), DIMENSION(achunk), INTENT(INOUT) :: Ppll
REAL(rp), DIMENSION(achunk), INTENT(INOUT) :: Pmu
REAL(rp), DIMENSION(achunk) :: Bmag
REAL(rp), DIMENSION(achunk) :: B_R,B_PHI,B_Z
REAL(rp), DIMENSION(achunk) :: curlb_R,curlb_PHI,curlb_Z
REAL(rp), DIMENSION(achunk) :: gradB_R,gradB_PHI,gradB_Z,ntot
REAL(rp), DIMENSION(achunk) :: E_R,E_Z,E_PHI_LAC
REAL(rp), DIMENSION(achunk), INTENT(OUT) :: E_PHI,ne,Te,PSIp
REAL(rp), DIMENSION(achunk), INTENT(IN) :: Y_R,Y_PHI,Y_Z
INTEGER(is), DIMENSION(achunk), INTENT(INOUT) :: flagCol
INTEGER(is), DIMENSION(achunk), INTENT(INOUT) :: flagCon
REAL(rp), INTENT(IN) :: me
REAL(rp), DIMENSION(achunk) :: Zeff
REAL(rp), DIMENSION(achunk) :: nAr0,nAr1,nAr2,nAr3
REAL(rp), DIMENSION(achunk) :: nD,nD1
REAL(rp), DIMENSION(achunk,2) :: dW
REAL(rp), DIMENSION(achunk,2) :: rnd1
REAL(rp) :: dt,time
REAL(rp), DIMENSION(achunk) :: pm,pm0
REAL(rp), DIMENSION(achunk) :: dp
REAL(rp), DIMENSION(achunk) :: xi,xi0
REAL(rp), DIMENSION(achunk) :: dxi
REAL(rp), DIMENSION(achunk) :: v,gam
!! speed of particle
REAL(rp), DIMENSION(achunk) :: CAL
REAL(rp) , DIMENSION(achunk) :: dCAL
REAL(rp), DIMENSION(achunk) :: CFL
REAL(rp), DIMENSION(achunk) :: CBL
REAL(rp), DIMENSION(achunk) :: SC_p,SC_xi,BREM_p
REAL(rp) :: kappa
integer :: cc,ii
integer(ip),INTENT(IN) :: tt
REAL(rp), DIMENSION(achunk,params%num_impurity_species) :: nimp
dt=cparams_ss%coll_per_dump_dt
time=params%init_time+(params%it-1)*params%dt+ &
tt*cparams_ss%coll_per_dump_dt
if (params%field_eval.eq.'eqn') then
call analytical_fields_GC_p(achunk,F,Y_R,Y_PHI, &
Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,PSIp)
else if (params%field_eval.eq.'interp') then
#ifdef PSPLINE
if (F%axisymmetric_fields) then
if (F%Bflux) then
if (.not.params%SC_E) then
call calculate_GCfields_p(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp)
else
call calculate_GCfields_p_FS(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp)
end if
else if (F%ReInterp_2x1t) then
call calculate_GCfieldswE_p(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp)
else if (F%Bflux3D) then
call calculate_GCfields_2x1t_p(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp,time)
else if (F%dBfield) then
call calculate_2DBdBfields_p(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp)
else
call interp_fields_p(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon)
end if
else
if (F%dBfield) then
call calculate_2DBdBfields_p(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp)
else
call interp_fields_3D_p(achunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, &
E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
gradB_R,gradB_PHI,gradB_Z,flagCon)
end if
endif
#endif
if(.not.F%ReInterp_2x1t) call add_analytical_E_p(params,tt,F,E_PHI,Y_R,Y_Z)
end if
if (params%profile_model(1:10).eq.'ANALYTICAL') then
call analytical_profiles_p(achunk,time,params,Y_R,Y_Z,P,F, &
ne,Te,Zeff,PSIp)
else if (params%profile_model(1:8).eq.'EXTERNAL') then
#ifdef PSPLINE
if (params%profile_model(10:10).eq.'H') then
call interp_Hcollision_p(achunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, &
nAr0,nAr1,nAr2,nAr3,nD,nD1,flagCon)
do cc=1_idef,achunk
nimp(cc,1)=nAr0(cc)
nimp(cc,2)=nAr1(cc)
nimp(cc,3)=nAr2(cc)
nimp(cc,4)=nAr3(cc)
nimp(cc,5)=nD(cc)
nimp(cc,6)=nD1(cc)
end do
else
call interp_FOcollision_p(achunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff,flagCon)
endif
#endif
end if
E_PHI_LAC=E_PHI
if (.not.params%FokPlan) E_PHI=0._rp
!$OMP SIMD
! !$OMP& aligned (pm,xi,v,Ppll,Bmag,Pmu)
do cc=1_idef,achunk
Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+B_Z(cc)*B_Z(cc))
! Transform p_pll,mu to P,eta
pm(cc) = SQRT(Ppll(cc)*Ppll(cc)+2*me*Bmag(cc)*Pmu(cc))
pm0(cc)=pm(cc)
xi(cc) = Ppll(cc)/pm(cc)
xi0(cc)=xi(cc)
gam(cc) = sqrt(1+pm(cc)*pm(cc))
v(cc) = pm(cc)/gam(cc)
! normalized speed (v_K=v_P/c)
end do
!$OMP END SIMD
! write(output_unit_write,'("ne: "E17.10)') ne
! write(output_unit_write,'("Te: "E17.10)') Te
! write(output_unit_write,'("Bmag: "E17.10)') Bmag
! write(output_unit_write,'("v: ",E17.10)') v
! write(output_unit_write,'("xi: ",E17.10)') xi
! write(output_unit_write,'("size(E_PHI_GC): ",I16)') size(E_PHI)
!$OMP SIMD
! !$OMP& aligned(rnd1,dW,CAL,dCAL,CFL,CBL,v,ne,Te,Zeff,dp, &
! !$OMP& flagCon,flagCol,dxi,xi,pm,Ppll,Pmu,Bmag)
do cc=1_idef,achunk
#ifdef PARALLEL_RANDOM
rnd1(cc,1) = get_random()
rnd1(cc,2) = get_random()
! rnd1(:,1) = get_random_mkl()
! rnd1(:,2) = get_random_mkl()
#else
call RANDOM_NUMBER(rnd1)
#endif
dW(cc,1) = SQRT(3*dt)*(-1+2*rnd1(cc,1))
dW(cc,2) = SQRT(3*dt)*(-1+2*rnd1(cc,2))
! write(output_unit_write,'("dW1: ",E17.10)') dW(cc,1)
! write(output_unit_write,'("dW2: ",E17.10)') dW(cc,2)
if (params%profile_model(10:10).eq.'H') then
CAL(cc) = CA_SD(v(cc),ne(cc),Te(cc))
dCAL(cc)= dCA_SD(v(cc),me,ne(cc),Te(cc))
CFL(cc) = CF_SD_FIO(params,v(cc),ne(cc),Te(cc),nimp(cc,:))
CBL(cc) = (CB_ee_SD(v(cc),ne(cc),Te(cc),Zeff(cc))+ &
CB_ei_SD_FIO(params,v(cc),ne(cc),Te(cc),nimp(cc,:),Zeff(cc)))
! using *_FIO routine to include additions due to Hollmann species
! profile
else
CAL(cc) = CA_SD(v(cc),ne(cc),Te(cc))
dCAL(cc)= dCA_SD(v(cc),me,ne(cc),Te(cc))
CFL(cc) = CF_SD(params,v(cc),ne(cc),Te(cc))
CBL(cc) = (CB_ee_SD(v(cc),ne(cc),Te(cc),Zeff(cc))+ &
CB_ei_SD(params,v(cc),ne(cc),Te(cc),Zeff(cc)))
endif
if (.not.energy_diffusion) then
CAL(cc)=0._rp
dCAL(cc)=0._rp
endif
if (.not.pitch_diffusion) then
CBL(cc)=0._rp
endif
dp(cc)=REAL(flagCol(cc))*REAL(flagCon(cc))* &
((-CFL(cc)+dCAL(cc)+E_PHI(cc)*xi(cc))*dt+ &
sqrt(2.0_rp*CAL(cc))*dW(cc,1))
dxi(cc)=REAL(flagCol(cc))*REAL(flagCon(cc))* &
((-2*xi(cc)*CBL(cc)/(pm(cc)*pm(cc))+ &
E_PHI(cc)*(1-xi(cc)*xi(cc))/pm(cc))*dt- &
sqrt(2.0_rp*CBL(cc)*(1-xi(cc)*xi(cc)))/pm(cc)*dW(cc,2))
if(cparams_ss%sample_test) then
dp(cc)=0._rp
dxi(cc)=0._rp
end if
! write(output_unit_write,'("dp: ",E17.10)') dp(cc)
! write(output_unit_write,'("dxi: ",E17.10)') dxi(cc)
!write(6,*) 'gam,xi',gam(cc),xi(cc)
!write(6,*) 'dpE',E_PHI(cc)*xi(cc)*dt
!write(6,*) 'dpCF',CFL(cc)*dt
!write(6,*) 'dpCA',sqrt(2*CAL(cc)*dt)
!write(6,*) 'dxiE',E_PHI(cc)*(1-xi(cc)*xi(cc))/pm(cc)*dt
!write(6,*) 'dxiCB',2*xi(cc)*CBL(cc)/(pm(cc)*pm(cc))*dt
end do
!$OMP END SIMD
if (params%FokPlan.and.params%radiation) then
if(params%GC_rad_model.eq.'SDE') then
!$OMP SIMD
do cc=1_idef,achunk
SC_p(cc)=-gam(cc)*pm(cc)*(1-xi(cc)*xi(cc))/ &
(cparams_ss%taur/Bmag(cc)**2)
SC_xi(cc)=xi(cc)*(1-xi(cc)*xi(cc))/ &
((cparams_ss%taur/Bmag(cc)**2)*gam(cc))
kappa=2._rp*C_PI*C_RE**2._rp*C_ME*C_C**2._rp/ &
(params%cpp%length**2._rp*params%cpp%energy)
BREM_p(cc)=-2._rp*ne(cc)*kappa*Zeff(cc)*(Zeff(cc)+1._rp)* &
C_a/C_PI*(gam(cc)-1._rp)*(log(2._rp*gam(cc))-1._rp/3._rp)
!write(6,*) 'dpR',SC_p(cc)*dt
!write(6,*) 'dpB',BREM_p(cc)*dt
!write(6,*) 'dxiR',SC_xi(cc)*dt
if (.not.FP_bremsstrahlung) BREM_p(cc)=0._rp
dp(cc)=dp(cc)+(SC_p(cc)+BREM_p(cc))*dt* &
REAL(flagCol(cc))*REAL(flagCon(cc))
dxi(cc)=dxi(cc)+(SC_xi(cc))*dt* &
REAL(flagCol(cc))*REAL(flagCon(cc))
end do
!$OMP END SIMD
end if
end if
!$OMP SIMD
do cc=1_idef,achunk
pm(cc)=pm(cc)+dp(cc)
xi(cc)=xi(cc)+dxi(cc)
end do
!$OMP END SIMD
do cc=1_idef,achunk
! if (pm(cc)<0) pm(cc)=-pm(cc)
! Keep xi between [-1,1]
if (xi(cc)>1) then
! write(output_unit_write,'("High xi at: ",E17.10," with dxi: ",E17.10)') &
! time*params%cpp%time, dxi(cc)
xi(cc)=1-mod(xi(cc),1._rp)
else if (xi(cc)<-1) then
xi(cc)=-1-mod(xi(cc),-1._rp)
! write(output_unit_write,'("Low xi at: ",E17.10," with dxi: ",E17.10)') &
! time*params%cpp%time, dxi(cc)
endif
if ((pm(cc).lt.min(cparams_ss%p_min*cparams_ss%pmin_scale, &
p_therm)).and.flagCol(cc).eq.1_ip) then
! write(output_unit_write,'("Momentum less than zero")')
! stop
! write(output_unit_write,'("Particle not tracked at: ",E17.10," &
! & with xi: ",E17.10)') time*params%cpp%time, xi(cc)
flagCol(cc)=0_ip
end if
end do
! write(output_unit_write,'("rnd1: ",E17.10)') rnd1
! write(output_unit_write,'("flag: ",I16)') flag
! write(output_unit_write,'("CA: ",E17.10)') CAL
! write(output_unit_write,'("dCA: ",E17.10)') dCAL
! write(output_unit_write,'("CF ",E17.10)') CFL
! write(output_unit_write,'("CB: ",E17.10)') CBL
! write(output_unit_write,'("dp: ",E17.10)') dp
! write(output_unit_write,'("dxi: ",E17.10)') dxi
! write(output_unit_write,'("Ppll: ",E17.10)') Ppll
! write(output_unit_write,'("Pmu: ",E17.10)') Pmu
! write(output_unit_write,'("E_PHI_COL: ",E17.10)') E_PHI
! if (tt .EQ. 1_ip) then
! write(output_unit_write,'("dp_rad: ",E17.10)') &
! -gam(1)*pm(1)*(1-xi(1)*xi(1))/ &
! (cparams_ss%taur/Bmag(1)**2)*dt
! write(output_unit_write,'("dxi_rad: ",E17.10)') &
! xi(1)*(1-xi(1)*xi(1))/ &
! ((cparams_ss%taur/Bmag(1)**2)*gam(1))*dt
! end if
! if (tt .EQ. 1_ip) then
! write(output_unit_write,'("CA: ",E17.10)') CAL(1)
! write(output_unit_write,'("dCA: ",E17.10)') dCAL(1)
! write(output_unit_write,'("CF ",E17.10)') CFL(1)
! write(output_unit_write,'("CB: ",E17.10)') CBL(1)
! end if
#if DBG_CHECK
do cc=1_idef,achunk
if (ISNAN(xi(cc)).or.(abs(xi(cc)).gt.1._rp)) then
write(6,*) 'xi is NaN or >1 before LAC'
write(6,*) 'p0,xi0',pm0(cc),xi0(cc)
write(6,*) 'p,xi',pm(cc),xi(cc)
write(6,*) 'dp,dxi',dp(cc),dxi(cc)
write(6,*) 'CBL',CBL(cc)
write(6,*) 'v,ne,Te,Zeff',v(cc),ne(cc)*params%cpp%density,Te(cc)*params%cpp%temperature,Zeff(cc)
write(6,*) 'ppll,pmu,Bmag',Ppll(cc),Pmu(cc),Bmag(cc)
call korc_abort(24)
end if
end do
# endif
if (cparams_ss%avalanche) then
!$OMP SIMD
do cc=1_idef,achunk
ntot(cc)=ne(cc)
if (params%bound_electron_model.eq.'HESSLOW') then
if ((cparams_ms%Zj(1).eq.0.0).and. &
(neut_prof.eq.'UNIFORM')) then
ntot(cc)=ntot(cc)+cparams_ms%nz(1)* &
(cparams_ms%Zo(1)-cparams_ms%Zj(1))
else if ((cparams_ms%Zj(1).eq.0.0).and. &
(neut_prof.eq.'HOLLOW')) then
ntot(cc)=ntot(cc)+max(cparams_ms%nz(1)-ne(cc),0._rp)* &
(cparams_ms%Zo(1)-cparams_ms%Zj(1))
else
ntot(cc)=ntot(cc)+ne(cc)*cparams_ms%nz(1)/cparams_ms%ne* &
(cparams_ms%Zo(1)-cparams_ms%Zj(1))
endif
do ii=2,cparams_ms%num_impurity_species
ntot(cc)=ntot(cc)+ne(cc)*cparams_ms%nz(ii)/cparams_ms%ne* &
(cparams_ms%Zo(ii)-cparams_ms%Zj(ii))
end do
end if
end do
!$OMP END SIMD
!write(6,*) 'ntot',ntot*params%cpp%density
call large_angle_source(spp,params,achunk,F,Y_R,Y_PHI,Y_Z, &
pm,xi,ne,ntot,Te,Bmag,E_PHI_LAC,me,flagCol,flagCon)
end if
#if DBG_CHECK
do cc=1_idef,achunk
if (ISNAN(xi(cc)).or.(abs(xi(cc)).gt.1._rp)) then
write(6,*) 'xi is NaN or >1 after LAC'
write(6,*) 'p0,xi0',pm0(cc),xi0(cc)
write(6,*) 'p,xi',pm(cc),xi(cc)
write(6,*) 'dp,dxi',dp(cc),dxi(cc)
write(6,*) 'CBL',CBL(cc)
write(6,*) 'v,ne,Te,Zeff',v(cc),ne(cc)*params%cpp%density,Te(cc)*params%cpp%temperature,Zeff(cc)
call korc_abort(24)
end if
end do
# endif
!$OMP SIMD
do cc=1_idef,achunk
! Transform P,xi to p_pll,mu
Ppll(cc)=pm(cc)*xi(cc)
Pmu(cc)=(pm(cc)*pm(cc)-Ppll(cc)*Ppll(cc))/(2*me*Bmag(cc))
end do
!$OMP END SIMD
#if DBG_CHECK
do cc=1_idef,achunk
if (Pmu(cc).lt.0._rp) then
write(6,*) 'mu is negative'
write(6,*) 'p,xi',pm(cc),xi(cc)
write(6,*) 'V_PLL,V_MU',Ppll(cc),Pmu(cc)
call korc_abort(24)
endif
end do
#endif
E_PHI=E_PHI_LAC
end subroutine include_CoulombCollisionsLA_GC_p