scale_random.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_random
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: random_setup
  public :: random_uniform
  public :: random_normal

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: random_reset

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical, private :: random_fix = .false.

  integer, private, allocatable :: random_seedvar(:)

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine random_setup
    use scale_prc, only: &
       prc_abort
    implicit none

    namelist / param_random / &
       random_fix

    integer :: nseeds, ierr
    !---------------------------------------------------------------------------

    log_newline
    log_info("RANDOM_setup",*) 'Setup'

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_random,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("RANDOM_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("RANDOM_setup",*) 'Not appropriate names in namelist PARAM_RANDOM. Check!'
       call prc_abort
    endif
    log_nml(param_random)

    call random_seed
    call random_seed(size=nseeds)

    allocate( random_seedvar(nseeds))

    log_newline
    log_info("RANDOM_setup",*) 'Array size for random seed:', nseeds
    if ( random_fix ) then
       log_info("RANDOM_setup",*) 'random seed is fixed.'
    endif

    call random_reset

    return
  end subroutine random_setup

  !-----------------------------------------------------------------------------
  subroutine random_reset
    use scale_prc, only: &
       prc_myrank
    implicit none

    integer  :: time1(8) ! date and time information
    real(RP) :: time2    ! CPU time
    !---------------------------------------------------------------------------

    if ( random_fix ) then
       ! birthday of SCALE
       time1(1) = 2011 ! The year
       time1(2) = 12   ! The month
       time1(3) = 5    ! The day of the month
       time1(4) = 9    ! Time difference with UTC in minutes
       time1(5) = 10   ! The hour of the day
       time1(6) = 20   ! The minutes of the hour
       time1(7) = 41   ! The seconds of the minute
       time1(8) = 0    ! The milliseconds of the second
       time2    = 0.0_rp
    else
       call date_and_time(values=time1)
       call cpu_time(time2)
    endif

    random_seedvar(:) = prc_myrank                     &
                      + ( time1(1) - 1970 ) * 32140800 &
                      + time1(2) * 2678400             &
                      + time1(3) * 86400               &
                      + time1(5) * 60                  &
                      + time1(6) * 3600                &
                      + time1(7) * 60                  &
                      + time1(8)                       &
                      + int(time2*1.e6_rp)

    call random_seed(put=random_seedvar)

    return
  end subroutine random_reset

  !-----------------------------------------------------------------------------
  subroutine random_uniform( var )
    implicit none

    real(RP), intent(out) :: var(:,:,:)
    !---------------------------------------------------------------------------

    call random_number(var)

    return
  end subroutine random_uniform

  !-----------------------------------------------------------------------------
  subroutine random_normal( var )
    implicit none
    real(RP), intent(out) :: var(:,:,:)
    integer :: n

    n = size(var)
    call get_normal( n, var(:,:,:) )

    return
  end subroutine random_normal

  ! private

  subroutine get_normal( n, var )
    use scale_const, only: &
       pi => const_pi
    implicit none
    integer,  intent(in)  :: n
    real(RP), intent(out) :: var(n)

    real(RP) :: rnd(n+1)
    real(RP) :: fact
    real(RP) :: theta
    integer :: i
    !---------------------------------------------------------------------------

    call random_number(rnd)

    !$omp parallel do &
    !$omp private(fact,theta)
    do i = 1, n/2
       fact = sqrt(-2.0_rp * log( rnd(i*2-1) ) )
       theta = 2.0_rp * pi * rnd(i*2)
       var(i*2-1) = fact * cos(theta)
       var(i*2  ) = fact * sin(theta)
    end do
    if ( mod(n,2) == 1 ) then
       fact = sqrt(-2.0_rp * log( rnd(n) ) )
       theta = 2.0_rp * pi * rnd(n+1)
       var(n) = fact * cos(theta)
    end if

    return
  end subroutine get_normal

end module scale_random