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

Last change on this file since 1561 was 1426, checked in by knoop, 10 years ago

last commit documented

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