diff --git a/AUTHORS.txt b/AUTHORS.txt
index 38bd3331b9e9cac3fa2a3e48c8de094e8811d565..338c3351b580ecd411fd8c749f9962b1126efbca 100644
--- a/AUTHORS.txt
+++ b/AUTHORS.txt
@@ -10,6 +10,7 @@
# ---------------------------------------------------------------
# A more complete author list is provided in the ICON repository.
+Daniel Hupp
Daniel Reinert
Daniel Rieger
Florian Prill
diff --git a/README.md b/README.md
index 09e56b86f28e88c7509d132ca232180d880a6f3d..7baf4ef60422f5f191e233c4201f14225f9ccbf3 100644
--- a/README.md
+++ b/README.md
@@ -50,6 +50,7 @@ The `libfortran-support` library includes some general Fortran supporting module
- `mo_io_units`: io unit definitions
- `mo_namelist`: open/close namelist files
- `mo_octree`: octree data structure and operations
+- `mo_random_number_generators`: generators for pseudo random numbers (should be moved to math-support once available)
- `mo_simple_dump`: array value dumping
- `mo_util_backtrace`: function backtrace
- `mo_util_file`: file operations
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7db2325be8cb497271dee192f7aaabade11733cc..a1427110ef9a262fc78b3f9c7f516f017b5351c1 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -32,6 +32,7 @@ add_library(fortran-support
mo_io_units.F90
mo_namelist.F90
mo_octree.F90
+ mo_random_number_generators.F90
mo_simple_dump.F90
mo_util_backtrace.F90
mo_util_file.F90
diff --git a/src/mo_random_number_generators.F90 b/src/mo_random_number_generators.F90
new file mode 100644
index 0000000000000000000000000000000000000000..e77e04acbde284a696883202b50a4af8d9df76e0
--- /dev/null
+++ b/src/mo_random_number_generators.F90
@@ -0,0 +1,152 @@
+! ICON
+!
+! ---------------------------------------------------------------
+! Copyright (C) 2004-2024, 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
+! ---------------------------------------------------------------
+
+MODULE mo_random_number_generators
+
+ USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+ USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: int32, int64, real64
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC :: random_normal_values, clip, initialize_random_number_generator, generate_uniform_random_number
+
+ INTERFACE generate_uniform_random_number
+ MODULE PROCEDURE generate_uniform_random_number_wp
+ END INTERFACE generate_uniform_random_number
+
+ INTERFACE random_normal_values
+ MODULE PROCEDURE random_normal_values_wp
+ END INTERFACE random_normal_values
+
+ INTERFACE clip
+ MODULE PROCEDURE clip_wp
+ END INTERFACE clip
+
+ ! integer parameters for using revised minstd parameters as in
+ ! https://en.wikipedia.org/wiki/Lehmer_random_number_generator#Parameters_in_common_use
+ INTEGER(KIND=int64), PARAMETER :: minstd_A0 = 16807
+ INTEGER(KIND=int64), PARAMETER :: minstd_A = 48271
+
+ INTEGER(KIND=int64), PARAMETER :: minstd_M = 2147483647
+ REAL(KIND=real64), PARAMETER :: minstd_Scale = 1.0_real64/REAL(minstd_M, real64)
+
+CONTAINS
+
+ ! initializes Lehmer random number generator
+ PURE FUNCTION initialize_random_number_generator(iseed, jseed) RESULT(minstd_state)
+
+ INTEGER(KIND=int32), INTENT(IN) :: iseed
+ INTEGER(KIND=int32), INTENT(IN) :: jseed
+
+ INTEGER(KIND=int32) :: minstd_state
+
+ ! doing kind of an xorshift32 see https://en.wikipedia.org/wiki/Xorshift
+ minstd_state = iseed
+ minstd_state = IEOR(minstd_state, jseed*2**13)
+ minstd_state = IEOR(minstd_state, jseed/2**17)
+ minstd_state = IEOR(minstd_state, jseed*2**5)
+
+ minstd_state = MOD(minstd_A0*minstd_state, minstd_M)
+
+ minstd_state = MOD(minstd_A*minstd_state, minstd_M)
+
+ END FUNCTION initialize_random_number_generator
+
+ ! generated Lehmer random number generator https://en.wikipedia.org/wiki/Lehmer_random_number_generator
+ FUNCTION generate_uniform_random_number_wp(minstd_state) RESULT(random_number)
+
+ INTEGER(KIND=int32), INTENT(INOUT) :: minstd_state
+ REAL(KIND=real64) :: random_number
+
+ minstd_state = MOD(minstd_A*minstd_state, minstd_M)
+
+ random_number = ABS(minstd_Scale*minstd_state)
+
+ END FUNCTION generate_uniform_random_number_wp
+
+ SUBROUTINE clip_wp(values_range, values)
+
+ ! Subroutine arguments (in/out/inout)
+ REAL(wp), INTENT(IN) :: values_range ! allowed range of random numbers
+ REAL(wp), DIMENSION(:), INTENT(INOUT) :: values ! array of value to be filled with random numbers
+
+ ! Limit random number range
+ WHERE (values(:) > values_range)
+ values(:) = values_range
+ ELSEWHERE(values(:) < -values_range)
+ values(:) = -values_range
+ END WHERE
+
+ END SUBROUTINE clip_wp
+
+ SUBROUTINE random_normal_values_wp(seed, values_range, values, mean, stdv)
+
+ ! Subroutine arguments (in/out/inout)
+ INTEGER, INTENT(IN) :: seed ! seed to be used
+ REAL(wp), INTENT(IN) :: values_range ! allowed range of random numbers
+ REAL(wp), DIMENSION(:), INTENT(INOUT) :: values ! array of value to be filled with random numbers
+ REAL(wp), INTENT(IN), OPTIONAL :: mean
+ REAL(wp), INTENT(IN), OPTIONAL :: stdv
+
+ REAL(wp), PARAMETER :: pi = 3.141592653589793238_wp ! replace with math constants from math-support once available
+
+ ! Local variables
+ INTEGER(KIND=int32) :: i, rng_state, n
+ REAL(KIND=wp) :: u1, u2, z0, z1, magnitude, theta, mu, sigma
+
+ IF (PRESENT(mean)) THEN
+ mu = mean
+ ELSE
+ mu = 0.0_wp
+ END IF
+
+ IF (PRESENT(stdv)) THEN
+ sigma = stdv
+ ELSE
+ sigma = 1.0_wp
+ END IF
+
+ n = SIZE(values)
+
+ ! -----------------------------------------------------------------------
+ ! Initialize a random number generator, by using the MINSTD linear
+ ! congruential generator (LCG). "nmaxstreams" indicates that random
+ ! numbers will be requested in blocks of this length. The generator
+ ! is seeded with "seed".
+ !-------------------------------------------------------------------------
+
+ DO i = 1, n/2 + 1
+ rng_state = initialize_random_number_generator(seed, i)
+
+ u1 = generate_uniform_random_number(rng_state)
+ DO WHILE (u1 <= EPSILON(0.0_wp))
+ u1 = generate_uniform_random_number(rng_state)
+ END DO
+
+ u2 = generate_uniform_random_number(rng_state)
+
+ ! using the Box-Muller transform https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
+ magnitude = sigma*SQRT(-2.0_wp*LOG(u1));
+ theta = 2.0_wp*pi*u2
+
+ z0 = magnitude*COS(theta) + mu;
+ z1 = magnitude*SIN(theta) + mu;
+ values(2*i - 1) = z0
+ IF (2*i <= n) values(2*i) = z1
+
+ END DO
+
+ CALL clip(values_range, values)
+
+ END SUBROUTINE random_normal_values_wp
+
+END MODULE mo_random_number_generators
diff --git a/test/fortran/helpers.f90 b/test/fortran/helpers.f90
index fc2c604e7364f1ba59a5e9f5102eb29b1dd53936..c30ead76b568b328f16ed64f0ce01be35b234adc 100644
--- a/test/fortran/helpers.f90
+++ b/test/fortran/helpers.f90
@@ -13,6 +13,7 @@ MODULE helpers
USE ISO_C_BINDING, ONLY: c_double
USE mo_io_units, ONLY: find_next_free_unit
USE mo_exception, ONLY: message
+ USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
CONTAINS
@@ -52,4 +53,57 @@ CONTAINS
OPEN (unit=nerr, file=file, status='old', action='read')
END SUBROUTINE open_logfile
+
+ FUNCTION calculate_mean_wp(array) RESULT(mean)
+ REAL(wp), INTENT(IN) :: array(:)
+ REAL(wp) :: mean
+ INTEGER :: i, n
+
+ mean = 0.0_wp
+ variance = 0.0_wp
+
+ n = SIZE(array)
+
+ DO i = 1, SIZE(array)
+ mean = mean + array(i)
+ END DO
+ mean = mean/n
+
+ END FUNCTION
+
+ FUNCTION calculate_variance_wp(array) RESULT(variance)
+ REAL(wp), INTENT(IN) :: array(:)
+ REAL(wp) :: mean, variance
+ INTEGER :: i, n
+
+ mean = 0.0_wp
+ variance = 0.0_wp
+
+ n = SIZE(array)
+ mean = calculate_mean_wp(array)
+
+ DO i = 1, SIZE(array)
+ variance = variance + (array(i) - mean)**2
+ END DO
+ variance = variance/n
+
+ END FUNCTION
+
+ SUBROUTINE assert_statistics(test_name, array, mean, variance, tol, max, min)
+ USE FORTUTF
+ CHARACTER(LEN=*), INTENT(IN) :: test_name
+ REAL(wp), INTENT(IN) :: array(:)
+ REAL(wp), INTENT(IN) :: mean, variance, tol, max, min
+
+ CALL TAG_TEST(test_name//"__mean")
+ CALL ASSERT_ALMOST_EQUAL(calculate_mean_wp(array), mean, tol)
+ CALL TAG_TEST(test_name//"__variance")
+ CALL ASSERT_ALMOST_EQUAL(calculate_variance_wp(array), variance, tol)
+ CALL TAG_TEST(test_name//"__max")
+ CALL ASSERT_LESS_THAN_EQUAL(MAXVAL(array), max)
+ CALL TAG_TEST(test_name//"__min")
+ CALL ASSERT_GREATER_THAN_EQUAL(MINVAL(array), min)
+
+ END SUBROUTINE
+
END MODULE helpers
diff --git a/test/fortran/test_random_number_generators.f90 b/test/fortran/test_random_number_generators.f90
new file mode 100644
index 0000000000000000000000000000000000000000..78f4932561920490df80d85f4390f6ad34c4e636
--- /dev/null
+++ b/test/fortran/test_random_number_generators.f90
@@ -0,0 +1,101 @@
+! ICON
+!
+! ---------------------------------------------------------------
+! Copyright (C) 2004-2024, 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
+! ---------------------------------------------------------------
+
+MODULE TEST_mo_random_number_generators
+ USE FORTUTF
+ USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
+ USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: int32
+ USE helpers, ONLY: assert_statistics, calculate_mean_wp, calculate_variance_wp
+
+CONTAINS
+ SUBROUTINE TEST_CALCULATE_STATISTICS_WP
+ REAL(wp) :: array(8) = (/2.0_wp, 4.0_wp, 4.0_wp, 4.0_wp, 5.0_wp, 5.0_wp, 7.0_wp, 9.0_wp/)
+ CALL TAG_TEST("TEST_calculate_mean_wp")
+ CALL ASSERT_ALMOST_EQUAL(calculate_mean_wp(array), 5.0_wp)
+ CALL TAG_TEST("TEST_calculate_variance_wp")
+ CALL ASSERT_ALMOST_EQUAL(calculate_variance_wp(array), 4.0_wp)
+ END SUBROUTINE
+
+ SUBROUTINE TEST_CLIP_WP
+ USE mo_random_number_generators
+ REAL(wp) :: array(6) = (/4.1_wp, -4.1_wp, 2.0_wp, -2.0_wp, 1.0_wp, -1.0_wp/)
+ CALL TAG_TEST("TEST_clip_wp")
+ CALL clip(2.0_wp, array)
+ CALL ASSERT_EQUAL(array, (/2.0_wp, -2.0_wp, 2.0_wp, -2.0_wp, 1.0_wp, -1.0_wp/))
+ END SUBROUTINE
+
+ SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_BASIC
+ USE mo_random_number_generators
+ INTEGER, PARAMETER :: n = 1e5
+ INTEGER :: i
+ INTEGER(int32) :: rng_state
+ REAL(wp) :: array(n)
+
+ rng_state = initialize_random_number_generator(11, 17)
+ DO i = 1, n
+ array(i) = generate_uniform_random_number(rng_state)
+ END DO
+
+ CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_BASIC", &
+ array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp)
+
+ END SUBROUTINE
+
+ SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_ISEED
+ USE mo_random_number_generators
+ INTEGER, PARAMETER :: n = 1e5
+ INTEGER :: i
+ INTEGER(int32) :: rng_state
+ REAL(wp) :: array(n)
+
+ DO i = 1, n
+ rng_state = initialize_random_number_generator(i, 11)
+ array(i) = generate_uniform_random_number(rng_state)
+ END DO
+
+ CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_ISEED", &
+ array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp)
+
+ END SUBROUTINE
+
+ SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_JSEED
+ USE mo_random_number_generators
+ INTEGER, PARAMETER :: n = 1e5
+ INTEGER :: i
+ INTEGER(int32) :: rng_state
+ REAL(wp) :: array(n)
+
+ DO i = 1, n
+ rng_state = initialize_random_number_generator(11, i)
+ array(i) = generate_uniform_random_number(rng_state)
+ END DO
+
+ CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_JSEED", &
+ array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp)
+
+ END SUBROUTINE
+
+ SUBROUTINE TEST_RANDOM_NORMAL_VALUES_WP
+ USE mo_random_number_generators
+ INTEGER, PARAMETER :: n = 1e5 + 1
+ REAL(wp) :: array(n)
+ REAL(wp) :: mean = -5.0_wp
+ REAL(wp) :: variance = 7.1_wp
+ REAL(wp) :: value_range = 1.0e32_wp
+
+ CALL random_normal_values(seed=11, values_range=value_range, values=array, mean=mean, stdv=SQRT(variance))
+
+ CALL assert_statistics("TEST_RANDOM_NORMAL_VALUES_WP", &
+ array, mean=mean, variance=variance, tol=1.0e-2_wp, max=value_range, min=-value_range)
+
+ END SUBROUTINE
+
+END MODULE TEST_mo_random_number_generators