diff --git a/CMakeLists.txt b/CMakeLists.txt
index dd7dd2beed83fe5a263e6c9e7347a76860700274..5406c2ad6d94d2b54f2821d44a533bc9924e7419 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -23,6 +23,7 @@ option(FS_ENABLE_BACKTRACE_TEST "Test backtrace function" ON)
 option(FS_ENABLE_OMP "Build with OpenMP support" OFF)
 option(FS_ENABLE_OPENACC "Build with OpenACC support" OFF)
 option(FS_ENABLE_MIXED_PRECISION "Use mixed precision" OFF)
+option(FS_ENABLE_SINGLE_PRECISION "Use single precision" OFF)
 
 include(GNUInstallDirs)
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index f25f2227e33eb52b999425406aea92c58370f835..3162cb3ebd1b86c01f439e897d6c8a686f4603d4 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -27,6 +27,7 @@ add_library(
   fortran-support
   ${CMAKE_CURRENT_BINARY_DIR}/config.h
   fortran_support.F90
+  mo_iconlib_kind.F90
   mo_exception.F90
   mo_expression.F90
   mo_fortran_tools.F90
@@ -75,7 +76,15 @@ set_target_properties(
              EXPORT_NAME ${PROJECT_NAME}::fortran-support)
 
 if(FS_ENABLE_MIXED_PRECISION)
-  target_compile_definitions(fortran-support PRIVATE __MIXED_PRECISION)
+  if(FS_ENABLE_SINGLE_PRECISION)
+    # Do not enable MIXED_PRECISION if SINGLE_PRECISION enabled
+  else()
+    target_compile_definitions(fortran-support PRIVATE __MIXED_PRECISION)
+  endif()
+endif()
+
+if(FS_ENABLE_SINGLE_PRECISION)
+  target_compile_definitions(fortran-support PRIVATE __SINGLE_PRECISION)
 endif()
 
 include("${PROJECT_SOURCE_DIR}/cmake/check_macro.cmake")
diff --git a/src/fortran_support.F90 b/src/fortran_support.F90
index 932bd0af4d11ba16e4c86aa66e51fb01c15a59bc..04a570fef0c6b031b96a0861d06aa4766955c690 100644
--- a/src/fortran_support.F90
+++ b/src/fortran_support.F90
@@ -16,16 +16,41 @@ MODULE fortran_support
     & set_msg_timestamp, enable_logging, message, warning, finish, &
     & print_value, message_to_own_unit, message_text
   USE mo_expression, ONLY: expression, parse_expression_string
-  USE mo_fortran_tools, ONLY: assign_if_present, t_ptr_2d3d, t_ptr_2d3d_vp, &
-    & assign_if_present_allocatable, if_associated, t_ptr_1d, t_ptr_1d_sp, &
-    & t_ptr_1d_int, t_ptr_1d_ptr_1d, t_ptr_2d, t_ptr_2d_sp, t_ptr_2d_int, &
-    & t_ptr_3d, t_ptr_3d_sp, t_ptr_3d_int, t_ptr_i2d3d, t_ptr_4d, t_ptr_4d_sp, &
-    & t_ptr_4d_int, t_ptr_tracer, copy, init, swap, negative2zero, var_scale, &
-    & var_add, init_zero_contiguous_dp, init_zero_contiguous_sp, &
+  USE mo_fortran_tools, ONLY: assign_if_present, assign_if_present_allocatable, if_associated, &
+    & t_ptr_1d_dp, t_ptr_1d_sp, t_ptr_1d_int, &
+    & t_ptr_2d_dp, t_ptr_2d_sp, t_ptr_2d_int, &
+    & t_ptr_3d_dp, t_ptr_3d_sp, t_ptr_3d_int, t_ptr_i2d3d, &
+    & t_ptr_4d_dp, t_ptr_4d_sp, t_ptr_4d_int, &
+    & t_ptr_2d3d, t_ptr_2d3d_vp, t_ptr_tracer, &
+    & copy, init, swap, negative2zero, var_scale, &
+    & var_add, init_zero_contiguous, init_contiguous, &
+    & init_zero_contiguous_dp, init_zero_contiguous_sp, &
     & init_contiguous_dp, init_contiguous_sp, init_contiguous_i4, &
     & init_contiguous_l, minval_1d, minval_2d, resize_arr_c1d, DO_DEALLOCATE, &
     & DO_PTR_DEALLOCATE, insert_dimension, assert_acc_host_only, &
     & assert_acc_device_only, set_acc_host_or_device
+#ifdef __SINGLE_PRECISION
+  USE mo_fortran_tools, ONLY: &
+    & t_ptr_1d_wp => t_ptr_1d_sp, &
+    & t_ptr_2d_wp => t_ptr_2d_sp, &
+    & t_ptr_3d_wp => t_ptr_3d_sp, &
+    & t_ptr_4d_wp => t_ptr_4d_sp, &
+    & t_ptr_1d => t_ptr_1d_sp, &
+    & t_ptr_2d => t_ptr_2d_sp, &
+    & t_ptr_3d => t_ptr_3d_sp, &
+    & t_ptr_4d => t_ptr_4d_sp
+#else
+  USE mo_fortran_tools, ONLY: &
+    & t_ptr_1d_wp => t_ptr_1d_dp, &
+    & t_ptr_2d_wp => t_ptr_2d_dp, &
+    & t_ptr_3d_wp => t_ptr_3d_dp, &
+    & t_ptr_4d_wp => t_ptr_4d_dp, &
+    & t_ptr_1d => t_ptr_1d_dp, &
+    & t_ptr_2d => t_ptr_2d_dp, &
+    & t_ptr_3d => t_ptr_3d_dp, &
+    & t_ptr_4d => t_ptr_4d_dp
+#endif
+
   USE mo_hash_table, ONLY: t_HashTable, hashTable_make, t_HashIterator
   USE mo_io_units, ONLY: filename_max, nerr, nlog, nnml, nstat, ngmt, nin, &
     & nout, nnml_output, find_next_free_unit
@@ -82,12 +107,15 @@ MODULE fortran_support
   PUBLIC :: expression, parse_expression_string
 
   ! From mo_fortran_tools
-  PUBLIC :: assign_if_present, t_ptr_2d3d, t_ptr_2d3d_vp, &
-    & assign_if_present_allocatable, if_associated, t_ptr_1d, t_ptr_1d_sp, &
-    & t_ptr_1d_int, t_ptr_1d_ptr_1d, t_ptr_2d, t_ptr_2d_sp, t_ptr_2d_int, &
-    & t_ptr_3d, t_ptr_3d_sp, t_ptr_3d_int, t_ptr_i2d3d, t_ptr_4d, t_ptr_4d_sp, &
-    & t_ptr_4d_int, t_ptr_tracer, copy, init, swap, negative2zero, var_scale, &
-    & var_add, init_zero_contiguous_dp, init_zero_contiguous_sp, &
+  PUBLIC :: assign_if_present, assign_if_present_allocatable, if_associated, &
+    & t_ptr_1d, t_ptr_1d_wp, t_ptr_1d_dp, t_ptr_1d_sp, t_ptr_1d_int, &
+    & t_ptr_2d, t_ptr_2d_wp, t_ptr_2d_dp, t_ptr_2d_sp, t_ptr_2d_int, &
+    & t_ptr_3d, t_ptr_3d_wp, t_ptr_3d_dp, t_ptr_3d_sp, t_ptr_3d_int, t_ptr_i2d3d, &
+    & t_ptr_4d, t_ptr_4d_wp, t_ptr_4d_dp, t_ptr_4d_sp, t_ptr_4d_int, &
+    & t_ptr_2d3d, t_ptr_2d3d_vp, t_ptr_tracer, &
+    & copy, init, swap, negative2zero, var_scale, &
+    & var_add, init_zero_contiguous, init_contiguous, &
+    & init_zero_contiguous_dp, init_zero_contiguous_sp, &
     & init_contiguous_dp, init_contiguous_sp, init_contiguous_i4, &
     & init_contiguous_l, minval_1d, minval_2d, resize_arr_c1d, DO_DEALLOCATE, &
     & DO_PTR_DEALLOCATE, insert_dimension, assert_acc_host_only, &
diff --git a/src/mo_exception.F90 b/src/mo_exception.F90
index a5574a18c6b5b3e9836cfeb92e1a6bd02f8cc3c4..7f6df9cf646d269ea1859025ee92f8b1de25135d 100644
--- a/src/mo_exception.F90
+++ b/src/mo_exception.F90
@@ -38,8 +38,7 @@ MODULE mo_exception
   !
   !       abort: routine to end the code in case finish is called
 
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64, &
-    &                                      i8 => int64
+  USE mo_iconlib_kind, ONLY: wp, dp, sp, i8
   USE mo_io_units, ONLY: filename_max
   USE mo_util_backtrace, ONLY: ftn_util_backtrace
 
@@ -72,7 +71,8 @@ MODULE mo_exception
     MODULE PROCEDURE print_lvalue !< logical
     MODULE PROCEDURE print_ivalue !< integer
     MODULE PROCEDURE print_i8value !< integer(i8)
-    MODULE PROCEDURE print_rvalue !< real
+    MODULE PROCEDURE print_rvalue_dp !< real64
+    MODULE PROCEDURE print_rvalue_sp !< real32
   END INTERFACE
 
   INTERFACE
@@ -398,10 +398,10 @@ CONTAINS
 
   END SUBROUTINE print_i8value
 
-  SUBROUTINE print_rvalue(mstring, rvalue, routine)
+  SUBROUTINE print_rvalue_dp(mstring, rvalue, routine)
 
     CHARACTER(len=*), INTENT(IN)   :: mstring
-    REAL(wp), INTENT(IN)           :: rvalue
+    REAL(dp), INTENT(IN)           :: rvalue
     CHARACTER(len=*), TARGET, OPTIONAL, INTENT(IN) :: routine
     CHARACTER(len=:), POINTER :: rtn
     CHARACTER(len=1), TARGET :: dummy
@@ -416,6 +416,26 @@ CONTAINS
     WRITE (tmp, '(a60,1x,":",g12.5)') mstring, rvalue
     CALL param(rtn, tmp)
 
-  END SUBROUTINE print_rvalue
+  END SUBROUTINE print_rvalue_dp
+
+  SUBROUTINE print_rvalue_sp(mstring, rvalue, routine)
+
+    CHARACTER(len=*), INTENT(IN)   :: mstring
+    REAL(sp), INTENT(IN)           :: rvalue
+    CHARACTER(len=*), TARGET, OPTIONAL, INTENT(IN) :: routine
+    CHARACTER(len=:), POINTER :: rtn
+    CHARACTER(len=1), TARGET :: dummy
+    CHARACTER(len=filename_max) :: tmp = ""
+
+    IF (PRESENT(routine)) THEN
+      rtn => routine
+    ELSE
+      dummy = ' '
+      rtn => dummy(1:0)
+    END IF
+    WRITE (tmp, '(a60,1x,":",g12.5)') mstring, rvalue
+    CALL param(rtn, tmp)
+
+  END SUBROUTINE print_rvalue_sp
 
 END MODULE mo_exception
diff --git a/src/mo_expression.F90 b/src/mo_expression.F90
index e83b5af489161ebeaa12205fc9d6e134b3d4b6f6..f73f6d6178a92690e5fefd006f8910cd07f75e36 100644
--- a/src/mo_expression.F90
+++ b/src/mo_expression.F90
@@ -33,12 +33,10 @@
 
 MODULE mo_expression
   USE, INTRINSIC :: ISO_C_BINDING
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: real64
+  USE mo_iconlib_kind, ONLY: wp
   IMPLICIT NONE
   PRIVATE
 
-  INTEGER, PARAMETER :: real_kind = real64
-
   INTEGER, PARAMETER :: MAX_BUF_LEN = 1024 ! max length of expression operator list
   INTEGER, PARAMETER :: RESULT_STACKSIZE = 1024 ! max size of result stack
   INTEGER, PARAMETER :: MAX_NAME_LEN = 32 ! max length of variable name
@@ -89,9 +87,9 @@ MODULE mo_expression
 
   ! single field on result stack: scalar, 2D or 3D results
   TYPE t_result_item
-    REAL(real_kind), POINTER :: ptr_0D => NULL()
-    REAL(real_kind), POINTER :: ptr_2D(:, :) => NULL()
-    REAL(real_kind), POINTER :: ptr_3D(:, :, :) => NULL()
+    REAL(wp), POINTER :: ptr_0D => NULL()
+    REAL(wp), POINTER :: ptr_2D(:, :) => NULL()
+    REAL(wp), POINTER :: ptr_3D(:, :, :) => NULL()
   END TYPE t_result_item
 
   ! stack for result computation, base class
@@ -219,7 +217,7 @@ MODULE mo_expression
   END INTERFACE
 
   TYPE, PUBLIC, EXTENDS(t_stack_op) :: t_op_value ! add value to result stack
-    REAL(real_kind) :: val
+    REAL(wp) :: val
   CONTAINS
     PROCEDURE, PUBLIC :: eval_0D => stack_op_value_eval_0D
     PROCEDURE, PUBLIC :: eval_2D => stack_op_value_eval_2D
@@ -365,7 +363,7 @@ CONTAINS
 
   SUBROUTINE result_stack_pop_0D(this, res)
     CLASS(t_result_stack_0D), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER :: res
+    REAL(wp), POINTER :: res
     res => this%rstack(this%ridx)%ptr_0D
     NULLIFY (this%rstack(this%ridx)%ptr_0D)
     this%ridx = this%ridx - 1
@@ -373,7 +371,7 @@ CONTAINS
 
   SUBROUTINE result_stack_pop_2D(this, res)
     CLASS(t_result_stack_2D), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER :: res(:, :)
+    REAL(wp), POINTER :: res(:, :)
     res => this%rstack(this%ridx)%ptr_2D
     NULLIFY (this%rstack(this%ridx)%ptr_2D)
     this%ridx = this%ridx - 1
@@ -381,7 +379,7 @@ CONTAINS
 
   SUBROUTINE result_stack_pop_3D(this, res)
     CLASS(t_result_stack_3D), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER :: res(:, :, :)
+    REAL(wp), POINTER :: res(:, :, :)
     res => this%rstack(this%ridx)%ptr_3D
     NULLIFY (this%rstack(this%ridx)%ptr_3D)
     this%ridx = this%ridx - 1
@@ -389,21 +387,21 @@ CONTAINS
 
   SUBROUTINE result_stack_push_0D(this, val)
     CLASS(t_result_stack_0D), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER :: val
+    REAL(wp), POINTER :: val
     this%ridx = this%ridx + 1
     this%rstack(this%ridx)%ptr_0D => val
   END SUBROUTINE result_stack_push_0D
 
   SUBROUTINE result_stack_push_2D(this, val)
     CLASS(t_result_stack_2D), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER                :: val(:, :)
+    REAL(wp), POINTER                :: val(:, :)
     this%ridx = this%ridx + 1
     this%rstack(this%ridx)%ptr_2D => val
   END SUBROUTINE result_stack_push_2D
 
   SUBROUTINE result_stack_push_3D(this, val)
     CLASS(t_result_stack_3D), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER                :: val(:, :, :)
+    REAL(wp), POINTER                :: val(:, :, :)
     this%ridx = this%ridx + 1
     this%rstack(this%ridx)%ptr_3D => val
   END SUBROUTINE result_stack_push_3D
@@ -442,7 +440,7 @@ CONTAINS
       SELECT CASE (cqueue%list(i)%etype)
       CASE (VALUE)
         ALLOCATE (op_value)
-        op_value%val = REAL(cqueue%list(i)%val, KIND=real_kind)
+        op_value%val = REAL(cqueue%list(i)%val, KIND=wp)
         expr_list%list(i)%p => op_value
       CASE (VARIABLE)
         ALLOCATE (op_variable)
@@ -575,12 +573,12 @@ CONTAINS
   ! evaluate an expression, i.e. evaluate its list of operators.
   SUBROUTINE expr_list_evaluate_0D(this, evaluate)
     CLASS(expression), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER :: evaluate
+    REAL(wp), POINTER :: evaluate
     ! local variables
     INTEGER                    :: i, dim1, dim2, dim3
     TYPE(t_result_stack_0D)    :: result_stack
     TYPE(t_var_ptr)            :: var_shape
-    REAL(real_kind), POINTER   :: ptr
+    REAL(wp), POINTER   :: ptr
 
     IF (this%isize <= 0) THEN
       this%err_no = ERR_EXPRESSION_EMPTY
@@ -616,12 +614,12 @@ CONTAINS
 
   SUBROUTINE expr_list_evaluate_2D(this, evaluate)
     CLASS(expression), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER :: evaluate(:, :)
+    REAL(wp), POINTER :: evaluate(:, :)
     ! local variables
     INTEGER                    :: i, dim1, dim2, dim3, size1, size2
     TYPE(t_result_stack_2D)    :: result_stack
     TYPE(t_var_ptr)            :: var_shape
-    REAL(real_kind), DIMENSION(:, :), POINTER :: ptr
+    REAL(wp), DIMENSION(:, :), POINTER :: ptr
 
     IF (this%isize <= 0) THEN
       this%err_no = ERR_EXPRESSION_EMPTY
@@ -658,7 +656,7 @@ CONTAINS
             evaluate = result_stack%rstack(1)%ptr_2D
           ELSE
             this%err_no = ERR_DIMENSION_MISMATCH
-            evaluate(:, :) = 0._real_kind
+            evaluate(:, :) = 0._wp
           END IF
         END IF
       END IF
@@ -669,12 +667,12 @@ CONTAINS
 
   SUBROUTINE expr_list_evaluate_3D(this, evaluate)
     CLASS(expression), INTENT(INOUT) :: this
-    REAL(real_kind), POINTER :: evaluate(:, :, :)
+    REAL(wp), POINTER :: evaluate(:, :, :)
     ! local variables
     INTEGER                    :: i, dim1, dim2, dim3, size1, size2, size3
     TYPE(t_result_stack_3D)    :: result_stack
     TYPE(t_var_ptr)            :: var_shape
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: ptr
+    REAL(wp), DIMENSION(:, :, :), POINTER :: ptr
 
     IF (this%isize <= 0) THEN
       this%err_no = ERR_EXPRESSION_EMPTY
@@ -711,7 +709,7 @@ CONTAINS
           evaluate = result_stack%rstack(1)%ptr_3D
         ELSE
           this%err_no = ERR_DIMENSION_MISMATCH
-          evaluate(:, :, :) = 0._real_kind
+          evaluate(:, :, :) = 0._wp
         END IF
       END IF
       CALL result_stack%pop(ptr)
@@ -741,7 +739,7 @@ CONTAINS
   SUBROUTINE expr_list_link_0D(this, string, variable)
     CLASS(expression), INTENT(INOUT) :: this
     CHARACTER(len=*), INTENT(IN)    :: string
-    REAL(real_kind), TARGET          :: variable
+    REAL(wp), TARGET          :: variable
     ! local variables
     INTEGER :: i
 
@@ -765,7 +763,7 @@ CONTAINS
   SUBROUTINE expr_list_link_2D(this, string, variable)
     CLASS(expression), INTENT(INOUT) :: this
     CHARACTER(len=*), INTENT(IN)    :: string
-    REAL(real_kind), TARGET          :: variable(:, :)
+    REAL(wp), TARGET          :: variable(:, :)
     ! local variables
     INTEGER :: i
 
@@ -789,7 +787,7 @@ CONTAINS
   SUBROUTINE expr_list_link_3D(this, string, variable)
     CLASS(expression), INTENT(INOUT) :: this
     CHARACTER(len=*), INTENT(IN)    :: string
-    REAL(real_kind), TARGET          :: variable(:, :, :)
+    REAL(wp), TARGET          :: variable(:, :, :)
     ! local variables
     INTEGER :: i
 
@@ -812,7 +810,7 @@ CONTAINS
   SUBROUTINE stack_op_value_eval_0D(this, result_stack)
     CLASS(t_op_value), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: ptr
+    REAL(wp), POINTER :: ptr
     ALLOCATE (ptr)
     ptr = this%val
     CALL result_stack%push(ptr)
@@ -821,7 +819,7 @@ CONTAINS
   SUBROUTINE stack_op_value_eval_2D(this, result_stack)
     CLASS(t_op_value), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: ptr
+    REAL(wp), DIMENSION(:, :), POINTER :: ptr
     ALLOCATE (ptr(result_stack%dim1, result_stack%dim2))
     ptr(:, :) = this%val
     CALL result_stack%push(ptr)
@@ -830,7 +828,7 @@ CONTAINS
   SUBROUTINE stack_op_value_eval_3D(this, result_stack)
     CLASS(t_op_value), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: ptr
+    REAL(wp), DIMENSION(:, :, :), POINTER :: ptr
     ALLOCATE (ptr(result_stack%dim1, result_stack%dim2, result_stack%dim3))
     ptr(:, :, :) = this%val
     CALL result_stack%push(ptr)
@@ -839,7 +837,7 @@ CONTAINS
   SUBROUTINE stack_op_variable_eval_0D(this, result_stack)
     CLASS(t_op_variable), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: ptr
+    REAL(wp), POINTER :: ptr
     IF (ASSOCIATED(this%var%p%ptr_0D)) THEN
       ALLOCATE (ptr)
       ptr = this%var%p%ptr_0D
@@ -852,7 +850,7 @@ CONTAINS
   SUBROUTINE stack_op_variable_eval_2D(this, result_stack)
     CLASS(t_op_variable), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: ptr
+    REAL(wp), DIMENSION(:, :), POINTER :: ptr
     INTEGER :: pdim1, pdim2
     ALLOCATE (ptr(result_stack%dim1, result_stack%dim2))
     IF (ASSOCIATED(this%var%p%ptr_0D)) THEN
@@ -876,7 +874,7 @@ CONTAINS
   SUBROUTINE stack_op_variable_eval_3D(this, result_stack)
     CLASS(t_op_variable), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: ptr
+    REAL(wp), DIMENSION(:, :, :), POINTER :: ptr
     INTEGER :: pdim1, pdim2, pdim3, i
     ALLOCATE (ptr(result_stack%dim1, result_stack%dim2, result_stack%dim3))
     IF (ASSOCIATED(this%var%p%ptr_0D)) THEN
@@ -917,7 +915,7 @@ CONTAINS
   SUBROUTINE stack_op_plus_eval_0D(this, result_stack)
     CLASS(t_op_plus), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1 + arg2
@@ -928,7 +926,7 @@ CONTAINS
   SUBROUTINE stack_op_plus_eval_2D(this, result_stack)
     CLASS(t_op_plus), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1 + arg2
@@ -939,7 +937,7 @@ CONTAINS
   SUBROUTINE stack_op_plus_eval_3D(this, result_stack)
     CLASS(t_op_plus), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1 + arg2
@@ -950,7 +948,7 @@ CONTAINS
   SUBROUTINE stack_op_minus_eval_0D(this, result_stack)
     CLASS(t_op_minus), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1 - arg2
@@ -961,7 +959,7 @@ CONTAINS
   SUBROUTINE stack_op_minus_eval_2D(this, result_stack)
     CLASS(t_op_minus), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1 - arg2
@@ -972,7 +970,7 @@ CONTAINS
   SUBROUTINE stack_op_minus_eval_3D(this, result_stack)
     CLASS(t_op_minus), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1 - arg2
@@ -983,7 +981,7 @@ CONTAINS
   SUBROUTINE stack_op_mul_eval_0D(this, result_stack)
     CLASS(t_op_mul), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1*arg2
@@ -994,7 +992,7 @@ CONTAINS
   SUBROUTINE stack_op_mul_eval_2D(this, result_stack)
     CLASS(t_op_mul), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1*arg2
@@ -1005,7 +1003,7 @@ CONTAINS
   SUBROUTINE stack_op_mul_eval_3D(this, result_stack)
     CLASS(t_op_mul), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1*arg2
@@ -1016,7 +1014,7 @@ CONTAINS
   SUBROUTINE stack_op_div_eval_0D(this, result_stack)
     CLASS(t_op_div), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1/arg2
@@ -1027,7 +1025,7 @@ CONTAINS
   SUBROUTINE stack_op_div_eval_2D(this, result_stack)
     CLASS(t_op_div), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1/arg2
@@ -1038,7 +1036,7 @@ CONTAINS
   SUBROUTINE stack_op_div_eval_3D(this, result_stack)
     CLASS(t_op_div), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = arg1/arg2
@@ -1049,12 +1047,12 @@ CONTAINS
   SUBROUTINE stack_op_pow_eval_0D(this, result_stack)
     CLASS(t_op_pow), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     INTEGER :: exponent
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     exponent = NINT(arg2)
-    IF (arg2 == REAL(exponent, real_kind)) THEN
+    IF (arg2 == REAL(exponent, wp)) THEN
       arg1 = arg1**exponent
     ELSE
       arg1 = arg1**arg2
@@ -1066,12 +1064,12 @@ CONTAINS
   SUBROUTINE stack_op_pow_eval_2D(this, result_stack)
     CLASS(t_op_pow), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     INTEGER, DIMENSION(:, :), ALLOCATABLE :: exponent
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     exponent = NINT(arg2)
-    IF (ALL(arg2 == REAL(exponent, real_kind))) THEN
+    IF (ALL(arg2 == REAL(exponent, wp))) THEN
       arg1 = arg1**exponent
     ELSE
       arg1 = arg1**arg2
@@ -1083,12 +1081,12 @@ CONTAINS
   SUBROUTINE stack_op_pow_eval_3D(this, result_stack)
     CLASS(t_op_pow), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: exponent
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     exponent = NINT(arg2)
-    IF (ALL(arg2 == REAL(exponent, real_kind))) THEN
+    IF (ALL(arg2 == REAL(exponent, wp))) THEN
       arg1 = arg1**exponent
     ELSE
       arg1 = arg1**arg2
@@ -1100,13 +1098,13 @@ CONTAINS
   SUBROUTINE stack_op_gt_eval_0D(this, result_stack)
     CLASS(t_op_gt), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     IF (arg1 > arg2) THEN
-      arg1 = 1._real_kind
+      arg1 = 1._wp
     ELSE
-      arg1 = 0._real_kind
+      arg1 = 0._wp
     END IF
     CALL result_stack%push(arg1)
     DEALLOCATE (arg2)
@@ -1115,13 +1113,13 @@ CONTAINS
   SUBROUTINE stack_op_gt_eval_2D(this, result_stack)
     CLASS(t_op_gt), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     WHERE (arg1(:, :) > arg2(:, :))
-      arg1 = 1._real_kind
+      arg1 = 1._wp
     ELSEWHERE
-      arg1 = 0._real_kind
+      arg1 = 0._wp
     END WHERE
     CALL result_stack%push(arg1)
     DEALLOCATE (arg2)
@@ -1130,13 +1128,13 @@ CONTAINS
   SUBROUTINE stack_op_gt_eval_3D(this, result_stack)
     CLASS(t_op_gt), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     WHERE (arg1(:, :, :) > arg2(:, :, :))
-      arg1 = 1._real_kind
+      arg1 = 1._wp
     ELSEWHERE
-      arg1 = 0._real_kind
+      arg1 = 0._wp
     END WHERE
     CALL result_stack%push(arg1)
     DEALLOCATE (arg2)
@@ -1145,13 +1143,13 @@ CONTAINS
   SUBROUTINE stack_op_lt_eval_0D(this, result_stack)
     CLASS(t_op_lt), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     IF (arg1 < arg2) THEN
-      arg1 = 1._real_kind
+      arg1 = 1._wp
     ELSE
-      arg1 = 0._real_kind
+      arg1 = 0._wp
     END IF
     CALL result_stack%push(arg1)
     DEALLOCATE (arg2)
@@ -1160,13 +1158,13 @@ CONTAINS
   SUBROUTINE stack_op_lt_eval_2D(this, result_stack)
     CLASS(t_op_lt), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     WHERE (arg1(:, :) < arg2)
-      arg1(:, :) = 1._real_kind
+      arg1(:, :) = 1._wp
     ELSEWHERE
-      arg1(:, :) = 0._real_kind
+      arg1(:, :) = 0._wp
     END WHERE
     CALL result_stack%push(arg1)
     DEALLOCATE (arg2)
@@ -1175,13 +1173,13 @@ CONTAINS
   SUBROUTINE stack_op_lt_eval_3D(this, result_stack)
     CLASS(t_op_lt), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     WHERE (arg1(:, :, :) < arg2(:, :, :))
-      arg1 = 1._real_kind
+      arg1 = 1._wp
     ELSEWHERE
-      arg1 = 0._real_kind
+      arg1 = 0._wp
     END WHERE
     CALL result_stack%push(arg1)
     DEALLOCATE (arg2)
@@ -1190,7 +1188,7 @@ CONTAINS
   SUBROUTINE stack_op_exp_eval_0D(this, result_stack)
     CLASS(t_op_exp), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1
+    REAL(wp), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = EXP(arg1)
     CALL result_stack%push(arg1)
@@ -1199,7 +1197,7 @@ CONTAINS
   SUBROUTINE stack_op_exp_eval_2D(this, result_stack)
     CLASS(t_op_exp), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = EXP(arg1)
     CALL result_stack%push(arg1)
@@ -1208,7 +1206,7 @@ CONTAINS
   SUBROUTINE stack_op_exp_eval_3D(this, result_stack)
     CLASS(t_op_exp), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = EXP(arg1)
     CALL result_stack%push(arg1)
@@ -1217,7 +1215,7 @@ CONTAINS
   SUBROUTINE stack_op_log_eval_0D(this, result_stack)
     CLASS(t_op_log), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1
+    REAL(wp), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = LOG(arg1)
     CALL result_stack%push(arg1)
@@ -1226,7 +1224,7 @@ CONTAINS
   SUBROUTINE stack_op_log_eval_2D(this, result_stack)
     CLASS(t_op_log), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = LOG(arg1)
     CALL result_stack%push(arg1)
@@ -1235,7 +1233,7 @@ CONTAINS
   SUBROUTINE stack_op_log_eval_3D(this, result_stack)
     CLASS(t_op_log), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = LOG(arg1)
     CALL result_stack%push(arg1)
@@ -1244,7 +1242,7 @@ CONTAINS
   SUBROUTINE stack_op_sin_eval_0D(this, result_stack)
     CLASS(t_op_sin), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1
+    REAL(wp), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = SIN(arg1)
     CALL result_stack%push(arg1)
@@ -1253,7 +1251,7 @@ CONTAINS
   SUBROUTINE stack_op_sin_eval_2D(this, result_stack)
     CLASS(t_op_sin), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = SIN(arg1)
     CALL result_stack%push(arg1)
@@ -1262,7 +1260,7 @@ CONTAINS
   SUBROUTINE stack_op_sin_eval_3D(this, result_stack)
     CLASS(t_op_sin), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = SIN(arg1)
     CALL result_stack%push(arg1)
@@ -1271,7 +1269,7 @@ CONTAINS
   SUBROUTINE stack_op_cos_eval_0D(this, result_stack)
     CLASS(t_op_cos), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1
+    REAL(wp), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = COS(arg1)
     CALL result_stack%push(arg1)
@@ -1280,7 +1278,7 @@ CONTAINS
   SUBROUTINE stack_op_cos_eval_2D(this, result_stack)
     CLASS(t_op_cos), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = COS(arg1)
     CALL result_stack%push(arg1)
@@ -1289,7 +1287,7 @@ CONTAINS
   SUBROUTINE stack_op_cos_eval_3D(this, result_stack)
     CLASS(t_op_cos), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = COS(arg1)
     CALL result_stack%push(arg1)
@@ -1298,7 +1296,7 @@ CONTAINS
   SUBROUTINE stack_op_min_eval_0D(this, result_stack)
     CLASS(t_op_min), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = MIN(arg1, arg2)
@@ -1308,7 +1306,7 @@ CONTAINS
   SUBROUTINE stack_op_min_eval_2D(this, result_stack)
     CLASS(t_op_min), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = MIN(arg1, arg2)
@@ -1318,7 +1316,7 @@ CONTAINS
   SUBROUTINE stack_op_min_eval_3D(this, result_stack)
     CLASS(t_op_min), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = MIN(arg1, arg2)
@@ -1328,7 +1326,7 @@ CONTAINS
   SUBROUTINE stack_op_max_eval_0D(this, result_stack)
     CLASS(t_op_max), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2
+    REAL(wp), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = MAX(arg1, arg2)
@@ -1339,7 +1337,7 @@ CONTAINS
   SUBROUTINE stack_op_max_eval_2D(this, result_stack)
     CLASS(t_op_max), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = MAX(arg1, arg2)
@@ -1350,7 +1348,7 @@ CONTAINS
   SUBROUTINE stack_op_max_eval_3D(this, result_stack)
     CLASS(t_op_max), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
     arg1 = MAX(arg1, arg2)
@@ -1361,7 +1359,7 @@ CONTAINS
   SUBROUTINE stack_op_if_eval_0D(this, result_stack)
     CLASS(t_op_if), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1, arg2, arg3
+    REAL(wp), POINTER :: arg1, arg2, arg3
     CALL result_stack%pop(arg3)
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
@@ -1378,11 +1376,11 @@ CONTAINS
   SUBROUTINE stack_op_if_eval_2D(this, result_stack)
     CLASS(t_op_if), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1, arg2, arg3
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1, arg2, arg3
     CALL result_stack%pop(arg3)
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
-    WHERE (arg1(:, :) <= 0._real_kind)
+    WHERE (arg1(:, :) <= 0._wp)
       arg2(:, :) = arg3(:, :)
     END WHERE
     CALL result_stack%push(arg2)
@@ -1392,11 +1390,11 @@ CONTAINS
   SUBROUTINE stack_op_if_eval_3D(this, result_stack)
     CLASS(t_op_if), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1, arg2, arg3
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1, arg2, arg3
     CALL result_stack%pop(arg3)
     CALL result_stack%pop(arg2)
     CALL result_stack%pop(arg1)
-    WHERE (arg1(:, :, :) <= 0._real_kind)
+    WHERE (arg1(:, :, :) <= 0._wp)
       arg2(:, :, :) = arg3(:, :, :)
     END WHERE
     CALL result_stack%push(arg2)
@@ -1406,7 +1404,7 @@ CONTAINS
   SUBROUTINE stack_op_sqrt_eval_0D(this, result_stack)
     CLASS(t_op_sqrt), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: p
+    REAL(wp), POINTER :: p
     CALL result_stack%pop(p)
     p = SQRT(p)
     CALL result_stack%push(p)
@@ -1415,7 +1413,7 @@ CONTAINS
   SUBROUTINE stack_op_sqrt_eval_2D(this, result_stack)
     CLASS(t_op_sqrt), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: p(:, :)
+    REAL(wp), POINTER :: p(:, :)
     CALL result_stack%pop(p)
     p = SQRT(p)
     CALL result_stack%push(p)
@@ -1424,7 +1422,7 @@ CONTAINS
   SUBROUTINE stack_op_sqrt_eval_3D(this, result_stack)
     CLASS(t_op_sqrt), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: p(:, :, :)
+    REAL(wp), POINTER :: p(:, :, :)
     CALL result_stack%pop(p)
     p = SQRT(p)
     CALL result_stack%push(p)
@@ -1433,7 +1431,7 @@ CONTAINS
   SUBROUTINE stack_op_erf_eval_0D(this, result_stack)
     CLASS(t_op_erf), INTENT(INOUT) :: this
     CLASS(t_result_stack_0D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), POINTER :: arg1
+    REAL(wp), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = ERF(arg1)
     CALL result_stack%push(arg1)
@@ -1442,7 +1440,7 @@ CONTAINS
   SUBROUTINE stack_op_erf_eval_2D(this, result_stack)
     CLASS(t_op_erf), INTENT(INOUT) :: this
     CLASS(t_result_stack_2D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = ERF(arg1)
     CALL result_stack%push(arg1)
@@ -1451,7 +1449,7 @@ CONTAINS
   SUBROUTINE stack_op_erf_eval_3D(this, result_stack)
     CLASS(t_op_erf), INTENT(INOUT) :: this
     CLASS(t_result_stack_3D), INTENT(INOUT) :: result_stack
-    REAL(real_kind), DIMENSION(:, :, :), POINTER :: arg1
+    REAL(wp), DIMENSION(:, :, :), POINTER :: arg1
     CALL result_stack%pop(arg1)
     arg1 = ERF(arg1)
     CALL result_stack%push(arg1)
diff --git a/src/mo_fortran_tools.F90 b/src/mo_fortran_tools.F90
index 708296f6be9f2258417571dd56489527e9b162bf..4adddc68ed48ab81501e4947129fa022ebfdde95 100644
--- a/src/mo_fortran_tools.F90
+++ b/src/mo_fortran_tools.F90
@@ -16,10 +16,7 @@
 ! subroutines.
 MODULE mo_fortran_tools
 
-  USE ISO_FORTRAN_ENV, ONLY: wp => real64, &
-                             sp => real32, &
-                             dp => real64, &
-                             ik4 => int32
+  USE mo_iconlib_kind, ONLY: wp, vp, dp, sp, ik4 => i4 ! i4 used for indexing in many subroutines
   USE mo_exception, ONLY: finish
 #ifdef _OPENACC
   USE openacc
@@ -33,18 +30,18 @@ MODULE mo_fortran_tools
   PUBLIC :: t_ptr_2d3d, t_ptr_2d3d_vp
   PUBLIC :: assign_if_present_allocatable
   PUBLIC :: if_associated
-  PUBLIC :: t_ptr_1d, t_ptr_1d_sp, t_ptr_1d_int
-  PUBLIC :: t_ptr_1d_ptr_1d
-  PUBLIC :: t_ptr_2d, t_ptr_2d_sp, t_ptr_2d_int
-  PUBLIC :: t_ptr_3d, t_ptr_3d_sp, t_ptr_3d_int
+  PUBLIC :: t_ptr_1d_sp, t_ptr_1d_dp, t_ptr_1d_int
+  PUBLIC :: t_ptr_2d_sp, t_ptr_2d_dp, t_ptr_2d_int
+  PUBLIC :: t_ptr_3d_sp, t_ptr_3d_dp, t_ptr_3d_int
+  PUBLIC :: t_ptr_4d_sp, t_ptr_4d_dp, t_ptr_4d_int
   PUBLIC :: t_ptr_i2d3d
-  PUBLIC :: t_ptr_4d, t_ptr_4d_sp, t_ptr_4d_int
   PUBLIC :: t_ptr_tracer
   PUBLIC :: copy, init, swap, negative2zero
   PUBLIC :: var_scale, var_add
   PUBLIC :: init_zero_contiguous_dp, init_zero_contiguous_sp
   PUBLIC :: init_contiguous_dp, init_contiguous_sp
   PUBLIC :: init_contiguous_i4, init_contiguous_l
+  PUBLIC :: init_contiguous, init_zero_contiguous
   PUBLIC :: minval_1d
   PUBLIC :: minval_2d
   PUBLIC :: resize_arr_c1d
@@ -57,82 +54,72 @@ MODULE mo_fortran_tools
 
   PRIVATE
 
-#ifdef __MIXED_PRECISION
-  INTEGER, PARAMETER :: vp = sp
-#else
-  INTEGER, PARAMETER :: vp = wp
-#endif
-
-  TYPE t_ptr_1d
-    REAL(wp), POINTER :: p(:) ! pointer to 1D (spatial) array
-  END TYPE t_ptr_1d
+  TYPE t_ptr_1d_dp
+    REAL(dp), POINTER :: p(:) => NULL() ! pointer to 1D (spatial) array
+  END TYPE t_ptr_1d_dp
 
   TYPE t_ptr_1d_sp
-    REAL(sp), POINTER :: p(:) ! pointer to 1D (spatial) array
+    REAL(sp), POINTER :: p(:) => NULL() ! pointer to 1D (spatial) array
   END TYPE t_ptr_1d_sp
 
   TYPE t_ptr_1d_int
-    INTEGER, POINTER :: p(:) ! pointer to 1D (spatial) array
+    INTEGER, POINTER :: p(:) => NULL() ! pointer to 1D (spatial) array
   END TYPE t_ptr_1d_int
 
-  TYPE t_ptr_1d_ptr_1d
-    TYPE(t_ptr_1d), POINTER :: p(:) ! pointer to a 1D array of pointers to 1D (spatial) arrays
-  END TYPE t_ptr_1d_ptr_1d
-
-  TYPE t_ptr_2d
-    REAL(dp), POINTER :: p(:, :) ! pointer to 2D (spatial) array
-  END TYPE t_ptr_2d
+  TYPE t_ptr_2d_dp
+    REAL(dp), POINTER :: p(:, :) => NULL() ! pointer to 2D (spatial) array
+  END TYPE t_ptr_2d_dp
 
   TYPE t_ptr_2d_sp
-    REAL(sp), POINTER :: p(:, :) ! pointer to 2D (spatial) array
+    REAL(sp), POINTER :: p(:, :) => NULL() ! pointer to 2D (spatial) array
   END TYPE t_ptr_2d_sp
 
   TYPE t_ptr_2d_int
-    INTEGER, POINTER :: p(:, :) ! pointer to 2D (spatial) array
+    INTEGER, POINTER :: p(:, :) => NULL() ! pointer to 2D (spatial) array
   END TYPE t_ptr_2d_int
 
-  TYPE t_ptr_3d
-    REAL(dp), POINTER :: p(:, :, :) ! pointer to 3D (spatial) array
-  END TYPE t_ptr_3d
+  TYPE t_ptr_3d_dp
+    REAL(dp), POINTER :: p(:, :, :) => NULL() ! pointer to 3D (spatial) array
+  END TYPE t_ptr_3d_dp
 
   TYPE t_ptr_3d_sp
-    REAL(sp), POINTER :: p(:, :, :) ! pointer to 3D (spatial) array
+    REAL(sp), POINTER :: p(:, :, :) => NULL() ! pointer to 3D (spatial) array
   END TYPE t_ptr_3d_sp
 
   TYPE t_ptr_3d_int
-    INTEGER, POINTER :: p(:, :, :) ! pointer to 3D (spatial) array
+    INTEGER, POINTER :: p(:, :, :) => NULL() ! pointer to 3D (spatial) array
   END TYPE t_ptr_3d_int
 
-  TYPE t_ptr_4d
-    REAL(dp), POINTER :: p(:, :, :, :) ! pointer to 3D (spatial) array
-  END TYPE t_ptr_4d
-
   TYPE t_ptr_4d_sp
-    REAL(sp), POINTER :: p(:, :, :, :) ! pointer to 3D (spatial) array
+    REAL(sp), POINTER :: p(:, :, :, :) => NULL() ! pointer to 4D (spatial) array
   END TYPE t_ptr_4d_sp
 
+  TYPE t_ptr_4d_dp
+    REAL(dp), POINTER :: p(:, :, :, :) => NULL() ! pointer to 3D (spatial) array
+  END TYPE t_ptr_4d_dp
+
   TYPE t_ptr_4d_int
-    INTEGER, POINTER :: p(:, :, :, :) ! pointer to 3D (spatial) array
+    INTEGER, POINTER :: p(:, :, :, :) => NULL() ! pointer to 3D (spatial) array
   END TYPE t_ptr_4d_int
 
   TYPE t_ptr_2d3d
-    REAL(wp), POINTER :: p_3d(:, :, :) ! REAL pointer to 3D (spatial) array
-    REAL(wp), POINTER :: p_2d(:, :) ! REAL pointer to 2D (spatial) array
+    REAL(wp), POINTER :: p_3d(:, :, :) => NULL() ! REAL pointer to 3D (spatial) array
+    REAL(wp), POINTER :: p_2d(:, :) => NULL() ! REAL pointer to 2D (spatial) array
   END TYPE t_ptr_2d3d
 
   TYPE t_ptr_2d3d_vp
-    REAL(vp), POINTER :: p_3d(:, :, :) ! REAL pointer to 3D (spatial) array
-    REAL(vp), POINTER :: p_2d(:, :) ! REAL pointer to 2D (spatial) array
+    REAL(vp), POINTER :: p_3d(:, :, :) => NULL() ! REAL pointer to 3D (spatial) array
+    REAL(vp), POINTER :: p_2d(:, :) => NULL() ! REAL pointer to 2D (spatial) array
   END TYPE t_ptr_2d3d_vp
 
   TYPE t_ptr_i2d3d
-    INTEGER, POINTER :: p_3d(:, :, :) ! INTEGER pointer to 3D (spatial) array
-    INTEGER, POINTER :: p_2d(:, :) ! INTEGER pointer to 2D (spatial) array
+    INTEGER, POINTER :: p_3d(:, :, :) => NULL() ! INTEGER pointer to 3D (spatial) array
+    INTEGER, POINTER :: p_2d(:, :) => NULL() ! INTEGER pointer to 2D (spatial) array
   END TYPE t_ptr_i2d3d
 
   ! Type to pass pointer arrays to convection and turbulent diffusion subroutines
   TYPE t_ptr_tracer
-    REAL(wp), POINTER :: ptr(:, :)
+    REAL(wp), POINTER :: ptr(:, :) => NULL()
     INTEGER           :: idx_tracer
   END TYPE t_ptr_tracer
 
@@ -173,7 +160,12 @@ MODULE mo_fortran_tools
     MODULE PROCEDURE copy_3d_dp
     MODULE PROCEDURE copy_4d_dp
     MODULE PROCEDURE copy_5d_dp
+    MODULE PROCEDURE copy_1d_sp
+    MODULE PROCEDURE copy_2d_sp
+    MODULE PROCEDURE copy_3d_sp
+    MODULE PROCEDURE copy_4d_sp
     MODULE PROCEDURE copy_5d_sp
+    MODULE PROCEDURE copy_3d_dpsp
     MODULE PROCEDURE copy_2d_spdp
     MODULE PROCEDURE copy_3d_spdp
     MODULE PROCEDURE copy_4d_spdp
@@ -189,6 +181,7 @@ MODULE mo_fortran_tools
     MODULE PROCEDURE init_zero_1d_dp
     MODULE PROCEDURE init_zero_1d_sp
     MODULE PROCEDURE init_zero_2d_dp
+    MODULE PROCEDURE init_zero_2d_sp
     MODULE PROCEDURE init_zero_2d_i4
     MODULE PROCEDURE init_zero_3d_dp
     MODULE PROCEDURE init_zero_3d_sp
@@ -197,8 +190,11 @@ MODULE mo_fortran_tools
     MODULE PROCEDURE init_zero_4d_sp
     MODULE PROCEDURE init_zero_4d_i4
     MODULE PROCEDURE init_1d_dp
+    MODULE PROCEDURE init_1d_sp
     MODULE PROCEDURE init_2d_dp
+    MODULE PROCEDURE init_2d_sp
     MODULE PROCEDURE init_3d_dp
+    MODULE PROCEDURE init_3d_sp
     MODULE PROCEDURE init_3d_spdp
     MODULE PROCEDURE init_5d_dp
     MODULE PROCEDURE init_5d_sp
@@ -208,10 +204,12 @@ MODULE mo_fortran_tools
 
   INTERFACE negative2zero
     MODULE PROCEDURE negative2zero_4d_dp
+    MODULE PROCEDURE negative2zero_4d_sp
   END INTERFACE negative2zero
 
   INTERFACE var_scale
     MODULE PROCEDURE var_scale_3d_dp
+    MODULE PROCEDURE var_scale_3d_sp
   END INTERFACE var_scale
 
   INTERFACE var_add
@@ -228,6 +226,10 @@ MODULE mo_fortran_tools
     MODULE PROCEDURE DO_DEALLOCATE_r3D
     MODULE PROCEDURE DO_DEALLOCATE_r2D
     MODULE PROCEDURE DO_DEALLOCATE_r1D
+    MODULE PROCEDURE DO_DEALLOCATE_s4D
+    MODULE PROCEDURE DO_DEALLOCATE_s3D
+    MODULE PROCEDURE DO_DEALLOCATE_s2D
+    MODULE PROCEDURE DO_DEALLOCATE_s1D
     MODULE PROCEDURE DO_DEALLOCATE_i3D
     MODULE PROCEDURE DO_DEALLOCATE_i2D
     MODULE PROCEDURE DO_DEALLOCATE_i1D
@@ -253,6 +255,18 @@ MODULE mo_fortran_tools
     MODULE PROCEDURE insert_dimension_i4_6_5, insert_dimension_i4_6_5_s
   END INTERFACE insert_dimension
 
+  INTERFACE init_contiguous
+    MODULE PROCEDURE init_contiguous_sp
+    MODULE PROCEDURE init_contiguous_dp
+    MODULE PROCEDURE init_contiguous_i4
+    MODULE PROCEDURE init_contiguous_l
+  END INTERFACE init_contiguous
+
+  INTERFACE init_zero_contiguous
+    MODULE PROCEDURE init_zero_contiguous_dp
+    MODULE PROCEDURE init_zero_contiguous_sp
+  END INTERFACE init_zero_contiguous
+
   INTEGER, PARAMETER :: SUCCESS = 0
 
 CONTAINS
@@ -321,8 +335,8 @@ CONTAINS
   END SUBROUTINE assign_if_present_integers
 
   SUBROUTINE assign_if_present_real(y, x)
-    REAL(wp), INTENT(INOUT)        :: y
-    REAL(wp), INTENT(IN), OPTIONAL :: x
+    REAL(dp), INTENT(INOUT)        :: y
+    REAL(dp), INTENT(IN), OPTIONAL :: x
     IF (.NOT. PRESENT(x)) RETURN
     IF (x == -HUGE(x)) RETURN
     y = x
@@ -672,6 +686,129 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE copy_5d_dp
 
+  !> copy state, omp parallel, does not wait for other threads to complete
+  SUBROUTINE copy_1d_sp(src, dest, lacc, opt_acc_async)
+    REAL(sp), INTENT(IN) :: src(:)
+    REAL(sp), INTENT(OUT) :: dest(:)
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+    INTEGER :: i1, m1
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(dest, 1)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) IF(lzacc)
+    !$omp do private(i1)
+    DO i1 = 1, m1
+      dest(i1) = src(i1)
+    END DO
+    !$omp end do nowait
+    !$ACC END PARALLEL LOOP
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE copy_1d_sp
+
+  !> copy state, omp parallel, does not wait for other threads to complete
+  SUBROUTINE copy_2d_sp(src, dest, lacc, opt_acc_async)
+    REAL(sp), INTENT(IN) :: src(:, :)
+    REAL(sp), INTENT(OUT) :: dest(:, :)
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+    INTEGER :: i1, i2, m1, m2
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(dest, 1)
+    m2 = SIZE(dest, 2)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(2) IF(lzacc)
+#ifdef __INTEL_COMPILER
+!$omp do private(i1,i2)
+#else
+!$omp do collapse(2)
+#endif
+    DO i2 = 1, m2
+      DO i1 = 1, m1
+        dest(i1, i2) = src(i1, i2)
+      END DO
+    END DO
+!$omp end do nowait
+    CALL acc_wait_if_requested(1, opt_acc_async)
+
+  END SUBROUTINE copy_2d_sp
+
+  !> copy state, omp parallel, does not wait for other threads to complete
+  SUBROUTINE copy_3d_sp(src, dest, lacc, opt_acc_async)
+    REAL(sp), INTENT(IN) :: src(:, :, :)
+    REAL(sp), INTENT(OUT) :: dest(:, :, :)
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+    INTEGER :: i1, i2, i3, m1, m2, m3
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(dest, 1)
+    m2 = SIZE(dest, 2)
+    m3 = SIZE(dest, 3)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(3) IF(lzacc)
+#if (defined(_CRAYFTN) || defined(__INTEL_COMPILER))
+!$omp do private(i1,i2,i3)
+#else
+!$omp do collapse(3)
+#endif
+    DO i3 = 1, m3
+      DO i2 = 1, m2
+        DO i1 = 1, m1
+          dest(i1, i2, i3) = src(i1, i2, i3)
+        END DO
+      END DO
+    END DO
+!$omp end do nowait
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE copy_3d_sp
+
+  !> copy state, omp parallel, does not wait for other threads to complete
+  SUBROUTINE copy_4d_sp(src, dest, lacc, opt_acc_async)
+    REAL(sp), INTENT(IN) :: src(:, :, :, :)
+    REAL(sp), INTENT(OUT) :: dest(:, :, :, :)
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+    INTEGER :: i1, i2, i3, i4, m1, m2, m3, m4
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(dest, 1)
+    m2 = SIZE(dest, 2)
+    m3 = SIZE(dest, 3)
+    m4 = SIZE(dest, 4)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(4) IF(lzacc)
+#if (defined(_CRAYFTN) || defined(__INTEL_COMPILER))
+!$omp do private(i1,i2,i3,i4)
+#else
+!$omp do collapse(4)
+#endif
+    DO i4 = 1, m4
+      DO i3 = 1, m3
+        DO i2 = 1, m2
+          DO i1 = 1, m1
+            dest(i1, i2, i3, i4) = src(i1, i2, i3, i4)
+          END DO
+        END DO
+      END DO
+    END DO
+!$omp end do nowait
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE copy_4d_sp
+
   !> copy state, omp parallel, does not wait for other threads to complete
   SUBROUTINE copy_5d_sp(src, dest, lacc, opt_acc_async)
     REAL(sp), INTENT(IN) :: src(:, :, :, :, :)
@@ -711,6 +848,39 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE copy_5d_sp
 
+  !> copy state, omp parallel, does not wait for other threads to complete
+  SUBROUTINE copy_3d_dpsp(src, dest, lacc, opt_acc_async)
+    REAL(dp), INTENT(IN) :: src(:, :, :)
+    REAL(sp), INTENT(OUT) :: dest(:, :, :)
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+    INTEGER :: i1, i2, i3, m1, m2, m3
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(dest, 1)
+    m2 = SIZE(dest, 2)
+    m3 = SIZE(dest, 3)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(3) IF(lzacc)
+#if (defined(__INTEL_COMPILER))
+!$omp do private(i1,i2,i3)
+#else
+!$omp do collapse(2)
+#endif
+    DO i3 = 1, m3
+      DO i2 = 1, m2
+        DO i1 = 1, m1
+          dest(i1, i2, i3) = REAL(src(i1, i2, i3), KIND=sp)
+        END DO
+      END DO
+    END DO
+!$omp end do nowait
+    CALL acc_wait_if_requested(1, opt_acc_async)
+
+  END SUBROUTINE copy_3d_dpsp
+
   !> copy state, omp parallel, does not wait for other threads to complete
   SUBROUTINE copy_2d_spdp(src, dest, lacc, opt_acc_async)
     REAL(sp), INTENT(IN) :: src(:, :)
@@ -1058,6 +1228,34 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE init_zero_2d_dp
 
+  SUBROUTINE init_zero_2d_sp(init_var, lacc, opt_acc_async)
+    REAL(sp), INTENT(OUT) :: init_var(:, :)
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+    INTEGER :: i1, i2, m1, m2
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(init_var, 1)
+    m2 = SIZE(init_var, 2)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(2) IF(lzacc)
+#if (defined(__INTEL_COMPILER))
+!$omp do private(i1,i2)
+#else
+!$omp do collapse(2)
+#endif
+    DO i2 = 1, m2
+      DO i1 = 1, m1
+        init_var(i1, i2) = 0.0_sp
+      END DO
+    END DO
+!$omp end do nowait
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE init_zero_2d_sp
+
   SUBROUTINE init_zero_2d_i4(init_var, lacc, opt_acc_async)
     INTEGER(ik4), INTENT(OUT) :: init_var(:, :)
     LOGICAL, INTENT(IN) :: lacc
@@ -1307,6 +1505,28 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE init_1d_dp
 
+  SUBROUTINE init_1d_sp(init_var, init_val, lacc, opt_acc_async)
+    REAL(sp), INTENT(OUT) :: init_var(:)
+    REAL(sp), INTENT(IN) :: init_val
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+    INTEGER :: i1, m1
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(init_var, 1)
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) IF(lzacc)
+    !$omp do private(i1)
+    DO i1 = 1, m1
+      init_var(i1) = init_val
+    END DO
+    !$omp end do nowait
+    !$ACC END PARALLEL LOOP
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE init_1d_sp
+
   SUBROUTINE init_2d_dp(init_var, init_val, lacc, opt_acc_async)
     REAL(dp), INTENT(OUT) :: init_var(:, :)
     REAL(dp), INTENT(IN) :: init_val
@@ -1337,6 +1557,36 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE init_2d_dp
 
+  SUBROUTINE init_2d_sp(init_var, init_val, lacc, opt_acc_async)
+    REAL(sp), INTENT(OUT) :: init_var(:, :)
+    REAL(sp), INTENT(IN) :: init_val
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+
+    INTEGER :: i1, i2, m1, m2
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(init_var, 1)
+    m2 = SIZE(init_var, 2)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(2) IF(lzacc)
+#if (defined(__INTEL_COMPILER))
+!$omp do private(i1,i2)
+#else
+!$omp do collapse(2)
+#endif
+    DO i2 = 1, m2
+      DO i1 = 1, m1
+        init_var(i1, i2) = init_val
+      END DO
+    END DO
+!$omp end do nowait
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE init_2d_sp
+
   SUBROUTINE init_3d_dp(init_var, init_val, lacc, opt_acc_async)
     REAL(dp), INTENT(OUT) :: init_var(:, :, :)
     REAL(dp), INTENT(IN) :: init_val
@@ -1370,6 +1620,39 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE init_3d_dp
 
+  SUBROUTINE init_3d_sp(init_var, init_val, lacc, opt_acc_async)
+    REAL(sp), INTENT(OUT) :: init_var(:, :, :)
+    REAL(sp), INTENT(IN) :: init_val
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+
+    INTEGER :: i1, i2, i3, m1, m2, m3
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(init_var, 1)
+    m2 = SIZE(init_var, 2)
+    m3 = SIZE(init_var, 3)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(3) IF(lzacc)
+#if (defined(__INTEL_COMPILER))
+!$omp do private(i1,i2,i3)
+#else
+!$omp do collapse(3)
+#endif
+    DO i3 = 1, m3
+      DO i2 = 1, m2
+        DO i1 = 1, m1
+          init_var(i1, i2, i3) = init_val
+        END DO
+      END DO
+    END DO
+!$omp end do nowait
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE init_3d_sp
+
   SUBROUTINE init_3d_spdp(init_var, init_val, lacc, opt_acc_async)
     REAL(sp), INTENT(OUT) :: init_var(:, :, :)
     REAL(dp), INTENT(IN) :: init_val
@@ -1592,6 +1875,39 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE var_scale_3d_dp
 
+  SUBROUTINE var_scale_3d_sp(var, scale_val, lacc, opt_acc_async)
+    REAL(sp), INTENT(inout) :: var(:, :, :)
+    REAL(sp), INTENT(in) :: scale_val
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+
+    INTEGER :: i1, i2, i3, m1, m2, m3
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(var, 1)
+    m2 = SIZE(var, 2)
+    m3 = SIZE(var, 3)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) ASYNC(1) COLLAPSE(3) IF(lzacc)
+#if (defined(__INTEL_COMPILER))
+!$omp do private(i1,i2,i3)
+#else
+!$omp do collapse(3)
+#endif
+    DO i3 = 1, m3
+      DO i2 = 1, m2
+        DO i1 = 1, m1
+          var(i1, i2, i3) = var(i1, i2, i3)*scale_val
+        END DO
+      END DO
+    END DO
+!$omp end do nowait
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE var_scale_3d_sp
+
   ! add a constant value to a 3D field
   SUBROUTINE var_addc_3d_dp(var, add_val, lacc, opt_acc_async)
     REAL(dp), INTENT(inout) :: var(:, :, :)
@@ -1663,6 +1979,43 @@ CONTAINS
     CALL acc_wait_if_requested(1, opt_acc_async)
   END SUBROUTINE negative2zero_4d_dp
 
+  SUBROUTINE negative2zero_4d_sp(var, lacc, opt_acc_async)
+    REAL(sp), INTENT(inout) :: var(:, :, :, :)
+    LOGICAL, INTENT(IN) :: lacc
+    LOGICAL, INTENT(IN), OPTIONAL :: opt_acc_async
+
+    INTEGER :: i1, i2, i3, i4, m1, m2, m3, m4
+    REAL(sp) :: v
+    LOGICAL :: lzacc
+
+    CALL set_acc_host_or_device(lzacc, lacc)
+
+    m1 = SIZE(var, 1)
+    m2 = SIZE(var, 2)
+    m3 = SIZE(var, 3)
+    m4 = SIZE(var, 4)
+
+    !$ACC PARALLEL LOOP DEFAULT(PRESENT) PRIVATE(v) ASYNC(1) COLLAPSE(4) IF(lzacc)
+#if (defined(__INTEL_COMPILER))
+!$omp do private(i1,i2,i3,i4)
+#else
+!$omp do collapse(4)
+#endif
+    DO i4 = 1, m4
+      DO i3 = 1, m3
+        DO i2 = 1, m2
+          DO i1 = 1, m1
+            v = var(i1, i2, i3, i4)
+            var(i1, i2, i3, i4) = (ABS(v) + v)*0.5_sp
+          END DO
+        END DO
+      END DO
+    END DO
+!$omp end do nowait
+
+    CALL acc_wait_if_requested(1, opt_acc_async)
+  END SUBROUTINE negative2zero_4d_sp
+
   SUBROUTINE init_contiguous_dp(var, n, v, lacc, opt_acc_async)
     INTEGER, INTENT(IN) :: n
     REAL(dp), INTENT(OUT) :: var(n)
@@ -2208,7 +2561,7 @@ CONTAINS
 
   ! AUXILIARY ROUTINES FOR DEALLOCATION
   SUBROUTINE DO_DEALLOCATE_r4D(object)
-    REAL(wp), ALLOCATABLE, INTENT(INOUT) :: object(:, :, :, :)
+    REAL(dp), ALLOCATABLE, INTENT(INOUT) :: object(:, :, :, :)
     INTEGER :: ierrstat
     IF (ALLOCATED(object)) THEN
       DEALLOCATE (object, STAT=ierrstat)
@@ -2217,7 +2570,7 @@ CONTAINS
   END SUBROUTINE DO_DEALLOCATE_R4D
 
   SUBROUTINE DO_DEALLOCATE_r3D(object)
-    REAL(wp), ALLOCATABLE, INTENT(INOUT) :: object(:, :, :)
+    REAL(dp), ALLOCATABLE, INTENT(INOUT) :: object(:, :, :)
     INTEGER :: ierrstat
     IF (ALLOCATED(object)) THEN
       DEALLOCATE (object, STAT=ierrstat)
@@ -2226,7 +2579,7 @@ CONTAINS
   END SUBROUTINE DO_DEALLOCATE_R3D
 
   SUBROUTINE DO_DEALLOCATE_r2D(object)
-    REAL(wp), ALLOCATABLE, INTENT(INOUT) :: object(:, :)
+    REAL(dp), ALLOCATABLE, INTENT(INOUT) :: object(:, :)
     INTEGER :: ierrstat
     IF (ALLOCATED(object)) THEN
       DEALLOCATE (object, STAT=ierrstat)
@@ -2235,7 +2588,7 @@ CONTAINS
   END SUBROUTINE DO_DEALLOCATE_R2D
 
   SUBROUTINE DO_DEALLOCATE_r1D(object)
-    REAL(wp), ALLOCATABLE, INTENT(INOUT) :: object(:)
+    REAL(dp), ALLOCATABLE, INTENT(INOUT) :: object(:)
     INTEGER :: ierrstat
     IF (ALLOCATED(object)) THEN
       DEALLOCATE (object, STAT=ierrstat)
@@ -2243,6 +2596,42 @@ CONTAINS
     END IF
   END SUBROUTINE DO_DEALLOCATE_R1D
 
+  SUBROUTINE DO_DEALLOCATE_s4D(object)
+    REAL(sp), ALLOCATABLE, INTENT(INOUT) :: object(:, :, :, :)
+    INTEGER :: ierrstat
+    IF (ALLOCATED(object)) THEN
+      DEALLOCATE (object, STAT=ierrstat)
+      IF (ierrstat /= SUCCESS) CALL finish("DO_DEALLOCATE_s4D", "DEALLOCATE failed!")
+    END IF
+  END SUBROUTINE DO_DEALLOCATE_s4D
+
+  SUBROUTINE DO_DEALLOCATE_s3D(object)
+    REAL(sp), ALLOCATABLE, INTENT(INOUT) :: object(:, :, :)
+    INTEGER :: ierrstat
+    IF (ALLOCATED(object)) THEN
+      DEALLOCATE (object, STAT=ierrstat)
+      IF (ierrstat /= SUCCESS) CALL finish("DO_DEALLOCATE_s3D", "DEALLOCATE failed!")
+    END IF
+  END SUBROUTINE DO_DEALLOCATE_s3D
+
+  SUBROUTINE DO_DEALLOCATE_s2D(object)
+    REAL(sp), ALLOCATABLE, INTENT(INOUT) :: object(:, :)
+    INTEGER :: ierrstat
+    IF (ALLOCATED(object)) THEN
+      DEALLOCATE (object, STAT=ierrstat)
+      IF (ierrstat /= SUCCESS) CALL finish("DO_DEALLOCATE_s2D", "DEALLOCATE failed!")
+    END IF
+  END SUBROUTINE DO_DEALLOCATE_s2D
+
+  SUBROUTINE DO_DEALLOCATE_s1D(object)
+    REAL(sp), ALLOCATABLE, INTENT(INOUT) :: object(:)
+    INTEGER :: ierrstat
+    IF (ALLOCATED(object)) THEN
+      DEALLOCATE (object, STAT=ierrstat)
+      IF (ierrstat /= SUCCESS) CALL finish("DO_DEALLOCATE_s1D", "DEALLOCATE failed!")
+    END IF
+  END SUBROUTINE DO_DEALLOCATE_s1D
+
   SUBROUTINE DO_DEALLOCATE_i3D(object)
     INTEGER, ALLOCATABLE, INTENT(INOUT) :: object(:, :, :)
     INTEGER :: ierrstat
diff --git a/src/mo_iconlib_kind.F90 b/src/mo_iconlib_kind.F90
new file mode 100644
index 0000000000000000000000000000000000000000..8f470d7b9d27e30adb10163de0bef2485680cd99
--- /dev/null
+++ b/src/mo_iconlib_kind.F90
@@ -0,0 +1,47 @@
+! ICON
+!
+! ---------------------------------------------------------------
+! Copyright (C) 2004-2025, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
+! Contact information: icon-model.org
+!
+! See AUTHORS.TXT for a list of authors
+! See LICENSES/ for license information
+! SPDX-License-Identifier: BSD-3-Clause
+! ---------------------------------------------------------------
+
+!>
+!!   Contains real and integer kinds for libfortran-support and libiconmath
+!!
+MODULE mo_iconlib_kind
+
+  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: real64, real32, int32, int64
+
+  IMPLICIT NONE
+
+  PRIVATE
+
+  PUBLIC :: wp !< Selected working precision
+  PUBLIC :: vp !< Selected variable precision
+  PUBLIC :: dp !< 8 byte real (double precision)
+  PUBLIC :: sp !< 4 byte real (single precision)
+  PUBLIC :: i8 !< 8 byte integer
+  PUBLIC :: i4 !< 4 byte integer
+
+  INTEGER, PARAMETER :: dp = real64
+  INTEGER, PARAMETER :: sp = real32
+  INTEGER, PARAMETER :: i8 = int64
+  INTEGER, PARAMETER :: i4 = int32
+
+#ifdef __SINGLE_PRECISION
+  INTEGER, PARAMETER :: wp = sp
+#else
+  INTEGER, PARAMETER :: wp = dp
+#endif
+
+#ifdef __MIXED_PRECISION
+  INTEGER, PARAMETER :: vp = sp
+#else
+  INTEGER, PARAMETER :: vp = wp
+#endif
+
+END MODULE mo_iconlib_kind
diff --git a/src/mo_octree.F90 b/src/mo_octree.F90
index 33b81ae4c7d4876cb6df742b3ad8216821361323..7cc9c547c76f2a4a5aca0188b60fd41c3f82904b 100644
--- a/src/mo_octree.F90
+++ b/src/mo_octree.F90
@@ -14,7 +14,7 @@
 
 MODULE mo_octree
 
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+  USE mo_iconlib_kind, ONLY: wp
   USE mo_exception, ONLY: finish
   USE mo_util_sort, ONLY: quicksort
 #ifdef __SX__
diff --git a/src/mo_simple_dump.F90 b/src/mo_simple_dump.F90
index 87f0441fe05192eea635066f45f848860faea595..dbbb3a53780708f4c5ced0b907efc1fd6f06283c 100644
--- a/src/mo_simple_dump.F90
+++ b/src/mo_simple_dump.F90
@@ -11,7 +11,7 @@
 
 MODULE mo_simple_dump
 
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+  USE mo_iconlib_kind, ONLY: wp
 
   IMPLICIT NONE
 
diff --git a/src/mo_util_file.F90 b/src/mo_util_file.F90
index b0d0fa08ccfb0f7a58ac50da540fbbe6263999ac..3797476716ccd80685451545ecef608fabd377ca 100644
--- a/src/mo_util_file.F90
+++ b/src/mo_util_file.F90
@@ -12,7 +12,7 @@
 MODULE mo_util_file
 
   USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_int, c_char, c_null_char, c_long
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: i8 => int64
+  USE mo_iconlib_kind, ONLY: i8
   USE mo_exception, ONLY: finish
   USE mo_io_units, ONLY: filename_max
 
diff --git a/src/mo_util_sort.F90 b/src/mo_util_sort.F90
index e5cf130b2ab80a862f693d5f7558086d3b50d65f..c0e45df996925039e7544e10e7f9c1b8924c39db 100644
--- a/src/mo_util_sort.F90
+++ b/src/mo_util_sort.F90
@@ -14,7 +14,7 @@
 
 MODULE mo_util_sort
 
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+  USE mo_iconlib_kind, ONLY: wp, dp, sp
 
   IMPLICIT NONE
 
@@ -27,14 +27,16 @@ MODULE mo_util_sort
 
   ! generic interface for in-situ QuickSort sorting routine
   INTERFACE quicksort
-    MODULE PROCEDURE quicksort_real
+    MODULE PROCEDURE quicksort_real_sp
+    MODULE PROCEDURE quicksort_real_dp
     MODULE PROCEDURE quicksort_int
     MODULE PROCEDURE quicksort_permutation_int
     MODULE PROCEDURE quicksort_string
   END INTERFACE quicksort
 #ifdef __SX__
   INTERFACE radixsort
-    MODULE PROCEDURE radixsort_real
+    MODULE PROCEDURE radixsort_real_dp
+    MODULE PROCEDURE radixsort_real_sp
     MODULE PROCEDURE radixsort_int
   END INTERFACE
 #endif
@@ -58,13 +60,13 @@ CONTAINS
   !
   !  Ordering after the sorting process: smallest...largest.
   !
-  RECURSIVE SUBROUTINE quicksort_real(a, permutation, l_in, r_in)
-    REAL(wp), INTENT(INOUT)           :: a(:) !< array for in-situ sorting
+  RECURSIVE SUBROUTINE quicksort_real_sp(a, permutation, l_in, r_in)
+    REAL(sp), INTENT(INOUT)           :: a(:) !< array for in-situ sorting
     INTEGER, INTENT(INOUT), OPTIONAL :: permutation(:) !< (optional) permutation of indices
     INTEGER, INTENT(IN), OPTIONAL :: l_in, r_in !< left, right partition indices
     ! local variables
     INTEGER  :: i, j, l, r, t_p
-    REAL(wp) :: t, v
+    REAL(sp) :: t, v
 
     IF (PRESENT(l_in)) THEN
       l = l_in
@@ -110,7 +112,61 @@ CONTAINS
       CALL quicksort(a, permutation, l, i - 1)
       CALL quicksort(a, permutation, i + 1, r)
     END IF
-  END SUBROUTINE quicksort_real
+  END SUBROUTINE quicksort_real_sp
+
+  RECURSIVE SUBROUTINE quicksort_real_dp(a, permutation, l_in, r_in)
+    REAL(dp), INTENT(INOUT)           :: a(:) !< array for in-situ sorting
+    INTEGER, INTENT(INOUT), OPTIONAL :: permutation(:) !< (optional) permutation of indices
+    INTEGER, INTENT(IN), OPTIONAL :: l_in, r_in !< left, right partition indices
+    ! local variables
+    INTEGER  :: i, j, l, r, t_p
+    REAL(dp) :: t, v
+
+    IF (PRESENT(l_in)) THEN
+      l = l_in
+    ELSE
+      l = 1
+    END IF
+    IF (PRESENT(r_in)) THEN
+      r = r_in
+    ELSE
+      r = SIZE(a, 1)
+    END IF
+    IF (r > l) THEN
+      v = a(r)
+      i = l - 1
+      j = r
+      LOOP: DO
+        CNTLOOP1: DO
+          i = i + 1
+          IF (a(i) >= v) EXIT CNTLOOP1
+        END DO CNTLOOP1
+        CNTLOOP2: DO
+          j = j - 1
+          IF ((a(j) <= v) .OR. (j == 1)) EXIT CNTLOOP2
+        END DO CNTLOOP2
+        t = a(i)
+        a(i) = a(j)
+        a(j) = t
+        IF (PRESENT(permutation)) THEN
+          t_p = permutation(i)
+          permutation(i) = permutation(j)
+          permutation(j) = t_p
+        END IF
+        IF (j <= i) EXIT LOOP
+      END DO LOOP
+      a(j) = a(i)
+      a(i) = a(r)
+      a(r) = t
+      IF (PRESENT(permutation)) THEN
+        permutation(j) = permutation(i)
+        permutation(i) = permutation(r)
+        permutation(r) = t_p
+      END IF
+      CALL quicksort(a, permutation, l, i - 1)
+      CALL quicksort(a, permutation, i + 1, r)
+    END IF
+  END SUBROUTINE quicksort_real_dp
 
 #ifdef __SX__
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -125,24 +181,24 @@ CONTAINS
   !!    perm: permutation state to record the sorting procedure
   !!          Must have the same size as array!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  SUBROUTINE radixsort_real(array, perm)
+  SUBROUTINE radixsort_real_dp(array, perm)
 
     IMPLICIT NONE
 
-    REAL(KIND=wp), DIMENSION(:), INTENT(INOUT) :: array
+    REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: array
     INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: perm
-    INTEGER(KIND=wp), DIMENSION(SIZE(array)) :: iarray
+    INTEGER(KIND=dp), DIMENSION(SIZE(array)) :: iarray
 
-    INTEGER(KIND=wp), DIMENSION(SIZE(array)) :: bucket00, bucket01, bucket10, bucket11
+    INTEGER(KIND=dp), DIMENSION(SIZE(array)) :: bucket00, bucket01, bucket10, bucket11
     INTEGER, DIMENSION(SIZE(array)) :: pbucket00, pbucket01, pbucket10, pbucket11
     INTEGER :: idx00, idx01, idx10, idx11
 
-    INTEGER(KIND=wp), PARAMETER :: bitmaskr2i = -1
-    ! REAL(KIND=wp), PARAMETER :: bitmaski2r = Z'FFFFFFFFFFFFFFFF'
+    INTEGER(KIND=dp), PARAMETER :: bitmaskr2i = -1
+    ! REAL(KIND=dp), PARAMETER :: bitmaski2r = Z'FFFFFFFFFFFFFFFF'
 
     INTEGER :: i, j
     INTEGER :: n, offs
-    INTEGER(KIND=wp) :: tmp
+    INTEGER(KIND=dp) :: tmp
 
     n = SIZE(array)
 
@@ -286,7 +342,183 @@ CONTAINS
     array = TRANSFER(iarray, Z'FFFFFFFFFFFFFFFF', n)
 
     RETURN
-  END SUBROUTINE radixsort_real
+  END SUBROUTINE radixsort_real_dp
+
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !! Description:
+  !!    Sorts an array of type sp with a radix sort using two bits at once.
+  !!    Doubles the required memory,
+  !!    Performance of sp version not tested
+  !!
+  !!    Scales linearly with the number of elements for large number of elements.
+  !!    The difference between initial and final state is recorded by a permutation state perm
+  !! Variables:
+  !!    array: The sp-array to be sorted
+  !!    perm: permutation state to record the sorting procedure
+  !!          Must have the same size as array!
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  SUBROUTINE radixsort_real_sp(array, perm)
+
+    IMPLICIT NONE
+
+    REAL(KIND=sp), DIMENSION(:), INTENT(INOUT) :: array
+    INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: perm
+    INTEGER(KIND=sp), DIMENSION(SIZE(array)) :: iarray
+
+    INTEGER(KIND=sp), DIMENSION(SIZE(array)) :: bucket00, bucket01, bucket10, bucket11
+    INTEGER, DIMENSION(SIZE(array)) :: pbucket00, pbucket01, pbucket10, pbucket11
+    INTEGER :: idx00, idx01, idx10, idx11
+
+    INTEGER(KIND=sp), PARAMETER :: bitmaskr2i = -1
+    ! REAL(KIND=sp), PARAMETER :: bitmaski2r = Z'FFFFFFFF'
+
+    INTEGER :: i, j
+    INTEGER :: n, offs
+    INTEGER(KIND=sp) :: tmp
+
+    n = SIZE(array)
+
+    ! Transfer the bit pattern of the real array to an i8 array
+    iarray = TRANSFER(array, bitmaskr2i, n)
+
+    ! Loop over the bits
+    DO i = 0, STORAGE_SIZE(tmp) - 2, 2
+      idx00 = 0
+      idx01 = 0
+      idx10 = 0
+      idx11 = 0
+      ! sort numbers into buckets based on their bit i
+      DO j = 1, n
+        ! extract the i-th bit
+        tmp = IBITS(iarray(j), i, 2)
+        IF (tmp < 2) THEN
+          IF (tmp == 0) THEN
+            !increase the counter of the selected bucket
+            idx00 = idx00 + 1
+            ! add number to bucket
+            bucket00(idx00) = iarray(j)
+            ! same operation on the permutation
+            IF (PRESENT(perm)) pbucket00(idx00) = perm(j)
+          ELSE
+            !increase the counter of the selected bucket
+            idx01 = idx01 + 1
+            ! add number to bucket
+            bucket01(idx01) = iarray(j)
+            ! same operation on the permutation
+            IF (PRESENT(perm)) pbucket01(idx01) = perm(j)
+          END IF
+        ELSE
+          IF (tmp == 2) THEN
+            !increase the counter of the selected bucket
+            idx10 = idx10 + 1
+            ! add number to bucket
+            bucket10(idx10) = iarray(j)
+            ! same operation on the permutation
+            IF (PRESENT(perm)) pbucket10(idx10) = perm(j)
+          ELSE
+            !increase the counter of the selected bucket
+            idx11 = idx11 + 1
+            ! add number to bucket
+            bucket11(idx11) = iarray(j)
+            ! same operation on the permutation
+            IF (PRESENT(perm)) pbucket11(idx11) = perm(j)
+          END IF
+        END IF
+      END DO
+
+      ! copy the presorted numbers back onto the array
+      DO j = 1, idx00
+        iarray(j) = bucket00(j)
+        IF (PRESENT(perm)) perm(j) = pbucket00(j)
+      END DO
+      offs = idx00
+      DO j = 1, idx01
+        iarray(offs + j) = bucket01(j)
+        IF (PRESENT(perm)) perm(offs + j) = pbucket01(j)
+      END DO
+      offs = offs + idx01
+      DO j = 1, idx10
+        iarray(offs + j) = bucket10(j)
+        IF (PRESENT(perm)) perm(offs + j) = pbucket10(j)
+      END DO
+      offs = offs + idx10
+      DO j = 1, idx11
+        iarray(offs + j) = bucket11(j)
+        IF (PRESENT(perm)) perm(offs + j) = pbucket11(j)
+      END DO
+    END DO
+
+    ! sort by sign
+    idx00 = 0
+    idx01 = 0
+    idx10 = 0
+    idx11 = 0
+    i = STORAGE_SIZE(tmp) - 2
+    ! sort numbers into buckets based on their first bit and
+    DO j = 1, n
+      tmp = IBITS(iarray(j), i, 2)
+      IF (tmp < 2) THEN
+        IF (tmp == 0) THEN
+          !increase the counter of the selected bucket
+          idx00 = idx00 + 1
+          ! add number to bucket
+          bucket00(idx00) = iarray(j)
+          ! same operation on the permutation
+          IF (PRESENT(perm)) pbucket00(idx00) = perm(j)
+        ELSE
+          !increase the counter of the selected bucket
+          idx01 = idx01 + 1
+          ! add number to bucket
+          bucket01(idx01) = iarray(j)
+          ! same operation on the permutation
+          IF (PRESENT(perm)) pbucket01(idx01) = perm(j)
+        END IF
+      ELSE
+        IF (tmp == 2) THEN
+          !increase the counter of the selected bucket
+          idx10 = idx10 + 1
+          ! add number to bucket
+          bucket10(idx10) = iarray(j)
+          ! same operation on the permutation
+          IF (PRESENT(perm)) pbucket10(idx10) = perm(j)
+        ELSE
+          !increase the counter of the selected bucket
+          idx11 = idx11 + 1
+          ! add number to bucket
+          bucket11(idx11) = iarray(j)
+          ! same operation on the permutation
+          IF (PRESENT(perm)) pbucket11(idx11) = perm(j)
+        END IF
+      END IF
+    END DO
+
+    ! copy the presorted numbers back onto the array
+    ! the half inverse order is a result of the ieee 745 float sign bit standard
+    DO j = 1, idx11
+      iarray(idx11 - j + 1) = bucket11(j)
+      IF (PRESENT(perm)) perm(idx11 - j + 1) = pbucket11(j)
+    END DO
+    offs = idx11
+    DO j = 1, idx10
+      iarray(offs + idx10 - j + 1) = bucket10(j)
+      IF (PRESENT(perm)) perm(offs + idx10 - j + 1) = pbucket10(j)
+    END DO
+    offs = offs + idx10
+    DO j = 1, idx00
+      iarray(offs + j) = bucket00(j)
+      IF (PRESENT(perm)) perm(offs + j) = pbucket00(j)
+    END DO
+    offs = offs + idx00
+    DO j = 1, idx01
+      iarray(offs + j) = bucket01(j)
+      IF (PRESENT(perm)) perm(offs + j) = pbucket01(j)
+    END DO
+
+    ! Transfer the bit pattern of the i4 array to the r4 array
+    array = TRANSFER(iarray, Z'FFFFFFFF', n)
+
+    RETURN
+  END SUBROUTINE radixsort_real_sp
 #endif
   SUBROUTINE swap_int(a, i, j)
     !> array for in-situ sorting
diff --git a/test/fortran/CMakeLists.txt b/test/fortran/CMakeLists.txt
index a95e5d0bd0d8bc2d2d33d06258ed25222b651cf9..8ec99d582220ecbcb31f4b94636a7b7d9365e4c6 100644
--- a/test/fortran/CMakeLists.txt
+++ b/test/fortran/CMakeLists.txt
@@ -34,3 +34,9 @@ FortUTF_Find_Tests()
 set(PROJECT_NAME ${name_holder})
 add_test(NAME FortUTF_UnitTest COMMAND fortran_support_Tests)
 set_property(TEST FortUTF_UnitTest PROPERTY LABELS Fortran)
+
+# preprocessing is required for precision logic in testfile
+set_source_files_properties(
+  test_fortran_tools.f90
+  PROPERTIES Fortran_PREPROCESS ON
+)
diff --git a/test/fortran/helpers.f90 b/test/fortran/helpers.f90
index 5a65dbab0a3d3c05d9de170e482d1ffadb3f8a14..f77f27c270b15550dc654bf4d38b3f6d13719ff8 100644
--- a/test/fortran/helpers.f90
+++ b/test/fortran/helpers.f90
@@ -12,7 +12,6 @@
 MODULE helpers
   USE ISO_C_BINDING, ONLY: c_double
   USE fortran_support, ONLY: find_next_free_unit, message
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
 
 CONTAINS
 
diff --git a/test/fortran/test_exception.f90 b/test/fortran/test_exception.f90
index 544b0867add518cc271e6ca6d7f7922767864161..27f74e88b868d0a3c6b104432c4a314ad3270e1e 100644
--- a/test/fortran/test_exception.f90
+++ b/test/fortran/test_exception.f90
@@ -11,8 +11,7 @@
 
 MODULE test_mo_exception
   USE FORTUTF
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64, &
-    &                                      i8 => int64
+  USE mo_iconlib_kind, ONLY: dp, sp, i8
   ! From mo_exception
   USE fortran_support, ONLY: init_logger, message, message_to_own_unit, &
     & debug_on, debug_off, warning, print_value, finish, set_msg_timestamp, &
@@ -160,7 +159,8 @@ CONTAINS
     CALL print_value('-', .TRUE.)
     CALL print_value('-', 6)
     CALL print_value('-', 6_i8)
-    CALL print_value('-', 6.0_wp)
+    CALL print_value('-', 6.0_dp)
+    CALL print_value('-', 6.0_sp)
 
     CLOSE (nerr)
 
@@ -179,7 +179,11 @@ CONTAINS
     CALL STRING_CONTAINS(log_in_file, '- :         6')
 
     READ (nerr, '(A)') log_in_file
-    CALL TAG_TEST("TEST_param_real")
+    CALL TAG_TEST("TEST_param_dp")
+    CALL STRING_CONTAINS(log_in_file, '- :  6.0000')
+
+    READ (nerr, '(A)') log_in_file
+    CALL TAG_TEST("TEST_param_sp")
     CALL STRING_CONTAINS(log_in_file, '- :  6.0000')
 
     CLOSE (nerr)
diff --git a/test/fortran/test_expression.f90 b/test/fortran/test_expression.f90
index 9f3f81ca401ca8fa4f2a69e53704dc11927d4540..7d4b20bb0f6ae2509f74c45e995ff7f02a7caa28 100644
--- a/test/fortran/test_expression.f90
+++ b/test/fortran/test_expression.f90
@@ -12,7 +12,7 @@
 MODULE test_mo_expression
   USE FORTUTF
   USE fortran_support, ONLY: expression
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+  USE mo_iconlib_kind, ONLY: wp
 
 CONTAINS
 
diff --git a/test/fortran/test_fortran_tools.f90 b/test/fortran/test_fortran_tools.f90
index 8a83902c2a2333bf92750170ed7e93b70d77a4a7..36012a270f65223a849a45f2ca105bed57bb5d3d 100644
--- a/test/fortran/test_fortran_tools.f90
+++ b/test/fortran/test_fortran_tools.f90
@@ -11,16 +11,18 @@
 
 MODULE test_fortran_tools
   USE FORTUTF
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: dp => real64, &
-    &                                      sp => real32, &
-    &                                      i4 => int32
+  USE mo_iconlib_kind, ONLY: wp, dp, sp, i4
+
   ! From mo_fortran_tools
   USE fortran_support, ONLY: assign_if_present, assign_if_present_allocatable, &
     & if_associated, copy, init, swap, negative2zero, var_scale, var_add, &
     & init_zero_contiguous_dp, init_zero_contiguous_sp, init_contiguous_dp, &
     & init_contiguous_sp, init_contiguous_i4, init_contiguous_l, minval_1d, &
     & minval_2d, resize_arr_c1d, DO_DEALLOCATE, DO_PTR_DEALLOCATE, &
-    & insert_dimension
+    & insert_dimension, &
+    & t_ptr_2d_wp, t_ptr_2d_dp, t_ptr_2d_sp, &
+    & t_ptr_3d_wp, t_ptr_3d_dp, t_ptr_3d_sp, &
+    & t_ptr_4d_wp, t_ptr_4d_dp, t_ptr_4d_sp
   ! From other modules
   USE fortran_support, ONLY: init_logger, finish
   USE helpers, ONLY: open_logfile, open_new_logfile, custom_exit_dummy
@@ -347,6 +349,58 @@ CONTAINS
     CALL ASSERT_EQUAL(assert_real_5d_array(src, dest), .TRUE.)
   END SUBROUTINE
 
+  SUBROUTINE Test_copy_1d_sp
+    REAL(sp) :: src(10) = 1.0, dest(10)
+
+    CALL TAG_TEST("Test_copy_1d_sp_ones")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_array(src, dest), .TRUE.)
+
+    CALL RANDOM_NUMBER(src)
+    CALL TAG_TEST("Test_copy_1d_sp_random")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_array(src, dest), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_copy_2d_sp
+    REAL(sp) :: src(10, 10) = 1.0, dest(10, 10)
+
+    CALL TAG_TEST("Test_copy_2d_sp_ones")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_2d_array(src, dest), .TRUE.)
+
+    CALL RANDOM_NUMBER(src)
+    CALL TAG_TEST("Test_copy_2d_sp_random")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_2d_array(src, dest), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_copy_3d_sp
+    REAL(sp) :: src(10, 10, 10) = 1.0, dest(10, 10, 10)
+
+    CALL TAG_TEST("Test_copy_3d_sp_ones")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_3d_array(src, dest), .TRUE.)
+
+    CALL RANDOM_NUMBER(src)
+    CALL TAG_TEST("Test_copy_3d_sp_random")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_3d_array(src, dest), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_copy_4d_sp
+    REAL(sp) :: src(5, 5, 5, 5) = 1.0, dest(5, 5, 5, 5)
+
+    CALL TAG_TEST("Test_copy_4d_sp_ones")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_4d_array(src, dest), .TRUE.)
+
+    CALL RANDOM_NUMBER(src)
+    CALL TAG_TEST("Test_copy_4d_sp_random")
+    CALL copy(src, dest, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_4d_array(src, dest), .TRUE.)
+  END SUBROUTINE
+
   SUBROUTINE Test_copy_5d_sp
     REAL(sp) :: src(5, 5, 5, 5, 5) = 1.0, dest(5, 5, 5, 5, 5)
 
@@ -500,6 +554,14 @@ CONTAINS
     CALL ASSERT_EQUAL(assert_real_2d_array(arr, zeros), .TRUE.)
   END SUBROUTINE
 
+  SUBROUTINE Test_init_zero_2d_sp
+    REAL(sp) :: arr(10, 10), zeros(10, 10) = 0.0
+
+    CALL TAG_TEST("Test_init_zero_2d_sp")
+    CALL init(arr, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_2d_array(arr, zeros), .TRUE.)
+  END SUBROUTINE
+
   SUBROUTINE Test_init_zero_2d_i4
     INTEGER(i4) :: arr(10, 10), zeros(10, 10) = 0
 
@@ -589,6 +651,30 @@ CONTAINS
     CALL ASSERT_EQUAL(assert_real_spdp_3d_array(arr, ones), .TRUE.)
   END SUBROUTINE
 
+  SUBROUTINE Test_init_1d_sp
+    REAL(sp) :: arr(10), ones(10) = 1.0
+
+    CALL TAG_TEST("Test_init_1d_sp")
+    CALL init(arr, 1.0_sp, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_array(arr, ones), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_init_2d_sp
+    REAL(sp) :: arr(10, 10), ones(10, 10) = 1.0
+
+    CALL TAG_TEST("Test_init_2d_sp")
+    CALL init(arr, 1.0_sp, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_2d_array(arr, ones), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_init_3d_sp
+    REAL(sp) :: arr(10, 10, 10), ones(10, 10, 10) = 1.0
+
+    CALL TAG_TEST("Test_init_3d_sp")
+    CALL init(arr, 1.0_sp, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_3d_array(arr, ones), .TRUE.)
+  END SUBROUTINE
+
   SUBROUTINE Test_init_5d_dp
     REAL(dp) :: arr(5, 5, 5, 5, 5), ones(5, 5, 5, 5, 5) = 1.0
 
@@ -601,7 +687,7 @@ CONTAINS
     REAL(sp) :: arr(5, 5, 5, 5, 5), ones(5, 5, 5, 5, 5) = 1.0
 
     CALL TAG_TEST("Test_init_5d_sp")
-    CALL init(arr, 1.0, .FALSE.)
+    CALL init(arr, 1.0_sp, .FALSE.)
     CALL ASSERT_EQUAL(assert_real_sp_5d_array(arr, ones), .TRUE.)
   END SUBROUTINE
 
@@ -621,15 +707,24 @@ CONTAINS
     CALL ASSERT_EQUAL(assert_logical_5d_array(arr, trues), .TRUE.)
   END SUBROUTINE
 
-  SUBROUTINE Test_var_scale_3d
+  SUBROUTINE Test_var_scale_3d_dp
     REAL(dp) :: arr(10, 10, 10) = 1.0, scale = 5.0
     REAL(dp) :: ans(10, 10, 10) = 5.0
 
-    CALL TAG_TEST("Test_var_scale_3d")
+    CALL TAG_TEST("Test_var_scale_3d_dp")
     CALL var_scale(arr, scale, .FALSE.)
     CALL ASSERT_EQUAL(assert_real_3d_array(arr, ans), .TRUE.)
   END SUBROUTINE
 
+  SUBROUTINE Test_var_scale_3d_sp
+    REAL(sp) :: arr(10, 10, 10) = 1.0, scale = 5.0
+    REAL(sp) :: ans(10, 10, 10) = 5.0
+
+    CALL TAG_TEST("Test_var_scale_3d_sp")
+    CALL var_scale(arr, scale, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_3d_array(arr, ans), .TRUE.)
+  END SUBROUTINE
+
   SUBROUTINE Test_var_addc_3d_dp
     REAL(dp) :: arr(10, 10, 10) = 2.0, const = 3.0
     REAL(dp) :: ans(10, 10, 10) = 5.0
@@ -664,6 +759,31 @@ CONTAINS
     CALL ASSERT_EQUAL(assert_real_4d_array(arr, ans), .TRUE.)
   END SUBROUTINE
 
+  SUBROUTINE Test_negative2zero_4d_sp
+    REAL(sp) :: arr(5, 5, 5, 5), ans(5, 5, 5, 5)
+    INTEGER :: i, j, k, l
+
+    CALL RANDOM_NUMBER(arr)
+    DO i = 1, 5
+      DO j = 1, 5
+        DO k = 1, 5
+          DO l = 1, 5
+            arr(i, j, k, l) = arr(i, j, k, l) - 0.5
+            IF (arr(i, j, k, l) < 0.0) THEN
+              ans(i, j, k, l) = 0.0
+            ELSE
+              ans(i, j, k, l) = arr(i, j, k, l)
+            END IF
+          END DO
+        END DO
+      END DO
+    END DO
+
+    CALL TAG_TEST("Test_negative2zero_4d_sp")
+    CALL negative2zero(arr, .FALSE.)
+    CALL ASSERT_EQUAL(assert_real_sp_4d_array(arr, ans), .TRUE.)
+  END SUBROUTINE
+
   SUBROUTINE Test_init_contiguous_dp
     REAL(dp) :: arr(10), ans(10) = 7.0
 
@@ -1406,6 +1526,42 @@ CONTAINS
     CALL ASSERT_EQUAL(ALLOCATED(arr), .FALSE.)
   END SUBROUTINE
 
+  SUBROUTINE Test_DO_DEALLOCATE_s4D
+    REAL(sp), ALLOCATABLE :: arr(:, :, :, :)
+
+    ALLOCATE (arr(5, 5, 5, 5))
+    CALL TAG_TEST("Test_DO_DEALLOCATE_s4D")
+    CALL DO_DEALLOCATE(arr)
+    CALL ASSERT_EQUAL(ALLOCATED(arr), .FALSE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_DO_DEALLOCATE_s3D
+    REAL(sp), ALLOCATABLE :: arr(:, :, :)
+
+    ALLOCATE (arr(10, 10, 10))
+    CALL TAG_TEST("Test_DO_DEALLOCATE_s3D")
+    CALL DO_DEALLOCATE(arr)
+    CALL ASSERT_EQUAL(ALLOCATED(arr), .FALSE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_DO_DEALLOCATE_s2D
+    REAL(sp), ALLOCATABLE :: arr(:, :)
+
+    ALLOCATE (arr(10, 10))
+    CALL TAG_TEST("Test_DO_DEALLOCATE_s2D")
+    CALL DO_DEALLOCATE(arr)
+    CALL ASSERT_EQUAL(ALLOCATED(arr), .FALSE.)
+  END SUBROUTINE
+
+  SUBROUTINE Test_DO_DEALLOCATE_s1D
+    REAL(sp), ALLOCATABLE :: arr(:)
+
+    ALLOCATE (arr(10))
+    CALL TAG_TEST("Test_DO_DEALLOCATE_s1D")
+    CALL DO_DEALLOCATE(arr)
+    CALL ASSERT_EQUAL(ALLOCATED(arr), .FALSE.)
+  END SUBROUTINE
+
   SUBROUTINE Test_DO_DEALLOCATE_i3D
     INTEGER, ALLOCATABLE :: arr(:, :, :)
 
@@ -1483,7 +1639,7 @@ CONTAINS
     CALL ASSERT_EQUAL(ASSOCIATED(ptr), .FALSE.)
   END SUBROUTINE
 
-! Support functions for testing
+  ! Support functions for testing
 
   LOGICAL FUNCTION assert_logical_array(array1, array2)
     LOGICAL, INTENT(IN) :: array1(:), array2(:)
@@ -1609,6 +1765,21 @@ CONTAINS
     END DO
   END FUNCTION assert_real_sp_array
 
+  LOGICAL FUNCTION assert_real_sp_2d_array(array1, array2)
+    REAL(sp), INTENT(IN) :: array1(:, :), array2(:, :)
+    INTEGER :: i, j
+
+    assert_real_sp_2d_array = .TRUE.
+    DO i = 1, SIZE(array1, 1)
+      DO j = 1, SIZE(array1, 2)
+        IF (array1(i, j) /= array2(i, j)) THEN
+          assert_real_sp_2d_array = .FALSE.
+          EXIT
+        END IF
+      END DO
+    END DO
+  END FUNCTION assert_real_sp_2d_array
+
   LOGICAL FUNCTION assert_real_sp_3d_array(array1, array2)
     REAL(sp), INTENT(IN) :: array1(:, :, :), array2(:, :, :)
     INTEGER :: i, j, k
@@ -1835,4 +2006,121 @@ CONTAINS
     END DO
   END FUNCTION assert_logical_5d_array
 
+  ! Intentionally reverse t_ptr_2d_dp and t_ptr_2d_wp when passing inputs to test compiler
+  LOGICAL FUNCTION assert_t_ptr_2d_tmp(ptr1, ptr2)
+    ! See: https://stackoverflow.com/a/33304488
+#ifdef __SINGLE_PRECISION
+    TYPE(t_ptr_2d_sp), INTENT(IN) :: ptr1
+#else
+    TYPE(t_ptr_2d_dp), INTENT(IN) :: ptr1
+#endif
+    TYPE(t_ptr_2d_wp), INTENT(IN) :: ptr2
+    INTEGER :: i, j
+
+    assert_t_ptr_2d_tmp = .TRUE.
+    IF (ASSOCIATED(ptr1%p) .OR. ASSOCIATED(ptr2%p)) THEN
+      assert_t_ptr_2d_tmp = ASSOCIATED(ptr1%p, ptr2%p)
+    END IF
+  END FUNCTION assert_t_ptr_2d_tmp
+
+  LOGICAL FUNCTION assert_t_ptr_2d(ptr1, ptr2)
+    ! See: https://stackoverflow.com/a/33304488
+#ifdef __SINGLE_PRECISION
+    TYPE(t_ptr_2d_sp), INTENT(IN) :: ptr1
+#else
+    TYPE(t_ptr_2d_dp), INTENT(IN) :: ptr1
+#endif
+    TYPE(t_ptr_2d_wp), INTENT(IN) :: ptr2
+    INTEGER :: i, j
+
+    assert_t_ptr_2d = .TRUE.
+    IF (ASSOCIATED(ptr1%p) .OR. ASSOCIATED(ptr2%p)) THEN
+      assert_t_ptr_2d = ASSOCIATED(ptr1%p, ptr2%p)
+    END IF
+  END FUNCTION assert_t_ptr_2d
+
+  LOGICAL FUNCTION assert_t_ptr_2d_arr(ptr1, ptr2)
+#ifdef __SINGLE_PRECISION
+    TYPE(t_ptr_2d_sp), INTENT(IN) :: ptr1(:)
+#else
+    TYPE(t_ptr_2d_dp), INTENT(IN) :: ptr1(:)
+#endif
+    TYPE(t_ptr_2d_wp), INTENT(IN) :: ptr2(:)
+    INTEGER :: i
+
+    assert_t_ptr_2d_arr = .TRUE.
+    DO i = 1, SIZE(ptr1)
+      IF (assert_t_ptr_2d(ptr1(i), ptr2(i)) .NEQV. .TRUE.) THEN
+        assert_t_ptr_2d_arr = .FALSE.
+        EXIT
+      END IF
+    END DO
+  END FUNCTION assert_t_ptr_2d_arr
+
+  LOGICAL FUNCTION assert_t_ptr_3d(ptr1, ptr2)
+#ifdef __SINGLE_PRECISION
+    TYPE(t_ptr_3d_sp), INTENT(IN) :: ptr1
+#else
+    TYPE(t_ptr_3d_dp), INTENT(IN) :: ptr1
+#endif
+    TYPE(t_ptr_3d_wp), INTENT(IN) :: ptr2
+    INTEGER :: i, j, k
+
+    assert_t_ptr_3d = .TRUE.
+    IF (ASSOCIATED(ptr1%p) .OR. ASSOCIATED(ptr2%p)) THEN
+      assert_t_ptr_3d = ASSOCIATED(ptr1%p, ptr2%p)
+    END IF
+  END FUNCTION assert_t_ptr_3d
+
+  LOGICAL FUNCTION assert_t_ptr_3d_arr(ptr1, ptr2)
+#ifdef __SINGLE_PRECISION
+    TYPE(t_ptr_3d_sp), INTENT(IN) :: ptr1(:)
+#else
+    TYPE(t_ptr_3d_dp), INTENT(IN) :: ptr1(:)
+#endif
+    TYPE(t_ptr_3d_wp), INTENT(IN) :: ptr2(:)
+    INTEGER :: i
+
+    assert_t_ptr_3d_arr = .TRUE.
+    DO i = 1, SIZE(ptr1)
+      IF (assert_t_ptr_3d(ptr1(i), ptr2(i)) .NEQV. .TRUE.) THEN
+        assert_t_ptr_3d_arr = .FALSE.
+        EXIT
+      END IF
+    END DO
+  END FUNCTION assert_t_ptr_3d_arr
+
+  LOGICAL FUNCTION assert_t_ptr_4d(ptr1, ptr2)
+#ifdef __SINGLE_PRECISION
+    TYPE(t_ptr_4d_sp), INTENT(IN) :: ptr1
+#else
+    TYPE(t_ptr_4d_dp), INTENT(IN) :: ptr1
+#endif
+    TYPE(t_ptr_4d_wp), INTENT(IN) :: ptr2
+    INTEGER :: i, j, k
+
+    assert_t_ptr_4d = .TRUE.
+    IF (ASSOCIATED(ptr1%p) .OR. ASSOCIATED(ptr2%p)) THEN
+      assert_t_ptr_4d = ASSOCIATED(ptr1%p, ptr2%p)
+    END IF
+  END FUNCTION assert_t_ptr_4d
+
+  LOGICAL FUNCTION assert_t_ptr_4d_arr(ptr1, ptr2)
+#ifdef __SINGLE_PRECISION
+    TYPE(t_ptr_4d_sp), INTENT(IN) :: ptr1(:)
+#else
+    TYPE(t_ptr_4d_dp), INTENT(IN) :: ptr1(:)
+#endif
+    TYPE(t_ptr_4d_wp), INTENT(IN) :: ptr2(:)
+    INTEGER :: i
+
+    assert_t_ptr_4d_arr = .TRUE.
+    DO i = 1, SIZE(ptr1)
+      IF (assert_t_ptr_4d(ptr1(i), ptr2(i)) .NEQV. .TRUE.) THEN
+        assert_t_ptr_4d_arr = .FALSE.
+        EXIT
+      END IF
+    END DO
+  END FUNCTION assert_t_ptr_4d_arr
+
 END MODULE
diff --git a/test/fortran/test_sort.f90 b/test/fortran/test_sort.f90
index d1d62ad9d13e3781e2fb9f4f1d0ecee45a25b211..1d03a24a7e35ceb6f6429e83513386e1ee56034c 100644
--- a/test/fortran/test_sort.f90
+++ b/test/fortran/test_sort.f90
@@ -11,55 +11,103 @@
 
 MODULE TEST_mo_util_sort
   USE FORTUTF
-  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+  USE mo_iconlib_kind, ONLY: wp, dp, sp
   USE fortran_support, ONLY: quicksort, insertion_sort
 
 CONTAINS
-  SUBROUTINE TEST_quicksort_real
-    REAL(wp) :: to_sort(6) = (/144.4, 58.6, 4.3, 7.8, 10.0, 11.0/)
-    CALL TAG_TEST("TEST_quicksort_real_before")
-    CALL ASSERT_EQUAL(is_sorted_real(to_sort), .FALSE.)
+  SUBROUTINE TEST_quicksort_real_dp
+    REAL(dp) :: to_sort(6) = (/144.4, 58.6, 4.3, 7.8, 10.0, 11.0/)
+    CALL TAG_TEST("TEST_quicksort_real_dp_before")
+    CALL ASSERT_EQUAL(is_sorted_real_dp(to_sort), .FALSE.)
 
     CALL quicksort(to_sort)
 
-    CALL TAG_TEST("TEST_quicksort_real_after")
-    CALL ASSERT_EQUAL(is_sorted_real(to_sort), .TRUE.)
+    CALL TAG_TEST("TEST_quicksort_real_dp_after")
+    CALL ASSERT_EQUAL(is_sorted_real_dp(to_sort), .TRUE.)
   END SUBROUTINE
 
-  SUBROUTINE TEST_quicksort_real2
-    REAL(wp) :: to_sort(6) = (/144.4, 11.0, 4.3, 58.6, 10.0, 7.8/)
-    CALL TAG_TEST("TEST_quicksort_real2_before")
-    CALL ASSERT_EQUAL(is_sorted_real(to_sort), .FALSE.)
+  SUBROUTINE TEST_quicksort_real_dp2
+    REAL(dp) :: to_sort(6) = (/144.4, 11.0, 4.3, 58.6, 10.0, 7.8/)
+    CALL TAG_TEST("TEST_quicksort_real_dp2_before")
+    CALL ASSERT_EQUAL(is_sorted_real_dp(to_sort), .FALSE.)
 
     CALL quicksort(to_sort)
 
-    CALL TAG_TEST("TEST_quicksort_real2_after")
-    CALL ASSERT_EQUAL(is_sorted_real(to_sort), .TRUE.)
+    CALL TAG_TEST("TEST_quicksort_real_dp2_after")
+    CALL ASSERT_EQUAL(is_sorted_real_dp(to_sort), .TRUE.)
   END SUBROUTINE
 
-  SUBROUTINE TEST_quicksort_real_random
-    REAL(wp) :: to_sort(6)
+  SUBROUTINE TEST_quicksort_real_dp_random
+    REAL(dp) :: to_sort(6)
     ! Generate random numbers geq 0.0 and < 256.0
     CALL RANDOM_NUMBER(to_sort)
     to_sort = to_sort*256
 
     CALL quicksort(to_sort)
 
-    CALL TAG_TEST("TEST_quicksort_real_random")
-    CALL ASSERT_EQUAL(is_sorted_real(to_sort), .TRUE.)
+    CALL TAG_TEST("TEST_quicksort_real_dp_random")
+    CALL ASSERT_EQUAL(is_sorted_real_dp(to_sort), .TRUE.)
   END SUBROUTINE
 
-  SUBROUTINE TEST_quicksort_permutation_real
-    REAL(wp) :: to_sort(6) = (/144.4, 58.6, 4.3, 7.8, 10.0, 11.0/)
+  SUBROUTINE TEST_quicksort_permutation_real_dp
+    REAL(dp) :: to_sort(6) = (/144.4, 58.6, 4.3, 7.8, 10.0, 11.0/)
     INTEGER :: idx_permutation(6) = (/1, 2, 3, 4, 5, 6/)
     CALL TAG_TEST("TEST_quicksort_permutation_real_before")
-    CALL ASSERT_EQUAL(is_sorted_real(to_sort), .FALSE.)
+    CALL ASSERT_EQUAL(is_sorted_real_dp(to_sort), .FALSE.)
 
     CALL quicksort(to_sort, idx_permutation)
 
-    CALL TAG_TEST("TEST_quicksort_permutation_real_after")
-    CALL ASSERT_EQUAL(is_sorted_real(to_sort), .TRUE.)
-    CALL TAG_TEST("TEST_quicksort_permutation_real_permutation")
+    CALL TAG_TEST("TEST_quicksort_permutation_real_dp_after")
+    CALL ASSERT_EQUAL(is_sorted_real_dp(to_sort), .TRUE.)
+    CALL TAG_TEST("TEST_quicksort_permutation_real_dp_permutation")
+    CALL ASSERT_EQUAL(has_same_values_int(idx_permutation, (/3, 4, 5, 6, 2, 1/)), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE TEST_quicksort_real_sp
+    REAL(sp) :: to_sort(6) = (/144.4, 58.6, 4.3, 7.8, 10.0, 11.0/)
+    CALL TAG_TEST("TEST_quicksort_real_sp_before")
+    CALL ASSERT_EQUAL(is_sorted_real_sp(to_sort), .FALSE.)
+
+    CALL quicksort(to_sort)
+
+    CALL TAG_TEST("TEST_quicksort_real_sp_after")
+    CALL ASSERT_EQUAL(is_sorted_real_sp(to_sort), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE TEST_quicksort_real_sp2
+    REAL(sp) :: to_sort(6) = (/144.4, 11.0, 4.3, 58.6, 10.0, 7.8/)
+    CALL TAG_TEST("TEST_quicksort_real_sp2_before")
+    CALL ASSERT_EQUAL(is_sorted_real_sp(to_sort), .FALSE.)
+
+    CALL quicksort(to_sort)
+
+    CALL TAG_TEST("TEST_quicksort_real_sp2_after")
+    CALL ASSERT_EQUAL(is_sorted_real_sp(to_sort), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE TEST_quicksort_real_sp_random
+    REAL(sp) :: to_sort(6)
+    ! Generate random numbers geq 0.0 and < 256.0
+    CALL RANDOM_NUMBER(to_sort)
+    to_sort = to_sort*256
+
+    CALL quicksort(to_sort)
+
+    CALL TAG_TEST("TEST_quicksort_real_sp_random")
+    CALL ASSERT_EQUAL(is_sorted_real_sp(to_sort), .TRUE.)
+  END SUBROUTINE
+
+  SUBROUTINE TEST_quicksort_permutation_real_sp
+    REAL(sp) :: to_sort(6) = (/144.4, 58.6, 4.3, 7.8, 10.0, 11.0/)
+    INTEGER :: idx_permutation(6) = (/1, 2, 3, 4, 5, 6/)
+    CALL TAG_TEST("TEST_quicksort_permutation_real_before")
+    CALL ASSERT_EQUAL(is_sorted_real_sp(to_sort), .FALSE.)
+
+    CALL quicksort(to_sort, idx_permutation)
+
+    CALL TAG_TEST("TEST_quicksort_permutation_real_sp_after")
+    CALL ASSERT_EQUAL(is_sorted_real_sp(to_sort), .TRUE.)
+    CALL TAG_TEST("TEST_quicksort_permutation_real_sp_permutation")
     CALL ASSERT_EQUAL(has_same_values_int(idx_permutation, (/3, 4, 5, 6, 2, 1/)), .TRUE.)
   END SUBROUTINE
 
@@ -194,19 +242,33 @@ CONTAINS
 
   END FUNCTION has_same_values_int
 
-  LOGICAL FUNCTION is_sorted_real(array)
-    REAL(wp), INTENT(IN) :: array(:)
+  LOGICAL FUNCTION is_sorted_real_dp(array)
+    REAL(dp), INTENT(IN) :: array(:)
+    INTEGER :: i
+
+    is_sorted_real_dp = .TRUE.
+    DO i = 1, SIZE(array) - 1
+      IF (array(i) > array(i + 1)) THEN
+        is_sorted_real_dp = .FALSE.
+        EXIT
+      END IF
+    END DO
+
+  END FUNCTION is_sorted_real_dp
+
+  LOGICAL FUNCTION is_sorted_real_sp(array)
+    REAL(sp), INTENT(IN) :: array(:)
     INTEGER :: i
 
-    is_sorted_real = .TRUE.
+    is_sorted_real_sp = .TRUE.
     DO i = 1, SIZE(array) - 1
       IF (array(i) > array(i + 1)) THEN
-        is_sorted_real = .FALSE.
+        is_sorted_real_sp = .FALSE.
         EXIT
       END IF
     END DO
 
-  END FUNCTION is_sorted_real
+  END FUNCTION is_sorted_real_sp
 
   LOGICAL FUNCTION is_sorted_int(array)
     INTEGER, INTENT(IN) :: array(:)