source: palm/trunk/SOURCE/lpm_set_attributes.f90 @ 850

Last change on this file since 850 was 850, checked in by raasch, 9 years ago

last commit documented

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