source: palm/trunk/SOURCE/set_particle_attributes.f90 @ 754

Last change on this file since 754 was 623, checked in by raasch, 14 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 8.9 KB
Line 
1 SUBROUTINE set_particle_attributes
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: set_particle_attributes.f90 623 2010-12-10 08:52:17Z raasch $
11!
12! 622 2010-12-10 08:08:13Z raasch
13! optional barriers included in order to speed up collective operations
14!
15! 271 2009-03-26 00:47:14Z raasch
16! Initial version
17!
18! Description:
19! ------------
20! This routine sets certain particle attributes depending on the values that
21! other PALM variables have at the current particle position.
22!------------------------------------------------------------------------------!
23
24    USE arrays_3d
25    USE control_parameters
26    USE cpulog
27    USE dvrp_variables
28    USE grid_variables
29    USE indices
30    USE interfaces
31    USE particle_attributes
32    USE pegrid
33    USE statistics
34
35    IMPLICIT NONE
36
37    INTEGER ::  i, j, k, n
38    REAL    ::  aa, absuv, bb, cc, dd, gg, height, pt_int, pt_int_l, pt_int_u, &
39                u_int, u_int_l, u_int_u, v_int, v_int_l, v_int_u, w_int,       &
40                w_int_l, w_int_u, x, y
41
42
43    CALL cpu_log( log_point_s(49), 'set_particle_attrib', 'start' )
44
45!
46!-- Set particle color
47    IF ( particle_color == 'absuv' )  THEN
48
49!
50!--    Set particle color depending on the absolute value of the horizontal
51!--    velocity
52       DO  n = 1, number_of_particles
53!
54!--       Interpolate u velocity-component, determine left, front, bottom
55!--       index of u-array
56          i = ( particles(n)%x + 0.5 * dx ) * ddx
57          j =   particles(n)%y * ddy
58          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
59              + offset_ocean_nzt                     ! only exact if equidistant
60
61!
62!--       Interpolation of the velocity components in the xy-plane
63          x  = particles(n)%x + ( 0.5 - i ) * dx
64          y  = particles(n)%y - j * dy
65          aa = x**2          + y**2
66          bb = ( dx - x )**2 + y**2
67          cc = x**2          + ( dy - y )**2
68          dd = ( dx - x )**2 + ( dy - y )**2
69          gg = aa + bb + cc + dd
70
71          u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
72                    + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
73                    ) / ( 3.0 * gg ) - u_gtrans
74          IF ( k+1 == nzt+1 )  THEN
75             u_int = u_int_l
76          ELSE
77             u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
78                       + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
79                       ) / ( 3.0 * gg ) - u_gtrans
80             u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
81                               ( u_int_u - u_int_l )
82          ENDIF
83
84!
85!--       Same procedure for interpolation of the v velocity-component (adopt
86!--       index k from u velocity-component)
87          i =   particles(n)%x * ddx
88          j = ( particles(n)%y + 0.5 * dy ) * ddy
89
90          x  = particles(n)%x - i * dx
91          y  = particles(n)%y + ( 0.5 - j ) * dy
92          aa = x**2          + y**2
93          bb = ( dx - x )**2 + y**2
94          cc = x**2          + ( dy - y )**2
95          dd = ( dx - x )**2 + ( dy - y )**2
96          gg = aa + bb + cc + dd
97
98          v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
99                 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
100                 ) / ( 3.0 * gg ) - v_gtrans
101          IF ( k+1 == nzt+1 )  THEN
102             v_int = v_int_l
103          ELSE
104             v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
105                       + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
106                       ) / ( 3.0 * gg ) - v_gtrans
107             v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
108                               ( v_int_u - v_int_l )
109          ENDIF
110
111          absuv = SQRT( u_int**2 + v_int**2 )
112
113!
114!--       Limit values by the given interval and normalize to interval [0,1]
115          absuv = MIN( absuv, color_interval(2) )
116          absuv = MAX( absuv, color_interval(1) )
117
118          absuv = ( absuv - color_interval(1) ) / &
119                  ( color_interval(2) - color_interval(1) )
120
121!
122!--       Number of available colors is defined in init_dvrp
123          particles(n)%color = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
124
125       ENDDO
126
127    ELSEIF ( particle_color == 'pt*' )  THEN
128!
129!--    Set particle color depending on the resolved scale temperature
130!--    fluctuation.
131!--    First, calculate the horizontal average of the potential temperature
132!--    (This is also done in flow_statistics, but flow_statistics is called
133!--    after this routine.)
134       sums_l(:,4,0) = 0.0
135       DO  i = nxl, nxr
136          DO  j =  nys, nyn
137             DO  k = nzb, nzt+1
138                sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i)
139             ENDDO
140          ENDDO
141       ENDDO
142
143#if defined( __parallel )
144
145       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
146       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
147                           MPI_REAL, MPI_SUM, comm2d, ierr )
148
149#else
150
151       sums(:,4) = sums_l(:,4,0)
152
153#endif
154
155       sums(:,4) = sums(:,4) / ngp_2dh(0)
156
157       DO  n = 1, number_of_particles
158!
159!--       Interpolate temperature to the current particle position
160          i = particles(n)%x * ddx
161          j = particles(n)%y * ddy
162          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
163              + offset_ocean_nzt                     ! only exact if equidistant
164
165          x  = particles(n)%x - i * dx
166          y  = particles(n)%y - j * dy
167          aa = x**2          + y**2
168          bb = ( dx - x )**2 + y**2
169          cc = x**2          + ( dy - y )**2
170          dd = ( dx - x )**2 + ( dy - y )**2
171          gg = aa + bb + cc + dd
172
173          pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
174                     + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
175                     ) / ( 3.0 * gg ) - sums(k,4)
176
177          pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
178                     + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
179                     ) / ( 3.0 * gg ) - sums(k,4)
180
181          pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
182                              ( pt_int_u - pt_int_l )
183
184!
185!--       Limit values by the given interval and normalize to interval [0,1]
186          pt_int = MIN( pt_int, color_interval(2) )
187          pt_int = MAX( pt_int, color_interval(1) )
188
189          pt_int = ( pt_int - color_interval(1) ) / &
190                   ( color_interval(2) - color_interval(1) )
191
192!
193!--       Number of available colors is defined in init_dvrp
194          particles(n)%color = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
195
196       ENDDO
197
198    ELSEIF ( particle_color == 'z' )  THEN
199!
200!--    Set particle color depending on the height above the bottom
201!--    boundary (z=0)
202       DO  n = 1, number_of_particles
203
204          height = particles(n)%z
205!
206!--       Limit values by the given interval and normalize to interval [0,1]
207          height = MIN( height, color_interval(2) )
208          height = MAX( height, color_interval(1) )
209
210          height = ( height - color_interval(1) ) / &
211                   ( color_interval(2) - color_interval(1) )
212
213!
214!--       Number of available colors is defined in init_dvrp
215          particles(n)%color = 1 + height * ( dvrp_colortable_entries_prt - 1 )
216
217       ENDDO
218
219    ENDIF
220
221!
222!-- Set particle size for dvrp graphics
223    IF ( particle_dvrpsize == 'absw' )  THEN
224
225       DO  n = 1, number_of_particles
226!
227!--       Interpolate w-component to the current particle position
228          i = particles(n)%x * ddx
229          j = particles(n)%y * ddy
230          k = particles(n)%z / dz
231
232          x  = particles(n)%x - i * dx
233          y  = particles(n)%y - j * dy
234          aa = x**2          + y**2
235          bb = ( dx - x )**2 + y**2
236          cc = x**2          + ( dy - y )**2
237          dd = ( dx - x )**2 + ( dy - y )**2
238          gg = aa + bb + cc + dd
239
240          w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
241                    + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
242                    ) / ( 3.0 * gg )
243
244          IF ( k+1 == nzt+1 )  THEN
245             w_int = w_int_l
246          ELSE
247             w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + ( gg-bb ) * w(k+1,j,i+1)   &
248                       + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) &
249                       ) / ( 3.0 * gg )
250             w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
251                               ( w_int_u - w_int_l )
252          ENDIF
253
254!
255!--       Limit values by the given interval and normalize to interval [0,1]
256          w_int = ABS( w_int )
257          w_int = MIN( w_int, dvrpsize_interval(2) )
258          w_int = MAX( w_int, dvrpsize_interval(1) )
259
260          w_int = ( w_int - dvrpsize_interval(1) ) / &
261                  ( dvrpsize_interval(2) - dvrpsize_interval(1) )
262
263          particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx
264
265       ENDDO
266
267    ENDIF
268
269    CALL cpu_log( log_point_s(49), 'set_particle_attrib', 'stop' )
270
271
272 END SUBROUTINE set_particle_attributes
Note: See TracBrowser for help on using the repository browser.