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

Last change on this file since 264 was 264, checked in by raasch, 14 years ago

new dvrp features added

  • Property svn:keywords set to Id
File size: 4.8 KB
Line 
1 SUBROUTINE set_particle_attributes
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: set_particle_attributes.f90 264 2009-03-21 08:14:44Z raasch $
11!
12! Initial version (2009/03/21)
13!
14! Description:
15! ------------
16! This routine sets certain particle attributes depending on the values that
17! other PALM variables have at the current particle position.
18!------------------------------------------------------------------------------!
19
20    USE arrays_3d
21    USE control_parameters
22    USE cpulog
23    USE dvrp_variables
24    USE grid_variables
25    USE indices
26    USE interfaces
27    USE particle_attributes
28    USE pegrid
29    USE statistics
30
31    IMPLICIT NONE
32
33    INTEGER ::  i, j, k, n
34    REAL    ::  aa, bb, cc, dd, gg, pt_int, pt_int_l, pt_int_u, w_int, &
35                w_int_l, w_int_u, x, y
36
37
38    CALL cpu_log( log_point_s(49), 'set_particle_attributes', 'start' )
39
40!
41!-- Set particle color
42    IF ( particle_color == 'absuv' )  THEN
43
44    ELSEIF ( particle_color == 'pt*' )  THEN
45!
46!--    Set particle color depending on the resolved scale temperature
47!--    fluctuation.
48!--    First, calculate the horizontal average of the potential temperature
49!--    (This is also done in flow_statistics, but flow_statistics is called
50!--    after this routine.)
51       sums_l(:,4,0) = 0.0
52       DO  i = nxl, nxr
53          DO  j =  nys, nyn
54             DO  k = nzb, nzt+1
55                sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i)
56             ENDDO
57          ENDDO
58       ENDDO
59
60#if defined( __parallel )
61
62       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
63                           MPI_REAL, MPI_SUM, comm2d, ierr )
64
65#else
66
67       sums(:,4) = sums_l(:,4,0)
68
69#endif
70
71       sums(:,4) = sums(:,4) / ngp_2dh(0)
72
73       DO  n = 1, number_of_particles
74!
75!--       Interpolate temperature to the current particle position
76          i = particles(n)%x * ddx
77          j = particles(n)%y * ddy
78          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
79              + offset_ocean_nzt                     ! only exact if equidistant
80
81          x  = particles(n)%x - i * dx
82          y  = particles(n)%y - j * dy
83          aa = x**2          + y**2
84          bb = ( dx - x )**2 + y**2
85          cc = x**2          + ( dy - y )**2
86          dd = ( dx - x )**2 + ( dy - y )**2
87          gg = aa + bb + cc + dd
88
89          pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
90                     + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
91                     ) / ( 3.0 * gg ) - sums(k,4)
92
93          pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
94                     + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
95                     ) / ( 3.0 * gg ) - sums(k,4)
96
97          pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
98                              ( pt_int_u - pt_int_l )
99
100!
101!--       Limit values by the given interval and normalize to interval [0,1]
102          pt_int = MIN( pt_int, color_interval(2) )
103          pt_int = MAX( pt_int, color_interval(1) )
104
105          pt_int = ( pt_int - color_interval(1) ) / &
106                   ( color_interval(2) - color_interval(1) )
107
108!
109!--       Number of available colors is defined in init_dvrp
110          particles(n)%color = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
111
112    ENDDO
113
114    ELSEIF ( particle_color == 'z' )  THEN
115
116    ENDIF
117
118!
119!-- Set particle size for dvrp graphics
120    IF ( particle_dvrpsize == 'absw' )  THEN
121
122       DO  n = 1, number_of_particles
123!
124!--       Interpolate w-component to the current particle position
125          i = particles(n)%x * ddx
126          j = particles(n)%y * ddy
127          k = particles(n)%z / dz
128
129          x  = particles(n)%x - i * dx
130          y  = particles(n)%y - j * dy
131          aa = x**2          + y**2
132          bb = ( dx - x )**2 + y**2
133          cc = x**2          + ( dy - y )**2
134          dd = ( dx - x )**2 + ( dy - y )**2
135          gg = aa + bb + cc + dd
136
137          w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
138                    + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
139                    ) / ( 3.0 * gg )
140
141          IF ( k+1 == nzt+1 )  THEN
142             w_int = w_int_l
143          ELSE
144             w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + ( gg-bb ) * w(k+1,j,i+1)   &
145                       + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) &
146                       ) / ( 3.0 * gg )
147             w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
148                               ( w_int_u - w_int_l )
149          ENDIF
150
151          particles(n)%dvrp_psize = ( 0.2 + MIN( 0.6, ABS( w_int ) / 6.0 ) ) * &
152                                    dx * 1.4
153
154       ENDDO
155
156    ENDIF
157
158    CALL cpu_log( log_point_s(49), 'set_particle_attributes', 'stop' )
159
160
161 END SUBROUTINE set_particle_attributes
Note: See TracBrowser for help on using the repository browser.