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

Last change on this file since 1319 was 1319, checked in by raasch, 8 years ago

last commit documented

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