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

Last change on this file since 4338 was 4329, checked in by motisi, 5 years ago

Renamed wall_flags_0 to wall_flags_static_0

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