source: LMDZ6/trunk/libf/phylmd/ecrad/random_numbers_mix.F90 @ 4582

!**** *RANDOM_NUMBERS_MIX*  - Portable Random Number Generator

! (C) Copyright 2002- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     Purpose.
!     --------
!           Generate machine-independent pseudo-random numbers

!**   Interface.
!     ----------
!        CALL initialize_random_numbers (kseed, yd_stream)
!        CALL uniform_distribution      (px   , yd_stream)
!        CALL gaussian_distribution     (px   , yd_stream)
!        CALL random_number_restartfile (fname, action)
!        CALL wr_rangen_state           (nunit)

!        Explicit arguments :
!        --------------------
!        kseed  (input)    : integer seed in the range [0,HUGE(kseed)]
!        yd_stream (optional) : the state of the random number generator
!        px     (output)   : array to receive random numbers in the range

!        In the case of uniform_distribution, px has values in the range [0.0,1.0)

!        Implicit arguments :
!        --------------------
!        None

!     Method.
!     -------
!        Based loosly on ZUFALL (Petersen, 1994).

!        The main difference between this generator and ZUFALL is that integer arithmetic
!        is used. This ensures portability to vector machines that implement different
!        real arithmetic. In particular, vector machines often implement non-IEEE
!        arithmetic for their vector pipes. This routine will give identical results for
!        any integer type with at least 32 bits.

!        The generator is a lagged-Fibonacci generator: x(i) = x(i-p) + x(i-q) mod 2**m.
!        Lagged-Fibonacci generators have very long repeat periods: (2**q -1) * 2**(m-1)
!        (i.e about 2.85E191 for q=607, m=30). They pass most tests for randomness.

!        p and q must be chosen carefully. Values from the following table are OK.
!        Larger values give better random numbers, but smaller values are more
!        cache-friendly.

!          q         p
!        9689      4187
!        4423      2098
!        2281      1029
!        1279       418
!         607       273
!         521       168
!         250       103
!         127        63
!          97        33
!          55        24
!          43        22
!          31        13
!          24        10

!        The initial q values of x are set using the binary shirt register method of
!        Burns and Pryor 1999.

!        Mascagni et al (1995) show how to choose different sets of initial values that
!        are guaranteed to be drawn from different maximal-length cycles. This requires
!        the initial values of x(1)...x(q) to be in "canonical form". Specifically,
!        x(1) must be zero and all but a particular one or two values of x must be
!        even. For q=607 and p=273, only one element (jpq-jps) must be odd.

!     Externals.
!     ----------
!        None

!     Reference.
!     ----------
!        Burns P.J. and Pryor D.V. 1999,
!                             Surface Radiative Transport at Large Scale via Monte Carlo.
!                             Annual Review of Heat Transfer, Vol 9.
!
!        Petersen W.P., 1994, Lagged Fibonacci Series Random Number Generator
!                             for the NEC SX-3. International Journal of High Speed Computing
!                             Vol. 6, No. 3, pp387-398.
!
!        Mascagni M., Cuccaro S.A., Pryor D.V., Robinson M.L., 1995,
!                             A Fast, High Quality and Reproducible Parallel Lagged-Fibonacci
!                             Pseudorandom Number Generator. Journal of Computational Physics
!                             Vol 119. pp211-219.

!     Author.
!     -------
!        Mike Fisher *ECMWF*

!     Modifications.
!     --------------
!        Original : 2002-09-25
!        Made parallel friendly: 2003-08-11 Robert Pincus
!        M Leutbecher: 2004-05-10 restart capability
!        M Fisher:     2005-03-30 replaced LCG initialization with shift register
!     ------------------------------------------------------------------

#ifdef RS6K
@PROCESS HOT(NOVECTOR) NOSTRICT
#endif
MODULE RANDOM_NUMBERS_MIX
USE YOMHOOK,  ONLY : LHOOK, DR_HOOK
USE PARKIND1, ONLY : JPIM, JPRB

IMPLICIT NONE

SAVE

PRIVATE
PUBLIC RANDOMNUMBERSTREAM,WR_RANGEN_STATE, &
     & INITIALIZE_RANDOM_NUMBERS, UNIFORM_DISTRIBUTION, GAUSSIAN_DISTRIBUTION ! ,&
!     & RANDOM_NUMBER_RESTARTFILE, 

INTEGER(KIND=JPIM), PARAMETER      :: JPP=273, JPQ=607, JPS=105
INTEGER(KIND=JPIM), PARAMETER      :: JPMM=30
INTEGER(KIND=JPIM), PARAMETER      :: JPM=2**JPMM
INTEGER(KIND=JPIM), PARAMETER      :: JPNUMSPLIT=(JPQ-2)/(JPP-1)
INTEGER(KIND=JPIM), PARAMETER      :: JPLENSPLIT=(JPQ-JPP+JPNUMSPLIT-1)/JPNUMSPLIT
INTEGER(KIND=JPIM), PARAMETER      :: INITVALUE = 12345678

TYPE RANDOMNUMBERSTREAM
  PRIVATE
  INTEGER(KIND=JPIM)                 :: IUSED
  INTEGER(KIND=JPIM)                 :: INITTEST ! Should initialize to zero, but can't in F90
  INTEGER(KIND=JPIM), DIMENSION(JPQ) :: IX 
  REAL(KIND=JPRB)                    :: ZRM
END TYPE RANDOMNUMBERSTREAM

CONTAINS
!-------------------------------------------------------------------------------
SUBROUTINE INITIALIZE_RANDOM_NUMBERS (KSEED, YD_STREAM) 
  !-------------------------------------------------------------------------------
  ! Initialize fibgen
  !-------------------------------------------------------------------------------
  INTEGER(KIND=JPIM),                INTENT(IN   ) :: KSEED
  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
  
  INTEGER, PARAMETER :: JPMASK=123459876
  INTEGER(KIND=JPIM), PARAMETER     :: JPWARMUP_SHFT=64, JPWARMUP_LFG=999
  INTEGER(KIND=JPIM)                :: IDUM,JK,JJ,JBIT
  REAL(KIND=JPRB), DIMENSION(JPWARMUP_LFG)   :: ZWARMUP

  !-------------------------------------------------------------------------------
  ! Initialize the buffer using a binary shift register (Burns and Pryor, 1999).
  ! The Galois representation is used for the shift register as it is more
  ! efficient than the Fibonacci representation. The magic numbers 31 and 87
  ! define the shift register primitive polynomial=(32,7,5,3,2,1,0).
  !
  ! To ensure that different seeds produce distinct initial buffer states in
  ! canonical form, bits 0...jpmm-2 of the initial seed (after XORing with jpmask
  ! and spinning up using the linear congruential generator) are used to construct
  ! x(2), and the remaining bits are used to construct x(jpq).
  !-------------------------------------------------------------------------------
  
  REAL(KIND=JPRB) :: ZHOOK_HANDLE
  IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:INITIALIZE_RANDOM_NUMBERS',0,ZHOOK_HANDLE)
  IDUM = ABS(IEOR(KSEED,JPMASK))
  IF (IDUM==0) IDUM=JPMASK

  DO JJ=1,JPWARMUP_SHFT
    IF (BTEST(IDUM,31)) THEN
      IDUM=IBSET(ISHFT(IEOR(IDUM,87),1),0)
    ELSE
      IDUM=IBCLR(ISHFT(IDUM,1),0)
    ENDIF
  ENDDO

  YD_STREAM%IX(1:JPQ-1)= 0
  YD_STREAM%IX(2)      = ISHFT(IBITS(IDUM,0,JPMM-1),1)
  YD_STREAM%IX(JPQ)    = IBITS(IDUM,JPMM-1,BIT_SIZE(IDUM)+1-JPMM)

  DO JBIT=1,JPMM-1
    DO JJ=3,JPQ-1
      IF (BTEST(IDUM,31)) THEN
        IDUM=IBSET(ISHFT(IEOR(IDUM,87),1),0)
        YD_STREAM%IX(JJ)=IBSET(YD_STREAM%IX(JJ),JBIT)
      ELSE
        IDUM=IBCLR(ISHFT(IDUM,1),0)
      ENDIF
    ENDDO
  ENDDO

  YD_STREAM%IX(JPQ-JPS) = IBSET(YD_STREAM%IX(JPQ-JPS),0)
  
  !-------------------------------------------------------------------------------
  ! Initialize some constants
  !-------------------------------------------------------------------------------
  
  YD_STREAM%IUSED=JPQ
  YD_STREAM%ZRM=1.0_JPRB/REAL(JPM,JPRB)
  
  !-------------------------------------------------------------------------------
  ! Check the calculation of jpnumsplit and jplensplit.
  !-------------------------------------------------------------------------------
  
  IF (JPP+JPNUMSPLIT*JPLENSPLIT < JPQ) THEN
    CALL ABOR1 ('initialize_random_numbers: upper limit of last loop < jpq')
  ENDIF
  
  IF (JPLENSPLIT >=JPP) THEN
    CALL ABOR1 ('initialize_random_numbers: loop length > jpp')
  ENDIF
  
  IF (JPNUMSPLIT>1) THEN
    IF ((JPQ-JPP+JPNUMSPLIT-2)/(JPNUMSPLIT-1) < JPP) THEN
      CALL ABOR1 ('initialize_random_numbers: jpnumsplit is bigger than necessary')
    ENDIF
  ENDIF

  !-------------------------------------------------------------------------------
  ! Set initTest to show that the stream is initialized.
  !-------------------------------------------------------------------------------

  YD_STREAM%INITTEST = INITVALUE
  
  !-------------------------------------------------------------------------------
  ! Warm up the generator.
  !-------------------------------------------------------------------------------

  CALL UNIFORM_DISTRIBUTION (ZWARMUP, YD_STREAM)

IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:INITIALIZE_RANDOM_NUMBERS',1,ZHOOK_HANDLE)
END SUBROUTINE INITIALIZE_RANDOM_NUMBERS

!@PROCESS HOT NOSTRICT
SUBROUTINE UNIFORM_DISTRIBUTION (PX,YD_STREAM)
  !--------------------------------------------------------------------------------
  ! Generate uniformly distributed random numbers in the range 0.0<= px < 1.0
  !--------------------------------------------------------------------------------
  INTEGER(KIND=JPIM), PARAMETER :: IVAR = INT(Z"3FFFFFFF",JPIM)
  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
  REAL(KIND=JPRB), DIMENSION(:),     INTENT(  OUT) :: PX

  INTEGER(KIND=JPIM)                :: JJ, JK, IN, IFILLED
  
  ! This test is a little dirty but Fortran 90 doesn't allow for the initialization
  !   of components of derived types. 
  REAL(KIND=JPRB) :: ZHOOK_HANDLE
! DR_HOOK removed to reduce overhead
! IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',0,ZHOOK_HANDLE)
  IF(YD_STREAM%INITTEST /= INITVALUE) &
    & CALL ABOR1 ('uniform_distribution called before initialize_random_numbers')

  !--------------------------------------------------------------------------------
  ! Copy numbers that were generated during the last call, but not used.
  !--------------------------------------------------------------------------------
  
  IN=SIZE(PX)
  IFILLED=0
  
  DO JJ=YD_STREAM%IUSED+1,MIN(JPQ,IN+YD_STREAM%IUSED)
    PX(JJ-YD_STREAM%IUSED) = YD_STREAM%IX(JJ)*YD_STREAM%ZRM
    IFILLED=IFILLED+1
  ENDDO
  
  YD_STREAM%IUSED=YD_STREAM%IUSED+IFILLED
  
  IF (IFILLED==IN)  THEN 
!   IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE)
! DR_HOOK removed to reduce overhead
    RETURN
  ENDIF
  
  !--------------------------------------------------------------------------------
  ! Generate batches of jpq numbers until px has been filled
  !--------------------------------------------------------------------------------
  
  DO WHILE (IFILLED<IN)
  
  !--------------------------------------------------------------------------------
  ! Generate jpq numbers in vectorizable loops. The first loop is length jpp. The
  ! remaining jpq-jpp elements are calculated in loops of length shorter than jpp.
  !--------------------------------------------------------------------------------
  
  !OCL NOVREC
    DO JJ=1,JPP
!     yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp+jpq)
!     if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm
      YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP+JPQ))
    ENDDO
  
    DO JK=1,JPNUMSPLIT
  !OCL NOVREC
      DO JJ=1+JPP+(JK-1)*JPLENSPLIT,MIN(JPQ,JPP+JK*JPLENSPLIT)
!       yd_stream%ix(jj) = yd_stream%ix(jj) + yd_stream%ix(jj-jpp)
!       if (yd_stream%ix(jj)>=jpm) yd_stream%ix(jj) = yd_stream%ix(jj)-jpm
        YD_STREAM%IX(JJ) = IAND(IVAR,YD_STREAM%IX(JJ) + YD_STREAM%IX(JJ-JPP))
      ENDDO
    ENDDO
  
    YD_STREAM%IUSED = MIN(JPQ,IN-IFILLED)
    PX(IFILLED+1:IFILLED+YD_STREAM%IUSED) = YD_STREAM%IX(1:YD_STREAM%IUSED)*YD_STREAM%ZRM
    IFILLED = IFILLED+YD_STREAM%IUSED
  ENDDO
  
!IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:UNIFORM_DISTRIBUTION',1,ZHOOK_HANDLE)
! DR_HOOK removed to reduce overhead
END SUBROUTINE UNIFORM_DISTRIBUTION
!-------------------------------------------------------------------------------
SUBROUTINE GAUSSIAN_DISTRIBUTION (PX, YD_STREAM)
  TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
  REAL(KIND=JPRB),                   INTENT(  OUT) :: PX(:)
  !--------------------------------------------------------------------------------
  ! Generate normally-distributed random numbers using the Box-Muller method.
  !
  ! NB: this routine does not use buffering. This means that the following calls:
  !     call gaussian_distribution (zx(1:k))
  !     call gaussian_distribution (zx(k+1:n))
  ! will produce different numbers for elements k+1 onwards than the single call:
  !     call gaussian_distribution (zx(1:n))
  !--------------------------------------------------------------------------------
  
  INTEGER(KIND=JPIM) :: ILEN, J
  REAL(KIND=JPRB) :: ZFAC, ZTWOPI
  REAL(KIND=JPRB) :: ZX(SIZE(PX)+1)
  
  !--------------------------------------------------------------------------------
  ! Generate uniform random points in the range [0,1)
  !--------------------------------------------------------------------------------

    REAL(KIND=JPRB) :: ZHOOK_HANDLE
    IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',0,ZHOOK_HANDLE)
    CALL UNIFORM_DISTRIBUTION (ZX, YD_STREAM)

  !--------------------------------------------------------------------------------
  ! Generate gaussian deviates using Box-Muller method
  !--------------------------------------------------------------------------------
  
  ZTWOPI = 8.0_JPRB*ATAN(1.0_JPRB)
  ILEN=SIZE(PX)
  
  DO J=1,ILEN-1,2
    ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(J)))
    PX(J  ) = ZFAC*COS(ZTWOPI*ZX(J+1))
    PX(J+1) = ZFAC*SIN(ZTWOPI*ZX(J+1))
  ENDDO
  
  !--------------------------------------------------------------------------------
  ! Generate the last point if ilen is odd
  !--------------------------------------------------------------------------------
  
  IF (MOD(ILEN,2) /= 0) THEN
    ZFAC = SQRT(-2.0_JPRB*LOG(1.0_JPRB-ZX(ILEN)))
    PX(ILEN) = ZFAC*COS(ZTWOPI*ZX(ILEN+1))
  ENDIF
  
IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:GAUSSIAN_DISTRIBUTION',1,ZHOOK_HANDLE)
END SUBROUTINE GAUSSIAN_DISTRIBUTION
!-------------------------------------------------------------------------------
!!$SUBROUTINE RANDOM_NUMBER_RESTARTFILE( CDFNAME, CDACTION,YD_STREAM )
!!$!--------------------------------------------------------------------------------
!!$!
!!$! read (cdaction='r') or write (cdaction='w') restart file
!!$! for random number generator
!!$!
!!$!--------------------------------------------------------------------------------
!!$CHARACTER (LEN=*),   INTENT(IN) :: CDFNAME
!!$CHARACTER (LEN=1  ), INTENT(IN) :: CDACTION
!!$TYPE(RANDOMNUMBERSTREAM), INTENT(INOUT) :: YD_STREAM
!!$  
!!$INTEGER(KIND=JPIM) :: IUNIT, IRET, IBYTES_IN_JPIM
!!$
!!$REAL(KIND=JPRB) :: ZHOOK_HANDLE
!!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',0,ZHOOK_HANDLE)
!!$IBYTES_IN_JPIM= CEILING(REAL(BIT_SIZE(YD_STREAM%IUSED))/8.0_JPRB - TINY(1.0_JPRB))
!!$
!!$IF (IBYTES_IN_JPIM /= 4) THEN
!!$  CALL ABOR1('random_number_restartfile: number of bytes for JPIM is not 4 ')        
!!$ENDIF
!!$
!!$CALL PBOPEN(IUNIT, CDFNAME, CDACTION, IRET)
!!$IF (IRET /= 0) THEN
!!$  CALL ABOR1('random_number_restartfile: PBOPEN FAILED opening '//CDFNAME)    
!!$ENDIF
!!$
!!$
!!$IF (CDACTION=='r' .OR. CDACTION=='R') THEN
!!$  CALL PBREAD(IUNIT, YD_STREAM%IX,    IBYTES_IN_JPIM*JPQ, IRET)
!!$  IF (IRET < 0) THEN
!!$    CALL ABOR1('random_number_restartfile: PBREAD could not read ix from '//CDFNAME)    
!!$  ENDIF
!!$  CALL PBREAD(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM    , IRET)
!!$  IF (IRET < 0) THEN
!!$    CALL ABOR1('random_number_restartfile: PBREAD could not read iused from '//CDFNAME)    
!!$  ENDIF
!!$
!!$!  l_initialized = .TRUE.
!!$  YD_STREAM%INITTEST = INITVALUE
!!$  YD_STREAM%ZRM=1.0_JPRB/REAL(JPM,JPRB)
!!$ELSEIF(CDACTION=='w' .OR. CDACTION=='W') THEN
!!$  CALL PBWRITE(IUNIT, YD_STREAM%IX, IBYTES_IN_JPIM*JPQ, IRET)
!!$  IF (IRET < 0) THEN
!!$    CALL ABOR1('random_number_restartfile: PBWRITE could not write ix on '//CDFNAME)    
!!$  ENDIF
!!$  CALL PBWRITE(IUNIT, YD_STREAM%IUSED, IBYTES_IN_JPIM , IRET)
!!$  IF (IRET < 0) THEN
!!$    CALL ABOR1('random_number_restartfile: PBWRITE could not write iused on '//CDFNAME)    
!!$  ENDIF
!!$
!!$ELSE
!!$  CALL ABOR1 ('random_number_restartfile: cdaction = '//CDACTION//' is undefined.')
!!$ENDIF
!!$
!!$CALL PBCLOSE(IUNIT, IRET)
!!$IF (IRET /= 0) THEN
!!$  CALL ABOR1('random_number_restartfile: PBCLOSE FAILED closing '//CDFNAME)    
!!$ENDIF
!!$
!!$IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:RANDOM_NUMBER_RESTARTFILE',1,ZHOOK_HANDLE)
!!$END SUBROUTINE RANDOM_NUMBER_RESTARTFILE


SUBROUTINE WR_RANGEN_STATE( KUNIT, YD_STREAM )
!--------------------------------------------------------------------------------
! write state of random number generator to unit kunit
!--------------------------------------------------------------------------------
INTEGER(KIND=JPIM), INTENT(IN) :: KUNIT
TYPE(RANDOMNUMBERSTREAM), INTENT(IN) :: YD_STREAM

REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',0,ZHOOK_HANDLE)
WRITE( KUNIT, * ) 'module random_numbers_mix, generator state is'
WRITE( KUNIT, '(8I10)') YD_STREAM%IX
WRITE( KUNIT, '(I10)')  YD_STREAM%IUSED

IF (LHOOK) CALL DR_HOOK('RANDOM_NUMBERS_MIX:WR_RANGEN_STATE',1,ZHOOK_HANDLE)
END SUBROUTINE WR_RANGEN_STATE

END MODULE RANDOM_NUMBERS_MIX