source: palm/trunk/SOURCE/lpm.f90 @ 1319

Last change on this file since 1319 was 1319, checked in by raasch, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 9.8 KB
Line 
1 SUBROUTINE lpm
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!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm.f90 1319 2014-03-17 15:08:44Z raasch $
27!
28! 1318 2014-03-17 13:35:16Z raasch
29! module interfaces removed
30!
31! 1036 2012-10-22 13:43:42Z raasch
32! code put under GPL (PALM 3.9)
33!
34! 851 2012-03-15 14:32:58Z raasch
35! Bugfix: resetting of particle_mask and tail mask moved from routine
36! lpm_exchange_horiz to here (end of sub-timestep loop)
37!
38! 849 2012-03-15 10:35:09Z raasch
39! original routine advec_particles split into several subroutines and renamed
40! lpm
41!
42! 831 2012-02-22 00:29:39Z raasch
43! thermal_conductivity_l and diff_coeff_l now depend on temperature and
44! pressure
45!
46! 828 2012-02-21 12:00:36Z raasch
47! fast hall/wang kernels with fixed radius/dissipation classes added,
48! particle feature color renamed class, routine colker renamed
49! recalculate_kernel,
50! lower limit for droplet radius changed from 1E-7 to 1E-8
51!
52! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
53!
54! 825 2012-02-19 03:03:44Z raasch
55! droplet growth by condensation may include curvature and solution effects,
56! initialisation of temporary particle array for resorting removed,
57! particle attributes speed_x|y|z_sgs renamed rvar1|2|3,
58! module wang_kernel_mod renamed lpm_collision_kernels_mod,
59! wang_collision_kernel renamed wang_kernel
60!
61! 799 2011-12-21 17:48:03Z franke
62! Implementation of Wang collision kernel and corresponding new parameter
63! wang_collision_kernel
64!
65! 792 2011-12-01 raasch
66! particle arrays (particles, particles_temp) implemented as pointers in
67! order to speed up sorting (see routine sort_particles)
68!
69! 759 2011-09-15 13:58:31Z raasch
70! Splitting of parallel I/O (routine write_particles)
71!
72! Revision 1.1  1999/11/25 16:16:06  raasch
73! Initial revision
74!
75!
76! Description:
77! ------------
78! Particle advection
79!------------------------------------------------------------------------------!
80
81    USE arrays_3d
82    USE control_parameters
83    USE cpulog
84    USE particle_attributes
85    USE pegrid
86    USE statistics
87
88    IMPLICIT NONE
89
90    INTEGER ::  m
91
92
93    CALL cpu_log( log_point(25), 'lpm', 'start' )
94
95!
96!-- Write particle data at current time on file.
97!-- This has to be done here, before particles are further processed,
98!-- because they may be deleted within this timestep (in case that
99!-- dt_write_particle_data = dt_prel = particle_maximum_age).
100    time_write_particle_data = time_write_particle_data + dt_3d
101    IF ( time_write_particle_data >= dt_write_particle_data )  THEN
102
103       CALL lpm_data_output_particles
104!
105!--    The MOD function allows for changes in the output interval with restart
106!--    runs.
107       time_write_particle_data = MOD( time_write_particle_data, &
108                                  MAX( dt_write_particle_data, dt_3d ) )
109    ENDIF
110
111
112!
113!-- Initialize arrays for marking those particles/tails to be deleted after the
114!-- (sub-) timestep
115    particle_mask     = .TRUE.
116    deleted_particles = 0
117
118    IF ( use_particle_tails )  THEN
119       tail_mask = .TRUE.
120    ENDIF
121    deleted_tails = 0
122
123
124!
125!-- Initialize variables used for accumulating the number of particles
126!-- exchanged between the subdomains during all sub-timesteps (if sgs
127!-- velocities are included). These data are output further below on the
128!-- particle statistics file.
129    trlp_count_sum      = 0
130    trlp_count_recv_sum = 0
131    trrp_count_sum      = 0
132    trrp_count_recv_sum = 0
133    trsp_count_sum      = 0
134    trsp_count_recv_sum = 0
135    trnp_count_sum      = 0
136    trnp_count_recv_sum = 0
137
138
139!
140!-- Calculate exponential term used in case of particle inertia for each
141!-- of the particle groups
142    DO  m = 1, number_of_particle_groups
143       IF ( particle_groups(m)%density_ratio /= 0.0 )  THEN
144          particle_groups(m)%exp_arg  =                                        &
145                    4.5 * particle_groups(m)%density_ratio *                   &
146                    molecular_viscosity / ( particle_groups(m)%radius )**2
147          particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
148                                             dt_3d )
149       ENDIF
150    ENDDO
151
152
153!
154!-- Particle (droplet) growth by condensation/evaporation and collision
155    IF ( cloud_droplets )  THEN
156
157!
158!--    Reset summation arrays
159       ql_c = 0.0;  ql_v = 0.0;  ql_vp = 0.0
160
161!
162!--    Droplet growth by condensation / evaporation
163       CALL lpm_droplet_condensation
164
165!
166!--    Particle growth by collision
167       IF ( collision_kernel /= 'none' )  CALL lpm_droplet_collision
168
169    ENDIF
170
171
172!
173!-- If particle advection includes SGS velocity components, calculate the
174!-- required SGS quantities (i.e. gradients of the TKE, as well as horizontally
175!-- averaged profiles of the SGS TKE and the resolved-scale velocity variances)
176    IF ( use_sgs_for_particles )  CALL lpm_init_sgs_tke
177
178
179!
180!-- Initialize the variable storing the total time that a particle has advanced
181!-- within the timestep procedure
182    particles(1:number_of_particles)%dt_sum = 0.0
183
184
185!
186!-- Timestep loop for particle advection.
187!-- This loop has to be repeated until the advection time of every particle
188!-- (within the total domain!) has reached the LES timestep (dt_3d).
189!-- In case of including the SGS velocities, the particle timestep may be
190!-- smaller than the LES timestep (because of the Lagrangian timescale restric-
191!-- tion) and particles may require to undergo several particle timesteps,
192!-- before the LES timestep is reached. Because the number of these particle
193!-- timesteps to be carried out is unknown at first, these steps are carried
194!-- out in the following infinite loop with exit condition.
195    DO
196
197       CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' )
198
199!
200!--    Initialize the switch used for the loop exit condition checked at the
201!--    end of this loop.
202!--    If at least one particle has failed to reach the LES timestep, this
203!--    switch will be set false.
204       dt_3d_reached_l = .TRUE.
205
206!
207!--    Particle advection
208       CALL lpm_advec
209
210!
211!--    Particle reflection from walls
212       CALL lpm_boundary_conds( 'walls' )
213
214!
215!--    User-defined actions after the calculation of the new particle position
216       CALL user_lpm_advec
217
218!
219!--    Find out, if all particles on every PE have completed the LES timestep
220!--    and set the switch corespondingly
221#if defined( __parallel )
222       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
223       CALL MPI_ALLREDUCE( dt_3d_reached_l, dt_3d_reached, 1, MPI_LOGICAL, &
224                           MPI_LAND, comm2d, ierr )
225#else
226       dt_3d_reached = dt_3d_reached_l
227#endif
228
229       CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' )
230
231!
232!--    Increment time since last release
233       IF ( dt_3d_reached )  time_prel = time_prel + dt_3d
234
235!
236!--    If necessary, release new set of particles
237       IF ( time_prel >= dt_prel  .AND.  end_time_prel > simulated_time  .AND. &
238            dt_3d_reached )  THEN
239
240          CALL lpm_release_set
241
242!
243!--       The MOD function allows for changes in the output interval with
244!--       restart runs.
245          time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
246
247       ENDIF
248
249!
250!--    Horizontal boundary conditions including exchange between subdmains
251       CALL lpm_exchange_horiz
252
253!
254!--    Apply boundary conditions to those particles that have crossed the top or
255!--    bottom boundary and delete those particles, which are older than allowed
256       CALL lpm_boundary_conds( 'bottom/top' )
257
258!
259!--    Pack particles (eliminate those marked for deletion),
260!--    determine new number of particles
261       IF ( number_of_particles > 0  .AND.  deleted_particles > 0 )  THEN
262          CALL lpm_pack_arrays
263       ENDIF
264
265!
266!--    Initialize variables for the next (sub-) timestep, i.e. for marking those
267!--    particles to be deleted after the timestep
268       particle_mask     = .TRUE.
269       deleted_particles = 0
270
271       IF ( use_particle_tails )  THEN
272          tail_mask     = .TRUE.
273       ENDIF
274       deleted_tails = 0
275
276
277       IF ( dt_3d_reached )  EXIT
278
279    ENDDO   ! timestep loop
280
281
282!
283!-- Sort particles in the sequence the gridboxes are stored in the memory
284    time_sort_particles = time_sort_particles + dt_3d
285    IF ( time_sort_particles >= dt_sort_particles )  THEN
286       CALL lpm_sort_arrays
287       time_sort_particles = MOD( time_sort_particles, &
288                                  MAX( dt_sort_particles, dt_3d ) )
289    ENDIF
290
291
292!
293!-- Calculate the new liquid water content for each grid box
294    IF ( cloud_droplets )  CALL lpm_calc_liquid_water_content
295
296
297!
298!-- Set particle attributes.
299!-- Feature is not available if collision is activated, because the respective
300!-- particle attribute (class) is then used for storing the particle radius
301!-- class.
302    IF ( collision_kernel == 'none' )  CALL lpm_set_attributes
303
304
305!
306!-- Set particle attributes defined by the user
307    CALL user_lpm_set_attributes
308
309
310!
311!-- If required, add the current particle positions to the particle tails
312    IF ( use_particle_tails )  CALL lpm_extend_tails
313
314
315!
316!-- Write particle statistics (in particular the number of particles
317!-- exchanged between the subdomains) on file
318    IF ( write_particle_statistics )  CALL lpm_write_exchange_statistics
319
320    CALL cpu_log( log_point(25), 'lpm', 'stop' )
321
322
323 END SUBROUTINE lpm
Note: See TracBrowser for help on using the repository browser.