diff --git a/src/Makefile.manual b/src/Makefile.manual
index 8e912741b..7cc8816ea 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -37,6 +37,7 @@ clean:
 
 # Fortran module dependencies
 f18estop.o: stdlib_experimental_error.o
+stdlib_experimental_error.o: stdlib_experimental_optval.o
 stdlib_experimental_io.o: \
 	stdlib_experimental_error.o \
 	stdlib_experimental_optval.o \
diff --git a/src/stdlib_experimental_error.f90 b/src/stdlib_experimental_error.f90
index 3d932d6c9..21c6f501e 100644
--- a/src/stdlib_experimental_error.f90
+++ b/src/stdlib_experimental_error.f90
@@ -1,34 +1,64 @@
 module stdlib_experimental_error
-use, intrinsic :: iso_fortran_env, only: stderr=>error_unit
+use, intrinsic :: iso_fortran_env, only: stderr => error_unit
+use stdlib_experimental_optval, only: optval
 implicit none
 private
 
 interface ! f{08,18}estop.f90
-module subroutine error_stop(msg, code)
-character(*), intent(in) :: msg
-integer, intent(in), optional :: code
-end subroutine error_stop
+    module subroutine error_stop(msg, code)
+        character(*), intent(in) :: msg
+        integer, intent(in), optional :: code
+    end subroutine error_stop
 end interface
 
-public :: assert, error_stop
+public :: check, error_stop
 
 contains
 
-subroutine assert(condition, code)
-! If condition == .false., it aborts the program.
-!
-! Arguments
-! ---------
-!
-logical, intent(in) :: condition
-integer, intent(in), optional :: code
-!
-! Example
-! -------
-!
-! call assert(a == 5)
-
-if (.not. condition) call error_stop("Assert failed.", code)
-end subroutine
-
-end module
+subroutine check(condition, msg, code, warn)
+
+    ! Checks the value of a logical condition. If condition == .false. and:
+    !
+    !   * No other arguments are provided, it stops the program with the default
+    !     message and exit code 1;
+    !   * msg is provided, it prints the value of msg;
+    !   * code is provided, it stops the program with the given exit code;
+    !   * warn is provided and .true., it doesn't stop the program and prints
+    !   * the message.
+    !
+    ! Arguments
+    ! ---------
+
+    logical, intent(in) :: condition
+    character(*), intent(in), optional :: msg
+    integer, intent(in), optional :: code
+    logical, intent(in), optional :: warn
+    character(*), parameter :: msg_default = 'Check failed.'
+
+    ! Examples
+    ! --------
+    !
+    ! ! If a /= 5, stops the program with exit code 1
+    ! ! and prints 'Check failed.'
+    ! call check(a == 5)
+    !
+    ! ! As above, but prints 'a == 5 failed.'
+    ! call check(a == 5, msg='a == 5 failed.')
+    !
+    ! ! As above, but doesn't stop the program.
+    ! call check(a == 5, msg='a == 5 failed.', warn=.true.)
+    !
+    ! ! As example #2, but stops the program with exit code 77
+    ! call check(a == 5, msg='a == 5 failed.', code=77)
+
+    if (.not. condition) then
+        if (optval(warn, .false.)) then
+            write(stderr,*) optval(msg, msg_default)
+        else
+            call error_stop(optval(msg, msg_default), optval(code, 1))
+        end if
+    end if
+
+end subroutine check
+
+end module stdlib_experimental_error
diff --git a/src/tests/ascii/test_ascii.f90 b/src/tests/ascii/test_ascii.f90
index add28c698..3ff788fe0 100644
--- a/src/tests/ascii/test_ascii.f90
+++ b/src/tests/ascii/test_ascii.f90
@@ -1,6 +1,6 @@
 program test_ascii
 
-    use stdlib_experimental_error, only: assert
+    use stdlib_experimental_error, only: check
     use stdlib_experimental_ascii, only: lowercase, uppercase, digits, &
         octal_digits, fullhex_digits, hex_digits, lowerhex_digits, &
         whitespace, letters, is_alphanum, is_alpha, is_lower, is_upper, &
@@ -72,12 +72,12 @@ program test_ascii
 
     subroutine test_is_alphanum_short
         print *, "test_is_alphanum_short"
-        call assert(is_alphanum('A'))
-        call assert(is_alphanum('1'))
-        call assert(.not. is_alphanum('#'))
+        call check(is_alphanum('A'))
+        call check(is_alphanum('1'))
+        call check(.not. is_alphanum('#'))
 
         ! N.B.: does not return true for non-ASCII Unicode alphanumerics
-        call assert(.not. is_alphanum('á'))
+        call check(.not. is_alphanum('á'))
     end subroutine
 
     subroutine test_is_alphanum_long
@@ -88,23 +88,23 @@ subroutine test_is_alphanum_long
 
         clist = digits//octal_digits//fullhex_digits//letters//lowercase//uppercase
         do i = 1, len(clist)
-            call assert(is_alphanum(clist(i:i)))
+            call check(is_alphanum(clist(i:i)))
         end do
 
         clist = whitespace
         do i = 1, len(clist)
-            call assert(.not. is_alphanum(clist(i:i)))
+            call check(.not. is_alphanum(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_alpha_short
         print *, "test_is_alpha_short"
-        call assert(is_alpha('A'))
-        call assert(.not. is_alpha('1'))
-        call assert(.not. is_alpha('#'))
+        call check(is_alpha('A'))
+        call check(.not. is_alpha('1'))
+        call check(.not. is_alpha('#'))
 
         ! N.B.: does not return true for non-ASCII Unicode alphabetic characters
-        call assert(.not. is_alpha('á'))
+        call check(.not. is_alpha('á'))
     end subroutine
 
     subroutine test_is_alpha_long
@@ -115,24 +115,24 @@ subroutine test_is_alpha_long
 
         clist = letters//lowercase//uppercase
         do i = 1, len(clist)
-            call assert(is_alpha(clist(i:i)))
+            call check(is_alpha(clist(i:i)))
         end do
 
         clist = digits//octal_digits//whitespace
         do i = 1, len(clist)
-            call assert(.not. is_alpha(clist(i:i)))
+            call check(.not. is_alpha(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_lower_short
         print *, "test_is_lower_short"
-        call assert(is_lower('a'))
-        call assert(.not. is_lower('A'))
-        call assert(.not. is_lower('#'))
+        call check(is_lower('a'))
+        call check(.not. is_lower('A'))
+        call check(.not. is_lower('#'))
 
         ! N.B.: does not return true for non-ASCII Unicode lowercase letters
-        call assert(.not. is_lower('á'))
-        call assert(.not. is_lower('Á'))
+        call check(.not. is_lower('á'))
+        call check(.not. is_lower('Á'))
     end subroutine
 
     subroutine test_is_lower_long
@@ -141,24 +141,24 @@ subroutine test_is_lower_long
 
         print *, "test_is_lower_long"
         do i = 1, len(lowercase)
-            call assert(is_lower(lowercase(i:i)))
+            call check(is_lower(lowercase(i:i)))
         end do
 
         clist = digits//uppercase//whitespace
         do i = 1, len(clist)
-            call assert(.not. is_lower(clist(i:i)))
+            call check(.not. is_lower(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_upper_short
         print *, "test_is_upper_short"
-        call assert(is_upper('A'))
-        call assert(.not. is_upper('a'))
-        call assert(.not. is_upper('#'))
+        call check(is_upper('A'))
+        call check(.not. is_upper('a'))
+        call check(.not. is_upper('#'))
 
         ! N.B.: does not return true for non-ASCII Unicode uppercase letters
-        call assert(.not. is_upper('á'))
-        call assert(.not. is_upper('Á'))
+        call check(.not. is_upper('á'))
+        call check(.not. is_upper('Á'))
     end subroutine
 
     subroutine test_is_upper_long
@@ -166,26 +166,26 @@ subroutine test_is_upper_long
         character(len=:), allocatable :: clist
         print *, "test_is_upper_long"
         do i = 1, len(uppercase)
-            call assert(is_upper(uppercase(i:i)))
+            call check(is_upper(uppercase(i:i)))
         end do
 
         clist = digits//lowercase//whitespace
         do i = 1, len(clist)
-            call assert(.not. is_upper(clist(i:i)))
+            call check(.not. is_upper(clist(i:i)))
         end do
     end subroutine
 
 
     subroutine test_is_digit_short
         print *, "test_is_digit_short"
-        call assert(is_digit('3'))
-        call assert(is_digit('8'))
-        call assert(.not. is_digit('B'))
-        call assert(.not. is_digit('#'))
+        call check(is_digit('3'))
+        call check(is_digit('8'))
+        call check(.not. is_digit('B'))
+        call check(.not. is_digit('#'))
 
         ! N.B.: does not return true for non-ASCII Unicode numbers
-        call assert(.not. is_digit('０')) ! full-width digit zero (U+FF10)
-        call assert(.not. is_digit('４')) ! full-width digit four (U+FF14))
+        call check(.not. is_digit('０')) ! full-width digit zero (U+FF10)
+        call check(.not. is_digit('４')) ! full-width digit four (U+FF14))
     end subroutine
 
     subroutine test_is_digit_long
@@ -193,22 +193,22 @@ subroutine test_is_digit_long
         character(len=:), allocatable :: clist
         print *, "test_is_digit_long"
         do i = 1, len(digits)
-            call assert(is_digit(digits(i:i)))
+            call check(is_digit(digits(i:i)))
         end do
 
         clist = letters//whitespace
         do i = 1, len(clist)
-            call assert(.not. is_digit(clist(i:i)))
+            call check(.not. is_digit(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_octal_digit_short
         print *, "test_is_octal_digit_short"
-        call assert(is_octal_digit('0'))
-        call assert(is_octal_digit('7'))
-        call assert(.not. is_octal_digit('8'))
-        call assert(.not. is_octal_digit('A'))
-        call assert(.not. is_octal_digit('#'))
+        call check(is_octal_digit('0'))
+        call check(is_octal_digit('7'))
+        call check(.not. is_octal_digit('8'))
+        call check(.not. is_octal_digit('A'))
+        call check(.not. is_octal_digit('#'))
     end subroutine
 
     subroutine test_is_octal_digit_long
@@ -216,22 +216,22 @@ subroutine test_is_octal_digit_long
         character(len=:), allocatable :: clist
         print *, "test_is_octal_digit_long"
         do i = 1, len(octal_digits)
-            call assert(is_octal_digit(octal_digits(i:i)))
+            call check(is_octal_digit(octal_digits(i:i)))
         end do
         clist = letters//'89'//whitespace
         do i = 1, len(clist)
-            call assert(.not. is_octal_digit(clist(i:i)))
+            call check(.not. is_octal_digit(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_hex_digit_short
         print *, "test_is_hex_digit_short"
-        call assert(is_hex_digit('0'))
-        call assert(is_hex_digit('A'))
-        call assert(is_hex_digit('f')) !! lowercase hex digits are accepted
-        call assert(.not. is_hex_digit('g'))
-        call assert(.not. is_hex_digit('G'))
-        call assert(.not. is_hex_digit('#'))
+        call check(is_hex_digit('0'))
+        call check(is_hex_digit('A'))
+        call check(is_hex_digit('f')) !! lowercase hex digits are accepted
+        call check(.not. is_hex_digit('g'))
+        call check(.not. is_hex_digit('G'))
+        call check(.not. is_hex_digit('#'))
     end subroutine
 
     subroutine test_is_hex_digit_long
@@ -239,22 +239,22 @@ subroutine test_is_hex_digit_long
         character(len=:), allocatable :: clist
         print *, "test_is_hex_digit_long"
         do i = 1, len(fullhex_digits)
-            call assert(is_hex_digit(fullhex_digits(i:i)))
+            call check(is_hex_digit(fullhex_digits(i:i)))
         end do
         clist = lowercase(7:)//uppercase(7:)//whitespace
         do i = 1, len(clist)
-            call assert(.not. is_hex_digit(clist(i:i)))
+            call check(.not. is_hex_digit(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_white_short
         print *, "test_is_white_short"
-        call assert(is_white(' '))
-        call assert(is_white(TAB))
-        call assert(is_white(LF))
-        call assert(.not. is_white('1'))
-        call assert(.not. is_white('a'))
-        call assert(.not. is_white('#'))
+        call check(is_white(' '))
+        call check(is_white(TAB))
+        call check(is_white(LF))
+        call check(.not. is_white('1'))
+        call check(.not. is_white('a'))
+        call check(.not. is_white('#'))
     end subroutine
 
     subroutine test_is_white_long
@@ -262,21 +262,21 @@ subroutine test_is_white_long
         character(len=:), allocatable :: clist
         print *, "test_is_white_long"
         do i = 1, len(whitespace)
-            call assert(is_white(whitespace(i:i)))
+            call check(is_white(whitespace(i:i)))
         end do
         clist = digits//letters
         do i = 1, len(clist)
-            call assert(.not. is_white(clist(i:i)))
+            call check(.not. is_white(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_blank_short
         print *, "test_is_blank_short"
-        call assert(is_blank(' '))
-        call assert(is_blank(TAB))
-        call assert(.not. is_blank('1'))
-        call assert(.not. is_blank('a'))
-        call assert(.not. is_blank('#'))
+        call check(is_blank(' '))
+        call check(is_blank(TAB))
+        call check(.not. is_blank('1'))
+        call check(.not. is_blank('a'))
+        call check(.not. is_blank('#'))
     end subroutine
 
     subroutine test_is_blank_long
@@ -285,14 +285,14 @@ subroutine test_is_blank_long
         print *, "test_is_blank_long"
         do i = 1, len(whitespace)
             if (whitespace(i:i) == ' ' .or. whitespace(i:i) == TAB) then
-                call assert(is_blank(whitespace(i:i)))
+                call check(is_blank(whitespace(i:i)))
             else
-                call assert(.not. is_blank(whitespace(i:i)))
+                call check(.not. is_blank(whitespace(i:i)))
             end if
         end do
         clist = digits//letters
         do i = 1, len(clist)
-            call assert(.not. is_blank(clist(i:i)))
+            call check(.not. is_blank(clist(i:i)))
         end do
     end subroutine
 
@@ -300,11 +300,11 @@ subroutine test_is_control_short
         print *, "test_is_control_short"
         ! print *, is_control('\0')
         ! print *, is_control('\022')
-        call assert(is_control(new_line('a'))) ! newline is both whitespace and control
-        call assert(.not. is_control(' '))
-        call assert(.not. is_control('1'))
-        call assert(.not. is_control('a'))
-        call assert(.not. is_control('#'))
+        call check(is_control(new_line('a'))) ! newline is both whitespace and control
+        call check(.not. is_control(' '))
+        call check(.not. is_control('1'))
+        call check(.not. is_control('a'))
+        call check(.not. is_control('#'))
 
         ! N.B.: non-ASCII Unicode control characters are not recognized:
         ! print *, .not. is_control('\u0080')
@@ -317,32 +317,32 @@ subroutine test_is_control_long
         character(len=:), allocatable :: clist
         print *, "test_is_control_long"
         do i = 0, 31
-            call assert(is_control(achar(i)))
+            call check(is_control(achar(i)))
         end do
-        call assert(is_control(DEL))
+        call check(is_control(DEL))
 
         clist = digits//letters//' '
         do i = 1, len(clist)
-            call assert(.not. is_control(clist(i:i)))
+            call check(.not. is_control(clist(i:i)))
         end do
     end subroutine
 
     subroutine test_is_punctuation_short
         print *, "test_is_punctuation_short"
-        call assert(is_punctuation('.'))
-        call assert(is_punctuation(','))
-        call assert(is_punctuation(':'))
-        call assert(is_punctuation('!'))
-        call assert(is_punctuation('#'))
-        call assert(is_punctuation('~'))
-        call assert(is_punctuation('+'))
-        call assert(is_punctuation('_'))
-
-        call assert(.not. is_punctuation('1'))
-        call assert(.not. is_punctuation('a'))
-        call assert(.not. is_punctuation(' '))
-        call assert(.not. is_punctuation(LF)) ! new line character
-        call assert(.not. is_punctuation(NUL))
+        call check(is_punctuation('.'))
+        call check(is_punctuation(','))
+        call check(is_punctuation(':'))
+        call check(is_punctuation('!'))
+        call check(is_punctuation('#'))
+        call check(is_punctuation('~'))
+        call check(is_punctuation('+'))
+        call check(is_punctuation('_'))
+
+        call check(.not. is_punctuation('1'))
+        call check(.not. is_punctuation('a'))
+        call check(.not. is_punctuation(' '))
+        call check(.not. is_punctuation(LF)) ! new line character
+        call check(.not. is_punctuation(NUL))
 
         ! N.B.: Non-ASCII Unicode punctuation characters are not recognized.
         ! print *, is_punctuation('\u2012') ! (U+2012 = en-dash)
@@ -355,24 +355,24 @@ subroutine test_is_punctuation_long
         do i = 0, 127
             c = achar(i)
             if (is_control(c) .or. is_alphanum(c) .or. c == ' ') then
-                call assert(.not. is_punctuation(c))
+                call check(.not. is_punctuation(c))
             else
-                call assert(is_punctuation(c))
+                call check(is_punctuation(c))
             end if
         end do
     end subroutine
 
     subroutine test_is_graphical_short
         print *, "test_is_graphical"
-        call assert(is_graphical('1'))
-        call assert(is_graphical('a'))
-        call assert(is_graphical('#'))
-        call assert(.not. is_graphical(' ')) ! whitespace is not graphical
-        call assert(.not. is_graphical(LF))
-        call assert(.not. is_graphical(NUL))
+        call check(is_graphical('1'))
+        call check(is_graphical('a'))
+        call check(is_graphical('#'))
+        call check(.not. is_graphical(' ')) ! whitespace is not graphical
+        call check(.not. is_graphical(LF))
+        call check(.not. is_graphical(NUL))
 
         ! N.B.: Unicode graphical characters are not regarded as such.
-        call assert(.not. is_graphical('ä'))
+        call check(.not. is_graphical('ä'))
     end subroutine
 
     subroutine test_is_graphical_long
@@ -382,23 +382,23 @@ subroutine test_is_graphical_long
         do i = 0, 127
             c = achar(i)
             if (is_control(c) .or. c == ' ') then
-                call assert(.not. is_graphical(c))
+                call check(.not. is_graphical(c))
             else
-                call assert(is_graphical(c))
+                call check(is_graphical(c))
             end if
         end do
     end subroutine
 
     subroutine test_is_printable_short
         print *, "test_is_printable_short"
-        call assert(is_printable(' ')) ! whitespace is printable
-        call assert(is_printable('1'))
-        call assert(is_printable('a'))
-        call assert(is_printable('#'))
-        call assert(.not. is_printable(NUL)) ! control characters are not printable
+        call check(is_printable(' ')) ! whitespace is printable
+        call check(is_printable('1'))
+        call check(is_printable('a'))
+        call check(is_printable('#'))
+        call check(.not. is_printable(NUL)) ! control characters are not printable
 
         ! N.B.: Printable non-ASCII Unicode characters are not recognized.
-        call assert(.not. is_printable('ä'))
+        call check(.not. is_printable('ä'))
     end subroutine
 
     subroutine test_is_printable_long
@@ -408,34 +408,34 @@ subroutine test_is_printable_long
         do i = 0, 127
             c = achar(i)
             if (is_control(c)) then
-                call assert(.not. is_printable(c))
+                call check(.not. is_printable(c))
             else
-                call assert(is_printable(c))
+                call check(is_printable(c))
             end if
         end do
     end subroutine
 
     subroutine test_is_ascii_short()
         print *, "test_is_ascii_short"
-        call assert(is_ascii('a'))
-        call assert(.not. is_ascii('ä'))
+        call check(is_ascii('a'))
+        call check(.not. is_ascii('ä'))
     end subroutine
 
     subroutine test_is_ascii_long()
         integer :: i
         print *, "test_is_ascii_long"
         do i = 0, 127
-            call assert(is_ascii(achar(i)))
+            call check(is_ascii(achar(i)))
         end do
-        call assert(.not. is_ascii(achar(128))) ! raises compiler warning
+        call check(.not. is_ascii(achar(128))) ! raises compiler warning
 
     end subroutine
 
     subroutine test_to_lower_short()
         print *, "test_to_lower_short"
-        call assert(to_lower('a') == 'a')
-        call assert(to_lower('A') == 'a')
-        call assert(to_lower('#') == '#')
+        call check(to_lower('a') == 'a')
+        call check(to_lower('A') == 'a')
+        call check(to_lower('#') == '#')
     end subroutine
 
     subroutine test_to_lower_long()
@@ -443,23 +443,23 @@ subroutine test_to_lower_long()
         character(len=1) :: c
         print *, "test_to_lower_long"
         do i = 1, len(uppercase)
-            call assert(to_lower(uppercase(i:i)) == lowercase(i:i))
+            call check(to_lower(uppercase(i:i)) == lowercase(i:i))
         end do
         do i = 0, 127
             c = achar(i)
             if (c < 'A' .or. c > 'Z') then
-                call assert(to_lower(c) == c)
+                call check(to_lower(c) == c)
             else
-                call assert(to_lower(c) /= c)
+                call check(to_lower(c) /= c)
             end if
         end do
     end subroutine
 
     subroutine test_to_upper_short()
         print *, "test_to_upper_short"
-        call assert(to_upper('a') == 'A')
-        call assert(to_upper('A') == 'A')
-        call assert(to_upper('#') == '#')
+        call check(to_upper('a') == 'A')
+        call check(to_upper('A') == 'A')
+        call check(to_upper('#') == '#')
     end subroutine
 
     subroutine test_to_upper_long()
@@ -467,15 +467,15 @@ subroutine test_to_upper_long()
         character(len=1) :: c
         print *, "test_to_upper_long"
         do i = 1, len(lowercase)
-            call assert(to_upper(lowercase(i:i)) == uppercase(i:i))
+            call check(to_upper(lowercase(i:i)) == uppercase(i:i))
         end do
 
         do i = 0, 127
             c = achar(i)
             if (c < 'a' .or. c > 'z') then
-                call assert(to_upper(c) == c)
+                call check(to_upper(c) == c)
             else
-                call assert(to_upper(c) /= c)
+                call check(to_upper(c) /= c)
             end if
         end do
     end subroutine
@@ -539,6 +539,6 @@ pure logical function validation_func_interface(c)
             write(*,'(I3,2X,12(L4),2X,I3)') j, (table(j,i),i=1,12), count(table(j,:))
         end do
         write(*,'(5X,12(I4))') (count(table(:,i)),i=1,12)
-    end subroutine
+    end subroutine test_ascii_table
 
-end program
\ No newline at end of file
+end program test_ascii
diff --git a/src/tests/io/test_open.f90 b/src/tests/io/test_open.f90
index deaee3593..2f7805652 100644
--- a/src/tests/io/test_open.f90
+++ b/src/tests/io/test_open.f90
@@ -1,6 +1,6 @@
 program test_open
 use stdlib_experimental_io, only: open
-use stdlib_experimental_error, only: assert
+use stdlib_experimental_error, only: check
 implicit none
 
 character(:), allocatable :: filename
@@ -17,7 +17,7 @@ program test_open
 ! Test mode "r"
 u = open(filename, "r")
 read(u, *) a
-call assert(all(a == [1, 2, 3]))
+call check(all(a == [1, 2, 3]))
 close(u)
 
 ! Test mode "a"
@@ -26,9 +26,9 @@ program test_open
 close(u)
 u = open(filename, "r")
 read(u, *) a
-call assert(all(a == [1, 2, 3]))
+call check(all(a == [1, 2, 3]))
 read(u, *) a
-call assert(all(a == [4, 5, 6]))
+call check(all(a == [4, 5, 6]))
 close(u)
 
 
@@ -44,7 +44,7 @@ program test_open
 ! Test mode "r"
 u = open(filename, "rb")
 read(u) a
-call assert(all(a == [1, 2, 3]))
+call check(all(a == [1, 2, 3]))
 close(u)
 
 ! Test mode "a"
@@ -53,9 +53,9 @@ program test_open
 close(u)
 u = open(filename, "rb")
 read(u) a
-call assert(all(a == [1, 2, 3]))
+call check(all(a == [1, 2, 3]))
 read(u) a
-call assert(all(a == [4, 5, 6]))
+call check(all(a == [4, 5, 6]))
 close(u)
 
 
@@ -64,21 +64,21 @@ program test_open
 filename = get_outpath() // "/io_open.stream"
 
 u = open(filename, "rb", io)
-call assert(io == 0)
+call check(io == 0)
 if (io == 0) close(u)
 
 u = open(filename, "ab", io)
-call assert(io == 0)
+call check(io == 0)
 if (io == 0) close(u)
 
 
 filename = get_outpath() // "/does_not_exist.error"
 
 u = open(filename, "a", io)
-call assert(io /= 0)
+call check(io /= 0)
 
 u = open(filename, "r", io)
-call assert(io /= 0)
+call check(io /= 0)
 
 
 contains
diff --git a/src/tests/io/test_parse_mode.f90 b/src/tests/io/test_parse_mode.f90
index 9d7d68b63..add7fcfd7 100644
--- a/src/tests/io/test_parse_mode.f90
+++ b/src/tests/io/test_parse_mode.f90
@@ -1,6 +1,6 @@
 program test_parse_mode
 use stdlib_experimental_io, only: parse_mode
-use stdlib_experimental_error, only: assert
+use stdlib_experimental_error, only: check
 implicit none
 
 call test_parse_mode_expected_order()
@@ -16,156 +16,156 @@ program test_parse_mode
     subroutine test_parse_mode_expected_order()
     character(3) :: m
     m = parse_mode("")
-    call assert(m == "r t")
+    call check(m == "r t")
 
     m = parse_mode("r")
-    call assert(m == "r t")
+    call check(m == "r t")
     m = parse_mode("w")
-    call assert(m == "w t")
+    call check(m == "w t")
     m = parse_mode("a")
-    call assert(m == "a t")
+    call check(m == "a t")
     m = parse_mode("x")
-    call assert(m == "x t")
+    call check(m == "x t")
 
     m = parse_mode("rt")
-    call assert(m == "r t")
+    call check(m == "r t")
     m = parse_mode("wt")
-    call assert(m == "w t")
+    call check(m == "w t")
     m = parse_mode("at")
-    call assert(m == "a t")
+    call check(m == "a t")
     m = parse_mode("xt")
-    call assert(m == "x t")
+    call check(m == "x t")
 
     m = parse_mode("rb")
-    call assert(m == "r b")
+    call check(m == "r b")
     m = parse_mode("wb")
-    call assert(m == "w b")
+    call check(m == "w b")
     m = parse_mode("ab")
-    call assert(m == "a b")
+    call check(m == "a b")
     m = parse_mode("xb")
-    call assert(m == "x b")
+    call check(m == "x b")
 
     m = parse_mode("r+")
-    call assert(m == "r+t")
+    call check(m == "r+t")
     m = parse_mode("w+")
-    call assert(m == "w+t")
+    call check(m == "w+t")
     m = parse_mode("a+")
-    call assert(m == "a+t")
+    call check(m == "a+t")
     m = parse_mode("x+")
-    call assert(m == "x+t")
+    call check(m == "x+t")
 
     m = parse_mode("r+t")
-    call assert(m == "r+t")
+    call check(m == "r+t")
     m = parse_mode("w+t")
-    call assert(m == "w+t")
+    call check(m == "w+t")
     m = parse_mode("a+t")
-    call assert(m == "a+t")
+    call check(m == "a+t")
     m = parse_mode("x+t")
-    call assert(m == "x+t")
+    call check(m == "x+t")
 
     m = parse_mode("r+b")
-    call assert(m == "r+b")
+    call check(m == "r+b")
     m = parse_mode("w+b")
-    call assert(m == "w+b")
+    call check(m == "w+b")
     m = parse_mode("a+b")
-    call assert(m == "a+b")
+    call check(m == "a+b")
     m = parse_mode("x+b")
-    call assert(m == "x+b")
+    call check(m == "x+b")
 
     end subroutine
 
     subroutine test_parse_mode_reverse_order()
     character(3) :: m
     m = parse_mode("")
-    call assert(m == "r t")
+    call check(m == "r t")
 
     m = parse_mode("tr")
-    call assert(m == "r t")
+    call check(m == "r t")
     m = parse_mode("tw")
-    call assert(m == "w t")
+    call check(m == "w t")
     m = parse_mode("ta")
-    call assert(m == "a t")
+    call check(m == "a t")
     m = parse_mode("tx")
-    call assert(m == "x t")
+    call check(m == "x t")
 
     m = parse_mode("br")
-    call assert(m == "r b")
+    call check(m == "r b")
     m = parse_mode("bw")
-    call assert(m == "w b")
+    call check(m == "w b")
     m = parse_mode("ba")
-    call assert(m == "a b")
+    call check(m == "a b")
     m = parse_mode("bx")
-    call assert(m == "x b")
+    call check(m == "x b")
 
     m = parse_mode("+r")
-    call assert(m == "r+t")
+    call check(m == "r+t")
     m = parse_mode("+w")
-    call assert(m == "w+t")
+    call check(m == "w+t")
     m = parse_mode("+a")
-    call assert(m == "a+t")
+    call check(m == "a+t")
     m = parse_mode("+x")
-    call assert(m == "x+t")
+    call check(m == "x+t")
 
     m = parse_mode("t+r")
-    call assert(m == "r+t")
+    call check(m == "r+t")
     m = parse_mode("t+w")
-    call assert(m == "w+t")
+    call check(m == "w+t")
     m = parse_mode("t+a")
-    call assert(m == "a+t")
+    call check(m == "a+t")
     m = parse_mode("t+x")
-    call assert(m == "x+t")
+    call check(m == "x+t")
 
     m = parse_mode("b+r")
-    call assert(m == "r+b")
+    call check(m == "r+b")
     m = parse_mode("b+w")
-    call assert(m == "w+b")
+    call check(m == "w+b")
     m = parse_mode("b+a")
-    call assert(m == "a+b")
+    call check(m == "a+b")
     m = parse_mode("x+b")
-    call assert(m == "x+b")
+    call check(m == "x+b")
 
     end subroutine
 
     subroutine test_parse_mode_random_order()
     character(3) :: m
     m = parse_mode("")
-    call assert(m == "r t")
+    call check(m == "r t")
 
     m = parse_mode("t r")
-    call assert(m == "r t")
+    call check(m == "r t")
     m = parse_mode(" tw ")
-    call assert(m == "w t")
+    call check(m == "w t")
     m = parse_mode("ta  ")
-    call assert(m == "a t")
+    call check(m == "a t")
     m = parse_mode("  t   x   ")
-    call assert(m == "x t")
+    call check(m == "x t")
 
     m = parse_mode("+ r ")
-    call assert(m == "r+t")
+    call check(m == "r+t")
     m = parse_mode("w   +")
-    call assert(m == "w+t")
+    call check(m == "w+t")
     m = parse_mode(" a+")
-    call assert(m == "a+t")
+    call check(m == "a+t")
     m = parse_mode(" x+   t  ")
-    call assert(m == "x+t")
+    call check(m == "x+t")
 
     m = parse_mode("tr+ ")
-    call assert(m == "r+t")
+    call check(m == "r+t")
     m = parse_mode("wt + ")
-    call assert(m == "w+t")
+    call check(m == "w+t")
     m = parse_mode("a + t")
-    call assert(m == "a+t")
+    call check(m == "a+t")
     m = parse_mode(" xt + ")
-    call assert(m == "x+t")
+    call check(m == "x+t")
 
     m = parse_mode(" + t")
-    call assert(m == "r+t")
+    call check(m == "r+t")
     m = parse_mode(" +w  b")
-    call assert(m == "w+b")
+    call check(m == "w+b")
     m = parse_mode("a + b")
-    call assert(m == "a+b")
+    call check(m == "a+b")
     m = parse_mode(" b + x  ")
-    call assert(m == "x+b")
+    call check(m == "x+b")
 
     end subroutine
 
@@ -173,13 +173,13 @@ subroutine test_parse_mode_always_fail()
     character(3) :: m
 
     m = parse_mode("r+w")
-    call assert(m /= "r t")
+    call check(m /= "r t")
 
     m = parse_mode("tt")
-    call assert(m /= "r t")
+    call check(m /= "r t")
 
     m = parse_mode("bt")
-    call assert(m /= "r t")
+    call check(m /= "r t")
 
     end subroutine
 
diff --git a/src/tests/io/test_savetxt.f90 b/src/tests/io/test_savetxt.f90
index 8b7d13fcc..8ffd4d6be 100644
--- a/src/tests/io/test_savetxt.f90
+++ b/src/tests/io/test_savetxt.f90
@@ -1,7 +1,7 @@
 program test_savetxt
 use stdlib_experimental_kinds, only: int32, sp, dp
 use stdlib_experimental_io, only: loadtxt, savetxt
-use stdlib_experimental_error, only: assert
+use stdlib_experimental_error, only: check
 implicit none
 
 character(:), allocatable :: outpath
@@ -36,14 +36,14 @@ subroutine test_iint32(outpath)
     d = reshape([1, 2, 3, 4, 5, 6], [3, 2])
     call savetxt(outpath, d)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [3, 2]))
-    call assert(all(abs(d-d2) == 0))
+    call check(all(shape(d2) == [3, 2]))
+    call check(all(abs(d-d2) == 0))
 
     e = reshape([1, 2, 3, 4, 5, 6], [2, 3])
     call savetxt(outpath, e)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [2, 3]))
-    call assert(all(abs(e-d2) == 0))
+    call check(all(shape(d2) == [2, 3]))
+    call check(all(abs(e-d2) == 0))
     end subroutine
 
 
@@ -54,14 +54,14 @@ subroutine test_rsp(outpath)
     d = reshape([1, 2, 3, 4, 5, 6], [3, 2])
     call savetxt(outpath, d)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [3, 2]))
-    call assert(all(abs(d-d2) < epsilon(1._sp)))
+    call check(all(shape(d2) == [3, 2]))
+    call check(all(abs(d-d2) < epsilon(1._sp)))
 
     e = reshape([1, 2, 3, 4, 5, 6], [2, 3])
     call savetxt(outpath, e)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [2, 3]))
-    call assert(all(abs(e-d2) < epsilon(1._sp)))
+    call check(all(shape(d2) == [2, 3]))
+    call check(all(abs(e-d2) < epsilon(1._sp)))
     end subroutine test_rsp
 
 
@@ -72,14 +72,14 @@ subroutine test_rdp(outpath)
     d = reshape([1, 2, 3, 4, 5, 6], [3, 2])
     call savetxt(outpath, d)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [3, 2]))
-    call assert(all(abs(d-d2) < epsilon(1._dp)))
+    call check(all(shape(d2) == [3, 2]))
+    call check(all(abs(d-d2) < epsilon(1._dp)))
 
     e = reshape([1, 2, 3, 4, 5, 6], [2, 3])
     call savetxt(outpath, e)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [2, 3]))
-    call assert(all(abs(e-d2) < epsilon(1._dp)))
+    call check(all(shape(d2) == [2, 3]))
+    call check(all(abs(e-d2) < epsilon(1._dp)))
     end subroutine test_rdp
 
     subroutine test_csp(outpath)
@@ -89,14 +89,14 @@ subroutine test_csp(outpath)
     d = cmplx(1, 1)* reshape([1, 2, 3, 4, 5, 6], [3, 2])
     call savetxt(outpath, d)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [3, 2]))
-    call assert(all(abs(d-d2) < epsilon(1._sp)))
+    call check(all(shape(d2) == [3, 2]))
+    call check(all(abs(d-d2) < epsilon(1._sp)))
 
     e = cmplx(1, 1)* reshape([1, 2, 3, 4, 5, 6], [2, 3])
     call savetxt(outpath, e)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [2, 3]))
-    call assert(all(abs(e-d2) < epsilon(1._sp)))
+    call check(all(shape(d2) == [2, 3]))
+    call check(all(abs(e-d2) < epsilon(1._sp)))
     end subroutine test_csp
 
     subroutine test_cdp(outpath)
@@ -106,14 +106,14 @@ subroutine test_cdp(outpath)
     d = cmplx(1._dp, 1._dp)* reshape([1, 2, 3, 4, 5, 6], [3, 2])
     call savetxt(outpath, d)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [3, 2]))
-    call assert(all(abs(d-d2) < epsilon(1._dp)))
+    call check(all(shape(d2) == [3, 2]))
+    call check(all(abs(d-d2) < epsilon(1._dp)))
 
     e = cmplx(1, 1)* reshape([1, 2, 3, 4, 5, 6], [2, 3])
     call savetxt(outpath, e)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [2, 3]))
-    call assert(all(abs(e-d2) < epsilon(1._dp)))
+    call check(all(shape(d2) == [2, 3]))
+    call check(all(abs(e-d2) < epsilon(1._dp)))
     end subroutine test_cdp
 
 end program test_savetxt
diff --git a/src/tests/io/test_savetxt_qp.f90 b/src/tests/io/test_savetxt_qp.f90
index 8ebb7151a..2c93f7093 100644
--- a/src/tests/io/test_savetxt_qp.f90
+++ b/src/tests/io/test_savetxt_qp.f90
@@ -1,7 +1,7 @@
 program test_savetxt_qp
 use stdlib_experimental_kinds, only: qp
 use stdlib_experimental_io, only: loadtxt, savetxt
-use stdlib_experimental_error, only: assert
+use stdlib_experimental_error, only: check
 implicit none
 
 character(:), allocatable :: outpath
@@ -33,14 +33,14 @@ subroutine test_rqp(outpath)
     d = reshape([1, 2, 3, 4, 5, 6], [3, 2])
     call savetxt(outpath, d)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [3, 2]))
-    call assert(all(abs(d-d2) < epsilon(1._qp)))
+    call check(all(shape(d2) == [3, 2]))
+    call check(all(abs(d-d2) < epsilon(1._qp)))
 
     e = reshape([1, 2, 3, 4, 5, 6], [2, 3])
     call savetxt(outpath, e)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [2, 3]))
-    call assert(all(abs(e-d2) < epsilon(1._qp)))
+    call check(all(shape(d2) == [2, 3]))
+    call check(all(abs(e-d2) < epsilon(1._qp)))
     end subroutine test_rqp
 
     subroutine test_cqp(outpath)
@@ -50,14 +50,14 @@ subroutine test_cqp(outpath)
     d = reshape([1, 2, 3, 4, 5, 6], [3, 2])
     call savetxt(outpath, d)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [3, 2]))
-    call assert(all(abs(d-d2) < epsilon(1._qp)))
+    call check(all(shape(d2) == [3, 2]))
+    call check(all(abs(d-d2) < epsilon(1._qp)))
 
     e = reshape([1, 2, 3, 4, 5, 6], [2, 3])
     call savetxt(outpath, e)
     call loadtxt(outpath, d2)
-    call assert(all(shape(d2) == [2, 3]))
-    call assert(all(abs(e-d2) < epsilon(1._qp)))
+    call check(all(shape(d2) == [2, 3]))
+    call check(all(abs(e-d2) < epsilon(1._qp)))
     end subroutine test_cqp
 
 end program test_savetxt_qp
diff --git a/src/tests/optval/test_optval.f90 b/src/tests/optval/test_optval.f90
index 00a40a8fa..d9a8f3a23 100644
--- a/src/tests/optval/test_optval.f90
+++ b/src/tests/optval/test_optval.f90
@@ -2,7 +2,7 @@ program test_optval
   use, intrinsic :: iso_fortran_env, only: &
        sp => real32, dp => real64, qp => real128, &
        int8, int16, int32, int64
-  use stdlib_experimental_error, only: assert
+  use stdlib_experimental_error, only: check
   use stdlib_experimental_optval, only: optval
 
   implicit none
@@ -42,8 +42,8 @@ program test_optval
 
   subroutine test_optval_rsp
     print *, "test_optval_rsp"
-    call assert(foo_sp(1.0_sp) == 1.0_sp)
-    call assert(foo_sp() == 2.0_sp)
+    call check(foo_sp(1.0_sp) == 1.0_sp)
+    call check(foo_sp() == 2.0_sp)
   end subroutine test_optval_rsp
 
 
@@ -56,8 +56,8 @@ function foo_sp(x) result(z)
 
   subroutine test_optval_rdp
     print *, "test_optval_rdp"
-    call assert(foo_dp(1.0_dp) == 1.0_dp)
-    call assert(foo_dp() == 2.0_dp)
+    call check(foo_dp(1.0_dp) == 1.0_dp)
+    call check(foo_dp() == 2.0_dp)
   end subroutine test_optval_rdp
 
 
@@ -70,8 +70,8 @@ function foo_dp(x) result(z)
 
   subroutine test_optval_rqp
     print *, "test_optval_rqp"
-    call assert(foo_qp(1.0_qp) == 1.0_qp)
-    call assert(foo_qp() == 2.0_qp)
+    call check(foo_qp(1.0_qp) == 1.0_qp)
+    call check(foo_qp() == 2.0_qp)
   end subroutine test_optval_rqp
 
 
@@ -86,8 +86,8 @@ subroutine test_optval_csp
     complex(sp) :: z1
     print *, "test_optval_csp"
     z1 = cmplx(1.0_sp, 2.0_sp, kind=sp)
-    call assert(foo_csp(z1) == z1)
-    call assert(foo_csp() == z1)
+    call check(foo_csp(z1) == z1)
+    call check(foo_csp() == z1)
   end subroutine test_optval_csp
 
   function foo_csp(x) result(z)
@@ -101,8 +101,8 @@ subroutine test_optval_cdp
     complex(dp) :: z1
     print *, "test_optval_cdp"
     z1 = cmplx(1.0_dp, 2.0_dp)
-    call assert(foo_cdp(z1) == z1)
-    call assert(foo_cdp() == z1)
+    call check(foo_cdp(z1) == z1)
+    call check(foo_cdp() == z1)
   end subroutine test_optval_cdp
 
   function foo_cdp(x) result(z)
@@ -116,8 +116,8 @@ subroutine test_optval_cqp
     complex(qp) :: z1
     print *, "test_optval_cqp"
     z1 = cmplx(1.0_qp, 2.0_qp, kind=qp)
-    call assert(foo_cqp(z1) == z1)
-    call assert(foo_cqp() == z1)
+    call check(foo_cqp(z1) == z1)
+    call check(foo_cqp() == z1)
   end subroutine test_optval_cqp
 
   function foo_cqp(x) result(z)
@@ -129,8 +129,8 @@ function foo_cqp(x) result(z)
 
   subroutine test_optval_iint8
     print *, "test_optval_iint8"
-    call assert(foo_int8(1_int8) == 1_int8)
-    call assert(foo_int8() == 2_int8)
+    call check(foo_int8(1_int8) == 1_int8)
+    call check(foo_int8() == 2_int8)
   end subroutine test_optval_iint8
 
 
@@ -143,8 +143,8 @@ function foo_int8(x) result(z)
 
   subroutine test_optval_iint16
     print *, "test_optval_iint16"
-    call assert(foo_int16(1_int16) == 1_int16)
-    call assert(foo_int16() == 2_int16)
+    call check(foo_int16(1_int16) == 1_int16)
+    call check(foo_int16() == 2_int16)
   end subroutine test_optval_iint16
 
 
@@ -157,8 +157,8 @@ function foo_int16(x) result(z)
 
   subroutine test_optval_iint32
     print *, "test_optval_iint32"
-    call assert(foo_int32(1_int32) == 1_int32)
-    call assert(foo_int32() == 2_int32)
+    call check(foo_int32(1_int32) == 1_int32)
+    call check(foo_int32() == 2_int32)
   end subroutine test_optval_iint32
 
 
@@ -171,8 +171,8 @@ function foo_int32(x) result(z)
 
   subroutine test_optval_iint64
     print *, "test_optval_int64"
-    call assert(foo_int64(1_int64) == 1_int64)
-    call assert(foo_int64() == 2_int64)
+    call check(foo_int64(1_int64) == 1_int64)
+    call check(foo_int64() == 2_int64)
   end subroutine test_optval_iint64
 
 
@@ -185,8 +185,8 @@ function foo_int64(x) result(z)
 
   subroutine test_optval_logical
     print *, "test_optval_logical"
-    call assert(foo_logical(.true.))
-    call assert(.not.foo_logical())
+    call check(foo_logical(.true.))
+    call check(.not.foo_logical())
   end subroutine test_optval_logical
 
 
@@ -199,8 +199,8 @@ function foo_logical(x) result(z)
 
   subroutine test_optval_character
     print *, "test_optval_character"
-    call assert(foo_character("x") == "x")
-    call assert(foo_character() == "y")
+    call check(foo_character("x") == "x")
+    call check(foo_character() == "y")
   end subroutine test_optval_character
 
 
@@ -213,8 +213,8 @@ function foo_character(x) result(z)
 
   subroutine test_optval_rsp_arr
     print *, "test_optval_rsp_arr"
-    call assert(all(foo_sp_arr([1.0_sp, -1.0_sp]) == [1.0_sp, -1.0_sp]))
-    call assert(all(foo_sp_arr() == [2.0_sp, -2.0_sp]))
+    call check(all(foo_sp_arr([1.0_sp, -1.0_sp]) == [1.0_sp, -1.0_sp]))
+    call check(all(foo_sp_arr() == [2.0_sp, -2.0_sp]))
   end subroutine test_optval_rsp_arr
 
 
@@ -227,8 +227,8 @@ end function foo_sp_arr
 
   subroutine test_optval_rdp_arr
     print *, "test_optval_rdp_arr"
-    call assert(all(foo_dp_arr([1.0_dp, -1.0_dp]) == [1.0_dp, -1.0_dp]))
-    call assert(all(foo_dp_arr() == [2.0_dp, -2.0_dp]))
+    call check(all(foo_dp_arr([1.0_dp, -1.0_dp]) == [1.0_dp, -1.0_dp]))
+    call check(all(foo_dp_arr() == [2.0_dp, -2.0_dp]))
   end subroutine test_optval_rdp_arr
 
 
@@ -241,8 +241,8 @@ end function foo_dp_arr
 
   subroutine test_optval_rqp_arr
     print *, "test_optval_qp_arr"
-    call assert(all(foo_qp_arr([1.0_qp, -1.0_qp]) == [1.0_qp, -1.0_qp]))
-    call assert(all(foo_qp_arr() == [2.0_qp, -2.0_qp]))
+    call check(all(foo_qp_arr([1.0_qp, -1.0_qp]) == [1.0_qp, -1.0_qp]))
+    call check(all(foo_qp_arr() == [2.0_qp, -2.0_qp]))
   end subroutine test_optval_rqp_arr
 
 
@@ -258,8 +258,8 @@ subroutine test_optval_csp_arr
     print *, "test_optval_csp_arr"
     z1 = cmplx(1.0_sp, 2.0_sp, kind=sp)*[1.0_sp, -1.0_sp]
     z2 = cmplx(2.0_sp, 2.0_sp, kind=sp)*[1.0_sp, -1.0_sp]
-    call assert(all(foo_csp_arr(z1) == z1))
-    call assert(all(foo_csp_arr() == z2))
+    call check(all(foo_csp_arr(z1) == z1))
+    call check(all(foo_csp_arr() == z2))
   end subroutine test_optval_csp_arr
 
 
@@ -275,8 +275,8 @@ subroutine test_optval_cdp_arr
     print *, "test_optval_cdp_arr"
     z1 = cmplx(1.0_dp, 2.0_dp, kind=dp)*[1.0_dp, -1.0_dp]
     z2 = cmplx(2.0_dp, 2.0_dp, kind=dp)*[1.0_dp, -1.0_dp]
-    call assert(all(foo_cdp_arr(z1) == z1))
-    call assert(all(foo_cdp_arr() == z2))
+    call check(all(foo_cdp_arr(z1) == z1))
+    call check(all(foo_cdp_arr() == z2))
   end subroutine test_optval_cdp_arr
 
 
@@ -292,8 +292,8 @@ subroutine test_optval_cqp_arr
     print *, "test_optval_cqp_arr"
     z1 = cmplx(1.0_qp, 2.0_qp, kind=qp)*[1.0_qp, -1.0_qp]
     z2 = cmplx(2.0_qp, 2.0_qp, kind=qp)*[1.0_qp, -1.0_qp]
-    call assert(all(foo_cqp_arr(z1) == z1))
-    call assert(all(foo_cqp_arr() == z2))
+    call check(all(foo_cqp_arr(z1) == z1))
+    call check(all(foo_cqp_arr() == z2))
   end subroutine test_optval_cqp_arr
 
 
@@ -306,8 +306,8 @@ end function foo_cqp_arr
 
   subroutine test_optval_iint8_arr
     print *, "test_optval_int8_arr"
-    call assert(all(foo_int8_arr([1_int8, -1_int8]) == [1_int8, -1_int8]))
-    call assert(all(foo_int8_arr() == [2_int8, -2_int8]))
+    call check(all(foo_int8_arr([1_int8, -1_int8]) == [1_int8, -1_int8]))
+    call check(all(foo_int8_arr() == [2_int8, -2_int8]))
   end subroutine test_optval_iint8_arr
 
 
@@ -320,8 +320,8 @@ end function foo_int8_arr
 
   subroutine test_optval_iint16_arr
     print *, "test_optval_int16_arr"
-    call assert(all(foo_int16_arr([1_int16, -1_int16]) == [1_int16, -1_int16]))
-    call assert(all(foo_int16_arr() == [2_int16, -2_int16]))
+    call check(all(foo_int16_arr([1_int16, -1_int16]) == [1_int16, -1_int16]))
+    call check(all(foo_int16_arr() == [2_int16, -2_int16]))
   end subroutine test_optval_iint16_arr
 
 
@@ -334,8 +334,8 @@ end function foo_int16_arr
 
   subroutine test_optval_iint32_arr
     print *, "test_optval_int32_arr"
-    call assert(all(foo_int32_arr([1_int32, -1_int32]) == [1_int32, -1_int32]))
-    call assert(all(foo_int32_arr() == [2_int32, -2_int32]))
+    call check(all(foo_int32_arr([1_int32, -1_int32]) == [1_int32, -1_int32]))
+    call check(all(foo_int32_arr() == [2_int32, -2_int32]))
   end subroutine test_optval_iint32_arr
 
 
@@ -348,8 +348,8 @@ end function foo_int32_arr
 
   subroutine test_optval_iint64_arr
     print *, "test_optval_int64_arr"
-    call assert(all(foo_int64_arr([1_int64, -1_int64]) == [1_int64, -1_int64]))
-    call assert(all(foo_int64_arr() == [2_int64, -2_int64]))
+    call check(all(foo_int64_arr([1_int64, -1_int64]) == [1_int64, -1_int64]))
+    call check(all(foo_int64_arr() == [2_int64, -2_int64]))
   end subroutine test_optval_iint64_arr
 
 
@@ -362,8 +362,8 @@ end function foo_int64_arr
 
   subroutine test_optval_logical_arr
     print *, "test_optval_logical_arr"
-    call assert(all(foo_logical_arr()))
-    call assert(all(.not.foo_logical_arr()))
+    call check(all(foo_logical_arr()))
+    call check(all(.not.foo_logical_arr()))
   end subroutine test_optval_logical_arr
 
 
diff --git a/src/tests/quadrature/test_trapz.f90 b/src/tests/quadrature/test_trapz.f90
index 1265e4b6b..0cf58ddea 100644
--- a/src/tests/quadrature/test_trapz.f90
+++ b/src/tests/quadrature/test_trapz.f90
@@ -1,6 +1,6 @@
 program test_trapz
     use stdlib_experimental_kinds, only: sp, dp, qp
-    use stdlib_experimental_error, only: assert
+    use stdlib_experimental_error, only: check
     use stdlib_experimental_quadrature, only: trapz, trapz_weights
 
     implicit none
@@ -33,21 +33,21 @@ subroutine test_trapz_sp
 
         val = trapz(y, 1.0_sp)
         ans = 128.0_sp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         val = trapz(y, 0.5_sp)
         ans = 64.0_sp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = [((i-1)*4.0_sp/real(n-1, sp), i = 1, n)]
         val = trapz(y, x)
         ans = 32.0_sp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = y**2
         val = trapz(y, x)
         ans = 2728.0_sp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
     end subroutine test_trapz_sp
 
     subroutine test_trapz_dp
@@ -64,21 +64,21 @@ subroutine test_trapz_dp
 
         val = trapz(y, 1.0_dp)
         ans = 128.0_dp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         val = trapz(y, 0.5_dp)
         ans = 64.0_dp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = [((i-1)*4.0_dp/real(n-1, dp), i = 1, n)]
         val = trapz(y, x)
         ans = 32.0_dp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = y**2
         val = trapz(y, x)
         ans = 2728.0_sp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
     end subroutine test_trapz_dp
 
 
@@ -96,21 +96,21 @@ subroutine test_trapz_qp
 
         val = trapz(y, 1.0_qp)
         ans = 128.0_qp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         val = trapz(y, 0.5_qp)
         ans = 64.0_qp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = [((i-1)*4.0_qp/real(n-1, qp), i = 1, n)]
         val = trapz(y, x)
         ans = 32.0_qp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = y**2
         val = trapz(y, x)
         ans = 2728.0_qp
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
     end subroutine test_trapz_qp
 
 
@@ -131,13 +131,13 @@ subroutine test_trapz_weights_sp
         w = trapz_weights(x)
         val = dot_product(w, y)
         ans = trapz(y, x)
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = y**2
         w = trapz_weights(x)
         val = dot_product(w, y)
         ans = trapz(y, x)
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
     end subroutine test_trapz_weights_sp
 
@@ -159,13 +159,13 @@ subroutine test_trapz_weights_dp
         w = trapz_weights(x)
         val = dot_product(w, y)
         ans = trapz(y, x)
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = y**2
         w = trapz_weights(x)
         val = dot_product(w, y)
         ans = trapz(y, x)
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
     end subroutine test_trapz_weights_dp
 
@@ -187,13 +187,13 @@ subroutine test_trapz_weights_qp
         w = trapz_weights(x)
         val = dot_product(w, y)
         ans = trapz(y, x)
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
         x = y**2
         w = trapz_weights(x)
         val = dot_product(w, y)
         ans = trapz(y, x)
-        call assert(abs(val - ans) < epsilon(ans))
+        call check(abs(val - ans) < epsilon(ans))
 
     end subroutine test_trapz_weights_qp
 
@@ -203,10 +203,10 @@ subroutine test_trapz_zero_sp
 
         print *, "test_trapz_zero_sp"
 
-        call assert(abs(trapz(a, 1.0_sp)) < epsilon(0.0_sp))
-        call assert(abs(trapz([1.0_sp], 1.0_sp)) < epsilon(0.0_sp))
-        call assert(abs(trapz(a, a)) < epsilon(0.0_sp))
-        call assert(abs(trapz([1.0_sp], [1.0_sp])) < epsilon(0.0_sp))
+        call check(abs(trapz(a, 1.0_sp)) < epsilon(0.0_sp))
+        call check(abs(trapz([1.0_sp], 1.0_sp)) < epsilon(0.0_sp))
+        call check(abs(trapz(a, a)) < epsilon(0.0_sp))
+        call check(abs(trapz([1.0_sp], [1.0_sp])) < epsilon(0.0_sp))
     end subroutine test_trapz_zero_sp
 
 
@@ -215,10 +215,10 @@ subroutine test_trapz_zero_dp
 
         print *, "test_trapz_zero_dp"
 
-        call assert(abs(trapz(a, 1.0_dp)) < epsilon(0.0_dp))
-        call assert(abs(trapz([1.0_dp], 1.0_dp)) < epsilon(0.0_dp))
-        call assert(abs(trapz(a, a)) < epsilon(0.0_dp))
-        call assert(abs(trapz([1.0_dp], [1.0_dp])) < epsilon(0.0_dp))
+        call check(abs(trapz(a, 1.0_dp)) < epsilon(0.0_dp))
+        call check(abs(trapz([1.0_dp], 1.0_dp)) < epsilon(0.0_dp))
+        call check(abs(trapz(a, a)) < epsilon(0.0_dp))
+        call check(abs(trapz([1.0_dp], [1.0_dp])) < epsilon(0.0_dp))
     end subroutine test_trapz_zero_dp
 
 
@@ -227,10 +227,10 @@ subroutine test_trapz_zero_qp
 
         print *, "test_trapz_zero_qp"
 
-        call assert(abs(trapz(a, 1.0_qp)) < epsilon(0.0_qp))
-        call assert(abs(trapz([1.0_qp], 1.0_qp)) < epsilon(0.0_qp))
-        call assert(abs(trapz(a, a)) < epsilon(0.0_qp))
-        call assert(abs(trapz([1.0_qp], [1.0_qp])) < epsilon(0.0_qp))
+        call check(abs(trapz(a, 1.0_qp)) < epsilon(0.0_qp))
+        call check(abs(trapz([1.0_qp], 1.0_qp)) < epsilon(0.0_qp))
+        call check(abs(trapz(a, a)) < epsilon(0.0_qp))
+        call check(abs(trapz([1.0_qp], [1.0_qp])) < epsilon(0.0_qp))
     end subroutine test_trapz_zero_qp
 
 end program test_trapz
diff --git a/src/tests/stats/test_mean.f90 b/src/tests/stats/test_mean.f90
index 4b6f62d1e..2721123d3 100644
--- a/src/tests/stats/test_mean.f90
+++ b/src/tests/stats/test_mean.f90
@@ -1,5 +1,5 @@
 program test_mean
-use stdlib_experimental_error, only: assert
+use stdlib_experimental_error, only: check
 use stdlib_experimental_kinds, only: sp, dp, int32, int64
 use stdlib_experimental_io, only: loadtxt
 use stdlib_experimental_stats, only: mean
@@ -25,65 +25,65 @@ program test_mean
 !sp
 call loadtxt("array3.dat", s)
 
-call assert( abs(mean(s) - sum(s)/real(size(s), sp)) < sptol)
-call assert( sum( abs( mean(s,1) - sum(s,1)/real(size(s,1), sp) )) < sptol)
-call assert( sum( abs( mean(s,2) - sum(s,2)/real(size(s,2), sp) )) < sptol)
+call check( abs(mean(s) - sum(s)/real(size(s), sp)) < sptol)
+call check( sum( abs( mean(s,1) - sum(s,1)/real(size(s,1), sp) )) < sptol)
+call check( sum( abs( mean(s,2) - sum(s,2)/real(size(s,2), sp) )) < sptol)
 
 ! check reduction of rank one array to scalar
-call assert(abs(mean(s1) - sum(s1) / real(size(s1), sp)) < sptol)
-call assert(abs(mean(s1, dim=1) - sum(s1, dim=1) / real(size(s1, dim=1), sp)) < sptol)
+call check(abs(mean(s1) - sum(s1) / real(size(s1), sp)) < sptol)
+call check(abs(mean(s1, dim=1) - sum(s1, dim=1) / real(size(s1, dim=1), sp)) < sptol)
 
 
 !dp
 call loadtxt("array3.dat", d)
 
-call assert( abs(mean(d) - sum(d)/real(size(d), dp)) < dptol)
-call assert( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) < dptol)
-call assert( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) < dptol)
+call check( abs(mean(d) - sum(d)/real(size(d), dp)) < dptol)
+call check( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) < dptol)
+call check( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) < dptol)
 
 !csp
 
 call loadtxt("array3.dat", d)
 cs = cmplx(1._sp, 1._sp)*d
 
-call assert( abs(mean(cs) - sum(cs)/real(size(cs), sp)) < sptol)
-call assert( sum( abs( mean(cs,1) - sum(cs,1)/real(size(cs,1), sp) )) < sptol)
-call assert( sum( abs( mean(cs,2) - sum(cs,2)/real(size(cs,2), sp) )) < sptol)
+call check( abs(mean(cs) - sum(cs)/real(size(cs), sp)) < sptol)
+call check( sum( abs( mean(cs,1) - sum(cs,1)/real(size(cs,1), sp) )) < sptol)
+call check( sum( abs( mean(cs,2) - sum(cs,2)/real(size(cs,2), sp) )) < sptol)
 
 !cdp
 
 call loadtxt("array3.dat", d)
 cd = cmplx(1._dp, 1._dp)*d
 
-call assert( abs(mean(cd) - sum(cd)/real(size(cd), dp)) < dptol)
-call assert( sum( abs( mean(cd,1) - sum(cd,1)/real(size(cd,1), dp) )) < dptol)
-call assert( sum( abs( mean(cd,2) - sum(cd,2)/real(size(cd,2), dp) )) < dptol)
+call check( abs(mean(cd) - sum(cd)/real(size(cd), dp)) < dptol)
+call check( sum( abs( mean(cd,1) - sum(cd,1)/real(size(cd,1), dp) )) < dptol)
+call check( sum( abs( mean(cd,2) - sum(cd,2)/real(size(cd,2), dp) )) < dptol)
 
 ! check mask = .false.
 
-call assert( ieee_is_nan(mean(d, .false.)))
-call assert( any(ieee_is_nan(mean(d, 1, .false.))))
-call assert( any(ieee_is_nan(mean(d, 2, .false.))))
+call check( ieee_is_nan(mean(d, .false.)))
+call check( any(ieee_is_nan(mean(d, 1, .false.))))
+call check( any(ieee_is_nan(mean(d, 2, .false.))))
 
 ! check mask of the same shape as input
-call assert( abs(mean(d, d > 0) - sum(d, d > 0)/real(count(d > 0), dp)) < dptol)
-call assert( sum(abs(mean(d, 1, d > 0) - sum(d, 1, d > 0)/real(count(d > 0, 1), dp))) < dptol)
-call assert( sum(abs(mean(d, 2, d > 0) - sum(d, 2, d > 0)/real(count(d > 0, 2), dp))) < dptol)
+call check( abs(mean(d, d > 0) - sum(d, d > 0)/real(count(d > 0), dp)) < dptol)
+call check( sum(abs(mean(d, 1, d > 0) - sum(d, 1, d > 0)/real(count(d > 0, 1), dp))) < dptol)
+call check( sum(abs(mean(d, 2, d > 0) - sum(d, 2, d > 0)/real(count(d > 0, 2), dp))) < dptol)
 
 !int32
 call loadtxt("array3.dat", d)
 
-call assert( abs(mean(int(d, int32)) - sum(real(int(d, int32),dp))/real(size(d), dp)) < dptol)
-call assert( sum(abs( mean(int(d, int32),1) - sum(real(int(d, int32),dp),1)/real(size(d,1), dp) )) < dptol)
-call assert( sum(abs( mean(int(d, int32),2) - sum(real(int(d, int32),dp),2)/real(size(d,2), dp) )) < dptol)
+call check( abs(mean(int(d, int32)) - sum(real(int(d, int32),dp))/real(size(d), dp)) < dptol)
+call check( sum(abs( mean(int(d, int32),1) - sum(real(int(d, int32),dp),1)/real(size(d,1), dp) )) < dptol)
+call check( sum(abs( mean(int(d, int32),2) - sum(real(int(d, int32),dp),2)/real(size(d,2), dp) )) < dptol)
 
 
 !int64
 call loadtxt("array3.dat", d)
 
-call assert( abs(mean(int(d, int64)) - sum(real(int(d, int64),dp))/real(size(d), dp)) < dptol)
-call assert( sum(abs( mean(int(d, int64),1) - sum(real(int(d, int64),dp),1)/real(size(d,1), dp) )) < dptol)
-call assert( sum(abs( mean(int(d, int64),2) - sum(real(int(d, int64),dp),2)/real(size(d,2), dp) )) < dptol)
+call check( abs(mean(int(d, int64)) - sum(real(int(d, int64),dp))/real(size(d), dp)) < dptol)
+call check( sum(abs( mean(int(d, int64),1) - sum(real(int(d, int64),dp),1)/real(size(d,1), dp) )) < dptol)
+call check( sum(abs( mean(int(d, int64),2) - sum(real(int(d, int64),dp),2)/real(size(d,2), dp) )) < dptol)
 
 
 !dp rank 3
@@ -92,10 +92,10 @@ program test_mean
 d3(:,:,2)=d*1.5;
 d3(:,:,3)=d*4;
 
-call assert( abs(mean(d3) - sum(d3)/real(size(d3), dp)) < dptol)
-call assert( sum( abs( mean(d3,1) - sum(d3,1)/real(size(d3,1), dp) )) < dptol)
-call assert( sum( abs( mean(d3,2) - sum(d3,2)/real(size(d3,2), dp) )) < dptol)
-call assert( sum( abs( mean(d3,3) - sum(d3,3)/real(size(d3,3), dp) )) < dptol)
+call check( abs(mean(d3) - sum(d3)/real(size(d3), dp)) < dptol)
+call check( sum( abs( mean(d3,1) - sum(d3,1)/real(size(d3,1), dp) )) < dptol)
+call check( sum( abs( mean(d3,2) - sum(d3,2)/real(size(d3,2), dp) )) < dptol)
+call check( sum( abs( mean(d3,3) - sum(d3,3)/real(size(d3,3), dp) )) < dptol)
 
 !cdp rank 3
 allocate(cd3(size(d,1),size(d,2),3))
@@ -104,10 +104,10 @@ program test_mean
 cd3(:,:,3)=d*4;
 cd3 = cmplx(1._sp, 1._sp)*cd3
 
-call assert( abs(mean(cd3) - sum(cd3)/real(size(cd3), dp)) < dptol)
-call assert( sum( abs( mean(cd3,1) - sum(cd3,1)/real(size(cd3,1), dp) )) < dptol)
-call assert( sum( abs( mean(cd3,2) - sum(cd3,2)/real(size(cd3,2), dp) )) < dptol)
-call assert( sum( abs( mean(cd3,3) - sum(cd3,3)/real(size(cd3,3), dp) )) < dptol)
+call check( abs(mean(cd3) - sum(cd3)/real(size(cd3), dp)) < dptol)
+call check( sum( abs( mean(cd3,1) - sum(cd3,1)/real(size(cd3,1), dp) )) < dptol)
+call check( sum( abs( mean(cd3,2) - sum(cd3,2)/real(size(cd3,2), dp) )) < dptol)
+call check( sum( abs( mean(cd3,3) - sum(cd3,3)/real(size(cd3,3), dp) )) < dptol)
 
 
 !dp rank 4
@@ -118,26 +118,26 @@ program test_mean
 d4(:,:,3,1)=d*4;
 d4(:,:,3,9)=d*4;
 
-call assert( abs(mean(d4) - sum(d4)/real(size(d4), dp)) < dptol)
-call assert( sum( abs( mean(d4,1) - sum(d4,1)/real(size(d4,1), dp) )) < dptol)
-call assert( sum( abs( mean(d4,2) - sum(d4,2)/real(size(d4,2), dp) )) < dptol)
-call assert( sum( abs( mean(d4,3) - sum(d4,3)/real(size(d4,3), dp) )) < dptol)
-call assert( sum( abs( mean(d4,4) - sum(d4,4)/real(size(d4,4), dp) )) < dptol)
+call check( abs(mean(d4) - sum(d4)/real(size(d4), dp)) < dptol)
+call check( sum( abs( mean(d4,1) - sum(d4,1)/real(size(d4,1), dp) )) < dptol)
+call check( sum( abs( mean(d4,2) - sum(d4,2)/real(size(d4,2), dp) )) < dptol)
+call check( sum( abs( mean(d4,3) - sum(d4,3)/real(size(d4,3), dp) )) < dptol)
+call check( sum( abs( mean(d4,4) - sum(d4,4)/real(size(d4,4), dp) )) < dptol)
 
 ! check mask = .false.
 
-call assert( ieee_is_nan(mean(d4, .false.)))
-call assert( any(ieee_is_nan(mean(d4, 1, .false.))))
-call assert( any(ieee_is_nan(mean(d4, 2, .false.))))
-call assert( any(ieee_is_nan(mean(d4, 3, .false.))))
-call assert( any(ieee_is_nan(mean(d4, 4, .false.))))
+call check( ieee_is_nan(mean(d4, .false.)))
+call check( any(ieee_is_nan(mean(d4, 1, .false.))))
+call check( any(ieee_is_nan(mean(d4, 2, .false.))))
+call check( any(ieee_is_nan(mean(d4, 3, .false.))))
+call check( any(ieee_is_nan(mean(d4, 4, .false.))))
 
 
 ! check mask of the same shape as input
-call assert( abs(mean(d4, d4 > 0) - sum(d4, d4 > 0)/real(count(d4 > 0), dp)) < dptol)
-call assert( any(ieee_is_nan(mean(d4, 1, d4 > 0))) )
-call assert( any(ieee_is_nan(mean(d4, 2, d4 > 0))) )
-call assert( any(ieee_is_nan(mean(d4, 3, d4 > 0))) )
-call assert( sum(abs(mean(d4, 4, d4 > 0) - sum(d4, 4, d4 > 0)/real(count(d4 > 0, 4), dp))) < dptol)
+call check( abs(mean(d4, d4 > 0) - sum(d4, d4 > 0)/real(count(d4 > 0), dp)) < dptol)
+call check( any(ieee_is_nan(mean(d4, 1, d4 > 0))) )
+call check( any(ieee_is_nan(mean(d4, 2, d4 > 0))) )
+call check( any(ieee_is_nan(mean(d4, 3, d4 > 0))) )
+call check( sum(abs(mean(d4, 4, d4 > 0) - sum(d4, 4, d4 > 0)/real(count(d4 > 0, 4), dp))) < dptol)
 
 end program
diff --git a/src/tests/stats/test_mean_f03.f90 b/src/tests/stats/test_mean_f03.f90
index c1d9cd2e1..cf05a4f39 100644
--- a/src/tests/stats/test_mean_f03.f90
+++ b/src/tests/stats/test_mean_f03.f90
@@ -1,5 +1,5 @@
 program test_mean
-use stdlib_experimental_error, only: assert
+use stdlib_experimental_error, only: check
 use stdlib_experimental_kinds, only: sp, dp, int32, int64
 use stdlib_experimental_io, only: loadtxt
 use stdlib_experimental_stats, only: mean
@@ -15,9 +15,9 @@ program test_mean
 !dp
 call loadtxt("array3.dat", d)
 
-call assert( mean(d) - sum(d)/real(size(d), dp) < dptol)
-call assert( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) < dptol)
-call assert( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) < dptol)
+call check( mean(d) - sum(d)/real(size(d), dp) < dptol)
+call check( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) < dptol)
+call check( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) < dptol)
 
 !dp rank 8
 allocate(d8(size(d,1), size(d,2), 3, 4, 5, 6, 7, 8))
@@ -25,31 +25,31 @@ program test_mean
 d8(:, :, 2, 4, 5 ,6 ,7 ,8)=d * 1.5_dp;
 d8(:, :, 3, 4, 5 ,6 ,7 ,8)=d * 4._dp;
 
-call assert( mean(d8) - sum(d8)/real(size(d8), dp) < dptol)
+call check( mean(d8) - sum(d8)/real(size(d8), dp) < dptol)
 
-call assert( sum( abs( mean(d8,1) - sum(d8,1)/real(size(d8,1), dp) )) < dptol)
-call assert( sum( abs( mean(d8,2) - sum(d8,2)/real(size(d8,2), dp) )) < dptol)
-call assert( sum( abs( mean(d8,3) - sum(d8,3)/real(size(d8,3), dp) )) < dptol)
-call assert( sum( abs( mean(d8,4) - sum(d8,4)/real(size(d8,4), dp) )) < dptol)
-call assert( sum( abs( mean(d8,5) - sum(d8,5)/real(size(d8,5), dp) )) < dptol)
-call assert( sum( abs( mean(d8,6) - sum(d8,6)/real(size(d8,6), dp) )) < dptol)
-call assert( sum( abs( mean(d8,7) - sum(d8,7)/real(size(d8,7), dp) )) < dptol)
-call assert( sum( abs( mean(d8,8) - sum(d8,8)/real(size(d8,8), dp) )) < dptol)
+call check( sum( abs( mean(d8,1) - sum(d8,1)/real(size(d8,1), dp) )) < dptol)
+call check( sum( abs( mean(d8,2) - sum(d8,2)/real(size(d8,2), dp) )) < dptol)
+call check( sum( abs( mean(d8,3) - sum(d8,3)/real(size(d8,3), dp) )) < dptol)
+call check( sum( abs( mean(d8,4) - sum(d8,4)/real(size(d8,4), dp) )) < dptol)
+call check( sum( abs( mean(d8,5) - sum(d8,5)/real(size(d8,5), dp) )) < dptol)
+call check( sum( abs( mean(d8,6) - sum(d8,6)/real(size(d8,6), dp) )) < dptol)
+call check( sum( abs( mean(d8,7) - sum(d8,7)/real(size(d8,7), dp) )) < dptol)
+call check( sum( abs( mean(d8,8) - sum(d8,8)/real(size(d8,8), dp) )) < dptol)
 
 !cdp rank 8
 allocate(cd8(size(d,1), size(d,2), 3, 4, 5, 6, 7, 8))
 cd8 = cmplx(1._dp, 1._dp, kind=dp)*d8
 
-call assert( abs(mean(cd8) - sum(cd8)/real(size(cd8), dp)) < dptol)
+call check( abs(mean(cd8) - sum(cd8)/real(size(cd8), dp)) < dptol)
 
-call assert( sum( abs( mean(cd8,1) - sum(cd8,1)/real(size(cd8,1), dp) )) < dptol)
-call assert( sum( abs( mean(cd8,2) - sum(cd8,2)/real(size(cd8,2), dp) )) < dptol)
-call assert( sum( abs( mean(cd8,3) - sum(cd8,3)/real(size(cd8,3), dp) )) < dptol)
-call assert( sum( abs( mean(cd8,4) - sum(cd8,4)/real(size(cd8,4), dp) )) < dptol)
-call assert( sum( abs( mean(cd8,5) - sum(cd8,5)/real(size(cd8,5), dp) )) < dptol)
-call assert( sum( abs( mean(cd8,6) - sum(cd8,6)/real(size(cd8,6), dp) )) < dptol)
-call assert( sum( abs( mean(cd8,7) - sum(cd8,7)/real(size(cd8,7), dp) )) < dptol)
-call assert( sum( abs( mean(cd8,8) - sum(cd8,8)/real(size(cd8,8), dp) )) < dptol)
+call check( sum( abs( mean(cd8,1) - sum(cd8,1)/real(size(cd8,1), dp) )) < dptol)
+call check( sum( abs( mean(cd8,2) - sum(cd8,2)/real(size(cd8,2), dp) )) < dptol)
+call check( sum( abs( mean(cd8,3) - sum(cd8,3)/real(size(cd8,3), dp) )) < dptol)
+call check( sum( abs( mean(cd8,4) - sum(cd8,4)/real(size(cd8,4), dp) )) < dptol)
+call check( sum( abs( mean(cd8,5) - sum(cd8,5)/real(size(cd8,5), dp) )) < dptol)
+call check( sum( abs( mean(cd8,6) - sum(cd8,6)/real(size(cd8,6), dp) )) < dptol)
+call check( sum( abs( mean(cd8,7) - sum(cd8,7)/real(size(cd8,7), dp) )) < dptol)
+call check( sum( abs( mean(cd8,8) - sum(cd8,8)/real(size(cd8,8), dp) )) < dptol)
 contains
 
 end program
diff --git a/src/tests/stats/test_moment.f90 b/src/tests/stats/test_moment.f90
index 6e7457994..8eea4529d 100644
--- a/src/tests/stats/test_moment.f90
+++ b/src/tests/stats/test_moment.f90
@@ -1,5 +1,5 @@
 program test_moment
-    use stdlib_experimental_error, only: assert
+    use stdlib_experimental_error, only: check
     use stdlib_experimental_kinds, only: sp, dp, int32, int64
     use stdlib_experimental_stats, only: moment
     use,intrinsic :: ieee_arithmetic, only : ieee_is_nan
@@ -45,32 +45,32 @@ subroutine test_sp(x1, x2)
 
         !1dim
         print*,' test_sp_1dim', order
-        call assert( abs(moment(x1, order)) < sptol)
-        call assert( abs(moment(x1, order, dim=1)) < sptol)
+        call check( abs(moment(x1, order)) < sptol)
+        call check( abs(moment(x1, order, dim=1)) < sptol)
 
         print*,' test_sp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_sp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5)) < sptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5)) < sptol)
+        call check( abs(moment(x1, order, x1 < 5)) < sptol)
+        call check( abs(moment(x1, order, 1, x1 < 5)) < sptol)
 
         !2dim
         print*,' test_sp_2dim', order
-        call assert( abs(moment(x2, order)) < sptol)
-        call assert( all( abs( moment(x2, order, 1)) < sptol))
-        call assert( all( abs( moment(x2, order, 2)) < sptol))
+        call check( abs(moment(x2, order)) < sptol)
+        call check( all( abs( moment(x2, order, 1)) < sptol))
+        call check( all( abs( moment(x2, order, 2)) < sptol))
 
         print*,' test_sp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_sp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)) < sptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11)) < sptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11)) < sptol))
+        call check( abs(moment(x2, order, x2 < 11)) < sptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11)) < sptol))
+        call check( all( abs( moment(x2, order, 2, x2 < 11)) < sptol))
 
         !3dim
         allocate(x3(size(x2,1),size(x2,2),3))
@@ -79,73 +79,73 @@ subroutine test_sp(x1, x2)
         x3(:,:,3)=x2*4;
     
         print*,' test_sp_3dim', order
-        call assert( abs(moment(x3, order)) < sptol)
-        call assert( all( abs( moment(x3, order, 1)) < sptol))
-        call assert( all( abs( moment(x3, order, 2)) < sptol))
-        call assert( all( abs( moment(x3, order, 3)) < sptol))
+        call check( abs(moment(x3, order)) < sptol)
+        call check( all( abs( moment(x3, order, 1)) < sptol))
+        call check( all( abs( moment(x3, order, 2)) < sptol))
+        call check( all( abs( moment(x3, order, 3)) < sptol))
     
         print*,' test_sp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_sp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) ) < sptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45)) < sptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45)) < sptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45)) < sptol ))
+        call check( abs(moment(x3, order, x3 < 11) ) < sptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45)) < sptol ))
+        call check( all( abs( moment(x3, order, 2, x3 < 45)) < sptol ))
+        call check( all( abs( moment(x3, order, 3, x3 < 45)) < sptol ))
  
 
         order = 2
 
         !1dim
         print*,' test_sp_1dim', order
-        call assert( abs(moment(x1, order) - 2._sp) < sptol)
-        call assert( abs(moment(x1, order, dim=1) - 2._sp) < sptol)
+        call check( abs(moment(x1, order) - 2._sp) < sptol)
+        call check( abs(moment(x1, order, dim=1) - 2._sp) < sptol)
 
         print*,' test_sp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_sp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5) - 1.25_sp) < sptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_sp) < sptol)
+        call check( abs(moment(x1, order, x1 < 5) - 1.25_sp) < sptol)
+        call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_sp) < sptol)
 
         !2dim
         print*,' test_sp_2dim', order
-        call assert( abs(moment(x2, order) - 107.25_sp/9.) < sptol)
-        call assert( all( abs( moment(x2, order, 1) - [5._sp, 5._sp, 1.25_sp]) < sptol))
-        call assert( all( abs( moment(x2, order, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]*2./3.) < sptol))
+        call check( abs(moment(x2, order) - 107.25_sp/9.) < sptol)
+        call check( all( abs( moment(x2, order, 1) - [5._sp, 5._sp, 1.25_sp]) < sptol))
+        call check( all( abs( moment(x2, order, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]*2./3.) < sptol))
 
         print*,' test_sp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_sp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)- 2.75_sp*3.) < sptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11) -&
+        call check( abs(moment(x2, order, x2 < 11)- 2.75_sp*3.) < sptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11) -&
                       [5._sp, 5._sp, 0.25_sp]) < sptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11) -&
+        call check( all( abs( moment(x2, order, 2, x2 < 11) -&
                       [19._sp*2./3., 43._sp/9.*2., 0.25_sp , 0.25_sp]) < sptol))
 
         !3dim
         print*,' test_sp_3dim', order
-        call assert( abs(moment(x3, order) - 153.4_sp*35./36.) < sptol)
-        call assert( all( abs( moment(x3, order, 1) -&
+        call check( abs(moment(x3, order) - 153.4_sp*35./36.) < sptol)
+        call check( all( abs( moment(x3, order, 1) -&
                  reshape([20._sp / 3., 20._sp / 3., 5._sp / 3.,&
                           4* 20._sp / 3., 4* 20._sp / 3., 4* 5._sp / 3.,&
                           16* 20._sp / 3., 16* 20._sp / 3., 16* 5._sp / 3.],&
                           [size(x3,2), size(x3,3)])*3._sp/4.)&
                  < sptol))
-        call assert( all( abs( moment(x3, order, 2) -&
+        call check( all( abs( moment(x3, order, 2) -&
                  reshape([19._sp, 43._sp / 3., 31._sp / 3. , 7.0_sp,&
                           4* 19.0_sp, 4* 43._sp / 3., 4* 31._sp / 3. , 4* 7.0_sp,&
                           16* 19.0_sp, 16* 43._sp / 3., 16* 31._sp / 3. , 16* 7.0_sp],&
                           [size(x3,1), size(x3,3)] )*2._sp/3.)&
                  < sptol))
-        call assert( all( abs( moment(x3, order, 3) -&
+        call check( all( abs( moment(x3, order, 3) -&
                  reshape([ 7._sp/3., 21._sp, 175._sp/3.,&
                            343._sp/3., 28._sp/3., 112._sp/3.,&
                            84._sp, 448._sp/3., 189._sp,&
@@ -154,24 +154,24 @@ subroutine test_sp(x1, x2)
                  < sptol))
     
         print*,' test_sp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_sp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < sptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45) -&
+        call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < sptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45) -&
                       reshape([5._sp, 5._sp, 1.25_sp,  20._sp, 20._sp, 5._sp,&
                                80._sp, 80._sp, 32._sp/3.],&
                                [size(x3, 2), size(x3, 3)])) < sptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 2, x3 < 45) -&
                       reshape([ 38._sp/3., 86._sp/9., 62._sp/9., 14._sp/3., 152._sp/3.,&
                                 344._sp/9., 248._sp/9., 168._sp/9., 1824._sp/9.,&
                                 1376._sp/9., 992._sp/9., 4._sp&
                                ],&
                       [size(x3, 1), size(x3, 3)])) < sptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 3, x3 < 45) -&
                      reshape([14._sp/9., 14._sp, 350._sp/9., 686._sp/9., 56._sp/9.,&
                               224._sp/9., 56._sp, 896._sp/9., 126._sp, 1400._sp/9.,&
                               1694._sp/9., 36._sp&
@@ -190,32 +190,32 @@ subroutine test_dp(x1, x2)
 
         !1dim
         print*,' test_dp_1dim', order
-        call assert( abs(moment(x1, order)) < dptol)
-        call assert( abs(moment(x1, order, dim=1)) < dptol)
+        call check( abs(moment(x1, order)) < dptol)
+        call check( abs(moment(x1, order, dim=1)) < dptol)
 
         print*,' test_dp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_dp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5)) < dptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5)) < dptol)
+        call check( abs(moment(x1, order, x1 < 5)) < dptol)
+        call check( abs(moment(x1, order, 1, x1 < 5)) < dptol)
 
         !2dim
         print*,' test_dp_2dim', order
-        call assert( abs(moment(x2, order)) < dptol)
-        call assert( all( abs( moment(x2, order, 1)) < dptol))
-        call assert( all( abs( moment(x2, order, 2)) < dptol))
+        call check( abs(moment(x2, order)) < dptol)
+        call check( all( abs( moment(x2, order, 1)) < dptol))
+        call check( all( abs( moment(x2, order, 2)) < dptol))
 
         print*,' test_dp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_dp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)) < dptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11)) < dptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11)) < dptol))
+        call check( abs(moment(x2, order, x2 < 11)) < dptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11)) < dptol))
+        call check( all( abs( moment(x2, order, 2, x2 < 11)) < dptol))
 
         !3dim
         allocate(x3(size(x2,1),size(x2,2),3))
@@ -224,74 +224,74 @@ subroutine test_dp(x1, x2)
         x3(:,:,3)=x2*4;
     
         print*,' test_dp_3dim', order
-        call assert( abs(moment(x3, order)) < dptol)
-        call assert( all( abs( moment(x3, order, 1)) < dptol))
-        call assert( all( abs( moment(x3, order, 2)) < dptol))
-        call assert( all( abs( moment(x3, order, 3)) < dptol))
+        call check( abs(moment(x3, order)) < dptol)
+        call check( all( abs( moment(x3, order, 1)) < dptol))
+        call check( all( abs( moment(x3, order, 2)) < dptol))
+        call check( all( abs( moment(x3, order, 3)) < dptol))
     
         print*,' test_dp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_dp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) ) < dptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45)) < dptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45)) < dptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45)) < dptol ))
+        call check( abs(moment(x3, order, x3 < 11) ) < dptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45)) < dptol ))
+        call check( all( abs( moment(x3, order, 2, x3 < 45)) < dptol ))
+        call check( all( abs( moment(x3, order, 3, x3 < 45)) < dptol ))
  
 
         order = 2
 
         !1dim
         print*,' test_dp_1dim', order
-        call assert( abs(moment(x1, order) - 2._dp) < dptol)
-        call assert( abs(moment(x1, order, dim=1) - 2._dp) < dptol)
+        call check( abs(moment(x1, order) - 2._dp) < dptol)
+        call check( abs(moment(x1, order, dim=1) - 2._dp) < dptol)
 
         print*,' test_dp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_dp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol)
+        call check( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol)
+        call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol)
 
         !2dim
         print*,' test_dp_2dim', order
-        call assert( abs(moment(x2, order) - 107.25_dp/9.) < dptol)
-        call assert( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol))
-        call assert( all( abs( moment(x2, order, 2) -&
+        call check( abs(moment(x2, order) - 107.25_dp/9.) < dptol)
+        call check( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol))
+        call check( all( abs( moment(x2, order, 2) -&
                       [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol))
 
         print*,' test_dp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_dp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11) -&
+        call check( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11) -&
                       [5._dp, 5._dp, 0.25_dp]) < dptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11) -&
+        call check( all( abs( moment(x2, order, 2, x2 < 11) -&
                       [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol))
 
         !3dim
         print*,' test_dp_3dim', order
-        call assert( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol)
-        call assert( all( abs( moment(x3, order, 1) -&
+        call check( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol)
+        call check( all( abs( moment(x3, order, 1) -&
                  reshape([20._dp / 3., 20._dp / 3., 5._dp / 3.,&
                           4* 20._dp / 3., 4* 20._dp / 3., 4* 5._dp / 3.,&
                           16* 20._dp / 3., 16* 20._dp / 3., 16* 5._dp / 3.],&
                           [size(x3,2), size(x3,3)])*3._dp/4.)&
                  < dptol))
-        call assert( all( abs( moment(x3, order, 2) -&
+        call check( all( abs( moment(x3, order, 2) -&
                  reshape([19._dp, 43._dp / 3., 31._dp / 3. , 7.0_dp,&
                           4* 19.0_dp, 4* 43._dp / 3., 4* 31._dp / 3. , 4* 7.0_dp,&
                           16* 19.0_dp, 16* 43._dp / 3., 16* 31._dp / 3. , 16* 7.0_dp],&
                           [size(x3,1), size(x3,3)] )*2._dp/3.)&
                  < dptol))
-        call assert( all( abs( moment(x3, order, 3) -&
+        call check( all( abs( moment(x3, order, 3) -&
                  reshape([ 7._dp/3., 21._dp, 175._dp/3.,&
                            343._dp/3., 28._dp/3., 112._dp/3.,&
                            84._dp, 448._dp/3., 189._dp,&
@@ -300,24 +300,24 @@ subroutine test_dp(x1, x2)
                  < dptol))
     
         print*,' test_dp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_dp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45) -&
+        call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45) -&
                       reshape([5._dp, 5._dp, 1.25_dp,  20._dp, 20._dp, 5._dp,&
                                80._dp, 80._dp, 32._dp/3.],&
                                [size(x3, 2), size(x3, 3)])) < dptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 2, x3 < 45) -&
                       reshape([ 38._dp/3., 86._dp/9., 62._dp/9., 14._dp/3., 152._dp/3.,&
                                 344._dp/9., 248._dp/9., 168._dp/9., 1824._dp/9.,&
                                 1376._dp/9., 992._dp/9., 4._dp&
                                ],&
                       [size(x3, 1), size(x3, 3)])) < dptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 3, x3 < 45) -&
                      reshape([14._dp/9., 14._dp, 350._dp/9., 686._dp/9., 56._dp/9.,&
                               224._dp/9., 56._dp, 896._dp/9., 126._dp, 1400._dp/9.,&
                               1694._dp/9., 36._dp&
@@ -336,32 +336,32 @@ subroutine test_int32(x1, x2)
 
         !1dim
         print*,' test_dp_1dim', order
-        call assert( abs(moment(x1, order)) < dptol)
-        call assert( abs(moment(x1, order, dim=1)) < dptol)
+        call check( abs(moment(x1, order)) < dptol)
+        call check( abs(moment(x1, order, dim=1)) < dptol)
 
         print*,' test_dp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_dp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5)) < dptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5)) < dptol)
+        call check( abs(moment(x1, order, x1 < 5)) < dptol)
+        call check( abs(moment(x1, order, 1, x1 < 5)) < dptol)
 
         !2dim
         print*,' test_dp_2dim', order
-        call assert( abs(moment(x2, order)) < dptol)
-        call assert( all( abs( moment(x2, order, 1)) < dptol))
-        call assert( all( abs( moment(x2, order, 2)) < dptol))
+        call check( abs(moment(x2, order)) < dptol)
+        call check( all( abs( moment(x2, order, 1)) < dptol))
+        call check( all( abs( moment(x2, order, 2)) < dptol))
 
         print*,' test_dp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_dp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)) < dptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11)) < dptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11)) < dptol))
+        call check( abs(moment(x2, order, x2 < 11)) < dptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11)) < dptol))
+        call check( all( abs( moment(x2, order, 2, x2 < 11)) < dptol))
 
         !3dim
         allocate(x3(size(x2,1),size(x2,2),3))
@@ -370,74 +370,74 @@ subroutine test_int32(x1, x2)
         x3(:,:,3)=x2*4;
     
         print*,' test_dp_3dim', order
-        call assert( abs(moment(x3, order)) < dptol)
-        call assert( all( abs( moment(x3, order, 1)) < dptol))
-        call assert( all( abs( moment(x3, order, 2)) < dptol))
-        call assert( all( abs( moment(x3, order, 3)) < dptol))
+        call check( abs(moment(x3, order)) < dptol)
+        call check( all( abs( moment(x3, order, 1)) < dptol))
+        call check( all( abs( moment(x3, order, 2)) < dptol))
+        call check( all( abs( moment(x3, order, 3)) < dptol))
     
         print*,' test_dp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_dp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) ) < dptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45)) < dptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45)) < dptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45)) < dptol ))
+        call check( abs(moment(x3, order, x3 < 11) ) < dptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45)) < dptol ))
+        call check( all( abs( moment(x3, order, 2, x3 < 45)) < dptol ))
+        call check( all( abs( moment(x3, order, 3, x3 < 45)) < dptol ))
  
 
         order = 2
 
         !1dim
         print*,' test_dp_1dim', order
-        call assert( abs(moment(x1, order) - 2._dp) < dptol)
-        call assert( abs(moment(x1, order, dim=1) - 2._dp) < dptol)
+        call check( abs(moment(x1, order) - 2._dp) < dptol)
+        call check( abs(moment(x1, order, dim=1) - 2._dp) < dptol)
 
         print*,' test_dp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_dp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol)
+        call check( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol)
+        call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol)
 
         !2dim
         print*,' test_dp_2dim', order
-        call assert( abs(moment(x2, order) - 107.25_dp/9.) < dptol)
-        call assert( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol))
-        call assert( all( abs( moment(x2, order, 2) -&
+        call check( abs(moment(x2, order) - 107.25_dp/9.) < dptol)
+        call check( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol))
+        call check( all( abs( moment(x2, order, 2) -&
                       [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol))
 
         print*,' test_dp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_dp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11) -&
+        call check( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11) -&
                       [5._dp, 5._dp, 0.25_dp]) < dptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11) -&
+        call check( all( abs( moment(x2, order, 2, x2 < 11) -&
                       [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol))
 
         !3dim
         print*,' test_dp_3dim', order
-        call assert( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol)
-        call assert( all( abs( moment(x3, order, 1) -&
+        call check( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol)
+        call check( all( abs( moment(x3, order, 1) -&
                  reshape([20._dp / 3., 20._dp / 3., 5._dp / 3.,&
                           4* 20._dp / 3., 4* 20._dp / 3., 4* 5._dp / 3.,&
                           16* 20._dp / 3., 16* 20._dp / 3., 16* 5._dp / 3.],&
                           [size(x3,2), size(x3,3)])*3._dp/4.)&
                  < dptol))
-        call assert( all( abs( moment(x3, order, 2) -&
+        call check( all( abs( moment(x3, order, 2) -&
                  reshape([19._dp, 43._dp / 3., 31._dp / 3. , 7.0_dp,&
                           4* 19.0_dp, 4* 43._dp / 3., 4* 31._dp / 3. , 4* 7.0_dp,&
                           16* 19.0_dp, 16* 43._dp / 3., 16* 31._dp / 3. , 16* 7.0_dp],&
                           [size(x3,1), size(x3,3)] )*2._dp/3.)&
                  < dptol))
-        call assert( all( abs( moment(x3, order, 3) -&
+        call check( all( abs( moment(x3, order, 3) -&
                  reshape([ 7._dp/3., 21._dp, 175._dp/3.,&
                            343._dp/3., 28._dp/3., 112._dp/3.,&
                            84._dp, 448._dp/3., 189._dp,&
@@ -446,24 +446,24 @@ subroutine test_int32(x1, x2)
                  < dptol))
     
         print*,' test_dp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_dp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45) -&
+        call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45) -&
                       reshape([5._dp, 5._dp, 1.25_dp,  20._dp, 20._dp, 5._dp,&
                                80._dp, 80._dp, 32._dp/3.],&
                                [size(x3, 2), size(x3, 3)])) < dptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 2, x3 < 45) -&
                       reshape([ 38._dp/3., 86._dp/9., 62._dp/9., 14._dp/3., 152._dp/3.,&
                                 344._dp/9., 248._dp/9., 168._dp/9., 1824._dp/9.,&
                                 1376._dp/9., 992._dp/9., 4._dp&
                                ],&
                       [size(x3, 1), size(x3, 3)])) < dptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 3, x3 < 45) -&
                      reshape([14._dp/9., 14._dp, 350._dp/9., 686._dp/9., 56._dp/9.,&
                               224._dp/9., 56._dp, 896._dp/9., 126._dp, 1400._dp/9.,&
                               1694._dp/9., 36._dp&
@@ -482,32 +482,32 @@ subroutine test_int64(x1, x2)
 
         !1dim
         print*,' test_dp_1dim', order
-        call assert( abs(moment(x1, order)) < dptol)
-        call assert( abs(moment(x1, order, dim=1)) < dptol)
+        call check( abs(moment(x1, order)) < dptol)
+        call check( abs(moment(x1, order, dim=1)) < dptol)
 
         print*,' test_dp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_dp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5)) < dptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5)) < dptol)
+        call check( abs(moment(x1, order, x1 < 5)) < dptol)
+        call check( abs(moment(x1, order, 1, x1 < 5)) < dptol)
 
         !2dim
         print*,' test_dp_2dim', order
-        call assert( abs(moment(x2, order)) < dptol)
-        call assert( all( abs( moment(x2, order, 1)) < dptol))
-        call assert( all( abs( moment(x2, order, 2)) < dptol))
+        call check( abs(moment(x2, order)) < dptol)
+        call check( all( abs( moment(x2, order, 1)) < dptol))
+        call check( all( abs( moment(x2, order, 2)) < dptol))
 
         print*,' test_dp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_dp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)) < dptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11)) < dptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11)) < dptol))
+        call check( abs(moment(x2, order, x2 < 11)) < dptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11)) < dptol))
+        call check( all( abs( moment(x2, order, 2, x2 < 11)) < dptol))
 
         !3dim
         allocate(x3(size(x2,1),size(x2,2),3))
@@ -516,74 +516,74 @@ subroutine test_int64(x1, x2)
         x3(:,:,3)=x2*4;
     
         print*,' test_dp_3dim', order
-        call assert( abs(moment(x3, order)) < dptol)
-        call assert( all( abs( moment(x3, order, 1)) < dptol))
-        call assert( all( abs( moment(x3, order, 2)) < dptol))
-        call assert( all( abs( moment(x3, order, 3)) < dptol))
+        call check( abs(moment(x3, order)) < dptol)
+        call check( all( abs( moment(x3, order, 1)) < dptol))
+        call check( all( abs( moment(x3, order, 2)) < dptol))
+        call check( all( abs( moment(x3, order, 3)) < dptol))
     
         print*,' test_dp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_dp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) ) < dptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45)) < dptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45)) < dptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45)) < dptol ))
+        call check( abs(moment(x3, order, x3 < 11) ) < dptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45)) < dptol ))
+        call check( all( abs( moment(x3, order, 2, x3 < 45)) < dptol ))
+        call check( all( abs( moment(x3, order, 3, x3 < 45)) < dptol ))
  
 
         order = 2
 
         !1dim
         print*,' test_dp_1dim', order
-        call assert( abs(moment(x1, order) - 2._dp) < dptol)
-        call assert( abs(moment(x1, order, dim=1) - 2._dp) < dptol)
+        call check( abs(moment(x1, order) - 2._dp) < dptol)
+        call check( abs(moment(x1, order, dim=1) - 2._dp) < dptol)
 
         print*,' test_dp_1dim_mask', order
-        call assert( ieee_is_nan(moment(x1, order, .false.)))
-        call assert( ieee_is_nan(moment(x1, order, 1, .false.)))
+        call check( ieee_is_nan(moment(x1, order, .false.)))
+        call check( ieee_is_nan(moment(x1, order, 1, .false.)))
 
         print*,' test_dp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol)
-        call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol)
+        call check( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol)
+        call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol)
 
         !2dim
         print*,' test_dp_2dim', order
-        call assert( abs(moment(x2, order) - 107.25_dp/9.) < dptol)
-        call assert( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol))
-        call assert( all( abs( moment(x2, order, 2) -&
+        call check( abs(moment(x2, order) - 107.25_dp/9.) < dptol)
+        call check( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol))
+        call check( all( abs( moment(x2, order, 2) -&
                       [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol))
 
         print*,' test_dp_2dim_mask', order
-        call assert( ieee_is_nan(moment(x2, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x2, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x2, order, 2, .false.))))
+        call check( ieee_is_nan(moment(x2, order, .false.)))
+        call check( any(ieee_is_nan(moment(x2, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x2, order, 2, .false.))))
 
         print*,' test_dp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol)
-        call assert( all( abs( moment(x2, order, 1, x2 < 11) -&
+        call check( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol)
+        call check( all( abs( moment(x2, order, 1, x2 < 11) -&
                       [5._dp, 5._dp, 0.25_dp]) < dptol))
-        call assert( all( abs( moment(x2, order, 2, x2 < 11) -&
+        call check( all( abs( moment(x2, order, 2, x2 < 11) -&
                       [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol))
 
         !3dim
         print*,' test_dp_3dim', order
-        call assert( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol)
-        call assert( all( abs( moment(x3, order, 1) -&
+        call check( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol)
+        call check( all( abs( moment(x3, order, 1) -&
                  reshape([20._dp / 3., 20._dp / 3., 5._dp / 3.,&
                           4* 20._dp / 3., 4* 20._dp / 3., 4* 5._dp / 3.,&
                           16* 20._dp / 3., 16* 20._dp / 3., 16* 5._dp / 3.],&
                           [size(x3,2), size(x3,3)])*3._dp/4.)&
                  < dptol))
-        call assert( all( abs( moment(x3, order, 2) -&
+        call check( all( abs( moment(x3, order, 2) -&
                  reshape([19._dp, 43._dp / 3., 31._dp / 3. , 7.0_dp,&
                           4* 19.0_dp, 4* 43._dp / 3., 4* 31._dp / 3. , 4* 7.0_dp,&
                           16* 19.0_dp, 16* 43._dp / 3., 16* 31._dp / 3. , 16* 7.0_dp],&
                           [size(x3,1), size(x3,3)] )*2._dp/3.)&
                  < dptol))
-        call assert( all( abs( moment(x3, order, 3) -&
+        call check( all( abs( moment(x3, order, 3) -&
                  reshape([ 7._dp/3., 21._dp, 175._dp/3.,&
                            343._dp/3., 28._dp/3., 112._dp/3.,&
                            84._dp, 448._dp/3., 189._dp,&
@@ -592,24 +592,24 @@ subroutine test_int64(x1, x2)
                  < dptol))
     
         print*,' test_dp_3dim_mask', order
-        call assert( ieee_is_nan(moment(x3, order, .false.)))
-        call assert( any(ieee_is_nan(moment(x3, order, 1, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 2, .false.))))
-        call assert( any(ieee_is_nan(moment(x3, order, 3, .false.))))
+        call check( ieee_is_nan(moment(x3, order, .false.)))
+        call check( any(ieee_is_nan(moment(x3, order, 1, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 2, .false.))))
+        call check( any(ieee_is_nan(moment(x3, order, 3, .false.))))
     
         print*,' test_dp_3dim_mask_array', order
-        call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol)
-        call assert( all( abs( moment(x3, order, 1, x3 < 45) -&
+        call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol)
+        call check( all( abs( moment(x3, order, 1, x3 < 45) -&
                       reshape([5._dp, 5._dp, 1.25_dp,  20._dp, 20._dp, 5._dp,&
                                80._dp, 80._dp, 32._dp/3.],&
                                [size(x3, 2), size(x3, 3)])) < dptol ))
-        call assert( all( abs( moment(x3, order, 2, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 2, x3 < 45) -&
                       reshape([ 38._dp/3., 86._dp/9., 62._dp/9., 14._dp/3., 152._dp/3.,&
                                 344._dp/9., 248._dp/9., 168._dp/9., 1824._dp/9.,&
                                 1376._dp/9., 992._dp/9., 4._dp&
                                ],&
                       [size(x3, 1), size(x3, 3)])) < dptol ))
-        call assert( all( abs( moment(x3, order, 3, x3 < 45) -&
+        call check( all( abs( moment(x3, order, 3, x3 < 45) -&
                      reshape([14._dp/9., 14._dp, 350._dp/9., 686._dp/9., 56._dp/9.,&
                               224._dp/9., 56._dp, 896._dp/9., 126._dp, 1400._dp/9.,&
                               1694._dp/9., 36._dp&
@@ -628,60 +628,60 @@ subroutine test_csp(x1, x2)
 
         !1dim
         print*,' test_sp_1dim', order
-        call assert( abs(moment(x1, order)) < sptol)
-        call assert( abs(moment(x1, order, dim=1)) < sptol)
+        call check( abs(moment(x1, order)) < sptol)
+        call check( abs(moment(x1, order, dim=1)) < sptol)
 
         print*,' test_sp_1dim_mask', order
-        call assert( ieee_is_nan(abs(moment(x1, order, .false.))))
-        call assert( ieee_is_nan(abs(moment(x1, order, 1, .false.))))
+        call check( ieee_is_nan(abs(moment(x1, order, .false.))))
+        call check( ieee_is_nan(abs(moment(x1, order, 1, .false.))))
 
         print*,' test_sp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, aimag(x1) == 0)) < sptol)
-        call assert( abs(moment(x1, order, 1, aimag(x1) == 0)) < sptol)
+        call check( abs(moment(x1, order, aimag(x1) == 0)) < sptol)
+        call check( abs(moment(x1, order, 1, aimag(x1) == 0)) < sptol)
 
         !2dim
         print*,' test_sp_2dim', order
-        call assert( abs(moment(x2, order)) < sptol)
-        call assert( all( abs( moment(x2, order, 1)) < sptol))
-        call assert( all( abs( moment(x2, order, 2)) < sptol))
+        call check( abs(moment(x2, order)) < sptol)
+        call check( all( abs( moment(x2, order, 1)) < sptol))
+        call check( all( abs( moment(x2, order, 2)) < sptol))
 
         print*,' test_sp_2dim_mask', order
-        call assert( ieee_is_nan(abs(moment(x2, order, .false.))))
-        call assert( any(ieee_is_nan(abs(moment(x2, order, 1, .false.)))))
-        call assert( any(ieee_is_nan(abs(moment(x2, order, 2, .false.)))))
+        call check( ieee_is_nan(abs(moment(x2, order, .false.))))
+        call check( any(ieee_is_nan(abs(moment(x2, order, 1, .false.)))))
+        call check( any(ieee_is_nan(abs(moment(x2, order, 2, .false.)))))
 
         print*,' test_sp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, aimag(x2) == 0)) < sptol)
-        call assert( all( abs( moment(x2, order, 1, aimag(x2) == 0)) < sptol))
-        call assert( any(ieee_is_nan( abs( moment(x2, order, 2, aimag(x2) == 0)))))
+        call check( abs(moment(x2, order, aimag(x2) == 0)) < sptol)
+        call check( all( abs( moment(x2, order, 1, aimag(x2) == 0)) < sptol))
+        call check( any(ieee_is_nan( abs( moment(x2, order, 2, aimag(x2) == 0)))))
 
         order = 2
 
         !1dim
         print*,' test_sp_1dim', order
-        call assert( abs(moment(x1, order) - (-6.459422410E-02,-0.556084037)) < sptol)
-        call assert( abs(moment(x1, order, dim=1) -&
+        call check( abs(moment(x1, order) - (-6.459422410E-02,-0.556084037)) < sptol)
+        call check( abs(moment(x1, order, dim=1) -&
                         (-6.459422410E-02,-0.556084037)) < sptol)
 
         print*,' test_sp_1dim_mask', order
-        call assert( ieee_is_nan(abs(moment(x1, order, .false.))))
-        call assert( ieee_is_nan(abs(moment(x1, order, 1, .false.))))
+        call check( ieee_is_nan(abs(moment(x1, order, .false.))))
+        call check( ieee_is_nan(abs(moment(x1, order, 1, .false.))))
 
         print*,' test_sp_1dim_mask_array', order
-        call assert( abs(moment(x1, order, aimag(x1) == 0) -&
+        call check( abs(moment(x1, order, aimag(x1) == 0) -&
                           (8.969944715E-02,0.00000000)) < sptol)
-        call assert( abs(moment(x1, order, 1, aimag(x1) == 0) -&
+        call check( abs(moment(x1, order, 1, aimag(x1) == 0) -&
                           (8.969944715E-02,0.00000000)) < sptol)
 
         !2dim
         print*,' test_sp_2dim', order
-        call assert( abs(moment(x2, order) - (-0.163121477,-1.86906016)) < sptol)
-        call assert( all( abs( moment(x2, order, 1) -&
+        call check( abs(moment(x2, order) - (-0.163121477,-1.86906016)) < sptol)
+        call check( all( abs( moment(x2, order, 1) -&
                      [(-6.459422410E-02,-0.556084037),&
                       (-0.581347823,-5.00475645),&
                       (-0.145336956,-1.25118911)]&
                      ) < sptol))
-        call assert( all( abs( moment(x2, order, 2) -&
+        call check( all( abs( moment(x2, order, 2) -&
                      [(0.240498722,0.00000000),&
                      (-1.49895227,0.00000000),&
                      (1.15390968,0.00000000),&
@@ -690,14 +690,14 @@ subroutine test_csp(x1, x2)
                      ) < sptol))
 
         print*,' test_sp_2dim_mask', order
-        call assert( ieee_is_nan(abs(moment(x2, order, .false.))))
-        call assert( any(ieee_is_nan(abs(moment(x2, order, 1, .false.)))))
-        call assert( any(ieee_is_nan(abs(moment(x2, order, 2, .false.)))))
+        call check( ieee_is_nan(abs(moment(x2, order, .false.))))
+        call check( any(ieee_is_nan(abs(moment(x2, order, 1, .false.)))))
+        call check( any(ieee_is_nan(abs(moment(x2, order, 2, .false.)))))
 
         print*,' test_sp_2dim_mask_array', order
-        call assert( abs(moment(x2, order, aimag(x2) == 0)-&
+        call check( abs(moment(x2, order, aimag(x2) == 0)-&
                      (1.08109438,0.00000000)) < sptol)
-        call assert( all( abs( moment(x2, order, 1, aimag(x2)==0) -&
+        call check( all( abs( moment(x2, order, 1, aimag(x2)==0) -&
                       [(8.969944715E-02,0.00000000),&
                        (0.807295084,0.00000000),&
                        (0.201823771,0.00000000)]&
diff --git a/src/tests/stats/test_var.f90 b/src/tests/stats/test_var.f90
index 36d7ac74b..892a718c4 100644
--- a/src/tests/stats/test_var.f90
+++ b/src/tests/stats/test_var.f90
@@ -1,5 +1,5 @@
 program test_var
-    use stdlib_experimental_error, only: assert
+    use stdlib_experimental_error, only: check
     use stdlib_experimental_kinds, only: sp, dp, int32, int64
     use stdlib_experimental_stats, only: var
     use,intrinsic :: ieee_arithmetic, only : ieee_is_nan
@@ -42,35 +42,35 @@ program test_var
     !sp
     !1dim
     print*,' test_sp_1dim'
-    call assert( abs(var(s1) - 2.5) < sptol)
-    call assert( abs(var(s1, dim=1) - 2.5) < sptol)
+    call check( abs(var(s1) - 2.5) < sptol)
+    call check( abs(var(s1, dim=1) - 2.5) < sptol)
 
     print*,' test_sp_1dim_mask'
-    call assert( ieee_is_nan(var(s1, .false.)))
-    call assert( ieee_is_nan(var(s1, 1, .false.)))
+    call check( ieee_is_nan(var(s1, .false.)))
+    call check( ieee_is_nan(var(s1, 1, .false.)))
 
     print*,' test_sp_1dim_mask_array'
-    call assert( abs(var(s1, s1 < 5) - 5./3.) < sptol)
-    call assert( ieee_is_nan((var(s1, s1 < 0.))))
-    call assert( ieee_is_nan((var(s1, s1 == 1.))))
-    call assert( abs(var(s1, 1, s1 < 5) - 5./3.) < sptol)
+    call check( abs(var(s1, s1 < 5) - 5./3.) < sptol)
+    call check( ieee_is_nan((var(s1, s1 < 0.))))
+    call check( ieee_is_nan((var(s1, s1 == 1.))))
+    call check( abs(var(s1, 1, s1 < 5) - 5./3.) < sptol)
 
     !2dim
     print*,' test_sp_2dim'
     s = d
-    call assert( abs(var(s) - 13) < sptol)
-    call assert( all( abs( var(s, 1) - [20. / 3., 20. / 3., 5. / 3.]) < sptol))
-    call assert( all( abs( var(s, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]) < sptol))
+    call check( abs(var(s) - 13) < sptol)
+    call check( all( abs( var(s, 1) - [20. / 3., 20. / 3., 5. / 3.]) < sptol))
+    call check( all( abs( var(s, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]) < sptol))
 
     print*,' test_sp_2dim_mask'
-    call assert( ieee_is_nan(var(s, .false.)))
-    call assert( any(ieee_is_nan(var(s, 1, .false.))))
-    call assert( any(ieee_is_nan(var(s, 2, .false.))))
+    call check( ieee_is_nan(var(s, .false.)))
+    call check( any(ieee_is_nan(var(s, 1, .false.))))
+    call check( any(ieee_is_nan(var(s, 2, .false.))))
 
     print*,' test_sp_2dim_mask_array'
-    call assert( abs(var(s, s < 11) - 27.5 / 3.) < sptol)
-    call assert( all( abs( var(s, 1, s < 11) - [20. / 3., 20. / 3., 0.5]) < sptol))
-    call assert( all( abs( var(s, 2, s < 11) - [19.0, 43. / 3., 0.5 , 0.5]) < sptol))
+    call check( abs(var(s, s < 11) - 27.5 / 3.) < sptol)
+    call check( all( abs( var(s, 1, s < 11) - [20. / 3., 20. / 3., 0.5]) < sptol))
+    call check( all( abs( var(s, 2, s < 11) - [19.0, 43. / 3., 0.5 , 0.5]) < sptol))
 
 
     !3dim
@@ -80,20 +80,20 @@ program test_var
     s3(:,:,3)=s*4;
 
     print*,' test_sp_3dim'
-    call assert( abs(var(s3) - 153.4) < sptol)
-    call assert( all( abs( var(s3, 1) -&
+    call check( abs(var(s3) - 153.4) < sptol)
+    call check( all( abs( var(s3, 1) -&
                  reshape([20. / 3., 20. / 3., 5. / 3.,&
                           4* 20. / 3., 4* 20. / 3., 4* 5. / 3.,&
                           16* 20. / 3., 16* 20. / 3., 16* 5. / 3.],&
                           [size(s3,2), size(s3,3)]))&
                  < sptol))
-    call assert( all( abs( var(s3, 2) -&
+    call check( all( abs( var(s3, 2) -&
                  reshape([19.0, 43. / 3., 31. / 3. , 7.0,&
                           4* 19.0, 4* 43. / 3., 4* 31. / 3. , 4* 7.0,&
                           16* 19.0, 16* 43. / 3., 16* 31. / 3. , 16* 7.0],&
                           [size(s3,1), size(s3,3)] ))&
                  < sptol))
-    call assert( all(abs( var(s3, 3) -&
+    call check( all(abs( var(s3, 3) -&
                  reshape([ 7./3., 21., 175./3.,&
                            343./3., 28./3., 112./3.,&
                            84., 448./3., 189.,&
@@ -101,19 +101,19 @@ program test_var
                  < sptol))
 
     print*,' test_sp_3dim_mask'
-    call assert( ieee_is_nan(var(s3, .false.)))
-    call assert( any(ieee_is_nan(var(s3, 1, .false.))))
-    call assert( any(ieee_is_nan(var(s3, 2, .false.))))
-    call assert( any(ieee_is_nan(var(s3, 3, .false.))))
+    call check( ieee_is_nan(var(s3, .false.)))
+    call check( any(ieee_is_nan(var(s3, 1, .false.))))
+    call check( any(ieee_is_nan(var(s3, 2, .false.))))
+    call check( any(ieee_is_nan(var(s3, 3, .false.))))
 
     print*,' test_sp_3dim_mask_array'
-    call assert( abs(var(s3, s3 < 11) - 8.2205877_sp) < sptol)
-    call assert( all( abs( var(s3, 1, s3 < 45) -&
+    call check( abs(var(s3, s3 < 11) - 8.2205877_sp) < sptol)
+    call check( all( abs( var(s3, 1, s3 < 45) -&
                   reshape([20./3., 20./3., 5./3.,  80./3., 80./3., 20./3.,&
                            320./3., 320./3., 16.],&
                   [size(s3, 2), size(s3, 3)])) < sptol ))
-    call assert( any( ieee_is_nan( var(s3, 2, s3 < 25))))
-    call assert( all( abs( var(s3, 3, s3 < 25) -&
+    call check( any( ieee_is_nan( var(s3, 2, s3 < 25))))
+    call check( all( abs( var(s3, 3, s3 < 25) -&
                  reshape([ 7./3., 21., 175./3.,&
                            24.5, 28./3., 112./3.,&
                            84., 32., 40.5,&
@@ -124,37 +124,37 @@ program test_var
     !dp
     !1dim
     print*,' test_dp_1dim'
-    call assert( abs(var(d1) - 2.5) < dptol)
-    call assert( abs(var(d1, 1) - 2.5) < dptol)
+    call check( abs(var(d1) - 2.5) < dptol)
+    call check( abs(var(d1, 1) - 2.5) < dptol)
 
     print*,' test_dp_1dim_mask'
-    call assert( ieee_is_nan(var(d1, .false.)))
-    call assert( ieee_is_nan(var(d1, 1, .false.)))
+    call check( ieee_is_nan(var(d1, .false.)))
+    call check( ieee_is_nan(var(d1, 1, .false.)))
 
     print*,' test_dp_1dim_mask_array'
-    call assert( abs(var(d1, d1 < 5) - 5._dp/3._dp) < dptol)
-    call assert( ieee_is_nan((var(d1, d1 < 0.))))
-    call assert( ieee_is_nan((var(d1, d1 == 1.))))
-    call assert( abs(var(d1, 1, d1 < 5) - 5._dp/3._dp) < dptol)
+    call check( abs(var(d1, d1 < 5) - 5._dp/3._dp) < dptol)
+    call check( ieee_is_nan((var(d1, d1 < 0.))))
+    call check( ieee_is_nan((var(d1, d1 == 1.))))
+    call check( abs(var(d1, 1, d1 < 5) - 5._dp/3._dp) < dptol)
 
     !2dim
     print*,' test_dp_2dim'
-    call assert( abs(var(d) - 13) < dptol)
-    call assert( all( abs( var(d,1) -&
+    call check( abs(var(d) - 13) < dptol)
+    call check( all( abs( var(d,1) -&
                  [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol))
-    call assert( all( abs( var(d,2) -&
+    call check( all( abs( var(d,2) -&
                  [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol))
 
     print*,' test_dp_2dim_mask'
-    call assert( ieee_is_nan(var(d, .false.)))
-    call assert( any(ieee_is_nan(var(d, 1, .false.))))
-    call assert( any(ieee_is_nan(var(d, 2, .false.))))
+    call check( ieee_is_nan(var(d, .false.)))
+    call check( any(ieee_is_nan(var(d, 1, .false.))))
+    call check( any(ieee_is_nan(var(d, 2, .false.))))
 
     print*,' test_dp_2dim_mask_array'
-    call assert( abs(var(d, d < 11) - 27.5_dp / 3._dp) < dptol)
-    call assert( all( abs( var(d, 1, d < 11) -&
+    call check( abs(var(d, d < 11) - 27.5_dp / 3._dp) < dptol)
+    call check( all( abs( var(d, 1, d < 11) -&
                  [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol))
-    call assert( all( abs( var(d, 2, d < 11) -&
+    call check( all( abs( var(d, 2, d < 11) -&
                  [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol))
 
     !3dim
@@ -164,15 +164,15 @@ program test_var
     d3(:,:,3)=d*4;
 
     print*,' test_dp_3dim'
-    call assert( abs(var(d3) - 153.4_dp) < dptol)
-    call assert( all( abs( var(d3, 1) -&
+    call check( abs(var(d3) - 153.4_dp) < dptol)
+    call check( all( abs( var(d3, 1) -&
                  reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,&
                           4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,&
                           16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],&
                           [size(d3,2), size(d3,3)]))&
                  < dptol))
     print*,' test_dp_3dim'
-    call assert( all( abs( var(d3, 2) -&
+    call check( all( abs( var(d3, 2) -&
                  reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,&
                           4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,&
                           16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,&
@@ -180,7 +180,7 @@ program test_var
                           [size(d3,1), size(d3,3)] ))&
                  < dptol))
     print*,' test_dp_3dim'
-    call assert( all(abs( var(d3, 3) -&
+    call check( all(abs( var(d3, 3) -&
                  reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                            343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,&
                            84._dp, 448._dp/3._dp, 189._dp,&
@@ -189,21 +189,21 @@ program test_var
                  < dptol))
 
     print*,' test_dp_3dim_mask'
-    call assert( ieee_is_nan(var(d3, .false.)))
-    call assert( any(ieee_is_nan(var(d3, 1, .false.))))
-    call assert( any(ieee_is_nan(var(d3, 2, .false.))))
-    call assert( any(ieee_is_nan(var(d3, 3, .false.))))
+    call check( ieee_is_nan(var(d3, .false.)))
+    call check( any(ieee_is_nan(var(d3, 1, .false.))))
+    call check( any(ieee_is_nan(var(d3, 2, .false.))))
+    call check( any(ieee_is_nan(var(d3, 3, .false.))))
 
     print*,' test_dp_3dim_mask_array'
-    call assert( abs(var(d3, d3 < 25) - 46.041379310344823_dp) < dptol)
-    call assert( all( abs( var(d3, 1, d3 < 45) -&
+    call check( abs(var(d3, d3 < 25) - 46.041379310344823_dp) < dptol)
+    call check( all( abs( var(d3, 1, d3 < 45) -&
                   reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,&
                            80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,&
                            320._dp/3._dp, 320._dp/3._dp, 16._dp],&
                            [size(d3, 2), size(d3, 3)]))&
                   < dptol ))
-    call assert( any( ieee_is_nan( var(d3, 2, d3 < 25))))
-    call assert( all( abs( var(d3, 3, d3 < 25) -&
+    call check( any( ieee_is_nan( var(d3, 2, d3 < 25))))
+    call check( all( abs( var(d3, 3, d3 < 25) -&
                  reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                            24.5_dp, 28._dp/3._dp, 112._dp/3._dp,&
                            84._dp, 32._dp, 40.5_dp,&
@@ -216,38 +216,38 @@ program test_var
     !int32
     !1dim
     print*,' test_int32_1dim'
-    call assert( abs(var(i321) - 2.5) < dptol)
-    call assert( abs(var(i321, 1) - 2.5) < dptol)
+    call check( abs(var(i321) - 2.5) < dptol)
+    call check( abs(var(i321, 1) - 2.5) < dptol)
 
     print*,' test_int32_1dim_mask'
-    call assert( ieee_is_nan(var(i321, .false.)))
-    call assert( ieee_is_nan(var(i321, 1, .false.)))
+    call check( ieee_is_nan(var(i321, .false.)))
+    call check( ieee_is_nan(var(i321, 1, .false.)))
 
     print*,' test_int32_1dim_mask_array'
-    call assert( abs(var(i321, i321 < 5) - 5._dp/3._dp) < dptol)
-    call assert( ieee_is_nan((var(i321, i321 < 0))))
-    call assert( ieee_is_nan((var(i321, i321 == 1))))
-    call assert( abs(var(i321, 1, i321 < 5) - 5._dp/3._dp) < dptol)
+    call check( abs(var(i321, i321 < 5) - 5._dp/3._dp) < dptol)
+    call check( ieee_is_nan((var(i321, i321 < 0))))
+    call check( ieee_is_nan((var(i321, i321 == 1))))
+    call check( abs(var(i321, 1, i321 < 5) - 5._dp/3._dp) < dptol)
 
     !2dim
     print*,' test_int32_2dim'
     i32 = d
-    call assert( abs(var(i32) - 13) < dptol)
-    call assert( all( abs( var(i32,1) -&
+    call check( abs(var(i32) - 13) < dptol)
+    call check( all( abs( var(i32,1) -&
                  [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol))
-    call assert( all( abs( var(i32,2) -&
+    call check( all( abs( var(i32,2) -&
                  [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol))
 
     print*,' test_int32_2dim_mask'
-    call assert( ieee_is_nan(var(i32, .false.)))
-    call assert( any(ieee_is_nan(var(i32, 1, .false.))))
-    call assert( any(ieee_is_nan(var(i32, 2, .false.))))
+    call check( ieee_is_nan(var(i32, .false.)))
+    call check( any(ieee_is_nan(var(i32, 1, .false.))))
+    call check( any(ieee_is_nan(var(i32, 2, .false.))))
 
     print*,' test_int32_2dim_mask_array'
-    call assert( abs(var(i32, i32 < 11) - 27.5_dp / 3._dp) < dptol)
-    call assert( all( abs( var(i32, 1, i32 < 11) -&
+    call check( abs(var(i32, i32 < 11) - 27.5_dp / 3._dp) < dptol)
+    call check( all( abs( var(i32, 1, i32 < 11) -&
                  [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol))
-    call assert( all( abs( var(i32, 2, i32 < 11) -&
+    call check( all( abs( var(i32, 2, i32 < 11) -&
                  [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol))
 
     !3dim
@@ -257,21 +257,21 @@ program test_var
     i323(:,:,3)=d*4;
 
     print*,' test_int32_3dim'
-    call assert( abs(var(i323) - 153.4_dp) < dptol)
-    call assert( all( abs( var(i323, 1) -&
+    call check( abs(var(i323) - 153.4_dp) < dptol)
+    call check( all( abs( var(i323, 1) -&
                  reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,&
                           4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,&
                           16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],&
                           [size(i323,2), size(i323,3)]))&
                  < dptol))
-    call assert( all( abs( var(i323, 2) -&
+    call check( all( abs( var(i323, 2) -&
                  reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,&
                           4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,&
                           16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,&
                           16* 7.0_dp],&
                           [size(i323,1), size(i323,3)] ))&
                  < dptol))
-    call assert( all(abs( var(i323, 3) -&
+    call check( all(abs( var(i323, 3) -&
                  reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                            343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,&
                            84._dp, 448._dp/3._dp, 189._dp,&
@@ -280,21 +280,21 @@ program test_var
                  < dptol))
 
     print*,' test_int32_3dim_mask'
-    call assert( ieee_is_nan(var(i323, .false.)))
-    call assert( any(ieee_is_nan(var(i323, 1, .false.))))
-    call assert( any(ieee_is_nan(var(i323, 2, .false.))))
-    call assert( any(ieee_is_nan(var(i323, 3, .false.))))
+    call check( ieee_is_nan(var(i323, .false.)))
+    call check( any(ieee_is_nan(var(i323, 1, .false.))))
+    call check( any(ieee_is_nan(var(i323, 2, .false.))))
+    call check( any(ieee_is_nan(var(i323, 3, .false.))))
 
     print*,' test_int32_3dim_mask_array'
-    call assert( abs(var(i323, i323 < 25) - 46.041379310344823_dp) < dptol)
-    call assert( all( abs( var(i323, 1, i323 < 45) -&
+    call check( abs(var(i323, i323 < 25) - 46.041379310344823_dp) < dptol)
+    call check( all( abs( var(i323, 1, i323 < 45) -&
                   reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,&
                            80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,&
                            320._dp/3._dp, 320._dp/3._dp, 16._dp],&
                            [size(i323, 2), size(i323, 3)]))&
                   < dptol ))
-    call assert( any( ieee_is_nan( var(i323, 2, i323 < 25))))
-    call assert( all( abs( var(i323, 3, i323 < 25) -&
+    call check( any( ieee_is_nan( var(i323, 2, i323 < 25))))
+    call check( all( abs( var(i323, 3, i323 < 25) -&
                  reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                            24.5_dp, 28._dp/3._dp, 112._dp/3._dp,&
                            84._dp, 32._dp, 40.5_dp,&
@@ -306,38 +306,38 @@ program test_var
     !int64
     !1dim
     print*,' test_int64_1dim'
-    call assert( abs(var(i641) - 2.5) < dptol)
-    call assert( abs(var(i641, 1) - 2.5) < dptol)
+    call check( abs(var(i641) - 2.5) < dptol)
+    call check( abs(var(i641, 1) - 2.5) < dptol)
 
     print*,' test_int64_1dim_mask'
-    call assert( ieee_is_nan(var(i641, .false.)))
-    call assert( ieee_is_nan(var(i641, 1, .false.)))
+    call check( ieee_is_nan(var(i641, .false.)))
+    call check( ieee_is_nan(var(i641, 1, .false.)))
 
     print*,' test_int641_1dim_mask_array'
-    call assert( abs(var(i641, i641 < 5) - 5._dp/3._dp) < dptol)
-    call assert( ieee_is_nan((var(i641, i641 < 0))))
-    call assert( ieee_is_nan((var(i641, i641 == 1))))
-    call assert( abs(var(i641, 1, i641 < 5) - 5._dp/3._dp) < dptol)
+    call check( abs(var(i641, i641 < 5) - 5._dp/3._dp) < dptol)
+    call check( ieee_is_nan((var(i641, i641 < 0))))
+    call check( ieee_is_nan((var(i641, i641 == 1))))
+    call check( abs(var(i641, 1, i641 < 5) - 5._dp/3._dp) < dptol)
 
     !2dim
     print*,' test_int64_2dim'
     i64 = d
-    call assert( abs(var(i64) - 13) < dptol)
-    call assert( all( abs( var(i64,1) -&
+    call check( abs(var(i64) - 13) < dptol)
+    call check( all( abs( var(i64,1) -&
                  [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol))
-    call assert( all( abs( var(i64,2) -&
+    call check( all( abs( var(i64,2) -&
                  [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol))
 
     print*,' test_int64_2dim_mask'
-    call assert( ieee_is_nan(var(i64, .false.)))
-    call assert( any(ieee_is_nan(var(i64, 1, .false.))))
-    call assert( any(ieee_is_nan(var(i64, 2, .false.))))
+    call check( ieee_is_nan(var(i64, .false.)))
+    call check( any(ieee_is_nan(var(i64, 1, .false.))))
+    call check( any(ieee_is_nan(var(i64, 2, .false.))))
 
     print*,' test_int64_2dim_mask_array'
-    call assert( abs(var(i64, i64 < 11) - 27.5_dp / 3._dp) < dptol)
-    call assert( all( abs( var(i64, 1, i64 < 11) -&
+    call check( abs(var(i64, i64 < 11) - 27.5_dp / 3._dp) < dptol)
+    call check( all( abs( var(i64, 1, i64 < 11) -&
                  [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol))
-    call assert( all( abs( var(i64, 2, i64 < 11) -&
+    call check( all( abs( var(i64, 2, i64 < 11) -&
                  [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol))
 
     !3dim
@@ -347,21 +347,21 @@ program test_var
     i643(:,:,3)=d*4;
 
     print*,' test_int32_3dim'
-    call assert( abs(var(i643) - 153.4_dp) < dptol)
-    call assert( all( abs( var(i643, 1) -&
+    call check( abs(var(i643) - 153.4_dp) < dptol)
+    call check( all( abs( var(i643, 1) -&
                  reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,&
                           4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,&
                           16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],&
                           [size(i643,2), size(i643,3)]))&
                  < dptol))
-    call assert( all( abs( var(i643, 2) -&
+    call check( all( abs( var(i643, 2) -&
                  reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,&
                           4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,&
                           16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,&
                           16* 7.0_dp],&
                           [size(i643,1), size(i643,3)] ))&
                  < dptol))
-    call assert( all(abs( var(i643, 3) -&
+    call check( all(abs( var(i643, 3) -&
                  reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                            343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,&
                            84._dp, 448._dp/3._dp, 189._dp,&
@@ -370,21 +370,21 @@ program test_var
                  < dptol))
 
     print*,' test_int32_3dim_mask'
-    call assert( ieee_is_nan(var(i643, .false.)))
-    call assert( any(ieee_is_nan(var(i643, 1, .false.))))
-    call assert( any(ieee_is_nan(var(i643, 2, .false.))))
-    call assert( any(ieee_is_nan(var(i643, 3, .false.))))
+    call check( ieee_is_nan(var(i643, .false.)))
+    call check( any(ieee_is_nan(var(i643, 1, .false.))))
+    call check( any(ieee_is_nan(var(i643, 2, .false.))))
+    call check( any(ieee_is_nan(var(i643, 3, .false.))))
 
     print*,' test_int64_3dim_mask_array'
-    call assert( abs(var(i643, i643 < 25) - 46.041379310344823_dp) < dptol)
-    call assert( all( abs( var(i643, 1, i643 < 45) -&
+    call check( abs(var(i643, i643 < 25) - 46.041379310344823_dp) < dptol)
+    call check( all( abs( var(i643, 1, i643 < 45) -&
                   reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,&
                            80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,&
                            320._dp/3._dp, 320._dp/3._dp, 16._dp],&
                            [size(i643, 2), size(i643, 3)]))&
                   < dptol ))
-    call assert( any( ieee_is_nan( var(i643, 2, i643 < 25))))
-    call assert( all( abs( var(i643, 3, i643 < 25) -&
+    call check( any( ieee_is_nan( var(i643, 2, i643 < 25))))
+    call check( all( abs( var(i643, 3, i643 < 25) -&
                  reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,&
                            24.5_dp, 28._dp/3._dp, 112._dp/3._dp,&
                            84._dp, 32._dp, 40.5_dp,&
@@ -395,16 +395,16 @@ program test_var
     !csp
     !1dim
     print*,' test_csp_1dim'
-    call assert( abs(var(cs1) - (var(real(cs1)) + var(aimag(cs1)))) < sptol)
-    call assert( abs(var(cs1, dim=1) - (var(real(cs1),1) + var(aimag(cs1), 1)) ) < sptol)
+    call check( abs(var(cs1) - (var(real(cs1)) + var(aimag(cs1)))) < sptol)
+    call check( abs(var(cs1, dim=1) - (var(real(cs1),1) + var(aimag(cs1), 1)) ) < sptol)
 
     print*,' test_csp_1dim_mask'
-    call assert( ieee_is_nan(var(cs1, .false.)))
-    call assert( ieee_is_nan(var(cs1, 1, .false.)))
+    call check( ieee_is_nan(var(cs1, .false.)))
+    call check( ieee_is_nan(var(cs1, 1, .false.)))
 
     print*,' test_csp_1dim_mask_array'
-    call assert( abs(var(cs1, aimag(cs1) == 0) - var(real(cs1), aimag(cs1) == 0)) < sptol)
-    call assert( abs(var(cs1, 1, aimag(cs1) == 0) - var(real(cs1), 1, aimag(cs1) == 0)) < sptol)
+    call check( abs(var(cs1, aimag(cs1) == 0) - var(real(cs1), aimag(cs1) == 0)) < sptol)
+    call check( abs(var(cs1, 1, aimag(cs1) == 0) - var(real(cs1), 1, aimag(cs1) == 0)) < sptol)
 
     !2dim
     cs(:,1) = cs1
@@ -412,33 +412,33 @@ program test_var
     cs(:,3) = cs1*1.5_sp
 
     print*,' test_csp_2dim'
-    call assert( abs(var(cs) - (var(real(cs)) + var(aimag(cs)))) < sptol)
-    call assert( all( abs( var(cs, 1) - (var(real(cs), 1) + var(aimag(cs), 1))) < sptol))
-    call assert( all( abs( var(cs, 2) - (var(real(cs), 2) + var(aimag(cs), 2))) < sptol))
+    call check( abs(var(cs) - (var(real(cs)) + var(aimag(cs)))) < sptol)
+    call check( all( abs( var(cs, 1) - (var(real(cs), 1) + var(aimag(cs), 1))) < sptol))
+    call check( all( abs( var(cs, 2) - (var(real(cs), 2) + var(aimag(cs), 2))) < sptol))
 
     print*,' test_csp_2dim_mask'
-    call assert( ieee_is_nan(var(cs, .false.)))
-    call assert( any(ieee_is_nan(var(cs, 1, .false.))))
-    call assert( any(ieee_is_nan(var(cs, 2, .false.))))
+    call check( ieee_is_nan(var(cs, .false.)))
+    call check( any(ieee_is_nan(var(cs, 1, .false.))))
+    call check( any(ieee_is_nan(var(cs, 2, .false.))))
 
     print*,' test_csp_2dim_mask_array'
-    call assert( abs(var(cs, aimag(cs) == 0) - var(real(cs), aimag(cs) == 0)) < sptol)
-    call assert( all( abs( var(cs, 1, aimag(cs) == 0) - var(real(cs), 1, aimag(cs) == 0)) < sptol))
-    call assert( any( ieee_is_nan( var(cs, 2, aimag(cs) == 0))))
+    call check( abs(var(cs, aimag(cs) == 0) - var(real(cs), aimag(cs) == 0)) < sptol)
+    call check( all( abs( var(cs, 1, aimag(cs) == 0) - var(real(cs), 1, aimag(cs) == 0)) < sptol))
+    call check( any( ieee_is_nan( var(cs, 2, aimag(cs) == 0))))
 
     !cdp
     !1dim
     print*,' test_cdp_1dim'
-    call assert( abs(var(cd1) - (var(real(cd1)) + var(aimag(cd1)))) < dptol)
-    call assert( abs(var(cd1, dim=1) - (var(real(cd1),1) + var(aimag(cd1), 1)) ) < dptol)
+    call check( abs(var(cd1) - (var(real(cd1)) + var(aimag(cd1)))) < dptol)
+    call check( abs(var(cd1, dim=1) - (var(real(cd1),1) + var(aimag(cd1), 1)) ) < dptol)
 
     print*,' test_cdp_1dim_mask'
-    call assert( ieee_is_nan(var(cd1, .false.)))
-    call assert( ieee_is_nan(var(cd1, 1, .false.)))
+    call check( ieee_is_nan(var(cd1, .false.)))
+    call check( ieee_is_nan(var(cd1, 1, .false.)))
 
     print*,' test_cdp_1dim_mask_array'
-    call assert( abs(var(cd1, aimag(cd1) == 0) - var(real(cd1), aimag(cd1) == 0)) < dptol)
-    call assert( abs(var(cd1, 1, aimag(cd1) == 0) - var(real(cd1), 1, aimag(cd1) == 0)) < dptol)
+    call check( abs(var(cd1, aimag(cd1) == 0) - var(real(cd1), aimag(cd1) == 0)) < dptol)
+    call check( abs(var(cd1, 1, aimag(cd1) == 0) - var(real(cd1), 1, aimag(cd1) == 0)) < dptol)
 
     !2dim
     cd(:,1) = cd1
@@ -446,18 +446,18 @@ program test_var
     cd(:,3) = cd1*1.5_sp
 
     print*,' test_cdp_2dim'
-    call assert( abs(var(cd) - (var(real(cd)) + var(aimag(cd)))) < dptol)
-    call assert( all( abs( var(cd, 1) - (var(real(cd), 1) + var(aimag(cd), 1))) < dptol))
-    call assert( all( abs( var(cd, 2) - (var(real(cd), 2) + var(aimag(cd), 2))) < dptol))
+    call check( abs(var(cd) - (var(real(cd)) + var(aimag(cd)))) < dptol)
+    call check( all( abs( var(cd, 1) - (var(real(cd), 1) + var(aimag(cd), 1))) < dptol))
+    call check( all( abs( var(cd, 2) - (var(real(cd), 2) + var(aimag(cd), 2))) < dptol))
 
     print*,' test_cdp_2dim_mask'
-    call assert( ieee_is_nan(var(cd, .false.)))
-    call assert( any(ieee_is_nan(var(cd, 1, .false.))))
-    call assert( any(ieee_is_nan(var(cd, 2, .false.))))
+    call check( ieee_is_nan(var(cd, .false.)))
+    call check( any(ieee_is_nan(var(cd, 1, .false.))))
+    call check( any(ieee_is_nan(var(cd, 2, .false.))))
 
     print*,' test_cdp_2dim_mask_array'
-    call assert( abs(var(cd, aimag(cd) == 0) - var(real(cd), aimag(cd) == 0)) < dptol)
-    call assert( all( abs( var(cd, 1, aimag(cd) == 0) - var(real(cd), 1, aimag(cd) == 0)) < dptol))
-    call assert( any( ieee_is_nan( var(cd, 2, aimag(cd) == 0))))
+    call check( abs(var(cd, aimag(cd) == 0) - var(real(cd), aimag(cd) == 0)) < dptol)
+    call check( all( abs( var(cd, 1, aimag(cd) == 0) - var(real(cd), 1, aimag(cd) == 0)) < dptol))
+    call check( any( ieee_is_nan( var(cd, 2, aimag(cd) == 0))))
 
 end program
diff --git a/src/tests/stats/test_varn.f90 b/src/tests/stats/test_varn.f90
index 9e33b4e67..af731e704 100644
--- a/src/tests/stats/test_varn.f90
+++ b/src/tests/stats/test_varn.f90
@@ -1,5 +1,5 @@
 program test_varn
-    use stdlib_experimental_error, only: assert
+    use stdlib_experimental_error, only: check
     use stdlib_experimental_kinds, only: sp, dp, int32, int64
     use stdlib_experimental_stats, only: var
     use,intrinsic :: ieee_arithmetic, only : ieee_is_nan
@@ -34,37 +34,37 @@ program test_varn
     !sp
     !1dim
     print*,' test_sp_1dim'
-    call assert( abs(var(s1, corrected=.false.) - 2.5*(4./5.)) < sptol)
-    call assert( abs(var(s1, dim=1, corrected=.false.) - 2.5*(4./5.)) < sptol)
+    call check( abs(var(s1, corrected=.false.) - 2.5*(4./5.)) < sptol)
+    call check( abs(var(s1, dim=1, corrected=.false.) - 2.5*(4./5.)) < sptol)
 
     print*,' test_sp_1dim_mask'
-    call assert( ieee_is_nan(var(s1, .false., corrected=.false.)))
-    call assert( ieee_is_nan(var(s1, 1, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(s1, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(s1, 1, .false., corrected=.false.)))
 
     print*,' test_sp_1dim_mask_array'
-    call assert( abs(var(s1, s1 < 5, corrected=.false.) - 5./4.) < sptol)
-    call assert( ieee_is_nan((var(s1, s1 < 0., corrected=.false.))))
-    call assert( abs(var(s1, s1 == 1., corrected=.false.)) < sptol)
-    call assert( abs(var(s1, 1, s1 < 5, corrected=.false.) - 5./4.) < sptol)
+    call check( abs(var(s1, s1 < 5, corrected=.false.) - 5./4.) < sptol)
+    call check( ieee_is_nan((var(s1, s1 < 0., corrected=.false.))))
+    call check( abs(var(s1, s1 == 1., corrected=.false.)) < sptol)
+    call check( abs(var(s1, 1, s1 < 5, corrected=.false.) - 5./4.) < sptol)
 
     !2dim
     print*,' test_sp_2dim'
     s = d
-    call assert( abs(var(s, corrected=.false.) - 13.*11./12.) < sptol)
-    call assert( all( abs( var(s, 1, corrected=.false.) - [20., 20., 5.]/4.) < sptol))
-    call assert( all( abs( var(s, 2, corrected=.false.) -&
+    call check( abs(var(s, corrected=.false.) - 13.*11./12.) < sptol)
+    call check( all( abs( var(s, 1, corrected=.false.) - [20., 20., 5.]/4.) < sptol))
+    call check( all( abs( var(s, 2, corrected=.false.) -&
                       [19.0, 43. / 3., 31. / 3. , 7.0]*2./3.) < sptol))
 
     print*,' test_sp_2dim_mask'
-    call assert( ieee_is_nan(var(s, .false., corrected=.false.)))
-    call assert( any(ieee_is_nan(var(s, 1, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(s, 2, .false., corrected=.false.))))
+    call check( ieee_is_nan(var(s, .false., corrected=.false.)))
+    call check( any(ieee_is_nan(var(s, 1, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(s, 2, .false., corrected=.false.))))
 
     print*,' test_sp_2dim_mask_array'
-    call assert( abs(var(s, s < 11, corrected=.false.) - 2.75*3.) < sptol)
-    call assert( all( abs( var(s, 1, s < 11, corrected=.false.) -&
+    call check( abs(var(s, s < 11, corrected=.false.) - 2.75*3.) < sptol)
+    call check( all( abs( var(s, 1, s < 11, corrected=.false.) -&
                       [5., 5., 0.25]) < sptol))
-    call assert( all( abs( var(s, 2, s < 11, corrected=.false.) -&
+    call check( all( abs( var(s, 2, s < 11, corrected=.false.) -&
                       [19.0*2./3., 43./9.*2., 0.25 , 0.25]) < sptol))
 
 
@@ -75,20 +75,20 @@ program test_varn
     s3(:,:,3)=s*4;
 
     print*,' test_sp_3dim'
-    call assert( abs(var(s3, corrected=.false.) - 153.4*35./36.) < sptol)
-    call assert( all( abs( var(s3, 1, corrected=.false.) -&
+    call check( abs(var(s3, corrected=.false.) - 153.4*35./36.) < sptol)
+    call check( all( abs( var(s3, 1, corrected=.false.) -&
                  reshape([20. / 3., 20. / 3., 5. / 3.,&
                           4* 20. / 3., 4* 20. / 3., 4* 5. / 3.,&
                           16* 20. / 3., 16* 20. / 3., 16* 5. / 3.],&
                           [size(s3,2), size(s3,3)])*3./4.)&
                  < sptol))
-    call assert( all( abs( var(s3, 2, corrected=.false.) -&
+    call check( all( abs( var(s3, 2, corrected=.false.) -&
                  reshape([19.0, 43. / 3., 31. / 3. , 7.0,&
                           4* 19.0, 4* 43. / 3., 4* 31. / 3. , 4* 7.0,&
                           16* 19.0, 16* 43. / 3., 16* 31. / 3. , 16* 7.0],&
                           [size(s3,1), size(s3,3)] )*2./3.)&
                  < sptol))
-    call assert( all(abs( var(s3, 3, corrected=.false.) -&
+    call check( all(abs( var(s3, 3, corrected=.false.) -&
                  reshape([ 7./3., 21., 175./3.,&
                            343./3., 28./3., 112./3.,&
                            84., 448./3., 189.,&
@@ -96,23 +96,23 @@ program test_varn
                  < sptol))
 
     print*,' test_sp_3dim_mask'
-    call assert( ieee_is_nan(var(s3, .false., corrected=.false.)))
-    call assert( any(ieee_is_nan(var(s3, 1, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(s3, 2, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(s3, 3, .false., corrected=.false.))))
+    call check( ieee_is_nan(var(s3, .false., corrected=.false.)))
+    call check( any(ieee_is_nan(var(s3, 1, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(s3, 2, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(s3, 3, .false., corrected=.false.))))
 
     print*,' test_sp_3dim_mask_array'
-    call assert( abs(var(s3, s3 < 11, corrected=.false.) - 7.73702383_sp) < sptol)
-    call assert( all( abs( var(s3, 1, s3 < 45, corrected=.false.) -&
+    call check( abs(var(s3, s3 < 11, corrected=.false.) - 7.73702383_sp) < sptol)
+    call check( all( abs( var(s3, 1, s3 < 45, corrected=.false.) -&
                   reshape([5., 5., 1.25,  20., 20., 5., 80., 80., 32./3.],&
                   [size(s3, 2), size(s3, 3)])) < sptol ))
-    call assert( all( abs( var(s3, 2, s3 < 45, corrected=.false.) -&
+    call check( all( abs( var(s3, 2, s3 < 45, corrected=.false.) -&
                   reshape([ 38./3., 86./9., 6.88888931, 14./3., 152./3.,&
                             38.2222214, 27.5555573, 18.6666660, 202.666672,&
                             152.888885, 110.222229, 4.&
                            ],&
                   [size(s3, 1), size(s3, 3)])) < sptol ))
-    call assert( all( abs( var(s3, 3, s3 < 45, corrected=.false.) -&
+    call check( all( abs( var(s3, 3, s3 < 45, corrected=.false.) -&
                  reshape([1.555555, 14., 38.888888, 76.222222, 6.2222222,&
                           24.888888, 56., 99.5555, 126., 155.555555, 188.22222, 36.&
                            ], [size(s3,1), size(s3,2)] ))&
@@ -121,36 +121,36 @@ program test_varn
     !dp
     !1dim
     print*,' test_dp_1dim'
-    call assert( abs(var(d1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
-    call assert( abs(var(d1, dim=1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
+    call check( abs(var(d1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
+    call check( abs(var(d1, dim=1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
 
     print*,' test_dp_1dim_mask'
-    call assert( ieee_is_nan(var(d1, .false., corrected=.false.)))
-    call assert( ieee_is_nan(var(d1, 1, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(d1, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(d1, 1, .false., corrected=.false.)))
 
     print*,' test_dp_1dim_mask_array'
-    call assert( abs(var(d1, d1 < 5, corrected=.false.) - 5._dp/4.) < dptol)
-    call assert( ieee_is_nan((var(d1, d1 < 0, corrected=.false.))))
-    call assert( abs(var(d1, d1 == 1, corrected=.false.)) < dptol)
-    call assert( abs(var(d1, 1, d1 < 5, corrected=.false.) - 5._dp/4.) < dptol)
+    call check( abs(var(d1, d1 < 5, corrected=.false.) - 5._dp/4.) < dptol)
+    call check( ieee_is_nan((var(d1, d1 < 0, corrected=.false.))))
+    call check( abs(var(d1, d1 == 1, corrected=.false.)) < dptol)
+    call check( abs(var(d1, 1, d1 < 5, corrected=.false.) - 5._dp/4.) < dptol)
 
     !2dim
     print*,' test_dp_2dim'
-    call assert( abs(var(d, corrected=.false.) - 13._dp*11./12.) < dptol)
-    call assert( all( abs( var(d, 1, corrected=.false.) - [20., 20., 5.]/4._dp) < dptol))
-    call assert( all( abs( var(d, 2, corrected=.false.) -&
+    call check( abs(var(d, corrected=.false.) - 13._dp*11./12.) < dptol)
+    call check( all( abs( var(d, 1, corrected=.false.) - [20., 20., 5.]/4._dp) < dptol))
+    call check( all( abs( var(d, 2, corrected=.false.) -&
                       [38._dp, 86._dp / 3._dp, 62._dp / 3._dp , 14._dp]/3._dp) < dptol))
 
     print*,' test_dp_2dim_mask'
-    call assert( ieee_is_nan(var(d, .false., corrected=.false.)))
-    call assert( any(ieee_is_nan(var(d, 1, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(d, 2, .false., corrected=.false.))))
+    call check( ieee_is_nan(var(d, .false., corrected=.false.)))
+    call check( any(ieee_is_nan(var(d, 1, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(d, 2, .false., corrected=.false.))))
 
     print*,' test_dp_2dim_mask_array'
-    call assert( abs(var(d, d < 11, corrected=.false.) - 2.75_dp*3._dp) < dptol)
-    call assert( all( abs( var(d, 1, d < 11, corrected=.false.) -&
+    call check( abs(var(d, d < 11, corrected=.false.) - 2.75_dp*3._dp) < dptol)
+    call check( all( abs( var(d, 1, d < 11, corrected=.false.) -&
                       [5._dp, 5._dp, 0.25_dp]) < dptol))
-    call assert( all( abs( var(d, 2, d < 11, corrected=.false.) -&
+    call check( all( abs( var(d, 2, d < 11, corrected=.false.) -&
                       [38._dp/3., 86._dp/9., 0.25_dp , 0.25_dp]) < dptol))
 
 
@@ -161,21 +161,21 @@ program test_varn
     d3(:,:,3)=d*4;
 
     print*,' test_dp_3dim'
-    call assert( abs(var(d3, corrected=.false.) - 153.4_dp*35._dp/36._dp) < dptol)
-    call assert( all( abs( var(d3, 1, corrected=.false.) -&
+    call check( abs(var(d3, corrected=.false.) - 153.4_dp*35._dp/36._dp) < dptol)
+    call check( all( abs( var(d3, 1, corrected=.false.) -&
                  reshape([20._dp , 20._dp , 5._dp ,&
                           4* 20._dp , 4* 20._dp , 4* 5._dp ,&
                           16* 20._dp , 16* 20._dp , 16* 5._dp ],&
                           [size(d3,2), size(d3,3)])/4._dp)&
                  < dptol))
-    call assert( all( abs( var(d3, 2, corrected=.false.) -&
+    call check( all( abs( var(d3, 2, corrected=.false.) -&
                  reshape([38._dp, 86. / 3._dp, 62. / 3._dp , 14._dp,&
                           8* 19._dp, 8* 43. / 3._dp, 8* 31. / 3._dp, 8* 7._dp,&
                           32* 19._dp, 32* 43. / 3._dp, 32* 31. / 3._dp, 32* 7._dp],&
                           [size(d3,1), size(d3,3)] )/3._dp)&
                  < dptol))
     print*,' test_dp_3dim'
-    call assert( all(abs( var(d3, 3, corrected=.false.) -&
+    call check( all(abs( var(d3, 3, corrected=.false.) -&
                  reshape([ 7._dp/3., 21._dp, 175._dp/3.,&
                            343._dp/3., 28._dp/3., 112._dp/3.,&
                            84._dp, 448._dp/3., 189._dp,&
@@ -184,26 +184,26 @@ program test_varn
                  < dptol))
 
     print*,' test_dp_3dim_mask'
-    call assert( ieee_is_nan(var(d3, .false., corrected=.false.)))
-    call assert( any(ieee_is_nan(var(d3, 1, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(d3, 2, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(d3, 3, .false., corrected=.false.))))
+    call check( ieee_is_nan(var(d3, .false., corrected=.false.)))
+    call check( any(ieee_is_nan(var(d3, 1, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(d3, 2, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(d3, 3, .false., corrected=.false.))))
 
     print*,' test_dp_3dim_mask_array'
-    call assert( abs(var(d3, d3 < 11, corrected=.false.) -&
+    call check( abs(var(d3, d3 < 11, corrected=.false.) -&
                      7.7370242214532876_dp) < dptol)
-    call assert( all( abs( var(d3, 1, d3 < 45, corrected=.false.) -&
+    call check( all( abs( var(d3, 1, d3 < 45, corrected=.false.) -&
                   reshape([5._dp, 5._dp, 1.25_dp,  20._dp, 20._dp, 5._dp,&
                            80._dp, 80._dp, 32._dp/3.],&
                            [size(d3, 2), size(d3, 3)])) < dptol ))
-    call assert( all( abs( var(d3, 2, d3 < 45, corrected=.false.) -&
+    call check( all( abs( var(d3, 2, d3 < 45, corrected=.false.) -&
                   reshape([ 38._dp/3., 86._dp/9., 62._dp/9., 14._dp/3., 152._dp/3.,&
                             344._dp/9., 248._dp/9., 168._dp/9., 1824._dp/9.,&
                             1376._dp/9., 992._dp/9., 4._dp&
                            ],&
                   [size(d3, 1), size(d3, 3)])) < dptol ))
     print*,' test_dp_3dim_mask_array'
-    call assert( all( abs( var(d3, 3, d3 < 45, corrected=.false.) -&
+    call check( all( abs( var(d3, 3, d3 < 45, corrected=.false.) -&
                  reshape([14._dp/9., 14._dp, 350._dp/9., 686._dp/9., 56._dp/9.,&
                           224._dp/9., 56._dp, 896._dp/9., 126._dp, 1400._dp/9.,&
                           1694._dp/9., 36._dp&
@@ -213,38 +213,38 @@ program test_varn
     !int32
     !1dim
     print*,' test_int32_1dim'
-    call assert( abs(var(i321, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
-    call assert( abs(var(i321, dim=1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
+    call check( abs(var(i321, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
+    call check( abs(var(i321, dim=1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol)
 
     print*,' test_int32_1dim_mask'
-    call assert( ieee_is_nan(var(i321, .false., corrected=.false.)))
-    call assert( ieee_is_nan(var(i321, 1, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(i321, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(i321, 1, .false., corrected=.false.)))
 
     print*,' test_int32_1dim_mask_array'
-    call assert( abs(var(i321, i321 < 5, corrected=.false.) - 5._dp/4.) < dptol)
-    call assert( ieee_is_nan((var(i321, i321 < 0, corrected=.false.))))
-    call assert( abs(var(i321, i321 == 1, corrected=.false.)) < dptol)
-    call assert( abs(var(i321, 1, i321 < 5, corrected=.false.) - 5._dp/4.) < dptol)
+    call check( abs(var(i321, i321 < 5, corrected=.false.) - 5._dp/4.) < dptol)
+    call check( ieee_is_nan((var(i321, i321 < 0, corrected=.false.))))
+    call check( abs(var(i321, i321 == 1, corrected=.false.)) < dptol)
+    call check( abs(var(i321, 1, i321 < 5, corrected=.false.) - 5._dp/4.) < dptol)
 
     !2dim
     i32 = d
     print*,' test_int32_2dim'
-    call assert( abs(var(i32, corrected=.false.) - 13._dp*11./12.) < dptol)
-    call assert( all( abs( var(i32, 1, corrected=.false.) -&
+    call check( abs(var(i32, corrected=.false.) - 13._dp*11./12.) < dptol)
+    call check( all( abs( var(i32, 1, corrected=.false.) -&
                       [20., 20., 5.]/4._dp) < dptol))
-    call assert( all( abs( var(i32, 2, corrected=.false.) -&
+    call check( all( abs( var(i32, 2, corrected=.false.) -&
                       [38._dp, 86._dp / 3._dp, 62._dp / 3._dp , 14._dp]/3._dp) < dptol))
 
     print*,' test_int32_2dim_mask'
-    call assert( ieee_is_nan(var(i32, .false., corrected=.false.)))
-    call assert( any(ieee_is_nan(var(i32, 1, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(i32, 2, .false., corrected=.false.))))
+    call check( ieee_is_nan(var(i32, .false., corrected=.false.)))
+    call check( any(ieee_is_nan(var(i32, 1, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(i32, 2, .false., corrected=.false.))))
 
     print*,' test_int32_2dim_mask_array'
-    call assert( abs(var(i32, i32 < 11, corrected=.false.) - 2.75_dp*3._dp) < dptol)
-    call assert( all( abs( var(i32, 1, i32 < 11, corrected=.false.) -&
+    call check( abs(var(i32, i32 < 11, corrected=.false.) - 2.75_dp*3._dp) < dptol)
+    call check( all( abs( var(i32, 1, i32 < 11, corrected=.false.) -&
                       [5._dp, 5._dp, 0.25_dp]) < dptol))
-    call assert( all( abs( var(i32, 2, i32 < 11, corrected=.false.) -&
+    call check( all( abs( var(i32, 2, i32 < 11, corrected=.false.) -&
                       [38._dp/3., 86._dp/9., 0.25_dp , 0.25_dp]) < dptol))
 
 
@@ -255,21 +255,21 @@ program test_varn
     i323(:,:,3)=i32*4;
 
     print*,' test_int32_3dim'
-    call assert( abs(var(i323, corrected=.false.) - 153.4_dp*35._dp/36._dp) < dptol)
-    call assert( all( abs( var(i323, 1, corrected=.false.) -&
+    call check( abs(var(i323, corrected=.false.) - 153.4_dp*35._dp/36._dp) < dptol)
+    call check( all( abs( var(i323, 1, corrected=.false.) -&
                  reshape([20._dp , 20._dp , 5._dp ,&
                           4* 20._dp , 4* 20._dp , 4* 5._dp ,&
                           16* 20._dp , 16* 20._dp , 16* 5._dp ],&
                           [size(i323,2), size(i323,3)])/4._dp)&
                  < dptol))
-    call assert( all( abs( var(i323, 2, corrected=.false.) -&
+    call check( all( abs( var(i323, 2, corrected=.false.) -&
                  reshape([38._dp, 86. / 3._dp, 62. / 3._dp , 14._dp,&
                           8* 19._dp, 8* 43. / 3._dp, 8* 31. / 3._dp, 8* 7._dp,&
                           32* 19._dp, 32* 43. / 3._dp, 32* 31. / 3._dp, 32* 7._dp],&
                           [size(i323,1), size(i323,3)] )/3._dp)&
                  < dptol))
     print*,' test_int32_3dim'
-    call assert( all(abs( var(i323, 3, corrected=.false.) -&
+    call check( all(abs( var(i323, 3, corrected=.false.) -&
                  reshape([ 7._dp/3., 21._dp, 175._dp/3.,&
                            343._dp/3., 28._dp/3., 112._dp/3.,&
                            84._dp, 448._dp/3., 189._dp,&
@@ -278,26 +278,26 @@ program test_varn
                  < dptol))
 
     print*,' test_int32_3dim_mask'
-    call assert( ieee_is_nan(var(i323, .false., corrected=.false.)))
-    call assert( any(ieee_is_nan(var(i323, 1, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(i323, 2, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(i323, 3, .false., corrected=.false.))))
+    call check( ieee_is_nan(var(i323, .false., corrected=.false.)))
+    call check( any(ieee_is_nan(var(i323, 1, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(i323, 2, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(i323, 3, .false., corrected=.false.))))
 
     print*,' test_int32_3dim_mask_array'
-    call assert( abs(var(i323, i323 < 11, corrected=.false.) -&
+    call check( abs(var(i323, i323 < 11, corrected=.false.) -&
                      7.7370242214532876_dp) < dptol)
-    call assert( all( abs( var(i323, 1, i323 < 45, corrected=.false.) -&
+    call check( all( abs( var(i323, 1, i323 < 45, corrected=.false.) -&
                   reshape([5._dp, 5._dp, 1.25_dp,  20._dp, 20._dp, 5._dp,&
                            80._dp, 80._dp, 32._dp/3.],&
                            [size(i323, 2), size(i323, 3)])) < dptol ))
-    call assert( all( abs( var(i323, 2, i323 < 45, corrected=.false.) -&
+    call check( all( abs( var(i323, 2, i323 < 45, corrected=.false.) -&
                   reshape([ 38._dp/3., 86._dp/9., 62._dp/9., 14._dp/3., 152._dp/3.,&
                             344._dp/9., 248._dp/9., 168._dp/9., 1824._dp/9.,&
                             1376._dp/9., 992._dp/9., 4._dp&
                            ],&
                   [size(i323, 1), size(i323, 3)])) < dptol ))
     print*,' test_int32_3dim_mask_array'
-    call assert( all( abs( var(i323, 3, i323 < 45, corrected=.false.) -&
+    call check( all( abs( var(i323, 3, i323 < 45, corrected=.false.) -&
                  reshape([14._dp/9., 14._dp, 350._dp/9., 686._dp/9., 56._dp/9.,&
                           224._dp/9., 56._dp, 896._dp/9., 126._dp, 1400._dp/9.,&
                           1694._dp/9., 36._dp&
@@ -307,25 +307,25 @@ program test_varn
     !cdp
     !1dim
     print*,' test_cdp_1dim'
-    call assert( abs(var(cd1, corrected=.false.) -&
+    call check( abs(var(cd1, corrected=.false.) -&
                  (var(real(cd1), corrected=.false.) +&
                   var(aimag(cd1), corrected=.false.))) < dptol)
-    call assert( abs(var(cd1, dim=1, corrected=.false.) -&
+    call check( abs(var(cd1, dim=1, corrected=.false.) -&
                  (var(real(cd1), dim=1, corrected=.false.) +&
                   var(aimag(cd1), dim=1, corrected=.false.))) < dptol)
 
     print*,' test_cdp_1dim_mask'
-    call assert( ieee_is_nan(var(cd1, .false., corrected=.false.)))
-    call assert( ieee_is_nan(var(cd1, 1, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(cd1, .false., corrected=.false.)))
+    call check( ieee_is_nan(var(cd1, 1, .false., corrected=.false.)))
 
     print*,' test_cdp_1dim_mask_array'
-    call assert( abs(var(cd1, aimag(cd1) == 0, corrected=.false.) -&
+    call check( abs(var(cd1, aimag(cd1) == 0, corrected=.false.) -&
                      var(real(cd1), aimag(cd1) == 0, corrected=.false.)) < dptol)
-    call assert( abs(var(cd1, 1, aimag(cd1) == 0, corrected=.false.) -&
+    call check( abs(var(cd1, 1, aimag(cd1) == 0, corrected=.false.) -&
                      var(real(cd1), 1, aimag(cd1) == 0, corrected=.false.)) < dptol)
-    call assert( ieee_is_nan((var(cd1, (real(cd1) == 0 .and. aimag(cd1) == 0),&
+    call check( ieee_is_nan((var(cd1, (real(cd1) == 0 .and. aimag(cd1) == 0),&
                             corrected=.false.))))
-    call assert( abs(var(cd1, (real(cd1) > 1.2 .and. aimag(cd1) == 0),&
+    call check( abs(var(cd1, (real(cd1) > 1.2 .and. aimag(cd1) == 0),&
                      corrected=.false.)) < dptol)
 
     !2dim
@@ -334,26 +334,26 @@ program test_varn
     cd(:,3) = cd1*1.5_sp
 
     print*,' test_cdp_2dim'
-    call assert( abs(var(cd, corrected=.false.) -&
+    call check( abs(var(cd, corrected=.false.) -&
                     (var(real(cd), corrected=.false.) +&
                      var(aimag(cd), corrected=.false.))) < dptol)
-    call assert( all( abs(var(cd, 1, corrected=.false.) -&
+    call check( all( abs(var(cd, 1, corrected=.false.) -&
                          (var(real(cd), 1, corrected=.false.) +&
                           var(aimag(cd), 1, corrected=.false.))) < dptol))
-    call assert( all( abs(var(cd, 2, corrected=.false.) -&
+    call check( all( abs(var(cd, 2, corrected=.false.) -&
                          (var(real(cd), 2, corrected=.false.) +&
                           var(aimag(cd), 2, corrected=.false.))) < dptol))
 
     print*,' test_cdp_2dim_mask'
-    call assert( ieee_is_nan(var(cd, .false., corrected=.false.)))
-    call assert( any(ieee_is_nan(var(cd, 1, .false., corrected=.false.))))
-    call assert( any(ieee_is_nan(var(cd, 2, .false., corrected=.false.))))
+    call check( ieee_is_nan(var(cd, .false., corrected=.false.)))
+    call check( any(ieee_is_nan(var(cd, 1, .false., corrected=.false.))))
+    call check( any(ieee_is_nan(var(cd, 2, .false., corrected=.false.))))
 
     print*,' test_cdp_2dim_mask_array'
-    call assert( abs(var(cd, aimag(cd) == 0, corrected=.false.) -&
+    call check( abs(var(cd, aimag(cd) == 0, corrected=.false.) -&
                      var(real(cd), aimag(cd) == 0, corrected=.false.)) < dptol)
-    call assert( all( abs( var(cd, 1, aimag(cd) == 0, corrected=.false.) -&
+    call check( all( abs( var(cd, 1, aimag(cd) == 0, corrected=.false.) -&
                            var(real(cd), 1, aimag(cd) == 0, corrected=.false.)) < dptol))
-    call assert( any( ieee_is_nan( var(cd, 2, aimag(cd) == 0, corrected=.false.))))
+    call check( any( ieee_is_nan( var(cd, 2, aimag(cd) == 0, corrected=.false.))))
 
 end program
diff --git a/src/tests/test_always_fail.f90 b/src/tests/test_always_fail.f90
index 48ab13ed1..64d0fa393 100644
--- a/src/tests/test_always_fail.f90
+++ b/src/tests/test_always_fail.f90
@@ -1,8 +1,8 @@
 program test_always_fail
 
-use stdlib_experimental_error, only : assert
+use stdlib_experimental_error, only: check
 implicit none
 
-call assert(.false.)
+call check(.false.)
 
 end program
diff --git a/src/tests/test_always_skip.f90 b/src/tests/test_always_skip.f90
index 84092d2fc..1cd6bbb96 100644
--- a/src/tests/test_always_skip.f90
+++ b/src/tests/test_always_skip.f90
@@ -1,8 +1,8 @@
 program test_always_skip
 
-use stdlib_experimental_error, only : assert
+use stdlib_experimental_error, only: check
 implicit none
 
-call assert(.false., 77)
+call check(.false., code=77)
 
 end program