source: LMDZ6/trunk/libf/phylmd/cosp2/mo_rng.F90 @ 5069

Last change on this file since 5069 was 3358, checked in by idelkadi, 6 years ago

Implementation de la nouvelle version COSPv2 dans LMDZ.
Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"

File size: 5.4 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2015, Regents of the University of Colorado
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History:
30! May 2015- D. Swales - Original version
31!
32! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33MODULE mod_rng
34  USE cosp_kinds, ONLY: dp, sp, wp
35  IMPLICIT NONE
36
37  INTEGER, parameter :: ki9    = selected_int_kind(R=9)
38  integer :: testInt
39 
40  TYPE rng_state
41     INTEGER(ki9) :: seed ! 32 bit integer
42  END TYPE rng_state
43 
44  INTERFACE init_rng
45     MODULE PROCEDURE init_rng_1, init_rng_n
46  END INTERFACE init_rng
47 
48  INTERFACE get_rng
49     MODULE PROCEDURE get_rng_1, get_rng_n, get_rng_v
50  END INTERFACE get_rng
51 
52CONTAINS
53 
54  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55  ! Set single seed
56  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57  SUBROUTINE init_rng_1(s, seed_in)
58    TYPE(rng_state), INTENT(INOUT) :: s
59    INTEGER,         INTENT(IN   ) :: seed_in
60    s%seed = seed_in     
61  END SUBROUTINE init_rng_1
62 
63  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64  ! Set vector of seeds
65  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66  SUBROUTINE init_rng_n(s, seed_in)
67    TYPE(rng_state), DIMENSION(:), INTENT(INOUT) :: s
68    INTEGER,         DIMENSION(:), INTENT(IN   ) :: seed_in
69   
70    INTEGER :: i
71    DO i = 1, SIZE(seed_in)
72       s(i)%seed = seed_in(i)
73    END DO
74  END SUBROUTINE init_rng_n
75 
76  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77  ! Create single random number from seed
78  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79  FUNCTION get_rng_1(s) 
80    TYPE(rng_state), INTENT(INOUT) :: s
81    REAL(WP)                       :: get_rng_1
82    REAL(SP)                       :: r
83
84    ! Return the next random numbers     
85    ! Marsaglia CONG algorithm
86    s%seed=69069*s%seed+1234567
87    ! mod 32 bit overflow
88    s%seed=mod(s%seed,2_ki9**30_ki9)   
89    r = s%seed*0.931322574615479E-09
90
91    ! convert to range 0-1 (32 bit only)
92    ! DJS2016: What is being done here is an intentional integer overflow and a test to
93    !          see if this occured. Some compilers check for integer overflows during
94    !          compilation (ie. gfortan), while others do not (ie. pgi and ifort). When
95    !          using gfortran, you cannot use the overflow and test for overflow method,
96    !          so we use sizeof(someInt) to determine wheter it is on 32 bit.
97    !if ( i2_16*i2_16 .le. huge32 ) then
98    if (digits(testInt) .le. 31) then
99    !if (sizeof(testInt) .eq. 4) then
100       r=r+1
101       r=r-int(r)
102    endif
103    get_rng_1 = REAL(r, KIND = WP)
104   
105  END FUNCTION get_rng_1
106 
107  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108  ! Create single random number N times
109  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110  FUNCTION get_rng_n(s, n) RESULT (r)
111    integer,intent(inout) :: n
112    TYPE(rng_state),INTENT(INOUT) :: s
113    ! Return the next N random numbers
114    REAL(WP), DIMENSION (n)        :: r
115   
116    INTEGER :: i
117   
118    DO i = 1, N
119       r(i) = get_rng_1(s)
120    END DO
121  END FUNCTION get_rng_n
122 
123  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124  ! Create a vector of random numbers from a vector of input seeds
125  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126  FUNCTION get_rng_v(s) RESULT (r)
127    ! Return the next N random numbers
128    TYPE(rng_state), DIMENSION(:), INTENT(INOUT) :: s
129    REAL(WP),        DIMENSION (SIZE(S))         :: r
130   
131    INTEGER :: i
132   
133    DO i = 1, size(s)
134       r(i) = get_rng_1(s(i))
135    END DO
136  END FUNCTION get_rng_v
137 
138END MODULE mod_rng
Note: See TracBrowser for help on using the repository browser.