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