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

Last change on this file since 1321 was 1321, checked in by raasch, 7 years ago

last commit documented

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