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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

  • Property svn:keywords set to Id
File size: 11.7 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! ONLY-attribute added to USE-statements,
23! kind-parameters added to all INTEGER and REAL declaration statements,
24! kinds are defined in new module kinds,
25! revision history before 2012 removed,
26! comment fields (!:) to be used for variable explanations added to
27! all variable declaration statements
28!
29! Former revisions:
30! -----------------
31! $Id: lpm_set_attributes.f90 1320 2014-03-20 08:40:49Z raasch $
32!
33! 1318 2014-03-17 13:35:16Z raasch
34! module interfaces removed
35!
36! 1036 2012-10-22 13:43:42Z raasch
37! code put under GPL (PALM 3.9)
38!
39! 849 2012-03-15 10:35:09Z raasch
40! routine renamed: set_particle_attributes -> lpm_set_attributes
41!
42! 828 2012-02-21 12:00:36Z raasch
43! particle feature color renamed class
44!
45! 271 2009-03-26 00:47:14Z raasch
46! Initial version
47!
48! Description:
49! ------------
50! This routine sets certain particle attributes depending on the values that
51! other PALM variables have at the current particle position.
52!------------------------------------------------------------------------------!
53
54    USE arrays_3d,                                                             &
55        ONLY:  pt, u, v, w, zu, zw
56
57    USE control_parameters,                                                    &
58        ONLY:  atmos_ocean_sign, u_gtrans, v_gtrans, dz
59
60    USE cpulog,                                                                &
61        ONLY:  cpu_log, log_point_s
62
63    USE dvrp_variables,                                                        &
64        ONLY:  color_interval, dvrp_colortable_entries_prt, dvrpsize_interval, &
65               particle_color, particle_dvrpsize
66
67    USE grid_variables,                                                        &
68        ONLY:  ddx, dx, ddy, dy
69
70    USE indices,                                                               &
71        ONLY:  ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
72
73    USE kinds
74
75    USE particle_attributes,                                                   &
76        ONLY:  number_of_particles, offset_ocean_nzt, particles
77
78    USE pegrid
79
80    USE statistics,                                                            &
81        ONLY:  sums, sums_l
82
83    IMPLICIT NONE
84
85    INTEGER(iwp) ::  i        !:
86    INTEGER(iwp) ::  j        !:
87    INTEGER(iwp) ::  k        !:
88    INTEGER(iwp) ::  n        !:
89
90    REAL(wp)    ::  aa        !:
91    REAL(wp)    ::  absuv     !:
92    REAL(wp)    ::  bb        !:
93    REAL(wp)    ::  cc        !:
94    REAL(wp)    ::  dd        !:
95    REAL(wp)    ::  gg        !:
96    REAL(wp)    ::  height    !:
97    REAL(wp)    ::  pt_int    !:
98    REAL(wp)    ::  pt_int_l  !:
99    REAL(wp)    ::  pt_int_u  !:
100    REAL(wp)    ::  u_int     !:
101    REAL(wp)    ::  u_int_l   !:
102    REAL(wp)    ::  u_int_u   !:
103    REAL(wp)    ::  v_int     !:
104    REAL(wp)    ::  v_int_l   !:
105    REAL(wp)    ::  v_int_u   !:
106    REAL(wp)    ::  w_int     !:
107    REAL(wp)    ::  w_int_l   !:
108    REAL(wp)    ::  w_int_u   !:
109    REAL(wp)    ::  x         !:
110    REAL(wp)    ::  y         !:
111
112    CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'start' )
113
114!
115!-- Set particle color
116    IF ( particle_color == 'absuv' )  THEN
117
118!
119!--    Set particle color depending on the absolute value of the horizontal
120!--    velocity
121       DO  n = 1, number_of_particles
122!
123!--       Interpolate u velocity-component, determine left, front, bottom
124!--       index of u-array
125          i = ( particles(n)%x + 0.5 * dx ) * ddx
126          j =   particles(n)%y * ddy
127          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
128              + offset_ocean_nzt                     ! only exact if equidistant
129
130!
131!--       Interpolation of the velocity components in the xy-plane
132          x  = particles(n)%x + ( 0.5 - i ) * dx
133          y  = particles(n)%y - j * dy
134          aa = x**2          + y**2
135          bb = ( dx - x )**2 + y**2
136          cc = x**2          + ( dy - y )**2
137          dd = ( dx - x )**2 + ( dy - y )**2
138          gg = aa + bb + cc + dd
139
140          u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
141                    + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
142                    ) / ( 3.0 * gg ) - u_gtrans
143          IF ( k+1 == nzt+1 )  THEN
144             u_int = u_int_l
145          ELSE
146             u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
147                       + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
148                       ) / ( 3.0 * gg ) - u_gtrans
149             u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
150                               ( u_int_u - u_int_l )
151          ENDIF
152
153!
154!--       Same procedure for interpolation of the v velocity-component (adopt
155!--       index k from u velocity-component)
156          i =   particles(n)%x * ddx
157          j = ( particles(n)%y + 0.5 * dy ) * ddy
158
159          x  = particles(n)%x - i * dx
160          y  = particles(n)%y + ( 0.5 - j ) * dy
161          aa = x**2          + y**2
162          bb = ( dx - x )**2 + y**2
163          cc = x**2          + ( dy - y )**2
164          dd = ( dx - x )**2 + ( dy - y )**2
165          gg = aa + bb + cc + dd
166
167          v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
168                 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
169                 ) / ( 3.0 * gg ) - v_gtrans
170          IF ( k+1 == nzt+1 )  THEN
171             v_int = v_int_l
172          ELSE
173             v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
174                       + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
175                       ) / ( 3.0 * gg ) - v_gtrans
176             v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
177                               ( v_int_u - v_int_l )
178          ENDIF
179
180          absuv = SQRT( u_int**2 + v_int**2 )
181
182!
183!--       Limit values by the given interval and normalize to interval [0,1]
184          absuv = MIN( absuv, color_interval(2) )
185          absuv = MAX( absuv, color_interval(1) )
186
187          absuv = ( absuv - color_interval(1) ) / &
188                  ( color_interval(2) - color_interval(1) )
189
190!
191!--       Number of available colors is defined in init_dvrp
192          particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
193
194       ENDDO
195
196    ELSEIF ( particle_color == 'pt*' )  THEN
197!
198!--    Set particle color depending on the resolved scale temperature
199!--    fluctuation.
200!--    First, calculate the horizontal average of the potential temperature
201!--    (This is also done in flow_statistics, but flow_statistics is called
202!--    after this routine.)
203       sums_l(:,4,0) = 0.0
204       DO  i = nxl, nxr
205          DO  j =  nys, nyn
206             DO  k = nzb, nzt+1
207                sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i)
208             ENDDO
209          ENDDO
210       ENDDO
211
212#if defined( __parallel )
213
214       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
215       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
216                           MPI_REAL, MPI_SUM, comm2d, ierr )
217
218#else
219
220       sums(:,4) = sums_l(:,4,0)
221
222#endif
223
224       sums(:,4) = sums(:,4) / ngp_2dh(0)
225
226       DO  n = 1, number_of_particles
227!
228!--       Interpolate temperature to the current particle position
229          i = particles(n)%x * ddx
230          j = particles(n)%y * ddy
231          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
232              + offset_ocean_nzt                     ! only exact if equidistant
233
234          x  = particles(n)%x - i * dx
235          y  = particles(n)%y - j * dy
236          aa = x**2          + y**2
237          bb = ( dx - x )**2 + y**2
238          cc = x**2          + ( dy - y )**2
239          dd = ( dx - x )**2 + ( dy - y )**2
240          gg = aa + bb + cc + dd
241
242          pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
243                     + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
244                     ) / ( 3.0 * gg ) - sums(k,4)
245
246          pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
247                     + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
248                     ) / ( 3.0 * gg ) - sums(k,4)
249
250          pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
251                              ( pt_int_u - pt_int_l )
252
253!
254!--       Limit values by the given interval and normalize to interval [0,1]
255          pt_int = MIN( pt_int, color_interval(2) )
256          pt_int = MAX( pt_int, color_interval(1) )
257
258          pt_int = ( pt_int - color_interval(1) ) / &
259                   ( color_interval(2) - color_interval(1) )
260
261!
262!--       Number of available colors is defined in init_dvrp
263          particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
264
265       ENDDO
266
267    ELSEIF ( particle_color == 'z' )  THEN
268!
269!--    Set particle color depending on the height above the bottom
270!--    boundary (z=0)
271       DO  n = 1, number_of_particles
272
273          height = particles(n)%z
274!
275!--       Limit values by the given interval and normalize to interval [0,1]
276          height = MIN( height, color_interval(2) )
277          height = MAX( height, color_interval(1) )
278
279          height = ( height - color_interval(1) ) / &
280                   ( color_interval(2) - color_interval(1) )
281
282!
283!--       Number of available colors is defined in init_dvrp
284          particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 )
285
286       ENDDO
287
288    ENDIF
289
290!
291!-- Set particle size for dvrp graphics
292    IF ( particle_dvrpsize == 'absw' )  THEN
293
294       DO  n = 1, number_of_particles
295!
296!--       Interpolate w-component to the current particle position
297          i = particles(n)%x * ddx
298          j = particles(n)%y * ddy
299          k = particles(n)%z / dz
300
301          x  = particles(n)%x - i * dx
302          y  = particles(n)%y - j * dy
303          aa = x**2          + y**2
304          bb = ( dx - x )**2 + y**2
305          cc = x**2          + ( dy - y )**2
306          dd = ( dx - x )**2 + ( dy - y )**2
307          gg = aa + bb + cc + dd
308
309          w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
310                    + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
311                    ) / ( 3.0 * gg )
312
313          IF ( k+1 == nzt+1 )  THEN
314             w_int = w_int_l
315          ELSE
316             w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + ( gg-bb ) * w(k+1,j,i+1)   &
317                       + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) &
318                       ) / ( 3.0 * gg )
319             w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
320                               ( w_int_u - w_int_l )
321          ENDIF
322
323!
324!--       Limit values by the given interval and normalize to interval [0,1]
325          w_int = ABS( w_int )
326          w_int = MIN( w_int, dvrpsize_interval(2) )
327          w_int = MAX( w_int, dvrpsize_interval(1) )
328
329          w_int = ( w_int - dvrpsize_interval(1) ) / &
330                  ( dvrpsize_interval(2) - dvrpsize_interval(1) )
331
332          particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx
333
334       ENDDO
335
336    ENDIF
337
338    CALL cpu_log( log_point_s(49), 'lpm_set_attributes', 'stop' )
339
340
341 END SUBROUTINE lpm_set_attributes
Note: See TracBrowser for help on using the repository browser.