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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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