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

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

code has been put under the GNU General Public License (v3)

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