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

Last change on this file since 1279 was 1037, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 6.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!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
[484]20! Current revisions:
[1]21! -----------------
22!
23! Former revisions:
24! -----------------
[3]25! $Id: disturb_field.f90 1037 2012-10-22 14:10:22Z raasch $
[77]26!
[1037]27! 1036 2012-10-22 13:43:42Z raasch
28! code put under GPL (PALM 3.9)
29!
[668]30! 667 2010-12-23 12:06:00Z suehring/gryschka
31! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
32! Calls of exchange_horiz are modified.
33!
[482]34! 420 2010-01-13 15:10:53Z franke
35! Loop was split to make runs reproducible when using ifort compiler
36!
[77]37! 75 2007-03-22 09:54:05Z raasch
38! xrp, ynp eliminated, 2nd+3rd argument removed from exchange horiz
39!
[3]40! RCS Log replace by Id keyword, revision history cleaned up
41!
[1]42! Revision 1.11  2006/08/04 14:31:59  raasch
43! izuf renamed iran
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    USE cpulog
60    USE grid_variables
61    USE indices
62    USE interfaces
63    USE random_function_mod
64
65    IMPLICIT NONE
66
[75]67    INTEGER ::  i, j, k
[667]68    INTEGER ::  nzb_uv_inner(nysg:nyng,nxlg:nxrg)
[1]69
70    REAL    ::  randomnumber,                             &
[667]71                dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
72                field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
[1]73    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  dist2
74
75
76    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
77
78!
79!-- Create an additional temporary array and initialize the arrays needed
80!-- to store the disturbance
[667]81    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]82    dist1 = 0.0
83    dist2 = 0.0
84
85!
86!-- Create the random perturbation and store it on temporary array
87    IF ( random_generator == 'numerical-recipes' )  THEN
88       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
89          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
90             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
91                randomnumber = 3.0 * disturbance_amplitude * &
92                               ( random_function( iran ) - 0.5 )
93                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  &
94                     nyn >= j ) &
95                THEN
96                   dist1(k,j,i) = randomnumber
97                ENDIF
98             ENDDO
99          ENDDO
100       ENDDO
101    ELSEIF ( random_generator == 'system-specific' )  THEN
102       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
103          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
104             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
105#if defined( __nec )
106                randomnumber = 3.0 * disturbance_amplitude * &
107                               ( RANDOM( 0 ) - 0.5 )
108#else
109                CALL RANDOM_NUMBER( randomnumber )
110                randomnumber = 3.0 * disturbance_amplitude * &
111                                ( randomnumber - 0.5 )
112#endif
113                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j ) &
114                THEN
115                   dist1(k,j,i) = randomnumber
116                ENDIF
117             ENDDO
118          ENDDO
119       ENDDO
120
121    ENDIF
122
123!
124!-- Exchange of ghost points for the random perturbation
[420]125
[667]126    CALL exchange_horiz( dist1, nbgp )
[1]127!
128!-- Applying the Shuman filter in order to smooth the perturbations.
129!-- Neighboured grid points in all three directions are used for the
130!-- filter operation.
[420]131!-- Loop has been splitted to make runs reproducible on HLRN systems using
132!-- compiler option -O3
133     DO  i = nxl, nxr
134        DO  j = nys, nyn
[1]135          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
[420]136             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) &
137                            + dist1(k,j+1,i) + dist1(k+1,j,i) &
138                            ) / 12.0
139          ENDDO
140          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
141              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
[1]142                            + 6.0 * dist1(k,j,i)                               &
143                            ) / 12.0
144          ENDDO
[420]145        ENDDO
146     ENDDO
[1]147
148!
149!-- Exchange of ghost points for the filtered perturbation.
150!-- Afterwards, filter operation and exchange of ghost points are repeated.
[667]151    CALL exchange_horiz( dist2, nbgp )
[420]152
[1]153    DO  i = nxl, nxr
154       DO  j = nys, nyn
155          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
156             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
157                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
158                            + 6.0 * dist2(k,j,i)                               &
159                            ) / 12.0
160          ENDDO
161       ENDDO
162    ENDDO
[420]163
[667]164    CALL exchange_horiz( dist1, nbgp )
[1]165
166!
167!-- Remove perturbations below topography (including one gridpoint above it
168!-- in order to allow for larger timesteps at the beginning of the simulation
169!-- (diffusion criterion))
170    IF ( TRIM( topography ) /= 'flat' )  THEN
[667]171       DO  i = nxlg, nxrg
172          DO  j = nysg, nyng
[1]173             dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0
174          ENDDO
175       ENDDO
176    ENDIF
177
178!
179!-- Random perturbation is added to the array to be disturbed.
[667]180    DO  i = nxlg, nxrg
181       DO  j = nysg, nyng
[1]182          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
183             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
184          ENDDO
185       ENDDO
186    ENDDO
187
188!
189!-- Deallocate the temporary array
190    DEALLOCATE( dist2 )
191
192!
193!-- Set a flag, which indicates that a random perturbation is imposed
194    disturbance_created = .TRUE.
195
196
197    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )
198
199
200 END SUBROUTINE disturb_field
Note: See TracBrowser for help on using the repository browser.