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

Last change on this file since 3783 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

  • Property svn:keywords set to Id
File size: 3.4 KB
RevLine 
[1682]1!> @file random_gauss.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1036]9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[3655]17! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1321]22!
[2001]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: random_gauss.f90 3655 2019-01-07 16:51:22Z forkel $
[2716]27! Corrected "Former revisions" section
28!
29! 2696 2017-12-14 17:12:51Z kanani
30! Change in file header (GPL part)
[1321]31!
[2716]32! 2101 2017-01-05 16:42:31Z suehring
33!
[2001]34! 2000 2016-08-20 18:09:15Z knoop
35! Forced header and separation lines into 80 columns
36!
[1683]37! 1682 2015-10-07 23:56:08Z knoop
38! Code annotations made doxygen readable
39!
[1343]40! 1342 2014-03-26 17:04:47Z kanani
41! REAL constants defined as wp-kind
42!
[1321]43! 1320 2014-03-20 08:40:49Z raasch
[1320]44! ONLY-attribute added to USE-statements,
45! kind-parameters added to all INTEGER and REAL declaration statements,
46! kinds are defined in new module kinds,
47! old module precision_kind is removed,
48! revision history before 2012 removed,
49! comment fields (!:) to be used for variable explanations added to
50! all variable declaration statements
[1037]51!
52! 1036 2012-10-22 13:43:42Z raasch
53! code put under GPL (PALM 3.9)
54!
[3]55! RCS Log replace by Id keyword, revision history cleaned up
56!
[1]57! Revision 1.1  1998/03/25 20:09:47  raasch
58! Initial revision
59!
60!
61! Description:
62! ------------
[1682]63!> Generates a gaussian distributed random number (mean value 1, sigma = 1)
64!> This routine is taken from the "numerical recipies".
[1]65!------------------------------------------------------------------------------!
[1682]66 FUNCTION random_gauss( idum, upper_limit )
67 
[1]68
[1320]69    USE kinds
[1]70
[1320]71    USE random_function_mod,                                                   &
72        ONLY:  random_function
73
[1]74    IMPLICIT NONE
75
[1682]76    INTEGER(iwp) ::  idum          !<
77    INTEGER(iwp) ::  iset          !<
[1]78
[1682]79    REAL(wp)     ::  fac           !<
80    REAL(wp)     ::  gset          !<
81    REAL(wp)     ::  random_gauss  !<
82    REAL(wp)     ::  rsq           !<
83    REAL(wp)     ::  upper_limit   !<
84    REAL(wp)     ::  v1            !<
85    REAL(wp)     ::  v2            !<
[1320]86
[1]87    SAVE  iset, gset
88
89    DATA  iset /0/
90
91!
92!-- Random numbers are created as long as they do not fall below the given
93!-- upper limit
94    DO
95
96       IF ( iset == 0 )  THEN
[1342]97          rsq = 0.0_wp
98          DO  WHILE ( rsq >= 1.0_wp  .OR.  rsq == 0.0_wp )
99             v1  = 2.0_wp * random_function( idum ) - 1.0_wp
100             v2  = 2.0_wp * random_function( idum ) - 1.0_wp
[1]101             rsq = v1**2 + v2**2
102          ENDDO
[1342]103          fac          = SQRT( -2.0_wp * LOG( rsq ) / rsq )
[1]104          gset         = v1 * fac
[1342]105          random_gauss = v2 * fac + 1.0_wp
[1]106          iset         = 1
107       ELSE
[1342]108          random_gauss = gset + 1.0_wp
[1]109          iset         = 0
110       ENDIF
111
[1342]112       IF ( ABS( random_gauss - 1.0_wp ) <= upper_limit )  EXIT
[1]113
114    ENDDO
115
116 END FUNCTION random_gauss
Note: See TracBrowser for help on using the repository browser.