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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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