all subroutines for CD computation grouped in eccd module
This commit is contained in:
parent
e92ff7cee1
commit
87de4c9cc2
11
Makefile
11
Makefile
@ -3,8 +3,8 @@ EXE=gray
|
|||||||
|
|
||||||
# Objects list
|
# Objects list
|
||||||
MAINOBJ=gray.o
|
MAINOBJ=gray.o
|
||||||
OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eierf.o \
|
OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eccd.o eierf.o \
|
||||||
graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
|
graydata_anequil.o graydata_flags.o graydata_par.o \
|
||||||
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
||||||
reflections.o simplespline.o utils.o beamdata.o
|
reflections.o simplespline.o utils.o beamdata.o
|
||||||
|
|
||||||
@ -25,17 +25,18 @@ $(EXE): $(MAINOBJ) $(OTHOBJ)
|
|||||||
$(FC) $(FFLAGS) -o $@ $^
|
$(FC) $(FFLAGS) -o $@ $^
|
||||||
|
|
||||||
# Dependencies on modules
|
# Dependencies on modules
|
||||||
gray.o: const_and_precisions.o conical.o dierckx.o dispersion.o \
|
gray.o: const_and_precisions.o dierckx.o dispersion.o eccd.o \
|
||||||
graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
|
graydata_anequil.o graydata_flags.o graydata_par.o \
|
||||||
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
|
||||||
reflections.o simplespline.o utils.o beamdata.o
|
reflections.o simplespline.o utils.o beamdata.o
|
||||||
conical.o: const_and_precisions.o
|
conical.o: const_and_precisions.o
|
||||||
dierckx.o: const_and_precisions.o
|
dierckx.o: const_and_precisions.o
|
||||||
dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
|
dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
|
||||||
|
eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.o
|
||||||
|
eierf.o: const_and_precisions.o
|
||||||
graydata_anequil.o: const_and_precisions.o
|
graydata_anequil.o: const_and_precisions.o
|
||||||
graydata_flags.o: const_and_precisions.o
|
graydata_flags.o: const_and_precisions.o
|
||||||
graydata_par.o: const_and_precisions.o
|
graydata_par.o: const_and_precisions.o
|
||||||
green_func_p.o: const_and_precisions.o numint.o
|
|
||||||
interp_eqprof.o: const_and_precisions.o
|
interp_eqprof.o: const_and_precisions.o
|
||||||
magsurf_data.o: const_and_precisions.o
|
magsurf_data.o: const_and_precisions.o
|
||||||
math.o: const_and_precisions.o
|
math.o: const_and_precisions.o
|
||||||
|
|||||||
889
src/eccd.f90
Normal file
889
src/eccd.f90
Normal file
@ -0,0 +1,889 @@
|
|||||||
|
module eccd
|
||||||
|
use const_and_precisions, only : wp_
|
||||||
|
implicit none
|
||||||
|
real(wp_), parameter, private :: cst2min=1.0e-6_wp_ ! min width of trap. cone
|
||||||
|
integer, parameter, private :: nfpp=13, & ! number of extra parameters passed
|
||||||
|
nfpp1=nfpp+ 1, nfpp2=nfpp+ 2, & ! to the integrand function fpp
|
||||||
|
nfpp3=nfpp+ 3, nfpp4=nfpp+ 4, &
|
||||||
|
nfpp5=nfpp+ 5
|
||||||
|
!########################################################################
|
||||||
|
! the following parameters are used by N.M. subroutines:
|
||||||
|
! The module contains few subroutines which are requested to calculate
|
||||||
|
! the current drive value by adjoint approach
|
||||||
|
!########################################################################
|
||||||
|
CHARACTER, PRIVATE, PARAMETER :: adj_appr(1:6) = & ! adj. approach switcher
|
||||||
|
(/ 'l', & ! (1)='l': collisionless limit
|
||||||
|
! (1)='c': collisional (classical) limit,
|
||||||
|
! w/o trap. part.
|
||||||
|
'm', & ! (2)='m': momentum conservation
|
||||||
|
! (2)='h': high-speed limit
|
||||||
|
!---
|
||||||
|
'l', & ! DO NOT CHANGE!
|
||||||
|
'r', & ! DO NOT CHANGE!
|
||||||
|
'v', & ! DO NOT CHANGE!
|
||||||
|
'i' /) ! DO NOT CHANGE!
|
||||||
|
!-------
|
||||||
|
REAL(wp_), PRIVATE :: r2,q2,gp1 ! coefficients for HSL integrand function
|
||||||
|
!-------
|
||||||
|
REAL(wp_), PRIVATE, PARAMETER :: delta = 1e-4 ! border for recalculation
|
||||||
|
!------- for N.M. subroutines (variational principle) -------
|
||||||
|
REAL(wp_), PRIVATE :: sfd(1:4) ! polyn. exp. of the "Spitzer"-function
|
||||||
|
INTEGER, PRIVATE, PARAMETER :: nre = 2 ! order of rel. correct.
|
||||||
|
REAL(wp_), PRIVATE, PARAMETER :: vp_mee(0:4,0:4,0:2) = &
|
||||||
|
RESHAPE((/0.0, 0.0, 0.0, 0.0, 0.0, &
|
||||||
|
0.0, 0.184875, 0.484304, 1.06069, 2.26175, &
|
||||||
|
0.0, 0.484304, 1.41421, 3.38514, 7.77817, &
|
||||||
|
0.0, 1.06069, 3.38514, 8.73232, 21.4005, &
|
||||||
|
0.0, 2.26175, 7.77817, 21.4005, 55.5079, &
|
||||||
|
! &
|
||||||
|
0.0, -1.33059,-2.57431, -5.07771, -10.3884, &
|
||||||
|
-0.846284,-1.46337, -1.4941, -0.799288, 2.57505, &
|
||||||
|
-1.1601, -1.4941, 2.25114, 14.159, 50.0534, &
|
||||||
|
-1.69257, -0.799288, 14.159, 61.4168, 204.389, &
|
||||||
|
-2.61022, 2.57505, 50.0534, 204.389, 683.756, &
|
||||||
|
! &
|
||||||
|
0.0, 2.62498, 0.985392,-5.57449, -27.683, &
|
||||||
|
0.0, 3.45785, 5.10096, 9.34463, 22.9831, &
|
||||||
|
-0.652555, 5.10096, 20.5135, 75.8022, 268.944, &
|
||||||
|
-2.11571, 9.34463, 75.8022, 330.42, 1248.69, &
|
||||||
|
-5.38358, 22.9831, 268.944, 1248.69, 4876.48/),&
|
||||||
|
(/5,5,3/))
|
||||||
|
REAL(wp_), PRIVATE, PARAMETER :: vp_mei(0:4,0:4,0:2) = &
|
||||||
|
RESHAPE((/0.0, 0.886227, 1.0, 1.32934, 2.0, &
|
||||||
|
0.886227,1.0, 1.32934, 2.0, 3.32335, &
|
||||||
|
1.0, 1.32934, 2.0, 3.32335, 6.0, &
|
||||||
|
1.32934, 2.0, 3.32335, 6.0, 11.6317, &
|
||||||
|
2.0, 3.32335, 6.0, 11.6317, 24.0, &
|
||||||
|
! &
|
||||||
|
0.0, 0.332335, 1.0, 2.49251, 6.0, &
|
||||||
|
1.66168, 1.0, 2.49251, 6.0, 14.5397, &
|
||||||
|
3.0, 2.49251, 6.0, 14.5397, 36.0, &
|
||||||
|
5.81586, 6.0, 14.5397, 36.0, 91.5999, &
|
||||||
|
12.0, 14.5397, 36.0, 91.5999, 240.0, &
|
||||||
|
! &
|
||||||
|
0.0, -0.103855, 0.0, 1.09047, 6.0, &
|
||||||
|
0.726983,0.0, 1.09047, 6.0, 24.5357, &
|
||||||
|
3.0, 1.09047, 6.0, 24.5357, 90.0, &
|
||||||
|
9.81427, 6.0, 24.5357, 90.0, 314.875, &
|
||||||
|
30.0, 24.5357, 90.0, 314.875, 1080.0 /), &
|
||||||
|
(/5,5,3/))
|
||||||
|
REAL(wp_), PRIVATE, PARAMETER :: vp_oee(0:4,0:4,0:2) = &
|
||||||
|
RESHAPE((/0.0, 0.56419, 0.707107, 1.0073, 1.59099, &
|
||||||
|
0.56419, 0.707107, 1.0073, 1.59099, 2.73981, &
|
||||||
|
0.707107,1.0073, 1.59099, 2.73981, 5.08233, &
|
||||||
|
1.0073, 1.59099, 2.73981, 5.08233, 10.0627, &
|
||||||
|
1.59099, 2.73981, 5.08233, 10.0627, 21.1138, &
|
||||||
|
! &
|
||||||
|
0.0, 1.16832, 1.90035, 3.5758, 7.41357, &
|
||||||
|
2.17562, 1.90035, 3.5758, 7.41357, 16.4891, &
|
||||||
|
3.49134, 3.5758, 7.41357, 16.4891, 38.7611, &
|
||||||
|
6.31562, 7.41357, 16.4891, 38.7611, 95.4472, &
|
||||||
|
12.4959, 16.4891, 38.7611, 95.4472, 244.803, &
|
||||||
|
! &
|
||||||
|
0.0, 2.65931, 4.64177, 9.6032, 22.6941, &
|
||||||
|
4.8652, 4.64177, 9.6032, 22.6941, 59.1437, &
|
||||||
|
9.51418, 9.6032, 22.6941, 59.1437, 165.282, &
|
||||||
|
21.061, 22.6941, 59.1437, 165.282, 485.785, &
|
||||||
|
50.8982, 59.1437, 165.282, 485.785, 1483.22/), &
|
||||||
|
(/5,5,3/))
|
||||||
|
REAL(wp_), PRIVATE, PARAMETER :: vp_g(0:4,0:2) = &
|
||||||
|
RESHAPE((/1.32934, 2.0, 3.32335, 6.0, 11.6317, &
|
||||||
|
2.49251, 0.0, 2.90793, 12.0, 39.2571, &
|
||||||
|
1.09047, 6.0, 11.45, 30.0, 98.9606/), &
|
||||||
|
(/5,3/))
|
||||||
|
!########################################################################
|
||||||
|
|
||||||
|
interface setcdcoeff
|
||||||
|
module procedure setcdcoeff_notrap,setcdcoeff_cohen,setcdcoeff_ncl
|
||||||
|
end interface setcdcoeff
|
||||||
|
|
||||||
|
contains
|
||||||
|
|
||||||
|
subroutine initeccd(ieccd)
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: ieccd
|
||||||
|
end subroutine initeccd
|
||||||
|
|
||||||
|
subroutine setcdcoeff_notrap(zeff,cst2,eccdpar)
|
||||||
|
implicit none
|
||||||
|
real(wp_), intent(in) :: zeff
|
||||||
|
real(wp_), intent(out) :: cst2
|
||||||
|
real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
|
||||||
|
|
||||||
|
cst2=0.0_wp_
|
||||||
|
allocate(eccdpar(1))
|
||||||
|
eccdpar(1)=zeff
|
||||||
|
end subroutine setcdcoeff_notrap
|
||||||
|
|
||||||
|
subroutine setcdcoeff_cohen(zeff,rbn,rbx,cst2,eccdpar)
|
||||||
|
! cohen model
|
||||||
|
! rbn=B/B_min
|
||||||
|
! rbx=B/B_max
|
||||||
|
! cst2=1.0_wp_-B/B_max
|
||||||
|
! alams=sqrt(1-B_min/B_max)
|
||||||
|
! Zeff < 31 !!!
|
||||||
|
! fp0s= P_a (alams)
|
||||||
|
use conical, only : fconic
|
||||||
|
implicit none
|
||||||
|
real(wp_), intent(in) :: zeff,rbn,rbx
|
||||||
|
real(wp_), intent(out) :: cst2
|
||||||
|
real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
|
||||||
|
real(wp_) :: alams,pa,fp0s
|
||||||
|
|
||||||
|
cst2=1.0_wp_-rbx
|
||||||
|
if(cst2<cst2min) cst2=0.0_wp_
|
||||||
|
alams=sqrt(1.0_wp_-rbx/rbn)
|
||||||
|
pa=sqrt(32.0_wp_/(Zeff+1.0_wp_)-1.0_wp_)/2.0_wp_
|
||||||
|
fp0s=fconic(alams,pa,0)
|
||||||
|
allocate(eccdpar(5))
|
||||||
|
eccdpar(1)=zeff
|
||||||
|
eccdpar(2)=rbn
|
||||||
|
eccdpar(3)=alams
|
||||||
|
eccdpar(4)=pa
|
||||||
|
eccdpar(5)=fp0s
|
||||||
|
end subroutine setcdcoeff_cohen
|
||||||
|
|
||||||
|
subroutine setcdcoeff_ncl(zeff,rbx,fc,amu,rhop,cst2,eccdpar)
|
||||||
|
use magsurf_data, only : ch,tjp,tlm,njpt,nlmt
|
||||||
|
use dierckx, only : profil
|
||||||
|
implicit none
|
||||||
|
integer, parameter :: ksp=3
|
||||||
|
real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop
|
||||||
|
real(wp_), intent(out) :: cst2
|
||||||
|
real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
|
||||||
|
real(wp_) :: alams,pa,fp0s
|
||||||
|
real(wp_), dimension(nlmt) :: chlm
|
||||||
|
integer :: nlm,ierr,npar
|
||||||
|
|
||||||
|
cst2=1.0_wp_-rbx
|
||||||
|
if(cst2<cst2min) cst2=0.0_wp_
|
||||||
|
call SpitzFuncCoeff(amu,Zeff,fc)
|
||||||
|
nlm=nlmt
|
||||||
|
call profil(0,tjp,njpt,tlm,nlmt,ch,ksp,ksp,rhop,nlm,chlm,ierr)
|
||||||
|
if(ierr>0) write(*,*) ' Hlambda profil =',ierr
|
||||||
|
npar=3+2*nlm
|
||||||
|
allocate(eccdpar(npar))
|
||||||
|
eccdpar(1)=zeff
|
||||||
|
eccdpar(2) = fc
|
||||||
|
eccdpar(3) = rbx
|
||||||
|
eccdpar(4:3+nlm) = tlm
|
||||||
|
eccdpar(4+nlm:npar) = chlm
|
||||||
|
end subroutine setcdcoeff_ncl
|
||||||
|
|
||||||
|
subroutine eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmn,nhmx, &
|
||||||
|
ithn,cst2,fcur,eccdpar,effjcd,iokhawa,ierr)
|
||||||
|
use const_and_precisions, only : pi,qesi=>e_,mesi=>me_, &
|
||||||
|
vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_
|
||||||
|
use quadpack, only : dqagsmv
|
||||||
|
implicit none
|
||||||
|
! local constants
|
||||||
|
real(wp_), parameter :: mc2m2=1.0_wp_/mc2**2, &
|
||||||
|
canucc=2.0e13_wp_*pi*qe**4/(me**2*vc**3),ceff=qesi/(mesi*vcsi)
|
||||||
|
real(wp_), parameter :: epsa=0.0_wp_,epsr=1.0e-2_wp_,xxcr=16.0_wp_
|
||||||
|
real(wp_), parameter :: dumin=1.0e-6_wp_
|
||||||
|
integer, parameter :: lw=5000,liw=lw/4
|
||||||
|
! arguments
|
||||||
|
integer :: i,nhmn,nhmx,ithn,iokhawa,ierr
|
||||||
|
real(wp_) :: yg,anpl,anprre,dens,amu,cst2,effjcd
|
||||||
|
real(wp_), dimension(:) :: eccdpar
|
||||||
|
complex(wp_) :: ex,ey,ez
|
||||||
|
! local variables
|
||||||
|
integer :: nhn,neval,ier,last,npar
|
||||||
|
integer, dimension(liw) :: iw
|
||||||
|
real(wp_) :: anpl2,dnl,ygn,ygn2,resji,rdu2,upltp,upltm,uplp,uplm, &
|
||||||
|
rdu,rdu2t,duu,uu1,uu2,xx1,xx2,resj,resp,eji,epp,anum,denom, &
|
||||||
|
cstrdut,anucc
|
||||||
|
real(wp_), dimension(lw) :: w
|
||||||
|
real(wp_), dimension(nfpp+size(eccdpar)) :: apar
|
||||||
|
real(wp_), dimension(0:1) :: uleft,uright
|
||||||
|
! common/external functions/variables
|
||||||
|
real(wp_), external :: fcur
|
||||||
|
!
|
||||||
|
! effjpl = <J_parallel>/<p_d> /(B_min/<B>) [A m /W]
|
||||||
|
!
|
||||||
|
apar(1) = yg
|
||||||
|
apar(2) = anpl
|
||||||
|
apar(3) = amu
|
||||||
|
apar(4) = anprre
|
||||||
|
apar(5) = dble(ex)
|
||||||
|
apar(6) = dimag(ex)
|
||||||
|
apar(7) = dble(ey)
|
||||||
|
apar(8) = dimag(ey)
|
||||||
|
apar(9) = dble(ez)
|
||||||
|
apar(10) = dimag(ez)
|
||||||
|
apar(11) = dble(ithn)
|
||||||
|
|
||||||
|
npar=size(apar)
|
||||||
|
apar(nfpp+1:npar) = eccdpar
|
||||||
|
|
||||||
|
anpl2=anpl*anpl
|
||||||
|
|
||||||
|
effjcd=0.0_wp_
|
||||||
|
anum=0.0_wp_
|
||||||
|
denom=0.0_wp_
|
||||||
|
iokhawa=0
|
||||||
|
ierr=0
|
||||||
|
do nhn=nhmn,nhmx
|
||||||
|
ygn=nhn*yg
|
||||||
|
ygn2=ygn*ygn
|
||||||
|
|
||||||
|
rdu2=anpl2+ygn2-1.0_wp_
|
||||||
|
|
||||||
|
if (rdu2.lt.0.0_wp_) cycle
|
||||||
|
rdu=sqrt(rdu2)
|
||||||
|
dnl=1.0_wp_-anpl2
|
||||||
|
uplp=(anpl*ygn+rdu)/dnl
|
||||||
|
uplm=(anpl*ygn-rdu)/dnl
|
||||||
|
|
||||||
|
uu1=uplm
|
||||||
|
uu2=uplp
|
||||||
|
xx1=amu*(anpl*uu1+ygn-1.0_wp_)
|
||||||
|
xx2=amu*(anpl*uu2+ygn-1.0_wp_)
|
||||||
|
|
||||||
|
if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0_wp_)/anpl
|
||||||
|
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl
|
||||||
|
duu=abs(uu1-uu2)
|
||||||
|
|
||||||
|
if(duu.le.dumin) cycle
|
||||||
|
|
||||||
|
apar(12) = dble(nhn)
|
||||||
|
apar(13) = ygn
|
||||||
|
|
||||||
|
call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp, &
|
||||||
|
epp,neval,ier,liw,lw,last,iw,w)
|
||||||
|
if (ier.gt.0) then
|
||||||
|
ierr=90
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
|
||||||
|
rdu2t=cst2*anpl2+ygn2-1.0_wp_
|
||||||
|
|
||||||
|
if (rdu2t.gt.0.0_wp_.and.cst2.gt.0.0_wp_) then
|
||||||
|
!
|
||||||
|
! resonance curve crosses the trapping region
|
||||||
|
!
|
||||||
|
iokhawa=1
|
||||||
|
cstrdut=sqrt(cst2*rdu2t)
|
||||||
|
upltm=(cst2*anpl*ygn-cstrdut)/(1.0_wp_-cst2*anpl2)
|
||||||
|
upltp=(cst2*anpl*ygn+cstrdut)/(1.0_wp_-cst2*anpl2)
|
||||||
|
uleft(0)=uplm
|
||||||
|
uright(0)=upltm
|
||||||
|
uleft(1)=upltp
|
||||||
|
uright(1)=uplp
|
||||||
|
else
|
||||||
|
!
|
||||||
|
! resonance curve does not cross the trapping region
|
||||||
|
!
|
||||||
|
iokhawa=0
|
||||||
|
uleft(0)=uplm
|
||||||
|
uright(0)=uplp
|
||||||
|
end if
|
||||||
|
|
||||||
|
resj=0.0_wp_
|
||||||
|
do i=0,1
|
||||||
|
xx1=amu*(anpl*uleft(i)+ygn-1.0_wp_)
|
||||||
|
xx2=amu*(anpl*uright(i)+ygn-1.0_wp_)
|
||||||
|
if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
|
||||||
|
if(xx2.gt.xxcr) uright(i)=(xxcr/amu-ygn+1.0_wp_)/anpl
|
||||||
|
if(xx1.gt.xxcr) uleft(i)=(xxcr/amu-ygn+1.0_wp_)/anpl
|
||||||
|
duu=abs(uleft(i)-uright(i))
|
||||||
|
if(duu.gt.dumin) then
|
||||||
|
call dqagsmv(fcur,uleft(i),uright(i),apar,npar,epsa,epsr, &
|
||||||
|
resji,eji,neval,ier,liw,lw,last,iw,w)
|
||||||
|
if (ier.gt.0) then
|
||||||
|
if (abs(resji).lt.1.0e-10_wp_) then
|
||||||
|
resji=0.0_wp_
|
||||||
|
else
|
||||||
|
ierr=91+iokhawa+i
|
||||||
|
return
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
resj=resj+resji
|
||||||
|
if(iokhawa.eq.0) exit
|
||||||
|
end do
|
||||||
|
anum=anum+resj
|
||||||
|
denom=denom+resp
|
||||||
|
end do
|
||||||
|
|
||||||
|
if(denom.gt.0.0_wp_) then
|
||||||
|
anucc=canucc*dens*(48.0_wp_-log(1.0e7_wp_*dens*mc2m2*amu**2))
|
||||||
|
effjcd=-ceff*anum/(anucc*denom)
|
||||||
|
end if
|
||||||
|
|
||||||
|
end subroutine eccdeff
|
||||||
|
|
||||||
|
function fpp(upl,extrapar,npar)
|
||||||
|
!
|
||||||
|
! computation of integral for power density, integrand function fpp
|
||||||
|
!
|
||||||
|
! ith=0 : polarization term = const
|
||||||
|
! ith=1 : polarization term Larmor radius expansion to lowest order
|
||||||
|
! ith=2 : full polarization term (J Bessel)
|
||||||
|
!
|
||||||
|
! integration variable upl passed explicitly. other variables passed
|
||||||
|
! as array of extra parameters of length npar=size(extrapar)
|
||||||
|
!
|
||||||
|
! extrapar(1) = yg
|
||||||
|
! extrapar(2) = anpl
|
||||||
|
! extrapar(3) = amu
|
||||||
|
! extrapar(4) = Re(anprw)
|
||||||
|
! extrapar(5) = Re(ex)
|
||||||
|
! extrapar(6) = Im(ex)
|
||||||
|
! extrapar(7) = Re(ey)
|
||||||
|
! extrapar(8) = Im(ey)
|
||||||
|
! extrapar(9) = Re(ez)
|
||||||
|
! extrapar(10) = Im(ez)
|
||||||
|
! extrapar(11) = double(ithn)
|
||||||
|
! extrapar(12) = double(nhn)
|
||||||
|
! extrapar(13) = ygn
|
||||||
|
!
|
||||||
|
use const_and_precisions, only : ui=>im
|
||||||
|
use math, only : fact
|
||||||
|
implicit none
|
||||||
|
! arguments
|
||||||
|
integer :: npar
|
||||||
|
real(wp_) :: upl,fpp
|
||||||
|
real(wp_), dimension(npar) :: extrapar
|
||||||
|
! local variables
|
||||||
|
integer :: ithn,nhn,nm,np
|
||||||
|
real(wp_) :: yg,anpl,amu,anprre,ygn,upr,upr2,gam,ee,thn2,thn2u,bb,cth, &
|
||||||
|
ajbnm,ajbnp,ajbn
|
||||||
|
complex(wp_) :: ex,ey,ez,emxy,epxy
|
||||||
|
|
||||||
|
yg=extrapar(1)
|
||||||
|
anpl=extrapar(2)
|
||||||
|
amu=extrapar(3)
|
||||||
|
anprre=extrapar(4)
|
||||||
|
ex=cmplx(extrapar(5),extrapar(6),wp_)
|
||||||
|
ey=cmplx(extrapar(7),extrapar(8),wp_)
|
||||||
|
ez=cmplx(extrapar(9),extrapar(10),wp_)
|
||||||
|
ithn=int(extrapar(11))
|
||||||
|
nhn=int(extrapar(12))
|
||||||
|
ygn=extrapar(13)
|
||||||
|
|
||||||
|
gam=anpl*upl+ygn
|
||||||
|
upr2=gam*gam-1.0_wp_-upl*upl
|
||||||
|
ee=exp(-amu*(gam-1))
|
||||||
|
|
||||||
|
! thn2=1.0_wp_
|
||||||
|
thn2u=upr2 !*thn2
|
||||||
|
if(ithn.gt.0) then
|
||||||
|
emxy=ex-ui*ey
|
||||||
|
epxy=ex+ui*ey
|
||||||
|
if(upr2.gt.0.0_wp_) then
|
||||||
|
upr=sqrt(upr2)
|
||||||
|
bb=anprre*upr/yg
|
||||||
|
if(ithn.eq.1) then
|
||||||
|
! Larmor radius expansion polarization term at lowest order
|
||||||
|
cth=1.0_wp_
|
||||||
|
if(nhn.gt.1) cth=(0.5_wp_*bb)**(nhn-1)*nhn/fact(nhn)
|
||||||
|
thn2=(0.5_wp_*cth*abs(emxy+ez*anprre*upl/ygn))**2
|
||||||
|
thn2u=upr2*thn2
|
||||||
|
else
|
||||||
|
! Full polarization term
|
||||||
|
nm=nhn-1
|
||||||
|
np=nhn+1
|
||||||
|
ajbnm=dbesjn(nm, bb)
|
||||||
|
ajbnp=dbesjn(np, bb)
|
||||||
|
ajbn=dbesjn(nhn, bb)
|
||||||
|
thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0_wp_))**2
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
|
||||||
|
fpp=ee*thn2u
|
||||||
|
end function fpp
|
||||||
|
|
||||||
|
function fjch(upl,extrapar,npar)
|
||||||
|
!
|
||||||
|
! computation of integral for current density
|
||||||
|
! integrand for Cohen model with trapping
|
||||||
|
!
|
||||||
|
! integration variable upl passed explicitly. Other variables passed
|
||||||
|
! as array of extra parameters of length npar=size(extrapar).
|
||||||
|
! variables with index 1..nfpp must be passed to fpp
|
||||||
|
! variable with index nfpp+1 is zeff
|
||||||
|
! variables with index gt nfpp+1 are specific of the cd model
|
||||||
|
!
|
||||||
|
! extrapar(2) = anpl
|
||||||
|
! extrapar(4) = Re(anprw)
|
||||||
|
! extrapar(13) = ygn
|
||||||
|
!
|
||||||
|
! extrapar(14) = zeff
|
||||||
|
! extrapar(15) = rb
|
||||||
|
! extrapar(16) = alams
|
||||||
|
! extrapar(17) = pa
|
||||||
|
! extrapar(18) = fp0s
|
||||||
|
!
|
||||||
|
use conical, only : fconic
|
||||||
|
implicit none
|
||||||
|
! arguments
|
||||||
|
integer :: npar
|
||||||
|
real(wp_) :: upl,fjch
|
||||||
|
real(wp_), dimension(npar) :: extrapar
|
||||||
|
! local variables
|
||||||
|
integer :: nhn
|
||||||
|
real(wp_) :: anpl,anprre,ygn,zeff,rb,alams,pa,fp0s, &
|
||||||
|
upr2,gam,u2,u,z5,xi,xib,xibi,fu2b,fu2,gu,gg,dgg,alam,fp0, &
|
||||||
|
dfp0,fh,dfhl,eta
|
||||||
|
|
||||||
|
anpl=extrapar(2)
|
||||||
|
anprre=extrapar(4)
|
||||||
|
ygn=extrapar(13)
|
||||||
|
zeff=extrapar(nfpp1)
|
||||||
|
rb=extrapar(nfpp2)
|
||||||
|
alams=extrapar(nfpp3)
|
||||||
|
pa=extrapar(nfpp4)
|
||||||
|
fp0s=extrapar(nfpp5)
|
||||||
|
|
||||||
|
gam=anpl*upl+ygn
|
||||||
|
u2=gam*gam-1.0_wp_
|
||||||
|
upr2=u2-upl*upl
|
||||||
|
u=sqrt(u2)
|
||||||
|
z5=Zeff+5.0_wp_
|
||||||
|
xi=1.0_wp_/z5**2
|
||||||
|
xib=1.0_wp_-xi
|
||||||
|
xibi=1.0_wp_/xib
|
||||||
|
fu2b=1.0_wp_+xib*u2
|
||||||
|
fu2=1.0_wp_+xi*u2
|
||||||
|
gu=(1.0_wp_-1.0_wp_/fu2b**xibi)/sqrt(fu2)
|
||||||
|
gg=u*gu/z5
|
||||||
|
dgg=(gu+u2*(2.0_wp_/fu2b**(1.0_wp_+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
|
||||||
|
|
||||||
|
alam=sqrt(1.0_wp_-upr2/u2/rb)
|
||||||
|
fp0=fconic(alam,pa,0)
|
||||||
|
dfp0=-(pa*pa/2.0_wp_+0.125_wp_)
|
||||||
|
if (alam.lt.1.0_wp_) dfp0=-fconic(alam,pa,1)/sqrt(1.0_wp_-alam**2)
|
||||||
|
fh=alam*(1.0_wp_-alams*fp0/(alam*fp0s))
|
||||||
|
dfhl=1.0_wp_-alams*dfp0/fp0s
|
||||||
|
|
||||||
|
eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam)
|
||||||
|
|
||||||
|
if(upl.lt.0.0_wp_) eta=-eta
|
||||||
|
fjch=eta*fpp(upl,extrapar(1:nfpp),nfpp)
|
||||||
|
|
||||||
|
end function fjch
|
||||||
|
|
||||||
|
function fjch0(upl,extrapar,npar)
|
||||||
|
!
|
||||||
|
! computation of integral for current density
|
||||||
|
! integrand for Cohen model without trapping
|
||||||
|
!
|
||||||
|
! integration variable upl passed explicitly. Other variables passed
|
||||||
|
! as array of extra parameters of length npar=size(extrapar).
|
||||||
|
! variables with index 1..nfpp must be passed to fpp
|
||||||
|
! variable with index nfpp+1 is zeff
|
||||||
|
! variables with index gt nfpp+1 are specific of the cd model
|
||||||
|
!
|
||||||
|
! extrapar(2) = anpl
|
||||||
|
! extrapar(13) = ygn
|
||||||
|
!
|
||||||
|
! extrapar(14) = zeff
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
! arguments
|
||||||
|
real(wp_) :: upl,fjch0
|
||||||
|
integer :: npar
|
||||||
|
real(wp_), dimension(npar) :: extrapar
|
||||||
|
! local variables
|
||||||
|
real(wp_) :: anpl,ygn,zeff,gam,u2,u,z5,xi,xib,xibi,fu2b,fu2,gu,gg,dgg,eta
|
||||||
|
!
|
||||||
|
anpl=extrapar(2)
|
||||||
|
ygn=extrapar(13)
|
||||||
|
zeff=extrapar(nfpp1)
|
||||||
|
|
||||||
|
gam=anpl*upl+ygn
|
||||||
|
u2=gam*gam-1.0_wp_
|
||||||
|
u=sqrt(u2)
|
||||||
|
z5=Zeff+5.0_wp_
|
||||||
|
xi=1.0_wp_/z5**2
|
||||||
|
xib=1.0_wp_-xi
|
||||||
|
xibi=1.0_wp_/xib
|
||||||
|
fu2b=1.0_wp_+xib*u2
|
||||||
|
fu2=1.0_wp_+xi*u2
|
||||||
|
gu=(1.0_wp_-1.0_wp_/fu2b**xibi)/sqrt(fu2)
|
||||||
|
gg=u*gu/z5
|
||||||
|
dgg=(gu+u2*(2.0_wp_/fu2b**(1.0_wp_+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
|
||||||
|
eta=anpl*gg+gam*upl*dgg/u
|
||||||
|
fjch0=eta*fpp(upl,extrapar(1:nfpp),nfpp)
|
||||||
|
|
||||||
|
end function fjch0
|
||||||
|
|
||||||
|
function fjncl(upl,extrapar,npar)
|
||||||
|
!
|
||||||
|
! computation of integral for current density
|
||||||
|
! integrand for momentum conserv. model K(u) from Maruschenko
|
||||||
|
! gg=F(u)/u with F(u) as in Cohen paper
|
||||||
|
!
|
||||||
|
! integration variable upl passed explicitly. Other variables passed
|
||||||
|
! as array of extra parameters of length npar=size(extrapar).
|
||||||
|
! variables with index 1..nfpp must be passed to fpp
|
||||||
|
! variable with index nfpp+1 is zeff
|
||||||
|
! variables with index gt nfpp+1 are specific of the cd model
|
||||||
|
!
|
||||||
|
! extrapar(2) = anpl
|
||||||
|
! extrapar(3) = amu
|
||||||
|
! extrapar(13) = ygn
|
||||||
|
!
|
||||||
|
! extrapar(14) = zeff
|
||||||
|
! extrapar(15) = fc
|
||||||
|
! extrapar(16) = rbx
|
||||||
|
! extrapar(17:16+(npar-16)/2) = tlm
|
||||||
|
! extrapar(17+(npar-16)/2:npar) = chlm
|
||||||
|
!
|
||||||
|
use dierckx, only : splev,splder
|
||||||
|
implicit none
|
||||||
|
! arguments
|
||||||
|
integer :: npar
|
||||||
|
real(wp_) :: upl,fjncl
|
||||||
|
real(wp_), dimension(npar) :: extrapar
|
||||||
|
! local variables
|
||||||
|
integer :: nlm
|
||||||
|
real(wp_) :: anpl,amu,ygn,zeff,fc,rbx,gam,u2,u,upr2, &
|
||||||
|
bth,uth,fk,dfk,alam,fu,dfu,eta
|
||||||
|
! local variables
|
||||||
|
integer :: ier
|
||||||
|
real(wp_), dimension((npar-nfpp3)/2) :: wrk
|
||||||
|
real(wp_), dimension(1) :: xs,ys
|
||||||
|
!
|
||||||
|
anpl=extrapar(2)
|
||||||
|
amu=extrapar(3)
|
||||||
|
ygn=extrapar(13)
|
||||||
|
zeff=extrapar(nfpp1)
|
||||||
|
fc=extrapar(nfpp2)
|
||||||
|
rbx=extrapar(nfpp3)
|
||||||
|
|
||||||
|
gam=anpl*upl+ygn
|
||||||
|
u2=gam*gam-1.0_wp_
|
||||||
|
u=sqrt(u2)
|
||||||
|
upr2=u2-upl*upl
|
||||||
|
bth=sqrt(2.0_wp_/amu)
|
||||||
|
uth=u/bth
|
||||||
|
call GenSpitzFunc(Zeff,fc,uth,u,gam,fk,dfk)
|
||||||
|
fk=fk*(4.0_wp_/amu**2)
|
||||||
|
dfk=dfk*(2.0_wp_/amu)*bth
|
||||||
|
|
||||||
|
alam=upr2/u2/rbx
|
||||||
|
xs(1)=alam
|
||||||
|
nlm=(npar-nfpp3)/2
|
||||||
|
!
|
||||||
|
! extrapar(17:16+(npar-16)/2) = tlm
|
||||||
|
! extrapar(17+(npar-16)/2:npar) = chlm
|
||||||
|
!
|
||||||
|
call splev(extrapar(nfpp4:nfpp3+nlm),nlm,extrapar(nfpp4+nlm:npar),3, &
|
||||||
|
xs(1),ys(1),1,ier)
|
||||||
|
fu=ys(1)
|
||||||
|
call splder(extrapar(nfpp4:nfpp3+nlm),nlm,extrapar(nfpp4+nlm:npar),3,1, &
|
||||||
|
xs(1),ys(1),1,wrk,ier)
|
||||||
|
dfu=ys(1)
|
||||||
|
|
||||||
|
eta=gam*fu*dfk/u-2.0_wp_*(anpl-gam*upl/u2)*fk*dfu*upl/u2/rbx
|
||||||
|
if(upl.lt.0) eta=-eta
|
||||||
|
fjncl=eta*fpp(upl,extrapar(1:nfpp),nfpp)
|
||||||
|
end function fjncl
|
||||||
|
|
||||||
|
SUBROUTINE GenSpitzFunc(Zeff,fc,u,q,gam, K,dKdu)
|
||||||
|
!=======================================================================
|
||||||
|
! Author: N.B.Marushchenko
|
||||||
|
! June 2005: as start point the subroutine of Ugo Gasparino (198?)
|
||||||
|
! SpitzFunc() is taken and modified.
|
||||||
|
! 1. adapted to the Fortran-95
|
||||||
|
! 2. derivative of Spitzer function is added
|
||||||
|
! 3. separation for 2 brunches is done:
|
||||||
|
! 1st is referenced as 'with conservation of the moment',
|
||||||
|
! 2nd - as 'high speed limit'.
|
||||||
|
! The last one is taken from the Lin-Liu formulation
|
||||||
|
! (Phys.Plasmas 10 (2003) 4064) with K = F*fc.
|
||||||
|
! The asymptotical high speed limit (Taguchi-Fisch model)
|
||||||
|
! is also included as the reference case.
|
||||||
|
! Feb. 2008: non-relativ. version is replaced by the relativistic one;
|
||||||
|
! the method is the the same, but the trial-function is
|
||||||
|
! based on the relativistic formulation.
|
||||||
|
! The relativistic corrections for the collisional operator
|
||||||
|
! up to the second order, i.e. (1/mu)**2, are applied.
|
||||||
|
! Sep. 2008: generalized Spitzer function for arbitrary collisionality
|
||||||
|
! is implemented. The model is based on the concept of
|
||||||
|
! the "effective trapped particles fraction".
|
||||||
|
! The different.-integral kinetic equation for the generalized
|
||||||
|
! Spitzer function is produced with help of subroutines
|
||||||
|
! ArbColl_TrappFract_Array and ArbColl_SpitzFunc_Array,
|
||||||
|
! where the subroutines of H. Maassberg are called).
|
||||||
|
!========================================================================
|
||||||
|
! Spitzer function with & w/o trapped particle effects is given by:
|
||||||
|
!
|
||||||
|
! K(x) = x/gamma*(d1*x+d2*x^2+d4*x^3+d4*x^4),
|
||||||
|
!
|
||||||
|
! where x = v/v_th and gamma=1 for non-relativistic version (Ugo),
|
||||||
|
! or x = p/p_th for relativistic version (N.M., February 2008).
|
||||||
|
! Note, that somewhere the function F(x) instead of K(x) is applied,
|
||||||
|
!
|
||||||
|
! F(x) = K(x)/fc.
|
||||||
|
!
|
||||||
|
! Numerical inversion of the 5x5 symmetric matrix obtained from the
|
||||||
|
! generalized Spitzer problem (see paper of Taguchi for the equation
|
||||||
|
! and paper of Hirshman for the variational approach bringing to the
|
||||||
|
! matrix to be inverted).
|
||||||
|
!
|
||||||
|
! The numerical method used is an improved elimination scheme
|
||||||
|
! (Banachiewiczs-Cholesky-Crout method).
|
||||||
|
! This method is particularly simple for symmetric matrix.
|
||||||
|
! As a reference see "Mathematical Handbook" by Korn & Korn, p.635-636.
|
||||||
|
!
|
||||||
|
! Refs.: 1. S.P. Hirshman, Phys. Fluids 23 (1980) 1238
|
||||||
|
! 2. M. Rome' et al., Plasma Phys. Contr. Fus. 40 (1998) 511
|
||||||
|
! 3. N.B. Marushchenko et al., Fusion Sci. Technol. 55 (2009) 180
|
||||||
|
!========================================================================
|
||||||
|
! INPUTS:
|
||||||
|
! u - p/sqrt(2mT)
|
||||||
|
! q - p/mc;
|
||||||
|
! gam - relativistic factor;
|
||||||
|
! Zeff - effective charge;
|
||||||
|
! fc - fraction of circulating particles.
|
||||||
|
!
|
||||||
|
! OUTPUTS:
|
||||||
|
! K - Spitzer's function
|
||||||
|
! dKdu = dK/du, i.e. its derivative over normalized momentum
|
||||||
|
!=======================================================================
|
||||||
|
use const_and_precisions, only : comp_eps
|
||||||
|
IMPLICIT NONE
|
||||||
|
REAL(wp_), INTENT(in) :: Zeff,fc,u,q,gam
|
||||||
|
REAL(wp_), INTENT(out) :: K,dKdu
|
||||||
|
REAL(wp_) :: gam1,gam2,gam3
|
||||||
|
|
||||||
|
K = 0
|
||||||
|
dKdu = 0
|
||||||
|
IF (u < comp_eps) RETURN
|
||||||
|
|
||||||
|
SELECT CASE(adj_appr(2))
|
||||||
|
CASE('m') !--------------- momentum conservation ------------------!
|
||||||
|
gam1 = gam !
|
||||||
|
IF (adj_appr(4) == 'n') gam1 = 1 !
|
||||||
|
gam2 = gam1*gam1 !
|
||||||
|
gam3 = gam1*gam2 !
|
||||||
|
K = u/gam1*u*(sfd(1)+u*(sfd(2)+u*(sfd(3)+u*sfd(4)))) !
|
||||||
|
dKdu = u/gam3* (sfd(1)*(1+ gam2)+u*(sfd(2)*(1+2*gam2)+ & !
|
||||||
|
u*(sfd(3)*(1+3*gam2)+u* sfd(4)*(1+4*gam2)))) !
|
||||||
|
!--------------------- end momentum conservation -------------------!
|
||||||
|
CASE('h') !---------------- high-speed-limit ----------------------!
|
||||||
|
IF (adj_appr(4) == 'n') THEN !- non-relativ. asymptotic form -!
|
||||||
|
K = u**4 *fc/(Zeff+1+4*fc) !- (Taguchi-Fisch model) -!
|
||||||
|
dKdu = 4*u**3 *fc/(Zeff+1+4*fc) !
|
||||||
|
ELSEIF (adj_appr(4) == 'r') THEN !- relativistic, Lin-Liu form. -!
|
||||||
|
CALL SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu) !
|
||||||
|
ENDIF !
|
||||||
|
CASE default !----------------------------------------------------!
|
||||||
|
PRINT*,'GenSpitzFunc: WARNING! Spitzer function is not defined.'
|
||||||
|
RETURN
|
||||||
|
END SELECT
|
||||||
|
END SUBROUTINE GenSpitzFunc
|
||||||
|
|
||||||
|
SUBROUTINE SpitzFuncCoeff(mu,Zeff,fc)
|
||||||
|
!=======================================================================
|
||||||
|
! Calculates the matrix coefficients required for the subroutine
|
||||||
|
! "GenSpitzFunc", where the Spitzer function is defined through the
|
||||||
|
! variational principle.
|
||||||
|
!
|
||||||
|
! Weakly relativistic (upgraded) version (10.09.2008).
|
||||||
|
! Apart of the non-relativistic matrix coefficients, taken from the
|
||||||
|
! old subroutine of Ugo Gasparino, the relativistic correction written
|
||||||
|
! as series in 1/mu^n (mu=mc2/T) powers is added. Two orders are taken
|
||||||
|
! into account, i.e. n=0,1,2.
|
||||||
|
!
|
||||||
|
! In this version, the coefficients "oee", i.e. Omega_ij, are formulated
|
||||||
|
! for arbitrary collisionality.
|
||||||
|
!
|
||||||
|
! INPUT VARIABLES:
|
||||||
|
! rho = sqrt(SS) with SS - flux-surface label (norm. magn. flux)
|
||||||
|
! ne - density, 1/m^3
|
||||||
|
! mu - mc2/Te
|
||||||
|
! Zeff - effective charge
|
||||||
|
! fc - fraction of circulating particles
|
||||||
|
!
|
||||||
|
! OUTPUT VARIABLES (defined as a global ones):
|
||||||
|
! sfd(1),...,sfd(4) - coefficients of the polynomial expansion of the
|
||||||
|
! "Spitzer"-function (the same as in the Hirshman paper)
|
||||||
|
!=======================================================================
|
||||||
|
use const_and_precisions, only : mc2_
|
||||||
|
IMPLICIT NONE
|
||||||
|
REAL(wp_), INTENT(in) :: mu,Zeff,fc
|
||||||
|
INTEGER :: n,i,j
|
||||||
|
REAL(wp_) :: rtc,rtc1,y,tn(1:nre)
|
||||||
|
REAL(wp_) :: m(0:4,0:4),g(0:4)
|
||||||
|
REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, &
|
||||||
|
gam22,gam32,gam42,gam02, &
|
||||||
|
gam33,gam43,gam03, &
|
||||||
|
gam44,gam04,gam00
|
||||||
|
REAL(wp_) :: alp12,alp13,alp14,alp10, &
|
||||||
|
alp23,alp24,alp20, &
|
||||||
|
alp34,alp30,alp40
|
||||||
|
REAL(wp_) :: bet0,bet1,bet2,bet3,bet4,d0
|
||||||
|
LOGICAL :: renew,rel,newmu,newZ,newfc
|
||||||
|
REAL(wp_), SAVE :: sfdx(1:4) = 0
|
||||||
|
REAL(wp_), SAVE :: mu_old =-1, Zeff_old =-1, fc_old =-1
|
||||||
|
|
||||||
|
rel = mu < mc2_
|
||||||
|
newmu = abs(mu -mu_old ) > delta*mu
|
||||||
|
newZ = abs(Zeff-Zeff_old) > delta*Zeff
|
||||||
|
newfc = abs(fc -fc_old ) > delta*fc
|
||||||
|
SELECT CASE(adj_appr(1))
|
||||||
|
CASE ('l','c')
|
||||||
|
renew = (newmu .and. rel) .OR. newZ .OR. newfc
|
||||||
|
END SELECT
|
||||||
|
IF (.not.renew) THEN
|
||||||
|
sfd(:) = sfdx(:)
|
||||||
|
RETURN
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
tn(:) = 0
|
||||||
|
IF (adj_appr(4) == 'r') THEN
|
||||||
|
IF (nre > 0) THEN
|
||||||
|
!mu = min(mu,1.e3*mc2_)
|
||||||
|
tn(1) = 1/mu
|
||||||
|
DO n=2,min(2,nre)
|
||||||
|
tn(n) = tn(n-1)/mu
|
||||||
|
ENDDO
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
SELECT CASE(adj_appr(1))
|
||||||
|
CASE ('l','c') !---- both classical & collisionless limits ----!
|
||||||
|
rtc = (1-fc)/fc; rtc1 = rtc+1 !
|
||||||
|
!--- !
|
||||||
|
DO i=0,4 !
|
||||||
|
g(i) = vp_g(i,0) !
|
||||||
|
DO n=1,min(2,nre) !
|
||||||
|
g(i) = g(i) + tn(n)*vp_g(i,n) !
|
||||||
|
ENDDO !
|
||||||
|
!--- !
|
||||||
|
DO j=0,4 !
|
||||||
|
IF (i == 0 .or. j == 0 .or. j >= i) THEN !
|
||||||
|
y = vp_mee(i,j,0) + rtc *vp_oee(i,j,0) + & !
|
||||||
|
Zeff*rtc1*vp_mei(i,j,0) !
|
||||||
|
DO n=1,min(2,nre) !
|
||||||
|
y = y + (vp_mee(i,j,n) + rtc *vp_oee(i,j,n) + & !
|
||||||
|
Zeff*rtc1*vp_mei(i,j,n))*tn(n) !
|
||||||
|
ENDDO !
|
||||||
|
m(i,j) = y !
|
||||||
|
ENDIF !
|
||||||
|
ENDDO !
|
||||||
|
ENDDO !
|
||||||
|
DO i=2,4 !
|
||||||
|
DO j=1,i-1 !
|
||||||
|
m(i,j) = m(j,i) !
|
||||||
|
ENDDO !
|
||||||
|
ENDDO !
|
||||||
|
m(0,0) = 0 !
|
||||||
|
CASE default !------------------------------------------------!
|
||||||
|
PRINT*,'Green_Func: WARNING! Adjoint approach is not defined.'
|
||||||
|
RETURN
|
||||||
|
END SELECT
|
||||||
|
|
||||||
|
gam11 = m(1,1)
|
||||||
|
gam21 = m(2,1)
|
||||||
|
gam31 = m(3,1)
|
||||||
|
gam41 = m(4,1)
|
||||||
|
gam01 = m(0,1)
|
||||||
|
|
||||||
|
alp12 = m(1,2)/m(1,1)
|
||||||
|
alp13 = m(1,3)/m(1,1)
|
||||||
|
alp14 = m(1,4)/m(1,1)
|
||||||
|
alp10 = m(1,0)/m(1,1)
|
||||||
|
|
||||||
|
gam22 = m(2,2)-gam21*alp12
|
||||||
|
gam32 = m(3,2)-gam31*alp12
|
||||||
|
gam42 = m(4,2)-gam41*alp12
|
||||||
|
gam02 = m(0,2)-gam01*alp12
|
||||||
|
|
||||||
|
alp23 = gam32/gam22
|
||||||
|
alp24 = gam42/gam22
|
||||||
|
alp20 = gam02/gam22
|
||||||
|
|
||||||
|
gam33 = m(3,3)-gam31*alp13-gam32*alp23
|
||||||
|
gam43 = m(4,3)-gam41*alp13-gam42*alp23
|
||||||
|
gam03 = m(0,3)-gam01*alp13-gam02*alp23
|
||||||
|
|
||||||
|
alp34 = gam43/gam33
|
||||||
|
alp30 = gam03/gam33
|
||||||
|
|
||||||
|
gam44 = m(4,4)-gam41*alp14-gam42*alp24-gam43*alp34
|
||||||
|
gam04 = m(0,4)-gam01*alp14-gam02*alp24-gam03*alp34
|
||||||
|
|
||||||
|
alp40 = gam04/gam44
|
||||||
|
|
||||||
|
gam00 = m(0,0)-gam01*alp10-gam02*alp20-gam03*alp30-gam04*alp40
|
||||||
|
|
||||||
|
bet1 = g(1)/m(1,1)
|
||||||
|
bet2 = (g(2)-gam21*bet1)/gam22
|
||||||
|
bet3 = (g(3)-gam31*bet1-gam32*bet2)/gam33
|
||||||
|
bet4 = (g(4)-gam41*bet1-gam42*bet2-gam43*bet3)/gam44
|
||||||
|
bet0 = (g(0)-gam01*bet1-gam02*bet2-gam03*bet3-gam04*bet4)/gam00
|
||||||
|
|
||||||
|
d0 = bet0
|
||||||
|
sfd(4) = bet4-alp40*d0
|
||||||
|
sfd(3) = bet3-alp30*d0-alp34*sfd(4)
|
||||||
|
sfd(2) = bet2-alp20*d0-alp24*sfd(4)-alp23*sfd(3)
|
||||||
|
sfd(1) = bet1-alp10*d0-alp14*sfd(4)-alp13*sfd(3)-alp12*sfd(2)
|
||||||
|
|
||||||
|
fc_old = fc
|
||||||
|
mu_old = mu
|
||||||
|
Zeff_old = Zeff
|
||||||
|
sfdx(1:4) = sfd(1:4)
|
||||||
|
|
||||||
|
END SUBROUTINE SpitzFuncCoeff
|
||||||
|
|
||||||
|
SUBROUTINE SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu)
|
||||||
|
!=======================================================================
|
||||||
|
! Calculates the "Spitzer function" in high velocity limit, relativistic
|
||||||
|
! formulation: Lin-Liu et al., Phys.Pl. (2003),v10, 4064, Eq.(33).
|
||||||
|
!
|
||||||
|
! Inputs:
|
||||||
|
! Zeff - effective charge
|
||||||
|
! fc - fraction of circulating electrons
|
||||||
|
! u - p/(m*vte)
|
||||||
|
! q - p/mc
|
||||||
|
! gam - relativ. factor
|
||||||
|
!
|
||||||
|
! Outputs:
|
||||||
|
! K - Spitzer function
|
||||||
|
! dKdu - its derivative
|
||||||
|
!=======================================================================
|
||||||
|
use const_and_precisions, only : zero,one
|
||||||
|
use numint, only : quanc8
|
||||||
|
IMPLICIT NONE
|
||||||
|
REAL(wp_), INTENT(in) :: Zeff,fc,u,q,gam
|
||||||
|
REAL(wp_), INTENT(out) :: K,dKdu
|
||||||
|
INTEGER :: nfun
|
||||||
|
REAL(wp_) :: gam2,err,flag,Integr
|
||||||
|
REAL(wp_), PARAMETER :: a = zero, b = one, rtol = 1e-4_wp_, atol = 1e-12_wp_
|
||||||
|
|
||||||
|
r2 = (1+Zeff)/fc ! global parameter needed for integrand, HSL_f(t)
|
||||||
|
|
||||||
|
IF (u < 1e-2) THEN
|
||||||
|
K = u**4/(r2+4)
|
||||||
|
dKdu = 4*u**3/(r2+4)
|
||||||
|
RETURN
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
q2 = q*q ! for the integrand, HSL_f
|
||||||
|
gp1 = gam+1 ! ..
|
||||||
|
CALL quanc8(HSL_f,zero,one,atol,rtol,Integr,err,nfun,flag)
|
||||||
|
|
||||||
|
gam2 = gam*gam
|
||||||
|
K = u**4 * Integr
|
||||||
|
dKdu = (u/gam)**3 * (1-r2*gam2*Integr)
|
||||||
|
END SUBROUTINE SpitzFunc_HighSpeedLimit
|
||||||
|
|
||||||
|
FUNCTION HSL_f(t) RESULT(f)
|
||||||
|
!=======================================================================
|
||||||
|
! Integrand for the high-speed limit approach (Lin-Liu's formulation)
|
||||||
|
!=======================================================================
|
||||||
|
IMPLICIT NONE
|
||||||
|
REAL(wp_), INTENT(in) :: t
|
||||||
|
REAL(wp_) :: f,g
|
||||||
|
g = sqrt(1+t*t*q2)
|
||||||
|
f = t**(3+r2)/g**3 * (gp1/(g+1))**r2
|
||||||
|
END FUNCTION HSL_f
|
||||||
|
|
||||||
|
end module eccd
|
||||||
884
src/gray.f
884
src/gray.f
File diff suppressed because it is too large
Load Diff
@ -1,437 +0,0 @@
|
|||||||
!########################################################################
|
|
||||||
|
|
||||||
MODULE green_func_p
|
|
||||||
|
|
||||||
!########################################################################
|
|
||||||
!
|
|
||||||
! The module contains few subroutines which are requested to calculate
|
|
||||||
! the current drive value by adjoint approach
|
|
||||||
!
|
|
||||||
!########################################################################
|
|
||||||
USE const_and_precisions, only : wp_
|
|
||||||
!-------
|
|
||||||
IMPLICIT NONE
|
|
||||||
CHARACTER(Len=1), PRIVATE :: adj_appr(6) ! adjoint approach switcher
|
|
||||||
!-------
|
|
||||||
REAL(wp_), PRIVATE :: r2,q2,gp1
|
|
||||||
!-------
|
|
||||||
REAL(wp_), PRIVATE, PARAMETER :: delta = 1e-4 ! border for recalculation
|
|
||||||
!------- for N.M. subroutines (variational principle) -------
|
|
||||||
REAL(wp_), PRIVATE :: sfd(1:4)
|
|
||||||
INTEGER, PRIVATE, PARAMETER :: nre = 2 ! order of rel. correct.
|
|
||||||
REAL(wp_), PRIVATE, PARAMETER :: vp_mee(0:4,0:4,0:2) = &
|
|
||||||
RESHAPE((/0.0, 0.0, 0.0, 0.0, 0.0, &
|
|
||||||
0.0, 0.184875, 0.484304, 1.06069, 2.26175, &
|
|
||||||
0.0, 0.484304, 1.41421, 3.38514, 7.77817, &
|
|
||||||
0.0, 1.06069, 3.38514, 8.73232, 21.4005, &
|
|
||||||
0.0, 2.26175, 7.77817, 21.4005, 55.5079, &
|
|
||||||
! &
|
|
||||||
0.0, -1.33059,-2.57431, -5.07771, -10.3884, &
|
|
||||||
-0.846284,-1.46337, -1.4941, -0.799288, 2.57505, &
|
|
||||||
-1.1601, -1.4941, 2.25114, 14.159, 50.0534, &
|
|
||||||
-1.69257, -0.799288, 14.159, 61.4168, 204.389, &
|
|
||||||
-2.61022, 2.57505, 50.0534, 204.389, 683.756, &
|
|
||||||
! &
|
|
||||||
0.0, 2.62498, 0.985392,-5.57449, -27.683, &
|
|
||||||
0.0, 3.45785, 5.10096, 9.34463, 22.9831, &
|
|
||||||
-0.652555, 5.10096, 20.5135, 75.8022, 268.944, &
|
|
||||||
-2.11571, 9.34463, 75.8022, 330.42, 1248.69, &
|
|
||||||
-5.38358, 22.9831, 268.944, 1248.69, 4876.48/),&
|
|
||||||
(/5,5,3/))
|
|
||||||
REAL(wp_), PRIVATE, PARAMETER :: vp_mei(0:4,0:4,0:2) = &
|
|
||||||
RESHAPE((/0.0, 0.886227, 1.0, 1.32934, 2.0, &
|
|
||||||
0.886227,1.0, 1.32934, 2.0, 3.32335, &
|
|
||||||
1.0, 1.32934, 2.0, 3.32335, 6.0, &
|
|
||||||
1.32934, 2.0, 3.32335, 6.0, 11.6317, &
|
|
||||||
2.0, 3.32335, 6.0, 11.6317, 24.0, &
|
|
||||||
! &
|
|
||||||
0.0, 0.332335, 1.0, 2.49251, 6.0, &
|
|
||||||
1.66168, 1.0, 2.49251, 6.0, 14.5397, &
|
|
||||||
3.0, 2.49251, 6.0, 14.5397, 36.0, &
|
|
||||||
5.81586, 6.0, 14.5397, 36.0, 91.5999, &
|
|
||||||
12.0, 14.5397, 36.0, 91.5999, 240.0, &
|
|
||||||
! &
|
|
||||||
0.0, -0.103855, 0.0, 1.09047, 6.0, &
|
|
||||||
0.726983,0.0, 1.09047, 6.0, 24.5357, &
|
|
||||||
3.0, 1.09047, 6.0, 24.5357, 90.0, &
|
|
||||||
9.81427, 6.0, 24.5357, 90.0, 314.875, &
|
|
||||||
30.0, 24.5357, 90.0, 314.875, 1080.0 /), &
|
|
||||||
(/5,5,3/))
|
|
||||||
REAL(wp_), PRIVATE, PARAMETER :: vp_oee(0:4,0:4,0:2) = &
|
|
||||||
RESHAPE((/0.0, 0.56419, 0.707107, 1.0073, 1.59099, &
|
|
||||||
0.56419, 0.707107, 1.0073, 1.59099, 2.73981, &
|
|
||||||
0.707107,1.0073, 1.59099, 2.73981, 5.08233, &
|
|
||||||
1.0073, 1.59099, 2.73981, 5.08233, 10.0627, &
|
|
||||||
1.59099, 2.73981, 5.08233, 10.0627, 21.1138, &
|
|
||||||
! &
|
|
||||||
0.0, 1.16832, 1.90035, 3.5758, 7.41357, &
|
|
||||||
2.17562, 1.90035, 3.5758, 7.41357, 16.4891, &
|
|
||||||
3.49134, 3.5758, 7.41357, 16.4891, 38.7611, &
|
|
||||||
6.31562, 7.41357, 16.4891, 38.7611, 95.4472, &
|
|
||||||
12.4959, 16.4891, 38.7611, 95.4472, 244.803, &
|
|
||||||
! &
|
|
||||||
0.0, 2.65931, 4.64177, 9.6032, 22.6941, &
|
|
||||||
4.8652, 4.64177, 9.6032, 22.6941, 59.1437, &
|
|
||||||
9.51418, 9.6032, 22.6941, 59.1437, 165.282, &
|
|
||||||
21.061, 22.6941, 59.1437, 165.282, 485.785, &
|
|
||||||
50.8982, 59.1437, 165.282, 485.785, 1483.22/), &
|
|
||||||
(/5,5,3/))
|
|
||||||
REAL(wp_), PRIVATE, PARAMETER :: vp_g(0:4,0:2) = &
|
|
||||||
RESHAPE((/1.32934, 2.0, 3.32335, 6.0, 11.6317, &
|
|
||||||
2.49251, 0.0, 2.90793, 12.0, 39.2571, &
|
|
||||||
1.09047, 6.0, 11.45, 30.0, 98.9606/), &
|
|
||||||
(/5,3/))
|
|
||||||
!########################################################################
|
|
||||||
|
|
||||||
CONTAINS
|
|
||||||
|
|
||||||
!#######################################################################
|
|
||||||
|
|
||||||
SUBROUTINE Setup_SpitzFunc
|
|
||||||
!=======================================================================
|
|
||||||
IMPLICIT NONE
|
|
||||||
!=======================================================================
|
|
||||||
adj_appr(1) = 'l' ! collisionless limit
|
|
||||||
! adj_appr(1) = 'c' ! collisional (classical) limit, w/o trap. part.
|
|
||||||
adj_appr(2) = 'm' ! momentum conservation
|
|
||||||
! adj_appr(2) = 'h' ! high-speed limit
|
|
||||||
!---
|
|
||||||
adj_appr(3) = 'l' ! DO NOT CHANGE!
|
|
||||||
adj_appr(4) = 'r' ! DO NOT CHANGE!
|
|
||||||
adj_appr(5) = 'v' ! DO NOT CHANGE!
|
|
||||||
adj_appr(6) = 'i' ! DO NOT CHANGE!
|
|
||||||
!=======================================================================
|
|
||||||
!.....
|
|
||||||
!=======================================================================
|
|
||||||
RETURN
|
|
||||||
END SUBROUTINE Setup_SpitzFunc
|
|
||||||
|
|
||||||
|
|
||||||
SUBROUTINE GenSpitzFunc(Zeff,fc,u,q,gam, K,dKdu)
|
|
||||||
!=======================================================================
|
|
||||||
! Author: N.B.Marushchenko
|
|
||||||
! June 2005: as start point the subroutine of Ugo Gasparino (198?)
|
|
||||||
! SpitzFunc() is taken and modified.
|
|
||||||
! 1. adapted to the Fortran-95
|
|
||||||
! 2. derivative of Spitzer function is added
|
|
||||||
! 3. separation for 2 brunches is done:
|
|
||||||
! 1st is referenced as 'with conservation of the moment',
|
|
||||||
! 2nd - as 'high speed limit'.
|
|
||||||
! The last one is taken from the Lin-Liu formulation
|
|
||||||
! (Phys.Plasmas 10 (2003) 4064) with K = F*fc.
|
|
||||||
! The asymptotical high speed limit (Taguchi-Fisch model)
|
|
||||||
! is also included as the reference case.
|
|
||||||
! Feb. 2008: non-relativ. version is replaced by the relativistic one;
|
|
||||||
! the method is the the same, but the trial-function is
|
|
||||||
! based on the relativistic formulation.
|
|
||||||
! The relativistic corrections for the collisional operator
|
|
||||||
! up to the second order, i.e. (1/mu)**2, are applied.
|
|
||||||
! Sep. 2008: generalized Spitzer function for arbitrary collisionality
|
|
||||||
! is implemented. The model is based on the concept of
|
|
||||||
! the "effective trapped particles fraction".
|
|
||||||
! The different.-integral kinetic equation for the generalized
|
|
||||||
! Spitzer function is produced with help of subroutines
|
|
||||||
! ArbColl_TrappFract_Array and ArbColl_SpitzFunc_Array,
|
|
||||||
! where the subroutines of H. Maassberg are called).
|
|
||||||
!========================================================================
|
|
||||||
! Spitzer function with & w/o trapped particle effects is given by:
|
|
||||||
!
|
|
||||||
! K(x) = x/gamma*(d1*x+d2*x^2+d4*x^3+d4*x^4),
|
|
||||||
!
|
|
||||||
! where x = v/v_th and gamma=1 for non-relativistic version (Ugo),
|
|
||||||
! or x = p/p_th for relativistic version (N.M., February 2008).
|
|
||||||
! Note, that somewhere the function F(x) instead of K(x) is applied,
|
|
||||||
!
|
|
||||||
! F(x) = K(x)/fc.
|
|
||||||
!
|
|
||||||
! Numerical inversion of the 5x5 symmetric matrix obtained from the
|
|
||||||
! generalized Spitzer problem (see paper of Taguchi for the equation
|
|
||||||
! and paper of Hirshman for the variational approach bringing to the
|
|
||||||
! matrix to be inverted).
|
|
||||||
!
|
|
||||||
! The numerical method used is an improved elimination scheme
|
|
||||||
! (Banachiewiczs-Cholesky-Crout method).
|
|
||||||
! This method is particularly simple for symmetric matrix.
|
|
||||||
! As a reference see "Mathematical Handbook" by Korn & Korn, p.635-636.
|
|
||||||
!
|
|
||||||
! Refs.: 1. S.P. Hirshman, Phys. Fluids 23 (1980) 1238
|
|
||||||
! 2. M. Rome' et al., Plasma Phys. Contr. Fus. 40 (1998) 511
|
|
||||||
! 3. N.B. Marushchenko et al., Fusion Sci. Technol. 55 (2009) 180
|
|
||||||
!========================================================================
|
|
||||||
! INPUTS:
|
|
||||||
! u - p/sqrt(2mT)
|
|
||||||
! q - p/mc;
|
|
||||||
! gam - relativistic factor;
|
|
||||||
! Zeff - effective charge;
|
|
||||||
! fc - fraction of circulating particles.
|
|
||||||
!
|
|
||||||
! OUTPUTS:
|
|
||||||
! K - Spitzer's function
|
|
||||||
! dKdu = dK/du, i.e. its derivative over normalized momentum
|
|
||||||
!=======================================================================
|
|
||||||
use const_and_precisions, only : comp_eps
|
|
||||||
IMPLICIT NONE
|
|
||||||
REAL(wp_), INTENT(in) :: Zeff,fc,u,q,gam
|
|
||||||
REAL(wp_), INTENT(out) :: K,dKdu
|
|
||||||
REAL(wp_) :: gam1,gam2,gam3
|
|
||||||
!=======================================================================
|
|
||||||
K = 0
|
|
||||||
dKdu = 0
|
|
||||||
IF (u < comp_eps) RETURN
|
|
||||||
!---
|
|
||||||
SELECT CASE(adj_appr(2))
|
|
||||||
CASE('m') !--------------- momentum conservation ------------------!
|
|
||||||
gam1 = gam !
|
|
||||||
IF (adj_appr(4) == 'n') gam1 = 1 !
|
|
||||||
gam2 = gam1*gam1 !
|
|
||||||
gam3 = gam1*gam2 !
|
|
||||||
K = u/gam1*u*(sfd(1)+u*(sfd(2)+u*(sfd(3)+u*sfd(4)))) !
|
|
||||||
dKdu = u/gam3* (sfd(1)*(1+ gam2)+u*(sfd(2)*(1+2*gam2)+ & !
|
|
||||||
u*(sfd(3)*(1+3*gam2)+u* sfd(4)*(1+4*gam2)))) !
|
|
||||||
!--------------------- end momentum conservation -------------------!
|
|
||||||
CASE('h') !---------------- high-speed-limit ----------------------!
|
|
||||||
IF (adj_appr(4) == 'n') THEN !- non-relativ. asymptotic form -!
|
|
||||||
K = u**4 *fc/(Zeff+1+4*fc) !- (Taguchi-Fisch model) -!
|
|
||||||
dKdu = 4*u**3 *fc/(Zeff+1+4*fc) !
|
|
||||||
ELSEIF (adj_appr(4) == 'r') THEN !- relativistic, Lin-Liu form. -!
|
|
||||||
CALL SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu) !
|
|
||||||
ENDIF !
|
|
||||||
CASE default !----------------------------------------------------!
|
|
||||||
PRINT*,'GenSpitzFunc: WARNING! Spitzer function is not defined.'
|
|
||||||
RETURN
|
|
||||||
END SELECT
|
|
||||||
!=======================================================================
|
|
||||||
RETURN
|
|
||||||
END SUBROUTINE GenSpitzFunc
|
|
||||||
|
|
||||||
!#######################################################################
|
|
||||||
!#######################################################################
|
|
||||||
!#######################################################################
|
|
||||||
|
|
||||||
SUBROUTINE SpitzFuncCoeff(mu,Zeff,fc)
|
|
||||||
!=======================================================================
|
|
||||||
! Calculates the matrix coefficients required for the subroutine
|
|
||||||
! "GenSpitzFunc", where the Spitzer function is defined through the
|
|
||||||
! variational principle.
|
|
||||||
!
|
|
||||||
! Weakly relativistic (upgraded) version (10.09.2008).
|
|
||||||
! Apart of the non-relativistic matrix coefficients, taken from the
|
|
||||||
! old subroutine of Ugo Gasparino, the relativistic correction written
|
|
||||||
! as series in 1/mu^n (mu=mc2/T) powers is added. Two orders are taken
|
|
||||||
! into account, i.e. n=0,1,2.
|
|
||||||
!
|
|
||||||
! In this version, the coefficients "oee", i.e. Omega_ij, are formulated
|
|
||||||
! for arbitrary collisionality.
|
|
||||||
!
|
|
||||||
! INPUT VARIABLES:
|
|
||||||
! rho = sqrt(SS) with SS - flux-surface label (norm. magn. flux)
|
|
||||||
! ne - density, 1/m^3
|
|
||||||
! mu - mc2/Te
|
|
||||||
! Zeff - effective charge
|
|
||||||
! fc - fraction of circulating particles
|
|
||||||
!
|
|
||||||
! OUTPUT VARIABLES (defined as a global ones):
|
|
||||||
! sfd(1),...,sfd(4) - coefficients of the polynomial expansion of the
|
|
||||||
! "Spitzer"-function (the same as in the Hirshman paper)
|
|
||||||
!=======================================================================
|
|
||||||
use const_and_precisions, only : mc2_
|
|
||||||
IMPLICIT NONE
|
|
||||||
REAL(wp_), INTENT(in) :: mu,Zeff,fc
|
|
||||||
INTEGER :: n,i,j
|
|
||||||
REAL(wp_) :: rtc,rtc1,y,tn(1:nre)
|
|
||||||
REAL(wp_) :: m(0:4,0:4),g(0:4)
|
|
||||||
REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, &
|
|
||||||
gam22,gam32,gam42,gam02, &
|
|
||||||
gam33,gam43,gam03, &
|
|
||||||
gam44,gam04,gam00
|
|
||||||
REAL(wp_) :: alp12,alp13,alp14,alp10, &
|
|
||||||
alp23,alp24,alp20, &
|
|
||||||
alp34,alp30,alp40
|
|
||||||
REAL(wp_) :: bet0,bet1,bet2,bet3,bet4,d0
|
|
||||||
LOGICAL :: renew,rel,newmu,newZ,newfc
|
|
||||||
REAL(wp_), SAVE :: sfdx(1:4) = 0
|
|
||||||
REAL(wp_), SAVE :: mu_old =-1, Zeff_old =-1, fc_old =-1
|
|
||||||
!=======================================================================
|
|
||||||
rel = mu < mc2_
|
|
||||||
newmu = abs(mu -mu_old ) > delta*mu
|
|
||||||
newZ = abs(Zeff-Zeff_old) > delta*Zeff
|
|
||||||
newfc = abs(fc -fc_old ) > delta*fc
|
|
||||||
SELECT CASE(adj_appr(1))
|
|
||||||
CASE ('l','c')
|
|
||||||
renew = (newmu .and. rel) .OR. newZ .OR. newfc
|
|
||||||
END SELECT
|
|
||||||
!---
|
|
||||||
IF (.not.renew) THEN
|
|
||||||
sfd(:) = sfdx(:)
|
|
||||||
RETURN
|
|
||||||
ENDIF
|
|
||||||
!=======================================================================
|
|
||||||
tn(:) = 0
|
|
||||||
IF (adj_appr(4) == 'r') THEN
|
|
||||||
IF (nre > 0) THEN
|
|
||||||
!mu = min(mu,1.e3*mc2_)
|
|
||||||
tn(1) = 1/mu
|
|
||||||
DO n=2,min(2,nre)
|
|
||||||
tn(n) = tn(n-1)/mu
|
|
||||||
ENDDO
|
|
||||||
ENDIF
|
|
||||||
ENDIF
|
|
||||||
!---
|
|
||||||
SELECT CASE(adj_appr(1))
|
|
||||||
CASE ('l','c') !---- both classical & collisionless limits ----!
|
|
||||||
rtc = (1-fc)/fc; rtc1 = rtc+1 !
|
|
||||||
!--- !
|
|
||||||
DO i=0,4 !
|
|
||||||
g(i) = vp_g(i,0) !
|
|
||||||
DO n=1,min(2,nre) !
|
|
||||||
g(i) = g(i) + tn(n)*vp_g(i,n) !
|
|
||||||
ENDDO !
|
|
||||||
!--- !
|
|
||||||
DO j=0,4 !
|
|
||||||
IF (i == 0 .or. j == 0 .or. j >= i) THEN !
|
|
||||||
y = vp_mee(i,j,0) + rtc *vp_oee(i,j,0) + & !
|
|
||||||
Zeff*rtc1*vp_mei(i,j,0) !
|
|
||||||
DO n=1,min(2,nre) !
|
|
||||||
y = y + (vp_mee(i,j,n) + rtc *vp_oee(i,j,n) + & !
|
|
||||||
Zeff*rtc1*vp_mei(i,j,n))*tn(n) !
|
|
||||||
ENDDO !
|
|
||||||
m(i,j) = y !
|
|
||||||
ENDIF !
|
|
||||||
ENDDO !
|
|
||||||
ENDDO !
|
|
||||||
DO i=2,4 !
|
|
||||||
DO j=1,i-1 !
|
|
||||||
m(i,j) = m(j,i) !
|
|
||||||
ENDDO !
|
|
||||||
ENDDO !
|
|
||||||
m(0,0) = 0 !
|
|
||||||
CASE default !------------------------------------------------!
|
|
||||||
PRINT*,'Green_Func: WARNING! Adjoint approach is not defined.'
|
|
||||||
RETURN
|
|
||||||
END SELECT
|
|
||||||
!=======================================================================
|
|
||||||
gam11 = m(1,1)
|
|
||||||
gam21 = m(2,1)
|
|
||||||
gam31 = m(3,1)
|
|
||||||
gam41 = m(4,1)
|
|
||||||
gam01 = m(0,1)
|
|
||||||
!
|
|
||||||
alp12 = m(1,2)/m(1,1)
|
|
||||||
alp13 = m(1,3)/m(1,1)
|
|
||||||
alp14 = m(1,4)/m(1,1)
|
|
||||||
alp10 = m(1,0)/m(1,1)
|
|
||||||
!
|
|
||||||
gam22 = m(2,2)-gam21*alp12
|
|
||||||
gam32 = m(3,2)-gam31*alp12
|
|
||||||
gam42 = m(4,2)-gam41*alp12
|
|
||||||
gam02 = m(0,2)-gam01*alp12
|
|
||||||
!
|
|
||||||
alp23 = gam32/gam22
|
|
||||||
alp24 = gam42/gam22
|
|
||||||
alp20 = gam02/gam22
|
|
||||||
!
|
|
||||||
gam33 = m(3,3)-gam31*alp13-gam32*alp23
|
|
||||||
gam43 = m(4,3)-gam41*alp13-gam42*alp23
|
|
||||||
gam03 = m(0,3)-gam01*alp13-gam02*alp23
|
|
||||||
!
|
|
||||||
alp34 = gam43/gam33
|
|
||||||
alp30 = gam03/gam33
|
|
||||||
!
|
|
||||||
gam44 = m(4,4)-gam41*alp14-gam42*alp24-gam43*alp34
|
|
||||||
gam04 = m(0,4)-gam01*alp14-gam02*alp24-gam03*alp34
|
|
||||||
!
|
|
||||||
alp40 = gam04/gam44
|
|
||||||
!
|
|
||||||
gam00 = m(0,0)-gam01*alp10-gam02*alp20-gam03*alp30-gam04*alp40
|
|
||||||
!
|
|
||||||
bet1 = g(1)/m(1,1)
|
|
||||||
bet2 = (g(2)-gam21*bet1)/gam22
|
|
||||||
bet3 = (g(3)-gam31*bet1-gam32*bet2)/gam33
|
|
||||||
bet4 = (g(4)-gam41*bet1-gam42*bet2-gam43*bet3)/gam44
|
|
||||||
bet0 = (g(0)-gam01*bet1-gam02*bet2-gam03*bet3-gam04*bet4)/gam00
|
|
||||||
!
|
|
||||||
d0 = bet0
|
|
||||||
sfd(4) = bet4-alp40*d0
|
|
||||||
sfd(3) = bet3-alp30*d0-alp34*sfd(4)
|
|
||||||
sfd(2) = bet2-alp20*d0-alp24*sfd(4)-alp23*sfd(3)
|
|
||||||
sfd(1) = bet1-alp10*d0-alp14*sfd(4)-alp13*sfd(3)-alp12*sfd(2)
|
|
||||||
!=======================================================================
|
|
||||||
fc_old = fc
|
|
||||||
mu_old = mu
|
|
||||||
Zeff_old = Zeff
|
|
||||||
!---
|
|
||||||
sfdx(1:4) = sfd(1:4)
|
|
||||||
!=======================================================================
|
|
||||||
RETURN
|
|
||||||
END SUBROUTINE SpitzFuncCoeff
|
|
||||||
|
|
||||||
!#######################################################################
|
|
||||||
!#######################################################################
|
|
||||||
!#######################################################################
|
|
||||||
|
|
||||||
SUBROUTINE SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu)
|
|
||||||
!=======================================================================
|
|
||||||
! Calculates the "Spitzer function" in high velocity limit, relativistic
|
|
||||||
! formulation: Lin-Liu et al., Phys.Pl. (2003),v10, 4064, Eq.(33).
|
|
||||||
!
|
|
||||||
! Inputs:
|
|
||||||
! Zeff - effective charge
|
|
||||||
! fc - fraction of circulating electrons
|
|
||||||
! u - p/(m*vte)
|
|
||||||
! q - p/mc
|
|
||||||
! gam - relativ. factor
|
|
||||||
!
|
|
||||||
! Outputs:
|
|
||||||
! K - Spitzer function
|
|
||||||
! dKdu - its derivative
|
|
||||||
!=======================================================================
|
|
||||||
use const_and_precisions, only : zero,one
|
|
||||||
use numint, only : quanc8
|
|
||||||
IMPLICIT NONE
|
|
||||||
REAL(wp_), INTENT(in) :: Zeff,fc,u,q,gam
|
|
||||||
REAL(wp_), INTENT(out) :: K,dKdu
|
|
||||||
INTEGER :: nfun
|
|
||||||
REAL(8) :: gam2,err,flag,Integr
|
|
||||||
REAL(8), PARAMETER :: a = 0d0, b = 1d0, rtol = 1d-4, atol = 1d-12
|
|
||||||
!=======================================================================
|
|
||||||
r2 = (1+Zeff)/fc ! global parameter needed for integrand, HSL_f(t)
|
|
||||||
!------------------
|
|
||||||
IF (u < 1e-2) THEN
|
|
||||||
K = u**4/(r2+4)
|
|
||||||
dKdu = 4*u**3/(r2+4)
|
|
||||||
RETURN
|
|
||||||
ENDIF
|
|
||||||
!=======================================================================
|
|
||||||
q2 = q*q ! for the integrand, HSL_f
|
|
||||||
gp1 = gam+1 ! ..
|
|
||||||
!---
|
|
||||||
CALL quanc8(HSL_f,zero,one,atol,rtol,Integr,err,nfun,flag)
|
|
||||||
!=======================================================================
|
|
||||||
gam2 = gam*gam
|
|
||||||
!---
|
|
||||||
K = u**4 * Integr
|
|
||||||
dKdu = (u/gam)**3 * (1-r2*gam2*Integr)
|
|
||||||
!=======================================================================
|
|
||||||
RETURN
|
|
||||||
END SUBROUTINE SpitzFunc_HighSpeedLimit
|
|
||||||
|
|
||||||
!#######################################################################
|
|
||||||
!#######################################################################
|
|
||||||
!#######################################################################
|
|
||||||
|
|
||||||
FUNCTION HSL_f(t) RESULT(f)
|
|
||||||
!=======================================================================
|
|
||||||
! Integrand for the high-speed limit approach (Lin-Liu's formulation)
|
|
||||||
!=======================================================================
|
|
||||||
IMPLICIT NONE
|
|
||||||
REAL(8), INTENT(in) :: t
|
|
||||||
REAL(8) :: f,g
|
|
||||||
g = sqrt(1+t*t*q2)
|
|
||||||
f = t**(3+r2)/g**3 * (gp1/(g+1))**r2
|
|
||||||
END FUNCTION HSL_f
|
|
||||||
|
|
||||||
!#######################################################################
|
|
||||||
|
|
||||||
END MODULE green_func_p
|
|
||||||
|
|
||||||
!#######################################################################
|
|
||||||
Loading…
Reference in New Issue
Block a user