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

Last change on this file since 1321 was 1321, checked in by raasch, 10 years ago

last commit documented

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