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

Last change on this file since 1319 was 1319, checked in by raasch, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 6.7 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 1319 2014-03-17 15:08:44Z raasch $
27!
28! 1318 2014-03-17 13:35:16Z raasch
29! module interfaces removed
30!
31! 1036 2012-10-22 13:43:42Z raasch
32! code put under GPL (PALM 3.9)
33!
34! 667 2010-12-23 12:06:00Z suehring/gryschka
35! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
36! Calls of exchange_horiz are modified.
37!
38! 420 2010-01-13 15:10:53Z franke
39! Loop was split to make runs reproducible when using ifort compiler
40!
41! 75 2007-03-22 09:54:05Z raasch
42! xrp, ynp eliminated, 2nd+3rd argument removed from exchange horiz
43!
44! RCS Log replace by Id keyword, revision history cleaned up
45!
46! Revision 1.11  2006/08/04 14:31:59  raasch
47! izuf renamed iran
48!
49! Revision 1.1  1998/02/04 15:40:45  raasch
50! Initial revision
51!
52!
53! Description:
54! ------------
55! Imposing a random perturbation on a 3D-array.
56! On parallel computers, the random number generator is as well called for all
57! gridpoints of the total domain to ensure, regardless of the number of PEs
58! used, that the elements of the array have the same values in the same
59! order in every case. The perturbation range is steered by dist_range.
60!------------------------------------------------------------------------------!
61
62    USE control_parameters
63    USE cpulog
64    USE grid_variables
65    USE indices
66    USE random_function_mod
67
68    IMPLICIT NONE
69
70    INTEGER ::  i, j, k
71    INTEGER ::  nzb_uv_inner(nysg:nyng,nxlg:nxrg)
72
73    REAL    ::  randomnumber,                             &
74                dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
75                field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
76    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  dist2
77
78
79    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
80
81!
82!-- Create an additional temporary array and initialize the arrays needed
83!-- to store the disturbance
84    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
85    dist1 = 0.0
86    dist2 = 0.0
87
88!
89!-- Create the random perturbation and store it on temporary array
90    IF ( random_generator == 'numerical-recipes' )  THEN
91       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
92          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
93             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
94                randomnumber = 3.0 * disturbance_amplitude * &
95                               ( random_function( iran ) - 0.5 )
96                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  &
97                     nyn >= j ) &
98                THEN
99                   dist1(k,j,i) = randomnumber
100                ENDIF
101             ENDDO
102          ENDDO
103       ENDDO
104    ELSEIF ( random_generator == 'system-specific' )  THEN
105       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
106          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
107             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
108#if defined( __nec )
109                randomnumber = 3.0 * disturbance_amplitude * &
110                               ( RANDOM( 0 ) - 0.5 )
111#else
112                CALL RANDOM_NUMBER( randomnumber )
113                randomnumber = 3.0 * disturbance_amplitude * &
114                                ( randomnumber - 0.5 )
115#endif
116                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j ) &
117                THEN
118                   dist1(k,j,i) = randomnumber
119                ENDIF
120             ENDDO
121          ENDDO
122       ENDDO
123
124    ENDIF
125
126!
127!-- Exchange of ghost points for the random perturbation
128
129    CALL exchange_horiz( dist1, nbgp )
130!
131!-- Applying the Shuman filter in order to smooth the perturbations.
132!-- Neighboured grid points in all three directions are used for the
133!-- filter operation.
134!-- Loop has been splitted to make runs reproducible on HLRN systems using
135!-- compiler option -O3
136     DO  i = nxl, nxr
137        DO  j = nys, nyn
138          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
139             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) &
140                            + dist1(k,j+1,i) + dist1(k+1,j,i) &
141                            ) / 12.0
142          ENDDO
143          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
144              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
145                            + 6.0 * dist1(k,j,i)                               &
146                            ) / 12.0
147          ENDDO
148        ENDDO
149     ENDDO
150
151!
152!-- Exchange of ghost points for the filtered perturbation.
153!-- Afterwards, filter operation and exchange of ghost points are repeated.
154    CALL exchange_horiz( dist2, nbgp )
155
156    DO  i = nxl, nxr
157       DO  j = nys, nyn
158          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
159             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
160                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
161                            + 6.0 * dist2(k,j,i)                               &
162                            ) / 12.0
163          ENDDO
164       ENDDO
165    ENDDO
166
167    CALL exchange_horiz( dist1, nbgp )
168
169!
170!-- Remove perturbations below topography (including one gridpoint above it
171!-- in order to allow for larger timesteps at the beginning of the simulation
172!-- (diffusion criterion))
173    IF ( TRIM( topography ) /= 'flat' )  THEN
174       DO  i = nxlg, nxrg
175          DO  j = nysg, nyng
176             dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0
177          ENDDO
178       ENDDO
179    ENDIF
180
181!
182!-- Random perturbation is added to the array to be disturbed.
183    DO  i = nxlg, nxrg
184       DO  j = nysg, nyng
185          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
186             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
187          ENDDO
188       ENDDO
189    ENDDO
190
191!
192!-- Deallocate the temporary array
193    DEALLOCATE( dist2 )
194
195!
196!-- Set a flag, which indicates that a random perturbation is imposed
197    disturbance_created = .TRUE.
198
199
200    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )
201
202
203 END SUBROUTINE disturb_field
Note: See TracBrowser for help on using the repository browser.