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

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

Code annotations made doxygen readable

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