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

Last change on this file since 3065 was 3065, checked in by Giersch, 3 years ago

New vertical stretching procedure has been introduced

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