add INI configuration file

This adds a new configuration file based on the INI format.
The new format will allow adding GRAY parameters without breaking
compatibility with existing configurations, unlike as of the old
gray_params.data.

input/gray.ini

@@ -0,0 +1,190 @@
+[raytracing]
+; Number of rays in the radial/angular direction
+nrayr  = 1
+nrayth = 16
+; Normalized maximum radius of the beam power
+rwmax = 1
+
+; Switch between simple raytracing (0) and beamtracing (1)
+igrad = 0
+
+; Max number of passes inside the plasma  [multipass module]
+; When positive reflections occur on the plasma limiter provided by the
+; G-EQDSK file; when negative on a simple limiter at R=`rwall`, see below.
+ipass = 1
+
+; Whether to compute the wave polarisation
+ipol = 0
+; Step size (cm) for the numerical integration
+dst = 0.1
+; Max number of integration steps
+nstep = 12000
+
+; Choice of the integration variable
+;  0: path length (s)
+;  1: "time" (actually c⋅t)
+;  2: real part of the eikonal (S_r)
+idst = 0
+
+
+[ecrh_cd]
+; Choice of the power absorption model
+;  0: no absorption at all
+;  1: weakly relativistic
+;  2: fully relativistic (faster variant)
+;  3: fully relativistic (slower variant)
+;  4: tenuous plasma (very fast)
+; Note: iwarm>0 is required for current driver
+iwarm = 2
+
+; Order of the electron Larmor radius expansion
+; (used by some absorption models)
+ilarm = 5
+
+; Max number of iterations for the solution of the dispersion relation.
+; (used by some absorption models)
+; Note: if negative the result of the first iteration will be used in
+; case the result doesn't converge within |imx| iterations.
+imx = -20
+
+; Current drive model
+;   0: no current drive at all
+;   1: Cohen
+;   2: no trapping
+;   3: Neoclassical
+ieccd = 3
+
+
+[antenna]
+; Wave launch angles (deg)
+alpha = 45 ; Poloidal angle (positive → up)
+beta  = 0  ; Toroidal angle (positive → right)
+
+; Injected power (MW)
+power = 1
+
+; Polarisation mode
+;   1: ordinary (O)
+;   2: extraordinary (X)
+iox = 1
+; Alternatively, parameters of the polarisation ellipse
+;   χ: angle between the principal axes and the (x,y) axes
+;   ψ: atan(ε), where ε is the ellipticity
+chi = 0
+psi = 0
+
+; Beam description kind
+;  0: simple beam shape
+;  1: 1D table
+;  2: 2D table
+ibeam = 0
+
+; Filepath of the beam data (relative to this file)
+filenm = "beamdata.txt"
+
+
+[equilibrium]
+; MHD equilibrium kind
+;   0: vacuum (i.e. no plasma at all)
+;   1: analytical
+;   2: G-EQDSK format - data valid on the whole domain
+;   3: G-EQDSK format - data valid only inside the LCFS
+iequil = 3
+
+; Filepath of the equilibrium data (relative to this file)
+filenm = "magneticdata.eqdsk"
+; COCOS index
+icocos = 0
+
+; Normalisation of the poloidal function
+;   0: G-EQDSK, ψ → |ψ - ψ(edge)|/|ψ(axis) - ψ(edge)|
+;   1: no normalisation
+ipsinorm = 0
+
+; G-EQDSK format parameters
+; Whether header starts with a description, a.k.a identification string
+idesc = 1
+; Whether the records have variable length
+; Note: some non-compliant programs output numbers formatted with variable length
+; instead of using the single (5e16.9) specifier.
+ifreefmt = 0
+
+; Position of the X point
+;   -1: bottom
+;    0: no X point
+;   +1: top
+ixp = 0
+
+; Tension of splines
+; Note: 0 means perfect interpolation
+ssplps = 0.005  ; for ψ(R,Z), normalised poloidal flux
+ssplf = 0.01    ; for I(ψ)=R⋅B_T, poloidal current function
+
+; Sign of toroidal field/current (used when COCOS index is 0,10)
+; When viewing from above: +1 → counter-clockwise, -1 → clockwise
+sgnb = -1
+sgni = +1
+; Rescaling factor for the magnetic field
+factb = 1
+
+
+[profiles]
+; (input) plasma profiles parameters
+
+; Profiles kind
+;   0: analytical
+;   1: numerical
+iprof = 1
+
+; Profile radial coordinate
+;  0: ρ_t = √Φ (where Φ is the normalised toroidal flux)
+;  1: ρ_p = √ψ (where ψ is the normalised poloidal flux)
+;  2: normalised poloidal flux ψ
+irho = 0
+
+; Filepath of the equilibrium (relative to this file)
+filenm = "profiles.txt"
+
+; Value of ψ at the plasma boundary  [multipass module]
+; Notes:
+;   1. boundary means zero density;
+;   2. determines when a ray is considered inside the plasma
+psnbnd = 1.007
+
+; Tension of the density spline
+; Note: 0 means perfect interpolation
+sspld = 0.1
+
+; Rescaling factor for electron temperature/density
+factte = 1
+factne = 1
+
+; Choice of model for rescaling the temperature/density with the
+; magnetic field
+;   1: costant Greenwald density (ab=1)
+;   2: no rescaling (ab=0)
+iscal = 2
+
+
+[output]
+; Output data parameters
+
+; ECRH&CD profiles grid:
+; Radial coordinate
+;  1: ρ_p = √ψ (where ψ is the normalised poloidal flux)
+;  2: ρ_t = √Φ (where Φ is the normalised toroidal flux)
+ipec = 1
+; Number of points
+nrho = 501
+
+; Subsampling factors:
+istpr = 5  ; beam cross section (units 8, 12)
+istpl = 5  ; outer rays data (unit 33)
+
+
+[misc]
+; Other parameters
+
+; Radius of the inner wall (m)  [multipass module]
+; (when ipass<0, used to build a simple limiter for reflections)
+rwall = 1.36

src/gray_cli.f90

@@ -15,6 +15,7 @@ module gray_cli
     ! Files
     character(len=:), allocatable :: output_dir
     character(len=:), allocatable :: params_file
+    character(len=:), allocatable :: config_file
     character(len=:), allocatable :: sum_filelist
     ! others
     integer :: verbose
@@ -47,13 +48,15 @@ contains
     print '(a)', '  -q, --quiet               suppress all non-critical messages'
     print '(a)', '  -o, --output-dir DIR      specify where to write the output files'
     print '(a)', '                            (default: current directory)'
-    print '(a)', '  -p, --params-file FILE    set the parameters file'
+    print '(a)', '  -p, --params-file FILE    set the (legacy) parameters file'
     print '(a)', '                            (default: gray_params.data)'
+    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)', '                            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;'
+    print '(a)', '                            specified via --params-file/--config-file;'
     print '(a)', '                            the ID is GROUP.NAME, ex. antenna.fghz;'
     print '(a)', '                            see the manual for available parameters.'
     print '(a)', ''
@@ -155,6 +158,10 @@ contains
           call get_command_string(i + 1, opts%params_file)
           skip_next = .true.
 
+        case ('-c', '--config-file')
+          call get_command_string(i + 1, opts%config_file)
+          skip_next = .true.
+
         case ('-s', '--sum')
           call get_command_string(i + 1, opts%sum_filelist)
           skip_next = .true.
@@ -201,6 +208,7 @@ contains
     ! Reads GRAY parameters from CLI and overrides `params` accordingly
 
     use gray_params, only : gray_parameters, update_parameter
+    use ini_parser,  only : ERR_VALUE, ERR_UNKNOWN
 
     implicit none
 
@@ -211,7 +219,7 @@ contains
     character(len=:), allocatable :: argument, temp, id, val
     logical :: skip_next = .false.
     integer :: i, nargs
-    integer :: error, sep
+    integer :: sep
 
     nargs = command_argument_count()
     do i = 1, nargs
@@ -242,12 +250,12 @@ contains
 
           ! match the name string to a parameter
           select case (update_parameter(params, id, val))
-            case (1)
+            case (ERR_VALUE)
               print '(4a)', 'invalid value for ', id, ': ', val
               deallocate(temp)
               call exit(1)
 
-            case (2)
+            case (ERR_UNKNOWN)
               print '(a,a)', 'unknown GRAY parameter: ', id
               deallocate(temp)
               call exit(1)
@@ -297,6 +305,7 @@ contains
 
     if (allocated(opts%output_dir))   deallocate(opts%output_dir)
     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)
   end subroutine deinit_cli_options

src/gray_jetto1beam.f90

@@ -31,8 +31,8 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
 
     ! Read parameters from external file
     init_params: block
-      use gray_params, only : read_parameters, set_globals
+      use gray_params, only : read_gray_params, set_globals
-      call read_parameters('gray.data', params)
+      call read_gray_params('gray.data', params)
 
       ! Override some parameters
       params%misc%rwall         = r(1)

src/gray_params.f90

@@ -25,31 +25,46 @@ module gray_params
 
   ! MHD equilibrium parameters
   type equilibrium_parameters
-    real(wp_) :: ssplps, ssplf, factb
+    real(wp_) :: ssplps, ssplf      ! Spline tensions for ψ(R,Z), I(ψ)
-    integer :: sgnb, sgni, ixp
+    real(wp_) :: factb              ! Rescaling factor for B
-    integer :: iequil, icocos, ipsinorm, idesc, ifreefmt
+    integer :: sgnb, sgni           ! Sign of B, I
-    character(len=lenfnm) :: filenm
+    integer :: ixp                  ! Position of X point
+    integer :: iequil               ! Equilibrium kind
+    integer :: icocos               ! COCOS index
+    integer :: ipsinorm             ! Normalisation of ψ
+    integer :: idesc, ifreefmt      ! G-EQDSK quirks
+    character(len=lenfnm) :: filenm ! Equilibrium data filepath
   end type
 
   ! Kinetic plasma profiles
   type profiles_parameters
-    real(wp_) :: psnbnd ! plasma boundary (ψ: ne(ψ)=0)
+    real(wp_) :: psnbnd             ! Plasma boundary (ψ: ne(ψ)=0)
-    real(wp_) :: sspld, factne, factte
+    real(wp_) :: sspld              ! Density spline tension
-    integer :: iscal, irho
+    real(wp_) :: factne, factte     ! Rescale factors for n_e, T_e
-    integer :: iprof
+    integer :: iscal                ! Rescaling model
-    character(len=lenfnm) :: filenm
+    integer :: irho                 ! Radial coordinate
+    integer :: iprof                ! Profile kind
+    character(len=lenfnm) :: filenm ! Profiles data filepath
   end type
 
   ! Raytracing parameters
   type raytracing_parameters
-    real(wp_) :: rwmax, dst
+    real(wp_) :: dst                ! Integration step size
-    integer :: nrayr, nrayth, nstep
+    real(wp_) :: rwmax              ! Normalized maximum radius of beam power
-    integer :: igrad, idst, ipass, ipol
+    integer :: nrayr, nrayth        ! Radial/angular number of rays
+    integer :: igrad                ! Complex eikonal switch
+    integer :: nstep                ! Max number of integration steps
+    integer :: idst                 ! Choice of the integration variable
+    integer :: ipass                ! Number of plasma passes
+    integer :: ipol                 ! Whether to compute wave polarisation
   end type
 
   ! EC resonant heating & Current Drive parameters
   type ecrh_cd_parameters
-    integer :: iwarm, ilarm, imx, ieccd
+    integer :: iwarm  ! Choice of power absorption model
+    integer :: ieccd  ! Choice of current drive model
+    integer :: ilarm  ! Larmor-radius expansion order
+    integer :: imx    ! Max iterations for solving the dispersion relation
   end type
 
   ! Output data parameters
@@ -231,28 +246,71 @@ contains
   function update_parameter(params, name, value) result(error)
     ! Updates the value of a parameter, addressed by a string
     ! The return error is:
-    !   0 on success;
+    !   ERR_SUCCESS on success;
-    !   1 on invalid parameter value;
+    !   ERR_VALUE   on invalid parameter value;
-    !   2 on unknown parameter name.
+    !   ERR_UNKNOWN on unknown parameter name.
     !
     ! Ex. update_parameter(params, 'raytracing.nrayr', '10')
+    use ini_parser, only : ini_error, ERR_SUCCESS, ERR_VALUE, ERR_UNKNOWN
+
     implicit none
 
     ! function arguments
     type(gray_parameters), intent(inout) :: params
-    character(len=:), allocatable, intent(in)    :: name, value
+    character(*),          intent(in)    :: name, value
-    integer :: error
+    integer(kind(ini_error)) :: error
+
+    ! hope for the best...
+    error = ERR_SUCCESS
 
     select case (name)
 #include "gray_params.inc"
       case default
-        error = 2
+        ! unknown parameter
+        error = ERR_UNKNOWN
     end select
 
   end function update_parameter
 
 
-  subroutine read_parameters(filename, params, unit)
+  subroutine read_gray_config(filename, params)
+    ! Reads the GRAY parameters from the gray.ini configuration file
+    use ini_parser, only : parse_ini, property_handler, ini_error, ERR_SUCCESS
+
+    implicit none
+
+    ! subroutine arguments
+    character(len=*),       intent(in) :: filename
+    type(gray_parameters), intent(out) :: params
+
+    ! local variables
+    procedure(property_handler), pointer :: p
+    integer :: error
+
+    p => ini_handler
+    call parse_ini(filename, p, error)
+    if (error /= ERR_SUCCESS) call exit(1)
+
+    contains
+
+      function ini_handler(section, name, value) result(error)
+        ! This function handles a single INI property and updates
+        ! the `params` structure
+
+        implicit none
+
+        ! function arguments
+        character(*), intent(in) :: section, name, value
+        integer(kind(ini_error)) :: error
+
+        error = update_parameter(params, section // "." // name, value)
+      end function ini_handler
+
+  end subroutine read_gray_config
+
+
+  subroutine read_gray_params(filename, params)
+    ! Reads the GRAY parameters from the legacy gray_params.data file
     use utils,  only : get_free_unit
     use logger, only : log_error
 
@@ -261,125 +319,57 @@ contains
     ! subrouting arguments
     character(len=*), intent(in) :: filename
     type(gray_parameters), intent(out) :: params
-    integer, intent(in), optional :: unit
 
     ! local variables
     integer :: u, err
 
-    u = get_free_unit(unit)
+    u = get_free_unit()
 
     open(u, file=filename, status='old', action='read', iostat=err)
     if (err /= 0) then
       call log_error('opening gray_params file ('//filename//') failed!', &
-                     mod='gray_params', proc='read_parameters')
+                     mod='gray_params', proc='read_gray_params')
       call exit(1)
     end if
 
-    ! Raytracing
+    read(u, *) params%raytracing%nrayr, params%raytracing%nrayth, &
-    ! ========================================================================
+               params%raytracing%rwmax
-    ! nrayr      :number of rays in radial direction
+    read(u, *) params%raytracing%igrad, params%raytracing%ipass, &
-    ! nrayth     :number of rays in angular direction
+               params%raytracing%ipol
-    ! rwmax      :normalized maximum radius of beam power
+    read(u, *) params%raytracing%dst, params%raytracing%nstep, &
-    ! rwmax=1 -> :last ray at radius = waist
+               params%raytracing%idst
-    read(u, *) params%raytracing%nrayr, params%raytracing%nrayth, params%raytracing%rwmax
-    ! igrad=0    :optical ray-tracing, initial conditions as for beam
-    ! igrad=1    :quasi-optical ray-tracing
-    ! igrad=-1   :ray-tracing, init. condit.
-    ! from center of mirror and with angular spread
-    ! ipass=1/2  :1 or 2 passes into plasma
-    ! ipol=0     :compute mode polarization at antenna, ipol=1 use polariz angles
-    read(u, *) params%raytracing%igrad, params%raytracing%ipass, params%raytracing%ipol
-    ! dst        :integration step
-    ! nstep      :maximum number of integration steps
-    ! idst=0/1/2 :0 integration in s, 1 integr. in ct, 2 integr. in Sr
-    read(u, *) params%raytracing%dst, params%raytracing%nstep, params%raytracing%idst
 
-    ! Heating & Current drive
+    read(u, *) params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, params%ecrh_cd%imx
-    ! ========================================================================
-    ! iwarm=0  :no absorption and cd
-    ! iwarm=1  :weakly relativistic absorption
-    ! iwarm=2  :relativistic absorption, n<1 asymptotic expansion
-    ! iwarm=3  :relativistic absorption, numerical integration
-    ! ilarm    :order of larmor expansion
-    ! imx      :max n of iterations in dispersion, imx<0 uses 1st
-    !           iteration in case of failure after |imx| iterations
-    read(u, *) params%ecrh_cd%iwarm,params%ecrh_cd%ilarm,params%ecrh_cd%imx
-    ! ieccd 0/1 NO/YES  ECCD calculation ieccd>0 different CD models
     read(u, *) params%ecrh_cd%ieccd
 
-    ! Wave launcher
-    ! ========================================================================
-    ! alpha0, beta0 (cartesian) launching angles
     read(u, *) params%antenna%alpha, params%antenna%beta
-    ! p0mw      injected power (MW)
     read(u, *) params%antenna%power
-    ! abs(iox)=1/2    OM/XM
-    ! psipol0,chipol0 polarization angles at the antenna (if iox<0)
     read(u, *) params%antenna%iox, params%antenna%psi, params%antenna%chi
     read(u, *) params%antenna%ibeam
     read(u, *) params%antenna%filenm
 
-    ! MHD equilibrium
-    ! ========================================================================
-    ! iequil=0  :vacuum (no plasma)
-    ! iequil=1  :analytical equilibrium
-    ! iequil=2  :read eqdsk
-    ! iequil=2  :read eqdsk, data only valid inside last closed flux surface
     read(u, *) params%equilibrium%iequil
     read(u, *) params%equilibrium%filenm
-    ! icocos    :index for equilibrium from COCOS - O. Sauter Feb 2012
+    read(u, *) params%equilibrium%icocos, params%equilibrium%ipsinorm, &
-    ! ipsinorm  :0 standard EQDSK format, 1 format Portone summer 2004
+               params%equilibrium%idesc, params%equilibrium%ifreefmt
-    read(u, *) params%equilibrium%icocos, params%equilibrium%ipsinorm, params%equilibrium%idesc, params%equilibrium%ifreefmt
+    read(u, *) params%equilibrium%ixp, params%equilibrium%ssplps, &
-    ! ixp=0,-1,+1  : no X point , bottom/up X point
+               params%equilibrium%ssplf
-    ! ssplps       : spline parameter for psi interpolation
+    read(u, *) params%equilibrium%sgnb, params%equilibrium%sgni, &
-    read(u, *) params%equilibrium%ixp, params%equilibrium%ssplps !, params%equilibrium%ssplf
+               params%equilibrium%factb
-    params%equilibrium%ssplf=0.01_wp_
-    ! signum of toroidal B and I
-    ! factb     factor for magnetic field (only for numerical equil)
-    ! scaling adopted: beta=const, qpsi=const, nustar=const
-    read(u, *) params%equilibrium%sgnb, params%equilibrium%sgni, params%equilibrium%factb
 
-    ! Wall
-    ! ========================================================================
     read(u, *) params%misc%rwall
 
-    ! Profiles
+    read(u, *) params%profiles%iprof, params%profiles%irho
-    ! ========================================================================
-    ! iprof=0  :analytical density and temp. profiles
-    ! iprof>0  :numerical density and temp. profiles
-    read(u, *) params%profiles%iprof, params%profiles%irho  ! irho=0,1,2 -> num profiles vs rhot,rhop,psin
     read(u, *) params%profiles%filenm
-    ! psbnd     value of psi ( > 1 ) of density boundary
     read(u, *) params%profiles%psnbnd, params%profiles%sspld
-    ! prfparam%sspld=0.001_wp_
+    read(u, *) params%profiles%factte, params%profiles%factne, &
-    ! iscal       :ne Te scaling 0: nustar=const,  1: n_greenw=const; 2 no rescaling
+               params%profiles%iscal
-    ! factT factn :factor for Te&ne scaling
-    read(u, *) params%profiles%factte, params%profiles%factne, params%profiles%iscal
 
-    ! Output
-    ! ========================================================================
-    ! ipec=0/1  :pec profiles grid in psi/rhop
-    ! nrho      :number of grid steps for pec profiles +1
     read(u, *) params%output%ipec, params%output%nrho
-    ! istpr0    :projection step = dsdt*istprj
-    ! istpl0    :plot step = dsdt*istpl
     read(u, *) params%output%istpr, params%output%istpl
 
     close(u)
-
+  end subroutine read_gray_params
-    ! Convert all filenames to absolute paths
-    block
-      use utils, only : dirname, isrelative
-      if (isrelative(params%antenna%filenm))     &
-        params%antenna%filenm     = dirname(filename) // params%antenna%filenm
-
-      if (isrelative(params%equilibrium%filenm)) &
-        params%equilibrium%filenm = dirname(filename) // params%equilibrium%filenm
-
-      if (isrelative(params%profiles%filenm))    &
-        params%profiles%filenm    = dirname(filename) // params%profiles%filenm
-    end block
-  end subroutine read_parameters
 
 
   subroutine set_globals(params)

src/gray_params.sh

@@ -24,7 +24,7 @@ for set in $sets; do
     cat <<EOF
       case ('$set.$param')
         read (value, *, iostat=error) params%$set%$param
-        if (error > 0) error = 1
+        if (error > 0) error = ERR_VALUE
 EOF
   done
 done

src/ini_parser.f90

@@ -0,0 +1,193 @@
+! This module provides a minimal INI parser
+!
+! The format is:
+!
+!   ; comment
+!   [section-name]
+!   property-name = property-value ; comment
+!
+! The `parse_ini` subroutine takes a file and a handler
+! function that is called with the section, name and value
+! of each property in the INI file.
+!
+module ini_parser
+
+  use logger, only : log_error
+
+  ! INI syntax constants
+  character, parameter :: comment_sign  = ';'
+  character, parameter :: property_sep  = '='
+  character, parameter :: section_start = '['
+  character, parameter :: section_stop  = ']'
+
+  ! Errors
+  enum, bind(C)
+    enumerator :: ini_error   = -1
+    enumerator :: ERR_SUCCESS = 0  ! no errors
+    enumerator :: ERR_SYNTAX  = 1  ! syntax error in the INI file
+    enumerator :: ERR_VALUE   = 2  ! invalid value for a property
+    enumerator :: ERR_UNKNOWN = 3  ! unknown property name
+    enumerator :: ERR_IO      = 4  ! I/O error
+  end enum
+
+  abstract interface
+    function property_handler(section, name, value) result(error)
+      character(*), intent(in) :: section, name, value
+      integer(kind(ini_error)) :: error
+    end function
+  end interface
+
+  private
+  public parse_ini
+  public property_handler
+  public ini_error, ERR_SUCCESS, ERR_SYNTAX, ERR_VALUE, ERR_UNKNOWN, ERR_IO
+
+contains
+
+  subroutine parse_ini(filepath, handler, error)
+    ! Parses a INI file
+    !   filepath: path of the INI file to pase
+    !   handler:  handler function
+    !
+    ! The handler must have the following signature:
+    !
+    !   function handler(section, name, value) result(error)
+    !
+    ! where the error should be:
+    !   ERR_SUCCESS on success;
+    !   ERR_VALUE   on invalid values for this property;
+    !   ERR_UNKNOWN on unknown property.
+    !
+    use utils, only : get_free_unit
+
+    implicit none
+
+    ! function argument
+    character(*),                         intent(in)  :: filepath
+    procedure(property_handler), pointer, intent(in)  :: handler
+    integer(kind(ini_error)),             intent(out) :: error
+
+    ! local variables
+    integer :: ini, sep, n
+    character(256) :: msg
+    character(len=:), allocatable :: line
+    character(len=:), allocatable :: section, name, value
+
+    ! open the INI file
+    ini = get_free_unit()
+    open(unit=ini, file=filepath, iostat=error)
+    if (error /= 0) then
+      write (msg, '("failed to open INI file: ", a)') filepath
+      call log_error(msg, proc='parse_ini', mod='ini_parser')
+      error = ERR_IO
+      return
+    end if
+
+    n = 1 ! line number
+    do
+      ! get one line
+      call getline(ini, line, error)
+      if (error /= 0) exit
+
+      ! skip empty lines
+      if (len(line) == 0) cycle
+
+      ! skip comments
+      if (line(1:1) == comment_sign) cycle
+
+      ! parse section header
+      if (line(1:1) == section_start) then
+        ! split at section stop (ex. [section<here>])
+        sep = str_index(line, section_stop)
+
+        if (sep == 0) then
+          write (msg, '("invalid section header at line ",g0,": ",a)') n, line
+          call log_error(msg, proc='parse_ini', mod='ini_parser')
+          error = ERR_SYNTAX
+          exit
+        end if
+
+        ! update the current section
+        section = line(2:sep - 1)
+        cycle
+      end if
+
+      ! split line at separator (ex. name<here>=value)
+      sep = str_index(line, property_sep)
+      if (sep == 0) then
+        write (msg, '("invalid property definition at line ",g0,": ",a)') n, line
+        call log_error(msg, proc='parse_ini', mod='ini_parser')
+        error = ERR_SYNTAX
+        exit
+      end if
+      name  = trim(line(1:sep - 1))
+      value = trim(line(sep + 1:))
+
+      ! call the handler
+      select case (handler(section, name, value))
+        case (ERR_SUCCESS)
+          cycle
+
+        case (ERR_VALUE)
+          write (msg, '("invalid value for property ",a,": ", a)') name, value
+          call log_error(msg, proc='parse_ini', mod='ini_parser')
+          deallocate(line)
+          error = ERR_VALUE
+          exit
+
+        case (ERR_UNKNOWN)
+          write (msg, '("unknown property ",a)') name
+          call log_error(msg, proc='parse_ini', mod='ini_parser')
+          error = ERR_UNKNOWN
+          exit
+      end select
+    end do
+
+    ! parsed the whole file
+    if (error < 0) error = ERR_SUCCESS
+    close(ini)
+  end subroutine parse_ini
+
+
+  subroutine getline(unit, line, error)
+    ! Reads a line into a deferred length string
+
+    ! subroutine arguments
+    integer,                       intent(in)  :: unit
+    character(len=:), allocatable, intent(out) :: line
+    integer,                       intent(out) :: error
+
+    integer, parameter     :: bufsize = 512
+    character(len=bufsize) :: buffer
+    integer :: chunk
+
+    allocate(character(len=0) :: line)
+    do
+      read(unit, '(a)', advance='no', iostat=error, size=chunk) buffer
+      if (error > 0) exit
+      line = line // buffer(:chunk)
+      if (error < 0) then
+        if (is_iostat_eor(error)) error = 0
+        exit
+      end if
+    end do
+  end subroutine getline
+
+
+  pure function str_index(str, char) result(n)
+    ! Returns the index of the first occurence of `char` within `str`
+    ! If not `char` is not found returns 0
+
+    implicit none
+
+    ! function arguments
+    character(*), intent(in) :: str
+    character,    intent(in) :: char
+
+    ! local variables
+    integer i, n
+
+    n = findloc([(str(i:i) == char, i = 1, len(str))], .true., 1)
+  end function str_index
+
+end module ini_parser

src/main.f90

@@ -2,11 +2,14 @@ program main
   use const_and_precisions, only : wp_, one, zero
   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
   use gray_core,            only : gray_main
   use gray_params,          only : gray_parameters, gray_data, gray_results, &
-                                   read_parameters, params_set_globals => set_globals
+                                   read_gray_params, read_gray_config, &
+                                   params_set_globals => set_globals
+
   implicit none
 
   ! CLI options
@@ -18,6 +21,10 @@ program main
   type(gray_results)    :: results ! Outputs
   integer               :: err     ! Exit code
 
+  ! Store the original working directory
+  character(len=256) :: cwd
+  call getcwd(cwd)
+
   ! Parse the command-line options
   call parse_cli_options(opts)
 
@@ -28,10 +35,28 @@ program main
   ! Activate the given output units
   call set_active_units(opts%units)
 
-  ! Load the parameters from file, apply CLI
+  ! Load the parameters from file and move to its directory
-  ! overrides, and copy them into global
+  ! (all other filepaths are assumed relative to it)
-  ! variables exported by the gray_params
+  if (allocated(opts%config_file)) then
-  call read_parameters(opts%params_file, params)
+    ! Use the newer gray.ini configuration file
+    call log_message(level=INFO, mod='main', msg='using GRAY INI config')
+    call read_gray_config(opts%config_file, params)
+    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')
+    call read_gray_params(opts%params_file, params)
+    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
+
+  ! Apply CLI overrides and copy the parameters into
+  ! global variables exported by the gray_params module
   call parse_param_overrides(params)
   call params_set_globals(params)
 
@@ -42,9 +67,9 @@ program main
   call init_antenna(params%antenna)
   call init_misc(params, data)
 
-  ! Change the current directory to output files here
+  ! Change the current directory to output files there
   if (allocated(opts%output_dir)) then
-    if (chdir(opts%output_dir) /= 0) then
+    if (chdir(trim(cwd) // '/' // opts%output_dir) /= 0) then
       call log_message(level=ERROR, mod='main', &
         msg='chdir to output_dir ('//opts%output_dir//') failed!')
       call exit(1)
@@ -254,9 +279,9 @@ contains
 
       ! Convert psrad to ψ
       select case (params%irho)
-        case (0) ! psrad is ρ_t
+        case (0) ! psrad is ρ_t = √Φ (toroidal flux)
           data%psrad = frhopolv(data%psrad)**2
-        case (1) ! psrad is ρ_p
+        case (1) ! psrad is ρ_p = √ψ (poloidal flux)
           data%psrad = data%psrad**2
         case default ! psrad is already ψ
       end select