source: palm/trunk/SOURCE/random_generator_parallel_mod.f90 @ 3436

Last change on this file since 3436 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

  • Property svn:keywords set to Id
File size: 19.0 KB
RevLine 
[1850]1!> @file random_generator_parallel_mod.f90
[1400]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1400]4!
[2000]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.
[1400]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/>.
16!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[1400]18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
[2173]23!
[1400]24! Former revisions:
25! -----------------
26! $Id: random_generator_parallel_mod.f90 3241 2018-09-12 15:02:00Z gronemeier $
[3241]27! unused variable removed
28!
29! 2718 2018-01-02 08:49:38Z maronga
[2716]30! Corrected "Former revisions" section
31!
32! 2696 2017-12-14 17:12:51Z kanani
33! Change in file header (GPL part)
[2696]34! Formatting (MS)
35!
36! 2255 2017-06-07 11:58:09Z knoop
[2255]37! Implemented PALM specific error message handling
38!
39! 2173 2017-03-08 15:56:54Z knoop
[1400]40!
[2173]41! 2172 2017-03-08 15:55:25Z knoop
42! Bugfix: added global initialization routine and removed global array
43!
[2145]44! 2144 2017-02-07 15:32:45Z knoop
45! Bugfix: compile time parameter replacement prevented
46!
[2001]47! 2000 2016-08-20 18:09:15Z knoop
48! Forced header and separation lines into 80 columns
49!
[1851]50! 1850 2016-04-08 13:29:27Z maronga
51! Module renamed
52!
53!
[1683]54! 1682 2015-10-07 23:56:08Z knoop
55! Code annotations made doxygen readable
56!
[1401]57! 1400 2014-05-09 14:03:54Z knoop
58! Initial revision
59!
[1400]60!
61! Description:
62! ------------
[1682]63!> This module contains and supports the random number generating routine ran_parallel.
64!> ran_parallel returns a uniform random deviate between 0.0 and 1.0
65!> (exclusive of the end point values).
66!> Additionally it provides the generator with five integer for use as initial state space.
67!> The first tree integers (iran, jran, kran) are maintained as non negative values,
68!> while the last two (mran, nran) have 32-bit nonzero values.
69!> Also provided by this module is support for initializing or reinitializing
70!> the state space to a desired standard sequence number, hashing the initial
71!> values to random values, and allocating and deallocating the internal workspace
72!> Random number generator, produces numbers equally distributed in interval
73!>
74!> This routine is taken from the "numerical recipies vol. 2"
[1400]75!------------------------------------------------------------------------------!
[1682]76MODULE random_generator_parallel
77 
[1400]78
79   USE kinds
80   
81   IMPLICIT NONE
82   
83   PRIVATE
84   PUBLIC random_number_parallel, random_seed_parallel, random_dummy,          &
[2172]85          id_random_array, seq_random_array, init_parallel_random_generator
[1400]86   
[1682]87   INTEGER(isp), SAVE :: lenran=0             !<
88   INTEGER(isp), SAVE :: seq=0                !<
89   INTEGER(isp), SAVE :: iran0                !<
90   INTEGER(isp), SAVE :: jran0                !<
91   INTEGER(isp), SAVE :: kran0                !<
92   INTEGER(isp), SAVE :: mran0                !<
93   INTEGER(isp), SAVE :: nran0                !<
94   INTEGER(isp), SAVE :: rans                 !<
[1400]95   
[1682]96   INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds   !<
[1400]97   
[1682]98   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: iran   !<
99   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: jran   !<
100   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: kran   !<
101   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: mran   !<
102   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: nran   !<
103   INTEGER(isp), DIMENSION(:), POINTER, SAVE :: ranv   !<
[1400]104   
105   
106   
[1682]107   INTEGER(isp), DIMENSION(:,:), ALLOCATABLE   ::  id_random_array    !<
108   INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE ::  seq_random_array   !<
[1400]109   
[1682]110   REAL(wp), SAVE :: amm   !<
[1400]111   
[1682]112   REAL(wp) :: random_dummy=0.0   !<
[1400]113   
[2172]114   INTERFACE init_parallel_random_generator
115      MODULE PROCEDURE init_parallel_random_generator
116   END INTERFACE
117   
[1400]118   INTERFACE random_number_parallel
119      MODULE PROCEDURE ran0_s
120   END INTERFACE
121   
122   INTERFACE random_seed_parallel
123      MODULE PROCEDURE random_seed_parallel
124   END INTERFACE
125   
126   INTERFACE ran_hash
127      MODULE PROCEDURE ran_hash_v
128   END INTERFACE
129   
130   INTERFACE reallocate
131      MODULE PROCEDURE reallocate_iv,reallocate_im
132   END INTERFACE
133   
134   INTERFACE arth
135      MODULE PROCEDURE arth_i
136   END INTERFACE
137
138 CONTAINS
139 
140!------------------------------------------------------------------------------!
141! Description:
142! ------------
[2172]143!> Initialize the parallel random number generator for a specific subdomain
144!------------------------------------------------------------------------------!
[3241]145   SUBROUTINE init_parallel_random_generator( nx, nys, nyn, nxl, nxr )
[2172]146
147      USE kinds
148
149      USE control_parameters,                                                  &
150          ONLY: ensemble_member_nr
151
152      IMPLICIT NONE
153
154      INTEGER(isp), INTENT(IN) ::  nx    !<
155      INTEGER(isp), INTENT(IN) ::  nys   !<
156      INTEGER(isp), INTENT(IN) ::  nyn   !<
157      INTEGER(isp), INTENT(IN) ::  nxl   !<
158      INTEGER(isp), INTENT(IN) ::  nxr   !<
159
160      INTEGER(iwp) ::  i   !<
161      INTEGER(iwp) ::  j   !<
162
163!--   Allocate ID-array and state-space-array
164      ALLOCATE ( seq_random_array(5,nys:nyn,nxl:nxr) )
165      ALLOCATE ( id_random_array(nys:nyn,nxl:nxr) )
166      seq_random_array = 0
167      id_random_array  = 0
168
[2696]169!--   Asigning an ID to every vertical gridpoint column
170!--   dependig on the ensemble run number.
171      DO i=nxl, nxr
172         DO j=nys, nyn
173            id_random_array(j,i) = j*(nx+1.0_wp) + i + 1.0_wp + 1E6 * &
174                                   ( ensemble_member_nr - 1000 )
175         ENDDO
176      ENDDO
177!--   Initializing with random_seed_parallel for every vertical
178!--   gridpoint column.
179      random_dummy=0
180      DO i = nxl, nxr
181         DO j = nys, nyn
182            CALL random_seed_parallel (random_sequence=id_random_array(j, i))
183            CALL random_number_parallel (random_dummy)
184            CALL random_seed_parallel (get=seq_random_array(:, j, i))
185         ENDDO
186      ENDDO
[2172]187
188   END SUBROUTINE init_parallel_random_generator
189 
190!------------------------------------------------------------------------------!
191! Description:
192! ------------
[1682]193!> Lagged Fibonacci generator combined with a Marsaglia shiftsequence.
194!> Returns as harvest a uniform random deviate between 0.0 and 1.0 (exclusive of the end point values).
195!> This generator has the same calling and initialization conventions as Fortran 90's random_number routine.
196!> Use random_seed_parallel to initialize or reinitialize to a particular sequence.
197!> The period of this generator is about 2.0 x 10^28, and it fully vectorizes.
198!> Validity of the integer model assumed by this generator is tested at initialization.
[1400]199!------------------------------------------------------------------------------!
200   SUBROUTINE ran0_s(harvest)
201
202      USE kinds
203     
204      IMPLICIT NONE
205     
[1682]206      REAL(wp), INTENT(OUT) :: harvest   !<
[1400]207     
208      IF  (lenran < 1) CALL ran_init(1)  !- Initialization routine in ran_state.
209     
210      !- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69.
211      rans = iran0 - kran0   
212     
213      IF  (rans < 0) rans = rans + 2147483579_isp
214     
215      iran0 = jran0
216      jran0 = kran0
217      kran0 = rans
218     
219      nran0 = ieor( nran0, ishft (nran0, 13) ) !- Update Marsaglia shift sequence with period 2^32 - 1.
220      nran0 = ieor( nran0, ishft (nran0, -17) )
221      nran0 = ieor( nran0, ishft (nran0, 5) )
222     
223      rans  = ieor( nran0, rans )   !- Combine the generators.
224     
225      harvest = amm * merge( rans, not(rans), rans < 0 ) !- Make the result positive definite (note that amm is negative).
226     
227   END SUBROUTINE ran0_s
228
229!------------------------------------------------------------------------------!
230! Description:
231! ------------
[1682]232!> Initialize or reinitialize the random generator state space to vectors of size length.
233!> The saved variable seq is hashed (via calls to the module routine ran_hash)
234!> to create unique starting seeds, different for each vector component.
[1400]235!------------------------------------------------------------------------------!
236   SUBROUTINE ran_init( length )
237   
238      USE kinds
239     
240      IMPLICIT NONE
241     
[1682]242      INTEGER(isp), INTENT(IN) ::  length   !<
[1400]243   
[2144]244      INTEGER(isp) ::  hg    !<
245      INTEGER(isp) ::  hgm   !<
246      INTEGER(isp) ::  hgng  !<
[1400]247     
[1682]248      INTEGER(isp) ::  new   !<
249      INTEGER(isp) ::  j     !<
250      INTEGER(isp) ::  hgt   !<
[1400]251     
252      IF ( length < lenran ) RETURN !- Simply return if enough space is already allocated.
253     
[2144]254      hg=huge(1_isp)
255      hgm=-hg
256      hgng=hgm-1
[1400]257      hgt = hg
258     
259      !- The following lines check that kind value isp is in fact a 32-bit integer
260      !- with the usual properties that we expect it to have (under negation and wrap-around addition).
261      !- If all of these tests are satisfied, then the routines that use this module are portable,
262      !- even though they go beyond Fortran 90's integer model.
263     
[2255]264      IF  ( hg /= 2147483647 ) CALL ran_error('arithmetic assumption 1 failed')
265      IF  ( hgng >= 0 )        CALL ran_error('arithmetic assumption 2 failed')
266      IF  ( hgt+1 /= hgng )    CALL ran_error('arithmetic assumption 3 failed')
267      IF  ( not(hg) >= 0 )     CALL ran_error('arithmetic assumption 4 failed')
268      IF  ( not(hgng) < 0 )    CALL ran_error('arithmetic assumption 5 failed')
269      IF  ( hg+hgng >= 0 )     CALL ran_error('arithmetic assumption 6 failed')
270      IF  ( not(-1_isp) < 0 )  CALL ran_error('arithmetic assumption 7 failed')
271      IF  ( not(0_isp) >= 0 )  CALL ran_error('arithmetic assumption 8 failed')
272      IF  ( not(1_isp) >= 0 )  CALL ran_error('arithmetic assumption 9 failed')
[1400]273     
274      IF  ( lenran > 0) THEN                          !- Reallocate space, or ...
275     
276         ranseeds => reallocate( ranseeds, length, 5)
277         ranv => reallocate( ranv, length-1)
278         new = lenran+1
279         
280      ELSE                                            !- allocate space.
281     
282         ALLOCATE(ranseeds(length,5))
283         ALLOCATE(ranv(length-1))
284         new = 1   !- Index of first location not yet initialized.
285         amm = nearest(1.0_wp,-1.0_wp)/hgng
286         !- Use of nearest is to ensure that returned random deviates are strictly lessthan 1.0.
287         IF  (amm*hgng >= 1.0 .or. amm*hgng <= 0.0)                            &
[2255]288            CALL ran_error('arithmetic assumption 10 failed')
[1400]289           
290      END IF 
291     
292      !- Set starting values, unique by seq and vector component.
293      ranseeds(new:,1) = seq
294      ranseeds(new:,2:5)=spread(arth(new,1,size(ranseeds(new:,1))),2,4)
295     
296      DO j=1,4   !- Hash them.
297         CALL ran_hash(ranseeds(new:,j), ranseeds(new:,j+1))
298      END DO
299     
300      WHERE (ranseeds (new: ,1:3) < 0)                                         & 
301         ranseeds(new: ,1:3)=not(ranseeds(new: ,1:3))  !- Enforce nonnegativity.
302         
303      WHERE (ranseeds(new: ,4:5) == 0) ranseeds(new: ,4:5)=1 !- Enforce nonzero.
304     
305      IF  (new == 1) THEN !- Set scalar seeds.
306     
307         iran0 = ranseeds(1,1)
308         jran0 = ranseeds(1,2)
309         kran0 = ranseeds(1,3)
310         mran0 = ranseeds(1,4)
311         nran0 = ranseeds(1,5)
312         rans  = nran0
313         
314      END IF
315     
316      IF  (length > 1) THEN   !- Point to vector seeds.
317     
318         iran => ranseeds(2:,1)
319         jran => ranseeds(2:,2)
320         kran => ranseeds(2:,3)
321         mran => ranseeds(2:,4)
322         nran => ranseeds(2:,5)
323         ranv = nran
324         
325      END IF
326     
327      lenran = length
328     
329   END SUBROUTINE ran_init
330
331!------------------------------------------------------------------------------!
332! Description:
333! ------------
[1682]334!> User interface to release the workspace used by the random number routines.
[1400]335!------------------------------------------------------------------------------!
336   SUBROUTINE ran_deallocate
337   
338      IF  ( lenran > 0 ) THEN
339     
340         DEALLOCATE(ranseeds, ranv)
341         NULLIFY(ranseeds, ranv, iran, jran, kran, mran, nran)
342         lenran = 0
343         
344      END IF
345     
346   END SUBROUTINE ran_deallocate
347
348!------------------------------------------------------------------------------!
349! Description:
350! ------------
[1682]351!> User interface for seeding the random number routines.
352!> Syntax is exactly like Fortran 90's random_seed routine,
353!> with one additional argument keyword: random_sequence, set to any integer
354!> value, causes an immediate new initialization, seeded by that integer.
[1400]355!------------------------------------------------------------------------------!
356   SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
357   
358      IMPLICIT NONE
359     
[1682]360      INTEGER(isp), OPTIONAL, INTENT(IN)  ::  random_sequence   !<
361      INTEGER(isp), OPTIONAL, INTENT(OUT) ::  state_size        !<
[1400]362     
[1682]363      INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN)  ::  put   !<
364      INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) ::  get   !<
[1400]365     
366      IF  ( present(state_size) ) THEN
367     
368         state_size = 5 * lenran
369         
370      ELSE IF  ( present(put) ) THEN
371     
372         IF  ( lenran == 0 ) RETURN
373         
374         ranseeds = reshape( put,shape(ranseeds) )
375         
376         WHERE (ranseeds(:,1:3) < 0) ranseeds(: ,1:3) = not(ranseeds(: ,1:3))
377         !- Enforce nonnegativity and nonzero conditions on any user-supplied seeds.
378         
379         WHERE (ranseeds(:,4:5) == 0) ranseeds(:,4:5) = 1
380         
381         iran0 = ranseeds(1,1)
382         jran0 = ranseeds(1,2)
383         kran0 = ranseeds(1,3)
384         mran0 = ranseeds(1,4)
385         nran0 = ranseeds(1,5)
386         
387      ELSE IF  ( present(get) ) THEN
388     
389         IF  (lenran == 0) RETURN
390         
391         ranseeds(1,1:5) = (/ iran0,jran0,kran0,mran0,nran0 /)
392         get = reshape( ranseeds, shape(get) )
393         
394      ELSE IF  ( present(random_sequence) ) THEN
395     
396         CALL ran_deallocate
397         seq = random_sequence
398         
399      END IF
400     
401   END SUBROUTINE random_seed_parallel
402
403!------------------------------------------------------------------------------!
404! Description:
405! ------------
[1682]406!> DES-like hashing of two 32-bit integers, using shifts,
407!> xor's, and adds to make the internal nonlinear function.
[1400]408!------------------------------------------------------------------------------!
409   SUBROUTINE ran_hash_v( il, ir )
410   
411      IMPLICIT NONE
412     
[1682]413      INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  il   !<
414      INTEGER(isp), DIMENSION(:), INTENT(INOUT) ::  ir   !<
[1400]415     
[1682]416      INTEGER(isp), DIMENSION(size(il)) ::  is   !<
[1400]417     
[1682]418      INTEGER(isp) :: j   !<
[1400]419     
420      DO j=1,4
421     
422         is = ir
423         ir = ieor( ir, ishft(ir,5) ) + 1422217823
424         ir = ieor( ir, ishft(ir,-16) ) + 1842055030
425         ir = ieor( ir, ishft(ir,9) ) + 80567781
426         ir = ieor( il, ir )
427         il = is
428      END DO
429     
430   END SUBROUTINE ran_hash_v
431
432!------------------------------------------------------------------------------!
433! Description:
434! ------------
[1682]435!> User interface to process error-messages
436!> produced by the random_number_parallel module
[1400]437!------------------------------------------------------------------------------!
438   SUBROUTINE ran_error(string)
[2255]439
440      USE control_parameters,                                                &
441        ONLY: message_string
442
443      CHARACTER(LEN=*), INTENT(IN) ::  string   !< Error message string
444
445      message_string = 'incompatible integer arithmetic: ' // TRIM( string )
446      CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 )
447
[1400]448   END SUBROUTINE ran_error
449
450!------------------------------------------------------------------------------!
451! Description:
452! ------------
[1682]453!> Reallocates the generators state space "ranseeds" to vectors of size length.
[1400]454!------------------------------------------------------------------------------!
455   FUNCTION reallocate_iv( p, n )
[2255]456
457      USE control_parameters,                                                &
458        ONLY: message_string
[1400]459   
[1682]460      INTEGER(isp), DIMENSION(:), POINTER ::  p               !<
461      INTEGER(isp), DIMENSION(:), POINTER ::  reallocate_iv   !<
[1400]462     
[1682]463      INTEGER(isp), INTENT(IN) ::  n   !<
[1400]464     
[1682]465      INTEGER(isp) ::  nold   !<
466      INTEGER(isp) ::  ierr   !<
[1400]467     
468      ALLOCATE(reallocate_iv(n),stat=ierr)
469     
[2255]470      IF (ierr /= 0) THEN
471         message_string = 'problem in attempt to allocate memory'
472         CALL message( 'reallocate_iv', 'PA0454', 1, 2, 0, 6, 0 )
473      END IF
474
[1400]475      IF (.not. associated(p)) RETURN
476     
477      nold = size(p)
478     
479      reallocate_iv(1:min(nold,n)) = p(1:min(nold,n))
480     
481      DEALLOCATE(p)
482     
483   END FUNCTION reallocate_iv
484   
485   FUNCTION reallocate_im( p, n, m )
[2255]486
487      USE control_parameters,                                                &
488        ONLY: message_string
[1400]489   
[1682]490      INTEGER(isp), DIMENSION(:,:), POINTER ::  p               !<
491      INTEGER(isp), DIMENSION(:,:), POINTER ::  reallocate_im   !<
[1400]492     
[1682]493      INTEGER(isp), INTENT(IN) ::  m   !<
494      INTEGER(isp), INTENT(IN) ::  n   !<
[1400]495     
[1682]496      INTEGER(isp) ::  mold   !<
497      INTEGER(isp) ::  nold   !<
498      INTEGER(isp) ::  ierr   !<
[1400]499     
500      ALLOCATE(reallocate_im(n,m),stat=ierr)
501     
[2255]502      IF (ierr /= 0) THEN
503         message_string = 'problem in attempt to allocate memory'
504         CALL message( 'reallocate_im', 'PA0454', 1, 2, 0, 6, 0 )
505      END IF
506
[1400]507      IF (.not. associated(p)) RETURN
508     
509      nold = size(p,1)
510      mold = size(p,2)
511     
512      reallocate_im(1:min(nold,n),1:min(mold,m)) =                             &
513         p(1:min(nold,n),1:min(mold,m))
514         
515      DEALLOCATE(p)
516     
517   END FUNCTION reallocate_im
518
519!------------------------------------------------------------------------------!
520! Description:
521! ------------
[1682]522!> Reallocates the generators state space "ranseeds" to vectors of size length.
[1400]523!------------------------------------------------------------------------------!
524   FUNCTION arth_i(first,increment,n)
525   
[1682]526      INTEGER(isp), INTENT(IN) ::  first       !<
527      INTEGER(isp), INTENT(IN) ::  increment   !<
528      INTEGER(isp), INTENT(IN) ::  n           !<
[1400]529     
[1682]530      INTEGER(isp), DIMENSION(n) ::  arth_i    !<
[1400]531     
[1682]532      INTEGER(isp) ::  k      !<
533      INTEGER(isp) ::  k2     !<
534      INTEGER(isp) ::  temp   !<
[1400]535     
[1682]536      INTEGER(isp), PARAMETER ::  npar_arth=16   !<
537      INTEGER(isp), PARAMETER ::  npar2_arth=8   !<
[1400]538     
539      IF (n > 0) arth_i(1) = first
540     
541      IF (n <= npar_arth) THEN
542     
543         DO k=2,n
544            arth_i(k) = arth_i(k-1)+increment
545         END DO
546         
547      ELSE
548     
549         DO k=2,npar2_arth
550            arth_i(k) = arth_i(k-1) + increment
551         END DO
552         
553         temp = increment * npar2_arth
554         k = npar2_arth
555         
556         DO
557            IF (k >= n) EXIT
558            k2 = k + k
559            arth_i(k+1:min(k2,n)) = temp + arth_i(1:min(k,n-k))
560            temp = temp + temp
561            k = k2
562         END DO
563         
564      END IF
565     
566   END FUNCTION arth_i
567
568END MODULE random_generator_parallel
Note: See TracBrowser for help on using the repository browser.