diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90
index 2bad8d9..615cc11 100644
--- a/src/coreprofiles.f90
+++ b/src/coreprofiles.f90
@@ -1,5 +1,5 @@
 module coreprofiles
-  use const_and_precisions, only : wp_,zero,one
+  use const_and_precisions, only : wp_, zero, one
   implicit none
 
   integer, save :: npp, nsfd
@@ -185,25 +185,33 @@ contains
   end subroutine read_profiles
 
 
-  subroutine read_profiles_an(filenm,te,ne,zeff,unit)
-    use utils,  only : get_free_unit
-    use logger, only : log_error
+  subroutine read_profiles_an(filenm, data, unit)
+    ! Reads the plasma profiles `data` in the analytical format
+    ! from params%filenm.
+    ! If given, the file is opened in the `unit` number.
+    !
+    ! TODO: add format description
+    use gray_params, only : profiles_data
+    use utils,       only : get_free_unit
+    use logger,      only : log_error
 
     implicit none
-! arguments
-    character(len=*), intent(in) :: filenm
-    real(wp_), dimension(:), allocatable, intent(out) :: te,ne,zeff
-    integer, optional, intent(in) :: unit
-! local variables
+
+    ! subroutine arguments
+    character(len=*),    intent(in)  :: filenm
+    type(profiles_data), intent(out) :: data
+    integer, optional,   intent(in)  :: unit
+
+    ! local variables
     integer :: u
     integer :: err
         
     u = get_free_unit(unit)
 
-    if(allocated(te)) deallocate(te)
-    if(allocated(ne)) deallocate(ne)
-    if(allocated(zeff)) deallocate(zeff)
-    allocate(te(4),ne(3),zeff(1))
+    if (allocated(data%terad)) deallocate(data%terad)
+    if (allocated(data%derad)) deallocate(data%derad)
+    if (allocated(data%zfc))   deallocate(data%zfc)
+    allocate(data%terad(4), data%derad(3), data%zfc(1))
 
     open(file=filenm, status='old', action='read', unit=u, iostat=err)
     if (err /= 0) then
@@ -212,9 +220,9 @@ contains
       call exit(1)
     end if
 
-    read(u,*) ne(1:3) ! dens0,aln1,aln2
-    read(u,*) te(1:4) ! te0,dte0,alt1,alt2
-    read(u,*) zeff(1) ! zeffan
+    read (u,*) data%derad(1:3)  ! dens0, aln1, aln2
+    read (u,*) data%terad(1:4)  ! te0, dte0, alt1, alt2
+    read (u,*) data%zfc(1)      ! zeffan
     close(u)
   end subroutine read_profiles_an
 
@@ -264,7 +272,7 @@ contains
   end subroutine scale_profiles
 
 
-  subroutine set_prfspl(params, data)
+  subroutine set_profiles_spline(params, data)
     ! Computes splines for the plasma profiles data and stores them
     ! in their respective global variables, see the top of this file.
     use simplespline, only : difcs
@@ -361,37 +369,43 @@ contains
         psdbnd=min(psdbnd,xxp)
       end if
       write (msg, '(a,g0.3)') 'density boundary: ψ=', psdbnd
-      call log_info(msg, mod="coreprofiles", proc="set_prfspl")
+      call log_info(msg, mod="coreprofiles", proc="set_profiles_spline")
     end if
 
     deallocate(iwrkf,wrkf,wf)
-  end subroutine set_prfspl
+  end subroutine set_profiles_spline
 
 
-  subroutine unset_prfspl
+  subroutine unset_profiles_spline
     implicit none
 
-    if(allocated(psrad)) deallocate(psrad)
-    if(allocated(ct)) deallocate(ct)
-    if(allocated(cz)) deallocate(cz)
-    if(allocated(tfn)) deallocate(tfn)
-    if(allocated(cfn)) deallocate(cfn)
-  end subroutine unset_prfspl
+    if (allocated(psrad)) deallocate(psrad)
+    if (allocated(ct))    deallocate(ct)
+    if (allocated(cz))    deallocate(cz)
+    if (allocated(tfn))   deallocate(tfn)
+    if (allocated(cfn))   deallocate(cfn)
+  end subroutine unset_profiles_spline
+
+
+  subroutine set_profiles_an(data)
+    ! Stores the analytical profiles data in their respective
+    ! global variables, see the top of this file.
+    use gray_params, only : profiles_data
 
-  subroutine set_prfan(te,ne,zeff)
     implicit none
-    REAL(wp_), dimension(:), intent(in) :: te,ne,zeff
 
-    te0=te(1)
-    dte0=te(2)
-    alt1=te(3)
-    alt2=te(4)
-    dens0=ne(1)
-    aln1=ne(2)
-    aln2=ne(3)
-    zeffan=zeff(1)
+    ! subroutine arguments
+    type(profiles_data), intent(in) :: data
 
-    psdbnd=1.0_wp_
-  end subroutine set_prfan
+    te0    = data%terad(1)
+    dte0   = data%terad(2)
+    alt1   = data%terad(3)
+    alt2   = data%terad(4)
+    dens0  = data%derad(1)
+    aln1   = data%derad(2)
+    aln2   = data%derad(3)
+    zeffan = data%zfc(1)
+    psdbnd = one
+  end subroutine set_profiles_an
   
 end module coreprofiles
diff --git a/src/equilibrium.f90 b/src/equilibrium.f90
index 21fbcb1..3bae324 100644
--- a/src/equilibrium.f90
+++ b/src/equilibrium.f90
@@ -160,20 +160,28 @@ contains
   end subroutine read_eqdsk
 
 
-  subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit)
-    use utils,  only : get_free_unit
-    use logger, only : log_error
+  subroutine read_equil_an(filenm, ipass, data, unit)
+    ! Reads the MHD equilibrium `data` in the analytical format
+    ! from params%filenm.
+    ! If given, the file is opened in the `unit` number.
+    !
+    ! TODO: add format description
+    use gray_params, only : equilibrium_data
+    use utils,       only : get_free_unit
+    use logger,      only : log_error
 
     implicit none
-! arguments
-    character(len=*), intent(in) :: filenm
-    integer, intent(in) :: ipass
-    integer, optional, intent(in) :: unit
-    real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q,rlim,zlim
-! local variables
+
+    ! subroutine arguments
+    character(len=*),       intent(in)  :: filenm
+    integer,                intent(in)  :: ipass
+    type(equilibrium_data), intent(out) :: data
+    integer, optional,      intent(in)  :: unit
+
+    ! local variables
     integer :: i, u, nlim
     integer :: err
-    real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen
+    real(wp_) :: rr0m, zr0m, rpam, b0, q0, qa, alq
 
     u = get_free_unit(unit)
 
@@ -184,48 +192,37 @@ contains
       call exit(1)
     end if
 
-    read(u,*) rr0m,zr0m,rpam
-    read(u,*) b0
-    read(u,*) q0,qa,alq
-!    rcen=rr0m
-!    btrcen=b0
-!    b0=isgnbphi*b0*factb
-!    rvac=rr0m
-!    rax=rr0m
-!    zax=zr0m
-!    zbmin=(zr0-rpa)/1.0e2_wp_
-!    zbmax=(zr0+rpa)/1.0e2_wp_
-    if(allocated(rv)) deallocate(rv)
-    if(allocated(zv)) deallocate(zv)
-    if(allocated(fpol)) deallocate(fpol)
-    if(allocated(q)) deallocate(q)
-    allocate(rv(2),zv(1),fpol(1),q(3))
-    rv(1)=rr0m
-    rv(2)=rpam
-    zv(1)=zr0m
-    fpol(1)=b0*rr0m
-    q(1)=q0
-    q(2)=qa
-    q(3)=alq
+    read(u, *) rr0m, zr0m, rpam
+    read(u, *) b0
+    read(u, *) q0, qa, alq
 
-  if(ipass.ge.2) then
+    if(allocated(data%rv))   deallocate(data%rv)
+    if(allocated(data%zv))   deallocate(data%zv)
+    if(allocated(data%fpol)) deallocate(data%fpol)
+    if(allocated(data%qpsi)) deallocate(data%qpsi)
+    allocate(data%rv(2), data%zv(1), data%fpol(1), data%qpsi(3))
 
-  ! get size of boundary and limiter arrays and allocate them
+    data%rv   = [rr0m, rpam]
+    data%zv   = [zr0m]
+    data%fpol = [b0*rr0m]
+    data%qpsi = [q0, qa, alq]
+
+    if(ipass >= 2) then
+      ! get size of boundary and limiter arrays and allocate them
       read (u,*) nlim
-      if (allocated(rlim)) deallocate(rlim)
-      if (allocated(zlim)) deallocate(zlim)
-  
-  ! store boundary and limiter data
-      if(nlim>0) then
-        allocate(rlim(nlim),zlim(nlim))
-        read(u,*) (rlim(i),zlim(i),i=1,nlim)
-      end if
+      if (allocated(data%rlim)) deallocate(data%rlim)
+      if (allocated(data%zlim)) deallocate(data%zlim)
 
-  end if
+      ! store boundary and limiter data
+      if(nlim > 0) then
+        allocate(data%rlim(nlim), data%zlim(nlim))
+        read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim)
+      end if
+    end if
 
     close(u)
   end subroutine read_equil_an
-  
+
 
   subroutine change_cocos(data, cocosin, cocosout, error)
     ! Convert the MHD equilibrium data from one coordinate convention
@@ -345,7 +342,7 @@ contains
   end subroutine eq_scal
 
 
-  subroutine set_eqspl(params, data)
+  subroutine set_equil_spline(params, data)
     ! Computes splines for the MHD equilibrium data and stores them
     ! in their respective global variables, see the top of this file.
     use const_and_precisions, only : zero, one
@@ -547,7 +544,7 @@ contains
 
     write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
       'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp
-    call log_info(msg, mod='equilibrium', proc='set_eqspl')
+    call log_info(msg, mod='equilibrium', proc='set_equil_spline')
 
     ! search for X-point if params%ixp /= 0
 
@@ -558,7 +555,7 @@ contains
         if(psinxptmp/=-1.0_wp_) then
           write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
             'r', r1, 'z', z1, 'ψ', psinxptmp
-          call log_info(msg, mod='equilibrium', proc='set_eqspl')
+          call log_info(msg, mod='equilibrium', proc='set_equil_spline')
 
           zbinf=z1
           psinop=psinoptmp
@@ -573,7 +570,7 @@ contains
         if(psinxptmp.ne.-1.0_wp_) then
           write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
             'r', r1, 'z', z1, 'ψ', psinxptmp
-          call log_info(msg, mod='equilibrium', proc='set_eqspl')
+          call log_info(msg, mod='equilibrium', proc='set_equil_spline')
 
           zbsup=z1
           psinop=psinoptmp
@@ -601,7 +598,7 @@ contains
       rbinf=r1
       write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
         'r', rbinf, rbsup, 'z', zbinf, zbsup
-      call log_info(msg, mod='equilibrium', proc='set_eqspl')
+      call log_info(msg, mod='equilibrium', proc='set_equil_spline')
     end if
 
     ! Save Bt value on axis (required in flux_average and used in Jcd def)
@@ -613,13 +610,13 @@ contains
     rcen   = data%rvac
 
     write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis
-    call log_info(msg, mod='equilibrium', proc='set_eqspl')
+    call log_info(msg, mod='equilibrium', proc='set_equil_spline')
 
     ! Compute rho_pol/rho_tor mapping based on input q profile
     call setqphi_num(data%psinr,abs(data%qpsi),abs(psia),rhotn)
     call set_rhospl(sqrt(data%psinr),rhotn)
 
-  end subroutine set_eqspl
+  end subroutine set_equil_spline
 
 
   subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
@@ -706,19 +703,33 @@ contains
   end subroutine setqphi_num
 
 
+  subroutine set_equil_an(data, n)
+    ! Computes the analytical equilibrium data and stores them
+    ! in their respective global variables, see the top of this file.
+    use const_and_precisions, only : pi, zero, one
+    use gray_params,          only : equilibrium_data
 
-  subroutine set_equian(rax,zax,a,bax,qax,q1,qexp,n)
-    use const_and_precisions, only : pi,zero,one
     implicit none
-! arguments
-    real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp
-    integer, intent(in), optional :: n
-! local variables
+
+    ! subroutine arguments
+    type(equilibrium_data), intent(in) :: data
+    integer, optional,      intent(in) :: n
+
+    ! local variables
     integer, parameter :: nqdef=101
     integer :: i
-    real(wp_) :: dr,fq0,fq1,qq,res,rn
+    real(wp_) :: dr, fq0, fq1, qq, res, rn
+    real(wp_) :: rax, zax, a, bax, qax, q1, qexp
     real(wp_), dimension(:), allocatable :: rhotn,rhopn
 
+    rax  = data%rv(1)
+    zax  = data%zv(1)
+    a    = data%rv(2)
+    bax  = data%fpol(1) / rax
+    qax  = data%qpsi(1)
+    q1   = data%qpsi(2)
+    qexp = data%qpsi(3)
+
     btaxis=bax
     rmaxis=rax
     zmaxis=zax
@@ -767,8 +778,7 @@ contains
     rhopn=sqrt(psinr)
     
     call set_rhospl(rhopn,rhotn)
-  end subroutine set_equian
-
+  end subroutine set_equil_an
 
 
   subroutine set_rhospl(rhop,rhot)
@@ -1305,36 +1315,37 @@ contains
 
 
 
-  subroutine unset_eqspl
+  subroutine unset_equil_spline
     implicit none
-    if(allocated(tr)) deallocate(tr)
-    if(allocated(tz)) deallocate(tz)
-    if(allocated(tfp)) deallocate(tfp)
-    if(allocated(cfp)) deallocate(cfp)
-    if(allocated(cceq)) deallocate(cceq)
-    if(allocated(cceq01)) deallocate(cceq01)
-    if(allocated(cceq10)) deallocate(cceq10)
-    if(allocated(cceq02)) deallocate(cceq02)
-    if(allocated(cceq20)) deallocate(cceq20)
-    if(allocated(cceq11)) deallocate(cceq11)
-    nsr=0
-    nsz=0
-    nsf=0
-  end subroutine unset_eqspl
+
+    if (allocated(tr))     deallocate(tr)
+    if (allocated(tz))     deallocate(tz)
+    if (allocated(tfp))    deallocate(tfp)
+    if (allocated(cfp))    deallocate(cfp)
+    if (allocated(cceq))   deallocate(cceq)
+    if (allocated(cceq01)) deallocate(cceq01)
+    if (allocated(cceq10)) deallocate(cceq10)
+    if (allocated(cceq02)) deallocate(cceq02)
+    if (allocated(cceq20)) deallocate(cceq20)
+    if (allocated(cceq11)) deallocate(cceq11)
+    nsr = 0
+    nsz = 0
+    nsf = 0
+  end subroutine unset_equil_spline
 
 
 
   subroutine unset_q
     implicit none
     
-    if(allocated(psinr)) deallocate(psinr)
-    if(allocated(cq)) deallocate(cq)
-    nq=0
+    if (allocated(psinr)) deallocate(psinr)
+    if (allocated(cq))    deallocate(cq)
+    nq = 0
   end subroutine unset_q
 
 
 
-  subroutine unset_rhospl
+  subroutine unset_rho_spline
     implicit none
     
     if(allocated(rhopr)) deallocate(rhopr)
@@ -1342,6 +1353,6 @@ contains
     if(allocated(crhop)) deallocate(crhop)
     if(allocated(crhot)) deallocate(crhot)
     nrho=0
-  end subroutine unset_rhospl
+  end subroutine unset_rho_spline
 
 end module equilibrium
diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90
index 999ba6e..8c4b706 100644
--- a/src/gray_jetto1beam.f90
+++ b/src/gray_jetto1beam.f90
@@ -64,7 +64,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
 
   ! Set MHD equilibrium data
   init_equilibrium: block
-    use equilibrium, only : set_eqspl
+    use equilibrium, only : set_equil_spline
 
     ! Copy argument arrays
     data%equilibrium%rv    = r
@@ -81,12 +81,12 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
     data%equilibrium%qpsi  = qpsi
 
     ! Compute splines
-    call set_eqspl(params%equilibrium, data%equilibrium)
+    call set_equil_spline(params%equilibrium, data%equilibrium)
   end block init_equilibrium
 
   ! Set plasma kinetic profiles
   init_profiles: block
-    use coreprofiles, only : set_prfspl
+    use coreprofiles, only : set_profiles_spline
 
     ! Copy argument arrays
     data%profiles%derad = dne
@@ -94,7 +94,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
     data%profiles%zfc   = zeff
 
     ! Compute splines
-    call set_prfspl(params%profiles, data%profiles)
+    call set_profiles_spline(params%profiles, data%profiles)
   end block init_profiles
 
   ! Set wave launcher parameters
@@ -115,16 +115,16 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
 
   ! Free memory
   free_memory: block
-    use equilibrium,  only : unset_eqspl, unset_rhospl, unset_q
-    use coreprofiles, only : unset_prfspl
+    use equilibrium,  only : unset_equil_spline, unset_rho_spline, unset_q
+    use coreprofiles, only : unset_profiles_spline
 
     ! Unset global variables of the `equilibrium` module
-    call unset_eqspl
-    call unset_rhospl
+    call unset_equil_spline
+    call unset_rho_spline
     call unset_q
 
     ! Unset global variables of the `coreprofiles` module
-    call unset_prfspl
+    call unset_profiles_spline
   end block free_memory
 
   ! Copy over the results
diff --git a/src/main.f90 b/src/main.f90
index c1a9102..e46e96d 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -182,7 +182,7 @@ contains
     ! 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
+                            set_equil_an, set_equil_spline, eq_scal
 
     implicit none
 
@@ -190,17 +190,11 @@ contains
     type(gray_parameters), intent(inout) :: params
     type(gray_data),       intent(out)   :: data
 
-    if(params%equilibrium%iequil < 2) then
+    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)
+                         data%equilibrium)
 
       ! Set psia sign to give the correct sign to Iphi
       ! (COCOS=3: psia<0 for Iphi>0)
@@ -216,17 +210,10 @@ contains
     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))
+    if (params%equilibrium%iequil < 2) then
+      call set_equil_an(data%equilibrium)
     else
-      call set_eqspl(params%equilibrium, data%equilibrium)
+      call set_equil_spline(params%equilibrium, data%equilibrium)
     end if
   end subroutine init_equilibrium
 
@@ -234,7 +221,7 @@ contains
   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
+    use equilibrium, only : unset_equil_spline, unset_rho_spline, unset_q
 
     implicit none
 
@@ -248,8 +235,8 @@ contains
     if (allocated(data%rlim)) deallocate(data%rlim, data%zlim)
 
     ! Unset global variables of the `equilibrium` module
-    call unset_eqspl
-    call unset_rhospl
+    call unset_equil_spline
+    call unset_rho_spline
     call unset_q
   end subroutine deinit_equilibrium
 
@@ -261,7 +248,8 @@ contains
     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
+                             scale_profiles, set_profiles_an, &
+                             set_profiles_spline
     implicit none
 
     ! subroutine arguments
@@ -269,10 +257,9 @@ contains
     real(wp_),                 intent(in)  :: factb
     type(profiles_data),       intent(out) :: data
 
-    if(params%iprof == 0) then
+    if (params%iprof == 0) then
       ! Analytical profiles
-      ! TODO: rewrite using derived type
-      call read_profiles_an(params%filenm, data%terad, data%derad, data%zfc)
+      call read_profiles_an(params%filenm, data)
     else
       ! Numerical profiles
       call read_profiles(params%filenm, data)
@@ -291,13 +278,12 @@ contains
     call scale_profiles(params, factb, data)
 
     ! Set global variables
-    if(params%iprof == 0) then
+    if (params%iprof == 0) then
       ! Analytical profiles
-      ! TODO: rewrite using derived type
-      call set_prfan(data%terad, data%derad, data%zfc)
+      call set_profiles_an(data)
     else
       ! Numerical profiles
-      call set_prfspl(params, data)
+      call set_profiles_spline(params, data)
     end if
   end subroutine init_profiles
 
@@ -305,7 +291,7 @@ contains
   subroutine deinit_profiles(data)
     ! Free all memory allocated by the init_profiles subroutine.
     use gray_params,  only : profiles_data
-    use coreprofiles, only : unset_prfspl
+    use coreprofiles, only : unset_profiles_spline
 
     implicit none
 
@@ -317,7 +303,7 @@ contains
     if (allocated(data%terad)) deallocate(data%terad, data%derad, data%zfc)
 
     ! Unset global variables of the `coreprofiles` module
-    call unset_prfspl
+    call unset_profiles_spline
   end subroutine deinit_profiles
 
 
@@ -375,7 +361,7 @@ contains
         allocate(data%equilibrium%rlim(5))
         allocate(data%equilibrium%zlim(5))
         params%raytracing%ipass = abs(params%raytracing%ipass)
-        if(params%equilibrium%iequil < 2) then
+        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))
@@ -406,7 +392,8 @@ contains
                           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 coreprofiles,         only : set_profiles_an, set_profiles_spline, &
+                                     temp, fzeff
     use dispersion,           only : expinit
     use gray_params,          only : gray_parameters, print_parameters
     use beams,                only : launchangles2n, xgygcoeff
@@ -457,7 +444,7 @@ contains
                   p0jk, ext, eyt, iiv)
 
     ! Initialise the dispersion module
-    if(params%ecrh_cd%iwarm > 1) call expinit
+    if (params%ecrh_cd%iwarm > 1) call expinit
 
     ! Initialise the magsurf_data module
     call flux_average ! requires frhotor for dadrhot,dvdrhot