main subroutine restructured; single indexing of rays; all subroutines updated to Fortran 90 (and all commons removed); two passes feature temporarily disabled.
This commit is contained in:
parent
198feb7a1f
commit
c36ffbc6b6
26
Makefile
26
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,12 +62,6 @@ 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 $<
|
||||
|
||||
.PHONY: clean install
|
||||
|
||||
251
src/beamdata.f90
251
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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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_
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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_)
|
||||
|
||||
4424
src/gray-externals.f
4424
src/gray-externals.f
File diff suppressed because it is too large
Load Diff
895
src/gray-externals.f90
Normal file
895
src/gray-externals.f90
Normal file
@ -0,0 +1,895 @@
|
||||
! program gray
|
||||
! use gray_params, only : ipass,igrad
|
||||
! implicit none
|
||||
!! local variables
|
||||
! real(wp_) :: p0mw1
|
||||
!! common/external functions/variables
|
||||
! integer :: ierr,index_rt
|
||||
! real(wp_) :: sox,p0mw,powrfl,taumn,taumx,pabstot,currtot,
|
||||
!!
|
||||
! common/ierr/ierr
|
||||
! common/mode/sox
|
||||
! common/p0/p0mw
|
||||
! common/powrfl/powrfl
|
||||
! common/index_rt/index_rt
|
||||
! common/taumnx/taumn,taumx,pabstot,currtot
|
||||
!!
|
||||
! if (ipass.gt.1) then
|
||||
!! second pass into plasma
|
||||
! p0mw1=p0mw
|
||||
! igrad=0
|
||||
!!
|
||||
! index_rt=2
|
||||
! p0mw=p0mw1*powrfl
|
||||
! call prfile
|
||||
! call vectinit2
|
||||
! call paraminit
|
||||
! call ic_rt2
|
||||
! call gray_integration
|
||||
! call after_gray_integration
|
||||
! pabstott=pabstott+pabstot
|
||||
! currtott=currtott+currtot
|
||||
!!
|
||||
! index_rt=3
|
||||
! sox=-sox
|
||||
! p0mw=p0mw1*(1.0_wp_-powrfl)
|
||||
! call prfile
|
||||
! call vectinit2
|
||||
! call paraminit
|
||||
! call ic_rt2
|
||||
! call gray_integration
|
||||
! call after_gray_integration
|
||||
! pabstott=pabstott+pabstot
|
||||
! currtott=currtott+currtot
|
||||
! end if
|
||||
!!
|
||||
! end program gray
|
||||
!
|
||||
!
|
||||
!
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
! 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
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!
|
||||
!! single pass is stopped when all the rays have crossed the plasma
|
||||
!! or complete absorption has occurred
|
||||
!! same for successive passes of multi-pass simulations (here exit
|
||||
!! from vessel is detected too
|
||||
!! first pass in multi-pass simulation is stopped when at least one
|
||||
!! ray has reflected and all rays are directed away from
|
||||
!! reflection point, or when no reflection has occurred and
|
||||
!! 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
|
||||
!! 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
|
||||
!! 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
|
||||
!! 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)
|
||||
!! look for first ray hitting the plasma, starting from the central
|
||||
!! 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
|
||||
!! look for first ray which completed the first pass in the plasma
|
||||
! if (iop(j1,k1).gt.1) then
|
||||
!! 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
|
||||
!! if no ray completed a first pass in the plasma, use central ray
|
||||
!! 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
|
||||
!
|
||||
!
|
||||
!
|
||||
! subroutine ic_rt(x00,y00,z00,anx0c,any0c,anz0c,ak0,xgcn,bres,
|
||||
! . wcsi,weta,rcicsi,rcieta,phiw,phir,sox,psipol0,chipol0)
|
||||
!! ray tracing initial conditions igrad=0
|
||||
!!
|
||||
! 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
|
||||
!! 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
|
||||
!! local constants
|
||||
! integer, parameter :: ndim=6,ndimm=3
|
||||
!! 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
|
||||
!! 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
|
||||
!
|
||||
!!
|
||||
! 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
|
||||
!!
|
||||
! 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
|
||||
!!
|
||||
! phiwrad=phiw*cvdr
|
||||
! csphiw=cos(phiwrad)
|
||||
! snphiw=sin(phiwrad)
|
||||
!!
|
||||
! dr=1.0_wp_
|
||||
! if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
|
||||
! da=2.0_wp_*pi/dble(nrayth)
|
||||
! z0t=0.0_wp_
|
||||
!!
|
||||
! 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
|
||||
!!
|
||||
!! csiw=u*dcsiw
|
||||
!! etaw=u*detaw
|
||||
!! csir=csiw
|
||||
!! etar=etaw
|
||||
!!
|
||||
! 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
|
||||
!!
|
||||
! ppcsi=u*dr*cos(alfak)*rcicsi
|
||||
! ppeta=u*dr*sin(alfak)*rcieta
|
||||
!!
|
||||
! anzt=1.0_wp_/sqrt(1.0_wp_+ppcsi**2+ppeta**2)
|
||||
! ancsi=ppcsi*anzt
|
||||
! aneta=ppeta*anzt
|
||||
!!
|
||||
! anxt=ancsi*csphiw-aneta*snphiw
|
||||
! anyt=ancsi*snphiw+aneta*csphiw
|
||||
!!
|
||||
! anx= anxt*csps+snps*(anyt*csth+anzt*snth)
|
||||
! any=-anxt*snps+csps*(anyt*csth+anzt*snth)
|
||||
! anz= anzt*csth-anyt*snth
|
||||
!!
|
||||
! an20=1.0_wp_
|
||||
! an0=sqrt(an20)
|
||||
! anx0=anx
|
||||
! any0=any
|
||||
! anz0=anz
|
||||
!!
|
||||
! xc0(1,j,k)=x0
|
||||
! xc0(2,j,k)=y0
|
||||
! xc0(3,j,k)=z0
|
||||
!!
|
||||
! 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
|
||||
!!
|
||||
! 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_
|
||||
!!
|
||||
! 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
|
||||
!! 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
|
||||
!!
|
||||
! 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_
|
||||
!!
|
||||
! dd=anx0**2+any0**2+anz0**2-an20
|
||||
! vgradi=0.0_wp_
|
||||
! ddi=2.0_wp_*vgradi
|
||||
!!
|
||||
! 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
|
||||
!!
|
||||
! if(nrayr.gt.1) then
|
||||
! iproj=0
|
||||
! nfilp=8
|
||||
! call projxyzt(iproj,nfilp)
|
||||
! end if
|
||||
!!
|
||||
! return
|
||||
!99 format(24(1x,e16.8e3))
|
||||
!111 format(3i5,20(1x,e16.8e3))
|
||||
! end
|
||||
|
||||
|
||||
|
||||
subroutine prfile
|
||||
implicit none
|
||||
write(4,*)' #sst R z phi psi rhot ne Te Btot '// &
|
||||
'Nperp Npl ki alpha tau Pt dIds nh iohkw index_rt ddr'
|
||||
write(8,*) ' #istep j k xt yt zt rt psin'
|
||||
write(9,*) ' #istep j k xt yt zt rt psin'
|
||||
write(17,*) ' #sst Dr_Nr1 Di_Nr1'
|
||||
write(33,*) ' #i jk 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"
|
||||
end subroutine prfile
|
||||
|
||||
|
||||
|
||||
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
|
||||
! local constants
|
||||
real(wp_), parameter :: eps=1.e-4_wp_
|
||||
! local variables
|
||||
integer :: i
|
||||
real(wp_) :: psin,rhop,rhot,ajphi,te,qq
|
||||
real(wp_) :: dens,ddens
|
||||
|
||||
write(55,*) ' #psi rhot ne Te q Jphi'
|
||||
do i=1,nq
|
||||
psin=psinr(i)
|
||||
rhop=sqrt(psin)
|
||||
|
||||
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
|
||||
end subroutine print_prof
|
||||
|
||||
subroutine print_prof_an
|
||||
use const_and_precisions, only : wp_
|
||||
use coreprofiles, only : density, temp
|
||||
use equilibrium, only : frhotor
|
||||
implicit none
|
||||
! local constants
|
||||
integer, parameter :: nst=51
|
||||
! local variables
|
||||
integer :: i
|
||||
real(wp_) :: psin,rhop,rhot,te
|
||||
real(wp_) :: dens,ddens
|
||||
|
||||
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,"(12(1x,e12.5))") psin,rhot,dens,te
|
||||
end do
|
||||
end subroutine print_prof_an
|
||||
|
||||
|
||||
|
||||
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,contours_psi
|
||||
use utils, only : locate, intlin
|
||||
implicit none
|
||||
! arguments
|
||||
integer, intent(in) :: nq
|
||||
real(wp_), dimension(nq), intent(in) :: psinq,qpsi
|
||||
real(wp_) :: qval
|
||||
! local variables
|
||||
integer :: ncnt,i1,ipr
|
||||
real(wp_) :: rup,zup,rlw,zlw,rhot,psival
|
||||
real(wp_), dimension(npoints) :: rcn,zcn
|
||||
|
||||
ncnt=(npoints-1)/2
|
||||
|
||||
! locate psi surface for q=qval
|
||||
call locate(abs(qpsi),nq,qval,i1)
|
||||
if (i1>0.and.i1<nq) then
|
||||
call intlin(abs(qpsi(i1)),psinq(i1),abs(qpsi(i1+1)),psinq(i1+1),qval,psival)
|
||||
rup=rmaxis
|
||||
rlw=rmaxis
|
||||
zup=(zbsup+zmaxis)/2.0_wp_
|
||||
zlw=(zmaxis+zbinf)/2.0_wp_
|
||||
ipr=1
|
||||
call contours_psi(psival,rup,zup,rlw,zlw,rcn,zcn,ipr)
|
||||
rhot=frhotor(sqrt(psival))
|
||||
print'(4(a,f8.5))','q = ',qval, ' psi = ',psival, &
|
||||
' rhop = ',sqrt(psival),' rhot = ',rhot
|
||||
end if
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine bfield_res(rv,zv,nr,nz,bres)
|
||||
use const_and_precisions, only : wp_
|
||||
use equilibrium, only : bfield
|
||||
implicit none
|
||||
! arguments
|
||||
integer, intent(in) :: nr, nz
|
||||
real(wp_), intent(in) :: rv(nr), zv(nz), bres
|
||||
! local constants
|
||||
integer, parameter :: icmx=2002
|
||||
! local variables
|
||||
integer :: j,k,n,nconts,inc,nctot
|
||||
integer, dimension(10) :: ncpts
|
||||
real(wp_) :: btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
|
||||
real(wp_), dimension(icmx) :: rrcb,zzcb
|
||||
real(wp_), dimension(nr,nz) :: btotal
|
||||
|
||||
! Btotal on psi grid
|
||||
btmx=-1.0e30_wp_
|
||||
btmn=1.0e30_wp_
|
||||
do j=1,nr
|
||||
rrj=rv(j)
|
||||
do k=1,nz
|
||||
zzk=zv(k)
|
||||
call bfield(rrj,zzk,bbphi,bbr,bbz)
|
||||
btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2)
|
||||
if(btotal(j,k).ge.btmx) btmx=btotal(j,k)
|
||||
if(btotal(j,k).le.btmn) btmn=btotal(j,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! compute Btot=Bres/n with n=1,5
|
||||
write(70,*)'#i Btot R z'
|
||||
do n=1,5
|
||||
bbb=bres/dble(n)
|
||||
if (bbb.ge.btmn.and.bbb.le.btmx) then
|
||||
call cniteq(rv,zv,btotal,nr,nz,bbb,nconts,ncpts,nctot,rrcb,zzcb)
|
||||
do inc=1,nctot
|
||||
write(70,'(i6,12(1x,e12.5))') inc,bbb,rrcb(inc),zzcb(inc)
|
||||
end do
|
||||
end if
|
||||
write(70,*)
|
||||
end do
|
||||
end subroutine bfield_res
|
||||
|
||||
subroutine bres_anal(bres)
|
||||
use const_and_precisions, only : wp_,pi
|
||||
use equilibrium, only : aminor,rmaxis,zmaxis
|
||||
implicit none
|
||||
! arguments
|
||||
real(wp_) :: bres
|
||||
! local variables
|
||||
integer :: i
|
||||
integer, parameter :: ngrid=51
|
||||
real(wp_) :: dxgrid
|
||||
real(wp_), dimension(ngrid) :: rv,zv
|
||||
|
||||
dxgrid=2.0_wp_*aminor/dble(ngrid-1)
|
||||
do i=1,ngrid
|
||||
rv(i) = rmaxis - aminor + dxgrid*(i-1)
|
||||
zv(i) = zmaxis - aminor + dxgrid*(i-1)
|
||||
end do
|
||||
call bfield_res(rv,zv,ngrid,ngrid,bres)
|
||||
end subroutine bres_anal
|
||||
|
||||
|
||||
|
||||
subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
|
||||
use const_and_precisions, only : wp_
|
||||
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
|
||||
! (based on an older code)
|
||||
use const_and_precisions, only : wp_
|
||||
implicit none
|
||||
! local constants
|
||||
integer, parameter :: icmx=2002
|
||||
! arguments
|
||||
integer :: nr,nz
|
||||
real(wp_), dimension(nr) :: rqgrid
|
||||
real(wp_), dimension(nz) :: zqgrid
|
||||
real(wp_), dimension(nr,nz) :: matr2dgrid
|
||||
integer :: ncon,icount
|
||||
integer, dimension(10) :: npts
|
||||
real(wp_) :: h
|
||||
real(wp_), dimension(icmx) :: rcon,zcon
|
||||
! local variables
|
||||
integer :: i,j,k,l,ico,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb
|
||||
integer :: jabs,jnb,kx,ikx,itm,inext,in
|
||||
integer, dimension(3,2) :: ja
|
||||
integer, dimension(1000) :: lx
|
||||
real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y
|
||||
real(wp_), dimension(nr*nz) :: a
|
||||
logical :: flag1, flag2
|
||||
|
||||
px=0.5_wp_
|
||||
|
||||
a=reshape(matr2dgrid,(/nr*nz/))
|
||||
|
||||
do ico=1,icmx
|
||||
rcon(ico)=0.0_wp_
|
||||
zcon(ico)=0.0_wp_
|
||||
enddo
|
||||
|
||||
nrqmax=nr
|
||||
nr=nr
|
||||
nz=nz
|
||||
drgrd=rqgrid(2)-rqgrid(1)
|
||||
dzgrd=zqgrid(2)-zqgrid(1)
|
||||
|
||||
ncon = 0
|
||||
|
||||
npts = 0
|
||||
|
||||
iclast = 0
|
||||
icount = 0
|
||||
mpl=0
|
||||
ix=0
|
||||
mxr = nrqmax * (nz - 1)
|
||||
n1 = nr - 1
|
||||
|
||||
do jx=2,n1
|
||||
do jm=jx,mxr,nrqmax
|
||||
j = jm + nrqmax
|
||||
ah=a(j)-h
|
||||
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
|
||||
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)<tiny(walln)) then
|
||||
! wall never hit
|
||||
dstvac=0.0_wp_
|
||||
xvend=xv0
|
||||
ivac=-1
|
||||
return
|
||||
end if
|
||||
! search second wall interface (inside-outside)
|
||||
st=smax
|
||||
xvend=xv0+st*anv0
|
||||
call inters_linewall(xvend/1.0e2_wp_,anv0,rlim(1:nlim), &
|
||||
zlim(1:nlim),nlim,smax,walln)
|
||||
smax=smax*1.0e2_wp_+st
|
||||
end if
|
||||
i=0
|
||||
do
|
||||
st=i*dst
|
||||
xvend=xv0+st*anv0
|
||||
rrm=1.0e-2_wp_*sqrt(xvend(1)**2+xvend(2)**2)
|
||||
zzm=1.0e-2_wp_*xvend(3)
|
||||
plfound=inside_plasma(rrm,zzm)
|
||||
if (st.ge.smax.or.plfound) exit
|
||||
i=i+1
|
||||
end do
|
||||
|
||||
if (plfound) then
|
||||
ivac=1
|
||||
dstvac=st
|
||||
else
|
||||
ivac=2
|
||||
dstvac=smax
|
||||
xvend=xv0+smax*anv0
|
||||
end if
|
||||
end subroutine vacuum_rt
|
||||
1434
src/graycore.f90
1434
src/graycore.f90
File diff suppressed because it is too large
Load Diff
@ -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
|
||||
|
||||
@ -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)
|
||||
! ========================================================================
|
||||
|
||||
|
||||
385
src/pec.f90
Normal file
385
src/pec.f90
Normal file
@ -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
|
||||
152
src/polarization.f90
Normal file
152
src/polarization.f90
Normal file
@ -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
|
||||
Loading…
Reference in New Issue
Block a user