New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
ran_num.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/ran_num.F90 @ 2587

Last change on this file since 2587 was 2587, checked in by vidard, 13 years ago

refer to ticket #798
File size: 7.9 KB
Line 
1MODULE ran_num
2   !!======================================================================
3   !!                       *** MODULE ran_num ***
4   !! NEMOVAR: Random number routines
5   !!======================================================================
6   !!----------------------------------------------------------------------
7   !! gaustb    : returns a gaussian randum number with mean and std.
8   !! gaustb_2d : returns a gaussian randum number with mean and std for a
9   !!             given horizontal grid point
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE par_kind       ! Kind variables
13   USE dom_oce        ! Domain variables
14   USE in_out_manager ! I/O stuff
15   USE mt19937ar, ONLY : &
16    & init_mtrand,       &
17    & mtrand_real1
18
19   IMPLICIT NONE
20
21   !! * Routine accessibility
22
23   PRIVATE
24
25   PUBLIC &
26      & gaustb, &
27      & gaustb_2d
28
29CONTAINS
30
31   FUNCTION gaustb( kseed, pamp, pmean  )
32      !!----------------------------------------------------------------------
33      !!               ***  ROUTINE gaustb ***
34      !!         
35      !! ** Purpose : Returns a gaussian randum number with mean and std.
36      !!
37      !! ** Method  : Generate Gaussian random variables.
38      !!              The standard deviation and mean of the variables are
39      !!              specified by the variables pamp and pmean.
40      !!
41      !! ** Action  :
42      !!
43      !! History :
44      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
45      !!----------------------------------------------------------------------
46      !! * Function return
47      REAL(wp) :: &
48         & gaustb
49      !! * Arguments
50      INTEGER, INTENT(INOUT) :: &
51         & kseed        ! Seed
52      REAL(wp), INTENT(IN) :: &
53         & pamp, &      ! Amplitude
54         & pmean        ! Mean value
55      !! * Local declarations
56     
57      gaustb = pamp * gausva( kseed ) + pmean
58
59   END FUNCTION gaustb
60
61
62   FUNCTION gausva( kdum )
63      !!----------------------------------------------------------------------
64      !!               ***  ROUTINE gausva ***
65      !!         
66      !! ** Purpose : Returns a normally distributed deviate with 0 mean
67      !!              and unit variance using the unifva(kdum) as the
68      !!              source of uniform deviates.
69      !!
70      !! ** Method  :
71      !!
72      !! ** Action  :
73      !!
74      !! History :
75      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
76      !!----------------------------------------------------------------------
77      !! * Function return
78      REAL(wp) :: &
79         & gausva
80      !! * Arguments
81      INTEGER, INTENT(INOUT) :: &
82         & kdum          ! Seed
83      !! * Local declarations
84      REAL(wp), SAVE :: &
85         & gset
86      INTEGER, SAVE :: &
87         & niset = 0
88      REAL(wp) :: &
89         & zfac, &
90         & zrsq, &
91         & zv1,  &
92         & zv2
93
94      ! Begin main
95
96      IF ( niset == 0 ) THEN
97
98         zv1   = 2.0_wp * psrandom( kdum ) - 1.0_wp
99         zv2   = 2.0_wp * psrandom( kdum ) - 1.0_wp
100         zrsq  = zv1**2 + zv2**2
101         
102         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 
103
104            zv1   = 2.0_wp * psrandom( kdum ) - 1.0_wp
105            zv2   = 2.0_wp * psrandom( kdum ) - 1.0_wp
106            zrsq  = zv1**2 + zv2**2
107
108         END DO
109
110         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
111         gset   = zv1 * zfac
112         gausva = zv2 * zfac
113         niset  = 1
114
115      ELSE
116
117         gausva = gset
118         niset  = 0
119
120      ENDIF
121       
122   END FUNCTION gausva
123
124   FUNCTION gaustb_2d( ki, kj, kseed, pamp, pmean  )
125      !!----------------------------------------------------------------------
126      !!               ***  ROUTINE gaustb_2d ***
127      !!         
128      !! ** Purpose : Returns a Gaussian randum number with mean and std
129      !!              for a given horizontal grid point.
130      !!
131      !! ** Method  : Generate Gaussian random variables.
132      !!              The standard deviation and mean of the variables are
133      !!              specified by the variables pamp and pmean.
134      !!
135      !! ** Action  :
136      !!
137      !! History :
138      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
139      !!----------------------------------------------------------------------
140      !! * Function return
141      REAL(wp) :: &
142         & gaustb_2d
143      !! * Arguments
144      INTEGER, INTENT(IN) :: &
145         & ki, &        ! Indices in seed array
146         & kj
147      INTEGER, INTENT(INOUT), DIMENSION(jpi,jpj) :: &
148         & kseed        ! Seed
149      REAL(wp), INTENT(IN) :: &
150         & pamp, &      ! Amplitude
151         & pmean        ! Mean value
152      !! * Local declarations
153     
154      gaustb_2d = pamp * gausva_2d( ki, kj, kseed ) + pmean
155
156   END FUNCTION gaustb_2d
157
158   FUNCTION gausva_2d( ki, kj, kdum )
159      !!----------------------------------------------------------------------
160      !!               ***  ROUTINE gausva_2d ***
161      !!         
162      !! ** Purpose : Returns a normally distributed deviate with 0 mean
163      !!              and unit variance.
164      !!
165      !! ** Method  :
166      !!
167      !! ** Action  :
168      !!
169      !! History :
170      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
171      !!----------------------------------------------------------------------
172      !! * Function return
173      REAL(wp) :: &
174         & gausva_2d
175      !! * Arguments
176      INTEGER, INTENT(IN) :: &
177         & ki, &         ! Indices in seed array
178         & kj
179      INTEGER, INTENT(INOUT), DIMENSION(jpi,jpj) :: &
180         & kdum          ! Seed
181      !! * Local declarations
182      REAL(wp), SAVE, DIMENSION(jpi,jpj) :: &
183         & gset
184      INTEGER, SAVE, DIMENSION(jpi,jpj) :: &
185         & niset
186      LOGICAL, SAVE :: &
187         & llinit = .FALSE.
188      REAL(wp) :: &
189         & zfac, &
190         & zrsq, &
191         & zv1,  &
192         & zv2
193
194      ! Initialization
195
196      IF ( .NOT. llinit ) THEN
197
198         niset(:,:) = 0
199         llinit     = .TRUE.
200
201      ENDIF
202
203      ! Begin main
204
205      IF ( niset(ki,kj) == 0 ) THEN
206
207         zv1   = 2.0_wp * psrandom( kdum(ki,kj) ) - 1.0_wp
208         zv2   = 2.0_wp * psrandom( kdum(ki,kj) ) - 1.0_wp
209         zrsq  = zv1**2 + zv2**2
210
211         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 
212
213            zv1   = 2.0_wp * psrandom( kdum(ki,kj) ) - 1.0_wp
214            zv2   = 2.0_wp * psrandom( kdum(ki,kj) ) - 1.0_wp
215            zrsq  = zv1**2 + zv2**2
216
217         END DO
218
219         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
220         gset(ki,kj) = zv1 * zfac
221         gausva_2d   = zv2 * zfac
222         niset(ki,kj) = 1
223
224      ELSE
225
226         gausva_2d = gset(ki,kj)
227         niset(ki,kj)  = 0
228
229      ENDIF
230       
231   END FUNCTION gausva_2d
232
233   FUNCTION psrandom( kdum )
234      !!----------------------------------------------------------------------
235      !!               ***  ROUTINE psrandom ***
236      !!         
237      !! ** Purpose : Pseudo-Random number generator.
238      !!
239      !! ** Method  : Returns a pseudo-random number from a uniform distribution
240      !!              between 0 and 1
241      !!              Call with kdum a negative integer to initialize.
242      !!              Thereafter, do not alter kdum between successive deviates
243      !!              in sequence.
244      !!
245      !! ** Action  :
246      !!
247      !! History :
248      !!        !  10-02  (F. Vigilant)  Original code
249      !!----------------------------------------------------------------------
250      !! * Function return
251      REAL(wp) ::  &
252         & psrandom
253      !! * Arguments
254      INTEGER, INTENT(INOUT) :: &
255         & kdum          ! Seed
256      LOGICAL, SAVE :: &
257         & llinit = .FALSE.
258      INTEGER :: &
259   & kdum1, &
260   & kdum2
261
262      ! Initialization
263      IF ( .NOT. llinit ) THEN
264    kdum2 = 596035
265    kdum1 = kdum + nproc * kdum2
266         CALL init_mtrand(kdum)
267         llinit     = .TRUE.
268      ENDIF
269
270      psrandom = mtrand_real1()
271     
272   END FUNCTION psrandom
273
274
275END MODULE ran_num
Note: See TracBrowser for help on using the repository browser.