diff --git a/doc/specs/stdlib_io.md b/doc/specs/stdlib_io.md
index 46befe2ea..68e85a027 100644
--- a/doc/specs/stdlib_io.md
+++ b/doc/specs/stdlib_io.md
@@ -17,25 +17,30 @@ Loads a rank-2 `array` from a text file.
 
 ### Syntax
 
-`call ` [[stdlib_io(module):loadtxt(interface)]] `(filename, array [, skiprows] [, max_rows] [, fmt] [, delimiter])`
+`call ` [[stdlib_io(module):loadtxt(interface)]] `(filename, array  [, comments] [, delimiter] [, skiplines] [, max_rows]  [, usecols])`
+
+`call ` [[stdlib_io(module):loadtxt(interface)]] `(unit, array  [, comments] [, delimiter] [, skiplines] [, max_rows]  [, usecols])`
 
 ### Arguments
 
-`filename`: Shall be  a character expression containing the file name from which to load the rank-2 `array`.
+`filename or unit`: Shall be  a character expression containing the file name or an integer containing the unit of an already open file from which to load the rank-2 `array`.
 
 `array`: Shall be an allocatable rank-2 array of type `real`, `complex` or `integer`.
 
-`skiprows` (optional): Skip the first `skiprows` lines. If skipping more rows than present, a 0-sized array will be returned. The default is 0.
+`comments` (optional): Shall  be a character expression of any length used to indicate the start of a    comment. Default: `#`.
 
-`max_rows` (optional): Read `max_rows` lines of content after `skiprows` lines. A negative value results in reading all lines. A value of zero results in no lines to be read. The default value is -1.
+`delimiter` (optional): Shall be a character expression of length 1 that contains the delimiter used to separate the columns. The default is an empty string `''` indicating that any number of whitespace will be considered a delimiter.
 
-`fmt` (optional): Fortran format specifier for the text read.  Defaults to the write format for the data type.  Setting fmt='*' will specify list directed read.   
+`skiplines` (optional): Skip the first `skiplines` lines from file, including comments. If skipping more lines than present, a 0-sized array will be returned. The default is 0.
+
+`max_rows` (optional): Shall be an integer indicating that `max_rows` **rows of data** after `skiprows` will be read. A negative value results in reading all data. The default is to read all lines of data.
+
+`usecols` (optional): Shall be an integer array indicating what columns will be read. For example, ``usecols = (1,3,5)`` will extract the first, third and fifth columns. The default is to read all columns.
 
-`delimiter` (optional): Shall be a character expression of length 1 that contains the delimiter used to separate the columns. The default is `' '`.
 
 ### Return value
 
-Returns an allocated rank-2 `array` with the content of `filename`.
+Returns an allocated rank-2 `array` with the content of the file.
 
 ### Example
 
diff --git a/example/io/example_loadtxt.f90 b/example/io/example_loadtxt.f90
index bd20c93f0..6c4946fb3 100644
--- a/example/io/example_loadtxt.f90
+++ b/example/io/example_loadtxt.f90
@@ -1,12 +1,11 @@
 program example_loadtxt
-  use stdlib_io, only: loadtxt
-  implicit none
-  real, allocatable :: x(:, :)
-  call loadtxt('example.dat', x)
-  
-  ! Can also use list directed format if the default read fails.
-  call loadtxt('example.dat', x, fmt='*')
+    use stdlib_kinds, only: dp
+    use stdlib_io, only: loadtxt
+    implicit none
+    real(dp), allocatable :: x(:, :)
 
-  call loadtxt('example.csv', x, delimiter=',')
+    call loadtxt('example.dat', x)
+
+    call loadtxt('example.csv', x, delimiter=',')
 
 end program example_loadtxt
diff --git a/src/io/stdlib_io.fypp b/src/io/stdlib_io.fypp
index 92d4306c0..77bb5b5c3 100644
--- a/src/io/stdlib_io.fypp
+++ b/src/io/stdlib_io.fypp
@@ -13,8 +13,11 @@ module stdlib_io
                                  FMT_COMPLEX_SP, FMT_COMPLEX_DP, FMT_COMPLEX_XDP, FMT_COMPLEX_QP
   use stdlib_error, only: error_stop, state_type, STDLIB_IO_ERROR
   use stdlib_optval, only: optval
-  use stdlib_ascii, only: is_blank
+  use stdlib_ascii, only: is_blank, CR, LF, TAB
   use stdlib_string_type, only : string_type, assignment(=), move
+  use stdlib_strings, only: starts_with
+  use stdlib_str2num, only: to_num_from_stream
+
   implicit none
   private
   ! Public API
@@ -46,7 +49,10 @@ module stdlib_io
   public :: parse_mode
 
   !> Default delimiter for loadtxt, savetxt and number_of_columns
-  character(len=1), parameter :: delimiter_default = " "
+  character(len=1), parameter :: delimiter_default = ""
+  character(len=1), parameter :: comment_default = "#"
+  character(len=2), parameter :: nl = CR//LF
+  character(len=*), parameter :: blanks = " "//TAB
 
   public :: FMT_INT, FMT_REAL_SP, FMT_REAL_DP, FMT_REAL_XDP, FMT_REAL_QP
   public :: FMT_COMPLEX_SP, FMT_COMPLEX_DP, FMT_COMPLEX_XDP, FMT_COMPLEX_QP
@@ -66,8 +72,10 @@ module stdlib_io
     !!
     !! Loads a 2D array from a text file
     !! ([Specification](../page/specs/stdlib_io.html#description))
+    #:for arg1 in ['filename', 'unit']
     #:for k1, t1 in KINDS_TYPES
-      module procedure loadtxt_${t1[0]}$${k1}$
+      module procedure loadtxt_${t1[0]}$${k1}$${arg1[0]}$
+    #:endfor
     #:endfor
   end interface loadtxt
 
@@ -83,153 +91,319 @@ module stdlib_io
 
 contains
 
+  #:for arg1 in ['filename', 'unit']
   #:for k1, t1 in KINDS_TYPES
-    subroutine  loadtxt_${t1[0]}$${k1}$(filename, d, skiprows, max_rows, fmt, delimiter)
-      !! version: experimental
-      !!
-      !! Loads a 2D array from a text file.
-      !!
-      !! Arguments
-      !! ---------
-      !!
-      !! Filename to load the array from
-      character(len=*), intent(in) :: filename
-      !! The array 'd' will be automatically allocated with the correct dimensions
-      ${t1}$, allocatable, intent(out) :: d(:,:)
-      !! Skip the first `skiprows` lines. If skipping more rows than present, a 0-sized array will be returned. The default is 0.
-      integer, intent(in), optional :: skiprows
-      !! Read `max_rows` lines of content after `skiprows` lines.
-      !! A negative value results in reading all lines.
-      !! A value of zero results in no lines to be read.
-      !! The default value is -1.
-      integer, intent(in), optional :: max_rows
-      character(len=*), intent(in), optional :: fmt
-      character(len=1), intent(in), optional :: delimiter
-      character(len=:), allocatable :: fmt_
-      character(len=1) :: delimiter_
-      !!
-      !! Example
-      !! -------
-      !!
-      !!```fortran
-      !! ${t1}$, allocatable :: data(:, :)
-      !! call loadtxt("log.txt", data)  ! 'data' will be automatically allocated
-      !!```
-      !!
-      !! Where 'log.txt' contains for example::
-      !!
-      !!     1 2 3
-      !!     2 4 6
-      !!     8 9 10
-      !!     11 12 13
-      !!     ...
-      !!
-      integer :: s
-      integer :: nrow, ncol, i, j, ios, skiprows_, max_rows_, istart, iend
-      character(len=:), allocatable :: line, iomsg_
-      character(len=1024) :: iomsg, msgout
-
-      skiprows_ = max(optval(skiprows, 0), 0)
-      max_rows_ = optval(max_rows, -1)
-      delimiter_ = optval(delimiter, delimiter_default)
+    subroutine loadtxt_${t1[0]}$${k1}$${arg1[0]}$ (${arg1}$, d, comments, delimiter, skiplines, max_rows, usecols)
+        !! version: experimental
+        !!
+        !! Loads a 2D array from a text file.
+        !!
+        !! Arguments
+        !! ---------
+        !!
+      #:if 'filename' in arg1
+        !! Filename with the array to load
+        character(len=*), intent(in) :: filename  ! File to save the array to
+      #:elif 'unit' in arg1
+        !! unit of an open file from where to load the array
+        integer, intent(in) :: unit
+      #:endif
+        !! The array 'd' will be automatically allocated with the correct dimensions
+        ${t1}$, allocatable, intent(out) :: d(:, :) 
+        integer, intent(in), optional :: skiplines !! Skip the first `skiplines` lines. If skipping more rows than present, a 0-sized array will be returned. The default is 0.
+        integer, intent(in), optional :: max_rows !! Read `max_rows` lines of content after `skiplines` lines. A negative value results in reading all lines. A value of zero results in no lines to be read. The default is to read all rows.
+        character(len=*), intent(in), optional :: comments !! from comments symbol until line end everything else will be ignored. The default is '#'.
+        character(len=*), intent(in), optional :: delimiter !! Character used to separate values in a line. The default is any number of spaces and or tabs.
+        integer, intent(in), optional :: usecols(:) !! Array of column indices to read. If not provided, all columns are read.
+        !!
+        !! Example
+        !! -------
+        !!
+        !!```fortran
+        !! ${t1}$, allocatable :: data(:, :)
+        !! call loadtxt("log.txt", data)  ! 'data' will be automatically allocated
+        !!```
+        !!
+        !! Where 'log.txt' contains for example::
+        !!
+        !!     1 2 3
+        !!     2 4 6
+        !!     8 9 10
+        !!     11 12 13
+        !!     ...
+        !!
+        integer :: iostat
+        integer :: skiplines_, max_rows_
+        integer :: fsze, nrows, nrows_effective, ncols, j, start_effective
+        character(:), allocatable, target :: ff
+        character(len=:), pointer :: ffp
+        integer :: line_start, line_end
+        character(len=:), allocatable :: delim_
+        character(len=:), allocatable :: comment_
+        integer, allocatable ::  usecols_(:)
+        character(len=1024) :: iomsg, msgout, fout
+        character(len=16) :: readable
+        logical :: opened
+        !
+        integer :: row, row_effective, err
+        ${t1}$, allocatable :: cols(:)
+      #:if 'complex' in t1
+        real(${k1}$) :: reval, imval
+      #:endif
+      #:if 'filename' in arg1
+        integer :: unit
+      #:endif
 
-      s = open(filename)
+        comment_ = optval(comments, comment_default)
+        delim_ = optval(delimiter, delimiter_default)
+        skiplines_ = optval(skiplines, 0)
+        ! max_rows will be set later, after determining number of rows
+
+        !----------------------------------------- Check file
+      #:if 'filename' in arg1
+        unit = open (filename, "rb", iostat=iostat)
+        if (iostat /= 0) then
+            write(msgout,'(a)') "loadtxt: error opening file "//trim(filename)
+            call error_stop(msg=trim(msgout))
+        end if
+        fout = filename           ! fout is used for unified error message later
+      #:else 
+        ! first argument is unit
+        inquire (unit=unit, opened=opened, action=readable) 
+        if((.not. opened) .or. (readable(1:1) /= 'R')) then
+            write (msgout,'(a,i0,a)') 'loadtxt error: unit ',unit,' not open for reading'
+            call error_stop(msg=trim(msgout))
+        end if
+        write(fout,'(i0)') unit
+        fout = adjustl(fout)  ! fout is used for unified error message later
+      #:endif
+            
+        !----------------------------------------- Load file in a single string
+        inquire (unit=unit, size=fsze)
+        if(fsze == 0) then
+            write (msgout,'(a,i0,a)') 'loadtxt error: file empty'
+            call error_stop(msg=trim(msgout))
+        else
+            allocate (character(fsze) :: ff)
+            read (unit, iostat=iostat, iomsg=iomsg) ff
+            if (iostat /= 0) then
+                write(msgout,'(a)') "loadtxt: error reading file "//trim(fout)//"("//trim(iomsg)//")"
+                call error_stop(msg=trim(msgout))
+            end if
+        end if
+      #:if 'filename' in arg1
+        close (unit)
+      #:endif
 
-      ! determine number or rows
-      nrow = number_of_rows(s)
-      skiprows_ = min(skiprows_, nrow)
-      if ( max_rows_ < 0 .or. max_rows_ > (nrow - skiprows_) ) max_rows_ = nrow - skiprows_
+        ffp => ff
+        start_effective = 1     ! Start after skiplines
+
+        !----------------------------------------- Count lines and columns
+        nrows = 0               ! Total number of rows (including empty and commented lines)
+        nrows_effective = 0     ! rows with data
+        ncols = 0
+        do while (len(ffp) > 0)
+            line_end = shift_to_eol(ffp)
+            if (line_end > len(ffp)) exit  ! No more lines
+            line_start = shift_to_nonwhitespace(ffp(:line_end)) ! Skip initial blanks in line
+            nrows = nrows + 1
+            if (nrows <= skiplines_) then
+                start_effective = start_effective + line_end ! Remember position to use aS starting point when reading
+                ffp => ffp(line_end + 1:) ! Skip the line
+                cycle
+            end if
 
-      ! determine number of columns
-      ncol = 0
-      if ( skiprows_ < nrow ) ncol = number_of_columns(s, skiprows=skiprows_, delimiter=delimiter_)
+            if (starts_with(ffp(line_start:), comment_) .or. (line_start == line_end)) then
+                ffp => ffp(line_end + 1:) ! Skip comment lines and blank lines
+                cycle
+            end if
+            nrows_effective = nrows_effective + 1
+            !
+            ! if ncols is not set yet, determine the number of columns
+            if (ncols == 0) ncols = number_cols_line(ffp(line_start:line_end), delim_, comment_)
+            ffp => ffp(line_end + 1:) ! go to next line
+        end do
+
+        !----------------------------------------- Allocate and read data
       #:if 'complex' in t1
-      ncol = ncol / 2
+        ncols = ncols / 2
       #:endif
 
-      allocate(d(max_rows_, ncol))
-      if (max_rows_ == 0 .or. ncol == 0) return
+        max_rows_ = min(optval(max_rows, nrows_effective), nrows_effective)
+        ! If there is no data we will return an empty array
+        if ((max_rows_ <= 0) .or. (ncols == 0)) then
+            allocate (d(0, 0))
+            return
+        end if
 
-      do i = 1, skiprows_
-        read(s, *, iostat=ios, iomsg=iomsg)
-        
-        if (ios/=0) then 
-           write(msgout,1) trim(iomsg),i,trim(filename) 
-           1 format('loadtxt: error <',a,'> skipping line ',i0,' of ',a,'.')
-           call error_stop(msg=trim(msgout))
+        if (present(usecols)) then ! user set columns to extract
+            usecols_ = usecols
+        else                    ! extract all columns
+            usecols_ = [(j, j=1, ncols)]
         end if
-        
-      end do
-      
-      ! Default to format used for savetxt if fmt not specified.
-      #:if 'real' in t1
-      fmt_ = optval(fmt, "(*"//FMT_REAL_${k1}$(1:len(FMT_REAL_${k1}$)-1)//",:,1x))")
-      #:elif 'complex' in t1
-      fmt_ = optval(fmt, "(*"//FMT_COMPLEX_${k1}$(1:len(FMT_COMPLEX_${k1}$)-1)//",:,1x))")
-      #:else
-      fmt_ = optval(fmt, "*")
-      #:endif      
-
-      if ( fmt_ == '*' ) then
-        ! Use list directed read if user has specified fmt='*'
-        if (is_blank(delimiter_) .or. delimiter_ == ",") then
-          do i = 1, max_rows_
-            read (s,*,iostat=ios,iomsg=iomsg) d(i, :)
-            
-            if (ios/=0) then 
-              write(msgout,2) trim(iomsg),size(d,2),i,trim(filename)
-              call error_stop(msg=trim(msgout))
-            end if          
-            
-          enddo
-        ! Otherwise read each value separately
-        else
-          do i = 1, max_rows_
-            call get_line(s, line, ios, iomsg_)
-            if (ios/=0) then 
-               write(msgout,2) trim(iomsg_),size(d,2),i,trim(filename)
-               call error_stop(msg=trim(msgout))
+        allocate (d(max_rows_, size(usecols_)))
+        allocate (cols(ncols))  ! Used to hold each row
+
+        row_effective = 0
+        ffp => ff(start_effective:) ! Reset pointer to the beginning of the file after skiplines
+        nrows = nrows - skiplines_
+
+        do row = 1, nrows
+            line_end = shift_to_eol(ffp)
+            line_start = shift_to_nonwhitespace(ffp(:line_end)) ! Avoid initial blanks in line
+
+            if (starts_with(ffp(line_start:), comment_) .or. &
+                (line_start == line_end)) then
+                ffp => ffp(line_end + 1:) ! Skip comment lines and blank lines
+                cycle
             end if
-  
-            istart = 0
-            do j = 1, ncol - 1
-              iend = index(line(istart+1:), delimiter_)
-              read (line(istart+1:istart+iend-1),*,iostat=ios,iomsg=iomsg) d(i, j)
-              if (ios/=0) then 
-                 write(msgout,2) trim(iomsg),size(d,2),i,trim(filename)
-                 call error_stop(msg=trim(msgout))
-              end if
-              istart = istart + iend
+
+            row_effective = row_effective + 1
+            do j = 1, ncols     ! Read a row
+              #:if 'complex' in t1
+                reval = to_num_from_stream(ffp, reval)
+                if (in_delim(ffp, delim_)) then
+                    ffp => ffp(shift_to_nondelim(ffp, delim_):)
+                end if
+                imval = to_num_from_stream(ffp, imval)
+                if (in_delim(ffp, delim_)) then
+                    ffp => ffp(shift_to_nondelim(ffp, delim_):)
+                end if
+                cols(j) = cmplx(reval, imval, kind(cols(j)))
+              #:else  
+                cols(j) = to_num_from_stream(ffp, cols(j))
+                if (in_delim(ffp, delim_)) then
+                    ffp => ffp(shift_to_nondelim(ffp, delim_):)
+                end if
+              #:endif
+
+                if (scan(ffp(1:1), nl) /= 0) then ! If EOL => no more cols
+                    exit
+                end if
+            end do
+            ! Copy the columns of the current row to d(row_effective,:)
+            do j = 1, size(usecols_)
+                d(row_effective, j) = cols(usecols_(j))
             end do
-  
-            read (line(istart+1:),*,iostat=ios,iomsg=iomsg) d(i, ncol)
-            if (ios/=0) then 
-               write(msgout,2) trim(iomsg),size(d,2),i,trim(filename)
-               call error_stop(msg=trim(msgout))
-            end if          
-            
-          enddo
-        end if
-      else
-        ! Otherwise pass default or user specified fmt string.
-        do i = 1, max_rows_
-          read (s,fmt_,iostat=ios,iomsg=iomsg) d(i, :)
-          
-          if (ios/=0) then 
-             write(msgout,2) trim(iomsg),size(d,2),i,trim(filename)
-             call error_stop(msg=trim(msgout))
-          end if             
-          
-        enddo
-      endif
 
-      close(s)
-      
-      2 format('loadtxt: error <',a,'> reading ',i0,' values from line ',i0,' of ',a,'.')
+            if (row_effective >= max_rows_) return
+            line_end = shift_to_eol(ffp)
+            ffp => ffp(line_end + 1:)
+        end do
 
-    end subroutine loadtxt_${t1[0]}$${k1}$
+    end subroutine loadtxt_${t1[0]}$${k1}$${arg1[0]}$
   #:endfor
+  #:endfor
+
+    function number_cols_line(row, delimiter, comment) result(ncols)
+        implicit none
+        character(len=*), intent(in) :: row !<
+        character(len=*), intent(in) :: delimiter !<
+        character(len=*), intent(in) :: comment !<
+        character(len=:), allocatable :: line
+        integer :: ncols
+        integer :: pos, p
+        ncols = 0
+        pos = index(row, comment)
+        if (pos == 0) pos = len(row)
+        line = trim(adjustl(row(:pos))) ! Line with no blanks around and no comments
+        pos = 1
+        do
+            p = shift_after_next_delim(line(pos:), delimiter)
+            pos = pos + p - 1 ! Find delimiter
+            if (pos >= len(line)) exit
+            ncols = ncols + 1
+        end do
+        ncols = ncols + 1
+    end function number_cols_line
+
+    elemental function in_delim(s, delim) result(m)
+        !! Check if current position is the init of a delimiter
+        character(len=*), intent(in) :: s !! character chain
+        character(len=*), intent(in) :: delim !! character chain
+        logical :: m !! True or False
+        !----------------------------------------------
+        if (delim == delimiter_default) then
+            m = (scan(s(1:1), blanks) /= 0)
+        else
+            m = starts_with(s(shift_to_nonwhitespace(s):), delim)
+        end if
+    end function in_delim
+
+    elemental function shift_to_eol(s) result(p)
+        !! move string to position of the next end-of-line character
+        character(len=*), intent(in) :: s !! character chain
+        integer :: p !! position
+        !----------------------------------------------
+        p = scan(s, nl)
+        if (p < len(s)) then ! If CRLF, move to LF
+            if (s(p:p + 1) == nl) p = p + 1
+        end if
+
+    end function shift_to_eol
+
+    function shift_after_next_delim(s, delim) result(p)
+      !! move string to position of the next non delimiter character
+        character(len=*), intent(in) :: s !! character chain
+        character(len=*), intent(in) :: delim !! character chain
+        integer :: p !! position
+        !----------------------------------------------
+        if (delim == delimiter_default) then
+            p = 1
+            if (.not. is_blank(s(p:p))) p = shift_to_whitespace(s)
+            p = p + shift_to_nonwhitespace(s(p:)) - 1
+        else
+            p = index(s, delim)
+            if (p == 0) then
+                p = len(s)
+            else
+                p = p + len(delim)
+            end if
+        end if
+        if (p > len(s)) p = len(s)
+    end function shift_after_next_delim
+
+    elemental function shift_to_nondelim(s, delim) result(p)
+      !! move string to position of the next non delimiter character
+      !! Assumes that it is in a delim
+        character(len=*), intent(in) :: s !! character chain
+        character(len=*), intent(in) :: delim !! character chain
+        integer :: p !! position
+        !----------------------------------------------
+        if (delim == delimiter_default) then
+            p = shift_to_nonwhitespace(s)
+        else
+            ! Check first if we are at the beginning of a delimiter
+            p = index(s, delim)
+            if (p == 0) then
+                p = len(s)
+            else
+                p = p + len(delim)
+                p = p + shift_to_nonwhitespace(s(p:)) - 1 ! Extra-spaces make to_num_from_stream fail
+                if (p > len(s)) p = len(s)
+            end if
+        end if
 
+    end function shift_to_nondelim
+
+    elemental function shift_to_nonwhitespace(s) result(p)
+    !! move string to position of the next non white space character
+        character(len=*), intent(in) :: s !! character chain
+        integer :: p !! position
+        !----------------------------------------------
+        ! p = verify(s, blanks//nl)
+        p = verify(s, blanks)
+        if (p == 0) p = len(s)
+    end function shift_to_nonwhitespace
+
+    elemental function shift_to_whitespace(s) result(p)
+    !! move string to position of the next white space character
+        character(len=*), intent(in) :: s !! character chain
+        integer :: p !! position
+        !----------------------------------------------
+        p = scan(s, blanks)
+        if (p == 0) p = len(s)
+    end function shift_to_whitespace
 
   #:for k1, t1 in KINDS_TYPES
     subroutine savetxt_${t1[0]}$${k1}$(filename, d, delimiter)
@@ -291,75 +465,6 @@ contains
   #:endfor
 
 
-  integer function number_of_columns(s, skiprows, delimiter)
-    !! version: experimental
-    !!
-    !! determine number of columns
-    integer,intent(in) :: s
-    integer, intent(in), optional :: skiprows
-    character(len=1), intent(in), optional :: delimiter
-
-    integer :: ios, skiprows_, i
-    character :: c
-    character(len=:), allocatable :: line
-    character(len=1) :: delimiter_
-    logical :: last_delim
-
-    skiprows_ = optval(skiprows, 0)
-    delimiter_ = optval(delimiter, delimiter_default)
-
-    rewind(s)
-
-    do i = 1, skiprows_
-      read(s, *)
-    end do
-    number_of_columns = 0
-    
-    ! Read first non-skipped line as a whole
-    call get_line(s, line, ios)
-    if (ios/=0 .or. .not.allocated(line)) return
-
-    last_delim = .true.
-    if (delimiter_ == delimiter_default) then
-      do i = 1,len(line)
-        c = line(i:i)
-        if (last_delim .and. .not. is_blank(c)) number_of_columns = number_of_columns + 1
-        last_delim = is_blank(c)
-      end do
-    else
-      do i = 1,len(line)
-        if (line(i:i) == delimiter_) number_of_columns = number_of_columns + 1
-      end do
-      if (number_of_columns == 0) then
-        if (len_trim(line) /= 0) number_of_columns = 1
-      else
-        number_of_columns = number_of_columns + 1
-      end if
-    end if
-    rewind(s)
-
-  end function number_of_columns
-
-
-  integer function number_of_rows(s) result(nrows)
-    !! version: experimental
-    !!
-    !! Determine the number or rows in a file
-    integer, intent(in)::s
-    integer :: ios
-
-    rewind(s)
-    nrows = 0
-    do
-      read(s, *, iostat=ios)
-      if (ios /= 0) exit
-      nrows = nrows + 1
-    end do
-
-    rewind(s)
-
-  end function number_of_rows
-
 
   integer function open(filename, mode, iostat) result(u)
     !! version: experimental
diff --git a/test/io/test_loadtxt.f90 b/test/io/test_loadtxt.f90
index 88277fcf3..3c692dd41 100644
--- a/test/io/test_loadtxt.f90
+++ b/test/io/test_loadtxt.f90
@@ -14,273 +14,282 @@ subroutine collect_loadtxt(testsuite)
         type(unittest_type), allocatable, intent(out) :: testsuite(:)
 
         testsuite = [ &
-            new_unittest("loadtxt_int32", test_loadtxt_int32), &
-            new_unittest("loadtxt_sp", test_loadtxt_sp), &
-            new_unittest("loadtxt_sp_huge", test_loadtxt_sp_huge), &
-            new_unittest("loadtxt_sp_tiny", test_loadtxt_sp_tiny), &
-            new_unittest("loadtxt_dp", test_loadtxt_dp), &
-            new_unittest("loadtxt_dp_max_skip", test_loadtxt_dp_max_skip), &
-            new_unittest("loadtxt_dp_huge", test_loadtxt_dp_huge), &
-            new_unittest("loadtxt_dp_tiny", test_loadtxt_dp_tiny), &
-            new_unittest("loadtxt_complex", test_loadtxt_complex) &
-        ]
+                    new_unittest("loadtxt_int32", test_loadtxt_int32), &
+                    new_unittest("loadtxt_sp", test_loadtxt_sp), &
+                    new_unittest("loadtxt_sp_huge", test_loadtxt_sp_huge), &
+                    new_unittest("loadtxt_sp_tiny", test_loadtxt_sp_tiny), &
+                    new_unittest("loadtxt_dp", test_loadtxt_dp), &
+                    new_unittest("loadtxt_dp_max_skip", test_loadtxt_dp_max_skip), &
+                    new_unittest("loadtxt_dp_huge", test_loadtxt_dp_huge), &
+                    new_unittest("loadtxt_dp_tiny", test_loadtxt_dp_tiny), &
+                    new_unittest("loadtxt_complex", test_loadtxt_complex), &
+                    new_unittest("loadtxt_dp_comm_skip_max", test_loadtxt_dp_comm_skip_max) &
+                    ]
 
     end subroutine collect_loadtxt
 
-
     subroutine test_loadtxt_int32(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        integer(int32), allocatable :: input(:,:), expected(:,:)
-        real(sp), allocatable :: harvest(:,:)
+        integer(int32), allocatable :: input(:, :), expected(:, :)
+        real(sp), allocatable :: harvest(:, :)
         integer :: n
-        allocate(harvest(10,10))
-        allocate(input(10,10))
-        allocate(expected(10,10))
+        allocate (harvest(10, 10))
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
         do n = 1, 10
             call random_number(harvest)
             input = int(harvest * 100)
             call savetxt('test_int32.txt', input)
             call loadtxt('test_int32.txt', expected)
-            call check(error, all(input == expected),'Default list directed read failed')
-            if (allocated(error)) return
-            call loadtxt('test_int32.txt', expected, fmt='*')
-            call check(error, all(input == expected),'User specified list directed read faile')
+            call check(error, all(input == expected), 'Default list directed read failed')
             if (allocated(error)) return
+            !
             call savetxt('test_int32.txt', input, delimiter=',')
             call loadtxt('test_int32.txt', expected, delimiter=',')
-            call check(error, all(input == expected),'User specified delimiter `,` read failed')
+            call check(error, all(input == expected), 'User specified delimiter `,` read failed')
             if (allocated(error)) return
+            !
             call savetxt('test_int32.txt', input, delimiter='-')
             call loadtxt('test_int32.txt', expected, delimiter='-')
-            call check(error, all(input == expected),'User specified delimiter `-` read failed')
+            call check(error, all(input == expected), 'User specified delimiter `-` read failed')
             if (allocated(error)) return
         end do
 
     end subroutine test_loadtxt_int32
 
-
     subroutine test_loadtxt_sp(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        real(sp), allocatable :: input(:,:), expected(:,:)
-        character(len=*), parameter :: FMT_REAL_SP = '(es15.8e2)'
+        real(sp), allocatable :: input(:, :), expected(:, :)
         integer :: n
 
-        allocate(input(10,10))
-        allocate(expected(10,10))
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
         do n = 1, 10
             call random_number(input)
             input = input - 0.5
             call savetxt('test_sp.txt', input)
             call loadtxt('test_sp.txt', expected)
-            call check(error, all(input == expected),'Default format read failed')
-            if (allocated(error)) return
-            call loadtxt('test_sp.txt', expected, fmt='*')
-            call check(error, all(input == expected),'List directed read failed')
-            if (allocated(error)) return
-            call loadtxt('test_sp.txt', expected, fmt="(*"//FMT_REAL_sp(1:len(FMT_REAL_sp)-1)//",1x))")
-            call check(error, all(input == expected),'User specified format failed')
+            call check(error, all(input == expected), 'Default format read failed')
             if (allocated(error)) return
+            !
             call savetxt('test_sp.txt', input, delimiter=',')
             call loadtxt('test_sp.txt', expected, delimiter=',')
-            call check(error, all(input == expected),'User specified delimiter `,` read failed')
+            call check(error, all(input == expected), 'User specified delimiter `,` read failed')
             if (allocated(error)) return
+            !
             call savetxt('test_sp.txt', input, delimiter=';')
             call loadtxt('test_sp.txt', expected, delimiter=';')
-            call check(error, all(input == expected),'User specified delimiter `;` read failed')
+            call check(error, all(input == expected), 'User specified delimiter `;` read failed')
             if (allocated(error)) return
         end do
 
     end subroutine test_loadtxt_sp
 
-
     subroutine test_loadtxt_sp_huge(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        real(sp), allocatable :: input(:,:), expected(:,:)
+        real(sp), allocatable :: input(:, :), expected(:, :)
         integer :: n
-        character(len=*), parameter :: FMT_REAL_SP = '(es15.8e2)'
-        
-        allocate(input(10,10))
-        allocate(expected(10,10))
+
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
 
         do n = 1, 10
             call random_number(input)
             input = (input - 0.5) * huge(input)
             call savetxt('test_sp_huge.txt', input)
             call loadtxt('test_sp_huge.txt', expected)
-            call check(error, all(input == expected),'Default format read failed')
-            if (allocated(error)) return
-            call loadtxt('test_sp_huge.txt', expected, fmt='*')
-            call check(error, all(input == expected),'List directed read failed')
-            if (allocated(error)) return
-            call loadtxt('test_sp_huge.txt', expected, fmt="(*"//FMT_REAL_sp(1:len(FMT_REAL_sp)-1)//",1x))")
-            call check(error, all(input == expected),'User specified format failed')
+            call check(error, all(input == expected), 'Default format read failed')
             if (allocated(error)) return
         end do
 
     end subroutine test_loadtxt_sp_huge
 
-
     subroutine test_loadtxt_sp_tiny(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        real(sp), allocatable :: input(:,:), expected(:,:)
+        real(sp), allocatable :: input(:, :), expected(:, :)
         integer :: n
-        character(len=*), parameter :: FMT_REAL_SP = '(es15.8e2)'
 
-        allocate(input(10,10))
-        allocate(expected(10,10))
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
 
         do n = 1, 10
             call random_number(input)
             input = (input - 0.5) * tiny(input)
             call savetxt('test_sp_tiny.txt', input)
             call loadtxt('test_sp_tiny.txt', expected)
-            call check(error, all(input == expected),'Default format read failed')
-            if (allocated(error)) return
-            call loadtxt('test_sp_tiny.txt', expected, fmt='*')
-            call check(error, all(input == expected),'List directed read failed')
-            if (allocated(error)) return
-            call loadtxt('test_sp_tiny.txt', expected, fmt="(*"//FMT_REAL_sp(1:len(FMT_REAL_sp)-1)//",1x))")
-            call check(error, all(input == expected),'User specified format failed')
+            call check(error, all(input == expected), 'Default format read failed')
             if (allocated(error)) return
         end do
 
     end subroutine test_loadtxt_sp_tiny
 
-
     subroutine test_loadtxt_dp(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        real(dp), allocatable :: input(:,:), expected(:,:)
+        real(dp), allocatable :: input(:, :), expected(:, :)
         integer :: n
-        character(len=*), parameter :: FMT_REAL_DP = '(es24.16e3)'
-
-        allocate(input(10,10))
-        allocate(expected(10,10))
-
+        real(dp) :: eps
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
+        eps = epsilon(1._dp)
         do n = 1, 10
             call random_number(input)
             input = input - 0.5
             call savetxt('test_dp.txt', input)
             call loadtxt('test_dp.txt', expected)
-            call check(error, all(input == expected),'Default format read failed')
-            if (allocated(error)) return
-            call loadtxt('test_dp.txt', expected, fmt='*')
-            call check(error, all(input == expected),'List directed read failed')
-            if (allocated(error)) return
-            call loadtxt('test_dp.txt', expected, fmt="(*"//FMT_REAL_dp(1:len(FMT_REAL_dp)-1)//",1x))")
-            call check(error, all(input == expected),'User specified format failed')
+            call check(error, all(abs(input - expected) < eps), 'Default format read failed')
             if (allocated(error)) return
             call savetxt('test_dp.txt', input, delimiter=',')
             call loadtxt('test_dp.txt', expected, delimiter=',')
-            call check(error, all(input == expected),'User specified delimiter read failed')
+            call check(error, all(abs(input - expected) < eps), 'User specified delimiter read failed')
             if (allocated(error)) return
         end do
 
     end subroutine test_loadtxt_dp
 
+    subroutine test_loadtxt_dp_comm_skip_max(error)
+        !> Error handling
+        type(error_type), allocatable, intent(out) :: error
+        real(dp), allocatable :: expected(:, :)
+        integer :: unit
+        real(dp) :: eps
+        character, parameter :: nl = new_line('a')
+        character(len=*), parameter :: s = "# x(arb. units)     double (a.u.)           x²/5"//nl//&
+        &"#"//nl//&
+        &"# Multiline header"//nl//&
+        &"5.0000000000000000E-001    1.0000000000000000E+000  5.0000000000000003E-002"//nl//&
+        &nl//&
+        &"1.5555555555555556E+000  3.1111111111111112E+000  4.8395061728395061E-001"//nl//&
+        &"  2.6111111111111112E+000  5.2222222222222223E+000  1.3635802469135803E+000 # leading spaces"//nl//&
+        &"3.6666666666666670E+000 7.3333333333333339E+000  2.6888888888888891E+000 "//nl//&
+        &"4.7222222222222223E+000  9.4444444444444446E+000  4.4598765432098766E+000"//nl//&
+        &"# footer"
+        character(len=*), parameter :: fname = 'test_dp_comm_skip_max.dat'
+        real(dp) :: input(5, 3)
+        input = reshape([5.0000000000000000e-1_dp, 1.0000000000000000e+0_dp, 5.0000000000000000e-2_dp, &
+      & 1.5555555555555556e+0_dp, 3.1111111111111112e+0_dp, 4.8395061728395061e-1_dp,&
+      & 2.6111111111111112e+0_dp, 5.2222222222222223e+0_dp, 1.3635802469135803e+0_dp,&
+      & 3.6666666666666670e+0_dp, 7.3333333333333339e+0_dp, 2.6888888888888891e+0_dp,&
+      & 4.7222222222222223e+0_dp, 9.4444444444444446e+0_dp, 4.4598765432098766e+0_dp], &
+      & shape(input), order=[2, 1])
+
+        allocate (expected(10, 10))
+        eps = epsilon(1._dp)
+        open (newunit=unit, file=fname)
+        write (unit, '(a)') s
+        close (unit)
+
+        ! Test default values
+        call loadtxt(fname, expected) ! Read all data
+        call check(error, all(abs(input - expected) < eps), 'Default read failed')
+        if (allocated(error)) return
+        ! Test skiplines option
+        call loadtxt(fname, expected, skiplines=4) ! Skip comment and first line
+        call check(error, all(abs(input(2:, :) - expected) < eps), 'skiplines read failed')
+        if (allocated(error)) return
+        ! Test max_rows option
+        call loadtxt(fname, expected, max_rows=4) ! Skip comment and first line
+        call check(error, all(abs(input(:4, :) - expected) < eps), 'max_rows read failed')
+        if (allocated(error)) return
+        ! Test usecols option
+        call loadtxt(fname, expected, usecols=[3, 1, 1, 2]) ! Skip comment and first line
+        call check(error, all(abs(input(:, 3) - expected(:, 1)) < eps) &
+                   .or. all(abs(input(:, 1) - expected(:, 2)) < eps) &
+                   .or. all(abs(input(:, 1) - expected(:, 3)) < eps) &
+                   .or. all(abs(input(:, 2) - expected(:, 4)) < eps), 'usecols read failed')
+        if (allocated(error)) return
+
+    end subroutine test_loadtxt_dp_comm_skip_max
 
     subroutine test_loadtxt_dp_max_skip(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        real(dp), allocatable :: input(:,:), expected(:,:)
+        real(dp), allocatable :: input(:, :), expected(:, :)
         integer :: n, m
-        character(len=*), parameter :: FMT_REAL_DP = '(es24.16e3)'
-
-        allocate(input(10,10))
+        real(dp) :: eps
+        eps = epsilon(1._dp)
+        allocate (input(10, 10))
 
         do m = 0, 5
             do n = 1, 11
                 call random_number(input)
                 input = input - 0.5
                 call savetxt('test_dp_max_skip.txt', input)
-                call loadtxt('test_dp_max_skip.txt', expected, skiprows=m, max_rows=n)
-                call check(error, all(input(m+1:min(n+m,10),:) == expected),'Default format read failed')
-                if (allocated(error)) return
-                call loadtxt('test_dp_max_skip.txt', expected, skiprows=m, max_rows=n, fmt='*')
-                call check(error, all(input(m+1:min(n+m,10),:) == expected),'List directed read failed')
-                if (allocated(error)) return
-                call loadtxt('test_dp_max_skip.txt', expected, fmt="(*"//FMT_REAL_dp(1:len(FMT_REAL_dp)-1)//",1x))")
-                call check(error, all(input == expected),'User specified format failed')
-                deallocate(expected)
+                call loadtxt('test_dp_max_skip.txt', expected, skiplines=m, max_rows=n)
+                call check(error, all(abs(input(m + 1:min(n + m, 10), :) - expected) < eps),&
+                    &'max_rows and skiplines read failed')
                 if (allocated(error)) return
             end do
         end do
 
     end subroutine test_loadtxt_dp_max_skip
 
-
     subroutine test_loadtxt_dp_huge(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        real(dp), allocatable :: input(:,:), expected(:,:)
+        real(dp), allocatable :: input(:, :), expected(:, :)
         integer :: n
-        character(len=*), parameter :: FMT_REAL_DP = '(es24.16e3)'
+        integer, allocatable :: j(:)
+        real(dp) :: eps
+        eps = 10 * epsilon(1._dp) * huge(1._dp)
 
-        allocate(input(10,10))
-        allocate(expected(10,10))
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
 
         do n = 1, 10
             call random_number(input)
             input = (input - 0.5) * huge(input)
             call savetxt('test_dp_huge.txt', input)
             call loadtxt('test_dp_huge.txt', expected)
-            call check(error, all(input == expected),'Default format read failed')
-            if (allocated(error)) return
-            call loadtxt('test_dp_huge.txt', expected, fmt='*')
-            call check(error, all(input == expected),'List directed read failed')
-            if (allocated(error)) return
-            call loadtxt('test_dp_huge.txt', expected, fmt="(*"//FMT_REAL_dp(1:len(FMT_REAL_dp)-1)//",1x))")
-            call check(error, all(input == expected),'User specified format failed')
-            if (allocated(error)) return
+            call check(error, all(abs(input - expected) < eps), 'Huge read failed')
+            if (allocated(error)) then
+                ! j = maxloc(abs(input - expected))
+                ! print *, maxval(abs(input - expected)), input(j(1), j(2)), expected(j(1), j(2)), eps
+                return
+            end if
         end do
 
     end subroutine test_loadtxt_dp_huge
 
-
     subroutine test_loadtxt_dp_tiny(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        real(dp), allocatable :: input(:,:), expected(:,:)
+        real(dp), allocatable :: input(:, :), expected(:, :)
         integer :: n
-        character(len=*), parameter :: FMT_REAL_DP = '(es24.16e3)'
-        
-        allocate(input(10,10))
-        allocate(expected(10,10))
+        real(dp) :: eps
+        eps = epsilon(1._dp)
+
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
 
         do n = 1, 10
             call random_number(input)
             input = (input - 0.5) * tiny(input)
             call savetxt('test_dp_tiny.txt', input)
             call loadtxt('test_dp_tiny.txt', expected)
-            call check(error, all(input == expected),'Default format read failed')
-            if (allocated(error)) return
-            call loadtxt('test_dp_tiny.txt', expected, fmt='*')
-            call check(error, all(input == expected),'List directed read failed')
-            if (allocated(error)) return
-            call loadtxt('test_dp_tiny.txt', expected, fmt="(*"//FMT_REAL_dp(1:len(FMT_REAL_dp)-1)//",1x))")
-            call check(error, all(input == expected),'User specified format failed')
+            call check(error, all(abs(input - expected) < eps), 'Default format read failed')
             if (allocated(error)) return
         end do
 
     end subroutine test_loadtxt_dp_tiny
 
-
     subroutine test_loadtxt_complex(error)
         !> Error handling
         type(error_type), allocatable, intent(out) :: error
-        complex(dp), allocatable :: input(:,:), expected(:,:)
-        real(dp), allocatable :: re(:,:), im(:,:)
+        complex(dp), allocatable :: input(:, :), expected(:, :)
+        real(dp), allocatable :: re(:, :), im(:, :)
         integer :: n
-        character(len=*), parameter :: FMT_COMPLEX_DP = '(es24.16e3,1x,es24.16e3)'
+        real(dp) :: eps
+        eps = epsilon(1._dp)
 
-        allocate(re(10,10))
-        allocate(im(10,10))
-        allocate(input(10,10))
-        allocate(expected(10,10))
+        allocate (re(10, 10))
+        allocate (im(10, 10))
+        allocate (input(10, 10))
+        allocate (expected(10, 10))
 
         do n = 1, 10
             call random_number(re)
@@ -288,17 +297,17 @@ subroutine test_loadtxt_complex(error)
             input = cmplx(re, im)
             call savetxt('test_complex.txt', input)
             call loadtxt('test_complex.txt', expected)
-            call check(error, all(input == expected))
-            call loadtxt('test_complex.txt', expected, fmt="(*"//FMT_COMPLEX_dp(1:len(FMT_COMPLEX_dp)-1)//",1x))")
-            call check(error, all(input == expected))
+            call check(error, all(abs(input - expected) < eps))
             if (allocated(error)) return
+            !
             call savetxt('test_complex.txt', input, delimiter=',')
             call loadtxt('test_complex.txt', expected, delimiter=',')
-            call check(error, all(input == expected))
+            call check(error, all(abs(input - expected) < eps))
             if (allocated(error)) return
+            !
             call savetxt('test_complex.txt', input, delimiter=';')
             call loadtxt('test_complex.txt', expected, delimiter=';')
-            call check(error, all(input == expected))
+            call check(error, all(abs(input - expected) < eps))
             if (allocated(error)) return
         end do
 
@@ -306,11 +315,10 @@ end subroutine test_loadtxt_complex
 
 end module test_loadtxt
 
-
 program tester
-    use, intrinsic :: iso_fortran_env, only : error_unit
-    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
-    use test_loadtxt, only : collect_loadtxt
+    use, intrinsic :: iso_fortran_env, only: error_unit
+    use testdrive, only: run_testsuite, new_testsuite, testsuite_type
+    use test_loadtxt, only: collect_loadtxt
     implicit none
     integer :: stat, is
     type(testsuite_type), allocatable :: testsuites(:)
@@ -318,16 +326,16 @@ program tester
 
     stat = 0
     testsuites = [ &
-        new_testsuite("loadtxt", collect_loadtxt) &
-        ]
+                 new_testsuite("loadtxt", collect_loadtxt) &
+                 ]
 
     do is = 1, size(testsuites)
-        write(error_unit, fmt) "Testing:", testsuites(is)%name
+        write (error_unit, fmt) "Testing:", testsuites(is)%name
         call run_testsuite(testsuites(is)%collect, error_unit, stat)
     end do
 
     if (stat > 0) then
-        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+        write (error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
         error stop
     end if
 
diff --git a/test/io/test_loadtxt_qp.fypp b/test/io/test_loadtxt_qp.fypp
index fc041a60a..0e004db63 100644
--- a/test/io/test_loadtxt_qp.fypp
+++ b/test/io/test_loadtxt_qp.fypp
@@ -30,6 +30,8 @@ contains
 #:if WITH_QP
         real(qp), allocatable :: input(:,:), expected(:,:)
         integer :: n
+        real(qp) :: eps
+        eps = epsilon(1._qp)
 
         allocate(input(10,10))
         allocate(expected(10,10))
@@ -39,7 +41,7 @@ contains
             input = input - 0.5
             call savetxt('test_qp.txt', input)
             call loadtxt('test_qp.txt', expected)
-            call check(error, all(input == expected))
+            call check(error, all(abs(input - expected) < eps))
             if (allocated(error)) return
         end do
 #:else
@@ -55,6 +57,8 @@ contains
 #:if WITH_QP
         real(qp), allocatable :: input(:,:), expected(:,:)
         integer :: n
+        real(qp) :: eps
+        eps = 10 * epsilon(1._qp) * huge(1._qp)
 
         allocate(input(10,10))
         allocate(expected(10,10))
@@ -64,7 +68,7 @@ contains
             input = (input - 0.5) * huge(input)
             call savetxt('test_qp_huge.txt', input)
             call loadtxt('test_qp_huge.txt', expected)
-            call check(error, all(input == expected))
+            call check(error, all(abs(input - expected) < eps), 'Huge number read failed')
             if (allocated(error)) return
         end do
 #:else
@@ -80,6 +84,8 @@ contains
 #:if WITH_QP
         real(qp), allocatable :: input(:,:), expected(:,:)
         integer :: n
+        real(qp) :: eps
+        eps = epsilon(1._qp)
 
         allocate(input(10,10))
         allocate(expected(10,10))
@@ -89,7 +95,7 @@ contains
             input = (input - 0.5) * tiny(input)
             call savetxt('test_qp_tiny.txt', input)
             call loadtxt('test_qp_tiny.txt', expected)
-            call check(error, all(input == expected))
+            call check(error, all(abs(input - expected) < eps), 'Tiny numbers read failed')
             if (allocated(error)) return
         end do
 #:else