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

Last change on this file since 1401 was 1401, checked in by knoop, 7 years ago

last commit documented

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