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

Last change on this file since 600 was 484, checked in by raasch, 15 years ago

typo in file headers removed

  • Property svn:keywords set to Id
File size: 5.7 KB
Line 
1 SUBROUTINE disturb_field( nzb_uv_inner, dist1, field )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: disturb_field.f90 484 2010-02-05 07:36:54Z raasch $
11!
12! 420 2010-01-13 15:10:53Z franke
13! Loop was split to make runs reproducible when using ifort compiler
14!
15! 75 2007-03-22 09:54:05Z raasch
16! xrp, ynp eliminated, 2nd+3rd argument removed from exchange horiz
17!
18! RCS Log replace by Id keyword, revision history cleaned up
19!
20! Revision 1.11  2006/08/04 14:31:59  raasch
21! izuf renamed iran
22!
23! Revision 1.1  1998/02/04 15:40:45  raasch
24! Initial revision
25!
26!
27! Description:
28! ------------
29! Imposing a random perturbation on a 3D-array.
30! On parallel computers, the random number generator is as well called for all
31! gridpoints of the total domain to ensure, regardless of the number of PEs
32! used, that the elements of the array have the same values in the same
33! order in every case. The perturbation range is steered by dist_range.
34!------------------------------------------------------------------------------!
35
36    USE control_parameters
37    USE cpulog
38    USE grid_variables
39    USE indices
40    USE interfaces
41    USE random_function_mod
42
43    IMPLICIT NONE
44
45    INTEGER ::  i, j, k
46    INTEGER ::  nzb_uv_inner(nys-1:nyn+1,nxl-1:nxr+1)
47
48    REAL    ::  randomnumber,                             &
49                dist1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
50                field(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
51    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  dist2
52
53
54    CALL cpu_log( log_point(20), 'disturb_field', 'start' )
55
56!
57!-- Create an additional temporary array and initialize the arrays needed
58!-- to store the disturbance
59    ALLOCATE( dist2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
60    dist1 = 0.0
61    dist2 = 0.0
62
63!
64!-- Create the random perturbation and store it on temporary array
65    IF ( random_generator == 'numerical-recipes' )  THEN
66       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
67          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
68             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
69                randomnumber = 3.0 * disturbance_amplitude * &
70                               ( random_function( iran ) - 0.5 )
71                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.  &
72                     nyn >= j ) &
73                THEN
74                   dist1(k,j,i) = randomnumber
75                ENDIF
76             ENDDO
77          ENDDO
78       ENDDO
79    ELSEIF ( random_generator == 'system-specific' )  THEN
80       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
81          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
82             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
83#if defined( __nec )
84                randomnumber = 3.0 * disturbance_amplitude * &
85                               ( RANDOM( 0 ) - 0.5 )
86#else
87                CALL RANDOM_NUMBER( randomnumber )
88                randomnumber = 3.0 * disturbance_amplitude * &
89                                ( randomnumber - 0.5 )
90#endif
91                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j ) &
92                THEN
93                   dist1(k,j,i) = randomnumber
94                ENDIF
95             ENDDO
96          ENDDO
97       ENDDO
98
99    ENDIF
100
101!
102!-- Exchange of ghost points for the random perturbation
103
104    CALL exchange_horiz( dist1 )
105
106!
107!-- Applying the Shuman filter in order to smooth the perturbations.
108!-- Neighboured grid points in all three directions are used for the
109!-- filter operation.
110!-- Loop has been splitted to make runs reproducible on HLRN systems using
111!-- compiler option -O3
112     DO  i = nxl, nxr
113        DO  j = nys, nyn
114          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
115             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) &
116                            + dist1(k,j+1,i) + dist1(k+1,j,i) &
117                            ) / 12.0
118          ENDDO
119          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
120              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
121                            + 6.0 * dist1(k,j,i)                               &
122                            ) / 12.0
123          ENDDO
124        ENDDO
125     ENDDO
126
127!
128!-- Exchange of ghost points for the filtered perturbation.
129!-- Afterwards, filter operation and exchange of ghost points are repeated.
130    CALL exchange_horiz( dist2 )
131
132    DO  i = nxl, nxr
133       DO  j = nys, nyn
134          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
135             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
136                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
137                            + 6.0 * dist2(k,j,i)                               &
138                            ) / 12.0
139          ENDDO
140       ENDDO
141    ENDDO
142
143    CALL exchange_horiz( dist1 )
144
145!
146!-- Remove perturbations below topography (including one gridpoint above it
147!-- in order to allow for larger timesteps at the beginning of the simulation
148!-- (diffusion criterion))
149    IF ( TRIM( topography ) /= 'flat' )  THEN
150       DO  i = nxl-1, nxr+1
151          DO  j = nys-1, nyn+1
152             dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0
153          ENDDO
154       ENDDO
155    ENDIF
156
157!
158!-- Random perturbation is added to the array to be disturbed.
159    DO  i = nxl-1, nxr+1
160       DO  j = nys-1, nyn+1
161          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
162             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
163          ENDDO
164       ENDDO
165    ENDDO
166
167!
168!-- Deallocate the temporary array
169    DEALLOCATE( dist2 )
170
171!
172!-- Set a flag, which indicates that a random perturbation is imposed
173    disturbance_created = .TRUE.
174
175
176    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )
177
178
179 END SUBROUTINE disturb_field
Note: See TracBrowser for help on using the repository browser.