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

Last change on this file since 1267 was 1037, checked in by raasch, 12 years ago

last commit documented

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