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

Last change on this file since 2089 was 2001, checked in by knoop, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 17.7 KB
Line 
1!> @file lpm_set_attributes.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2016 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_set_attributes.f90 2001 2016-08-20 18:41:22Z suehring $
27!
28! 2000 2016-08-20 18:09:15Z knoop
29! Forced header and separation lines into 80 columns
30!
31! 1822 2016-04-07 07:49:42Z hoffmann
32! Unused variables removed.
33!
34! 1682 2015-10-07 23:56:08Z knoop
35! Code annotations made doxygen readable
36!
37! 1359 2014-04-11 17:15:14Z hoffmann
38! New particle structure integrated.
39! Kind definition added to all floating point numbers.
40!
41! 1320 2014-03-20 08:40:49Z raasch
42! ONLY-attribute added to USE-statements,
43! kind-parameters added to all INTEGER and REAL declaration statements,
44! kinds are defined in new module kinds,
45! revision history before 2012 removed,
46! comment fields (!:) to be used for variable explanations added to
47! all variable declaration statements
48!
49! 1318 2014-03-17 13:35:16Z raasch
50! module interfaces removed
51!
52! 1036 2012-10-22 13:43:42Z raasch
53! code put under GPL (PALM 3.9)
54!
55! 849 2012-03-15 10:35:09Z raasch
56! routine renamed: set_particle_attributes -> lpm_set_attributes
57!
58! 828 2012-02-21 12:00:36Z raasch
59! particle feature color renamed class
60!
61! 271 2009-03-26 00:47:14Z raasch
62! Initial version
63!
64! Description:
65! ------------
66!> This routine sets certain particle attributes depending on the values that
67!> other PALM variables have at the current particle position.
68!------------------------------------------------------------------------------!
69 SUBROUTINE lpm_set_attributes
70 
71
72    USE arrays_3d,                                                             &
73        ONLY:  pt, u, v, w, zu, zw
74
75    USE control_parameters,                                                    &
76        ONLY:  u_gtrans, v_gtrans, dz
77
78    USE cpulog,                                                                &
79        ONLY:  cpu_log, log_point_s
80
81    USE dvrp_variables,                                                        &
82        ONLY:  color_interval, dvrp_colortable_entries_prt, dvrpsize_interval, &
83               particle_color, particle_dvrpsize
84
85    USE grid_variables,                                                        &
86        ONLY:  dx, dy
87
88    USE indices,                                                               &
89        ONLY:  ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
90
91    USE kinds
92
93    USE particle_attributes,                                                   &
94        ONLY:  block_offset, grid_particles, number_of_particles, particles,   &
95               prt_count
96
97    USE pegrid
98
99    USE statistics,                                                            &
100        ONLY:  sums, sums_l
101
102    IMPLICIT NONE
103
104    INTEGER(iwp) ::  i        !<
105    INTEGER(iwp) ::  ip       !<
106    INTEGER(iwp) ::  j        !<
107    INTEGER(iwp) ::  jp       !<
108    INTEGER(iwp) ::  k        !<
109    INTEGER(iwp) ::  kp       !<
110    INTEGER(iwp) ::  n        !<
111    INTEGER(iwp) ::  nb       !<
112
113    INTEGER(iwp), DIMENSION(0:7) ::  start_index !<
114    INTEGER(iwp), DIMENSION(0:7) ::  end_index   !<
115
116    REAL(wp)    ::  aa        !<
117    REAL(wp)    ::  absuv     !<
118    REAL(wp)    ::  bb        !<
119    REAL(wp)    ::  cc        !<
120    REAL(wp)    ::  dd        !<
121    REAL(wp)    ::  gg        !<
122    REAL(wp)    ::  height    !<
123    REAL(wp)    ::  pt_int    !<
124    REAL(wp)    ::  pt_int_l  !<
125    REAL(wp)    ::  pt_int_u  !<
126    REAL(wp)    ::  u_int_l   !<
127    REAL(wp)    ::  u_int_u   !<
128    REAL(wp)    ::  v_int_l   !<
129    REAL(wp)    ::  v_int_u   !<
130    REAL(wp)    ::  w_int     !<
131    REAL(wp)    ::  w_int_l   !<
132    REAL(wp)    ::  w_int_u   !<
133    REAL(wp)    ::  x         !<
134    REAL(wp)    ::  y         !<
135
136    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_int                  !<
137    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_int                  !<
138    REAL(wp), DIMENSION(:), ALLOCATABLE ::  xv                     !<
139    REAL(wp), DIMENSION(:), ALLOCATABLE ::  yv                     !<
140    REAL(wp), DIMENSION(:), ALLOCATABLE ::  zv                     !<
141
142    CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'start' )
143
144!
145!-- Set particle color
146    IF ( particle_color == 'absuv' )  THEN
147
148!
149!--    Set particle color depending on the absolute value of the horizontal
150!--    velocity
151       DO  ip = nxl, nxr
152          DO  jp = nys, nyn
153             DO  kp = nzb+1, nzt
154
155                number_of_particles = prt_count(kp,jp,ip)
156                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
157                IF ( number_of_particles <= 0 )  CYCLE
158                start_index = grid_particles(kp,jp,ip)%start_index
159                end_index   = grid_particles(kp,jp,ip)%end_index
160
161                ALLOCATE( u_int(1:number_of_particles), &
162                          v_int(1:number_of_particles), &
163                          xv(1:number_of_particles),    &
164                          yv(1:number_of_particles),    &
165                          zv(1:number_of_particles) )
166
167                xv = particles(1:number_of_particles)%x
168                yv = particles(1:number_of_particles)%y
169                zv = particles(1:number_of_particles)%z
170
171                DO  nb = 0,7
172
173                   i = ip
174                   j = jp + block_offset(nb)%j_off
175                   k = kp + block_offset(nb)%k_off
176
177                   DO  n = start_index(nb), end_index(nb)
178!
179!--                   Interpolation of the velocity components in the xy-plane
180                      x  = xv(n) + ( 0.5_wp - i ) * dx
181                      y  = yv(n) - j * dy
182                      aa = x**2          + y**2
183                      bb = ( dx - x )**2 + y**2
184                      cc = x**2          + ( dy - y )**2
185                      dd = ( dx - x )**2 + ( dy - y )**2
186                      gg = aa + bb + cc + dd
187
188                      u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) *     &
189                                  u(k,j,i+1) + ( gg - cc ) * u(k,j+1,i) +      &
190                                  ( gg - dd ) * u(k,j+1,i+1)                   &
191                                ) / ( 3.0_wp * gg ) - u_gtrans
192
193                      IF ( k+1 == nzt+1 )  THEN
194                         u_int(n) = u_int_l
195                      ELSE
196                         u_int_u = ( ( gg - aa ) * u(k+1,j,i)   + ( gg - bb ) *  &
197                                     u(k+1,j,i+1) + ( gg - cc ) * u(k+1,j+1,i) + &
198                                     ( gg - dd ) * u(k+1,j+1,i+1)                &
199                                   ) / ( 3.0_wp * gg ) - u_gtrans
200                         u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *         &
201                                           ( u_int_u - u_int_l )
202                      ENDIF
203
204                   ENDDO
205
206                   i = ip + block_offset(nb)%i_off
207                   j = jp
208                   k = kp + block_offset(nb)%k_off
209
210                   DO  n = start_index(nb), end_index(nb)
211!
212!--                   Same procedure for interpolation of the v velocity-component
213                      x  = xv(n) - i * dx
214                      y  = yv(n) + ( 0.5_wp - j ) * dy
215                      aa = x**2          + y**2
216                      bb = ( dx - x )**2 + y**2
217                      cc = x**2          + ( dy - y )**2
218                      dd = ( dx - x )**2 + ( dy - y )**2
219                      gg = aa + bb + cc + dd
220
221                      v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) *     &
222                                  v(k,j,i+1) + ( gg - cc ) * v(k,j+1,i) +      &
223                                  ( gg - dd ) * v(k,j+1,i+1)                   &
224                                ) / ( 3.0_wp * gg ) - v_gtrans
225
226                      IF ( k+1 == nzt+1 )  THEN
227                         v_int(n) = v_int_l
228                      ELSE
229                         v_int_u  = ( ( gg - aa ) * v(k+1,j,i) + ( gg - bb ) *    &
230                                      v(k+1,j,i+1) + ( gg - cc ) * v(k+1,j+1,i) + &
231                                      ( gg - dd ) * v(k+1,j+1,i+1)                &
232                                    ) / ( 3.0_wp * gg ) - v_gtrans
233                         v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz *         &
234                                           ( v_int_u - v_int_l )
235                      ENDIF
236
237                   ENDDO
238
239                ENDDO
240
241                DO  n = 1, number_of_particles
242
243                   absuv = SQRT( u_int(n)**2 + v_int(n)**2 )
244
245!
246!--                Limit values by the given interval and normalize to
247!--                interval [0,1]
248                   absuv = MIN( absuv, color_interval(2) )
249                   absuv = MAX( absuv, color_interval(1) )
250
251                   absuv = ( absuv - color_interval(1) ) / &
252                           ( color_interval(2) - color_interval(1) )
253
254!
255!--                Number of available colors is defined in init_dvrp
256                   particles(n)%class = 1 + absuv *                            &
257                                            ( dvrp_colortable_entries_prt - 1 )
258
259                ENDDO
260
261                DEALLOCATE( u_int, v_int, xv, yv, zv )
262
263             ENDDO
264          ENDDO
265       ENDDO
266
267    ELSEIF ( particle_color == 'pt*' )  THEN
268!
269!--    Set particle color depending on the resolved scale temperature
270!--    fluctuation.
271!--    First, calculate the horizontal average of the potential temperature
272!--    (This is also done in flow_statistics, but flow_statistics is called
273!--    after this routine.)
274       sums_l(:,4,0) = 0.0_wp
275       DO  i = nxl, nxr
276          DO  j =  nys, nyn
277             DO  k = nzb, nzt+1
278                sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i)
279             ENDDO
280          ENDDO
281       ENDDO
282
283#if defined( __parallel )
284       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
285       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
286                           MPI_REAL, MPI_SUM, comm2d, ierr )
287#else
288       sums(:,4) = sums_l(:,4,0)
289#endif
290       sums(:,4) = sums(:,4) / ngp_2dh(0)
291
292       DO  ip = nxl, nxr
293          DO  jp = nys, nyn
294             DO  kp = nzb+1, nzt
295
296                number_of_particles = prt_count(kp,jp,ip)
297                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
298                IF ( number_of_particles <= 0 )  CYCLE
299                start_index = grid_particles(kp,jp,ip)%start_index
300                end_index   = grid_particles(kp,jp,ip)%end_index
301
302                ALLOCATE( xv(1:number_of_particles), &
303                          yv(1:number_of_particles), &
304                          zv(1:number_of_particles) )
305
306                xv = particles(1:number_of_particles)%x
307                yv = particles(1:number_of_particles)%y
308                zv = particles(1:number_of_particles)%z
309
310                DO  nb = 0,7
311
312                   i = ip + block_offset(nb)%i_off
313                   j = jp + block_offset(nb)%j_off
314                   k = kp + block_offset(nb)%k_off
315
316                   DO  n = start_index(nb), end_index(nb)
317!
318!--                   Interpolate temperature to the current particle position
319                      x  = xv(n) - i * dx
320                      y  = yv(n) - j * dy
321                      aa = x**2          + y**2
322                      bb = ( dx - x )**2 + y**2
323                      cc = x**2          + ( dy - y )**2
324                      dd = ( dx - x )**2 + ( dy - y )**2
325                      gg = aa + bb + cc + dd
326
327                      pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) *   &
328                                   pt(k,j,i+1) + ( gg - cc ) * pt(k,j+1,i) +   &
329                                   ( gg - dd ) * pt(k,j+1,i+1)                 &
330                                 ) / ( 3.0_wp * gg ) - sums(k,4)
331
332                      pt_int_u = ( ( gg - aa ) * pt(k+1,j,i) + ( gg - bb ) *     &
333                                   pt(k+1,j,i+1) + ( gg - cc ) * pt(k+1,j+1,i) + &
334                                   ( gg - dd ) * pt(k+1,j+1,i+1)                 &
335                                 ) / ( 3.0_wp * gg ) - sums(k,4)
336
337                      pt_int = pt_int_l + ( zv(n) - zu(k) ) / dz *    &
338                                          ( pt_int_u - pt_int_l )
339
340!
341!--                   Limit values by the given interval and normalize to
342!--                   interval [0,1]
343                      pt_int = MIN( pt_int, color_interval(2) )
344                      pt_int = MAX( pt_int, color_interval(1) )
345
346                      pt_int = ( pt_int - color_interval(1) ) /                &
347                               ( color_interval(2) - color_interval(1) )
348
349!
350!--                   Number of available colors is defined in init_dvrp
351                      particles(n)%class = 1 + pt_int *                        &
352                                           ( dvrp_colortable_entries_prt - 1 )
353
354                   ENDDO
355                ENDDO
356
357                DEALLOCATE( xv, yv, zv )
358
359             ENDDO
360          ENDDO
361       ENDDO
362
363    ELSEIF ( particle_color == 'z' )  THEN
364!
365!--    Set particle color depending on the height above the bottom
366!--    boundary (z=0)
367       DO  ip = nxl, nxr
368          DO  jp = nys, nyn
369             DO  kp = nzb+1, nzt
370
371                number_of_particles = prt_count(kp,jp,ip)
372                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
373                IF ( number_of_particles <= 0 )  CYCLE
374                DO  n = 1, number_of_particles
375
376                   height = particles(n)%z
377!
378!--                Limit values by the given interval and normalize to
379!--                interval [0,1]
380                   height = MIN( height, color_interval(2) )
381                   height = MAX( height, color_interval(1) )
382
383                   height = ( height - color_interval(1) ) / &
384                            ( color_interval(2) - color_interval(1) )
385
386!
387!--                Number of available colors is defined in init_dvrp
388                   particles(n)%class = 1 + height *                           &
389                                            ( dvrp_colortable_entries_prt - 1 )
390
391                ENDDO
392
393             ENDDO
394          ENDDO
395       ENDDO
396
397    ENDIF
398
399!
400!-- Set particle size for dvrp graphics
401    IF ( particle_dvrpsize == 'absw' )  THEN
402
403       DO  ip = nxl, nxr
404          DO  jp = nys, nyn
405             DO  kp = nzb+1, nzt
406
407                number_of_particles = prt_count(kp,jp,ip)
408                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
409                IF ( number_of_particles <= 0 )  CYCLE
410                start_index = grid_particles(kp,jp,ip)%start_index
411                end_index   = grid_particles(kp,jp,ip)%end_index
412
413                ALLOCATE( xv(1:number_of_particles), &
414                          yv(1:number_of_particles) )
415
416                xv = particles(1:number_of_particles)%x
417                yv = particles(1:number_of_particles)%y
418                zv = particles(1:number_of_particles)%z
419
420                DO  nb = 0,7
421
422                   i = ip + block_offset(nb)%i_off
423                   j = jp + block_offset(nb)%j_off
424                   k = kp-1
425
426                   DO  n = start_index(nb), end_index(nb)
427!
428!--                   Interpolate w-component to the current particle position
429                      x  = xv(n) - i * dx
430                      y  = yv(n) - j * dy
431                      aa = x**2          + y**2
432                      bb = ( dx - x )**2 + y**2
433                      cc = x**2          + ( dy - y )**2
434                      dd = ( dx - x )**2 + ( dy - y )**2
435                      gg = aa + bb + cc + dd
436
437                      w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) *     &
438                                  w(k,j,i+1) + ( gg - cc ) * w(k,j+1,i) +      &
439                                  ( gg - dd ) * w(k,j+1,i+1)                   &
440                                ) / ( 3.0_wp * gg )
441
442                      IF ( k+1 == nzt+1 )  THEN
443                         w_int = w_int_l
444                      ELSE
445                         w_int_u = ( ( gg - aa ) * w(k+1,j,i)   + ( gg - bb ) *  &
446                                     w(k+1,j,i+1) + ( gg - cc ) * w(k+1,j+1,i) + &
447                                     ( gg - dd ) * w(k+1,j+1,i+1)                &
448                                   ) / ( 3.0_wp * gg )
449                         w_int = w_int_l + ( zv(n) - zw(k) ) / dz *     &
450                                           ( w_int_u - w_int_l )
451                      ENDIF
452
453!
454!--                   Limit values by the given interval and normalize to
455!--                   interval [0,1]
456                      w_int = ABS( w_int )
457                      w_int = MIN( w_int, dvrpsize_interval(2) )
458                      w_int = MAX( w_int, dvrpsize_interval(1) )
459
460                      w_int = ( w_int - dvrpsize_interval(1) ) / &
461                              ( dvrpsize_interval(2) - dvrpsize_interval(1) )
462
463                      particles(n)%dvrp_psize = ( 0.25_wp + w_int * 0.6_wp ) * &
464                                                dx
465
466                   ENDDO
467                ENDDO
468
469                DEALLOCATE( xv, yv, zv )
470
471             ENDDO
472          ENDDO
473       ENDDO
474
475    ENDIF
476
477    CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'stop' )
478
479
480 END SUBROUTINE lpm_set_attributes
Note: See TracBrowser for help on using the repository browser.