source: palm/trunk/SOURCE/random_gauss.f90 @ 1

Last change on this file since 1 was 1, checked in by raasch, 14 years ago

Initial repository layout and content

File size: 1.7 KB
Line 
1 FUNCTION random_gauss( idum, upper_limit )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: random_gauss.f90,v $
11! Revision 1.4  2006/08/04 15:01:48  raasch
12! Range of random number is limited by an upper limit (new second parameter)
13!
14! Revision 1.3  2003/10/29 09:08:04  raasch
15! Module random_function_mod is used
16!
17! Revision 1.2  2001/01/22 08:03:31  raasch
18! Comments translated into English
19!
20! Revision 1.1  1998/03/25 20:09:47  raasch
21! Initial revision
22!
23!
24! Description:
25! ------------
26! Generates a gaussian distributed random number (mean value 1, sigma = 1)
27! This routine is taken from the "numerical recipies".
28!------------------------------------------------------------------------------!
29
30    USE random_function_mod
31
32    IMPLICIT NONE
33
34    INTEGER :: idum, iset
35    REAL    :: fac, gset, random_gauss, rsq, upper_limit, v1, v2
36
37    SAVE  iset, gset
38
39    DATA  iset /0/
40
41!
42!-- Random numbers are created as long as they do not fall below the given
43!-- upper limit
44    DO
45
46       IF ( iset == 0 )  THEN
47          rsq = 0.0
48          DO  WHILE ( rsq >= 1.0  .OR.  rsq == 0.0 )
49             v1  = 2.0 * random_function( idum ) - 1.0
50             v2  = 2.0 * random_function( idum ) - 1.0
51             rsq = v1**2 + v2**2
52          ENDDO
53          fac          = SQRT( -2.0 * LOG( rsq ) / rsq )
54          gset         = v1 * fac
55          random_gauss = v2 * fac + 1.0
56          iset         = 1
57       ELSE
58          random_gauss = gset + 1.0
59          iset         = 0
60       ENDIF
61
62       IF ( ABS( random_gauss - 1.0 ) <= upper_limit )  EXIT
63
64    ENDDO
65
66 END FUNCTION random_gauss
Note: See TracBrowser for help on using the repository browser.