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

Last change on this file since 4178 was 3849, checked in by knoop, 2 years ago

Bugfix: added proper OpenACC support to disturb_field

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