pkg/darwin/darwin_random.F

8fbfd1f382 Oliv*0001 #include "DARWIN_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: DARWIN_RANDOM_INIT
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE DARWIN_RANDOM_INIT(seed, myThid)
                0009
                0010 C     !DESCRIPTION:
                0011 C     Initializes the random number generator.
                0012 C     - seed must be positive.
                0013 C     - uses portable random number generator of Knuth [see Numerical
                0014 C       Recipes, Ch.7.1: ran3].  We use the floating-point version.
                0015 C       In order to obtain unique sequences of random numbers, the seed
                0016 C       should be between 1 and 1618032.
                0017 C     - NOTE: not thread-safe (yet)!!!
                0018
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021 #include "EEPARAMS.h"
                0022
                0023 C     !INPUT/OUTPUT PARAMETERS:
                0024 C     seed   :: seed for random number generator; must be positive
                0025 C     myThid :: thread number
                0026       INTEGER seed
                0027       INTEGER myThid
                0028 CEOP
                0029
                0030 #ifdef ALLOW_DARWIN
                0031
                0032 C     !FUNCTIONS:
                0033       real*8 port_rand
                0034       external port_rand
                0035
                0036 C     !LOCAL VARIABLES:
                0037 C     msgBuf :: informational/error meesage buffer
                0038       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0039       _RL RandNo
                0040       REAL*8 Dseed
                0041       CHARACTER*16 random_name
                0042
                0043       IF (myThid .GT. 1) THEN
                0044         CALL PRINT_ERROR('DARWIN_RANDOM_INIT: threading no supported',
                0045      &      myThid)
                0046         STOP 'ABNORMAL END: S/R DARWIN_RANDOM_INIT'
                0047       END IF
                0048
                0049       IF (seed .LE. 0) THEN
                0050         CALL PRINT_ERROR('DARWIN_RANDOM_INIT: seed must be positive'
                0051      &                  , myThid)
                0052       END IF
                0053
                0054       Dseed = float(seed)
                0055       RandNo = port_rand(Dseed)
                0056 C     need to call again to get a non-zero random number
                0057       Dseed = -1.D0
                0058       RandNo = port_rand(Dseed)
                0059       random_name = 'port_rand'
                0060
                0061       WRITE(msgbuf,'(A,A,I10,X,F20.16)')
                0062      &   'DARWIN_RANDOM_INIT: ', random_name, seed, RandNo
                0063       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0064      &                    SQUEEZE_RIGHT , myThid)
                0065
                0066 #endif
                0067
                0068       RETURN
                0069       END
                0070
                0071 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0072 CBOP
                0073 C     !ROUTINE: DARWIN_RANDOM
                0074
                0075 C     !INTERFACE:
                0076       FUNCTION DARWIN_RANDOM(myThid)
                0077
                0078 C     !DESCRIPTION:
                0079 C     Returns a uniform random number between 0 and 1
                0080
                0081 C     !USES:
                0082       IMPLICIT NONE
                0083 #include "EEPARAMS.h"
                0084
                0085 C     !INPUT/OUTPUT PARAMETERS:
                0086 C     DARWIN_RANDOM :: uniform random number
                0087 C     myThid     :: thread number
                0088       _RL DARWIN_RANDOM
                0089       INTEGER myThid
                0090 CEOP
                0091
                0092 #ifdef ALLOW_DARWIN
                0093
                0094 C     !FUNCTIONS:
                0095       real*8 port_rand
                0096       external port_rand
                0097
                0098 C     !LOCAL VARIABLES:
                0099       real*8 Dseed
                0100
                0101       IF (myThid .GT. 1) THEN
                0102         CALL PRINT_ERROR('DARWIN_RANDOM: threading no supported',
                0103      &      myThid)
                0104         STOP 'ABNORMAL END: S/R DARWIN_RANDOM'
                0105       END IF
                0106
                0107       Dseed = -1.d0
                0108       darwin_random = port_rand(Dseed)
                0109
                0110 #endif
                0111
                0112       RETURN
                0113       END
                0114
                0115 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0116 CBOP
                0117 C     !ROUTINE: DARWIN_RANDOM_NORMAL
                0118
                0119 C     !INTERFACE:
                0120       FUNCTION DARWIN_RANDOM_NORMAL(myThid)
                0121
                0122 C     !DESCRIPTION:
                0123 C     Returns a normally distributed random number with mean 0
                0124 C     and stddev 1
                0125
                0126 C     !USES:
                0127       IMPLICIT NONE
                0128 #include "EEPARAMS.h"
                0129
                0130 C     !INPUT/OUTPUT PARAMETERS:
                0131 C     DARWIN_RANDOM_NORMAL :: normally distributed random number
                0132 C     myThid            :: thread number
                0133       _RL DARWIN_RANDOM_NORMAL
                0134       INTEGER myThid
                0135 CEOP
                0136
                0137 #ifdef ALLOW_DARWIN
                0138
                0139 C     !FUNCTIONS:
                0140       real*8 port_rand_norm
                0141       external port_rand_norm
                0142
                0143 C     !LOCAL VARIABLES:
                0144
                0145       IF (myThid .GT. 1) THEN
                0146         CALL PRINT_ERROR('DARWIN_RANDOM: threading no supported',
                0147      &      myThid)
                0148         STOP 'ABNORMAL END: S/R DARWIN_RANDOM'
                0149       END IF
                0150
                0151       darwin_random_normal = port_rand_norm()
                0152
                0153 #endif
                0154
                0155       RETURN
                0156       END
                0157 c ==========================================================
                0158
                0159 c Inverse normal distribution
                0160 c returns inverse normal cumulative distribution
                0161 c from p:[0,1] -> y:[-inf,+inf] centered on mu with stdev of sigma
                0162 c p is RandNo passed in, y is return variable for deviate
                0163 c
                0164 c  Scott Grant, Spring 2006
                0165
                0166       SUBROUTINE DARWIN_INVNORMAL(y,p,mean,sigma)
                0167       implicit none
                0168
                0169 c local variables
                0170       real*8 mean
                0171       real*8 sigma
                0172       real*8 q
                0173       real*8 r
                0174       real*8 x
                0175       real*8 p
                0176       real*8 plow
                0177       real*8 phigh
                0178       real*8 y
                0179       real*8 a(6)
                0180       real*8 b(5)
                0181       real*8 c(6)
                0182       real*8 d(4)
                0183
                0184
                0185 c Create random variable from -inf to +inf
                0186 c Coefficients in rational approximations.
                0187       a(1) = -3.969683028665376d+01
                0188       a(2) =  2.209460984245205d+02
                0189       a(3) = -2.759285104469687d+02
                0190       a(4) =  1.383577518672690d+02
                0191       a(5) = -3.066479806614716d+01
                0192       a(6) =  2.506628277459239d+00
                0193
                0194       b(1) = -5.447609879822406d+01
                0195       b(2) =  1.615858368580409d+02
                0196       b(3) = -1.556989798598866d+02
                0197       b(4) =  6.680131188771972d+01
                0198       b(5) = -1.328068155288572d+01
                0199
                0200       c(1) = -7.784894002430293d-03
                0201       c(2) = -3.223964580411365d-01
                0202       c(3) = -2.400758277161838d+00
                0203       c(4) = -2.549732539343734d+00
                0204       c(5) =  4.374664141464968d+00
                0205       c(6) =  2.938163982698783d+00
                0206
                0207       d(1) =  7.784695709041462d-03
                0208       d(2) =  3.224671290700398d-01
                0209       d(3) =  2.445134137142996d+00
                0210       d(4) =  3.754408661907416d+00
                0211
                0212 c  Define break-points.
                0213
                0214       plow  = 0.02425d0
                0215       phigh = 1.d0 - plow
                0216
                0217 c  Rational approximation for lower region.
                0218
                0219       if ((0.d0 .lt. p) .and. (p .lt. plow))then
                0220          q = sqrt(-2.0d0*log(p))
                0221          x = (((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
                0222      &      ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1.d0)
                0223       endif
                0224
                0225 c  Rational approximation for central region.
                0226
                0227       if ((plow .le. p).and.(p .le. phigh))then
                0228          q = p - 0.5d0
                0229          r = q*q
                0230          x = (((((a(1)*r+a(2))*r+a(3))*r+a(4))*r+a(5))*r+a(6))*q /
                0231      &     (((((b(1)*r+b(2))*r+b(3))*r+b(4))*r+b(5))*r+1.d0)
                0232       endif
                0233
                0234 c  Rational approximation for upper region.
                0235
                0236       if ((phigh .lt. p).and.(p .lt. 1.d0))then
                0237          q = sqrt(-2.0d0*log(1.d0-p))
                0238          x = -(((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
                0239      &      ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1.d0)
                0240       endif
                0241
                0242 c Normal Deviate about mean
                0243 c        write(6,*)'DEVIATE',x
                0244       y = sigma*sqrt(2.0d0)*x + mean
                0245 c        write(6,*)'Normal PDF Value INSIDE:',y
                0246 c        write(6,*)'MEAN:',mean
                0247 c        write(6,*)'SIGMA:',sigma
                0248
                0249       RETURN
                0250       END