Home

Context Navigation

random_generator_parallel_mod.f90 @ 4875

Last change on this file since 4875 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 22.4 KB

Rev	Line
[1850]	1	!> @file random_generator_parallel_mod.f90
[4592]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1400]	4	!
[4592]	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.
[1400]	8	!
[4592]	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.
[1400]	12	!
[4592]	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/>.
[1400]	15	!
[4828]	16	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4592]	17	!--------------------------------------------------------------------------------------------------!
[1400]	18	!
[4592]	19	!
[1400]	20	! Current revisions:
	21	! -----------------
[4592]	22	!
	23	!
[1400]	24	! Former revisions:
	25	! -----------------
	26	! $Id: random_generator_parallel_mod.f90 4828 2021-01-05 11:21:41Z scharf $
[4592]	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	!
[4545]	32	! 4438 2020-03-03 20:49:28Z raasch
[4438]	33	! - Rename variables to avoid confusion with subdomain grid indices
	34	! - Some formatting adjustments
[4592]	35	! - New routine to initialize spatial 1D-arrays independent on the 2D random-number array
	36	!
[4438]	37	! 4360 2020-01-07 11:25:50Z suehring
[4182]	38	! Corrected "Former revisions" section
[4592]	39	!
[4182]	40	! 3655 2019-01-07 16:51:22Z knoop
[4592]	41	! Unused variable removed
[4182]	42	!
	43	! 1400 2014-05-09 14:03:54Z knoop
	44	! Initial revision
[1400]	45	!
[4592]	46	!
	47	!--------------------------------------------------------------------------------------------------!
[1400]	48	! Description:
	49	! ------------
[1682]	50	!> This module contains and supports the random number generating routine ran_parallel.
[4592]	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.
[1682]	59	!>
	60	!> This routine is taken from the "numerical recipies vol. 2"
[4592]	61	!--------------------------------------------------------------------------------------------------!
	62	MODULE random_generator_parallel
[1400]	63
[4592]	64	USE kinds
[4438]	65
[4592]	66	IMPLICIT NONE
[4438]	67
[4592]	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 !<
[4438]	76
[4592]	77	INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds !<
[4438]	78
[4592]	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 !<
[4438]	85
[4592]	86	INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: id_random_array !<
	87	INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE :: seq_random_array !<
[4438]	88
[4592]	89	REAL(wp), SAVE :: amm !<
[4438]	90
[4592]	91	REAL(wp) :: random_dummy = 0.0 !<
[4438]	92
[4592]	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
[4438]	97
[4592]	98	INTERFACE random_number_parallel
	99	MODULE PROCEDURE ran0_s
	100	END INTERFACE
[4438]	101
[4592]	102	INTERFACE random_number_parallel_gauss
	103	MODULE PROCEDURE gasdev_s
	104	END INTERFACE
[4545]	105
[4592]	106	INTERFACE random_seed_parallel
	107	MODULE PROCEDURE random_seed_parallel
	108	END INTERFACE
[4438]	109
[4592]	110	INTERFACE ran_hash
	111	MODULE PROCEDURE ran_hash_v
	112	END INTERFACE
[4438]	113
[4592]	114	INTERFACE reallocate
	115	MODULE PROCEDURE reallocate_iv,reallocate_im
	116	END INTERFACE
[4438]	117
[4592]	118	INTERFACE arth
	119	MODULE PROCEDURE arth_i
	120	END INTERFACE
[1400]	121
[4592]	122	PRIVATE
[4438]	123
[4592]	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
[4438]	131
[1400]	132	CONTAINS
[4438]	133
[4592]	134	!--------------------------------------------------------------------------------------------------!
[1400]	135	! Description:
	136	! ------------
[4438]	137	!> Initialize the parallel random number generator for a 1-dimensional array.
[4592]	138	!--------------------------------------------------------------------------------------------------!
	139	SUBROUTINE init_parallel_random_generator_1d( nxy, ns, ne, id_rand, seq_rand )
[4438]	140
[4592]	141	USE control_parameters, &
	142	ONLY: ensemble_member_nr
[4438]	143
[4592]	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
[4438]	147
[4592]	148	INTEGER(iwp) :: j !< loop index
[4438]	149
[4592]	150	INTEGER(isp), DIMENSION(ns:ne) :: id_rand !< initial IDs
	151	INTEGER(isp), DIMENSION(5,ns:ne) :: seq_rand !< initial random seeds
[4438]	152
	153	!
[4592]	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
[4438]	158	!
[4592]	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
[4438]	166
[4592]	167	END SUBROUTINE init_parallel_random_generator_1d
[4438]	168
[4592]	169	!--------------------------------------------------------------------------------------------------!
[4438]	170	! Description:
	171	! ------------
[4592]	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 )
[2172]	175
[4592]	176	USE kinds
[2172]	177
[4592]	178	USE control_parameters, &
	179	ONLY: ensemble_member_nr
[2172]	180
[4592]	181	IMPLICIT NONE
[2172]	182
[4592]	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
[2172]	188
[4592]	189	INTEGER(iwp) :: i !< grid index x-direction
	190	INTEGER(iwp) :: j !< grid index y-direction
[2172]	191
[4592]	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
[2172]	197
[4592]	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
[2172]	214
[4592]	215	END SUBROUTINE init_parallel_random_generator_2d
	216
	217	!--------------------------------------------------------------------------------------------------!
[2172]	218	! Description:
	219	! ------------
[4592]	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.
[1682]	224	!> The period of this generator is about 2.0 x 10^28, and it fully vectorizes.
	225	!> Validity of the integer model assumed by this generator is tested at initialization.
[4592]	226	!--------------------------------------------------------------------------------------------------!
	227	SUBROUTINE ran0_s( harvest )
[1400]	228
[4592]	229	USE kinds
[1400]	230
[4592]	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	!--------------------------------------------------------------------------------------------------!
[1400]	261	! Description:
	262	! ------------
[4592]	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 )
[4545]	268
[4592]	269	REAL(wp), INTENT(OUT) :: harvest !<
[4545]	270
[4592]	271	REAL(wp) :: rsq !<
	272	REAL(wp) :: v1 !<
	273	REAL(wp) :: v2 !<
	274	REAL(wp), SAVE :: g !<
[4545]	275
[4592]	276	LOGICAL, SAVE :: gaus_stored = .FALSE. !<
[4545]	277	!
[4592]	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.
[4545]	282	!
[4592]	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
[4545]	291	!
[4592]	292	!-- See if they are in the unit circle
	293	rsq = v12 + v22
[4545]	294	!
[4592]	295	!-- Otherwise try again.
	296	IF ( rsq > 0.0 .AND. rsq < 1.0 ) EXIT
	297	ENDDO
[4545]	298	!
[4592]	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
[4545]	306
[4592]	307	END SUBROUTINE gasdev_s
[4545]	308
[4592]	309	!--------------------------------------------------------------------------------------------------!
[4545]	310	! Description:
	311	! ------------
[4592]	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 )
[1400]	317
[4592]	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	!--------------------------------------------------------------------------------------------------!
[1400]	414	! Description:
	415	! ------------
[1682]	416	!> User interface to release the workspace used by the random number routines.
[4592]	417	!--------------------------------------------------------------------------------------------------!
	418	SUBROUTINE ran_deallocate
[1400]	419
[4592]	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	!--------------------------------------------------------------------------------------------------!
[1400]	431	! Description:
	432	! ------------
[4592]	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 )
[1400]	439
[4592]	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	!--------------------------------------------------------------------------------------------------!
[1400]	486	! Description:
	487	! ------------
[4592]	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 )
[1400]	492
[4592]	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	!--------------------------------------------------------------------------------------------------!
[1400]	515	! Description:
	516	! ------------
[4592]	517	!> User interface to process error-messages produced by the random_number_parallel module.
	518	!--------------------------------------------------------------------------------------------------!
	519	SUBROUTINE ran_error( string )
[2255]	520
[4592]	521	USE control_parameters, &
	522	ONLY: message_string
[2255]	523
[4592]	524	CHARACTER(LEN=*), INTENT(IN) :: string !< Error message string
[2255]	525
[4592]	526	message_string = 'incompatible integer arithmetic: ' // TRIM( string )
	527	CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 )
[2255]	528
[4592]	529	END SUBROUTINE ran_error
[1400]	530
[4592]	531	!--------------------------------------------------------------------------------------------------!
[1400]	532	! Description:
	533	! ------------
[4592]	534	!> Reallocates the generators state space "ranseeds" to vectors of size length.
	535	!--------------------------------------------------------------------------------------------------!
	536	FUNCTION reallocate_iv( p, n )
[2255]	537
[4592]	538	USE control_parameters, &
	539	ONLY: message_string
[2255]	540
[4592]	541	INTEGER(isp), DIMENSION(:), POINTER :: p !<
	542	INTEGER(isp), DIMENSION(:), POINTER :: reallocate_iv !<
[2255]	543
[4592]	544	INTEGER(isp), INTENT(IN) :: n !<
[2255]	545
[4592]	546	INTEGER(isp) :: nold !<
	547	INTEGER(isp) :: ierr !<
[1400]	548
[4592]	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	!--------------------------------------------------------------------------------------------------!
[1400]	600	! Description:
	601	! ------------
[4592]	602	!> Reallocates the generators state space "ranseeds" to vectors of size length.
	603	!--------------------------------------------------------------------------------------------------!
	604	FUNCTION arth_i( first, increment, n )
[1400]	605
[4592]	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
	647
[1400]	648	END MODULE random_generator_parallel

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: