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

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

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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