MODULE probMod

IMPLICIT NONE
SAVE
PRIVATE
PUBLIC :: uniRan, poissonRan, normalRan, getFrequency1D
CONTAINS
REAL(4) FUNCTION ranGen (idum)
INTEGER :: idum, im1, im2, imm1, ia1, ia2, iq1, iq2, ir1, ir2, ntab, ndiv, j, k
REAL(4) :: am, eps, rnmx
PARAMETER (im1=2147483563, im2=2147483399, am=1.0/im1, imm1=im1-1, &
ia1=40014, ia2=40692, iq1=53668, iq2=52774, &
ir1=12211, ir2=3791, &
ntab=32, ndiv=1+imm1/ntab, eps=1.2e-7, rnmx=1.0-eps)
INTEGER, SAVE :: iv(32), idum2=123456789, iy=0
DATA iv/ntab*0/
!------
! Initialise
IF (idum <= 0) THEN
idum = MAX(-idum,1) ! Be sure to prevent idum=0
idum2 = idum
DO j = ntab+8, 1, -1 ! Load the shuffle table after 8 warm-ups
k = idum/iq1
idum = ia1*(idum-k*iq1) - k*ir1
IF (idum < 0) idum = idum + im1
IF (j <= ntab) iv(j) = idum
END DO
iy = iv(1)
END IF
!
! Start here when not initialising
k = idum/iq1
!
! Compute idum=mod(ia2*idum2,im1) without overflow by Schragge's method
idum = ia1*(idum-k*iq1) - k*ir1
IF (idum < 0) idum = idum + im1
k = idum2/iq2
!
! Compute idum=mod(ia2*idum2,im2), likewise
idum2 = ia2*(idum2-k*iq2) - k*ir2
IF (idum2 < 0) idum2 = idum2 + im2
j = 1 + iy/ndiv ! Will be in the range 1:ntab
!
! Here idum is shuffled. idum and idum2 are combined to generate output
iy = iv(j) - idum2
iv(j) = idum
IF (iy < 1) iy = iy + imm1
!
! Because users don't expect endpoint values
ranGen = MIN (am*iy, rnmx)
END FUNCTION ranGen
!
! Uniform random number generator
REAL(4) FUNCTION uniRan ()
!------
! Uniform 0,1 random number generator
!
! Based on code in 'Numerical Recipes' by Press et al. (1986)
!----
INTEGER, SAVE :: idum
LOGical, SAVE :: initSeed=.true.
INTEGER, PARAMETER :: initialDum=-2*153351
!------
! Initialize on first CALL
IF (initSeed) THEN
initSeed = .false.
idum = initialDum
IF (idum > 0) idum = -idum
uniRan = ranGen (idum)
END IF
!------
! Generate a ranDOm number
uniRan = ranGen (idum)
END FUNCTION uniRan
REAL(4) FUNCTION poissonRan (xm)
!------
! Poisson random number generator
!
! Based on code in 'Numerical Recipes' by Press et al. (1986)
!----
REAL(4), INTENT(in) :: xm
REAL(4), PARAMETER :: pi=3.141592654
REAL(4) :: em, t, y
REAL(4), SAVE :: alxm, g, sq, oldm=-1.0
!------
IF (xm < 12.0) THEN
IF (xm /= oldm) THEN
oldm = xm
g = EXP(-xm)
END IF
em = -1
t = 1.0
DO
em = em + 1.0
t = t*uniRan()
IF (t <= g) EXIT
END DO
ELSE
IF (xm /= oldm) THEN
oldm = xm
sq = SQRT(2.0*xm)
alxm = LOG(xm)
g = xm*alxm - gammln(xm+1.0)
END IF
DO
y = TAN(pi*uniRan())
em = sq*y + xm
IF (em < 0.0) CYCLE
em = INT(em)
t=0.9*(1.0+y**2)*EXP(em*alxm-gammln(em+1.0)-g)
IF (uniRan() <= t) EXIT
END DO
END IF
poissonRan = em
END FUNCTION poissonRan
REAL(4) FUNCTION gammln (xx)
REAL(4), INTENT(in) :: xx
INTEGER :: j
REAL(8) :: ser, tmp, x, y
REAL(8), PARAMETER :: cof(6)=(/76.18009172947146d0,-86.50532032941677d0, &
24.01409824083091d0, -1.231739572450155d0, &
.1208650973866179d-2, -.5395239384953d-5/), &
stp=2.5066282746310005d0
!-----
x = xx
y = x
tmp = x + 5.5d0
tmp = (x+0.5d0)*LOG(tmp) - tmp
ser = 1.000000000190015d0
DO j = 1, 6
y = y + 1.d0
ser = ser + cof(j)/y
END DO
gammln = tmp + LOG(stp*ser/x)
END FUNCTION gammln
! Function produces normal deviate corresponding to lower tail area of p
! ifault = 0 means ok
! = 1 means p is not in range (0,1)
REAL(8) FUNCTION ppnd (p, ifault)
!-----
! Algorithm AS 111 Applied Statistics (1977) vol 26, no.1
! "The percentage points of the normal distribution" by J.D. Beasley and S.G. Sp
ringer
IMPLICIT NONE
INTEGER :: ifault
REAL(8),PARAMETER :: zero=0.0d0, split=0.42d0, half=0.5d0, one=1.0d0
REAL(8) :: p, q, r
REAL(8), PARAMETER :: a0 = 2.50662823884d0, a1 = -18.61500062529d0, &
a2 = 41.39119773534d0, a3 = -25.44106049637d0, &
b1 = -8.47351093090d0, b2 = 23.08336743743d0, &
b3 = -21.06224101826d0, b4 = 3.13082909833d0, &
c0 = -2.78718931138d0, c1 = -2.29796479134d0, &
c2 = 4.85014127135d0, c3 = 2.32121276858d0, &
d1 = 3.54388924762d0, d2 = 1.63706781897d0
!-----
ifault = 0
q = p - half
IF (ABS(q).le.split) THEN
r = q*q
ppnd = q*(((a3*r+a2)*r+a1)*r+a0) / ((((b4*r+b3)*r+b2)*r+b1)*r+one)
ELSE
r = p
IF (q > zero) r = one - p
IF (r <= zero) THEN
ifault = 1
ppnd = zero
RETURN
END IF
r = SQRT (-DLOG(r))
ppnd = (((c3*r+c2)*r+c1)*r+c0) / ((d2*r+d1)*r+one)
IF (q < zero) ppnd = -ppnd
END IF
END FUNCTION ppnd
!
! Function produces normal deviate with mean and stdDev
REAL(4) FUNCTION normalRan (mean, stdDev)
IMPLICIT NONE
REAL(4), INTENT(in) :: mean, stdDev
REAL(8) :: prob
INTEGER :: fault
!----
prob = DBLE(uniRan())
IF (stdDev <= 0.0) THEN
WRITE(*,*) 'f-normal/Negative standard deviation'
STOP
ELSE
normalRan = mean + stdDev*ppnd (prob,fault)
END IF
END FUNCTION normalRan
!------------------------------
pure subroutine getFrequency1D(x,xminIn,xmaxIn,xfreq,freq,nOutsideRange,err,mess
age)
! Purpose: Given (x) samples, computes frequency on a specified 1D mesh.
! Programmer: Dmitri Kavetski, University of Newcastle, 2002.
implicit none
! Dummies
real,intent(in)::x(:)
real,intent(in),optional::xminIn,xmaxIn
real,intent(out)::xfreq(:)
integer,intent(out)::freq(:)
integer,intent(out),optional::nOutsideRange ! number of samples outside of range
integer,intent(out)::err
character(*),intent(out)::message
! Locals
integer::nXbins,i,n,k,nout_loc
real::dx,xmin,xmax
! Start procedure here
n=size(x); nXbins=size(freq)
if(present(xminIn))then
xmin=xminIn
else
xmin=minval(x)
endif
if(present(xmaxIn))then
xmax=xmaxIn
else
xmax=maxval(x)
endif
if(xmin>=xmax)then
err=100; message="f-getFrequency1D/xmin>=xmax"
return
elseif(nXbins<1)then
err=200; message="f-getFrequency1D/nXbin<1"
return
elseif(nXbins/=size(xfreq))then
err=300; message="f-getFrequency1D/arrayDimError"
return
endif
! get frequency within bins of width dx and frequency outside plot area
freq=0; nout_loc=0; dx=(xMax-xMin)/real(nXbins)
do i=1,nXbins
xfreq(i)=xmin+0.5*dx+real(i-1)*dx
enddo
do k=1,n
if(xMin<=x(k).and.x(k)<=xMax)then
i=min(nXbins,int((x(k)-xMin)/dx)+1)
freq(i)=freq(i)+1
else
nout_loc=nout_loc+1
endif
enddo
if(present(nOutsideRange))nOutsideRange=nout_loc
err=0; message="getFrequency1D/ok"
! End procedure here
endsubroutine getFrequency1D
!------------------------------
END MODULE probMod
!=========================================================================
program tute10
use probMod
implicit none
integer,parameter::nSamp=1000000
integer::i,nDamage
integer::nSamp
real,dimension(nSamp)::damages

real,parameter::meanFlood=42000
real,parameter::stdFlood=20000
real,parameter::FloodThres=80000
real::lambda,zeta
real::Q
integer,parameter::unt=1,nHisto=50
real::xfreq(nHisto)
integer::freq(nHisto)
integer::err
character(64)::message
!write(*,*)"please enter number of samples to be trialled"
!read(*,*) nSamp
!write(*,*)"the number of samples chosen is",nSamp
!write(*,*)"press enter to terminate"
!read(*,*)
zeta = sqrt(log(1.0+(stdFlood/meanFlood)**2))
lambda = LOG(meanFlood)-zeta*zeta/2.0
nDamage=0
do i=1,nSamp
Q=exp(normalRan(lambda,zeta))
if(Q>FloodThres)then
nDamage=nDamage+1
damages(nDamage)=(10000*(Q-FloodThres)**2)
endif
enddo
open(unit=123,file="flooddamagethistogram.out",status="unknown",action="write",i
ostat=err)
if(err/=0)then
write(*,*)"Output file open error"
stop
endif
write(123,'("Probability of Damage: ",f)')real(nDamage)/nSamp
damages=damages/1e9
call getFrequency1D(x=damages(1:nDamage),xminIn=0.0,xmaxIn=50000.0,xfreq=xfreq,f
req=freq,err=err,message=message)
write(123,*)
write(123,'(a,i0,a)')"Distribution of 'damage' after ",nSamp," samples"
write(123,'(a4,2x,a10,2x,a6)')"i","damage","count"
do i=1,nHisto
write(123,'(i4,2x,es10.3,2x,i6)')i,xfreq(i),freq(i)
enddo
write(123,*)"outside=",nDamage-sum(freq)
close(123)
end program