diff --git a/src/main.f90 b/src/main.f90
index 5bfcc1f..8bd982d 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -206,6 +206,7 @@ contains
     ! Reads the MHD equilibrium file (either in the G-EQDSK format
     ! or an analytical description) and initialises the respective
     ! 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, &
                             set_equil_an, set_equil_spline, scale_equil
     use logger, only : log_debug
@@ -215,35 +216,41 @@ contains
     type(gray_data),       intent(out)   :: data
     integer,               intent(out)   :: err
 
-    if (params%equilibrium%iequil < 2) then
-      ! Analytical equilibrium
-      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
-    else
-      ! Numerical equilibrium
-      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 if
+    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
 
     ! Rescale B, I and/or force their signs
     call scale_equil(params%equilibrium, data%equilibrium)
 
-    ! Set global variables (for splines)
-    if (params%equilibrium%iequil < 2) then
-      call set_equil_an
-    else
-      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 if
+    ! Set global variables
+    select case (params%equilibrium%iequil)
+      case (EQ_VACUUM, EQ_ANALYTICAL)
+        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
   end subroutine init_equilibrium
 
 
@@ -270,7 +277,9 @@ contains
     ! 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 gray_params,  only : profiles_parameters, profiles_data, &
+                             RHO_TOR, RHO_POL, RHO_PSI, &
+                             PROF_ANALYTIC, PROF_NUMERIC
     use equilibrium,  only : frhopol
     use coreprofiles, only : read_profiles_an, read_profiles, &
                              scale_profiles, set_profiles_an, &
@@ -287,42 +296,42 @@ contains
     ! local variables
     integer :: i
 
-    if (params%iprof == 0) then
-      ! Analytical profiles
-      call log_debug('loading analytical file',  &
-                     mod='main', proc='init_profiles')
-      call read_profiles_an(params%filenm, data, err)
-      if (err /= 0) return
-    else
-      ! Numerical profiles
-      call log_debug('loading numerical file', &
-                     mod='main', proc='init_profiles')
-      call read_profiles(params%filenm, data, err)
-      if (err /= 0) return
+    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
 
-      ! Convert psrad to ψ
-      select case (params%irho)
-        case (0) ! psrad is ρ_t = √Φ (toroidal flux)
-          data%psrad = [(frhopol(data%psrad(i))**2, i=1,size(data%psrad))]
-        case (1) ! psrad is ρ_p = √ψ (poloidal flux)
-          data%psrad = data%psrad**2
-        case default ! psrad is already ψ
-      end select
-    end if
+      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
 
     ! Rescale input data
     call scale_profiles(params, factb, data)
 
     ! Set global variables
-    if (params%iprof == 0) then
-      ! Analytical profiles
-      call set_profiles_an(params, data)
-    else
-      ! Numerical profiles
-      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 if
+    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
+
   end subroutine init_profiles
 
 
@@ -348,7 +357,7 @@ contains
     ! 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
+    use gray_params, only : antenna_parameters, BEAM_0D, BEAM_1D, BEAM_2D
 
     ! subroutine arguments
     type(antenna_parameters), intent(inout) :: params
@@ -356,16 +365,16 @@ contains
 
     ! Note: α, β are loaded from gray_params.data
     select case (params%ibeam)
-      case (2)
+      case (BEAM_2D)
         ! 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, err=err)
-      case (1)
+      case (BEAM_1D)
         ! 1 degree of freedom
         ! w(z, α), 1/R(z, α)
         call read_beam1(params, err)
-      case default
+      case (BEAM_0D)
         ! fixed w(z), 1/R(z)
         call read_beam0(params, err)
     end select
@@ -374,6 +383,8 @@ contains
 
   subroutine init_misc(params, data)
     ! 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 limiter,              only : limiter_set_globals=>set_globals
     use const_and_precisions, only : cm
@@ -403,16 +414,18 @@ contains
                   params%raytracing%dst * params%raytracing%nstep * cm
 
         ! Max radius, either due to the plasma extent or equilibrium grid
-        if (params%equilibrium%iequil < 2) then
-          ! analytic equilibrium, use R₀+a
-          block
-            use equilibrium, only : model
-            R_max = model%R0 + model%a
-          end block
-        else
-          ! numeric equilibrium, use max R of the grid
-          R_max = data%equilibrium%rv(size(data%equilibrium%rv))
-        end if
+        select case (params%equilibrium%iequil)
+          case (EQ_VACUUM, EQ_ANALYTICAL)
+            ! 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
 
         ! Avoid clipping out the launcher
         R_max = max(R_launch, R_max)