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 NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC – NEMO

source: NEMO/branches/NERC/dev_release-3.4_NEMOTAM_consolidated/NEMOGCM/NEMO/OPATAM_SRC/ran_num.F90 @ 11774

Last change on this file since 11774 was 4598, checked in by pabouttier, 10 years ago

Add routines to allow restartability of PRNG, see Ticket #1284
  • Property svn:executable set to *
File size: 8.8 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    & mtrand_seeddump,   &
19    & mtrand_seedread
20 
21   IMPLICIT NONE
22
23   !! * Routine accessibility
24
25   PRIVATE
26
27   PUBLIC &
28      & gaustb, &
29      & gaustb_2d, &
30      & random_seeddump, &
31      & random_seedread
32
33CONTAINS
34
35   FUNCTION gaustb( pamp, pmean  )
36      !!----------------------------------------------------------------------
37      !!               ***  ROUTINE gaustb ***
38      !!         
39      !! ** Purpose : Returns a gaussian randum number with mean and std.
40      !!
41      !! ** Method  : Generate Gaussian random variables.
42      !!              The standard deviation and mean of the variables are
43      !!              specified by the variables pamp and pmean.
44      !!
45      !! ** Action  :
46      !!
47      !! History :
48      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
49      !!----------------------------------------------------------------------
50      !! * Function return
51      REAL(wp) :: &
52         & gaustb
53      !! * Arguments
54      REAL(wp), INTENT(IN) :: &
55         & pamp, &      ! Amplitude
56         & pmean        ! Mean value
57      !! * Local declarations
58     
59      gaustb = pamp * gausva( ) + pmean
60
61   END FUNCTION gaustb
62
63
64   FUNCTION gausva( )
65      !!----------------------------------------------------------------------
66      !!               ***  ROUTINE gausva ***
67      !!         
68      !! ** Purpose : Returns a normally distributed deviate with 0 mean
69      !!              and unit variance using the unifva(kdum) as the
70      !!              source of uniform deviates.
71      !!
72      !! ** Method  :
73      !!
74      !! ** Action  :
75      !!
76      !! History :
77      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
78      !!----------------------------------------------------------------------
79      !! * Function return
80      REAL(wp) :: &
81         & gausva
82      !! * Local declarations
83      REAL(wp), SAVE :: &
84         & gset
85      INTEGER, SAVE :: &
86         & niset = 0
87      REAL(wp) :: &
88         & zfac, &
89         & zrsq, &
90         & zv1,  &
91         & zv2
92
93      ! Begin main
94
95      IF ( niset == 0 ) THEN
96
97         zv1   = 2.0_wp * psrandom( ) - 1.0_wp
98         zv2   = 2.0_wp * psrandom( ) - 1.0_wp
99         zrsq  = zv1**2 + zv2**2
100         
101         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 
102
103            zv1   = 2.0_wp * psrandom( ) - 1.0_wp
104            zv2   = 2.0_wp * psrandom( ) - 1.0_wp
105            zrsq  = zv1**2 + zv2**2
106
107         END DO
108
109         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
110         gset   = zv1 * zfac
111         gausva = zv2 * zfac
112         niset  = 1
113
114      ELSE
115
116         gausva = gset
117         niset  = 0
118
119      ENDIF
120       
121   END FUNCTION gausva
122
123   FUNCTION gaustb_2d( ki, kj, pamp, pmean  )
124      !!----------------------------------------------------------------------
125      !!               ***  ROUTINE gaustb_2d ***
126      !!         
127      !! ** Purpose : Returns a Gaussian randum number with mean and std
128      !!              for a given horizontal grid point.
129      !!
130      !! ** Method  : Generate Gaussian random variables.
131      !!              The standard deviation and mean of the variables are
132      !!              specified by the variables pamp and pmean.
133      !!
134      !! ** Action  :
135      !!
136      !! History :
137      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
138      !!----------------------------------------------------------------------
139      !! * Function return
140      REAL(wp) :: &
141         & gaustb_2d
142      !! * Arguments
143      INTEGER, INTENT(IN) :: &
144         & ki, &        ! Indices in seed array
145         & kj
146      REAL(wp), INTENT(IN) :: &
147         & pamp, &      ! Amplitude
148         & pmean        ! Mean value
149      !! * Local declarations
150     
151      gaustb_2d = pamp * gausva_2d( ki, kj ) + pmean
152
153   END FUNCTION gaustb_2d
154
155   FUNCTION gausva_2d( ki, kj )
156      !!----------------------------------------------------------------------
157      !!               ***  ROUTINE gausva_2d ***
158      !!         
159      !! ** Purpose : Returns a normally distributed deviate with 0 mean
160      !!              and unit variance.
161      !!
162      !! ** Method  :
163      !!
164      !! ** Action  :
165      !!
166      !! History :
167      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
168      !!----------------------------------------------------------------------
169      !! * Function return
170      REAL(wp) :: &
171         & gausva_2d
172      !! * Arguments
173      INTEGER, INTENT(IN) :: &
174         & ki, &         ! Indices in seed array
175         & kj
176      !! * Local declarations
177      REAL(wp), DIMENSION(jpi,jpj) :: &
178         & gset
179      INTEGER, DIMENSION(jpi,jpj) :: &
180         & niset
181      LOGICAL, SAVE :: &
182         & llinit = .FALSE.
183      REAL(wp) :: &
184         & zfac, &
185         & zrsq, &
186         & zv1,  &
187         & zv2
188
189      ! Initialization
190
191      IF ( .NOT. llinit ) THEN
192
193         niset(:,:) = 0
194         llinit     = .TRUE.
195
196      ENDIF
197
198      ! Begin main
199
200      IF ( niset(ki,kj) == 0 ) THEN
201
202         zv1   = 2.0_wp * psrandom( ) - 1.0_wp
203         zv2   = 2.0_wp * psrandom( ) - 1.0_wp
204         zrsq  = zv1**2 + zv2**2
205
206         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 
207
208            zv1   = 2.0_wp * psrandom( ) - 1.0_wp
209            zv2   = 2.0_wp * psrandom( ) - 1.0_wp
210            zrsq  = zv1**2 + zv2**2
211
212         END DO
213
214         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
215         gset(ki,kj) = zv1 * zfac
216         gausva_2d   = zv2 * zfac
217         niset(ki,kj) = 1
218
219      ELSE
220
221         gausva_2d = gset(ki,kj)
222         niset(ki,kj)  = 0
223
224      ENDIF
225       
226   END FUNCTION gausva_2d
227
228   FUNCTION psrandom()
229      !!----------------------------------------------------------------------
230      !!               ***  ROUTINE psrandom ***
231      !!
232      !! ** Purpose : Pseudo-Random number generator.
233      !!
234      !! ** Method  : Returns a pseudo-random number from a uniform distribution
235      !!              between 0 and 1
236      !!              Call with kdum a negative integer to initialize.
237      !!              Thereafter, do not alter kdum between successive deviates
238      !!              in sequence.
239      !!
240      !! ** Action  :
241      !!
242      !! History :
243      !!        !  10-02  (F. Vigilant)  Original code
244      !!        ! 2011-08 (D. Lea)  initialize with kdum1
245      !!----------------------------------------------------------------------
246      !! * Function return
247      REAL(wp) ::  &
248         & psrandom
249      !! * Arguments
250      INTEGER  :: &
251         & kdum          ! Seed
252      LOGICAL, SAVE :: &
253         & llinit = .TRUE.    ! Initialize on first call
254
255      ! Initialization
256      IF ( llinit ) THEN
257         kdum = 456953 + nproc * 596035
258         CALL init_mtrand(kdum)
259         llinit     = .FALSE.    ! On subsequent calls do not reinitialize
260      ENDIF
261
262      psrandom = mtrand_real1()
263
264   END FUNCTION psrandom
265
266   
267   SUBROUTINE random_seeddump(cdfile)
268      !!----------------------------------------------------------------------
269      !!               ***  ROUTINE psrandom ***
270      !!
271      !! ** Purpose : Pseudo-Random number generator.
272      !!
273      !! ** Method  : Dump the seed
274      !!
275      !! ** Action  :
276      !!
277      !! History :
278      !!        ! 2013-03 (K. Mogensen) Original code
279      !!----------------------------------------------------------------------
280      !! * Arguments
281      CHARACTER(len=128) :: cdfile
282      CALL mtrand_seeddump(cdfile)
283
284   END SUBROUTINE random_seeddump
285
286   SUBROUTINE random_seedread(cdfile)
287      !!----------------------------------------------------------------------
288      !!               ***  ROUTINE psrandom ***
289      !!
290      !! ** Purpose : Pseudo-Random number generator.
291      !!
292      !! ** Method  : Read the seed
293      !!
294      !! ** Action  :
295      !!
296      !! History :
297      !!        ! 2013-03 (K. Mogensen) Original code
298      !!----------------------------------------------------------------------
299      !! * Arguments
300      CHARACTER(len=128) :: cdfile
301      CALL mtrand_seedread(cdfile)
302
303   END SUBROUTINE random_seedread
304
305END MODULE ran_num
Note: See TracBrowser for help on using the repository browser.