# [Dart-dev] [5679] DART/branches/development/random_seq/random_seq_mod.f90: add some bookkeeping in the random gauss routine.

nancy at ucar.edu nancy at ucar.edu
Mon Apr 9 15:06:24 MDT 2012

```Revision: 5679
Author:   nancy
Date:     2012-04-09 15:06:24 -0600 (Mon, 09 Apr 2012)
Log Message:
-----------
add some bookkeeping in the random gauss routine.  if you loop inside that routine
for more than 100 times without successfully generating a good gaussian value, then
most likely you've called the random number generator without initializing it.
call the error handler instead of sitting in an infinite loop.

Modified Paths:
--------------
DART/branches/development/random_seq/random_seq_mod.f90

-------------- next part --------------
Modified: DART/branches/development/random_seq/random_seq_mod.f90
===================================================================
--- DART/branches/development/random_seq/random_seq_mod.f90	2012-04-09 21:02:26 UTC (rev 5678)
+++ DART/branches/development/random_seq/random_seq_mod.f90	2012-04-09 21:06:24 UTC (rev 5679)
@@ -59,6 +59,7 @@
end type random_seq_type

logical, save :: module_initialized = .false.
+character(len=128) :: errstring

contains
@@ -314,6 +315,7 @@
real(r8) :: ran_gauss

real(digits12) :: x, y, r2, t
+integer :: lc

if ( .not. module_initialized ) call initialize_module

@@ -326,7 +328,9 @@
s%gset = .false.
else

-   10 continue
+   lc = 0
+10 continue
+   lc = lc + 1

! choose x,y in uniform square (-1,-1) to (+1,+1)
x = -1.0_digits12 + 2.0_digits12 * ran_unif(s)
@@ -334,7 +338,16 @@

! repeat if it is outside the unit circle or at origin
r2 = x*x + y*y
-   if (r2 >= 1.0_digits12 .or. r2 == 0.0_digits12) goto 10
+   if (r2 >= 1.0_digits12 .or. r2 == 0.0_digits12) then
+      if (lc > 100) then
+         write(errstring, *) 'x, y = ', x, y
+         call error_handler(E_ERR, 'ran_gauss', &
+                            'if both x and y are -1, random number generator probably not initialized', &
+                            source, revision, revdate, &
+                            text2 = errstring);
+      endif
+      goto 10
+   endif

t = sqrt(-2.0_digits12 * log(r2) / r2)
s%lastg   = real(x * t, r8)
```