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

Last change on this file since 2173 was 2173, checked in by knoop, 4 years ago

last commit documented

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