From c36ffbc6b63f47cc19273db892cc3590cdc00179 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 12 Nov 2015 11:57:43 +0000 Subject: [PATCH] main subroutine restructured; single indexing of rays; all subroutines updated to Fortran 90 (and all commons removed); two passes feature temporarily disabled. --- Makefile | 28 +- src/beamdata.f90 | 251 +- src/beams.f90 | 34 +- src/const_and_precisions.f90 | 2 + src/dispersion.f90 | 3 +- src/eccd.f90 | 2 +- src/gray-externals.f | 4424 ---------------------------------- src/gray-externals.f90 | 895 +++++++ src/graycore.f90 | 1436 ++++++++++- src/magsurf_data.f90 | 92 + src/main.f90 | 2 +- src/pec.f90 | 385 +++ src/polarization.f90 | 152 ++ 13 files changed, 3064 insertions(+), 4642 deletions(-) delete mode 100644 src/gray-externals.f create mode 100644 src/gray-externals.f90 create mode 100644 src/pec.f90 create mode 100644 src/polarization.f90 diff --git a/Makefile b/Makefile index a375ee5..4933240 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ MAINOBJ=main.o OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \ dierckx.o dispersion.o eccd.o eierf.o graycore.o gray-externals.o \ gray_params.o equilibrium.o magsurf_data.o math.o minpack.o numint.o \ - quadpack.o reflections.o simplespline.o utils.o + pec.o polarization.o quadpack.o reflections.o simplespline.o utils.o # Alternative search paths vpath %.f90 src @@ -28,17 +28,18 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) # Dependencies on modules main.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \ graycore.o gray_params.o reflections.o -graycore.o: const_and_precisions.o beamdata.o beams.o equilibrium.o \ - gray-externals.o gray_params.o reflections.o +graycore.o: const_and_precisions.o beamdata.o beams.o coreprofiles.o \ + dispersion.o equilibrium.o gray-externals.o gray_params.o \ + pec.o polarization.o reflections.o utils.o gray-externals.o: const_and_precisions.o beams.o coreprofiles.o dierckx.o \ - dispersion.o eccd.o gray_params.o \ - equilibrium.o magsurf_data.o math.o numint.o quadpack.o \ - reflections.o simplespline.o utils.o beamdata.o + dispersion.o eccd.o gray_params.o \ + equilibrium.o magsurf_data.o math.o numint.o quadpack.o \ + reflections.o simplespline.o utils.o beamdata.o beams.o: const_and_precisions.o simplespline.o utils.o beamdata.o: const_and_precisions.o gray_params.o conical.o: const_and_precisions.o coreprofiles.o: const_and_precisions.o dierckx.o gray_params.o simplespline.o \ - utils.o + utils.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 @@ -47,10 +48,13 @@ gray_params.o: const_and_precisions.o utils.o equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o \ utils.o gray_params.o magsurf_data.o: const_and_precisions.o gray_params.o equilibrium.o dierckx.o \ - simplespline.o utils.o + reflections.o simplespline.o utils.o math.o: const_and_precisions.o minpack.o: const_and_precisions.o numint.o: const_and_precisions.o +pec.o: const_and_precisions.o beamdata.o equilibrium.o gray_params.o \ + magsurf_data.o utils.o +polarization.o: const_and_precisions.o quadpack.o: const_and_precisions.o reflections.o: const_and_precisions.o utils.o simplespline.o: const_and_precisions.o @@ -58,13 +62,7 @@ utils.o: const_and_precisions.o # General object compilation command %.o: %.f90 - $(FC) $(FFLAGS) -c $< - -%.o: %.f - $(FC) $(FFLAGS) -c $< - -gray-externals.o:gray-externals.f - $(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< + $(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< .PHONY: clean install # Remove output files diff --git a/src/beamdata.f90 b/src/beamdata.f90 index 83221ca..7e89d3d 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -2,39 +2,46 @@ module beamdata use const_and_precisions, only : wp_ implicit none - integer, parameter :: jmx=31,kmx=36,nmx=8000 - integer, save :: nray,nrayr,nrayth,nstep - real(wp_) :: dst,rwmax - real(wp_), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav - real(wp_), dimension(:,:,:), allocatable, save :: pdjki,currj,didst - integer, dimension(:,:), allocatable, save :: iiv,iop,iow,ihcd,istore - real(wp_), dimension(:,:), allocatable, save :: tau1v - real(wp_), dimension(:), allocatable, save :: q - real(wp_), dimension(:,:,:), allocatable, save :: yyrfl !(6,:,:) - real(wp_), dimension(:,:,:), allocatable, save :: ywrk,ypwrk !(6,:,:) - real(wp_), dimension(:,:,:), allocatable, save :: xc,xco,du1,du1o !(3,:,:) - real(wp_), dimension(:,:,:), allocatable, save :: gri,dgrad2v !(3,:,:) - real(wp_), dimension(:,:,:,:), allocatable, save :: ggri !(3,3,:,:) - real(wp_), dimension(:,:), allocatable, save :: grad2 - real(wp_), dimension(:), allocatable, save :: dffiu,ddffiu - complex(wp_), dimension(:,:,:), allocatable, save :: ext,eyt + integer, save :: nray,nrayr,nrayth,nstep,jray1 + real(wp_), save :: dst,h,hh,h6,rwmax,twodr2 + integer, parameter :: nfileproj0 = 8, nfilew = 12 contains - subroutine init_rtr(rtrparam) + subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri, & + psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) use gray_params, only : rtrparam_type + use const_and_precisions, only : zero,half,two implicit none type(rtrparam_type), intent(in) :: rtrparam + real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & + gri,psjki,tauv,alphav,ppabs,didst,ccci + real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri + real(wp_), dimension(:), intent(out), allocatable :: p0jk + complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt + integer, dimension(:), intent(out), allocatable :: iiv dst=rtrparam%dst - rwmax=rtrparam%rwmax + h=dst + hh=h*half + h6=h/6.0_wp_ + nrayr=rtrparam%nrayr nrayth=rtrparam%nrayth if(nrayr==1) nrayth=1 nray=(nrayr-1)*nrayth+1 + + rwmax=rtrparam%rwmax + if(nrayr>1) then + twodr2 = two*(rwmax/(nrayr-1))**2 + else + twodr2 = two + end if + nstep=rtrparam%nstep - call alloc_beam -! call alloc_beam1 + + call alloc_beam(ywork,ypwork,xc,du1,gri,ggri, & + psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) end subroutine init_rtr function rayi2jk(i) result(jk) @@ -45,10 +52,10 @@ contains if (i>1) then ioff = i - 2 - jk(1) = ioff/nrayth ! jr-2 - jk(2) = ioff - (jk(1))*nrayth + 1 ! kt -! jk(2) = mod(ioff,nrayth) + 1 ! kt - jk(1) = jk(1) + 2 ! jr + jk(1) = ioff/nrayth ! jr-2 + jk(2) = ioff - jk(1)*nrayth + 1 ! kt +! jk(2) = mod(ioff,nrayth) + 1 ! kt + jk(1) = jk(1) + 2 ! jr else jk = 1 end if @@ -93,77 +100,149 @@ contains end if end function rayjk2i - subroutine alloc_beam + subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri, & + psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) implicit none + real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & + gri,psjki,tauv,alphav,ppabs,didst,ccci + real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri + real(wp_), dimension(:), intent(out), allocatable :: p0jk + complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt + integer, dimension(:), intent(out), allocatable :: iiv - call dealloc_beam - allocate(psjki(nrayr,nrayth,nstep), ppabs(nrayr,nrayth,nstep), & - pdjki(nrayr,nrayth,nstep), ccci(nrayr,nrayth,nstep), & - currj(nrayr,nrayth,nstep), didst(nrayr,nrayth,nstep), & - tauv(nrayr,nrayth,nstep), alphav(nrayr,nrayth,nstep), & - iiv(nrayr,nrayth), iop(nrayr,nrayth), & - iow(nrayr,nrayth), tau1v(nrayr,nrayth), & - ihcd(nrayr,nrayth), istore(nrayr,nrayth), & - q(nrayr), yyrfl(6,nrayr,nrayth), & - ywrk(6,nrayr,nrayth), ypwrk(6,nrayr,nrayth), & - xc(3,nrayr,nrayth), xco(3,nrayr,nrayth), & - du1(3,nrayr,nrayth), du1o(3,nrayr,nrayth), & - gri(3,nrayr,nrayth), dgrad2v(3,nrayr,nrayth), & - ggri(3,3,nrayr,nrayth), grad2(nrayr,nrayth), & - dffiu(nrayr), ddffiu(nrayr), & - ext(nrayr,nrayth,0:4), eyt(nrayr,nrayth,0:4)) + call dealloc_beam(ywork,ypwork,xc,du1,gri,ggri, & + psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + + allocate(ywork(6,nray), ypwork(6,nray), gri(3,nray), ggri(3,3,nray), & + xc(3,nrayth,nrayr), du1(3,nrayth,nrayr), & + psjki(nray,nstep), tauv(nray,nstep), alphav(nray,nstep), & + ppabs(nray,nstep), didst(nray,nstep), ccci(nray,nstep), & + p0jk(nray), ext(nray), eyt(nray), iiv(nray)) end subroutine alloc_beam -! subroutine alloc_beam1 -! implicit none -! integer, intent(out) :: ierr -! -! call dealloc_beam -! allocate(psjki(nray,nstep), ppabs(nray,nstep), pdjki(nray,nstep), & -! ccci(nray,nstep), currj(nray,nstep), didst(nray,nstep), & -! tauv(nray,nstep), alphav(nray,nstep), & -! ywrk(6,nray), ypwrk(6,nray), & -! xc(3,nray), xco(3,nray), & -! du1(3,nray), du1o(3,nray), & -! grad2(nray), dgrad2v(3,nray), & -! gri(3,nray), ggri(3,3,nray), & -! dffiu(nrayr), ddffiu(nrayr), q(nrayr), & -! iiv(nray), iop(nray), iow(nray), & -! ihcd(nray), istore(nray), tau1v(nray), & -! yyrfl(6,nray), ext(nray,0:4), eyt(nray,0:4)) -! end subroutine alloc_beam1 -! - subroutine dealloc_beam + + subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri, & + psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) implicit none + real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & + gri,psjki,tauv,alphav,ppabs,didst,ccci + real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri + real(wp_), dimension(:), intent(out), allocatable :: p0jk + complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt + integer, dimension(:), intent(out), allocatable :: iiv + + if (allocated(ywork)) deallocate(ywork) + if (allocated(ypwork)) deallocate(ypwork) + if (allocated(xc)) deallocate(xc) + if (allocated(du1)) deallocate(du1) + if (allocated(gri)) deallocate(gri) + if (allocated(ggri)) deallocate(ggri) if (allocated(psjki)) deallocate(psjki) - if (allocated(ppabs)) deallocate(ppabs) - if (allocated(pdjki)) deallocate(pdjki) - if (allocated(ccci)) deallocate(ccci) - if (allocated(currj)) deallocate(currj) - if (allocated(didst)) deallocate(didst) if (allocated(tauv)) deallocate(tauv) if (allocated(alphav)) deallocate(alphav) - if (allocated(iiv)) deallocate(iiv) - if (allocated(iop)) deallocate(iop) - if (allocated(iow)) deallocate(iow) - if (allocated(ihcd)) deallocate(ihcd) - if (allocated(istore)) deallocate(istore) - if (allocated(tau1v)) deallocate(tau1v) - if (allocated(q)) deallocate(q) - if (allocated(yyrfl)) deallocate(yyrfl) - if (allocated(ywrk)) deallocate(ywrk) - if (allocated(ypwrk)) deallocate(ypwrk) - if (allocated(xc)) deallocate(xc) - if (allocated(xco)) deallocate(xco) - if (allocated(du1)) deallocate(du1) - if (allocated(du1o)) deallocate(du1o) - if (allocated(gri)) deallocate(gri) - if (allocated(dgrad2v)) deallocate(dgrad2v) - if (allocated(ggri)) deallocate(ggri) - if (allocated(grad2)) deallocate(grad2) - if (allocated(dffiu)) deallocate(dffiu) - if (allocated(ddffiu)) deallocate(ddffiu) + if (allocated(ppabs)) deallocate(ppabs) + if (allocated(didst)) deallocate(didst) + if (allocated(ccci)) deallocate(ccci) + if (allocated(p0jk)) deallocate(p0jk) if (allocated(ext)) deallocate(ext) if (allocated(eyt)) deallocate(eyt) + if (allocated(iiv)) deallocate(iiv) end subroutine dealloc_beam + + subroutine pweight(p0,p0jk) +! power associated to jk-th ray p0jk(j) for total beam power p0 + use const_and_precisions, only : wp_, zero, one, half, two + implicit none +! arguments + real(wp_), intent(in) :: p0 + real(wp_), dimension(:), intent(out) :: p0jk +! local variables + integer :: j,jk,jkn + real(wp_) :: dr,r,w,r0,w0,summ + real(wp_), dimension(nrayr) :: q + + if(nray==1) then + q(1) = one + else + dr = rwmax/dble(nrayr - 1) + summ = zero + w0 = one + do j = 1, nrayr + r = (dble(j) - half)*dr + w = exp(-two*r**2) + q(j) = w - w0 + summ = summ + q(j) + r0 = r + w0 = w + end do + q = q/summ + q(2:) = q(2:)/nrayth + end if + + p0jk(1)=q(1)*p0 + jk=2 + do j=2,nrayr + jkn=jk+nrayth + p0jk(jk:jkn-1)=q(j)*p0 + jk=jkn + end do + end subroutine pweight + + subroutine print_projxyzt(st,ywrk,iproj) + use const_and_precisions, only : wp_, comp_huge, zero, one + implicit none +! arguments + real(wp_), intent(in) :: st + real(wp_), dimension(:,:), intent(in) :: ywrk + integer, intent(in) :: iproj +! local variables + integer :: jk,jkz,nfile + integer, dimension(2) ::jkv + real(wp_), dimension(3) :: xv1,dir,dxv + real(wp_) :: dirm,rtimn,rtimx,csth1,snth1,csps1,snps1,xti,yti,zti,rti +! common/external functions/variables + + nfile = nfileproj0 + iproj + + xv1 = ywrk(1:3,1) + dir = ywrk(4:6,1) + dirm = sqrt(dir(1)**2 + dir(2)**2 + dir(3)**2) + dir = dir/dirm + csth1 = dir(3) + snth1 = sqrt(one - csth1**2) + if(snth1 > zero) then + csps1=dir(2)/snth1 + snps1=dir(1)/snth1 + else + csps1=one + snps1=zero + end if + + if(iproj==0) then + jkz = nray - nrayth + 1 + else + jkz = 1 + end if + + rtimn = comp_huge + rtimx = zero + do jk = jkz, nray + dxv = ywrk(1:3,jk) - xv1 + xti = dxv(1)*csps1 - dxv(2)*snps1 + yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1 + zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1 + rti = sqrt(xti**2 + yti**2) + + jkv=rayi2jk(jk) + if(.not.(iproj==0 .and. jk==1)) & + write(nfile,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') st,jkv,xti,yti,zti,rti + if(iproj==1 .and. jkv(2)==nrayth) write(nfile,*) ' ' + + if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti + if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti + end do + write(nfile,*) ' ' + write(nfilew,'(3(1x,e16.8e3))') st,rtimn,rtimx + end subroutine print_projxyzt + end module beamdata diff --git a/src/beams.f90 b/src/beams.f90 index 8aef627..7b71504 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -164,34 +164,36 @@ contains end subroutine read_beam1 - subroutine launchangles2n(alpha,beta,x,y,z,anx,any,anz) + subroutine launchangles2n(alpha,beta,xv,anv) use const_and_precisions, only : degree implicit none ! arguments - real(wp_), intent(in) :: alpha,beta,x,y,z - real(wp_), intent(out) :: anx,any,anz + real(wp_), intent(in) :: alpha,beta,xv(3) + real(wp_), intent(out) :: anv(3) ! local variables - real(wp_) :: r,anr,anphi + real(wp_) :: r,anr,anphi,a,b - r=sqrt(x**2+y**2) + r=sqrt(xv(1)**2+xv(2)**2) ! phi=atan2(y,x) print'(a,2f8.3)','alpha0, beta0 = ',alpha,beta - print'(a,4f8.3)','x00, y00, R00, z00 = ',x,y,r,z + print'(a,4f8.3)','x00, y00, R00, z00 = ',xv(1:2),r,xv(3) print*,' ' + a = degree*alpha + b = degree*beta ! ! angles alpha, beta in a local reference system as proposed by Gribov et al ! - anr = -cos(degree*beta)*cos(degree*alpha) - anphi = sin(degree*beta) -! anx = -cos(degree*beta)*cos(degree*alpha) -! any = sin(degree*beta) + anr = -cos(b)*cos(a) + anphi = sin(b) +! anx = -cos(b)*cos(a) +! any = sin(b) - anx = (anr*x - anphi*y)/r - any = (anr*y + anphi*x)/r -! anr = (anx*x + any*y)/r -! anphi = (any*x - anx*y)/r + anv(1) = (anr*xv(1) - anphi*xv(2))/r ! = anx + anv(2) = (anr*xv(2) + anphi*xv(1))/r ! = any +! anr = (anx*xv(1) + any*xv(2))/r +! anphi = (any*xv(1) - anx*xv(2))/r - anz =-cos(degree*beta)*sin(degree*alpha) + anv(3) =-cos(b)*sin(a) ! = anz end subroutine launchangles2n subroutine xgygcoeff(fghz,ak0,bres,xgcn) @@ -212,6 +214,6 @@ contains ! ! xg=xgcn*dens19 ! - xgcn=1.0e-5_wp_*qe**2/(pi*me*fghz**2) ! [10^-19 m^3] + xgcn=4.0e13_wp_*pi*qe**2/(me*omega**2) ! [10^-19 m^3] end subroutine xgygcoeff end module beams diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90 index 0789d2e..326ae57 100755 --- a/src/const_and_precisions.f90 +++ b/src/const_and_precisions.f90 @@ -32,7 +32,9 @@ !======================================================================== integer, parameter :: izero = 0 REAL(wp_), PARAMETER :: zero = 0.0_wp_ + REAL(wp_), PARAMETER :: half = 0.5_wp_ REAL(wp_), PARAMETER :: one = 1.0_wp_ + REAL(wp_), PARAMETER :: two = 2.0_wp_ real(wp_), parameter :: pi = 3.141592653589793_wp_ ! 3.141592653589793238462643383280 real(wp_), parameter :: pihalf = 1.57079632679489661923_wp_ REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_ diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 3818d54..f1010d3 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -200,10 +200,11 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez) ! ! if(i.gt.imx) print*,' i>imx ',yg,errnpr,i ! - if(sqrt(dble(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).eq.huge(one).or. & + if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).eq.huge(one).or. & abs(npr2).le.tiny(one)) then write(*,"(' X =',f7.4,' Y =',f7.4,' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(npr2)) npr2=czero + err=99 end if ! if(dble(npr2).lt.zero) then ! npr2=zero diff --git a/src/eccd.f90 b/src/eccd.f90 index 0b39c07..a8d5e89 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -275,7 +275,7 @@ contains end if resj=0.0_wp_ - do i=0,1 + do i=0,iokhawa resji=0.0_wp_ xx1=amu*(anpl*uleft(i)+ygn-1.0_wp_) xx2=amu*(anpl*uright(i)+ygn-1.0_wp_) diff --git a/src/gray-externals.f b/src/gray-externals.f deleted file mode 100644 index e95dc42..0000000 --- a/src/gray-externals.f +++ /dev/null @@ -1,4424 +0,0 @@ -C program gray -C use gray_params, only : ipass,igrad -C implicit none -Cc local variables -C real(wp_) :: p0mw1 -Cc common/external functions/variables -C integer :: ierr,index_rt -C real(wp_) :: sox,p0mw,powrfl,taumn,taumx,pabstot,currtot, -Cc -C common/ierr/ierr -C common/mode/sox -C common/p0/p0mw -C common/powrfl/powrfl -C common/index_rt/index_rt -C common/taumnx/taumn,taumx,pabstot,currtot -Cc -C if (ipass.gt.1) then -Cc second pass into plasma -C p0mw1=p0mw -C igrad=0 -Cc -C index_rt=2 -C p0mw=p0mw1*powrfl -C call prfile -C call vectinit2 -C call paraminit -C call ic_rt2 -C call gray_integration -C call after_gray_integration -C pabstott=pabstott+pabstot -C currtott=currtott+currtot -Cc -C index_rt=3 -C sox=-sox -C p0mw=p0mw1*(1.0_wp_-powrfl) -C call prfile -C call vectinit2 -C call paraminit -C call ic_rt2 -C call gray_integration -C call after_gray_integration -C pabstott=pabstott+pabstot -C currtott=currtott+currtot -C end if -Cc -C end program gray -c -c -c - subroutine after_gray_integration(pabs,currt,dpdv,ajphiv) - use const_and_precisions, only : wp_,zero - use gray_params, only : iwarm,nnd,ipec - use beamdata, only : nrayr - implicit none -c arguments - real(wp_) :: pabs,currt - real(wp_), dimension(nnd) :: dpdv,ajphiv -c local variables - integer :: iproj,nfilp,i - real(wp_) :: currtka - real(wp_), dimension(nnd) :: pins,currins - real(wp_), dimension(nnd) :: dvol,darea,ratjav,ratjbv,ratjplv - real(wp_), dimension(nnd) :: rt_in,rhop_tab,rhot_tab - real(wp_), dimension(0:nnd) :: rtabpsi1 - real(wp_) :: rhotpav,drhotpav - real(wp_) :: rhotjava,drhotjava -c common/external functions/variables - integer :: index_rt - real(wp_) :: st,taumn,taumx,pabstot,currtot - -! arguments -c - common/ss/st - common/index_rt/index_rt - common/taumnx/taumn,taumx,pabstot,currtot -c -c print all ray positions in local reference system -c - if(nrayr.gt.1) then - iproj=1 - nfilp=9 - call projxyzt(iproj,nfilp) - end if -c -c print final results on screen -c - write(*,*) - write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st - if(iwarm.gt.0) then - write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx - else - write(*,'(a,2f9.4)') 'taumn, taumx = ', zero,zero - end if -c - write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot - currtka =currtot*1.0e3_wp_ - write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka -c -c compute power and current density profiles for all rays -c - pabs=pabstot - currt=currtot - - call pec_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1, - . dvol,darea,ratjav,ratjbv,ratjplv) - - call spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv,pins,currins) - - call postproc_profiles(pabs,currt,rhot_tab,dvol,darea,dpdv,ajphiv, - . rhotpav,drhotpav,rhotjava,drhotjava) - - write(*,*) 'rhotpav,drhotpav,rhotjava,drhotjava' - write(*,99) rhotpav,drhotpav,rhotjava,drhotjava - - - do i=1,nnd - write(48,99) rhop_tab(i),rhot_tab(i),ajphiv(i), - . ajphiv(i)*ratjbv(i),dpdv(i),currins(i),pins(i) - end do - - return -99 format(16(1x,e16.8e3)) - end -c -c -c - subroutine after_onestep(p0,sox,ak0,xgcn,bres,i,istop,ier) - use const_and_precisions, only : wp_,pi - use reflections, only : inside,rlim,zlim,nlim,wall_refl - use gray_params, only : iwarm,istpr0,istpl0,ipass,igrad - use equilibrium, only : zbinf,zbsup - use beamdata, only : nrayr,nrayth,dst,psjki,ppabs,ccci,iiv,tauv, - . ihcd,istore,iop,iow,tau1v,yyrfl,ywrk,ext,eyt - implicit none -c local constants - real(wp_), parameter :: taucr=12.0_wp_ -c arguments - real(wp_), intent(in) :: p0,sox,ak0,xgcn,bres - integer, intent(in) :: i - integer, intent(out) :: istop, ier -c local variables - integer :: iproj,nfilp,irfl,j,k,kkk,iopmin,iowmax,iowmin,ivac,j1, - . k1,kkkk - real(wp_) :: qqout,uuout,vvout,qqin2,uuin2,vvin2,rrm,zzm,anvjkr, - . aknmin,akdotn,anwclr - real(wp_), dimension(6) :: y,dery - real(wp_), dimension(3) :: xvrfl,anvrfl,xvvac,anw - real(wp_), dimension(3,nrayr,nrayth) :: xvjk,anvjk - complex(wp_) :: extr,eytr,exin2,eyin2 - logical :: ins_pl -c common/external functions/variables - integer :: istpr,istpl,jclosest,ierr,index_rt - real(wp_) :: psinv,psinv11,taumn,taumx,pabstot,currtot,psipol, - . chipol,powrfl,strfl11 - real(wp_), dimension(3) :: anwcl,xwcl,xv,anv - logical :: inside_plasma -c - external inside_plasma -c - common/istgr/istpr,istpl - common/refln/anwcl,xwcl,jclosest - common/psival/psinv - common/psinv11/psinv11 - common/ierr/ierr - common/taumnx/taumn,taumx,pabstot,currtot - common/xv/xv - common/anv/anv - common/powrfl/powrfl - common/strfl11/strfl11 - common/index_rt/index_rt -c - pabstot=0.0_wp_ - currtot=0.0_wp_ - taumn=1.0e+30_wp_ - taumx=-1.0e+30_wp_ - psinv11=1.0_wp_ - iopmin=100 - iowmin=100 - iowmax=0 - istop=0 - ier=0 - -c - do j=1,nrayr - kkk=nrayth - if(j.eq.1) kkk=1 - do k=1,kkk - call gwork(sox,xgcn,bres,j,k) -c - if(ierr.gt.0) then - write(*,*) ' IERR = ', ierr - ier=ierr - if(ierr.eq.97) then -c igrad=0 -c ierr=0 - else - istop=1 -c exit - end if - end if - - psjki(j,k,i)=psinv - rrm=sqrt(xv(1)**2+xv(2)**2)/100.0_wp_ - zzm=xv(3)/100.0_wp_ - if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then - if(psinv.ge.0.and.psinv.le.1.0_wp_.and. - . zzm.ge.zbinf.and.zzm.le.zbsup) then - call pabs_curr(p0,sox,ak0,i,j,k) - iiv(j,k)=i - else - if(iiv(j,k).gt.1) tauv(j,k,i)=tauv(j,k,i-1) - end if - if(tauv(j,k,i).le.taumn) taumn=tauv(j,k,i) - if(tauv(j,k,i).ge.taumx) taumx=tauv(j,k,i) - pabstot=pabstot+ppabs(j,k,iiv(j,k)) - currtot=currtot+ccci(j,k,iiv(j,k)) - end if - call print_output(p0,i,j,k) - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - ins_pl=inside_plasma(rrm,zzm) - if (mod(iop(j,k),2).eq.0 .and. ins_pl) then - iop(j,k)=iop(j,k)+1 - call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) - - if (ipass.gt.1 .and. index_rt.eq.1 .and. - . iowmax.gt.1 .and. istore(j,k).eq.0) then - istore(j,k)=istore(j,k)+1 - yyrfl(j,k,1:3)=xv - yyrfl(j,k,4:6)=anv - ihcd(j,k)=0 - end if - else if (mod(iop(j,k),2).eq.1.and. - . .not.ins_pl) then - iop(j,k)=iop(j,k)+1 - call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) - end if - - if (ipass.gt.1) then - if (iow(j,k).eq.0 .and. inside(rlim,zlim,nlim,rrm,zzm)) then - iow(j,k)=1 - else if (iow(j,k).eq.1 .and. - . .not.inside(rlim,zlim,nlim,rrm,zzm)) then - iow(j,k)=2 - if (ins_pl) then - iop(j,k)=iop(j,k)+1 - call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) - end if - call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)), - . eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl) - istore(j,k)=istore(j,k)+1 - yyrfl(j,k,1:3)=xvrfl - yyrfl(j,k,4:6)=anvrfl - tau1v(j,k)=tauv(j,k,iiv(j,k)) - ext(j,k,iop(j,k))=extr - eyt(j,k,iop(j,k))=eytr - if (j.lt.jclosest) then - jclosest=j - anwcl=anw - xwcl=xvrfl - end if - xv=xvrfl - anv=anvrfl - rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2) - zzm=1.0e-2_wp_*xv(3) - ywrk(1:3,j,k)=xv - ywrk(4:6,j,k)=anv - igrad=0 - call gwork(sox,xgcn,bres,j,k) - if (ins_pl) then - iop(j,k)=iop(j,k)+1 - call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) - if (index_rt.eq.1) ihcd(j,k)=0 - end if - end if - end if - - if(index_rt.eq.1 .and. j.eq.1) psinv11=psinv - if(iop(j,k).lt.iopmin) iopmin=iop(j,k) - if(iow(j,k).lt.iowmin) iowmin=iow(j,k) - if(iow(j,k).gt.iowmax) iowmax=iow(j,k) - - xvjk(:,j,k)=xv - anvjk(:,j,k)=anv - - end do - end do - if(jclosest.le.nrayr) then - aknmin=1.0_wp_ - do j=1,nrayr - kkk=nrayth - if(j.eq.1) kkk=1 - do k=1,kkk - print*,i,j,k - print*,anwcl,xwcl,anvjk(1:2,j,k) - anwclr=(anwcl(1)*xwcl(1)+anwcl(2)*xwcl(2)) - . /sqrt(xwcl(1)**2+xwcl(2)**2) - anvjkr=(anvjk(1,j,k)*xvjk(1,j,k)+anvjk(2,j,k)*xvjk(2,j,k)) - . /sqrt(xvjk(1,j,k)**2+xvjk(2,j,k)**2) - akdotn=anwclr*anvjkr+anwcl(3)*anvjk(3,j,k) - if(akdotn.lt.aknmin) aknmin=akdotn - end do - end do - else - aknmin=-1.0_wp_ - end if - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c print ray positions for j=nrayr in local reference system - - istpr=istpr+1 - if (istpr.eq.istpr0) then -c print*,'istep = ',i - if(nrayr.gt.1) then - iproj=0 - nfilp=8 - call projxyzt(iproj,nfilp) - end if - istpr=0 - end if -c - if (istpl.eq.istpl0) istpl=0 - istpl=istpl+1 - - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -c single pass is stopped when all the rays have crossed the plasma -c or complete absorption has occurred -c same for successive passes of multi-pass simulations (here exit -c from vessel is detected too -c first pass in multi-pass simulation is stopped when at least one -c ray has reflected and all rays are directed away from -c reflection point, or when no reflection has occurred and -c central ray re-enters the plasma - - if((ipass.eq.1 .and. ((iopmin.gt.1) .or. - . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr))) - . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or. - . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr)))) then - istop=1 - else if(ipass.gt.1 .and. index_rt.eq.1 .and. - . ((iowmin.gt.1 .and. aknmin.gt.0) .or. - . (iowmax.le.1 .and. iop(1,1).gt.2))) then -c flag second pass mode coupling as unset - powrfl=-1.0_wp_ - qqout=0.0_wp_ - uuout=0.0_wp_ - vvout=0.0_wp_ - do j=1,nrayr - kkk=nrayth - if(j.eq.1) kkk=1 - do k=1,kkk -c store missing initial conditions for the second pass - if (istore(j,k).eq.0) then - istore(j,k)=istore(j,k)+1 - yyrfl(j,k,1:3)=xvjk(:,j,k) - yyrfl(j,k,4:6)=anvjk(:,j,k) - tau1v(j,k)=tauv(j,k,iiv(j,k)) - end if -c determine mode coupling at the plasma boundary - if (powrfl.lt.0.0_wp_) then - call vacuum_rt(xvjk(:,j,k),anvjk(:,j,k),xvvac,ivac) -c look for first ray hitting the plasma, starting from the central -c and evaluate polarization - if (ivac.eq.1) then - y(1:3)=xvjk(:,j,k) - y(4:6)=anvjk(:,j,k) - call fwork(sox,xgcn,bres,y,dery) - call pol_limit(sox,exin2,eyin2) - call stokes(exin2,eyin2,qqin2,uuin2,vvin2) - powloop: do j1=1,nrayr - kkkk=nrayth - if(j1.eq.1) kkkk=1 - do k1=1,kkkk -c look for first ray which completed the first pass in the plasma - if (iop(j1,k1).gt.1) then -c if found, use its polarization state to compute mode coupling - call stokes(ext(j1,k1,2),eyt(j1,k1,2), - . qqout,uuout,vvout) - exit powloop - end if - end do - end do powloop -c if no ray completed a first pass in the plasma, use central ray -c initial polarization (possibly reflected) - if (qqout.le.0.0_wp_) then - call stokes(ext(1,1,0),eyt(1,1,0),qqout,uuout,vvout) - end if - powrfl=0.5_wp_*(1.0_wp_+vvout*vvin2+ - . uuout*uuin2+qqout*qqin2) - end if - end if - end do - end do - strfl11=i*dst - write(6,*) ' ' - write(6,*) 'Reflected power fraction =',powrfl - write(66,*) psipol,chipol,powrfl - istop=1 - end if - - return - end -c -c -c - subroutine print_output(p0mw,i,j,k) - use const_and_precisions, only : wp_,pi - use gray_params, only : istpl0 - use coreprofiles, only : psdbnd - use beamdata, only : ywrk,psjki,tauv,alphav,pdjki, - . currj,didst,q,tau1v,nrayr !,nrayth - use equilibrium, only : frhotor - implicit none -c arguments - real(wp_), intent(in) :: p0mw -c local constants - integer, parameter :: ndim=6 - real(wp_), parameter :: taucr=12.0_wp_ -c arguments - integer :: i,j,k -c local variables - real(wp_) :: x,y,z,rr,rrm,zzm,stm,xxm,yym,anx,any,anz,anr,anphi, - . phi,phideg,psi,rhot,bbr,bbz,bpol,dens11,dids11,dpdv11,ajcd11, - . alpha11,taujk,tau11,pt11 -c common/external functions/variables - integer :: nharm,nhf,iohkawa,index_rt,istpr,istpl - real(wp_) :: st,brr,bphi,bzz,ajphi,btot,xg,yg,dens,ddens, - . tekev,alpha,effjcd,akim,tau0,anpl,anpr,ddr,an2s,an2,fdia,bdotgr, - . ddi,ddr11,anprr,anpri - complex(wp_) :: ex,ey,ez -c - common/nharm/nharm,nhf - common/iokh/iohkawa - common/index_rt/index_rt - common/ss/st - common/istgr/istpr,istpl - common/parpl/brr,bphi,bzz,ajphi - common/btot/btot - common/xgxg/xg - common/ygyg/yg - common/dens/dens,ddens - common/tete/tekev - common/absor/alpha,effjcd,akim,tau0 - common/epolar/ex,ey,ez - common/nplr/anpl,anpr - common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 - common/nprw/anprr,anpri -c - x=ywrk(1,j,k) - y=ywrk(2,j,k) - z=ywrk(3,j,k) - rr=sqrt(x*x+y*y) -c - anx=ywrk(4,j,k) - any=ywrk(5,j,k) - anz=ywrk(6,j,k) - anr=(anx*x+any*y)/rr - anphi=(any*x-anx*y)/rr -c cnphi=(any*x-anx*y) -c - rrm=rr*1.0e-2_wp_ - zzm=z*1.0e-2_wp_ - stm=st*1.0e-2_wp_ - xxm=x*1.0e-2_wp_ - yym=y*1.0e-2_wp_ -c - if(index_rt.gt.1) then - taujk=tauv(j,k,i)+tau1v(j,k) - else - taujk=tauv(j,k,i) - end if - -c central ray only begin - if(j.eq.1) then - phi=acos(x/sqrt(y*y+x*x)) - if(y.lt.0.0_wp_) phi=-phi - phideg=phi*180.0_wp_/pi - psi=psjki(j,k,i) - rhot=1.0_wp_ - bbr=0.0_wp_ - bbz=0.0_wp_ - bpol=0.0_wp_ - rhot=1.0_wp_ - dens11=0.0_wp_ - if(psi.ge.0.0_wp_) then - if (psi.le.1.0_wp_) rhot=frhotor(psi) - bbr=brr - bbz=bzz - bpol=sqrt(bbr**2+bbz**2) - else - tekev=0.0_wp_ - akim=0.0_wp_ - end if - ddr11=ddr -c cutoff=xg-(1-yg)*(1-anpl**2) - -c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 - - if(psi.le.psdbnd.and.psi.ge.0.0_wp_) dens11=dens - dids11=didst(j,k,i)*1.0e2_wp_/(p0mw*q(j)) - dpdv11=pdjki(j,k,i)/(p0mw*q(j)) - ajcd11=currj(j,k,i)/(p0mw*q(j)) - alpha11=alphav(j,k,i) - tau11=taujk - pt11=exp(-taujk) - - write(4,99) stm,rrm,zzm,phideg, - . psi,rhot,dens11,tekev, - . bbr,bphi,bbz,ajphi*1.0e-6_wp_,sqrt(an2),anpl,akim*100, - . alpha11,tau11,pt11,ajcd11,dids11, - . dble(nhf),dble(iohkawa),dble(index_rt) - -c call polarcold(sox,exf,eyif,ezf,elf,etf) - end if -c central ray only end - - if(k.eq.1.and.j.eq.nrayr) write(17,99) st,ddr11,ddr,ddi - -c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3 - if(istpl.eq.istpl0) then - if(j.eq.nrayr) then -c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then - write(33,111) i,j,k,stm,xxm,yym,rrm,zzm, - . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) - end if -c if(k.eq.nrayth) write(33,*) ' ' - end if -c - return - 99 format(30(1x,e16.8e3)) -111 format(3i5,16(1x,e16.8e3)) - end -c -c -c - subroutine prfile - implicit none -c common/external functions/variables - integer :: index_rt -c - write(4,*)' #sst R z phi psi rhot ne Te BR Bphi Bz Jphi '// - . 'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' - write(8,*) ' #istep j k xt yt zt rt psin' - write(9,*) ' #istep j k xt yt zt rt psin' - write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' - write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' - write(12,*) ' #i sst psi w1 w2' - write(7,*)'#Icd Pa Jphip dPdVp '// - . 'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// - . 'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj '// - . 'drhotp' - write(48,*) '#rhop rhot Jphi Jcdb dPdV Icdins Pins' - write(66,*) "# psipol0 chipol0 powrfl" -c - return - end -c -c -c - subroutine bres_anal(bres) - use const_and_precisions, only : wp_,pi - use equilibrium, only : btaxis,aminor,rmaxis,zmaxis - implicit none -c arguments - real(wp_) :: bres -c local variables - integer :: i - real(wp_) :: rres,zmx,zmn,zzres -c -c print resonant B iequil=1 - write(70,*)'#i Btot R z' - rres=btaxis*rmaxis/bres - zmx=zmaxis+aminor - zmn=zmaxis-aminor - do i=1,21 - zzres=zmn+(i-1)*(zmx-zmn)/2.0e1_wp_ - write(70,111) i,bres,rres,zzres - end do - - return -111 format(i4,20(1x,e16.8e3)) - end -c -c -c - subroutine print_prof - use const_and_precisions, only : wp_ - use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi - use coreprofiles, only : density, temp - implicit none -c local constants - real(wp_), parameter :: eps=1.e-4_wp_ -c local variables - integer :: i - real(wp_) :: psin,rhop,rhot,ajphi,te,qq - real(wp_) :: dens,ddens -c - write(55,*) ' #psi rhot ne Te q Jphi' - do i=1,nq - psin=psinr(i) - rhop=sqrt(psin) -c - call density(psin,dens,ddens) - te=temp(psin) - qq=fq(psin) - rhot=frhotor(rhop) - call tor_curr_psi(max(eps,psin),ajphi) - write(55,"(12(1x,e12.5))") psin,rhot,dens,te,qq,ajphi*1.e-6_wp_ - end do -c - return - end -c -c -c - subroutine print_prof_an - use const_and_precisions, only : wp_ - use coreprofiles, only : density, temp - use equilibrium, only : frhotor - implicit none -c local constants - integer, parameter :: nst=51 -c local variables - integer :: i - real(wp_) :: psin,rhop,rhot,te - real(wp_) :: dens,ddens -c - write(55,*) ' #psi rhot ne Te' - do i=1,nst - psin=dble(i-1)/dble(nst-1) - rhop=sqrt(psin) - rhot=frhotor(rhop) - call density(psin,dens,ddens) - te=temp(psin) - write(55,111) psin,rhot,dens,te - end do -c - return - 111 format(12(1x,e12.5)) - end -c -c -c - subroutine surfq(psinq,qpsi,nq,qval) - use const_and_precisions, only : wp_ - use equilibrium, only : rmaxis,zmaxis,zbinf,zbsup,frhotor - use magsurf_data, only : npoints - use utils, only : locate, intlin - implicit none -c arguments - integer, intent(in) :: nq - real(wp_), dimension(nq), intent(in) :: psinq,qpsi - real(wp_) :: qval -c local variables - integer :: ncnt,i1,ipr - real(wp_) :: rup,zup,rlw,zlw,rhot,psival - real(wp_), dimension(npoints) :: rcn,zcn -c - ncnt=(npoints-1)/2 -c -c locate psi surface for q=qval -c - call locate(abs(qpsi),nq,qval,i1) - if (i1>0.and.i1cceq,tr,tz,kspl,points_tgo - use dierckx, only : profil,sproota - implicit none -c local constants - integer, parameter :: mest=4 -c arguments - integer, intent(in) :: ipr - real(wp_), intent(in) :: h - real(wp_), intent(inout) :: rup,zup,rlw,zlw - real(wp_), dimension(npoints), intent(out) :: rcn,zcn -c local variables - integer :: np,info,ic,ier,ii,iopt,m - real(wp_) :: ra,rb,za,zb,th,zc,val - real(wp_), dimension(mest) :: zeroc - real(wp_), dimension(nsr) :: czc -c - np=(npoints-1)/2 -c - ra=rup - rb=rlw - za=zup - zb=zlw - call points_tgo(ra,za,rup,zup,h,info) - call points_tgo(rb,zb,rlw,zlw,h,info) -c - th=pi/dble(np) - rcn(1)=rlw - zcn(1)=zlw - rcn(npoints)=rlw - zcn(npoints)=zlw - rcn(np+1)=rup - zcn(np+1)=zup - do ic=2,np - zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ - iopt=1 - call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) - if(ier.gt.0) print*,' profil =',ier - val=h*psiant+psinop - call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) - if (zeroc(1).gt.rwallm) then - rcn(ic)=zeroc(1) - rcn(npoints+1-ic)=zeroc(2) - else - rcn(ic)=zeroc(2) - rcn(npoints+1-ic)=zeroc(3) - end if - zcn(ic)=zc - zcn(npoints+1-ic)=zc - end do - if (ipr.gt.0) then - do ii=1,npoints - write(71,111) ii,h,rcn(ii),zcn(ii) - end do - write(71,*) - write(71,*) - end if -c - return -111 format(i6,12(1x,e12.5)) - end -c -c -c -C subroutine flux_average -C use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv -C use magsurf_data -C use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge, -C . equinum_fpol,bfield,fq,frhotor,zbsup,zbinf -C use simplespline, only : difcs -C use dierckx, only : regrid,coeff_parder -C implicit none -Cc local constants -C integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, -C . njest=nnintp+ksp+1,nlest=nlam+ksp+1, -C . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, -C . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam -Cc local variables -C integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr -C integer, dimension(kwrk) :: iwrk -C real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, -C . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height,height2,r2iav,currp, -C . area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, -C . dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, -C . shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, -C . bphi,brr,bzz,riav,fp,rup,rlw,zup,zlw -C real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl -C real(wp_), dimension(2*ncnt) :: dlpv -C real(wp_), dimension(2*ncnt+1) :: bv,bpv -C real(wp_), dimension(nlam) :: alam,weights -C real(wp_), dimension(nnintp,nlam) :: fhlam -C real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam -C real(wp_), dimension(lwrk) :: wrk -C real(wp_), dimension(:), allocatable :: rctemp,zctemp -Cc common/external functions/variables -C real(wp_) :: fpolv -Cc -C npsi=nnintp -C ninpr=(npsi-1)/10 -C npoints = 2*ncnt+1 -Cc -C call alloc_surfvec(ierr) -C if(allocated(tjp)) deallocate(tjp) -C if(allocated(tlm)) deallocate(tlm) -C if(allocated(ch)) deallocate(ch) -C allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), -C . rctemp(npoints),zctemp(npoints),stat=ierr) -C if (ierr.ne.0) return -C -Cc computation of flux surface averaged quantities -C -C write(71,*)' #i psin R z' -C -C dlam=1.0_wp_/dble(nlam-1) -C do l=1,nlam-1 -C alam(l)=dble(l-1)*dlam -C fhlam(1,l)=sqrt(1.0_wp_-alam(l)) -C ffhlam(l)=fhlam(1,l) -C dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) -C weights(l)=1.0_wp_ -C end do -C weights(1)=0.5_wp_ -C weights(nlam)=0.5_wp_ -C alam(nlam)=1.0_wp_ -C fhlam(1,nlam)=0.0_wp_ -C ffhlam(nlam)=0.0_wp_ -C dffhlam(nlam)=-99999.0_wp_ -C -C jp=1 -C anorm=2.0_wp_*pi*rmaxis/abs(btaxis) -C b2av=btaxis**2 -C dvdpsi=2.0_wp_*pi*anorm -C dadpsi=2.0_wp_*pi/abs(btaxis) -C ratio_cdator=abs(btaxis/btrcen) -C ratio_cdbtor=1.0_wp_ -C ratio_pltor=1.0_wp_ -C qq=1.0_wp_ -C fc=1.0_wp_ -C -C psicon(1)=0.0_wp_ -C rcon(1,:)=rmaxis -C zcon(1,:)=zmaxis -C pstab(1)=0.0_wp_ -C rhot_eq(1)=0.0_wp_ -C rpstab(1)=0.0_wp_ -C rhotqv(1)=0.0_wp_ -C vcurrp(1)=0.0_wp_ -C vajphiav(1)=0.0_wp_ -C bmxpsi(1)=abs(btaxis) -C bmnpsi(1)=abs(btaxis) -C bav(1)=abs(btaxis) -C rbav(1)=1.0_wp_ -C rri(1)=rmaxis -C varea(1)=0.0_wp_ -C vvol(1)=0.0_wp_ -C vratjpl(1)=ratio_pltor -C vratja(1)=ratio_cdator -C vratjb(1)=ratio_cdbtor -C ffc(1)=fc -C rhot2q=0.0_wp_ -C qqv(1)=1.0_wp_ -C -C dadrhotv(1)=0.0_wp_ -C dvdrhotv(1)=0.0_wp_ -C -C rup=rmaxis -C rlw=rmaxis -C zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ -C zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ -C -C do jp=2,npsi -C height=dble(jp-1)/dble(npsi-1) -C psicon(jp)=height -C if(jp.eq.npsi) height=0.9999_wp_ -C ipr=0 -C jpr=mod(jp,ninpr) -C if(jpr.eq.1) ipr=1 -C height2=height*height -C call contours_psi(height2,rup,zup,rlw,zlw,rctemp,zctemp,ipr) -C rcon(jp,:) = rctemp -C zcon(jp,:) = zctemp -Cc -C r2iav=0.0_wp_ -C anorm=0.0_wp_ -C dadpsi=0.0_wp_ -C currp=0.0_wp_ -C b2av=0.0_wp_ -C area=0.0_wp_ -C volume=0.0_wp_ -C ajphiav=0.0_wp_ -C bbav=0.0_wp_ -C bmmx=-1.0e+30_wp_ -C bmmn=1.0e+30_wp_ -C -C call tor_curr(rctemp(1),zctemp(1),ajphi0) -C call bfield(rctemp(1),zctemp(1),br=brr,bz=bzz) -C call equinum_fpol(height2,fpolv) -C bphi=fpolv/rctemp(1) -C btot0=sqrt(bphi**2+brr**2+bzz**2) -C bpoloid0=sqrt(brr**2+bzz**2) -C bv(1)=btot0 -C bpv(1)=bpoloid0 -C rpsim0=rctemp(1) -C -C do inc=1,npoints-1 -C inc1=inc+1 -C dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) -C dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) -C dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ -C . (zctemp(inc1)-zctemp(inc))**2) -C drc=(rctemp(inc1)-rctemp(inc)) -C -Cc compute length, area and volume defined by psi=height^2 -C -C ph=0.5_wp_*(dla+dlb+dlp) -C area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) -C area=area+sqrt(area2) -C rzp=rctemp(inc1)*zctemp(inc1) -C rz=rctemp(inc)*zctemp(inc) -C volume=pi*(rzp+rz)*drc+volume -C -Cc compute line integral on the contour psi=height^2 -C -C rpsim=rctemp(inc1) -C zpsim=zctemp(inc1) -C call bfield(rpsim,zpsim,br=brr,bz=bzz) -C call tor_curr(rpsim,zpsim,ajphi) -C! call equinum_fpol(height2,fpolv) -C bphi=fpolv/rpsim -C btot=sqrt(bphi**2+brr**2+bzz**2) -C bpoloid=sqrt(brr**2+bzz**2) -C dlpv(inc)=dlp -C bv(inc1)=btot -C bpv(inc1)=bpoloid -C -C dlph=0.5_wp_*dlp -C anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) -C dadpsi=dadpsi+dlph* -C . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) -C currp=currp+dlph*(bpoloid+bpoloid0) -C b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) -C bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) -C r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) -C . +1.0_wp_/(bpoloid0*rpsim0**2)) -C ajphiav=ajphiav+dlph* -C . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) -C -C ajphi0=ajphi -C rpsim0=rpsim -C bpoloid0=bpoloid -C btot0=btot -C -Cc computation maximum/minimum B values on given flux surface -C -C if(btot.le.bmmn) bmmn=btot -C if(btot.ge.bmmx) bmmx=btot -C end do -C -Cc bav= [T] , b2av= [T^2] , rbav=/b_min -Cc anorm = int d l_p/B_p = dV/dpsi/(2pi) -Cc r2iav=<1/R^2> [m^-2] , -Cc riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), -Cc rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] -Cc currp = plasma current within psi=const -C -C bbav=bbav/anorm -C r2iav=r2iav/anorm -C dvdpsi=2.0_wp_*pi*anorm -C riav=dadpsi/anorm -C b2av=b2av/anorm -C vcurrp(jp)=ccj*currp -C vajphiav(jp)=ajphiav/dadpsi -C -Cc area == varea, volume == vvol -Cc flux surface minor radius == (area/pi)^1/2 -Cc ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 -Cc ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / -Cc ratio_pltor = Jcd_||/J_phi Jcd_|| = -C -C pstab(jp)=height2 -C rpstab(jp)=height -C vvol(jp)=abs(volume) -C varea(jp)=area -C bav(jp)=bbav -C rbav(jp)=bbav/bmmn -C bmxpsi(jp)=bmmx -C bmnpsi(jp)=bmmn -C rri(jp)=bav(jp)/abs(fpolv*r2iav) -C ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) -C ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) -C ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) -C vratjpl(jp)=ratio_pltor -C vratja(jp)=ratio_cdator -C vratjb(jp)=ratio_cdbtor -C qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) -C qqv(jp)=qq -C -C dadrhotv(jp)=phitedge*frhotor(height)/fq(height2) -C . *dadpsi/pi -C dvdrhotv(jp)=phitedge*frhotor(height)/fq(height2) -C . *dvdpsi/pi -C -Cc computation of rhot from calculated q profile -C -C rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) -C . /dble(npsi-1) -C rhotqv(jp)=sqrt(rhot2q) -C -Cc computation of fraction of circulating/trapped fraction fc, ft -Cc and of function H(lambda,rhop) -Cc ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ -C -C fc=0.0_wp_ -C shlam=0.0_wp_ -C do l=nlam,1,-1 -C lam=alam(l) -C srl=0.0_wp_ -C rl2=1.0_wp_-lam*bv(1)/bmmx -C rl0=0.0_wp_ -C if(rl2.gt.0) rl0=sqrt(rl2) -C do inc=1,npoints-1 -C rl2=1.0_wp_-lam*bv(inc+1)/bmmx -C rl=0.0_wp_ -C if(rl2.gt.0) rl=sqrt(rl2) -C srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) -C rl0=rl -C end do -C srl=srl/anorm -C dhlam=0.5_wp_/srl -C fc=fc+lam/srl*weights(l) -C if(l.eq.nlam) then -C fhlam(jp,l)=0.0_wp_ -C ffhlam(nlam*(jp-1)+l)=0.0_wp_ -C dffhlam(nlam*(jp-1)+l)=-dhlam -C dhlam0=dhlam -C else -C shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam -C fhlam(jp,l)=shlam -C dffhlam(nlam*(jp-1)+l)=-dhlam -C dhlam0=dhlam -C end if -C end do -C fc=0.75_wp_*b2av/bmmx**2*fc*dlam -C ffc(jp)=fc -C -C ccfh=bmmn/bmmx/fc -C do l=1,nlam -C ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) -C dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) -C end do -C end do -C -C write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// -C .' Area Vol |I_pl| qq fc ratioJa ratioJb' -C -C qqv(1)=qqv(2) -C vajphiav(1)=vajphiav(2) -C do jp=1,npsi -C rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) -C if(jp.eq.npsi) then -C rhotqv(jp)=1.0_wp_ -C rpstab(jp)=1.0_wp_ -C pstab(jp)=1.0_wp_ -C end if -C rhot_eq(jp)=frhotor(sqrt(pstab(jp))) -C write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), -C . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), -C . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), -C . vratja(jp),vratjb(jp) -C end do -C -Cc rarea=sqrt(varea(npsi)/pi) -C -Cc spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi -Cc used for computations of dP/dV and J_cd -Cc spline coefficients of rhot -C -C iopt=0 -C call difcs(rpstab,vvol,npsi,iopt,cvol,ier) -C iopt=0 -C call difcs(rpstab,rbav,npsi,iopt,crbav,ier) -C iopt=0 -C call difcs(rpstab,rri,npsi,iopt,crri,ier) -C iopt=0 -C call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) -C iopt=0 -C call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) -C iopt=0 -C call difcs(rpstab,vratja,npsi,iopt,cratja,ier) -C iopt=0 -C call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) -C iopt=0 -C call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) -C iopt=0 -C call difcs(rpstab,varea,npsi,iopt,carea,ier) -C iopt=0 -C call difcs(rpstab,ffc,npsi,iopt,cfc,ier) -C iopt=0 -C call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) -C iopt=0 -C call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) -C iopt=0 -C call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) -C -Cc spline interpolation of H(lambda,rhop) and dH/dlambda -C -C iopt=0 -C s=0.0_wp_ -C call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, -C . zero,one,zero,one,ksp,ksp,s, -C . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) -C njpt=njp -C nlmt=nlm -C -C return -C 99 format(20(1x,e12.5)) -C end -Cc -Cc -Cc -C function fdadrhot(rpsi) -C use const_and_precisions, only : wp_ -C use magsurf_data, only : rpstab,cdadrhot,npsi -C use simplespline, only :spli -C implicit none -Cc arguments -C real(wp_) :: rpsi,fdadrhot -Cc local variables -C integer :: ip -C real(wp_) :: dps -Cc -C ip=int((npsi-1)*rpsi+1) -Cc if(ip.eq.0) ip=1 -Cc if(ip.eq.npsi) ip=npsi-1 -C ip=min(max(1,ip),npsi-1) -C dps=rpsi-rpstab(ip) -C fdadrhot=spli(cdadrhot,npsi,ip,dps) -C return -C end -Cc -Cc -Cc -C function fdvdrhot(rpsi) -C use const_and_precisions, only : wp_ -C use magsurf_data, only : rpstab,cdvdrhot,npsi -C use simplespline, only :spli -C implicit none -Cc arguments -C real(wp_) :: rpsi,fdvdrhot -Cc local variables -C integer :: ip -C real(wp_) :: dps -Cc -C ip=int((npsi-1)*rpsi+1) -C ip=min(max(1,ip),npsi-1) -C dps=rpsi-rpstab(ip) -C fdvdrhot=spli(cdvdrhot,npsi,ip,dps) -Cc -C return -C end -c -c -c -C subroutine flux_average_an -C use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv -C use magsurf_data -C use equilibrium, only : btrcen,frhotor,btaxis,rmaxis,zmaxis, -C . fq,aminor,equian -C use simplespline, only : difcs -C use dierckx, only : regrid,coeff_parder -C implicit none -Cc local constants -C integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, -C . njest=nnintp+ksp+1,nlest=nlam+ksp+1, -C . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, -C . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam -Cc local variables -C integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr -C integer, dimension(kwrk) :: iwrk -C real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, -C . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height2,r2iav,currp,area, -C . volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb,dlp, -C . drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0,shlam, -C . srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi,bphi,brr,bzz,riav,fp, -C . dadr,dbvcdc13,dbvcdc31,dvdr,height,phitedge,qqan,rn -C real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl -C real(wp_), dimension(2*ncnt) :: dlpv -C real(wp_), dimension(2*ncnt+1) :: bv,bpv -C real(wp_), dimension(nlam) :: alam,weights -C real(wp_), dimension(nnintp,nlam) :: fhlam -C real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam -C real(wp_), dimension(lwrk) :: wrk -C real(wp_), dimension(:), allocatable :: rctemp,zctemp -C real(wp_) :: psinv !unused. can be removed when equian is modified -C !to accept optional arguments -Cc common/external functions/variables -C real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv -Cc -C common/fpol/fpolv,dfpv -C common/derip1/dpsidr,dpsidz -C common/derip2/ddpsidrr,ddpsidzz,ddpsidrz -Cc -C npsi=nnintp -C ninpr=(npsi-1)/10 -C npoints = 2*ncnt+1 -Cc -C call alloc_surfvec(ierr) -C if(allocated(tjp)) deallocate(tjp) -C if(allocated(tlm)) deallocate(tlm) -C if(allocated(ch)) deallocate(ch) -C allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), -C . rctemp(npoints),zctemp(npoints),stat=ierr) -C if (ierr.ne.0) return -Cc -Cc computation of flux surface averaged quantities -C -C phitedge=pi*aminor**2*btaxis -C -C write(71,*)' #i psin R z' -C -C dlam=1.0_wp_/dble(nlam-1) -C do l=1,nlam-1 -C alam(l)=dble(l-1)*dlam -C fhlam(1,l)=sqrt(1.0_wp_-alam(l)) -C ffhlam(l)=fhlam(1,l) -C dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) -C weights(l)=1.0_wp_ -C end do -C weights(1)=0.5_wp_ -C weights(nlam)=0.5_wp_ -C alam(nlam)=1.0_wp_ -C fhlam(1,nlam)=0.0_wp_ -C ffhlam(nlam)=0.0_wp_ -C dffhlam(nlam)=-99999.0_wp_ -C -C jp=1 -C anorm=2.0_wp_*pi*rmaxis/abs(btaxis) -C b2av=btaxis**2 -C dvdpsi=2.0_wp_*pi*anorm -C dadpsi=2.0_wp_*pi/abs(btaxis) -C ratio_cdator=abs(btaxis/btrcen) -C ratio_cdbtor=1.0_wp_ -C ratio_pltor=1.0_wp_ -C qq=1.0_wp_ -C fc=1.0_wp_ -C -C psicon(1)=0.0_wp_ -C rcon(1,:)=rmaxis -C zcon(1,:)=zmaxis -C pstab(1)=0.0_wp_ -C rhot_eq(1)=0.0_wp_ -C rpstab(1)=0.0_wp_ -C rhotqv(1)=0.0_wp_ -C vcurrp(1)=0.0_wp_ -C vajphiav(1)=0.0_wp_ -C bmxpsi(1)=abs(btaxis) -C bmnpsi(1)=abs(btaxis) -C bav(1)=abs(btaxis) -C rbav(1)=1.0_wp_ -C rri(1)=rmaxis -C varea(1)=0.0_wp_ -C vvol(1)=0.0_wp_ -C vratjpl(1)=ratio_pltor -C vratja(1)=ratio_cdator -C vratjb(1)=ratio_cdbtor -C ffc(1)=fc -C rhot2q=0.0_wp_ -C dadrhotv(1)=0.0_wp_ -C dvdrhotv(1)=0.0_wp_ -C qqv(1)=1.0_wp_ -C -Cc rup=rmaxis -Cc rlw=rmaxis -Cc zup=(zbmax+zmaxis)/2.0_wp_ -Cc zlw=(zmaxis+zbmin)/2.0_wp_ -C -C do jp=2,npsi -C height=dble(jp-1)/dble(npsi-1) -C psicon(jp)=height -C if(jp.eq.npsi) height=0.9999_wp_ -C height2=height*height -C ipr=0 -C jpr=mod(jp,ninpr) -C if(jpr.eq.1) ipr=1 -C call contours_psi_an(height2,rctemp,zctemp,ipr) -C rcon(jp,:) = rctemp -C zcon(jp,:) = zctemp -Cc -C r2iav=0.0_wp_ -C anorm=0.0_wp_ -C dadpsi=0.0_wp_ -C currp=0.0_wp_ -C b2av=0.0_wp_ -C area=0.0_wp_ -C volume=0.0_wp_ -C ajphiav=0.0_wp_ -C bbav=0.0_wp_ -C bmmx=-1.0e+30_wp_ -C bmmn=1.0e+30_wp_ -C -C call equian(rctemp(1),zctemp(1),psinv) -C dbvcdc13=-ddpsidzz/rctemp(1) -C dbvcdc31= ddpsidrr/rctemp(1)-bzz/rctemp(1) -C ajphi=ccj*(dbvcdc13-dbvcdc31) -C brr=-dpsidz/rctemp(1) -C bzz= dpsidr/rctemp(1) -C bphi=fpolv/rctemp(1) -C btot0=sqrt(bphi**2+brr**2+bzz**2) -C bpoloid0=sqrt(brr**2+bzz**2) -C bv(1)=btot0 -C bpv(1)=bpoloid0 -C rpsim0=rctemp(1) -C -C do inc=1,npoints-1 -C inc1=inc+1 -C dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) -C dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) -C dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ -C . (zctemp(inc1)-zctemp(inc))**2) -C drc=(rctemp(inc1)-rctemp(inc)) -C -Cc compute length, area and volume defined by psi=height^2 -C -C ph=0.5_wp_*(dla+dlb+dlp) -C area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) -C area=area+sqrt(area2) -C rzp=rctemp(inc1)*zctemp(inc1) -C rz=rctemp(inc)*zctemp(inc) -C volume=pi*(rzp+rz)*drc+volume -C -Cc compute line integral on the contour psi=height^2 -C -C rpsim=rctemp(inc1) -C zpsim=zctemp(inc1) -C call equian(rpsim,zpsim,psinv) -C brr=-dpsidz/rpsim -C bzz= dpsidr/rpsim -C dbvcdc13=-ddpsidzz/rpsim -C dbvcdc31= ddpsidrr/rpsim-bzz/rpsim -C ajphi=ccj*(dbvcdc13-dbvcdc31) -C bphi=fpolv/rpsim -C btot=sqrt(bphi**2+brr**2+bzz**2) -C bpoloid=sqrt(brr**2+bzz**2) -C dlpv(inc)=dlp -C bv(inc1)=btot -C bpv(inc1)=bpoloid -C -C dlph=0.5_wp_*dlp -C anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) -C dadpsi=dadpsi+dlph* -C . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) -C currp=currp+dlph*(bpoloid+bpoloid0) -C b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) -C bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) -C r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) -C . +1.0_wp_/(bpoloid0*rpsim0**2)) -C ajphiav=ajphiav+dlph* -C . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) -C -C ajphi0=ajphi -C rpsim0=rpsim -C bpoloid0=bpoloid -C btot0=btot -Cc computation maximum/minimum B values on given flux surface -C -C if(btot.le.bmmn) bmmn=btot -C if(btot.ge.bmmx) bmmx=btot -C end do -C -Cc bav= [T] , b2av= [T^2] , rbav=/b_min -Cc anorm = int d l_p/B_p = dV/dpsi/(2pi) -Cc r2iav=<1/R^2> [m^-2] , -Cc riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), -Cc rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] -Cc currp = plasma current within psi=const -C -C bbav=bbav/anorm -C r2iav=r2iav/anorm -C dvdpsi=2.0_wp_*pi*anorm -C riav=dadpsi/anorm -C b2av=b2av/anorm -C vcurrp(jp)=ccj*currp -C vajphiav(jp)=ajphiav/dadpsi -C -Cc area == varea, volume == vvol -Cc flux surface minor radius == (area/pi)^1/2 -Cc ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 -Cc ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / -Cc ratio_pltor = Jcd_||/J_phi Jcd_|| = -C -C pstab(jp)=height2 -C rpstab(jp)=height -C vvol(jp)=abs(volume) -C varea(jp)=area -C bav(jp)=bbav -C rbav(jp)=bbav/bmmn -C bmxpsi(jp)=bmmx -C bmnpsi(jp)=bmmn -C rri(jp)=bav(jp)/abs(fpolv*r2iav) -C ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) -C ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) -C ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) -C vratjpl(jp)=ratio_pltor -C vratja(jp)=ratio_cdator -C vratjb(jp)=ratio_cdbtor -C qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) -C qqv(jp)=qq -C -C rn=frhotor(sqrt(pstab(jp))) -C qqan=fq(pstab(jp)) -C -C dadr=2.0_wp_*pi*rn*aminor**2 -C dvdr=4.0_wp_*pi**2*rn*rmaxis*aminor**2 -C -Cc dadrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dadpsi -Cc dvdrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dvdpsi -C dadrhotv(jp)=phitedge*rn*dadpsi/pi/qqan -C dvdrhotv(jp)=phitedge*rn*dvdpsi/pi/qqan -C -Cc computation of rhot from calculated q profile -C rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) -C . /dble(npsi-1) -Cc print*,jp,rhot2q,qqv(jp),rpstab(jp),qqv(jp-1),rpstab(jp-1) -C rhotqv(jp)=sqrt(rhot2q) -Cc rhotqv(jp)=rn -Cc -C write(57,99) rpstab(jp),rn,rhotqv(jp),qqv(jp),qqan,r2iav, -C . dadr,dadrhotv(jp),dadpsi,dvdpsi,fpolv -C -Cc computation of fraction of circulating/trapped fraction fc, ft -Cc and of function H(lambda,rhop) -Cc ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ -C -C fc=0.0_wp_ -C shlam=0.0_wp_ -C do l=nlam,1,-1 -C lam=alam(l) -C srl=0.0_wp_ -C rl2=1.0_wp_-lam*bv(1)/bmmx -C rl0=0.0_wp_ -C if(rl2.gt.0) rl0=sqrt(rl2) -C do inc=1,npoints-1 -C rl2=1.0_wp_-lam*bv(inc+1)/bmmx -C rl=0.0_wp_ -C if(rl2.gt.0) rl=sqrt(rl2) -C srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) -C rl0=rl -C end do -C srl=srl/anorm -C dhlam=0.5_wp_/srl -C fc=fc+lam/srl*weights(l) -C if(l.eq.nlam) then -C fhlam(jp,l)=0.0_wp_ -C ffhlam(nlam*(jp-1)+l)=0.0_wp_ -C dffhlam(nlam*(jp-1)+l)=-dhlam -C dhlam0=dhlam -C else -C shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam -C fhlam(jp,l)=shlam -C dffhlam(nlam*(jp-1)+l)=-dhlam -C dhlam0=dhlam -C end if -C end do -C fc=0.75_wp_*b2av/bmmx**2*fc*dlam -C ffc(jp)=fc -C -C ccfh=bmmn/bmmx/fc -C do l=1,nlam -C ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) -C dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) -C end do -C end do -C -C write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// -C .' Area Vol |I_pl| qq fc ratioJa ratioJb dadr dvdr' -C -C qqv(1)=qqv(2) -C vajphiav(1)=vajphiav(2) -C do jp=1,npsi -C rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) -C if(jp.eq.npsi) then -C rhotqv(jp)=1.0_wp_ -C rpstab(jp)=1.0_wp_ -C pstab(jp)=1.0_wp_ -C end if -C rhot_eq(jp)=frhotor(sqrt(pstab(jp))) -C write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), -C . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), -C . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), -C . vratja(jp),vratjb(jp),dadrhotv(jp),dvdrhotv(jp) -C end do -C -Cc rarea=sqrt(varea(npsi)/pi) -C -Cc spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi -Cc used for computations of dP/dV and J_cd -Cc spline coefficients of rhot -C -C iopt=0 -C call difcs(rpstab,vvol,npsi,iopt,cvol,ier) -C iopt=0 -C call difcs(rpstab,rbav,npsi,iopt,crbav,ier) -C iopt=0 -C call difcs(rpstab,rri,npsi,iopt,crri,ier) -C iopt=0 -C call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) -C iopt=0 -C call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) -C iopt=0 -C call difcs(rpstab,vratja,npsi,iopt,cratja,ier) -C iopt=0 -C call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) -C iopt=0 -C call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) -C iopt=0 -C call difcs(rpstab,varea,npsi,iopt,carea,ier) -C iopt=0 -C call difcs(rpstab,ffc,npsi,iopt,cfc,ier) -C iopt=0 -C call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) -C iopt=0 -C call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) -C iopt=0 -C call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) -C -Cc spline interpolation of H(lambda,rhop) -C -C iopt=0 -C s=0.0_wp_ -C call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, -C . zero,one,zero,one,ksp,ksp,s, -C . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) -C njpt=njp -C nlmt=nlm -C -C return -C 99 format(20(1x,e12.5)) -C end -c -c -c - subroutine contours_psi_an(h,rcn,zcn,ipr) - use const_and_precisions, only : wp_,pi - use magsurf_data, only : npoints - use equilibrium, only : frhotor,aminor,rmaxis,zmaxis - implicit none -c arguments - integer :: ipr - real(wp_) :: h - real(wp_), dimension(npoints) :: rcn,zcn -c local variables - integer :: np,ic - real(wp_) :: rn,th -c - np=(npoints-1)/2 - th=pi/dble(np) - rn=frhotor(sqrt(h)) - do ic=1,npoints - zcn(ic)=zmaxis+aminor*rn*sin(th*(ic-1)) - rcn(ic)=rmaxis+aminor*rn*cos(th*(ic-1)) -c - if (ipr.gt.0) write(71,111) ic,h,rcn(ic),zcn(ic) - end do - write(71,*) -c - return -111 format(i6,12(1x,e12.5)) - end -c -c -c - subroutine vectinit - use const_and_precisions, only : wp_ - use dispersion, only: expinit - use gray_params, only : iwarm - use beamdata, only : psjki,tauv,alphav,pdjki,ppabs, - . currj,didst,ccci,iiv,iop,iow,ihcd,istore,tau1v,nrayr,nrayth, - . nstep - implicit none -c local variables - integer :: i,j,k -c common/external functions/variables - integer :: jclosest - real(wp_), dimension(3) :: anwcl,xwcl -c - common/refln/anwcl,xwcl,jclosest -c - jclosest=nrayr+1 - anwcl(1:3)=0.0_wp_ - xwcl(1:3)=0.0_wp_ -c - do i=1,nstep - do k=1,nrayth - do j=1,nrayr - psjki(j,k,i)=0.0_wp_ - tauv(j,k,i)=0.0_wp_ - alphav(j,k,i)=0.0_wp_ - pdjki(j,k,i)=0.0_wp_ - ppabs(j,k,i)=0.0_wp_ - didst(j,k,i)=0.0_wp_ - ccci(j,k,i)=0.0_wp_ - currj(j,k,i)=0.0_wp_ - iiv(j,k)=1 - iop(j,k)=0 - iow(j,k)=0 - ihcd(j,k)=1 - istore(j,k)=0 - tau1v(j,k)=0.0_wp_ - end do - end do - end do -c - if(iwarm.gt.1) call expinit -c - return - end -c -c -c - subroutine vectinit2 - use const_and_precisions, only : wp_ - use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj, - . didst,ccci,iiv,iop,iow,ihcd,nrayr,nrayth,nstep - implicit none -c local variables - integer :: i,j,k -c - do i=1,nstep - do k=1,nrayth - do j=1,nrayr - psjki(j,k,i)=0.0_wp_ - tauv(j,k,i)=0.0_wp_ - alphav(j,k,i)=0.0_wp_ - pdjki(j,k,i)=0.0_wp_ - ppabs(j,k,i)=0.0_wp_ - didst(j,k,i)=0.0_wp_ - ccci(j,k,i)=0.0_wp_ - currj(j,k,i)=0.0_wp_ - iiv(j,k)=1 - iop(j,k)=0 - iow(j,k)=0 - ihcd(j,k)=1 - end do - end do - end do - - return - end -c -c -c - subroutine paraminit - use const_and_precisions, only : wp_ - implicit none -c common/external functions/variables - integer :: istep,istpr,istpl,ierr,istop,index_rt -c - common/istep/istep - common/istgr/istpr,istpl - common/ierr/ierr - common/istop/istop - common/index_rt/index_rt -c - istpr=0 - istpl=1 - ierr=0 - istep=0 - istop=0 - index_rt=1 -c - return - end -c -c -c - subroutine updatepos - use const_and_precisions, only : wp_ - use beamdata, only : xc,xco,du1,du1o,ywrk,nrayr,nrayth - implicit none -c local variables - integer :: j,k -c - do j=1,nrayr - do k=1,nrayth - if(j.eq.1.and.k.gt.1) then - xco(1,j,k)=xco(1,j,1) - xco(2,j,k)=xco(2,j,1) - xco(3,j,k)=xco(3,j,1) - xc(1,j,k)=xc(1,j,1) - xc(2,j,k)=xc(2,j,1) - xc(3,j,k)=xc(3,j,1) - else - xco(1,j,k)=xc(1,j,k) - xco(2,j,k)=xc(2,j,k) - xco(3,j,k)=xc(3,j,k) - xc(1,j,k)=ywrk(1,j,k) - xc(2,j,k)=ywrk(2,j,k) - xc(3,j,k)=ywrk(3,j,k) - end if - du1o(1,j,k)=du1(1,j,k) - du1o(2,j,k)=du1(2,j,k) - du1o(3,j,k)=du1(3,j,k) - end do - end do -c - return - end -c -c -c - subroutine gradi - use const_and_precisions, only : wp_ - use beamdata, only : dffiu,ddffiu,xc,xco,du1,du1o,gri,ggri, - . dgrad2v,grad2,nrayr,nrayth - implicit none -c local variables - integer :: iv,j,jm,jp,k,km,kp - real(wp_) :: ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz, - . dfu,dfuu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz - real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu, - . dgg1,dgg2,dgg3,df1,df2,df3 -c -c compute grad u and grad(S_I) -c - do k=1,nrayth - do j=1,nrayr - if(j.eq.1) then - gri(1,j,k)=0.0_wp_ - gri(2,j,k)=0.0_wp_ - gri(3,j,k)=0.0_wp_ - jp=j+1 - km=k-1 - if(k.eq.1) km=nrayth - kp=k+1 - if(k.eq.nrayth) kp=1 - do iv=1,3 - dxv1(iv)=xc(iv,jp,k)-xc(iv,j,k) - dxv2(iv)=xc(iv,jp,kp)-xc(iv,jp,km) - dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) - end do - call solg0(dxv1,dxv2,dxv3,dgu) - else - jm=j-1 - km=k-1 - if(k.eq.1) km=nrayth - kp=k+1 - if(k.eq.nrayth) kp=1 - do iv=1,3 - dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) - dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) - dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) - end do - call solg0(dxv1,dxv2,dxv3,dgu) - end if - du1(1,j,k)=dgu(1) - du1(2,j,k)=dgu(2) - du1(3,j,k)=dgu(3) - gri(1,j,k)=dgu(1)*dffiu(j) - gri(2,j,k)=dgu(2)*dffiu(j) - gri(3,j,k)=dgu(3)*dffiu(j) - grad2(j,k)=gri(1,j,k)**2+gri(2,j,k)**2+gri(3,j,k)**2 - end do - end do -c -c compute derivatives of grad u and grad(S_I) -c - do k=1,nrayth - do j=1,nrayr - if(j.eq.1) then - jp=j+1 - km=k-1 - if(k.eq.1) km=nrayth - kp=k+1 - if(k.eq.nrayth) kp=1 - do iv=1,3 - dxv1(iv)=xc(iv,jp,kp)-xc(iv,jp,km) - dxv2(iv)=xc(iv,jp,k)-xc(iv,j,k) - dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) - end do - df1(1)=du1(1,jp,kp)-du1(1,jp,km) - df1(2)=du1(1,jp,k)-du1(1,j,k) - df1(3)=du1(1,j,k)-du1o(1,j,k) - df2(1)=du1(2,jp,kp)-du1(2,jp,km) - df2(2)=du1(2,jp,k)-du1(2,j,k) - df2(3)=du1(2,j,k)-du1o(2,j,k) - df3(1)=du1(3,jp,kp)-du1(3,jp,km) - df3(2)=du1(3,jp,k)-du1(3,j,k) - df3(3)=du1(3,j,k)-du1o(3,j,k) - call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3) - else - jm=j-1 - km=k-1 - if(k.eq.1) km=nrayth - kp=k+1 - if(k.eq.nrayth) kp=1 - do iv=1,3 - dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) - dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) - dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) - end do - df1(1)=du1(1,j,k)-du1(1,jm,k) - df1(2)=du1(1,j,kp)-du1(1,j,km) - df1(3)=du1(1,j,k)-du1o(1,j,k) - df2(1)=du1(2,j,k)-du1(2,jm,k) - df2(2)=du1(2,j,kp)-du1(2,j,km) - df2(3)=du1(2,j,k)-du1o(2,j,k) - df3(1)=du1(3,j,k)-du1(3,jm,k) - df3(2)=du1(3,j,kp)-du1(3,j,km) - df3(3)=du1(3,j,k)-du1o(3,j,k) - call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3) - end if -c -c derivatives of u -c - ux=du1(1,j,k) - uy=du1(2,j,k) - uz=du1(3,j,k) - uxx=dgg1(1) - uyy=dgg2(2) - uzz=dgg3(3) - uxy=(dgg1(2)+dgg2(1))/2.0_wp_ - uxz=(dgg1(3)+dgg3(1))/2.0_wp_ - uyz=(dgg2(3)+dgg3(2))/2.0_wp_ -c -c derivatives of S_I and Grad(S_I) -c - dfu=dffiu(j) - dfuu=ddffiu(j) - gx=ux*dfu - gy=uy*dfu - gz=uz*dfu - gxx=dfuu*ux*ux+dfu*uxx - gyy=dfuu*uy*uy+dfu*uyy - gzz=dfuu*uz*uz+dfu*uzz - gxy=dfuu*ux*uy+dfu*uxy - gxz=dfuu*ux*uz+dfu*uxz - gyz=dfuu*uy*uz+dfu*uyz -c - ggri(1,1,j,k)=gxx - ggri(2,2,j,k)=gyy - ggri(3,3,j,k)=gzz - ggri(1,2,j,k)=gxy - ggri(2,1,j,k)=gxy - ggri(1,3,j,k)=gxz - ggri(3,1,j,k)=gxz - ggri(2,3,j,k)=gyz - ggri(3,2,j,k)=gyz -c -c derivatives of |Grad(S_I)|^2 -c - dgrad2v(1,j,k)=2.0_wp_*(gx*gxx+gy*gxy+gz*gxz) - dgrad2v(2,j,k)=2.0_wp_*(gx*gxy+gy*gyy+gz*gyz) - dgrad2v(3,j,k)=2.0_wp_*(gx*gxz+gy*gyz+gz*gzz) - end do - end do -c - return - end -c -c -c - subroutine solg0(dxv1,dxv2,dxv3,dgg) -c solution of the linear system of 3 eqs : dgg . dxv = dff -c input vectors : dxv1, dxv2, dxv3, dff -c output vector : dgg -c dff=(1,0,0) -c - use const_and_precisions, only : wp_ - implicit none -c arguments - real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgg -c local variables - real(wp_) :: denom,aa1,aa2,aa3 -c - aa1=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) - aa2=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) - aa3=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) - denom = dxv1(1)*aa1-dxv2(1)*aa2+dxv3(1)*aa3 - dgg(1) = aa1/denom - dgg(2) = -(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))/denom - dgg(3) = (dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))/denom -c - return - end -c -c -c - subroutine solg3(dxv1,dxv2,dxv3,df1,df2,df3,dg1,dg2,dg3) -c three rhs vectors df1, df2, df3 -c - use const_and_precisions, only : wp_ - implicit none -c arguments - real(wp_), dimension(3) :: dxv1,dxv2,dxv3,df1,df2,df3, - . dg1,dg2,dg3 -c local variables - real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 -c - a11=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) - a21=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) - a31=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) - a12=(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3)) - a22=(dxv1(1)*dxv3(3)-dxv1(3)*dxv3(1)) - a32=(dxv1(1)*dxv2(3)-dxv1(3)*dxv2(1)) - a13=(dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2)) - a23=(dxv1(1)*dxv3(2)-dxv1(2)*dxv3(1)) - a33=(dxv1(1)*dxv2(2)-dxv1(2)*dxv2(1)) - denom = dxv1(1)*a11-dxv2(1)*a21+dxv3(1)*a31 - dg1(1) = (df1(1)*a11-df1(2)*a21+df1(3)*a31)/denom - dg1(2) = -(df1(1)*a12-df1(2)*a22+df1(3)*a32)/denom - dg1(3) = (df1(1)*a13-df1(2)*a23+df1(3)*a33)/denom - dg2(1) = (df2(1)*a11-df2(2)*a21+df2(3)*a31)/denom - dg2(2) = -(df2(1)*a12-df2(2)*a22+df2(3)*a32)/denom - dg2(3) = (df2(1)*a13-df2(2)*a23+df2(3)*a33)/denom - dg3(1) = (df3(1)*a11-df3(2)*a21+df3(3)*a31)/denom - dg3(2) = -(df3(1)*a12-df3(2)*a22+df3(3)*a32)/denom - dg3(3) = (df3(1)*a13-df3(2)*a23+df3(3)*a33)/denom -c - return - end -c -c -c - subroutine rkint4(sox,xgcn,bres) -c Runge-Kutta integrator -c - use const_and_precisions, only : wp_ - use gray_params, only : igrad - use beamdata, only : dst,nrayr,nrayth,ywrk,ypwrk,grad2,dgrad2v, - . gri,ggri - implicit none -c arguments - real(wp_), intent(in) :: sox,xgcn,bres -c parameter - integer, parameter :: ndim=6 -c local variables - integer :: ieq,iv,j,jv,k,kkk - real(wp_) :: h,hh,h6 - real(wp_), dimension(ndim) :: y,yy,fk1,fk2,fk3,fk4 -c common/external functions/variables - real(wp_) :: gr2 - real(wp_), dimension(3) :: dgr2,dgr - real(wp_), dimension(3,3) :: ddgr -c - common/gr/gr2 - common/dgr/dgr2,dgr,ddgr -c - h=dst - hh=h*0.5_wp_ - h6=h/6.0_wp_ -c - do j=1,nrayr - kkk=nrayth - if(j.eq.1) kkk=1 - do k=1,kkk - do ieq=1,ndim - y(ieq)=ywrk(ieq,j,k) - fk1(ieq)=ypwrk(ieq,j,k) - yy(ieq)=y(ieq)+fk1(ieq)*hh - end do - - if(igrad.eq.1) then - gr2=grad2(j,k) - do iv=1,3 - dgr2(iv)=dgrad2v(iv,j,k) - dgr(iv)=gri(iv,j,k) - do jv=1,3 - ddgr(iv,jv)=ggri(iv,jv,j,k) - end do - end do - end if - - call fwork(sox,xgcn,bres,yy,fk2) -c - do ieq=1,ndim - yy(ieq)=y(ieq)+fk2(ieq)*hh - end do - call fwork(sox,xgcn,bres,yy,fk3) -c - do ieq=1,ndim - yy(ieq)=y(ieq)+fk3(ieq)*h - end do - call fwork(sox,xgcn,bres,yy,fk4) -c - do ieq=1,ndim - ywrk(ieq,j,k)=y(ieq) - . +h6*(fk1(ieq)+2.0_wp_*fk2(ieq)+2.0_wp_*fk3(ieq)+fk4(ieq)) - end do - end do - end do -c - call updatepos -c - if(igrad.eq.1) call gradi -c - return - end -c -c -c - subroutine gwork(sox,xgcn,bres,j,k) - use const_and_precisions, only : wp_ - use gray_params, only : igrad - use beamdata, only : ywrk,ypwrk,grad2,dgrad2v,gri,ggri - implicit none -c local constants - integer, parameter :: ndim=6 -c arguments - real(wp_), intent(in) :: sox,xgcn,bres - integer :: j,k -c local variables - integer :: ieq,iv,jv - real(wp_), dimension(ndim) :: yy,yyp -c common/external functions/variables - real(wp_) :: gr2 - real(wp_), dimension(3) :: dgr2,dgr - real(wp_), dimension(3,3) :: ddgr -c - common/gr/gr2 - common/dgr/dgr2,dgr,ddgr -c -c begin --- update vector yy -c - do ieq=1,ndim - yy(ieq)=ywrk(ieq,j,k) - end do -c - if(igrad.eq.1) then - gr2=grad2(j,k) - do iv=1,3 - dgr2(iv)=dgrad2v(iv,j,k) - dgr(iv)=gri(iv,j,k) - do jv=1,3 - ddgr(iv,jv)=ggri(iv,jv,j,k) - end do - end do - end if -c - call fwork(sox,xgcn,bres,yy,yyp) -c - do ieq=1,ndim - ypwrk(ieq,j,k)=yyp(ieq) - end do -c -c end --- update vector yy -c - return - end -c -c -c - subroutine fwork(sox,xgcn,bres,y,dery) - use const_and_precisions, only : wp_ - use gray_params, only : idst,igrad - implicit none -c local constants - integer, parameter :: ndim=6 -c arguments - real(wp_), dimension(ndim) :: y,dery - real(wp_), intent(in) :: sox,xgcn,bres -c local variables - integer :: iv,jv - real(wp_) :: xx,yy,zz,yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl, - . dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel, - . derdom,dfdiadnpl,dfdiadxg,dfdiadyg - real(wp_), dimension(3) :: vgv,derdxv,danpldxv,derdnv,dbgr -c common/external functions/variables - integer :: ierr - real(wp_) :: gr2,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,anpl, - . anpr,xg,yg,vgm,derdnm,dersdst - real(wp_), dimension(3) :: dgr2,dgr,xv,anv,bv,derxg,deryg - real(wp_), dimension(3,3) :: ddgr,derbv -c - common/gr/gr2 - common/dgr/dgr2,dgr,ddgr - common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 - common/nplr/anpl,anpr - common/bb/bv - common/dbb/derbv - common/xgxg/xg - common/ygyg/yg - common/dxgyg/derxg,deryg - common/vgv/vgm,derdnm - common/dersdst/dersdst - common/ierr/ierr - common/anv/anv - common/xv/xv -c - xx=y(1) - yy=y(2) - zz=y(3) - - xv(1)=y(1) - xv(2)=y(2) - xv(3)=y(3) -c - anv(1) = y(4) - anv(2) = y(5) - anv(3) = y(6) -c -c rr=sqrt(xx**2+yy**2) -c phi=acos(xx/rr) -c if (yy.lt.0.0_wp_) then -c phi=-phi -c end if -c - call plas_deriv(xgcn,bres,xx,yy,zz) -c - an2=anv(1)*anv(1)+anv(2)*anv(2)+anv(3)*anv(3) - anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3) -c - if(abs(anpl).gt.0.99_wp_) then - if(abs(anpl).le.1.05_wp_) then - ierr=97 - else - ierr=98 - end if - end if -c - anpl2=anpl*anpl - dnl=1.0_wp_-anpl2 - anpr2=an2-anpl2 - anpr=0.0_wp_ - if(anpr2.gt.0.0_wp_) anpr=sqrt(anpr2) - yg2=yg**2 -c - an2s=1.0_wp_ - dan2sdxg=0.0_wp_ - dan2sdyg=0.0_wp_ - dan2sdnpl=0.0_wp_ - del=0.0_wp_ - fdia=0.0_wp_ - dfdiadnpl=0.0_wp_ - dfdiadxg=0.0_wp_ - dfdiadyg=0.0_wp_ -c - duh=1.0_wp_-xg-yg2 - if(xg.gt.0.0_wp_) then - del=sqrt(dnl*dnl+4.0_wp_*anpl*anpl*(1.0_wp_-xg)/yg2) - an2s=1.0_wp_-xg -xg*yg2*(1.0_wp_+anpl2+sox*del)/duh/2.0_wp_ -c - dan2sdxg=-yg2*(1.0_wp_-yg2)*(1.0_wp_+anpl2+sox*del)/ - . duh**2/2.0_wp_+sox*xg*anpl2/(del*duh)-1.0_wp_ - dan2sdyg=-xg*yg*(1.0_wp_-xg)*(1.0_wp_+anpl2+sox*del)/duh**2 - . +2.0_wp_*sox*xg*(1.0_wp_-xg)*anpl2/(yg*del*duh) - dan2sdnpl=-xg*yg2*anpl/duh - . -sox*xg*anpl*(2.0_wp_*(1.0_wp_-xg)-yg2*dnl)/(del*duh) -c - if(igrad.gt.0) then - ddelnpl2=2.0_wp_*(2.0_wp_*(1.0_wp_-xg) - . *(1.0_wp_+3.0_wp_*anpl2**2)-yg2*dnl**3)/yg2/del**3 - fdia=-xg*yg2*(1.0_wp_+sox*ddelnpl2/2.0_wp_)/duh - derdel=2.0_wp_*(1.0_wp_-xg)*anpl2*(1.0_wp_+3.0_wp_*anpl2**2) - . -dnl**2*(1.0_wp_+3.0_wp_*anpl2)*yg2 - derdel=4.0_wp_*derdel/(yg**5*del**5) - ddelnpl2y=2.0_wp_*(1.0_wp_-xg)*derdel - ddelnpl2x=yg*derdel - dfdiadnpl=24.0_wp_*sox*xg*(1.0_wp_-xg)*anpl*(1.0_wp_-anpl**4) - . /(yg2*del**5) - dfdiadxg=-yg2*(1.0_wp_-yg2)/duh**2-sox*yg2*((1.0_wp_-yg2) - . *ddelnpl2+xg*duh*ddelnpl2x)/(2.0_wp_*duh**2) - dfdiadyg=-2.0_wp_*yg*xg*(1.0_wp_-xg)/duh**2 - . -sox*xg*yg*(2.0_wp_*(1.0_wp_-xg)*ddelnpl2 - . +yg*duh*ddelnpl2y)/(2.0_wp_*duh**2) -c - end if - end if -c - bdotgr=0.0_wp_ - do iv=1,3 - bdotgr=bdotgr+bv(iv)*dgr(iv) - dbgr(iv)=0.0_wp_ - do jv=1,3 - dbgr(iv)=dbgr(iv)+dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv) - end do - end do -c - derdnm=0.0_wp_ -c - do iv=1,3 - danpldxv(iv)=anv(1)*derbv(1,iv)+anv(2)*derbv(2,iv) - . +anv(3)*derbv(3,iv) - derdxv(iv)=-(derxg(iv)*dan2sdxg - . +deryg(iv)*dan2sdyg+danpldxv(iv)*dan2sdnpl) - derdxv(iv)=derdxv(iv)-igrad*dgr2(iv) -c - derdxv(iv)=derdxv(iv)+fdia*bdotgr*dbgr(iv)+0.5_wp_*bdotgr**2* - . (derxg(iv)*dfdiadxg+deryg(iv)*dfdiadyg+danpldxv(iv)*dfdiadnpl) -c - derdnv(iv)=2.0_wp_*anv(iv)-bv(iv)*dan2sdnpl - . +0.5_wp_*bdotgr**2*bv(iv)*dfdiadnpl -c - derdnm=derdnm+derdnv(iv)**2 - end do -c - derdnm=sqrt(derdnm) -c - derdom=-2.0_wp_*an2+2.0_wp_*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl - . +2.0_wp_*igrad*gr2 - . -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0_wp_ - . +anpl*dfdiadnpl/2.0_wp_) -c - if (idst.eq.0) then -c integration variable: s - denom=-derdnm - else if (idst.eq.1) then -c integration variable: c*t - denom=derdom - else -c integration variable: Sr - denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3)) - end if -c -c coefficient for integration in s -c ds/dst, where st is the integration variable - dersdst=-derdnm/denom -c - dery(1) = -derdnv(1)/denom - dery(2) = -derdnv(2)/denom - dery(3) = -derdnv(3)/denom - dery(4) = derdxv(1)/denom - dery(5) = derdxv(2)/denom - dery(6) = derdxv(3)/denom -c -c vgv : ~ group velocity -c - vgm=0 - do iv=1,3 - vgv(iv)=-derdnv(iv)/derdom - vgm=vgm+vgv(iv)**2 - end do - vgm=sqrt(vgm) -c -c dd : dispersion relation (real part) -c ddi : dispersion relation (imaginary part) -c - dd=an2-an2s-igrad*(gr2-0.5_wp_*bdotgr**2*fdia) - ddi=derdnv(1)*dgr(1)+derdnv(2)*dgr(2)+derdnv(3)*dgr(3) -c - return - end -c -c -c - subroutine plas_deriv(xgcn,bres,xx,yy,zz) - use const_and_precisions, only : wp_,pi,ccj=>mu0inv - use gray_params, only : iequil - use equilibrium, only : equinum_fpol, equinum_psi, sgnbphi, equian - implicit none -c arguments - real(wp_) :: xgcn,bres,xx,yy,zz -c local variables - integer :: iv,jv - real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy, - . rr,rr2,rrm,snphi,zzm - real(wp_), dimension(3) :: dbtot,bvc - real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv -c common/external functions/variables - real(wp_) :: brr,bphi,bzz,ajphi,btot,xg,yg,dxgdpsi,dpsidr, - . dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,psinv - real(wp_), dimension(3) :: bv,derxg,deryg - real(wp_), dimension(3,3) :: derbv -c - common/parpl/brr,bphi,bzz,ajphi - common/btot/btot - common/bb/bv - common/dbb/derbv - common/xgxg/xg - common/dxgdps/dxgdpsi - common/ygyg/yg - common/dxgyg/derxg,deryg - common/derip1/dpsidr,dpsidz - common/derip2/ddpsidrr,ddpsidzz,ddpsidrz - common/fpol/fpolv,dfpv - common/psival/psinv - - xg=0.0_wp_ - yg=9.9e1_wp_ -c - do iv=1,3 - derxg(iv)=0.0_wp_ - deryg(iv)=0.0_wp_ - bv(iv)=0.0_wp_ - dbtot(iv)=0.0_wp_ - do jv=1,3 - dbv(iv,jv)=0.0_wp_ - derbv(iv,jv)=0.0_wp_ - dbvcdc(iv,jv)=0.0_wp_ - dbvcdc(iv,jv)=0.0_wp_ - dbvdc(iv,jv)=0.0_wp_ - end do - end do -c - if(iequil.eq.0) return -c -c cylindrical coordinates -c - rr2=xx**2+yy**2 - rr=sqrt(rr2) - csphi=xx/rr - snphi=yy/rr -c - bv(1)=-snphi*sgnbphi - bv(2)=csphi*sgnbphi -c -c convert from cm to meters -c - zzm=1.0e-2_wp_*zz - rrm=1.0e-2_wp_*rr -c - if(iequil.eq.1) then - call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, - . ddpsidrr,ddpsidzz,ddpsidrz) - end if -c - if(iequil.eq.2) then - call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz, - . ddpsidrz) - call equinum_fpol(psinv,fpolv,dfpv) - end if - - call sub_xg_derxg(xgcn,psinv) - yg=fpolv/(rrm*bres) - bphi=fpolv/rrm - btot=abs(bphi) - if(psinv.lt.0.0_wp_) return -c -c B = f(psi)/R e_phi+ grad psi x e_phi/R -c -c bvc(i) = B_i in cylindrical coordinates -c bv(i) = B_i in cartesian coordinates -c - bphi=fpolv/rrm - brr=-dpsidz/rrm - bzz= dpsidr/rrm -c - bvc(1)=brr - bvc(2)=bphi - bvc(3)=bzz -c - bv(1)=bvc(1)*csphi-bvc(2)*snphi - bv(2)=bvc(1)*snphi+bvc(2)*csphi - bv(3)=bvc(3) -c - b2tot=bv(1)**2+bv(2)**2+bv(3)**2 - btot=sqrt(b2tot) -c -c dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) -c - dbvcdc(1,1)=-ddpsidrz/rrm-brr/rrm - dbvcdc(1,3)=-ddpsidzz/rrm - dbvcdc(2,1)= dfpv*dpsidr/rrm-bphi/rrm - dbvcdc(2,3)= dfpv*dpsidz/rrm - dbvcdc(3,1)= ddpsidrr/rrm-bzz/rrm - dbvcdc(3,3)= ddpsidrz/rrm -c -c dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) -c - dbvdc(1,1) = dbvcdc(1,1)*csphi-dbvcdc(2,1)*snphi - dbvdc(1,2) = -bv(2) - dbvdc(1,3) = dbvcdc(1,3)*csphi-dbvcdc(2,3)*snphi - dbvdc(2,1) = dbvcdc(1,1)*snphi+dbvcdc(2,1)*csphi - dbvdc(2,2) = bv(1) - dbvdc(2,3) = dbvcdc(1,3)*snphi+dbvcdc(2,3)*csphi - dbvdc(3,1) = dbvcdc(3,1) - dbvdc(3,2) = dbvcdc(3,2) - dbvdc(3,3) = dbvcdc(3,3) -c - drrdx=csphi - drrdy=snphi - dphidx=-snphi/rrm - dphidy=csphi/rrm -c -c dbv(iv,jv) = d Bcart(iv) / dxvcart(jv) -c - dbv(1,1)=drrdx*dbvdc(1,1)+dphidx*dbvdc(1,2) - dbv(1,2)=drrdy*dbvdc(1,1)+dphidy*dbvdc(1,2) - dbv(1,3)=dbvdc(1,3) - dbv(2,1)=drrdx*dbvdc(2,1)+dphidx*dbvdc(2,2) - dbv(2,2)=drrdy*dbvdc(2,1)+dphidy*dbvdc(2,2) - dbv(2,3)=dbvdc(2,3) - dbv(3,1)=drrdx*dbvdc(3,1)+dphidx*dbvdc(3,2) - dbv(3,2)=drrdy*dbvdc(3,1)+dphidy*dbvdc(3,2) - dbv(3,3)=dbvdc(3,3) -c - dbtot(1)=(bv(1)*dbv(1,1)+bv(2)*dbv(2,1)+bv(3)*dbv(3,1))/btot - dbtot(2)=(bv(1)*dbv(1,2)+bv(2)*dbv(2,2)+bv(3)*dbv(3,2))/btot - dbtot(3)=(bv(1)*dbv(1,3)+bv(2)*dbv(2,3)+bv(3)*dbv(3,3))/btot -c - yg=btot/Bres -c -c convert spatial derivatives from dummy/m -> dummy/cm -c to be used in fwork -c -c bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z -c - do iv=1,3 - deryg(iv)=1.0e-2_wp_*dbtot(iv)/Bres - bv(iv)=bv(iv)/btot - do jv=1,3 - derbv(iv,jv)=1.0e-2_wp_*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot - end do - end do -c - derxg(1)=1.0e-2_wp_*drrdx*dpsidr*dxgdpsi - derxg(2)=1.0e-2_wp_*drrdy*dpsidr*dxgdpsi - derxg(3)=1.0e-2_wp_*dpsidz*dxgdpsi -c -c current density computation in Ampere/m^2 -c ccj==1/mu_0 -c - ajphi=ccj*(dbvcdc(1,3)-dbvcdc(3,1)) -c ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3)) -c ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2)) -c - return - end -c -c -c -C subroutine tor_curr_psi(h,ajphi) -C use const_and_precisions, only : wp_ -C use equilibrium, only : zmaxis -C implicit none -Cc arguments -C real(wp_) :: h,ajphi -Cc local variables -C real(wp_) :: r1,r2 -Cc -C call psi_raxis(h,r1,r2) -C call tor_curr(r2,zmaxis,ajphi) -Cc -C return -C end -Cc -Cc -Cc -C subroutine tor_curr(rpsim,zpsim,ajphi) -C use const_and_precisions, only : wp_,ccj=>mu0inv -C use equilibrium, only : equinum_psi -C implicit none -Cc arguments -C real(wp_) :: rpsim,zpsim,ajphi -Cc local variables -C real(wp_) :: bzz,dbvcdc13,dbvcdc31 -C real(wp_) :: dpsidr,ddpsidrr,ddpsidzz -Cc -C call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,ddpsidrr=ddpsidrr, -C . ddpsidzz=ddpsidzz) -C bzz= dpsidr/rpsim -C dbvcdc13=-ddpsidzz/rpsim -C dbvcdc31= ddpsidrr/rpsim-bzz/rpsim -C ajphi=ccj*(dbvcdc13-dbvcdc31) -Cc -C return -C end -Cc -Cc -Cc -C subroutine psi_raxis(h,r1,r2) -C use const_and_precisions, only : wp_ -C use equilibrium, only : psiant,psinop,zmaxis,nsr, -C . nsz,cc=>cceq,tr,tz,kspl -C use dierckx, only : profil,sproota -C implicit none -Cc local constants -C integer, parameter :: mest=4 -Cc arguments -C real(wp_) :: h,r1,r2 -Cc local variables -C integer :: iopt,ier,m -C real(wp_) :: zc,val -C real(wp_), dimension(mest) :: zeroc -C real(wp_), dimension(nsr) :: czc -Cc -C iopt=1 -C zc=zmaxis -C call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) -C if(ier.gt.0) print*,' profil =',ier -C val=h*psiant+psinop -C call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) -C r1=zeroc(1) -C r2=zeroc(2) -Cc -C return -C end -c -c -c - subroutine sub_xg_derxg(xgcn,psinv) - use const_and_precisions, only : wp_ - use equilibrium, only : psia - use coreprofiles, only : density - implicit none -c arguments - real(wp_) :: xgcn,psinv -c common/external functions/variables - real(wp_) :: xg,dxgdpsi - real(wp_) :: dens,ddenspsin -c - common/xgxg/xg - common/dxgdps/dxgdpsi - common/dens/dens,ddenspsin -c - xg=0.0_wp_ - dxgdpsi=0.0_wp_ -c if(psinv.le.psdbnd.and.psinv.ge.0) then - call density(psinv,dens,ddenspsin) - xg=xgcn*dens - dxgdpsi=xgcn*ddenspsin/psia -c end if -c - return - end -c -c -c - subroutine ic_gb(x00,y00,z00,anx0c,any0c,anz0c,ak0,xgcn,bres, - . wcsi,weta,rcicsi,rcieta,phiw,phir,sox,psipol0,chipol0) -c beam tracing initial conditions igrad=1 -c - use const_and_precisions, only : wp_,izero,zero,one,pi, - . cvdr=>degree,ui=>im - use math, only : catand - use gray_params, only : ipol - use beamdata, only : nrayr,nrayth,rwmax,ywrk0=>ywrk,ypwrk0=>ypwrk, - . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt - implicit none -c arguments - real(wp_), intent(in) :: x00,y00,z00,anx0c,any0c,anz0c - real(wp_), intent(in) :: ak0,xgcn,bres - real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir - real(wp_), intent(in) :: sox,psipol0,chipol0 -c local constants - integer, parameter :: ndim=6,ndimm=3 -c local variables - integer :: j,k,iproj,nfilp - real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw, - . wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy,rcixx, - . rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy,drcixx,drciyy, - . drcixy,dr,da,ddfu,u,alfak,dcsiw,detaw,dx0t,dy0t,x0t,y0t,z0t,dx0, - . dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2,gxxt,gyyt,gzzt,gxyt,gxzt,gyzt, - . dgr2xt,dgr2yt,dgr2zt,dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy, - . anzt,anxt,anyt,anx,any,anz,an20,an0,anx0,any0,anz0,du1tx,du1ty, - . du1tz,vgradi,r0,x0m,y0m,r0m,z0m,ddr110,deltapol,qq,uu,vv - real(wp_), dimension(ndim) :: ytmp,yptmp - complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1, - . dqi2,dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy -c common/external functions/variables - real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi, - . ddr11,psinv,dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi, - . psipol,chipol,psinv11 -c - common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 - common/nplr/anpl,anpr - common/psival/psinv - common/parpl/brr,bphi,bzz,ajphi - common/dens/dens,ddens - common/tete/tekev - common/polcof/psipol,chipol - common/psinv11/psinv11 -c - csth=anz0c - snth=sqrt(1.0_wp_-csth**2) - csps=1.0_wp_ - snps=0.0_wp_ - if(snth.gt.0.0_wp_) then - csps=any0c/snth - snps=anx0c/snth - end if -c - phiwrad=phiw*cvdr - phirrad=phir*cvdr - csphiw=cos(phiwrad) - snphiw=sin(phiwrad) -c csphir=cos(phirrad) -c snphir=sin(phirrad) -c - wwcsi=2.0_wp_/(ak0*wcsi**2) - wweta=2.0_wp_/(ak0*weta**2) -c - if(phir.ne.phiw) then - sk=(rcicsi+rcieta) - sw=(wwcsi+wweta) - dk=(rcicsi-rcieta) - dw=(wwcsi-wweta) - ts=-(dk*sin(2.0_wp_*phirrad)-ui*dw*sin(2.0_wp_*phiwrad)) - tc=(dk*cos(2.0_wp_*phirrad)-ui*dw*cos(2.0_wp_*phiwrad)) - phic=0.5_wp_*catand(ts/tc) - ddd=dk*cos(2*(phirrad+phic))-ui*dw*cos(2*(phiwrad+phic)) - sss=sk-ui*sw - qi1=0.5_wp_*(sss+ddd) - qi2=0.5_wp_*(sss-ddd) - rci1=dble(qi1) - ww1=-dimag(qi1) - rci2=dble(qi2) - ww2=-dimag(qi2) - else - rci1=rcicsi - rci2=rcieta - ww1=wwcsi - ww2=wweta - phic=-phiwrad - qi1=rci1-ui*ww1 - qi2=rci2-ui*ww2 - end if - -c w01=sqrt(2.0_wp_/(ak0*ww1)) -c z01=-rci1/(rci1**2+ww1**2) -c w02=sqrt(2.0_wp_/(ak0*ww2)) -c z02=-rci2/(rci2**2+ww2**2) - - qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2 - qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2 - qqxy=-(qi1-qi2)*sin(2.0_wp_*phic) - wwxx=-dimag(qqxx) - wwyy=-dimag(qqyy) - wwxy=-dimag(qqxy)/2.0_wp_ - rcixx=dble(qqxx) - rciyy=dble(qqyy) - rcixy=dble(qqxy)/2.0_wp_ - - dqi1=-qi1**2 - dqi2=-qi2**2 - d2qi1=2*qi1**3 - d2qi2=2*qi2**3 - dqqxx=dqi1*cos(phic)**2+dqi2*sin(phic)**2 - dqqyy=dqi1*sin(phic)**2+dqi2*cos(phic)**2 - dqqxy=-(dqi1-dqi2)*sin(2.0_wp_*phic) - d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2 - d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2 - d2qqxy=-(d2qi1-d2qi2)*sin(2.0_wp_*phic) - dwwxx=-dimag(dqqxx) - dwwyy=-dimag(dqqyy) - dwwxy=-dimag(dqqxy)/2.0_wp_ - d2wwxx=-dimag(d2qqxx) - d2wwyy=-dimag(d2qqyy) - d2wwxy=-dimag(d2qqxy)/2.0_wp_ - drcixx=dble(dqqxx) - drciyy=dble(dqqyy) - drcixy=dble(dqqxy)/2.0_wp_ - - dr=1.0_wp_ - if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) - da=2.0_wp_*pi/dble(nrayth) -c - ddfu=2.0_wp_*dr**2/ak0 - do j=1,nrayr - u=dble(j-1) -c ffi=u**2*ddfu/2.0_wp_ - dffiu(j)=u*ddfu - ddffiu(j)=ddfu - do k=1,nrayth - alfak=(k-1)*da - dcsiw=dr*cos(alfak)*wcsi - detaw=dr*sin(alfak)*weta - dx0t=dcsiw*csphiw-detaw*snphiw - dy0t=dcsiw*snphiw+detaw*csphiw - x0t=u*dx0t - y0t=u*dy0t - z0t=-0.5_wp_*(rcixx*x0t**2+rciyy*y0t**2+2*rcixy*x0t*y0t) -c -c csiw=u*dcsiw -c etaw=u*detaw -c csir=x0t*csphir+y0t*snphir -c etar=-x0t*snphir+y0t*csphir -c - dx0= x0t*csps+snps*(y0t*csth+z0t*snth) - dy0=-x0t*snps+csps*(y0t*csth+z0t*snth) - dz0= z0t*csth-y0t*snth - x0=x00+dx0 - y0=y00+dy0 - z0=z00+dz0 - - gxt=x0t*wwxx+y0t*wwxy - gyt=x0t*wwxy+y0t*wwyy - gzt=0.5_wp_*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy - gr2=gxt*gxt+gyt*gyt+gzt*gzt - gxxt=wwxx - gyyt=wwyy - gzzt=0.5_wp_*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy - gxyt=wwxy - gxzt=x0t*dwwxx+y0t*dwwxy - gyzt=x0t*dwwxy+y0t*dwwyy - dgr2xt=2.0_wp_*(gxt*gxxt+gyt*gxyt+gzt*gxzt) - dgr2yt=2.0_wp_*(gxt*gxyt+gyt*gyyt+gzt*gyzt) - dgr2zt=2.0_wp_*(gxt*gxzt+gyt*gyzt+gzt*gzzt) - dgr2x= dgr2xt*csps+snps*(dgr2yt*csth+dgr2zt*snth) - dgr2y=-dgr2xt*snps+csps*(dgr2yt*csth+dgr2zt*snth) - dgr2z= dgr2zt*csth-dgr2yt*snth - gri(1,j,k)=gxt*csps+snps*(gyt*csth+gzt*snth) - gri(2,j,k)=-gxt*snps+csps*(gyt*csth+gzt*snth) - gri(3,j,k)=gzt*csth-gyt*snth - ggri(1,1,j,k)=gxxt*csps**2+snps**2 - . *(gyyt*csth**2+gzzt*snth**2 - . +2.0_wp_*snth*csth*gyzt) - . +2.0_wp_*snps*csps*(gxyt*csth+gxzt*snth) - ggri(2,2,j,k)=gxxt*snps**2+csps**2 - . *(gyyt*csth**2+gzzt*snth**2 - . +2.0_wp_*snth*csth*gyzt) - . -2.0_wp_*snps*csps*(gxyt*csth+gxzt*snth) - ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2 - . -2.0_wp_*csth*snth*gyzt - ggri(1,2,j,k)=csps*snps*(-gxxt+csth**2*gyyt+snth**2*gzzt - . +2.0_wp_*csth*snth*gyzt) - . +(csps**2-snps**2)*(snth*gxzt+csth*gxyt) - ggri(1,3,j,k)=csth*snth*snps*(gzzt-gyyt)+(csth**2-snth**2) - . *snps*gyzt+csps*(csth*gxzt-snth*gxyt) - ggri(2,3,j,k)=csth*snth*csps*(gzzt-gyyt)+(csth**2-snth**2) - . *csps*gyzt+snps*(snth*gxyt-csth*gxzt) - ggri(2,1,j,k)=ggri(1,2,j,k) - ggri(3,1,j,k)=ggri(1,3,j,k) - ggri(3,2,j,k)=ggri(2,3,j,k) -c - du1tx=(dx0t*wwxx+dy0t*wwxy)/ddfu - du1ty=(dx0t*wwxy+dy0t*wwyy)/ddfu - du1tz=0.5_wp_*u*(dx0t**2*dwwxx+dy0t**2*dwwyy - . +2.0_wp_*dx0t*dy0t*dwwxy)/ddfu -c - pppx=x0t*rcixx+y0t*rcixy - pppy=x0t*rcixy+y0t*rciyy - denpp=pppx*gxt+pppy*gyt - if (denpp.ne.0.0_wp_) then - ppx=-pppx*gzt/denpp - ppy=-pppy*gzt/denpp - else - ppx=0.0_wp_ - ppy=0.0_wp_ - end if -c - anzt=sqrt((1.0_wp_+gr2)/(1.0_wp_+ppx**2+ppy**2)) - anxt=ppx*anzt - anyt=ppy*anzt -c - anx= anxt*csps+snps*(anyt*csth+anzt*snth) - any=-anxt*snps+csps*(anyt*csth+anzt*snth) - anz= anzt*csth-anyt*snth -c - an20=1.0_wp_+gr2 - an0=sqrt(an20) - anx0=anx - any0=any - anz0=anz -c - xc0(1,j,k)=x0 - xc0(2,j,k)=y0 - xc0(3,j,k)=z0 -c - ywrk0(1,j,k)=x0 - ywrk0(2,j,k)=y0 - ywrk0(3,j,k)=z0 - ywrk0(4,j,k)=anx0 - ywrk0(5,j,k)=any0 - ywrk0(6,j,k)=anz0 -c - ypwrk0(1,j,k) = anx0/an0 - ypwrk0(2,j,k) = any0/an0 - ypwrk0(3,j,k) = anz0/an0 - ypwrk0(4,j,k) = dgr2x/an0/2.0_wp_ - ypwrk0(5,j,k) = dgr2y/an0/2.0_wp_ - ypwrk0(6,j,k) = dgr2z/an0/2.0_wp_ -c - ytmp=ywrk0(:,j,k) - yptmp=ypwrk0(:,j,k) - call fwork(sox,xgcn,bres,ytmp,yptmp) - - if(ipol.eq.0) then - call pol_limit(sox,ext(j,k,0),eyt(j,k,0)) - qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 - uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) - vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) - call polellipse(qq,uu,vv,psipol0,chipol0) - else - qq=cos(2.0_wp_*chipol0*cvdr)*cos(2.0_wp_*psipol0*cvdr) - uu=cos(2.0_wp_*chipol0*cvdr)*sin(2.0_wp_*psipol0*cvdr) - vv=sin(2.0_wp_*chipol0*cvdr) - if(qq**2.lt.1) then - deltapol=asin(vv/sqrt(1.0_wp_-qq**2)) - ext(j,k,0)= sqrt((1.0_wp_+qq)/2) - eyt(j,k,0)= sqrt((1.0_wp_-qq)/2)*exp(-ui*deltapol) - else - ext(j,k,0)= 1.0_wp_ - eyt(j,k,0)= 0.0_wp_ - end if - endif - psipol=psipol0 - chipol=chipol0 -c - grad2(j,k)=gr2 - dgrad2v(1,j,k)=dgr2x - dgrad2v(2,j,k)=dgr2y - dgrad2v(3,j,k)=dgr2z -c - du10(1,j,k)= du1tx*csps+snps*(du1ty*csth+du1tz*snth) - du10(2,j,k)=-du1tx*snps+csps*(du1ty*csth+du1tz*snth) - du10(3,j,k)= du1tz*csth-du1ty*snth -c - dd=anx0**2+any0**2+anz0**2-an20 - vgradi=anxt*gxt+anyt*gyt+anzt*gzt - ddi=2.0_wp_*vgradi -c - r0=sqrt(x0**2+y0**2) - x0m=x0/1.0e2_wp_ - y0m=y0/1.0e2_wp_ - r0m=r0/1.0e2_wp_ - z0m=z0/1.0e2_wp_ - if(j.eq.nrayr) then - write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . psinv,zero,anpl,zero,one - end if - if(j.eq.1.and.k.eq.1) then - write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, - . psinv,one,dens,tekev,brr,bphi,bzz, - . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, - . zero,zero,zero,zero,zero,zero,zero,one - - psinv11=psinv - ddr110=dd - end if - if(j.eq.nrayr.and.k.eq.1) then - write(17,99) zero,ddr110,dd,ddi - end if - end do - end do - - call pweigth -c - if(nrayr.gt.1) then - iproj=0 - nfilp=8 - call projxyzt(iproj,nfilp) - end if -c - return - 99 format(24(1x,e16.8e3)) -111 format(3i5,20(1x,e16.8e3)) - end -c -c -c - subroutine ic_rt(x00,y00,z00,anx0c,any0c,anz0c,ak0,xgcn,bres, - . wcsi,weta,rcicsi,rcieta,phiw,phir,sox,psipol0,chipol0) -c ray tracing initial conditions igrad=0 -c - use const_and_precisions, only : wp_,izero,zero,one,pi, - . cvdr=>degree,ui=>im - use gray_params, only : ipol - use beamdata, only : nrayr,nrayth,rwmax,ywrk0=>ywrk,ypwrk0=>ypwrk, - . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt - implicit none -c arguments - real(wp_), intent(in) :: x00,y00,z00,anx0c,any0c,anz0c - real(wp_), intent(in) :: ak0,xgcn,bres - real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir - real(wp_), intent(in) :: sox,psipol0,chipol0 -c local constants - integer, parameter :: ndim=6,ndimm=3 -c local variables - integer :: j,k,iv,jv,iproj,nfilp - real(wp_) :: csth,snth,csps,snps,phiwrad,csphiw,snphiw,dr,da,u, - . alfak,dcsiw,detaw,dx0t,dy0t,x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0, - . anzt,anxt,anyt,anx,any,anz,an20,an0,anx0,any0,anz0,vgradi,r0, - . x0m,y0m,r0m,z0m,ancsi,aneta,ppcsi,ppeta,deltapol,qq,uu,vv - real(wp_), dimension(ndim) :: ytmp,yptmp -c common/external functions/variables - real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv,dens,ddens, - . tekev,anpl,anpr,brr,bphi,bzz,ajphi,psipol,chipol,psinv11 - -c - common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 - common/nplr/anpl,anpr - common/psival/psinv - common/parpl/brr,bphi,bzz,ajphi - common/dens/dens,ddens - common/tete/tekev - common/polcof/psipol,chipol - common/psinv11/psinv11 -c - csth=anz0c - snth=sqrt(1.0_wp_-csth**2) - csps=1.0_wp_ - snps=0.0_wp_ - if(snth.gt.0.0_wp_) then - csps=any0c/snth - snps=anx0c/snth - end if -c - phiwrad=phiw*cvdr - csphiw=cos(phiwrad) - snphiw=sin(phiwrad) -c - dr=1.0_wp_ - if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) - da=2.0_wp_*pi/dble(nrayth) - z0t=0.0_wp_ -c - do j=1,nrayr - u=dble(j-1) - dffiu(j)=0.0_wp_ - ddffiu(j)=0.0_wp_ - do k=1,nrayth - alfak=(k-1)*da - dcsiw=dr*cos(alfak)*wcsi - detaw=dr*sin(alfak)*weta - dx0t=dcsiw*csphiw-detaw*snphiw - dy0t=dcsiw*snphiw+detaw*csphiw - x0t=u*dx0t - y0t=u*dy0t -c -c csiw=u*dcsiw -c etaw=u*detaw -c csir=csiw -c etar=etaw -c - dx0= x0t*csps+snps*(y0t*csth+z0t*snth) - dy0=-x0t*snps+csps*(y0t*csth+z0t*snth) - dz0= z0t*csth-y0t*snth -c - x0=x00+dx0 - y0=y00+dy0 - z0=z00+dz0 -c - ppcsi=u*dr*cos(alfak)*rcicsi - ppeta=u*dr*sin(alfak)*rcieta -c - anzt=1.0_wp_/sqrt(1.0_wp_+ppcsi**2+ppeta**2) - ancsi=ppcsi*anzt - aneta=ppeta*anzt -c - anxt=ancsi*csphiw-aneta*snphiw - anyt=ancsi*snphiw+aneta*csphiw -c - anx= anxt*csps+snps*(anyt*csth+anzt*snth) - any=-anxt*snps+csps*(anyt*csth+anzt*snth) - anz= anzt*csth-anyt*snth -c - an20=1.0_wp_ - an0=sqrt(an20) - anx0=anx - any0=any - anz0=anz -c - xc0(1,j,k)=x0 - xc0(2,j,k)=y0 - xc0(3,j,k)=z0 -c - ywrk0(1,j,k)=x0 - ywrk0(2,j,k)=y0 - ywrk0(3,j,k)=z0 - ywrk0(4,j,k)=anx0 - ywrk0(5,j,k)=any0 - ywrk0(6,j,k)=anz0 -c - ypwrk0(1,j,k) = anx0/an0 - ypwrk0(2,j,k) = any0/an0 - ypwrk0(3,j,k) = anz0/an0 - ypwrk0(4,j,k) = 0.0_wp_ - ypwrk0(5,j,k) = 0.0_wp_ - ypwrk0(6,j,k) = 0.0_wp_ -c - ytmp=ywrk0(:,j,k) - yptmp=ypwrk0(:,j,k) - call fwork(sox,xgcn,bres,ytmp,yptmp) - - if(ipol.eq.0) then - call pol_limit(sox,ext(j,k,0),eyt(j,k,0)) - qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 - uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) - vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) - call polellipse(qq,uu,vv,psipol0,chipol0) - else - qq=cos(2.0_wp_*chipol0*cvdr)*cos(2.0_wp_*psipol0*cvdr) - uu=cos(2.0_wp_*chipol0*cvdr)*sin(2.0_wp_*psipol0*cvdr) - vv=sin(2.0_wp_*chipol0*cvdr) - if(qq**2.lt.1.0_wp_) then -c deltapol=phix-phiy, phix =0 - deltapol=atan2(vv,uu) - ext(j,k,0)= sqrt((1.0_wp_+qq)/2) - eyt(j,k,0)= sqrt((1.0_wp_-qq)/2)*exp(-ui*deltapol) - else - if(qq.gt.0.0_wp_) then - ext(j,k,0)= 1.0_wp_ - eyt(j,k,0)= 0.0_wp_ - else - eyt(j,k,0)= 1.0_wp_ - ext(j,k,0)= 0.0_wp_ - end if - end if - endif - psipol=psipol0 - chipol=chipol0 -c - do iv=1,3 - gri(iv,j,k)=0.0_wp_ - dgrad2v(iv,j,k)=0.0_wp_ - du10(iv,j,k)=0.0_wp_ - do jv=1,3 - ggri(iv,jv,j,k)=0.0_wp_ - end do - end do - grad2(j,k)=0.0_wp_ -c - dd=anx0**2+any0**2+anz0**2-an20 - vgradi=0.0_wp_ - ddi=2.0_wp_*vgradi -c - r0=sqrt(x0**2+y0**2) - x0m=x0/1.0e2_wp_ - y0m=y0/1.0e2_wp_ - r0m=r0/1.0e2_wp_ - z0m=z0/1.0e2_wp_ - if(j.eq.nrayr) then - write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . psinv,zero,anpl,zero,one - end if - if(j.eq.1.and.k.eq.1) then - psinv11=psinv - write(17,99) zero,zero,zero,zero - write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, - . psinv,one,dens,tekev,brr,bphi,bzz, - . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, - . zero,zero,zero,zero,zero,zero,zero,one - end if - end do - end do - - call pweigth -c - if(nrayr.gt.1) then - iproj=0 - nfilp=8 - call projxyzt(iproj,nfilp) - end if -c - return - 99 format(24(1x,e16.8e3)) -111 format(3i5,20(1x,e16.8e3)) - end -c -c -c - subroutine ic_rt2(sox,xgcn,bres) - use const_and_precisions, only : wp_,izero,zero,one,pi, - . cvdr=>degree - use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk, - . xc0=>xc,du10=>du1,grad2,dgrad2v,gri,ggri,yyrfl,ext,eyt - implicit none -c arguments - real(wp_), intent(in) :: sox,xgcn,bres -c local constants - integer, parameter :: ndim=6,ndimm=3 -c local variables - integer :: j,k,iv,jv,iproj,nfilp - real(wp_) :: x0,y0,z0,an20,an0,anx0,any0,anz0,vgradi, - . r0,x0m,y0m,r0m,z0m,strfl11,qq,uu,vv - real(wp_), dimension(ndim) :: ytmp,yptmp -c common/external functions/variables - real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv, - . dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi,psipol,chipol, - . psinv11 -c - common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 - common/nplr/anpl,anpr - common/psival/psinv - common/parpl/brr,bphi,bzz,ajphi - common/dens/dens,ddens - common/tete/tekev - common/polcof/psipol,chipol - common/psinv11/psinv11 -c - do j=1,nrayr - do k=1,nrayth - x0=yyrfl(j,k,1) - y0=yyrfl(j,k,2) - z0=yyrfl(j,k,3) - anx0=yyrfl(j,k,4) - any0=yyrfl(j,k,5) - anz0=yyrfl(j,k,6) - an20=anx0*anx0+any0*any0+anz0*anz0 - an0=sqrt(an20) -c - xc0(1,j,k)=x0 - xc0(2,j,k)=y0 - xc0(3,j,k)=z0 -c - ywrk0(1,j,k)=x0 - ywrk0(2,j,k)=y0 - ywrk0(3,j,k)=z0 - ywrk0(4,j,k)=anx0/an0 - ywrk0(5,j,k)=any0/an0 - ywrk0(6,j,k)=anz0/an0 -c - ypwrk0(1,j,k) = anx0/an0 - ypwrk0(2,j,k) = any0/an0 - ypwrk0(3,j,k) = anz0/an0 - ypwrk0(4,j,k) = 0.0_wp_ - ypwrk0(5,j,k) = 0.0_wp_ - ypwrk0(6,j,k) = 0.0_wp_ -c - ytmp=ywrk0(:,j,k) - yptmp=ypwrk0(:,j,k) - call fwork(sox,xgcn,bres,ytmp,yptmp) - - call pol_limit(sox,ext(j,k,0),eyt(j,k,0)) - qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 - uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) - vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) - call polellipse(qq,uu,vv,psipol,chipol) -c - do iv=1,3 - gri(iv,j,k)=0.0_wp_ - dgrad2v(iv,j,k)=0.0_wp_ - du10(iv,j,k)=0.0_wp_ - do jv=1,3 - ggri(iv,jv,j,k)=0.0_wp_ - end do - end do - grad2(j,k)=0.0_wp_ -c - dd=anx0**2+any0**2+anz0**2-an20 - vgradi=0.0_wp_ - ddi=2.0_wp_*vgradi -c - r0=sqrt(x0**2+y0**2) - x0m=x0/1.0e2_wp_ - y0m=y0/1.0e2_wp_ - r0m=r0/1.0e2_wp_ - z0m=z0/1.0e2_wp_ - if(j.eq.nrayr) then - write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . psinv,zero,anpl,zero,one - end if - if(j.eq.1.and.k.eq.1) then - psinv11=psinv - write(17,99) zero,zero,zero,zero - write(4,99) strfl11,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, - . psinv,one,dens,tekev,brr,bphi,bzz, - . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, - . zero,zero,zero,zero,zero,zero,zero,one - end if - end do - end do -c - call pweigth -c - if(nrayr.gt.1) then - iproj=0 - nfilp=8 - call projxyzt(iproj,nfilp) - end if -c - return - 99 format(24(1x,e16.8e3)) -111 format(3i5,20(1x,e16.8e3)) - end -c -c -c - subroutine pweigth -c ray power weigth coefficient q(j) -c - use const_and_precisions, only : wp_ - use beamdata, only : nrayr,nrayth,rwmax,q - implicit none -c local variables - integer :: j - real(wp_) :: dr,r0,r,w,w0,summ -c - q(1)=1.0_wp_ - 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 -Cc -C subroutine valpsisplpec(rhop,voli,areai) -C use const_and_precisions, only : wp_ -C use utils, only : locate -C use simplespline, only :spli -C use magsurf_data, only : rpstab,cvol,carea,npsi -C implicit none -Cc arguments -C real(wp_), intent(in) :: rhop -C real(wp_), intent(out) :: voli,areai -Cc local variables -C integer :: ip -C real(wp_) :: drh -Cc -C call locate(rpstab,npsi,rhop,ip) -C ip=min(max(1,ip),npsi-1) -C drh=rhop-rpstab(ip) -Cc -C voli=spli(cvol,npsi,ip,drh) -C areai=spli(carea,npsi,ip,drh) -Cc -C return -C end -c -C subroutine valpsisplpec(rhop,voli,areai) -C use const_and_precisions, only : wp_ -C use utils, only : locate -C use simplespline, only :spli -C use magsurf_data, only : rpstab,cvol,carea,npsi -C implicit none -Cc arguments -C real(wp_), intent(in) :: rhop -C real(wp_), intent(out) :: voli,areai -Cc local variables -C integer :: ip -C real(wp_) :: drh -Cc -C call locate(rpstab,npsi,rhop,ip) -C ip=min(max(1,ip),npsi-1) -C drh=rhop-rpstab(ip) -Cc -C voli=spli(cvol,npsi,ip,drh) -C areai=spli(carea,npsi,ip,drh) -Cc -C return -C end -Cc -Cc -Cc -C subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) -C use const_and_precisions, only : wp_ -C use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi -C use utils, only : locate -C use simplespline, only :spli -C implicit none -Cc arguments -C real(wp_) :: rpsi,ratjai,ratjbi,ratjpli -Cc local variables -C integer :: ip -C real(wp_) :: dps -Cc -C call locate(rpstab,npsi,rpsi,ip) -C ip=min(max(1,ip),npsi-1) -C dps=rpsi-rpstab(ip) -C ratjai=spli(cratja,npsi,ip,dps) -C ratjbi=spli(cratjb,npsi,ip,dps) -C ratjpli=spli(cratjpl,npsi,ip,dps) -Cc -C return -C end -c -c -c - subroutine pabs_curr(p0mw,sox,ak0,i,j,k) - use const_and_precisions, only : wp_,pi,mc2=>mc2_ - use gray_params, only : iwarm,ilarm,ieccd,imx - use beamdata, only : dst,psjki,tauv,alphav,pdjki,ppabs,currj, - . didst,ccci,q,tau1v - use coreprofiles, only : temp, fzeff - use dispersion, only : harmnumber, warmdisp - use equilibrium, only : rmaxis,sgnbphi - use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl - use magsurf_data, only : fluxval - - implicit none -c arguments - integer, intent(in) :: i,j,k - real(wp_), intent(in) :: p0mw, sox, ak0 -c local constants - real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ -c local variables - real(wp_) :: dvvl,rbavi,rrii,rhop0,rhop,dervoli, - . tau,effjcdav,dpdst,p0,pow,alpha0,didst0,ccci0 - integer :: lrm,ithn - real(wp_) :: amu,ratiovgr,rbn,rbx,zeff - real(wp_) :: cst2,bmxi,bmni,fci - real(wp_), dimension(:), allocatable :: eccdpar -c common/external functions/variables - integer :: nhmin,nhmax,iokhawa,ierr - real(wp_) :: dersdst - real(wp_) :: psinv,xg,yg,tekev,dens,ddens,btot - real(wp_) :: anpl,anpr,vgm,derdnm,anprre,anprim,alpha,effjcd, - . akim,tau0 - complex(wp_) :: ex,ey,ez -c - 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 - common/nplr/anpl,anpr - common/vgv/vgm,derdnm - common/nprw/anprre,anprim - common/epolar/ex,ey,ez - common/iokh/iokhawa - common/ierr/ierr -c - dvvl=1.0_wp_ - rbavi=1.0_wp_ - rrii=rmaxis - 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 -c -c absorption computation -c - alpha=0.0_wp_ - akim=0.0_wp_ - effjcd=0.0_wp_ -c - tekev=temp(psinv) - call fluxval(rhop,dervol=dervoli,rri=rrii,rbav=rbavi, - . bmn=bmni,bmx=bmxi,fc=fci) -c - 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.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) -c - 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 -c - if(ieccd.gt.0) then - ithn=1 - if(lrm.gt.nhmin) ithn=2 - zeff=fzeff(psinv) - rbn=btot/bmni - rbx=btot/bmxi - - 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 -c - alphav(j,k,i)=alpha - - tau=tau0+0.5_wp_*(alpha+alpha0)*dersdst*dst - tauv(j,k,i)=tau - - p0=p0mw*q(j)*exp(-tau1v(j,k)) - pow=p0*exp(-tau) - ppabs(j,k,i)=p0-pow - - 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 - - return - end -c -c -c -C subroutine valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) -C use const_and_precisions, only : wp_ -C use utils, only : locate -C use simplespline, only :spli,splid -C use magsurf_data, only : rpstab,crri,crbav,cvol,cbmn,cbmx,cfc,npsi -C implicit none -Cc arguments -C real(wp_), intent(in) :: rhop -C real(wp_), intent(out) :: rrii,rbavi,dervoli,bmni,bmxi,fci -Cc local variables -C integer :: ip -C real(wp_) :: drh -Cc -C call locate(rpstab,npsi,rhop,ip) -C ip=min(max(1,ip),npsi-1) -C drh=rhop-rpstab(ip) -C rrii=spli(crri,npsi,ip,drh) -C rbavi=spli(crbav,npsi,ip,drh) -C dervoli=splid(cvol,npsi,ip,drh) -C bmni=spli(cbmn,npsi,ip,drh) -C bmxi=spli(cbmx,npsi,ip,drh) -C fci=spli(cfc,npsi,ip,drh) -Cc -C return -C end -c -c -c - subroutine projxyzt(iproj,nfile) - use const_and_precisions, only : wp_ - use beamdata, only : ywrk,nrayr,nrayth - implicit none -c arguments - integer :: iproj,nfile -c local variables - integer :: jd,j,kkk,k - real(wp_) :: dir,rtimn,rtimx,dx,dy,dz,dirx,diry,dirz, - . csth1,snth1,csps1,snps1,xti,yti,zti,rti -c common/external functions/variables - integer :: istep - real(wp_) :: psinv11,st -c - common/psinv11/psinv11 - common/istep/istep - common/ss/st -c - rtimn=1.0e+30_wp_ - rtimx=-1.0e-30_wp_ -c - jd=1 - if(iproj.eq.0) jd=nrayr-1 - do j=1,nrayr,jd - kkk=nrayth - if(j.eq.1) kkk=1 - do k=1,kkk -c - dx=ywrk(1,j,k)-ywrk(1,1,1) - dy=ywrk(2,j,k)-ywrk(2,1,1) - dz=ywrk(3,j,k)-ywrk(3,1,1) -c - dirx=ywrk(4,j,k) - diry=ywrk(5,j,k) - dirz=ywrk(6,j,k) - dir=sqrt(dirx*dirx+diry*diry+dirz*dirz) -c - if(j.eq.1.and.k.eq.1) then - csth1=dirz/dir - snth1=sqrt(1.0_wp_-csth1**2) - csps1=1.0_wp_ - snps1=0.0_wp_ - if(snth1.gt.0.0_wp_) then - csps1=diry/(dir*snth1) - snps1=dirx/(dir*snth1) - end if - end if - xti= dx*csps1-dy*snps1 - yti=(dx*snps1+dy*csps1)*csth1-dz*snth1 - zti=(dx*snps1+dy*csps1)*snth1+dz*csth1 - rti=sqrt(xti**2+yti**2) -Cc -c dirxt= (dirx*csps1-diry*snps1)/dir -c diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir -c dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir -Cc -c if(k.eq.1) then -c xti1=xti -c yti1=yti -c zti1=zti -c rti1=rti -c end if - - if(.not.(iproj.eq.0.and.j.eq.1)) - . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 -c - if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti - if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti -c - end do -c -c if(.not.(iproj.eq.0.and.j.eq.1)) -c . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 - if(iproj.eq.1) write(nfile,*) ' ' - end do -c - write(nfile,*) ' ' -c - write(12,99) istep,st,psinv11,rtimn,rtimx - return - 99 format(i5,12(1x,e16.8e3)) -111 format(3i5,12(1x,e16.8e3)) - end - - - - subroutine pec_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1, - . dvol,darea,ratjav,ratjbv,ratjplv) - use const_and_precisions, only : wp_,zero,one,izero - use gray_params, only : nnd - use equilibrium, only : frhotor,frhopol - use magsurf_data, only : fluxval - implicit none - integer :: it,ipec - real(wp_) :: drt,rt,rt1,rhop1 - real(wp_) :: ratjai,ratjbi,ratjpli - real(wp_) :: voli0,voli1,areai0,areai1 - real(wp_), dimension(nnd), intent(in) :: rt_in - real(wp_), dimension(nnd), intent(out) :: rhop_tab,rhot_tab - real(wp_), dimension(0:nnd), intent(out) :: rtabpsi1 - real(wp_), dimension(nnd), intent(out) :: dvol,darea - real(wp_), dimension(nnd), intent(out) :: ratjav,ratjbv,ratjplv - -c ipec positive build equidistant grid dimension nnd -c ipec negative read input grid -c -c ipec=+/-1 rho_pol grid -c ipec=+/-2 rho_tor grid - - voli0=zero - areai0=zero - rtabpsi1(0)=zero - - do it=1,nnd - if(ipec > 0) then -c build equidistant radial grid - drt=one/dble(nnd-1) - rt=dble(it-1)*drt - else -c read radial grid from input - rt=rt_in(it) - drt=(rt_in(it+1)-rt)/2.0_wp_ - end if -c radial coordinate of i-(i+1) interval mid point - if(it < nnd) then - rt1=rt+drt/2.0_wp_ - else - rt1=one - end if - if (abs(ipec) == 1) then - rhop_tab(it)=rt - rhot_tab(it)=frhotor(rt) - rhop1=rt1 - else - rhot_tab(it)=rt - rhop_tab(it)=frhopol(rt) - rhop1=frhopol(rt1) - end if -c psi grid at mid points, dimension nnd+1, for use in pec_tab - rtabpsi1(it)=rhop1**2 - - call fluxval(rhop1,area=areai1,vol=voli1) - dvol(it)=abs(voli1-voli0) - darea(it)=abs(areai1-areai0) - voli0=voli1 - areai0=areai1 - - call fluxval(rhop_tab(it),ratja=ratjai,ratjb=ratjbi, - . ratjpl=ratjpli) - ratjav(it)=ratjai - ratjbv(it)=ratjbi - ratjplv(it)=ratjpli - end do - end subroutine pec_init - - - - subroutine spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv, - . pins,currins) - use const_and_precisions, only : wp_,one,zero,izero - use gray_params, only : nnd - use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv - implicit none -c local constants - real(wp_), parameter :: rtbc=one -c arguments - real(wp_), intent(in) :: pabs,currt - real(wp_), dimension(nstep):: xxi,ypt,yamp - real(wp_), dimension(nnd), intent(out) :: dpdv,ajphiv - real(wp_), dimension(nnd), intent(out) :: pins,currins - real(wp_), dimension(nnd), intent(in) :: dvol,darea - real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1 -c local variables - integer :: i,ii,j,k,kkk - real(wp_) :: spds,sccs,facpds,facjs - real(wp_), dimension(nnd) :: wdpdv,wajphiv - -c calculation of dP and dI over radial grid - dpdv=zero - ajphiv=zero - kkk=1 - do j=1,nrayr - if(j > 1) kkk=nrayth - do k=1,kkk - ii=iiv(j,k) - if (ii < nstep ) then - if(psjki(j,k,ii+1) /= zero) ii=ii+1 - end if - xxi=zero - ypt=zero - yamp=zero - do i=1,ii - xxi(i)=abs(psjki(j,k,i)) - if(xxi(i) <= one) then - ypt(i)=ppabs(j,k,i) - yamp(i)=ccci(j,k,i) - end if - end do - call pec_tab(xxi,ypt,yamp,ii,rtabpsi1, - . wdpdv,wajphiv) - dpdv=dpdv+wdpdv - ajphiv=ajphiv+wajphiv - end do - end do - -c here dpdv is still dP (not divided yet by dV) -c here ajphiv is still dI (not divided yet by dA) - - spds=zero - sccs=zero - do i=1,nnd - spds=spds+dpdv(i) - sccs=sccs+ajphiv(i) - pins(i)=spds - currins(i)=sccs - dpdv(i)=dpdv(i)/dvol(i) - ajphiv(i)=ajphiv(i)/darea(i) - end do - - facpds=one - facjs=one - if(spds > zero) facpds=pabs/spds - if(sccs /= zero) facjs=currt/sccs - - do i=1,nnd - dpdv(i)=facpds*dpdv(i) - ajphiv(i)=facjs*ajphiv(i) - end do - -c now dpdv is dP/dV [MW/m^3] -c now ajphiv is J_phi=dI/dA [MA/m^2] - - end subroutine spec - - - - subroutine pec_tab(xxi,ypt,yamp,ii,rtabpsi1,wdpdv,wajphiv) -c Power and current projected on psi grid - mid points - use const_and_precisions, only : wp_,one,zero - use gray_params, only : nnd - use beamdata, only : nstep - use utils, only : locatex,intlin -c arguments - integer, intent(in) :: ii - integer, parameter :: llmx = 21 - integer, dimension(llmx) ::isev - real(wp_), dimension(nstep), intent(in) :: xxi,ypt,yamp - real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1 - real(wp_), dimension(nnd), intent(out) :: wdpdv,wajphiv -c local variables - real(wp_) :: ppa1,ppa2,cci1,cci2,dppa,didst,rt1 - integer :: i,is,ise0,idecr,iise0,iise,iis,iis1 - integer :: ind1,ind2,iind,ind,indi,itb1 - - isev=0 - ise0=0 - idecr=-1 - is=1 - wdpdv=zero - wajphiv=zero - do i=1,ii - if(ise0 == 0) then - if(xxi(i) < one) then - ise0=i - isev(is)=i-1 - is=is+1 - end if - else - if (idecr == -1) then - if(xxi(i) > xxi(i-1)) then - isev(is)=i-1 - is=is+1 - idecr=1 - end if - else - if(xxi(i) > one) exit - if(xxi(i) < xxi(i-1)) then - isev(is)=i-1 - is=is+1 - idecr=-1 - end if - end if - end if - end do - - isev(is)=i-1 - ppa1=zero - cci1=zero - - do iis=1,is-1 - iis1=iis+1 - iise0=isev(iis) - iise=isev(iis1) - if (mod(iis,2) /= 0) then - idecr=-1 - ind1=nnd - ind2=2 - iind=-1 - else - idecr=1 - ind1=1 - ind2=nnd - iind=1 - end if - do ind=ind1,ind2,iind - indi=ind - if (idecr == -1) indi=ind-1 - rt1=rtabpsi1(indi) - call locatex(xxi,iise,iise0,iise,rt1,itb1) - if(itb1 >= iise0 .and. itb1 < iise) then - call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1),ypt(itb1+1), - . rt1,ppa2) - call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1), - . rt1,cci2) - dppa=ppa2-ppa1 - wdpdv(ind)=wdpdv(ind)+dppa - didst=cci2-cci1 - wajphiv(ind)=wajphiv(ind)+didst - ppa1=ppa2 - cci1=cci2 - end if - end do - end do - - end subroutine pec_tab - - - - subroutine postproc_profiles(pabs,currt,rhot_tab,dvol,darea, - . dpdv,ajphiv,rhotpav,drhotpav,rhotjava,drhotjava) -c radial average values over power and current density profile - use const_and_precisions, only : wp_,one,zero,izero,pi - use gray_params, only : nnd,ieccd - use equilibrium, only : frhopol - use magsurf_data, only : fluxval - implicit none - real(wp_), intent(in) :: pabs,currt - real(wp_), intent(out) :: rhotpav,rhotjava - real(wp_), intent(out) :: drhotpav,drhotjava - real(wp_) :: rhopjava,rhoppav - real(wp_) :: dpdvp,dpdvmx,rhotp,drhotp - real(wp_) :: ajphip,ajmxfi,rhotjfi,drhotjfi - real(wp_) :: ratjamx,ratjbmx,ratjplmx - real(wp_), dimension(nnd), intent(in) :: rhot_tab,dvol,darea - real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv - - real(wp_) :: sccsa - real(wp_) :: rhotjav,rhot2pav,rhot2java,dvdrhotav,dadrhotava - - rhotpav=zero - rhot2pav=zero - rhotjav=zero - rhotjava=zero - rhot2java=zero - - if (pabs > zero) then - rhotpav=sum(rhot_tab(1:nnd)*dpdv(1:nnd)*dvol(1:nnd))/pabs - rhot2pav=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)* - . dpdv(1:nnd)*dvol(1:nnd))/pabs - end if - - if (ieccd /= 0) then - if (abs(currt) > zero) then - rhotjav=sum(rhot_tab(1:nnd)*ajphiv(1:nnd)*darea(1:nnd))/currt - end if - sccsa=sum(abs(ajphiv(1:nnd))*darea(1:nnd)) - if (sccsa > zero) then - rhotjava=sum(rhot_tab(1:nnd)*abs(ajphiv(1:nnd)) - . *darea(1:nnd))/sccsa - rhot2java=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)* - . abs(ajphiv(1:nnd))*darea(1:nnd))/sccsa - end if - end if - -c factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile - drhotpav=sqrt(8._wp_*(rhot2pav-rhotpav**2)) - drhotjava=sqrt(8._wp_*(rhot2java-rhotjava**2)) - - rhoppav=frhopol(rhotpav) - call fluxval(rhoppav,dvdrhot=dvdrhotav,ratja=ratjamx, - . ratjb=ratjbmx,ratjpl=ratjplmx) - - rhopjava=frhopol(rhotjava) - call fluxval(rhopjava,dadrhot=dadrhotava,ratja=ratjamx, - . ratjb=ratjbmx,ratjpl=ratjplmx) - - dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav) - ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava) - - call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp) - if(ieccd.ne.0) then - call profwidth(nnd,rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi) - else - rhotjfi=0.0d0 - ajmxfi=0.0d0 - ajphip=0.0d0 - drhotjfi=0.0d0 - ratjamx=1.0d0 - ratjbmx=1.0d0 - ratjplmx=1.0d0 - end if - - end subroutine postproc_profiles - - - - subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe) - use const_and_precisions, only : wp_,emn1 - use utils, only : locatex, locate, intlin, vmaxmini - implicit none -c arguments - integer :: nd - real(wp_), dimension(nd) :: xx,yy - real(wp_), intent(out) :: xpk,ypk,dxxe -c local variables - integer :: imn,imx,ipk,ie - real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2 - real(wp_) :: ypkp,ypkm -c - call vmaxmini(yy,nd,ymn,ymx,imn,imx) - ypk=0.0_wp_ - xmx=xx(imx) - xmn=xx(imn) - if (abs(ymx).gt.abs(ymn)) then - ipk=imx - ypkp=ymx - xpkp=xmx - if(abs(ymn/ymx).lt.1.0e-2_wp_) ymn=0.0_wp_ - ypkm=ymn - xpkm=xmn - else - ipk=imn - ypkp=ymn - xpkp=xmn - if(abs(ymx/ymn).lt.1.0e-2_wp_) ymx=0.0_wp_ - ypkm=ymx - xpkm=xmx - end if - if(xpkp.gt.0.0_wp_) then - xpk=xpkp - ypk=ypkp - yye=ypk*emn1 - call locatex(yy,nd,1,ipk,yye,ie) - if(ie.gt.0.and.ie.lt.nd) then - call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1) - else - rte1=0.0_wp_ - end if - call locatex(yy,nd,ipk,nd,yye,ie) - if(ie.gt.0.and.ie.lt.nd) then - call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) - else - rte2=0.0_wp_ - end if - else - ipk=2 - xpk=xx(2) - ypk=yy(2) - rte1=0.0_wp_ - yye=ypk*emn1 - call locate(yy,nd,yye,ie) - if(ie.gt.0.and.ie.lt.nd) then - call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) - else - rte2=0.0_wp_ - end if - end if - dxxe=rte2-rte1 - if(ymx.ne.0.and.ymn.ne.0) dxxe=-dxxe - - return - end -c -c -c - subroutine polarcold(sox,exf,eyif,ezf,elf,etf) - use const_and_precisions, only : wp_ - implicit none -c arguments - real(wp_) :: sox,exf,eyif,ezf,elf,etf -c local variables - real(wp_) :: anpl2,anpr2,an2,yg2,dy2,aa,e3,qq,p -c common/external functions/variables - real(wp_) :: anpl,anpr,xg,yg -c - common/nplr/anpl,anpr - common/xgxg/xg - common/ygyg/yg -c -c dcold dispersion -c dielectric tensor (transposed) -c -c exf=0.0_wp_ -c eyif=0.0_wp_ -c ezf=0.0_wp_ -c if(xg.le.0) return -c - anpl2=anpl*anpl - anpr2=anpr*anpr - an2=anpl2+anpr2 - yg2=yg**2 - dy2=1.0_wp_-yg2 - aa=1.0_wp_-xg-yg2 - e3=1.0_wp_-xg -c - if(xg.gt.0.0_wp_) then - if (anpl.ne.0.0_wp_) then - qq=xg*yg/(an2*dy2-aa) - p=(anpr2-e3)/(anpl*anpr) - ezf=1.0_wp_/sqrt(1.0_wp_+p*p*(1.0_wp_+qq*qq)) - exf=p*ezf - eyif=qq*exf - else - if(sox.lt.0.0_wp_) then - ezf=1 - exf=0 - eyif=0 - else - ezf=0 - qq=-aa/(xg*yg) - exf=1.0_wp_/sqrt(1.0_wp_+qq*qq) - eyif=qq*exf - end if - end if - elf=(anpl*ezf+anpr*exf)/sqrt(an2) - etf=sqrt(1.0_wp_-elf*elf) - else - if(sox.lt.0.0_wp_) then - ezf=1 - exf=0 - eyif=0 - else - ezf=0 - exf=0.0_wp_ - eyif=1.0_wp_ - end if - elf=0 - etf=1 - end if -c - return - end -c -c -c - subroutine pol_limit(sox,ext,eyt) - use const_and_precisions, only : wp_,ui=>im,pi - implicit none -c arguments - real(wp_), intent(in) :: sox - complex(wp_) :: ext,eyt -c local variables - real(wp_) :: anx,any,anz,xe2om,ye2om,xe2xm,ye2xm,an2,an,anxy, - . sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2,deno,denx,anpl2,dnl, - . del0,gam - complex(wp_) :: exom,eyom,exxm,eyxm -c common/external functions/variables - real(wp_) :: anpl,anpr,yg - real(wp_), dimension(3) :: bv,anv -c - common/anv/anv - common/nplr/anpl,anpr - common/ygyg/yg - common/bb/bv -c - anx=anv(1) - any=anv(2) - anz=anv(3) - anpl2=anpl*anpl - dnl=1.0_wp_-anpl2 - del0=sqrt(dnl**2+4.0_wp_*anpl2/yg**2) - ffo=0.5_wp_*yg*(dnl+del0) - ffx=0.5_wp_*yg*(dnl-del0) - an2=anx*anx+any*any+anz*anz - an=sqrt(an2) - anxy=sqrt(anx*anx+any*any) - sngam=(anz*anpl-an2*bv(3))/(an*anxy*anpr) - csgam=-(any*bv(1)-anx*bv(2))/anxy/anpr - csg2=csgam**2 - sng2=sngam**2 - ffo2=ffo*ffo - ffx2=ffx*ffx - deno=ffo2+anpl2 - denx=ffx2+anpl2 - xe2om=(ffo2*csg2+anpl2*sng2)/deno - ye2om=(ffo2*sng2+anpl2*csg2)/deno - xe2xm=(ffx2*csg2+anpl2*sng2)/denx - ye2xm=(ffx2*sng2+anpl2*csg2)/denx -c - exom=(ffo*csgam-ui*anpl*sngam)/sqrt(deno) - eyom=(-ffo*sngam-ui*anpl*csgam)/sqrt(deno) -c -c if anpl=0 the expressions for exxm and eyxm are 0/0 - if (anpl.eq.0.0_wp_) then - exxm=-ui*sngam - eyxm=-ui*csgam - else - exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx) - eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx) - end if -c - if (sox.lt.0.0_wp_) then - ext=exom - eyt=eyom - else - ext=exxm - eyt=eyxm - endif - - gam=atan(sngam/csgam)*180.0_wp_/pi - - return - end -c -c -c - subroutine stokes(ext,eyt,qq,uu,vv) - use const_and_precisions, only : wp_ - implicit none -c arguments - complex(wp_) :: ext,eyt - real(wp_) :: qq,uu,vv -c - qq=abs(ext)**2-abs(eyt)**2 - uu=2.0_wp_*dble(ext*dconjg(eyt)) - vv=2.0_wp_*dimag(ext*dconjg(eyt)) -c - end subroutine stokes -c -c -c - subroutine polellipse(qq,uu,vv,psipol,chipol) - use const_and_precisions, only : wp_,pi - implicit none -c arguments - real(wp_) :: qq,uu,vv,psipol,chipol -c -c real*8 llm,aa,bb,ell -c llm=sqrt(qq**2+uu**2) -c aa=sqrt((1+llm)/2.0_wp_) -c bb=sqrt((1-llm)/2.0_wp_) -c ell=bb/aa - psipol=0.5_wp_*atan2(uu,qq)*180.0_wp_/pi - chipol=0.5_wp_*asin(vv)*180.0_wp_/pi -c - end subroutine polellipse -c -c -c - logical function inside_plasma(rrm,zzm) - use const_and_precisions, only : wp_ - use gray_params, only : iequil - use coreprofiles, only : psdbnd - use equilibrium, only : zbinf,zbsup,equinum_psi,equian - implicit none -c arguments - real(wp_) :: rrm,zzm -c local variables - real(wp_) :: psinv -c - if(iequil.eq.1) then - call equian(rrm,zzm,psinv) - else - call equinum_psi(rrm,zzm,psinv) - end if -c - if (psinv.ge.0.0_wp_.and.psinv.lt.psdbnd) then - if (psinv.lt.1.0_wp_.and.zzm.lt.zbinf.or.zzm.gt.zbsup) then - inside_plasma=.false. - else - inside_plasma=.true. - end if - else - inside_plasma=.false. - end if -c - end function inside_plasma -c -c -c - subroutine vacuum_rt(xvstart,anv,xvend,ivac) - use const_and_precisions, only : wp_ - use reflections, only : inters_linewall,inside,rlim,zlim,nlim - use beamdata, only : dst - implicit none -c arguments - integer :: ivac - real(wp_), dimension(3) :: xvstart,anv,xvend -c local variables - integer :: i - real(wp_) :: st,rrm,zzm,smax - real(wp_), dimension(3) :: walln,xv0,anv0 - logical :: plfound -c common/external functions/variables - real(wp_) :: dstvac - logical :: inside_plasma -c - external inside_plasma -c - common/dstvac/dstvac -c -c ivac=1 plasma hit before wall reflection -c ivac=2 wall hit before plasma -c ivac=-1 vessel (and thus plasma) never crossed - -c the real argument associated to xvstart is in a common block -c used by fwork!!! - xv0=xvstart - anv0=anv - call inters_linewall(xv0/1.0e2_wp_,anv,rlim(1:nlim),zlim(1:nlim), - . nlim,smax,walln) - smax=smax*1.0e2_wp_ - rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2) - zzm=1.0e-2_wp_*xv0(3) - if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then - ! first wall interface is outside-inside - if (dot_product(walln,walln)0.and.i1 h .or. & + ah > 0.0_wp_ .and. a(jm) <= h) then + ix=ix+1 + lx(ix)=-j + end if + if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & + ah > 0.0_wp_ .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + end do + end do + + do jm=nr,mxr,nrqmax + j = jm + nrqmax + ah=a(j)-h + if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & + ah > 0.0_wp_ .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + if (ah <= 0.0_wp_ .and. a(jm) > h .or. & + ah > 0.0_wp_ .and. a(jm) <= h) then + ix=ix+1 + lx(ix)=-j + end if + end do + + do jm=1,mxr,nrqmax + j = jm + nrqmax + if (a(j) <= h .and. a(jm) > h .or. & + a(j) > h .and. a(jm) <= h) then + ix=ix+1 + lx(ix) =-j + end if + end do + + do j=2,nr + if (a(j) <= h .and. a(j-1) > h .or. & + a(j) > h .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + end do + + if(ix<=0) return + +bb: do + in=ix + jx=lx(in) + jfor=0 + lda=1 + ldb=2 + + do + if(jx<0) then + jabs=-jx + jnb = jabs - nrqmax + else + jabs=jx + jnb=jabs-1 + end if + + adn=a(jabs)-a(jnb) + if(adn/=0) px=(a(jabs)-h)/adn + kx = (jabs - 1) / nrqmax + ikx = jabs - nrqmax * kx - 1 + + if(jx<0) then + x = drgrd * ikx + y = dzgrd * (kx - px) + else + x = drgrd * (ikx - px) + y = dzgrd * kx + end if + + icount = icount + 1 + rcon(icount) = x + rqgrid(1) + zcon(icount) = y + zqgrid(1) + mpl= icount + itm = 1 + ja(1,1) = jabs + nrqmax + j=1 + + if(jx<=0) then + ja(1,1) = -jabs-1 + j=2 + end if + + ja(2,1) = -ja(1,1) + ja(3,1) = -jx + 1 - nrqmax + ja(3,2) = -jx + ja(j,2) = jabs - nrqmax + k= 3-j + ja(k,2) = 1-jabs + + if (kx<=0 .or. ikx<=0) then + lda=1 + ldb=lda + else if (ikx + 1 - nr >= 0 .and. jx <= 0) then + lda=2 + ldb=lda + else if(jfor/=0) then + lda=2 + do i=1,3 + if(jfor==ja(i,2)) then + lda=1 + exit + end if + end do + ldb=lda + end if + + flag1=.false. + aa: do k=1,3 + do l=lda,ldb + do i=1,ix + if(lx(i)==ja(k,l)) then + itm=itm+1 + inext= i + if(jfor/=0) exit aa + if(itm .gt. 3) then + flag1=.true. + exit aa + end if + end if + end do + end do + end do aa + + if(.not.flag1) then + lx(in)=0 + if(itm .eq. 1) exit + end if + + jfor=jx + jx=lx(inext) + in = inext + end do + + do + if(lx(ix)/=0) then + if(mpl>=4) then + ncon = ncon + 1 + npts(ncon) = icount - iclast + iclast = icount + end if + exit + end if + ix= ix-1 + if(ix<=0) exit bb + end do + + end do bb + + if(mpl >= 4) then + ncon = ncon + 1 + npts(ncon) = icount - iclast + iclast = icount + end if + end subroutine cniteq + + + + logical function inside_plasma(rrm,zzm) + use const_and_precisions, only : wp_, zero, one + use gray_params, only : iequil + use coreprofiles, only : psdbnd + use equilibrium, only : zbinf,zbsup,equinum_psi,equian + implicit none +! arguments + real(wp_), intent(in) :: rrm,zzm +! local variables + real(wp_) :: psinv + + if(iequil.eq.1) then + call equian(rrm,zzm,psinv) + else + call equinum_psi(rrm,zzm,psinv) + end if + + inside_plasma = (psinv >= zero .and. psinv < psdbnd) .and. & + (psinv >= one .or. (zzm >= zbinf .and. zzm <= zbsup)) + end function inside_plasma + + + + subroutine vacuum_rt(xv0,anv0,xvend,dstvac,ivac) + use const_and_precisions, only : wp_ + use reflections, only : inters_linewall,inside,rlim,zlim,nlim + use beamdata, only : dst + implicit none +! arguments + real(wp_), dimension(3), intent(in) :: xv0,anv0 + real(wp_), dimension(3), intent(out) :: xvend + real(wp_), intent(out) :: dstvac + integer, intent(out) :: ivac +! local variables + integer :: i + real(wp_) :: st,rrm,zzm,smax + real(wp_), dimension(3) :: walln + logical :: plfound +! common/external functions/variables + logical, external :: inside_plasma + +! ivac=1 plasma hit before wall reflection +! ivac=2 wall hit before plasma +! ivac=-1 vessel (and thus plasma) never crossed + + call inters_linewall(xv0/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), & + nlim,smax,walln) + smax=smax*1.0e2_wp_ + rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2) + zzm=1.0e-2_wp_*xv0(3) + if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then + ! first wall interface is outside-inside + if (dot_product(walln,walln) 1) call expinit ! ======= set environment END ====== ! ======= pre-proc prints BEGIN ====== ! print Btot=Bres ! print ne, Te, q, Jphi versus psi, rhop, rhot - if (iequil==1) then - call bres_anal(bres) - call print_prof_an - else - call bfield_res(rv,zv,size(rv),size(zv),bres) - call print_prof - end if + if (iequil<2) then + call bres_anal(bres) + call print_prof_an + else + call bfield_res(rv,zv,size(rv),size(zv),bres) + call print_prof + end if + call prfile ! ======= pre-proc prints END ====== ! ======= main loop BEGIN ====== - iox=iox0 - sox=-1.0_wp_ - if(iox.eq.2) sox=1.0_wp_ - index_rt=1 - call prfile - call paraminit - call vectinit + iox=iox0 + sox=-1.0_wp_ + if(iox==2) sox=1.0_wp_ + call vectinit(psjki,tauv,alphav,ppabs,didst,ccci,iiv) + call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri) + + psipol=psipol0 + chipol=chipol0 + call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) + call pweight(p0,p0jk) - if(igrad==0) then - call ic_rt(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, & - w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0) - else - call ic_gb(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, & - w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0) - end if + st0 = zero + if(nray>1) call print_projxyzt(st0,yw,0) ! iproj=0 ==> nfilp=8 +! test if at least part of the beam has entered the plsama + somein = .false. + istop = 0 ! beam/ray propagation + do i=1,nstep -! st0=0.0_wp_ -! if(index_rt.gt.1) st0=strfl11 - do i=1,nstep - istep=i - st=i*dst !+st0 -! advance one step - call rkint4(sox,xgcn,bres) -! calculations after one step: - call after_onestep(p0,sox,ak0,xgcn,bres,i,istop,ierr) - if(istop.eq.1) exit - end do +! advance one step with "frozen" grad(S_I) + st=i*dst+st0 + do jk=1,nray + call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk)) + end do -! postprocessing - call after_gray_integration(pec,icd,dpdv,jcd) +! update position and grad +! if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) + call gradi_upd(yw,ak0,xc,du1,gri,ggri) +! test if the beam is completely out of the plsama + allout = .true. + do jk=1,nray +! compute derivatives with updated gradient and local plasma values + xv = yw(1:3,jk) + anv = yw(4:6,jk) + call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), & + psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm) + + if( abs(anpl) > anplth1) then ! anplth1=0.99_wp_ + if(abs(anpl) <= anplth2) then ! anplth2=1.05_wp_ + ierr=97 +! igrad=0 + else + ierr=98 + istop=1 + end if + else + ierr=0 + end if + + tekev=zero + zzm = xv(3)*0.01_wp_ + ins_pl = (psinv>=zero .and. psinv=zbinf .and. zzm<=zbsup) + allout = allout .and. .not.ins_pl + somein = somein .or. ins_pl + +! compute ECRH&CD + if(ierr==0 .and. iwarm>0 .and. ins_pl) then +! print*,i,jk,rayi2jk(jk),psinv,zzm,anpl + tekev=temp(psinv) + if(tekev>zero) then + if (ieccd> 0) zeff=fzeff(psinv) + call alpha_effj(psinv,xg,yg,dens,tekev,zeff,ak0,bres,derdnm, & + anpl,anpr,sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierr) + iiv(jk)=i + end if + else + alpha=zero + anprim=zero + anprre=anpr + didp=zero + end if + + psjki(jk,i) = psinv + +! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s) + tau=tauv(jk,i-1)+0.5_wp_*(alpha+alphav(jk,i-1))*dersdst*dst + tauv(jk,i)=tau + alphav(jk,i)=alpha + pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk)) + ppabs(jk,i)=p0jk(jk)-pow + + dpdst=pow*alpha*dersdst + didst(jk,i)=didp*dpdst + ccci(jk,i)=ccci(jk,i-1)+0.5_wp_*(didst(jk,i)+didst(jk,i-1))*dst + + call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,ak0,anpl,anpr, & + anprim,dens,tekev,alphav(jk,i),tauv(jk,i),didst(jk,i),nhf,iokhawa, & + index_rt,ddr,ddi) + +! print error code + select case (ierr) + case(97) !+1 + print*,i,jk,' IERR = ', ierr,' N// = ',anpl + case(98) !+2 + print*,i,jk,' IERR = ', ierr,' N// = ',anpl + case(99) !+1*4 + print*,i,jk,' IERR = ', ierr,' Nwarm = ',anprre,anprim + case(94) !+2*4 + print*,i,jk,' IERR = ', ierr,' alpha < 0' + case(90) !+1*16 + print*,i,jk,' IERR = ', ierr,' fpp integration error' + case(91:93) !+2..4*16 + print*,i,jk,' IERR = ', ierr,' fcur integration error' + end select + + end do + +! print ray positions for j=nrayr in local reference system + if (mod(i,istpr0) == 0) then + if(nray > 1) call print_projxyzt(st,yw,0) + end if + +! test if trajectory integration must be stopped + call vmaxmin(tauv(:,i),nray,taumn,taumx) + if ((taumn > taucr) .or. (somein .and. allout)) then + pabs = sum(ppabs(:,i)) + icd = sum(ccci(:,i)) + istop = 1 + end if + if(istop == 1) exit + end do ! ======= main loop END ====== +! ======= post-proc BEGIN ====== +! print all ray positions in local reference system + if(nray > 1) call print_projxyzt(st,yw,1) +! print final results on screen + write(*,*) + write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st + write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx + write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs + write(*,'(a,f9.4)') 'I_tot (kA) = ',icd*1.0e3_wp_ + +! compute power and current density profiles for all rays + call pec_init(ipec,sqrt(psinr)) + nnd=size(rhop_tab) + allocate(jphi(nnd),pins(nnd),currins(nnd)) + call spec(psjki,ppabs,ccci,iiv,pabs,icd,dpdv,jphi,jcd,pins,currins) + + call postproc_profiles(pabs,icd,rhot_tab,dpdv,jphi, & + rhotpav,drhotpav,rhotjava,drhotjava) + +! print power and current density profiles + do i=1,nnd + write(48,'(7(1x,e16.8e3))') rhop_tab(i),rhot_tab(i), & + jphi(i),jcd(i),dpdv(i),currins(i),pins(i) + end do +! ======= post-proc END ====== + ! ======= free memory BEGIN ====== - call dealloc_beam -! call unset... + call dealloc_beam(yw,ypw,xc,du1,gri,ggri, & + psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) +! call unset_eqspl +! call unset_q +! call unset_rhospl +! call unset_prfspl + call dealloc_pec + deallocate(jphi,pins,currins) ! ======= free memory END ====== -end subroutine gray + end subroutine gray + + + subroutine vectinit(psjki,tauv,alphav,ppabs,didst,ccci,iiv) + use const_and_precisions, only : wp_, zero, one + use dispersion, only: expinit + use gray_params, only : iwarm + implicit none +! arguments + real(wp_), dimension(:,:), intent(out) :: psjki,tauv,alphav,ppabs,didst,ccci + integer, dimension(:), intent(out) :: iiv +!! common/external functions/variables +! integer :: jclosest +! real(wp_), dimension(3) :: anwcl,xwcl +! +! common/refln/anwcl,xwcl,jclosest +! +! jclosest=nrayr+1 +! anwcl(1:3)=0.0_wp_ +! xwcl(1:3)=0.0_wp_ + + psjki = zero + tauv = zero + alphav = zero + ppabs = zero + didst = zero + ccci = zero + iiv = one + + end subroutine vectinit + + + subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, & + ywrk0,ypwrk0,xc0,du10,gri,ggri) +! beam tracing initial conditions igrad=1 +! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! + use const_and_precisions, only : wp_,izero,zero,one,pi,half,two,degree,ui=>im + use math, only : catand + use gray_params, only : ipol,idst + use beamdata, only : nray,nrayr,nrayth,rwmax + implicit none +! arguments + real(wp_), dimension(3), intent(in) :: xv0c,anv0c + real(wp_), intent(in) :: ak0 + real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir + real(wp_), dimension(6,nray), intent(out) :: ywrk0,ypwrk0 + real(wp_), dimension(3,nray), intent(out) :: gri + real(wp_), dimension(3,3,nray), intent(out) :: ggri + real(wp_), dimension(3,nrayth,nrayr), intent(out) :: xc0,du10 + +! local variables + integer :: j,k,jk + real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw,alfak + real(wp_) :: wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy + real(wp_) :: rcixx,rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy + real(wp_) :: drcixx,drciyy,drcixy,dr,da,ddfu,dcsiw,detaw,dx0t,dy0t + real(wp_) :: x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2 + real(wp_) :: gxxt,gyyt,gzzt,gxyt,gxzt,gyzt,dgr2xt,dgr2yt,dgr2zt + real(wp_) :: dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy + real(wp_) :: anzt,anxt,anyt,anx,any,anz,an20,an0 + real(wp_) :: du1tx,du1ty,du1tz,denom,ddr,ddi + real(wp_), dimension(nrayr) :: uj + real(wp_), dimension(nrayth) :: sna,csa + complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1,dqi2 + complex(wp_) :: dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy + + csth=anv0c(3) + snth=sqrt(one-csth**2) + if(snth > zero) then + csps=anv0c(2)/snth + snps=anv0c(1)/snth + else + csps=one + snps=zero + end if + + phiwrad = phiw*degree + phirrad = phir*degree + csphiw = cos(phiwrad) + snphiw = sin(phiwrad) +! csphir = cos(phirrad) +! snphir = sin(phirrad) + + wwcsi = two/(ak0*wcsi**2) + wweta = two/(ak0*weta**2) + + if(phir/=phiw) then + sk = rcicsi + rcieta + sw = wwcsi + wweta + dk = rcicsi - rcieta + dw = wwcsi - wweta + ts = -(dk*sin(2*phirrad) - ui*dw*sin(2*phiwrad)) + tc = (dk*cos(2*phirrad) - ui*dw*cos(2*phiwrad)) + phic = half*catand(ts/tc) + ddd = dk*cos(2*(phirrad+phic)) - ui*dw*cos(2*(phiwrad+phic)) + sss = sk - ui*sw + qi1 = half*(sss + ddd) + qi2 = half*(sss - ddd) + rci1 = dble(qi1) + rci2 = dble(qi2) + ww1 =-dimag(qi1) + ww2 =-dimag(qi2) + else + rci1 = rcicsi + rci2 = rcieta + ww1 = wwcsi + ww2 = wweta + phic = -phiwrad + qi1 = rci1 - ui*ww1 + qi2 = rci2 - ui*ww2 + end if + +! w01=sqrt(2.0_wp_/(ak0*ww1)) +! z01=-rci1/(rci1**2+ww1**2) +! w02=sqrt(2.0_wp_/(ak0*ww2)) +! z02=-rci2/(rci2**2+ww2**2) + + qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2 + qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2 + qqxy = -(qi1 - qi2)*sin(2*phic) + wwxx = -dimag(qqxx) + wwyy = -dimag(qqyy) + wwxy = -half*dimag(qqxy) + rcixx = dble(qqxx) + rciyy = dble(qqyy) + rcixy = half* dble(qqxy) + + dqi1 = -qi1**2 + dqi2 = -qi2**2 + d2qi1 = 2*qi1**3 + d2qi2 = 2*qi2**3 + dqqxx = dqi1*cos(phic)**2 + dqi2*sin(phic)**2 + dqqyy = dqi1*sin(phic)**2 + dqi2*cos(phic)**2 + dqqxy = -(dqi1 - dqi2)*sin(2*phic) + d2qqxx = d2qi1*cos(phic)**2 + d2qi2*sin(phic)**2 + d2qqyy = d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2 + d2qqxy = -(d2qi1 - d2qi2)*sin(2*phic) + + dwwxx = -dimag(dqqxx) + dwwyy = -dimag(dqqyy) + dwwxy = -half*dimag(dqqxy) + d2wwxx = -dimag(d2qqxx) + d2wwyy = -dimag(d2qqyy) + d2wwxy = -half*dimag(d2qqxy) + drcixx = dble(dqqxx) + drciyy = dble(dqqyy) + drcixy = half* dble(dqqxy) + + if(nrayr > 1) then + dr = rwmax/dble(nrayr-1) + else + dr = one + end if + ddfu = two*dr**2/ak0 ! twodr2 = 2*dr**2 = 2*rwmax/dble(nrayr-1) + do j = 1, nrayr + uj(j) = dble(j-1) + end do + + da=2*pi/dble(nrayth) + do k=1,nrayth + alfak = (k-1)*da + sna(k) = sin(alfak) + csa(k) = cos(alfak) + end do + +! central ray + jk=1 + gri(:,1) = zero + ggri(:,:,1) = zero + + ywrk0(1:3,1) = xv0c + ywrk0(4:6,1) = anv0c + ypwrk0(1:3,1) = anv0c + ypwrk0(4:6,1) = zero + + do k=1,nrayth + dcsiw = dr*csa(k)*wcsi + detaw = dr*sna(k)*weta + dx0t = dcsiw*csphiw - detaw*snphiw + dy0t = dcsiw*snphiw + detaw*csphiw + du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu + du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu + + xc0(:,k,1) = xv0c + du10(1,k,1) = du1tx*csps + snps*du1ty*csth + du10(2,k,1) = -du1tx*snps + csps*du1ty*csth + du10(3,k,1) = -du1ty*snth + end do + ddr = zero + ddi = zero + +! loop on rays jk>1 + j=2 + k=0 + do jk=2,nray + k=k+1 + if(k > nrayth) then + j=j+1 + k=1 + end if + +! csiw=u*dcsiw +! etaw=u*detaw +! csir=x0t*csphir+y0t*snphir +! etar=-x0t*snphir+y0t*csphir + dcsiw = dr*csa(k)*wcsi + detaw = dr*sna(k)*weta + dx0t = dcsiw*csphiw - detaw*snphiw + dy0t = dcsiw*snphiw + detaw*csphiw + x0t = uj(j)*dx0t + y0t = uj(j)*dy0t + z0t = -half*(rcixx*x0t**2 + rciyy*y0t**2 + 2*rcixy*x0t*y0t) + + dx0 = x0t*csps + snps*(y0t*csth + z0t*snth) + dy0 = -x0t*snps + csps*(y0t*csth + z0t*snth) + dz0 = z0t*csth - y0t*snth + x0 = xv0c(1) + dx0 + y0 = xv0c(2) + dy0 + z0 = xv0c(3) + dz0 + + gxt = x0t*wwxx + y0t*wwxy + gyt = x0t*wwxy + y0t*wwyy + gzt = half*(x0t**2*dwwxx + y0t**2*dwwyy ) + x0t*y0t*dwwxy + gr2 = gxt*gxt + gyt*gyt + gzt*gzt + gxxt = wwxx + gyyt = wwyy + gzzt = half*(x0t**2*d2wwxx + y0t**2*d2wwyy) + x0t*y0t*d2wwxy + gxyt = wwxy + gxzt = x0t*dwwxx + y0t*dwwxy + gyzt = x0t*dwwxy + y0t*dwwyy + dgr2xt = 2*(gxt*gxxt + gyt*gxyt + gzt*gxzt) + dgr2yt = 2*(gxt*gxyt + gyt*gyyt + gzt*gyzt) + dgr2zt = 2*(gxt*gxzt + gyt*gyzt + gzt*gzzt) + dgr2x = dgr2xt*csps + snps*(dgr2yt*csth + dgr2zt*snth) + dgr2y = -dgr2xt*snps + csps*(dgr2yt*csth + dgr2zt*snth) + dgr2z = dgr2zt*csth - dgr2yt*snth + + gri(1,jk) = gxt*csps + snps*(gyt*csth + gzt*snth) + gri(2,jk) = -gxt*snps + csps*(gyt*csth + gzt*snth) + gri(3,jk) = gzt*csth - gyt*snth + ggri(1,1,jk) = gxxt*csps**2 & + + snps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) & + +2*snps*csps*(gxyt*csth + gxzt*snth) + ggri(2,1,jk) = csps*snps & + *(-gxxt+csth**2*gyyt + snth**2*gzzt + 2*csth*snth*gyzt) & + +(csps**2 - snps**2)*(snth*gxzt + csth*gxyt) + ggri(3,1,jk) = csth*snth*snps*(gzzt - gyyt) + (csth**2 - snth**2) & + *snps*gyzt + csps*(csth*gxzt - snth*gxyt) + ggri(1,2,jk) = ggri(2,1,jk) + ggri(2,2,jk) = gxxt*snps**2 & + + csps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) & + -2*snps*csps*(gxyt*csth + gxzt*snth) + ggri(3,2,jk) = csth*snth*csps*(gzzt - gyyt) + (csth**2-snth**2) & + *csps*gyzt + snps*(snth*gxyt - csth*gxzt) + ggri(1,3,jk) = ggri(3,1,jk) + ggri(2,3,jk) = ggri(3,2,jk) + ggri(3,3,jk) = gzzt*csth**2 + gyyt*snth**2 - 2*csth*snth*gyzt + + du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu + du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu + du1tz = half*uj(j)*(dx0t**2*dwwxx + dy0t**2*dwwyy + 2*dx0t*dy0t*dwwxy)/ddfu + + du10(1,k,j) = du1tx*csps + snps*(du1ty*csth + du1tz*snth) + du10(2,k,j) = -du1tx*snps + csps*(du1ty*csth + du1tz*snth) + du10(3,k,j) = du1tz*csth - du1ty*snth + + pppx = x0t*rcixx + y0t*rcixy + pppy = x0t*rcixy + y0t*rciyy + denpp = pppx*gxt + pppy*gyt + if (denpp/=zero) then + ppx = -pppx*gzt/denpp + ppy = -pppy*gzt/denpp + else + ppx = zero + ppy = zero + end if + + anzt = sqrt((one + gr2)/(one + ppx**2 + ppy**2)) + anxt = ppx*anzt + anyt = ppy*anzt + + anx = anxt*csps + snps*(anyt*csth + anzt*snth) + any =-anxt*snps + csps*(anyt*csth + anzt*snth) + anz = anzt*csth - anyt*snth + + an20 = one + gr2 + an0 = sqrt(an20) + + xc0(1,k,j) = x0 + xc0(2,k,j) = y0 + xc0(3,k,j) = z0 + + ywrk0(1,jk) = x0 + ywrk0(2,jk) = y0 + ywrk0(3,jk) = z0 + ywrk0(4,jk) = anx + ywrk0(5,jk) = any + ywrk0(6,jk) = anz + + select case(idst) + case(1) +! integration variable: c*t + denom = one + case(2) +! integration variable: Sr + denom = an20 + case default ! idst=0 +! integration variable: s + denom = an0 + end select + ypwrk0(1,jk) = anx/denom + ypwrk0(2,jk) = any/denom + ypwrk0(3,jk) = anz/denom + ypwrk0(4,jk) = dgr2x/(2*denom) + ypwrk0(5,jk) = dgr2y/(2*denom) + ypwrk0(6,jk) = dgr2z/(2*denom) + + ddr = anx**2 + any**2 + anz**2 - an20 + ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt) + end do + write(17,'(3(1x,e16.8e3))') zero,ddr,ddi + end subroutine ic_gb + + subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr) +! Runge-Kutta integrator + use const_and_precisions, only : wp_ + use gray_params, only : igrad + use beamdata, only : h,hh,h6 + implicit none + real(wp_), intent(in) :: sox,bres,xgcn + real(wp_), dimension(6), intent(inout) :: y + real(wp_), dimension(6), intent(in) :: yp + real(wp_), dimension(3), intent(in) :: dgr + real(wp_), dimension(3,3), intent(in) :: ddgr + + real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4 + real(wp_) :: gr2 + real(wp_), dimension(3) :: dgr2 + +! if(igrad.eq.1) then + gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2 + dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3)) +! end if + fk1 = yp + + yy = y + fk1*hh + call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2) + yy = y + fk2*hh + call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3) + yy = y + fk3*h + call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4) + + y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4) + end subroutine rkstep + + subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery) +! Compute right-hand side terms of the ray equations (dery) +! used in R-K integrator + use const_and_precisions, only : wp_ + use gray_params, only : idst,igrad + implicit none +! arguments + real(wp_), dimension(6), intent(in) :: y + real(wp_), intent(in) :: sox,bres,xgcn,gr2 + real(wp_), dimension(3), intent(in) :: dgr2,dgr + real(wp_), dimension(3,3), intent(in) :: ddgr + real(wp_), dimension(6), intent(out) :: dery +! local variables + real(wp_) :: psinv,dens,btot,xg,yg,anpl,anpr,ajphi + real(wp_) :: ddr,ddi,dersdst,derdnm + real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg + real(wp_), dimension(3,3) :: derbv + + xv = y(1:3) + call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg, & + ajphi) + + anv = y(4:6) + call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & + dery,anpl,anpr,ddr,ddi,dersdst,derdnm) + end subroutine rhs + + + subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & + xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm) +! Compute right-hand side terms of the ray equations (dery) +! used after full R-K step and grad(S_I) update + use gray_params, only : igrad + implicit none +! arguments + real(wp_), dimension(3), intent(in) :: xv,anv + real(wp_), dimension(3), intent(in) :: dgr + real(wp_), dimension(3,3), intent(in) :: ddgr + real(wp_), intent(in) :: sox,bres,xgcn + real(wp_), dimension(6), intent(out) :: dery + real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr + real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm +! local variables + real(wp_) :: gr2,ajphi + real(wp_), dimension(3) :: dgr2,bv,derxg,deryg + real(wp_), dimension(3,3) :: derbv + +! if(igrad == 1) then + gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2 + dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3)) +! end if + call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi) + call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & + dery,anpl,anpr,ddr,ddi,dersdst,derdnm) + end subroutine ywppla_upd + + + subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri) + use const_and_precisions, only : wp_,zero,half + use beamdata, only : nray,nrayr,nrayth,twodr2 + implicit none + real(wp_), intent(in) :: ak0 + real(wp_), dimension(6,nray), intent(in) :: ywrk + real(wp_), dimension(3,nrayth,nrayr), intent(inout) :: xc,du1 + real(wp_), dimension(3,nray), intent(out) :: gri + real(wp_), dimension(3,3,nray), intent(out) :: ggri +! local variables + real(wp_), dimension(3,nrayth,nrayr) :: xco,du1o + integer :: jk,j,jm,jp,k,km,kp + real(wp_) :: ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz + real(wp_) :: dfuu,dffiu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz + real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu + real(wp_), dimension(3,3) :: dgg,dff + +! update position and du1 vectors + xco = xc + du1o = du1 + + jk = 1 + do j=1,nrayr + do k=1,nrayth + if(j>1) jk=jk+1 + xc(1:3,k,j)=ywrk(1:3,jk) + end do + end do + +! compute grad u1 for central ray + j = 1 + jp = 2 + do k=1,nrayth + if(k == 1) then + km = nrayth + else + km = k-1 + end if + if(k == nrayth) then + kp = 1 + else + kp = k+1 + end if + dxv1 = xc(:,k ,jp) - xc(:,k ,j) + dxv2 = xc(:,kp,jp) - xc(:,km,jp) + dxv3 = xc(:,k ,j) - xco(:,k ,j) + call solg0(dxv1,dxv2,dxv3,dgu) + du1(:,k,j) = dgu + end do + gri(:,1) = zero + +! compute grad u1 and grad(S_I) for all the other rays + dfuu=twodr2/ak0 ! twodr2 = 2*dr**2 = 2*(rwmax/(nrayr-1))**2 + jm=1 + j=2 + k=0 + dffiu = dfuu + do jk=2,nray + k=k+1 + if(k > nrayth) then + jm = j + j = j+1 + k = 1 + dffiu = dfuu*jm + end if + kp = k+1 + km = k-1 + if (k == 1) then + km=nrayth + else if (k == nrayth) then + kp=1 + end if + dxv1 = xc(:,k ,j) - xc(:,k ,jm) + dxv2 = xc(:,kp,j) - xc(:,km,j) + dxv3 = xc(:,k ,j) - xco(:,k ,j) + call solg0(dxv1,dxv2,dxv3,dgu) + du1(:,k,j) = dgu + gri(:,jk) = dgu(:)*dffiu + end do + +! compute derivatives of grad u and grad(S_I) for rays jk>1 + ggri(:,:,1) = zero + jm=1 + j=2 + k=0 + dffiu = dfuu + do jk=2,nray + k=k+1 + if(k > nrayth) then + jm=j + j=j+1 + k=1 + dffiu = dfuu*jm + end if + kp=k+1 + km=k-1 + if (k == 1) then + km=nrayth + else if (k == nrayth) then + kp=1 + end if + dxv1 = xc(:,k ,j) - xc(:,k ,jm) + dxv2 = xc(:,kp,j) - xc(:,km,j) + dxv3 = xc(:,k ,j) - xco(:,k ,j) + dff(:,1) = du1(:,k ,j) - du1(:,k ,jm) + dff(:,2) = du1(:,kp,j) - du1(:,km,j) + dff(:,3) = du1(:,k ,j) - du1o(:,k ,j) + call solg3(dxv1,dxv2,dxv3,dff,dgg) + +! derivatives of u + ux = du1(1,k,j) + uy = du1(2,k,j) + uz = du1(3,k,j) + uxx = dgg(1,1) + uyy = dgg(2,2) + uzz = dgg(3,3) + uxy = (dgg(1,2) + dgg(2,1))*half + uxz = (dgg(1,3) + dgg(3,1))*half + uyz = (dgg(2,3) + dgg(3,2))*half + +! derivatives of S_I and Grad(S_I) + gx = ux*dffiu + gy = uy*dffiu + gz = uz*dffiu + gxx = dfuu*ux*ux + dffiu*uxx + gyy = dfuu*uy*uy + dffiu*uyy + gzz = dfuu*uz*uz + dffiu*uzz + gxy = dfuu*ux*uy + dffiu*uxy + gxz = dfuu*ux*uz + dffiu*uxz + gyz = dfuu*uy*uz + dffiu*uyz + + ggri(1,1,jk)=gxx + ggri(2,1,jk)=gxy + ggri(3,1,jk)=gxz + ggri(1,2,jk)=gxy + ggri(2,2,jk)=gyy + ggri(3,2,jk)=gyz + ggri(1,3,jk)=gxz + ggri(2,3,jk)=gyz + ggri(3,3,jk)=gzz + end do + + end subroutine gradi_upd + + subroutine solg0(dxv1,dxv2,dxv3,dgg) +! solution of the linear system of 3 eqs : dgg . dxv = dff +! input vectors : dxv1, dxv2, dxv3, dff +! output vector : dgg +! dff=(1,0,0) + use const_and_precisions, only : wp_ + implicit none +! arguments + real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 + real(wp_), dimension(3), intent(out) :: dgg +! local variables + real(wp_) :: denom,aa1,aa2,aa3 + + aa1 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3)) + aa2 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3)) + aa3 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3)) + + denom = dxv1(1)*aa1 - dxv2(1)*aa2 + dxv3(1)*aa3 + + dgg(1) = aa1/denom + dgg(2) = -(dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3))/denom + dgg(3) = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2))/denom + end subroutine solg0 + + subroutine solg3(dxv1,dxv2,dxv3,dff,dgg) +! rhs "matrix" dff, result in dgg + use const_and_precisions, only : wp_ + implicit none +! arguments + real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 + real(wp_), dimension(3,3), intent(in) :: dff + real(wp_), dimension(3,3), intent(out) :: dgg +! local variables + real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 + + a11 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3)) + a21 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3)) + a31 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3)) + + a12 = (dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3)) + a22 = (dxv1(1)*dxv3(3) - dxv3(1)*dxv1(3)) + a32 = (dxv1(1)*dxv2(3) - dxv2(1)*dxv1(3)) + + a13 = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2)) + a23 = (dxv1(1)*dxv3(2) - dxv3(1)*dxv1(2)) + a33 = (dxv1(1)*dxv2(2) - dxv2(1)*dxv1(2)) + + denom = dxv1(1)*a11 - dxv2(1)*a21 + dxv3(1)*a31 + + dgg(:,1) = ( dff(:,1)*a11 - dff(:,2)*a21 + dff(:,3)*a31)/denom + dgg(:,2) = (-dff(:,1)*a12 + dff(:,2)*a22 - dff(:,3)*a32)/denom + dgg(:,3) = ( dff(:,1)*a13 - dff(:,2)*a23 + dff(:,3)*a33)/denom + end subroutine solg3 + + + subroutine plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv, & + xg,yg,derxg,deryg,ajphi) + use const_and_precisions, only : wp_,zero,pi,ccj=>mu0inv + use gray_params, only : iequil + use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi + use coreprofiles, only : density + implicit none +! arguments + real(wp_), dimension(3), intent(in) :: xv + real(wp_), intent(in) :: xgcn,bres + real(wp_), intent(out) :: psinv,dens,btot,xg,yg + real(wp_), dimension(3), intent(out) :: bv,derxg,deryg + real(wp_), dimension(3,3), intent(out) :: derbv +! local variables + integer :: iv,jv + real(wp_) :: xx,yy,zz + real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm + real(wp_), dimension(3) :: dbtot,bvc + real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv + real(wp_) :: brr,bphi,bzz,ajphi,dxgdpsi + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddenspsin + + xg = zero + yg = 99._wp_ + psinv = -1._wp_ + dens = zero + btot = zero + ajphi = zero + derxg = zero + deryg = zero + bv = zero + derbv = zero + + if(iequil==0) return + + dbtot = zero + dbv = zero + dbvcdc = zero + dbvcdc = zero + dbvdc = zero + + xx = xv(1) + yy = xv(2) + zz = xv(3) + +! cylindrical coordinates + rr2 = xx**2 + yy**2 + rr = sqrt(rr2) + csphi = xx/rr + snphi = yy/rr + + bv(1) = -snphi*sgnbphi + bv(2) = csphi*sgnbphi + +! convert from cm to meters + zzm = 1.0e-2_wp_*zz + rrm = 1.0e-2_wp_*rr + + if(iequil==1) then + call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, & + ddpsidrr,ddpsidzz,ddpsidrz) + else + call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz) + call equinum_fpol(psinv,fpolv,dfpv) + end if + +! compute yg and derivative + if(psinv < zero) then + bphi = fpolv/rrm + btot = abs(bphi) + yg = btot/bres + return + end if + +! compute xg and derivative + call density(psinv,dens,ddenspsin) + xg = xgcn*dens + dxgdpsi = xgcn*ddenspsin/psia + +! B = f(psi)/R e_phi+ grad psi x e_phi/R + bphi = fpolv/rrm + brr =-dpsidz/rrm + bzz = dpsidr/rrm + +! bvc(i) = B_i in cylindrical coordinates + bvc(1) = brr + bvc(2) = bphi + bvc(3) = bzz + +! bv(i) = B_i in cartesian coordinates + bv(1)=bvc(1)*csphi - bvc(2)*snphi + bv(2)=bvc(1)*snphi + bvc(2)*csphi + bv(3)=bvc(3) + +! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) + dbvcdc(1,1) = -ddpsidrz/rrm - brr/rrm + dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm + dbvcdc(3,1) = ddpsidrr/rrm - bzz/rrm + dbvcdc(1,3) = -ddpsidzz/rrm + dbvcdc(2,3) = dfpv*dpsidz/rrm + dbvcdc(3,3) = ddpsidrz/rrm + +! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) + dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi + dbvdc(2,1) = dbvcdc(1,1)*snphi + dbvcdc(2,1)*csphi + dbvdc(3,1) = dbvcdc(3,1) + dbvdc(1,2) = -bv(2) + dbvdc(2,2) = bv(1) + dbvdc(3,2) = dbvcdc(3,2) + dbvdc(1,3) = dbvcdc(1,3)*csphi - dbvcdc(2,3)*snphi + dbvdc(2,3) = dbvcdc(1,3)*snphi + dbvcdc(2,3)*csphi + dbvdc(3,3) = dbvcdc(3,3) + + drrdx = csphi + drrdy = snphi + dphidx = -snphi/rrm + dphidy = csphi/rrm + +! dbv(iv,jv) = d Bcart(iv) / dxvcart(jv) + dbv(:,1) = drrdx*dbvdc(:,1) + dphidx*dbvdc(:,2) + dbv(:,2) = drrdy*dbvdc(:,1) + dphidy*dbvdc(:,2) + dbv(:,3) = dbvdc(:,3) + +! B magnitude and derivatives + b2tot = bv(1)**2 + bv(2)**2 + bv(3)**2 + btot = sqrt(b2tot) + + dbtot = (bv(1)*dbv(1,:) + bv(2)*dbv(2,:) + bv(3)*dbv(3,:))/btot + + yg = btot/Bres + +! convert spatial derivatives from dummy/m -> dummy/cm +! to be used in rhs + +! bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z + deryg = 1.0e-2_wp_*dbtot/Bres + bv = bv/btot + do jv=1,3 + derbv(:,jv) = 1.0e-2_wp_*(dbv(:,jv) - bv(:)*dbtot(jv))/btot + end do + + derxg(1) = 1.0e-2_wp_*drrdx*dpsidr*dxgdpsi + derxg(2) = 1.0e-2_wp_*drrdy*dpsidr*dxgdpsi + derxg(3) = 1.0e-2_wp_*dpsidz *dxgdpsi + +! current density computation in Ampere/m^2, ccj==1/mu_0 + ajphi = ccj*(dbvcdc(1,3) - dbvcdc(3,1)) +! ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3)) +! ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2)) + + end subroutine plas_deriv + + + subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & + dery,anpl,anpr,ddr,ddi,dersdst,derdnm) + use const_and_precisions, only : wp_,zero,one,half,two + use gray_params, only : idst,igrad + implicit none +! arguments + real(wp_), intent(in) :: xg,yg,gr2,sox + real(wp_), intent(out) :: anpl,anpr,ddr,ddi,derdnm,dersdst + real(wp_), dimension(3), intent(in) :: anv,bv,derxg,deryg + real(wp_), dimension(3), intent(in) :: dgr2,dgr + real(wp_), dimension(3,3), intent(in) :: ddgr,derbv + real(wp_), dimension(6), intent(out) :: dery +! local variables + integer :: iv,jv + real(wp_) :: yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,an2,an2s + real(wp_) :: dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel + real(wp_) :: derdom,dfdiadnpl,dfdiadxg,dfdiadyg,fdia,bdotgr !,vgm + real(wp_), dimension(3) :: derdxv,danpldxv,derdnv,dbgr !,vgv + + an2 = anv(1)*anv(1) + anv(2)*anv(2) + anv(3)*anv(3) + anpl = anv(1)*bv(1) + anv(2)*bv(2) + anv(3)*bv(3) + + anpl2 = anpl**2 + dnl = one - anpl2 + anpr2 = max(an2-anpl2,zero) + anpr = sqrt(anpr2) + yg2 = yg**2 + + an2s = one + dan2sdxg = zero + dan2sdyg = zero + dan2sdnpl = zero + del = zero + fdia = zero + dfdiadnpl = zero + dfdiadxg = zero + dfdiadyg = zero + + duh = one - xg - yg2 + if(xg > zero) then + del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2) + an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh + + dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 & + + sox*xg*anpl2/(del*duh) - one + dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 & + + two*sox*xg*(one - xg)*anpl2/(yg*del*duh) + dan2sdnpl = - xg*yg2*anpl/duh & + - sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh) + + if(igrad > 0) then + ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) & + - yg2*dnl**3)/yg2/del**3 + fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh + derdel = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2) & + - dnl**2*(one + 3.0_wp_*anpl2)*yg2 + derdel = 4.0_wp_*derdel/(yg*del)**5 + ddelnpl2y = two*(one - xg)*derdel + ddelnpl2x = yg*derdel + dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) & + /(yg2*del**5) + dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) & + *ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2) + dfdiadyg = - two*yg*xg*(one - xg)/duh**2 & + - sox*xg*yg*(two*(one - xg)*ddelnpl2 & + + yg*duh*ddelnpl2y)/(two*duh**2) + end if + end if + + bdotgr = bv(1)*dgr(1) + bv(2)*dgr(2) + bv(3)*dgr(3) + do iv=1,3 + dbgr(iv) = dgr(1)*derbv(1,iv) + bv(1)*ddgr(1,iv) & + + dgr(2)*derbv(2,iv) + bv(2)*ddgr(2,iv) & + + dgr(3)*derbv(3,iv) + bv(3)*ddgr(3,iv) + danpldxv(iv) = anv(1)*derbv(1,iv) + anv(2)*derbv(2,iv) + anv(3)*derbv(3,iv) + end do + + derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl + & + igrad*dgr2) & + + fdia*bdotgr*dbgr + half*bdotgr**2 & + *(derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl) + derdnv = two*anv + (half*bdotgr**2*dfdiadnpl - dan2sdnpl)*bv + + derdnm = sqrt(derdnv(1)**2 + derdnv(2)**2 + derdnv(3)**2) + + derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl & + + two*igrad*gr2 - bdotgr**2*(fdia + xg*dfdiadxg & + + half*yg*dfdiadyg & + + half*anpl*dfdiadnpl) + + if (idst == 0) then +! integration variable: s + denom = derdnm + else if (idst == 1) then +! integration variable: c*t + denom = -derdom + else +! integration variable: Sr + denom = anv(1)*derdnv(1) + anv(2)*derdnv(2) + anv(3)*derdnv(3) + end if + +! coefficient for integration in s +! ds/dst, where st is the integration variable + dersdst = derdnm/denom + +! rhs vector + dery(1:3) = derdnv(:)/denom + dery(4:6) = -derdxv(:)/denom + +! vgv : ~ group velocity +! vgm=0 +! do iv=1,3 +! vgv(iv)=-derdnv(iv)/derdom +! vgm=vgm+vgv(iv)**2 +! end do +! vgm=sqrt(vgm) + +! ddr : dispersion relation (real part) +! ddi : dispersion relation (imaginary part) + ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia) + ddi = derdnv(1)*dgr(1) + derdnv(2)*dgr(2) + derdnv(3)*dgr(3) + + end subroutine disp_deriv + + + subroutine alpha_effj(psinv,xg,yg,dens,tekev,zeff,ak0,bres,derdnm,anpl,anpr, & + sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,ierr) + use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_ + use gray_params, only : iwarm,ilarm,ieccd,imx + use equilibrium, only : rmaxis,sgnbphi +! use beamdata, only : psjki,tauv,alphav,ppabs,didst,ccci,tau1v,pdjki,currj + use dispersion, only : harmnumber, warmdisp + use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl + use magsurf_data, only : fluxval + implicit none +! arguments + real(wp_),intent(in) ::psinv,ak0,bres + real(wp_),intent(in) :: xg,yg,tekev,dens,zeff,anpl,anpr,derdnm,sox + real(wp_),intent(out) :: anprre,anprim,alpha,didp + integer, intent(out) :: nhmin,nhmax,iokhawa + integer, intent(out) :: ierr +! local constants + real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ +! local variables + real(wp_) :: rbavi,rrii,rhop + integer :: lrm,ithn + real(wp_) :: amu,ratiovgr,rbn,rbx + real(wp_) :: cst2,bmxi,bmni,fci + real(wp_), dimension(:), allocatable :: eccdpar + real(wp_) :: effjcd,effjcdav,akim,btot + complex(wp_) :: ex,ey,ez + + alpha=zero + anprim=zero + anprre=zero + didp=zero + + if(tekev>zero) then +! absorption computation + amu=mc2/tekev + call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) + if(nhmin.gt.0) then + lrm=max(ilarm,nhmax) + call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, & + iwarm,imx,ex,ey,ez) + akim=ak0*anprim + ratiovgr=2.0_wp_*anpr/derdnm!*vgm + alpha=2.0_wp_*akim*ratiovgr + if(alpha: effjcdav [A m/W ] + if(ieccd>0) then +! current drive computation + ithn=1 + if(lrm>nhmin) ithn=2 + rhop=sqrt(psinv) + call fluxval(rhop,rri=rrii,rbav=rbavi,bmn=bmni,bmx=bmxi,fc=fci) + btot=yg*bres + rbn=btot/bmni + rbx=btot/bmxi + + select case(ieccd) + case(1) +! 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) +! 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 +! 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) + effjcdav=rbavi*effjcd + didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii) + end if + end if + end if + end subroutine alpha_effj + + + + subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0) + use const_and_precisions, only : wp_,degree,zero,one,half,im + use beamdata, only : nray,nrayth + use equilibrium, only : bfield + use gray_params, only : ipol + use polarization, only : pol_limit, polellipse, stokes_ce, stokes_ell + implicit none +! arguments + real(wp_), dimension(6,nray), intent(in) :: ywrk0 + real(wp_), intent(in) :: sox,bres + real(wp_), intent(inout) :: psipol0, chipol0 + complex(wp_), dimension(nray), intent(out) :: ext0, eyt0 +! local variables + integer :: j,k,jk + real(wp_), dimension(3) :: xmv, anv, bv + real(wp_) :: rm, csphi, snphi, bphi, br, bz, qq, uu, vv, deltapol + + j=1 + k=0 + do jk=1,nray + k=k+1 + if(jk == 2 .or. k > nrayth) then + j=j+1 + k=1 + end if + + if(ipol == 0) then + xmv=ywrk0(1:3,jk)*0.01_wp_ ! convert from cm to m + anv=ywrk0(4:6,jk) + rm=sqrt(xmv(1)**2+xmv(2)**2) + csphi=xmv(1)/rm + snphi=xmv(2)/rm + call bfield(rm,xmv(3),bphi,br,bz) +! bv(i) = B_i in cartesian coordinates + bv(1)=br*csphi-bphi*snphi + bv(2)=br*snphi+bphi*csphi + bv(3)=bz + call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk)) + + if (jk == 1) then + call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv) + call polellipse(qq,uu,vv,psipol0,chipol0) + psipol0=psipol0/degree ! convert from rad to degree + chipol0=chipol0/degree + end if + else + call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv) + if(qq**2 < one) then + deltapol=asin(vv/sqrt(one - qq**2)) + ext0(jk)= sqrt(half*(one + qq)) + eyt0(jk)= sqrt(half*(one - qq))*exp(-im*deltapol) + else + ext0(jk)= one + eyt0(jk)= zero + end if + endif + end do + end subroutine set_pol + + subroutine print_output(i,jk,st,qj,xv,psinv,btot,ak0,anpl,anpr,anprim, & + dens,tekev,alpha,tau,didst,nhf,iokhawa,index_rt,ddr,ddi) + use const_and_precisions, only : degree,zero,one + use equilibrium, only : frhotor + use gray_params, only : istpl0 + use beamdata, only : nray,nrayth + implicit none +! arguments + integer, intent(in) :: i,jk,nhf,iokhawa,index_rt + real(wp_), dimension(3), intent(in) :: xv + real(wp_), intent(in) :: st,qj,psinv,btot,ak0,anpl,anpr,anprim + real(wp_), intent(in) :: dens,tekev,alpha,tau,didst,ddr,ddi +! local variables + real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,dids + integer :: k + + stm=st*1.0e-2_wp_ + xxm=xv(1)*1.0e-2_wp_ + yym=xv(2)*1.0e-2_wp_ + zzm=xv(3)*1.0e-2_wp_ + rrm=sqrt(xxm**2 + yym**2) + +! central ray only begin +! print dIds in A/m/W, ki in m^-1 + if(jk.eq.1) then + phideg=atan2(yym,xxm)/degree + if(psinv>=zero .and. psinv<=one) then + rhot=frhotor(psinv) + else + rhot=1.0_wp_ + end if + akim=anprim*ak0*1.0e2_wp_ + pt=exp(-tau) + dids=didst*1.0e2_wp_/qj + + write(4,'(30(1x,e16.8e3))') stm,rrm,zzm,phideg,psinv,rhot,dens,tekev, & + btot,anpr,anpl,akim,alpha,tau,pt,dids,dble(nhf),dble(iokhawa), & + dble(index_rt),ddr + end if +! central ray only end + +! print conservation of dispersion relation + if(jk==nray) write(17,'(30(1x,e16.8e3))') st,ddr,ddi + +! print outer trajectories + if(mod(i,istpl0)==0) then + k = jk + nrayth - nray + if(k>0) then + write(33,'(2i5,16(1x,e16.8e3))') i,k,stm,xxm,yym,rrm,zzm, & + psinv,tau,anpl,alpha,dble(index_rt) + end if + end if + end subroutine print_output end module graycore diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index fa69d23..926eed6 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -105,6 +105,98 @@ contains + subroutine contours_psi(h,rup,zup,rlw,zlw,rcn,zcn,ipr) + use const_and_precisions, only : wp_,pi + use equilibrium, only : psiant,psinop,nsr,nsz,cc=>cceq,tr,tz,kspl, & + points_tgo + use dierckx, only : profil,sproota + use reflections, only : rwallm + implicit none +! local constants + integer, parameter :: mest=4 +! arguments + integer, intent(in) :: ipr + real(wp_), intent(in) :: h + real(wp_), intent(inout) :: rup,zup,rlw,zlw + real(wp_), dimension(npoints), intent(out) :: rcn,zcn +! local variables + integer :: np,info,ic,ier,ii,iopt,m + real(wp_) :: ra,rb,za,zb,th,zc,val + real(wp_), dimension(mest) :: zeroc + real(wp_), dimension(nsr) :: czc + + np=(npoints-1)/2 + + ra=rup + rb=rlw + za=zup + zb=zlw + call points_tgo(ra,za,rup,zup,h,info) + call points_tgo(rb,zb,rlw,zlw,h,info) + + th=pi/dble(np) + rcn(1)=rlw + zcn(1)=zlw + rcn(npoints)=rlw + zcn(npoints)=zlw + rcn(np+1)=rup + zcn(np+1)=zup + do ic=2,np + zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ + iopt=1 + call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) + if(ier.gt.0) print*,' profil =',ier + val=h*psiant+psinop + call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) + if (zeroc(1).gt.rwallm) then + rcn(ic)=zeroc(1) + rcn(npoints+1-ic)=zeroc(2) + else + rcn(ic)=zeroc(2) + rcn(npoints+1-ic)=zeroc(3) + end if + zcn(ic)=zc + zcn(npoints+1-ic)=zc + end do + + if (ipr.gt.0) then + do ii=1,npoints + write(71,'(i6,12(1x,e12.5))') ii,h,rcn(ii),zcn(ii) + end do + write(71,*) + write(71,*) + end if + end subroutine contours_psi + + + + subroutine contours_psi_an(h,rcn,zcn,ipr) + use const_and_precisions, only : wp_,pi + use equilibrium, only : frhotor,aminor,rmaxis,zmaxis + implicit none +! arguments + integer :: ipr + real(wp_) :: h + real(wp_), dimension(npoints) :: rcn,zcn +! local variables + integer :: np,ic + real(wp_) :: rn,th + + np=(npoints-1)/2 + th=pi/dble(np) + rn=frhotor(sqrt(h)) + + do ic=1,npoints + zcn(ic)=zmaxis+aminor*rn*sin(th*(ic-1)) + rcn(ic)=rmaxis+aminor*rn*cos(th*(ic-1)) + + if (ipr.gt.0) write(71,'(i6,12(1x,e12.5))') ic,h,rcn(ic),zcn(ic) + end do + if (ipr.gt.0) write(71,*) + end subroutine contours_psi_an + + + subroutine flux_average use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use gray_params, only : iequil diff --git a/src/main.f90 b/src/main.f90 index ccbab13..0c697e1 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -104,7 +104,7 @@ program gray_main allocate(dpdv(outp%nrho),jcd(outp%nrho)) call gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & psrad,terad,derad,zfc,prfp, rlim,zlim, & - p0mw,fghz,alpha0,beta0,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir,iox0, & + p0mw,fghz,alpha0,beta0,(/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,iox0, & psipol0,chipol0, dpdv,jcd,pec,icd, outp,rtrp,hcdp,ierr) ! ======================================================================== diff --git a/src/pec.f90 b/src/pec.f90 new file mode 100644 index 0000000..428b133 --- /dev/null +++ b/src/pec.f90 @@ -0,0 +1,385 @@ +module pec + use const_and_precisions, only : wp_,zero,one + implicit none + real(wp_), dimension(:), allocatable, save :: rhop_tab,rhot_tab + real(wp_), dimension(:), allocatable, save :: rtabpsi1 + real(wp_), dimension(:), allocatable, save :: dvol,darea + real(wp_), dimension(:), allocatable, save :: ratjav,ratjbv,ratjplv + +contains + + subroutine pec_init(ipec,rt_in) + use equilibrium, only : frhotor,frhopol + use gray_params, only : nnd + use magsurf_data, only : fluxval + implicit none +! arguments + integer, intent(in) :: ipec + real(wp_), dimension(nnd), intent(in) :: rt_in +! local variables + integer :: it + real(wp_) :: drt,rt,rt1,rhop1 + real(wp_) :: ratjai,ratjbi,ratjpli + real(wp_) :: voli0,voli1,areai0,areai1 + +! ipec positive build equidistant grid dimension nnd +! ipec negative read input grid + +! ipec=+/-1 rho_pol grid +! ipec=+/-2 rho_tor grid + call dealloc_pec + allocate(rhop_tab(nnd),rhot_tab(nnd),rtabpsi1(0:nnd),dvol(nnd),darea(nnd), & + ratjav(nnd),ratjbv(nnd),ratjplv(nnd)) + + voli0 = zero + areai0 = zero + rtabpsi1(0) = zero + + do it=1,nnd + if(ipec > 0) then +! build equidistant radial grid + drt = one/dble(nnd-1) + rt = dble(it-1)*drt + else +! read radial grid from input + rt = rt_in(it) + drt = (rt_in(it+1)-rt)/2.0_wp_ !!!!! WARNING !!!!! non funziona per it==nnd + end if +! radial coordinate of i-(i+1) interval mid point + if(it < nnd) then + rt1 = rt + drt/2.0_wp_ + else + rt1 = one + end if + if (abs(ipec) == 1) then + rhop_tab(it) = rt + rhot_tab(it) = frhotor(rt) + rhop1 = rt1 + else + rhot_tab(it) = rt + rhop_tab(it) = frhopol(rt) + rhop1 = frhopol(rt1) + end if +! psi grid at mid points, dimension nnd+1, for use in pec_tab + rtabpsi1(it) = rhop1**2 + + call fluxval(rhop1,area=areai1,vol=voli1) + dvol(it) = abs(voli1 - voli0) + darea(it) = abs(areai1 - areai0) + voli0 = voli1 + areai0 = areai1 + + call fluxval(rhop_tab(it),ratja=ratjai,ratjb=ratjbi,ratjpl=ratjpli) + ratjav(it) = ratjai + ratjbv(it) = ratjbi + ratjplv(it) = ratjpli + end do + end subroutine pec_init + + + subroutine spec(psjki,ppabs,ccci,iiv,pabs,currt,dpdv,ajphiv,ajcd,pins,currins) + use gray_params, only : nnd + use beamdata, only : nray,nstep + implicit none +! local constants + real(wp_), parameter :: rtbc=one +! arguments + real(wp_), dimension(nray,nstep), intent(in) :: psjki,ppabs,ccci + integer, dimension(nray), intent(in) :: iiv + real(wp_), intent(in) :: pabs,currt + real(wp_), dimension(nnd), intent(out) :: dpdv,ajphiv,ajcd,pins,currins +! local variables + integer :: i,ii,jk + real(wp_) :: spds,sccs,facpds,facjs + real(wp_), dimension(nstep):: xxi,ypt,yamp + real(wp_), dimension(nnd) :: wdpdv,wajphiv + +! calculation of dP and dI over radial grid + dpdv=zero + ajphiv=zero + do jk=1,nray + ii=iiv(jk) + if (ii < nstep ) then + if(psjki(jk,ii+1) /= zero) ii=ii+1 + end if + xxi=zero + ypt=zero + yamp=zero + do i=1,ii + xxi(i)=abs(psjki(jk,i)) + if(xxi(i) <= one) then + ypt(i)=ppabs(jk,i) + yamp(i)=ccci(jk,i) + end if + end do + call pec_tab(xxi,ypt,yamp,ii,rtabpsi1,wdpdv,wajphiv) + dpdv = dpdv + wdpdv + ajphiv = ajphiv + wajphiv + end do + +! here dpdv is still dP (not divided yet by dV) +! here ajphiv is still dI (not divided yet by dA) + spds=zero + sccs=zero + do i=1,nnd + spds=spds+dpdv(i) + sccs=sccs+ajphiv(i) + pins(i)=spds + currins(i)=sccs + end do + + facpds=one + facjs=one + if(spds > zero) facpds=pabs/spds + if(sccs /= zero) facjs=currt/sccs + + dpdv=facpds*(dpdv/dvol) + ajphiv=facjs*(ajphiv/darea) + ajcd=ajphiv*ratjbv + +! now dpdv is dP/dV [MW/m^3] +! now ajphiv is J_phi=dI/dA [MA/m^2] + end subroutine spec + + + + subroutine pec_tab(xxi,ypt,yamp,ii,xtab1,wdpdv,wajphiv) +! Power and current projected on psi grid - mid points + use const_and_precisions, only : wp_,one,zero + use gray_params, only : nnd + use utils, only : locatex,intlin +! arguments + integer, intent(in) :: ii + real(wp_), dimension(ii), intent(in) :: xxi,ypt,yamp + real(wp_), dimension(0:nnd), intent(in) :: xtab1 + real(wp_), dimension(nnd), intent(out) :: wdpdv,wajphiv +! local variables + integer, parameter :: llmx = 21 + integer, dimension(llmx) ::isev + real(wp_) :: ppa1,ppa2,cci1,cci2,dppa,didst,rt1 + integer :: i,is,ise0,idecr,iise0,iise,iis,iis1 + integer :: ind1,ind2,iind,ind,indi,itb1 + + isev = 0 + ise0 = 0 + idecr = -1 + is = 1 + wdpdv = zero + wajphiv = zero + do i=1,ii + if(ise0 == 0) then + if(xxi(i) < one) then + ise0 = i + isev(is) = i - 1 + is = is + 1 + end if + else + if (idecr == -1) then + if(xxi(i) > xxi(i-1)) then + isev(is) = i - 1 + is = is + 1 + idecr = 1 + end if + else + if(xxi(i) > one) exit + if(xxi(i) < xxi(i-1)) then + isev(is) = i - 1 + is = is + 1 + idecr = -1 + end if + end if + end if + end do + + isev(is) = i-1 + ppa1 = zero + cci1 = zero + + do iis=1,is-1 + iis1 = iis + 1 + iise0 = isev(iis) + iise = isev(iis1) + if (mod(iis,2) /= 0) then + idecr = -1 + ind1 = nnd + ind2 = 2 + iind = -1 + else + idecr = 1 + ind1 = 1 + ind2 = nnd + iind = 1 + end if + do ind=ind1,ind2,iind + indi = ind + if (idecr == -1) indi = ind - 1 + rt1 = xtab1(indi) + call locatex(xxi,iise,iise0,iise,rt1,itb1) + if(itb1 >= iise0 .and. itb1 < iise) then + call intlin(xxi(itb1), ypt(itb1),xxi(itb1+1), ypt(itb1+1),rt1,ppa2) + call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1),rt1,cci2) + dppa = ppa2 - ppa1 + didst = cci2 - cci1 + wdpdv(ind) = wdpdv(ind) + dppa + wajphiv(ind) = wajphiv(ind) + didst + ppa1 = ppa2 + cci1 = cci2 + end if + end do + end do + end subroutine pec_tab + + + subroutine postproc_profiles(pabs,currt,rhot_tab,dpdv,ajphiv, & + rhotpav,drhotpav,rhotjava,drhotjava) +! radial average values over power and current density profile + use const_and_precisions, only : pi + use gray_params, only : nnd + use equilibrium, only : frhopol + use magsurf_data, only : fluxval + implicit none + real(wp_), intent(in) :: pabs,currt + real(wp_), dimension(nnd), intent(in) :: rhot_tab + real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv + real(wp_), intent(out) :: rhotpav,rhotjava + real(wp_), intent(out) :: drhotpav,drhotjava + real(wp_) :: rhopjava,rhoppav + real(wp_) :: dpdvp,dpdvmx,rhotp,drhotp + real(wp_) :: ajphip,ajmxfi,rhotjfi,drhotjfi + real(wp_) :: ratjamx,ratjbmx,ratjplmx + + real(wp_) :: sccsa + real(wp_) :: rhotjav,rhot2pav,rhot2java,dvdrhotav,dadrhotava + + rhotpav=zero + rhot2pav=zero + rhotjav=zero + rhotjava=zero + rhot2java=zero + + if (pabs > zero) then + rhotpav = sum(rhot_tab *dpdv*dvol)/pabs + rhot2pav = sum(rhot_tab**2*dpdv*dvol)/pabs + end if + + if (abs(currt) > zero) then + rhotjav = sum(rhot_tab*ajphiv*darea)/currt + end if + sccsa = sum(abs(ajphiv)*darea) + if (sccsa > zero) then + rhotjava = sum(rhot_tab *abs(ajphiv)*darea)/sccsa + rhot2java = sum(rhot_tab**2*abs(ajphiv)*darea)/sccsa + end if + +! factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile + drhotpav = sqrt(8._wp_*(rhot2pav -rhotpav**2)) + drhotjava = sqrt(8._wp_*(rhot2java-rhotjava**2)) + + rhoppav = frhopol(rhotpav) + rhopjava = frhopol(rhotjava) + + if (pabs > zero) then + call fluxval(rhoppav,dvdrhot=dvdrhotav) + dpdvp = pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav) + call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp) + else + dpdvp = zero + rhotp = zero + dpdvmx = zero + drhotp = zero + end if + + if (sccsa > zero) then + call fluxval(rhopjava,dadrhot=dadrhotava,ratja=ratjamx,ratjb=ratjbmx, & + ratjpl=ratjplmx) + ajphip = currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava) + call profwidth(nnd,rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi) + else + ajphip = zero + rhotjfi = zero + ajmxfi = zero + drhotjfi = zero + end if + end subroutine postproc_profiles + + + + subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe) + use const_and_precisions, only : wp_,emn1 + use utils, only : locatex, locate, intlin, vmaxmini + implicit none +! arguments + integer :: nd + real(wp_), dimension(nd) :: xx,yy + real(wp_), intent(out) :: xpk,ypk,dxxe +! local variables + integer :: imn,imx,ipk,ie + real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2 + real(wp_) :: ypkp,ypkm + + call vmaxmini(yy,nd,ymn,ymx,imn,imx) + ypk = zero + xmx = xx(imx) + xmn = xx(imn) + if (abs(ymx) > abs(ymn)) then + ipk = imx + ypkp = ymx + xpkp = xmx + if(abs(ymn/ymx) < 1.0e-2_wp_) ymn = 0.0_wp_ + ypkm = ymn + xpkm = xmn + else + ipk = imn + ypkp = ymn + xpkp = xmn + if(abs(ymx/ymn) < 1.0e-2_wp_) ymx = 0.0_wp_ + ypkm = ymx + xpkm = xmx + end if + if(xpkp > zero) then + xpk = xpkp + ypk = ypkp + yye = ypk*emn1 + call locatex(yy,nd,1,ipk,yye,ie) + if(ie > 0 .and. ie < nd) then + call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1) + else + rte1 = zero + end if + call locatex(yy,nd,ipk,nd,yye,ie) + if(ie > 0 .and. ie < nd) then + call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) + else + rte2 = zero + end if + else + ipk=2 + xpk=xx(2) + ypk=yy(2) + rte1=0.0_wp_ + yye=ypk*emn1 + call locate(yy,nd,yye,ie) + if(ie > 0 .and. ie < nd) then + call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) + else + rte2 = zero + end if + end if + dxxe = rte2 - rte1 + if(ymx /= zero .and. ymn /= zero) dxxe = -dxxe + end subroutine profwidth + + subroutine dealloc_pec + implicit none + + if (allocated(rhop_tab)) deallocate(rhop_tab) + if (allocated(rhot_tab)) deallocate(rhot_tab) + if (allocated(rtabpsi1)) deallocate(rtabpsi1) + if (allocated(dvol)) deallocate(dvol) + if (allocated(darea)) deallocate(darea) + if (allocated(ratjav)) deallocate(ratjav) + if (allocated(ratjbv)) deallocate(ratjbv) + if (allocated(ratjplv)) deallocate(ratjplv) + end subroutine dealloc_pec + +end module pec diff --git a/src/polarization.f90 b/src/polarization.f90 new file mode 100644 index 0000000..75c0390 --- /dev/null +++ b/src/polarization.f90 @@ -0,0 +1,152 @@ +module polarization + interface stokes + module procedure stokes_ce,stokes_ell + end interface + +contains + subroutine stokes_ce(ext,eyt,qq,uu,vv) + use const_and_precisions, only : wp_,two + implicit none +! arguments + complex(wp_), intent(in) :: ext,eyt + real(wp_), intent(out) :: qq,uu,vv + + qq = abs(ext)**2 - abs(eyt)**2 + uu = two* dble(ext*dconjg(eyt)) + vv = two*dimag(ext*dconjg(eyt)) + end subroutine stokes_ce + + + subroutine stokes_ell(chi,psi,qq,uu,vv) + use const_and_precisions, only : wp_,two + implicit none +! arguments + real(wp_), intent(in) :: chi,psi + real(wp_), intent(out) :: qq,uu,vv + + qq=cos(two*chi)*cos(two*psi) + uu=cos(two*chi)*sin(two*psi) + vv=sin(two*chi) + end subroutine stokes_ell + + + subroutine polellipse(qq,uu,vv,psi,chi) + use const_and_precisions, only : wp_,half + implicit none +! arguments + real(wp_), intent(in) :: qq,uu,vv + real(wp_), intent(out) :: psi,chi +! real(wp_) :: ll,aa,bb,ell + +! ll = sqrt(qq**2 + uu**2) +! aa = sqrt(half*(1 + ll)) +! bb = sqrt(half*(1 - ll)) +! ell = bb/aa + psi = half*atan2(uu,qq) + chi = half*asin(vv) + end subroutine polellipse + + subroutine pol_limit(anv,bv,bres,sox,ext,eyt) !,gam) + use const_and_precisions, only : wp_,ui=>im,pi,zero,one + implicit none +! arguments + real(wp_), dimension(3), intent(in) :: anv,bv + real(wp_), intent(in) :: bres,sox + complex(wp_), intent(out) :: ext,eyt +! real(wp_), optional, intent(out) :: gam +! local variables + real(wp_), dimension(3) :: bnv + real(wp_) :: anx,any,anz,an2,an,anpl2,anpl,anpr,anxy, & + btot,yg,den,dnl,del0,ff,ff2,sngam,csgam +! + btot = sqrt(bv(1)**2+bv(2)**2+bv(3)**2) + bnv = bv/btot + yg = btot/bres + + anx = anv(1) + any = anv(2) + anz = anv(3) + an2 = anx**2 + any**2 + anz**2 + an = sqrt(an2) + anxy = sqrt(anx**2 + any**2) + + anpl = (anv(1)*bnv(1) + anv(2)*bnv(2) + anv(3)*bnv(3)) + anpl2= anpl**2 + anpr = sqrt(an2 - anpl2) + + dnl = one - anpl2 + del0 = sqrt(dnl**2 + 4.0_wp_*anpl2/yg**2) + + sngam = (anz*anpl - an2*bnv(3))/(an*anxy*anpr) + csgam = -(any*bnv(1) - anx*bnv(2))/ (anxy*anpr) + + ff = 0.5_wp_*yg*(dnl - sox*del0) + ff2 = ff**2 + den = ff2 + anpl2 + if (den>zero) then + ext = (ff*csgam - ui*anpl*sngam)/sqrt(den) + eyt = (-ff*sngam - ui*anpl*csgam)/sqrt(den) + else ! only for XM (sox=+1) when N//=0 + ext = -ui*sngam + eyt = -ui*csgam + end if + +! gam = atan2(sngam,csgam)/degree + end subroutine pol_limit + + subroutine polarcold(anpl,anpr,xg,yg,sox,exf,eyif,ezf,elf,etf) + use const_and_precisions, only : wp_,zero,one + implicit none +! arguments + real(wp_), intent(in) :: anpl,anpr,xg,yg,sox + real(wp_), intent(out) :: exf,eyif,ezf,elf,etf +! local variables + real(wp_) :: anpl2,anpr2,an2,yg2,dy2,aa,e3,qq,p + + if(xg <= zero) then + exf = zero + if(sox < zero) then + ezf = one + eyif = zero + else + ezf = zero + eyif = one + end if + elf = zero + etf = one + else + anpl2 = anpl**2 + anpr2 = anpr**2 + an2 = anpl2 + anpr2 + + yg2=yg**2 + aa=1.0_wp_-xg-yg2 + + dy2 = one - yg2 + qq = xg*yg/(an2*dy2 - aa) + + if (anpl == zero) then + if(sox < zero) then + exf = zero + eyif = zero + ezf = one + else + qq = -aa/(xg*yg) + exf = one/sqrt(one + qq**2) + eyif = qq*exf + ezf = zero + end if + else + e3 = one - xg + p = (anpr2 - e3)/(anpl*anpr) ! undef for anpr==0 + exf = p*ezf + eyif = qq*exf + ezf = one/sqrt(one + p**2*(one + qq**2)) + end if + + elf = (anpl*ezf + anpr*exf)/sqrt(an2) + etf = sqrt(one - elf**2) + end if + end subroutine polarcold + +end module polarization