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

src: use derived type arguments (work in progress)

Browse Source 
This change structures the arguments of most functions, in particular
gray_main, into well-defined categories using derived types.

All types are defined in the gray_params.f90 (location subject to
change) and are organised as follows:

  gray_parameters (statically allocated data)
  ├── antenna_parameters
  ├── ecrh_cd_parameters
  ├── equilibrium_parameters
  ├── misc_parameters
  ├── output_parameters
  ├── profiles_parameters
  └── raytracing_parameters

  gray_data - inputs of gray_main (dynamically-allocated arrays)
  ├── equilibrium_data
  └── profiles_data

  gray_results - outputs of gray_main (dynamically-allocated arrays)
This commit is contained in:
Michele Guerini Rocco 2021-12-15 02:31:09 +01:00
parent 4f867bad14
commit 948a512254
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
9 changed files with 1400 additions and 1314 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/beamdata.f90
View File

@ -9,10 +9,10 @@ contains
subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
use gray_params, only : rtrparam_type use gray_params, only : raytracing_parameters
use const_and_precisions, only : zero,half,two use const_and_precisions, only : zero,half,two
implicit none implicit none
type(rtrparam_type), intent(in) :: rtrparam type(raytracing_parameters), intent(in) :: rtrparam
real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, & real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, &
gri,psjki,ppabs,ccci gri,psjki,ppabs,ccci
real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
@ -50,6 +50,10 @@ contains
nstep=rtrparam%nstep nstep=rtrparam%nstep
! Allocate for each ray:
! y = (x, k),
! yp = dy/dt,
! etc.
call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
end subroutine init_btr end subroutine init_btr
@ -97,6 +101,9 @@ contains
! Functions to map radial (j), poloidal (k) ray
! indices to a single global index (i)
function rayi2jk(i) result(jk) function rayi2jk(i) result(jk)
implicit none implicit none
integer, intent(in) :: i integer, intent(in) :: i

487
src/beams.f90
View File

@ -4,166 +4,175 @@ module beams
contains contains
subroutine read_beam0(file_beam,fghz,x00,y00,z00, & subroutine read_beam0(params, unit)
wcsi,weta,rcicsi,rcieta,phiw,phir,unit) ! Reads the wave launcher parameters for the simple case
use const_and_precisions, only : pi,vc=>ccgs_ ! where w(z) and 1/R(z) are fixed.
use utils, only : get_free_unit
use const_and_precisions, only : pi, vc=>ccgs_
use gray_params, only : antenna_parameters
use utils, only : get_free_unit
implicit none implicit none
! arguments
character(len=*), intent(in) :: file_beam ! subroutine arguments
real(wp_), intent(out) :: fGHz,x00,y00,z00 type(antenna_parameters), intent(inout) :: params
real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw
integer, intent(in), optional :: unit integer, intent(in), optional :: unit
! local variables
! local variables
integer :: u integer :: u
real(wp_) :: ak0,zrcsi,zreta,w0csi,w0eta,d0csi,d0eta real(wp_) :: ak0,zrcsi,zreta
if (present(unit)) then if (present(unit)) then
u=unit u = unit
else else
u = get_free_unit() u = get_free_unit()
end if end if
open(unit=u,file=trim(file_beam),status='OLD',action='READ')
! fghz wave frequency (GHz) open(unit=u, file=trim(params%filenm), status='OLD', action='READ')
read(u,*) fGHz read(u, *) params%fghz
! x00,y00,z00 coordinates of launching point in cm read(u, *) params%pos
read(u,*) x00, y00, z00 read(u, *) params%w, params%ri, params%phi(1)
! beams parameters in local reference system
! w0 -> waist, d0 -> waist distance from launching point
! phiw angle of spot ellipse
read(u,*) w0csi,w0eta,d0csi,d0eta,phiw
close(u) close(u)
ak0=2.0e9_wp_*pi*fghz/vc ak0 = 2.0e9_wp_* pi * params%fghz / vc
zrcsi=0.5_wp_*ak0*w0csi**2 zrcsi = 0.5_wp_ * ak0 * params%w(1)**2
zreta=0.5_wp_*ak0*w0eta**2 zreta = 0.5_wp_ * ak0 * params%w(2)**2
wcsi=w0csi*sqrt(1.0_wp_+(d0csi/zrcsi)**2)
weta=w0eta*sqrt(1.0_wp_+(d0eta/zreta)**2)
rcicsi=-d0csi/(d0csi**2+zrcsi**2)
rcieta=-d0eta/(d0eta**2+zreta**2)
phir=phiw
params%w(1) = params%w(1) * sqrt(1.0_wp_ + (params%ri(1)/zrcsi)**2)
params%w(2) = params%w(2) * sqrt(1.0_wp_ + (params%ri(2)/zreta)**2)
params%ri(1) = -params%ri(1) / (params%ri(1)**2 + zrcsi**2)
params%ri(2) = -params%ri(2) / (params%ri(2)**2 + zreta**2)
params%phi(2) = params%phi(1)
end subroutine read_beam0 end subroutine read_beam0
subroutine read_beam1(params, unit)
! Reads the wave launcher parameters for the case
! where w(z, α) and 1/R(z, α) depend on the launcher angle α.
subroutine read_beam1(file_beam,alpha0,beta0,fghz,x00,y00,z00, &
wcsi,weta,rcicsi,rcieta,phiw,phir,unit)
use const_and_precisions, only : pi,vc=>ccgs_ use const_and_precisions, only : pi,vc=>ccgs_
use simplespline, only : spli, difcs use gray_params, only : antenna_parameters
use utils, only : get_free_unit,locate use simplespline, only : spli, difcs
use utils, only : get_free_unit,locate
implicit none implicit none
! arguments
character(len=*), intent(in) :: file_beam ! subroutine arguments
real(wp_), intent(inout) :: alpha0 type(antenna_parameters), intent(inout) :: params
real(wp_), intent(out) :: fghz,x00,y00,z00,beta0
real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw
integer, intent(in), optional :: unit integer, intent(in), optional :: unit
! local variables
integer :: u,iopt,ier,nisteer,i,k,ii ! local variables
real(wp_) :: steer,dal integer :: u, iopt, ier, nisteer, i, k, ii
real(wp_), dimension(:), allocatable :: alphastv,betastv,x00v,y00v, & real(wp_) :: steer, dal
z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v, & real(wp_), dimension(:), allocatable :: &
cbeta,cx0,cy0,cz0,cwaist1,cwaist2, & alphastv, betastv, x00v, y00v, &
crci1,crci2,cphi1,cphi2 z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, &
cbeta, cx0, cy0, cz0, cwaist1, cwaist2, &
crci1, crci2, cphi1, cphi2
if (present(unit)) then if (present(unit)) then
u=unit u = unit
else else
u = get_free_unit() u = get_free_unit()
end if end if
open(unit=u,file=file_beam,status='OLD',action='READ')
read(u,*) fghz open(unit=u, file=params%filenm, status='OLD', action='READ')
read(u,*) params%fghz
read(u,*) nisteer read(u,*) nisteer
allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), & allocate(alphastv(nisteer), betastv(nisteer), waist1v(nisteer), &
waist2v(nisteer),rci1v(nisteer),rci2v(nisteer), & waist2v(nisteer), rci1v(nisteer), rci2v(nisteer), &
phi1v(nisteer),phi2v(nisteer),x00v(nisteer), & phi1v(nisteer), phi2v(nisteer), x00v(nisteer), &
y00v(nisteer),z00v(nisteer),cbeta(4*nisteer), & y00v(nisteer), z00v(nisteer), cbeta(4*nisteer), &
cx0(4*nisteer),cy0(4*nisteer),cz0(4*nisteer), & cx0(4*nisteer), cy0(4*nisteer), cz0(4*nisteer), &
cwaist1(4*nisteer),cwaist2(4*nisteer),crci1(4*nisteer), & cwaist1(4*nisteer), cwaist2(4*nisteer), crci1(4*nisteer), &
crci2(4*nisteer),cphi1(4*nisteer),cphi2(4*nisteer)) crci2(4*nisteer), cphi1(4*nisteer), cphi2(4*nisteer))
do i=1,nisteer do i=1,nisteer
read(u,*) steer,alphastv(i),betastv(i),x00v(i),y00v(i),z00v(i), & read(u, *) steer, alphastv(i), betastv(i), &
waist1v(i),waist2v(i),rci1v(i),rci2v(i),phi1v(i),phi2v(i) x00v(i), y00v(i), z00v(i), &
waist1v(i), waist2v(i), &
rci1v(i), rci2v(i), &
phi1v(i), phi2v(i)
end do end do
close(u) 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)
call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier)
call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier)
call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier)
call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier)
call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier)
call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier)
call difcs(alphastv,x00v,nisteer,iopt,cx0,ier)
call difcs(alphastv,y00v,nisteer,iopt,cy0,ier)
call difcs(alphastv,z00v,nisteer,iopt,cz0,ier)
if((alpha0 > alphastv(1)).and.(alpha0 < alphastv(nisteer))) then ! initial beam data measured in mm -> transformed to cm
call locate(alphastv,nisteer,alpha0,k) x00v = 0.1_wp_ * x00v
dal=alpha0-alphastv(k) y00v = 0.1_wp_ * y00v
beta0=spli(cbeta,nisteer,k,dal) z00v = 0.1_wp_ * z00v
x00=spli(cx0,nisteer,k,dal) waist1v = 0.1_wp_ * waist1v
y00=spli(cy0,nisteer,k,dal) waist2v = 0.1_wp_ * waist2v
z00=spli(cz0,nisteer,k,dal) rci1v = 10._wp_ * rci1v
wcsi=spli(cwaist1,nisteer,k,dal) rci2v = 10._wp_ * rci2v
weta=spli(cwaist2,nisteer,k,dal)
rcicsi=spli(crci1,nisteer,k,dal) iopt = 0
rcieta=spli(crci2,nisteer,k,dal) call difcs(alphastv, betastv, nisteer, iopt, cbeta, ier)
phiw=spli(cphi1,nisteer,k,dal) call difcs(alphastv, waist1v, nisteer, iopt, cwaist1, ier)
phir=spli(cphi2,nisteer,k,dal) call difcs(alphastv, rci1v, nisteer, iopt, crci1, ier)
call difcs(alphastv, waist2v, nisteer, iopt, cwaist2, ier)
call difcs(alphastv, rci2v, nisteer, iopt, crci2, ier)
call difcs(alphastv, phi1v, nisteer, iopt, cphi1, ier)
call difcs(alphastv, phi2v, nisteer, iopt, cphi2, ier)
call difcs(alphastv, x00v, nisteer, iopt, cx0, ier)
call difcs(alphastv, y00v, nisteer, iopt, cy0, ier)
call difcs(alphastv, z00v, nisteer, iopt, cz0, ier)
if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then
call locate(alphastv, nisteer, params%alpha , k)
dal = params%alpha - alphastv(k)
params%beta = spli(cbeta, nisteer, k, dal)
params%pos(1) = spli(cx0, nisteer, k, dal)
params%pos(2) = spli(cy0, nisteer, k, dal)
params%pos(3) = spli(cz0, nisteer, k, dal)
params%w(1) = spli(cwaist1, nisteer, k, dal)
params%w(2) = spli(cwaist2, nisteer, k, dal)
params%ri(1) = spli(crci1, nisteer, k, dal)
params%ri(2) = spli(crci2, nisteer, k, dal)
params%phi(1) = spli(cphi1, nisteer, k, dal)
params%phi(2) = spli(cphi2, nisteer, k, dal)
else else
! alpha0 outside table range ! params%alpha outside table range
if(alpha0 >= alphastv(nisteer)) ii=nisteer if(params%alpha >= alphastv(nisteer)) ii=nisteer
if(alpha0 <= alphastv(1)) ii=1 if(params%alpha <= alphastv(1)) ii=1
alpha0=alphastv(ii) params%alpha = alphastv(ii)
beta0=betastv(ii) params%beta = betastv(ii)
x00=x00v(ii) params%pos(1) = x00v(ii)
y00=y00v(ii) params%pos(2) = y00v(ii)
z00=z00v(ii) params%pos(3) = z00v(ii)
wcsi=waist1v(ii) params%w(1) = waist1v(ii)
weta=waist2v(ii) params%w(2) = waist2v(ii)
rcicsi=rci1v(ii) params%ri(1) = rci1v(ii)
rcieta=rci2v(ii) params%ri(2) = rci2v(ii)
phiw=phi1v(ii) params%phi(1) = phi1v(ii)
phir=phi2v(ii) params%phi(2) = phi2v(ii)
end if end if
deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v, & 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) cx0, cy0, cz0, cwaist1, cwaist2, &
crci1, crci2, cphi1, cphi2)
end subroutine read_beam1 end subroutine read_beam1
subroutine read_beam2(file_beam,beamid,alpha0,beta0,fghz,iox,x00,y00,z00, & subroutine read_beam2(params, beamid, unit)
wcsi,weta,rcicsi,rcieta,phiw,phir,unit) ! Reads the wave launcher parameters for the general case
use utils, only : get_free_unit, intlin, locate ! where w(z, α, β) and 1/R(z, α, β) depend on the launcher angles α, β.
use gray_params, only : antenna_parameters
use utils, only : get_free_unit, intlin, locate
use reflections, only : inside use reflections, only : inside
use dierckx, only : curfit, splev, surfit, bispev use dierckx, only : curfit, splev, surfit, bispev
implicit none implicit none
character(len=*), intent(in) :: file_beam
! subroutine arguments
type(antenna_parameters), intent(inout) :: params
integer, intent(in) :: beamid 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 integer, intent(in), optional :: unit
! local variables
character(len=20) :: beamname character(len=20) :: beamname
integer :: u integer :: u
integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta
@ -189,11 +198,12 @@ contains
real(wp_), parameter :: sspl=0.01_wp_ real(wp_), parameter :: sspl=0.01_wp_
if (present(unit)) then if (present(unit)) then
u=unit u = unit
else else
u = get_free_unit() u = get_free_unit()
end if end if
open(unit=u,file=file_beam,status='OLD',action='READ')
open(unit=u, file=params%filenm, status='OLD', action='READ')
!======================================================================================= !=======================================================================================
! # of beams ! # of beams
read(u,*) nbeam read(u,*) nbeam
@ -202,13 +212,13 @@ contains
jumprow=0 jumprow=0
! c==================================================================================== ! c====================================================================================
do i=1,beamid-1 do i=1,beamid-1
read(u,*) beamname, iox, fghz, nalpha, nbeta read(u,*) beamname, params%iox, params%fghz, nalpha, nbeta
jumprow = jumprow+nalpha*nbeta jumprow = jumprow+nalpha*nbeta
end do end do
! c==================================================================================== ! c====================================================================================
! !
! beam of interest ! beam of interest
read(u,*) beamname, iox, fghz, nalpha, nbeta read(u,*) beamname, params%iox, params%fghz, nalpha, nbeta
! !
! c==================================================================================== ! c====================================================================================
! unused beams' data grids ! unused beams' data grids
@ -216,7 +226,7 @@ contains
read(u,*) beamname read(u,*) beamname
end do end do
do i=1,jumprow do i=1,jumprow
read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 read(u,*) alphast,betast,params%pos(1),params%pos(2),params%pos(3),waist1,waist2,rci1,rci2,phi1,phi2
end do end do
! c==================================================================================== ! c====================================================================================
! !
@ -231,12 +241,12 @@ contains
! c==================================================================================== ! c====================================================================================
! beam data grid reading ! beam data grid reading
do i=1,nisteer do i=1,nisteer
read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 read(u,*) alphast,betast,params%pos(1),params%pos(2),params%pos(3),waist1,waist2,rci1,rci2,phi1,phi2
! !
! initial beam data (x00, y00, z00) are measured in mm -> transformed to cm ! initial beam data (params%pos(1), params%pos(2), params%pos(3)) are measured in mm -> transformed to cm
x00v(i)=0.1d0*x00 x00v(i)=0.1d0*params%pos(1)
y00v(i)=0.1d0*y00 y00v(i)=0.1d0*params%pos(2)
z00v(i)=0.1d0*z00 z00v(i)=0.1d0*params%pos(3)
alphastv(i)=alphast alphastv(i)=alphast
betastv(i)=betast betastv(i)=betast
waist1v(i)=0.1d0*waist1 waist1v(i)=0.1d0*waist1
@ -259,17 +269,17 @@ contains
! !
! no free variables ! no free variables
if(fdeg.eq.3) then if(fdeg.eq.3) then
alpha0=alphastv(1) params%alpha=alphastv(1)
beta0=betastv(1) params%beta=betastv(1)
x00=x00v(1) params%pos(1)=x00v(1)
y00=y00v(1) params%pos(2)=y00v(1)
z00=z00v(1) params%pos(3)=z00v(1)
wcsi=waist1v(1) params%w(1)=waist1v(1)
weta=waist2v(1) params%w(2)=waist2v(1)
rcicsi=rci1v(1) params%ri(1)=rci1v(1)
rcieta=rci2v(1) params%ri(2)=rci2v(1)
phiw=phi1v(1) params%phi(2)=phi1v(1)
phir=phi2v(1) params%phi(1)=phi2v(1)
deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, &
phi2v,x00v,y00v,z00v,xcoord,ycoord) phi2v,x00v,y00v,z00v,xcoord,ycoord)
return return
@ -283,8 +293,8 @@ contains
! alpha = dependent variable ! alpha = dependent variable
xcoord = betastv xcoord = betastv
ycoord = alphastv ycoord = alphastv
xcoord0 = beta0 xcoord0 = params%beta
ycoord0 = alpha0 ycoord0 = params%alpha
kx=min(nbeta-1,kspl) kx=min(nbeta-1,kspl)
! c==================================================================================== ! c====================================================================================
else else
@ -293,8 +303,8 @@ contains
! beta = dependent/independent (fdeg = 1/0) ! beta = dependent/independent (fdeg = 1/0)
xcoord = alphastv xcoord = alphastv
ycoord = betastv ycoord = betastv
xcoord0 = alpha0 xcoord0 = params%alpha
ycoord0 = beta0 ycoord0 = params%beta
nxcoord = nalpha nxcoord = nalpha
nycoord = nbeta nycoord = nbeta
kx=min(nalpha-1,kspl) kx=min(nalpha-1,kspl)
@ -461,23 +471,23 @@ contains
call splev(txycoord,nxycoord,cycoord,kx,(/xcoord0/),fi,1,ier) call splev(txycoord,nxycoord,cycoord,kx,(/xcoord0/),fi,1,ier)
ycoord0=fi(1) ycoord0=fi(1)
call splev(txwaist1,nxwaist1,cwaist1,kx,(/xcoord0/),fi,1,ier) call splev(txwaist1,nxwaist1,cwaist1,kx,(/xcoord0/),fi,1,ier)
wcsi=fi(1) params%w(1)=fi(1)
call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier) call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier)
weta=fi(1) params%w(2)=fi(1)
call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier) call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier)
rcicsi=fi(1) params%ri(1)=fi(1)
call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier) call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier)
rcieta=fi(1) params%ri(2)=fi(1)
call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier) call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier)
phiw=fi(1) params%phi(2)=fi(1)
call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier) call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier)
phir=fi(1) params%phi(1)=fi(1)
call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier) call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier)
x00=fi(1) params%pos(1)=fi(1)
call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier) call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier)
y00=fi(1) params%pos(2)=fi(1)
call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier) call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier)
z00=fi(1) params%pos(3)=fi(1)
! c---------------------------------------------------------------------------------- ! c----------------------------------------------------------------------------------
else else
! c---------------------------------------------------------------------------------- ! c----------------------------------------------------------------------------------
@ -486,15 +496,15 @@ contains
! !
xcoord0=xcoord(ii) xcoord0=xcoord(ii)
ycoord0=ycoord(ii) ycoord0=ycoord(ii)
x00=x00v(ii) params%pos(1)=x00v(ii)
y00=y00v(ii) params%pos(2)=y00v(ii)
z00=z00v(ii) params%pos(3)=z00v(ii)
wcsi=waist1v(ii) params%w(1)=waist1v(ii)
weta=waist2v(ii) params%w(2)=waist2v(ii)
rcicsi=rci1v(ii) params%ri(1)=rci1v(ii)
rcieta=rci2v(ii) params%ri(2)=rci2v(ii)
phiw=phi1v(ii) params%phi(2)=phi1v(ii)
phir=phi2v(ii) params%phi(1)=phi2v(ii)
end if end if
! c==================================================================================== ! c====================================================================================
else else
@ -586,53 +596,53 @@ contains
! 5: xcoord0 unchanged, ycoord0 moved on side C ! 5: xcoord0 unchanged, ycoord0 moved on side C
! 7: xcoord0 moved on side D, ycoord0 unchanged ! 7: xcoord0 moved on side D, ycoord0 unchanged
! 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates ! 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 ! in 1,3,5,7 incheck is set back to 1 to evaluate params%pos(1),params%pos(2),params%pos(3),waist,rci,phi in
! new (xcoord0,ycoord0) ! new (xcoord0,ycoord0)
! in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the ! in 2,4,6,8 incheck remains 0 and params%pos(1),params%pos(2),params%pos(3),waist,rci,phi values at the
! (xcoord0,ycoord0) vertex are used ! (xcoord0,ycoord0) vertex are used
alpha0 = xcoord0 params%alpha = xcoord0
beta0 = ycoord0 params%beta = ycoord0
SELECT CASE (in) SELECT CASE (in)
CASE (1) CASE (1)
! beta0 outside table range ! params%beta outside table range
! locate position of xcoord0 with respect to x coordinates of side A ! locate position of xcoord0 with respect to x coordinates of side A
call locate(xpolygA,nxcoord,xcoord0,ii) call locate(xpolygA,nxcoord,xcoord0,ii)
! find corresponding y value on side A for xcoord position ! find corresponding y value on side A for xcoord position
call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0) call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0)
incheck = 1 incheck = 1
CASE (2) CASE (2)
! alpha0 and beta0 outside table range ! params%alpha and params%beta outside table range
! xcoord0, ycoord0 set ! xcoord0, ycoord0 set
xcoord0 = xvert(2) xcoord0 = xvert(2)
ycoord0 = yvert(2) ycoord0 = yvert(2)
ii = nxcoord !indice per assegnare valori waist, rci, phi ii = nxcoord !indice per assegnare valori waist, rci, phi
CASE (3) CASE (3)
! alpha0 outside table range ! params%alpha outside table range
call locate(ypolygB,nycoord,ycoord0,ii) call locate(ypolygB,nycoord,ycoord0,ii)
call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0) call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0)
incheck = 1 incheck = 1
CASE (4) CASE (4)
! alpha0 and beta0 outside table range ! params%alpha and params%beta outside table range
xcoord0 = xvert(3) xcoord0 = xvert(3)
ycoord0 = yvert(3) ycoord0 = yvert(3)
ii = nxcoord+nycoord-1 ii = nxcoord+nycoord-1
CASE (5) CASE (5)
! beta0 outside table range ! params%beta outside table range
call locate(xpolygC,nxcoord,xcoord0,ii) call locate(xpolygC,nxcoord,xcoord0,ii)
call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0) call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0)
incheck = 1 incheck = 1
CASE (6) CASE (6)
! alpha0 and beta0 outside table range ! params%alpha and params%beta outside table range
xcoord0 = xvert(4) xcoord0 = xvert(4)
ycoord0 = yvert(4) ycoord0 = yvert(4)
ii = 2*nxcoord+nycoord-2 ii = 2*nxcoord+nycoord-2
CASE (7) CASE (7)
! alpha0 outside table range ! params%alpha outside table range
call locate(ypolygD,nycoord,ycoord0,ii) call locate(ypolygD,nycoord,ycoord0,ii)
call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0) call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0)
incheck = 1 incheck = 1
CASE (8) CASE (8)
! alpha0 and beta0 outside table range ! params%alpha and params%beta outside table range
xcoord0 = xvert(1) xcoord0 = xvert(1)
ycoord0 = yvert(1) ycoord0 = yvert(1)
ii = 1 ii = 1
@ -651,44 +661,44 @@ contains
allocate(wrk(lwrk),iwrk(kwrk)) allocate(wrk(lwrk),iwrk(kwrk))
call bispev(txwaist1,nxwaist1,tywaist1,nywaist1,cwaist1, & call bispev(txwaist1,nxwaist1,tywaist1,nywaist1,cwaist1, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
wcsi=fi(1) params%w(1)=fi(1)
call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, & call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
weta=fi(1) params%w(2)=fi(1)
call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1, & call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
rcicsi=fi(1) params%ri(1)=fi(1)
call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2, & call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
rcieta=fi(1) params%ri(2)=fi(1)
call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1, & call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
phiw=fi(1) params%phi(2)=fi(1)
call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2, & call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
phir=fi(1) params%phi(1)=fi(1)
call bispev(txx0,nxx0,tyx0,nyx0,cx0, & call bispev(txx0,nxx0,tyx0,nyx0,cx0, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
x00=fi(1) params%pos(1)=fi(1)
call bispev(txy0,nxy0,tyy0,nyy0,cy0, & call bispev(txy0,nxy0,tyy0,nyy0,cy0, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
y00=fi(1) params%pos(2)=fi(1)
call bispev(txz0,nxz0,tyz0,nyz0,cz0, & call bispev(txz0,nxz0,tyz0,nyz0,cz0, &
kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier)
z00=fi(1) params%pos(3)=fi(1)
deallocate(wrk,iwrk) deallocate(wrk,iwrk)
! c---------------------------------------------------------------------------------- ! c----------------------------------------------------------------------------------
else else
! c---------------------------------------------------------------------------------- ! c----------------------------------------------------------------------------------
x00=x00v(ii) params%pos(1)=x00v(ii)
y00=y00v(ii) params%pos(2)=y00v(ii)
z00=z00v(ii) params%pos(3)=z00v(ii)
wcsi=waist1v(ii) params%w(1)=waist1v(ii)
weta=waist2v(ii) params%w(2)=waist2v(ii)
rcicsi=rci1v(ii) params%ri(1)=rci1v(ii)
rcieta=rci2v(ii) params%ri(2)=rci2v(ii)
phiw=phi1v(ii) params%phi(2)=phi1v(ii)
phir=phi2v(ii) params%phi(1)=phi2v(ii)
end if end if
! c==================================================================================== ! c====================================================================================
end if end if
@ -709,11 +719,11 @@ contains
!####################################################################################### !#######################################################################################
! set correct values for alpha, beta ! set correct values for alpha, beta
if(fdeg.eq.2) then if(fdeg.eq.2) then
alpha0 = ycoord0 params%alpha = ycoord0
beta0 = xcoord0 params%beta = xcoord0
else else
alpha0 = xcoord0 params%alpha = xcoord0
beta0 = ycoord0 params%beta = ycoord0
end if end if
!####################################################################################### !#######################################################################################
deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, &
@ -722,53 +732,62 @@ contains
end subroutine read_beam2 end subroutine read_beam2
subroutine launchangles2n(alpha,beta,xv,anv) subroutine launchangles2n(params, anv)
use const_and_precisions, only : degree ! Given the wave launcher `params` computes the initial
implicit none ! wavevector `anv`, defined by n̅ = ck̅/ω, in cartesian coordinates.
! arguments
real(wp_), intent(in) :: alpha,beta,xv(3)
real(wp_), intent(out) :: anv(3)
! local variables
real(wp_) :: r,anr,anphi,a,b
r=sqrt(xv(1)**2+xv(2)**2) use const_and_precisions, only : degree
! phi=atan2(y,x) use gray_params, only : antenna_parameters
a = degree*alpha
b = degree*beta implicit none
!
! angles alpha, beta in a local reference system as proposed by Gribov et al ! subroutine arguments
! type(antenna_parameters), intent(in) :: params
real(wp_), intent(out) :: anv(3)
! local variables
real(wp_) :: r, anr, anphi, a, b
r = sqrt(params%pos(1)**2 + params%pos(2)**2)
a = degree*params%alpha
b = degree*params%beta
! Angles α, β in a local reference system
! as proposed by Gribov et al.
anr = -cos(b)*cos(a) anr = -cos(b)*cos(a)
anphi = sin(b) anphi = sin(b)
! anx = -cos(b)*cos(a)
! any = sin(b)
anv(1) = (anr*xv(1) - anphi*xv(2))/r ! = anx anv(1) = (anr*params%pos(1) - anphi*params%pos(2))/r ! = anx
anv(2) = (anr*xv(2) + anphi*xv(1))/r ! = any anv(2) = (anr*params%pos(2) + anphi*params%pos(1))/r ! = any
! anr = (anx*xv(1) + any*xv(2))/r anv(3) = -cos(b)*sin(a) ! = anz
! anphi = (any*xv(1) - anx*xv(2))/r
anv(3) =-cos(b)*sin(a) ! = anz
end subroutine launchangles2n end subroutine launchangles2n
subroutine xgygcoeff(fghz,ak0,bres,xgcn)
use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_,pi,wce1_ subroutine xgygcoeff(fghz, ak0, bres, xgcn)
! Given the EC wave frequency computes:
!
! 1. vacuum wavevector `k0` (k₀ = ω/c),
! 2. resonant magnetic field `bres` (qB/m = ω),
! 3. adimensional `xgcn` parameter (X = ω_p²/ω² = nq²/ε₀mω²).
use const_and_precisions, only : qe=>ecgs_, me=>mecgs_, &
vc=>ccgs_, pi, wce1_
implicit none implicit none
! arguments
real(wp_), intent(in) :: fghz ! subroutine arguments
real(wp_), intent(out) :: ak0,bres,xgcn real(wp_), intent(in) :: fghz
! local variables real(wp_), intent(out) :: ak0, bres, xgcn
! local variables
real(wp_) :: omega real(wp_) :: omega
omega=2.0e9_wp_*pi*fghz ! [rad/s] omega = 2.0e9_wp_*pi*fghz ! [rad/s]
ak0=omega/vc ! [rad/cm] ak0 = omega/vc ! [rad/cm]
!
! yg=btot/bres ! yg = btot/bres
! bres = omega/wce1_ ! [T]
bres=omega/wce1_ ! [T]
! ! xg = xgcn*dens19
! xg=xgcn*dens19 xgcn = 4.0e13_wp_ * pi * qe**2/(me * omega**2) ! [10^-19 m^3]
!
xgcn=4.0e13_wp_*pi*qe**2/(me*omega**2) ! [10^-19 m^3]
end subroutine xgygcoeff end subroutine xgygcoeff
end module beams end module beams

568
src/equilibrium.f90
View File

@ -13,7 +13,7 @@ module equilibrium
! === 2d spline psi(r,z), normalization and derivatives ========== ! === 2d spline psi(r,z), normalization and derivatives ==========
integer, save :: nsr, nsz integer, save :: nsr, nsz
real(wp_), save :: psia, psiant, psinop, psnbd real(wp_), save :: psia, psiant, psinop, psnbnd
real(wp_), dimension(:), allocatable, save :: tr,tz real(wp_), dimension(:), allocatable, save :: tr,tz
real(wp_), dimension(:), allocatable, save :: cceq, cceq01, cceq10, & real(wp_), dimension(:), allocatable, save :: cceq, cceq01, cceq10, &
cceq20, cceq02, cceq11 cceq20, cceq02, cceq11
@ -34,121 +34,127 @@ module equilibrium
contains contains
subroutine read_eqdsk(filenm,rv,zv,psin,psia,psinr,fpol,q,rvac,rax,zax, & subroutine read_eqdsk(params, data, unit)
rbnd,zbnd,rlim,zlim,ipsinorm,idesc,ifreefmt,unit) ! Reads the MHD equilibrium `data` from a G-EQDSK file (params%filenm).
! If given, the file is opened in the `unit` number.
! For a description of the G-EQDSK, see the GRAY user manual.
use const_and_precisions, only : one use const_and_precisions, only : one
use utils, only : get_free_unit use gray_params, only : equilibrium_parameters, equilibrium_data
use utils, only : get_free_unit
implicit none implicit none
! arguments
character(len=*), intent(in) :: filenm ! subroutine arguments
real(wp_), intent(out) :: psia,rvac,rax,zax type(equilibrium_parameters), intent(in) :: params
real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,psinr,fpol,q type(equilibrium_data), intent(out) :: data
real(wp_), dimension(:), allocatable, intent(out) :: rbnd,zbnd,rlim,zlim integer, optional, intent(in) :: unit
real(wp_), dimension(:,:), allocatable, intent(out) :: psin
integer, optional, intent(in) :: ipsinorm,idesc,ifreefmt,unit ! local variables
! local variables integer :: u, idum, i, j, nr, nz, nbnd, nlim
integer, parameter :: indef=0,iddef=1,iffdef=0
integer :: in,id,iffmt,u,idum,i,j,nr,nz,nbnd,nlim
character(len=48) :: string character(len=48) :: string
real(wp_) :: dr,dz,dps,rleft,zmid,zleft,xdum,psiedge,psiaxis real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis
real(wp_) :: xdum ! dummy variable, used to discard data
! set default values if optional arguments are absent
in=indef; id=iddef; iffmt=iffdef if(present(unit)) then
if(present(ipsinorm)) in=ipsinorm u = unit
if(present(idesc)) id=idesc
if(present(ifreefmt)) iffmt=ifreefmt
if (present(unit)) then
u=unit
else else
u=get_free_unit() u = get_free_unit()
end if end if
! open G EQDSK file (see http://fusion.gat.com/efit/g_eqdsk.html) ! Open the G-EQDSK file
open(file=trim(filenm),status='old',action='read',unit=u) open(file=trim(params%filenm), status='old', action='read', unit=u)
! get size of main arrays and allocate them ! get size of main arrays and allocate them
if (id==1) then if (params%idesc == 1) then
read (u,'(a48,3i4)') string,idum,nr,nz read (u,'(a48,3i4)') string,idum,nr,nz
else else
read (u,*) nr,nz read (u,*) nr, nz
end if end if
if (allocated(rv)) deallocate(rv) if (allocated(data%rv)) deallocate(data%rv)
if (allocated(zv)) deallocate(zv) if (allocated(data%zv)) deallocate(data%zv)
if (allocated(psin)) deallocate(psin) if (allocated(data%psin)) deallocate(data%psin)
if (allocated(psinr)) deallocate(psinr) if (allocated(data%psinr)) deallocate(data%psinr)
if (allocated(fpol)) deallocate(fpol) if (allocated(data%fpol)) deallocate(data%fpol)
if (allocated(q)) deallocate(q) if (allocated(data%qpsi)) deallocate(data%qpsi)
allocate(rv(nr),zv(nz),psin(nr,nz),psinr(nr),fpol(nr),q(nr)) allocate(data%rv(nr), data%zv(nz), &
data%psin(nr, nz), &
data%psinr(nr), &
data%fpol(nr), &
data%qpsi(nr))
! store 0D data and main arrays ! Store 0D data and main arrays
if (iffmt==1) then if (params%ifreefmt==1) then
read (u,*) dr,dz,rvac,rleft,zmid read (u, *) dr, dz, data%rvac, rleft, zmid
read (u,*) rax,zax,psiaxis,psiedge,xdum read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum
read (u,*) xdum,xdum,xdum,xdum,xdum read (u, *) xdum, xdum, xdum, xdum, xdum
read (u,*) xdum,xdum,xdum,xdum,xdum read (u, *) xdum, xdum, xdum, xdum, xdum
read (u,*) (fpol(i),i=1,nr) read (u, *) (data%fpol(i), i=1,nr)
read (u,*) (xdum,i=1,nr) read (u, *) (xdum,i=1, nr)
read (u,*) (xdum,i=1,nr) read (u, *) (xdum,i=1, nr)
read (u,*) (xdum,i=1,nr) read (u, *) (xdum,i=1, nr)
read (u,*) ((psin(i,j),i=1,nr),j=1,nz) read (u, *) ((data%psin(i,j), i=1,nr), j=1,nz)
read (u,*) (q(i),i=1,nr) read (u, *) (data%qpsi(i), i=1,nr)
else else
read (u,'(5e16.9)') dr,dz,rvac,rleft,zmid read (u, '(5e16.9)') dr,dz,data%rvac,rleft,zmid
read (u,'(5e16.9)') rax,zax,psiaxis,psiedge,xdum read (u, '(5e16.9)') data%rax,data%zax,psiaxis,psiedge,xdum
read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
read (u,'(5e16.9)') (fpol(i),i=1,nr) read (u, '(5e16.9)') (data%fpol(i),i=1,nr)
read (u,'(5e16.9)') (xdum,i=1,nr) read (u, '(5e16.9)') (xdum,i=1,nr)
read (u,'(5e16.9)') (xdum,i=1,nr) read (u, '(5e16.9)') (xdum,i=1,nr)
read (u,'(5e16.9)') (xdum,i=1,nr) read (u, '(5e16.9)') (xdum,i=1,nr)
read (u,'(5e16.9)') ((psin(i,j),i=1,nr),j=1,nz) read (u, '(5e16.9)') ((data%psin(i,j),i=1,nr),j=1,nz)
read (u,'(5e16.9)') (q(i),i=1,nr) read (u, '(5e16.9)') (data%qpsi(i),i=1,nr)
end if end if
! get size of boundary and limiter arrays and allocate them ! Get size of boundary and limiter arrays and allocate them
read (u,*) nbnd,nlim read (u,*) nbnd, nlim
if (allocated(rbnd)) deallocate(rbnd) if (allocated(data%rbnd)) deallocate(data%rbnd)
if (allocated(zbnd)) deallocate(zbnd) if (allocated(data%zbnd)) deallocate(data%zbnd)
if (allocated(rlim)) deallocate(rlim) if (allocated(data%rlim)) deallocate(data%rlim)
if (allocated(zlim)) deallocate(zlim) if (allocated(data%zlim)) deallocate(data%zlim)
! store boundary and limiter data ! Load plasma boundary data
if(nbnd>0) then if(nbnd > 0) then
allocate(rbnd(nbnd),zbnd(nbnd)) allocate(data%rbnd(nbnd), data%zbnd(nbnd))
if (iffmt==1) then if (params%ifreefmt == 1) then
read(u,*) (rbnd(i),zbnd(i),i=1,nbnd) read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd)
else else
read(u,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd) read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd)
end if
end if
if(nlim>0) then
allocate(rlim(nlim),zlim(nlim))
if (iffmt==1) then
read(u,*) (rlim(i),zlim(i),i=1,nlim)
else
read(u,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim)
end if end if
end if end if
! reading of G EQDSK file completed ! Load limiter data
if(nlim > 0) then
allocate(data%rlim(nlim), data%zlim(nlim))
if (params%ifreefmt == 1) then
read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim)
else
read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim)
end if
end if
! End of G-EQDSK file
close(u) close(u)
! build rv,zv,psinr arrays and normalize psin ! Build rv,zv,psinr arrays
zleft=zmid-0.5_wp_*dz zleft = zmid-0.5_wp_*dz
dr=dr/(nr-1) dr = dr/(nr-1)
dz=dz/(nz-1) dz = dz/(nz-1)
dps=one/(nr-1) dps = one/(nr-1)
do i=1,nr do i=1,nr
psinr(i)=(i-1)*dps data%psinr(i) = (i-1)*dps
rv(i)=rleft+(i-1)*dr data%rv(i) = rleft + (i-1)*dr
end do end do
do i=1,nz do i=1,nz
zv(i)=zleft+(i-1)*dz data%zv(i) = zleft + (i-1)*dz
end do end do
psia=psiedge-psiaxis
if(in==0) psin=(psin-psiaxis)/psia
end subroutine read_eqdsk
! Normalize psin
data%psia = psiedge - psiaxis
if(params%ipsinorm == 0) data%psin = (data%psin - psiaxis)/data%psia
end subroutine read_eqdsk
subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit) subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit)
@ -211,191 +217,228 @@ contains
close(u) close(u)
end subroutine read_equil_an end subroutine read_equil_an
subroutine change_cocos(psia,fpol,q,cocosin,cocosout,ierr)
use const_and_precisions, only : zero,one,pi subroutine change_cocos(data, cocosin, cocosout, error)
! Convert the MHD equilibrium data from one coordinate convention
! (COCOS) to another. These are specified by `cocosin` and
! `cocosout`, respectively.
!
! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010
use const_and_precisions, only : zero, one, pi
use gray_params, only : equilibrium_data
implicit none implicit none
! arguments
real(wp_), intent(inout) :: psia ! subroutine arguments
real(wp_), dimension(:), intent(inout) :: fpol,q type(equilibrium_data), intent(inout) :: data
integer, intent(in) :: cocosin, cocosout integer, intent(in) :: cocosin, cocosout
integer, intent(out), optional :: ierr integer, intent(out), optional :: error
! local variables
real(wp_) :: isign,bsign
integer :: exp2pi,exp2piout
logical :: phiccw,psiincr,qpos,phiccwout,psiincrout,qposout
call decode_cocos(cocosin,exp2pi,phiccw,psiincr,qpos) ! local variables
call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout) real(wp_) :: isign, bsign
integer :: exp2pi, exp2piout
logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout
! check sign consistency call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos)
isign=sign(one,psia) call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout)
if (.not.psiincr) isign=-isign
bsign=sign(one,fpol(size(fpol))) ! Check sign consistency
if (qpos.neqv.isign*bsign*q(size(q))>zero) then isign = sign(one, data%psia)
! warning: sign inconsistency found among q, Ipla and Bref if (.not.psiincr) isign = -isign
q=-q bsign = sign(one, data%fpol(size(data%fpol)))
if(present(ierr)) ierr=1 if (qpos .neqv. isign * bsign * data%qpsi(size(data%qpsi)) > zero) then
! Warning: sign inconsistency found among q, Ipla and Bref
data%qpsi = -data%qpsi
if (present(error)) error = 1
else else
if(present(ierr)) ierr=0 if (present(error)) error = 0
end if end if
! convert cocosin to cocosout ! Convert cocosin to cocosout
! opposite direction of toroidal angle phi in cocosin and cocosout ! Opposite direction of toroidal angle phi in cocosin and cocosout
if (phiccw.neqv.phiccwout) fpol=-fpol if (phiccw .neqv. phiccwout) data%fpol = -data%fpol
! q has opposite sign for given sign of Bphi*Ip
if (qpos .neqv. qposout) q=-q ! q has opposite sign for given sign of Bphi*Ip
! psi and Ip signs don't change accordingly if (qpos .neqv. qposout) data%qpsi = -data%qpsi
if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) psia=-psia
! convert Wb to Wb/rad or viceversa ! psi and Ip signs don't change accordingly
psia=psia*(2.0_wp_*pi)**(exp2piout-exp2pi) if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) &
data%psia = -data%psia
! Convert Wb to Wb/rad or viceversa
data%psia = data%psia * (2.0_wp_*pi)**(exp2piout - exp2pi)
end subroutine change_cocos end subroutine change_cocos
subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
implicit none
integer, intent(in) :: cocos
integer, intent(out) :: exp2pi
logical, intent(out) :: phiccw,psiincr,qpos
integer :: cmod10,cmod4
cmod10=mod(cocos,10) subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos)
cmod4=mod(cmod10,4) implicit none
! cocos>10 psi in Wb, cocos<10 psi in Wb/rad
exp2pi=(cocos-cmod10)/10 ! subroutine arguments
! cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW integer, intent(in) :: cocos
phiccw=(mod(cmod10,2)==1) integer, intent(out) :: exp2pi
! cocos mod 10 = 1,2,5,6: psi increasing with positive Ip logical, intent(out) :: phiccw, psiincr, qpos
psiincr=(cmod4==1 .or. cmod4==2)
! cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip ! local variables
qpos=(cmod10<3 .or. cmod10>6) integer :: cmod10, cmod4
cmod10 = mod(cocos, 10)
cmod4 = mod(cmod10, 4)
! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad
exp2pi = (cocos - cmod10)/10
! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW
phiccw = (mod(cmod10, 2)== 1)
! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip
psiincr = (cmod4==1 .or. cmod4==2)
! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip
qpos = (cmod10<3 .or. cmod10>6)
end subroutine decode_cocos end subroutine decode_cocos
subroutine eq_scal(psia,fpol,isign,bsign,factb)
subroutine eq_scal(params, data)
! Rescale the magnetic field (B) and the plasma current (I)
! and/or force their signs.
!
! Notes:
! 1. signi and signb are ignored on input if equal to 0.
! They are used to assign the direction of Bphi and Ipla BEFORE scaling.
! 2. cocos=3 assumed: CCW direction is >0
! 3. Bphi and Ipla scaled by the same factor factb to keep q unchanged
! 4. factb<0 reverses the directions of Bphi and Ipla
use const_and_precisions, only : one use const_and_precisions, only : one
use gray_params, only : equilibrium_parameters, equilibrium_data
implicit none implicit none
real(wp_), intent(inout) :: psia
integer, intent(inout) :: isign,bsign ! subroutine arguments
real(wp_), dimension(:), intent(inout) :: fpol type(equilibrium_parameters), intent(inout) :: params
real(wp_), intent(in) :: factb type(equilibrium_data), intent(inout) :: data
! local variables
real(wp_) :: last_fpol
last_fpol = data%fpol(size(data%fpol))
! isign and bsign ignored on input if equal to 0 if (params%sgni /=0) &
! used to assign the direction of Bphi and Ipla BEFORE scaling data%psia = sign(data%psia, real(-params%sgni, wp_))
! cocos=3 assumed: CCW direction is >0 if (params%sgnb /=0 .and. params%sgnb * last_fpol < 0) &
! Bphi and Ipla scaled by the same factor factb to keep q unchanged data%fpol = -data%fpol
! factb<0 reverses the directions of Bphi and Ipla
if(isign/=0) psia=sign(psia,real(-isign,wp_)) data%psia = data%psia * params%factb
if(bsign/=0 .and. bsign*fpol(size(fpol))<0) fpol=-fpol data%fpol = data%fpol * params%factb
psia=psia*factb params%sgni = int(sign(one, -data%psia))
fpol=fpol*factb params%sgnb = int(sign(one, last_fpol))
isign=int(sign(one,-psia))
bsign=int(sign(one,fpol(size(fpol))))
end subroutine eq_scal end subroutine eq_scal
subroutine set_eqspl(rv,zv,psin,psiwbrad,psinr,fpol,qpsi,sspl,ssfp, &
r0,rax,zax,rbnd,zbnd,ixp) subroutine set_eqspl(params, data)
use const_and_precisions, only : zero,one ! Computes splines for the MHD equilibrium data and stores them
use dierckx, only : regrid,coeff_parder,curfit,splev ! in their respective global variables, see the top of this file.
use gray_params, only : iequil use const_and_precisions, only : zero, one
use reflections, only : inside use gray_params, only : equilibrium_parameters, equilibrium_data
use utils, only : vmaxmin,vmaxmini use dierckx, only : regrid, coeff_parder, curfit, splev
use gray_params, only : iequil
use reflections, only : inside
use utils, only : vmaxmin, vmaxmini
implicit none implicit none
! arguments
real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol,qpsi ! subroutine arguments
real(wp_), dimension(:,:), intent(in) :: psin type(equilibrium_parameters), intent(in) :: params
real(wp_), intent(in) :: psiwbrad type(equilibrium_data), intent(in) :: data
real(wp_), intent(in) :: sspl,ssfp
real(wp_), intent(in), optional :: r0,rax,zax ! local constants
real(wp_), dimension(:), intent(in), optional :: rbnd,zbnd
integer, intent(in), optional :: ixp
! local constants
integer, parameter :: iopt=0 integer, parameter :: iopt=0
! local variables
! local variables
integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf
integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup
real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp
real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1 real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1
real(wp_), dimension(size(psinr)) :: rhotn real(wp_), dimension(size(data%psinr)) :: rhotn
real(wp_), dimension(1) :: fpoli real(wp_), dimension(1) :: fpoli
real(wp_), dimension(:), allocatable :: rv1d,zv1d,fvpsi,wf,wrk real(wp_), dimension(:), allocatable :: rv1d,zv1d,fvpsi,wf,wrk
integer, dimension(:), allocatable :: iwrk integer, dimension(:), allocatable :: iwrk
integer :: ier,ixploc,info,i,j,ij integer :: ier,ixploc,info,i,j,ij
! compute array sizes and prepare working space arrays ! compute array sizes and prepare working space arrays
nr=size(rv) nr=size(data%rv)
nz=size(zv) nz=size(data%zv)
nrz=nr*nz nrz=nr*nz
nrest=nr+ksplp nrest=nr+ksplp
nzest=nz+ksplp nzest=nz+ksplp
lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest) lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest)
liwrk=nz+nr+nrest+nzest+kspl liwrk=nz+nr+nrest+nzest+kspl
npsi=size(psinr) npsi=size(data%psinr)
npsest=npsi+ksplp npsest=npsi+ksplp
lwrkf=npsi*ksplp+npsest*(7+3*kspl) lwrkf=npsi*ksplp+npsest*(7+3*kspl)
allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest))) allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest)))
! spline fitting/interpolation of psin(i,j) and derivatives ! spline fitting/interpolation of data%psin(i,j) and derivatives
! allocate knots and spline coefficients arrays ! allocate knots and spline coefficients arrays
if (allocated(tr)) deallocate(tr) if (allocated(tr)) deallocate(tr)
if (allocated(tz)) deallocate(tz) if (allocated(tz)) deallocate(tz)
if (allocated(cceq)) deallocate(cceq) if (allocated(cceq)) deallocate(cceq)
allocate(tr(nrest),tz(nzest),cceq(nrz)) allocate(tr(nrest),tz(nzest),cceq(nrz))
! length in m !!! ! length in m !!!
rmnm=data%rv(1)
rmnm=rv(1) rmxm=data%rv(nr)
rmxm=rv(nr) zmnm=data%zv(1)
zmnm=zv(1) zmxm=data%zv(nz)
zmxm=zv(nz)
if (iequil>2) then if (iequil>2) then
! data valid only inside boundary (psin=0 outside), e.g. source==ESCO ! data valid only inside boundary (data%psin=0 outside), e.g. source==ESCO
! presence of boundary anticipated here to filter invalid data ! presence of boundary anticipated here to filter invalid data
if(present(rbnd).and.present(zbnd)) then nbnd = min(size(data%rbnd), size(data%zbnd))
nbnd=min(size(rbnd),size(zbnd))
else ! determine number of valid grid points
nbnd=0
end if
! determine number of valid grid points
nrz=0 nrz=0
do j=1,nz do j=1,nz
do i=1,nr do i=1,nr
if (nbnd.gt.0) then if (nbnd.gt.0) then
if(.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle
else else
if(psin(i,j).le.0.0d0) cycle if(data%psin(i,j).le.0.0d0) cycle
end if end if
nrz=nrz+1 nrz=nrz+1
end do end do
end do end do
! store valid data
! store valid data
allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz)) allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz))
ij=0 ij=0
do j=1,nz do j=1,nz
do i=1,nr do i=1,nr
if (nbnd.gt.0) then if (nbnd.gt.0) then
if(.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle
else else
if(psin(i,j).le.0.0d0) cycle if(data%psin(i,j).le.0.0d0) cycle
end if end if
ij=ij+1 ij=ij+1
rv1d(ij)=rv(i) rv1d(ij)=data%rv(i)
zv1d(ij)=zv(j) zv1d(ij)=data%zv(j)
fvpsi(ij)=psin(i,j) fvpsi(ij)=data%psin(i,j)
wf(ij)=1.0d0 wf(ij)=1.0d0
end do end do
end do end do
! fit as a scattered set of points ! fit as a scattered set of points
! use reduced number of knots to limit memory comsumption ? ! use reduced number of knots to limit memory comsumption ?
nsr=nr/4+4 nsr=nr/4+4
nsz=nz/4+4 nsz=nz/4+4
sspln=sspl sspln=params%ssplps
call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln, & call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln, &
rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier)
! if ier=-1 data are fitted using sspl=0 ! if ier=-1 data are fitted using params%ssplps=0
if(ier.eq.-1) then if(ier.eq.-1) then
sspln=0.0_wp_ sspln=0.0_wp_
nsr=nr/4+4 nsr=nr/4+4
@ -404,29 +447,29 @@ contains
rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier)
end if end if
deallocate(rv1d,zv1d,wf,fvpsi) deallocate(rv1d,zv1d,wf,fvpsi)
! reset nrz to the total number of grid points for next allocations ! reset nrz to the total number of grid points for next allocations
nrz=nr*nz nrz=nr*nz
else else
! iequil==2: data are valid on the full R,z grid ! iequil==2: data are valid on the full R,z grid
! reshape 2D psi array to 1D (transposed) array and compute spline coeffs ! reshape 2D psi array to 1D (transposed) array and compute spline coeffs
allocate(fvpsi(nrz)) allocate(fvpsi(nrz))
fvpsi=reshape(transpose(psin),(/nrz/)) fvpsi=reshape(transpose(data%psin),(/nrz/))
sspln=sspl sspln=params%ssplps
call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm, &
kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, &
wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
! if ier=-1 data are re-fitted using sspl=0 ! if ier=-1 data are re-fitted using params%ssplps=0
if(ier==-1) then if(ier==-1) then
sspln=0.0_wp_ sspln=0.0_wp_
call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm, &
kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, &
wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
end if end if
deallocate(fvpsi) deallocate(fvpsi)
end if end if
! compute spline coefficients for psi partial derivatives ! compute spline coefficients for psi partial derivatives
lw10 = nr*(ksplp-1) + nz*ksplp + nrz lw10 = nr*(ksplp-1) + nz*ksplp + nrz
lw01 = nr*ksplp + nz*(ksplp-1) + nrz lw01 = nr*ksplp + nz*(ksplp-1) + nrz
lw20 = nr*(ksplp-2) + nz*ksplp + nrz lw20 = nr*(ksplp-2) + nz*ksplp + nrz
@ -444,74 +487,57 @@ contains
call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier)
call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier) call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier)
! spline interpolation of fpol(i) ! Spline interpolation of data%fpol(i)
! allocate knots and spline coefficients arrays ! Allocate knots and spline coefficients arrays
if (allocated(tfp)) deallocate(tfp) if (allocated(tfp)) deallocate(tfp)
if (allocated(cfp)) deallocate(cfp) if (allocated(cfp)) deallocate(cfp)
allocate(tfp(npsest),cfp(npsest)) allocate(tfp(npsest),cfp(npsest))
allocate(wf(npsi)) allocate(wf(npsi))
wf(1:npsi-1)=one wf(1:npsi-1)=one
wf(npsi)=1.0e2_wp_ wf(npsi)=1.0e2_wp_
call curfit(iopt,npsi,psinr,fpol,wf,zero,one,kspl,ssfp,npsest,nsf, & call curfit(iopt,npsi,data%psinr,data%fpol,wf,zero,one,kspl,params%ssplf,npsest,nsf, &
tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier) tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier)
call splev(tfp,nsf,cfp,kspl,psinr(npsi:npsi),fpoli,1,ier) call splev(tfp,nsf,cfp,kspl,data%psinr(npsi:npsi),fpoli,1,ier)
! set vacuum value used outside 0<=psin<=1 range ! Set vacuum value used outside 0<=data%psin<=1 range
fpolas=fpoli(1) fpolas=fpoli(1)
sgnbphi=sign(one,fpolas) sgnbphi=sign(one,fpolas)
! free temporary arrays ! Free temporary arrays
deallocate(wrk,iwrk,wf) deallocate(wrk,iwrk,wf)
! re-normalize psi after spline computation ! Re-normalize psi after spline computation
! start with un-corrected psi ! Start with un-corrected psi
psia=data%psia
psia=psiwbrad
psinop=0.0_wp_ psinop=0.0_wp_
psiant=1.0_wp_ psiant=1.0_wp_
! use provided boundary to set an initial guess for the search of O/X points ! Use provided boundary to set an initial guess
! for the search of O/X points
nbnd=0 nbnd=min(size(data%rbnd), size(data%zbnd))
if(present(rbnd).and.present(zbnd)) then
nbnd=min(size(rbnd),size(zbnd))
end if
if (nbnd>0) then if (nbnd>0) then
call vmaxmini(zbnd,nbnd,zbinf,zbsup,ibinf,ibsup) call vmaxmini(data%zbnd,nbnd,zbinf,zbsup,ibinf,ibsup)
rbinf=rbnd(ibinf) rbinf=data%rbnd(ibinf)
rbsup=rbnd(ibsup) rbsup=data%rbnd(ibsup)
call vmaxmin(rbnd,nbnd,rbmin,rbmax) call vmaxmin(data%rbnd,nbnd,rbmin,rbmax)
else else
zbinf=zv(2) zbinf=data%zv(2)
zbsup=zv(nz-1) zbsup=data%zv(nz-1)
rbinf=rv((nr+1)/2) rbinf=data%rv((nr+1)/2)
rbsup=rbinf rbsup=rbinf
rbmin=rv(2) rbmin=data%rv(2)
rbmax=rv(nr-1) rbmax=data%rv(nr-1)
end if end if
! search for exact location of the magnetic axis ! Search for exact location of the magnetic axis
rax0=data%rax
if(present(rax)) then zax0=data%zax
rax0=rax
else
rax0=0.5_wp_*(rbmin+rbmax)
end if
if(present(zax)) then
zax0=zax
else
zax0=0.5_wp_*(zbinf+zbsup)
end if
call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info) call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info)
print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp
! search for X-point if ixp not = 0 ! search for X-point if params%ixp /= 0
if(present(ixp)) then ixploc = params%ixp
ixploc=ixp
else
ixploc=0
end if
if(ixploc/=0) then if(ixploc/=0) then
if(ixploc<0) then if(ixploc<0) then
call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info) call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info)
@ -543,11 +569,13 @@ contains
if (ixploc==0) then if (ixploc==0) then
psinop=psinoptmp psinop=psinoptmp
psiant=one-psinop psiant=one-psinop
! find upper horizontal tangent point
! Find upper horizontal tangent point
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
zbsup=z1 zbsup=z1
rbsup=r1 rbsup=r1
! find lower horizontal tangent point
! Find lower horizontal tangent point
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
zbinf=z1 zbinf=z1
rbinf=r1 rbinf=r1
@ -555,29 +583,23 @@ contains
end if end if
print*,' ' print*,' '
! save Bt value on axis (required in flux_average and used in Jcd def) ! Save Bt value on axis (required in flux_average and used in Jcd def)
! and vacuum value B0 at ref. radius R0 (used in Jcd_astra def) ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def)
call equinum_fpol(0.0_wp_,btaxis) call equinum_fpol(0.0_wp_,btaxis)
btaxis=btaxis/rmaxis btaxis = btaxis/rmaxis
if(present(r0)) then btrcen = fpolas/data%rvac
btrcen=fpolas/r0 rcen = data%rvac
rcen=r0 print '(a,f8.4)', 'BT_centr=', btrcen
else print '(a,f8.4)', 'BT_axis =', btaxis
btrcen=fpolas/rmaxis
rcen=rmaxis
end if
print'(a,f8.4)','BT_centr= ',btrcen
print'(a,f8.4)','BT_axis = ',btaxis
! compute rho_pol/rho_tor mapping based on input q profile ! Compute rho_pol/rho_tor mapping based on input q profile
call setqphi_num(psinr,abs(qpsi),abs(psia),rhotn) call setqphi_num(data%psinr,abs(data%qpsi),abs(psia),rhotn)
call set_rhospl(sqrt(psinr),rhotn) call set_rhospl(sqrt(data%psinr),rhotn)
end subroutine set_eqspl end subroutine set_eqspl
subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, & subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
tx,nknt_x,ty,nknt_y,coeff,ierr) tx,nknt_x,ty,nknt_y,coeff,ierr)
use const_and_precisions, only : wp_, comp_eps use const_and_precisions, only : wp_, comp_eps

16
src/gray_jetto1beam.f90
View File

@ -4,8 +4,8 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use units, only : ucenr,usumm,uprj0,uwbm,udisp,ubres,ucnt,uoutr,ueq,uprfin, & use units, only : ucenr,usumm,uprj0,uwbm,udisp,ubres,ucnt,uoutr,ueq,uprfin, &
uflx,upec,uprm,ubeam uflx,upec,uprm,ubeam
use gray_params, only : read_params,rtrparam_type,hcdparam_type,antctrl_type,& use gray_params, only : read_params,raytracing,ecrh_cd,antenna,&
eqparam_type,prfparam_type,outparam_type equilibrium,profiles,output
use beams, only : read_beam2 use beams, only : read_beam2
use graycore, only : gray_main use graycore, only : gray_main
use reflections, only : range2rect use reflections, only : range2rect
@ -29,12 +29,12 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
real(wp_) :: fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir real(wp_) :: fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir
real(wp_), dimension(5) :: rlim,zlim real(wp_), dimension(5) :: rlim,zlim
logical, save :: firstcall=.true. logical, save :: firstcall=.true.
type(rtrparam_type) :: rtrp type(raytracing) :: rtrp
type(hcdparam_type) :: hcdp type(ecrh_cd) :: hcdp
type(antctrl_type) :: antp type(antenna) :: antp
type(eqparam_type) :: eqp type(equilibrium) :: eqp
type(prfparam_type) :: prfp type(profiles) :: prfp
type(outparam_type) :: outp type(output) :: outp
! if first call read parameters from external file ! if first call read parameters from external file
if (firstcall) then if (firstcall) then

395
src/gray_params.f90
View File

@ -1,181 +1,246 @@
module gray_params module gray_params
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
implicit none implicit none
integer, parameter :: lenfnm=256 integer, parameter :: lenfnm = 256
integer, parameter :: headw=132,headl=21 integer, parameter :: headw = 132, headl = 21
type antctrl_type ! Antenna/wave launcher parameters
real(wp_) :: alpha, beta, power type antenna_parameters
real(wp_) :: psi, chi ! From gray_params.data:
integer :: iox real(wp_) :: alpha, beta ! Launching angles
integer :: ibeam real(wp_) :: power ! Initial power
character(len=lenfnm) :: filenm real(wp_) :: psi, chi ! Initial polarisation angles
end type antctrl_type integer :: iox ! Initial wave mode
integer :: ibeam ! Beam kind
character(len=lenfnm) :: filenm ! beamdata.txt filename
type eqparam_type ! From beamdata.txt:
real(wp_) :: fghz ! EC frequency
real(wp_) :: pos(3) ! Launcher position (tokamak frame)
real(wp_) :: w(2) ! Beam waists w(z) for the amplitude (local frame)
real(wp_) :: ri(2) ! Beam inverse radii 1/R(z) for the phase (local frame)
real(wp_) :: phi(2) ! Axes orientation angles for amplitude, phase (local frame)
end type
! MHD equilibrium parameters
type equilibrium_parameters
real(wp_) :: ssplps, ssplf, factb real(wp_) :: ssplps, ssplf, factb
integer :: sgnb, sgni, ixp integer :: sgnb, sgni, ixp
integer :: iequil, icocos, ipsinorm, idesc, ifreefmt integer :: iequil, icocos, ipsinorm, idesc, ifreefmt
character(len=lenfnm) :: filenm character(len=lenfnm) :: filenm
end type eqparam_type end type
type prfparam_type ! Kinetic plasma profiles
real(wp_) :: psnbnd, sspld, factne, factte type profiles_parameters
integer :: iscal, irho !, icrho, icte, icne, iczf real(wp_) :: psnbnd ! plasma boundary (ψ: ne(ψ)=0)
real(wp_) :: sspld, factne, factte
integer :: iscal, irho
integer :: iprof integer :: iprof
character(len=lenfnm) :: filenm character(len=lenfnm) :: filenm
end type prfparam_type end type
type rtrparam_type ! Raytracing parameters
type raytracing_parameters
real(wp_) :: rwmax, dst real(wp_) :: rwmax, dst
integer :: nrayr, nrayth, nstep integer :: nrayr, nrayth, nstep
integer :: igrad, idst, ipass, ipol integer :: igrad, idst, ipass, ipol
end type rtrparam_type end type
type hcdparam_type ! EC resonant heating & Current Drive parameters
type ecrh_cd_parameters
integer :: iwarm, ilarm, imx, ieccd integer :: iwarm, ilarm, imx, ieccd
end type hcdparam_type end type
type outparam_type ! Output data parameters
integer :: ipec, nrho, istpr, istpl type output_parameters
end type outparam_type integer :: ipec ! Grid spacing for profiles (profili EC)
integer :: nrho ! Number of grid steps for EC profiles + 1
integer :: istpr ! Subsampling factor for beam cross-section (units 8, 12)
integer :: istpl ! Subsampling factor for outer rays (unit 33)
end type
integer, save :: iequil,iprof,ipol ! Other parameters
integer, save :: iwarm,ilarm,imx,ieccd type misc_parameters
integer, save :: igrad,idst,ipass real(wp_) :: rwall ! R of the limiter (fallback)
integer, save :: istpr0,istpl0 end type
integer, save :: ipec,nnd
! MHD equilibrium data
type equilibrium_data
real(wp_), allocatable :: rv(:) ! R of the uniform grid
real(wp_), allocatable :: zv(:) ! Z of the uniform grid
real(wp_), allocatable :: rlim(:) ! R of the limiter contour (wall)
real(wp_), allocatable :: zlim(:) ! Z of the limiter contour
real(wp_), allocatable :: rbnd(:) ! R of the boundary contour (plasma)
real(wp_), allocatable :: zbnd(:) ! Z of the boundary contour
real(wp_), allocatable :: fpol(:) ! Poloidal current function
real(wp_), allocatable :: qpsi(:) ! Safety factor on the flux grid
real(wp_), allocatable :: psin(:,:) ! Poloidal flux on a uniform grid
real(wp_), allocatable :: psinr(:) ! Poloidal flux
real(wp_) :: psia ! Poloidal flux at edge - flux at magnetic axis
real(wp_) :: rvac ! Reference R₀ (B = B₀R₀ without the plasma)
real(wp_) :: rax ! R of the magnetic axis
real(wp_) :: zax ! Z of the magnetic axis
end type
! Kinetic plasma profiles data
type profiles_data
real(wp_), allocatable :: psrad(:) ! Radial coordinate
real(wp_), allocatable :: terad(:) ! Electron temperature profile
real(wp_), allocatable :: derad(:) ! Electron density profile
real(wp_), allocatable :: zfc(:) ! Effective charge profile
end type
! All GRAY parameters
type gray_parameters
type(antenna_parameters) :: antenna
type(equilibrium_parameters) :: equilibrium
type(profiles_parameters) :: profiles
type(raytracing_parameters) :: raytracing
type(ecrh_cd_parameters) :: ecrh_cd
type(output_parameters) :: output
type(misc_parameters) :: misc
end type
! All GRAY input data
type gray_data
type(equilibrium_data) :: equilibrium
type(profiles_data) :: profiles
end type
! GRAY final results (only some)
type gray_results
real(wp_) :: pabs ! Total absorbed power
real(wp_) :: icd ! Total driven current
real(wp_), allocatable :: dpdv(:) ! Absorbed power density
real(wp_), allocatable :: jcd(:) ! Driven current density
end type
! Global variables exposed for graycore
integer, save :: iequil, iprof, ipol
integer, save :: iwarm, ilarm, imx, ieccd
integer, save :: igrad, idst, ipass
integer, save :: istpr0, istpl0
integer, save :: ipec, nnd
contains contains
subroutine print_params(rtrparam,hcdparam,antctrl,eqparam,rwall, & subroutine print_parameters(params, strout)
prfparam,outparam,strout)
implicit none implicit none
! arguments ! subroutine arguments
type(rtrparam_type), intent(in) :: rtrparam type(gray_parameters), intent(in) :: params
type(hcdparam_type), intent(in) :: hcdparam character(len=*), dimension(:), intent(out) :: strout ! min len=110, dimension(21)
type(antctrl_type), intent(in) :: antctrl
type(eqparam_type), intent(in) :: eqparam
real(wp_), intent(in) :: rwall
type(prfparam_type), intent(in) :: prfparam
type(outparam_type), intent(in) :: outparam
character(len=*), dimension(:), intent(out) :: strout ! min len=110, dimension(21)
! local variables ! local variables
character(len=8) :: rdat character(len=8) :: date
character(len=10) :: rtim character(len=10) :: time
#ifndef REVISION #ifndef REVISION
character(len=*), parameter :: REVISION="unknown" character(len=*), parameter :: REVISION="unknown"
#endif #endif
! date and time ! Date and time
call date_and_time(rdat,rtim) call date_and_time(date, time)
write(strout(1),'("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') & write(strout(1), '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') &
rdat(1:4),rdat(5:6),rdat(7:8),rtim(1:2),rtim(3:4),rtim(5:10) date(1:4), date(5:6), date(7:8), &
time(1:2), time(3:4), time(5:10)
! Git revision ! Git revision
write(strout(2),'("# GRAY Git revision: ",a)') REVISION write(strout(2), '("# GRAY Git revision: ",a)') REVISION
! equilibrium input data ! MHD equilibrium parameters
if (eqparam%iequil > 0) then if (params%equilibrium%iequil > 0) then
write(strout(3),'("# EQL input: ",a)') trim(eqparam%filenm) write(strout(3), '("# EQL input: ",a)') trim(params%equilibrium%filenm)
!!! missing values ! TODO add missing values
write(strout(7),'("# EQL B0 R0 aminor Rax zax:",5(1x,e12.5))') & write(strout(7), '("# EQL B0 R0 aminor Rax zax:",5(1x,e12.5))') &
0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_ 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_
else else
write(strout(3),'("# EQL input: N/A (vacuum)")') write(strout(3), '("# EQL input: N/A (vacuum)")')
write(strout(7),'("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")') write(strout(7), '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
end if end if
write(strout(4),'("# EQL iequil sgnb sgni factb:",3(1x,i4),1x,e12.5)') & write(strout(4), '("# EQL iequil sgnb sgni factb:",3(1x,i4),1x,e12.5)') &
eqparam%iequil, eqparam%sgnb, eqparam%sgni, eqparam%factb params%equilibrium%iequil, params%equilibrium%sgnb, params%equilibrium%sgni, params%equilibrium%factb
if (eqparam%iequil > 1) then if (params%equilibrium%iequil > 1) then
write(strout(5),'("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,i4))') & write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,i4))') &
eqparam%icocos, eqparam%ipsinorm, eqparam%idesc, eqparam%ifreefmt params%equilibrium%icocos, params%equilibrium%ipsinorm, params%equilibrium%idesc, params%equilibrium%ifreefmt
write(strout(6),'("# EQL ssplps ssplf ixp:",2(1x,e12.5),1x,i4)') & write(strout(6), '("# EQL ssplps ssplf ixp:",2(1x,e12.5),1x,i4)') &
eqparam%ssplps, eqparam%ssplf, eqparam%ixp params%equilibrium%ssplps, params%equilibrium%ssplf, params%equilibrium%ixp
else else
write(strout(5),'("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")') write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
write(strout(6),'("# EQL ssplps ssplf ixp: N/A (analytical)")') write(strout(6), '("# EQL ssplps ssplf ixp: N/A (analytical)")')
end if end if
! profiles input data ! Profiles parameters
if (eqparam%iequil > 0) then if (params%equilibrium%iequil > 0) then
write(strout(8),'("# PRF input: ",a)') trim(prfparam%filenm) write(strout(8), '("# PRF input: ",a)') trim(params%profiles%filenm)
write(strout(9),'("# PRF iprof iscal factne factte:",2(1x,i4),2(1x,e12.5))') & write(strout(9), '("# PRF iprof iscal factne factte:",2(1x,i4),2(1x,e12.5))') &
prfparam%iprof,prfparam%iscal,prfparam%factne,prfparam%factte params%profiles%iprof, params%profiles%iscal, params%profiles%factne, params%profiles%factte
if (prfparam%iprof > 0) then if (params%profiles%iprof > 0) then
write(strout(10),'("# PRF irho psnbnd sspld:",1x,i4,2(1x,e12.5))') & write(strout(10), '("# PRF irho psnbnd sspld:",1x,i4,2(1x,e12.5))') &
prfparam%irho,prfparam%psnbnd,prfparam%sspld params%profiles%irho,params%profiles%psnbnd,params%profiles%sspld
else else
write(strout(10),'("# PRF irho psnbnd sspld: N/A (analytical)")') write(strout(10), '("# PRF irho psnbnd sspld: N/A (analytical)")')
end if end if
!!! missing values ! TODO: add missing values
write(strout(11),'("# PRF Te0 ne0 Zeff0:",3(1x,e12.5))') & write(strout(11), '("# PRF Te0 ne0 Zeff0:",3(1x,e12.5))') &
0._wp_, 0._wp_, 0._wp_ 0._wp_, 0._wp_, 0._wp_
else else
write(strout(8),'("# PRF input: N/A (vacuum)")') write(strout(8), '("# PRF input: N/A (vacuum)")')
write(strout(9),'("# PRF iprof iscal factne factte: N/A (vacuum)")') write(strout(9), '("# PRF iprof iscal factne factte: N/A (vacuum)")')
write(strout(10),'("# PRF irho psnbnd sspld: N/A (vacuum)")') write(strout(10), '("# PRF irho psnbnd sspld: N/A (vacuum)")')
write(strout(11),'("# PRF Te0 ne0 Zeff0: N/A (vacuum)")') write(strout(11), '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")')
end if end if
! launch parameters ! Antenna parameters
write(strout(12),'("# ANT input: ",a)') trim(antctrl%filenm) write(strout(12), '("# ANT input: ",a)') trim(params%antenna%filenm)
write(strout(13),'("# ANT ibeam iox psi chi:",2(1x,i4),2(1x,e12.5))') & write(strout(13), '("# ANT ibeam iox psi chi:",2(1x,i4),2(1x,e12.5))') &
antctrl%ibeam, antctrl%iox, antctrl%psi, antctrl%chi params%antenna%ibeam, params%antenna%iox, params%antenna%psi, params%antenna%chi
write(strout(14),'("# ANT alpha beta power:",3(1x,e12.5))') & write(strout(14), '("# ANT alpha beta power:",3(1x,e12.5))') &
antctrl%alpha, antctrl%beta, antctrl%power params%antenna%alpha, params%antenna%beta, params%antenna%power
!!! missing values ! TODO: add missing values
write(strout(15),'("# ANT x0 y0 z0:",3(1x,e12.5))') & write(strout(15), '("# ANT x0 y0 z0:",3(1x,e12.5))') &
0._wp_, 0._wp_, 0._wp_ 0._wp_, 0._wp_, 0._wp_
!!! missing values ! TODO: add missing values
write(strout(16),'("# ANT wx wy Rcix Rciy psiw psir:",6(1x,e12.5))') & write(strout(16), '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,e12.5))') &
0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_ 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_
! wall parameters ! Other parameters
write(strout(17),'("# RFL rwall:",1x,e12.5)') rwall write(strout(17), '("# RFL rwall:",1x,e12.5)') params%misc%rwall
! code parameters ! code parameters
write(strout(18),'("# COD igrad idst ipass ipol:",4(1x,i4))') & write(strout(18), '("# COD igrad idst ipass ipol:",4(1x,i4))') &
rtrparam%igrad, rtrparam%idst, rtrparam%ipass, rtrparam%ipol params%raytracing%igrad, params%raytracing%idst, params%raytracing%ipass, params%raytracing%ipol
write(strout(19),'("# COD nrayr nrayth nstep rwmax dst:",3(1x,i4),2(1x,e12.5))') & write(strout(19), '("# COD nrayr nrayth nstep rwmax dst:",3(1x,i4),2(1x,e12.5))') &
rtrparam%nrayr, rtrparam%nrayth, rtrparam%nstep, rtrparam%rwmax, rtrparam%dst params%raytracing%nrayr, params%raytracing%nrayth, params%raytracing%nstep, params%raytracing%rwmax, params%raytracing%dst
write(strout(20),'("# COD iwarm ilarm imx ieccd:",4(1x,i4))') & write(strout(20), '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') &
hcdparam%iwarm, hcdparam%ilarm, hcdparam%imx, hcdparam%ieccd params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, params%ecrh_cd%imx, params%ecrh_cd%ieccd
write(strout(21),'("# COD ipec nrho istpr istpl:",4(1x,i4))') & write(strout(21), '("# COD ipec nrho istpr istpl:",4(1x,i4))') &
outparam%ipec, outparam%nrho, outparam%istpr, outparam%istpl params%output%ipec, params%output%nrho, params%output%istpr, params%output%istpl
end subroutine print_params end subroutine print_parameters
subroutine read_params(filenm,rtrparam,hcdparam,antctrl,eqparam,rwall, & subroutine read_parameters(filename, params, unit)
prfparam,outparam,unit)
use utils, only : get_free_unit use utils, only : get_free_unit
implicit none implicit none
! arguments ! subrouting arguments
character(len=*), intent(in) :: filenm character(len=*), intent(in) :: filename
type(rtrparam_type), intent(out) :: rtrparam type(gray_parameters), intent(out) :: params
type(hcdparam_type), intent(out) :: hcdparam
type(antctrl_type), intent(out) :: antctrl
type(eqparam_type), intent(out) :: eqparam
real(wp_), intent(out) :: rwall
type(prfparam_type), intent(out) :: prfparam
type(outparam_type), intent(out) :: outparam
integer, intent(in), optional :: unit integer, intent(in), optional :: unit
! local variables ! local variables
integer :: u, iostat integer :: u, iostat
if (present(unit)) then if (present(unit)) then
u=unit u = unit
else else
u = get_free_unit() u = get_free_unit()
end if end if
open(u,file=filenm,status='old',action='read',iostat=iostat)
if (iostat>0) then open(u, file=filename, status='old', action='read', iostat=iostat)
print'(a)', 'gray_params file not found!' if (iostat > 0) then
print '(a)', 'gray_params file not found!'
call EXIT(1) call EXIT(1)
end if end if
@ -185,18 +250,18 @@ contains
! nrayth :number of rays in angular direction ! nrayth :number of rays in angular direction
! rwmax :normalized maximum radius of beam power ! rwmax :normalized maximum radius of beam power
! rwmax=1 -> :last ray at radius = waist ! rwmax=1 -> :last ray at radius = waist
read(u,*) rtrparam%nrayr, rtrparam%nrayth, rtrparam%rwmax read(u, *) params%raytracing%nrayr, params%raytracing%nrayth, params%raytracing%rwmax
! igrad=0 :optical ray-tracing, initial conditions as for beam ! igrad=0 :optical ray-tracing, initial conditions as for beam
! igrad=1 :quasi-optical ray-tracing ! igrad=1 :quasi-optical ray-tracing
! igrad=-1 :ray-tracing, init. condit. ! igrad=-1 :ray-tracing, init. condit.
! from center of mirror and with angular spread ! from center of mirror and with angular spread
! ipass=1/2 :1 or 2 passes into plasma ! ipass=1/2 :1 or 2 passes into plasma
! ipol=0 :compute mode polarization at antenna, ipol=1 use polariz angles ! ipol=0 :compute mode polarization at antenna, ipol=1 use polariz angles
read(u,*) rtrparam%igrad, rtrparam%ipass, rtrparam%ipol read(u, *) params%raytracing%igrad, params%raytracing%ipass, params%raytracing%ipol
! dst :integration step ! dst :integration step
! nstep :maximum number of integration steps ! nstep :maximum number of integration steps
! idst=0/1/2 :0 integration in s, 1 integr. in ct, 2 integr. in Sr ! idst=0/1/2 :0 integration in s, 1 integr. in ct, 2 integr. in Sr
read(u,*) rtrparam%dst, rtrparam%nstep, rtrparam%idst read(u, *) params%raytracing%dst, params%raytracing%nstep, params%raytracing%idst
! Heating & Current drive ! Heating & Current drive
! ======================================================================== ! ========================================================================
@ -207,104 +272,102 @@ contains
! ilarm :order of larmor expansion ! ilarm :order of larmor expansion
! imx :max n of iterations in dispersion, imx<0 uses 1st ! imx :max n of iterations in dispersion, imx<0 uses 1st
! iteration in case of failure after |imx| iterations ! iteration in case of failure after |imx| iterations
read(u,*) hcdparam%iwarm,hcdparam%ilarm,hcdparam%imx read(u, *) params%ecrh_cd%iwarm,params%ecrh_cd%ilarm,params%ecrh_cd%imx
! ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models ! ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models
read(u,*) hcdparam%ieccd read(u, *) params%ecrh_cd%ieccd
! Wave launcher ! Wave launcher
! ======================================================================== ! ========================================================================
! alpha0, beta0 (cartesian) launching angles ! alpha0, beta0 (cartesian) launching angles
read(u,*) antctrl%alpha, antctrl%beta read(u, *) params%antenna%alpha, params%antenna%beta
! p0mw injected power (MW) ! p0mw injected power (MW)
read(u,*) antctrl%power read(u, *) params%antenna%power
! abs(iox)=1/2 OM/XM ! abs(iox)=1/2 OM/XM
! psipol0,chipol0 polarization angles at the antenna (if iox<0) ! psipol0,chipol0 polarization angles at the antenna (if iox<0)
read(u,*) antctrl%iox, antctrl%psi, antctrl%chi read(u, *) params%antenna%iox, params%antenna%psi, params%antenna%chi
! ibeam=0 :read data for beam as above read(u, *) params%antenna%ibeam
! ibeam=1 :read data from file simple astigmatic beam read(u, *) params%antenna%filenm
! ibeam=2 :read data from file general astigmatic beam
read(u,*) antctrl%ibeam
read(u,*) antctrl%filenm
! Magnetic equilibrium ! MHD equilibrium
! ======================================================================== ! ========================================================================
! iequil=0 :vacuum ! iequil=0 :vacuum (no plasma)
! iequil=1 :analytical equilibrium ! iequil=1 :analytical equilibrium
! iequil=2 :read eqdsk ! iequil=2 :read eqdsk
read(u,*) eqparam%iequil ! iequil=2 :read eqdsk, data only valid inside last closed flux surface
read(u,*) eqparam%filenm read(u, *) params%equilibrium%iequil
read(u, *) params%equilibrium%filenm
! icocos :index for equilibrium from COCOS - O. Sauter Feb 2012 ! icocos :index for equilibrium from COCOS - O. Sauter Feb 2012
! ipsinorm :0 standard EQDSK format, 1 format Portone summer 2004 ! ipsinorm :0 standard EQDSK format, 1 format Portone summer 2004
read(u,*) eqparam%icocos, eqparam%ipsinorm, eqparam%idesc, eqparam%ifreefmt read(u, *) params%equilibrium%icocos, params%equilibrium%ipsinorm, params%equilibrium%idesc, params%equilibrium%ifreefmt
! ixp=0,-1,+1 : no X point , bottom/up X point ! ixp=0,-1,+1 : no X point , bottom/up X point
! ssplps : spline parameter for psi interpolation ! ssplps : spline parameter for psi interpolation
read(u,*) eqparam%ixp, eqparam%ssplps !, eqparam%ssplf read(u, *) params%equilibrium%ixp, params%equilibrium%ssplps !, params%equilibrium%ssplf
eqparam%ssplf=0.01_wp_ params%equilibrium%ssplf=0.01_wp_
! signum of toroidal B and I ! signum of toroidal B and I
! factb factor for magnetic field (only for numerical equil) ! factb factor for magnetic field (only for numerical equil)
! scaling adopted: beta=const, qpsi=const, nustar=const ! scaling adopted: beta=const, qpsi=const, nustar=const
read(u,*) eqparam%sgnb, eqparam%sgni, eqparam%factb read(u, *) params%equilibrium%sgnb, params%equilibrium%sgni, params%equilibrium%factb
! Wall ! Wall
! ======================================================================== ! ========================================================================
read(u,*) rwall read(u, *) params%misc%rwall
! Profiles ! Profiles
! ======================================================================== ! ========================================================================
! iprof=0 :analytical density and temp. profiles ! iprof=0 :analytical density and temp. profiles
! iprof>0 :numerical density and temp. profiles ! iprof>0 :numerical density and temp. profiles
read(u,*) prfparam%iprof, prfparam%irho ! irho=0,1,2 -> num profiles vs rhot,rhop,psin read(u, *) params%profiles%iprof, params%profiles%irho ! irho=0,1,2 -> num profiles vs rhot,rhop,psin
read(u,*) prfparam%filenm read(u, *) params%profiles%filenm
! psbnd value of psi ( > 1 ) of density boundary ! psbnd value of psi ( > 1 ) of density boundary
read(u,*) prfparam%psnbnd, prfparam%sspld read(u, *) params%profiles%psnbnd, params%profiles%sspld
! prfparam%sspld=0.001_wp_ ! prfparam%sspld=0.001_wp_
! iscal :ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling ! iscal :ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling
! factT factn :factor for Te&ne scaling ! factT factn :factor for Te&ne scaling
read(u,*) prfparam%factte, prfparam%factne, prfparam%iscal read(u, *) params%profiles%factte, params%profiles%factne, params%profiles%iscal
! Output ! Output
! ======================================================================== ! ========================================================================
! ipec=0/1 :pec profiles grid in psi/rhop ! ipec=0/1 :pec profiles grid in psi/rhop
! nrho :number of grid steps for pec profiles +1 ! nrho :number of grid steps for pec profiles +1
read(u,*) outparam%ipec, outparam%nrho read(u, *) params%output%ipec, params%output%nrho
! istpr0 :projection step = dsdt*istprj ! istpr0 :projection step = dsdt*istprj
! istpl0 :plot step = dsdt*istpl ! istpl0 :plot step = dsdt*istpl
read(u,*) outparam%istpr, outparam%istpl read(u, *) params%output%istpr, params%output%istpl
close(u) close(u)
end subroutine read_params end subroutine read_parameters
subroutine set_codepar(eqparam,prfparam,outparam,rtrparam,hcdparam) subroutine set_globals(params)
implicit none implicit none
type(eqparam_type), intent(in) :: eqparam
type(prfparam_type), intent(in) :: prfparam
type(outparam_type), intent(in) :: outparam
type(rtrparam_type), intent(in) :: rtrparam
type(hcdparam_type), intent(in) :: hcdparam
iequil=eqparam%iequil ! subroutine arguments
iprof=prfparam%iprof type(gray_parameters), intent(in) :: params
ipec=outparam%ipec iequil = params%equilibrium%iequil
nnd=outparam%nrho iprof = params%profiles%iprof
istpr0=outparam%istpr
istpl0=outparam%istpl
ipol=rtrparam%ipol ipec = params%output%ipec
igrad=rtrparam%igrad nnd = params%output%nrho
idst=rtrparam%idst istpr0 = params%output%istpr
ipass=rtrparam%ipass istpl0 = params%output%istpl
if (rtrparam%nrayr<5) then
igrad=0 ipol = params%raytracing%ipol
print*,' nrayr < 5 ! => OPTICAL CASE ONLY' igrad = params%raytracing%igrad
print*,' ' idst = params%raytracing%idst
ipass = params%raytracing%ipass
if (params%raytracing%nrayr < 5) then
igrad = 0
print *, ' nrayr < 5 ! => OPTICAL CASE ONLY'
print *, ' '
end if end if
iwarm=hcdparam%iwarm iwarm = params%ecrh_cd%iwarm
ilarm=hcdparam%ilarm ilarm = params%ecrh_cd%ilarm
imx=hcdparam%imx imx = params%ecrh_cd%imx
ieccd=hcdparam%ieccd ieccd = params%ecrh_cd%ieccd
end subroutine set_codepar end subroutine set_globals
end module gray_params end module gray_params

542
src/graycore.f90
View File

@ -4,56 +4,41 @@ module graycore
contains contains
subroutine gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, & subroutine gray_main(params, data, results, error, rhout)
eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, &
p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, &
psipol0,chipol0, dpdv,jcd,pabs,icd, outp,rtrp,hcdp,ierr, rhout)
use const_and_precisions, only : zero, one, degree, comp_tiny use const_and_precisions, only : zero, one, degree, comp_tiny
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl
use dispersion, only : expinit use dispersion, only : expinit
use gray_params, only : eqparam_type, prfparam_type, outparam_type, & use gray_params, only : gray_parameters, gray_data, gray_results, print_parameters, &
rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params, & iwarm, ipec, istpr0, igrad, headw, headl, ipass
iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl, ipass use beams, only : xgygcoeff, launchangles2n
use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff use beamdata, only : pweight, rayi2jk
use beamdata, only : pweight, rayi2jk use equilibrium, only : unset_eqspl, unset_rhospl, unset_q
use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, & use errcodes, only : check_err, print_errn, print_errhcd
zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q, psnbd
use errcodes, only : check_err, print_errn, print_errhcd
use magsurf_data, only : flux_average, dealloc_surfvec use magsurf_data, only : flux_average, dealloc_surfvec
use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab rhop_tab, rhot_tab
use limiter, only : set_lim, unset_lim use limiter, only : limiter_unset_globals=>unset_globals
use utils, only : vmaxmin use utils, only : vmaxmin
use reflections, only : inside use reflections, only : inside
use multipass, only : alloc_multipass, dealloc_multipass, initbeam, & use multipass, only : alloc_multipass, dealloc_multipass, initbeam, &
initmultipass, turnoffray, plasma_in, plasma_out, wall_out initmultipass, turnoffray, plasma_in, plasma_out, wall_out
use units, only : ucenr use units, only : ucenr
implicit none implicit none
! arguments
type(eqparam_type), intent(in) :: eqp
type(prfparam_type), intent(in) :: prfp
type(outparam_type), intent(in) :: outp
type(rtrparam_type), intent(in) :: rtrp
type(hcdparam_type), intent(in) :: hcdp
real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc ! Subroutine arguments
real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi type(gray_parameters), intent(in) :: params
real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim type(gray_data), intent(in) :: data
real(wp_), dimension(:,:), intent(in) :: psin type(gray_results), intent(out) :: results
real(wp_), intent(in) :: psia, rvac, rax, zax
integer, intent(in) :: iox0
real(wp_), intent(in) :: p0, fghz, psipol0, chipol0
real(wp_), intent(in) :: alpha0,beta0, w1,w2, ri1,ri2, phiw,phir
real(wp_), dimension(3), intent(in) :: xv0
real(wp_), intent(out) :: pabs,icd ! Predefined grid for the output profiles (optional)
real(wp_), dimension(:), intent(out) :: dpdv,jcd
real(wp_), dimension(:), intent(in), optional :: rhout real(wp_), dimension(:), intent(in), optional :: rhout
integer, intent(out) :: ierr ! Exit code
! local variables integer, intent(out) :: error
! local variables
real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr)
character, dimension(2), parameter :: mode=(/'O','X'/) character, dimension(2), parameter :: mode=(/'O','X'/)
@ -66,11 +51,17 @@ contains
real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl,cpl0 real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl,cpl0
real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg
real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null()
! Ray variables
real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null()
real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null() real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null()
integer :: i,j,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt
! i: integration step, jk: global ray index
integer :: i, jk
integer :: iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt
integer :: ip,ib,iopmin,ipar,iO integer :: ip,ib,iopmin,ipar,iO
integer :: igrad_b,ipol,istop_pass,nbeam_pass,nlim integer :: igrad_b,istop_pass,nbeam_pass,nlim
logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff
real(wp_), dimension(:,:,:), pointer :: yynext=>null(),yypnext=>null() real(wp_), dimension(:,:,:), pointer :: yynext=>null(),yypnext=>null()
@ -88,101 +79,92 @@ contains
logical, dimension(:), pointer :: iwait=>null() logical, dimension(:), pointer :: iwait=>null()
logical, dimension(:,:), pointer :: iroff=>null() logical, dimension(:,:), pointer :: iroff=>null()
! parameters log in file headers ! parameters log in file headers
character(len=headw), dimension(headl) :: strheader character(len=headw), dimension(headl) :: strheader
type(antctrl_type) :: antp
real(wp_) :: rwall
! ======== set environment BEGIN ======== ! ======== set environment BEGIN ========
call set_codepar(eqp,prfp,outp,rtrp,hcdp) ! Number of limiter contourn points
nlim = size(data%equilibrium%zlim)
call set_lim(rlim,zlim) ! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1)
nlim = size(zlim) call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
if(iequil<2) then ! Compute the initial cartesian wavevector (anv0)
call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) ! from launch angles α,β and the position x₀:
else ! NR(α, β, x₀)
call set_eqspl(rv,zv,psin, psia, psinr,fpol, qpsi, eqp%ssplps,eqp%ssplf, & ! Nφ(α, β, x₀)
rvac, rax,zax, rbnd,zbnd, eqp%ixp) ! Nz(α, β, x₀)
! qpsi used for rho_pol/rho_tor mapping (initializes fq,frhotor,frhopol) call launchangles2n(params%antenna, anv0)
end if
! compute flux surface averaged quantities
call flux_average ! requires frhotor for dadrhot,dvdrhot
if(iprof==0) then ! Initialise the ray variables (beamtracing)
call set_prfan(terad,derad,zfc) call init_btr(params%raytracing, yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, &
else tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv)
call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd)
end if
call xgygcoeff(fghz,ak0,bres,xgcn)
call launchangles2n(alpha0,beta0,xv0,anv0)
call init_btr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
! Initialise the dispersion module
if(iwarm > 1) call expinit if(iwarm > 1) call expinit
! Initialise the magsurf_data module
call flux_average ! requires frhotor for dadrhot,dvdrhot
call pec_init(ipec,rhout) ! Initialise the output profiles
nnd=size(rhop_tab) call pec_init(ipec, rhout)
call alloc_multipass(nnd,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & nnd = size(rhop_tab) ! number of radial profile points
stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, &
pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv)
! ========= set environment END ========= call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, &
stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, &
pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv)
! ======== pre-proc prints BEGIN ======== ! Allocate memory for the results...
antp%alpha=alpha0 allocate(results%dpdv(params%output%nrho))
antp%beta=beta0 allocate(results%jcd(params%output%nrho))
antp%power=p0
antp%psi=psipol0
antp%chi=chipol0
antp%iox=iox0
!!!!! missing values
antp%ibeam=0
antp%filenm=''
rwall=0._wp_
psnbd=prfp%psnbnd ! plasma boundary
ipol=rtrp%ipol
pabs=zero ! gray_main output ! ...and initialise them
icd=zero results%pabs = zero
dpdv=zero results%icd = zero
jcd=zero results%dpdv = zero
results%jcd = zero
call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader) ! ========= set environment END =========
call print_headers(strheader,0)
! print psi surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout ! ======== pre-proc prints BEGIN ========
call print_surfq((/1.5_wp_,2.0_wp_/)) call print_parameters(params, strheader)
! print call print_headers(strheader)
print*,' '
print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0 ! print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout
print'(a,4f8.3)','x00, y00, z00 = ',xv0 call print_surfq([1.5_wp_, 2.0_wp_])
! print Btot=Bres
! print ne, Te, q, Jphi versus psi, rhop, rhot ! print initial position
print *, ''
print '(a,2f8.3)', 'alpha0, beta0 = ', params%antenna%alpha, params%antenna%beta
print '(a,4f8.3)', 'x00, y00, z00 = ', params%antenna%pos
! print Btot=Bres
! print ne, Te, q, Jphi versus psi, rhop, rhot
call print_bres(bres) call print_bres(bres)
call print_prof call print_prof
call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2),sin(beta0*degree)) call print_maps(bres, xgcn, &
0.01_wp_*sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2), &
sin(params%antenna%beta*degree))
! ========= pre-proc prints END =========
! ========= pre-proc prints END ========= ! =========== main loop BEGIN ===========
call initmultipass(params%raytracing%ipol, params%antenna%iox, &
! =========== main loop BEGIN =========== iroff,yynext,yypnext,yw0,ypw0, &
call initmultipass(ipol,iox0,iroff,yynext,yypnext,yw0,ypw0, & stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv)
stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv)
if(ipol.eq.0) then if(params%raytracing%ipol .eq. 0) then
if(iox0.eq.2) then ! only X mode on 1st pass if(params%antenna%iox .eq. 2) then ! only X mode on 1st pass
cpl0 = (/zero,one/) cpl0 = (/zero,one/)
else ! only O mode on 1st pass else ! only O mode on 1st pass
cpl0 = (/one,zero/) cpl0 = (/one,zero/)
end if end if
end if end if
sox=one ! mode inverted for each beam sox=one ! mode inverted for each beam
iox=2 ! start with O: sox=-1, iox=1 iox=2 ! start with O: sox=-1, iox=1
psipol=psipol0 psipol = params%antenna%psi
chipol=chipol0 chipol = params%antenna%chi
call pweight(p0,p0jk) call pweight(params%antenna%power, p0jk)
nbeam_pass=1 ! max n of beam per pass nbeam_pass=1 ! max n of beam per pass
index_rt=0 ! global beam index: 1,O 2,X 1st pass index_rt=0 ! global beam index: 1,O 2,X 1st pass
@ -219,7 +201,6 @@ contains
end if end if
call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
! call print_headers((/' '/),index_rt)
if(ip.eq.1) then ! 1st pass if(ip.eq.1) then ! 1st pass
igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass
@ -230,14 +211,19 @@ contains
lgcpl1 = zero lgcpl1 = zero
p0ray = p0jk ! * initial beam power p0ray = p0jk ! * initial beam power
call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions call ic_gb(params%antenna%pos, anv0, ak0, &
params%antenna%w(1),params%antenna%w(2), &
params%antenna%ri(1),params%antenna%ri(2), &
params%antenna%phi(1),params%antenna%phi(2), &
yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions
call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization
do jk=1,nray do jk=1,nray
zzm = yw(3,jk)*0.01_wp_ zzm = yw(3,jk)*0.01_wp_
rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_
if(inside(rlim,zlim,nlim,rrm,zzm)) then ! * start propagation in/outside vessel? if(inside(data%equilibrium%rlim, data%equilibrium%zlim, &
nlim, rrm, zzm)) then ! * start propagation in/outside vessel?
iow(jk) = 1 ! + inside iow(jk) = 1 ! + inside
else else
iow(jk) = 0 ! + outside iow(jk) = 0 ! + outside
@ -273,7 +259,7 @@ contains
! update position and grad ! update position and grad
if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
ierr = 0 error = 0
istop = 0 ! stop flag for current beam istop = 0 ! stop flag for current beam
iopmin = 10 iopmin = 10
@ -289,7 +275,7 @@ contains
ierrn,igrad_b) ierrn,igrad_b)
! update global error code and print message ! update global error code and print message
if(ierrn/=0) then if(ierrn/=0) then
ierr = ior(ierr,ierrn) error = ior(error,ierrn)
call print_errn(ierrn,i,anpl) call print_errn(ierrn,i,anpl)
end if end if
@ -297,8 +283,9 @@ contains
zzm = xv(3)*0.01_wp_ zzm = xv(3)*0.01_wp_
rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_ rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_
ins_pl = (psinv>=zero .and. psinv<psnbd) ! .and. zzm>=zbinf .and. zzm<=zbsup ! in/out plasma? ins_pl = (psinv>=zero .and. psinv<params%profiles%psnbnd) ! in/out plasma?
ins_wl = inside(rlim,zlim,nlim,rrm,zzm) ! in/out vessel? ins_wl = inside(data%equilibrium%rlim, data%equilibrium%zlim, &
nlim,rrm,zzm) ! in/out vessel?
ent_pl = (mod(iop(jk),2).eq.0 .and. ins_pl) ! enter plasma ent_pl = (mod(iop(jk),2).eq.0 .and. ins_pl) ! enter plasma
ext_pl = (mod(iop(jk),2).eq.1 .and. .not.ins_pl) ! exit plasma ext_pl = (mod(iop(jk),2).eq.1 .and. .not.ins_pl) ! exit plasma
ent_wl = (mod(iow(jk),2).eq.0 .and. ins_wl) ! enter vessel ent_wl = (mod(iow(jk),2).eq.0 .and. ins_wl) ! enter vessel
@ -309,7 +296,7 @@ contains
if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass
if(ipol.eq.0) then ! + IF single mode propagation if(params%raytracing%ipol .eq. 0) then ! + IF single mode propagation
cpl = cpl0 cpl = cpl0
p0ray(jk) = p0ray(jk)*cpl(iox) p0ray(jk) = p0ray(jk)*cpl(iox)
else if(cpl(iox).lt.etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays else if(cpl(iox).lt.etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays
@ -376,7 +363,7 @@ contains
psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, & psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, &
ierrn,igrad_b) ! * update derivatives after reflection ierrn,igrad_b) ! * update derivatives after reflection
if(ierrn/=0) then ! * update global error code and print message if(ierrn/=0) then ! * update global error code and print message
ierr = ior(ierr,ierrn) error = ior(error,ierrn)
call print_errn(ierrn,i,anpl) call print_errn(ierrn,i,anpl)
end if end if
@ -429,7 +416,7 @@ contains
call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd) sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd)
if(ierrhcd/=0) then if(ierrhcd/=0) then
ierr = ior(ierr,ierrhcd) error = ior(error,ierrhcd)
call print_errhcd(ierrhcd,i,anprre,anprim,alpha) call print_errhcd(ierrhcd,i,anprre,anprim,alpha)
end if end if
else else
@ -464,9 +451,6 @@ contains
call print_output(i,jk,stv(jk),p0ray(jk),xv,psinv, & call print_output(i,jk,stv(jk),p0ray(jk),xv,psinv, &
btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, & btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
! call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv, &
! btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
! nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
end if end if
end do end do
@ -484,7 +468,7 @@ contains
end if end if
! check for any error code and stop if necessary ! check for any error code and stop if necessary
call check_err(ierr,istop) call check_err(error,istop)
! test whether further trajectory integration is unnecessary ! test whether further trajectory integration is unnecessary
call vmaxmin(tau1+tau0+lgcpl1,nray,taumn,taumx) ! test on tau + coupling call vmaxmin(tau1+tau0+lgcpl1,nray,taumn,taumx) ! test on tau + coupling
! if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1 ! (residual power~0) or (no ray active) => stop beam ! if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1 ! (residual power~0) or (no ray active) => stop beam
@ -563,8 +547,8 @@ contains
! ============ beam loop END ============ ! ============ beam loop END ============
! ======= cumulative prints BEGIN ======= ! ======= cumulative prints BEGIN =======
pabs = pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] results%pabs = results%pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output]
icd = icd + sum(icd_pass) results%icd = results%icd + sum(icd_pass)
! print final results for pass on screen ! print final results for pass on screen
write(*,*) write(*,*)
@ -581,16 +565,11 @@ contains
! print final results on screen ! print final results on screen
write(*,*) write(*,*)
write(*,'(a)') '## Final results:' write(*,'(a)') '## Final results:'
write(*,'(a,f9.4)') '## Pabs_tot (MW) = ',pabs write(*,'(a,f9.4)') '## Pabs_tot (MW) = ', results%pabs
write(*,'(a,f9.4)') '## I_tot (kA) = ',icd*1.0e3_wp_ write(*,'(a,f9.4)') '## I_tot (kA) = ', results%icd*1.0e3_wp_
! ========== free memory BEGIN ========== ! ========== free memory BEGIN ==========
call unset_eqspl
call unset_q
call unset_rhospl
call unset_prfspl
call unset_lim
call dealloc_surfvec call dealloc_surfvec
call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, & call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, &
alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
@ -603,170 +582,6 @@ contains
subroutine sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, &
eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, &
p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, &
psipol0,chipol0, jphi,jcd,dpdv,currins,pins,pabs,icd, &
jphip,dpdvp, &
rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, &
outp,rtrp,hcdp,ierr, rhout)
use const_and_precisions, only : zero, one, degree
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl
use dispersion, only : expinit
use gray_params, only : eqparam_type, prfparam_type, outparam_type, &
rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params, &
iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl
use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff
use beamdata, only : pweight, rayi2jk
use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, &
zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q
use errcodes, only : check_err, print_errn, print_errhcd
use magsurf_data, only : flux_average, dealloc_surfvec
use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab
use limiter, only : set_lim, unset_lim
use utils, only : vmaxmin
implicit none
! arguments
type(eqparam_type), intent(in) :: eqp
type(prfparam_type), intent(in) :: prfp
type(outparam_type), intent(in) :: outp
type(rtrparam_type), intent(in) :: rtrp
type(hcdparam_type), intent(in) :: hcdp
real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc
real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi
real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim
real(wp_), dimension(:,:), intent(in) :: psin
real(wp_), intent(in) :: psia, rvac, rax, zax
integer, intent(in) :: iox0
real(wp_), intent(in) :: p0, fghz, psipol0, chipol0
real(wp_), intent(in) :: alpha0,beta0, w1,w2, ri1,ri2, phiw,phir
real(wp_), dimension(3), intent(in) :: xv0
real(wp_), intent(in) :: pabs,icd
real(wp_), dimension(:), intent(in) :: jphi,jcd,dpdv,currins,pins
real(wp_), intent(out) :: jphip,dpdvp, &
rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx
real(wp_), dimension(:), intent(in), optional :: rhout
integer, intent(out) :: ierr
! local variables
real(wp_), parameter :: taucr = 12._wp_
real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre
real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0
real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx
real(wp_) :: drhotp,drhotj,dpdvmx,jphimx
real(wp_), dimension(3) :: xv,anv0,anv,bv
real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null()
real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null()
integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1
logical :: ins_pl, somein, allout
real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null()
real(wp_), dimension(:), pointer :: tau0=>null(),alphaabs0=>null(),dids0=>null(), &
ccci0=>null()
real(wp_), dimension(:), pointer :: p0jk=>null()
complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
integer, dimension(:), pointer :: iiv=>null()
! parameters log in file headers
character(len=headw), dimension(headl) :: strheader
type(antctrl_type) :: antp
real(wp_) :: rwall
! ======== set environment BEGIN ========
call set_codepar(eqp,prfp,outp,rtrp,hcdp)
call set_lim(rlim,zlim)
if(iequil<2) then
call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3))
else
call set_eqspl(rv,zv,psin, psia, psinr,fpol, qpsi, eqp%ssplps,eqp%ssplf, &
rvac, rax,zax, rbnd,zbnd, eqp%ixp)
! qpsi used for rho_pol/rho_tor mapping (initializes fq,frhotor,frhopol)
end if
! compute flux surface averaged quantities
call flux_average ! requires frhotor for dadrhot,dvdrhot
if(iprof==0) then
call set_prfan(terad,derad,zfc)
else
call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd)
end if
call xgygcoeff(fghz,ak0,bres,xgcn)
call launchangles2n(alpha0,beta0,xv0,anv0)
call init_btr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
if(iwarm > 1) call expinit
! ======= set environment END ======
! ======= pre-proc prints BEGIN ======
antp%alpha=alpha0
antp%beta=beta0
antp%power=p0
antp%psi=psipol0
antp%chi=chipol0
antp%iox=iox0
!!!!! missing values
antp%ibeam=0
antp%filenm=''
rwall=0._wp_
call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader)
call print_headers(strheader, 0)
! print psi surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout
call print_surfq((/1.5_wp_,2.0_wp_/))
! print
print*,' '
print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0
print'(a,4f8.3)','x00, y00, z00 = ',xv0
! print Btot=Bres
! print ne, Te, q, Jphi versus psi, rhop, rhot
call print_bres(bres)
call print_prof
call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2), &
sin(beta0*degree))
! ======= pre-proc prints END ======
! ======= post-proc BEGIN ======
! compute power and current density profiles for all rays
call pec_init(ipec,rhout)
nnd=size(rhop_tab)
! print power and current density profiles
call print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt)
! compute profiles width
call postproc_profiles(pabs,icd,rhot_tab,dpdv,jphi, &
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, &
rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx)
! print 0D results
call print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, &
drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
st,psipol,chipol,index_rt,p0,zero,zero) ! cpl1=cpl2=0 for a single pass
! ======= post-proc END ======
! ======= free memory BEGIN ======
call unset_eqspl
call unset_q
call unset_rhospl
call unset_prfspl
call unset_lim
call dealloc_surfvec
call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
call dealloc_pec
! ======= free memory END ======
end subroutine sum_profiles
subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
use const_and_precisions, only : wp_, zero use const_and_precisions, only : wp_, zero
@ -1167,7 +982,7 @@ contains
subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, &
bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr,igrad) bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,error,igrad)
! Compute right-hand side terms of the ray equations (dery) ! Compute right-hand side terms of the ray equations (dery)
! used after full R-K step and grad(S_I) update ! used after full R-K step and grad(S_I) update
use errcodes, only : pnpl use errcodes, only : pnpl
@ -1181,7 +996,7 @@ contains
real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr
real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm
real(wp_), dimension(3), intent(out) :: bv real(wp_), dimension(3), intent(out) :: bv
integer, intent(out) :: ierr integer, intent(out) :: error
real(wp_), dimension(3), intent(out) :: derxg real(wp_), dimension(3), intent(out) :: derxg
integer, intent(in) :: igrad integer, intent(in) :: igrad
! local variables ! local variables
@ -1196,12 +1011,12 @@ contains
call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, &
dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad) dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad)
ierr=0 error=0
if( abs(anpl) > anplth1) then if( abs(anpl) > anplth1) then
if(abs(anpl) > anplth2) then if(abs(anpl) > anplth2) then
ierr=ibset(ierr,pnpl+1) error=ibset(error,pnpl+1)
else else
ierr=ibset(ierr,pnpl) error=ibset(error,pnpl)
end if end if
end if end if
end subroutine ywppla_upd end subroutine ywppla_upd
@ -1692,7 +1507,7 @@ contains
subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,ierr) sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,error)
use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_ use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_
use gray_params, only : iwarm,ilarm,ieccd,imx use gray_params, only : iwarm,ilarm,ieccd,imx
use coreprofiles, only : fzeff use coreprofiles, only : fzeff
@ -1707,7 +1522,7 @@ contains
real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox
real(wp_),intent(out) :: anprre,anprim,alpha,didp real(wp_),intent(out) :: anprre,anprim,alpha,didp
integer, intent(out) :: nhmin,nhmax,iokhawa integer, intent(out) :: nhmin,nhmax,iokhawa
integer, intent(out) :: ierr integer, intent(out) :: error
! local constants ! local constants
real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_
! local variables ! local variables
@ -1726,7 +1541,7 @@ contains
nhmin=0 nhmin=0
nhmax=0 nhmax=0
iokhawa=0 iokhawa=0
ierr=0 error=0
if(tekev>zero) then if(tekev>zero) then
! absorption computation ! absorption computation
@ -1734,13 +1549,13 @@ contains
call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm)
if(nhmin.gt.0) then if(nhmin.gt.0) then
lrm=max(ilarm,nhmax) lrm=max(ilarm,nhmax)
call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, & call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,error,anprre,anprim, &
iwarm,imx,ex,ey,ez) iwarm,imx,ex,ey,ez)
akim=ak0*anprim akim=ak0*anprim
ratiovgr=2.0_wp_*anpr/derdnm!*vgm ratiovgr=2.0_wp_*anpr/derdnm!*vgm
alpha=2.0_wp_*akim*ratiovgr alpha=2.0_wp_*akim*ratiovgr
if(alpha<zero) then if(alpha<zero) then
ierr=ibset(ierr,palph) error=ibset(error,palph)
return return
end if end if
@ -1773,7 +1588,7 @@ contains
call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, &
ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd) ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd)
end select end select
ierr=ierr+ierrcd error=error+ierrcd
if(associated(eccdpar)) deallocate(eccdpar) if(associated(eccdpar)) deallocate(eccdpar)
effjcdav=rbavi*effjcd effjcdav=rbavi*effjcd
didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii) didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii)
@ -1844,94 +1659,6 @@ contains
end subroutine set_pol end subroutine set_pol
! logical function inside_plasma(rrm,zzm)
! use const_and_precisions, only : wp_, zero, one
! use gray_params, only : iequil
! use equilibrium, only : equian,equinum_psi,zbinf,zbsup
! use coreprofiles, only : psdbnd
! 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 beamdata, only : dst
! use limiter, only : rlim,zlim,nlim
! 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
!
! ! 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
subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon) subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. ! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
@ -2146,13 +1873,12 @@ bb: do
subroutine print_headers(strheader,index_rt) subroutine print_headers(strheader)
use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm
implicit none implicit none
! arguments ! subroutine arguments
character(len=*), dimension(:), intent(in) :: strheader character(len=*), dimension(:), intent(in) :: strheader
integer, intent(in) :: index_rt ! local variables
! local variables
integer :: i,l integer :: i,l
l=size(strheader) l=size(strheader)

51
src/limiter.f90
View File

@ -1,33 +1,44 @@
module limiter module limiter
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
! === 1D array limiter Rlim_i, Zlim_i
integer, public, save :: nlim ! Inner wall radius
real(wp_), save :: rwallm real(wp_), save :: rwallm
real(wp_), dimension(:), allocatable, save :: rlim,zlim
! Limiter contourn
integer, public, save :: nlim
real(wp_), dimension(:), allocatable, save :: rlim, zlim
contains contains
subroutine set_lim(rv,zv) subroutine set_globals(data)
! Set global variables exposed by this module.
use gray_params, only : equilibrium_data
implicit none implicit none
real(wp_), intent(in), dimension(:) :: rv,zv
! subroutine arguments
type(equilibrium_data), intent(in) :: data
if (allocated(rlim)) deallocate(rlim) if (allocated(rlim)) deallocate(rlim)
if (allocated(zlim)) deallocate(zlim) if (allocated(zlim)) deallocate(zlim)
nlim=size(rv) nlim = size(data%rlim)
allocate(rlim(nlim),zlim(nlim)) allocate(rlim(nlim), zlim(nlim))
rlim=rv rlim = data%rlim
zlim=zv zlim = data%zlim
rwallm=minval(rlim) rwallm = minval(rlim)
end subroutine set_lim end subroutine set_globals
subroutine unset_globals
subroutine unset_lim ! Unset global variables exposed by this module.
use const_and_precisions, only : zero use const_and_precisions, only : zero
implicit none
if(allocated(rlim)) deallocate(rlim)
if(allocated(zlim)) deallocate(zlim)
nlim=0
rwallm=zero
end subroutine unset_lim
end module limiter implicit none
if(allocated(rlim)) deallocate(rlim)
if(allocated(zlim)) deallocate(zlim)
nlim = 0
rwallm = zero
end subroutine unset_globals
end module limiter

6
src/magsurf_data.f90
View File

@ -110,12 +110,14 @@ contains
use dierckx, only : regrid,coeff_parder use dierckx, only : regrid,coeff_parder
use utils, only : get_free_unit use utils, only : get_free_unit
implicit none implicit none
! local constants
! local constants
integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, & integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, &
njest=nnintp+ksp+1,nlest=nlam+ksp+1, & njest=nnintp+ksp+1,nlest=nlam+ksp+1, &
lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, & lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, &
kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam
! local variables
! local variables
integer :: ier,ierr,l,jp,inc,inc1,iopt,njp,nlm,ninpr integer :: ier,ierr,l,jp,inc,inc1,iopt,njp,nlm,ninpr
integer, dimension(kwrk) :: iwrk integer, dimension(kwrk) :: iwrk
real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, & real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, &

638
src/main.f90
View File

@ -1,215 +1,451 @@
program main_std program main
use const_and_precisions, only : wp_,one use const_and_precisions, only : wp_, one, zero
use graycore, only : gray_main, sum_profiles use graycore, only : gray_main
use gray_params, only : read_params, antctrl_type, eqparam_type, & use gray_params, only : gray_parameters, gray_data, gray_results, &
prfparam_type, outparam_type, rtrparam_type, & read_parameters, params_set_globals => set_globals
hcdparam_type
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
use reflections, only : range2rect
implicit none implicit none
type(antctrl_type) :: antp
type(eqparam_type) :: eqp
type(prfparam_type) :: prfp
type(outparam_type) :: outp
type(rtrparam_type) :: rtrp
type(hcdparam_type) :: hcdp
real(wp_), dimension(:), allocatable :: psrad, terad, derad, zfc ! gray_main subroutine arguments
real(wp_), dimension(:), allocatable :: rv, zv, psinr, fpol, qpsi type(gray_parameters) :: params ! Inputs
real(wp_), dimension(:), allocatable :: rbnd, zbnd, rlim, zlim type(gray_data) :: data !
real(wp_), dimension(:,:), allocatable :: psin type(gray_results) :: results ! Outputs
real(wp_) :: psia, rvac, rax, zax integer :: error ! Exit code
real(wp_) :: fghz
real(wp_) :: x0, y0, z0, w1, w2, ri1, ri2, phiw, phir
real(wp_) :: pec,icd logical :: sum_mode = .false.
integer :: ierr ! Load the parameters and also copy them into
real(wp_), dimension(:), allocatable :: xrad, rhot, dpdv, jcd ! global variables exported by the gray_params
real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx call read_parameters('gray_params.data', params)
call params_set_globals(params)
logical :: sum_mode = .FALSE. ! Read the input data into set the global variables
! of the respective module. Note: order matters.
call init_equilibrium(params, data)
call init_profiles(params, data)
call init_antenna(params%antenna)
call init_misc(params, data)
! ------- sum mode variables -------
real(wp_), dimension(:), allocatable :: jphi, currins, pins, rtin, rpin
integer :: i,j,k,n,ngam,irt
character(len=255) :: fn,dumstr
real(wp_), dimension(5) :: f48v
real(wp_) :: gam,alp,bet, jphip,dpdvp, &
rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx
! ----------------------------------
call read_params('gray_params.data', rtrp, hcdp, antp, eqp, rwallm, prfp, outp)
! ======= read input data BEGIN =======
!------------ equilibrium ------------
if(eqp%iequil<2) then
call read_equil_an(eqp%filenm, rtrp%ipass, rv, zv, fpol, qpsi, rlim, zlim)
! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0)
psia = sign(one,qpsi(2)*fpol(1))
else
call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac, &
rax,zax, rbnd,zbnd, rlim,zlim, &
eqp%ipsinorm,eqp%idesc,eqp%ifreefmt)
call change_cocos(psia, fpol, qpsi, eqp%icocos, 3)
end if
! re-scale B/I and/or force signs. If sgn=0 on input, set to fpol/-psia signs on output
call eq_scal(psia, fpol, eqp%sgni, eqp%sgnb, eqp%factb)
! ??? analytical only? change for numerical!
! qpsi(1) = sign(qpsi(1),qpsi(1)*qpsi(2)*psia*fpol(1))
! qpsi(2) = sign(qpsi(2),psia*fpol(1))
! ----------------------------------
!------------- profiles -------------
if(prfp%iprof == 0) then
call read_profiles_an(prfp%filenm, terad, derad, zfc)
else
call read_profiles(prfp%filenm, xrad, terad, derad, zfc)
allocate(psrad(size(xrad)))
if(prfp%irho == 0) then ! xrad==rhot
allocate(rhot(size(psinr)))
call setqphi_num(psinr,qpsi,psia,rhot)
call set_rhospl(sqrt(psinr),rhot)
deallocate(rhot)
psrad = frhopolv(xrad)**2
else if(prfp%irho == 1) then ! xrad==rhop
psrad = xrad**2
else
psrad = xrad
end if
deallocate(xrad)
end if
! re-scale input data
call tene_scal(terad, derad, prfp%factte, prfp%factne, &
eqp%factb, prfp%iscal, prfp%iprof)
! ----------------------------------
!------------- antenna --------------
! interpolate beam table if antctrl%ibeam>0
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
! ----------------------------------
!--------------- wall ---------------
! set simple limiter if not read from EQDSK
! need to clean up...
r0m=sqrt(x0**2+y0**2)*0.01_wp_
dzmx=abs(rtrp%ipass)*rtrp%dst*rtrp%nstep*0.01_wp_
z0m=z0*0.01_wp_
if (.not.allocated(rlim).or.rtrp%ipass<0) then
if (allocated(rlim)) deallocate(rlim)
if (allocated(zlim)) deallocate(zlim)
allocate(rlim(5))
allocate(zlim(5))
if (rtrp%ipass<0) rtrp%ipass = -rtrp%ipass
if(eqp%iequil<2) then
rmxm=(rv(1)+rv(2))*0.01_wp_
else
rmxm=rv(size(rv))
end if
call range2rect(rwallm,max(r0m,rmxm),z0m-dzmx,z0m+dzmx,rlim,zlim)
end if
! ----------------------------------
! ======= read input data END =======
! ========================= MAIN SUBROUTINE CALL =========================
allocate(dpdv(outp%nrho), jcd(outp%nrho))
if (sum_mode) then if (sum_mode) then
allocate(jphi(outp%nrho), currins(outp%nrho), & sum: block
pins(outp%nrho), rtin(outp%nrho), rpin(outp%nrho)) real(wp_) :: pabs, icd, pec
real(wp_), dimension(:), allocatable :: dpdv, jcd, jphi
real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin
integer :: i, j, k, n, ngam, irt
character(len=255) :: filename
real(wp_), dimension(5) :: f48v
real(wp_) :: gam,alp,bet, jphip,dpdvp, &
rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx
allocate(jphi(params%output%nrho), currins(params%output%nrho), &
pins(params%output%nrho), rtin(params%output%nrho), &
rpin(params%output%nrho))
open(100,file='filelist.txt',action='read',status='old') open(100, file='filelist.txt', action='read', status='old')
read(100,*) n,ngam read(100, *) n, ngam
do i=1,n do i=1,n
read(100,*) fn read(100, *) filename
open(100+i,file=fn,action='read',status='old') open(100 + i, file=filename, action='read', status='old')
do j=1,22 do j=1,22
read(100+i,*) dumstr read(100 + i, *)
end do
end do
close(100)
open(100+n+1,file='f48sum.txt',action='write',status='unknown')
open(100+n+2,file='f7sum.txt',action='write',status='unknown')
do k=1,ngam
jphi=0.0_wp_
jcd=0.0_wp_
dpdv=0.0_wp_
currins=0.0_wp_
pins=0.0_wp_
do j=1,outp%nrho
do i=1,n
read(100+i,*) gam,alp,bet,rpin(j),rtin(j),f48v(1:5),irt
jphi(j)=jphi(j)+f48v(1)
jcd(j)=jcd(j)+f48v(2)
dpdv(j)=dpdv(j)+f48v(3)
currins(j)=currins(j)+f48v(4)
pins(j)=pins(j)+f48v(5)
end do end do
write(100+n+1,'(10(1x,e16.8e3),i5)') gam,alp,bet,rpin(j),rtin(j), &
jphi(j),jcd(j),dpdv(j),currins(j),pins(j),irt
end do end do
pec=pins(outp%nrho) close(100)
icd=currins(outp%nrho)
write(100+n+1,*) open(100 + n+1, file='f48sum.txt', action='write', status='unknown')
call sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & open(100 + n+2, file='f7sum.txt', action='write', status='unknown')
psrad,terad,derad,zfc,prfp, rlim,zlim, &
antp%beta,fghz,antp%alpha,antp%beta, & do k=1,ngam
(/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox, & jphi = zero
antp%psi, antp%chi, jphi,jcd,dpdv,currins,pins,pec,icd, & jcd = zero
jphip,dpdvp, rhotj,rhotjava,rhotp,rhotpav, & dpdv = zero
drhotjava,drhotpav, ratjamx,ratjbmx, outp,rtrp,hcdp,ierr) currins = zero
write(100+n+2,'(15(1x,e12.5),i5,4(1x,e12.5))') gam,alp,bet,icd,pec, & pins = zero
jphip,dpdvp, & do j=1,params%output%nrho
rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx do i=1,n
end do read(100+i, *) gam, alp, bet, rpin(j), rtin(j), f48v(1:5), irt
do i=1,n+2 jphi(j) = f48v(1) + jphi(j)
close(100+i) jcd(j) = f48v(2) + jcd(j)
end do dpdv(j) = f48v(3) + dpdv(j)
currins(j) = f48v(4) + currins(j)
pins(j) = f48v(5) + pins(j)
end do
write(100 + n+1,'(10(1x,e16.8e3),i5)') &
gam, alp, bet, rpin(j), rtin(j), &
jphi(j), jcd(j), dpdv(j), currins(j), pins(j), irt
end do
pec = pins(params%output%nrho)
icd = currins(params%output%nrho)
write(100 + n+1, *)
call sum_profiles(params, jphi, jcd, dpdv, currins, &
pins, pabs, icd, jphip, dpdvp, rhotj, &
rhotjava, rhotp, rhotpav, drhotjava, &
drhotpav, ratjamx, ratjbmx)
write(100 + n+2, '(15(1x,e12.5),i5,4(1x,e12.5))') &
gam, alp, bet, icd, pabs, jphip, dpdvp, &
rhotj, rhotjava, rhotp, rhotpav, &
drhotjava, drhotpav, ratjamx, ratjbmx
end do
do i=1,n+2
close(100 + i)
end do
deallocate(dpdv, jcd, jphi, currins, pins, rtin, rpin)
end block sum
else else
call gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & call gray_main(params, data, results, error)
psrad,terad,derad,zfc,prfp, rlim,zlim, &
antp%power,fghz,antp%alpha,antp%beta, &
(/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox, &
antp%psi,antp%chi, dpdv,jcd,pec,icd, outp,rtrp,hcdp,ierr)
end if end if
! ========================================================================
! ======= control prints BEGIN ======= print '(a)'
if(ierr/=0) print*,' IERR = ', ierr print '(a,f9.4)', 'Pabs (MW)=', results%pabs
print*,' ' print '(a,f9.4)', 'Icd (kA)=', results%icd * 1.0e3_wp_
print*,'Pabs (MW) = ', pec
print*,'Icd (kA) = ', icd*1.0e3_wp_
! ======= control prints END =======
! ======= free memory BEGIN ======= ! Free memory
if(allocated(psrad)) deallocate(psrad) call deinit_equilibrium(data%equilibrium)
if(allocated(terad)) deallocate(terad, derad, zfc) call deinit_profiles(data%profiles)
if(allocated(rv)) deallocate(rv, zv, fpol, qpsi) call deinit_misc
if(allocated(psin)) deallocate(psin, psinr) deallocate(results%dpdv, results%jcd)
if(allocated(rbnd)) deallocate(rbnd, zbnd)
if(allocated(rlim)) deallocate(rlim, zlim) contains
if(allocated(dpdv)) deallocate(dpdv, jcd)
if(sum_mode) deallocate(jphi, currins, pins, rtin, rpin) subroutine init_equilibrium(params, data)
! ======= free memory END ====== ! Reads the MHD equilibrium file (either in the G-EQDSK format
end program main_std ! or an analytical description) and initialises the respective
! GRAY parameters and data.
use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, &
set_equian, set_eqspl, eq_scal
implicit none
! subroutine arguments
type(gray_parameters), intent(inout) :: params
type(gray_data), intent(out) :: data
if(params%equilibrium%iequil < 2) then
! Analytical equilibrium
! TODO: rewrite using derived type
call read_equil_an(params%equilibrium%filenm, &
params%raytracing%ipass, &
data%equilibrium%rv, &
data%equilibrium%zv, &
data%equilibrium%fpol, &
data%equilibrium%qpsi, &
data%equilibrium%rlim, &
data%equilibrium%zlim)
! Set psia sign to give the correct sign to Iphi
! (COCOS=3: psia<0 for Iphi>0)
data%equilibrium%psia = sign(one, data%equilibrium%qpsi(2) &
* data%equilibrium%fpol(1))
else
! Numerical equilibrium
call read_eqdsk(params%equilibrium, data%equilibrium)
call change_cocos(data%equilibrium, params%equilibrium%icocos, 3)
end if
! Rescale B, I and/or force their signs
call eq_scal(params%equilibrium, data%equilibrium)
! Set global variables (for splines)
if(params%equilibrium%iequil < 2) then
! TODO: rewrite using derived type
call set_equian(data%equilibrium%rv(1), &
data%equilibrium%zv(1), &
data%equilibrium%rv(2), &
data%equilibrium%fpol(1) / data%equilibrium%rv(1), &
data%equilibrium%qpsi(1), &
data%equilibrium%qpsi(2), &
data%equilibrium%qpsi(3))
else
call set_eqspl(params%equilibrium, data%equilibrium)
end if
end subroutine init_equilibrium
subroutine deinit_equilibrium(data)
! Free all memory allocated by the init_equilibrium subroutine.
use gray_params, only : equilibrium_data
use equilibrium, only : unset_eqspl, unset_rhospl, unset_q
implicit none
! subroutine arguments
type(equilibrium_data), intent(inout) :: data
! Free the MHD equilibrium arrays
if (allocated(data%rv)) deallocate(data%rv, data%zv, data%fpol, data%qpsi)
if (allocated(data%psin)) deallocate(data%psin, data%psinr)
if (allocated(data%rbnd)) deallocate(data%rbnd, data%zbnd)
if (allocated(data%rlim)) deallocate(data%rlim, data%zlim)
! Unset global variables of the `equilibrium` module
call unset_eqspl
call unset_rhospl
call unset_q
end subroutine deinit_equilibrium
subroutine init_profiles(params, data)
! Reads the plasma kinetic profiles file (containing the elecron
! temperature, density and plasma effective charge) and initialises
! the respective GRAY data structure.
use gray_params, only : profiles_parameters, profiles_data
use equilibrium, only : frhopolv
use coreprofiles, only : read_profiles_an, read_profiles, tene_scal, &
set_prfan, set_prfspl
implicit none
! subroutine arguments
type(gray_parameters), intent(in) :: params
type(gray_data), intent(inout), target :: data
! local variables
type(profiles_parameters) :: profp
type(profiles_data), pointer :: profd
! Radial coordinate (depending on profp%irho: ρ_t, ρ_p, or ψ)
real(wp_), allocatable :: xrad(:)
profp = params%profiles
profd => data%profiles
if(params%profiles%iprof == 0) then
! Analytical profiles
! TODO: rewrite using derived type
call read_profiles_an(profp%filenm, profd%terad, profd%derad, profd%zfc)
else
! Numerical profiles
call read_profiles(profp%filenm, xrad, profd%terad, profd%derad, profd%zfc)
allocate(profd%psrad(size(xrad)))
select case (profp%irho)
case (0) ! xrad is rhot
profd%psrad = frhopolv(xrad)**2
case (1) ! xrad is rhop
profd%psrad = xrad**2
case default ! xrad is psi
profd%psrad = xrad
end select
deallocate(xrad)
end if
! Rescale input data
! TODO: rewrite using derived type
call tene_scal(profd%terad, profd%derad, profp%factte, profp%factne, &
params%equilibrium%factb, profp%iscal, profp%iprof)
! Set global variables
! TODO: rewrite using derived type
if(params%profiles%iprof == 0) then
call set_prfan(profd%terad, profd%derad, profd%zfc)
else
call set_prfspl(profd%psrad, profd%terad, profd%derad, profd%zfc, &
profp%sspld, profp%psnbnd)
end if
end subroutine init_profiles
subroutine deinit_profiles(data)
! Free all memory allocated by the init_profiles subroutine.
use gray_params, only : profiles_data
use coreprofiles, only : unset_prfspl
implicit none
! subroutine arguments
type(profiles_data), intent(inout) :: data
! Free the plasma kinetic profiles arrays
if (allocated(data%psrad)) deallocate(data%psrad)
if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc)
! Unset global variables of the `coreprofiles` module
call unset_prfspl
end subroutine deinit_profiles
subroutine init_antenna(params)
! Reads the wave launcher file (containing the wave frequency, launcher
! position, direction and beam description) and initialises the respective
! GRAY parameters.
use beams, only : read_beam0, read_beam1, read_beam2
use gray_params, only : antenna_parameters
implicit none
! subroutine arguments
type(antenna_parameters), intent(inout) :: params
! Note: α, β are loaded from gray_params.data
select case (params%ibeam)
case (2)
! 2 degrees of freedom
! w(z, α, β), 1/R(z, α, β)
! FIXME: 1st beam is always selected, iox read from table
call read_beam2(params, beamid=1)
case (1)
! 1 degree of freedom
! w(z, α), 1/R(z, α)
call read_beam1(params)
case default
! fixed w(z), 1/R(z)
call read_beam0(params)
end select
end subroutine init_antenna
subroutine init_misc(params, data)
! Performs miscellanous initial tasks, before the gray_main subroutine.
use reflections, only : range2rect
use limiter, only : limiter_set_globals=>set_globals
implicit none
! subroutine arguments
type(gray_parameters), intent(inout) :: params
type(gray_data), intent(inout) :: data
! Build a basic limiter when one is not provided by the EQDSK
if (.not. allocated(data%equilibrium%rlim) &
.or. params%raytracing%ipass < 0) then
block
real(wp_) :: rmxm, r0m, z0m, dzmx
r0m = sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2)* 0.01_wp_
dzmx = abs(params%raytracing%ipass) * &
params%raytracing%dst * params%raytracing%nstep * 0.01_wp_
z0m = params%antenna%pos(3) * 0.01_wp_
allocate(data%equilibrium%rlim(5))
allocate(data%equilibrium%zlim(5))
params%raytracing%ipass = abs(params%raytracing%ipass)
if(params%equilibrium%iequil < 2) then
rmxm = (data%equilibrium%rv(1) + data%equilibrium%rv(2)) * 0.01_wp_
else
rmxm = data%equilibrium%rv(size(data%equilibrium%rv))
end if
call range2rect(params%misc%rwall, max(r0m, rmxm), &
z0m - dzmx, z0m + dzmx, &
data%equilibrium%rlim, data%equilibrium%zlim)
end block
end if
! Set the global variables of the `limiter` module
call limiter_set_globals(data%equilibrium)
end subroutine init_misc
subroutine deinit_misc
! Free all memory allocated by the init_misc subroutine.
use limiter, only : limiter_unset_globals=>unset_globals
implicit none
! Unset the global variables of the `limiter` module
call limiter_unset_globals
end subroutine deinit_misc
subroutine sum_profiles(params, jphi, jcd, dpdv, currins, pins, pabs, icd, &
jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
drhotjava, drhotpav, ratjamx, ratjbmx)
use const_and_precisions, only : zero, degree
use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
use dispersion, only : expinit
use gray_params, only : gray_parameters, print_parameters, &
headw, headl
use beams, only : launchangles2n, xgygcoeff
use magsurf_data, only : flux_average, dealloc_surfvec
use beamdata, only : init_btr, dealloc_beam
use pec, only : pec_init, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab
use graycore, only : print_headers, print_finals, print_pec, &
print_bres, print_prof, print_maps, &
print_surfq
implicit none
! subroutine arguments
type(gray_parameters), intent(in) :: params
real(wp_), intent(in) :: pabs, icd
real(wp_), dimension(:), intent(in) :: jphi, jcd, dpdv, currins, pins
real(wp_), intent(out) :: jphip, dpdvp, rhotj, rhotjava, &
rhotp, rhotpav, drhotjava, drhotpav, &
ratjamx,ratjbmx
! local variables
real(wp_) :: ak0, bres, xgcn
real(wp_) :: chipol, psipol, st
real(wp_) :: drhotp, drhotj, dpdvmx, jphimx
real(wp_), dimension(3) :: anv0
real(wp_), dimension(:, :), pointer :: yw=>null(), ypw=>null(), gri=>null()
real(wp_), dimension(:, :, :), pointer :: xc=>null(), du1=>null(), ggri=>null()
real(wp_), dimension(:, :), pointer :: psjki=>null(), ppabs=>null(), ccci=>null()
real(wp_), dimension(:), pointer :: tau0=>null(), alphaabs0=>null(), &
dids0=>null(), ccci0=>null()
real(wp_), dimension(:), pointer :: p0jk=>null()
complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
integer, dimension(:), pointer :: iiv=>null()
! ======== set environment BEGIN ========
! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1)
call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
! Compute the initial cartesian wavevector (anv0)
call launchangles2n(params%antenna, anv0)
! Initialise the ray variables (beamtracing)
call init_btr(params%raytracing, yw, ypw, xc, du1, &
gri, ggri, psjki, ppabs, ccci, &
tau0, alphaabs0, dids0, ccci0, &
p0jk, ext, eyt, iiv)
! Initialise the dispersion module
if(params%ecrh_cd%iwarm > 1) call expinit
! Initialise the magsurf_data module
call flux_average ! requires frhotor for dadrhot,dvdrhot
! Initialise the output profiles
call pec_init(params%output%ipec)
! ======= set environment END ======
! ======== pre-proc prints BEGIN ========
block
! Parameters log in file headers
character(len=headw), dimension(headl) :: strheader
call print_parameters(params, strheader)
call print_headers(strheader)
end block
! Print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout
call print_surfq([1.5_wp_, 2.0_wp_])
! Print ne, Te, q, Jphi versus psi, rhop, rhot
call print_bres(bres)
call print_prof
call print_maps(bres, xgcn, &
0.01_wp_*sqrt(params%antenna%pos(1)**2 + params%antenna%pos(2)**2), &
sin(params%antenna%beta*degree))
! ========= pre-proc prints END =========
! Print power and current density profiles
call print_pec(rhop_tab, rhot_tab, jphi, jcd, &
dpdv, currins, pins, index_rt=1)
! Compute profiles width
call postproc_profiles(pabs, icd, rhot_tab, dpdv, jphi, &
rhotpav, drhotpav, rhotjava, drhotjava, &
dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, &
dpdvmx, jphimx, ratjamx, ratjbmx)
! Print 0D results
call print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, &
drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp, &
drhotj, ratjamx, ratjbmx, st, psipol, chipol, &
1, params%antenna%power, cpl1=zero, cpl2=zero)
! Free memory
call dealloc_surfvec ! for fluxval
call dealloc_beam(yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, &
tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv)
call dealloc_pec
end subroutine sum_profiles
end program main