diff --git a/monin_obukhov/monin_obukhov.F90 b/monin_obukhov/monin_obukhov.F90
index 883e4cbe34..ac8a89075f 100644
--- a/monin_obukhov/monin_obukhov.F90
+++ b/monin_obukhov/monin_obukhov.F90
@@ -274,16 +274,18 @@ subroutine stable_mix_3d(rich, mix)
 
 real, intent(in) , dimension(:,:,:)  :: rich
 real, intent(out), dimension(:,:,:)  :: mix
+integer :: n2 !< Size of dimension 2 of mix and rich
+integer :: n3 !< Size of dimension 3 of mix and rich
+integer :: i, j !< Loop indices
 
-integer :: n, ier
-
-if(.not.module_is_initialized) call error_mesg('stable_mix_3d in monin_obukhov_mod', &
-     'monin_obukhov_init has not been called', FATAL)
-
-n = size(rich,1)*size(rich,2)*size(rich,3)
-call monin_obukhov_stable_mix(stable_option, rich_crit, zeta_trans, &
-     & n, rich, mix, ier)
+n2 = size(mix, 2)
+n3 = size(mix, 3)
 
+do j=1, n3
+  do i=1, n2
+    call stable_mix(rich(:, i, j), mix(:, i, j))
+  enddo
+enddo
 
 end subroutine stable_mix_3d
 
@@ -943,16 +945,15 @@ subroutine stable_mix_2d(rich, mix)
 
 real, intent(in) , dimension(:,:)  :: rich
 real, intent(out), dimension(:,:)  :: mix
+integer :: n2 !< Size of dimension 2 of mix and rich
+integer :: i !< Loop index
 
-real, dimension(size(rich,1),size(rich,2),1) :: rich_3d, mix_3d
-
-rich_3d(:,:,1) = rich
+n2 = size(mix, 2)
 
-call stable_mix_3d(rich_3d, mix_3d)
-
-mix = mix_3d(:,:,1)
+do i=1, n2
+  call stable_mix(rich(:, i), mix(:, i))
+enddo
 
-return
 end subroutine stable_mix_2d
 
 
@@ -962,16 +963,17 @@ subroutine stable_mix_1d(rich, mix)
 
 real, intent(in) , dimension(:)  :: rich
 real, intent(out), dimension(:)  :: mix
+integer :: n !< Size of mix and rich
+integer :: ierr !< Error code set by monin_obukhov_stable_mix
 
-real, dimension(size(rich),1,1) :: rich_3d, mix_3d
-
-rich_3d(:,1,1) = rich
+if (.not.module_is_initialized) call error_mesg('stable_mix in monin_obukhov_mod', &
+     'monin_obukhov_init has not been called', FATAL)
 
-call stable_mix_3d(rich_3d, mix_3d)
+n = size(mix)
 
-mix = mix_3d(:,1,1)
+call monin_obukhov_stable_mix(stable_option, rich_crit, zeta_trans, &
+     & n, rich, mix, ierr)
 
-return
 end subroutine stable_mix_1d
 
 !=======================================================================
@@ -981,15 +983,12 @@ subroutine stable_mix_0d(rich, mix)
 real, intent(in) :: rich
 real, intent(out) :: mix
 
-real, dimension(1,1,1) :: rich_3d, mix_3d
-
-rich_3d(1,1,1) = rich
+real, dimension(1) :: mix_1d !< Representation of mix as a dimension(1) array
 
-call stable_mix_3d(rich_3d, mix_3d)
+call stable_mix([rich], mix_1d)
 
-mix = mix_3d(1,1,1)
+mix = mix_1d(1)
 
-return
 end subroutine stable_mix_0d
 !=======================================================================