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

Last change on this file since 2233 was 2233, checked in by suehring, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 9.8 KB
RevLine 
[1682]1!> @file disturb_field.f90
[2000]2!------------------------------------------------------------------------------!
[1036]3! This file is part of PALM.
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!
[2101]17! Copyright 1997-2017 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 2233 2017-05-30 18:08:54Z suehring $
27!
[2233]28! 2232 2017-05-30 17:47:52Z suehring
29! Modify referenced parameter for disturb_field, instead of nzb_uv_inner, pass
30! character to identify the respective grid (u- or v-grid).
31! Set perturbations within topography to zero using flags.
32!
[2173]33! 2172 2017-03-08 15:55:25Z knoop
34! Bugfix removed id_random_array from USE list
35!
[2001]36! 2000 2016-08-20 18:09:15Z knoop
37! Forced header and separation lines into 80 columns
38!
[1683]39! 1682 2015-10-07 23:56:08Z knoop
40! Code annotations made doxygen readable
41!
[1426]42! 1425 2014-07-05 10:57:53Z knoop
43! bugfix: Parallel random number generator loop: print-statement no needed
44!
[1401]45! 1400 2014-05-09 14:03:54Z knoop
46! Parallel random number generator added
47!
[1354]48! 1353 2014-04-08 15:21:23Z heinze
49! REAL constants provided with KIND-attribute
50!
[1323]51! 1322 2014-03-20 16:38:49Z raasch
52! REAL constants defined as wp-kind
53!
[1321]54! 1320 2014-03-20 08:40:49Z raasch
[1320]55! ONLY-attribute added to USE-statements,
56! kind-parameters added to all INTEGER and REAL declaration statements,
57! kinds are defined in new module kinds,
58! revision history before 2012 removed,
59! comment fields (!:) to be used for variable explanations added to
60! all variable declaration statements
[1]61!
[1037]62! 1036 2012-10-22 13:43:42Z raasch
63! code put under GPL (PALM 3.9)
64!
[1]65! Revision 1.1  1998/02/04 15:40:45  raasch
66! Initial revision
67!
68!
69! Description:
70! ------------
[1682]71!> Imposing a random perturbation on a 3D-array.
72!> On parallel computers, the random number generator is as well called for all
73!> gridpoints of the total domain to ensure, regardless of the number of PEs
74!> used, that the elements of the array have the same values in the same
75!> order in every case. The perturbation range is steered by dist_range.
[1]76!------------------------------------------------------------------------------!
[2232]77 SUBROUTINE disturb_field( var_char, dist1, field )
[1682]78 
[1]79
[2232]80    USE control_parameters,                                                    &
[1320]81        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range,             &
82               disturbance_amplitude, disturbance_created,                     &
83               disturbance_level_ind_b, disturbance_level_ind_t, iran,         &
84               random_generator, topography
85               
86    USE cpulog,                                                                &
87        ONLY:  cpu_log, log_point
88       
89    USE indices,                                                               &
[2232]90        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
91               nzt, wall_flags_0
[1320]92       
93    USE kinds
94   
95    USE random_function_mod,                                                   &
96        ONLY: random_function
[1400]97       
98    USE random_generator_parallel,                                             &
99        ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
[2172]100               seq_random_array
[1]101
102    IMPLICIT NONE
103
[2232]104    CHARACTER (LEN = *) ::  var_char !< flag to distinguish betwenn u- and v-component
[1]105
[2232]106    INTEGER(iwp) ::  flag_nr !< number of respective flag for u- or v-grid
107    INTEGER(iwp) ::  i       !< index variable
108    INTEGER(iwp) ::  j       !< index variable
109    INTEGER(iwp) ::  k       !< index variable
110
[1682]111    REAL(wp) ::  randomnumber  !<
[1320]112   
[1682]113    REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
114    REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
[1320]115   
[1682]116    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !<
[1]117
118
119    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
120!
[2232]121!-- Set flag number, 20 for u-grid, 21 for v-grid, required to mask topography
122    flag_nr = MERGE( 20, 21, TRIM(var_char) == 'u' )
123!
[1]124!-- Create an additional temporary array and initialize the arrays needed
125!-- to store the disturbance
[667]126    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1353]127    dist1 = 0.0_wp
128    dist2 = 0.0_wp
[1]129
130!
131!-- Create the random perturbation and store it on temporary array
132    IF ( random_generator == 'numerical-recipes' )  THEN
133       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
134          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
135             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
[1322]136                randomnumber = 3.0_wp * disturbance_amplitude *                &
137                               ( random_function( iran ) - 0.5_wp )
[1320]138                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
139                     nyn >= j )                                                &
[1]140                THEN
141                   dist1(k,j,i) = randomnumber
142                ENDIF
143             ENDDO
144          ENDDO
145       ENDDO
[1400]146    ELSEIF ( random_generator == 'random-parallel' )  THEN
147       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
148          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
149             CALL random_seed_parallel( put=seq_random_array(:, j, i) )
150             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
151                CALL random_number_parallel( random_dummy )
152                randomnumber = 3.0_wp * disturbance_amplitude *                &
153                               ( random_dummy - 0.5_wp )
154                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
155                     nyn >= j )                                                &
156                THEN
157                   dist1(k,j,i) = randomnumber
158                ENDIF
159             ENDDO
160             CALL random_seed_parallel( get=seq_random_array(:, j, i) )
161          ENDDO
162       ENDDO
[1]163    ELSEIF ( random_generator == 'system-specific' )  THEN
164       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
165          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
166             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
167#if defined( __nec )
[1322]168                randomnumber = 3.0_wp * disturbance_amplitude *                &
169                               ( RANDOM( 0 ) - 0.5_wp )
[1]170#else
171                CALL RANDOM_NUMBER( randomnumber )
[1322]172                randomnumber = 3.0_wp * disturbance_amplitude *                &
173                                ( randomnumber - 0.5_wp )
[1]174#endif
[1320]175                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
[1]176                THEN
177                   dist1(k,j,i) = randomnumber
178                ENDIF
179             ENDDO
180          ENDDO
181       ENDDO
182
183    ENDIF
184
185!
186!-- Exchange of ghost points for the random perturbation
[420]187
[667]188    CALL exchange_horiz( dist1, nbgp )
[1]189!
190!-- Applying the Shuman filter in order to smooth the perturbations.
191!-- Neighboured grid points in all three directions are used for the
192!-- filter operation.
[420]193!-- Loop has been splitted to make runs reproducible on HLRN systems using
194!-- compiler option -O3
195     DO  i = nxl, nxr
196        DO  j = nys, nyn
[1]197          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
[1320]198             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
199                            + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
[1322]200                            ) / 12.0_wp
[420]201          ENDDO
202          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
203              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
[1322]204                            + 6.0_wp * dist1(k,j,i)                            &
205                            ) / 12.0_wp
[1]206          ENDDO
[420]207        ENDDO
208     ENDDO
[1]209
210!
211!-- Exchange of ghost points for the filtered perturbation.
212!-- Afterwards, filter operation and exchange of ghost points are repeated.
[667]213    CALL exchange_horiz( dist2, nbgp )
[420]214
[1]215    DO  i = nxl, nxr
216       DO  j = nys, nyn
217          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
218             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
219                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
[1322]220                            + 6.0_wp * dist2(k,j,i)                            &
221                            ) / 12.0_wp
[1]222          ENDDO
223       ENDDO
224    ENDDO
[420]225
[667]226    CALL exchange_horiz( dist1, nbgp )
[1]227
228!
229!-- Remove perturbations below topography (including one gridpoint above it
230!-- in order to allow for larger timesteps at the beginning of the simulation
231!-- (diffusion criterion))
232    IF ( TRIM( topography ) /= 'flat' )  THEN
[667]233       DO  i = nxlg, nxrg
234          DO  j = nysg, nyng
[2232]235             DO  k = nzb, nzb_max
236                dist1(k,j,i) = MERGE( dist1(k,j,i), 0.0_wp,                    &
237                                      BTEST( wall_flags_0(k,j,i), flag_nr )    &
238                                    )
239             ENDDO
[1]240          ENDDO
241       ENDDO
242    ENDIF
243
244!
245!-- Random perturbation is added to the array to be disturbed.
[667]246    DO  i = nxlg, nxrg
247       DO  j = nysg, nyng
[1]248          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
249             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
250          ENDDO
251       ENDDO
252    ENDDO
253
254!
255!-- Deallocate the temporary array
256    DEALLOCATE( dist2 )
257
258!
259!-- Set a flag, which indicates that a random perturbation is imposed
260    disturbance_created = .TRUE.
261
262
263    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )
264
265
266 END SUBROUTINE disturb_field
Note: See TracBrowser for help on using the repository browser.