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

Last change on this file since 850 was 850, checked in by raasch, 10 years ago

last commit documented

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