2021-12-15 02:31:09 +01:00
|
|
|
|
program main
|
rework the analytical model
This change modifies the analytical equilibrium in order to simplify the
computation of the poloidal flux normalization and the derivatives.
In the power law parametrisation of the safety factor, ρ_t is replaced
with ρ_p and, similarly, the normalised poloidal radius is now
identified with ρ_p, instead of ρ_t.
With the same parameters (q₀,q₁,α...), this choice slightly changes the
plasma current distribution, but enables us to obtain a closed form for
ψ_a = ψ(r=a) and the relation ρ_t(ρ_p). In fact, both expressions are
now obtained by integrating the q(ρ_p), instead of 1/q(ρ_t), which has
no elementary antiderivative.
As the normalisation is now computed exactly, the values of the
normalised flux ψ_n = ψ/ψ_a and the gradient ∇ψ (entering the raytracing
equations in X and ∇X, respectively) are computed to the same precision.
Previously, ψ_n was computed to a lower precision due to the use of a
simple trapezoid integration of 1/q(ρ_p) for ψ_a, while ∇ψ was computed
up to machine precision using an exact formula.
This error effectively caused a very slight decoupling between X=ω_p²/ω²
and ∇X that introduced a systematic error in the numerical solution of
the raytracing equations.
The error manifests itself as a bias with a weak dependency on X in the
values taken by the dispersion function Λ(r̅, n̅) on the phase-space
points generated by the integrator. More specifically,
lim h→0 Λ(r̅_i, n̅_i) = -kX(r̅_i)
where h is the integrator step size;
r̅_i is the position at the i-th step;
k ≈ -3.258⋅10⁻⁵ and depends only on the number of points used to
perform the trapedoid integral for ψ_a (as ~ 1/n²).
After this change Λ behaves consistently with being a conserved quantity
(zero) up to the cumulative integration error of the 4° order
Runge-Kutta method. In fact we now have that:
Λ(r̅_i, n̅_i) ∝ - h⁴ ‖∂⁴X(r̅_i)/∂r̅⁴‖
It must be said that within this model the relation ρ_p(ρ_t) can't be
computed analytically (inverting ρ_t(ρ_p) produces a trascendental
equation of the form b = x + c x^α). However, this relation is not
necessary for raytracing and is easily solved, up to machine
precision, using minpack.
In addition, this change also makes the model consistetly use the
cocos=3 and fully implements the ability to force the signs of I_p, B_φ
(via equilibrium.sgni,sgnb) and rescaling the field (via
equilibrium.factb).
2023-10-11 17:29:24 +02:00
|
|
|
|
use const_and_precisions, only : wp_, zero
|
2021-12-18 18:57:38 +01:00
|
|
|
|
use logger, only : INFO, ERROR, set_log_level, log_message
|
2021-12-22 23:45:30 +01:00
|
|
|
|
use units, only : set_active_units, close_units
|
2022-05-09 21:11:51 +02:00
|
|
|
|
use utils, only : dirname
|
2022-04-29 01:52:50 +02:00
|
|
|
|
use gray_cli, only : cli_options, parse_cli_options, &
|
|
|
|
|
deinit_cli_options, parse_param_overrides
|
2021-12-15 02:31:14 +01:00
|
|
|
|
use gray_core, only : gray_main
|
2021-12-15 02:31:09 +01:00
|
|
|
|
use gray_params, only : gray_parameters, gray_data, gray_results, &
|
2022-05-09 21:11:51 +02:00
|
|
|
|
read_gray_params, read_gray_config, &
|
|
|
|
|
params_set_globals => set_globals
|
2015-11-18 17:34:33 +01:00
|
|
|
|
implicit none
|
2021-12-15 02:30:55 +01:00
|
|
|
|
|
2021-12-15 02:31:16 +01:00
|
|
|
|
! CLI options
|
|
|
|
|
type(cli_options) :: opts
|
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! gray_main subroutine arguments
|
|
|
|
|
type(gray_parameters) :: params ! Inputs
|
|
|
|
|
type(gray_data) :: data !
|
|
|
|
|
type(gray_results) :: results ! Outputs
|
2021-12-18 18:57:38 +01:00
|
|
|
|
integer :: err ! Exit code
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
2022-05-09 21:11:51 +02:00
|
|
|
|
! Store the original working directory
|
|
|
|
|
character(len=256) :: cwd
|
|
|
|
|
call getcwd(cwd)
|
|
|
|
|
|
2021-12-15 02:31:16 +01:00
|
|
|
|
! Parse the command-line options
|
|
|
|
|
call parse_cli_options(opts)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
2021-12-18 18:57:38 +01:00
|
|
|
|
! Initialise logging
|
|
|
|
|
if (opts%quiet) opts%verbose = ERROR
|
|
|
|
|
call set_log_level(opts%verbose)
|
|
|
|
|
|
2021-12-19 16:39:19 +01:00
|
|
|
|
! Activate the given output units
|
|
|
|
|
call set_active_units(opts%units)
|
|
|
|
|
|
2022-05-09 21:11:51 +02:00
|
|
|
|
! Load the parameters from file and move to its directory
|
|
|
|
|
! (all other filepaths are assumed relative to it)
|
|
|
|
|
if (allocated(opts%config_file)) then
|
|
|
|
|
! Use the newer gray.ini configuration file
|
|
|
|
|
call log_message(level=INFO, mod='main', msg='using GRAY INI config')
|
2023-09-12 16:29:09 +02:00
|
|
|
|
call read_gray_config(opts%config_file, params, err)
|
|
|
|
|
if (err /= 0) call exit(1)
|
2022-05-09 21:11:51 +02:00
|
|
|
|
err = chdir(dirname(opts%config_file))
|
|
|
|
|
else
|
|
|
|
|
! Use the legacy gray_params.data file
|
|
|
|
|
call log_message(level=INFO, mod='main', msg='using GRAY params file')
|
2023-09-12 16:29:09 +02:00
|
|
|
|
call read_gray_params(opts%params_file, params, err)
|
|
|
|
|
if (err /= 0) call exit(1)
|
2022-05-09 21:11:51 +02:00
|
|
|
|
err = chdir(dirname(opts%params_file))
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (err /= 0) then
|
|
|
|
|
call log_message(level=ERROR, mod='main', &
|
|
|
|
|
msg='chdir to configuration file directory failed!')
|
|
|
|
|
call exit(1)
|
|
|
|
|
end if
|
|
|
|
|
|
2023-04-02 16:18:43 +02:00
|
|
|
|
! Initialise antenna parameters before parsing the
|
|
|
|
|
! cli arguments, so antenna%pos can be overriden
|
|
|
|
|
call init_antenna(params%antenna, err)
|
|
|
|
|
if (err /= 0) call exit(err)
|
|
|
|
|
|
|
|
|
|
! Apply CLI overrides to the parameters
|
2022-04-29 01:52:50 +02:00
|
|
|
|
call parse_param_overrides(params)
|
2023-04-02 16:18:43 +02:00
|
|
|
|
|
|
|
|
|
! Copy the parameters into global variables
|
|
|
|
|
! exported by the gray_params module
|
2021-12-15 02:31:09 +01:00
|
|
|
|
call params_set_globals(params)
|
|
|
|
|
|
2021-12-19 02:35:24 +01:00
|
|
|
|
! Read the input data and set the global variables
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! of the respective module. Note: order matters.
|
2023-09-12 16:29:09 +02:00
|
|
|
|
call init_equilibrium(params, data, err)
|
|
|
|
|
if (err /= 0) call exit(err)
|
|
|
|
|
|
2022-05-21 22:56:57 +02:00
|
|
|
|
call init_profiles(params%profiles, params%equilibrium%factb, &
|
2023-09-12 16:29:09 +02:00
|
|
|
|
params%antenna%pos, data%profiles, err)
|
|
|
|
|
if (err /= 0) call exit(err)
|
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
call init_misc(params, data)
|
2021-12-15 02:30:55 +01:00
|
|
|
|
|
2023-11-23 18:51:47 +01:00
|
|
|
|
! Get back to the initial directory
|
|
|
|
|
err = chdir(cwd)
|
|
|
|
|
|
2022-05-09 21:11:51 +02:00
|
|
|
|
! Change the current directory to output files there
|
2021-12-15 02:31:16 +01:00
|
|
|
|
if (allocated(opts%output_dir)) then
|
2023-11-23 18:51:47 +01:00
|
|
|
|
if (chdir(opts%output_dir) /= 0) then
|
2021-12-18 18:57:38 +01:00
|
|
|
|
call log_message(level=ERROR, mod='main', &
|
|
|
|
|
msg='chdir to output_dir ('//opts%output_dir//') failed!')
|
2021-12-15 02:31:16 +01:00
|
|
|
|
call exit(1)
|
|
|
|
|
end if
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
if (allocated(opts%sum_filelist)) then
|
2021-12-18 18:57:38 +01:00
|
|
|
|
call log_message(level=INFO, mod='main', msg='summing profiles')
|
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
sum: block
|
2022-11-24 23:37:29 +01:00
|
|
|
|
real(wp_), dimension(:), allocatable :: jphi
|
2021-12-15 02:31:09 +01:00
|
|
|
|
real(wp_), dimension(:), allocatable :: currins, pins, rtin, rpin
|
2022-11-24 23:37:29 +01:00
|
|
|
|
real(wp_), dimension(:), allocatable :: extracol
|
|
|
|
|
integer :: i, j, k, l, n, ngam, nextracol, irt
|
|
|
|
|
character(len=255) :: filename, headerline, fmtstr
|
2021-12-15 02:31:09 +01:00
|
|
|
|
real(wp_), dimension(5) :: f48v
|
2022-11-24 23:37:29 +01:00
|
|
|
|
real(wp_) :: jphip,dpdvp, &
|
2021-12-15 02:31:09 +01:00
|
|
|
|
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))
|
|
|
|
|
|
2021-12-18 18:57:38 +01:00
|
|
|
|
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
|
|
|
|
|
|
2022-11-24 23:37:29 +01:00
|
|
|
|
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))
|
2021-12-15 02:31:09 +01:00
|
|
|
|
do i=1,n
|
|
|
|
|
read(100, *) filename
|
2021-12-18 18:57:38 +01:00
|
|
|
|
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
|
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
do j=1,22
|
2022-11-24 23:37:29 +01:00
|
|
|
|
read(100 + i, '(a255)') headerline
|
|
|
|
|
if (i == 1) then
|
|
|
|
|
write(100 -1, '(a)') trim(headerline)
|
|
|
|
|
write(100 -2, '(a)') trim(headerline)
|
|
|
|
|
end if
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end do
|
2021-12-15 02:30:55 +01:00
|
|
|
|
end do
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
2022-11-24 23:37:29 +01:00
|
|
|
|
allocate(results%jcd(params%output%nrho), results%dpdv(params%output%nrho))
|
2021-12-15 02:31:09 +01:00
|
|
|
|
do k=1,ngam
|
2022-11-24 23:37:29 +01:00
|
|
|
|
jphi = zero
|
|
|
|
|
results%jcd = zero
|
|
|
|
|
results%dpdv = zero
|
|
|
|
|
currins = zero
|
|
|
|
|
pins = zero
|
2021-12-15 02:31:09 +01:00
|
|
|
|
do j=1,params%output%nrho
|
|
|
|
|
do i=1,n
|
2022-11-24 23:37:29 +01:00
|
|
|
|
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)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end do
|
2022-11-24 23:37:29 +01:00
|
|
|
|
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
|
2021-12-15 02:30:55 +01:00
|
|
|
|
end do
|
2022-11-24 23:37:29 +01:00
|
|
|
|
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, &
|
2021-12-15 02:31:09 +01:00
|
|
|
|
drhotjava, drhotpav, ratjamx, ratjbmx
|
2021-12-15 02:30:55 +01:00
|
|
|
|
end do
|
2022-11-24 23:37:29 +01:00
|
|
|
|
do i=-2,n
|
2021-12-15 02:31:09 +01:00
|
|
|
|
close(100 + i)
|
|
|
|
|
end do
|
2022-11-24 23:37:29 +01:00
|
|
|
|
deallocate(jphi, currins, pins, rtin, rpin)
|
|
|
|
|
deallocate(extracol)
|
2022-04-29 01:52:50 +02:00
|
|
|
|
deallocate(opts%params_file)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end block sum
|
2021-12-15 02:30:55 +01:00
|
|
|
|
else
|
2021-12-18 18:57:38 +01:00
|
|
|
|
call gray_main(params, data, results, err)
|
2021-12-15 02:30:55 +01:00
|
|
|
|
end if
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
2021-12-18 18:57:38 +01:00
|
|
|
|
print_res: block
|
|
|
|
|
character(256) :: msg
|
|
|
|
|
write(msg, '(a,g0.3," MW")') 'total absoption: P=', results%pabs
|
|
|
|
|
call log_message(msg, level=INFO, mod='main')
|
|
|
|
|
write(msg, '(a,g0.3," kA")') 'total current drive: I=', results%icd * 1.0e3_wp_
|
|
|
|
|
call log_message(msg, level=INFO, mod='main')
|
|
|
|
|
end block print_res
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! Free memory
|
|
|
|
|
call deinit_equilibrium(data%equilibrium)
|
|
|
|
|
call deinit_profiles(data%profiles)
|
|
|
|
|
call deinit_misc
|
2022-04-29 01:52:50 +02:00
|
|
|
|
call deinit_cli_options(opts)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
deallocate(results%dpdv, results%jcd)
|
2021-12-22 23:45:30 +01:00
|
|
|
|
call close_units
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
contains
|
|
|
|
|
|
2023-09-12 16:29:09 +02:00
|
|
|
|
subroutine init_equilibrium(params, data, err)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! Reads the MHD equilibrium file (either in the G-EQDSK format
|
|
|
|
|
! or an analytical description) and initialises the respective
|
|
|
|
|
! GRAY parameters and data.
|
2024-01-29 02:09:27 +01:00
|
|
|
|
use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
|
2021-12-15 02:31:09 +01:00
|
|
|
|
use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, &
|
2022-12-18 14:09:40 +01:00
|
|
|
|
set_equil_an, set_equil_spline, scale_equil
|
rework the analytical model
This change modifies the analytical equilibrium in order to simplify the
computation of the poloidal flux normalization and the derivatives.
In the power law parametrisation of the safety factor, ρ_t is replaced
with ρ_p and, similarly, the normalised poloidal radius is now
identified with ρ_p, instead of ρ_t.
With the same parameters (q₀,q₁,α...), this choice slightly changes the
plasma current distribution, but enables us to obtain a closed form for
ψ_a = ψ(r=a) and the relation ρ_t(ρ_p). In fact, both expressions are
now obtained by integrating the q(ρ_p), instead of 1/q(ρ_t), which has
no elementary antiderivative.
As the normalisation is now computed exactly, the values of the
normalised flux ψ_n = ψ/ψ_a and the gradient ∇ψ (entering the raytracing
equations in X and ∇X, respectively) are computed to the same precision.
Previously, ψ_n was computed to a lower precision due to the use of a
simple trapezoid integration of 1/q(ρ_p) for ψ_a, while ∇ψ was computed
up to machine precision using an exact formula.
This error effectively caused a very slight decoupling between X=ω_p²/ω²
and ∇X that introduced a systematic error in the numerical solution of
the raytracing equations.
The error manifests itself as a bias with a weak dependency on X in the
values taken by the dispersion function Λ(r̅, n̅) on the phase-space
points generated by the integrator. More specifically,
lim h→0 Λ(r̅_i, n̅_i) = -kX(r̅_i)
where h is the integrator step size;
r̅_i is the position at the i-th step;
k ≈ -3.258⋅10⁻⁵ and depends only on the number of points used to
perform the trapedoid integral for ψ_a (as ~ 1/n²).
After this change Λ behaves consistently with being a conserved quantity
(zero) up to the cumulative integration error of the 4° order
Runge-Kutta method. In fact we now have that:
Λ(r̅_i, n̅_i) ∝ - h⁴ ‖∂⁴X(r̅_i)/∂r̅⁴‖
It must be said that within this model the relation ρ_p(ρ_t) can't be
computed analytically (inverting ρ_t(ρ_p) produces a trascendental
equation of the form b = x + c x^α). However, this relation is not
necessary for raytracing and is easily solved, up to machine
precision, using minpack.
In addition, this change also makes the model consistetly use the
cocos=3 and fully implements the ability to force the signs of I_p, B_φ
(via equilibrium.sgni,sgnb) and rescaling the field (via
equilibrium.factb).
2023-10-11 17:29:24 +02:00
|
|
|
|
use logger, only : log_debug
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! subroutine arguments
|
|
|
|
|
type(gray_parameters), intent(inout) :: params
|
|
|
|
|
type(gray_data), intent(out) :: data
|
2023-09-12 16:29:09 +02:00
|
|
|
|
integer, intent(out) :: err
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
2024-01-29 02:09:27 +01:00
|
|
|
|
select case (params%equilibrium%iequil)
|
|
|
|
|
case (EQ_VACUUM)
|
|
|
|
|
call log_debug('vacuum, no MHD equilibrium', &
|
|
|
|
|
mod='main', proc='init_equilibrium')
|
|
|
|
|
|
|
|
|
|
case (EQ_ANALYTICAL)
|
|
|
|
|
call log_debug('loading analytical file', &
|
|
|
|
|
mod='main', proc='init_equilibrium')
|
|
|
|
|
call read_equil_an(params%equilibrium%filenm, &
|
|
|
|
|
params%raytracing%ipass, &
|
|
|
|
|
data%equilibrium, err)
|
|
|
|
|
if (err /= 0) return
|
|
|
|
|
|
|
|
|
|
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
|
|
|
|
|
call log_debug('loading G-EQDK file', &
|
|
|
|
|
mod='main', proc='init_equilibrium')
|
|
|
|
|
call read_eqdsk(params%equilibrium, data%equilibrium, err)
|
|
|
|
|
if (err /= 0) return
|
|
|
|
|
call change_cocos(data%equilibrium, params%equilibrium%icocos, 3)
|
|
|
|
|
end select
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! Rescale B, I and/or force their signs
|
2022-12-18 14:09:40 +01:00
|
|
|
|
call scale_equil(params%equilibrium, data%equilibrium)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
2024-01-29 02:09:27 +01:00
|
|
|
|
! Set global variables
|
|
|
|
|
select case (params%equilibrium%iequil)
|
2024-01-30 10:14:24 +01:00
|
|
|
|
case (EQ_ANALYTICAL)
|
2024-01-29 02:09:27 +01:00
|
|
|
|
call set_equil_an
|
|
|
|
|
|
|
|
|
|
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
|
|
|
|
|
call log_debug('computing splines...', mod='main', proc='init_equilibrium')
|
|
|
|
|
call set_equil_spline(params%equilibrium, data%equilibrium, err)
|
|
|
|
|
if (err /= 0) return
|
|
|
|
|
call log_debug('splines computed', mod='main', proc='init_equilibrium')
|
|
|
|
|
end select
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end subroutine init_equilibrium
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
subroutine deinit_equilibrium(data)
|
|
|
|
|
! Free all memory allocated by the init_equilibrium subroutine.
|
|
|
|
|
use gray_params, only : equilibrium_data
|
2022-12-18 14:09:40 +01:00
|
|
|
|
use equilibrium, only : unset_equil_spline
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! subroutine arguments
|
|
|
|
|
type(equilibrium_data), intent(inout) :: data
|
|
|
|
|
|
|
|
|
|
! Free the MHD equilibrium arrays
|
|
|
|
|
if (allocated(data%rv)) deallocate(data%rv, data%zv, data%fpol, data%qpsi)
|
|
|
|
|
if (allocated(data%psin)) deallocate(data%psin, data%psinr)
|
|
|
|
|
if (allocated(data%rbnd)) deallocate(data%rbnd, data%zbnd)
|
|
|
|
|
if (allocated(data%rlim)) deallocate(data%rlim, data%zlim)
|
|
|
|
|
|
|
|
|
|
! Unset global variables of the `equilibrium` module
|
2022-05-11 00:30:34 +02:00
|
|
|
|
call unset_equil_spline
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end subroutine deinit_equilibrium
|
|
|
|
|
|
|
|
|
|
|
2023-09-12 16:29:09 +02:00
|
|
|
|
subroutine init_profiles(params, factb, launch_pos, data, err)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! Reads the plasma kinetic profiles file (containing the elecron
|
|
|
|
|
! temperature, density and plasma effective charge) and initialises
|
|
|
|
|
! the respective GRAY data structure.
|
2024-01-29 02:09:27 +01:00
|
|
|
|
use gray_params, only : profiles_parameters, profiles_data, &
|
|
|
|
|
RHO_TOR, RHO_POL, RHO_PSI, &
|
|
|
|
|
PROF_ANALYTIC, PROF_NUMERIC
|
rework the analytical model
This change modifies the analytical equilibrium in order to simplify the
computation of the poloidal flux normalization and the derivatives.
In the power law parametrisation of the safety factor, ρ_t is replaced
with ρ_p and, similarly, the normalised poloidal radius is now
identified with ρ_p, instead of ρ_t.
With the same parameters (q₀,q₁,α...), this choice slightly changes the
plasma current distribution, but enables us to obtain a closed form for
ψ_a = ψ(r=a) and the relation ρ_t(ρ_p). In fact, both expressions are
now obtained by integrating the q(ρ_p), instead of 1/q(ρ_t), which has
no elementary antiderivative.
As the normalisation is now computed exactly, the values of the
normalised flux ψ_n = ψ/ψ_a and the gradient ∇ψ (entering the raytracing
equations in X and ∇X, respectively) are computed to the same precision.
Previously, ψ_n was computed to a lower precision due to the use of a
simple trapezoid integration of 1/q(ρ_p) for ψ_a, while ∇ψ was computed
up to machine precision using an exact formula.
This error effectively caused a very slight decoupling between X=ω_p²/ω²
and ∇X that introduced a systematic error in the numerical solution of
the raytracing equations.
The error manifests itself as a bias with a weak dependency on X in the
values taken by the dispersion function Λ(r̅, n̅) on the phase-space
points generated by the integrator. More specifically,
lim h→0 Λ(r̅_i, n̅_i) = -kX(r̅_i)
where h is the integrator step size;
r̅_i is the position at the i-th step;
k ≈ -3.258⋅10⁻⁵ and depends only on the number of points used to
perform the trapedoid integral for ψ_a (as ~ 1/n²).
After this change Λ behaves consistently with being a conserved quantity
(zero) up to the cumulative integration error of the 4° order
Runge-Kutta method. In fact we now have that:
Λ(r̅_i, n̅_i) ∝ - h⁴ ‖∂⁴X(r̅_i)/∂r̅⁴‖
It must be said that within this model the relation ρ_p(ρ_t) can't be
computed analytically (inverting ρ_t(ρ_p) produces a trascendental
equation of the form b = x + c x^α). However, this relation is not
necessary for raytracing and is easily solved, up to machine
precision, using minpack.
In addition, this change also makes the model consistetly use the
cocos=3 and fully implements the ability to force the signs of I_p, B_φ
(via equilibrium.sgni,sgnb) and rescaling the field (via
equilibrium.factb).
2023-10-11 17:29:24 +02:00
|
|
|
|
use equilibrium, only : frhopol
|
2021-12-15 02:31:10 +01:00
|
|
|
|
use coreprofiles, only : read_profiles_an, read_profiles, &
|
2022-05-11 00:30:34 +02:00
|
|
|
|
scale_profiles, set_profiles_an, &
|
|
|
|
|
set_profiles_spline
|
2022-05-26 02:33:55 +02:00
|
|
|
|
use logger, only : log_debug
|
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! subroutine arguments
|
2022-05-21 22:56:57 +02:00
|
|
|
|
type(profiles_parameters), intent(inout) :: params
|
|
|
|
|
real(wp_), intent(in) :: factb
|
|
|
|
|
real(wp_), intent(in) :: launch_pos(3)
|
|
|
|
|
type(profiles_data), intent(out) :: data
|
2023-09-12 16:29:09 +02:00
|
|
|
|
integer, intent(out) :: err
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
rework the analytical model
This change modifies the analytical equilibrium in order to simplify the
computation of the poloidal flux normalization and the derivatives.
In the power law parametrisation of the safety factor, ρ_t is replaced
with ρ_p and, similarly, the normalised poloidal radius is now
identified with ρ_p, instead of ρ_t.
With the same parameters (q₀,q₁,α...), this choice slightly changes the
plasma current distribution, but enables us to obtain a closed form for
ψ_a = ψ(r=a) and the relation ρ_t(ρ_p). In fact, both expressions are
now obtained by integrating the q(ρ_p), instead of 1/q(ρ_t), which has
no elementary antiderivative.
As the normalisation is now computed exactly, the values of the
normalised flux ψ_n = ψ/ψ_a and the gradient ∇ψ (entering the raytracing
equations in X and ∇X, respectively) are computed to the same precision.
Previously, ψ_n was computed to a lower precision due to the use of a
simple trapezoid integration of 1/q(ρ_p) for ψ_a, while ∇ψ was computed
up to machine precision using an exact formula.
This error effectively caused a very slight decoupling between X=ω_p²/ω²
and ∇X that introduced a systematic error in the numerical solution of
the raytracing equations.
The error manifests itself as a bias with a weak dependency on X in the
values taken by the dispersion function Λ(r̅, n̅) on the phase-space
points generated by the integrator. More specifically,
lim h→0 Λ(r̅_i, n̅_i) = -kX(r̅_i)
where h is the integrator step size;
r̅_i is the position at the i-th step;
k ≈ -3.258⋅10⁻⁵ and depends only on the number of points used to
perform the trapedoid integral for ψ_a (as ~ 1/n²).
After this change Λ behaves consistently with being a conserved quantity
(zero) up to the cumulative integration error of the 4° order
Runge-Kutta method. In fact we now have that:
Λ(r̅_i, n̅_i) ∝ - h⁴ ‖∂⁴X(r̅_i)/∂r̅⁴‖
It must be said that within this model the relation ρ_p(ρ_t) can't be
computed analytically (inverting ρ_t(ρ_p) produces a trascendental
equation of the form b = x + c x^α). However, this relation is not
necessary for raytracing and is easily solved, up to machine
precision, using minpack.
In addition, this change also makes the model consistetly use the
cocos=3 and fully implements the ability to force the signs of I_p, B_φ
(via equilibrium.sgni,sgnb) and rescaling the field (via
equilibrium.factb).
2023-10-11 17:29:24 +02:00
|
|
|
|
! local variables
|
|
|
|
|
integer :: i
|
|
|
|
|
|
2024-01-29 02:09:27 +01:00
|
|
|
|
select case (params%iprof)
|
|
|
|
|
case (PROF_ANALYTIC)
|
|
|
|
|
call log_debug('loading analytical file', &
|
|
|
|
|
mod='main', proc='init_profiles')
|
|
|
|
|
call read_profiles_an(params%filenm, data, err)
|
|
|
|
|
if (err /= 0) return
|
|
|
|
|
|
|
|
|
|
case (PROF_NUMERIC)
|
|
|
|
|
call log_debug('loading numerical file', &
|
|
|
|
|
mod='main', proc='init_profiles')
|
|
|
|
|
call read_profiles(params%filenm, data, err)
|
|
|
|
|
if (err /= 0) return
|
|
|
|
|
|
|
|
|
|
! Convert psrad to ψ
|
|
|
|
|
select case (params%irho)
|
|
|
|
|
case (RHO_TOR) ! psrad is ρ_t = √Φ (toroidal flux)
|
|
|
|
|
data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))]
|
|
|
|
|
case (RHO_POL) ! psrad is ρ_p = √ψ (poloidal flux)
|
|
|
|
|
data%psrad = data%psrad**2
|
|
|
|
|
case (RHO_PSI) ! psrad is already ψ
|
|
|
|
|
end select
|
|
|
|
|
end select
|
2021-12-15 02:31:10 +01:00
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! Rescale input data
|
2021-12-15 02:31:10 +01:00
|
|
|
|
call scale_profiles(params, factb, data)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! Set global variables
|
2024-01-29 02:09:27 +01:00
|
|
|
|
select case (params%iprof)
|
|
|
|
|
case (PROF_ANALYTIC)
|
|
|
|
|
call set_profiles_an(params, data)
|
|
|
|
|
case (PROF_NUMERIC)
|
|
|
|
|
call log_debug('computing splines...', mod='main', proc='init_profiles')
|
|
|
|
|
call set_profiles_spline(params, data, err, launch_pos)
|
|
|
|
|
call log_debug('splines computed', mod='main', proc='init_profiles')
|
|
|
|
|
end select
|
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end subroutine init_profiles
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
subroutine deinit_profiles(data)
|
|
|
|
|
! Free all memory allocated by the init_profiles subroutine.
|
|
|
|
|
use gray_params, only : profiles_data
|
2022-05-11 00:30:34 +02:00
|
|
|
|
use coreprofiles, only : unset_profiles_spline
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! subroutine arguments
|
|
|
|
|
type(profiles_data), intent(inout) :: data
|
|
|
|
|
|
|
|
|
|
! Free the plasma kinetic profiles arrays
|
|
|
|
|
if (allocated(data%psrad)) deallocate(data%psrad)
|
|
|
|
|
if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc)
|
|
|
|
|
|
|
|
|
|
! Unset global variables of the `coreprofiles` module
|
2022-05-11 00:30:34 +02:00
|
|
|
|
call unset_profiles_spline
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end subroutine deinit_profiles
|
|
|
|
|
|
|
|
|
|
|
2023-09-12 16:29:09 +02:00
|
|
|
|
subroutine init_antenna(params, err)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! Reads the wave launcher file (containing the wave frequency, launcher
|
|
|
|
|
! position, direction and beam description) and initialises the respective
|
|
|
|
|
! GRAY parameters.
|
|
|
|
|
use beams, only : read_beam0, read_beam1, read_beam2
|
2024-01-29 02:09:27 +01:00
|
|
|
|
use gray_params, only : antenna_parameters, BEAM_0D, BEAM_1D, BEAM_2D
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! subroutine arguments
|
|
|
|
|
type(antenna_parameters), intent(inout) :: params
|
2023-09-12 16:29:09 +02:00
|
|
|
|
integer, intent(out) :: err
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! Note: α, β are loaded from gray_params.data
|
|
|
|
|
select case (params%ibeam)
|
2024-01-29 02:09:27 +01:00
|
|
|
|
case (BEAM_2D)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! 2 degrees of freedom
|
|
|
|
|
! w(z, α, β), 1/R(z, α, β)
|
|
|
|
|
! FIXME: 1st beam is always selected, iox read from table
|
2023-09-12 16:29:09 +02:00
|
|
|
|
call read_beam2(params, beamid=1, err=err)
|
2024-01-29 02:09:27 +01:00
|
|
|
|
case (BEAM_1D)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! 1 degree of freedom
|
|
|
|
|
! w(z, α), 1/R(z, α)
|
2023-09-12 16:29:09 +02:00
|
|
|
|
call read_beam1(params, err)
|
2024-01-29 02:09:27 +01:00
|
|
|
|
case (BEAM_0D)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
! fixed w(z), 1/R(z)
|
2023-09-12 16:29:09 +02:00
|
|
|
|
call read_beam0(params, err)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end select
|
|
|
|
|
end subroutine init_antenna
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
subroutine init_misc(params, data)
|
|
|
|
|
! Performs miscellanous initial tasks, before the gray_main subroutine.
|
2024-01-29 02:09:27 +01:00
|
|
|
|
use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, &
|
|
|
|
|
EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
|
2023-10-04 17:00:45 +02:00
|
|
|
|
use utils, only : range2rect
|
2023-09-20 15:12:17 +02:00
|
|
|
|
use limiter, only : limiter_set_globals=>set_globals
|
2024-01-30 10:14:24 +01:00
|
|
|
|
use const_and_precisions, only : cm, comp_huge
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! subroutine arguments
|
|
|
|
|
type(gray_parameters), intent(inout) :: params
|
|
|
|
|
type(gray_data), intent(inout) :: data
|
|
|
|
|
|
2023-09-20 15:12:17 +02:00
|
|
|
|
! Build a basic limiter if one is not provided by equilibrium.txt
|
2021-12-15 02:31:09 +01:00
|
|
|
|
if (.not. allocated(data%equilibrium%rlim) &
|
|
|
|
|
.or. params%raytracing%ipass < 0) then
|
2023-09-20 15:12:17 +02:00
|
|
|
|
|
|
|
|
|
! Restore sign of ipass
|
|
|
|
|
params%raytracing%ipass = abs(params%raytracing%ipass)
|
|
|
|
|
|
|
|
|
|
allocate(data%equilibrium%rlim(5))
|
|
|
|
|
allocate(data%equilibrium%zlim(5))
|
2021-12-15 02:31:09 +01:00
|
|
|
|
block
|
2023-09-20 15:12:17 +02:00
|
|
|
|
real(wp_) :: R_launch, z_launch, R_max, delta_z
|
|
|
|
|
|
|
|
|
|
! Launcher coordinates
|
|
|
|
|
R_launch = norm2(params%antenna%pos(1:2)) * cm
|
|
|
|
|
z_launch = params%antenna%pos(3) * cm
|
|
|
|
|
|
|
|
|
|
! Max vertical distance a ray can travel
|
|
|
|
|
delta_z = abs(params%raytracing%ipass) * &
|
|
|
|
|
params%raytracing%dst * params%raytracing%nstep * cm
|
|
|
|
|
|
|
|
|
|
! Max radius, either due to the plasma extent or equilibrium grid
|
2024-01-29 02:09:27 +01:00
|
|
|
|
select case (params%equilibrium%iequil)
|
2024-01-30 10:14:24 +01:00
|
|
|
|
case (EQ_VACUUM)
|
|
|
|
|
! Use a very large R, ~ unbounded
|
|
|
|
|
R_max = comp_huge
|
|
|
|
|
|
|
|
|
|
case (EQ_ANALYTICAL)
|
2024-01-29 02:09:27 +01:00
|
|
|
|
! use R₀+a
|
|
|
|
|
block
|
|
|
|
|
use equilibrium, only : model
|
|
|
|
|
R_max = model%R0 + model%a
|
|
|
|
|
end block
|
|
|
|
|
|
|
|
|
|
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
|
|
|
|
|
! use max R of the grid
|
|
|
|
|
R_max = data%equilibrium%rv(size(data%equilibrium%rv))
|
|
|
|
|
end select
|
2023-09-20 15:12:17 +02:00
|
|
|
|
|
|
|
|
|
! Avoid clipping out the launcher
|
|
|
|
|
R_max = max(R_launch, R_max)
|
|
|
|
|
|
|
|
|
|
! Convert to a list of R,z:
|
|
|
|
|
!
|
|
|
|
|
! (R_wall, z_launch+Δz)╔═════╗(R_max, z_launch+Δz)
|
|
|
|
|
! ║4 3║
|
|
|
|
|
! ║ ║
|
|
|
|
|
! ║1 2║
|
|
|
|
|
! (R_wall, z_launch-Δz)╚═════╝(R_max, z_launch-Δz)
|
|
|
|
|
!
|
|
|
|
|
call range2rect(xmin=params%misc%rwall, xmax=R_max, &
|
|
|
|
|
ymin=z_launch - delta_z, ymax=z_launch + delta_z, &
|
2023-10-04 17:00:45 +02:00
|
|
|
|
x=data%equilibrium%rlim, y=data%equilibrium%zlim)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
end block
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
! Set the global variables of the `limiter` module
|
|
|
|
|
call limiter_set_globals(data%equilibrium)
|
|
|
|
|
end subroutine init_misc
|
|
|
|
|
|
2022-04-19 23:09:24 +02:00
|
|
|
|
|
2021-12-15 02:31:09 +01:00
|
|
|
|
subroutine deinit_misc
|
|
|
|
|
! Free all memory allocated by the init_misc subroutine.
|
|
|
|
|
use limiter, only : limiter_unset_globals=>unset_globals
|
|
|
|
|
|
|
|
|
|
! Unset the global variables of the `limiter` module
|
|
|
|
|
call limiter_unset_globals
|
|
|
|
|
end subroutine deinit_misc
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
subroutine sum_profiles(params, jphi, jcd, dpdv, currins, pins, pabs, icd, &
|
|
|
|
|
jphip, dpdvp, rhotj, rhotjava, rhotp, rhotpav, &
|
|
|
|
|
drhotjava, drhotpav, ratjamx, ratjbmx)
|
|
|
|
|
use const_and_precisions, only : zero, degree
|
2022-05-11 00:30:34 +02:00
|
|
|
|
use coreprofiles, only : set_profiles_an, set_profiles_spline, &
|
|
|
|
|
temp, fzeff
|
2021-12-15 02:31:09 +01:00
|
|
|
|
use dispersion, only : expinit
|
2022-05-10 12:36:37 +02:00
|
|
|
|
use gray_params, only : gray_parameters, print_parameters
|
2021-12-15 02:31:09 +01:00
|
|
|
|
use beams, only : launchangles2n, xgygcoeff
|
|
|
|
|
use magsurf_data, only : flux_average, dealloc_surfvec
|
2024-04-23 17:00:40 +02:00
|
|
|
|
use beamdata, only : init_btr
|
2021-12-15 02:31:09 +01:00
|
|
|
|
use pec, only : pec_init, postproc_profiles, dealloc_pec, &
|
|
|
|
|
rhop_tab, rhot_tab
|
2021-12-15 02:31:14 +01:00
|
|
|
|
use gray_core, only : print_headers, print_finals, print_pec, &
|
2021-12-15 02:31:09 +01:00
|
|
|
|
print_bres, print_prof, print_maps, &
|
|
|
|
|
print_surfq
|
|
|
|
|
|
|
|
|
|
! subroutine arguments
|
2022-05-10 12:36:37 +02:00
|
|
|
|
type(gray_parameters), intent(inout) :: params
|
|
|
|
|
real(wp_), intent(in) :: pabs, icd
|
|
|
|
|
real(wp_), dimension(:), intent(in) :: jphi, jcd, dpdv, currins, pins
|
2021-12-15 02:31:09 +01:00
|
|
|
|
real(wp_), intent(out) :: jphip, dpdvp, rhotj, rhotjava, &
|
|
|
|
|
rhotp, rhotpav, drhotjava, drhotpav, &
|
2022-05-10 12:36:37 +02:00
|
|
|
|
ratjamx, ratjbmx
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! local variables
|
|
|
|
|
real(wp_) :: ak0, bres, xgcn
|
|
|
|
|
real(wp_) :: chipol, psipol, st
|
|
|
|
|
real(wp_) :: drhotp, drhotj, dpdvmx, jphimx
|
|
|
|
|
|
|
|
|
|
! ======== set environment BEGIN ========
|
2022-04-19 23:09:24 +02:00
|
|
|
|
! Compute X=(ω_pe/ω)² and Y=ω_ce/ω (with B=1)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn)
|
|
|
|
|
|
|
|
|
|
! Initialise the ray variables (beamtracing)
|
2024-04-23 17:00:40 +02:00
|
|
|
|
call init_btr(params%raytracing)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! Initialise the dispersion module
|
2022-05-11 00:30:34 +02:00
|
|
|
|
if (params%ecrh_cd%iwarm > 1) call expinit
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! Initialise the magsurf_data module
|
|
|
|
|
call flux_average ! requires frhotor for dadrhot,dvdrhot
|
|
|
|
|
|
|
|
|
|
! Initialise the output profiles
|
|
|
|
|
call pec_init(params%output%ipec)
|
|
|
|
|
! ======= set environment END ======
|
|
|
|
|
|
|
|
|
|
! ======== pre-proc prints BEGIN ========
|
2022-05-03 23:16:21 +02:00
|
|
|
|
call print_headers(params)
|
2021-12-15 02:31:09 +01:00
|
|
|
|
|
|
|
|
|
! Print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout
|
|
|
|
|
call print_surfq([1.5_wp_, 2.0_wp_])
|
|
|
|
|
|
|
|
|
|
! Print ne, Te, q, Jphi versus psi, rhop, rhot
|
|
|
|
|
call print_bres(bres)
|
2023-05-23 17:49:51 +02:00
|
|
|
|
call print_prof(params%profiles)
|
2021-12-19 02:35:24 +01:00
|
|
|
|
call print_maps(bres, xgcn, &
|
|
|
|
|
norm2(params%antenna%pos(1:2)) *0.01_wp_ , &
|
2021-12-15 02:31:09 +01:00
|
|
|
|
sin(params%antenna%beta*degree))
|
|
|
|
|
! ========= pre-proc prints END =========
|
|
|
|
|
|
|
|
|
|
! Print power and current density profiles
|
|
|
|
|
call print_pec(rhop_tab, rhot_tab, jphi, jcd, &
|
|
|
|
|
dpdv, currins, pins, index_rt=1)
|
|
|
|
|
! Compute profiles width
|
|
|
|
|
call postproc_profiles(pabs, icd, rhot_tab, dpdv, jphi, &
|
|
|
|
|
rhotpav, drhotpav, rhotjava, drhotjava, &
|
|
|
|
|
dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, &
|
|
|
|
|
dpdvmx, jphimx, ratjamx, ratjbmx)
|
|
|
|
|
! Print 0D results
|
|
|
|
|
call print_finals(pabs, icd, dpdvp, jphip, rhotpav, rhotjava, drhotpav, &
|
|
|
|
|
drhotjava, dpdvmx, jphimx, rhotp, rhotj, drhotp, &
|
|
|
|
|
drhotj, ratjamx, ratjbmx, st, psipol, chipol, &
|
|
|
|
|
1, params%antenna%power, cpl1=zero, cpl2=zero)
|
|
|
|
|
|
|
|
|
|
! Free memory
|
|
|
|
|
call dealloc_surfvec ! for fluxval
|
|
|
|
|
call dealloc_pec
|
|
|
|
|
end subroutine sum_profiles
|
|
|
|
|
|
|
|
|
|
end program main
|