How to code something in Fortran. Monte Carlo error runs

[zombiecakes]

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