rand_normal.f90

SUBROUTINE RAND_NORMAL(NUM,Z,SEEDNUM)
!
! ----------------------------------------------------------------------
!
!  GENERATES 'NUM' RANDOM NORMAL NUMBERS WITH ZERO MEAN AND UNIT VARIANCE
!
!  ALGORITHM FROM LEVA, 1992, IN ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE 18(4)
!  FOR SUCCINT DESCRIPTION, SEE THE LIST OF STEPS ON PP. 451-452.
!
!  THIS CODE IS MODIFIED TO ALLOW FOR FUNCTIONALITY WITH MULTIPLE PROCESSORS.
!  THE MODIFICATION IS ESSENTIALLY ONE EXTRA INPUT WHICH IS USED TO SEED
!  THE RANDOM GENERATOR IN SUCH A WAY THAT IT WILL RETURN DIFFERENT RANDOM
!  VARIATES THAN EACH OTHER PROCESSOR.
!
!  CALLS: NONE
!
!  VERSION HISTORY: WRITTEN - MD - MAY 2005
!                   1.1     - MD - AUG 2005 MODIFIED FOR MULTIPLE PROCESSORS
!                   1.2     - MD & RK - MAR 2014 MODIFIED FOR DYNAMIC SEED SIZE
!
!  NUM - NUMBER OF VARIATES TO GENERATE
!  Z - OUTPUT ARRAY OF VARIATES
!  SEEDNUM - SINGLE NUMBER USED IN SEED FOR THE RANDOM NUMBER GENERATOR

IMPLICIT NONE

INTEGER,INTENT(IN)::NUM,SEEDNUM
REAL,INTENT(OUT)::Z(NUM)
INTEGER N,RSZ,sout(2),I
INTEGER, ALLOCATABLE :: SEED(:)
REAL U,V,S,T,A,B,X,Y,Q

! DEFINE CONSTANTS
S=0.449871
T=-0.386595
A=0.19600
B=0.25472

! RSZ IS THE NUBMER OF INTEGERS USED TO HOLD THE SEED.  THIS IS OBTAINED
! WITH THE CALL TO RANDOM_SEED USING THE 'SIZE' ARGUMENT.
CALL RANDOM_SEED(SIZE=RSZ)
ALLOCATE(SEED(RSZ))

! NOW THE SEED IS COMPUTED BY ADDING THE PROCESS # AND THE NUMBER I=1,RSZ
DO I=1,RSZ
  SEED(I)=SEEDNUM+I
END DO

! FINALLY THE SEED IS SET USING THE SEED VARIABLE
CALL RANDOM_SEED(PUT=SEED)

N=0
DO WHILE(N<NUM)
  ! 1. GENERATE A POINT P=(U,V) UNIFORM IN R
  CALL RANDOM_NUMBER(U)
  CALL RANDOM_NUMBER(V)
  V=1.7156*(V-0.5) !SCALE V TO BE UNIFORM IN (-R,R) WHERE R=2/e = 0.8578

  ! 2. EVALUATE QUADRATIC FORM OF Q
  X=U-S
  Y=ABS(V)-T
  Q=X**2+Y*(A*Y-B*X)

  ! 3. DECIDE WHETHER TO ACCEPT OR REJECT THIS POINT
  IF (Q<0.27597) THEN
    !ACCEPT P IF Q IS LESS THAN INNER BOUNDARY VALUE
    N=N+1
    Z(N)=V/U
  ELSEIF (Q>0.27846.OR.V**2>-4*U**2*LOG(U)) THEN
    !REJECT THIS ONE IF IF Q IS GREATER THAN OUTER BOUNDARY VALUE OR IF Q IS
    !  OUTSIDE THE ACCEPTANCE REGION A
  ELSE
    !ACCEPT
    N=N+1
    Z(N)=V/U
  END IF
END DO  


END SUBROUTINE RAND_NORMAL