diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md
index de475c60d..b8a4f7f78 100644
--- a/doc/specs/stdlib_stats.md
+++ b/doc/specs/stdlib_stats.md
@@ -26,6 +26,10 @@ The Pearson correlation between two rows (or columns), say `x` and `y`, of `arra
 
 `result = [[stdlib_stats(module):corr(interface)]](array, dim [, mask])`
 
+### Class
+
+Generic subroutine
+
 ### Arguments
 
 `array`: Shall be a rank-1 or a rank-2 array of type `integer`, `real`, or `complex`.
@@ -83,6 +87,10 @@ The scaling can be changed with the logical argument `corrected`. If `corrected`
 
 `result = [[stdlib_stats(module):cov(interface)]](array, dim [, mask [, corrected]])`
 
+### Class
+
+Generic subroutine
+
 ### Arguments
 
 `array`: Shall be a rank-1 or a rank-2 array of type `integer`, `real`, or `complex`.
@@ -134,6 +142,10 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a
 
 `result = [[stdlib_stats(module):mean(interface)]](array, dim [, mask])`
 
+### Class
+
+Generic subroutine
+
 ### Arguments
 
 `array`: Shall be an array of type `integer`, `real`, or `complex`.
@@ -166,6 +178,80 @@ program demo_mean
 end program demo_mean
 ```
 
+## `median` - median of array elements
+
+### Status
+
+Experimental
+
+### Description
+
+Returns the median of all the elements of `array`, or of the elements of `array`
+along dimension `dim` if provided, and if the corresponding element in `mask` is
+`true`.
+
+The median of the elements of `array` is defined as the "middle"
+element, after that the elements are sorted in an increasing order, e.g. `array_sorted =
+sort(array)`. If `n = size(array)` is an even number, the median is:
+
+```
+median(array) = array_sorted( floor( (n + 1) / 2.))
+```
+
+and if `n` is an odd number, the median is:
+
+```
+median(array) = mean( array_sorted( floor( (n + 1) / 2.):floor( (n + 1) / 2.) + 1 ) )
+```
+
+The current implementation is a quite naive implementation that relies on sorting
+the whole array, using the subroutine `[[stdlib_sorting(module):ord_sort(interface)]]`
+provided by the `[[stdlib_sorting(module)]]` module.
+
+### Syntax
+
+`result = [[stdlib_stats(module):median(interface)]](array [, mask])`
+
+`result = [[stdlib_stats(module):median(interface)]](array, dim [, mask])`
+
+### Class
+
+Generic subroutine
+
+### 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`.
+
+`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
+
+### Return value
+
+If `array` is of type `real`, the result is of type `real` with the same kind as `array`.
+If `array` is of type `real` and contains IEEE `NaN`, the result is IEEE `NaN`.
+If `array` is of type `integer`, the result is of type `real(dp)`.
+
+If `dim` is absent, a scalar with the median 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.
+
+If `mask` is specified, the result is the median of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`.
+
+
+### Example
+
+```fortran
+program demo_median
+    use stdlib_stats, only: median
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    real :: y(1:2, 1:3) = reshape([ 1., 2., 3., 4., 5., 6. ], [ 2, 3])
+    print *, median(x)                                  !returns 3.5
+    print *, median(y)                                  !returns 3.5
+    print *, median(y, 1)                               !returns [ 1.5, 3.5, 5.5 ]
+    print *, median(y, 1,y > 3.)                        !returns [ NaN, 4.0, 5.5 ]
+end program demo_median
+```
+
 ## `moment` - central moments of array elements
 
 ### Status
@@ -199,6 +285,10 @@ The _k_-th order moment about `center` is defined as :
 
 `result = [[stdlib_stats(module):moment(interface)]](array, order, dim [, center [, mask]])`
 
+### Class
+
+Generic subroutine
+
 ### Arguments
 
 `array`: Shall be an array of type `integer`, `real`, or `complex`.
@@ -264,6 +354,10 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin
 
 `result = [[stdlib_stats(module):var(interface)]](array, dim [, mask [, corrected]])`
 
+### Class
+
+Generic subroutine
+
 ### Arguments
 
 `array`: Shall be an array of type `integer`, `real`, or `complex`.
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a059c4a5e..6d7f8483f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -18,6 +18,7 @@ set(fppFiles
     stdlib_stats_corr.fypp
     stdlib_stats_cov.fypp
     stdlib_stats_mean.fypp
+    stdlib_stats_median.fypp
     stdlib_stats_moment.fypp
     stdlib_stats_moment_all.fypp
     stdlib_stats_moment_mask.fypp
diff --git a/src/Makefile.manual b/src/Makefile.manual
index 28d00a6cc..02ee62e83 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -18,6 +18,7 @@ SRCFYPP =\
         stdlib_stats_corr.fypp \
         stdlib_stats_cov.fypp \
         stdlib_stats_mean.fypp \
+        stdlib_stats_median.fypp \
         stdlib_stats_moment.fypp \
         stdlib_stats_moment_all.fypp \
         stdlib_stats_moment_mask.fypp \
@@ -108,6 +109,11 @@ stdlib_stats_mean.o: \
 	stdlib_optval.o \
 	stdlib_kinds.o \
 	stdlib_stats.o
+stdlib_stats_median.o: \
+	stdlib_optval.o \
+	stdlib_kinds.o \
+	stdlib_sorting.o \
+	stdlib_stats.o
 stdlib_stats_moment.o: \
 	stdlib_optval.o \
 	stdlib_kinds.o \
diff --git a/src/common.fypp b/src/common.fypp
index 4c4fbe18b..0c6e86840 100644
--- a/src/common.fypp
+++ b/src/common.fypp
@@ -156,4 +156,67 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
   #:endcall
 #:enddef
 
+#!
+#! Generates an array rank suffix for subarrays along a dimension
+#!
+#! Args:
+#!   varname (str): Name of the variable to be used as origin
+#!   rank (int): Rank of the original variable
+#!   dim (int): Dimension of the variable
+#!
+#! Returns:
+#!   Array rank suffix string enclosed in braces
+#!
+#! E.g.,
+#!  select_subvector('j', 5, 2)
+#!   -> (j1, :, j3, j4, j5)
+#!
+#! Used, e.g., in
+#!  stdlib_stats_median.fypp
+#!
+#:def select_subvector(varname, rank, idim)
+  #:assert rank > 0
+  #:call join_lines(joinstr=", ", prefix="(", suffix=")")
+    #:for i in range(1, idim)
+      ${varname}$${i}$
+    #:endfor
+    :
+    #:for i in range(idim + 1, rank + 1)
+      ${varname}$${i}$
+    #:endfor
+  #:endcall
+#:enddef
+
+#!
+#! Generates an array rank suffix for arrays
+#!
+#! Args:
+#!   varname (str): Name of the variable to be used as origin
+#!   rank (int): Rank of the original array variable
+#!   idim (int): Dimension of the variable dropped
+#!
+#! Returns:
+#!   Array rank suffix string enclosed in braces
+#!
+#! E.g.,
+#!  reduce_subvector('j', 5, 2)
+#!   -> (j1, j3, j4, j5)
+#!
+#! Used, e.g., in
+#!  stdlib_stats_median.fypp
+#!
+#:def reduce_subvector(varname, rank, idim)
+  #:assert rank > 0
+  #:if rank > 1
+    #:call join_lines(joinstr=", ", prefix="(", suffix=")")
+      #:for i in range(1, idim)
+        ${varname}$${i}$
+      #:endfor
+      #:for i in range(idim + 1, rank + 1)
+        ${varname}$${i}$
+      #:endfor
+    #:endcall
+  #:endif
+#:enddef
+
 #:endmute
diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp
index 4aa9edaa1..038f27752 100644
--- a/src/stdlib_stats.fypp
+++ b/src/stdlib_stats.fypp
@@ -2,6 +2,7 @@
 #:set RANKS = range(1, MAXRANK + 1)
 #:set REDRANKS = range(2, MAXRANK + 1)
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+#:set IR_KINDS_TYPES_OUTPUT = list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) + list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS))
 module stdlib_stats
   !! Provides support for various statistical methods. This includes currently
   !! descriptive statistics
@@ -11,14 +12,14 @@ module stdlib_stats
   implicit none
   private
   ! Public API
-  public :: corr, cov, mean, moment, var
+  public :: corr, cov, mean, median, moment, var
 
 
   interface corr
     !! version: experimental
     !!
     !! Pearson correlation of array elements
-    !! ([Specification](../page/specs/stdlib_stats.html#description))
+    !! ([Specification](../page/specs/stdlib_stats.html#corr-pearson-correlation-of-array-elements))
     #:for k1, t1 in RC_KINDS_TYPES
       #:set RName = rname("corr",1, t1, k1)
       module function ${RName}$(x, dim, mask) result(res)
@@ -110,7 +111,7 @@ module stdlib_stats
     !! version: experimental
     !!
     !! Covariance of array elements
-    !! ([Specification](../page/specs/stdlib_stats.html#description_1))
+    !! ([Specification](../page/specs/stdlib_stats.html#cov-covariance-of-array-elements))
     #:for k1, t1 in RC_KINDS_TYPES
       #:set RName = rname("cov",1, t1, k1)
       module function ${RName}$(x, dim, mask, corrected) result(res)
@@ -209,7 +210,7 @@ module stdlib_stats
     !! version: experimental
     !!
     !! Mean of array elements
-    !! ([Specification](../page/specs/stdlib_stats.html#description_2))
+    !! ([Specification](../page/specs/stdlib_stats.html#mean-mean-of-array-elements))
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("mean_all",rank, t1, k1)
@@ -305,11 +306,65 @@ module stdlib_stats
   end interface mean
 
 
+  interface median
+    !! version: experimental
+    !!
+    !! Median of array elements
+    !! ([Specification](../page/specs/stdlib_stats.html#median-median-of-array-elements))
+    #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+      #:for rank in RANKS
+        #:set name = rname("median_all",rank, t1, k1, o1)
+        module function ${name}$ (x, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          logical, intent(in), optional :: mask
+          real(${o1}$) :: res
+        end function ${name}$
+      #:endfor
+    #:endfor
+  
+    #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+      #:for rank in RANKS
+        #:set name = rname("median",rank, t1, k1, o1)
+        module function ${name}$(x, dim, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: dim
+          logical, intent(in), optional :: mask
+          real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$
+        end function ${name}$
+      #:endfor
+    #:endfor
+  
+    #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+      #:for rank in RANKS
+        #:set name = rname('median_all_mask',rank, t1, k1, o1)
+        module function ${name}$(x, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          logical, intent(in) :: mask${ranksuffix(rank)}$
+          real(${o1}$) :: res
+        end function ${name}$
+      #:endfor
+    #:endfor
+  
+    #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+      #:for rank in RANKS
+        #:set name = rname('median_mask',rank, t1, k1, o1)
+        module function  ${name}$(x, dim, mask) result(res)
+          ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+          integer, intent(in) :: dim
+          logical, intent(in) :: mask${ranksuffix(rank)}$
+          real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$
+        end function ${name}$
+      #:endfor
+    #:endfor
+
+  end interface
+
+
   interface var
     !! version: experimental
     !!
     !! Variance of array elements
-    !! ([Specification](../page/specs/stdlib_stats.html#description_4))
+    !! ([Specification](../page/specs/stdlib_stats.html#var-variance-of-array-elements))
 
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
@@ -418,7 +473,7 @@ module stdlib_stats
     !! version: experimental
     !!
     !! Central moment of array elements
-    !! ([Specification](../page/specs/stdlib_stats.html#description_3))
+    !! ([Specification](../page/specs/stdlib_stats.html#moment-central-moments-of-array-elements))
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("moment_all",rank, t1, k1)
diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp
new file mode 100644
index 000000000..5adb7b256
--- /dev/null
+++ b/src/stdlib_stats_median.fypp
@@ -0,0 +1,255 @@
+#:include "common.fypp"
+#:set RANKS = range(1, MAXRANK + 1)
+#:set IR_KINDS_TYPES_OUTPUT = list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) + list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS))
+
+submodule (stdlib_stats) stdlib_stats_median
+
+  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_nan
+  use stdlib_error, only: error_stop
+  use stdlib_optval, only: optval
+  ! Use "ord_sort" rather than "sort" because the former can be much faster for arrays
+  ! that are already partly sorted. While it is slightly slower for random arrays,
+  ! ord_sort seems a better overall choice.
+  use stdlib_sorting, only: sort => ord_sort
+  implicit none
+
+contains
+
+  #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+    #:for rank in RANKS
+      #:set name = rname("median_all",rank, t1, k1, o1)
+      module function ${name}$ (x, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        logical, intent(in), optional :: mask
+        real(${o1}$) :: res
+
+        integer(kind = int64) :: c, n
+        ${t1}$, allocatable :: x_tmp(:)
+
+        if (.not.optval(mask, .true.) .or. size(x) == 0) then
+          res = ieee_value(1._${o1}$, ieee_quiet_nan)
+          return
+        end if
+
+        #:if t1[0] == 'r'
+          if (any(ieee_is_nan(x))) then
+            res = ieee_value(1._${o1}$, ieee_quiet_nan)
+            return
+          end if
+        #:endif
+
+        n = size(x, kind=int64)
+        c = floor( (n + 1) / 2._${o1}$, kind=int64 )
+
+        x_tmp = reshape(x, [n])
+
+        call sort(x_tmp)
+
+        if (mod(n, 2_int64) == 0) then
+          #:if t1[0] == 'r'
+            res = sum(x_tmp(c:c+1)) / 2._${o1}$
+          #:else
+            res = sum( real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$
+          #:endif
+        else
+            res = x_tmp(c)
+        end if
+
+      end function ${name}$
+    #:endfor
+  #:endfor
+
+  #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+    #:for rank in RANKS
+      #:set name = rname("median",rank, t1, k1, o1)
+      module function ${name}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$
+
+        integer :: c, n
+        #:if rank > 1
+        #:for fj in range(1, rank+1)
+        integer :: j${fj}$
+        #:endfor
+        #:endif
+        ${t1}$, allocatable :: x_tmp(:)
+
+        if (.not.optval(mask, .true.) .or. size(x) == 0) then
+          res = ieee_value(1._${o1}$, ieee_quiet_nan)
+          return
+        end if
+
+        n = size(x, dim)
+        c = floor( (n + 1) / 2._${o1}$ )
+
+        allocate(x_tmp(n))
+
+        select case(dim)
+          #:for fi in range(1, rank+1)
+          case(${fi}$)
+            ! Loop over every dimension of the array except "dim"
+            #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1))
+              do j${fj}$ = 1, size(x, ${fj}$)
+            #:endfor
+                x_tmp(:) = x${select_subvector('j', rank, fi)}$
+
+                #:if t1[0] == 'r'
+                  if (any(ieee_is_nan(x_tmp))) then
+                    res${reduce_subvector('j', rank, fi)}$ = &
+                      ieee_value(1._${o1}$, ieee_quiet_nan)
+                    #:if fi == 1
+                      return
+                    #:else
+                      cycle
+                    #:endif
+                  end if
+                #:endif
+
+                call sort(x_tmp)
+
+                if (mod(n, 2) == 0) then
+                    res${reduce_subvector('j', rank, fi)}$ = &
+                    #:if t1[0] == 'r'
+                        sum(x_tmp(c:c+1)) / 2._${o1}$
+                    #:else
+                        sum(real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$
+                    #:endif
+                else
+                    res${reduce_subvector('j', rank, fi)}$ = x_tmp(c)
+                end if
+            #:for fj in range(1, rank)
+              end do
+            #:endfor
+          #:endfor
+          case default
+            call error_stop("ERROR (median): wrong dimension")
+        end select
+
+      end function ${name}$
+    #:endfor
+  #:endfor
+
+
+  #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+    #:for rank in RANKS
+      #:set name = rname('median_all_mask',rank, t1, k1, o1)
+      module function ${name}$(x, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        real(${o1}$) :: res
+
+        integer(kind = int64) :: c, n
+        ${t1}$, allocatable   :: x_tmp(:)
+
+        if (any(shape(x) .ne. shape(mask))) then
+            call error_stop("ERROR (median): shapes of x and mask are different")
+        end if
+
+        #:if t1[0] == 'r'
+          if (any(ieee_is_nan(x))) then
+            res = ieee_value(1._${o1}$, ieee_quiet_nan)
+            return
+          end if
+        #:endif
+
+        x_tmp = pack(x, mask)
+
+        call sort(x_tmp)
+
+        n = size(x_tmp, kind=int64)
+        c = floor( (n + 1) / 2._${o1}$, kind=int64)
+
+        if (n == 0) then
+            res = ieee_value(1._${o1}$, ieee_quiet_nan)
+        else if (mod(n, 2_int64) == 0) then
+          #:if t1[0] == 'r'
+            res = sum(x_tmp(c:c+1)) / 2._${o1}$
+          #:else
+            res = sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$
+          #:endif
+        else if (mod(n, 2_int64) == 1) then
+            res = x_tmp(c)
+        end if
+
+      end function ${name}$
+    #:endfor
+  #:endfor
+
+  #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT
+    #:for rank in RANKS
+      #:set name = rname('median_mask',rank, t1, k1, o1)
+      module function ${name}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$
+
+        integer(kind = int64) :: c, n
+        #:if rank > 1
+        #:for fj in range(1, rank+1)
+        integer :: j${fj}$
+        #:endfor
+        #:endif
+        ${t1}$, allocatable :: x_tmp(:)
+
+        if (any(shape(x) .ne. shape(mask))) then
+            call error_stop("ERROR (median): shapes of x and mask are different")
+        end if
+
+        select case(dim)
+          #:for fi in range(1, rank+1)
+          case(${fi}$)
+            ! Loop over every dimension of the array except "dim"
+            #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1))
+              do j${fj}$ = 1, size(x, ${fj}$)
+            #:endfor
+                x_tmp = pack(x${select_subvector('j', rank, fi)}$, &
+                              mask${select_subvector('j', rank, fi)}$)
+
+                #:if t1[0] == 'r'
+                  if (any(ieee_is_nan(x_tmp))) then
+                    res${reduce_subvector('j', rank, fi)}$ = &
+                      ieee_value(1._${o1}$, ieee_quiet_nan)
+                    #:if rank == 1
+                      return
+                    #:else
+                      cycle
+                    #:endif
+                  end if
+                #:endif
+
+                call sort(x_tmp)
+
+                n = size(x_tmp, kind=int64)
+                c = floor( (n + 1) / 2._${o1}$, kind=int64 )
+
+                if (n == 0) then
+                    res${reduce_subvector('j', rank, fi)}$ = &
+                        ieee_value(1._${o1}$, ieee_quiet_nan)
+                else if (mod(n, 2_int64) == 0) then
+                    res${reduce_subvector('j', rank, fi)}$ = &
+                    #:if t1[0] == 'r'
+                        sum(x_tmp(c:c+1)) / 2._${o1}$
+                    #:else
+                        sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$
+                    #:endif
+                else if (mod(n, 2_int64) == 1) then
+                    res${reduce_subvector('j', rank, fi)}$ = x_tmp(c)
+                end if
+
+                deallocate(x_tmp)
+            #:for fj in range(1, rank)
+              end do
+            #:endfor
+          #:endfor
+          case default
+            call error_stop("ERROR (median): wrong dimension")
+        end select
+
+      end function ${name}$
+    #:endfor
+  #:endfor
+
+end submodule
diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt
index 38f9bb84b..600e925e5 100644
--- a/src/tests/stats/CMakeLists.txt
+++ b/src/tests/stats/CMakeLists.txt
@@ -1,6 +1,25 @@
+### Pre-process: .fpp -> .f90 via Fypp
+
+# Create a list of the files to be preprocessed
+set(fppFiles
+    test_median.fypp
+)
+
+# Custom preprocessor flags
+if(DEFINED CMAKE_MAXIMUM_RANK)
+  set(fyppFlags "-DMAXRANK=${CMAKE_MAXIMUM_RANK}")
+elseif(f03rank)
+  set(fyppFlags)
+else()
+  set(fyppFlags "-DVERSION90")
+endif()
+
+fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
+
 ADDTEST(corr)
 ADDTEST(cov)
 ADDTEST(mean)
+ADDTEST(median)
 ADDTEST(moment)
 ADDTEST(rawmoment)
 ADDTEST(var)
diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual
index a0731e8f4..55a21c5f4 100644
--- a/src/tests/stats/Makefile.manual
+++ b/src/tests/stats/Makefile.manual
@@ -1,5 +1,12 @@
+SRCFYPP =\
+        test_median.fypp
 
-PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 \
+SRCGEN = $(SRCFYPP:.fypp=.f90)
+
+$(SRCGEN): %.f90: %.fypp common.fypp
+	fypp $(FYPPFLAGS) $< $@
+
+PROGS_SRC = $(SRCGEN) test_mean.f90 test_moment.f90 test_var.f90 \
             test_distribution_PRNG.f90
 
-include ../Makefile.manual.test.mk
\ No newline at end of file
+include ../Makefile.manual.test.mk
diff --git a/src/tests/stats/common.fypp b/src/tests/stats/common.fypp
new file mode 100644
index 000000000..9ebaac6b9
--- /dev/null
+++ b/src/tests/stats/common.fypp
@@ -0,0 +1,37 @@
+#:mute
+
+#! Real kinds to be considered during templating
+#:set REAL_KINDS = ["sp", "dp", "qp"]
+
+#! Real types to be considered during templating
+#:set REAL_TYPES = ["real({})".format(k) for k in REAL_KINDS]
+
+#! Collected (kind, type) tuples for real types
+#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES))
+
+#! Integer kinds to be considered during templating
+#:set INT_KINDS = ["int8", "int16", "int32", "int64"]
+
+#! Integer types to be considered during templating
+#:set INT_TYPES = ["integer({})".format(k) for k in INT_KINDS]
+
+#! Collected (kind, type) tuples for integer types
+#:set INT_KINDS_TYPES = list(zip(INT_KINDS, INT_TYPES))
+
+
+#! Whether Fortran 90 compatible code should be generated
+#:set VERSION90 = defined('VERSION90')
+
+
+#! Ranks to be generated when templates are created
+#:if not defined('MAXRANK')
+  #:if VERSION90
+    #:set MAXRANK = 7
+  #:else
+    #:set MAXRANK = 15
+  #:endif
+#:endif
+
+
+
+#:endmute
diff --git a/src/tests/stats/test_median.fypp b/src/tests/stats/test_median.fypp
new file mode 100644
index 000000000..09962da27
--- /dev/null
+++ b/src/tests/stats/test_median.fypp
@@ -0,0 +1,234 @@
+#:include "common.fypp"
+#:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES
+program test_median
+    use stdlib_error, only: check
+    use stdlib_kinds, only: sp, dp, qp, &
+        int8, int16, int32, int64
+    use stdlib_stats, only: median, mean
+    use, intrinsic :: ieee_arithmetic, only : ieee_is_nan, ieee_value, ieee_quiet_nan
+    implicit none
+    real(sp), parameter :: sptol = 1000 * epsilon(1._sp)
+    real(dp), parameter :: dptol = 2000 * epsilon(1._dp)
+    real(qp), parameter :: qptol = 2000 * epsilon(1._qp)
+
+    #:for k1,t1 in IR_KINDS_TYPES
+    call test_median_${k1}$()
+    #:endfor
+
+contains
+
+#:for k1,t1 in INT_KINDS_TYPES
+subroutine test_median_${k1}$()
+    ${t1}$, allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:)
+
+    allocate(d0(0))
+    !check just to be sure that the setup of d0 is correct
+    call check(size(d0) == 0, 'test_median_${k1}$: should be of size 0')
+
+    d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20]
+    call check(mod(size(d1), 2) == 0, 'test_median_${k1}$: should be an even number')
+
+    d2 = reshape(d1, [3, 4])
+    call check(mod(size(d2), 2) == 0, 'test_median_${k1}$: should be an even number')
+
+    d3 = reshape(d1, [2, 3, 2])
+    call check(mod(size(d3), 2) == 0, 'test_median_${k1}$: should be an even number')
+
+    !median_all
+    call check( ieee_is_nan(median(d0)), '${k1}$ median(d0)' )
+    call check( ieee_is_nan(median(d0, .false.)), '${k1}$ median(d0, .false.)' )
+    call check( ieee_is_nan(median(d1, .false.)), '${k1}$ median(d1, .false.)' )
+    call check( ieee_is_nan(median(d2, .false.)), '${k1}$ median(d2, .false.)' )
+    call check( ieee_is_nan(median(d3, .false.)), '${k1}$ median(d3, .false.)' )
+
+    call check( abs(median(d1) - 1.5_dp) < dptol, '${k1}$ median(d1), even')
+    call check( abs(median(d2) - 1.5_dp) < dptol, '${k1}$ median(d2), even')
+    call check( abs(median(d3) - 1.5_dp) < dptol, '${k1}$ median(d3), even')
+
+    !median
+    call check( ieee_is_nan(median(d0, 1)), '${k1}$ median(d0, 1)' )
+    call check( ieee_is_nan(median(d0, 1, .false.)), '${k1}$ median(d0, 1, .false.)' )
+    call check( ieee_is_nan(median(d1, 1, .false.)), '${k1}$ median(d1, 1, .false.)' )
+    call check( any(ieee_is_nan(median(d2, 1, .false.))), '${k1}$ median(d2, 1, .false.)' )
+    call check( any(ieee_is_nan(median(d2, 2, .false.))), '${k1}$ median(d2, 2, .false.)' )
+    call check( any(ieee_is_nan(median(d3, 1, .false.))), '${k1}$ median(d3, 1, .false.)' )
+    call check( any(ieee_is_nan(median(d3, 2, .false.))), '${k1}$ median(d3, 2, .false.)' )
+    call check( any(ieee_is_nan(median(d3, 3, .false.))), '${k1}$ median(d3, 3, .false.)' )
+
+    call check( abs(median(d1, 1) - 1.5_dp) < dptol, '${k1}$ median(d1, 1), even')
+    call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp])) < dptol, &
+                 '${k1}$ median(d2, 1)')
+    call check( sum(abs(median(d2, 2) - [3.5_dp, 1.5_dp, 3._dp])) < dptol,       &
+                 '${k1}$ median(d2, 2)')
+
+    !median_mask_all
+    call check( ieee_is_nan(median(d0, d0 > 0)), '${k1}$ median(d0, d0 > 0)' )
+    call check( ieee_is_nan(median(d1, d1 > huge(d1))), '${k1}$ median(d1, d1 > huge(d1))' )
+    call check( ieee_is_nan(median(d2, d2 > huge(d2))), '${k1}$ median(d2, d2 > huge(d2))' )
+    call check( ieee_is_nan(median(d3, d3 > huge(d3))), '${k1}$ median(d3, d3 > huge(d3))' )
+    call check( (median(d1, d1 > 0) - 7._dp) < dptol, '${k1}$ median(d1, d1 > 0)' ) 
+    call check( (median(d2, d2 > 0) - 7._dp) < dptol, '${k1}$ median(d2, d2 > 0)' ) 
+    call check( (median(d3, d3 > 0) - 7._dp) < dptol, '${k1}$ median(d3, d3 > 0)' ) 
+
+    !median mask
+    call check( ieee_is_nan(median(d0, 1, d0 > 0)), '${k1}$ median(d0, 1, d0 > 0)' )
+    call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), '${k1}$ median(d1, 1, d1 > huge(d1))' )
+    call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), '${k1}$ median(d2, 1, d2 > huge(d2))' )
+    call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), '${k1}$ median(d3, 1, d3 > huge(d3))' )
+    call check( (median(d1, 1, d1 > 0) - 7._dp) < dptol, '${k1}$ median(d1, 1, d1 >0)') 
+    call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._dp, 6._dp, 8._dp, 10.5_dp] )  )) &
+                 < dptol, '${k1}$ median(d2, 1, d2 > 0 )') 
+    call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_dp, 2._dp, 14.5_dp] )  ))       &
+                 < dptol, '${k1}$ median(d2, 2, d2 > 0 )') 
+    call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), '${k1}$ median(d3, 1, d3 > 0)')
+
+    !odd number
+    d1 = [d1, 20_${k1}$]
+    call check(mod(size(d1), 2) == 1, 'test_median_${k1}$: should be an odd number')
+
+    call check( abs(median(d1) - 2._dp) < dptol, '${k1}$ median(d1), odd')
+
+    d2 = reshape(d1, [3, 5], pad = [0_${k1}$])
+    call check(mod(size(d2), 2) == 1, 'test_median_${k1}$: should be an odd number')
+    call check( abs(median(d2) - 1._dp) < dptol, '${k1}$ median(d2), odd')
+
+    d3 = reshape(d1, [1, 3, 5], pad = [0_${k1}$])
+    call check(mod(size(d3), 2) == 1, 'test_median_${k1}$: should be an odd number')
+    call check( abs(median(d3) - 1._dp) < dptol, '${k1}$ median(d3), odd')
+
+    call check( abs(median(d1, 1) - 2._dp) < dptol, '${k1}$ median(d1, 1)')
+    call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp, 0._dp])) < dptol, &
+                 '${k1}$ median(d2, 1), odd')
+    call check( sum(abs(median(d2, 2) - [7._dp, 1._dp, 0._dp])) < dptol,                &
+                 '${k1}$ median(d2, 2), odd')
+
+end subroutine
+#:endfor
+
+#:for k1,t1 in REAL_KINDS_TYPES
+subroutine test_median_${k1}$()
+    ${t1}$, allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:)
+    ${t1}$, allocatable :: tmp1(:)
+
+    allocate(d0(0))
+    !check just to be sure that the setup of d0 is correct
+    call check(size(d0) == 0, 'test_median_${k1}$: should be of size 0')
+
+    d1 = [10., 2., -3., -4., 6., -6., 7., -8., 9., 0., 1., 20.]
+    call check(mod(size(d1), 2) == 0, 'test_median_${k1}$: should be an even number')
+
+    d2 = reshape(d1, [3, 4])
+    call check(mod(size(d2), 2) == 0, 'test_median_${k1}$: should be an even number')
+
+    d3 = reshape(d1, [2, 3, 2])
+    call check(mod(size(d3), 2) == 0, 'test_median_${k1}$: should be an even number')
+
+    !median_all
+    call check( ieee_is_nan(median(d0)), '${k1}$ median(d0)' )
+    call check( ieee_is_nan(median(d0, .false.)), '${k1}$ median(d0, .false.)' )
+    call check( ieee_is_nan(median(d1, .false.)), '${k1}$ median(d1, .false.)' )
+    call check( ieee_is_nan(median(d2, .false.)), '${k1}$ median(d2, .false.)' )
+    call check( ieee_is_nan(median(d3, .false.)), '${k1}$ median(d3, .false.)' )
+
+    call check( abs(median(d1) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d1), even')
+    call check( abs(median(d2) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d2), even')
+    call check( abs(median(d3) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d3), even')
+
+    !median
+    call check( ieee_is_nan(median(d0, 1)), '${k1}$ median(d0, 1)' )
+    call check( ieee_is_nan(median(d0, 1, .false.)), '${k1}$ median(d0, 1, .false.)' )
+    call check( ieee_is_nan(median(d1, 1, .false.)), '${k1}$ median(d1, 1, .false.)' )
+    call check( any(ieee_is_nan(median(d2, 1, .false.))), '${k1}$ median(d2, 1, .false.)' )
+    call check( any(ieee_is_nan(median(d2, 2, .false.))), '${k1}$ median(d2, 2, .false.)' )
+    call check( any(ieee_is_nan(median(d3, 1, .false.))), '${k1}$ median(d3, 1, .false.)' )
+    call check( any(ieee_is_nan(median(d3, 2, .false.))), '${k1}$ median(d3, 2, .false.)' )
+    call check( any(ieee_is_nan(median(d3, 3, .false.))), '${k1}$ median(d3, 3, .false.)' )
+
+    call check( abs(median(d1, 1) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d1, 1), even')
+    call check( sum(abs(median(d2, 1) - [2._${k1}$, -4._${k1}$, 7._${k1}$, 1._${k1}$])) < ${k1}$tol, &
+                 '${k1}$ median(d2, 1)')
+    call check( sum(abs(median(d2, 2) - [3.5_${k1}$, 1.5_${k1}$, 3._${k1}$])) < ${k1}$tol,       &
+                 '${k1}$ median(d2, 2)')
+
+    !median_mask_all
+    call check( ieee_is_nan(median(d0, d0 > 0)), '${k1}$ median(d0, d0 > 0)' )
+    call check( ieee_is_nan(median(d1, d1 > huge(d1))), '${k1}$ median(d1, d1 > huge(d1))' )
+    call check( ieee_is_nan(median(d2, d2 > huge(d2))), '${k1}$ median(d2, d2 > huge(d2))' )
+    call check( ieee_is_nan(median(d3, d3 > huge(d3))), '${k1}$ median(d3, d3 > huge(d3))' )
+    call check( (median(d1, d1 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d1, d1 > 0)' ) 
+    call check( (median(d2, d2 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d2, d2 > 0)' ) 
+    call check( (median(d3, d3 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d3, d3 > 0)' ) 
+
+    !median mask
+    call check( ieee_is_nan(median(d0, 1, d0 > 0)), '${k1}$ median(d0, 1, d0 > 0)' )
+    call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), '${k1}$ median(d1, 1, d1 > huge(d1))' )
+    call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), '${k1}$ median(d2, 1, d2 > huge(d2))' )
+    call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), '${k1}$ median(d3, 1, d3 > huge(d3))' )
+    call check( (median(d1, 1, d1 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d1, 1, d1 >0)') 
+    call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._${k1}$, 6._${k1}$, 8._${k1}$, 10.5_${k1}$] )  )) &
+                 < ${k1}$tol, '${k1}$ median(d2, 1, d2 > 0 )') 
+    call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_${k1}$, 2._${k1}$, 14.5_${k1}$] )  ))       &
+                 < ${k1}$tol, '${k1}$ median(d2, 2, d2 > 0 )') 
+    call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), '${k1}$ median(d3, 1, d3 > 0)')
+
+    !odd number
+    d1 = [d1, 20._${k1}$]
+    call check(mod(size(d1), 2) == 1, 'test_median_${k1}$: should be an odd number')
+
+    call check( abs(median(d1) - 2._${k1}$) < ${k1}$tol, '${k1}$ median(d1), odd')
+
+    d2 = reshape(d1, [3, 5], pad = [0._${k1}$])
+    call check(mod(size(d2), 2) == 1, 'test_median_${k1}$: should be an odd number')
+    call check( abs(median(d2) - 1._${k1}$) < ${k1}$tol, '${k1}$ median(d2), odd')
+
+    d3 = reshape(d1, [1, 3, 5], pad = [0._${k1}$])
+    call check(mod(size(d3), 2) == 1, 'test_median_${k1}$: should be an odd number')
+    call check( abs(median(d3) - 1._${k1}$) < ${k1}$tol, '${k1}$ median(d3), odd')
+
+    call check( abs(median(d1, 1) - 2._${k1}$) < ${k1}$tol, '${k1}$ median(d1, 1)')
+    call check( sum(abs(median(d2, 1) - [2._${k1}$, -4._${k1}$, 7._${k1}$, 1._${k1}$, 0._${k1}$])) < ${k1}$tol, &
+                 '${k1}$ median(d2, 1), odd')
+    call check( sum(abs(median(d2, 2) - [7._${k1}$, 1._${k1}$, 0._${k1}$])) < ${k1}$tol,                &
+                 '${k1}$ median(d2, 2), odd')
+
+    !check IEEE NaN values in array
+    d1(1) = ieee_value(1._${k1}$, ieee_quiet_nan)
+    d2(1, 1) =  ieee_value(1._${k1}$, ieee_quiet_nan)
+    d3(1, 1, 1) =  ieee_value(1._${k1}$, ieee_quiet_nan)
+
+    !median_all
+    call check( ieee_is_nan(median(d1)), '${k1}$ median(d1), should be NaN')
+    call check( ieee_is_nan(median(d2)), '${k1}$ median(d2), should be NaN')
+    call check( ieee_is_nan(median(d3)), '${k1}$ median(d3), should be NaN')
+
+    !median
+    call check( any(ieee_is_nan(median(d2, 1))), '${k1}$ median(d2, 1) should contain at least 1 NaN' )
+    call check( any(ieee_is_nan(median(d2, 2))), '${k1}$ median(d2, 2) should contain at least 1 NaN' )
+
+    call check( any(ieee_is_nan(median(d3, 1))), '${k1}$ median(d3, 1) should contain at least 1 NaN' )
+    call check( any(ieee_is_nan(median(d3, 2))), '${k1}$ median(d3, 2) should contain at least 1 NaN' )
+    call check( any(ieee_is_nan(median(d3, 3))), '${k1}$ median(d3, 3) should contain at least 1 NaN' )
+
+    !median_mask_all
+    call check( ieee_is_nan(median(d1, d1 > 0)), '${k1}$ median(d1, d1 > 0) should be NaN' )
+    call check( ieee_is_nan(median(d2, d2 > 0)), '${k1}$ median(d2, d2 > 0) should be NaN' )
+    call check( ieee_is_nan(median(d3, d3 > 0)), '${k1}$ median(d3, d3 > 0) should be NaN' )
+
+    !median mask
+    call check( ieee_is_nan(median(d1, 1, d1 .ne. 0)), '${k1}$ median(d1, 1, d1.ne.0) should return NaN')
+
+    tmp1 = median(d2, 1, d2 .ne. 0)
+    call check( any(ieee_is_nan(tmp1)), '${k1}$ median(d2, 1, d2 .ne. 0 ) should contain at least 1 NaN')
+    call check( tmp1(2) == -4._${k1}$, '${k1}$ tmp1(2) == -4')
+    call check( any(ieee_is_nan(median(d2, 2, d2 .ne. 0))),              &
+       '${k1}$ median(d2, 2, d2 .ne. 0 ) should contain at least 1 NaN')
+
+    call check( any(ieee_is_nan(median(d3, 1, d3 .ne. 0))),              &
+       '${k1}$ median(d3, 1, d3 .ne. 0 ) should contain at least 1 NaN')
+
+end subroutine
+#:endfor
+
+
+
+end program