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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

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