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

Last change on this file since 4185 was 4182, checked in by scharf, 5 years ago
  • corrected "Former revisions" section
  • minor formatting in "Former revisions" section
  • added "Author" section
  • Property svn:keywords set to Id
File size: 9.1 KB
Line 
1!> @file disturb_field.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: disturb_field.f90 4182 2019-08-22 15:20:23Z oliver.maas $
27! Corrected "Former revisions" section
28!
29! 3849 2019-04-01 16:35:16Z knoop
30! Corrected "Former revisions" section
31!
32! Revision 1.1  1998/02/04 15:40:45  raasch
33! Initial revision
34!
35!
36! Description:
37! ------------
38!> Imposing a random perturbation on a 3D-array.
39!> On parallel computers, the random number generator is as well called for all
40!> gridpoints of the total domain to ensure, regardless of the number of PEs
41!> used, that the elements of the array have the same values in the same
42!> order in every case. The perturbation range is steered by dist_range.
43!------------------------------------------------------------------------------!
44 SUBROUTINE disturb_field( var_char, dist1, field )
45 
46
47    USE control_parameters,                                                    &
48        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range,             &
49               disturbance_amplitude, disturbance_created,                     &
50               disturbance_level_ind_b, disturbance_level_ind_t, iran,         &
51               random_generator, topography
52               
53    USE cpulog,                                                                &
54        ONLY:  cpu_log, log_point
55       
56    USE indices,                                                               &
57        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
58               nzt, wall_flags_0
59       
60    USE kinds
61   
62    USE random_function_mod,                                                   &
63        ONLY: random_function
64       
65    USE random_generator_parallel,                                             &
66        ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
67               seq_random_array
68
69    IMPLICIT NONE
70
71    CHARACTER (LEN = *) ::  var_char !< flag to distinguish betwenn u- and v-component
72
73    INTEGER(iwp) ::  flag_nr !< number of respective flag for u- or v-grid
74    INTEGER(iwp) ::  i       !< index variable
75    INTEGER(iwp) ::  j       !< index variable
76    INTEGER(iwp) ::  k       !< index variable
77
78    REAL(wp) ::  randomnumber  !<
79   
80    REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
81    REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
82   
83    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !<
84
85
86    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
87!
88!-- Set flag number, 20 for u-grid, 21 for v-grid, required to mask topography
89    flag_nr = MERGE( 20, 21, TRIM(var_char) == 'u' )
90!
91!-- Create an additional temporary array and initialize the arrays needed
92!-- to store the disturbance
93    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
94!$ACC DATA CREATE(dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
95
96!
97!-- dist1 is initialized on the host (see below) and then updated on the device.
98    dist1 = 0.0_wp
99    !$ACC KERNELS PRESENT(dist2)
100    dist2 = 0.0_wp
101    !$ACC END KERNELS
102
103!
104!-- Create the random perturbation and store it on temporary array
105    IF ( random_generator == 'numerical-recipes' )  THEN
106       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
107          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
108             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
109                randomnumber = 3.0_wp * disturbance_amplitude *                &
110                               ( random_function( iran ) - 0.5_wp )
111                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
112                     nyn >= j )                                                &
113                THEN
114                   dist1(k,j,i) = randomnumber
115                ENDIF
116             ENDDO
117          ENDDO
118       ENDDO
119    ELSEIF ( random_generator == 'random-parallel' )  THEN
120       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
121          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
122             CALL random_seed_parallel( put=seq_random_array(:, j, i) )
123             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
124                CALL random_number_parallel( random_dummy )
125                randomnumber = 3.0_wp * disturbance_amplitude *                &
126                               ( random_dummy - 0.5_wp )
127                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
128                     nyn >= j )                                                &
129                THEN
130                   dist1(k,j,i) = randomnumber
131                ENDIF
132             ENDDO
133             CALL random_seed_parallel( get=seq_random_array(:, j, i) )
134          ENDDO
135       ENDDO
136    ELSEIF ( random_generator == 'system-specific' )  THEN
137       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
138          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
139             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
140                CALL RANDOM_NUMBER( randomnumber )
141                randomnumber = 3.0_wp * disturbance_amplitude *                &
142                                ( randomnumber - 0.5_wp )
143                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
144                THEN
145                   dist1(k,j,i) = randomnumber
146                ENDIF
147             ENDDO
148          ENDDO
149       ENDDO
150
151    ENDIF
152
153!
154!-- Update dist1 on the device, this is expected by exchange_horiz!
155    !$ACC UPDATE DEVICE(dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
156
157!
158!-- Exchange of ghost points for the random perturbation
159
160    CALL exchange_horiz( dist1, nbgp )
161!
162!-- Applying the Shuman filter in order to smooth the perturbations.
163!-- Neighboured grid points in all three directions are used for the
164!-- filter operation.
165!-- Loop has been splitted to make runs reproducible on HLRN systems using
166!-- compiler option -O3
167     !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2)
168     DO  i = nxl, nxr
169        DO  j = nys, nyn
170          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
171             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
172                            + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
173                            ) / 12.0_wp
174          ENDDO
175          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
176              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
177                            + 6.0_wp * dist1(k,j,i)                            &
178                            ) / 12.0_wp
179          ENDDO
180        ENDDO
181     ENDDO
182
183!
184!-- Exchange of ghost points for the filtered perturbation.
185!-- Afterwards, filter operation and exchange of ghost points are repeated.
186    CALL exchange_horiz( dist2, nbgp )
187
188    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2)
189    DO  i = nxl, nxr
190       DO  j = nys, nyn
191          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
192             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
193                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
194                            + 6.0_wp * dist2(k,j,i)                            &
195                            ) / 12.0_wp
196          ENDDO
197       ENDDO
198    ENDDO
199
200    CALL exchange_horiz( dist1, nbgp )
201
202!
203!-- Remove perturbations below topography (including one gridpoint above it
204!-- in order to allow for larger timesteps at the beginning of the simulation
205!-- (diffusion criterion))
206    IF ( TRIM( topography ) /= 'flat' )  THEN
207       DO  i = nxlg, nxrg
208          DO  j = nysg, nyng
209             DO  k = nzb, nzb_max
210                dist1(k,j,i) = MERGE( dist1(k,j,i), 0.0_wp,                    &
211                                      BTEST( wall_flags_0(k,j,i), flag_nr )    &
212                                    )
213             ENDDO
214          ENDDO
215       ENDDO
216    ENDIF
217
218!
219!-- Random perturbation is added to the array to be disturbed.
220    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) PRESENT(field, dist1)
221    DO  i = nxlg, nxrg
222       DO  j = nysg, nyng
223          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
224             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
225          ENDDO
226       ENDDO
227    ENDDO
228
229!$ACC END DATA
230
231!
232!-- Deallocate the temporary array
233    DEALLOCATE( dist2 )
234
235!
236!-- Set a flag, which indicates that a random perturbation is imposed
237    disturbance_created = .TRUE.
238
239
240    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )
241
242
243 END SUBROUTINE disturb_field
Note: See TracBrowser for help on using the repository browser.