source: palm/trunk/SOURCE/disturb_field.f90 @ 1425

Last change on this file since 1425 was 1425, checked in by knoop, 7 years ago

Bugfix for parallel random number generator.

  • Property svn:keywords set to Id
File size: 8.7 KB
Line 
1 SUBROUTINE disturb_field( nzb_uv_inner, dist1, field )
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! bugfix: Parallel random number generator loop: print-statement no needed
23!
24! Former revisions:
25! -----------------
26! $Id: disturb_field.f90 1425 2014-07-05 10:57:53Z knoop $
27!
28! 1400 2014-05-09 14:03:54Z knoop
29! Parallel random number generator added
30!
31! 1353 2014-04-08 15:21:23Z heinze
32! REAL constants provided with KIND-attribute
33!
34! 1322 2014-03-20 16:38:49Z raasch
35! REAL constants defined as wp-kind
36!
37! 1320 2014-03-20 08:40:49Z raasch
38! ONLY-attribute added to USE-statements,
39! kind-parameters added to all INTEGER and REAL declaration statements,
40! kinds are defined in new module kinds,
41! revision history before 2012 removed,
42! comment fields (!:) to be used for variable explanations added to
43! all variable declaration statements
44!
45! 1036 2012-10-22 13:43:42Z raasch
46! code put under GPL (PALM 3.9)
47!
48! Revision 1.1  1998/02/04 15:40:45  raasch
49! Initial revision
50!
51!
52! Description:
53! ------------
54! Imposing a random perturbation on a 3D-array.
55! On parallel computers, the random number generator is as well called for all
56! gridpoints of the total domain to ensure, regardless of the number of PEs
57! used, that the elements of the array have the same values in the same
58! order in every case. The perturbation range is steered by dist_range.
59!------------------------------------------------------------------------------!
60
61    USE control_parameters,   &
62        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range,             &
63               disturbance_amplitude, disturbance_created,                     &
64               disturbance_level_ind_b, disturbance_level_ind_t, iran,         &
65               random_generator, topography
66               
67    USE cpulog,                                                                &
68        ONLY:  cpu_log, log_point
69       
70    USE indices,                                                               &
71        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
72       
73    USE kinds
74   
75    USE random_function_mod,                                                   &
76        ONLY: random_function
77       
78    USE random_generator_parallel,                                             &
79        ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
80               id_random_array, seq_random_array
81
82    IMPLICIT NONE
83
84    INTEGER(iwp) ::  i  !:
85    INTEGER(iwp) ::  j  !:
86    INTEGER(iwp) ::  k  !:
87   
88    INTEGER(iwp) ::  nzb_uv_inner(nysg:nyng,nxlg:nxrg) !:
89
90    REAL(wp) ::  randomnumber  !:
91   
92    REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !:
93    REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !:
94   
95    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !:
96
97
98    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
99
100!
101!-- Create an additional temporary array and initialize the arrays needed
102!-- to store the disturbance
103    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
104    dist1 = 0.0_wp
105    dist2 = 0.0_wp
106
107!
108!-- Create the random perturbation and store it on temporary array
109    IF ( random_generator == 'numerical-recipes' )  THEN
110       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
111          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
112             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
113                randomnumber = 3.0_wp * disturbance_amplitude *                &
114                               ( random_function( iran ) - 0.5_wp )
115                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
116                     nyn >= j )                                                &
117                THEN
118                   dist1(k,j,i) = randomnumber
119                ENDIF
120             ENDDO
121          ENDDO
122       ENDDO
123    ELSEIF ( random_generator == 'random-parallel' )  THEN
124       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
125          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
126             CALL random_seed_parallel( put=seq_random_array(:, j, i) )
127             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
128                CALL random_number_parallel( random_dummy )
129                randomnumber = 3.0_wp * disturbance_amplitude *                &
130                               ( random_dummy - 0.5_wp )
131                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
132                     nyn >= j )                                                &
133                THEN
134                   dist1(k,j,i) = randomnumber
135                ENDIF
136             ENDDO
137             CALL random_seed_parallel( get=seq_random_array(:, j, i) )
138          ENDDO
139       ENDDO
140    ELSEIF ( random_generator == 'system-specific' )  THEN
141       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
142          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
143             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
144#if defined( __nec )
145                randomnumber = 3.0_wp * disturbance_amplitude *                &
146                               ( RANDOM( 0 ) - 0.5_wp )
147#else
148                CALL RANDOM_NUMBER( randomnumber )
149                randomnumber = 3.0_wp * disturbance_amplitude *                &
150                                ( randomnumber - 0.5_wp )
151#endif
152                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
153                THEN
154                   dist1(k,j,i) = randomnumber
155                ENDIF
156             ENDDO
157          ENDDO
158       ENDDO
159
160    ENDIF
161
162!
163!-- Exchange of ghost points for the random perturbation
164
165    CALL exchange_horiz( dist1, nbgp )
166!
167!-- Applying the Shuman filter in order to smooth the perturbations.
168!-- Neighboured grid points in all three directions are used for the
169!-- filter operation.
170!-- Loop has been splitted to make runs reproducible on HLRN systems using
171!-- compiler option -O3
172     DO  i = nxl, nxr
173        DO  j = nys, nyn
174          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
175             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
176                            + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
177                            ) / 12.0_wp
178          ENDDO
179          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
180              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
181                            + 6.0_wp * dist1(k,j,i)                            &
182                            ) / 12.0_wp
183          ENDDO
184        ENDDO
185     ENDDO
186
187!
188!-- Exchange of ghost points for the filtered perturbation.
189!-- Afterwards, filter operation and exchange of ghost points are repeated.
190    CALL exchange_horiz( dist2, nbgp )
191
192    DO  i = nxl, nxr
193       DO  j = nys, nyn
194          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
195             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
196                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
197                            + 6.0_wp * dist2(k,j,i)                            &
198                            ) / 12.0_wp
199          ENDDO
200       ENDDO
201    ENDDO
202
203    CALL exchange_horiz( dist1, nbgp )
204
205!
206!-- Remove perturbations below topography (including one gridpoint above it
207!-- in order to allow for larger timesteps at the beginning of the simulation
208!-- (diffusion criterion))
209    IF ( TRIM( topography ) /= 'flat' )  THEN
210       DO  i = nxlg, nxrg
211          DO  j = nysg, nyng
212             dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0_wp
213          ENDDO
214       ENDDO
215    ENDIF
216
217!
218!-- Random perturbation is added to the array to be disturbed.
219    DO  i = nxlg, nxrg
220       DO  j = nysg, nyng
221          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
222             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
223          ENDDO
224       ENDDO
225    ENDDO
226
227!
228!-- Deallocate the temporary array
229    DEALLOCATE( dist2 )
230
231!
232!-- Set a flag, which indicates that a random perturbation is imposed
233    disturbance_created = .TRUE.
234
235
236    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )
237
238
239 END SUBROUTINE disturb_field
Note: See TracBrowser for help on using the repository browser.