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

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

last commit documented

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