*IF DEF,SCMA S_RANDOM.2
C ******************************COPYRIGHT****************************** S_RANDOM.3
C (c) CROWN COPYRIGHT 1998, METEOROLOGICAL OFFICE, All Rights Reserved. S_RANDOM.4
C S_RANDOM.5
C Use, duplication or disclosure of this code is subject to the S_RANDOM.6
C restrictions as set forth in the contract. S_RANDOM.7
C S_RANDOM.8
C Meteorological Office S_RANDOM.9
C London Road S_RANDOM.10
C BRACKNELL S_RANDOM.11
C Berkshire UK S_RANDOM.12
C RG12 2SZ S_RANDOM.13
C S_RANDOM.14
C If no contract has been raised with this copy of the code, the use, S_RANDOM.15
C duplication or disclosure of it is strictly prohibited. Permission S_RANDOM.16
C to do so must first be obtained in writing from the Head of Numerical S_RANDOM.17
C Modelling at the above address. S_RANDOM.18
C ******************************COPYRIGHT****************************** S_RANDOM.19
!LL Purpose of this deck : provide a number of routines equivalent S_RANDOM.20
!LL to the random number generation routines of the NAG library. S_RANDOM.21
!LL Where possible, the interface of the routines/function S_RANDOM.22
!LL matches its equivalent in the NAG library. S_RANDOM.23
!LL S_RANDOM.24
!LL Written for T3E & HP platforms indifferently (then run at their S_RANDOM.25
!LL own default precision). S_RANDOM.26
!LL S_RANDOM.27
!LL Model Modification history : S_RANDOM.28
!LL version Date S_RANDOM.29
!LL 4.5 04.10.98 Deck introduced as part of the Single Column S_RANDOM.30
!LL Model S_RANDOM.31
!LL S_RANDOM.32
!LL Programming standards : Unified Model documentation paper no. 4 S_RANDOM.33
!LL S_RANDOM.34
!LL Logical components : S_RANDOM.35
!LL S_RANDOM.36
!LL System task : S_RANDOM.37
!LL S_RANDOM.38
!LL Documentation : Unified Model documentation paper no ?? S_RANDOM.39
!LL S_RANDOM.40
! S_RANDOM.41
!------------------------------------------------------------------ S_RANDOM.42
! Function inspired from the 'Numerical Recipies in fortran 77' S_RANDOM.43
! well known book. Generates a random number uniformly distributed S_RANDOM.44
! over [0,1]. S_RANDOM.45
! It should not be called direct externally from outside this S_RANDOM.46
! deck. Instead it is accessed through the G05xx S_RANDOM.47
! routines/functions when they need a random number. S_RANDOM.48
! S_RANDOM.49
! Arguments : nil S_RANDOM.50
! Returned value : a random number uniformly distributed on [0,1]. S_RANDOM.51
! Side effects : the status of the /ran_jc/ common block S_RANDOM.52
! is updated at each call. S_RANDOM.53
Function RAN1_JC() 2S_RANDOM.54
Implicit none S_RANDOM.55
Integer idum,ia,im, iq,ir ,ntab, ndiv S_RANDOM.56
Real ran1_jc, am, eps, rnmx S_RANDOM.57
Parameter (ia=16807, im=2147483647, am=1./im, iq=127773, S_RANDOM.58
& ir=2836, ntab=32, S_RANDOM.59
& ndiv = 1 + (im-1)/ntab, eps=1.2e-7, rnmx=1.0-eps) S_RANDOM.60
S_RANDOM.61
C *Minimal* random number generator of Park and Miller with S_RANDOM.62
C Bays-Durham shuffle and added safeguards. Returns a uniform S_RANDOM.63
C random deviate between 0.0 and 1.0 (exclusive of the endpoint S_RANDOM.64
C values). S_RANDOM.65
C Call with idum a negative integer to initialize; thereafter, do S_RANDOM.66
C not alter idum between successive deviates in a sequence. S_RANDOM.67
C RNMX should approximate the largest Floatingvalue that is less S_RANDOM.68
C than 1. S_RANDOM.69
S_RANDOM.70
Integer j, k, iv(ntab), iy S_RANDOM.71
S_RANDOM.72
Common /ran_jc/ idum,iv,iy ! This common block constitutes the S_RANDOM.73
! 'memory' of the timeserie ; S_RANDOM.74
! it is manipulated by other routines S_RANDOM.75
! for intializing, saving, restoring S_RANDOM.76
! the timeserie. S_RANDOM.77
S_RANDOM.78
If (idum.le.0.or.iy.eq.0) then ! Initialize. S_RANDOM.79
idum=max(-idum,1) ! Be sure to prevent idum = 0 S_RANDOM.80
Do j=ntab+8,1,-1 ! Load the shuffle table (after 8 S_RANDOM.81
! warm-ups). S_RANDOM.82
k=idum/iq S_RANDOM.83
idum=ia*(idum-k*iq)-ir*k S_RANDOM.84
If (idum.lt.0) idum=idum+im S_RANDOM.85
If (j.le.ntab) iv(j)=idum S_RANDOM.86
Enddo S_RANDOM.87
iy=iv(1) S_RANDOM.88
Endif S_RANDOM.89
S_RANDOM.90
k=idum/IQ ! Start here when not initializing. S_RANDOM.91
idum=ia*(idum-k*iq)-ir*k ! Compute idum=mod(IA*idum,IM) without S_RANDOM.92
! overflows by Schrage's method. S_RANDOM.93
If (idum.lt.0) idum=idum+im S_RANDOM.94
S_RANDOM.95
j=1+iy/ndiv ! Will be in the range 1:NTAB. S_RANDOM.96
iy=iv(j) ! Output previously stored value and S_RANDOM.97
! refill the shuffle table S_RANDOM.98
iv(j)=idum S_RANDOM.99
ran1_jc=min(am*iy,rnmx) ! Because users don't expect endpoint S_RANDOM.100
! values. S_RANDOM.101
S_RANDOM.102
Return S_RANDOM.103
End S_RANDOM.104
S_RANDOM.105
!---------------------------------------------------------------- S_RANDOM.106
! S_RANDOM.107
! Initialize the random number generator RAN1_JC (and therefore S_RANDOM.108
! also the functions like g05xxx which use RAN_JC) to give S_RANDOM.109
! repeatable sequence S_RANDOM.110
! S_RANDOM.111
! Argument : iseed input, not destroyed after the call S_RANDOM.112
! Side effect : initialize the /ran_jc/ common block S_RANDOM.113
! S_RANDOM.114
! This routine *needs* to be called first when using the S_RANDOM.115
! random routines G05xxE S_RANDOM.116
subroutine G05CBE(iseed) 1S_RANDOM.117
Implicit none S_RANDOM.118
! Argument S_RANDOM.119
Integer iseed S_RANDOM.120
! Local variables S_RANDOM.121
Integer ntab S_RANDOM.122
Parameter (ntab=32) S_RANDOM.123
Integer idum, iv(ntab), iy ! This common block constitutes the S_RANDOM.124
Common /ran_jc/ idum,iv,iy ! 'memory' of the random timeserie. S_RANDOM.125
S_RANDOM.126
Integer j ! index S_RANDOM.127
S_RANDOM.128
S_RANDOM.129
! ---------------------- S_RANDOM.130
! Zero the arrays of the ran1_jc function S_RANDOM.131
Do j = 1, ntab S_RANDOM.132
iv(j) = 0 S_RANDOM.133
Enddo S_RANDOM.134
iy = 0 S_RANDOM.135
! initialize the ran1_jc function with the user provided value: S_RANDOM.136
idum = - int(iseed) ! idum must be negative when S_RANDOM.137
! initializing RAN1_JC. S_RANDOM.138
S_RANDOM.139
Return S_RANDOM.140
End S_RANDOM.141
S_RANDOM.142
!---------------------------------------------------------------- S_RANDOM.143
! S_RANDOM.144
! Save the state of the random number generator S_RANDOM.145
! S_RANDOM.146
! Arguments : they are all output value. It is the S_RANDOM.147
! task of the user to store those values somewhere S_RANDOM.148
! until s/he wants to reinitialise to an identical S_RANDOM.149
! serie of random numbers with G05CGE, called S_RANDOM.150
! with the same arguments. S_RANDOM.151
! Side effect : nil S_RANDOM.152
Subroutine G05CFE(idum_out,iv_out,iy_out) 2S_RANDOM.153
Implicit none S_RANDOM.154
S_RANDOM.155
Integer ntab S_RANDOM.156
Parameter (ntab=32) S_RANDOM.157
! Arguments to extract from the memory of the random routines S_RANDOM.158
Integer idum_out,iv_out(ntab),iy_out S_RANDOM.159
! Local variables S_RANDOM.160
Integer idum,iv(ntab),iy S_RANDOM.161
Common /ran_jc/ idum,iv,iy ! This common block constitutes the S_RANDOM.162
! 'memory' of the timeserie. S_RANDOM.163
Integer j S_RANDOM.164
! S_RANDOM.165
S_RANDOM.166
! Copy the values straight from the common into the _out S_RANDOM.167
! variables. S_RANDOM.168
Do j = 1, ntab S_RANDOM.169
iv_out(j) = iv(j) S_RANDOM.170
Enddo S_RANDOM.171
idum_out = idum S_RANDOM.172
iy_out = iy S_RANDOM.173
S_RANDOM.174
Return S_RANDOM.175
End S_RANDOM.176
S_RANDOM.177
S_RANDOM.178
!---------------------------------------------------------------- S_RANDOM.179
! S_RANDOM.180
! Restore the state of the random number generator S_RANDOM.181
! S_RANDOM.182
! Arguments : they are all input values. S_RANDOM.183
! They are the values which should be extracted S_RANDOM.184
! from the memory of the random routines by G05CFE S_RANDOM.185
! The arguments are not destroyed after the call. S_RANDOM.186
! Side effect : the /ran_jc/ common block is S_RANDOM.187
! updated with the values provided as argument. S_RANDOM.188
Subroutine G05CGE(idum_in,iv_in,iy_in) 1S_RANDOM.189
Implicit none S_RANDOM.190
S_RANDOM.191
Integer ntab S_RANDOM.192
Parameter (ntab=32) S_RANDOM.193
S_RANDOM.194
! arguments to put into the memory of the random routines S_RANDOM.195
Integer idum_in,iv_in(ntab),iy_in S_RANDOM.196
! Local variables S_RANDOM.197
Integer idum,iv(ntab),iy S_RANDOM.198
common /ran_jc/ idum,iv,iy ! This common block constitutes the S_RANDOM.199
! 'memory' of the timeserie. S_RANDOM.200
Integer j ! index S_RANDOM.201
S_RANDOM.202
! Copy the values straight into the common S_RANDOM.203
Do j = 1, ntab S_RANDOM.204
iv(j) = iv_in(j) S_RANDOM.205
Enddo S_RANDOM.206
idum = idum_in S_RANDOM.207
iy = iy_in S_RANDOM.208
S_RANDOM.209
Return S_RANDOM.210
End S_RANDOM.211
S_RANDOM.212
!---------------------------------------------------------------- S_RANDOM.213
! S_RANDOM.214
! Pseudo-random real numbers, Gaussian distribution of S_RANDOM.215
! mean m, standard deviation s. S_RANDOM.216
! S_RANDOM.217
! Also inspired from the 'Numerical Recipies in fortran 77' book. S_RANDOM.218
! S_RANDOM.219
! Arguments : Input (m,s) ; are not destroyed after the call. S_RANDOM.220
! calls RAN_JC() in order to access a uniform random number S_RANDOM.221
! generator over [0,1]. S_RANDOM.222
! Side effect : changes the values of the /ran_jc/ common S_RANDOM.223
! block S_RANDOM.224
Function G05DDE(m,s) 9,2S_RANDOM.225
Implicit none S_RANDOM.226
! Function type S_RANDOM.227
Real G05DDE S_RANDOM.228
! Arguments S_RANDOM.229
Real m,s ! (mean, standard deviation) of S_RANDOM.230
! the distribution. S_RANDOM.231
S_RANDOM.232
! Uses RAN1_JC S_RANDOM.233
External RAN1_JC S_RANDOM.234
Real RAN1_JC S_RANDOM.235
S_RANDOM.236
! Returns a normally distributed deviate with m mean and s S_RANDOM.237
! variance. Uusing ran1_jc() as the source of uniform deviates. S_RANDOM.238
S_RANDOM.239
! Local variables S_RANDOM.240
Integer iset S_RANDOM.241
Real fac,gset,rsq,v1,v2 S_RANDOM.242
Save iset, gset S_RANDOM.243
Data iset/0/ S_RANDOM.244
! S_RANDOM.245
If (iset.eq.0) then ! We don't have an extra deviate S_RANDOM.246
1 v1=2.*ran1_jc
()-1. ! handy, so pick two uniform numbers S_RANDOM.247
v2=2.*ran1_jc()-1. ! in the square extending from -1 to S_RANDOM.248
rsq=v1**2+v2**2 ! +1 in each direction see if they S_RANDOM.249
If (rsq.ge.1. .or. rsq.eq.0.) goto 1 ! are in the unit circle. S_RANDOM.250
! and if they are not, try again. S_RANDOM.251
fac=sqrt(-2.*log(rsq)/rsq) ! Now make the Box-Muller S_RANDOM.252
! transformation to get two normal S_RANDOM.253
! deviates. Return one and save the S_RANDOM.254
! other for next time. S_RANDOM.255
gset=v1*fac S_RANDOM.256
G05DDE=m+s*(v2*fac) S_RANDOM.257
iset=1 ! Set flag. S_RANDOM.258
else ! We have an extra deviate handy. S_RANDOM.259
G05DDE=m+s*gset ! so return it. S_RANDOM.260
iset=0 ! and unset the flag. S_RANDOM.261
endif S_RANDOM.262
S_RANDOM.263
return S_RANDOM.264
end S_RANDOM.265
*ENDIF S_RANDOM.266