source: palm/trunk/SOURCE/set_particle_attributes.f90 @ 622

Last change on this file since 622 was 622, checked in by raasch, 11 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


  • Property svn:keywords set to Id
File size: 8.9 KB
Line 
1 SUBROUTINE set_particle_attributes
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! optional barriers included in order to speed up collective operations
7!
8! Former revisions:
9! -----------------
10! $Id: set_particle_attributes.f90 622 2010-12-10 08:08:13Z raasch $
11!
12! 271 2009-03-26 00:47:14Z raasch
13! Initial version
14!
15! Description:
16! ------------
17! This routine sets certain particle attributes depending on the values that
18! other PALM variables have at the current particle position.
19!------------------------------------------------------------------------------!
20
21    USE arrays_3d
22    USE control_parameters
23    USE cpulog
24    USE dvrp_variables
25    USE grid_variables
26    USE indices
27    USE interfaces
28    USE particle_attributes
29    USE pegrid
30    USE statistics
31
32    IMPLICIT NONE
33
34    INTEGER ::  i, j, k, n
35    REAL    ::  aa, absuv, bb, cc, dd, gg, height, pt_int, pt_int_l, pt_int_u, &
36                u_int, u_int_l, u_int_u, v_int, v_int_l, v_int_u, w_int,       &
37                w_int_l, w_int_u, x, y
38
39
40    CALL cpu_log( log_point_s(49), 'set_particle_attrib', 'start' )
41
42!
43!-- Set particle color
44    IF ( particle_color == 'absuv' )  THEN
45
46!
47!--    Set particle color depending on the absolute value of the horizontal
48!--    velocity
49       DO  n = 1, number_of_particles
50!
51!--       Interpolate u velocity-component, determine left, front, bottom
52!--       index of u-array
53          i = ( particles(n)%x + 0.5 * dx ) * ddx
54          j =   particles(n)%y * ddy
55          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
56              + offset_ocean_nzt                     ! only exact if equidistant
57
58!
59!--       Interpolation of the velocity components in the xy-plane
60          x  = particles(n)%x + ( 0.5 - i ) * dx
61          y  = particles(n)%y - j * dy
62          aa = x**2          + y**2
63          bb = ( dx - x )**2 + y**2
64          cc = x**2          + ( dy - y )**2
65          dd = ( dx - x )**2 + ( dy - y )**2
66          gg = aa + bb + cc + dd
67
68          u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
69                    + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
70                    ) / ( 3.0 * gg ) - u_gtrans
71          IF ( k+1 == nzt+1 )  THEN
72             u_int = u_int_l
73          ELSE
74             u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
75                       + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
76                       ) / ( 3.0 * gg ) - u_gtrans
77             u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
78                               ( u_int_u - u_int_l )
79          ENDIF
80
81!
82!--       Same procedure for interpolation of the v velocity-component (adopt
83!--       index k from u velocity-component)
84          i =   particles(n)%x * ddx
85          j = ( particles(n)%y + 0.5 * dy ) * ddy
86
87          x  = particles(n)%x - i * dx
88          y  = particles(n)%y + ( 0.5 - j ) * dy
89          aa = x**2          + y**2
90          bb = ( dx - x )**2 + y**2
91          cc = x**2          + ( dy - y )**2
92          dd = ( dx - x )**2 + ( dy - y )**2
93          gg = aa + bb + cc + dd
94
95          v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
96                 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
97                 ) / ( 3.0 * gg ) - v_gtrans
98          IF ( k+1 == nzt+1 )  THEN
99             v_int = v_int_l
100          ELSE
101             v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
102                       + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
103                       ) / ( 3.0 * gg ) - v_gtrans
104             v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
105                               ( v_int_u - v_int_l )
106          ENDIF
107
108          absuv = SQRT( u_int**2 + v_int**2 )
109
110!
111!--       Limit values by the given interval and normalize to interval [0,1]
112          absuv = MIN( absuv, color_interval(2) )
113          absuv = MAX( absuv, color_interval(1) )
114
115          absuv = ( absuv - color_interval(1) ) / &
116                  ( color_interval(2) - color_interval(1) )
117
118!
119!--       Number of available colors is defined in init_dvrp
120          particles(n)%color = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
121
122       ENDDO
123
124    ELSEIF ( particle_color == 'pt*' )  THEN
125!
126!--    Set particle color depending on the resolved scale temperature
127!--    fluctuation.
128!--    First, calculate the horizontal average of the potential temperature
129!--    (This is also done in flow_statistics, but flow_statistics is called
130!--    after this routine.)
131       sums_l(:,4,0) = 0.0
132       DO  i = nxl, nxr
133          DO  j =  nys, nyn
134             DO  k = nzb, nzt+1
135                sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i)
136             ENDDO
137          ENDDO
138       ENDDO
139
140#if defined( __parallel )
141
142       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
143       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
144                           MPI_REAL, MPI_SUM, comm2d, ierr )
145
146#else
147
148       sums(:,4) = sums_l(:,4,0)
149
150#endif
151
152       sums(:,4) = sums(:,4) / ngp_2dh(0)
153
154       DO  n = 1, number_of_particles
155!
156!--       Interpolate temperature to the current particle position
157          i = particles(n)%x * ddx
158          j = particles(n)%y * ddy
159          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
160              + offset_ocean_nzt                     ! only exact if equidistant
161
162          x  = particles(n)%x - i * dx
163          y  = particles(n)%y - j * dy
164          aa = x**2          + y**2
165          bb = ( dx - x )**2 + y**2
166          cc = x**2          + ( dy - y )**2
167          dd = ( dx - x )**2 + ( dy - y )**2
168          gg = aa + bb + cc + dd
169
170          pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
171                     + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
172                     ) / ( 3.0 * gg ) - sums(k,4)
173
174          pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
175                     + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
176                     ) / ( 3.0 * gg ) - sums(k,4)
177
178          pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
179                              ( pt_int_u - pt_int_l )
180
181!
182!--       Limit values by the given interval and normalize to interval [0,1]
183          pt_int = MIN( pt_int, color_interval(2) )
184          pt_int = MAX( pt_int, color_interval(1) )
185
186          pt_int = ( pt_int - color_interval(1) ) / &
187                   ( color_interval(2) - color_interval(1) )
188
189!
190!--       Number of available colors is defined in init_dvrp
191          particles(n)%color = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
192
193       ENDDO
194
195    ELSEIF ( particle_color == 'z' )  THEN
196!
197!--    Set particle color depending on the height above the bottom
198!--    boundary (z=0)
199       DO  n = 1, number_of_particles
200
201          height = particles(n)%z
202!
203!--       Limit values by the given interval and normalize to interval [0,1]
204          height = MIN( height, color_interval(2) )
205          height = MAX( height, color_interval(1) )
206
207          height = ( height - color_interval(1) ) / &
208                   ( color_interval(2) - color_interval(1) )
209
210!
211!--       Number of available colors is defined in init_dvrp
212          particles(n)%color = 1 + height * ( dvrp_colortable_entries_prt - 1 )
213
214       ENDDO
215
216    ENDIF
217
218!
219!-- Set particle size for dvrp graphics
220    IF ( particle_dvrpsize == 'absw' )  THEN
221
222       DO  n = 1, number_of_particles
223!
224!--       Interpolate w-component to the current particle position
225          i = particles(n)%x * ddx
226          j = particles(n)%y * ddy
227          k = particles(n)%z / dz
228
229          x  = particles(n)%x - i * dx
230          y  = particles(n)%y - j * dy
231          aa = x**2          + y**2
232          bb = ( dx - x )**2 + y**2
233          cc = x**2          + ( dy - y )**2
234          dd = ( dx - x )**2 + ( dy - y )**2
235          gg = aa + bb + cc + dd
236
237          w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
238                    + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
239                    ) / ( 3.0 * gg )
240
241          IF ( k+1 == nzt+1 )  THEN
242             w_int = w_int_l
243          ELSE
244             w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + ( gg-bb ) * w(k+1,j,i+1)   &
245                       + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) &
246                       ) / ( 3.0 * gg )
247             w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
248                               ( w_int_u - w_int_l )
249          ENDIF
250
251!
252!--       Limit values by the given interval and normalize to interval [0,1]
253          w_int = ABS( w_int )
254          w_int = MIN( w_int, dvrpsize_interval(2) )
255          w_int = MAX( w_int, dvrpsize_interval(1) )
256
257          w_int = ( w_int - dvrpsize_interval(1) ) / &
258                  ( dvrpsize_interval(2) - dvrpsize_interval(1) )
259
260          particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx
261
262       ENDDO
263
264    ENDIF
265
266    CALL cpu_log( log_point_s(49), 'set_particle_attrib', 'stop' )
267
268
269 END SUBROUTINE set_particle_attributes
Note: See TracBrowser for help on using the repository browser.