Stuck on how to code something in Fortran. Monte Carlo simulation

Hey,

I am doing an atmospheric chemistry PhD and I use Fortran to model the reactions that I do in a flow tube.
I'm not really the best programmer, I can alter the model to my needs, but I wouldn't be able to write anything from scratch.

I use Monte Carlo simulations on my reactions in order to work out the errors in the rate coefficients. At the moment i use the code below and it does a "top hat" shaped distribution. Obviously there is more to model, as it actually models the fit to my data as well etc.

Do p=1,N
! here are parameters

Call random_number(A)

k1=k1o+(0.5-A(1))*k1width
k12=k12o+(0.5-A(2))*k12width
k13=k13o+(0.5-A(3))*k13width


Zed=(k2-Avk2)

If (zed.lt.0) then

zedn=zed

zednsq=zedn**2
Sumn=Sumn+zednsq

Xn=Xn+1

else

zedp=zed

zedpsq=zedp**2
Sump=Sump+zedpsq

Xp=Xp+1

endif

Enddo

mcerrorn=(SUMn/(Xn-1))**0.5
mcerrorp=(SUMp/(Xp-1))**0.5

write(9,*) 'N = ',N

write(*,*) 'error n = ',mcerrorn
write(9,*) 'error n = ',mcerrorn
write(9,*) 'sum n = ',Xn

write(*,*) 'error p = ',mcerrorp
write(9,*) 'error p = ',mcerrorp
write(9,*) 'sum p = ',Xp


Basically I work with rate coefficients and concentrations so there is a variable for each parameter in my model that has an error associated with it. The model calculates the overall error


However I need a normal distribution. I have been given some code by my supervisor (obviously not in Fortran...) but I have no idea about how to put it into my model.

17 x = mean + 25*sd*(0.5-RND)*2
p = exp(-(x-mean)^2/2/sd^2)
if RND > p the goto 17

What my supervisor said about the above code: I'm to assume the parameter x has a best fit value "mean" and is normally distributed with a standard deviation (or standard error ("sd").
In line 17 x is chosen to be somewhere within +/- 25 standard deviations of the mean (RND is a random number between 0 and 1). The probability p of x with respect to the mean is then calculated according to the normal distribution. A second random number is the generated. If RND < p, then the resulting value of x is used. Otherwise the routine loops back to line 17 and repeats.

Sorry if this makes sense of I haven't included enough information. Please feel free to ask me any questions.

Thanks in advance for any help
Lottes

Sign In or Register to comment.

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Categories