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 @ 1885

Last change on this file since 1885 was 1885, checked in by rblod, 14 years ago

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