module my_type INTEGER, PARAMETER :: IB = SELECTED_INT_KIND(9) !integer*4 INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6,37) !real*4 INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,307) !real*8 end module my_type program planck ! MC integration of the dimensionless Planck-function ! note: for n>1e6, sum1 and sum2 must be double prec! ! otherwise the rounding errors become too large, & ! since the manissa has only 6.5 digits use my_type implicit none integer(ib) :: seed integer(ib) :: n,i real(dp) :: sum1,sum2 real(sp) :: xmax,sigma,exact ! calculate the exact value of the integral exact= seed=5 print*,'give in sample size N and upper integration limit xmax' read*,n,xmax sum1=0. sum2=0. do i=1,n ! draw random variable (account for range different from 0...1) ! calculate function ! compute sums for MC-estimator of integral (sum1) and ! variance (sum2) enddo ! calculate final values for MC-integral (sum1) and error (sigma) print*,'n = ',n,' MC-integral = ',sum1,' +/- ',sigma,' exact= ',exact end