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

Last change on this file since 1886 was 1823, checked in by hoffmann, 9 years ago

last commit documented

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