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/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/ran_num.F90 @ 3611

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007

  • Property svn:executable set to *
File size: 7.5 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( 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      REAL(wp), INTENT(IN) :: &
51         & pamp, &      ! Amplitude
52         & pmean        ! Mean value
53      !! * Local declarations
54     
55      gaustb = pamp * gausva( ) + pmean
56
57   END FUNCTION gaustb
58
59
60   FUNCTION gausva( )
61      !!----------------------------------------------------------------------
62      !!               ***  ROUTINE gausva ***
63      !!         
64      !! ** Purpose : Returns a normally distributed deviate with 0 mean
65      !!              and unit variance using the unifva(kdum) as the
66      !!              source of uniform deviates.
67      !!
68      !! ** Method  :
69      !!
70      !! ** Action  :
71      !!
72      !! History :
73      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
74      !!----------------------------------------------------------------------
75      !! * Function return
76      REAL(wp) :: &
77         & gausva
78      !! * Local declarations
79      REAL(wp), SAVE :: &
80         & gset
81      INTEGER, SAVE :: &
82         & niset = 0
83      REAL(wp) :: &
84         & zfac, &
85         & zrsq, &
86         & zv1,  &
87         & zv2
88
89      ! Begin main
90
91      IF ( niset == 0 ) THEN
92
93         zv1   = 2.0_wp * psrandom( ) - 1.0_wp
94         zv2   = 2.0_wp * psrandom( ) - 1.0_wp
95         zrsq  = zv1**2 + zv2**2
96         
97         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 
98
99            zv1   = 2.0_wp * psrandom( ) - 1.0_wp
100            zv2   = 2.0_wp * psrandom( ) - 1.0_wp
101            zrsq  = zv1**2 + zv2**2
102
103         END DO
104
105         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
106         gset   = zv1 * zfac
107         gausva = zv2 * zfac
108         niset  = 1
109
110      ELSE
111
112         gausva = gset
113         niset  = 0
114
115      ENDIF
116       
117   END FUNCTION gausva
118
119   FUNCTION gaustb_2d( ki, kj, pamp, pmean  )
120      !!----------------------------------------------------------------------
121      !!               ***  ROUTINE gaustb_2d ***
122      !!         
123      !! ** Purpose : Returns a Gaussian randum number with mean and std
124      !!              for a given horizontal grid point.
125      !!
126      !! ** Method  : Generate Gaussian random variables.
127      !!              The standard deviation and mean of the variables are
128      !!              specified by the variables pamp and pmean.
129      !!
130      !! ** Action  :
131      !!
132      !! History :
133      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
134      !!----------------------------------------------------------------------
135      !! * Function return
136      REAL(wp) :: &
137         & gaustb_2d
138      !! * Arguments
139      INTEGER, INTENT(IN) :: &
140         & ki, &        ! Indices in seed array
141         & kj
142      REAL(wp), INTENT(IN) :: &
143         & pamp, &      ! Amplitude
144         & pmean        ! Mean value
145      !! * Local declarations
146     
147      gaustb_2d = pamp * gausva_2d( ki, kj ) + pmean
148
149   END FUNCTION gaustb_2d
150
151   FUNCTION gausva_2d( ki, kj )
152      !!----------------------------------------------------------------------
153      !!               ***  ROUTINE gausva_2d ***
154      !!         
155      !! ** Purpose : Returns a normally distributed deviate with 0 mean
156      !!              and unit variance.
157      !!
158      !! ** Method  :
159      !!
160      !! ** Action  :
161      !!
162      !! History :
163      !!        !  07-07  (K. Mogensen)  Original code based on gaustb.F
164      !!----------------------------------------------------------------------
165      !! * Function return
166      REAL(wp) :: &
167         & gausva_2d
168      !! * Arguments
169      INTEGER, INTENT(IN) :: &
170         & ki, &         ! Indices in seed array
171         & kj
172      !! * Local declarations
173      REAL(wp), DIMENSION(jpi,jpj) :: &
174         & gset
175      INTEGER, DIMENSION(jpi,jpj) :: &
176         & niset
177      LOGICAL, SAVE :: &
178         & llinit = .FALSE.
179      REAL(wp) :: &
180         & zfac, &
181         & zrsq, &
182         & zv1,  &
183         & zv2
184
185      ! Initialization
186
187      IF ( .NOT. llinit ) THEN
188
189         niset(:,:) = 0
190         llinit     = .TRUE.
191
192      ENDIF
193
194      ! Begin main
195
196      IF ( niset(ki,kj) == 0 ) THEN
197
198         zv1   = 2.0_wp * psrandom( ) - 1.0_wp
199         zv2   = 2.0_wp * psrandom( ) - 1.0_wp
200         zrsq  = zv1**2 + zv2**2
201
202         DO WHILE ( ( zrsq >= 1.0_wp ) .OR. ( zrsq == 0.0_wp ) ) 
203
204            zv1   = 2.0_wp * psrandom( ) - 1.0_wp
205            zv2   = 2.0_wp * psrandom( ) - 1.0_wp
206            zrsq  = zv1**2 + zv2**2
207
208         END DO
209
210         zfac   = SQRT( -2.0_wp * LOG( zrsq ) / zrsq )
211         gset(ki,kj) = zv1 * zfac
212         gausva_2d   = zv2 * zfac
213         niset(ki,kj) = 1
214
215      ELSE
216
217         gausva_2d = gset(ki,kj)
218         niset(ki,kj)  = 0
219
220      ENDIF
221       
222   END FUNCTION gausva_2d
223
224   FUNCTION psrandom()
225      !!----------------------------------------------------------------------
226      !!               ***  ROUTINE psrandom ***
227      !!
228      !! ** Purpose : Pseudo-Random number generator.
229      !!
230      !! ** Method  : Returns a pseudo-random number from a uniform distribution
231      !!              between 0 and 1
232      !!              Call with kdum a negative integer to initialize.
233      !!              Thereafter, do not alter kdum between successive deviates
234      !!              in sequence.
235      !!
236      !! ** Action  :
237      !!
238      !! History :
239      !!        !  10-02  (F. Vigilant)  Original code
240      !!        ! 2011-08 (D. Lea)  initialize with kdum1
241      !!----------------------------------------------------------------------
242      !! * Function return
243      REAL(wp) ::  &
244         & psrandom
245      !! * Arguments
246      INTEGER  :: &
247         & kdum          ! Seed
248      LOGICAL, SAVE :: &
249         & llinit = .TRUE.    ! Initialize on first call
250
251      ! Initialization
252      IF ( llinit ) THEN
253         kdum = 456953 + nproc * 596035
254         CALL init_mtrand(kdum)
255         llinit     = .FALSE.    ! On subsequent calls do not reinitialize
256      ENDIF
257
258      psrandom = mtrand_real1()
259
260   END FUNCTION psrandom
261
262
263END MODULE ran_num
Note: See TracBrowser for help on using the repository browser.