further cleaning of eccd routines
This commit is contained in:
parent
2333f83914
commit
464af38310
494
src/gray.f
494
src/gray.f
@ -5686,7 +5686,11 @@ c local constants
|
||||
c local variables
|
||||
integer :: lrm,ithn
|
||||
real(wp_) :: amu,ratiovgr,rbn,rbx,zeff
|
||||
integer :: npar
|
||||
real(wp_) :: cst2
|
||||
real(wp_), dimension(:), allocatable :: eccdpar
|
||||
c common/external functions/variables
|
||||
real(wp_), external :: fjch,fjncl,fjch0
|
||||
integer :: nhmin,nhmax,iokhawa,ierr
|
||||
real(wp_) :: xg,yg,anpl,anpr,vgm,derdnm,ak0,akinv,fhz,sox,
|
||||
. anprre,anprim,alpha,effjcd,akim,tau,psinv,tekev,dens,ddens,
|
||||
@ -5712,6 +5716,16 @@ c
|
||||
common/btot/btot
|
||||
common/bmxmn/bmax,bmin
|
||||
common/fc/fc
|
||||
c
|
||||
interface
|
||||
subroutine setcdcoeff(amu,zeff,rbn,rbx,fc,psinv,ieccd,
|
||||
. cst2,eccdpar,npar)
|
||||
use const_and_precisions, only : wp_
|
||||
integer :: ieccd,npar
|
||||
real(wp_) :: amu,Zeff,rbn,rbx,fc,psinv,cst2
|
||||
real(wp_), dimension(:), allocatable :: eccdpar
|
||||
end subroutine setcdcoeff
|
||||
end interface
|
||||
c
|
||||
c absorption computation
|
||||
c
|
||||
@ -5746,52 +5760,56 @@ c
|
||||
zeff=fzeff(psinv)
|
||||
rbn=btot/bmin
|
||||
rbx=btot/bmax
|
||||
call eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez,
|
||||
* dens,psinv,ieccd,nhmin,nhmax,ithn,effjcd,iokhawa,ierr)
|
||||
end if
|
||||
|
||||
call setcdcoeff(amu,zeff,rbn,rbx,fc,psinv,ieccd,cst2,
|
||||
. eccdpar,npar)
|
||||
select case(ieccd)
|
||||
case(1)
|
||||
c cohen model
|
||||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax,
|
||||
. ithn,cst2,fjch,eccdpar(1:npar),npar,effjcd,iokhawa,ierr)
|
||||
case(2:9)
|
||||
c no trapping
|
||||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax,
|
||||
. ithn,cst2,fjch0,eccdpar(1:npar),npar,effjcd,iokhawa,ierr)
|
||||
case(10:11)
|
||||
c neoclassical model
|
||||
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax,
|
||||
. ithn,cst2,fjncl,eccdpar(1:npar),npar,effjcd,iokhawa,ierr)
|
||||
CASE DEFAULT
|
||||
effjcd=0.0_wp_
|
||||
print*,'ieccd undefined'
|
||||
ierr=89
|
||||
return
|
||||
end select
|
||||
c
|
||||
deallocate(eccdpar)
|
||||
end if
|
||||
c
|
||||
return
|
||||
end
|
||||
c
|
||||
c
|
||||
c
|
||||
subroutine eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez,
|
||||
. dens,psinv,ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr)
|
||||
use const_and_precisions, only : wp_,pi,qesi=>e_,mesi=>me_,
|
||||
. vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_
|
||||
|
||||
subroutine setcdcoeff(amu,zeff,rbn,rbx,fc,psinv,ieccd,
|
||||
. cst2,eccdpar,npar)
|
||||
use const_and_precisions, only : wp_
|
||||
use green_func_p, only: SpitzFuncCoeff
|
||||
use conical, only : fconic
|
||||
use magsurf_data, only : ch,tjp,tlm,njpt,nlmt
|
||||
use dierckx, only : profil
|
||||
implicit none
|
||||
c 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 :: cst2min=1.0e-6_wp_
|
||||
integer, parameter :: ksp=3
|
||||
c arguments
|
||||
integer :: ieccd,nhmn,nhmx,ithn,iokhawa,ierr
|
||||
real(wp_) :: yg,anpl,anprre,amu,Zeff,rbn,rbx,
|
||||
. fc,dens,psinv,effjcd
|
||||
complex(wp_) :: ex,ey,ez
|
||||
integer :: ieccd,npar
|
||||
real(wp_) :: amu,Zeff,rbn,rbx,fc,psinv,cst2
|
||||
real(wp_), dimension(:), allocatable :: eccdpar ! dimension(max(5,3+nlmt))
|
||||
c local variables
|
||||
integer :: nhn,njp,nlm,npar
|
||||
real(wp_) :: anum,denom,resp,resj,coullog,anucc,alams,
|
||||
. fp0s,pa,cst2
|
||||
integer :: njp,nlm
|
||||
real(wp_) :: alams,fp0s,pa
|
||||
real(wp_) :: chlm(nlmt)
|
||||
real(wp_), dimension(3+2*nlmt) :: eccdpar ! dimension(max(5,3+nlmt))
|
||||
c common/external functions/variables
|
||||
real(wp_), external :: fjch,fjncl,fjch0
|
||||
c
|
||||
anum=0.0_wp_
|
||||
denom=0.0_wp_
|
||||
effjcd=0.0_wp_
|
||||
c
|
||||
coullog=48.0_wp_-log(1.0e7_wp_*dens*mc2m2*amu**2)
|
||||
anucc=canucc*dens*coullog
|
||||
c
|
||||
c nhmx=nhmn
|
||||
c
|
||||
eccdpar(1)=zeff
|
||||
c
|
||||
select case(ieccd)
|
||||
c
|
||||
@ -5802,114 +5820,115 @@ c rbx=B/B_max
|
||||
c cst2=1.0_wp_-B/B_max
|
||||
c alams=sqrt(1-B_min/B_max)
|
||||
c Zeff < 31 !!!
|
||||
c fp0s= P_a (alams)
|
||||
c fp0s= P_a (alams)
|
||||
cst2=1.0_wp_-rbx
|
||||
if(cst2.lt.1.0e-6_wp_) cst2=0.0_wp_
|
||||
if(cst2.lt.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)
|
||||
npar=5
|
||||
allocate(eccdpar(npar))
|
||||
eccdpar(1)=zeff
|
||||
eccdpar(2)=rbn
|
||||
eccdpar(3)=alams
|
||||
eccdpar(4)=pa
|
||||
eccdpar(5)=fp0s
|
||||
do nhn=nhmn,nhmx
|
||||
call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
|
||||
. fjch,eccdpar,5,resj,resp,iokhawa,ierr)
|
||||
anum=anum+resj
|
||||
denom=denom+resp
|
||||
end do
|
||||
|
||||
c
|
||||
case(2:9)
|
||||
cst2=0.0_wp_
|
||||
do nhn=nhmn,nhmx
|
||||
call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
|
||||
. fjch0,eccdpar,1,resj,resp,iokhawa,ierr)
|
||||
anum=anum+resj
|
||||
denom=denom+resp
|
||||
end do
|
||||
|
||||
npar=1
|
||||
allocate(eccdpar(npar))
|
||||
eccdpar(1)=zeff
|
||||
c
|
||||
case(10:11)
|
||||
c neoclassical model:
|
||||
c ft=1-fc trapped particle fraction
|
||||
c rzfc=(1+Zeff)/fc
|
||||
cst2=1.0_wp_-rbx
|
||||
if(cst2.lt.1.0e-6_wp_) cst2=0.0_wp_
|
||||
if(cst2.lt.cst2min) cst2=0.0_wp_
|
||||
call SpitzFuncCoeff(amu,Zeff,fc)
|
||||
eccdpar(2) = fc
|
||||
eccdpar(3) = rbx
|
||||
|
||||
njp=njpt
|
||||
nlm=nlmt
|
||||
call profil(0,tjp,njp,tlm,nlm,ch,ksp,ksp,sqrt(psinv),nlm,chlm,
|
||||
. ierr)
|
||||
if(ierr.gt.0) print*,' 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
|
||||
do nhn=nhmn,nhmx
|
||||
call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
|
||||
. fjncl,eccdpar,npar,resj,resp,iokhawa,ierr)
|
||||
anum=anum+resj
|
||||
denom=denom+resp
|
||||
end do
|
||||
|
||||
CASE DEFAULT
|
||||
print*,'ieccd undefined'
|
||||
|
||||
c
|
||||
end select
|
||||
c
|
||||
c effjpl = <J_parallel>/<p_d> /(B_min/<B>) [A m /W]
|
||||
c
|
||||
if(denom.gt.0.0_wp_) effjcd=-ceff*anum/(anucc*denom)
|
||||
return
|
||||
end
|
||||
c
|
||||
c
|
||||
c
|
||||
subroutine curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
|
||||
* fcur,eccdpar,necp,resj,resp,iokhawa,ierr)
|
||||
use const_and_precisions, only : wp_
|
||||
subroutine eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmn,nhmx,
|
||||
* ithn,cst2,fcur,eccdpar,necp,effjcd,iokhawa,ierr)
|
||||
use const_and_precisions, only : wp_,pi,qesi=>e_,mesi=>me_,
|
||||
. vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_
|
||||
use quadpack, only : dqagsmv
|
||||
implicit none
|
||||
c 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_
|
||||
integer, parameter :: lw=5000,liw=lw/4,nfpp=15
|
||||
real(wp_), parameter :: dumin=1.0e-6_wp_
|
||||
integer, parameter :: lw=5000,liw=lw/4,nfpp=13
|
||||
c arguments
|
||||
integer :: i,nhn,ithn,necp,iokhawa,ierr
|
||||
real(wp_) :: yg,anpl,anprre,amu,cst2,resj,resp
|
||||
integer :: i,nhmn,nhmx,ithn,necp,iokhawa,ierr
|
||||
real(wp_) :: yg,anpl,anprre,dens,amu,cst2,effjcd
|
||||
real(wp_), dimension(necp) :: eccdpar
|
||||
complex(wp_) :: ex,ey,ez
|
||||
c local variables
|
||||
integer :: neval,ier,last,npar
|
||||
integer :: nhn,neval,ier,last,npar
|
||||
integer, dimension(liw) :: iw
|
||||
real(wp_) :: anpl2,dnl,ygn,ygn2,resj1,resj2,rdu2,upltp,upltm,
|
||||
. rdu,rdut,rdu2t,duu,uu1,uu2,xx1,xx2,ej,ej1,ej2,epp,uplp,uplm,
|
||||
. cst
|
||||
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+necp) :: apar
|
||||
real(wp_), dimension(0:1) :: uleft,uright
|
||||
c common/external functions/variables
|
||||
real(wp_), external :: fcur,fpp
|
||||
c
|
||||
c EC power and current densities
|
||||
c effjpl = <J_parallel>/<p_d> /(B_min/<B>) [A m /W]
|
||||
c
|
||||
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)
|
||||
c
|
||||
npar=nfpp+necp
|
||||
apar(nfpp+1:npar) = eccdpar(1:necp)
|
||||
c
|
||||
anpl2=anpl*anpl
|
||||
dnl=1.0_wp_-anpl2
|
||||
ygn=nhn*yg
|
||||
ygn2=ygn*ygn
|
||||
|
||||
resj=0.0_wp_
|
||||
resj1=0.0_wp_
|
||||
resj2=0.0_wp_
|
||||
resp=0.0_wp_
|
||||
c
|
||||
rdu2=anpl2+ygn2-1.0_wp_
|
||||
uplp=0.0_wp_
|
||||
uplm=0.0_wp_
|
||||
upltp=0.0_wp_
|
||||
upltm=0.0_wp_
|
||||
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
|
||||
c
|
||||
if (rdu2.ge.0.0_wp_) then
|
||||
rdu2=anpl2+ygn2-1.0_wp_
|
||||
c
|
||||
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
|
||||
c
|
||||
@ -5922,92 +5941,72 @@ c
|
||||
if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl
|
||||
duu=abs(uu1-uu2)
|
||||
c
|
||||
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)
|
||||
if(duu.le.dumin) cycle
|
||||
c
|
||||
apar(12) = dble(nhn)
|
||||
apar(13) = uplp
|
||||
apar(14) = uplm
|
||||
apar(15) = ygn
|
||||
apar(13) = ygn
|
||||
c
|
||||
npar=nfpp+necp
|
||||
apar(nfpp+1:npar) = eccdpar(1:necp)
|
||||
c
|
||||
if(duu.gt.1.0e-6_wp_) then
|
||||
call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp,
|
||||
. epp,neval,ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) ierr=90
|
||||
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
|
||||
c
|
||||
rdu2t=cst2*anpl2+ygn2-1.0_wp_
|
||||
c
|
||||
if (rdu2t.lt.0.0_wp_.or.cst2.eq.0.0_wp_) then
|
||||
c
|
||||
c resonance curve does not cross the trapping region
|
||||
c
|
||||
iokhawa=0
|
||||
if(duu.gt.1.0e-4_wp_) then
|
||||
call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
|
||||
. resj,ej,neval,ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) ierr=91
|
||||
end if
|
||||
else
|
||||
if (rdu2t.gt.0.0_wp_.and.cst2.gt.0.0_wp_) then
|
||||
c
|
||||
c resonance curve crosses the trapping region
|
||||
c
|
||||
iokhawa=1
|
||||
rdut=sqrt(rdu2t)
|
||||
cst=sqrt(cst2)
|
||||
upltm=(cst2*anpl*ygn-cst*rdut)/(1.0_wp_-cst2*anpl2)
|
||||
upltp=(cst2*anpl*ygn+cst*rdut)/(1.0_wp_-cst2*anpl2)
|
||||
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
|
||||
c
|
||||
uu1=uplm
|
||||
uu2=upltm
|
||||
xx1=amu*(anpl*uu1+ygn-1.0_wp_)
|
||||
xx2=amu*(anpl*uu2+ygn-1.0_wp_)
|
||||
if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
|
||||
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.gt.1.0e-6_wp_) then
|
||||
call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
|
||||
. resj1,ej1,neval,ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) then
|
||||
if (abs(resj1).lt.1.0e-10_wp_) then
|
||||
resj1=0.0_wp_
|
||||
else
|
||||
ierr=92
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
c resonance curve does not cross the trapping region
|
||||
c
|
||||
uu1=upltp
|
||||
uu2=uplp
|
||||
xx1=amu*(anpl*uu1+ygn-1.0_wp_)
|
||||
xx2=amu*(anpl*uu2+ygn-1.0_wp_)
|
||||
if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
|
||||
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.gt.1.0e-6_wp_) then
|
||||
call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
|
||||
. resj2,ej2,neval,ier,liw,lw,last,iw,w)
|
||||
if (ier.gt.0) then
|
||||
if(ier.ne.2) ierr=93
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
resj=resj1+resj2
|
||||
iokhawa=0
|
||||
uleft(0)=uplm
|
||||
uright(0)=uplp
|
||||
end if
|
||||
c
|
||||
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
|
||||
c
|
||||
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
|
||||
|
||||
return
|
||||
@ -6035,9 +6034,7 @@ c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c extrapar(13) = ygn
|
||||
c
|
||||
use const_and_precisions, only : wp_,ui=>im
|
||||
use math, only : fact
|
||||
@ -6048,7 +6045,7 @@ c arguments
|
||||
real(wp_), dimension(npar) :: extrapar
|
||||
c local variables
|
||||
integer :: ithn,nhn,nm,np
|
||||
real(wp_) :: yg,anpl,amu,anprre,uplp,uplm,ygn,upr,upr2,gam,ee,
|
||||
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
|
||||
c
|
||||
@ -6061,23 +6058,21 @@ c
|
||||
ez=cmplx(extrapar(9),extrapar(10),wp_)
|
||||
ithn=int(extrapar(11))
|
||||
nhn=int(extrapar(12))
|
||||
uplp=extrapar(13)
|
||||
uplm=extrapar(14)
|
||||
ygn=extrapar(15)
|
||||
ygn=extrapar(13)
|
||||
|
||||
upr2=(1.0_wp_-anpl**2)*(uplp-upl)*(upl-uplm)
|
||||
gam=anpl*upl+ygn
|
||||
upr2=gam*gam-1.0_wp_-upl*upl
|
||||
ee=exp(-amu*(gam-1))
|
||||
|
||||
thn2=1.0_wp_
|
||||
thn2u=upr2*thn2
|
||||
! 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
|
||||
if(ithn.eq.1) then
|
||||
c 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)
|
||||
@ -6110,31 +6105,19 @@ c
|
||||
function fjch(upl,extrapar,npar)
|
||||
c integration variable upl passed explicitly. Other variables passed
|
||||
c as array of extra parameters of length npar=size(extrapar).
|
||||
c variables with index 1..15 must be passed to fpp
|
||||
c variable with index 16 is zeff
|
||||
c variables with index gt 16 are specific of the cd model
|
||||
c variables with index 1..13 must be passed to fpp
|
||||
c variable with index 14 is zeff
|
||||
c variables with index gt 14 are specific of the cd model
|
||||
c
|
||||
c extrapar(1) = yg
|
||||
c extrapar(2) = anpl
|
||||
c extrapar(3) = amu
|
||||
c extrapar(4) = Re(anprw)
|
||||
c extrapar(5) = Re(ex)
|
||||
c extrapar(6) = Im(ex)
|
||||
c extrapar(7) = Re(ey)
|
||||
c extrapar(8) = Im(ey)
|
||||
c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c extrapar(13) = ygn
|
||||
c
|
||||
c extrapar(16) = zeff
|
||||
c extrapar(17) = rb
|
||||
c extrapar(18) = alams
|
||||
c extrapar(19) = pa
|
||||
c extrapar(20) = fp0s
|
||||
c extrapar(14) = zeff
|
||||
c extrapar(15) = rb
|
||||
c extrapar(16) = alams
|
||||
c extrapar(17) = pa
|
||||
c extrapar(18) = fp0s
|
||||
c
|
||||
use const_and_precisions, only : wp_
|
||||
use conical, only : fconic
|
||||
@ -6145,25 +6128,23 @@ c arguments
|
||||
real(wp_), dimension(npar) :: extrapar
|
||||
c local variables
|
||||
integer :: nhn
|
||||
real(wp_) :: anpl,anprre,uplp,uplm,ygn,zeff,rb,alams,pa,fp0s,
|
||||
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,fpp
|
||||
complex(wp_) :: ex,ey,ez
|
||||
c
|
||||
anpl=extrapar(2)
|
||||
anprre=extrapar(4)
|
||||
uplp=extrapar(13)
|
||||
uplm=extrapar(14)
|
||||
ygn=extrapar(15)
|
||||
zeff=extrapar(16)
|
||||
rb=extrapar(17)
|
||||
alams=extrapar(18)
|
||||
pa=extrapar(19)
|
||||
fp0s=extrapar(20)
|
||||
ygn=extrapar(13)
|
||||
zeff=extrapar(14)
|
||||
rb=extrapar(15)
|
||||
alams=extrapar(16)
|
||||
pa=extrapar(17)
|
||||
fp0s=extrapar(18)
|
||||
|
||||
upr2=(1.0_wp_-anpl**2)*(uplp-upl)*(upl-uplm)
|
||||
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
|
||||
@ -6187,7 +6168,7 @@ c
|
||||
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,npar)
|
||||
fjch=eta*fpp(upl,extrapar(1:13),13)
|
||||
c
|
||||
return
|
||||
end
|
||||
@ -6201,23 +6182,10 @@ c variables with index 1..15 must be passed to fpp
|
||||
c variable with index 16 is zeff
|
||||
c variables with index gt 16 are specific of the cd model
|
||||
c
|
||||
c extrapar(1) = yg
|
||||
c extrapar(2) = anpl
|
||||
c extrapar(3) = amu
|
||||
c extrapar(4) = Re(anprw)
|
||||
c extrapar(5) = Re(ex)
|
||||
c extrapar(6) = Im(ex)
|
||||
c extrapar(7) = Re(ey)
|
||||
c extrapar(8) = Im(ey)
|
||||
c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c extrapar(13) = ygn
|
||||
c
|
||||
c extrapar(16) = zeff
|
||||
c extrapar(14) = zeff
|
||||
c
|
||||
use const_and_precisions, only : wp_
|
||||
implicit none
|
||||
@ -6232,8 +6200,8 @@ c local variables
|
||||
complex(wp_) :: ex,ey,ez
|
||||
c
|
||||
anpl=extrapar(2)
|
||||
ygn=extrapar(15)
|
||||
zeff=extrapar(16)
|
||||
ygn=extrapar(13)
|
||||
zeff=extrapar(14)
|
||||
|
||||
gam=anpl*upl+ygn
|
||||
u2=gam*gam-1.0_wp_
|
||||
@ -6248,7 +6216,7 @@ c
|
||||
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,npar)
|
||||
fjch0=eta*fpp(upl,extrapar(1:13),13)
|
||||
c
|
||||
return
|
||||
end
|
||||
@ -6258,31 +6226,19 @@ c
|
||||
function fjncl(upl,extrapar,npar)
|
||||
c integration variable upl passed explicitly. Other variables passed
|
||||
c as array of extra parameters of length npar=size(extrapar).
|
||||
c variables with index 1..15 must be passed to fpp
|
||||
c variable with index 16 is zeff
|
||||
c variables with index gt 16 are specific of the cd model
|
||||
c variables with index 1..13 must be passed to fpp
|
||||
c variable with index 14 is zeff
|
||||
c variables with index gt 14 are specific of the cd model
|
||||
c
|
||||
c extrapar(1) = yg
|
||||
c extrapar(2) = anpl
|
||||
c extrapar(3) = amu
|
||||
c extrapar(4) = Re(anprw)
|
||||
c extrapar(5) = Re(ex)
|
||||
c extrapar(6) = Im(ex)
|
||||
c extrapar(7) = Re(ey)
|
||||
c extrapar(8) = Im(ey)
|
||||
c extrapar(9) = Re(ez)
|
||||
c extrapar(10) = Im(ez)
|
||||
c extrapar(11) = double(ithn)
|
||||
c extrapar(12) = double(nhn)
|
||||
c extrapar(13) = uplp
|
||||
c extrapar(14) = uplm
|
||||
c extrapar(15) = ygn
|
||||
c extrapar(13) = ygn
|
||||
c
|
||||
c extrapar(16) = zeff
|
||||
c extrapar(17) = fc
|
||||
c extrapar(18) = rbx
|
||||
c extrapar(19:18+(npar-18)/2) = tlm
|
||||
c extrapar(19+(npar-18)/2:npar) = chlm
|
||||
c extrapar(14) = zeff
|
||||
c extrapar(15) = fc
|
||||
c extrapar(16) = rbx
|
||||
c extrapar(17:16+(npar-16)/2) = tlm
|
||||
c extrapar(17+(npar-16)/2:npar) = chlm
|
||||
c
|
||||
use const_and_precisions, only : wp_
|
||||
use green_func_p, only: GenSpitzFunc
|
||||
@ -6300,15 +6256,15 @@ c local variables
|
||||
. bth,uth,fk,dfk,alam,fu,dfu,eta,fpp
|
||||
c local variables
|
||||
integer :: ier
|
||||
real(wp_), dimension((npar-18)/2) :: wrk
|
||||
real(wp_), dimension((npar-16)/2) :: wrk
|
||||
real(wp_), dimension(1) :: xs,ys
|
||||
c
|
||||
anpl=extrapar(2)
|
||||
amu=extrapar(3)
|
||||
ygn=extrapar(15)
|
||||
zeff=extrapar(16)
|
||||
fc=extrapar(17)
|
||||
rbx=extrapar(18)
|
||||
ygn=extrapar(13)
|
||||
zeff=extrapar(14)
|
||||
fc=extrapar(15)
|
||||
rbx=extrapar(16)
|
||||
|
||||
gam=anpl*upl+ygn
|
||||
u2=gam*gam-1.0_wp_
|
||||
@ -6321,52 +6277,26 @@ c
|
||||
dfk=dfk*(2.0_wp_/amu)*bth
|
||||
|
||||
alam=upr2/u2/rbx
|
||||
nlm=(npar-18)/2
|
||||
xs(1)=alam
|
||||
nlm=(npar-16)/2
|
||||
c
|
||||
call splev(extrapar(19:18+nlm),nlm,extrapar(19+nlm:npar),3,
|
||||
c extrapar(17:16+(npar-16)/2) = tlm
|
||||
c extrapar(17+(npar-16)/2:npar) = chlm
|
||||
c
|
||||
call splev(extrapar(17:16+nlm),nlm,extrapar(17+nlm:npar),3,
|
||||
. xs(1),ys(1),1,ier)
|
||||
fu=ys(1)
|
||||
call splder(extrapar(19:18+nlm),nlm,extrapar(19+nlm:npar),3,1,
|
||||
call splder(extrapar(17:16+nlm),nlm,extrapar(17+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,npar)
|
||||
fjncl=eta*fpp(upl,extrapar(1:13),13)
|
||||
return
|
||||
end
|
||||
c
|
||||
c
|
||||
c
|
||||
C subroutine vlambda(alam,tlm,chlm,nlmt,fv,dfv)
|
||||
C use const_and_precisions, only : wp_
|
||||
C use dierckx, only : splev,splder
|
||||
C implicit none
|
||||
Cc local constants
|
||||
C integer, parameter :: ksp=3
|
||||
Cc arguments
|
||||
C integer :: nlmt
|
||||
C real(wp_) :: alam,fv,dfv
|
||||
C real(wp_), dimension(nlmt) :: tlm,chlm
|
||||
Cc local variables
|
||||
C integer :: nlm,ier
|
||||
C real(wp_), dimension(nlmt) :: wrk
|
||||
C real(wp_), dimension(1) :: xxs,ffs
|
||||
Cc
|
||||
C nlm=nlmt
|
||||
C xxs(1)=alam
|
||||
Cc
|
||||
C call splev(tlm,nlm,chlm,ksp,xxs(1),ffs(1),1,ier)
|
||||
C fv=ffs(1)
|
||||
Cc
|
||||
C call splder(tlm,nlm,chlm,ksp,1,xxs(1),ffs(1),1,wrk,ier)
|
||||
C dfv=ffs(1)
|
||||
Cc
|
||||
C return
|
||||
C end
|
||||
c
|
||||
c
|
||||
c
|
||||
subroutine projxyzt(iproj,nfile)
|
||||
use const_and_precisions, only : wp_
|
||||
|
||||
Loading…
Reference in New Issue
Block a user