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

src: implement toggling of output units

Browse Source 
This actually implements the --units option
This commit is contained in:
Michele Guerini Rocco 2021-12-19 16:39:19 +01:00
parent 96359bc3fd
commit 91a2e6cf07
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
5 changed files with 374 additions and 264 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

473
src/gray_core.f90
View File

@ -1,34 +1,36 @@
! This modules contains the core GRAY routines
module gray_core module gray_core
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
implicit none implicit none
contains contains
subroutine gray_main(params, data, results, error, rhout) subroutine gray_main(params, data, results, error, 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 : temp, fzeff
use dispersion, only : expinit use dispersion, only : expinit
use gray_params, only : gray_parameters, gray_data, gray_results, print_parameters, & use gray_params, only : gray_parameters, gray_data, gray_results, &
print_parameters, &
iwarm, ipec, istpr0, igrad, headw, headl, ipass iwarm, ipec, istpr0, igrad, headw, headl, ipass
use beams, only : xgygcoeff, launchangles2n use beams, only : xgygcoeff, launchangles2n
use beamdata, only : pweight, rayi2jk use beamdata, only : pweight, rayi2jk
use equilibrium, only : unset_eqspl, unset_rhospl, unset_q
use errcodes, only : check_err, print_errn, print_errhcd 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 : 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, &
use units, only : ucenr wall_out
use logger, only : log_info, log_debug use logger, only : log_info, log_debug
implicit none implicit none
! Subroutine arguments ! subroutine arguments
type(gray_parameters), intent(in) :: params type(gray_parameters), intent(in) :: params
type(gray_data), intent(in) :: data type(gray_data), intent(in) :: data
type(gray_results), intent(out) :: results type(gray_results), intent(out) :: results
@ -554,8 +556,6 @@ contains
end if end if
end if end if
write(ucenr,*) ''
call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, & call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, &
pins_beam,ip) ! *print power and current density profiles for current beam pins_beam,ip) ! *print power and current density profiles for current beam
@ -609,7 +609,7 @@ contains
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 : zero
implicit none implicit none
! arguments ! arguments
real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci
@ -642,7 +642,7 @@ contains
ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt) ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt)
! beam tracing initial conditions igrad=1 ! beam tracing initial conditions igrad=1
! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! ! !!!!!! check ray tracing initial conditions igrad=0 !!!!!!
use const_and_precisions, only : wp_,zero,one,pi,half,two,degree,ui=>im use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im
use math, only : catand use math, only : catand
use gray_params, only : idst use gray_params, only : idst
use beamdata, only : nray,nrayr,nrayth,rwmax use beamdata, only : nray,nrayr,nrayth,rwmax
@ -944,7 +944,6 @@ contains
subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad) subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad)
! Runge-Kutta integrator ! Runge-Kutta integrator
use const_and_precisions, only : wp_
! use gray_params, only : igrad ! use gray_params, only : igrad
use beamdata, only : h,hh,h6 use beamdata, only : h,hh,h6
implicit none implicit none
@ -980,7 +979,6 @@ contains
subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery,igrad) subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery,igrad)
! Compute right-hand side terms of the ray equations (dery) ! Compute right-hand side terms of the ray equations (dery)
! used in R-K integrator ! used in R-K integrator
use const_and_precisions, only : wp_
implicit none implicit none
! arguments ! arguments
real(wp_), dimension(6), intent(in) :: y real(wp_), dimension(6), intent(in) :: y
@ -1049,7 +1047,7 @@ contains
subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri) subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri)
use const_and_precisions, only : wp_,zero,half use const_and_precisions, only : zero,half
use beamdata, only : nray,nrayr,nrayth,twodr2 use beamdata, only : nray,nrayr,nrayth,twodr2
implicit none implicit none
real(wp_), intent(in) :: ak0 real(wp_), intent(in) :: ak0
@ -1199,7 +1197,6 @@ contains
! input vectors : dxv1, dxv2, dxv3, dff ! input vectors : dxv1, dxv2, dxv3, dff
! output vector : dgg ! output vector : dgg
! dff=(1,0,0) ! dff=(1,0,0)
use const_and_precisions, only : wp_
implicit none implicit none
! arguments ! arguments
real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
@ -1220,7 +1217,6 @@ contains
subroutine solg3(dxv1,dxv2,dxv3,dff,dgg) subroutine solg3(dxv1,dxv2,dxv3,dff,dgg)
! rhs "matrix" dff, result in dgg ! rhs "matrix" dff, result in dgg
use const_and_precisions, only : wp_
implicit none implicit none
! arguments ! arguments
real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3
@ -1252,7 +1248,7 @@ contains
subroutine plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv, & subroutine plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv, &
xg,yg,derxg,deryg,ajphi) xg,yg,derxg,deryg,ajphi)
use const_and_precisions, only : wp_,zero,ccj=>mu0inv use const_and_precisions, only : zero,ccj=>mu0inv
use gray_params, only : iequil use gray_params, only : iequil
use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi
use coreprofiles, only : density use coreprofiles, only : density
@ -1406,7 +1402,7 @@ contains
subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & subroutine 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)
use const_and_precisions, only : wp_,zero,one,half,two use const_and_precisions, only : zero,one,half,two
use gray_params, only : idst use gray_params, only : idst
implicit none implicit none
! arguments ! arguments
@ -1533,7 +1529,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,error) sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,error)
use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_ use const_and_precisions, only : 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
use equilibrium, only : sgnbphi use equilibrium, only : sgnbphi
@ -1623,7 +1619,7 @@ contains
subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0) subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0)
use const_and_precisions, only : wp_,degree,zero,one,half,im use const_and_precisions, only : degree,zero,one,half,im
use beamdata, only : nray,nrayth use beamdata, only : nray,nrayth
use equilibrium, only : bfield use equilibrium, only : bfield
use gray_params, only : ipol use gray_params, only : ipol
@ -1683,10 +1679,8 @@ contains
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_
! 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.
! (based on an older code) ! (based on an older code)
use const_and_precisions, only : wp_
implicit none implicit none
! arguments ! arguments
integer, intent(in) :: nr,nz integer, intent(in) :: nr,nz
@ -1897,75 +1891,95 @@ bb: do
subroutine print_headers(strheader) subroutine print_headers(strheader)
use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm ! Prints the headers of all gray_main output tables
use units, only : uprj0, uwbm, udisp, ucenr, uoutr, upec, usumm, &
unit_active
implicit none implicit none
! subroutine arguments ! subroutine arguments
character(len=*), dimension(:), intent(in) :: strheader character(len=*), dimension(:), intent(in) :: strheader
! local variables
integer :: i,l
l=size(strheader) ! local variables
do i=1,l integer :: i, j
write(uprj0,'(1x,a)') strheader(i) integer :: main_units(8)
write(uprj0+1,'(1x,a)') strheader(i) character(256) :: main_headers(8)
write(uwbm,'(1x,a)') strheader(i)
write(udisp,'(1x,a)') strheader(i) main_units = [uprj0, uprj0+1, uwbm, udisp, ucenr, uoutr, upec, usumm]
write(ucenr,'(1x,a)') strheader(i)
write(uoutr,'(1x,a)') strheader(i) main_headers(1) = '#sst j k xt yt zt rt'
write(upec,'(1x,a)') strheader(i) main_headers(2) = '#sst j k xt yt zt rt'
write(usumm,'(1x,a)') strheader(i) main_headers(3) = '#sst w1 w2'
main_headers(4) = '#sst Dr_Nr Di_Nr'
main_headers(5) = '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp ' &
// 'Npl Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax ' &
// 'iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz'
main_headers(6) = '#i k sst x y R z psin tau Npl alpha index_rt'
main_headers(7) = '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
main_headers(8) = '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' &
// 'drhotjava drhotpav ratjamx ratjbmx stmx psipol ' &
// 'chipol index_rt Jphimx dPdVmx drhotj drhotp P0 ' &
// 'cplO cplX'
do i=1,size(main_units)
if (unit_active(main_units(i))) then
do j=1,size(strheader)
write (main_units(i), '(1x,a)') strheader(j)
end do
write (main_units(i), '(1x,a)') trim(main_headers(i))
end if
end do end do
write(uprj0,'(1x,a)') '#sst j k xt yt zt rt'
write(uprj0+1,'(1x,a)') '#sst j k xt yt zt rt'
write(uwbm,'(1x,a)') '#sst w1 w2'
write(udisp,'(1x,a)') '#sst Dr_Nr Di_Nr'
write(ucenr,'(1x,a)') '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp Npl '// &
'Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz'
write(uoutr,'(1x,a)') '#i k sst x y R z psin tau Npl alpha index_rt'
write(upec,'(1x,a)') '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
write(usumm,'(1x,a)') '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // &
'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt ' // &
'Jphimx dPdVmx drhotj drhotp P0 cplO cplX'
end subroutine print_headers end subroutine print_headers
subroutine print_prof subroutine print_prof
use const_and_precisions, only : wp_ ! Prints the (input) plasma kinetic profiles (unit 98)
use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi
use equilibrium, only : psinr, nq, fq, frhotor, tor_curr_psi
use coreprofiles, only : density, temp use coreprofiles, only : density, temp
use units, only : uprfin use units, only : uprfin, unit_active
implicit none implicit none
! local constants
real(wp_), parameter :: eps=1.e-4_wp_ ! local constants
! local variables real(wp_), parameter :: eps = 1.e-4_wp_
! local variables
integer :: i integer :: i
real(wp_) :: psin,rhot,ajphi,dens,ddens real(wp_) :: psin, rhot, ajphi, dens, ddens
write(uprfin,*) ' #psi rhot ne Te q Jphi' if (.not. unit_active(uprfin)) return
write (uprfin, *) '#psi rhot ne Te q Jphi'
do i=1,nq do i=1,nq
psin=psinr(i) psin = psinr(i)
rhot=frhotor(sqrt(psin)) rhot = frhotor(sqrt(psin))
call density(psin,dens,ddens) call density(psin, dens, ddens)
call tor_curr_psi(max(eps,psin),ajphi) call tor_curr_psi(max(eps, psin), ajphi)
write(uprfin,"(12(1x,e12.5))") psin,rhot,dens,temp(psin),fq(psin),ajphi*1.e-6_wp_ write (uprfin, '(12(1x,e12.5))') &
psin, rhot, dens, temp(psin), fq(psin), ajphi*1.e-6_wp_
end do end do
end subroutine print_prof end subroutine print_prof
subroutine print_bres(bres) subroutine print_bres(bres)
use const_and_precisions, only : wp_ ! Prints the EC resonance surface table (unit 70)
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq
use units, only : ubres use units, only : ubres, unit_active
implicit none implicit none
! arguments
real(wp_) :: bres ! subroutine arguments
! local constants real(wp_), intent(in) :: bres
integer, parameter :: icmx=2002
! local variables ! local constants
integer, parameter :: icmx = 2002
! local variables
integer :: j,k,n,nconts,nctot integer :: j,k,n,nconts,nctot
integer, dimension(10) :: ncpts integer, dimension(10) :: ncpts
real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
@ -1973,98 +1987,109 @@ bb: do
real(wp_) :: rv(nq), zv(nq) real(wp_) :: rv(nq), zv(nq)
real(wp_), dimension(nq,nq) :: btotal real(wp_), dimension(nq,nq) :: btotal
if (.not. unit_active(ubres)) return
dr = (rmxm-rmnm)/(nq-1) dr = (rmxm - rmnm)/(nq - 1)
dz = (zmxm-zmnm)/(nq-1) dz = (zmxm - zmnm)/(nq - 1)
do j=1,nq do j=1,nq
rv(j) = rmnm + dr*(j-1) rv(j) = rmnm + dr*(j - 1)
zv(j) = zmnm + dz*(j-1) zv(j) = zmnm + dz*(j - 1)
end do end do
! Btotal on psi grid ! Btotal on psi grid
btmx=-1.0e30_wp_ btmx = -1.0e30_wp_
btmn=1.0e30_wp_ btmn = 1.0e30_wp_
do k=1,nq do k=1,nq
zzk=zv(k) zzk = zv(k)
do j=1,nq do j=1,nq
rrj=rv(j) rrj = rv(j)
call bfield(rrj,zzk,bbphi,bbr,bbz) call bfield(rrj, zzk, bbphi, bbr, bbz)
btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2) btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2)
if(btotal(j,k).ge.btmx) btmx=btotal(j,k) if(btotal(j,k) >= btmx) btmx = btotal(j,k)
if(btotal(j,k).le.btmn) btmn=btotal(j,k) if(btotal(j,k) <= btmn) btmn = btotal(j,k)
enddo end do
enddo end do
! compute Btot=Bres/n with n=1,5 ! compute Btot=Bres/n with n=1,5
write(ubres,*)'#i Btot R z' write (ubres, *) '#i Btot R z'
do n=1,5 do n=1,5
bbb=bres/dble(n) bbb = bres/dble(n)
if (bbb.ge.btmn.and.bbb.le.btmx) then if (bbb >= btmn .and. bbb <= btmx) then
nconts=size(ncpts) nconts = size(ncpts)
nctot=size(rrcb) nctot = size(rrcb)
call cniteq(rv,zv,btotal,nq,nq,bbb,nconts,ncpts,nctot,rrcb,zzcb) call cniteq(rv, zv, btotal, nq, nq, bbb, &
nconts, ncpts, nctot, rrcb, zzcb)
do j=1,nctot do j=1,nctot
write(ubres,'(i6,12(1x,e12.5))') j,bbb,rrcb(j),zzcb(j) write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j)
end do end do
end if end if
write(ubres,*) write (ubres, *)
end do end do
end subroutine print_bres end subroutine print_bres
subroutine print_maps(bres, xgcn, r0, anpl0)
! Prints several input quantities on a regular grid (unit 72)
subroutine print_maps(bres,xgcn,r0,anpl0)
use const_and_precisions, only : wp_
use gray_params, only : iequil use gray_params, only : iequil
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, & use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, &
equinum_fpol, nq equinum_fpol, nq
use coreprofiles, only : density, temp use coreprofiles, only : density, temp
use units, only : umaps use units, only : umaps, unit_active
implicit none implicit none
! arguments
! subroutine arguments
real(wp_), intent(in) :: bres,xgcn,r0,anpl0 real(wp_), intent(in) :: bres,xgcn,r0,anpl0
! local variables
! local variables
integer :: j,k integer :: j,k
real(wp_) :: dr,dz,zk,rj,bphi,br,bz,btot,psin,ne,dne,te,xg,yg,anpl real(wp_) :: dr, dz, zk, rj, bphi, br, bz, btot, &
psin, ne, dne, te, xg, yg, anpl
real(wp_), dimension(nq) :: r, z real(wp_), dimension(nq) :: r, z
if (.not. unit_active(umaps)) return
dr = (rmxm-rmnm)/(nq-1) dr = (rmxm-rmnm)/(nq - 1)
dz = (zmxm-zmnm)/(nq-1) dz = (zmxm-zmnm)/(nq - 1)
do j=1,nq do j=1,nq
r(j) = rmnm + dr*(j-1) r(j) = rmnm + dr*(j - 1)
z(j) = zmnm + dz*(j-1) z(j) = zmnm + dz*(j - 1)
end do end do
write(umaps,*)'#R z psin Br Bphi Bz Btot ne Te X Y Npl' write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl'
do j=1,nq do j=1,nq
rj=r(j) rj = r(j)
anpl=anpl0*r0/rj anpl = anpl0 * r0/rj
do k=1,nq do k=1,nq
zk=z(k) zk = z(k)
if (iequil < 2) then if (iequil < 2) then
call equian(rj,zk,psinv=psin,fpolv=bphi,dpsidr=bz,dpsidz=br) call equian(rj, zk, psinv=psin, fpolv=bphi, dpsidr=bz, dpsidz=br)
else else
call equinum_psi(rj,zk,psinv=psin,dpsidr=bz,dpsidz=br) call equinum_psi(rj, zk, psinv=psin, dpsidr=bz, dpsidz=br)
call equinum_fpol(psin,fpolv=bphi) call equinum_fpol(psin, fpolv=bphi)
end if end if
br = -br/rj br = -br/rj
bphi = bphi/rj bphi = bphi/rj
bz = bz/rj bz = bz/rj
btot = sqrt(br**2+bphi**2+bz**2) btot = sqrt(br**2 + bphi**2 + bz**2)
yg = btot/bres yg = btot/bres
te = temp(psin) te = temp(psin)
call density(psin,ne,dne) call density(psin, ne, dne)
xg = xgcn*ne xg = xgcn*ne
write(umaps,'(12(x,e12.5))') rj,zk,psin,br,bphi,bz,btot,ne,te,xg,yg,anpl write (umaps,'(12(x,e12.5))') &
enddo rj, zk, psin, br, bphi, bz, btot, ne, te, xg, yg, anpl
write(umaps,*) end do
enddo write (umaps,*)
end do
end subroutine print_maps end subroutine print_maps
subroutine print_surfq(qval) subroutine print_surfq(qval)
! Print constant ψ surfaces for a given `q` value
use equilibrium, only : psinr, nq, fq, frhotor, & use equilibrium, only : psinr, nq, fq, frhotor, &
rmaxis, zmaxis, zbsup, zbinf rmaxis, zmaxis, zbsup, zbinf
use magsurf_data, only : contours_psi, npoints, print_contour use magsurf_data, only : contours_psi, npoints, print_contour
@ -2094,16 +2119,16 @@ bb: do
do i=1,size(qval) do i=1,size(qval)
! FIXME: check for non monotonous q profile ! FIXME: check for non monotonous q profile
call locate(abs(qpsi),nq,qval(i),i1) call locate(abs(qpsi),nq,qval(i),i1)
if (i1>0.and.i1<nq) then if (i1>0 .and. i1<nq) then
call intlin(abs(qpsi(i1)),psinr(i1),abs(qpsi(i1+1)),psinr(i1+1), & call intlin(abs(qpsi(i1)), psinr(i1), abs(qpsi(i1+1)), psinr(i1+1), &
qval(i),psival) qval(i),psival)
rup=rmaxis rup = rmaxis
rlw=rmaxis rlw = rmaxis
zup=(zbsup+zmaxis)/2.0_wp_ zup = (zbsup + zmaxis)/2.0_wp_
zlw=(zmaxis+zbinf)/2.0_wp_ zlw = (zmaxis + zbinf)/2.0_wp_
call contours_psi(psival,rcn,zcn,rup,zup,rlw,zlw) call contours_psi(psival, rcn, zcn, rup, zup, rlw, zlw)
call print_contour(psival,rcn,zcn) call print_contour(psival, rcn, zcn)
rhot=frhotor(sqrt(psival)) rhot = frhotor(sqrt(psival))
write (msg, '(4(x,a,"=",g0.3))') & write (msg, '(4(x,a,"=",g0.3))') &
'q', qval(i), 'ψ', psival, 'rhop', sqrt(psival), 'rhot', rhot 'q', qval(i), 'ψ', psival, 'rhop', sqrt(psival), 'rhot', rhot
call log_info(msg, mod='gray_core', proc='print_surfq') call log_info(msg, mod='gray_core', proc='print_surfq')
@ -2113,36 +2138,44 @@ bb: do
end subroutine print_surfq end subroutine print_surfq
subroutine print_projxyzt(stv,ywrk,iproj) subroutine print_projxyzt(stv, ywrk, iproj)
use const_and_precisions, only : wp_, comp_huge, zero, one ! Prints the beam shape tables (unit 8, 12)
use const_and_precisions, only : comp_huge, zero, one
use beamdata, only : nray, nrayr, nrayth, rayi2jk use beamdata, only : nray, nrayr, nrayth, rayi2jk
use units, only : uprj0,uwbm use units, only : uprj0, uwbm, unit_active
implicit none implicit none
! arguments
! subroutine arguments
real(wp_), dimension(:), intent(in) :: stv real(wp_), dimension(:), intent(in) :: stv
real(wp_), dimension(:,:), intent(in) :: ywrk real(wp_), dimension(:,:), intent(in) :: ywrk
integer, intent(in) :: iproj integer, intent(in) :: iproj
! local variables
integer :: jk,jkz,uprj ! local variables
integer :: jk, jkz, uprj
integer, dimension(2) ::jkv integer, dimension(2) ::jkv
real(wp_), dimension(3) :: xv1,dir,dxv real(wp_), dimension(3) :: xv1, dir, dxv
real(wp_) :: dirm,rtimn,rtimx,csth1,snth1,csps1,snps1,xti,yti,zti,rti real(wp_) :: rtimn, rtimx, csth1, snth1, &
! common/external functions/variables csps1, snps1, xti, yti, zti, rti
if (.not.(unit_active(uprj0) .or. unit_active(uwbm))) return
! common/external functions/variables
uprj = uprj0 + iproj uprj = uprj0 + iproj
xv1 = ywrk(1:3,1) xv1 = ywrk(1:3, 1)
dir = ywrk(4:6,1) dir = ywrk(4:6, 1)
dirm = sqrt(dir(1)**2 + dir(2)**2 + dir(3)**2) dir = dir/norm2(dir)
dir = dir/dirm
csth1 = dir(3) csth1 = dir(3)
snth1 = sqrt(one - csth1**2) snth1 = sqrt(one - csth1**2)
if(snth1 > zero) then if(snth1 > zero) then
csps1=dir(2)/snth1 csps1 = dir(2)/snth1
snps1=dir(1)/snth1 snps1 = dir(1)/snth1
else else
csps1=one csps1 = one
snps1=zero snps1 = zero
end if end if
if(iproj==0) then if(iproj==0) then
@ -2158,116 +2191,148 @@ bb: do
xti = dxv(1)*csps1 - dxv(2)*snps1 xti = dxv(1)*csps1 - dxv(2)*snps1
yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1 yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1
zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1 zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1
rti = sqrt(xti**2 + yti**2) rti = hypot(xti, yti)
jkv=rayi2jk(jk) jkv=rayi2jk(jk)
if(.not.(iproj==0 .and. jk==1)) & if(.not.(iproj==0 .and. jk==1)) &
write(uprj,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') stv(jk),jkv,xti,yti,zti,rti write(uprj, '(1x,e16.8e3,2i5,4(1x,e16.8e3))') &
if(iproj==1 .and. jkv(2)==nrayth) write(uprj,*) stv(jk), jkv, xti, yti, zti, rti
if(iproj==1 .and. jkv(2)==nrayth) write (uprj, *)
if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti
if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti
end do end do
write(uprj,*) write (uprj, *)
write(uwbm,'(3(1x,e16.8e3))') stv(1),rtimn,rtimx write (uwbm, '(3(1x,e16.8e3))') stv(1), rtimn, rtimx
end subroutine print_projxyzt end subroutine print_projxyzt
subroutine print_output(i, jk, st, qj, xv, psinv, btot, bv, ak0, &
anpl, anpr, anv, anprim, dens, tekev, &
alpha, tau, dids, nhm, nhf, iokhawa, &
index_rt, ddr, ddi, xg, yg, derxg)
! Prints the ray variables and Hamiltonian error (units 4, 33, 17)
subroutine print_output(i,jk,st,qj,xv,psinv,btot,bv,ak0,anpl,anpr,anv, &
anprim,dens,tekev,alpha,tau,dids,nhm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg)
use const_and_precisions, only : degree,zero,one use const_and_precisions, only : degree,zero,one
use equilibrium, only : frhotor use equilibrium, only : frhotor
use gray_params, only : istpl0 use gray_params, only : istpl0
use beamdata, only : nray,nrayth,jkray1 use beamdata, only : nray,nrayth,jkray1
use units, only : ucenr,uoutr,udisp use units, only : ucenr, uoutr, udisp, unit_active
implicit none implicit none
! arguments
integer, intent(in) :: i,jk,nhm,nhf,iokhawa,index_rt ! subroutine arguments
real(wp_), dimension(3), intent(in) :: xv,bv,anv integer, intent(in) :: i, jk, nhm, nhf, iokhawa, index_rt
real(wp_), intent(in) :: st,qj,psinv,btot,ak0,anpl,anpr,anprim real(wp_), dimension(3), intent(in) :: xv, bv, anv
real(wp_), intent(in) :: dens,tekev,alpha,tau,dids,ddr,ddi real(wp_), intent(in) :: st, qj, psinv, btot, ak0, anpl, anpr, anprim
real(wp_), intent(in) :: xg,yg real(wp_), intent(in) :: dens, tekev, alpha, tau, dids, ddr, ddi
real(wp_), intent(in) :: xg, yg
real(wp_), dimension(3), intent(in) :: derxg real(wp_), dimension(3), intent(in) :: derxg
! local variables ! local variables
real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,didsn real(wp_) :: stm, xxm, yym, zzm, rrm, phideg, rhot, akim, pt, didsn
integer :: k integer :: k
stm=st*1.0e-2_wp_ stm = st * 1.0e-2_wp_
xxm=xv(1)*1.0e-2_wp_ xxm = xv(1) * 1.0e-2_wp_
yym=xv(2)*1.0e-2_wp_ yym = xv(2) * 1.0e-2_wp_
zzm=xv(3)*1.0e-2_wp_ zzm = xv(3) * 1.0e-2_wp_
rrm=sqrt(xxm**2 + yym**2) rrm = hypot(xxm, yym)
! print central ray trajectory. dIds in A/m/W, ki in m^-1 ! print central ray trajectory. dIds in A/m/W, ki in m⁻¹
if(jk.eq.1) then if(jk == 1) then
phideg=atan2(yym,xxm)/degree phideg = atan2(yym, xxm)/degree
if(psinv>=zero .and. psinv<=one) then if(psinv>=zero .and. psinv<=one) then
rhot=frhotor(sqrt(psinv)) rhot = frhotor(sqrt(psinv))
else else
rhot=1.0_wp_ rhot = 1.0_wp_
end if end if
akim=anprim*ak0*1.0e2_wp_ akim = anprim * ak0 * 1.0e2_wp_
pt=qj*exp(-tau) pt = qj*exp(-tau)
didsn=dids*1.0e2_wp_/qj didsn = dids * 1.0e2_wp_/qj
write(ucenr,'(22(1x,e16.8e3),4i5,6(1x,e16.8e3))') stm,rrm,zzm,phideg, & ! print central ray data
psinv,rhot,dens,tekev,btot,bv,anpr,anpl,anv,akim,alpha,tau,pt,didsn, & if (unit_active(ucenr)) then
nhm,nhf,iokhawa,index_rt,ddr,xg,yg,derxg write (ucenr,'(22(1x,e16.8e3),4i5,6(1x,e16.8e3))') &
stm, rrm, zzm, phideg, psinv, rhot, dens, tekev, &
btot, bv, anpr, anpl, anv, akim, alpha, tau, pt, &
didsn, nhm, nhf, iokhawa, index_rt, ddr, xg, yg, &
derxg
end if
end if end if
! print conservation of dispersion relation ! print conservation of dispersion relation (Hamiltonian error)
if(jk==nray) write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi if(unit_active(udisp) .and. jk==nray) then
write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi
end if
! print outer trajectories ! print outer ray trajectories
if(mod(i,istpl0)==0) then if(mod(i, istpl0) == 0) then
k = jk - jkray1 + 1 k = jk - jkray1 + 1
if(k>0 .and. k<=nrayth) then if(unit_active(uoutr) .and. k>0 .and. k<=nrayth) then
write(uoutr,'(2i5,9(1x,e16.8e3),i5)') i,k,stm,xxm,yym,rrm,zzm, & write(uoutr, '(2i5,9(1x,e16.8e3),i5)') &
psinv,tau,anpl,alpha,index_rt i, k, stm, xxm, yym, rrm, zzm, &
psinv, tau, anpl, alpha, index_rt
end if end if
end if end if
end subroutine print_output end subroutine print_output
subroutine print_pec(rhop_tab, rhot_tab, jphi, jcd, &
dpdv, currins, pins, index_rt)
! Prints the ECRH and CD profiles (unit 48)
use units, only : upec, unit_active
subroutine print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt)
use const_and_precisions, only : wp_
use units, only : upec
implicit none implicit none
! arguments
real(wp_), dimension(:), intent(in) :: rhop_tab,rhot_tab,jphi,jcd,dpdv, & ! subroutine arguments
currins,pins real(wp_), dimension(:), intent(in) :: rhop_tab, rhot_tab, jphi, jcd, &
dpdv, currins, pins
integer, intent(in) :: index_rt integer, intent(in) :: index_rt
! local variables
! local variables
integer :: i integer :: i
do i=1,size(rhop_tab) if (.not. unit_active(upec)) return
write(upec,'(7(1x,e16.8e3),i5)') rhop_tab(i),rhot_tab(i), &
jphi(i),jcd(i),dpdv(i),currins(i),pins(i),index_rt do i=1, size(rhop_tab)
write(upec, '(7(1x, e16.8e3), i5)') rhop_tab(i), rhot_tab(i), &
jphi(i), jcd(i), dpdv(i), currins(i), pins(i), index_rt
end do end do
write(upec,*) '' write(upec, *) ''
end subroutine print_pec end subroutine print_pec
subroutine print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, &
drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp, &
drhotj, ratjamx, ratjbmx, stmx, psipol, chipol, &
index_rt, p0, cpl1, cpl2)
! Prints global data and final results (unit 7)
use units, only : ucenr, usumm, unit_active
subroutine print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, &
drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, &
stmx,psipol,chipol,index_rt,p0,cpl1,cpl2)
use const_and_precisions, only : wp_
use units, only : usumm
implicit none implicit none
real(wp_), intent(in) :: pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, &
drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, & ! subroutine arguments
stmx,psipol,chipol,p0,cpl1,cpl2 real(wp_), intent(in) :: pabs, icd, dpdvp, jphip, rhotpav, rhotjava, &
drhotpav, drhotjava, dpdvmx, jphimx, rhotp, &
rhotj, drhotp, drhotj, ratjamx, ratjbmx, &
stmx, psipol, chipol, p0, cpl1, cpl2
integer, intent(in) :: index_rt integer, intent(in) :: index_rt
write(usumm,'(15(1x,e12.5),i5,7(1x,e12.5))') icd,pabs,jphip,dpdvp, & if (unit_active(ucenr)) then
rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, & write(ucenr, *) ''
stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp,p0, & end if
cpl1,cpl2
! write(usumm,*) '' if (unit_active(usumm)) then
write(usumm,'(15(1x,e12.5),i5,7(1x,e12.5))') &
icd, pabs, jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
drhotjava, drhotpav, ratjamx, ratjbmx, stmx, psipol, &
chipol, index_rt, jphimx, dpdvmx, drhotj, drhotp, &
p0, cpl1, cpl2
end if
end subroutine print_finals end subroutine print_finals
end module gray_core end module gray_core

13
src/gray_jetto1beam.f90
View File

@ -127,19 +127,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
call unset_prfspl call unset_prfspl
end block free_memory end block free_memory
! Close output units used for debugging
close_debug: block
use units, only : ucenr, usumm, uprj0, uwbm, udisp, ubres, &
ucnt, uoutr, ueq, uprfin, uflx, upec
integer :: i, debug_units(13)
debug_units = [ucenr, usumm, uprj0, uprj0+1, uwbm, udisp, &
ubres, ucnt, uoutr, ueq, uprfin, uflx, upec]
do i=1,size(debug_units)
close(debug_units(i))
end do
end block close_debug
! Copy over the results ! Copy over the results
pabs = res%pabs pabs = res%pabs
icd = res%icd icd = res%icd

51
src/magsurf_data.f90
View File

@ -535,44 +535,57 @@ contains
end subroutine contours_psi end subroutine contours_psi
subroutine print_contour(psin, rc, zc)
! Prints the flux surface contours (unit 71)
subroutine print_contour(psin,rc,zc)
use const_and_precisions, only : wp_, comp_tiny use const_and_precisions, only : wp_, comp_tiny
use units, only : ucnt use units, only : ucnt, unit_active
implicit none implicit none
! arguments
! subroutine arguments
real(wp_), intent(in) :: psin real(wp_), intent(in) :: psin
real(wp_), dimension(:), intent(in) :: rc, zc real(wp_), dimension(:), intent(in) :: rc, zc
! local variables
integer :: npoints,ic ! local variables
integer :: npoints, ic
if (.not. unit_active(ucnt)) return
if (psin < comp_tiny) then if (psin < comp_tiny) then
write(ucnt,*)' #i psin R z' write(ucnt, *)' #i psin R z'
else else
npoints=size(rc) npoints = size(rc)
do ic=1,npoints do ic=1,npoints
write(ucnt,'(i6,12(1x,e12.5))') ic,psin,rc(ic),zc(ic) write(ucnt, '(i6,12(1x,e12.5))') ic, psin, rc(ic), zc(ic)
end do end do
write(ucnt,*) write(ucnt, *)
write(ucnt,*) write(ucnt, *)
end if end if
end subroutine print_contour end subroutine print_contour
subroutine print_fluxav(psin, rhot, bav, bmx, bmn, area, vol, &
currp, ajphiav, ffc, ratja, ratjb, qq)
! Prints the flux-averaged quantities table (unit 56)
subroutine print_fluxav(psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, &
ffc,ratja,ratjb,qq)
use const_and_precisions, only : wp_, comp_tiny use const_and_precisions, only : wp_, comp_tiny
use units, only : uflx use units, only : uflx, unit_active
implicit none implicit none
! arguments
real(wp_), intent(in) :: psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, & ! subroutine arguments
ffc,ratja,ratjb,qq real(wp_), intent(in) :: psin, rhot, bav, bmx, bmn, area, vol, &
currp, ajphiav, ffc, ratja, ratjb, qq
if (.not. unit_active(uflx)) return
if (psin < comp_tiny) & if (psin < comp_tiny) &
write(uflx,*)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol |I_pl| <J_phi> fc ratJa ratJb qq' write(uflx, *)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol' // &
write(uflx,'(20(1x,e12.5))') psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, & ' |I_pl| <J_phi> fc ratJa ratJb qq'
ffc,ratja,ratjb,qq write(uflx, '(20(1x,e12.5))') &
psin, rhot, bav, bmx, bmn, area, vol, currp, ajphiav, &
ffc, ratja, ratjb, qq
end subroutine print_fluxav end subroutine print_fluxav
end module magsurf_data end module magsurf_data

4
src/main.f90
View File

@ -1,6 +1,7 @@
program main program main
use const_and_precisions, only : wp_, one, zero use const_and_precisions, only : wp_, one, zero
use logger, only : INFO, ERROR, set_log_level, log_message use logger, only : INFO, ERROR, set_log_level, log_message
use units, only : set_active_units
use gray_cli, only : cli_options, parse_cli_options use gray_cli, only : cli_options, parse_cli_options
use gray_core, only : gray_main use gray_core, only : gray_main
use gray_params, only : gray_parameters, gray_data, gray_results, & use gray_params, only : gray_parameters, gray_data, gray_results, &
@ -23,6 +24,9 @@ program main
if (opts%quiet) opts%verbose = ERROR if (opts%quiet) opts%verbose = ERROR
call set_log_level(opts%verbose) call set_log_level(opts%verbose)
! Activate the given output units
call set_active_units(opts%units)
! Load the parameters and also copy them into ! Load the parameters and also copy them into
! global variables exported by the gray_params ! global variables exported by the gray_params
call read_parameters(opts%params_file, params) call read_parameters(opts%params_file, params)

55
src/units.f90
View File

@ -1,15 +1,56 @@
! This modules contains the definitions of the unit numbers used for
! all the GRAY output ifles and a mechanism to toggle them
module units module units
! STANDARD
! Unit numbers
#ifdef JINTRAC
! JINTRAC
integer, parameter :: uprm =602, ubeam =603, uprf =644, ueq =644
integer, parameter :: uprfin =645, uflx =646, ubres =630, ucnt =631
integer, parameter :: uprj0 =608, uwbm =612, ucenr =604, uoutr =633
integer, parameter :: udisp =617, upec =648, usumm =607, umaps =632
integer, parameter :: uhead =638
#else
! STANDARD
integer, parameter :: uprm = 2, ubeam = 97, uprf = 98, ueq = 99 integer, parameter :: uprm = 2, ubeam = 97, uprf = 98, ueq = 99
integer, parameter :: uprfin = 55, uflx = 56, ubres = 70, ucnt = 71 integer, parameter :: uprfin = 55, uflx = 56, ubres = 70, ucnt = 71
integer, parameter :: uprj0 = 8, uwbm = 12, ucenr = 4, uoutr = 33 integer, parameter :: uprj0 = 8, uwbm = 12, ucenr = 4, uoutr = 33
integer, parameter :: udisp = 17, upec = 48, usumm = 7, umaps = 72 integer, parameter :: udisp = 17, upec = 48, usumm = 7, umaps = 72
integer, parameter :: uhead = 78 integer, parameter :: uhead = 78
! JINTRAC #endif
! integer, parameter :: uprm =602, ubeam =603, uprf =644, ueq =644
! integer, parameter :: uprfin =645, uflx =646, ubres =630, ucnt =631 ! List of active units
! integer, parameter :: uprj0 =608, uwbm =612, ucenr =604, uoutr =633 integer, allocatable, save :: active_units(:)
! integer, parameter :: udisp =617, upec =648, usumm =607, umaps =632
! integer, parameter :: uhead =638 contains
subroutine set_active_units(units)
! Activate the given units
!
! All units are inactive by default and no output
! will be directed to them.
implicit none
! subroutine arguments
integer, intent(in) :: units(:)
active_units = units
end subroutine set_active_units
function unit_active(unit) result(on)
! Checks whether the given `unit` is active
implicit none
! function arguments
integer, intent(in) :: unit
logical :: on
on = any(active_units == unit)
end function unit_active
end module units end module units