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

Last change on this file since 1036 was 1036, checked in by raasch, 11 years ago

code has been put under the GNU General Public License (v3)

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