istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

* Added case ibeam=2

Browse Source 
* Fixed intent in dierckx/bispev
* New use of bessel_jn function in eccd/fpp
* Some cleaning in cniteq
* Fixed wrong call to gradi_upd in case of single-ray runs
* Fixed call to pec_init for analytical equilibria
This commit is contained in:
Lorenzo Figini 2015-11-18 16:28:15 +00:00
parent 4226416c4a
commit c4a409f8c5
10 changed files with 691 additions and 144 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile
View File

@ -35,7 +35,7 @@ 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
beams.o: const_and_precisions.o simplespline.o utils.o
beams.o: const_and_precisions.o dierckx.o reflections.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 \

20
src/beamdata.f90
View File

@ -9,13 +9,13 @@ module beamdata
contains
subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri, &
psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv)
psjki,tauv,alphav,ppabs,dids,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
gri,psjki,tauv,alphav,ppabs,dids,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
@ -41,7 +41,7 @@ contains
nstep=rtrparam%nstep
call alloc_beam(ywork,ypwork,xc,du1,gri,ggri, &
psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv)
psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
end subroutine init_rtr
function rayi2jk(i) result(jk)
@ -101,31 +101,31 @@ contains
end function rayjk2i
subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri, &
psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv)
psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
implicit none
real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, &
gri,psjki,tauv,alphav,ppabs,didst,ccci
gri,psjki,tauv,alphav,ppabs,dids,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(ywork,ypwork,xc,du1,gri,ggri, &
psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv)
psjki,tauv,alphav,ppabs,dids,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), &
ppabs(nray,nstep), dids(nray,nstep), ccci(nray,nstep), &
p0jk(nray), ext(nray), eyt(nray), iiv(nray))
end subroutine alloc_beam
subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri, &
psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv)
psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
implicit none
real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, &
gri,psjki,tauv,alphav,ppabs,didst,ccci
gri,psjki,tauv,alphav,ppabs,dids,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
@ -141,7 +141,7 @@ contains
if (allocated(tauv)) deallocate(tauv)
if (allocated(alphav)) deallocate(alphav)
if (allocated(ppabs)) deallocate(ppabs)
if (allocated(didst)) deallocate(didst)
if (allocated(dids)) deallocate(dids)
if (allocated(ccci)) deallocate(ccci)
if (allocated(p0jk)) deallocate(p0jk)
if (allocated(ext)) deallocate(ext)

582
src/beams.f90
View File

@ -4,14 +4,13 @@ module beams
contains
subroutine read_beam0(file_beam,alpha0,beta0,fghz,x00,y00,z00, &
subroutine read_beam0(file_beam,fghz,x00,y00,z00, &
wcsi,weta,rcicsi,rcieta,phiw,phir,unit)
use const_and_precisions, only : pi,vc=>ccgs_
use utils, only : get_free_unit
implicit none
! arguments
character(len=*), intent(in) :: file_beam
real(wp_), intent(in) :: alpha0,beta0
real(wp_), intent(out) :: fGHz,x00,y00,z00
real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw
integer, intent(in), optional :: unit
@ -64,8 +63,7 @@ contains
integer, intent(in), optional :: unit
! local variables
integer :: u,ierr,iopt,ier,nisteer,i,k,ii
real(wp_) :: steer,alphast,betast,x00mm,y00mm,z00mm
real(wp_) :: w1,w2,rci1,rci2,phi1,phi2,dal
real(wp_) :: steer,dal
real(wp_), dimension(:), allocatable :: alphastv,betastv,x00v,y00v, &
z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v, &
cbeta,cx0,cy0,cz0,cwaist1,cwaist2, &
@ -93,29 +91,25 @@ contains
if (ierr/=0) then
close(u)
deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v, &
phi1v,phi2v,x00v,y00v,z00v,cbeta, &
phi1v,phi2v,x00v,y00v,z00v,cbeta, &
cx0,cy0,cz0,cwaist1,cwaist2,crci1,crci2,cphi1,cphi2)
write(*,*) 'cannot allocate arrays for beam data'
stop
end if
do i=1,nisteer
read(u,*) steer,alphast,betast,x00mm,y00mm,z00mm, &
w1,w2,rci1,rci2,phi1,phi2
! initial beam data measured in mm -> transformed to cm
x00v(i)=0.1_wp_*x00mm
y00v(i)=0.1_wp_*y00mm
z00v(i)=0.1_wp_*z00mm
alphastv(i)=alphast
betastv(i)=betast
waist1v(i)=0.1_wp_*w1
rci1v(i)=1.0e1_wp_*rci1
waist2v(i)=0.1_wp_*w2
rci2v(i)=1.0e1_wp_*rci2
phi1v(i)=phi1
phi2v(i)=phi2
read(u,*) steer,alphastv(i),betastv(i),x00v(i),y00v(i),z00v(i), &
waist1v(i),waist2v(i),rci1v(i),rci2v(i),phi1v(i),phi2v(i)
end do
close(u)
! initial beam data measured in mm -> transformed to cm
x00v = 0.1_wp_*x00v
y00v = 0.1_wp_*y00v
z00v = 0.1_wp_*z00v
waist1v = 0.1_wp_*waist1v
waist2v = 0.1_wp_*waist2v
rci1v = 10._wp_*rci1v
rci2v = 10._wp_*rci2v
iopt=0
call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier)
@ -164,6 +158,556 @@ contains
end subroutine read_beam1
subroutine read_beam2(file_beam,beamid,alpha0,beta0,fghz,iox,x00,y00,z00, &
wcsi,weta,rcicsi,rcieta,phiw,phir,unit)
use utils, only : get_free_unit, intlin, locate
use reflections, only : inside
use dierckx, only : curfit, splev, surfit, bispev
implicit none
character(len=*), intent(in) :: file_beam
integer, intent(in) :: beamid
real(wp_), intent(inout) :: alpha0,beta0
real(wp_), intent(out) :: fghz,x00,y00,z00, wcsi,weta,rcicsi,rcieta,phir,phiw
integer, intent(out) :: iox
integer, intent(in), optional :: unit
character(len=20) :: beamname
integer :: u
integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta
integer :: iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk
integer :: nxwaist1, nywaist1, nxwaist2, nywaist2, nxrci1, nyrci1, nxrci2
integer :: nyrci2, nxphi1, nyphi1, nxphi2, nyphi2, nxx0, nyx0, nxy0, nyy0
integer :: nxz0, nyz0, kx, ky, ii, npolyg, nmax, lwrk2, in
integer :: nxycoord
integer, DIMENSION(:), ALLOCATABLE :: iwrk
real(wp_) :: alphast,betast, waist1, waist2, rci1, rci2, phi1, phi2
real(wp_) :: fp, minx, maxx, miny, maxy, eps, xcoord0, ycoord0
real(wp_), DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv, &
betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord, &
ycoord, wrk, txwaist1, tywaist1, txwaist2, tywaist2, &
txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2, &
txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1, &
cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2, &
xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC, &
ypolygC, xpolygD, ypolygD
real(wp_), DIMENSION(4) :: xvert, yvert
real(wp_), dimension(1) :: fi
integer, parameter :: kspl=1
real(wp_), parameter :: sspl=0.01_wp_
if (present(unit)) then
u=unit
else
u = get_free_unit()
end if
open(unit=u,file=file_beam,status='OLD',action='READ')
!=======================================================================================
! # of beams
read(u,*) nbeam
!
! unused beams' data
jumprow=0
! c====================================================================================
do i=1,beamid-1
read(u,*) beamname, iox, fghz, nalpha, nbeta
jumprow = jumprow+nalpha*nbeta
end do
! c====================================================================================
!
! beam of interest
read(u,*) beamname, iox, fghz, nalpha, nbeta
!
! c====================================================================================
! unused beams' data grids
do i=1,(nbeam - beamid)
read(u,*) beamname
end do
do i=1,jumprow
read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2
end do
! c====================================================================================
!
! # of elements in beam data grid
nisteer = nalpha*nbeta
!
allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), &
waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), &
phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer), &
xcoord(nisteer),ycoord(nisteer))
!
! c====================================================================================
! beam data grid reading
do i=1,nisteer
read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2
!
! initial beam data (x00, y00, z00) are measured in mm -> transformed to cm
x00v(i)=0.1d0*x00
y00v(i)=0.1d0*y00
z00v(i)=0.1d0*z00
alphastv(i)=alphast
betastv(i)=betast
waist1v(i)=0.1d0*waist1
rci1v(i)=1.0d1*rci1
waist2v(i)=0.1d0*waist2
rci2v(i)=1.0d1*rci2
phi1v(i)=phi1
phi2v(i)=phi2
end do
close(u)
! c====================================================================================
!
! fdeg = 0 alpha, beta free variables
! 1 alpha free variable
! 2 beta free variable
! 3 no free variables
fdeg = 2*(1/nalpha) + 1/nbeta
!#######################################################################################
!
! no free variables
if(fdeg.eq.3) then
alpha0=alphastv(1)
beta0=betastv(1)
x00=x00v(1)
y00=y00v(1)
z00=z00v(1)
wcsi=waist1v(1)
weta=waist2v(1)
rcicsi=rci1v(1)
rcieta=rci2v(1)
phiw=phi1v(1)
phir=phi2v(1)
return
end if
!#######################################################################################
!
!
!#######################################################################################
if(fdeg.eq.2) then
! beta = independent variable
! alpha = dependent variable
xcoord = betastv
ycoord = alphastv
xcoord0 = beta0
ycoord0 = alpha0
kx=min(nbeta-1,kspl)
! c====================================================================================
else
! c====================================================================================
! alpha = independent variable
! beta = dependent/independent (fdeg = 1/0)
xcoord = alphastv
ycoord = betastv
xcoord0 = alpha0
ycoord0 = beta0
nxcoord = nalpha
nycoord = nbeta
kx=min(nalpha-1,kspl)
ky=min(nbeta-1,kspl)
end if
!#######################################################################################
!
iopt = 0
incheck = 0
!
!#######################################################################################
if(fdeg.ne.0) then
nxest = kx + nxcoord + 1
lwrk = (nxcoord*(kx+1)+nxest*(7+3*kx))
kwrk = nxest
allocate(cycoord(nxest), txycoord(nxest), cwaist1(nxest), &
txwaist1(nxest), cwaist2(nxest), txwaist2(nxest), &
crci1(nxest), txrci1(nxest), crci2(nxest), txrci2(nxest), &
cphi1(nxest), txphi1(nxest), cphi2(nxest), txphi2(nxest), &
cx0(nxest), txx0(nxest), cy0(nxest), txy0(nxest), &
cz0(nxest), txz0(nxest), w(nxcoord), wrk(lwrk), iwrk(kwrk))
!
w = 1.d0
!
! 2D interpolation
call curfit(iopt,nxcoord,xcoord,ycoord,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxycoord, &
txycoord,cycoord,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,waist1v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist1, &
txwaist1,cwaist1,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,waist2v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist2, &
txwaist2,cwaist2,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,rci1v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci1, &
txrci1,crci1,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,rci2v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci2, &
txrci2,crci2,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,phi1v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi1, &
txphi1,cphi1,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,phi2v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi2, &
txphi2,cphi2,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,x00v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxx0, &
txx0,cx0,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,y00v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxy0, &
txy0,cy0,fp,wrk,lwrk,iwrk,ier)
!
call curfit(iopt,nxcoord,xcoord,z00v,w, &
xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxz0, &
txz0,cz0,fp,wrk,lwrk,iwrk,ier)
!
! check if xcoord0 is out of table range
! incheck = 1 inside / 0 outside
if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer)) incheck=1
! c====================================================================================
else
! c====================================================================================
npolyg = 2*(nxcoord+nycoord-2)
minx = minval(xcoord)
maxx = maxval(xcoord)
miny = minval(ycoord)
maxy = maxval(ycoord)
nxest = kx + 1 + int(sqrt(nisteer/2.))
nyest = ky + 1 + int(sqrt(nisteer/2.))
nmax = max(nxest,nyest)
eps = 10.**(-8)
lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19
lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2
kwrk = nisteer+(nmax-3)*(nmax-3)
allocate(cwaist1(nxest*nyest), txwaist1(nmax), tywaist1(nmax), &
cwaist2(nxest*nyest), txwaist2(nmax), tywaist2(nmax), &
crci1(nxest*nyest), txrci1(nmax), tyrci1(nmax), &
crci2(nxest*nyest), txrci2(nmax), tyrci2(nmax), &
cphi1(nxest*nyest), txphi1(nmax), typhi1(nmax), &
cphi2(nxest*nyest), txphi2(nmax), typhi2(nmax), &
cx0(nxest*nyest), txx0(nmax), tyx0(nmax), &
cy0(nxest*nyest), txy0(nmax), tyy0(nmax), &
cz0(nxest*nyest), txz0(nmax), tyz0(nmax), &
wrk(lwrk), wrk2(lwrk2), iwrk(kwrk), &
xpolyg(npolyg), ypolyg(npolyg), w(nisteer))
!
w = 1.d0
!
! 3D interpolation
call surfit(iopt,nisteer,xcoord,ycoord,waist1v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxwaist1,txwaist1,nywaist1,tywaist1,cwaist1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,waist2v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxwaist2,txwaist2,nywaist2,tywaist2,cwaist2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,rci1v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxrci1,txrci1,nyrci1,tyrci1,crci1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,rci2v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxrci2,txrci2,nyrci2,tyrci2,crci2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,phi1v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxphi1,txphi1,nyphi1,typhi1,cphi1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,phi2v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxphi2,txphi2,nyphi2,typhi2,cphi2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,x00v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxx0,txx0,nyx0,tyx0,cx0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,y00v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxy0,txy0,nyy0,tyy0,cy0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
!
call surfit(iopt,nisteer,xcoord,ycoord,z00v,w, &
minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, &
nxz0,txz0,nyz0,tyz0,cz0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier)
! data range polygon
xpolyg(1:nxcoord) = xcoord(1:nxcoord)
ypolyg(1:nxcoord) = ycoord(1:nxcoord)
!
! c====================================================================================
do i=1,nycoord-2
xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord)
xpolyg(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1)
ypolyg(nxcoord+i) = ycoord((i+1)*nxcoord)
ypolyg(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1)
end do
! c====================================================================================
do i=1,nxcoord
xpolyg(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1)
ypolyg(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1)
end do
! c====================================================================================
!
! check if (xcoord0, ycoord0) is out of table range
! incheck = 1 inside / 0 outside
if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0)) incheck = 1
end if
deallocate(wrk,iwrk)
!#######################################################################################
!
!
!#######################################################################################
if(fdeg.ne.0) then
! c====================================================================================
if(incheck.eq.1) then
call splev(txycoord,nxycoord,cycoord,kx,(/xcoord0/),fi,1,ier)
ycoord0=fi(1)
call splev(txwaist1,nxwaist1,cwaist1,kx,(/xcoord0/),fi,1,ier)
wcsi=fi(1)
call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier)
weta=fi(1)
call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier)
rcicsi=fi(1)
call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier)
rcieta=fi(1)
call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier)
phiw=fi(1)
call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier)
phir=fi(1)
call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier)
x00=fi(1)
call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier)
y00=fi(1)
call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier)
z00=fi(1)
! c----------------------------------------------------------------------------------
else
! c----------------------------------------------------------------------------------
if(xcoord0.ge.xcoord(nisteer)) ii=nisteer
if(xcoord0.le.xcoord(1)) ii=1
!
xcoord0=xcoord(ii)
ycoord0=ycoord(ii)
x00=x00v(ii)
y00=y00v(ii)
z00=z00v(ii)
wcsi=waist1v(ii)
weta=waist2v(ii)
rcicsi=rci1v(ii)
rcieta=rci2v(ii)
phiw=phi1v(ii)
phir=phi2v(ii)
end if
! c====================================================================================
else
! c====================================================================================
if(incheck.eq.0) then
allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), &
ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), &
xpolygD(nycoord), ypolygD(nycoord))
! coordinates of vertices v1,v2,v3,v4
xvert(1) = xpolyg(1)
xvert(2) = xpolyg(nxcoord)
xvert(3) = xpolyg(nxcoord+nycoord-1)
xvert(4) = xpolyg(2*nxcoord+nycoord-2)
yvert(1) = ypolyg(1)
yvert(2) = ypolyg(nxcoord)
yvert(3) = ypolyg(nxcoord+nycoord-1)
yvert(4) = ypolyg(2*nxcoord+nycoord-2)
! coordinates of side A,B,C,D
xpolygA = xpolyg(1:nxcoord)
ypolygA = ypolyg(1:nxcoord)
xpolygB = xpolyg(nxcoord:nxcoord+nycoord-1)
ypolygB = ypolyg(nxcoord:nxcoord+nycoord-1)
xpolygC = xpolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
ypolygC = ypolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2)
xpolygD(1:nycoord-1) = xpolyg(2*nxcoord+nycoord-2:npolyg)
xpolygD(nycoord) = xpolyg(1)
ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg)
ypolygD(nycoord) = ypolyg(1)
! c----------------------------------------------------------------------------------
! search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid
!
! | |
! (6) (5) (4)
! | |
! _ _ _ v4 _________________v3_ _ _ _
! | C | (1)->(8) outside regions
! | |
! (7) D | | B (3) v1->v4 grid vertices
! | |
! _ _ _ _ |_________________| _ _ _ _ A-D grid sides
! v1 A v2
! | |
! (8) (1) (2)
! | |
!
if(xcoord0.gt.xvert(1).and.xcoord0.lt.xvert(2).and.ycoord0.le.maxval(ypolygA)) then
in=1
else if(ycoord0.gt.yvert(2).and.ycoord0.lt.yvert(3).and.xcoord0.ge.minval(xpolygB)) then
in=3
else if(xcoord0.lt.xvert(3).and.xcoord0.gt.xvert(4).and.ycoord0.ge.minval(ypolygC)) then
in=5
else if(ycoord0.lt.yvert(4).and.ycoord0.gt.yvert(1).and.xcoord0.le.maxval(xpolygD)) then
in=7
else if(xcoord0.ge.xvert(2).and.ycoord0.le.yvert(2)) then
in=2
else if(xcoord0.ge.xvert(3).and.ycoord0.ge.yvert(3)) then
in=4
else if(xcoord0.le.xvert(4).and.ycoord0.ge.yvert(4)) then
in=6
else if(xcoord0.le.xvert(1).and.ycoord0.le.yvert(1)) then
in=8
endif
! c----------------------------------------------------------------------------------
!
! c----------------------------------------------------------------------------------
! (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border
! depending on the region
! 1: xcoord0 unchanged, ycoord0 moved on side A
! 3: xcoord0 moved on side B, ycoord0 unchanged
! 5: xcoord0 unchanged, ycoord0 moved on side C
! 7: xcoord0 moved on side D, ycoord0 unchanged
! 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates
! in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi in
! new (xcoord0,ycoord0)
! in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the
! (xcoord0,ycoord0) vertex are used
alpha0 = xcoord0
beta0 = ycoord0
SELECT CASE (in)
CASE (1)
write(*,*) ' beta0 outside table range !!!'
! locate position of xcoord0 with respect to x coordinates of side A
call locate(xpolygA,nxcoord,xcoord0,ii)
! find corresponding y value on side A for xcoord position
call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0)
incheck = 1
CASE (2)
write(*,*) ' alpha0 and beta0 outside table range !!!'
! xcoord0, ycoord0 set
xcoord0 = xvert(2)
ycoord0 = yvert(2)
ii = nxcoord !indice per assegnare valori waist, rci, phi
CASE (3)
write(*,*) ' alpha0 outside table range !!!'
call locate(ypolygB,nycoord,ycoord0,ii)
call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0)
incheck = 1
CASE (4)
write(*,*) ' alpha0 and beta0 outside table range !!!'
xcoord0 = xvert(3)
ycoord0 = yvert(3)
ii = nxcoord+nycoord-1
CASE (5)
write(*,*) ' beta0 outside table range !!!'
call locate(xpolygC,nxcoord,xcoord0,ii)
call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0)
incheck = 1
CASE (6)
write(*,*) ' alpha0 and beta0 outside table range !!!'
xcoord0 = xvert(4)
ycoord0 = yvert(4)
ii = 2*nxcoord+nycoord-2
CASE (7)
write(*,*) ' alpha0 outside table range !!!'
call locate(ypolygD,nycoord,ycoord0,ii)
call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0)
incheck = 1
CASE (8)
write(*,*) ' alpha0 and beta0 outside table range !!!!'
xcoord0 = xvert(1)
ycoord0 = yvert(1)
ii = 1
END SELECT
! c----------------------------------------------------------------------------------
!
deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, xpolygD, ypolygD)
end if
! c====================================================================================
!
! c====================================================================================
if(incheck.eq.1) then
lwrk = 2*(kx+ky+2)
kwrk = 4
allocate(wrk(lwrk),iwrk(kwrk))
call bispev(txwaist1,nxwaist1,tywaist1,nywaist1,cwaist1, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
wcsi=fi(1)
call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
weta=fi(1)
call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
rcicsi=fi(1)
call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
rcieta=fi(1)
call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
phiw=fi(1)
call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
phir=fi(1)
call bispev(txx0,nxx0,tyx0,nyx0,cx0, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
x00=fi(1)
call bispev(txy0,nxy0,tyy0,nyy0,cy0, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
y00=fi(1)
call bispev(txz0,nxz0,tyz0,nyz0,cz0, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
z00=fi(1)
deallocate(wrk,iwrk)
! c----------------------------------------------------------------------------------
else
! c----------------------------------------------------------------------------------
x00=x00v(ii)
y00=y00v(ii)
z00=z00v(ii)
wcsi=waist1v(ii)
weta=waist2v(ii)
rcicsi=rci1v(ii)
rcieta=rci2v(ii)
phiw=phi1v(ii)
phir=phi2v(ii)
end if
! c====================================================================================
end if
!#######################################################################################
!
if(fdeg.ne.0) then
deallocate(cycoord, txycoord, cwaist1, txwaist1, cwaist2, &
txwaist2, crci1, txrci1, crci2, txrci2, cphi1, txphi1, &
cphi2, txphi2, cx0, txx0, cy0, txy0, cz0, txz0, w)
else
deallocate(cwaist1, txwaist1, tywaist1, cwaist2, txwaist2, tywaist2, &
crci1, txrci1, tyrci1, crci2, txrci2, tyrci2, &
cphi1, txphi1, typhi1, cphi2, txphi2, typhi2, &
cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0, &
wrk2, xpolyg, ypolyg, w)
end if
!
!#######################################################################################
! set correct values for alpha, beta
if(fdeg.eq.2) then
alpha0 = ycoord0
beta0 = xcoord0
else
alpha0 = xcoord0
beta0 = ycoord0
end if
!#######################################################################################
deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, &
phi2v,x00v,y00v,z00v,xcoord,ycoord)
!
end subroutine read_beam2
subroutine launchangles2n(alpha,beta,xv,anv)
use const_and_precisions, only : degree
implicit none

4
src/dierckx.f90
View File

@ -80,9 +80,9 @@ contains
integer, intent(in) :: nx, ny, kx, ky, mx, my, lwrk, kwrk
integer, intent(out) :: ier
integer, intent(inout) :: iwrk(kwrk)
real(wp_), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1))
real(wp_), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx), y(my)
real(wp_), intent(out) :: z(mx*my)
real(wp_), intent(inout) :: x(mx), y(my), wrk(lwrk)
real(wp_), intent(inout) :: wrk(lwrk)
! local variables
integer :: i, iw, lwest
! ..

24
src/eccd.f90
View File

@ -275,7 +275,8 @@ contains
end if
resj=0.0_wp_
do i=0,iokhawa
! do i=0,iokhawa
do i=0,1
resji=0.0_wp_
xx1=amu*(anpl*uleft(i)+ygn-1.0_wp_)
xx2=amu*(anpl*uright(i)+ygn-1.0_wp_)
@ -343,9 +344,10 @@ contains
real(wp_) :: upl,fpp
real(wp_), dimension(npar) :: extrapar
! local variables
integer :: ithn,nhn,nm,np
real(wp_) :: yg,anpl,amu,anprre,ygn,upr,upr2,gam,ee,thn2,thn2u,bb,cth, &
ajbnm,ajbnp,ajbn
integer :: ithn,nhn !,nm,np
real(wp_) :: yg,anpl,amu,anprre,ygn,upr,upr2,gam,ee,thn2,thn2u,bb,cth !, &
! ajbnm,ajbnp,ajbn
real(wp_), dimension(3) :: ajb
complex(wp_) :: ex,ey,ez,emxy,epxy
yg=extrapar(1)
@ -379,12 +381,14 @@ contains
thn2u=upr2*thn2
else
! Full polarization term
nm=nhn-1
np=nhn+1
ajbnm=dbesjn(nm, bb)
ajbnp=dbesjn(np, bb)
ajbn=dbesjn(nhn, bb)
thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0_wp_))**2
! nm=nhn-1
! np=nhn+1
! ajbnm=dbesjn(nm, bb)
! ajbnp=dbesjn(np, bb)
! ajbn=dbesjn(nhn, bb)
! thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0_wp_))**2
ajb=bessel_jn(nhn-1, nhn+1, bb)
thn2u=(abs(ez*ajb(2)*upl+upr*(ajb(3)*epxy+ajb(1)*emxy)/2.0_wp_))**2
end if
end if
end if

48
src/gray-externals.f90
View File

@ -556,6 +556,8 @@
do n=1,5
bbb=bres/dble(n)
if (bbb.ge.btmn.and.bbb.le.btmx) then
nconts=size(ncpts)
nctot=size(rrcb)
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)
@ -593,40 +595,34 @@
! (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
integer, intent(in) :: nr,nz
real(wp_), dimension(nr), intent(in) :: rqgrid
real(wp_), dimension(nz), intent(in) :: zqgrid
real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid
real(wp_), intent(in) :: h
integer, intent(inout) :: ncon, icount
integer, dimension(ncon), intent(out) :: npts
real(wp_), dimension(icount), intent(out) :: rcon,zcon
! local variables
integer :: i,j,k,l,ico,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb
integer :: i,j,k,l,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
integer, dimension(icount/2-1) :: lx
real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y
real(wp_), dimension(nr*nz) :: a
logical :: flag1, flag2
logical :: flag1
px=0.5_wp_
px = 0.5_wp_
a=reshape(matr2dgrid,(/nr*nz/))
a = reshape(matr2dgrid,(/nr*nz/))
do ico=1,icmx
rcon(ico)=0.0_wp_
zcon(ico)=0.0_wp_
enddo
rcon = 0.0_wp_
zcon = 0.0_wp_
nrqmax=nr
nr=nr
nz=nz
drgrd=rqgrid(2)-rqgrid(1)
dzgrd=zqgrid(2)-zqgrid(1)
nrqmax = nr
drgrd = rqgrid(2) - rqgrid(1)
dzgrd = zqgrid(2) - zqgrid(1)
ncon = 0
@ -634,8 +630,8 @@
iclast = 0
icount = 0
mpl=0
ix=0
mpl = 0
ix = 0
mxr = nrqmax * (nz - 1)
n1 = nr - 1

97
src/graycore.f90
View File

@ -12,8 +12,8 @@ contains
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
use dispersion, only : expinit
use gray_params, only : eqparam_type, prfparam_type, outparam_type, &
rtrparam_type, hcdparam_type, set_codepar, igrad, iequil, iprof, ieccd, &
iwarm, ipec, istpr0
rtrparam_type, hcdparam_type, set_codepar, iequil, iprof, ieccd, &
iwarm, ipec, istpr0, igrad
use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff
use beamdata, only : pweight, print_projxyzt, rayi2jk
use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, &
@ -54,8 +54,8 @@ contains
real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre
real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,zeff,dersdst,derdnm,st,st0
real(wp_) :: tau0,alphaabs0,didst0,ccci0
real(wp_) :: tau,pow,dpdst,ddr,ddi,taumn,taumx
real(wp_) :: tau0,alphaabs0,dids0,ccci0
real(wp_) :: tau,pow,ddr,ddi,taumn,taumx
real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava
real(wp_), dimension(3) :: xv,anv0,anv
real(wp_), dimension(:,:), allocatable :: yw,ypw,gri
@ -63,7 +63,7 @@ contains
integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,index_rt=1
logical :: ins_pl, somein, allout
real(wp_), dimension(:,:), allocatable :: psjki,tauv,alphav,ppabs,didst,ccci
real(wp_), dimension(:,:), allocatable :: psjki,tauv,alphav,ppabs,dids,ccci
real(wp_), dimension(:), allocatable :: p0jk
complex(wp_), dimension(:), allocatable :: ext, eyt
integer, dimension(:), allocatable :: iiv
@ -104,7 +104,7 @@ contains
call xgygcoeff(fghz,ak0,bres,xgcn)
call launchangles2n(alpha0,beta0,xv0,anv0)
call init_rtr(rtrp,yw,ypw,xc,du1,gri,ggri, &
psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv)
psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
if(iwarm > 1) call expinit
@ -127,7 +127,7 @@ contains
iox=iox0
sox=-1.0_wp_
if(iox==2) sox=1.0_wp_
call vectinit(psjki,tauv,alphav,ppabs,didst,ccci,iiv)
call vectinit(psjki,tauv,alphav,ppabs,dids,ccci,iiv)
call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri)
psipol=psipol0
@ -151,8 +151,7 @@ contains
end do
! 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)
if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
! test if the beam is completely out of the plsama
allout = .true.
@ -163,8 +162,8 @@ contains
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_
if( abs(anpl) > anplth1) then
if(abs(anpl) <= anplth2) then
ierr=97
! igrad=0
else
@ -175,16 +174,15 @@ contains
ierr=0
end if
tekev=zero
if(i==1) then
tau0=zero
alphaabs0=zero
didst0=zero
dids0=zero
ccci0=zero
else
tau0=tauv(jk,i-1)
alphaabs0=alphav(jk,i-1)
didst0=didst(jk,i-1)
dids0=dids(jk,i-1)
ccci0=ccci(jk,i-1)
end if
zzm = xv(3)*0.01_wp_
@ -196,13 +194,12 @@ contains
if(ierr==0 .and. iwarm>0 .and. ins_pl .and. tau0<=taucr) 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
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)
else
tekev=zero
zeff=zero
alpha=zero
didp=zero
anprim=zero
@ -211,9 +208,12 @@ contains
nhf=0
iokhawa=0
end if
if(nharm>0) iiv(jk)=i
! full storage required only for psjki,ppabs,ccci
! (jk,i) indexing can be removed from tauv,alphav,dids
! adding (jk) indexing to alphaabs0,tau0,dids0,ccci0
psjki(jk,i) = psinv
! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s)
tau=tau0+0.5_wp_*(alpha+alphaabs0)*dersdst*dst
tauv(jk,i)=tau
@ -221,12 +221,11 @@ contains
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)=ccci0+0.5_wp_*(didst0+didst(jk,i))*dst
dids(jk,i)=didp*pow*alpha
ccci(jk,i)=ccci0+0.5_wp_*(dids0+dids(jk,i))*dersdst*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, &
anprim,dens,tekev,alpha,tau,dids(jk,i),nhf,iokhawa, &
index_rt,ddr,ddi)
! print error code
@ -274,7 +273,7 @@ contains
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))
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)
@ -291,7 +290,7 @@ contains
! ======= free memory BEGIN ======
call dealloc_beam(yw,ypw,xc,du1,gri,ggri, &
psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv)
psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
! call unset_eqspl
! call unset_q
! call unset_rhospl
@ -302,13 +301,11 @@ contains
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
subroutine vectinit(psjki,tauv,alphav,ppabs,dids,ccci,iiv)
use const_and_precisions, only : wp_, zero
implicit none
! arguments
real(wp_), dimension(:,:), intent(out) :: psjki,tauv,alphav,ppabs,didst,ccci
real(wp_), dimension(:,:), intent(out) :: psjki,tauv,alphav,ppabs,dids,ccci
integer, dimension(:), intent(out) :: iiv
!! common/external functions/variables
! integer :: jclosest
@ -324,9 +321,9 @@ contains
tauv = zero
alphav = zero
ppabs = zero
didst = zero
dids = zero
ccci = zero
iiv = one
iiv = 1
end subroutine vectinit
@ -337,7 +334,7 @@ contains
! !!!!!! 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 gray_params, only : idst
use beamdata, only : nray,nrayr,nrayth,rwmax
implicit none
! arguments
@ -412,9 +409,9 @@ contains
end if
! w01=sqrt(2.0_wp_/(ak0*ww1))
! z01=-rci1/(rci1**2+ww1**2)
! d01=-rci1/(rci1**2+ww1**2)
! w02=sqrt(2.0_wp_/(ak0*ww2))
! z02=-rci2/(rci2**2+ww2**2)
! d02=-rci2/(rci2**2+ww2**2)
qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2
qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2
@ -625,7 +622,7 @@ contains
subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr)
! Runge-Kutta integrator
use const_and_precisions, only : wp_
use gray_params, only : igrad
! use gray_params, only : igrad
use beamdata, only : h,hh,h6
implicit none
real(wp_), intent(in) :: sox,bres,xgcn
@ -658,7 +655,6 @@ contains
! 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
@ -686,7 +682,7 @@ contains
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
! use gray_params, only : igrad
implicit none
! arguments
real(wp_), dimension(3), intent(in) :: xv,anv
@ -924,7 +920,7 @@ contains
real(wp_), dimension(3), intent(out) :: bv,derxg,deryg
real(wp_), dimension(3,3), intent(out) :: derbv
! local variables
integer :: iv,jv
integer :: 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
@ -1076,7 +1072,7 @@ contains
real(wp_), dimension(3,3), intent(in) :: ddgr,derbv
real(wp_), dimension(6), intent(out) :: dery
! local variables
integer :: iv,jv
integer :: iv
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
@ -1192,8 +1188,7 @@ contains
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 equilibrium, only : sgnbphi
use dispersion, only : harmnumber, warmdisp
use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl
use magsurf_data, only : fluxval
@ -1219,6 +1214,10 @@ contains
anprim=zero
anprre=zero
didp=zero
nhmin=0
nhmax=0
iokhawa=0
ierr=0
if(tekev>zero) then
! absorption computation
@ -1334,7 +1333,7 @@ contains
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)
dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)
use const_and_precisions, only : degree,zero,one
use equilibrium, only : frhotor
use gray_params, only : istpl0
@ -1344,9 +1343,9 @@ contains
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
real(wp_), intent(in) :: dens,tekev,alpha,tau,dids,ddr,ddi
! local variables
real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,dids
real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,didsn
integer :: k
stm=st*1.0e-2_wp_
@ -1366,10 +1365,10 @@ contains
end if
akim=anprim*ak0*1.0e2_wp_
pt=exp(-tau)
dids=didst*1.0e2_wp_/qj
didsn=dids*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), &
btot,anpr,anpl,akim,alpha,tau,pt,didsn,dble(nhf),dble(iokhawa), &
dble(index_rt),ddr
end if
! central ray only end

13
src/magsurf_data.f90
View File

@ -212,10 +212,10 @@ contains
lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, &
kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam
! local variables
integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr,u56,unit
integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr
integer, dimension(kwrk) :: iwrk
real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, &
ratio_cdbtor,ratio_pltor,fc,height,height2,r2iav,currp, &
ratio_cdbtor,ratio_pltor,fc,height,r2iav,currp, &
area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, &
dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, &
shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, &
@ -231,9 +231,6 @@ contains
! common/external functions/variables
real(wp_) :: fpolv,ddpsidrr,ddpsidzz
u56=get_free_unit()
open(file='f56.txt',unit=u56)
npsi=nnintp
ninpr=(npsi-1)/10
npoints = 2*ncnt+1
@ -478,7 +475,7 @@ contains
end do
end do
write(u56,*)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol |I_pl| <J_phi> fc ratJa ratJb'
write(56,*)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol |I_pl| <J_phi> fc ratJa ratJb'
do jp=1,npsi
if(jp.eq.npsi) then
@ -486,13 +483,11 @@ contains
pstab(jp)=1.0_wp_
end if
rhotjp=frhotor(rpstab(jp))
write(u56,99) rpstab(jp),rhotjp,bav(jp),bmxpsi(jp),bmnpsi(jp), &
write(56,99) rpstab(jp),rhotjp,bav(jp),bmxpsi(jp),bmnpsi(jp), &
varea(jp),vvol(jp),vcurrp(jp),vajphiav(jp),ffc(jp), &
vratja(jp),vratjb(jp)
end do
close(u56)
! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs
! used for computations of dP/dV and J_cd
iopt=0

21
src/main.f90
View File

@ -3,7 +3,7 @@ program gray_main
use graycore, only : gray
use gray_params, only : read_inputs,read_params, antctrl_type,eqparam_type, &
prfparam_type,outparam_type,rtrparam_type,hcdparam_type
use beams, only : read_beam0, read_beam1
use beams, only : read_beam0, read_beam1, read_beam2
use equilibrium, only : read_equil_an,read_eqdsk,change_cocos,eq_scal, &
set_rhospl,setqphi_num,frhopolv
use coreprofiles, only : read_profiles_an,read_profiles,tene_scal
@ -69,14 +69,21 @@ program gray_main
deallocate(xrad)
end if
! re-scale input data
call tene_scal(terad,derad,prfp%factte,prfp%factne,eqp%factb,prfp%iscal,prfp%iprof)
call tene_scal(terad,derad,prfp%factte,prfp%factne,eqp%factb,prfp%iscal, &
prfp%iprof)
!------------- antenna --------------
! interpolate beam table if antctrl%ibeam>0
if(antp%ibeam>0) then
call read_beam1(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir)
else
call read_beam0(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir)
end if
select case (antp%ibeam)
case (2)
! to be completed: now 1st beamd always selected, iox read from table
call read_beam2(antp%filenm,1,antp%alpha,antp%beta,fghz,antp%iox,x0,y0,z0, &
w1,w2,ri1,ri2,phiw,phir)
case (1)
call read_beam1(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0, &
w1,w2,ri1,ri2,phiw,phir)
case default
call read_beam0(antp%filenm,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir)
end select
alpha0=antp%alpha
beta0=antp%beta
p0mw=antp%power

24
src/pec.f90
View File

@ -15,18 +15,18 @@ contains
implicit none
! arguments
integer, intent(in) :: ipec
real(wp_), dimension(nnd), intent(in) :: rt_in
real(wp_), dimension(nnd), intent(in), optional :: 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
! rt_in present: read input grid
! else: build equidistant grid dimension nnd
! ipec=+/-1 rho_pol grid
! ipec=+/-2 rho_tor 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))
@ -36,14 +36,16 @@ contains
rtabpsi1(0) = zero
do it=1,nnd
if(ipec > 0) then
if(present(rt_in)) then
! read radial grid from input
rt = rt_in(it)
if(it<nnd) then
drt = rt_in(it+1)-rt
end if
else
! 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
@ -51,7 +53,7 @@ contains
else
rt1 = one
end if
if (abs(ipec) == 1) then
if (ipec == 1) then
rhop_tab(it) = rt
rhot_tab(it) = frhotor(rt)
rhop1 = rt1