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

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