ROMS
Loading...
Searching...
No Matches
gasdev.F File Reference
#include "cppdefs.h"
Include dependency graph for gasdev.F:

Go to the source code of this file.

Functions/Subroutines

subroutine gasdev_s (harvest)
 
subroutine gasdev_v (harvest)
 

Function/Subroutine Documentation

◆ gasdev_s()

subroutine gasdev_s ( real(r8), intent(out) harvest)

Definition at line 2 of file gasdev.F.

3!
4!git $Id$
5!================================================== Hernan G. Arango ===
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md !
9!=======================================================================
10! !
11! Return in harvest a normally distributed deviate with zero mean !
12! and unit variance, using RAN1 as the source of uniform deviates. !
13! !
14! Scalar version adapted from Numerical Recipes. !
15! !
16! Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, !
17! 1996: Numerical Recipes in Fortran 90, The Art of Parallel !
18! Scientific Computing, 2nd Edition, Cambridge Univ. Press. !
19! !
20!=======================================================================
21!
22 USE mod_kinds
23 USE nrutil, ONLY : ran1
24!
25! Imported variable declarations.
26!
27 real(r8), intent(out) :: harvest
28!
29! Local variable declarations.
30!
31 logical, save :: gaus_stored = .false.
32
33 real(r8) :: rsq, v1, v2
34 real(r8), save :: g
35!
36!-----------------------------------------------------------------------
37! Compute a normally distributed scalar deviate.
38!-----------------------------------------------------------------------
39!
40! We have an extra deviate handy, so return it, and unset the flag.
41!
42 IF (gaus_stored) THEN
43 harvest=g
44 gaus_stored=.false.
45!
46! We do not have an extra deviate handy, so pick two uniform numbers
47! in the square extending from -1 to +1 in each direction.
48!
49 ELSE
50 DO
51 CALL ran1 (v1)
52 CALL ran1 (v2)
53 v1=2.0_r8*v1-1.0_r8
54 v2=2.0_r8*v2-1.0_r8
55 rsq=v1*v1+v2*v2
56!
57! See if they are in the unit circle, and if they are not, try again.
58!
59 IF ((rsq.gt.0.0_r8).and.(rsq.lt.1.0_r8)) EXIT
60 END DO
61!
62! Now make the Box-Muller transformation to get two normal deviates.
63! Return one and save the other for next time.
64!
65 rsq=sqrt(-2.0_r8*log(rsq)/rsq)
66 harvest=v1*rsq
67 g=v2*rsq
68 gaus_stored=.true.
69 END IF
70
71 RETURN
Definition nrutil.F:1

◆ gasdev_v()

subroutine gasdev_v ( real(r8), dimension(:), intent(out) harvest)

Definition at line 74 of file gasdev.F.

75!
76!=======================================================================
77! !
78! Return in harvest a normally distributed deviate with zero mean !
79! and unit variance, using RAN1 as the source of uniform deviates. !
80! !
81! Vector version adapted from Numerical Recipes. !
82! !
83! Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, !
84! 1996: Numerical Recipes in Fortran 90, The Art of Parallel !
85! Scientific Computing, 2nd Edition, Cambridge Univ. Press. !
86! !
87!=======================================================================
88!
89 USE mod_kinds
90 USE nrutil, ONLY : array_copy
91 USE nrutil, ONLY : ran1
92!
93! Imported variable declarations.
94!
95 real(r8), dimension(:), intent(out) :: harvest
96!
97! Local variable declarations.
98!
99 logical, save :: gaus_stored = .false.
100
101 logical, dimension(SIZE(harvest)) :: mask
102
103 integer(i8b), save :: last_allocated = 0
104 integer(i8b) :: i, ii, m, mc, n, ng, nn
105
106 real(r8), dimension(SIZE(harvest)) :: rsq, v1, v2, v3
107 real(r8), allocatable, dimension(:), save :: g
108!
109!-----------------------------------------------------------------------
110! Compute a normally distributed vector deviate.
111!-----------------------------------------------------------------------
112!
113! We have an extra deviate handy, so return it, and unset the flag.
114!
115 n=SIZE(harvest)
116 IF (n.ne.last_allocated) THEN
117 IF (last_allocated.ne.0) DEALLOCATE (g)
118 ALLOCATE ( g(n) )
119 last_allocated=n
120 gaus_stored=.false.
121 END IF
122!
123! We do not have an extra deviate handy, so pick two uniform numbers
124! in the square extending from -1 to +1 in each direction.
125!
126 IF (gaus_stored) THEN
127 harvest=g
128 gaus_stored=.false.
129 ELSE
130 ng=1
131 DO
132 IF (ng.gt.n) EXIT
133 CALL ran1 (v1(ng:n))
134 CALL ran1 (v2(ng:n))
135 v1(ng:n)=2.0_r8*v1(ng:n)-1.0_r8
136 v2(ng:n)=2.0_r8*v2(ng:n)-1.0_r8
137!
138! See if they are in the unit circle, and if they are not, try again.
139! The original code was modified for portability with old F90 compilers
140! when using the PACK function (HGA/AMM).
141!
142 rsq(ng:n)=v1(ng:n)**2+v2(ng:n)**2
143 mask(ng:n)=((rsq(ng:n).gt.0.0_r8).and.(rsq(ng:n).lt.1.0_r8))
144 mc=count(mask(ng:n))
145
146 ii=0
147 DO i=ng,n
148 IF (mask(i)) THEN
149 ii=ii+1
150 v3(ii)=v1(i)
151 END IF
152 END DO
153 CALL array_copy (v3(1:mc), v1(ng:), nn, m)
154 ii=ng-1
155 DO i=ng,n
156 IF (mask(i)) THEN
157 ii=ii+1
158 v2(ii)=v2(i)
159 rsq(ii)=rsq(i)
160 END IF
161 END DO
162 ng=ng+nn
163 END DO
164!
165! Make the Box-Muller transformation to get two normal deviates.
166! Return one and save the other for next time.
167!
168 rsq=sqrt(-2.0_r8*log(rsq)/rsq)
169 harvest=v1*rsq
170 g=v2*rsq
171 gaus_stored=.true.
172 END IF
173
174 RETURN