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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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