source: palm/trunk/SOURCE/disturb_field.f90

Last change on this file was 4828, checked in by Giersch, 6 months ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

  • Property svn:keywords set to Id
File size: 9.8 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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: disturb_field.f90 4828 2021-01-05 11:21:41Z banzhafs $
26! file re-formatted to follow the PALM coding standard
27!
28! 4457 2020-03-11 14:20:43Z raasch
29! use statement for exchange horiz added
30!
31! 4360 2020-01-07 11:25:50Z suehring
32! Introduction of wall_flags_total_0, which currently sets bits based on static topography
33! information used in wall_flags_static_0
34!
35! 4329 2019-12-10 15:46:36Z motisi
36! Renamed wall_flags_0 to wall_flags_static_0
37!
38! 4237 2019-09-25 11:33:42Z knoop
39! Added missing OpenMP directives
40!
41! 4182 2019-08-22 15:20:23Z scharf
42! Corrected "Former revisions" section
43!
44! 3849 2019-04-01 16:35:16Z knoop
45! Corrected "Former revisions" section
46!
47! Revision 1.1  1998/02/04 15:40:45  raasch
48! Initial revision
49!
50!
51! Description:
52! ------------
53!> Imposing a random perturbation on a 3D-array.
54!> On parallel computers, the random number generator is as well called for all gridpoints of the
55!> total domain to ensure, regardless of the number of PEs used, that the elements of the array have
56!> the same values in the same order in every case. The perturbation range is steered by dist_range.
57!--------------------------------------------------------------------------------------------------!
58 SUBROUTINE disturb_field( var_char, dist1, field )
59
60
61    USE control_parameters,                                                                        &
62        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range, disturbance_amplitude,          &
63               disturbance_created, disturbance_level_ind_b, disturbance_level_ind_t, iran,        &
64               random_generator, topography
65
66    USE cpulog,                                                                                    &
67        ONLY:  cpu_log, log_point
68
69    USE exchange_horiz_mod,                                                                        &
70        ONLY:  exchange_horiz
71
72    USE indices,                                                                                   &
73        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, nzt,                &
74               wall_flags_total_0
75
76    USE kinds
77
78    USE random_function_mod,                                                                       &
79        ONLY: random_function
80
81    USE random_generator_parallel,                                                                 &
82        ONLY:  random_number_parallel, random_seed_parallel, random_dummy, seq_random_array
83
84    IMPLICIT NONE
85
86    CHARACTER (LEN = *) ::  var_char !< flag to distinguish betwenn u- and v-component
87
88    INTEGER(iwp) ::  flag_nr !< number of respective flag for u- or v-grid
89    INTEGER(iwp) ::  i       !< index variable
90    INTEGER(iwp) ::  j       !< index variable
91    INTEGER(iwp) ::  k       !< index variable
92
93    REAL(wp) ::  randomnumber  !<
94
95    REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
96    REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
97
98    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !<
99
100
101    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
102!
103!-- Set flag number, 20 for u-grid, 21 for v-grid, required to mask topography
104    flag_nr = MERGE( 20, 21, TRIM(var_char) == 'u' )
105!
106!-- Create an additional temporary array and initialize the arrays needed to store the disturbance
107    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
108!$ACC DATA CREATE(dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
109
110!
111!-- dist1 is initialized on the host (see below) and then updated on the device.
112    dist1 = 0.0_wp
113    !$ACC KERNELS PRESENT(dist2)
114    dist2 = 0.0_wp
115    !$ACC END KERNELS
116
117!
118!-- Create the random perturbation and store it on temporary array
119    IF ( random_generator == 'numerical-recipes' )  THEN
120       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
121          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
122             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
123                randomnumber = 3.0_wp * disturbance_amplitude * ( random_function( iran ) - 0.5_wp )
124                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND. nyn >= j )  THEN
125                   dist1(k,j,i) = randomnumber
126                ENDIF
127             ENDDO
128          ENDDO
129       ENDDO
130    ELSEIF ( random_generator == 'random-parallel' )  THEN
131       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
132          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
133             CALL random_seed_parallel( put=seq_random_array(:, j, i) )
134             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
135                CALL random_number_parallel( random_dummy )
136                randomnumber = 3.0_wp * disturbance_amplitude * ( random_dummy - 0.5_wp )
137                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  nyn >= j )  THEN
138                   dist1(k,j,i) = randomnumber
139                ENDIF
140             ENDDO
141             CALL random_seed_parallel( get=seq_random_array(:, j, i) )
142          ENDDO
143       ENDDO
144    ELSEIF ( random_generator == 'system-specific' )  THEN
145       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
146          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
147             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
148                CALL RANDOM_NUMBER( randomnumber )
149                randomnumber = 3.0_wp * disturbance_amplitude * ( randomnumber - 0.5_wp )
150                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  nyn >= j )  THEN
151                   dist1(k,j,i) = randomnumber
152                ENDIF
153             ENDDO
154          ENDDO
155       ENDDO
156
157    ENDIF
158
159!
160!-- Update dist1 on the device, this is expected by exchange_horiz!
161    !$ACC UPDATE DEVICE(dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
162
163!
164!-- Exchange of ghost points for the random perturbation
165
166    CALL exchange_horiz( dist1, nbgp )
167!
168!-- Applying the Shuman filter in order to smooth the perturbations.
169!-- Neighboured grid points in all three directions are used for the filter operation.
170!-- Loop has been splitted to make runs reproducible on HLRN systems using compiler option -O3
171     !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2)
172     !$OMP PARALLEL DO PRIVATE(i, j, k)
173     DO  i = nxl, nxr
174        DO  j = nys, nyn
175          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
176             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                                      &
177                            + dist1(k,j+1,i) + dist1(k+1,j,i)                                      &
178                            ) / 12.0_wp
179          ENDDO
180          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
181              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)                      &
182                            + 6.0_wp * dist1(k,j,i)                                                &
183                                            ) / 12.0_wp
184          ENDDO
185        ENDDO
186     ENDDO
187
188!
189!-- Exchange of ghost points for the filtered perturbation.
190!-- Afterwards, filter operation and exchange of ghost points are repeated.
191    CALL exchange_horiz( dist2, nbgp )
192
193    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) PRESENT(dist1, dist2)
194    !$OMP PARALLEL DO PRIVATE(i, j, k)
195    DO  i = nxl, nxr
196       DO  j = nys, nyn
197          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
198             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i)                     &
199                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i)                     &
200                            + 6.0_wp * dist2(k,j,i)                                                &
201                            ) / 12.0_wp
202          ENDDO
203       ENDDO
204    ENDDO
205
206    CALL exchange_horiz( dist1, nbgp )
207
208!
209!-- Remove perturbations below topography (including one gridpoint above it in order to allow for
210!-- larger timesteps at the beginning of the simulation (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_total_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.