diff --git a/src/stdlib_experimental_stat.f90 b/src/stdlib_experimental_stat.f90
index 73b4bf919..64cce049a 100644
--- a/src/stdlib_experimental_stat.f90
+++ b/src/stdlib_experimental_stat.f90
@@ -42,6 +42,7 @@ module function mean_1_int64_dp(x) result(res)
         real(dp) :: res
     end function mean_1_int64_dp
 
+
     module function mean_2_all_sp_sp(x) result(res)
         real(sp), intent(in) :: x(:,:)
         real(sp) :: res
@@ -109,6 +110,95 @@ module function mean_2_int64_dp(x, dim) result(res)
         real(dp) :: res(size(x)/size(x, dim))
     end function mean_2_int64_dp
 
+
+module function mean_3_all_sp_sp(x) result(res)
+    real(sp), intent(in) :: x(:,:,:)
+    real(sp) :: res
+end function mean_3_all_sp_sp
+module function mean_3_all_dp_dp(x) result(res)
+    real(dp), intent(in) :: x(:,:,:)
+    real(dp) :: res
+end function mean_3_all_dp_dp
+module function mean_3_all_qp_qp(x) result(res)
+    real(qp), intent(in) :: x(:,:,:)
+    real(qp) :: res
+end function mean_3_all_qp_qp
+
+module function mean_3_all_int8_dp(x) result(res)
+    integer(int8), intent(in) :: x(:,:,:)
+    real(dp) :: res
+end function mean_3_all_int8_dp
+module function mean_3_all_int16_dp(x) result(res)
+    integer(int16), intent(in) :: x(:,:,:)
+    real(dp) :: res
+end function mean_3_all_int16_dp
+module function mean_3_all_int32_dp(x) result(res)
+    integer(int32), intent(in) :: x(:,:,:)
+    real(dp) :: res
+end function mean_3_all_int32_dp
+module function mean_3_all_int64_dp(x) result(res)
+    integer(int64), intent(in) :: x(:,:,:)
+    real(dp) :: res
+end function mean_3_all_int64_dp
+
+module function mean_3_sp_sp(x, dim) result(res)
+    real(sp), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(sp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_sp_sp
+module function mean_3_dp_dp(x, dim) result(res)
+    real(dp), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_dp_dp
+module function mean_3_qp_qp(x, dim) result(res)
+    real(qp), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(qp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_qp_qp
+
+module function mean_3_int8_dp(x, dim) result(res)
+    integer(int8), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_int8_dp
+module function mean_3_int16_dp(x, dim) result(res)
+    integer(int16), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_int16_dp
+module function mean_3_int32_dp(x, dim) result(res)
+    integer(int32), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_int32_dp
+module function mean_3_int64_dp(x, dim) result(res)
+    integer(int64), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_int64_dp
+
 end interface
 
 end module
diff --git a/src/stdlib_experimental_stat.fypp.f90 b/src/stdlib_experimental_stat.fypp.f90
index 42aa5440d..8999eb5c5 100644
--- a/src/stdlib_experimental_stat.fypp.f90
+++ b/src/stdlib_experimental_stat.fypp.f90
@@ -32,6 +32,7 @@ module function mean_1_${k1}$_dp(x) result(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(:,:)
@@ -62,6 +63,43 @@ module function mean_2_${k1}$_dp(x, dim) result(res)
     end function mean_2_${k1}$_dp
 #:endfor
 
+
+#:for i1, k1, t1 in iktr
+module function mean_3_all_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    ${t1}$ :: res
+end function mean_3_all_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_3_all_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    real(dp) :: res
+end function mean_3_all_${k1}$_dp
+#:endfor
+
+#:for i1, k1, t1 in iktr
+module function mean_3_${k1}$_${k1}$(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    ${t1}$ :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_3_${k1}$_dp(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+end function mean_3_${k1}$_dp
+#:endfor
+
 end interface
 
 end module
diff --git a/src/stdlib_experimental_stat_mean.f90 b/src/stdlib_experimental_stat_mean.f90
index be6b5733a..578c470d1 100644
--- a/src/stdlib_experimental_stat_mean.f90
+++ b/src/stdlib_experimental_stat_mean.f90
@@ -61,7 +61,6 @@ module function mean_1_int64_dp(x) result(res)
 end function mean_1_int64_dp
 
 
-
 module function mean_2_all_sp_sp(x) result(res)
     real(sp), intent(in) :: x(:,:)
     real(sp) :: res
@@ -113,8 +112,6 @@ module function mean_2_all_int64_dp(x) result(res)
 
 end function mean_2_all_int64_dp
 
-
-
 module function mean_2_sp_sp(x, dim) result(res)
     real(sp), intent(in) :: x(:,:)
     integer, intent(in) :: dim
@@ -244,4 +241,288 @@ module function mean_2_int64_dp(x, dim) result(res)
 end function mean_2_int64_dp
 
 
+module function mean_3_all_sp_sp(x) result(res)
+    real(sp), intent(in) :: x(:,:,:)
+    real(sp) :: res
+
+    res = sum(x) / real(size(x), sp)
+
+end function mean_3_all_sp_sp
+module function mean_3_all_dp_dp(x) result(res)
+    real(dp), intent(in) :: x(:,:,:)
+    real(dp) :: res
+
+    res = sum(x) / real(size(x), dp)
+
+end function mean_3_all_dp_dp
+module function mean_3_all_qp_qp(x) result(res)
+    real(qp), intent(in) :: x(:,:,:)
+    real(qp) :: res
+
+    res = sum(x) / real(size(x), qp)
+
+end function mean_3_all_qp_qp
+
+module function mean_3_all_int8_dp(x) result(res)
+    integer(int8), intent(in) :: x(:,:,:)
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x), dp)
+
+end function mean_3_all_int8_dp
+module function mean_3_all_int16_dp(x) result(res)
+    integer(int16), intent(in) :: x(:,:,:)
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x), dp)
+
+end function mean_3_all_int16_dp
+module function mean_3_all_int32_dp(x) result(res)
+    integer(int32), intent(in) :: x(:,:,:)
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x), dp)
+
+end function mean_3_all_int32_dp
+module function mean_3_all_int64_dp(x) result(res)
+    integer(int64), intent(in) :: x(:,:,:)
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x), dp)
+
+end function mean_3_all_int64_dp
+
+module function mean_3_sp_sp(x, dim) result(res)
+    real(sp), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(sp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_sp_sp(x(:,i,j))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_sp_sp(x(i,:,j))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_sp_sp(x(i,j,:))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_sp_sp
+module function mean_3_dp_dp(x, dim) result(res)
+    real(dp), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(x(:,i,j))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(x(i,:,j))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(x(i,j,:))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_dp_dp
+module function mean_3_qp_qp(x, dim) result(res)
+    real(qp), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(qp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_qp_qp(x(:,i,j))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_qp_qp(x(i,:,j))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_qp_qp(x(i,j,:))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_qp_qp
+
+module function mean_3_int8_dp(x, dim) result(res)
+    integer(int8), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(:,i,j), dp))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,:,j), dp))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,j,:), dp))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_int8_dp
+module function mean_3_int16_dp(x, dim) result(res)
+    integer(int16), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(:,i,j), dp))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,:,j), dp))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,j,:), dp))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_int16_dp
+module function mean_3_int32_dp(x, dim) result(res)
+    integer(int32), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(:,i,j), dp))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,:,j), dp))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,j,:), dp))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_int32_dp
+module function mean_3_int64_dp(x, dim) result(res)
+    integer(int64), intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(:,i,j), dp))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,:,j), dp))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,j,:), dp))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_int64_dp
+
 end submodule
diff --git a/src/stdlib_experimental_stat_mean.fypp.f90 b/src/stdlib_experimental_stat_mean.fypp.f90
index 9f7b61ea2..9c7205e3e 100644
--- a/src/stdlib_experimental_stat_mean.fypp.f90
+++ b/src/stdlib_experimental_stat_mean.fypp.f90
@@ -36,7 +36,6 @@ 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(:,:)
@@ -57,8 +56,6 @@ module function mean_2_all_${k1}$_dp(x) result(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(:,:)
@@ -102,4 +99,96 @@ end function mean_2_${k1}$_dp
 #:endfor
 
 
+#:for i1, k1, t1 in iktr
+module function mean_3_all_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    ${t1}$ :: res
+
+    res = sum(x) / real(size(x), ${k1}$)
+
+end function mean_3_all_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_3_all_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    real(dp) :: res
+
+    res = sum(real(x, dp)) / real(size(x), dp)
+
+end function mean_3_all_${k1}$_dp
+#:endfor
+
+#:for i1, k1, t1 in iktr
+module function mean_3_${k1}$_${k1}$(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    ${t1}$ :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_${k1}$_${k1}$(x(:,i,j))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_${k1}$_${k1}$(x(i,:,j))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_${k1}$_${k1}$(x(i,j,:))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_3_${k1}$_dp(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:,:)
+    integer, intent(in) :: dim
+    integer :: j_
+    real(dp) :: res( &
+                  merge(size(x,2),size(x,1),mask = any((/(j_, j_ = dim, 3)/) == 1)), &
+                  merge(size(x,3),size(x,2),mask = any((/(j_, j_ = dim, 3)/) == 2)) )
+
+    integer :: i, j
+
+    select case(dim)
+     case(1)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(:,i,j), dp))
+            end do
+        end do
+     case(2)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,:,j), dp))
+            end do
+        end do
+     case(3)
+        do i=1, size(res, 1)
+            do j=1, size(res, 2)
+                res(i, j) = mean_1_dp_dp(real(x(i,j,:), dp))
+            end do
+        end do
+     case default
+    end select
+
+end function mean_3_${k1}$_dp
+#:endfor
+
 end submodule
diff --git a/src/tests/stat/test_mean.f90 b/src/tests/stat/test_mean.f90
index a4a4f832a..f50129b3e 100644
--- a/src/tests/stat/test_mean.f90
+++ b/src/tests/stat/test_mean.f90
@@ -8,6 +8,8 @@ program test_mean
 real(sp), allocatable :: s(:, :)
 real(dp), allocatable :: d(:, :)
 
+real(dp), allocatable :: d3(:, :, :)
+
 !sp
 call loadtxt("array1.dat", s)
 
@@ -33,6 +35,33 @@ program test_mean
 call assert(sum( mean(int(s, int32), dim = 2) - [1.5_dp, 3.5_dp, 5.5_dp, 7.5_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( sum( shape(mean(d3,1))-shape(sum(d3,1)) ) == 0)
+call assert( sum( shape(mean(d3,2))-shape(sum(d3,2)) ) == 0)
+call assert( sum( shape(mean(d3,3))-shape(sum(d3,3)) ) == 0)
+
+call assert( mean(d3) - sum(d3)/size(d3) == 0.0_dp)
+call assert( sum(abs( mean(d3,1) - &  
+             reshape([4.0_dp, 5.0_dp, 6.0_dp, 7.5_dp, 16.0_dp, 20.0_dp], shape(sum(d3,1))) ) ) &
+             == 0.0_dp)
+call assert( sum(abs( mean(d3,2) - &  
+             reshape([ 1.5_dp, 3.5_dp, 5.5_dp, 7.5_dp, &
+             2.25_dp, 5.25_dp, 8.25_dp, 11.25_dp, &
+             6._dp, 14._dp, 22._dp, 30._dp &
+             ], shape(sum(d3,2))) ) ) &
+             == 0.0_dp)
+call assert( sum(abs( mean(d3,3) - &  
+             reshape([2.1666666666666665_dp, 6.5_dp, 10.833333333333334_dp,&
+             15.166666666666666_dp, 4.333333333333333_dp, &
+             8.6666666666666661_dp, 13.0_dp, 17.333333333333332_dp &
+             ], shape(sum(d3,3))) ) ) &
+             == 0.0_dp)
+
 contains
 
 end program