diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index cf6507ed1..81d58d772 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -39,6 +39,9 @@ jobs:
       if: contains(matrix.os, 'ubuntu')
       run: ci/install_cmake.sh
 
+    - name: Install fypp
+      run: pip install --upgrade fypp
+
     - name: Install GFortran Linux
       if: contains( matrix.os, 'ubuntu')
       run: |
@@ -54,7 +57,7 @@ jobs:
       run: brew install gcc@${GCC_V} || brew upgrade gcc@${GCC_V} || true
 
     - name: Configure with CMake
-      run: cmake -Wdev -DCMAKE_BUILD_TYPE=Release -S . -B build
+      run: cmake -Wdev -DCMAKE_BUILD_TYPE=Release -DCMAKE_MAXIMUM_RANK=4 -S . -B build
 
     - name: Build and compile
       run: cmake --build build
@@ -77,6 +80,6 @@ jobs:
     - name: Test manual makefiles
       if: contains(matrix.os, 'ubuntu') && contains(matrix.gcc_v, '9')
       run: |
-        make -f Makefile.manual
+        make -f Makefile.manual FYPPFLAGS="-DMAXRANK=4"
         make -f Makefile.manual test
         make -f Makefile.manual clean
diff --git a/.github/workflows/ci_windows.yml b/.github/workflows/ci_windows.yml
index 0d34a9a5d..d8b7d4b8a 100644
--- a/.github/workflows/ci_windows.yml
+++ b/.github/workflows/ci_windows.yml
@@ -15,7 +15,10 @@ jobs:
     steps:
     - uses: actions/checkout@v1
 
-    - run: cmake -G "MinGW Makefiles" -DCMAKE_SH="CMAKE_SH-NOTFOUND" -Wdev -B build -DCMAKE_BUILD_TYPE=Debug -DCMAKE_Fortran_FLAGS_DEBUG="-Wall -Wextra -Wimplicit-interface -fPIC -g -fcheck=all -fbacktrace"
+    - name: Install fypp
+      run: pip install fypp
+
+    - run: cmake -G "MinGW Makefiles" -DCMAKE_SH="CMAKE_SH-NOTFOUND" -Wdev -B build -DCMAKE_BUILD_TYPE=Debug -DCMAKE_Fortran_FLAGS_DEBUG="-Wall -Wextra -Wimplicit-interface -fPIC -g -fcheck=all -fbacktrace" -DCMAKE_MAXIMUM_RANK=4
 
       env:
         FC: gfortran
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b858ebd6b..676bd09c2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,6 +15,11 @@ endif()
 include(CheckFortranSourceCompiles)
 include(CheckFortranSourceRuns)
 check_fortran_source_compiles("error stop i; end" f18errorstop SRC_EXT f90)
+check_fortran_source_compiles("real, allocatable :: array(:, :, :, :, :, :, :, :, :, :); end" f03rank SRC_EXT f90)
 check_fortran_source_runs("use, intrinsic :: iso_fortran_env, only : real128; real(real128) :: x; x = x+1; end" f03real128)
 
+if(DEFINED CMAKE_MAXIMUM_RANK)
+  set(CMAKE_MAXIMUM_RANK ${CMAKE_MAXIMUM_RANK})
+endif()
+
 add_subdirectory(src)
diff --git a/Makefile.manual b/Makefile.manual
index 3fd085ed6..b8280b102 100644
--- a/Makefile.manual
+++ b/Makefile.manual
@@ -2,9 +2,11 @@
 
 FC = gfortran
 FFLAGS = -Wall -Wextra -Wimplicit-interface -fPIC -g -fcheck=all
+FYPPFLAGS=
 
 export FC
 export FFLAGS
+export FYPPFLAGS
 
 .PHONY: all clean test
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4a8d3fe03..add1e9170 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,51 @@
+### Pre-process: .fpp -> .f90 via Fypp
+
+# Create a list of the files to be preprocessed
+set(fppFiles
+    stdlib_experimental_stats.fypp
+    stdlib_experimental_stats_mean.fypp
+)
+
+# Pre-process
+foreach(infileName IN LISTS fppFiles)
+
+    # Generate output file name
+    string(REGEX REPLACE ".fypp\$" ".f90" outfileName "${infileName}")
+
+    # Create the full path for the new file
+    set(outfile "${CMAKE_CURRENT_BINARY_DIR}/${outfileName}")
+
+    # Generate input file name
+    set(infile "${CMAKE_CURRENT_SOURCE_DIR}/${infileName}")
+
+    # Custom command to do the processing
+    if(DEFINED CMAKE_MAXIMUM_RANK)
+        add_custom_command(
+            OUTPUT "${outfile}"
+            COMMAND fypp -DMAXRANK=${CMAKE_MAXIMUM_RANK} "${infile}" "${outfile}"
+            MAIN_DEPENDENCY "${infile}"
+            VERBATIM)
+    elseif(f03rank)
+        add_custom_command(
+            OUTPUT "${outfile}"
+            COMMAND fypp "${infile}" "${outfile}"
+            MAIN_DEPENDENCY "${infile}"
+            VERBATIM)
+    else()
+        add_custom_command(
+            OUTPUT "${outfile}"
+            COMMAND fypp -DVERSION90 "${infile}" "${outfile}"
+            MAIN_DEPENDENCY "${infile}"
+            VERBATIM)
+    endif()
+
+    # Finally add output file to a list
+    set(outFiles ${outFiles} "${outfile}")
+
+endforeach(infileName)
+
+
+
 set(SRC
     stdlib_experimental_ascii.f90
     stdlib_experimental_io.f90
@@ -5,6 +53,7 @@ set(SRC
     stdlib_experimental_kinds.f90
     stdlib_experimental_optval.f90
     stdlib_experimental_system.F90
+    ${outFiles}
 )
 
 add_library(fortran_stdlib ${SRC})
diff --git a/src/Makefile.manual b/src/Makefile.manual
index 0807f3ff1..60e86874b 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -3,7 +3,9 @@ SRC = stdlib_experimental_ascii.f90 \
       stdlib_experimental_io.f90 \
       stdlib_experimental_optval.f90 \
       stdlib_experimental_kinds.f90 \
-      f18estop.f90
+      f18estop.f90 \
+      stdlib_experimental_stats.f90 \
+      stdlib_experimental_stats_mean.f90
 
 LIB = libstdlib.a
 
@@ -26,6 +28,8 @@ clean:
 %.o: %.f90
 	$(FC) $(FFLAGS) -c $<
 
+%.f90: %.fypp
+	fypp $(FYPPFLAGS) $< $@
 
 # Fortran module dependencies
 f18estop.o: stdlib_experimental_error.o
@@ -34,3 +38,9 @@ stdlib_experimental_io.o: \
 	stdlib_experimental_optval.o \
 	stdlib_experimental_kinds.o
 stdlib_experimental_optval.o: stdlib_experimental_kinds.o
+stdlib_experimental_stats_mean.o: \
+	stdlib_experimental_optval.o \
+	stdlib_experimental_kinds.o \
+	stdlib_experimental_stats.o
+stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp
+stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp
diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp
new file mode 100644
index 000000000..8cbf7f3df
--- /dev/null
+++ b/src/stdlib_experimental_stats.fypp
@@ -0,0 +1,126 @@
+module stdlib_experimental_stats
+
+#:set VERSION90 = defined('VERSION90')
+#:set REALKINDS = ["sp", "dp", "qp"]
+#:set INTKINDS = ["int8", "int16", "int32", "int64"]
+#:set REALTYPES = ["real({})".format(k) for k in REALKINDS]
+#:set INTTYPES = ["integer({})".format(k) for k in INTKINDS]
+#:set iktr = list(zip(range(len(REALKINDS)), REALKINDS, REALTYPES))
+#:set ikti = list(zip(range(len(INTKINDS)), INTKINDS, INTTYPES))
+
+use stdlib_experimental_kinds, only: sp, dp, qp, &
+    int8, int16, int32, int64
+implicit none
+private
+! Public API
+public :: mean
+
+interface mean
+#:for i1, k1, t1 in iktr
+    module function mean_1_${k1}$_${k1}$(x) result(res)
+        ${t1}$, intent(in) :: x(:)
+        ${t1}$ :: res
+    end function mean_1_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+    module function mean_1_${k1}$_dp(x) result(res)
+        ${t1}$, intent(in) :: x(:)
+        real(dp) :: res
+    end function mean_1_${k1}$_dp
+#:endfor
+
+
+#:for i1, k1, t1 in iktr
+    module function mean_2_all_${k1}$_${k1}$(x) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        ${t1}$ :: res
+    end function mean_2_all_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+    module function mean_2_all_${k1}$_dp(x) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        real(dp) :: res
+    end function mean_2_all_${k1}$_dp
+#:endfor
+
+#:for i1, k1, t1 in iktr
+    module function mean_2_${k1}$_${k1}$(x, dim) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        integer, intent(in) :: dim
+        ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1 < dim ))
+    end function mean_2_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+    module function mean_2_${k1}$_dp(x, dim) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        integer, intent(in) :: dim
+        real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1 < dim ))
+    end function mean_2_${k1}$_dp
+#:endfor
+
+
+#:def ranksuffix(rank)
+#{if rank > 0}#(${":" + ",:" * (rank - 1)}$)#{endif}#
+#:enddef
+
+#:if defined('MAXRANK')
+#:set ranks = range(3,MAXRANK+1)
+#:elif VERSION90
+#:set ranks = range(3,8)
+#:else
+#:set ranks = range(3,16)
+#:endif
+
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+    module function mean_${rank}$_all_${k1}$_${k1}$(x) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        ${t1}$ :: res
+    end function mean_${rank}$_all_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+    module function mean_${rank}$_all_${k1}$_dp(x) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        real(dp) :: res
+    end function mean_${rank}$_all_${k1}$_dp
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+    module function mean_${rank}$_${k1}$_${k1}$(x, dim) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: dim
+        ${t1}$ :: res( &
+#:for imerge in range(1,rank-1)
+                      merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                      merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+    end function mean_${rank}$_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+    module function mean_${rank}$_${k1}$_dp(x, dim) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: dim
+        real(dp) :: res( &
+#:for imerge in range(1,rank-1)
+                      merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                      merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+    end function mean_${rank}$_${k1}$_dp
+#:endfor
+#:endfor
+
+end interface
+
+end module
diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md
new file mode 100644
index 000000000..4cb850a5c
--- /dev/null
+++ b/src/stdlib_experimental_stats.md
@@ -0,0 +1,43 @@
+# Descriptive statistics
+
+## Implemented
+
+ * `mean`
+
+## `mean` - mean of array elements
+
+### Description
+
+Returns the mean of all the elements of `array`, or of the elements of `array` along dimension `dim` if provided.
+
+### Syntax
+
+`result = mean(array)`
+
+`result = mean(array, dim)`
+
+### Arguments
+
+`array`: Shall be an array of type `integer`, or `real`.
+
+`dim`: Shall be a scalar of type `integer` with a value in the range from 1 to n, where n is the rank of `array`.
+
+### Return value
+
+If `array` is of type `real`, the result is of the same type as `array`.
+If `array` is of type `integer`, the result is of type `double precision`.
+
+If `dim` is absent, a scalar with the mean of all elements in `array` is returned. Otherwise, an array of rank n-1, where n equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
+
+### Example
+
+```fortran
+program test
+    use stdlib_experimental_stats, only: mean
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    print *, mean(x)                            !returns 21.
+    print *, mean( reshape(x, [ 2, 3 ] ))       !returns 21.
+    print *, mean( reshape(x, [ 2, 3 ] ), 1)    !returns [ 3., 7., 11. ]
+end program
+```
diff --git a/src/stdlib_experimental_stats_mean.fypp b/src/stdlib_experimental_stats_mean.fypp
new file mode 100644
index 000000000..bb2daee4f
--- /dev/null
+++ b/src/stdlib_experimental_stats_mean.fypp
@@ -0,0 +1,178 @@
+submodule (stdlib_experimental_stats) stdlib_experimental_stats_mean
+
+#:set VERSION90 = defined('VERSION90')
+#:set REALKINDS = ["sp", "dp", "qp"]
+#:set INTKINDS = ["int8", "int16", "int32", "int64"]
+#:set REALTYPES = ["real({})".format(k) for k in REALKINDS]
+#:set INTTYPES = ["integer({})".format(k) for k in INTKINDS]
+#:set iktr = list(zip(range(len(REALKINDS)), REALKINDS, REALTYPES))
+#:set ikti = list(zip(range(len(INTKINDS)), INTKINDS, INTTYPES))
+
+use stdlib_experimental_error, only: error_stop
+implicit none
+
+contains
+
+#:for i1, k1, t1 in iktr
+module function mean_1_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x(:)
+    ${t1}$ :: res
+
+    res = sum(x) / real(size(x, kind = int64), ${k1}$)
+
+end function mean_1_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_1_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x(:)
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x, kind = int64), dp)
+
+end function mean_1_${k1}$_dp
+#:endfor
+
+
+#:for i1, k1, t1 in iktr
+module function mean_2_all_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    ${t1}$ :: res
+
+    res = sum(x) / real(size(x, kind = int64), ${k1}$)
+
+end function mean_2_all_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_2_all_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x, kind = int64), dp)
+
+end function mean_2_all_${k1}$_dp
+#:endfor
+
+#:for i1, k1, t1 in iktr
+module function mean_2_${k1}$_${k1}$(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    integer, intent(in) :: dim
+    ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1 < dim ))
+
+    select case(dim)
+     case(1)
+       res = sum(x, 1) / real(size(x, 1), ${k1}$)
+     case(2)
+       res = sum(x, 2) / real(size(x, 2), ${k1}$)
+     case default
+       call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+end function mean_2_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_2_${k1}$_dp(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    integer, intent(in) :: dim
+    real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1 < dim ))
+
+    select case(dim)
+     case(1)
+       res = sum(real(x, dp), 1) / real(size(x, 1), dp)
+     case(2)
+       res = sum(real(x, dp), 2) / real(size(x, 2), dp)
+     case default
+       call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+end function mean_2_${k1}$_dp
+#:endfor
+
+
+#:def ranksuffix(rank)
+#{if rank > 0}#(${":" + ",:" * (rank - 1)}$)#{endif}#
+#:enddef
+
+#:if defined('MAXRANK')
+#:set ranks = range(3,MAXRANK+1)
+#:elif VERSION90
+#:set ranks = range(3,8)
+#:else
+#:set ranks = range(3,16)
+#:endif
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+module function mean_${rank}$_all_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    ${t1}$ :: res
+
+    res = sum(x) / real(size(x, kind = int64), ${k1}$)
+
+end function mean_${rank}$_all_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+module function mean_${rank}$_all_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x, kind = int64), dp)
+
+end function mean_${rank}$_all_${k1}$_dp
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+module function mean_${rank}$_${k1}$_${k1}$(x, dim) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    integer, intent(in) :: dim
+    ${t1}$ :: res( &
+#:for imerge in range(1,rank-1)
+                  merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                  merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+
+    select case(dim)
+#:for fi in range(1,rank+1)
+     case(${fi}$)
+       res=sum(x, ${fi}$) / real(size(x, ${fi}$), ${k1}$)
+#:endfor
+     case default
+       call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+end function mean_${rank}$_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+module function mean_${rank}$_${k1}$_dp(x, dim) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    integer, intent(in) :: dim
+    real(dp) :: res( &
+#:for imerge in range(1,rank-1)
+                  merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                  merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+
+    select case(dim)
+#:for fi in range(1,rank+1)
+     case(${fi}$)
+       res=sum(real(x, dp), ${fi}$) / real(size(x, ${fi}$), dp)
+#:endfor
+     case default
+       call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+end function mean_${rank}$_${k1}$_dp
+#:endfor
+#:endfor
+
+end submodule
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index e72592579..df5bd0a09 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -9,6 +9,7 @@ endmacro(ADDTEST)
 add_subdirectory(ascii)
 add_subdirectory(io)
 add_subdirectory(optval)
+add_subdirectory(stats)
 add_subdirectory(system)
 
 ADDTEST(always_skip)
diff --git a/src/tests/Makefile.manual b/src/tests/Makefile.manual
index 5d2debb62..3b1e4dff2 100644
--- a/src/tests/Makefile.manual
+++ b/src/tests/Makefile.manual
@@ -4,13 +4,16 @@ all:
 	$(MAKE) -f Makefile.manual --directory=ascii
 	$(MAKE) -f Makefile.manual --directory=io
 	$(MAKE) -f Makefile.manual --directory=optval
+	$(MAKE) -f Makefile.manual --directory=stats
 
 test:
 	$(MAKE) -f Makefile.manual --directory=ascii test
 	$(MAKE) -f Makefile.manual --directory=io test
 	$(MAKE) -f Makefile.manual --directory=optval test
+	$(MAKE) -f Makefile.manual --directory=stats test
 
 clean:
 	$(MAKE) -f Makefile.manual --directory=ascii clean
 	$(MAKE) -f Makefile.manual --directory=io clean
 	$(MAKE) -f Makefile.manual --directory=optval clean
+	$(MAKE) -f Makefile.manual --directory=stats clean
diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt
new file mode 100644
index 000000000..e86fe23d7
--- /dev/null
+++ b/src/tests/stats/CMakeLists.txt
@@ -0,0 +1,9 @@
+ADDTEST(mean)
+
+if(DEFINED CMAKE_MAXIMUM_RANK)
+    if(${CMAKE_MAXIMUM_RANK} GREATER 7)
+        ADDTEST(mean_f03)
+    endif()
+elseif(f03rank)
+    ADDTEST(mean_f03)
+endif()
diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual
new file mode 100644
index 000000000..9faf154cb
--- /dev/null
+++ b/src/tests/stats/Makefile.manual
@@ -0,0 +1,3 @@
+PROGS_SRC = test_mean.f90
+
+include ../Makefile.manual.test.mk
diff --git a/src/tests/stats/array3.dat b/src/tests/stats/array3.dat
new file mode 100644
index 000000000..13b583f89
--- /dev/null
+++ b/src/tests/stats/array3.dat
@@ -0,0 +1,16 @@
+1.000000000000000021e-08 9.199998759392489944e+01
+1.024113254885563425e-08 9.199998731474968849e+01
+1.048233721895820948e-08 9.199998703587728244e+01
+1.072361403187881949e-08 9.199998675729767683e+01
+1.096496300919481796e-08 9.199998647900135040e+01
+1.120638417249036630e-08 9.199998620097916557e+01
+1.144787754335570897e-08 9.199998592322251056e+01
+1.168944314338753750e-08 9.199998564572304360e+01
+1.193108099418952317e-08 9.199998536847290609e+01
+1.217279111737088596e-08 9.199998509146449521e+01
+1.241457353454836993e-08 9.199998481469057765e+01
+1.265642826734443823e-08 9.199998453814424693e+01
+1.289835533738818635e-08 9.199998426181879552e+01
+1.314035476631514857e-08 9.199998398570787117e+01
+1.338242657576766519e-08 9.199998370980536322e+01
+1.362457078739434161e-08 9.199998343410533153e+01
diff --git a/src/tests/stats/test_mean.f90 b/src/tests/stats/test_mean.f90
new file mode 100644
index 000000000..d16071101
--- /dev/null
+++ b/src/tests/stats/test_mean.f90
@@ -0,0 +1,75 @@
+program test_mean
+use stdlib_experimental_error, only: assert
+use stdlib_experimental_kinds, only: sp, dp, int32, int64
+use stdlib_experimental_io, only: loadtxt
+use stdlib_experimental_stats, only: mean
+implicit none
+
+real(sp), allocatable :: s(:, :)
+real(dp), allocatable :: d(:, :)
+
+real(dp), allocatable :: d3(:, :, :)
+real(dp), allocatable :: d4(:, :, :, :)
+
+
+!sp
+call loadtxt("array3.dat", s)
+
+call assert( mean(s) - sum(s)/real(size(s), sp) == 0.0_sp)
+call assert( sum( abs( mean(s,1) - sum(s,1)/real(size(s,1), sp) )) == 0.0_sp)
+call assert( sum( abs( mean(s,2) - sum(s,2)/real(size(s,2), sp) )) == 0.0_sp)
+
+
+!dp
+call loadtxt("array3.dat", d)
+
+call assert( mean(d) - sum(d)/real(size(d), dp) == 0.0_dp)
+call assert( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) == 0.0_dp)
+
+
+!int32
+call loadtxt("array3.dat", d)
+
+call assert( mean(int(d, int32)) - sum(real(int(d, int32),dp))/real(size(d), dp) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int32),1) - sum(real(int(d, int32),dp),1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int32),2) - sum(real(int(d, int32),dp),2)/real(size(d,2), dp) )) == 0.0_dp)
+
+
+!int64
+call loadtxt("array3.dat", d)
+
+call assert( mean(int(d, int64)) - sum(real(int(d, int64),dp))/real(size(d), dp) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int64),1) - sum(real(int(d, int64),dp),1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int64),2) - sum(real(int(d, int64),dp),2)/real(size(d,2), dp) )) == 0.0_dp)
+
+
+!dp rank 3
+allocate(d3(size(d,1),size(d,2),3))
+d3(:,:,1)=d;
+d3(:,:,2)=d*1.5_dp;
+d3(:,:,3)=d*4._dp;
+
+call assert( mean(d3) - sum(d3)/real(size(d3), dp) == 0.0_dp)
+call assert( sum( abs( mean(d3,1) - sum(d3,1)/real(size(d3,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d3,2) - sum(d3,2)/real(size(d3,2), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d3,3) - sum(d3,3)/real(size(d3,3), dp) )) == 0.0_dp)
+
+
+!dp rank 4
+allocate(d4(size(d,1),size(d,2),3,9))
+d4 = 1.
+d4(:,:,1,1)=d;
+d4(:,:,2,1)=d*1.5_dp;
+d4(:,:,3,1)=d*4._dp;
+d4(:,:,3,9)=d*4._dp;
+
+call assert( mean(d4) - sum(d4)/real(size(d4), dp) == 0.0_dp)
+call assert( sum( abs( mean(d4,1) - sum(d4,1)/real(size(d4,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d4,2) - sum(d4,2)/real(size(d4,2), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d4,3) - sum(d4,3)/real(size(d4,3), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d4,4) - sum(d4,4)/real(size(d4,4), dp) )) == 0.0_dp)
+
+contains
+
+end program
diff --git a/src/tests/stats/test_mean_f03.f90 b/src/tests/stats/test_mean_f03.f90
new file mode 100644
index 000000000..fbb09b7ba
--- /dev/null
+++ b/src/tests/stats/test_mean_f03.f90
@@ -0,0 +1,38 @@
+program test_mean
+use stdlib_experimental_error, only: assert
+use stdlib_experimental_kinds, only: sp, dp, int32, int64
+use stdlib_experimental_io, only: loadtxt
+use stdlib_experimental_stats, only: mean
+implicit none
+
+real(dp), allocatable :: d(:, :)
+real(dp), allocatable :: d8(:, :, :, :, :, :, :, :)
+
+
+!dp
+call loadtxt("array3.dat", d)
+
+call assert( mean(d) - sum(d)/real(size(d), dp) == 0.0_dp)
+call assert( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) == 0.0_dp)
+
+!dp rank 8
+allocate(d8(size(d,1), size(d,2), 3, 4, 5, 6, 7, 8))
+d8(:, :, 1, 4, 5 ,6 ,7 ,8)=d;
+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) == 0.0_dp)
+
+call assert( sum( abs( mean(d8,1) - sum(d8,1)/real(size(d8,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,2) - sum(d8,2)/real(size(d8,2), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,3) - sum(d8,3)/real(size(d8,3), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,4) - sum(d8,4)/real(size(d8,4), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,5) - sum(d8,5)/real(size(d8,5), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,6) - sum(d8,6)/real(size(d8,6), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,7) - sum(d8,7)/real(size(d8,7), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,8) - sum(d8,8)/real(size(d8,8), dp) )) == 0.0_dp)
+
+contains
+
+end program