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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

  • Property svn:keywords set to Id
File size: 8.9 KB
Line 
1!> @file disturb_field.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2014 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21! Code annotations made doxygen readable
22!
23! Former revisions:
24! -----------------
25! $Id: disturb_field.f90 1682 2015-10-07 23:56:08Z knoop $
26!
27! 1425 2014-07-05 10:57:53Z knoop
28! bugfix: Parallel random number generator loop: print-statement no needed
29!
30! 1400 2014-05-09 14:03:54Z knoop
31! Parallel random number generator added
32!
33! 1353 2014-04-08 15:21:23Z heinze
34! REAL constants provided with KIND-attribute
35!
36! 1322 2014-03-20 16:38:49Z raasch
37! REAL constants defined as wp-kind
38!
39! 1320 2014-03-20 08:40:49Z raasch
40! ONLY-attribute added to USE-statements,
41! kind-parameters added to all INTEGER and REAL declaration statements,
42! kinds are defined in new module kinds,
43! revision history before 2012 removed,
44! comment fields (!:) to be used for variable explanations added to
45! all variable declaration statements
46!
47! 1036 2012-10-22 13:43:42Z raasch
48! code put under GPL (PALM 3.9)
49!
50! Revision 1.1  1998/02/04 15:40:45  raasch
51! Initial revision
52!
53!
54! Description:
55! ------------
56!> Imposing a random perturbation on a 3D-array.
57!> On parallel computers, the random number generator is as well called for all
58!> gridpoints of the total domain to ensure, regardless of the number of PEs
59!> used, that the elements of the array have the same values in the same
60!> order in every case. The perturbation range is steered by dist_range.
61!------------------------------------------------------------------------------!
62 SUBROUTINE disturb_field( nzb_uv_inner, dist1, field )
63 
64
65    USE control_parameters,   &
66        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range,             &
67               disturbance_amplitude, disturbance_created,                     &
68               disturbance_level_ind_b, disturbance_level_ind_t, iran,         &
69               random_generator, topography
70               
71    USE cpulog,                                                                &
72        ONLY:  cpu_log, log_point
73       
74    USE indices,                                                               &
75        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
76       
77    USE kinds
78   
79    USE random_function_mod,                                                   &
80        ONLY: random_function
81       
82    USE random_generator_parallel,                                             &
83        ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
84               id_random_array, seq_random_array
85
86    IMPLICIT NONE
87
88    INTEGER(iwp) ::  i  !<
89    INTEGER(iwp) ::  j  !<
90    INTEGER(iwp) ::  k  !<
91   
92    INTEGER(iwp) ::  nzb_uv_inner(nysg:nyng,nxlg:nxrg) !<
93
94    REAL(wp) ::  randomnumber  !<
95   
96    REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
97    REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
98   
99    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !<
100
101
102    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
103
104!
105!-- Create an additional temporary array and initialize the arrays needed
106!-- to store the disturbance
107    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
108    dist1 = 0.0_wp
109    dist2 = 0.0_wp
110
111!
112!-- Create the random perturbation and store it on temporary array
113    IF ( random_generator == 'numerical-recipes' )  THEN
114       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
115          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
116             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
117                randomnumber = 3.0_wp * disturbance_amplitude *                &
118                               ( random_function( iran ) - 0.5_wp )
119                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
120                     nyn >= j )                                                &
121                THEN
122                   dist1(k,j,i) = randomnumber
123                ENDIF
124             ENDDO
125          ENDDO
126       ENDDO
127    ELSEIF ( random_generator == 'random-parallel' )  THEN
128       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
129          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
130             CALL random_seed_parallel( put=seq_random_array(:, j, i) )
131             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
132                CALL random_number_parallel( random_dummy )
133                randomnumber = 3.0_wp * disturbance_amplitude *                &
134                               ( random_dummy - 0.5_wp )
135                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
136                     nyn >= j )                                                &
137                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#if defined( __nec )
149                randomnumber = 3.0_wp * disturbance_amplitude *                &
150                               ( RANDOM( 0 ) - 0.5_wp )
151#else
152                CALL RANDOM_NUMBER( randomnumber )
153                randomnumber = 3.0_wp * disturbance_amplitude *                &
154                                ( randomnumber - 0.5_wp )
155#endif
156                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
157                THEN
158                   dist1(k,j,i) = randomnumber
159                ENDIF
160             ENDDO
161          ENDDO
162       ENDDO
163
164    ENDIF
165
166!
167!-- Exchange of ghost points for the random perturbation
168
169    CALL exchange_horiz( dist1, nbgp )
170!
171!-- Applying the Shuman filter in order to smooth the perturbations.
172!-- Neighboured grid points in all three directions are used for the
173!-- filter operation.
174!-- Loop has been splitted to make runs reproducible on HLRN systems using
175!-- compiler option -O3
176     DO  i = nxl, nxr
177        DO  j = nys, nyn
178          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
179             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
180                            + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
181                            ) / 12.0_wp
182          ENDDO
183          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
184              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
185                            + 6.0_wp * dist1(k,j,i)                            &
186                            ) / 12.0_wp
187          ENDDO
188        ENDDO
189     ENDDO
190
191!
192!-- Exchange of ghost points for the filtered perturbation.
193!-- Afterwards, filter operation and exchange of ghost points are repeated.
194    CALL exchange_horiz( dist2, nbgp )
195
196    DO  i = nxl, nxr
197       DO  j = nys, nyn
198          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
199             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
200                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
201                            + 6.0_wp * dist2(k,j,i)                            &
202                            ) / 12.0_wp
203          ENDDO
204       ENDDO
205    ENDDO
206
207    CALL exchange_horiz( dist1, nbgp )
208
209!
210!-- Remove perturbations below topography (including one gridpoint above it
211!-- in order to allow for larger timesteps at the beginning of the simulation
212!-- (diffusion criterion))
213    IF ( TRIM( topography ) /= 'flat' )  THEN
214       DO  i = nxlg, nxrg
215          DO  j = nysg, nyng
216             dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0_wp
217          ENDDO
218       ENDDO
219    ENDIF
220
221!
222!-- Random perturbation is added to the array to be disturbed.
223    DO  i = nxlg, nxrg
224       DO  j = nysg, nyng
225          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
226             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
227          ENDDO
228       ENDDO
229    ENDDO
230
231!
232!-- Deallocate the temporary array
233    DEALLOCATE( dist2 )
234
235!
236!-- Set a flag, which indicates that a random perturbation is imposed
237    disturbance_created = .TRUE.
238
239
240    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )
241
242
243 END SUBROUTINE disturb_field
Note: See TracBrowser for help on using the repository browser.