diff --git a/src/beams.f90 b/src/beams.f90
index c046086..c6241ee 100644
--- a/src/beams.f90
+++ b/src/beams.f90
@@ -22,11 +22,7 @@ contains
     integer :: u
     real(wp_) :: ak0,zrcsi,zreta
 
-    if (present(unit)) then
-      u = unit
-    else
-      u = get_free_unit()
-    end if
+    u = get_free_unit(unit)
 
     open(unit=u, file=trim(params%filenm), status='OLD', action='READ')
     read(u, *) params%fghz
@@ -69,11 +65,7 @@ contains
       cbeta, cx0, cy0, cz0, cwaist1, cwaist2,              &
       crci1, crci2, cphi1, cphi2
 
-    if (present(unit)) then
-      u = unit
-    else
-      u = get_free_unit()
-    end if
+    u = get_free_unit(unit)
 
     open(unit=u, file=params%filenm, status='OLD', action='READ')
     read(u,*) params%fghz
@@ -196,11 +188,7 @@ contains
     integer, parameter :: kspl=1
     real(wp_), parameter :: sspl=0.01_wp_
 
-    if (present(unit)) then
-      u = unit
-    else
-      u = get_free_unit()
-    end if
+    u = get_free_unit(unit)
 
     open(unit=u, file=params%filenm, status='OLD', action='READ')
 !=======================================================================================
diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90
index 93bff77..a0a7cc6 100644
--- a/src/coreprofiles.f90
+++ b/src/coreprofiles.f90
@@ -151,11 +151,7 @@ contains
     if(allocated(data%derad)) deallocate(data%derad)
     if(allocated(data%zfc))   deallocate(data%zfc)
 
-    if (present(unit)) then
-      u = unit
-    else
-      u = get_free_unit()
-    end if
+    u = get_free_unit(unit)
 
     ! Read number of rows and allocate the arrays
     open(file=trim(filenm), status='old', action='read', unit=u)
@@ -185,11 +181,7 @@ contains
 ! local variables
     integer :: u
         
-    if (present(unit)) then
-      u=unit
-    else
-      u=get_free_unit()
-    end if
+    u = get_free_unit(unit)
 
     if(allocated(te)) deallocate(te)
     if(allocated(ne)) deallocate(ne)
diff --git a/src/equilibrium.f90 b/src/equilibrium.f90
index 1ee7085..2c2be51 100644
--- a/src/equilibrium.f90
+++ b/src/equilibrium.f90
@@ -55,11 +55,7 @@ contains
     real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis
     real(wp_) :: xdum ! dummy variable, used to discard data
 
-    if(present(unit)) then
-      u = unit
-    else
-      u = get_free_unit()
-    end if
+    u = get_free_unit(unit)
 
     ! Open the G-EQDSK file
     open(file=trim(params%filenm), status='old', action='read', unit=u)
@@ -169,11 +165,8 @@ contains
     integer :: i, u, nlim
     real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen
 
-    if (present(unit)) then
-      u=unit
-    else
-      u=get_free_unit()
-    end if
+    u = get_free_unit(unit)
+
     open(file=trim(filenm),status='old',action='read',unit=u)
     read(u,*) rr0m,zr0m,rpam
     read(u,*) b0
diff --git a/src/gray_params.f90 b/src/gray_params.f90
index 5af7fc8..9319f6a 100644
--- a/src/gray_params.f90
+++ b/src/gray_params.f90
@@ -232,11 +232,7 @@ contains
     ! local variables
     integer :: u, iostat
 
-    if (present(unit)) then
-      u = unit
-    else
-      u = get_free_unit()
-    end if
+    u = get_free_unit(unit)
 
     open(u, file=filename, status='old', action='read', iostat=iostat)
     if (iostat > 0) then
diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90
index 36adb32..1bbf45d 100644
--- a/src/magsurf_data.f90
+++ b/src/magsurf_data.f90
@@ -107,8 +107,7 @@ contains
     use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, &
                             equian,equinum_psi,bfield,frhotor,fq,tor_curr
     use simplespline, only : difcs
-    use dierckx, only : regrid,coeff_parder
-    use utils, only : get_free_unit
+    use dierckx,      only : regrid,coeff_parder
     implicit none
 
     ! local constants
diff --git a/src/utils.f90 b/src/utils.f90
index 477b388..7a93f93 100644
--- a/src/utils.f90
+++ b/src/utils.f90
@@ -246,33 +246,35 @@ contains
     end do
   end subroutine bubble
 
-  function get_free_unit(umin,umax) result(i)
+
+  function get_free_unit(unit) result(i)
+    ! Returns `unit` back or the first free unit
+    ! number `i` if `unit` is absent.
+    ! When no unit is available, returns -1.
+
     implicit none
+
+    ! function arguments
     integer :: i
-    integer, intent(in), optional :: umin, umax
+    integer, intent(in), optional :: unit
+
+    ! local variables
     integer, parameter :: max_allowed = 999
-    integer :: ierr, iend
+    integer :: error
     logical :: ex, op
 
-    if (present(umin)) then
-      i = max(0,umin)    ! start searching from unit min
-    else
-      i = 0
+    if (present(unit)) then
+      i = unit
+      return
     end if
-    if (present(umax)) then
-      iend = min(max(0,umax),max_allowed)
-    else
-      iend = max_allowed
-    end if
-    do
-      if (i>iend) then
-        i=-1  ! no free units found
-        exit
-      end if
-      inquire(unit=i,exist=ex,opened=op,iostat=ierr)
-      if (ierr==0.and.ex.and..not.op) exit  ! unit i exists and is not open
-      i = i + 1
+
+    do i=0,max_allowed
+      inquire(unit=i, exist=ex, opened=op, iostat=error)
+      ! if unit i exists and is free
+      if (error == 0 .and. ex .and. .not. op) return
     end do
+    i = -1
+
   end function get_free_unit