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

Last change on this file since 1386 was 1354, checked in by heinze, 11 years ago

last commit documented

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