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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

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