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

Last change on this file since 1359 was 1359, checked in by hoffmann, 7 years ago

new Lagrangian particle structure integrated

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