module number_generators implicit none integer, private :: idum, inext, inextp, inext1, inextp1 real*8, private, dimension(55) :: ma, ma2 logical :: switch CONTAINS !*********************************************************************** ! ran3: generates random values 0..1] ! seed (idum) is set to a negative integer on entry ! From Press et al.,(1992) Numerical recipes in FORTRAN (2nd ed.) CUP ! ! This routine runs through the ma vector pre-established when the ! generator is initialized !*********************************************************************** function ran3(idum) implicit none real*8 :: ran3, mj, mk real*8, parameter :: MBIG =4000000.0, MSEED =1618033.0, MZ =0.0, FAC = 1.0/MBIG integer :: i, ii, k, idum inext = inext + 1 if(inext == 56) inext = 1 inextp = inextp + 1 if(inextp == 56)inextp= 1 mj = ma(inext) - ma(inextp) if(mj < MZ) mj = mj + MBIG ma(inext) = mj ran3 = mj*FAC return end function ran3 !********************************************************************* ! rnorm: Unit Normal distribution !********************************************************************* function rnorm(idum) implicit none real*8 :: fac, r, v1, v2, rnorm real*8, save :: rnorm2 integer :: idum if (switch) then 10 v1 = 2.0 * ran3(idum) - 1.0 v2 = 2.0 * ran3(idum) - 1.0 r = v1**2 + v2**2 if ((r .ge. 1.0).or.(r .eq.0)) goto 10 fac = sqrt(-2.0 * log(r)/r) rnorm2 = v1 * fac switch = .false. rnorm = v2 * fac else switch = .true. rnorm = rnorm2 endif return end function rnorm !*********************************************************************** ! gaussian: returns a Normal deviate from N(mu, sd) !*********************************************************************** function gaussian(mu, sd) real*8, intent(in) :: mu, sd real*8 :: gaussian gaussian = (rnorm(idum) * sd) + mu ! Calculate Normal(mu, sigma) return end function gaussian !********************************************************************** ! FixSeed: Sets a fixed seed (-99) for the random-number routine ! and sets up the ma vector for the random number generator !********************************************************************** subroutine FixSeed real*8 :: mj, mk real*8, parameter :: MBIG =4000000.0, MSEED =1618033.0, MZ =0.0 integer :: i, ii, k idum = -99 write(*,'(A, I3)') ' Seed for random-number generator = ', idum mj= MSEED - iabs(idum) mj = mod(mj, MBIG) ma(55) = mj mk=1 do i=1,54 ii = mod(21*i, 55) ma(ii) = mk mk = mj - mk if(mk < MZ) mk = mk + MBIG mj = ma(ii) enddo do k = 1, 4 do i = 1, 55 ma(i) = ma(i) - ma(1 + mod(i+30, 55) ) if(ma(i) .lt. MZ) ma(i) = ma(i) + MBIG enddo enddo inext = 0 inextp = 31 return end subroutine FixSeed !********************************************************************** ! RandSeed: Sets a random seed for the random-number routines ! and sets up the ma vector for the random number generator !********************************************************************** subroutine RandSeed implicit none real*8 :: mj, mk real*8, parameter :: MBIG =4000000.0, MSEED =1618033.0, MZ =0.0 integer :: i, ii, k call system_clock(idum) 10 if (idum > 1000000) idum = idum / 10 if (idum > 1000000) go to 10 if (idum > 0) idum = idum *(-1) ! write(*,'(A, I8)') ' Seed for random-number generator = ', idum mj= MSEED - iabs(idum) mj = mod(mj, MBIG) ma(55) = mj mk=1 do i=1,54 ii = mod(21*i, 55) ma(ii) = mk mk = mj - mk if(mk < MZ) mk = mk + MBIG mj = ma(ii) enddo do k = 1, 4 do i = 1, 55 ma(i) = ma(i) - ma(1 + mod(i+30, 55) ) if(ma(i) .lt. MZ) ma(i) = ma(i) + MBIG enddo enddo inext = 0 inextp = 31 return end subroutine RandSeed !********************************************************************** ! get_ran_seed: gets current state of the random-number generator !********************************************************************** subroutine get_ran_seed(dummy, MAA) integer, intent(out), dimension(2) :: dummy real*8, dimension(55) :: MAA dummy(1) = inext dummy(2) = inextp MAA = ma return end subroutine get_ran_seed !********************************************************************** ! assign_seed: restores the state of the random-number generator !********************************************************************** subroutine assign_seed(iseed, MAA) integer, dimension(2), intent(in) :: iseed integer :: get_ran_seed real*8, dimension(55) :: MAA inext = iseed(1) inextp = iseed(2) ma = MAA switch = .true. return end subroutine assign_seed END MODULE Number_Generators