abstract the outputs units
This change replaces the output files (Fortran units) with a derived type called table, that hold the data in memory until further processing. The data stored in a table can be dumped to a file, as before, or processed in other ways, for example converted to other derived type.
This commit is contained in:
Michele Guerini Rocco
parent
624fbe3ec1
commit
15fc89110a
GPG Key ID:
BFBAF4C975F76450
@ -19,7 +19,7 @@ module gray_cli
|
||||
character(len=:), allocatable :: sum_filelist
|
||||
! others
|
||||
integer :: verbose
|
||||
integer, allocatable :: units(:)
|
||||
integer, allocatable :: tables(:)
|
||||
end type
|
||||
|
||||
private
|
||||
@ -54,8 +54,8 @@ contains
|
||||
print '(a)', ' -c, --config-file FILE set the (new) GRAY config file'
|
||||
print '(a)', ' (default: gray.ini)'
|
||||
print '(a)', ' -s, --sum FILE sum the output profiles from a list of files'
|
||||
print '(a)', ' -u, --units ID[,ID...] select which units to output (default: 4, 7);'
|
||||
print '(a)', ' use `all` to enable all, or `none` for no units.'
|
||||
print '(a)', ' -t, --tables ID[,ID...] select which tables to write (default: 4, 7);'
|
||||
print '(a)', ' use `all` to enable all, or `none` for no tables.'
|
||||
print '(a)', ' see the manual for all unit IDs.'
|
||||
print '(a)', ' -g, --gray-param ID=VAL set a GRAY parameter, overriding the value'
|
||||
print '(a)', ' specified via --params-file/--config-file;'
|
||||
@ -94,7 +94,7 @@ contains
|
||||
print '(a,a)' , ' param-file: ' , opts%params_file
|
||||
print '(a,a)' , ' sum: ' , opts%sum_filelist
|
||||
print '(a)' , 'others:'
|
||||
print '(a,20i3)' , ' - units: ' , opts%units
|
||||
print '(a,20i3)' , ' - tables: ' , opts%tables
|
||||
print '(a,l)' , ' - verbose: ' , opts%verbose
|
||||
end subroutine
|
||||
|
||||
@ -102,7 +102,6 @@ contains
|
||||
subroutine parse_cli_options(opts)
|
||||
! Parse the CLI arguments and initialise the options
|
||||
|
||||
use units, only : ucenr, usumm, all_enabled
|
||||
use logger, only : WARNING
|
||||
|
||||
! subroutine arguments
|
||||
@ -117,7 +116,7 @@ contains
|
||||
opts%verbose = WARNING
|
||||
opts%quiet = .false.
|
||||
opts%params_file = 'gray_params.data'
|
||||
opts%units = [ucenr, usumm]
|
||||
opts%tables = [4, 7]
|
||||
|
||||
nargs = command_argument_count()
|
||||
i = 1
|
||||
@ -156,24 +155,25 @@ contains
|
||||
case ('-s', '--sum')
|
||||
call get_next_command(i, opts%sum_filelist)
|
||||
|
||||
case ('-u', '--units')
|
||||
case ('-t', '--tables')
|
||||
call get_next_command(i, temp)
|
||||
|
||||
if (temp == 'none') then
|
||||
! disable all output units
|
||||
deallocate(opts%units)
|
||||
allocate(opts%units(0))
|
||||
! disable all output tables
|
||||
deallocate(opts%tables)
|
||||
allocate(opts%tables(0))
|
||||
elseif (temp == 'all') then
|
||||
! enable all output units
|
||||
all_enabled = .true.
|
||||
! enable all output tables
|
||||
deallocate(opts%tables)
|
||||
opts%tables = [-1]
|
||||
else
|
||||
! resize the array
|
||||
commas = count([(temp(i:i) == ',', i = 1, len(temp))])
|
||||
deallocate(opts%units)
|
||||
allocate(opts%units(commas + 1))
|
||||
deallocate(opts%tables)
|
||||
allocate(opts%tables(commas + 1))
|
||||
|
||||
! read the list of table IDs
|
||||
read (temp, *, iostat=error) opts%units
|
||||
read (temp, *, iostat=error) opts%tables
|
||||
if (error > 0) then
|
||||
print '(a,a)', 'invalid table IDs: ', temp
|
||||
deallocate(argument)
|
||||
@ -296,7 +296,7 @@ contains
|
||||
if (allocated(opts%params_file)) deallocate(opts%params_file)
|
||||
if (allocated(opts%config_file)) deallocate(opts%config_file)
|
||||
if (allocated(opts%sum_filelist)) deallocate(opts%sum_filelist)
|
||||
if (allocated(opts%units)) deallocate(opts%units)
|
||||
if (allocated(opts%tables)) deallocate(opts%tables)
|
||||
end subroutine deinit_cli_options
|
||||
|
||||
end module gray_cli
|
||||
|
||||
@ -8,16 +8,19 @@ module gray_core
|
||||
contains
|
||||
|
||||
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, comp_tiny
|
||||
use coreprofiles, only : temp, fzeff
|
||||
use dispersion, only : expinit
|
||||
use polarization, only : ellipse_to_field
|
||||
use gray_params, only : gray_parameters, gray_data, gray_results, &
|
||||
print_parameters, EQ_VACUUM
|
||||
use types, only : table, wrap
|
||||
use gray_params, only : gray_parameters, gray_data, gray_results, EQ_VACUUM
|
||||
use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, &
|
||||
store_beam_shape, find_flux_surfaces, &
|
||||
kinetic_profiles, ec_resonance, inputs_maps
|
||||
use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd
|
||||
use beams, only : xgygcoeff, launchangles2n
|
||||
use beamdata, only : pweight, rayi2jk
|
||||
use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd
|
||||
use magsurf_data, only : flux_average, dealloc_surfvec
|
||||
use magsurf_data, only : compute_flux_averages, dealloc_surfvec
|
||||
use beamdata, only : init_btr
|
||||
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
|
||||
rhop_tab, rhot_tab
|
||||
@ -106,6 +109,9 @@ contains
|
||||
! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1)
|
||||
call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
|
||||
|
||||
! Initialise the output tables
|
||||
call init_tables(params, results%tables)
|
||||
|
||||
! Compute the initial cartesian wavevector (anv0)
|
||||
! from launch angles α,β and the position
|
||||
call launchangles2n(params%antenna, anv0)
|
||||
@ -118,7 +124,7 @@ contains
|
||||
|
||||
if (params%equilibrium%iequil /= EQ_VACUUM) then
|
||||
! Initialise the magsurf_data module
|
||||
call flux_average ! requires frhotor for dadrhot,dvdrhot
|
||||
call compute_flux_averages(results%tables)
|
||||
|
||||
! Initialise the output profiles
|
||||
call pec_init(params%output%ipec, rhout)
|
||||
@ -136,11 +142,13 @@ contains
|
||||
results%jcd = zero
|
||||
! ========= set environment END =========
|
||||
|
||||
! ======== pre-proc prints BEGIN ========
|
||||
call print_headers(params)
|
||||
! Pre-determinted tables
|
||||
results%tables%kinetic_profiles = kinetic_profiles(params)
|
||||
results%tables%ec_resonance = ec_resonance(params, bres)
|
||||
results%tables%inputs_maps = inputs_maps(params, bres, xgcn)
|
||||
|
||||
! print ψ surface for q=1.5 and q=2 on file and log psi,rhot,rhop
|
||||
call print_surfq([1.5_wp_, 2.0_wp_])
|
||||
! print ψ surface for q=3/2 and q=2/1
|
||||
call find_flux_surfaces([1.5_wp_, 2.0_wp_], results%tables%flux_surfaces)
|
||||
|
||||
! print initial position
|
||||
write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos
|
||||
@ -150,17 +158,6 @@ contains
|
||||
'α', params%antenna%alpha, 'β', params%antenna%beta
|
||||
call log_info(msg, mod='gray_core', proc='gray_main')
|
||||
|
||||
! print Btot=Bres
|
||||
! print ne, Te, q, Jphi versus psi, rhop, rhot
|
||||
if (params%equilibrium%iequil /= EQ_VACUUM) then
|
||||
call print_bres(bres)
|
||||
call print_prof(params%profiles)
|
||||
call print_maps(bres, xgcn, &
|
||||
norm2(params%antenna%pos(1:2)) * 0.01_wp_, &
|
||||
sin(params%antenna%beta*degree))
|
||||
end if
|
||||
! ========= pre-proc prints END =========
|
||||
|
||||
! =========== main loop BEGIN ===========
|
||||
call initmultipass(params%raytracing%ipol, params%antenna%iox, &
|
||||
iroff,yynext,yypnext,yw0,ypw0, &
|
||||
@ -229,7 +226,8 @@ contains
|
||||
lgcpl1 = zero
|
||||
p0ray = p0jk ! * initial beam power
|
||||
|
||||
call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri, index_rt) ! * initial conditions
|
||||
call ic_gb(params, anv0, ak0, yw, ypw, stv, xc, du1, &
|
||||
gri, ggri, index_rt, results%tables) ! * initial conditions
|
||||
|
||||
do jk=1,params%raytracing%nray
|
||||
zzm = yw(3,jk)*0.01_wp_
|
||||
@ -259,7 +257,7 @@ contains
|
||||
iop = 0 ! start propagation outside plasma
|
||||
|
||||
if(params%raytracing%nray>1 .and. all(.not.iwait)) & ! iproj=0 ==> nfilp=8
|
||||
call print_projxyzt(stv,yw,0)
|
||||
call store_beam_shape(params%raytracing, results%tables, stv, yw)
|
||||
|
||||
! ======= propagation loop BEGIN =======
|
||||
call log_debug(' propagation loop start', mod='gray_core', proc='gray_main')
|
||||
@ -487,8 +485,9 @@ contains
|
||||
ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1)
|
||||
psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1)
|
||||
else
|
||||
call print_output(i,jk,stv(jk),p0ray(jk),xv,psinv, &
|
||||
btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
|
||||
call store_ray_data(results%tables, &
|
||||
i, jk, stv(jk), p0ray(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
|
||||
|
||||
@ -504,7 +503,7 @@ contains
|
||||
! print ray positions for j=nrayr in local reference system
|
||||
if(mod(i,params%output%istpr) == 0) then
|
||||
if(params%raytracing%nray > 1 .and. all(.not.iwait)) &
|
||||
call print_projxyzt(stv,yw,0)
|
||||
call store_beam_shape(params%raytracing, results%tables, stv, yw)
|
||||
end if
|
||||
|
||||
! test whether further trajectory integration is unnecessary
|
||||
@ -522,7 +521,9 @@ contains
|
||||
! ======== propagation loop END ========
|
||||
|
||||
! print all ray positions in local reference system
|
||||
if(params%raytracing%nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,1)
|
||||
if(params%raytracing%nray > 1 .and. all(.not.iwait)) &
|
||||
call store_beam_shape(params%raytracing, results%tables, &
|
||||
stv, yw, full=.true.)
|
||||
|
||||
! =========== post-proc BEGIN ===========
|
||||
! compute total absorbed power and driven current for current beam
|
||||
@ -581,17 +582,24 @@ contains
|
||||
end if
|
||||
|
||||
if (params%equilibrium%iequil /= EQ_VACUUM) then
|
||||
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
|
||||
call store_ec_profiles( &
|
||||
results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, &
|
||||
jcd_beam, dpdv_beam, currins_beam, pins_beam, ip)
|
||||
|
||||
call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam, &
|
||||
jphi_beam, rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, &
|
||||
rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam
|
||||
|
||||
call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav, &
|
||||
rhotjava,drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj, &
|
||||
drhotp,drhotj,ratjamx,ratjbmx,stv(1),psipv(index_rt), &
|
||||
chipv(index_rt),index_rt,sum(p0ray),cpl_beam1,cpl_beam2) ! *print 0D results for current beam
|
||||
if (results%tables%summary%active) then
|
||||
call results%tables%summary%append([ &
|
||||
wrap(icd_beam), wrap(pabs_beam), wrap(jphip), wrap(dpdvp), &
|
||||
wrap(rhotj), wrap(rhotjava), wrap(rhotp), wrap(rhotpav), &
|
||||
wrap(drhotjava), wrap(drhotpav), wrap(ratjamx), wrap(ratjbmx), &
|
||||
wrap(stv(1)), wrap(psipv(index_rt)), wrap(chipv(index_rt)), &
|
||||
wrap(index_rt), wrap(jphimx), wrap(dpdvmx), wrap(drhotj), &
|
||||
wrap(drhotp), wrap(sum(p0ray)), wrap(cpl_beam1), wrap(cpl_beam2)]) ! *print 0D results for current beam
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
! ============ post-proc END ============
|
||||
@ -621,10 +629,9 @@ contains
|
||||
call log_debug('pass loop end', mod='gray_core', proc='gray_main')
|
||||
! ============ main loop END ============
|
||||
|
||||
! ========== free memory BEGIN ==========
|
||||
! Free memory
|
||||
call dealloc_surfvec
|
||||
call dealloc_pec
|
||||
! =========== free memory END ===========
|
||||
|
||||
end block
|
||||
end associate
|
||||
@ -663,12 +670,15 @@ contains
|
||||
|
||||
|
||||
subroutine ic_gb(params, anv0c, ak0, ywrk0, ypwrk0, &
|
||||
stv, xc0, du10, gri, ggri, index_rt)
|
||||
stv, xc0, du10, gri, ggri, index_rt, &
|
||||
tables)
|
||||
! beam tracing initial conditions igrad=1
|
||||
! !!!!!! check ray tracing initial conditions igrad=0 !!!!!!
|
||||
use const_and_precisions, only : zero,one,pi,half,two,degree,ui=>im
|
||||
use gray_params, only : gray_parameters
|
||||
use gray_params, only : gray_parameters, gray_tables
|
||||
use beamdata, only : nray,nrayr,nrayth,rwmax
|
||||
use types, only : table
|
||||
use gray_tables, only : store_ray_data
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
@ -681,6 +691,9 @@ contains
|
||||
real(wp_), dimension(3, 3, nray), intent(out) :: ggri
|
||||
integer, intent(in) :: index_rt
|
||||
|
||||
! tables for storing initial rays data
|
||||
type(gray_tables), intent(inout), optional :: tables
|
||||
|
||||
! local variables
|
||||
real(wp_), dimension(3) :: xv0c
|
||||
real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir
|
||||
@ -973,10 +986,18 @@ contains
|
||||
|
||||
ddr = anx**2 + any**2 + anz**2 - an20
|
||||
ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt)
|
||||
! i=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0
|
||||
call print_output(0,jk,stv(jk),one,xc0(:,k,j),-one,zero,[zero,zero,zero], &
|
||||
ak0,zero,zero,[zero,zero,zero],zero,zero,zero,zero,zero,zero, &
|
||||
0,0,0,index_rt,ddr,ddi,zero,zero,[zero,zero,zero])
|
||||
|
||||
! save step "zero" data
|
||||
if (present(tables)) &
|
||||
call store_ray_data(tables, &
|
||||
i=0, jk=jk, s=stv(jk), P0=one, pos=xc0(:,k,j), &
|
||||
psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, &
|
||||
N_pl=zero, N_pr=zero, N=ywrk0(:, jk), N_pr_im=zero, &
|
||||
n_e=zero, T_e=zero, alpha=zero, tau=zero, dIds=zero, &
|
||||
nhm=0, nhf=0, iokhawa=0, index_rt=index_rt, &
|
||||
lambda_r=ddr, lambda_i=ddi, X=zero, Y=zero, &
|
||||
grad_X=[zero,zero,zero])
|
||||
|
||||
end do
|
||||
end subroutine ic_gb
|
||||
|
||||
@ -1848,700 +1869,4 @@ contains
|
||||
|
||||
end subroutine alpha_effj
|
||||
|
||||
|
||||
subroutine cniteq(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
|
||||
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
|
||||
! (based on an older code)
|
||||
! arguments
|
||||
integer, intent(in) :: nr,nz
|
||||
real(wp_), dimension(nr), intent(in) :: rqgrid
|
||||
real(wp_), dimension(nz), intent(in) :: zqgrid
|
||||
real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid
|
||||
real(wp_), intent(in) :: h
|
||||
integer, intent(inout) :: ncon, icount
|
||||
integer, dimension(ncon), intent(out) :: npts
|
||||
real(wp_), dimension(icount), intent(out) :: rcon,zcon
|
||||
! local variables
|
||||
integer :: i,j,k,l,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb
|
||||
integer :: jabs,jnb,kx,ikx,itm,inext,in
|
||||
integer, dimension(3,2) :: ja
|
||||
integer, dimension(icount/2-1) :: lx
|
||||
real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y
|
||||
real(wp_), dimension(nr*nz) :: a
|
||||
logical :: flag1
|
||||
|
||||
px = 0.5_wp_
|
||||
|
||||
a = reshape(matr2dgrid, [nr*nz])
|
||||
|
||||
rcon = 0.0_wp_
|
||||
zcon = 0.0_wp_
|
||||
|
||||
nrqmax = nr
|
||||
drgrd = rqgrid(2) - rqgrid(1)
|
||||
dzgrd = zqgrid(2) - zqgrid(1)
|
||||
|
||||
ncon = 0
|
||||
|
||||
npts = 0
|
||||
|
||||
iclast = 0
|
||||
icount = 0
|
||||
mpl = 0
|
||||
ix = 0
|
||||
mxr = nrqmax * (nz - 1)
|
||||
n1 = nr - 1
|
||||
|
||||
do jx=2,n1
|
||||
do jm=jx,mxr,nrqmax
|
||||
j = jm + nrqmax
|
||||
ah=a(j)-h
|
||||
if (ah <= 0.0_wp_ .and. a(jm) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(jm) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=-j
|
||||
end if
|
||||
if (ah <= 0.0_wp_ .and. a(j-1) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(j-1) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=j
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
do jm=nr,mxr,nrqmax
|
||||
j = jm + nrqmax
|
||||
ah=a(j)-h
|
||||
if (ah <= 0.0_wp_ .and. a(j-1) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(j-1) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=j
|
||||
end if
|
||||
if (ah <= 0.0_wp_ .and. a(jm) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(jm) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=-j
|
||||
end if
|
||||
end do
|
||||
|
||||
do jm=1,mxr,nrqmax
|
||||
j = jm + nrqmax
|
||||
if (a(j) <= h .and. a(jm) > h .or. &
|
||||
a(j) > h .and. a(jm) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix) =-j
|
||||
end if
|
||||
end do
|
||||
|
||||
do j=2,nr
|
||||
if (a(j) <= h .and. a(j-1) > h .or. &
|
||||
a(j) > h .and. a(j-1) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=j
|
||||
end if
|
||||
end do
|
||||
|
||||
if(ix<=0) return
|
||||
|
||||
bb: do
|
||||
in=ix
|
||||
jx=lx(in)
|
||||
jfor=0
|
||||
lda=1
|
||||
ldb=2
|
||||
|
||||
do
|
||||
if(jx<0) then
|
||||
jabs=-jx
|
||||
jnb = jabs - nrqmax
|
||||
else
|
||||
jabs=jx
|
||||
jnb=jabs-1
|
||||
end if
|
||||
|
||||
adn=a(jabs)-a(jnb)
|
||||
if(adn/=0) px=(a(jabs)-h)/adn
|
||||
kx = (jabs - 1) / nrqmax
|
||||
ikx = jabs - nrqmax * kx - 1
|
||||
|
||||
if(jx<0) then
|
||||
x = drgrd * ikx
|
||||
y = dzgrd * (kx - px)
|
||||
else
|
||||
x = drgrd * (ikx - px)
|
||||
y = dzgrd * kx
|
||||
end if
|
||||
|
||||
icount = icount + 1
|
||||
rcon(icount) = x + rqgrid(1)
|
||||
zcon(icount) = y + zqgrid(1)
|
||||
mpl= icount
|
||||
itm = 1
|
||||
ja(1,1) = jabs + nrqmax
|
||||
j=1
|
||||
|
||||
if(jx<=0) then
|
||||
ja(1,1) = -jabs-1
|
||||
j=2
|
||||
end if
|
||||
|
||||
ja(2,1) = -ja(1,1)
|
||||
ja(3,1) = -jx + 1 - nrqmax
|
||||
ja(3,2) = -jx
|
||||
ja(j,2) = jabs - nrqmax
|
||||
k= 3-j
|
||||
ja(k,2) = 1-jabs
|
||||
|
||||
if (kx<=0 .or. ikx<=0) then
|
||||
lda=1
|
||||
ldb=lda
|
||||
else if (ikx + 1 - nr >= 0 .and. jx <= 0) then
|
||||
lda=2
|
||||
ldb=lda
|
||||
else if(jfor/=0) then
|
||||
lda=2
|
||||
do i=1,3
|
||||
if(jfor==ja(i,2)) then
|
||||
lda=1
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
ldb=lda
|
||||
end if
|
||||
|
||||
flag1=.false.
|
||||
aa: do k=1,3
|
||||
do l=lda,ldb
|
||||
do i=1,ix
|
||||
if(lx(i)==ja(k,l)) then
|
||||
itm=itm+1
|
||||
inext= i
|
||||
if(jfor/=0) exit aa
|
||||
if(itm > 3) then
|
||||
flag1=.true.
|
||||
exit aa
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do aa
|
||||
|
||||
if(.not.flag1) then
|
||||
lx(in)=0
|
||||
if(itm == 1) exit
|
||||
end if
|
||||
|
||||
jfor=jx
|
||||
jx=lx(inext)
|
||||
in = inext
|
||||
end do
|
||||
|
||||
do
|
||||
if(lx(ix)/=0) then
|
||||
if(mpl>=4) then
|
||||
ncon = ncon + 1
|
||||
npts(ncon) = icount - iclast
|
||||
iclast = icount
|
||||
end if
|
||||
exit
|
||||
end if
|
||||
ix= ix-1
|
||||
if(ix<=0) exit bb
|
||||
end do
|
||||
|
||||
end do bb
|
||||
|
||||
if(mpl >= 4) then
|
||||
ncon = ncon + 1
|
||||
npts(ncon) = icount - iclast
|
||||
iclast = icount
|
||||
end if
|
||||
end subroutine cniteq
|
||||
|
||||
|
||||
|
||||
subroutine print_headers(params)
|
||||
! Prints the headers of all gray_main output tables
|
||||
|
||||
use units, only : uprj0, uwbm, udisp, ucenr, uoutr, upec, usumm, &
|
||||
unit_active, active_units
|
||||
use gray_params, only : gray_parameters, headw, headl, print_parameters
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
|
||||
! local variables
|
||||
integer :: i, j
|
||||
integer :: main_units(8)
|
||||
character(256) :: main_headers(8)
|
||||
character(len=headw), dimension(headl) :: header
|
||||
|
||||
main_units = [uprj0, uprj0+1, uwbm, udisp, ucenr, uoutr, upec, usumm]
|
||||
|
||||
main_headers(1) = '#sst j k xt yt zt rt'
|
||||
main_headers(2) = '#sst j k xt yt zt rt'
|
||||
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'
|
||||
|
||||
if (all(active_units == 0)) return
|
||||
|
||||
call print_parameters(params, header)
|
||||
|
||||
do i=1,size(main_units)
|
||||
if (unit_active(main_units(i))) then
|
||||
do j=1,size(header)
|
||||
write (main_units(i), '(1x,a)') header(j)
|
||||
end do
|
||||
write (main_units(i), '(1x,a)') trim(main_headers(i))
|
||||
end if
|
||||
end do
|
||||
|
||||
end subroutine print_headers
|
||||
|
||||
|
||||
subroutine print_prof(params)
|
||||
! Prints the (input) plasma kinetic profiles (unit 55)
|
||||
|
||||
use gray_params, only : profiles_parameters
|
||||
use equilibrium, only : fq, frhotor
|
||||
use coreprofiles, only : density, temp
|
||||
use units, only : uprfin, unit_active
|
||||
use magsurf_data, only : npsi, vajphiav
|
||||
|
||||
! suborutine arguments
|
||||
type(profiles_parameters), intent(in) :: params
|
||||
|
||||
! local variables
|
||||
integer :: i, ntail
|
||||
real(wp_) :: rho_t, J_phi, dens, ddens
|
||||
real(wp_) :: psi_n, rho_p, drho_p
|
||||
|
||||
if (.not. unit_active(uprfin)) return
|
||||
|
||||
! Δρ_p for the uniform ρ_p grid
|
||||
! Note: we don't use ψ_n because J_phi has been
|
||||
! sampled uniformly in ρ_p (see flux_average)
|
||||
drho_p = 1.0_wp_ / (npsi - 1)
|
||||
|
||||
! extra points to reach ψ=ψ_bnd
|
||||
ntail = ceiling((sqrt(params%psnbnd) - 1) / drho_p)
|
||||
|
||||
write (uprfin, *) '#psi rhot ne Te q Jphi'
|
||||
do i = 0, npsi + ntail
|
||||
rho_p = i * drho_p
|
||||
rho_t = frhotor(rho_p)
|
||||
psi_n = rho_p**2
|
||||
if (psi_n < 1) then
|
||||
J_phi = vajphiav(i+1) * 1.e-6_wp_
|
||||
else
|
||||
J_phi = 0
|
||||
end if
|
||||
call density(psi_n, dens, ddens)
|
||||
write (uprfin, '(12(1x,e12.5))') &
|
||||
psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi
|
||||
end do
|
||||
end subroutine print_prof
|
||||
|
||||
|
||||
subroutine print_bres(bres)
|
||||
! Prints the EC resonance surface table (unit 70)
|
||||
|
||||
use const_and_precisions, only : comp_eps
|
||||
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield
|
||||
use units, only : ubres, unit_active
|
||||
use magsurf_data, only : npsi
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: bres
|
||||
|
||||
! local constants
|
||||
integer, parameter :: icmx = 2002
|
||||
|
||||
! local variables
|
||||
integer :: j,k,n,nconts,nctot
|
||||
integer, dimension(10) :: ncpts
|
||||
real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
|
||||
real(wp_), dimension(icmx) :: rrcb,zzcb
|
||||
real(wp_) :: rv(npsi), zv(npsi)
|
||||
real(wp_), dimension(npsi,npsi) :: btotal
|
||||
|
||||
if (.not. unit_active(ubres)) return
|
||||
|
||||
! Build a regular (R, z) grid
|
||||
dr = (rmxm - rmnm - comp_eps)/(npsi - 1)
|
||||
dz = (zmxm - zmnm)/(npsi - 1)
|
||||
do j=1,npsi
|
||||
rv(j) = comp_eps + rmnm + dr*(j - 1)
|
||||
zv(j) = zmnm + dz*(j - 1)
|
||||
end do
|
||||
|
||||
! Btotal on psi grid
|
||||
btmx = -1.0e30_wp_
|
||||
btmn = 1.0e30_wp_
|
||||
do k = 1, npsi
|
||||
zzk = zv(k)
|
||||
do j = 1, npsi
|
||||
rrj = rv(j)
|
||||
call bfield(rrj, zzk, bbr, bbz, bbphi)
|
||||
btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2)
|
||||
if(btotal(j,k) >= btmx) btmx = btotal(j,k)
|
||||
if(btotal(j,k) <= btmn) btmn = btotal(j,k)
|
||||
end do
|
||||
end do
|
||||
|
||||
! compute Btot=Bres/n with n=1,5
|
||||
write (ubres, *) '#i Btot R z'
|
||||
do n = 1, 5
|
||||
bbb = bres/n
|
||||
if (bbb >= btmn .and. bbb <= btmx) then
|
||||
nconts = size(ncpts)
|
||||
nctot = size(rrcb)
|
||||
call cniteq(rv, zv, btotal, npsi, npsi, bbb, &
|
||||
nconts, ncpts, nctot, rrcb, zzcb)
|
||||
do j = 1, nctot
|
||||
write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j)
|
||||
end do
|
||||
end if
|
||||
write (ubres, *)
|
||||
end do
|
||||
|
||||
end subroutine print_bres
|
||||
|
||||
|
||||
subroutine print_maps(B_res, xgcn, R0, Npl0)
|
||||
! Prints several input quantities on a regular grid (unit 72)
|
||||
|
||||
use const_and_precisions, only : comp_eps
|
||||
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
|
||||
use coreprofiles, only : density, temp
|
||||
use units, only : umaps, unit_active
|
||||
use magsurf_data, only : npsi
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω
|
||||
real(wp_), intent(in) :: xgcn ! X normalisation, e²/ε₀m_eω²
|
||||
real(wp_), intent(in) :: R0 ! initial value of R
|
||||
real(wp_), intent(in) :: Npl0 ! initial value of N∥
|
||||
|
||||
! local variables
|
||||
integer :: j, k
|
||||
real(wp_) :: dR, dz, B_R, B_phi, B_z, B
|
||||
real(wp_) :: psi_n, ne, dne, Te, X, Y, Npl
|
||||
real(wp_), dimension(npsi) :: R, z
|
||||
|
||||
if (.not. unit_active(umaps)) return
|
||||
|
||||
! Build a regular (R, z) grid
|
||||
dR = (rmxm - rmnm - comp_eps)/(npsi - 1)
|
||||
dz = (zmxm - zmnm)/(npsi - 1)
|
||||
do j = 1, npsi
|
||||
R(j) = comp_eps + rmnm + dR*(j - 1)
|
||||
z(j) = zmnm + dz*(j - 1)
|
||||
end do
|
||||
|
||||
write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl'
|
||||
do j = 1, npsi
|
||||
Npl = Npl0 * R0/r(j)
|
||||
do k = 1, npsi
|
||||
call pol_flux(r(j), z(k), psi_n)
|
||||
call bfield(r(j), z(k), B_R, B_z, B_phi)
|
||||
call density(psi_n, ne, dne)
|
||||
B = sqrt(B_R**2 + B_phi**2 + B_z**2)
|
||||
X = xgcn*ne
|
||||
Y = B/B_res
|
||||
Te = temp(psi_n)
|
||||
write (umaps,'(12(x,e12.5))') &
|
||||
R(j), z(k), psi_n, B_R, B_phi, B_z, B, ne, Te, X, Y, Npl
|
||||
end do
|
||||
write (umaps,*)
|
||||
end do
|
||||
|
||||
end subroutine print_maps
|
||||
|
||||
|
||||
subroutine print_surfq(qvals)
|
||||
! Prints the constant ψ surfaces for a set of values of q
|
||||
|
||||
use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf
|
||||
use magsurf_data, only : contours_psi, npoints, print_contour
|
||||
use logger, only : log_info
|
||||
use minpack, only : hybrj1
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: qvals(:)
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
real(wp_) :: rho_t, rho_p, psi_n
|
||||
character(256) :: msg ! for log messages formatting
|
||||
|
||||
call log_info('constant ψ surfaces for:', &
|
||||
mod='gray_core', proc='print_surfq')
|
||||
do i = 1, size(qvals)
|
||||
! Find value of ψ_n for the current q
|
||||
block
|
||||
real(wp_) :: sol(1), fvec(1), fjac(1,1), wa(7)
|
||||
integer :: info
|
||||
|
||||
! Note: since we are mostly interested in q=3/2,2 we start
|
||||
! searching near the boundary in case q is not monotonic.
|
||||
sol = [0.8_wp_] ! first guess
|
||||
|
||||
! Solve fq(ψ_n) = qvals(i) for ψ_n
|
||||
call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, &
|
||||
ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7)
|
||||
! Solution
|
||||
psi_n = sol(1)
|
||||
|
||||
! Handle spurious or no solutions
|
||||
if (info /= 1 .or. psi_n < 0 .or. psi_n > 1) cycle
|
||||
end block
|
||||
|
||||
! Compute and print the ψ_n(R,z) = ψ_n₀ contour
|
||||
block
|
||||
real(wp_), dimension(npoints) :: R_cont, z_cont
|
||||
real(wp_) :: R_hi, R_lo, z_hi, z_lo
|
||||
|
||||
! Guesses for the high,low horzizontal-tangent points
|
||||
R_hi = rmaxis;
|
||||
z_hi = (zbsup + zmaxis)/2
|
||||
R_lo = rmaxis
|
||||
z_lo = (zbinf + zmaxis)/2
|
||||
|
||||
call contours_psi(psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo)
|
||||
call print_contour(psi_n, R_cont, z_cont)
|
||||
end block
|
||||
|
||||
! Log the values found for ψ_n, ρ_p, ρ_t
|
||||
rho_p = sqrt(psi_n)
|
||||
rho_t = frhotor(rho_p)
|
||||
write (msg, '(4(x,a,"=",g0.3))') &
|
||||
'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t
|
||||
call log_info(msg, mod='gray_core', proc='print_surfq')
|
||||
end do
|
||||
|
||||
contains
|
||||
|
||||
subroutine equation(n, x, f, df, ldf, flag)
|
||||
! The equation to solve: f(x) = q(x) - q₀ = 0
|
||||
use const_and_precisions, only : comp_eps
|
||||
|
||||
! subroutine arguments
|
||||
integer, intent(in) :: n, ldf, flag
|
||||
real(wp_), intent(in) :: x(n)
|
||||
real(wp_), intent(inout) :: f(n), df(ldf,n)
|
||||
|
||||
! optimal step size
|
||||
real(wp_), parameter :: e = comp_eps**(1/3.0_wp_)
|
||||
|
||||
if (flag == 1) then
|
||||
! return f(x)
|
||||
f(1) = fq(x(1)) - qvals(i)
|
||||
else
|
||||
! return f'(x), computed numerically
|
||||
df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e)
|
||||
end if
|
||||
end subroutine
|
||||
|
||||
end subroutine print_surfq
|
||||
|
||||
|
||||
subroutine print_projxyzt(stv, ywrk, iproj)
|
||||
! 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 units, only : uprj0, uwbm, unit_active
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), dimension(:), intent(in) :: stv
|
||||
real(wp_), dimension(:,:), intent(in) :: ywrk
|
||||
integer, intent(in) :: iproj
|
||||
|
||||
! local variables
|
||||
integer :: jk, jkz, uprj
|
||||
integer, dimension(2) ::jkv
|
||||
real(wp_), dimension(3) :: xv1, dir, dxv
|
||||
real(wp_) :: rtimn, rtimx, csth1, snth1, &
|
||||
csps1, snps1, xti, yti, zti, rti
|
||||
|
||||
if (.not.(unit_active(uprj0) .or. unit_active(uwbm))) return
|
||||
|
||||
! common/external functions/variables
|
||||
|
||||
uprj = uprj0 + iproj
|
||||
|
||||
xv1 = ywrk(1:3, 1)
|
||||
dir = ywrk(4:6, 1)
|
||||
dir = dir/norm2(dir)
|
||||
csth1 = dir(3)
|
||||
snth1 = sqrt(one - csth1**2)
|
||||
if(snth1 > zero) then
|
||||
csps1 = dir(2)/snth1
|
||||
snps1 = dir(1)/snth1
|
||||
else
|
||||
csps1 = one
|
||||
snps1 = zero
|
||||
end if
|
||||
|
||||
if(iproj==0) then
|
||||
jkz = nray - nrayth + 1
|
||||
else
|
||||
jkz = 1
|
||||
end if
|
||||
|
||||
rtimn = comp_huge
|
||||
rtimx = zero
|
||||
do jk = jkz, nray
|
||||
dxv = ywrk(1:3,jk) - xv1
|
||||
xti = dxv(1)*csps1 - dxv(2)*snps1
|
||||
yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1
|
||||
zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1
|
||||
rti = hypot(xti, yti)
|
||||
|
||||
jkv=rayi2jk(jk)
|
||||
if(.not.(iproj==0 .and. jk==1)) &
|
||||
write(uprj, '(1x,e16.8e3,2i5,4(1x,e16.8e3))') &
|
||||
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<=rtimn .and. jkv(1)==nrayr) rtimn = rti
|
||||
end do
|
||||
write (uprj, *)
|
||||
write (uwbm, '(3(1x,e16.8e3))') stv(1), rtimn, rtimx
|
||||
|
||||
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)
|
||||
|
||||
use const_and_precisions, only : degree,zero,one
|
||||
use equilibrium, only : frhotor
|
||||
use gray_params, only : istpl0
|
||||
use beamdata, only : nray,nrayth,jkray1
|
||||
use units, only : ucenr, uoutr, udisp, unit_active
|
||||
|
||||
! subroutine arguments
|
||||
integer, intent(in) :: i, jk, nhm, nhf, iokhawa, index_rt
|
||||
real(wp_), dimension(3), intent(in) :: xv, bv, anv
|
||||
real(wp_), intent(in) :: st, qj, psinv, btot, ak0, anpl, anpr, anprim
|
||||
real(wp_), intent(in) :: dens, tekev, alpha, tau, dids, ddr, ddi
|
||||
real(wp_), intent(in) :: xg, yg
|
||||
real(wp_), dimension(3), intent(in) :: derxg
|
||||
|
||||
! local variables
|
||||
real(wp_) :: stm, xxm, yym, zzm, rrm, phideg, rhot, akim, pt, didsn
|
||||
integer :: k
|
||||
|
||||
stm = st * 1.0e-2_wp_
|
||||
xxm = xv(1) * 1.0e-2_wp_
|
||||
yym = xv(2) * 1.0e-2_wp_
|
||||
zzm = xv(3) * 1.0e-2_wp_
|
||||
rrm = hypot(xxm, yym)
|
||||
|
||||
! print central ray trajectory. dIds in A/m/W, ki in m⁻¹
|
||||
if(jk == 1) then
|
||||
phideg = atan2(yym, xxm)/degree
|
||||
if(psinv>=zero .and. psinv<=one) then
|
||||
rhot = frhotor(sqrt(psinv))
|
||||
else
|
||||
rhot = 1.0_wp_
|
||||
end if
|
||||
akim = anprim * ak0 * 1.0e2_wp_
|
||||
pt = qj*exp(-tau)
|
||||
didsn = dids * 1.0e2_wp_/qj
|
||||
|
||||
! print central ray data
|
||||
if (unit_active(ucenr)) then
|
||||
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
|
||||
|
||||
! print conservation of dispersion relation (Hamiltonian error)
|
||||
if(unit_active(udisp) .and. jk==nray) then
|
||||
write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi
|
||||
end if
|
||||
|
||||
! print outer ray trajectories
|
||||
if(mod(i, istpl0) == 0) then
|
||||
k = jk - jkray1 + 1
|
||||
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, &
|
||||
psinv, tau, anpl, alpha, index_rt
|
||||
end if
|
||||
end if
|
||||
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 arguments
|
||||
real(wp_), dimension(:), intent(in) :: rhop_tab, rhot_tab, jphi, jcd, &
|
||||
dpdv, currins, pins
|
||||
integer, intent(in) :: index_rt
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
|
||||
if (.not. unit_active(upec)) return
|
||||
|
||||
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
|
||||
write(upec, *) ''
|
||||
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 arguments
|
||||
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
|
||||
|
||||
if (unit_active(ucenr)) then
|
||||
write(ucenr, *) ''
|
||||
end if
|
||||
|
||||
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 module gray_core
|
||||
|
||||
@ -2,7 +2,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam, &
|
||||
p0mw, alphain, betain, dpdv, jcd, pabs, icd, err)
|
||||
use const_and_precisions, only: wp_
|
||||
use units, only: set_active_units, close_units
|
||||
use gray_params, only: gray_parameters, gray_data, gray_results
|
||||
use gray_core, only: gray_main
|
||||
|
||||
@ -29,16 +28,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
first_call: if (firstcall) then
|
||||
firstcall = .false.
|
||||
|
||||
! Activate some debug units
|
||||
set_units: block
|
||||
use units
|
||||
call set_active_units([ &
|
||||
ucenr, usumm, uprj0, uprj0+1, &
|
||||
uwbm, udisp, ubres, ucnt, &
|
||||
uoutr, umaps, uprfin, uflx, &
|
||||
upec])
|
||||
end block set_units
|
||||
|
||||
! Read parameters from external file
|
||||
init_params: block
|
||||
use gray_params, only : read_gray_params, set_globals
|
||||
@ -47,6 +36,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
|
||||
! Override some parameters
|
||||
params%misc%rwall = r(1)
|
||||
params%misc%active_tables = [4, 7, 8, 9, 12, 17, 33, 70, 71, 55, 48]
|
||||
params%antenna%filenm = 'graybeam.data'
|
||||
params%equilibrium%filenm = 'JETTO'
|
||||
params%equilibrium%iequil = ijetto + 1
|
||||
@ -134,6 +124,16 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
! Call main subroutine for the ibeam-th beam
|
||||
call gray_main(params, data, res, err, rhout=sqrt(psrad))
|
||||
|
||||
write_debug_files: block
|
||||
integer :: i, err
|
||||
! Write results to file
|
||||
do i = 1, size(res%tables%iter)
|
||||
associate (tbl => res%tables%iter(i)%ptr)
|
||||
call tbl%save(err, filepath='gray_'//trim(tbl%title)//'.txt')
|
||||
end associate
|
||||
end do
|
||||
end block write_debug_files
|
||||
|
||||
! Free memory
|
||||
free_memory: block
|
||||
use equilibrium, only : unset_equil_spline
|
||||
@ -152,6 +152,4 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
|
||||
dpdv = res%dpdv
|
||||
jcd = res%jcd
|
||||
|
||||
call close_units
|
||||
|
||||
end subroutine gray_jetto1beam
|
||||
|
||||
@ -1,6 +1,7 @@
|
||||
module gray_params
|
||||
|
||||
use const_and_precisions, only : wp_
|
||||
use types, only : table
|
||||
|
||||
implicit none
|
||||
|
||||
@ -157,13 +158,14 @@ module gray_params
|
||||
type output_parameters
|
||||
integer :: ipec = RHO_TOR ! Radial coordinate of the EC profiles
|
||||
integer :: nrho = 501 ! Number of grid steps for EC profiles + 1
|
||||
integer :: istpr = 5 ! Subsampling factor for beam cross-section (units 8, 12)
|
||||
integer :: istpl = 5 ! Subsampling factor for outer rays (unit 33)
|
||||
integer :: istpr = 5 ! Subsampling factor for beam cross-section (tables 8, 12)
|
||||
integer :: istpl = 5 ! Subsampling factor for outer rays (table 33)
|
||||
end type
|
||||
|
||||
! Other parameters
|
||||
type misc_parameters
|
||||
real(wp_) :: rwall ! R of the limiter (fallback)
|
||||
integer, allocatable :: active_tables(:) ! IDs of output tables to fill in
|
||||
end type
|
||||
|
||||
! MHD equilibrium data
|
||||
@ -228,14 +230,33 @@ module gray_params
|
||||
type(profiles_data) :: profiles
|
||||
end type
|
||||
|
||||
! GRAY final results (only some)
|
||||
! Wrapper type for array of pointers
|
||||
type table_ptr
|
||||
type(table), pointer :: ptr => null()
|
||||
end type
|
||||
|
||||
! Extra outputs tables
|
||||
type gray_tables
|
||||
type(table) :: flux_averages, flux_surfaces
|
||||
type(table) :: summary, ec_profiles
|
||||
type(table) :: central_ray, outer_rays, dispersion
|
||||
type(table) :: beam_shape, beam_final, beam_size
|
||||
type(table) :: kinetic_profiles, ec_resonance, inputs_maps
|
||||
|
||||
! for iterating over the tables
|
||||
type(table_ptr) :: iter(13)
|
||||
end type
|
||||
|
||||
! GRAY final results
|
||||
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
|
||||
type(gray_tables) :: tables ! Extra outputs
|
||||
end type
|
||||
|
||||
|
||||
! Global variables exposed for gray_core
|
||||
integer, save :: iequil, iprof
|
||||
integer, save :: iwarm, ilarm, imx, ieccd
|
||||
@ -245,11 +266,12 @@ module gray_params
|
||||
|
||||
contains
|
||||
|
||||
subroutine print_parameters(params, strout)
|
||||
subroutine write_gray_header(params, file)
|
||||
! Writes the Gray output header to a file
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
character(len=*), dimension(:), intent(out) :: strout ! min len=110, dimension(21)
|
||||
integer :: file
|
||||
|
||||
! local variables
|
||||
character(len=8) :: date
|
||||
@ -260,91 +282,90 @@ contains
|
||||
|
||||
! Date and time
|
||||
call date_and_time(date, time)
|
||||
write(strout(1), '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') &
|
||||
write(file, '("# Run date/time: ",a4,2("/",a2),1x,2(a2,":"),a6)') &
|
||||
date(1:4), date(5:6), date(7:8), &
|
||||
time(1:2), time(3:4), time(5:10)
|
||||
|
||||
! Git revision
|
||||
write(strout(2), '("# GRAY Git revision: ",a)') REVISION
|
||||
write(file, '("# GRAY Git revision: ",a)') REVISION
|
||||
|
||||
! MHD equilibrium parameters
|
||||
if (params%equilibrium%iequil > 0) then
|
||||
write(strout(3), '("# EQL input: ",a)') trim(params%equilibrium%filenm(1:95))
|
||||
write(file, '("# EQL input: ",a)') trim(params%equilibrium%filenm)
|
||||
! TODO add missing values
|
||||
write(strout(7), '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') &
|
||||
0._wp_, 0._wp_, 0._wp_, 0._wp_, 0._wp_
|
||||
write(file, '("# EQL B0 R0 aminor Rax zax:",5(1x,g0.5))') &
|
||||
0, 0, 0, 0, 0
|
||||
else
|
||||
write(strout(3), '("# EQL input: N/A (vacuum)")')
|
||||
write(strout(7), '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
|
||||
write(file, '("# EQL input: N/A (vacuum)")')
|
||||
write(file, '("# EQL B0 R0 aminor Rax zax: N/A (vacuum)")')
|
||||
end if
|
||||
|
||||
write(strout(4), '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') &
|
||||
write(file, '("# EQL iequil sgnb sgni factb:",3(1x,g0),1x,g0.3)') &
|
||||
params%equilibrium%iequil, params%equilibrium%sgnb, &
|
||||
params%equilibrium%sgni, params%equilibrium%factb
|
||||
if (params%equilibrium%iequil > 1) then
|
||||
write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') &
|
||||
write(file, '("# EQL icocos ipsinorm idesc ifreefmt:",4(1x,g0))') &
|
||||
params%equilibrium%icocos, params%equilibrium%ipsinorm, &
|
||||
params%equilibrium%idesc, params%equilibrium%ifreefmt
|
||||
write(strout(6), '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') &
|
||||
write(file, '("# EQL ssplps ssplf ixp:",2(1x,g8.3e1),1x,g0)') &
|
||||
params%equilibrium%ssplps, params%equilibrium%ssplf, &
|
||||
params%equilibrium%ixp
|
||||
else
|
||||
write(strout(5), '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
|
||||
write(strout(6), '("# EQL ssplps ssplf ixp: N/A (analytical)")')
|
||||
write(file, '("# EQL icocos ipsinorm idesc ifreefmt: N/A (analytical)")')
|
||||
write(file, '("# EQL ssplps ssplf ixp: N/A (analytical)")')
|
||||
end if
|
||||
|
||||
! Profiles parameters
|
||||
if (params%equilibrium%iequil > 0) then
|
||||
write(strout(8), '("# PRF input: ",a)') trim(params%profiles%filenm(1:95))
|
||||
write(strout(9), '("# PRF iprof iscal factne factte:",4(1x,g0.4))') &
|
||||
write(file, '("# PRF input: ",a)') trim(params%profiles%filenm)
|
||||
write(file, '("# PRF iprof iscal factne factte:",4(1x,g0.4))') &
|
||||
params%profiles%iprof, params%profiles%iscal, &
|
||||
params%profiles%factne, params%profiles%factte
|
||||
if (params%profiles%iprof > 0) then
|
||||
write(strout(10), '("# PRF irho psnbnd sspld:",3(1x,g0.3))') &
|
||||
write(file, '("# PRF irho psnbnd sspld:",3(1x,g0.3))') &
|
||||
params%profiles%irho, params%profiles%psnbnd, params%profiles%sspld
|
||||
else
|
||||
write(strout(10), '("# PRF irho psnbnd sspld: N/A (analytical)")')
|
||||
write(file, '("# PRF irho psnbnd sspld: N/A (analytical)")')
|
||||
end if
|
||||
! TODO: add missing values
|
||||
write(strout(11), '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') &
|
||||
0._wp_, 0._wp_, 0._wp_
|
||||
write(file, '("# PRF Te0 ne0 Zeff0:",3(1x,g0.4))') 0, 0, 0
|
||||
else
|
||||
write(strout(8), '("# PRF input: 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(11), '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")')
|
||||
write(file, '("# PRF input: N/A (vacuum)")')
|
||||
write(file, '("# PRF iprof iscal factne factte: N/A (vacuum)")')
|
||||
write(file, '("# PRF irho psnbnd sspld: N/A (vacuum)")')
|
||||
write(file, '("# PRF Te0 ne0 Zeff0: N/A (vacuum)")')
|
||||
end if
|
||||
|
||||
! Antenna parameters
|
||||
write(strout(12), '("# ANT input: ",a)') trim(params%antenna%filenm(1:95))
|
||||
write(strout(13), '("# ANT ibeam iox psi chi:",4(1x,g0.4))') &
|
||||
write(file, '("# ANT input: ",a)') trim(params%antenna%filenm)
|
||||
write(file, '("# ANT ibeam iox psi chi:",4(1x,g0.4))') &
|
||||
params%antenna%ibeam, params%antenna%iox, &
|
||||
params%antenna%psi, params%antenna%chi
|
||||
write(strout(14), '("# ANT alpha beta power fghz:",4(1x,g0.3))') &
|
||||
write(file, '("# ANT alpha beta power fghz:",4(1x,g0.3))') &
|
||||
params%antenna%alpha, params%antenna%beta, &
|
||||
params%antenna%power, params%antenna%fghz
|
||||
write(strout(15), '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos
|
||||
write(strout(16), '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') &
|
||||
write(file, '("# ANT x0 y0 z0:",3(1x,g0.3))') params%antenna%pos
|
||||
write(file, '("# ANT wx wy Rcix Rciy psiw psir:",6(1x,g15.5e1))') &
|
||||
params%antenna%w, params%antenna%ri, params%antenna%phi
|
||||
|
||||
! Other parameters
|
||||
write(strout(17), '("# RFL rwall:",1x,e12.5)') params%misc%rwall
|
||||
write(file, '("# RFL rwall:",1x,e12.5)') params%misc%rwall
|
||||
|
||||
! code parameters
|
||||
write(strout(18), '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') &
|
||||
write(file, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') &
|
||||
params%raytracing%igrad, params%raytracing%idst, &
|
||||
params%raytracing%ipass, params%raytracing%ipol
|
||||
write(strout(19), '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') &
|
||||
write(file, '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') &
|
||||
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(file, '("# COD iwarm ilarm imx ieccd:",4(1x,i4))') &
|
||||
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(file, '("# COD ipec nrho istpr istpl:",4(1x,i4))') &
|
||||
params%output%ipec, params%output%nrho, &
|
||||
params%output%istpr, params%output%istpl
|
||||
end subroutine print_parameters
|
||||
end subroutine write_gray_header
|
||||
|
||||
|
||||
function update_parameter(params, name, value) result(err)
|
||||
|
||||
612
src/gray_tables.f90
Normal file
612
src/gray_tables.f90
Normal file
@ -0,0 +1,612 @@
|
||||
! This module contains the gray ouput tables and routines to operate on them
|
||||
module gray_tables
|
||||
|
||||
use const_and_precisions, only : wp_
|
||||
use types, only : queue, table, wrap
|
||||
|
||||
implicit none
|
||||
|
||||
private
|
||||
public init_tables ! Initialising main tables
|
||||
public find_flux_surfaces, store_flux_surface ! Filling the main tables
|
||||
public store_ray_data, store_ec_profiles !
|
||||
public store_beam_shape !
|
||||
public kinetic_profiles, ec_resonance, inputs_maps ! Extra tables
|
||||
|
||||
contains
|
||||
|
||||
pure function is_active(params, id)
|
||||
! Helper to check whether a table is active
|
||||
use gray_params, only : gray_parameters
|
||||
|
||||
! function arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
integer, intent(in) :: id
|
||||
logical :: is_active
|
||||
|
||||
associate (ids => params%misc%active_tables)
|
||||
is_active = any(ids == id)
|
||||
if (size(ids) == 1) then
|
||||
! special case, all enabled
|
||||
is_active = is_active .or. (ids(1) == -1)
|
||||
end if
|
||||
end associate
|
||||
end function
|
||||
|
||||
|
||||
subroutine init_tables(params, tbls)
|
||||
! Initialises the gray_main output tables
|
||||
|
||||
use gray_params, only : gray_parameters, gray_tables, ptr=>table_ptr
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
type(gray_tables), target, intent(out) :: tbls
|
||||
|
||||
call tbls%flux_averages%init( &
|
||||
title='flux-averages', id=56, active=is_active(params, 56), &
|
||||
labels=[character(64) :: 'ρ_p', 'ρ_t', 'B_avg', &
|
||||
'B_max', 'B_min', 'area', 'vol', 'I_pl', &
|
||||
'J_φ_avg', 'fc', 'ratio_Ja', 'ratio_Jb', 'q'])
|
||||
|
||||
call tbls%flux_surfaces%init( &
|
||||
title='flux-surfaces', id=71, active=is_active(params, 71), &
|
||||
labels=[character(64) :: 'i', 'ψ_n', 'R', 'z'])
|
||||
|
||||
call tbls%ec_profiles%init( &
|
||||
title='ec-profiles', id=48, active=is_active(params, 48), &
|
||||
labels=[character(64) :: 'ρ_p', 'ρ_t', 'J_φ', &
|
||||
'J_cdb', 'dPdV', 'I_cd_inside', &
|
||||
'P_inside', 'index_rt'])
|
||||
|
||||
call tbls%summary%init( &
|
||||
title='summary', id=7, active=is_active(params, 7), &
|
||||
labels=[character(64) :: 'I_cd', 'P_abs', 'J_φ_peak', &
|
||||
'dPdV_peak', 'ρ_max_J', 'ρ_avg_J', 'ρ_max_P', &
|
||||
'ρ_avg_P', 'Δρ_avg_J', 'Δρ_avg_P', &
|
||||
'ratio_Ja_max', 'ratio_Jb_max', 's_max', &
|
||||
'ψ', 'χ', 'index_rt', 'J_φ_max', 'dPdV_max', &
|
||||
'Δρ_J', 'Δρ_P', 'P0', 'cpl1', 'cpl2'])
|
||||
|
||||
call tbls%central_ray%init( &
|
||||
title='central-ray', id=4, active=is_active(params, 4), &
|
||||
labels=[character(64) :: 's', 'R', 'z', 'φ', &
|
||||
'ψ_n', 'ρ_t', 'n_e', 'T_e', 'B', 'b_x', &
|
||||
'b_y', 'b_z', 'N_⊥', 'N_∥', 'N_x', 'N_y', &
|
||||
'N_z', 'k_im', 'α', 'τ', 'P_t', 'dI/ds', &
|
||||
'n_harm_min', 'n_harm_max', 'i_okhawa', &
|
||||
'index_rt', 'Λ_r', 'X', 'Y', &
|
||||
'∇X_x', '∇X_y', '∇X_z'])
|
||||
|
||||
call tbls%outer_rays%init( &
|
||||
title='outer-rays', id=33, active=is_active(params, 33), &
|
||||
labels=[character(64) :: 'i', 'k', 's', 'x', 'y', &
|
||||
'R', 'z', 'ψ_n', 'τ', 'N_∥', 'α', 'index_rt'])
|
||||
|
||||
call tbls%dispersion%init( &
|
||||
title='dispersion', id=17, active=is_active(params, 17), &
|
||||
labels=[character(64) :: 's', 'Λ_r', 'Λ_i'])
|
||||
|
||||
call tbls%beam_shape%init( &
|
||||
title='beam-shape', id=8, active=is_active(params, 8), &
|
||||
labels=['s', 'j', 'k', 'x', 'y', 'z', 'r'])
|
||||
|
||||
call tbls%beam_final%init( &
|
||||
title='beam-shape-final', id=9, active=is_active(params, 9), &
|
||||
labels=['s', 'j', 'k', 'x', 'y', 'z', 'r'])
|
||||
|
||||
call tbls%beam_size%init( &
|
||||
title='beam-size', id=12, active=is_active(params, 12), &
|
||||
labels=[character(64) :: 's', 'r_min', 'r_max'])
|
||||
|
||||
! List of pointers to the tables for iterating over
|
||||
tbls%iter = [ &
|
||||
ptr(tbls%flux_averages), ptr(tbls%flux_surfaces), &
|
||||
ptr(tbls%summary), ptr(tbls%ec_profiles), &
|
||||
ptr(tbls%central_ray), ptr(tbls%outer_rays), ptr(tbls%dispersion), &
|
||||
ptr(tbls%beam_shape), ptr(tbls%beam_final), ptr(tbls%beam_size), &
|
||||
ptr(tbls%kinetic_profiles), ptr(tbls%ec_resonance), ptr(tbls%inputs_maps)]
|
||||
|
||||
end subroutine init_tables
|
||||
|
||||
|
||||
function kinetic_profiles(params) result(tbl)
|
||||
! Generates the plasma kinetic profiles
|
||||
|
||||
use gray_params, only : gray_parameters, EQ_VACUUM
|
||||
use equilibrium, only : fq, frhotor
|
||||
use coreprofiles, only : density, temp
|
||||
use magsurf_data, only : npsi, vajphiav
|
||||
|
||||
! function arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
type(table) :: tbl
|
||||
|
||||
! local variables
|
||||
integer :: i, ntail
|
||||
|
||||
real(wp_) :: rho_t, J_phi, dens, ddens
|
||||
real(wp_) :: psi_n, rho_p, drho_p
|
||||
|
||||
call tbl%init(title='kinetic-profiles', id=55, active=is_active(params, 25), &
|
||||
labels=[character(64) :: 'psi_n', 'rho_t', 'n_e', 'T_e', 'q', 'J_φ'])
|
||||
|
||||
if (.not. tbl%active) return
|
||||
if (params%equilibrium%iequil == EQ_VACUUM) return
|
||||
|
||||
! Δρ_p for the uniform ρ_p grid
|
||||
! Note: we don't use ψ_n because J_phi has been
|
||||
! sampled uniformly in ρ_p (see flux_average)
|
||||
drho_p = 1.0_wp_ / (npsi - 1)
|
||||
|
||||
! extra points to reach ψ=ψ_bnd
|
||||
ntail = ceiling((sqrt(params%profiles%psnbnd) - 1) / drho_p)
|
||||
|
||||
do i = 0, npsi + ntail
|
||||
rho_p = i * drho_p
|
||||
rho_t = frhotor(rho_p)
|
||||
psi_n = rho_p**2
|
||||
if (psi_n < 1) then
|
||||
J_phi = vajphiav(i+1) * 1.e-6_wp_
|
||||
else
|
||||
J_phi = 0
|
||||
end if
|
||||
call density(psi_n, dens, ddens)
|
||||
call tbl%append([psi_n, rho_t, dens, temp(psi_n), fq(psi_n), J_phi])
|
||||
end do
|
||||
|
||||
end function kinetic_profiles
|
||||
|
||||
|
||||
function ec_resonance(params, B_res) result(tbl)
|
||||
! Generates the EC resonance surface table
|
||||
|
||||
use const_and_precisions, only : comp_eps
|
||||
use gray_params, only : gray_parameters, EQ_VACUUM
|
||||
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield
|
||||
use magsurf_data, only : npsi
|
||||
use types, only : item
|
||||
use cniteq, only : cniteq_f
|
||||
|
||||
! function arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
real(wp_), intent(in) :: B_res
|
||||
type(table) :: tbl
|
||||
|
||||
call tbl%init(title='ec-resonance', id=70, active=is_active(params, 70), &
|
||||
labels=[character(64) :: 'i', 'B_tot', 'R', 'z'])
|
||||
|
||||
if (.not. tbl%active) return
|
||||
if (params%equilibrium%iequil == EQ_VACUUM) return
|
||||
|
||||
block
|
||||
! local variables
|
||||
integer :: j, k, n
|
||||
integer :: n_conts, n_points(10), n_tot
|
||||
real(wp_) :: dr, dz, B_min, B_max
|
||||
real(wp_) :: B, B_R, B_z, B_phi, B_tot(npsi, npsi)
|
||||
real(wp_), dimension(2002) :: R_cont, z_cont
|
||||
real(wp_), dimension(npsi) :: R, z
|
||||
|
||||
! Build a regular (R, z) grid
|
||||
dr = (rmxm - rmnm - comp_eps)/(npsi - 1)
|
||||
dz = (zmxm - zmnm)/(npsi - 1)
|
||||
do j=1,npsi
|
||||
R(j) = comp_eps + rmnm + dr*(j - 1)
|
||||
z(j) = zmnm + dz*(j - 1)
|
||||
end do
|
||||
|
||||
! B_tot on psi grid
|
||||
B_max = -1.0e30_wp_
|
||||
B_min = +1.0e30_wp_
|
||||
do k = 1, npsi
|
||||
do j = 1, npsi
|
||||
call bfield(R(j), z(k), B_R, B_z, B_phi)
|
||||
B_tot(j,k) = sqrt(B_R**2 + B_z**2 + B_phi**2)
|
||||
if(B_tot(j,k) >= B_max) B_max = B_tot(j,k)
|
||||
if(B_tot(j,k) <= B_min) B_min = B_tot(j,k)
|
||||
end do
|
||||
end do
|
||||
|
||||
! compute Btot=Bres/n with n=1,5
|
||||
do n = 1, 5
|
||||
B = B_res/n
|
||||
if (B >= B_min .and. B <= B_max) then
|
||||
n_conts = size(n_points)
|
||||
n_tot = size(R_cont)
|
||||
call cniteq_f(R, z, B_tot, size(R), size(z), B, &
|
||||
n_conts, n_points, n_tot, R_cont, z_cont)
|
||||
do j = 1, size(R_cont)
|
||||
call tbl%append([wrap(j), wrap(B), wrap(R_cont(j)), wrap(z_cont(j))])
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
end block
|
||||
|
||||
end function ec_resonance
|
||||
|
||||
|
||||
function inputs_maps(params, B_res, X_norm) result(tbl)
|
||||
! Generates 2D maps of several input quantities
|
||||
|
||||
use const_and_precisions, only : comp_eps, cm, degree
|
||||
use gray_params, only : gray_parameters, EQ_VACUUM
|
||||
use equilibrium, only : rmnm, rmxm, zmnm, zmxm, pol_flux, bfield
|
||||
use coreprofiles, only : density, temp
|
||||
use magsurf_data, only : npsi
|
||||
|
||||
! function arguments
|
||||
type(gray_parameters), intent(in) :: params
|
||||
real(wp_), intent(in) :: B_res ! resonant magnetic field, e/m_eω
|
||||
real(wp_), intent(in) :: X_norm ! X normalisation, e²/ε₀m_eω²
|
||||
type(table) :: tbl
|
||||
|
||||
! local variables
|
||||
integer :: j, k
|
||||
real(wp_) :: R0, Npl0
|
||||
real(wp_) :: dR, dz, B_R, B_phi, B_z, B
|
||||
real(wp_) :: psi_n, ne, dne, Te, X, Y, Npl
|
||||
real(wp_), dimension(npsi) :: R, z
|
||||
|
||||
call tbl%init(title='inputs-maps', id=72, active=is_active(params, 72), &
|
||||
labels=[character(64) :: 'R', 'z', 'ψ_n', 'B_R', 'B_z', 'B', &
|
||||
'n_e', 'T_e', 'X', 'Y', 'N_∥'])
|
||||
|
||||
if (.not. tbl%active) return
|
||||
if (params%equilibrium%iequil == EQ_VACUUM) return
|
||||
|
||||
R0 = norm2(params%antenna%pos(1:2)) * cm ! initial value of R
|
||||
Npl0 = sin(params%antenna%beta*degree) ! initial value of N∥
|
||||
|
||||
! Build a regular (R, z) grid
|
||||
dR = (rmxm - rmnm - comp_eps)/(npsi - 1)
|
||||
dz = (zmxm - zmnm)/(npsi - 1)
|
||||
do j = 1, npsi
|
||||
R(j) = comp_eps + rmnm + dR*(j - 1)
|
||||
z(j) = zmnm + dz*(j - 1)
|
||||
end do
|
||||
|
||||
do j = 1, npsi
|
||||
Npl = Npl0 * R0/r(j)
|
||||
do k = 1, npsi
|
||||
call pol_flux(r(j), z(k), psi_n)
|
||||
call bfield(r(j), z(k), B_R, B_z, B_phi)
|
||||
call density(psi_n, ne, dne)
|
||||
B = sqrt(B_R**2 + B_phi**2 + B_z**2)
|
||||
X = X_norm*ne
|
||||
Y = B/B_res
|
||||
Te = temp(psi_n)
|
||||
call tbl%append([R(j), z(k), psi_n, B_R, B_phi, &
|
||||
B_z, B, ne, Te, X, Y, Npl])
|
||||
end do
|
||||
end do
|
||||
|
||||
end function inputs_maps
|
||||
|
||||
|
||||
subroutine find_flux_surfaces(qvals, tbl)
|
||||
! Finds the ψ for a set of values of q and stores the
|
||||
! associated surfaces to the flux surface table
|
||||
|
||||
use equilibrium, only : fq, frhotor, rmaxis, zmaxis, zbsup, zbinf
|
||||
use magsurf_data, only : contours_psi, npoints
|
||||
use logger, only : log_info
|
||||
use minpack, only : hybrj1
|
||||
use types, only : table
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: qvals(:)
|
||||
type(table), intent(inout) :: tbl
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
real(wp_) :: rho_t, rho_p, psi_n
|
||||
character(256) :: msg ! for log messages formatting
|
||||
|
||||
call log_info('constant ψ surfaces for:', &
|
||||
mod='gray_tables', proc='find_flux_surfaces')
|
||||
do i = 1, size(qvals)
|
||||
! Find value of ψ_n for the current q
|
||||
block
|
||||
real(wp_) :: sol(1), fvec(1), fjac(1,1), wa(7)
|
||||
integer :: info
|
||||
|
||||
! Note: since we are mostly interested in q=3/2,2 we start
|
||||
! searching near the boundary in case q is not monotonic.
|
||||
sol = [0.8_wp_] ! first guess
|
||||
|
||||
! Solve fq(ψ_n) = qvals(i) for ψ_n
|
||||
call hybrj1(equation, n=1, x=sol, fvec=fvec, fjac=fjac, &
|
||||
ldfjac=1, tol=1e-3_wp_, info=info, wa=wa, lwa=7)
|
||||
! Solution
|
||||
psi_n = sol(1)
|
||||
|
||||
! Handle spurious or no solutions
|
||||
if (info /= 1 .or. psi_n < 0 .or. psi_n > 1) cycle
|
||||
end block
|
||||
|
||||
! Compute and print the ψ_n(R,z) = ψ_n₀ contour
|
||||
block
|
||||
real(wp_), dimension(npoints) :: R_cont, z_cont
|
||||
real(wp_) :: R_hi, R_lo, z_hi, z_lo
|
||||
|
||||
! Guesses for the high,low horzizontal-tangent points
|
||||
R_hi = rmaxis;
|
||||
z_hi = (zbsup + zmaxis)/2
|
||||
R_lo = rmaxis
|
||||
z_lo = (zbinf + zmaxis)/2
|
||||
|
||||
call contours_psi(psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo)
|
||||
call store_flux_surface(tbl, psi_n, R_cont, z_cont)
|
||||
end block
|
||||
|
||||
! Log the values found for ψ_n, ρ_p, ρ_t
|
||||
rho_p = sqrt(psi_n)
|
||||
rho_t = frhotor(rho_p)
|
||||
write (msg, '(4(x,a,"=",g0.3))') &
|
||||
'q', qvals(i), 'ψ_n', psi_n, 'ρ_p', rho_p, 'ρ_t', rho_t
|
||||
call log_info(msg, mod='gray_tables', proc='find_flux_surfaces')
|
||||
end do
|
||||
|
||||
contains
|
||||
|
||||
subroutine equation(n, x, f, df, ldf, flag)
|
||||
! The equation to solve: f(x) = q(x) - q₀ = 0
|
||||
use const_and_precisions, only : comp_eps
|
||||
|
||||
! subroutine arguments
|
||||
integer, intent(in) :: n, ldf, flag
|
||||
real(wp_), intent(in) :: x(n)
|
||||
real(wp_), intent(inout) :: f(n), df(ldf,n)
|
||||
|
||||
! optimal step size
|
||||
real(wp_), parameter :: e = comp_eps**(1/3.0_wp_)
|
||||
|
||||
if (flag == 1) then
|
||||
! return f(x)
|
||||
f(1) = fq(x(1)) - qvals(i)
|
||||
else
|
||||
! return f'(x), computed numerically
|
||||
df(1,1) = (fq(x(1) + e) - fq(x(1) - e)) / (2*e)
|
||||
end if
|
||||
end subroutine
|
||||
|
||||
end subroutine find_flux_surfaces
|
||||
|
||||
|
||||
subroutine store_flux_surface(tbl, psi_n, R, z)
|
||||
! Helper to add a contour to the flux surfaces table
|
||||
|
||||
use const_and_precisions, only : wp_, comp_tiny
|
||||
use types, only : table, wrap
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: psi_n
|
||||
real(wp_), intent(in) :: R(:)
|
||||
real(wp_), intent(in) :: z(:)
|
||||
type(table), intent(inout) :: tbl
|
||||
|
||||
! local variables
|
||||
integer :: i
|
||||
|
||||
if (.not. tbl%active) return
|
||||
if (psi_n < comp_tiny) return
|
||||
|
||||
do i = 1, size(R)
|
||||
call tbl%append([wrap(i), wrap(psi_n), wrap(R(i)), wrap(z(i))])
|
||||
end do
|
||||
end subroutine store_flux_surface
|
||||
|
||||
|
||||
subroutine store_ray_data(tables, &
|
||||
i, jk, s, P0, pos, psi_n, B, b_n, k0, &
|
||||
N_pl, N_pr, N, N_pr_im, n_e, T_e, &
|
||||
alpha, tau, dIds, nhm, nhf, iokhawa, &
|
||||
index_rt, lambda_r, lambda_i, X, Y, grad_X)
|
||||
! Stores some ray variables and local quantities
|
||||
|
||||
use const_and_precisions, only : degree, cm
|
||||
use equilibrium, only : frhotor
|
||||
use gray_params, only : gray_tables, istpl0
|
||||
use beamdata, only : nray,nrayth,jkray1
|
||||
|
||||
! subroutine arguments
|
||||
|
||||
! tables where to store the data
|
||||
type(gray_tables), intent(inout) :: tables
|
||||
|
||||
integer, intent(in) :: i, jk ! step, ray index
|
||||
real(wp_), intent(in) :: s ! arclength
|
||||
real(wp_), intent(in) :: P0 ! initial power
|
||||
real(wp_), intent(in) :: pos(3) ! position vector
|
||||
real(wp_), intent(in) :: B, b_n(3) ! magnetic field and unit vector
|
||||
real(wp_), intent(in) :: N(3) ! refractive index vector, N
|
||||
real(wp_), intent(in) :: psi_n ! normalised poloidal flux, ψ_n
|
||||
real(wp_), intent(in) :: N_pl, N_pr ! N_∥, N_⊥
|
||||
real(wp_), intent(in) :: N_pr_im ! Im(N_⊥) from warm dispersion
|
||||
real(wp_), intent(in) :: k0 ! vacuum wavenumber (k₀=ω/c)
|
||||
real(wp_), intent(in) :: n_e, T_e ! electron density and temperature
|
||||
real(wp_), intent(in) :: alpha, tau ! absorption coeff., optical depth
|
||||
real(wp_), intent(in) :: dIds ! ECCD current density
|
||||
integer, intent(in) :: nhm, nhf ! EC harmonic numbers
|
||||
integer, intent(in) :: iokhawa !
|
||||
integer, intent(in) :: index_rt !
|
||||
real(wp_), intent(in) :: lambda_r, lambda_i ! quasioptical dispersion relation
|
||||
real(wp_), intent(in) :: X, Y ! X=(ω_pe/ω)², Y=ω_ce/ω
|
||||
real(wp_), intent(in) :: grad_X(3) ! gradient of X
|
||||
|
||||
! local variables
|
||||
real(wp_) :: s_m, pos_m(3), R_m, phi_deg ! coordinates in SI units
|
||||
real(wp_) :: rho_t, k_im, Pt, dIds_n
|
||||
integer :: k
|
||||
|
||||
s_m = s * cm
|
||||
pos_m = pos * cm
|
||||
R_m = hypot(pos_m(1), pos_m(2))
|
||||
|
||||
! Store central ray data
|
||||
if (jk == 1 .and. tables%central_ray%active) then
|
||||
phi_deg = atan2(pos_m(2), pos_m(1)) / degree
|
||||
if(psi_n >= 0 .and. psi_n <= 1) then
|
||||
rho_t = frhotor(sqrt(psi_n))
|
||||
else
|
||||
rho_t = 1.0_wp_
|
||||
end if
|
||||
k_im = N_pr_im * k0 / cm
|
||||
Pt = P0 * exp(-tau)
|
||||
dIds_n = (1/P0 * dIds) / cm
|
||||
|
||||
call tables%central_ray%append([ &
|
||||
wrap(s_m), wrap(R_m), wrap(pos_m(3)), wrap(phi_deg), wrap(psi_n), &
|
||||
wrap(rho_t), wrap(n_e), wrap(T_e), wrap(B), wrap(b_n(1)), &
|
||||
wrap(b_n(2)), wrap(b_n(3)), wrap(N_pr), wrap(N_pl), wrap(N(1)), &
|
||||
wrap(N(2)), wrap(N(3)), wrap(k_im), wrap(alpha), wrap(tau), wrap(Pt), &
|
||||
wrap(dIds_n), wrap(nhm), wrap(nhf), wrap(iokhawa), wrap(index_rt), &
|
||||
wrap(lambda_r), wrap(X), wrap(Y), wrap(grad_X(1)), wrap(grad_X(2)), &
|
||||
wrap(grad_X(3))])
|
||||
end if
|
||||
|
||||
! Store outer rays data
|
||||
if (mod(i, istpl0) == 0) then
|
||||
k = jk - jkray1 + 1
|
||||
if(k > 0 .and. k <= nrayth .and. tables%outer_rays%active) then
|
||||
call tables%outer_rays%append([ &
|
||||
wrap(i), wrap(k), wrap(s_m), wrap(pos_m(1)), wrap(pos_m(2)), &
|
||||
wrap(R_m), wrap(pos_m(3)), wrap(psi_n), wrap(tau), wrap(N_pl), &
|
||||
wrap(alpha), wrap(index_rt)])
|
||||
end if
|
||||
end if
|
||||
|
||||
! Store dispersion relation data
|
||||
if (jk == nray .and. tables%dispersion%active) then
|
||||
call tables%dispersion%append([s, lambda_r, lambda_i])
|
||||
end if
|
||||
|
||||
end subroutine store_ray_data
|
||||
|
||||
|
||||
subroutine store_ec_profiles(tbl, rho_p, rho_t, J_phi, J_cd, &
|
||||
dPdV, I_ins, P_ins, index_rt)
|
||||
! Helper to add data to the ECRH&CD profiles table
|
||||
|
||||
use types, only : table, wrap
|
||||
|
||||
! subroutine arguments
|
||||
type(table), intent(inout) :: tbl
|
||||
real(wp_), dimension(:), intent(in) :: rho_p, rho_t, J_phi, J_cd, &
|
||||
dPdV, I_ins, P_ins
|
||||
integer, intent(in) :: index_rt
|
||||
|
||||
integer :: i
|
||||
|
||||
if (.not. tbl%active) return
|
||||
|
||||
do i = 1, size(rho_p)
|
||||
call tbl%append([wrap(rho_p(i)), wrap(rho_t(i)), wrap(J_phi(i)), &
|
||||
wrap(J_cd(i)), wrap(dPdV(i)), wrap(I_ins(i)), &
|
||||
wrap(P_ins(i)), wrap(index_rt)])
|
||||
end do
|
||||
end subroutine store_ec_profiles
|
||||
|
||||
|
||||
subroutine store_beam_shape(params, tables, s, y, full)
|
||||
! Computes the beam shape tables
|
||||
|
||||
use const_and_precisions, only : zero, one, comp_huge
|
||||
use types, only : table, wrap
|
||||
use gray_params, only : raytracing_parameters, gray_tables
|
||||
use beamdata, only : rayi2jk
|
||||
|
||||
! subroutine arguments
|
||||
type(raytracing_parameters), intent(in) :: params
|
||||
type(gray_tables), target, intent(inout) :: tables
|
||||
real(wp_), intent(in) :: s(:) ! arclength
|
||||
real(wp_), intent(in) :: y(:, :) ! (x̅, N̅) vector
|
||||
|
||||
! whether to show all rays or only the outer ring
|
||||
logical, intent(in), optional :: full
|
||||
|
||||
! local variables
|
||||
logical :: full_
|
||||
integer :: jk, jk0, jkv(2)
|
||||
real(wp_) :: n(3), M(3,3), dx(3), dx_loc(3), c
|
||||
real(wp_) :: r, r_min, r_max
|
||||
type(table), pointer :: beam_shape
|
||||
|
||||
if (.not. (tables%beam_shape%active .or. tables%beam_size%active)) return
|
||||
|
||||
full_ = .false.
|
||||
if (present(full)) full_ = full
|
||||
|
||||
if (full_) then
|
||||
beam_shape => tables%beam_final
|
||||
else
|
||||
beam_shape => tables%beam_shape
|
||||
end if
|
||||
|
||||
! Convert to the local beam basis
|
||||
!
|
||||
! If n̅ is the normalised direction of the central ray,
|
||||
! and the {e̅₁,e̅₂,e̅₃} the standard basis, the unit vectors of
|
||||
! the beam basis are given by:
|
||||
!
|
||||
! x̅ = n̅×e̅₃ / |n̅×e̅₃| = (n₂/c)e̅₁ + (-n₁/c)e̅₂
|
||||
! y̅ = n̅×(n̅×e̅₃) / |n̅×e̅₃| = (n₁n₃/c)e̅₁ + (n₂n₃/c)e̅₂ - ce̅₃
|
||||
! z̅ = n̅ = n₁e̅₁ + n₂e̅₂ + n₃e̅₃
|
||||
!
|
||||
! where c = |n̅×e̅₃| = √(n₁² + n₂²)
|
||||
!
|
||||
! Or in the case where n̅ ∥ e̅₃ by:
|
||||
!
|
||||
! x̅ = e̅₁
|
||||
! y̅ = (n₃e̅₃)×e̅₁ = n₃e̅₂
|
||||
! z̅ = n₃e̅₃
|
||||
!
|
||||
! So, the change of basis matrix is
|
||||
!
|
||||
! [n₂/c -n₁/c 0]
|
||||
! M = [n₁n₃/c n₂n₃/c -c], c > 0
|
||||
! [n₁ n₂ n₃]
|
||||
!
|
||||
! or M = diag(1, n₃, n₃) if c = 0.
|
||||
!
|
||||
! By definition then we have: x̅_loc = M x̅.
|
||||
!
|
||||
n = y(4:6, 1) / norm2(y(4:6, 1))
|
||||
c = norm2(n(1:2))
|
||||
if (c > 0) then
|
||||
M = reshape([ n(2)/c, n(1)*n(3)/c, n(1), &
|
||||
-n(1)/c, n(2)*n(3)/c, n(2), &
|
||||
zero, -c, n(3) ], [3, 3])
|
||||
else
|
||||
M = reshape([ one, zero, zero, &
|
||||
zero, n(3), zero, &
|
||||
zero, zero, n(3) ], [3, 3])
|
||||
end if
|
||||
|
||||
! start from the central ray or the last ring
|
||||
jk0 = 1 + merge(0, params%nray - params%nrayth, full_)
|
||||
|
||||
r_min = comp_huge
|
||||
r_max = 0
|
||||
do jk = jk0, params%nray
|
||||
dx = y(1:3, jk) - y(1:3, 1) ! distance from the central ray
|
||||
dx_loc = matmul(M, dx) ! local ray position
|
||||
r = norm2(dx_loc(1:2)) ! cross-section radius
|
||||
jkv = rayi2jk(jk) ! [j, k] vector
|
||||
|
||||
if (full_ .or. jk /= 1) &
|
||||
call beam_shape%append([ &
|
||||
wrap(s(jk)), wrap(jkv(1)), wrap(jkv(2)), &
|
||||
wrap(dx_loc(1)), wrap(dx_loc(2)), wrap(dx_loc(3)), wrap(r)])
|
||||
|
||||
! Compute min/max radius
|
||||
if (r>=r_max .and. jkv(1)==params%nrayr) r_max = r
|
||||
if (r<=r_min .and. jkv(1)==params%nrayr) r_min = r
|
||||
end do
|
||||
|
||||
call tables%beam_size%append([s(1), r_min, r_max])
|
||||
|
||||
end subroutine store_beam_shape
|
||||
|
||||
end module gray_tables
|
||||
@ -95,11 +95,16 @@ contains
|
||||
end subroutine dealloc_surfvec
|
||||
|
||||
|
||||
subroutine flux_average
|
||||
subroutine compute_flux_averages(tables)
|
||||
use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv
|
||||
use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, &
|
||||
bfield,frhotor,fq,tor_curr,psia,pol_flux
|
||||
use dierckx, only : regrid,coeff_parder
|
||||
use types, only : table, wrap
|
||||
use gray_params, only : gray_tables
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_tables), intent(inout), optional :: tables
|
||||
|
||||
! local constants
|
||||
integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3, &
|
||||
@ -384,21 +389,29 @@ contains
|
||||
njpt=njp
|
||||
nlmt=nlm
|
||||
|
||||
deallocate(rctemp, zctemp)
|
||||
|
||||
if (.not. present(tables)) return
|
||||
|
||||
do jp=1,npsi
|
||||
call print_fluxav(rpstab(jp),frhotor(rpstab(jp)),bav(jp),bmxpsi(jp), &
|
||||
if (.not. tables%flux_averages%active) exit
|
||||
call tables%flux_averages%append([ &
|
||||
rpstab(jp), frhotor(rpstab(jp)), bav(jp), bmxpsi(jp), &
|
||||
bmnpsi(jp), varea(jp), vvol(jp), vcurrp(jp), vajphiav(jp), &
|
||||
ffc(jp),vratja(jp),vratjb(jp),qqv(jp))
|
||||
ffc(jp), vratja(jp), vratjb(jp), qqv(jp)])
|
||||
end do
|
||||
|
||||
ninpr = (npsi-1)/10
|
||||
do jp = ninpr+1, npsi, ninpr
|
||||
call print_contour(psicon(jp),rcon(:,jp),zcon(:,jp))
|
||||
! Note: can't use store_flux_surface to avoid a circular dependency
|
||||
if (.not. tables%flux_surfaces%active) exit
|
||||
do l = 1, npoints
|
||||
call tables%flux_surfaces%append( &
|
||||
[wrap(l), wrap(psicon(jp)), wrap(rcon(l,jp)), wrap(zcon(l,jp))])
|
||||
end do
|
||||
end do
|
||||
|
||||
deallocate(rctemp, zctemp)
|
||||
end subroutine flux_average
|
||||
|
||||
end subroutine compute_flux_averages
|
||||
|
||||
|
||||
subroutine fluxval(rhop,area,vol,dervol,dadrhot,dvdrhot, &
|
||||
@ -514,54 +527,4 @@ contains
|
||||
end if
|
||||
end subroutine contours_psi
|
||||
|
||||
|
||||
subroutine print_contour(psin, rc, zc)
|
||||
! Prints the flux surface contours (unit 71)
|
||||
|
||||
use const_and_precisions, only : wp_, comp_tiny
|
||||
use units, only : ucnt, unit_active
|
||||
|
||||
! subroutine arguments
|
||||
real(wp_), intent(in) :: psin
|
||||
real(wp_), dimension(:), intent(in) :: rc, zc
|
||||
|
||||
! local variables
|
||||
integer :: npoints, ic
|
||||
|
||||
if (.not. unit_active(ucnt)) return
|
||||
|
||||
if (psin < comp_tiny) then
|
||||
write(ucnt, *)' #i psin R z'
|
||||
else
|
||||
npoints = size(rc)
|
||||
do ic=1,npoints
|
||||
write(ucnt, '(i6,12(1x,e12.5))') ic, psin, rc(ic), zc(ic)
|
||||
end do
|
||||
write(ucnt, *)
|
||||
write(ucnt, *)
|
||||
end if
|
||||
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)
|
||||
|
||||
use const_and_precisions, only : wp_, comp_tiny
|
||||
use units, only : uflx, unit_active
|
||||
|
||||
! subroutine arguments
|
||||
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) &
|
||||
write(uflx, *)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol' // &
|
||||
' |I_pl| <J_phi> fc 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 module magsurf_data
|
||||
|
||||
268
src/main.f90
268
src/main.f90
@ -1,7 +1,6 @@
|
||||
program main
|
||||
use const_and_precisions, only : wp_, zero
|
||||
use const_and_precisions, only : wp_
|
||||
use logger, only : INFO, ERROR, set_log_level, log_message
|
||||
use units, only : set_active_units, close_units
|
||||
use utils, only : dirname
|
||||
use gray_cli, only : cli_options, parse_cli_options, &
|
||||
deinit_cli_options, parse_param_overrides
|
||||
@ -31,9 +30,6 @@ program main
|
||||
if (opts%quiet) opts%verbose = ERROR
|
||||
call set_log_level(opts%verbose)
|
||||
|
||||
! Activate the given output units
|
||||
call set_active_units(opts%units)
|
||||
|
||||
! Load the parameters from file and move to its directory
|
||||
! (all other filepaths are assumed relative to it)
|
||||
if (allocated(opts%config_file)) then
|
||||
@ -79,6 +75,9 @@ program main
|
||||
|
||||
call init_misc(params, data)
|
||||
|
||||
! Activate the given output tables
|
||||
params%misc%active_tables = opts%tables
|
||||
|
||||
! Get back to the initial directory
|
||||
err = chdir(cwd)
|
||||
|
||||
@ -93,93 +92,7 @@ program main
|
||||
|
||||
if (allocated(opts%sum_filelist)) then
|
||||
call log_message(level=INFO, mod='main', msg='summing profiles')
|
||||
|
||||
sum: block
|
||||
real(wp_), dimension(:), allocatable :: jphi
|
||||
real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin
|
||||
real(wp_), dimension(:), allocatable :: extracol
|
||||
integer :: i, j, k, l, n, ngam, nextracol, irt
|
||||
character(len=255) :: filename, headerline, fmtstr
|
||||
real(wp_), dimension(5) :: f48v
|
||||
real(wp_) :: 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=opts%sum_filelist, action='read', status='old', iostat=err)
|
||||
if (err /= 0) then
|
||||
call log_message(level=ERROR, mod='main', &
|
||||
msg='opening file list ('//opts%sum_filelist//') failed!')
|
||||
call exit(1)
|
||||
end if
|
||||
|
||||
open(unit=100 -1, file='f48sum.txt', action='write', status='unknown')
|
||||
open(unit=100 -2, file='f7sum.txt', action='write', status='unknown')
|
||||
|
||||
read(100, *) n, ngam, nextracol
|
||||
allocate(extracol(nextracol))
|
||||
do i=1,n
|
||||
read(100, *) filename
|
||||
open(100 + i, file=filename, action='read', status='old', iostat=err)
|
||||
if (err /= 0) then
|
||||
call log_message(level=ERROR, mod='main', &
|
||||
msg='opening summand file ('//trim(filename)//') failed!')
|
||||
call exit(1)
|
||||
end if
|
||||
|
||||
do j=1,22
|
||||
read(100 + i, '(a255)') headerline
|
||||
if (i == 1) then
|
||||
write(100 -1, '(a)') trim(headerline)
|
||||
write(100 -2, '(a)') trim(headerline)
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
allocate(results%jcd(params%output%nrho), results%dpdv(params%output%nrho))
|
||||
do k=1,ngam
|
||||
jphi = zero
|
||||
results%jcd = zero
|
||||
results%dpdv = zero
|
||||
currins = zero
|
||||
pins = zero
|
||||
do j=1,params%output%nrho
|
||||
do i=1,n
|
||||
read(100+i, *) (extracol(l), l=1,nextracol), &
|
||||
rpin(j), rtin(j), f48v(1:5), irt
|
||||
jphi(j) = f48v(1) + jphi(j)
|
||||
results%jcd(j) = f48v(2) + results%jcd(j)
|
||||
results%dpdv(j) = f48v(3) + results%dpdv(j)
|
||||
currins(j) = f48v(4) + currins(j)
|
||||
pins(j) = f48v(5) + pins(j)
|
||||
end do
|
||||
write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracol+7
|
||||
write(100 -1,trim(fmtstr)) &
|
||||
(extracol(l), l=1,nextracol), rpin(j), rtin(j), &
|
||||
jphi(j), results%jcd(j), results%dpdv(j), &
|
||||
currins(j), pins(j), irt
|
||||
end do
|
||||
results%pabs = pins(params%output%nrho)
|
||||
results%icd = currins(params%output%nrho)
|
||||
write(100 -1, *)
|
||||
call sum_profiles(params, jphi, results%jcd, results%dpdv, &
|
||||
currins, pins, results%pabs, results%icd, &
|
||||
jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
|
||||
drhotjava, drhotpav, ratjamx, ratjbmx)
|
||||
write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracol+12
|
||||
write(100 -2, trim(fmtstr)) &
|
||||
(extracol(l), l=1,nextracol), results%icd, results%pabs, &
|
||||
jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
|
||||
drhotjava, drhotpav, ratjamx, ratjbmx
|
||||
end do
|
||||
do i=-2,n
|
||||
close(100 + i)
|
||||
end do
|
||||
deallocate(jphi, currins, pins, rtin, rpin)
|
||||
deallocate(extracol)
|
||||
deallocate(opts%params_file)
|
||||
end block sum
|
||||
call sum_profiles(params, results)
|
||||
else
|
||||
call gray_main(params, data, results, err)
|
||||
end if
|
||||
@ -192,16 +105,40 @@ program main
|
||||
call log_message(msg, level=INFO, mod='main')
|
||||
end block print_res
|
||||
|
||||
write_res: block
|
||||
integer :: i
|
||||
|
||||
! Write results to file
|
||||
do i = 1, size(results%tables%iter)
|
||||
if (.not. associated(results%tables%iter(i)%ptr)) cycle
|
||||
|
||||
associate (tbl => results%tables%iter(i)%ptr)
|
||||
call tbl%save(err, make_header=make_header)
|
||||
if (err /= 0) then
|
||||
call log_message('failed to save '//trim(tbl%title)//' table', &
|
||||
level=ERROR, mod='main')
|
||||
end if
|
||||
end associate
|
||||
end do
|
||||
end block write_res
|
||||
|
||||
! Free memory
|
||||
call deinit_equilibrium(data%equilibrium)
|
||||
call deinit_profiles(data%profiles)
|
||||
call deinit_misc
|
||||
call deinit_cli_options(opts)
|
||||
deallocate(results%dpdv, results%jcd)
|
||||
call close_units
|
||||
|
||||
contains
|
||||
|
||||
subroutine make_header(file)
|
||||
! Writes the gray output files header
|
||||
use gray_params, only : write_gray_header
|
||||
integer, intent(in) :: file
|
||||
call write_gray_header(params, file)
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine init_equilibrium(params, data, err)
|
||||
! Reads the MHD equilibrium file (either in the G-EQDSK format
|
||||
! or an analytical description) and initialises the respective
|
||||
@ -462,83 +399,120 @@ contains
|
||||
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_profiles_an, set_profiles_spline, &
|
||||
temp, fzeff
|
||||
use dispersion, only : expinit
|
||||
use gray_params, only : gray_parameters, print_parameters
|
||||
use beams, only : launchangles2n, xgygcoeff
|
||||
use magsurf_data, only : flux_average, dealloc_surfvec
|
||||
subroutine sum_profiles(params, results)
|
||||
use gray_params, only : gray_parameters
|
||||
use magsurf_data, only : compute_flux_averages, dealloc_surfvec
|
||||
use beamdata, only : init_btr
|
||||
use pec, only : pec_init, postproc_profiles, dealloc_pec, &
|
||||
rhop_tab, rhot_tab
|
||||
use gray_core, only : print_headers, print_finals, print_pec, &
|
||||
print_bres, print_prof, print_maps, &
|
||||
print_surfq
|
||||
rhot_tab
|
||||
|
||||
! subroutine arguments
|
||||
type(gray_parameters), intent(inout) :: 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
|
||||
type(gray_results), intent(inout) :: results
|
||||
|
||||
real(wp_), dimension(:), allocatable :: jphi
|
||||
real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin
|
||||
real(wp_), dimension(:), allocatable :: extracol
|
||||
integer :: i, j, k, l, n, ngam, nextracol, irt
|
||||
character(len=255) :: filename, headerline, fmtstr
|
||||
real(wp_), dimension(5) :: f48v
|
||||
real(wp_) :: 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
|
||||
|
||||
! ======== set environment BEGIN ========
|
||||
! Compute X=(ω_pe/ω)² and Y=ω_ce/ω (with B=1)
|
||||
call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
|
||||
|
||||
! Initialise the ray variables (beamtracing)
|
||||
call init_btr(params%raytracing)
|
||||
|
||||
! 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
|
||||
call compute_flux_averages
|
||||
|
||||
! Initialise the output profiles
|
||||
call pec_init(params%output%ipec)
|
||||
! ======= set environment END ======
|
||||
|
||||
! ======== pre-proc prints BEGIN ========
|
||||
call print_headers(params)
|
||||
associate(nrho =>params%output%nrho)
|
||||
allocate(jphi(nrho), currins(nrho), pins(nrho), rtin(nrho), rpin(nrho))
|
||||
allocate(results%jcd(nrho), results%dpdv(nrho))
|
||||
end associate
|
||||
|
||||
! 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_])
|
||||
open(100, file=opts%sum_filelist, action='read', status='old', iostat=err)
|
||||
if (err /= 0) then
|
||||
call log_message(level=ERROR, mod='main', &
|
||||
msg='opening file list ('//opts%sum_filelist//') failed!')
|
||||
call exit(1)
|
||||
end if
|
||||
|
||||
! Print ne, Te, q, Jphi versus psi, rhop, rhot
|
||||
call print_bres(bres)
|
||||
call print_prof(params%profiles)
|
||||
call print_maps(bres, xgcn, &
|
||||
norm2(params%antenna%pos(1:2)) *0.01_wp_ , &
|
||||
sin(params%antenna%beta*degree))
|
||||
! ========= pre-proc prints END =========
|
||||
open(unit=100 -1, file='f48sum.txt', action='write', status='unknown')
|
||||
open(unit=100 -2, file='f7sum.txt', action='write', status='unknown')
|
||||
|
||||
! 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, &
|
||||
read(100, *) n, ngam, nextracol
|
||||
allocate(extracol(nextracol))
|
||||
do i=1,n
|
||||
read(100, *) filename
|
||||
open(100 + i, file=filename, action='read', status='old', iostat=err)
|
||||
if (err /= 0) then
|
||||
call log_message(level=ERROR, mod='main', &
|
||||
msg='opening summand file ('//trim(filename)//') failed!')
|
||||
call exit(1)
|
||||
end if
|
||||
|
||||
do j=1,22
|
||||
read(100 + i, '(a255)') headerline
|
||||
if (i == 1) then
|
||||
write(100 -1, '(a)') trim(headerline)
|
||||
write(100 -2, '(a)') trim(headerline)
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
do k=1,ngam
|
||||
jphi = 0
|
||||
results%jcd = 0
|
||||
results%dpdv = 0
|
||||
currins = 0
|
||||
pins = 0
|
||||
do j=1,params%output%nrho
|
||||
do i=1,n
|
||||
read(100+i, *) (extracol(l), l=1,nextracol), &
|
||||
rpin(j), rtin(j), f48v(1:5), irt
|
||||
jphi(j) = f48v(1) + jphi(j)
|
||||
results%jcd(j) = f48v(2) + results%jcd(j)
|
||||
results%dpdv(j) = f48v(3) + results%dpdv(j)
|
||||
currins(j) = f48v(4) + currins(j)
|
||||
pins(j) = f48v(5) + pins(j)
|
||||
end do
|
||||
write(fmtstr, '("(",i0,"(1x,e16.8e3),i5)")') nextracol+7
|
||||
write(100 -1,trim(fmtstr)) &
|
||||
(extracol(l), l=1,nextracol), rpin(j), rtin(j), &
|
||||
jphi(j), results%jcd(j), results%dpdv(j), &
|
||||
currins(j), pins(j), irt
|
||||
end do
|
||||
results%pabs = pins(params%output%nrho)
|
||||
results%icd = currins(params%output%nrho)
|
||||
write(100 -1, *)
|
||||
! Recompute final results
|
||||
call postproc_profiles(results%pabs, results%icd, rhot_tab, &
|
||||
results%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)
|
||||
write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracol+12
|
||||
write(100 -2, trim(fmtstr)) &
|
||||
(extracol(l), l=1,nextracol), results%icd, results%pabs, &
|
||||
jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
|
||||
drhotjava, drhotpav, ratjamx, ratjbmx
|
||||
end do
|
||||
|
||||
do i=-2,n
|
||||
close(100 + i)
|
||||
end do
|
||||
deallocate(jphi, currins, pins, rtin, rpin)
|
||||
deallocate(extracol)
|
||||
deallocate(opts%params_file)
|
||||
|
||||
! Free memory
|
||||
call dealloc_surfvec ! for fluxval
|
||||
call dealloc_surfvec
|
||||
call dealloc_pec
|
||||
end subroutine sum_profiles
|
||||
|
||||
|
||||
354
src/types.f90
Normal file
354
src/types.f90
Normal file
@ -0,0 +1,354 @@
|
||||
module types
|
||||
|
||||
use const_and_precisions, only : wp_
|
||||
|
||||
implicit none
|
||||
|
||||
type item
|
||||
! An item is the building block of a queue
|
||||
class(*), allocatable :: value ! the actual content
|
||||
type(item), pointer :: next => null() ! pointer to the next item
|
||||
end type
|
||||
|
||||
type queue
|
||||
! A queue is a list of items with O(1) insertion and extraction.
|
||||
! The first item inserted (`put` operation) is the first to be
|
||||
! extracted (`get` operation).
|
||||
|
||||
! References to the first and last items in the queue
|
||||
type(item), pointer :: first => null(), last => null()
|
||||
|
||||
contains
|
||||
procedure :: put => queue_put ! inserts an item at the end of the queue
|
||||
procedure :: get => queue_get ! extracts the first item in the queue
|
||||
procedure :: empty => queue_empty ! checks whether the queue is empty
|
||||
end type
|
||||
|
||||
type table
|
||||
! A type for storing tabular data before serialisation
|
||||
integer :: id ! the table unique ID
|
||||
logical :: active ! whether to process the table
|
||||
character(64) :: title ! title of the table
|
||||
character(64), allocatable :: labels(:) ! labels for each column
|
||||
type(queue), allocatable :: columns(:) ! values for each column
|
||||
|
||||
contains
|
||||
procedure :: init => table_init ! creates an empty table
|
||||
procedure :: append => table_append ! appends an extra row to the table
|
||||
procedure :: empty => table_empty ! checks whether the queue is empty
|
||||
procedure :: save => table_save ! processes and writes the table to file
|
||||
end type
|
||||
|
||||
type wrap
|
||||
! A wrapper type for storing (references of) heterogeneous
|
||||
! values in an array or other homogeneuous container type
|
||||
class(*), pointer :: value
|
||||
end type
|
||||
|
||||
! Interface for custom type constructor
|
||||
interface wrap
|
||||
procedure :: wrap_init
|
||||
end interface
|
||||
|
||||
! Interface for subroutines that write to file
|
||||
interface
|
||||
subroutine writer(file)
|
||||
integer, intent(in) :: file
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine queue_put(self, val)
|
||||
! Inserts an item of value `val` at the end of the `self` queue
|
||||
class(queue), intent(inout) :: self
|
||||
class(*), intent(in) :: val
|
||||
|
||||
type(item), pointer :: i
|
||||
|
||||
! set up the item and copy the value
|
||||
allocate(i)
|
||||
allocate(i%value, source=val)
|
||||
|
||||
if (self%empty()) then
|
||||
! first item received, store it separately
|
||||
self%first => i
|
||||
else
|
||||
! otherwise, link the last item to the new one
|
||||
self%last%next => i
|
||||
end if
|
||||
|
||||
! update the last item
|
||||
self%last => i
|
||||
end subroutine queue_put
|
||||
|
||||
|
||||
subroutine queue_get(self, val)
|
||||
! Extracts the item from the front of the queue.
|
||||
! Notes:
|
||||
! - this assumes the queue is not empty.
|
||||
! - this should be converted to a function once
|
||||
! gcc bug 115072 is fixed.
|
||||
class(queue), intent(inout) :: self
|
||||
class(*), allocatable, intent(out) :: val
|
||||
|
||||
type(item), pointer :: tmp
|
||||
|
||||
! extract the value
|
||||
allocate(val, source=self%first%value)
|
||||
|
||||
tmp => self%first ! get a grip
|
||||
self%first => tmp%next ! shift the reference
|
||||
deallocate(tmp) ! clear the item
|
||||
end subroutine queue_get
|
||||
|
||||
|
||||
pure function queue_empty(self)
|
||||
! Whether the `self` queue is empty
|
||||
class(queue), intent(in) :: self
|
||||
logical :: queue_empty
|
||||
|
||||
queue_empty = .not. associated(self%first)
|
||||
end function
|
||||
|
||||
|
||||
function wrap_init(val)
|
||||
! Initialise a `wrap` value
|
||||
!
|
||||
! Note: this is needed because, apparently, gfortran can't handle
|
||||
! unlimited polymorphic arguments in the default constructor.
|
||||
class(*), target, intent(in) :: val
|
||||
type(wrap) :: wrap_init
|
||||
wrap_init%value => val
|
||||
end function wrap_init
|
||||
|
||||
|
||||
subroutine table_init(self, id, title, labels, active)
|
||||
! Initialises the table
|
||||
class(table), intent(out) :: self
|
||||
integer, intent(in) :: id
|
||||
character(*), intent(in) :: title
|
||||
character(*), intent(in) :: labels(:)
|
||||
|
||||
logical, intent(in), optional :: active
|
||||
|
||||
integer :: i
|
||||
|
||||
self%id = id
|
||||
self%title = title
|
||||
self%labels = labels
|
||||
self%columns = [(queue(), i=1,size(labels))]
|
||||
|
||||
if (present(active)) self%active = active
|
||||
end subroutine table_init
|
||||
|
||||
|
||||
pure function table_empty(self)
|
||||
! Checks whether the table is empty
|
||||
class(table), intent(in) :: self
|
||||
logical :: table_empty
|
||||
|
||||
table_empty = self%columns(1)%empty()
|
||||
end function table_empty
|
||||
|
||||
|
||||
subroutine table_append(self, row)
|
||||
! Appends an extra row to the table
|
||||
!
|
||||
! Heterogenous values are supported, provided they
|
||||
! are wrapped in a `wrap` type.
|
||||
|
||||
class(table), intent(inout) :: self
|
||||
class(*), intent(in) :: row(:)
|
||||
|
||||
integer :: i
|
||||
|
||||
do i = 1, size(self%columns)
|
||||
select type (row)
|
||||
type is (wrap)
|
||||
call self%columns(i)%put(row(i)%value)
|
||||
class default
|
||||
call self%columns(i)%put(row(i))
|
||||
end select
|
||||
end do
|
||||
|
||||
end subroutine table_append
|
||||
|
||||
|
||||
subroutine table_save(self, error, filepath, make_header)
|
||||
! Save the table to a file
|
||||
!
|
||||
! Note: this operation consumes the table
|
||||
use utils, only : get_free_unit, digits
|
||||
|
||||
! suborutine arguments
|
||||
class(table), intent(inout) :: self
|
||||
integer, intent(out) :: error
|
||||
character(*), intent(in), optional :: filepath
|
||||
procedure(writer), optional :: make_header
|
||||
|
||||
! local variables
|
||||
integer :: i, file
|
||||
class(*), allocatable :: val
|
||||
character(len=len(self%title) + 4) :: filepath_
|
||||
character(len=10) :: fmt
|
||||
|
||||
if (self%empty()) return
|
||||
|
||||
file = get_free_unit()
|
||||
|
||||
! set default file name
|
||||
if (present(filepath)) then
|
||||
filepath_ = filepath
|
||||
else
|
||||
write(filepath_, '(a,".",i0,".txt")') trim(self%title), self%id
|
||||
end if
|
||||
|
||||
! open output file
|
||||
open(unit=file, file=filepath_, action='write', iostat=error)
|
||||
if (error /= 0) return
|
||||
|
||||
! write a custom file header
|
||||
if (present(make_header)) call make_header(file)
|
||||
|
||||
! write table header
|
||||
write(file, '(a)', advance='no') '#'
|
||||
do i = 1, size(self%columns)
|
||||
write(file, '(x,a)', advance='no') trim(self%labels(i))
|
||||
end do
|
||||
write(file, '(x)')
|
||||
|
||||
! iterate the table
|
||||
do while (.not. self%empty())
|
||||
|
||||
! write a row
|
||||
do i = 1, size(self%columns)
|
||||
call self%columns(i)%get(val)
|
||||
|
||||
! format a record
|
||||
select type (val)
|
||||
type is (integer)
|
||||
write(fmt, '("(i",i0,")")') max(5, 1 + digits(val))
|
||||
write(file, fmt, advance='no') val
|
||||
type is (real(wp_))
|
||||
write(file, '(es16.8e3)', advance='no') val
|
||||
type is (logical)
|
||||
write(file, '(l1)', advance='no') val
|
||||
type is (character(len=*))
|
||||
write(file, '(a)', advance='no') trim(val)
|
||||
end select
|
||||
|
||||
! add space between values
|
||||
if (i < size(self%columns)) write(file, '(x)', advance='no')
|
||||
end do
|
||||
|
||||
! start newline
|
||||
write(file, '(x)')
|
||||
end do
|
||||
|
||||
end subroutine table_save
|
||||
|
||||
|
||||
subroutine test_queue()
|
||||
integer, target :: x=1, y=2
|
||||
real, target :: z = 1.231
|
||||
logical, target :: p=.false.
|
||||
character(len=4) :: s='ciao'
|
||||
type(queue) :: q
|
||||
class(*), allocatable :: c
|
||||
|
||||
call assert(q%empty(), 'queue is empty')
|
||||
|
||||
! put some items in
|
||||
call q%put(x)
|
||||
call q%put(y)
|
||||
call q%put(z)
|
||||
call q%put(p)
|
||||
call q%put(s)
|
||||
|
||||
call assert(.not. q%empty(), 'queue is not empty')
|
||||
|
||||
! get them out and check them
|
||||
call q%get(c)
|
||||
call assert(transfer(c, x) == x, 'got back an integer')
|
||||
call q%get(c)
|
||||
call assert(transfer(c, y) == y, 'got back another integer')
|
||||
call q%get(c)
|
||||
call assert(transfer(c, z) == z, 'got back a real')
|
||||
call q%get(c)
|
||||
call assert(transfer(c, p) .eqv. p, 'got back a logical')
|
||||
call q%get(c)
|
||||
call assert(cast_char(c, 4) == s, 'got back a string')
|
||||
|
||||
call assert(q%empty(), 'queue is empty again')
|
||||
|
||||
! this should be finalised
|
||||
call q%put(s)
|
||||
|
||||
contains
|
||||
|
||||
subroutine assert(p, msg)
|
||||
logical, intent(in) :: p
|
||||
character(len=*), intent(in) :: msg
|
||||
|
||||
if (.not. p) then
|
||||
print '(a, a)', 'assertion failed: ', msg
|
||||
stop
|
||||
end if
|
||||
end subroutine
|
||||
|
||||
pure function cast_char(x, n) result(y)
|
||||
class(*), intent(in) :: x
|
||||
integer, intent(in) :: n
|
||||
character(len=n) :: y
|
||||
|
||||
select type (x)
|
||||
type is (character(len=*))
|
||||
y = x
|
||||
class default
|
||||
error stop 'invalid cast to character'
|
||||
end select
|
||||
end function
|
||||
|
||||
end subroutine test_queue
|
||||
|
||||
|
||||
subroutine test_table()
|
||||
type(table) :: tbl
|
||||
integer :: i
|
||||
class(*), allocatable :: x, y
|
||||
|
||||
call tbl%init(title='testing table', id=0, labels=['x', 'y'])
|
||||
|
||||
! fill in some values
|
||||
call assert(tbl%empty(), 'table is empty')
|
||||
do i = 1, 4
|
||||
call tbl%append([0.2*i, (0.2*i)**2])
|
||||
end do
|
||||
call assert(.not. tbl%empty(), 'table is not empty')
|
||||
|
||||
! get them back
|
||||
do i = 1, 4
|
||||
call tbl%columns(1)%get(x)
|
||||
call tbl%columns(2)%get(y)
|
||||
call assert(transfer(x, 1.0) == (0.2*i), 'got correct x')
|
||||
call assert(transfer(y, 1.0) == (0.2*i)**2, 'got correct y')
|
||||
end do
|
||||
|
||||
call assert(tbl%empty(), 'table is empty again')
|
||||
|
||||
contains
|
||||
|
||||
subroutine assert(p, msg)
|
||||
logical, intent(in) :: p
|
||||
character(len=*), intent(in) :: msg
|
||||
|
||||
if (.not. p) then
|
||||
print '(a, a)', 'assertion failed: ', msg
|
||||
stop
|
||||
end if
|
||||
end subroutine
|
||||
|
||||
end subroutine test_table
|
||||
|
||||
end module types
|
||||
@ -1,70 +0,0 @@
|
||||
! This modules contains the definitions of the unit numbers used for
|
||||
! all the GRAY output ifles and a mechanism to toggle them
|
||||
module units
|
||||
|
||||
implicit none
|
||||
|
||||
! 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 :: uprfin = 55, uflx = 56, ubres = 70, ucnt = 71
|
||||
integer, parameter :: uprj0 = 8, uwbm = 12, ucenr = 4, uoutr = 33
|
||||
integer, parameter :: udisp = 17, upec = 48, usumm = 7, umaps = 72
|
||||
integer, parameter :: uhead = 78
|
||||
#endif
|
||||
|
||||
! List of active units
|
||||
logical, save :: all_enabled = .false.
|
||||
integer, allocatable, save :: active_units(:)
|
||||
|
||||
contains
|
||||
|
||||
subroutine set_active_units(units)
|
||||
! Activate the given units
|
||||
!
|
||||
! All units are inactive by default and no output
|
||||
! will be directed to them.
|
||||
|
||||
! subroutine arguments
|
||||
integer, intent(in) :: units(:)
|
||||
|
||||
active_units = units
|
||||
|
||||
end subroutine set_active_units
|
||||
|
||||
|
||||
subroutine close_units
|
||||
! Close all the active units to flush the buffer.
|
||||
|
||||
integer :: i
|
||||
|
||||
if (allocated(active_units)) then
|
||||
do i = 1, size(active_units)
|
||||
close(active_units(i))
|
||||
end do
|
||||
deallocate(active_units)
|
||||
end if
|
||||
|
||||
end subroutine close_units
|
||||
|
||||
|
||||
function unit_active(unit) result(on)
|
||||
! Checks whether the given `unit` is active
|
||||
|
||||
! function arguments
|
||||
integer, intent(in) :: unit
|
||||
logical :: on
|
||||
|
||||
on = all_enabled .or. any(active_units == unit)
|
||||
|
||||
end function unit_active
|
||||
|
||||
end module units
|
||||
@ -378,4 +378,20 @@ contains
|
||||
isrelative = (filepath(1:1) /= '/')
|
||||
end function isrelative
|
||||
|
||||
|
||||
pure function digits(n)
|
||||
! Returns the number of digits of an integer
|
||||
|
||||
! function arguments
|
||||
integer, intent(in) :: n
|
||||
integer :: digits
|
||||
|
||||
! How it works:
|
||||
! 1. `bit_size` returns the number of bits + 1 (sign)
|
||||
! 2. subtracting the number of leading zeros gives the significant bits
|
||||
! 3. multiplying by log₁₀(2) gives the approximate number of digits
|
||||
! 4. `ceiling` rounds to the next closest integer
|
||||
digits = ceiling((bit_size(abs(n)) - leadz(abs(n))) * log10(2.0_wp_))
|
||||
end function digits
|
||||
|
||||
end module utils
|
||||
|
||||
218
src/vendor/cniteq.f90
vendored
Normal file
218
src/vendor/cniteq.f90
vendored
Normal file
@ -0,0 +1,218 @@
|
||||
module cniteq
|
||||
|
||||
use const_and_precisions, only : wp_
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
subroutine cniteq_f(rqgrid,zqgrid,matr2dgrid,nr,nz,h,ncon,npts,icount,rcon,zcon)
|
||||
! v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking.
|
||||
! (based on an older code)
|
||||
! arguments
|
||||
integer, intent(in) :: nr,nz
|
||||
real(wp_), dimension(nr), intent(in) :: rqgrid
|
||||
real(wp_), dimension(nz), intent(in) :: zqgrid
|
||||
real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid
|
||||
real(wp_), intent(in) :: h
|
||||
integer, intent(inout) :: ncon, icount
|
||||
integer, dimension(ncon), intent(out) :: npts
|
||||
real(wp_), dimension(icount), intent(out) :: rcon,zcon
|
||||
! local variables
|
||||
integer :: i,j,k,l,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb
|
||||
integer :: jabs,jnb,kx,ikx,itm,inext,in
|
||||
integer, dimension(3,2) :: ja
|
||||
integer, dimension(icount/2-1) :: lx
|
||||
real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y
|
||||
real(wp_), dimension(nr*nz) :: a
|
||||
logical :: flag1
|
||||
|
||||
px = 0.5_wp_
|
||||
|
||||
a = reshape(matr2dgrid, [nr*nz])
|
||||
|
||||
rcon = 0.0_wp_
|
||||
zcon = 0.0_wp_
|
||||
|
||||
nrqmax = nr
|
||||
drgrd = rqgrid(2) - rqgrid(1)
|
||||
dzgrd = zqgrid(2) - zqgrid(1)
|
||||
|
||||
ncon = 0
|
||||
|
||||
npts = 0
|
||||
|
||||
iclast = 0
|
||||
icount = 0
|
||||
mpl = 0
|
||||
ix = 0
|
||||
mxr = nrqmax * (nz - 1)
|
||||
n1 = nr - 1
|
||||
|
||||
do jx=2,n1
|
||||
do jm=jx,mxr,nrqmax
|
||||
j = jm + nrqmax
|
||||
ah=a(j)-h
|
||||
if (ah <= 0.0_wp_ .and. a(jm) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(jm) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=-j
|
||||
end if
|
||||
if (ah <= 0.0_wp_ .and. a(j-1) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(j-1) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=j
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
do jm=nr,mxr,nrqmax
|
||||
j = jm + nrqmax
|
||||
ah=a(j)-h
|
||||
if (ah <= 0.0_wp_ .and. a(j-1) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(j-1) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=j
|
||||
end if
|
||||
if (ah <= 0.0_wp_ .and. a(jm) > h .or. &
|
||||
ah > 0.0_wp_ .and. a(jm) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=-j
|
||||
end if
|
||||
end do
|
||||
|
||||
do jm=1,mxr,nrqmax
|
||||
j = jm + nrqmax
|
||||
if (a(j) <= h .and. a(jm) > h .or. &
|
||||
a(j) > h .and. a(jm) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix) =-j
|
||||
end if
|
||||
end do
|
||||
|
||||
do j=2,nr
|
||||
if (a(j) <= h .and. a(j-1) > h .or. &
|
||||
a(j) > h .and. a(j-1) <= h) then
|
||||
ix=ix+1
|
||||
lx(ix)=j
|
||||
end if
|
||||
end do
|
||||
|
||||
if(ix<=0) return
|
||||
|
||||
bb: do
|
||||
in=ix
|
||||
jx=lx(in)
|
||||
jfor=0
|
||||
lda=1
|
||||
ldb=2
|
||||
|
||||
do
|
||||
if(jx<0) then
|
||||
jabs=-jx
|
||||
jnb = jabs - nrqmax
|
||||
else
|
||||
jabs=jx
|
||||
jnb=jabs-1
|
||||
end if
|
||||
|
||||
adn=a(jabs)-a(jnb)
|
||||
if(adn/=0) px=(a(jabs)-h)/adn
|
||||
kx = (jabs - 1) / nrqmax
|
||||
ikx = jabs - nrqmax * kx - 1
|
||||
|
||||
if(jx<0) then
|
||||
x = drgrd * ikx
|
||||
y = dzgrd * (kx - px)
|
||||
else
|
||||
x = drgrd * (ikx - px)
|
||||
y = dzgrd * kx
|
||||
end if
|
||||
|
||||
icount = icount + 1
|
||||
rcon(icount) = x + rqgrid(1)
|
||||
zcon(icount) = y + zqgrid(1)
|
||||
mpl= icount
|
||||
itm = 1
|
||||
ja(1,1) = jabs + nrqmax
|
||||
j=1
|
||||
|
||||
if(jx<=0) then
|
||||
ja(1,1) = -jabs-1
|
||||
j=2
|
||||
end if
|
||||
|
||||
ja(2,1) = -ja(1,1)
|
||||
ja(3,1) = -jx + 1 - nrqmax
|
||||
ja(3,2) = -jx
|
||||
ja(j,2) = jabs - nrqmax
|
||||
k= 3-j
|
||||
ja(k,2) = 1-jabs
|
||||
|
||||
if (kx<=0 .or. ikx<=0) then
|
||||
lda=1
|
||||
ldb=lda
|
||||
else if (ikx + 1 - nr >= 0 .and. jx <= 0) then
|
||||
lda=2
|
||||
ldb=lda
|
||||
else if(jfor/=0) then
|
||||
lda=2
|
||||
do i=1,3
|
||||
if(jfor==ja(i,2)) then
|
||||
lda=1
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
ldb=lda
|
||||
end if
|
||||
|
||||
flag1=.false.
|
||||
aa: do k=1,3
|
||||
do l=lda,ldb
|
||||
do i=1,ix
|
||||
if(lx(i)==ja(k,l)) then
|
||||
itm=itm+1
|
||||
inext= i
|
||||
if(jfor/=0) exit aa
|
||||
if(itm > 3) then
|
||||
flag1=.true.
|
||||
exit aa
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do aa
|
||||
|
||||
if(.not.flag1) then
|
||||
lx(in)=0
|
||||
if(itm == 1) exit
|
||||
end if
|
||||
|
||||
jfor=jx
|
||||
jx=lx(inext)
|
||||
in = inext
|
||||
end do
|
||||
|
||||
do
|
||||
if(lx(ix)/=0) then
|
||||
if(mpl>=4) then
|
||||
ncon = ncon + 1
|
||||
npts(ncon) = icount - iclast
|
||||
iclast = icount
|
||||
end if
|
||||
exit
|
||||
end if
|
||||
ix= ix-1
|
||||
if(ix<=0) exit bb
|
||||
end do
|
||||
|
||||
end do bb
|
||||
|
||||
if(mpl >= 4) then
|
||||
ncon = ncon + 1
|
||||
npts(ncon) = icount - iclast
|
||||
iclast = icount
|
||||
end if
|
||||
end subroutine cniteq_f
|
||||
|
||||
end module
|
||||
Loading…
Reference in New Issue
Block a user