Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Jul 8, 2020 9:56:29 AM (4 years ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/random_generator_parallel_mod.f90

    r4545 r4592  
    11!> @file random_generator_parallel_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
    27 ! Add generator (using parallel mode) returning gaussian distributed random
    28 ! number
    29 !
     27! File re-formatted to follow the PALM coding standard
     28!
     29! 4545 2020-05-22 13:17:57Z schwenkel
     30! Add generator (using parallel mode) returning gaussian distributed random number
     31!
    3032! 4438 2020-03-03 20:49:28Z raasch
    3133! - Rename variables to avoid confusion with subdomain grid indices
    3234! - Some formatting adjustments
    33 ! - New routine to initialize spatial 1D-arrays independent on the 2D random-
    34 !   number array
    35 !
     35! - New routine to initialize spatial 1D-arrays independent on the 2D random-number array
     36!
    3637! 4360 2020-01-07 11:25:50Z suehring
    3738! Corrected "Former revisions" section
    38 ! 
     39!
    3940! 3655 2019-01-07 16:51:22Z knoop
    40 ! unused variable removed
     41! Unused variable removed
    4142!
    4243! 1400 2014-05-09 14:03:54Z knoop
    4344! Initial revision
    44 !
    45 !
     45!
     46!
     47!--------------------------------------------------------------------------------------------------!
    4648! Description:
    4749! ------------
    4850!> This module contains and supports the random number generating routine ran_parallel.
    49 !> ran_parallel returns a uniform random deviate between 0.0 and 1.0
    50 !> (exclusive of the end point values).
    51 !> Moreover, it contains a routine returning a normally distributed random number
    52 !> with mean zero and unity standard deviation.
    53 !> Additionally it provides the generator with five integer for use as initial state space.
    54 !> The first tree integers (iran, jran, kran) are maintained as non negative values,
    55 !> while the last two (mran, nran) have 32-bit nonzero values.
    56 !> Also provided by this module is support for initializing or reinitializing
    57 !> the state space to a desired standard sequence number, hashing the initial
    58 !> values to random values, and allocating and deallocating the internal workspace
    59 !> Random number generator, produces numbers equally distributed in interval
     51!> ran_parallel returns a uniform random deviate between 0.0 and 1.0 (exclusive of the end point
     52!> values). Moreover, it contains a routine returning a normally distributed random number with mean
     53!> zero and unity standard deviation. Additionally, it provides the generator with five integers for
     54!> use as initial state space. The first tree integers (iran, jran, kran) are maintained as non
     55!> negative values, while the last two (mran, nran) have 32-bit nonzero values. Also provided by
     56!> this module is support for initializing or reinitializing the state space to a desired standard
     57!> sequence number, hashing the initial values to random values and allocating and deallocating the
     58!> internal workspace.
    6059!>
    6160!> This routine is taken from the "numerical recipies vol. 2"
    62 !------------------------------------------------------------------------------!
    63 MODULE random_generator_parallel
    64 
    65    USE kinds
    66 
    67    IMPLICIT NONE
    68 
    69    INTEGER(isp), SAVE :: lenran=0             !<
    70    INTEGER(isp), SAVE :: seq=0                !<
    71    INTEGER(isp), SAVE :: iran0                !<
    72    INTEGER(isp), SAVE :: jran0                !<
    73    INTEGER(isp), SAVE :: kran0                !<
    74    INTEGER(isp), SAVE :: mran0                !<
    75    INTEGER(isp), SAVE :: nran0                !<
    76    INTEGER(isp), SAVE :: rans                 !<
    77 
    78    INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds   !<
    79 
    80    INTEGER(isp), DIMENSION(:), POINTER, SAVE :: iran   !<
    81    INTEGER(isp), DIMENSION(:), POINTER, SAVE :: jran   !<
    82    INTEGER(isp), DIMENSION(:), POINTER, SAVE :: kran   !<
    83    INTEGER(isp), DIMENSION(:), POINTER, SAVE :: mran   !<
    84    INTEGER(isp), DIMENSION(:), POINTER, SAVE :: nran   !<
    85    INTEGER(isp), DIMENSION(:), POINTER, SAVE :: ranv   !<
    86 
    87    INTEGER(isp), DIMENSION(:,:), ALLOCATABLE   ::  id_random_array    !<
    88    INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE ::  seq_random_array   !<
    89 
    90    REAL(wp), SAVE :: amm   !<
    91 
    92    REAL(wp) :: random_dummy=0.0   !<
    93 
    94    INTERFACE init_parallel_random_generator
    95       MODULE PROCEDURE init_parallel_random_generator_1d
    96       MODULE PROCEDURE init_parallel_random_generator_2d
    97    END INTERFACE
    98 
    99    INTERFACE random_number_parallel
    100       MODULE PROCEDURE ran0_s
    101    END INTERFACE
    102 
    103    INTERFACE random_number_parallel_gauss
    104       MODULE PROCEDURE gasdev_s
    105    END INTERFACE
    106 
    107    INTERFACE random_seed_parallel
    108       MODULE PROCEDURE random_seed_parallel
    109    END INTERFACE
    110 
    111    INTERFACE ran_hash
    112       MODULE PROCEDURE ran_hash_v
    113    END INTERFACE
    114 
    115    INTERFACE reallocate
    116       MODULE PROCEDURE reallocate_iv,reallocate_im
    117    END INTERFACE
    118 
    119    INTERFACE arth
    120       MODULE PROCEDURE arth_i
    121    END INTERFACE
    122 
    123    PRIVATE
    124 
    125    PUBLIC random_number_parallel, random_seed_parallel, random_dummy,          &
    126           id_random_array, seq_random_array, init_parallel_random_generator,   &
    127           random_number_parallel_gauss
     61!--------------------------------------------------------------------------------------------------!
     62 MODULE random_generator_parallel
     63
     64    USE kinds
     65
     66    IMPLICIT NONE
     67
     68    INTEGER(isp), SAVE ::  lenran = 0  !<
     69    INTEGER(isp), SAVE ::  seq = 0     !<
     70    INTEGER(isp), SAVE ::  iran0       !<
     71    INTEGER(isp), SAVE ::  jran0       !<
     72    INTEGER(isp), SAVE ::  kran0       !<
     73    INTEGER(isp), SAVE ::  mran0       !<
     74    INTEGER(isp), SAVE ::  nran0       !<
     75    INTEGER(isp), SAVE ::  rans        !<
     76
     77    INTEGER(isp), DIMENSION(:, :), POINTER, SAVE ::  ranseeds  !<
     78
     79    INTEGER(isp), DIMENSION(:), POINTER, SAVE ::  iran  !<
     80    INTEGER(isp), DIMENSION(:), POINTER, SAVE ::  jran  !<
     81    INTEGER(isp), DIMENSION(:), POINTER, SAVE ::  kran  !<
     82    INTEGER(isp), DIMENSION(:), POINTER, SAVE ::  mran  !<
     83    INTEGER(isp), DIMENSION(:), POINTER, SAVE ::  nran  !<
     84    INTEGER(isp), DIMENSION(:), POINTER, SAVE ::  ranv  !<
     85
     86    INTEGER(isp), DIMENSION(:,:), ALLOCATABLE   ::  id_random_array   !<
     87    INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE ::  seq_random_array  !<
     88
     89    REAL(wp), SAVE ::  amm  !<
     90
     91    REAL(wp) ::  random_dummy = 0.0  !<
     92
     93    INTERFACE init_parallel_random_generator
     94       MODULE PROCEDURE init_parallel_random_generator_1d
     95       MODULE PROCEDURE init_parallel_random_generator_2d
     96    END INTERFACE
     97
     98    INTERFACE random_number_parallel
     99       MODULE PROCEDURE ran0_s
     100    END INTERFACE
     101
     102    INTERFACE random_number_parallel_gauss
     103       MODULE PROCEDURE gasdev_s
     104    END INTERFACE
     105
     106    INTERFACE random_seed_parallel
     107       MODULE PROCEDURE random_seed_parallel
     108    END INTERFACE
     109
     110    INTERFACE ran_hash
     111       MODULE PROCEDURE ran_hash_v
     112    END INTERFACE
     113
     114    INTERFACE reallocate
     115       MODULE PROCEDURE reallocate_iv,reallocate_im
     116    END INTERFACE
     117
     118    INTERFACE arth
     119       MODULE PROCEDURE arth_i
     120    END INTERFACE
     121
     122    PRIVATE
     123
     124    PUBLIC random_number_parallel,                                                                 &
     125           random_seed_parallel,                                                                   &
     126           random_dummy,                                                                           &
     127           id_random_array,                                                                        &
     128           seq_random_array,                                                                       &
     129           init_parallel_random_generator,                                                         &
     130           random_number_parallel_gauss
    128131
    129132 CONTAINS
    130133
    131 !------------------------------------------------------------------------------!
     134!--------------------------------------------------------------------------------------------------!
    132135! Description:
    133136! ------------
    134137!> Initialize the parallel random number generator for a 1-dimensional array.
    135 !------------------------------------------------------------------------------!
    136    SUBROUTINE init_parallel_random_generator_1d( nxy, ns, ne, id_rand, seq_rand )
    137 
    138       USE control_parameters,                                                  &
    139           ONLY: ensemble_member_nr
    140 
    141       INTEGER(isp), INTENT(IN) ::  nxy  !< constant scaling with grid dimensions
    142       INTEGER(isp), INTENT(IN) ::  ns   !< start index on subdomain
    143       INTEGER(isp), INTENT(IN) ::  ne   !< end index on subdomain
    144 
    145       INTEGER(iwp) ::  j   !< loop index
    146 
    147       INTEGER(isp), DIMENSION(ns:ne)   ::  id_rand  !< initial IDs
    148       INTEGER(isp), DIMENSION(5,ns:ne) ::  seq_rand !< initial random seeds
    149 
    150 !
    151 !--   Asigning an ID to every vertical gridpoint column
    152 !--   dependig on the ensemble run number.
    153       DO  j = ns, ne
    154          id_rand(j) = j * ( nxy + 1.0_wp ) + 1.0_wp + 1E6 *                    &
    155                                 ( ensemble_member_nr - 1000 )
    156       ENDDO
    157 !
    158 !--   Initializing with random_seed_parallel for every vertical
    159 !--   gridpoint column.
    160       random_dummy=0
    161       DO  j = ns, ne
    162          CALL random_seed_parallel( random_sequence=id_rand(j) )
    163          CALL random_number_parallel( random_dummy )
    164          CALL random_seed_parallel( get=seq_rand(:,j) )
    165       ENDDO
    166 
    167    END SUBROUTINE init_parallel_random_generator_1d
    168 
    169 !------------------------------------------------------------------------------!
    170 ! Description:
    171 ! ------------
    172 !> Initialize the parallel random number generator for a specific subdomain
    173 !------------------------------------------------------------------------------!
    174    SUBROUTINE init_parallel_random_generator_2d( nx_l, nys_l, nyn_l, nxl_l, nxr_l )
    175 
    176       USE kinds
    177 
    178       USE control_parameters,                                                  &
    179           ONLY: ensemble_member_nr
    180 
    181       IMPLICIT NONE
    182 
    183       INTEGER(isp), INTENT(IN) ::  nx_l    !< constant
    184       INTEGER(isp), INTENT(IN) ::  nys_l   !< local lower subdomain bound index in y-direction
    185       INTEGER(isp), INTENT(IN) ::  nyn_l   !< local upper subdomain bound index in y-direction
    186       INTEGER(isp), INTENT(IN) ::  nxl_l   !< local lower subdomain bound index in x-direction
    187       INTEGER(isp), INTENT(IN) ::  nxr_l   !< local upper subdomain bound index in x-direction
    188 
    189       INTEGER(iwp) ::  i   !< grid index x-direction
    190       INTEGER(iwp) ::  j   !< grid index y-direction
    191 
    192 !--   Allocate ID-array and state-space-array
    193       ALLOCATE ( seq_random_array(5,nys_l:nyn_l,nxl_l:nxr_l) )
    194       ALLOCATE ( id_random_array(nys_l:nyn_l,nxl_l:nxr_l) )
    195       seq_random_array = 0
    196       id_random_array  = 0
    197 
    198 !--   Asigning an ID to every vertical gridpoint column
    199 !--   dependig on the ensemble run number.
    200       DO  i = nxl_l, nxr_l
    201          DO  j = nys_l, nyn_l
    202             id_random_array(j,i) = j * ( nx_l + 1.0_wp ) + i + 1.0_wp + 1E6 *  &
    203                                    ( ensemble_member_nr - 1000 )
    204          ENDDO
    205       ENDDO
    206 !--   Initializing with random_seed_parallel for every vertical
    207 !--   gridpoint column.
    208       random_dummy=0
    209       DO  i = nxl_l, nxr_l
    210          DO  j = nys_l, nyn_l
    211             CALL random_seed_parallel (random_sequence=id_random_array(j, i))
    212             CALL random_number_parallel (random_dummy)
    213             CALL random_seed_parallel (get=seq_random_array(:, j, i))
    214          ENDDO
    215       ENDDO
    216 
    217    END SUBROUTINE init_parallel_random_generator_2d
    218  
    219 !------------------------------------------------------------------------------!
    220 ! Description:
    221 ! ------------
    222 !> Lagged Fibonacci generator combined with a Marsaglia shiftsequence.
    223 !> Returns as harvest a uniform random deviate between 0.0 and 1.0 (exclusive of the end point values).
    224 !> This generator has the same calling and initialization conventions as Fortran 90's random_number routine.
    225 !> Use random_seed_parallel to initialize or reinitialize to a particular sequence.
     138!--------------------------------------------------------------------------------------------------!
     139 SUBROUTINE init_parallel_random_generator_1d( nxy, ns, ne, id_rand, seq_rand )
     140
     141    USE control_parameters,                                                                        &
     142        ONLY: ensemble_member_nr
     143
     144    INTEGER(isp), INTENT(IN) ::  nxy  !< constant scaling with grid dimensions
     145    INTEGER(isp), INTENT(IN) ::  ns   !< start index on subdomain
     146    INTEGER(isp), INTENT(IN) ::  ne   !< end index on subdomain
     147
     148    INTEGER(iwp) ::  j  !< loop index
     149
     150    INTEGER(isp), DIMENSION(ns:ne)   ::  id_rand   !< initial IDs
     151    INTEGER(isp), DIMENSION(5,ns:ne) ::  seq_rand  !< initial random seeds
     152
     153!
     154!-- Asigning an ID to every vertical gridpoint column dependig on the ensemble run number.
     155    DO  j = ns, ne
     156       id_rand(j) = j * ( nxy + 1.0_wp ) + 1.0_wp + 1E6 * ( ensemble_member_nr - 1000 )
     157    ENDDO
     158!
     159!-- Initializing with random_seed_parallel for every vertical gridpoint column.
     160    random_dummy = 0
     161    DO  j = ns, ne
     162       CALL random_seed_parallel( random_sequence=id_rand(j) )
     163       CALL random_number_parallel( random_dummy )
     164       CALL random_seed_parallel( get=seq_rand(:,j) )
     165    ENDDO
     166
     167 END SUBROUTINE init_parallel_random_generator_1d
     168
     169!--------------------------------------------------------------------------------------------------!
     170! Description:
     171! ------------
     172!> Initialize the parallel random number generator for a specific subdomain.
     173!--------------------------------------------------------------------------------------------------!
     174 SUBROUTINE init_parallel_random_generator_2d( nx_l, nys_l, nyn_l, nxl_l, nxr_l )
     175
     176    USE kinds
     177
     178    USE control_parameters,                                                                        &
     179        ONLY: ensemble_member_nr
     180
     181    IMPLICIT NONE
     182
     183    INTEGER(isp), INTENT(IN) ::  nx_l   !< constant
     184    INTEGER(isp), INTENT(IN) ::  nys_l  !< local lower subdomain bound index in y-direction
     185    INTEGER(isp), INTENT(IN) ::  nyn_l  !< local upper subdomain bound index in y-direction
     186    INTEGER(isp), INTENT(IN) ::  nxl_l  !< local lower subdomain bound index in x-direction
     187    INTEGER(isp), INTENT(IN) ::  nxr_l  !< local upper subdomain bound index in x-direction
     188
     189    INTEGER(iwp) ::  i  !< grid index x-direction
     190    INTEGER(iwp) ::  j  !< grid index y-direction
     191
     192!-- Allocate ID-array and state-space-array
     193    ALLOCATE ( seq_random_array(5,nys_l:nyn_l,nxl_l:nxr_l) )
     194    ALLOCATE ( id_random_array(nys_l:nyn_l,nxl_l:nxr_l) )
     195    seq_random_array = 0
     196    id_random_array  = 0
     197
     198!-- Asigning an ID to every vertical gridpoint column dependig on the ensemble run number.
     199    DO  i = nxl_l, nxr_l
     200       DO  j = nys_l, nyn_l
     201          id_random_array(j,i) = j * ( nx_l + 1.0_wp ) + i + 1.0_wp + 1E6 *                        &
     202                                 ( ensemble_member_nr - 1000 )
     203       ENDDO
     204    ENDDO
     205!-- Initializing with random_seed_parallel for every vertical gridpoint column.
     206    random_dummy = 0
     207    DO  i = nxl_l, nxr_l
     208       DO  j = nys_l, nyn_l
     209          CALL random_seed_parallel( random_sequence = id_random_array(j,i) )
     210          CALL random_number_parallel( random_dummy )
     211          CALL random_seed_parallel( get = seq_random_array(:,j,i) )
     212       ENDDO
     213    ENDDO
     214
     215 END SUBROUTINE init_parallel_random_generator_2d
     216
     217!--------------------------------------------------------------------------------------------------!
     218! Description:
     219! ------------
     220!> Lagged Fibonacci generator combined with a Marsaglia shift sequence. Returns as harvest a uniform
     221!> random deviate between 0.0 and 1.0 (exclusive of the end point values). This generator has the
     222!> same calling and initialization conventions as Fortran 90's random_number routine.
     223!> Use random_seed_parallel to initialize or reinitialize to a particular sequence.
    226224!> The period of this generator is about 2.0 x 10^28, and it fully vectorizes.
    227225!> Validity of the integer model assumed by this generator is tested at initialization.
    228 !------------------------------------------------------------------------------!
    229    SUBROUTINE ran0_s(harvest)
    230 
    231       USE kinds
    232      
    233       IMPLICIT NONE
    234      
    235       REAL(wp), INTENT(OUT) :: harvest   !<
    236      
    237       IF  (lenran < 1) CALL ran_init(1)  !- Initialization routine in ran_state.
    238      
    239       !- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69.
    240       rans = iran0 - kran0   
    241      
    242       IF  (rans < 0) rans = rans + 2147483579_isp
    243      
    244       iran0 = jran0
    245       jran0 = kran0
    246       kran0 = rans
    247      
    248       nran0 = ieor( nran0, ishft (nran0, 13) ) !- Update Marsaglia shift sequence with period 2^32 - 1.
    249       nran0 = ieor( nran0, ishft (nran0, -17) )
    250       nran0 = ieor( nran0, ishft (nran0, 5) )
    251      
    252       rans  = ieor( nran0, rans )   !- Combine the generators.
    253      
    254       harvest = amm * merge( rans, not(rans), rans < 0 ) !- Make the result positive definite (note that amm is negative).
    255      
    256    END SUBROUTINE ran0_s
    257 
    258 !------------------------------------------------------------------------------!
    259 ! Description:
    260 ! ------------
    261 !> Returns in harvest a normally distributed deviate with zero mean and unit
    262 !> variance, using ran0_s as the source of uniform deviates. Following
    263 !> Numerical Recipes in Fortran90 (Press et al., 2nd edition, 1996, pp 1152ff).
    264 !> Note, instead of ran1_s, ran0_s is used.
    265 !------------------------------------------------------------------------------!
    266    SUBROUTINE gasdev_s(harvest)
    267 
    268       REAL(wp), INTENT(OUT) ::  harvest   !<
    269 
    270       REAL(wp) ::  rsq      !<
    271       REAL(wp) ::  v1       !<
    272       REAL(wp) ::  v2       !<
    273       REAL(wp), SAVE ::  g  !<
    274 
    275       LOGICAL, SAVE ::  gaus_stored = .FALSE. !<
    276 !
    277 !--   We have an extra deviate handy, so return it, and unset the flag.
    278       IF (gaus_stored)  THEN
    279          harvest = g
    280          gaus_stored = .FALSE.
    281 !
    282 !--   We don’t have an extra deviate handy, so pick two uniform numbers in the
    283 !--   square extending from -1 to +1 in each direction
    284       ELSE
    285          DO
    286             CALL ran0_s(v1)
    287             CALL ran0_s(v2)
    288             v1 = 2.0_wp * v1 - 1.0_wp
    289             v2 = 2.0_wp * v2 - 1.0_wp
    290 !
    291 !--         see if they are in the unit circle
    292             rsq = v1**2 + v2**2
    293 !
    294 !--         otherwise try again.
    295             IF (rsq > 0.0  .AND.  rsq < 1.0) EXIT
    296          ENDDO
    297 !
    298 !--      Now make the Box-Muller transformation to get two normal deviates.
    299 !--      Return one and save the other for next time. Set flag.
    300          rsq = SQRT(-2.0_sp * LOG(rsq)/rsq)
    301          harvest = v1 * rsq
    302          g = v2 * rsq
    303          gaus_stored = .TRUE.
    304       ENDIF
    305 
    306    END SUBROUTINE gasdev_s
    307 
    308 !------------------------------------------------------------------------------!
    309 ! Description:
    310 ! ------------
    311 !> Initialize or reinitialize the random generator state space to vectors of size length.
    312 !> The saved variable seq is hashed (via calls to the module routine ran_hash)
    313 !> to create unique starting seeds, different for each vector component.
    314 !------------------------------------------------------------------------------!
    315    SUBROUTINE ran_init( length )
    316    
    317       USE kinds
    318      
    319       IMPLICIT NONE
    320      
    321       INTEGER(isp), INTENT(IN) ::  length   !<
    322    
    323       INTEGER(isp) ::  hg    !<
    324       INTEGER(isp) ::  hgm   !<
    325       INTEGER(isp) ::  hgng  !<
    326      
    327       INTEGER(isp) ::  new   !<
    328       INTEGER(isp) ::  j     !<
    329       INTEGER(isp) ::  hgt   !<
    330      
    331       IF ( length < lenran ) RETURN !- Simply return if enough space is already allocated.
    332      
    333       hg=huge(1_isp)
    334       hgm=-hg
    335       hgng=hgm-1
    336       hgt = hg
    337      
    338       !- The following lines check that kind value isp is in fact a 32-bit integer
    339       !- with the usual properties that we expect it to have (under negation and wrap-around addition).
    340       !- If all of these tests are satisfied, then the routines that use this module are portable,
    341       !- even though they go beyond Fortran 90's integer model.
    342      
    343       IF  ( hg /= 2147483647 ) CALL ran_error('arithmetic assumption 1 failed')
    344       IF  ( hgng >= 0 )        CALL ran_error('arithmetic assumption 2 failed')
    345       IF  ( hgt+1 /= hgng )    CALL ran_error('arithmetic assumption 3 failed')
    346       IF  ( not(hg) >= 0 )     CALL ran_error('arithmetic assumption 4 failed')
    347       IF  ( not(hgng) < 0 )    CALL ran_error('arithmetic assumption 5 failed')
    348       IF  ( hg+hgng >= 0 )     CALL ran_error('arithmetic assumption 6 failed')
    349       IF  ( not(-1_isp) < 0 )  CALL ran_error('arithmetic assumption 7 failed')
    350       IF  ( not(0_isp) >= 0 )  CALL ran_error('arithmetic assumption 8 failed')
    351       IF  ( not(1_isp) >= 0 )  CALL ran_error('arithmetic assumption 9 failed')
    352      
    353       IF  ( lenran > 0) THEN                          !- Reallocate space, or ...
    354      
    355          ranseeds => reallocate( ranseeds, length, 5)
    356          ranv => reallocate( ranv, length-1)
    357          new = lenran+1
    358          
    359       ELSE                                            !- allocate space.
    360      
    361          ALLOCATE(ranseeds(length,5))
    362          ALLOCATE(ranv(length-1))
    363          new = 1   !- Index of first location not yet initialized.
    364          amm = nearest(1.0_wp,-1.0_wp)/hgng
    365          !- Use of nearest is to ensure that returned random deviates are strictly lessthan 1.0.
    366          IF  (amm*hgng >= 1.0 .or. amm*hgng <= 0.0)                            &
    367             CALL ran_error('arithmetic assumption 10 failed')
    368            
    369       END IF
    370      
    371       !- Set starting values, unique by seq and vector component.
    372       ranseeds(new:,1) = seq
    373       ranseeds(new:,2:5)=spread(arth(new,1,size(ranseeds(new:,1))),2,4)
    374      
    375       DO j=1,4   !- Hash them.
    376          CALL ran_hash(ranseeds(new:,j), ranseeds(new:,j+1))
    377       END DO
    378      
    379       WHERE (ranseeds (new: ,1:3) < 0)                                         &
    380          ranseeds(new: ,1:3)=not(ranseeds(new: ,1:3))  !- Enforce nonnegativity.
    381          
    382       WHERE (ranseeds(new: ,4:5) == 0) ranseeds(new: ,4:5)=1 !- Enforce nonzero.
    383      
    384       IF  (new == 1) THEN !- Set scalar seeds.
    385      
    386          iran0 = ranseeds(1,1)
    387          jran0 = ranseeds(1,2)
    388          kran0 = ranseeds(1,3)
    389          mran0 = ranseeds(1,4)
    390          nran0 = ranseeds(1,5)
    391          rans  = nran0
    392          
    393       END IF
    394      
    395       IF  (length > 1) THEN   !- Point to vector seeds.
    396      
    397          iran => ranseeds(2:,1)
    398          jran => ranseeds(2:,2)
    399          kran => ranseeds(2:,3)
    400          mran => ranseeds(2:,4)
    401          nran => ranseeds(2:,5)
    402          ranv = nran
    403          
    404       END IF
    405      
    406       lenran = length
    407      
    408    END SUBROUTINE ran_init
    409 
    410 !------------------------------------------------------------------------------!
     226!--------------------------------------------------------------------------------------------------!
     227 SUBROUTINE ran0_s( harvest )
     228
     229    USE kinds
     230
     231    IMPLICIT NONE
     232
     233    REAL(wp), INTENT(OUT) ::  harvest  !<
     234!
     235!-- Initialization routine in ran_state.
     236    IF ( lenran < 1 )  CALL ran_init( 1 )
     237!
     238!-- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69.
     239    rans = iran0 - kran0
     240
     241    IF ( rans < 0 )  rans = rans + 2147483579_isp
     242
     243    iran0 = jran0
     244    jran0 = kran0
     245    kran0 = rans
     246!
     247!-- Update Marsaglia shift sequence with period 2^32 - 1.
     248    nran0 = IEOR( nran0, ISHFT( nran0, 13 ) )
     249    nran0 = IEOR( nran0, ISHFT( nran0, -17 ) )
     250    nran0 = IEOR( nran0, ISHFT( nran0, 5 ) )
     251!
     252!-- Combine the generators.
     253    rans  = IEOR( nran0, rans )
     254!
     255!-- Make the result positive definite (note that amm is negative).
     256    harvest = amm * MERGE( rans, NOT( rans ), rans < 0 )
     257
     258 END SUBROUTINE ran0_s
     259
     260!--------------------------------------------------------------------------------------------------!
     261! Description:
     262! ------------
     263!> Returns in harvest a normally distributed deviate with zero mean and unit variance, using ran0_s
     264!> as the source of uniform deviates. Following Numerical Recipes in Fortran90 (Press et al., 2nd
     265!> edition, 1996, pp 1152ff). Note, instead of ran1_s, ran0_s is used.
     266!--------------------------------------------------------------------------------------------------!
     267 SUBROUTINE gasdev_s( harvest )
     268
     269    REAL(wp), INTENT(OUT) ::  harvest  !<
     270
     271    REAL(wp) ::  rsq  !<
     272    REAL(wp) ::  v1   !<
     273    REAL(wp) ::  v2   !<
     274    REAL(wp), SAVE ::  g  !<
     275
     276    LOGICAL, SAVE ::  gaus_stored = .FALSE.  !<
     277!
     278!-- We have an extra deviate handy, so return it, and unset the flag.
     279    IF ( gaus_stored )  THEN
     280       harvest = g
     281       gaus_stored = .FALSE.
     282!
     283!-- We don’t have an extra deviate handy, so pick two uniform numbers in the square extending
     284!-- from -1 to +1 in each direction
     285    ELSE
     286       DO
     287          CALL ran0_s( v1 )
     288          CALL ran0_s( v2 )
     289          v1 = 2.0_wp * v1 - 1.0_wp
     290          v2 = 2.0_wp * v2 - 1.0_wp
     291!
     292!--       See if they are in the unit circle
     293          rsq = v1**2 + v2**2
     294!
     295!--       Otherwise try again.
     296          IF ( rsq > 0.0  .AND.  rsq < 1.0 )  EXIT
     297       ENDDO
     298!
     299!--    Now make the Box-Muller transformation to get two normal deviates.
     300!--    Return one and save the other for next time. Set flag.
     301       rsq = SQRT( - 2.0_sp * LOG( rsq ) / rsq )
     302       harvest = v1 * rsq
     303       g = v2 * rsq
     304       gaus_stored = .TRUE.
     305    ENDIF
     306
     307 END SUBROUTINE gasdev_s
     308
     309!--------------------------------------------------------------------------------------------------!
     310! Description:
     311! ------------
     312!> Initialize or reinitialize the random generator state space to vectors of size length.
     313!> The saved variable seq is hashed (via calls to the module routine ran_hash) to create unique
     314!> starting seeds, different for each vector component.
     315!--------------------------------------------------------------------------------------------------!
     316 SUBROUTINE ran_init( length )
     317
     318    USE kinds
     319
     320    IMPLICIT NONE
     321
     322    INTEGER(isp), INTENT(IN) ::  length  !<
     323
     324    INTEGER(isp) ::  hg    !<
     325    INTEGER(isp) ::  hgm   !<
     326    INTEGER(isp) ::  hgng  !<
     327
     328    INTEGER(isp) ::  new   !<
     329    INTEGER(isp) ::  j     !<
     330    INTEGER(isp) ::  hgt   !<
     331!
     332!-- Simply return if enough space is already allocated.
     333    IF ( length < lenran )  RETURN
     334
     335    hg = HUGE( 1_isp )
     336    hgm = - hg
     337    hgng = hgm - 1
     338    hgt = hg
     339
     340!-- The following lines check that kind value isp is in fact a 32-bit integer with the usual
     341!-- properties that we expect it to have (under negation and wrap-around addition).
     342!-- If all of these tests are satisfied, then the routines that use this module are portable, even
     343!-- though they go beyond Fortran 90's integer model.
     344
     345    IF ( hg /= 2147483647   )  CALL ran_error( 'arithmetic assumption 1 failed' )
     346    IF ( hgng >= 0          )  CALL ran_error( 'arithmetic assumption 2 failed' )
     347    IF ( hgt + 1 /= hgng    )  CALL ran_error( 'arithmetic assumption 3 failed' )
     348    IF ( NOT( hg ) >= 0     )  CALL ran_error( 'arithmetic assumption 4 failed' )
     349    IF ( NOT( hgng ) < 0    )  CALL ran_error( 'arithmetic assumption 5 failed' )
     350    IF ( hg + hgng >= 0     )  CALL ran_error( 'arithmetic assumption 6 failed' )
     351    IF ( NOT( - 1_isp ) < 0 )  CALL ran_error( 'arithmetic assumption 7 failed' )
     352    IF ( NOT( 0_isp ) >= 0  )  CALL ran_error( 'arithmetic assumption 8 failed' )
     353    IF ( NOT( 1_isp ) >= 0  )  CALL ran_error( 'arithmetic assumption 9 failed' )
     354
     355    IF ( lenran > 0 )  THEN                          ! Reallocate space, or ...
     356
     357       ranseeds => reallocate( ranseeds, length, 5 )
     358       ranv => reallocate( ranv, length - 1 )
     359       new = lenran + 1
     360
     361    ELSE                                            ! Allocate space.
     362
     363       ALLOCATE( ranseeds(length,5) )
     364       ALLOCATE( ranv(length-1) )
     365       new = 1                               !- Index of first location not yet initialized.
     366       amm = NEAREST( 1.0_wp, - 1.0_wp ) / hgng
     367!
     368!--    Use of nearest is to ensure that returned random deviates are strictly less than 1.0.
     369       IF  ( amm * hgng >= 1.0  .OR.  amm * hgng <= 0.0 )                                          &
     370          CALL ran_error( 'arithmetic assumption 10 failed' )
     371
     372    END IF
     373!
     374!-- Set starting values, unique by seq and vector component.
     375    ranseeds(new:,1) = seq
     376    ranseeds(new:,2:5) = SPREAD( arth( new, 1, SIZE( ranseeds(new:,1) ) ), 2, 4 )
     377
     378    DO  j = 1, 4   ! Hash them.
     379       CALL ran_hash( ranseeds(new:,j), ranseeds(new:,j+1) )
     380    END DO
     381
     382    WHERE ( ranseeds(new:,1:3) < 0 )                                                              &
     383       ranseeds(new:,1:3) = NOT( ranseeds(new:,1:3) )  ! Enforce nonnegativity.
     384
     385    WHERE ( ranseeds(new:,4:5) == 0 ) ranseeds(new:,4:5) = 1 ! Enforce nonzero.
     386
     387    IF ( new == 1 )  THEN ! Set scalar seeds.
     388
     389       iran0 = ranseeds(1,1)
     390       jran0 = ranseeds(1,2)
     391       kran0 = ranseeds(1,3)
     392       mran0 = ranseeds(1,4)
     393       nran0 = ranseeds(1,5)
     394       rans  = nran0
     395
     396    END IF
     397
     398    IF ( length > 1 )  THEN   ! Point to vector seeds.
     399
     400       iran => ranseeds(2:,1)
     401       jran => ranseeds(2:,2)
     402       kran => ranseeds(2:,3)
     403       mran => ranseeds(2:,4)
     404       nran => ranseeds(2:,5)
     405       ranv = nran
     406
     407    END IF
     408
     409    lenran = length
     410
     411 END SUBROUTINE ran_init
     412
     413!--------------------------------------------------------------------------------------------------!
    411414! Description:
    412415! ------------
    413416!> User interface to release the workspace used by the random number routines.
    414 !------------------------------------------------------------------------------!
    415    SUBROUTINE ran_deallocate
    416    
    417       IF  ( lenran > 0 ) THEN
    418      
    419          DEALLOCATE(ranseeds, ranv)
    420          NULLIFY(ranseeds, ranv, iran, jran, kran, mran, nran)
    421          lenran = 0
    422          
    423       END IF
    424      
    425    END SUBROUTINE ran_deallocate
    426 
    427 !------------------------------------------------------------------------------!
    428 ! Description:
    429 ! ------------
    430 !> User interface for seeding the random number routines.
    431 !> Syntax is exactly like Fortran 90's random_seed routine,
    432 !> with one additional argument keyword: random_sequence, set to any integer
    433 !> value, causes an immediate new initialization, seeded by that integer.
    434 !------------------------------------------------------------------------------!
    435    SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
    436    
    437       IMPLICIT NONE
    438      
    439       INTEGER(isp), OPTIONAL, INTENT(IN)  ::  random_sequence   !<
    440       INTEGER(isp), OPTIONAL, INTENT(OUT) ::  state_size        !<
    441      
    442       INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN)  ::  put   !<
    443       INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) ::  get   !<
    444      
    445       IF  ( present(state_size) ) THEN
    446      
    447          state_size = 5 * lenran
    448          
    449       ELSE IF  ( present(put) ) THEN
    450      
    451          IF  ( lenran == 0 ) RETURN
    452          
    453          ranseeds = reshape( put,shape(ranseeds) )
    454          
    455          WHERE (ranseeds(:,1:3) < 0) ranseeds(: ,1:3) = not(ranseeds(: ,1:3))
    456          !- Enforce nonnegativity and nonzero conditions on any user-supplied seeds.
    457          
    458          WHERE (ranseeds(:,4:5) == 0) ranseeds(:,4:5) = 1
    459          
    460          iran0 = ranseeds(1,1)
    461          jran0 = ranseeds(1,2)
    462          kran0 = ranseeds(1,3)
    463          mran0 = ranseeds(1,4)
    464          nran0 = ranseeds(1,5)
    465          
    466       ELSE IF  ( present(get) ) THEN
    467      
    468          IF  (lenran == 0) RETURN
    469          
    470          ranseeds(1,1:5) = (/ iran0,jran0,kran0,mran0,nran0 /)
    471          get = reshape( ranseeds, shape(get) )
    472          
    473       ELSE IF  ( present(random_sequence) ) THEN
    474      
    475          CALL ran_deallocate
    476          seq = random_sequence
    477          
    478       END IF
    479      
    480    END SUBROUTINE random_seed_parallel
    481 
    482 !------------------------------------------------------------------------------!
    483 ! Description:
    484 ! ------------
    485 !> DES-like hashing of two 32-bit integers, using shifts,
    486 !> xor's, and adds to make the internal nonlinear function.
    487 !------------------------------------------------------------------------------!
    488    SUBROUTINE ran_hash_v( il, ir )
    489    
    490       IMPLICIT NONE
    491      
    492       INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  il   !<
    493       INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  ir   !<
    494      
    495       INTEGER(isp), DIMENSION(size(il)) ::  is   !<
    496      
    497       INTEGER(isp) :: j   !<
    498      
    499       DO j=1,4
    500      
    501          is = ir
    502          ir = ieor( ir, ishft(ir,5) ) + 1422217823
    503          ir = ieor( ir, ishft(ir,-16) ) + 1842055030
    504          ir = ieor( ir, ishft(ir,9) ) + 80567781
    505          ir = ieor( il, ir )
    506          il = is
    507       END DO
    508      
    509    END SUBROUTINE ran_hash_v
    510 
    511 !------------------------------------------------------------------------------!
    512 ! Description:
    513 ! ------------
    514 !> User interface to process error-messages
    515 !> produced by the random_number_parallel module
    516 !------------------------------------------------------------------------------!
    517    SUBROUTINE ran_error(string)
    518 
    519       USE control_parameters,                                                &
    520         ONLY: message_string
    521 
    522       CHARACTER(LEN=*), INTENT(IN) ::  string   !< Error message string
    523 
    524       message_string = 'incompatible integer arithmetic: ' // TRIM( string )
    525       CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 )
    526 
    527    END SUBROUTINE ran_error
    528 
    529 !------------------------------------------------------------------------------!
    530 ! Description:
    531 ! ------------
    532 !> Reallocates the generators state space "ranseeds" to vectors of size length.
    533 !------------------------------------------------------------------------------!
    534    FUNCTION reallocate_iv( p, n )
    535 
    536       USE control_parameters,                                                &
    537         ONLY: message_string
    538    
    539       INTEGER(isp), DIMENSION(:), POINTER ::  p               !<
    540       INTEGER(isp), DIMENSION(:), POINTER ::  reallocate_iv   !<
    541      
    542       INTEGER(isp), INTENT(IN) ::  n   !<
    543      
    544       INTEGER(isp) ::  nold   !<
    545       INTEGER(isp) ::  ierr   !<
    546      
    547       ALLOCATE(reallocate_iv(n),stat=ierr)
    548      
    549       IF (ierr /= 0) THEN
    550          message_string = 'problem in attempt to allocate memory'
    551          CALL message( 'reallocate_iv', 'PA0454', 1, 2, 0, 6, 0 )
    552       END IF
    553 
    554       IF (.not. associated(p)) RETURN
    555      
    556       nold = size(p)
    557      
    558       reallocate_iv(1:min(nold,n)) = p(1:min(nold,n))
    559      
    560       DEALLOCATE(p)
    561      
    562    END FUNCTION reallocate_iv
    563    
    564    FUNCTION reallocate_im( p, n, m )
    565 
    566       USE control_parameters,                                                &
    567         ONLY: message_string
    568    
    569       INTEGER(isp), DIMENSION(:,:), POINTER ::  p               !<
    570       INTEGER(isp), DIMENSION(:,:), POINTER ::  reallocate_im   !<
    571      
    572       INTEGER(isp), INTENT(IN) ::  m   !<
    573       INTEGER(isp), INTENT(IN) ::  n   !<
    574      
    575       INTEGER(isp) ::  mold   !<
    576       INTEGER(isp) ::  nold   !<
    577       INTEGER(isp) ::  ierr   !<
    578      
    579       ALLOCATE(reallocate_im(n,m),stat=ierr)
    580      
    581       IF (ierr /= 0) THEN
    582          message_string = 'problem in attempt to allocate memory'
    583          CALL message( 'reallocate_im', 'PA0454', 1, 2, 0, 6, 0 )
    584       END IF
    585 
    586       IF (.not. associated(p)) RETURN
    587      
    588       nold = size(p,1)
    589       mold = size(p,2)
    590      
    591       reallocate_im(1:min(nold,n),1:min(mold,m)) =                             &
    592          p(1:min(nold,n),1:min(mold,m))
    593          
    594       DEALLOCATE(p)
    595      
    596    END FUNCTION reallocate_im
    597 
    598 !------------------------------------------------------------------------------!
    599 ! Description:
    600 ! ------------
    601 !> Reallocates the generators state space "ranseeds" to vectors of size length.
    602 !------------------------------------------------------------------------------!
    603    FUNCTION arth_i(first,increment,n)
    604    
    605       INTEGER(isp), INTENT(IN) ::  first       !<
    606       INTEGER(isp), INTENT(IN) ::  increment   !<
    607       INTEGER(isp), INTENT(IN) ::  n           !<
    608      
    609       INTEGER(isp), DIMENSION(n) ::  arth_i    !<
    610      
    611       INTEGER(isp) ::  k      !<
    612       INTEGER(isp) ::  k2     !<
    613       INTEGER(isp) ::  temp   !<
    614      
    615       INTEGER(isp), PARAMETER ::  npar_arth=16   !<
    616       INTEGER(isp), PARAMETER ::  npar2_arth=8   !<
    617      
    618       IF (n > 0) arth_i(1) = first
    619      
    620       IF (n <= npar_arth) THEN
    621      
    622          DO k=2,n
    623             arth_i(k) = arth_i(k-1)+increment
    624          END DO
    625          
    626       ELSE
    627      
    628          DO k=2,npar2_arth
    629             arth_i(k) = arth_i(k-1) + increment
    630          END DO
    631          
    632          temp = increment * npar2_arth
    633          k = npar2_arth
    634          
    635          DO
    636             IF (k >= n) EXIT
    637             k2 = k + k
    638             arth_i(k+1:min(k2,n)) = temp + arth_i(1:min(k,n-k))
    639             temp = temp + temp
    640             k = k2
    641          END DO
    642          
    643       END IF
    644      
    645    END FUNCTION arth_i
     417!--------------------------------------------------------------------------------------------------!
     418 SUBROUTINE ran_deallocate
     419
     420    IF ( lenran > 0 )  THEN
     421
     422       DEALLOCATE( ranseeds, ranv )
     423       NULLIFY( ranseeds, ranv, iran, jran, kran, mran, nran )
     424       lenran = 0
     425
     426    END IF
     427
     428 END SUBROUTINE ran_deallocate
     429
     430!--------------------------------------------------------------------------------------------------!
     431! Description:
     432! ------------
     433!> User interface for seeding the random number routines.
     434!> Syntax is exactly like Fortran 90's random_seed routine, with one additional argument keyword:
     435!> random_sequence, set to any integer value, causes an immediate new initialization, seeded by that
     436!> integer.
     437!--------------------------------------------------------------------------------------------------!
     438 SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
     439
     440    IMPLICIT NONE
     441
     442    INTEGER(isp), OPTIONAL, INTENT(IN)  ::  random_sequence  !<
     443    INTEGER(isp), OPTIONAL, INTENT(OUT) ::  state_size       !<
     444
     445    INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN)  ::  put  !<
     446    INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) ::  get  !<
     447
     448    IF ( PRESENT( state_size ) )  THEN
     449
     450       state_size = 5 * lenran
     451
     452    ELSE IF ( PRESENT( put ) )  THEN
     453
     454       IF ( lenran == 0 )  RETURN
     455
     456       ranseeds = RESHAPE( put, SHAPE( ranseeds ) )
     457!
     458!--    Enforce nonnegativity and nonzero conditions on any user-supplied seeds.
     459       WHERE ( ranseeds(:,1:3) < 0 ) ranseeds(:,1:3) = NOT( ranseeds(:,1:3) )
     460
     461       WHERE ( ranseeds(:,4:5) == 0 ) ranseeds(:,4:5) = 1
     462
     463       iran0 = ranseeds(1,1)
     464       jran0 = ranseeds(1,2)
     465       kran0 = ranseeds(1,3)
     466       mran0 = ranseeds(1,4)
     467       nran0 = ranseeds(1,5)
     468
     469    ELSE IF ( PRESENT( get ) )  THEN
     470
     471       IF ( lenran == 0 )  RETURN
     472
     473       ranseeds(1,1:5) = (/ iran0, jran0, kran0, mran0, nran0 /)
     474       get = RESHAPE( ranseeds, SHAPE( get ) )
     475
     476    ELSE IF ( PRESENT( random_sequence ) ) THEN
     477
     478       CALL ran_deallocate
     479       seq = random_sequence
     480
     481    END IF
     482
     483 END SUBROUTINE random_seed_parallel
     484
     485!--------------------------------------------------------------------------------------------------!
     486! Description:
     487! ------------
     488!> DES-like hashing of two 32-bit integers, using shifts, xor's, and adds to make the internal
     489!> nonlinear function.
     490!--------------------------------------------------------------------------------------------------!
     491 SUBROUTINE ran_hash_v( il, ir )
     492
     493    IMPLICIT NONE
     494
     495    INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  il  !<
     496    INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  ir  !<
     497
     498    INTEGER(isp), DIMENSION(size(il)) ::  is  !<
     499
     500    INTEGER(isp) :: j  !<
     501
     502    DO  j = 1, 4
     503
     504       is = ir
     505       ir = IEOR( ir, ISHFT( ir, 5 ) ) + 1422217823
     506       ir = IEOR( ir, ISHFT( ir, - 16 ) ) + 1842055030
     507       ir = IEOR( ir, ISHFT( ir, 9 ) ) + 80567781
     508       ir = IEOR( il, ir )
     509       il = is
     510    END DO
     511
     512 END SUBROUTINE ran_hash_v
     513
     514!--------------------------------------------------------------------------------------------------!
     515! Description:
     516! ------------
     517!> User interface to process error-messages produced by the random_number_parallel module.
     518!--------------------------------------------------------------------------------------------------!
     519 SUBROUTINE ran_error( string )
     520
     521    USE control_parameters,                                                                        &
     522      ONLY: message_string
     523
     524    CHARACTER(LEN=*), INTENT(IN) ::  string  !< Error message string
     525
     526    message_string = 'incompatible integer arithmetic: ' // TRIM( string )
     527    CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 )
     528
     529 END SUBROUTINE ran_error
     530
     531!--------------------------------------------------------------------------------------------------!
     532! Description:
     533! ------------
     534!> Reallocates the generators state space "ranseeds" to vectors of size length.
     535!--------------------------------------------------------------------------------------------------!
     536 FUNCTION reallocate_iv( p, n )
     537
     538    USE control_parameters,                                                                        &
     539      ONLY: message_string
     540
     541    INTEGER(isp), DIMENSION(:), POINTER ::  p              !<
     542    INTEGER(isp), DIMENSION(:), POINTER ::  reallocate_iv  !<
     543
     544    INTEGER(isp), INTENT(IN) ::  n  !<
     545
     546    INTEGER(isp) ::  nold  !<
     547    INTEGER(isp) ::  ierr  !<
     548
     549    ALLOCATE( reallocate_iv(n), stat = ierr )
     550
     551    IF ( ierr /= 0 )  THEN
     552       message_string = 'problem in attempt to allocate memory'
     553       CALL message( 'reallocate_iv', 'PA0454', 1, 2, 0, 6, 0 )
     554    END IF
     555
     556    IF ( .NOT. ASSOCIATED( p ) )  RETURN
     557
     558    nold = SIZE( p )
     559
     560    reallocate_iv(1:MIN( nold, n )) = p(1:MIN( nold, n ))
     561
     562    DEALLOCATE( p )
     563
     564 END FUNCTION reallocate_iv
     565
     566 FUNCTION reallocate_im( p, n, m )
     567
     568    USE control_parameters,                                                                        &
     569      ONLY: message_string
     570
     571    INTEGER(isp), DIMENSION(:,:), POINTER ::  p              !<
     572    INTEGER(isp), DIMENSION(:,:), POINTER ::  reallocate_im  !<
     573
     574    INTEGER(isp), INTENT(IN) ::  m  !<
     575    INTEGER(isp), INTENT(IN) ::  n  !<
     576
     577    INTEGER(isp) ::  mold  !<
     578    INTEGER(isp) ::  nold  !<
     579    INTEGER(isp) ::  ierr  !<
     580
     581    ALLOCATE( reallocate_im(n,m), stat = ierr )
     582
     583    IF ( ierr /= 0 )  THEN
     584       message_string = 'problem in attempt to allocate memory'
     585       CALL message( 'reallocate_im', 'PA0454', 1, 2, 0, 6, 0 )
     586    END IF
     587
     588    IF ( .NOT. ASSOCIATED( p ) )  RETURN
     589
     590    nold = SIZE( p, 1 )
     591    mold = SIZE( p, 2 )
     592
     593    reallocate_im(1:MIN( nold, n ),1:MIN( mold, m )) = p(1:MIN( nold, n ),1:MIN( mold, m ))
     594
     595    DEALLOCATE(p)
     596
     597 END FUNCTION reallocate_im
     598
     599!--------------------------------------------------------------------------------------------------!
     600! Description:
     601! ------------
     602!> Reallocates the generators state space "ranseeds" to vectors of size length.
     603!--------------------------------------------------------------------------------------------------!
     604 FUNCTION arth_i( first, increment, n )
     605
     606    INTEGER(isp), INTENT(IN) ::  first      !<
     607    INTEGER(isp), INTENT(IN) ::  increment  !<
     608    INTEGER(isp), INTENT(IN) ::  n          !<
     609
     610    INTEGER(isp), DIMENSION(n) ::  arth_i  !<
     611
     612    INTEGER(isp) ::  k     !<
     613    INTEGER(isp) ::  k2    !<
     614    INTEGER(isp) ::  temp  !<
     615
     616    INTEGER(isp), PARAMETER ::  npar_arth = 16  !<
     617    INTEGER(isp), PARAMETER ::  npar2_arth = 8  !<
     618
     619    IF ( n > 0 )  arth_i(1) = first
     620
     621    IF ( n <= npar_arth )  THEN
     622
     623       DO  k = 2, n
     624          arth_i(k) = arth_i(k-1) + increment
     625       END DO
     626
     627    ELSE
     628
     629       DO  k = 2, npar2_arth
     630          arth_i(k) = arth_i(k-1) + increment
     631       END DO
     632
     633       temp = increment * npar2_arth
     634       k = npar2_arth
     635
     636       DO
     637          IF ( k >= n )  EXIT
     638          k2 = k + k
     639          arth_i(k+1:MIN( k2, n )) = temp + arth_i(1:MIN( k, n - k ))
     640          temp = temp + temp
     641          k = k2
     642       END DO
     643
     644    END IF
     645
     646 END FUNCTION arth_i
    646647
    647648END MODULE random_generator_parallel
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |