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

Last change on this file since 3709 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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