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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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