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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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