diff --git a/Makefile b/Makefile index 520296b..d9a88e5 100644 --- a/Makefile +++ b/Makefile @@ -3,8 +3,8 @@ EXE=gray # Objects list MAINOBJ=gray.o -OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eierf.o \ - graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.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 \ interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ reflections.o simplespline.o utils.o beamdata.o @@ -25,17 +25,18 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) $(FC) $(FFLAGS) -o $@ $^ # Dependencies on modules -gray.o: const_and_precisions.o conical.o dierckx.o dispersion.o \ - graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \ +gray.o: const_and_precisions.o dierckx.o dispersion.o eccd.o \ + graydata_anequil.o graydata_flags.o graydata_par.o \ interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ reflections.o simplespline.o utils.o beamdata.o conical.o: const_and_precisions.o dierckx.o: const_and_precisions.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_flags.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 magsurf_data.o: const_and_precisions.o math.o: const_and_precisions.o diff --git a/src/eccd.f90 b/src/eccd.f90 new file mode 100644 index 0000000..ea4b391 --- /dev/null +++ b/src/eccd.f90 @@ -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(cst20) 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 = / /(B_min/) [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 \ No newline at end of file diff --git a/src/gray.f b/src/gray.f index c106093..9f11f9d 100644 --- a/src/gray.f +++ b/src/gray.f @@ -655,7 +655,6 @@ c subroutine read_data use const_and_precisions, only : wp_,qe=>ecgs_,me=>mecgs_, . vc=>ccgs_,cvdr=>degree,pi - use green_func_p, only:Setup_SpitzFunc use graydata_flags use graydata_par use graydata_anequil @@ -743,7 +742,7 @@ c c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models c read(2,*) ieccd - if (ieccd.ge.10) call Setup_SpitzFunc + !if (ieccd.ge.10) call Setup_SpitzFunc c c ipec=0/1 :pec profiles grid in psi/rhop c nnd :number of grid steps for pec profiles +1 @@ -1639,13 +1638,13 @@ c locate psi surface for q=1.5 and q=2 call vmaxmini(qpsi,nr,qmin,qmax,iqmn,iqmx) if (q15.gt.qmin.and.q15.lt.qmax) then call surfq(q15,psi15) - rhot15=frhotor(psi15) + rhot15=frhotor(sqrt(psi15)) print'(3(a,f8.5))','psi_1.5 = ',psi15,' rhop_1.5 = ' . ,sqrt(psi15),' rhot_1.5 = ',rhot15 end if if (q2.gt.qmin.and.q2.lt.qmax) then call surfq(q2,psi2) - rhot2=frhotor(psi2) + rhot2=frhotor(sqrt(psi2)) print'(3(a,f8.5))','psi_2 = ',psi2,' rhop_2 = ' . ,sqrt(psi2),' rhot_2 = ',rhot2 end if @@ -1852,7 +1851,7 @@ c else qq=qpsi(nr) end if - rhot=frhotor(psin) + rhot=frhotor(rhop) call tor_curr_psi(psin,ajphi) write(55,111) psin,rhot,dens,te,qq,ajphi*1.0e-6_wp_ end do @@ -1881,7 +1880,7 @@ c do i=1,nst psin=dble(i-1)/dble(nst-1) rhop=sqrt(psin) - rhot=frhotor(psin) + rhot=frhotor(rhop) call density(psin) te=temp(psin) write(55,111) psin,rhot,dens,te @@ -2176,41 +2175,42 @@ c c c c - function frhotor(psi) + function frhotor(rhop) use const_and_precisions, only : wp_ use graydata_flags, only : iequil implicit none c arguments - real(wp_) :: psi,frhotor + real(wp_) :: rhop,frhotor c common/external functions/variables real(wp_) :: frhotor_eq,frhotor_an c - if(iequil.eq.2) frhotor=frhotor_eq(sqrt(psi)) - if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi)) + if(iequil.eq.1) then + frhotor=frhotor_an(rhop) + else + frhotor=frhotor_eq(rhop) + end if c return end c c c - function frhotor_av(psi) + function frhotor_av(rhop) use const_and_precisions, only : wp_ use magsurf_data, only : rpstab, crhotq, npsi - use simplespline, only :spli + use simplespline, only : spli + use utils, only : locate implicit none c arguments - real(wp_) :: psi,frhotor_av + real(wp_) :: rhop,frhotor_av c local variables integer :: ip - real(wp_) :: rpsi,dps + real(wp_) :: drh c - rpsi=sqrt(psi) - ip=int((npsi-1)*rpsi+1) -c if(ip.eq.0) ip=1 -c if(ip.eq.npsi) ip=npsi-1 + call locate(rpstab,npsi,rhop,ip) ip=min(max(1,ip),npsi-1) - dps=rpsi-rpstab(ip) - frhotor_av=spli(crhotq,npsi,ip,dps) + drh=rhop-rpstab(ip) + frhotor_av=spli(crhotq,npsi,ip,drh) c return end @@ -2244,7 +2244,7 @@ c do i=1,nnr rhop(i)=(i-1)*dr psin=rhop(i)*rhop(i) - rhot(i)=frhotor(psin) + rhot(i)=frhotor(rhop(i)) wp(i)=1.0_wp_ end do wp(1)=1.0e3_wp_ @@ -4841,7 +4841,7 @@ c c c function temp(arg) - use const_and_precisions, only : wp_ + use const_and_precisions, only : wp_,zero,one use graydata_flags, only : iprof use graydata_anequil, only : te0,dte0,alt1,alt2 use interp_eqprof, only : psrad,ct,npp @@ -4854,15 +4854,14 @@ c local variables integer :: k real(wp_) :: proft,dps c - temp=0.0_wp_ - if(arg.ge.1.0_wp_.or.arg.lt.0.0_wp_) return + temp=zero + if(arg.ge.one.or.arg.lt.zero) return if(iprof.eq.0) then proft=(1.0_wp_-arg**alt1)**alt2 temp=(te0-dte0)*proft+dte0 else call locate(psrad,npp,arg,k) - if(k.eq.0) k=1 - if(k.eq.npp) k=npp-1 + k=max(1,min(k,npp-1)) dps=arg-psrad(k) temp=spli(ct,npp,k,dps) endif @@ -4872,8 +4871,8 @@ c c c c - function fzeff(arg) - use const_and_precisions, only : wp_ + function fzeff(ps) + use const_and_precisions, only : wp_,zero,one use graydata_flags, only : iprof use graydata_anequil, only : zeffan use interp_eqprof, only : psrad,cz,npp @@ -4881,14 +4880,14 @@ c use simplespline, only :spli implicit none c arguments - real(wp_) :: arg,fzeff + real(wp_), intent(in) :: ps + real(wp_) :: fzeff c local variables integer :: k - real(wp_) :: ps,dps + real(wp_) :: dps c - fzeff=1 - ps=arg - if(ps.gt.1.0_wp_.and.ps.lt.0.0_wp_) return + fzeff=one + if(ps.gt.one.or.ps.lt.zero) return if(iprof.eq.0) then fzeff=zeffan else @@ -5535,70 +5534,50 @@ c implicit none c local variables integer :: j,k - real(wp_) :: dr,r1,r2,summ,sm + real(wp_) :: dr,r0,r,w,w0,summ c - dr=1.0_wp_ - if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) - r1=0.0_wp_ - summ=0.0_wp_ q(1)=1.0_wp_ - if(nrayr.gt.1) then - do j=1,nrayr - r2=(dble(j)-0.5_wp_)*dr - q(j)=(exp(-2.0_wp_*r1**2)-exp(-2.0_wp_*r2**2)) - r1=r2 - summ=summ+q(j) - end do -c - q(1)=q(1)/summ - sm=q(1) - j=1 - k=1 - do j=2,nrayr - q(j)=q(j)/nrayth/summ - do k=1,nrayth - sm=sm+q(j) - end do - end do - end if + if(nrayr.le.1) return + + dr=rwmax/dble(nrayr-1) + + summ=0.0_wp_ + w0=1.0_wp_ + do j=1,nrayr + r=(dble(j)-0.5_wp_)*dr + w=exp(-2.0_wp_*r**2) + q(j)=w-w0 + r0=r + w0=w + summ=summ+q(j) + end do + q=q/summ + q(2:)=q(2:)/nrayth c return end c c c - subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi, - . bmxi,bmni,fci,intp) + subroutine valpsisplpec(rhop,voli,areai) use const_and_precisions, only : wp_ - use simplespline, only :spli,splid - use magsurf_data, only : rpstab,cvol,crri,crbav,cbmx,cbmn, - . carea,cfc,npsi + use utils, only : locate + use simplespline, only :spli + use magsurf_data, only : rpstab,cvol,carea,npsi implicit none c arguments - integer :: intp - real(wp_) :: rpsi,voli,dervoli,areai,rrii,rbavi,bmxi,bmni,fci + real(wp_), intent(in) :: rhop + real(wp_), intent(out) :: voli,areai c local variables integer :: ip - real(wp_) :: dps + real(wp_) :: drh c - ip=int((npsi-1)*rpsi+1) -c if(ip.eq.0) ip=1 -c if(ip.eq.npsi) ip=npsi-1 + call locate(rpstab,npsi,rhop,ip) ip=min(max(1,ip),npsi-1) + drh=rhop-rpstab(ip) c - dps=rpsi-rpstab(ip) -c - areai=spli(carea,npsi,ip,dps) - voli=spli(cvol,npsi,ip,dps) - dervoli=splid(cvol,npsi,ip,dps) - rrii=spli(crri,npsi,ip,dps) -c - if(intp.eq.0) return -c - rbavi=spli(crbav,npsi,ip,dps) - bmxi=spli(cbmx,npsi,ip,dps) - bmni=spli(cbmn,npsi,ip,dps) - fci=spli(cfc,npsi,ip,dps) + voli=spli(cvol,npsi,ip,drh) + areai=spli(carea,npsi,ip,drh) c return end @@ -5608,6 +5587,7 @@ c subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) use const_and_precisions, only : wp_ use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi + use utils, only : locate use simplespline, only :spli implicit none c arguments @@ -5616,9 +5596,7 @@ c local variables integer :: ip real(wp_) :: dps c - ip=int((npsi-1)*rpsi+1) -c if(ip.eq.0) ip=1 -c if(ip.eq.npsi) ip=npsi-1 + call locate(rpstab,npsi,rpsi,ip) ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) ratjai=spli(cratja,npsi,ip,dps) @@ -5631,93 +5609,43 @@ c c c subroutine pabs_curr(i,j,k) - use const_and_precisions, only : wp_,pi - use graydata_flags, only : dst + use const_and_precisions, only : wp_,pi,mc2=>mc2_ + use graydata_flags, only : dst,iwarm,ilarm,ieccd,imx use graydata_par, only : sgnbphi use graydata_anequil, only : rr0m use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,didst, . ccci,q,tau1v + use dispersion, only : harmnumber, warmdisp + use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl implicit none c arguments - integer :: i,j,k -c local variables - integer :: intp - real(wp_) :: dvvl,rbavi,rrii,rhop0,rhop,voli,dervoli,area,adnm, - . dtau,effjcdav,dcidst,dpdst -c common/external functions/variables - real(wp_) :: p0mw,dersdst,alpha,effjcd,akim,tau0,psinv, - . bmxi,bmni,fci -c - common/p0/p0mw - common/dersdst/dersdst - common/absor/alpha,effjcd,akim,tau0 - common/psival/psinv - common/bmxmn/bmxi,bmni - common/fc/fci -c - dvvl=1.0_wp_ - rbavi=1.0_wp_ - rrii=rr0m - rhop0=sqrt(psjki(j,k,i-1)) - intp=1 - psinv=psjki(j,k,i) - rhop=sqrt(psinv) - call valpsispl(rhop,voli,dervoli,area,rrii, - . rbavi,bmxi,bmni,fci,intp) - dvvl=dervoli*abs(rhop-rhop0) - if(dvvl.le.0.0_wp_) dvvl=1.0_wp_ -c - adnm=2.0_wp_*pi*rrii -c -c calcolo della corrente cohen: currtot [MA] -c calcolo della densita' corrente cohen: currj [MA/m^2] -c calcolo della densita' potenza: pdjki [MW/m^3] -c calcolo della efficienza : effjcdav [A m/W ] -c - tau0=tauv(j,k,i-1) -c - call ecrh_cd -c - alphav(j,k,i)=alpha - dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0_wp_ - tauv(j,k,i)=tauv(j,k,i-1)+dtau - dpdst=p0mw*q(j)*exp(-tauv(j,k,i)-tau1v(j,k))* - . alphav(j,k,i)*dersdst - pdjki(j,k,i)=dpdst*dst/dvvl - ppabs(j,k,i)=p0mw*q(j)*exp(-tau1v(j,k))* - . (1.0_wp_-exp(-tauv(j,k,i))) - effjcdav=rbavi*effjcd - currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i) - didst(j,k,i)=sgnbphi*effjcdav*dpdst/adnm - dcidst=(didst(j,k,i)+didst(j,k,i-1))/2.0_wp_ - ccci(j,k,i)=ccci(j,k,i-1)+dcidst*dst - return - end -c -c -c - subroutine ecrh_cd - use const_and_precisions, only : wp_,mc2=>mc2_ - use dispersion, only : harmnumber, warmdisp - use graydata_flags, only : iwarm,ilarm,ieccd,imx - implicit none + integer, intent(in) :: i,j,k c local constants real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ c local variables + integer :: intp + real(wp_) :: dvvl,rbavi,rrii,rhop0,rhop,voli,dervoli,area, + . tau,effjcdav,dpdst,p0,pow,alpha0,didst0,ccci0 integer :: lrm,ithn real(wp_) :: amu,ratiovgr,rbn,rbx,zeff - integer :: npar - real(wp_) :: cst2 + real(wp_) :: cst2,bmxi,bmni,fci 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, - . btot,bmax,bmin,fc + real(wp_) :: p0mw,ak0,akinv,fhz,dersdst + real(wp_) :: psinv,xg,yg,tekev,dens,ddens,btot + real(wp_) :: anpl,anpr,vgm,derdnm,sox,anprre,anprim,alpha,effjcd, + . akim,tau0 complex(wp_) :: ex,ey,ez - real(wp_) :: fzeff,temp + real(wp_), external :: fzeff,temp c + common/p0/p0mw + common/dersdst/dersdst + common/absor/alpha,effjcd,akim,tau0 ! solo per print_output + common/psival/psinv + common/tete/tekev ! per print_output + common/dens/dens,ddens + common/btot/btot common/nharm/nhmin,nhmax common/xgxg/xg common/ygyg/yg @@ -5727,25 +5655,26 @@ c common/mode/sox common/nprw/anprre,anprim common/epolar/ex,ey,ez - common/absor/alpha,effjcd,akim,tau - common/ierr/ierr common/iokh/iokhawa - common/psival/psinv - common/tete/tekev - common/dens/dens,ddens - common/btot/btot - common/bmxmn/bmax,bmin - common/fc/fc + common/ierr/ierr +c + dvvl=1.0_wp_ + rbavi=1.0_wp_ + rrii=rr0m + rhop=sqrt(psinv) +c +c +c calcolo della corrente cohen: currtot [MA] +c calcolo della densita' corrente cohen: currj [MA/m^2] +c calcolo della densita' potenza: pdjki [MW/m^3] +c calcolo della efficienza : effjcdav [A m/W ] +c + rhop0=sqrt(psjki(j,k,i-1)) + alpha0=alphav(j,k,i-1) + tau0=tauv(j,k,i-1) + didst0=didst(j,k,i-1) + ccci0=ccci(j,k,i-1) 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 @@ -5754,569 +5683,109 @@ c effjcd=0.0_wp_ c tekev=temp(psinv) + call valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) c - if(tekev.le.0.0_wp_.or.tau.gt.taucr) return + if(tekev.gt.0.0_wp_.and.tau0.le.taucr) then c - amu=mc2/tekev - call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) - if(nhmin.eq.0) return + amu=mc2/tekev + call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) + if(nhmin.gt.0) then c - lrm=max(ilarm,nhmax) - call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, - * iwarm,imx,ex,ey,ez) + lrm=max(ilarm,nhmax) + call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, + * iwarm,imx,ex,ey,ez) c - akim=ak0*anprim - ratiovgr=2.0_wp_*anpr/derdnm -c ratiovgr=2.0_wp_*anpr/derdnm*vgm - alpha=2.0_wp_*akim*ratiovgr - if(alpha.lt.0.0_wp_) then - ierr=94 - print*,' IERR = ', ierr,' alpha negative' - end if -c - if(ieccd.gt.0) then - ithn=1 - if(lrm.gt.nhmin) ithn=2 - zeff=fzeff(psinv) - rbn=btot/bmin - rbx=btot/bmax - - 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 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 -c local constants - real(wp_), parameter :: cst2min=1.0e-6_wp_ - integer, parameter :: ksp=3 -c arguments - 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 :: njp,nlm - real(wp_) :: alams,fp0s,pa - real(wp_) :: chlm(nlmt) -c - select case(ieccd) -c - case(1) -c cohen model -c rbn=B/B_min -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) - cst2=1.0_wp_-rbx - 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 -c - case(2:9) - cst2=0.0_wp_ - 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.cst2min) cst2=0.0_wp_ - call SpitzFuncCoeff(amu,Zeff,fc) - 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 -c - end select -c - return - end -c -c -c - 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_ - real(wp_), parameter :: dumin=1.0e-6_wp_ - integer, parameter :: lw=5000,liw=lw/4,nfpp=13 -c arguments - 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 :: 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+necp) :: apar - real(wp_), dimension(0:1) :: uleft,uright -c common/external functions/variables - real(wp_), external :: fcur,fpp -c -c effjpl = / /(B_min/) [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 -c - 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 - 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 - uu1=uplm - uu2=uplp - xx1=amu*(anpl*uu1+ygn-1.0_wp_) - xx2=amu*(anpl*uu2+ygn-1.0_wp_) -c - 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) -c - if(duu.le.dumin) cycle -c - apar(12) = dble(nhn) - apar(13) = ygn -c - 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.gt.0.0_wp_.and.cst2.gt.0.0_wp_) then -c -c resonance curve crosses the trapping region -c - 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 -c -c resonance curve does not cross the trapping region -c - 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 + akim=ak0*anprim + ratiovgr=2.0_wp_*anpr/derdnm!*vgm + alpha=2.0_wp_*akim*ratiovgr + if(alpha.lt.0.0_wp_) then + ierr=94 + print*,' IERR = ', ierr,' alpha negative' 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 + if(ieccd.gt.0) then + ithn=1 + if(lrm.gt.nhmin) ithn=2 + zeff=fzeff(psinv) + rbn=btot/bmni + rbx=btot/bmxi - return - end -c -c computation of integral for power density, integrand function fpp -c -c ith=0 : polarization term = const -c ith=1 : polarization term Larmor radius expansion to lowest order -c ith=2 : full polarization term (J Bessel) -c - function fpp(upl,extrapar,npar) -c integration variable upl passed explicitly. other variables passed -c as array of extra parameters of length npar=size(extrapar) -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) = ygn -c - use const_and_precisions, only : wp_,ui=>im - use math, only : fact - implicit none -c arguments - integer :: npar - real(wp_) :: upl,fpp - real(wp_), dimension(npar) :: extrapar -c 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 -c - 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 -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) - thn2=(0.5_wp_*cth*abs(emxy+ez*anprre*upl/ygn))**2 - thn2u=upr2*thn2 - else -c 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 + select case(ieccd) + case(1) +c cohen model + call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar) + call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, + . ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierr) + case(2) +c no trapping + call setcdcoeff(zeff,cst2,eccdpar) + call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, + . ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierr) + case default +c neoclassical model + call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar) + call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, + . ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierr) + end select + !deallocate(eccdpar) end if end if end if - - fpp=ee*thn2u c - return - end -c -c computation of integral for current density -c fjch integrand for Cohen model with trapping -c fjch0 integrand for Cohen model without trapping -c fjncl integrand for momentum conserv. model K(u) from Maruschenko -c gg=F(u)/u with F(u) as in Cohen paper -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..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(2) = anpl -c extrapar(4) = Re(anprw) -c extrapar(13) = ygn -c -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 - implicit none -c arguments - integer :: npar - real(wp_) :: upl,fjch - real(wp_), dimension(npar) :: extrapar -c 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,fpp - complex(wp_) :: ex,ey,ez -c - anpl=extrapar(2) - anprre=extrapar(4) - ygn=extrapar(13) - zeff=extrapar(14) - rb=extrapar(15) - alams=extrapar(16) - pa=extrapar(17) - fp0s=extrapar(18) + alphav(j,k,i)=alpha - 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 + tau=tau0+0.5_wp_*(alpha+alpha0)*dersdst*dst + tauv(j,k,i)=tau - 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_) then - dfp0=-fconic(alam,pa,1)/sqrt(1.0_wp_-alam**2) - end if - fh=alam*(1.0_wp_-alams*fp0/(alam*fp0s)) - dfhl=1.0_wp_-alams*dfp0/fp0s + p0=p0mw*q(j)*exp(-tau1v(j,k)) + pow=p0*exp(-tau) + ppabs(j,k,i)=p0-pow - eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam) + dpdst=pow*alpha*dersdst + + dvvl=dervoli*abs(rhop-rhop0) + if(dvvl.le.0.0_wp_) dvvl=1.0_wp_ + pdjki(j,k,i)=dpdst*dst/dvvl + + effjcdav=rbavi*effjcd + currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i) + + didst(j,k,i)=sgnbphi*effjcdav*dpdst/(2.0_wp_*pi*rrii) + + ccci(j,k,i)=ccci0+0.5_wp_*(didst(j,k,i)+didst0)*dst - if(upl.lt.0.0_wp_) eta=-eta - fjch=eta*fpp(upl,extrapar(1:13),13) -c return end c c c - function fjch0(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 -c extrapar(2) = anpl -c extrapar(13) = ygn -c -c extrapar(14) = zeff -c + subroutine valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) use const_and_precisions, only : wp_ + use utils, only : locate + use simplespline, only :spli,splid + use magsurf_data, only : rpstab,crri,crbav,cvol,cbmn,cbmx,cfc,npsi implicit none c arguments - real(wp_) :: upl,fjch0 - integer :: npar - real(wp_), dimension(npar) :: extrapar + real(wp_), intent(in) :: rhop + real(wp_), intent(out) :: rrii,rbavi,dervoli,bmni,bmxi,fci c local variables - integer :: nhn - real(wp_) :: anpl,anprre,ygn,zeff,gam,u2,u,z5,xi,xib,xibi, - . fu2b,fu2,gu,gg,dgg,eta,fpp - complex(wp_) :: ex,ey,ez + integer :: ip + real(wp_) :: drh c - anpl=extrapar(2) - ygn=extrapar(13) - zeff=extrapar(14) - - 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:13),13) + call locate(rpstab,npsi,rhop,ip) + ip=min(max(1,ip),npsi-1) + drh=rhop-rpstab(ip) + rrii=spli(crri,npsi,ip,drh) + rbavi=spli(crbav,npsi,ip,drh) + dervoli=splid(cvol,npsi,ip,drh) + bmni=spli(cbmn,npsi,ip,drh) + bmxi=spli(cbmx,npsi,ip,drh) + fci=spli(cfc,npsi,ip,drh) c return end c c -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..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(2) = anpl -c extrapar(3) = amu -c extrapar(13) = ygn -c -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 - use dierckx, only : splev,splder - implicit none -c arguments - integer :: npar - real(wp_) :: upl,fjncl - real(wp_), dimension(npar) :: extrapar -c local constants - integer, parameter :: ksp=3 -c local variables - integer :: nlm - real(wp_) :: anpl,amu,ygn,zeff,fc,rbx,gam,u2,u,upr2, - . bth,uth,fk,dfk,alam,fu,dfu,eta,fpp -c local variables - integer :: ier - real(wp_), dimension((npar-16)/2) :: wrk - real(wp_), dimension(1) :: xs,ys -c - anpl=extrapar(2) - amu=extrapar(3) - ygn=extrapar(13) - zeff=extrapar(14) - fc=extrapar(15) - rbx=extrapar(16) - - 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-16)/2 -c -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(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(1:13),13) - return - end -c -c c subroutine projxyzt(iproj,nfile) use const_and_precisions, only : wp_ @@ -6449,7 +5918,7 @@ c radial coordinate of i-(i+1) interval mid point end if if (abs(ipec) == 1) then rhop_tab(it)=rt - rhot_tab(it)=frhotor(rt*rt) + rhot_tab(it)=frhotor(rt) rhop1=rt1 else rhot_tab(it)=rt @@ -6459,8 +5928,7 @@ c radial coordinate of i-(i+1) interval mid point c psi grid at mid points, dimension nnd+1, for use in pec_tab rtabpsi1(it)=rhop1**2 - call valpsispl(rhop1,voli1,dervoli,areai1,rrii, - . rbavi,bmxi,bmni,fci,izero) + call valpsisplpec(rhop1,voli1,areai1) dvol(it)=abs(voli1-voli0) darea(it)=abs(areai1-areai0) voli0=voli1 diff --git a/src/green_func_p.f90 b/src/green_func_p.f90 deleted file mode 100644 index 1a6f93a..0000000 --- a/src/green_func_p.f90 +++ /dev/null @@ -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 - -!#######################################################################