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

Last change on this file since 845 was 829, checked in by raasch, 12 years ago

last commit documented

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