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

src/main.f90: use enums

Browse Source
This commit is contained in:
Michele Guerini Rocco 2024-01-29 02:09:27 +01:00
parent c5a4b180bc
commit fac0c6ded8
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
1 changed files with 83 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

77
src/main.f90
View File

@ -206,6 +206,7 @@ contains
! Reads the MHD equilibrium file (either in the G-EQDSK format ! Reads the MHD equilibrium file (either in the G-EQDSK format
! or an analytical description) and initialises the respective ! or an analytical description) and initialises the respective
! GRAY parameters and data. ! GRAY parameters and data.
use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, & use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, &
set_equil_an, set_equil_spline, scale_equil set_equil_an, set_equil_spline, scale_equil
use logger, only : log_debug use logger, only : log_debug
@ -215,35 +216,41 @@ contains
type(gray_data), intent(out) :: data type(gray_data), intent(out) :: data
integer, intent(out) :: err integer, intent(out) :: err
if (params%equilibrium%iequil < 2) then select case (params%equilibrium%iequil)
! Analytical equilibrium case (EQ_VACUUM)
call log_debug('vacuum, no MHD equilibrium', &
mod='main', proc='init_equilibrium')
case (EQ_ANALYTICAL)
call log_debug('loading analytical file', & call log_debug('loading analytical file', &
mod='main', proc='init_equilibrium') mod='main', proc='init_equilibrium')
call read_equil_an(params%equilibrium%filenm, & call read_equil_an(params%equilibrium%filenm, &
params%raytracing%ipass, & params%raytracing%ipass, &
data%equilibrium, err) data%equilibrium, err)
if (err /= 0) return if (err /= 0) return
else
! Numerical equilibrium case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
call log_debug('loading G-EQDK file', & call log_debug('loading G-EQDK file', &
mod='main', proc='init_equilibrium') mod='main', proc='init_equilibrium')
call read_eqdsk(params%equilibrium, data%equilibrium, err) call read_eqdsk(params%equilibrium, data%equilibrium, err)
if (err /= 0) return if (err /= 0) return
call change_cocos(data%equilibrium, params%equilibrium%icocos, 3) call change_cocos(data%equilibrium, params%equilibrium%icocos, 3)
end if end select
! Rescale B, I and/or force their signs ! Rescale B, I and/or force their signs
call scale_equil(params%equilibrium, data%equilibrium) call scale_equil(params%equilibrium, data%equilibrium)
! Set global variables (for splines) ! Set global variables
if (params%equilibrium%iequil < 2) then select case (params%equilibrium%iequil)
case (EQ_VACUUM, EQ_ANALYTICAL)
call set_equil_an call set_equil_an
else
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
call log_debug('computing splines...', mod='main', proc='init_equilibrium') call log_debug('computing splines...', mod='main', proc='init_equilibrium')
call set_equil_spline(params%equilibrium, data%equilibrium, err) call set_equil_spline(params%equilibrium, data%equilibrium, err)
if (err /= 0) return if (err /= 0) return
call log_debug('splines computed', mod='main', proc='init_equilibrium') call log_debug('splines computed', mod='main', proc='init_equilibrium')
end if end select
end subroutine init_equilibrium end subroutine init_equilibrium
@ -270,7 +277,9 @@ contains
! Reads the plasma kinetic profiles file (containing the elecron ! Reads the plasma kinetic profiles file (containing the elecron
! temperature, density and plasma effective charge) and initialises ! temperature, density and plasma effective charge) and initialises
! the respective GRAY data structure. ! the respective GRAY data structure.
use gray_params, only : profiles_parameters, profiles_data use gray_params, only : profiles_parameters, profiles_data, &
RHO_TOR, RHO_POL, RHO_PSI, &
PROF_ANALYTIC, PROF_NUMERIC
use equilibrium, only : frhopol use equilibrium, only : frhopol
use coreprofiles, only : read_profiles_an, read_profiles, & use coreprofiles, only : read_profiles_an, read_profiles, &
scale_profiles, set_profiles_an, & scale_profiles, set_profiles_an, &
@ -287,14 +296,14 @@ contains
! local variables ! local variables
integer :: i integer :: i
if (params%iprof == 0) then select case (params%iprof)
! Analytical profiles case (PROF_ANALYTIC)
call log_debug('loading analytical file', & call log_debug('loading analytical file', &
mod='main', proc='init_profiles') mod='main', proc='init_profiles')
call read_profiles_an(params%filenm, data, err) call read_profiles_an(params%filenm, data, err)
if (err /= 0) return if (err /= 0) return
else
! Numerical profiles case (PROF_NUMERIC)
call log_debug('loading numerical file', & call log_debug('loading numerical file', &
mod='main', proc='init_profiles') mod='main', proc='init_profiles')
call read_profiles(params%filenm, data, err) call read_profiles(params%filenm, data, err)
@ -302,27 +311,27 @@ contains
! Convert psrad to ψ ! Convert psrad to ψ
select case (params%irho) select case (params%irho)
case (0) ! psrad is ρ_t = √Φ (toroidal flux) case (RHO_TOR) ! psrad is ρ_t = √Φ (toroidal flux)
data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))] data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))]
case (1) ! psrad is ρ_p = √ψ (poloidal flux) case (RHO_POL) ! psrad is ρ_p = √ψ (poloidal flux)
data%psrad = data%psrad**2 data%psrad = data%psrad**2
case default ! psrad is already ψ case (RHO_PSI) ! psrad is already ψ
end select
end select end select
end if
! Rescale input data ! Rescale input data
call scale_profiles(params, factb, data) call scale_profiles(params, factb, data)
! Set global variables ! Set global variables
if (params%iprof == 0) then select case (params%iprof)
! Analytical profiles case (PROF_ANALYTIC)
call set_profiles_an(params, data) call set_profiles_an(params, data)
else case (PROF_NUMERIC)
! Numerical profiles
call log_debug('computing splines...', mod='main', proc='init_profiles') call log_debug('computing splines...', mod='main', proc='init_profiles')
call set_profiles_spline(params, data, err, launch_pos) call set_profiles_spline(params, data, err, launch_pos)
call log_debug('splines computed', mod='main', proc='init_profiles') call log_debug('splines computed', mod='main', proc='init_profiles')
end if end select
end subroutine init_profiles end subroutine init_profiles
@ -348,7 +357,7 @@ contains
! position, direction and beam description) and initialises the respective ! position, direction and beam description) and initialises the respective
! GRAY parameters. ! GRAY parameters.
use beams, only : read_beam0, read_beam1, read_beam2 use beams, only : read_beam0, read_beam1, read_beam2
use gray_params, only : antenna_parameters use gray_params, only : antenna_parameters, BEAM_0D, BEAM_1D, BEAM_2D
! subroutine arguments ! subroutine arguments
type(antenna_parameters), intent(inout) :: params type(antenna_parameters), intent(inout) :: params
@ -356,16 +365,16 @@ contains
! Note: α, β are loaded from gray_params.data ! Note: α, β are loaded from gray_params.data
select case (params%ibeam) select case (params%ibeam)
case (2) case (BEAM_2D)
! 2 degrees of freedom ! 2 degrees of freedom
! w(z, α, β), 1/R(z, α, β) ! w(z, α, β), 1/R(z, α, β)
! FIXME: 1st beam is always selected, iox read from table ! FIXME: 1st beam is always selected, iox read from table
call read_beam2(params, beamid=1, err=err) call read_beam2(params, beamid=1, err=err)
case (1) case (BEAM_1D)
! 1 degree of freedom ! 1 degree of freedom
! w(z, α), 1/R(z, α) ! w(z, α), 1/R(z, α)
call read_beam1(params, err) call read_beam1(params, err)
case default case (BEAM_0D)
! fixed w(z), 1/R(z) ! fixed w(z), 1/R(z)
call read_beam0(params, err) call read_beam0(params, err)
end select end select
@ -374,6 +383,8 @@ contains
subroutine init_misc(params, data) subroutine init_misc(params, data)
! Performs miscellanous initial tasks, before the gray_main subroutine. ! Performs miscellanous initial tasks, before the gray_main subroutine.
use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, &
EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
use utils, only : range2rect use utils, only : range2rect
use limiter, only : limiter_set_globals=>set_globals use limiter, only : limiter_set_globals=>set_globals
use const_and_precisions, only : cm use const_and_precisions, only : cm
@ -403,16 +414,18 @@ contains
params%raytracing%dst * params%raytracing%nstep * cm params%raytracing%dst * params%raytracing%nstep * cm
! Max radius, either due to the plasma extent or equilibrium grid ! Max radius, either due to the plasma extent or equilibrium grid
if (params%equilibrium%iequil < 2) then select case (params%equilibrium%iequil)
! analytic equilibrium, use R₀+a case (EQ_VACUUM, EQ_ANALYTICAL)
! use R₀+a
block block
use equilibrium, only : model use equilibrium, only : model
R_max = model%R0 + model%a R_max = model%R0 + model%a
end block end block
else
! numeric equilibrium, use max R of the grid case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
! use max R of the grid
R_max = data%equilibrium%rv(size(data%equilibrium%rv)) R_max = data%equilibrium%rv(size(data%equilibrium%rv))
end if end select
! Avoid clipping out the launcher ! Avoid clipping out the launcher
R_max = max(R_launch, R_max) R_max = max(R_launch, R_max)