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

Last change on this file since 1574 was 1360, checked in by hoffmann, 11 years ago

last commit documented

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