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

Last change on this file since 4251 was 4237, checked in by knoop, 5 years ago

Added missing OpenMP directives within "disturb_field", "surface_layer_fluxes" and "timestep"

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