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

Last change on this file since 2675 was 2123, checked in by hoffmann, 8 years ago

last commit documented

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