random_numbers.F90

Go to the documentation of this file.

!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************

module random_numbers_mod 
  ! Generic module to wrap random number generators. 
  !   The module defines a type that identifies the particular stream of random 
  !   numbers, and has procedures for initializing it and getting real numbers 
  !   in the range 0 to 1. 
  ! This version uses the Mersenne Twister to generate random numbers on [0, 1]. 
  !
  use mersennetwister_mod, only: randomnumbersequence, & ! The random number engine.
                                 new_randomnumbersequence, getrandomreal 
  use time_manager_mod, only: time_type, get_date
  implicit none
  private
  
  type randomnumberstream
    type(randomnumbersequence) :: thenumbers
  end type randomnumberstream
  
  interface getrandomnumbers
    module procedure getrandomnumber_scalar, getrandomnumber_1d, getrandomnumber_2d
  end interface getrandomnumbers
  
  interface initializerandomnumberstream
    module procedure initializerandomnumberstream_s, initializerandomnumberstream_v
  end interface initializerandomnumberstream

  public :: randomnumberstream,                             &
            initializerandomnumberstream, getrandomnumbers, &
            constructseed
contains
  ! ---------------------------------------------------------
  ! Initialization
  ! ---------------------------------------------------------
  function initializerandomnumberstream_s(seed) result(new) 
    integer, intent( in)     :: seed
    type(randomnumberstream) :: new
    
    new%theNumbers = new_randomnumbersequence(seed)
    
  end function initializerandomnumberstream_s
  ! ---------------------------------------------------------
  function initializerandomnumberstream_v(seed) result(new) 
    integer, dimension(:), intent( in) :: seed
    type(randomnumberstream)           :: new
    
    new%theNumbers = new_randomnumbersequence(seed)
    
  end function initializerandomnumberstream_v
  ! ---------------------------------------------------------
  ! Procedures for drawing random numbers
  ! ---------------------------------------------------------
  subroutine getrandomnumber_scalar(stream, number)
    type(randomNumberStream), intent(inout) :: stream
    real,                     intent(  out) :: number
    
    number = getrandomreal(stream%theNumbers)
  end subroutine getrandomnumber_scalar
  ! ---------------------------------------------------------
  subroutine getrandomnumber_1d(stream, numbers)
    type(randomNumberStream), intent(inout) :: stream
    real, dimension(:),       intent(  out) :: numbers
    
    ! Local variables
    integer :: i
    
    do i = 1, size(numbers)
      numbers(i) = getrandomreal(stream%theNumbers)
    end do
  end subroutine getrandomnumber_1d
  ! ---------------------------------------------------------
  subroutine getrandomnumber_2d(stream, numbers)
    type(randomNumberStream), intent(inout) :: stream
    real, dimension(:, :),    intent(  out) :: numbers
    
    ! Local variables
    integer :: i
    
    do i = 1, size(numbers, 2)
      call getrandomnumber_1d(stream, numbers(:, i))
    end do
  end subroutine getrandomnumber_2d
  ! ---------------------------------------------------------
  ! Constructs a unique seed from grid cell index and model date/time
  !   The perm is supplied we generate a different seed by 
  !   circularly shifting the bits of the seed - this is useful 
  !   if we want to create more than one seed for a given 
  !   column and model date/time. 
  !   Note that abs(perm) must be <= the number of bits used 
  !   to represent the default integer (likely 32) 
  ! ---------------------------------------------------------
  function constructseed(i, j, time, perm) result(seed)
    integer,           intent( in)  :: i, j
    type(time_type),   intent( in) :: time
    integer, optional, intent( in) :: perm
    integer, dimension(8) :: seed
    
    ! Local variables
    integer :: year, month, day, hour, minute, second
    
    
    call get_date(time, year, month, day, hour, minute, second)
    seed = (/ i, j, year, month, day, hour, minute, second /)
    if(present(perm)) seed = ishftc(seed, perm)
  end function constructseed
end module random_numbers_mod