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

Last change on this file since 2803 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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