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

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

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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