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

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

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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