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

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

Parallel random number generator added (preliminary version).

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