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

Last change on this file since 1936 was 1936, checked in by suehring, 5 years ago

deallocation of unused particle memory, formatting adjustments

  • Property svn:keywords set to Id
File size: 14.9 KB
Line 
1!> @file lpm.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21! Call routine for deallocation of unused memory.
22! Formatting adjustments
23!
24! Former revisions:
25! -----------------
26! $Id: lpm.f90 1936 2016-06-13 13:37:44Z suehring $
27!
28! 1929 2016-06-09 16:25:25Z suehring
29! Call wall boundary conditions only if particles are in the vertical range of
30! topography.
31!
32! 1822 2016-04-07 07:49:42Z hoffmann
33! Tails removed.
34!
35! Initialization of sgs model not necessary for the use of cloud_droplets and
36! use_sgs_for_particles.
37!
38! lpm_release_set integrated.
39!
40! Unused variabled removed.
41!
42! 1682 2015-10-07 23:56:08Z knoop
43! Code annotations made doxygen readable
44!
45! 1416 2014-06-04 16:04:03Z suehring
46! user_lpm_advec is called for each gridpoint.
47! Bugfix: in order to prevent an infinite loop, time_loop_done is set .TRUE.
48! at the head of the do-loop. 
49!
50! 1359 2014-04-11 17:15:14Z hoffmann
51! New particle structure integrated.
52! Kind definition added to all floating point numbers.
53!
54! 1320 2014-03-20 08:40:49Z raasch
55! ONLY-attribute added to USE-statements,
56! kind-parameters added to all INTEGER and REAL declaration statements,
57! kinds are defined in new module kinds,
58! revision history before 2012 removed,
59! comment fields (!:) to be used for variable explanations added to
60! all variable declaration statements
61!
62! 1318 2014-03-17 13:35:16Z raasch
63! module interfaces removed
64!
65! 1036 2012-10-22 13:43:42Z raasch
66! code put under GPL (PALM 3.9)
67!
68! 851 2012-03-15 14:32:58Z raasch
69! Bugfix: resetting of particle_mask and tail mask moved from routine
70! lpm_exchange_horiz to here (end of sub-timestep loop)
71!
72! 849 2012-03-15 10:35:09Z raasch
73! original routine advec_particles split into several subroutines and renamed
74! lpm
75!
76! 831 2012-02-22 00:29:39Z raasch
77! thermal_conductivity_l and diff_coeff_l now depend on temperature and
78! pressure
79!
80! 828 2012-02-21 12:00:36Z raasch
81! fast hall/wang kernels with fixed radius/dissipation classes added,
82! particle feature color renamed class, routine colker renamed
83! recalculate_kernel,
84! lower limit for droplet radius changed from 1E-7 to 1E-8
85!
86! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
87!
88! 825 2012-02-19 03:03:44Z raasch
89! droplet growth by condensation may include curvature and solution effects,
90! initialisation of temporary particle array for resorting removed,
91! particle attributes speed_x|y|z_sgs renamed rvar1|2|3,
92! module wang_kernel_mod renamed lpm_collision_kernels_mod,
93! wang_collision_kernel renamed wang_kernel
94!
95!
96! Revision 1.1  1999/11/25 16:16:06  raasch
97! Initial revision
98!
99!
100! Description:
101! ------------
102!> Particle advection
103!------------------------------------------------------------------------------!
104 SUBROUTINE lpm
105 
106
107    USE arrays_3d,                                                             &
108        ONLY:  ql_c, ql_v, ql_vp
109
110    USE control_parameters,                                                    &
111        ONLY:  cloud_droplets, dt_3d, dt_3d_reached, dt_3d_reached_l,          &
112               molecular_viscosity, simulated_time, topography
113
114    USE cpulog,                                                                &
115        ONLY:  cpu_log, log_point, log_point_s
116
117    USE indices,                                                               &
118        ONLY: nxl, nxr, nys, nyn, nzb, nzb_max, nzt, nzb_w_inner
119
120    USE kinds
121
122    USE lpm_exchange_horiz_mod,                                                &
123        ONLY:  dealloc_particles_array, lpm_exchange_horiz, lpm_move_particle
124
125    USE lpm_init_mod,                                                          &
126        ONLY: lpm_create_particle, PHASE_RELEASE
127
128    USE lpm_pack_arrays_mod,                                                   &
129        ONLY:  lpm_pack_all_arrays
130
131    USE particle_attributes,                                                   &
132        ONLY:  collision_kernel, deleted_particles, deallocate_memory,         &
133               dt_write_particle_data, dt_prel, end_time_prel,                 &
134               grid_particles, number_of_particles, number_of_particle_groups, &
135               particles, particle_groups, prt_count, step_dealloc,            &
136               time_prel, time_write_particle_data, trlp_count_sum,            &
137               trlp_count_recv_sum, trnp_count_sum,                            &
138               trnp_count_recv_sum, trrp_count_sum, trrp_count_recv_sum,       &
139               trsp_count_sum, trsp_count_recv_sum,                            &
140               use_sgs_for_particles, write_particle_statistics
141
142    USE pegrid
143
144    IMPLICIT NONE
145
146    INTEGER(iwp)       ::  i                  !<
147    INTEGER(iwp)       ::  ie                 !<
148    INTEGER(iwp)       ::  is                 !<
149    INTEGER(iwp)       ::  j                  !<
150    INTEGER(iwp)       ::  je                 !<
151    INTEGER(iwp)       ::  js                 !<
152    INTEGER(iwp), SAVE ::  lpm_count = 0      !<
153    INTEGER(iwp)       ::  k                  !<
154    INTEGER(iwp)       ::  ke                 !<
155    INTEGER(iwp)       ::  ks                 !<
156    INTEGER(iwp)       ::  m                  !<
157    INTEGER(iwp), SAVE ::  steps = 0          !<
158
159    LOGICAL            ::  first_loop_stride  !<
160
161    CALL cpu_log( log_point(25), 'lpm', 'start' )
162
163!
164!-- Write particle data at current time on file.
165!-- This has to be done here, before particles are further processed,
166!-- because they may be deleted within this timestep (in case that
167!-- dt_write_particle_data = dt_prel = particle_maximum_age).
168    time_write_particle_data = time_write_particle_data + dt_3d
169    IF ( time_write_particle_data >= dt_write_particle_data )  THEN
170
171       CALL lpm_data_output_particles
172!
173!--    The MOD function allows for changes in the output interval with restart
174!--    runs.
175       time_write_particle_data = MOD( time_write_particle_data, &
176                                  MAX( dt_write_particle_data, dt_3d ) )
177    ENDIF
178
179!
180!-- Initialize arrays for marking those particles to be deleted after the
181!-- (sub-) timestep
182    deleted_particles = 0
183
184!
185!-- Initialize variables used for accumulating the number of particles
186!-- exchanged between the subdomains during all sub-timesteps (if sgs
187!-- velocities are included). These data are output further below on the
188!-- particle statistics file.
189    trlp_count_sum      = 0
190    trlp_count_recv_sum = 0
191    trrp_count_sum      = 0
192    trrp_count_recv_sum = 0
193    trsp_count_sum      = 0
194    trsp_count_recv_sum = 0
195    trnp_count_sum      = 0
196    trnp_count_recv_sum = 0
197
198
199!
200!-- Calculate exponential term used in case of particle inertia for each
201!-- of the particle groups
202    DO  m = 1, number_of_particle_groups
203       IF ( particle_groups(m)%density_ratio /= 0.0_wp )  THEN
204          particle_groups(m)%exp_arg  =                                        &
205                    4.5_wp * particle_groups(m)%density_ratio *                &
206                    molecular_viscosity / ( particle_groups(m)%radius )**2
207
208          particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg *     &
209                    dt_3d )
210       ENDIF
211    ENDDO
212
213!
214!-- If necessary, release new set of particles
215    IF ( time_prel >= dt_prel  .AND.  end_time_prel > simulated_time )  THEN
216
217       CALL lpm_create_particle(PHASE_RELEASE)
218!
219!--    The MOD function allows for changes in the output interval with
220!--    restart runs.
221       time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
222
223    ENDIF
224!
225!-- Reset summation arrays
226    IF ( cloud_droplets)  THEN
227       ql_c  = 0.0_wp
228       ql_v  = 0.0_wp
229       ql_vp = 0.0_wp
230    ENDIF
231
232    first_loop_stride = .TRUE.
233    grid_particles(:,:,:)%time_loop_done = .TRUE.
234!
235!-- Timestep loop for particle advection.
236!-- This loop has to be repeated until the advection time of every particle
237!-- (within the total domain!) has reached the LES timestep (dt_3d).
238!-- In case of including the SGS velocities, the particle timestep may be
239!-- smaller than the LES timestep (because of the Lagrangian timescale
240!-- restriction) and particles may require to undergo several particle
241!-- timesteps, before the LES timestep is reached. Because the number of these
242!-- particle timesteps to be carried out is unknown at first, these steps are
243!-- carried out in the following infinite loop with exit condition.
244    DO
245       CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' )
246       CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
247       
248       grid_particles(:,:,:)%time_loop_done = .TRUE.
249!
250!--    If particle advection includes SGS velocity components, calculate the
251!--    required SGS quantities (i.e. gradients of the TKE, as well as
252!--    horizontally averaged profiles of the SGS TKE and the resolved-scale
253!--    velocity variances)
254
255       IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
256          CALL lpm_init_sgs_tke
257       ENDIF
258
259       DO  i = nxl, nxr
260          DO  j = nys, nyn
261             DO  k = nzb+1, nzt
262
263                number_of_particles = prt_count(k,j,i)
264!
265!--             If grid cell gets empty, flag must be true
266                IF ( number_of_particles <= 0 )  THEN
267                   grid_particles(k,j,i)%time_loop_done = .TRUE.
268                   CYCLE
269                ENDIF
270
271                IF ( .NOT. first_loop_stride  .AND.  &
272                     grid_particles(k,j,i)%time_loop_done ) CYCLE
273
274                particles => grid_particles(k,j,i)%particles(1:number_of_particles)
275
276                particles(1:number_of_particles)%particle_mask = .TRUE.
277!
278!--             Initialize the variable storing the total time that a particle
279!--             has advanced within the timestep procedure
280                IF ( first_loop_stride )  THEN
281                   particles(1:number_of_particles)%dt_sum = 0.0_wp
282                ENDIF
283!
284!--             Particle (droplet) growth by condensation/evaporation and
285!--             collision
286                IF ( cloud_droplets  .AND.  first_loop_stride)  THEN
287!
288!--                Droplet growth by condensation / evaporation
289                   CALL lpm_droplet_condensation(i,j,k)
290!
291!--                Particle growth by collision
292                   IF ( collision_kernel /= 'none' )  THEN
293                      CALL lpm_droplet_collision(i,j,k)
294                   ENDIF
295
296                ENDIF
297!
298!--             Initialize the switch used for the loop exit condition checked
299!--             at the end of this loop. If at least one particle has failed to
300!--             reach the LES timestep, this switch will be set false in
301!--             lpm_advec.
302                dt_3d_reached_l = .TRUE.
303
304!
305!--             Particle advection
306                CALL lpm_advec(i,j,k)
307!
308!--             Particle reflection from walls. Only applied if the particles
309!--             are in the vertical range of the topography. (Here, some
310!--             optimization is still possible.)
311                IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
312                   CALL lpm_boundary_conds( 'walls' )
313                ENDIF
314!
315!--             User-defined actions after the calculation of the new particle
316!--             position
317                CALL user_lpm_advec(i,j,k)
318!
319!--             Apply boundary conditions to those particles that have crossed
320!--             the top or bottom boundary and delete those particles, which are
321!--             older than allowed
322                CALL lpm_boundary_conds( 'bottom/top' )
323!
324!---            If not all particles of the actual grid cell have reached the
325!--             LES timestep, this cell has to to another loop iteration. Due to
326!--             the fact that particles can move into neighboring grid cell,
327!--             these neighbor cells also have to perform another loop iteration
328                IF ( .NOT. dt_3d_reached_l )  THEN
329                   ks = MAX(nzb+1,k)
330                   ke = MIN(nzt,k)
331                   js = MAX(nys,j)
332                   je = MIN(nyn,j)
333                   is = MAX(nxl,i)
334                   ie = MIN(nxr,i)
335                   grid_particles(ks:ke,js:je,is:ie)%time_loop_done = .FALSE.
336                ENDIF
337
338             ENDDO
339          ENDDO
340       ENDDO
341
342       steps = steps + 1
343       dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done)
344!
345!--    Find out, if all particles on every PE have completed the LES timestep
346!--    and set the switch corespondingly
347#if defined( __parallel )
348       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
349       CALL MPI_ALLREDUCE( dt_3d_reached_l, dt_3d_reached, 1, MPI_LOGICAL, &
350                           MPI_LAND, comm2d, ierr )
351#else
352       dt_3d_reached = dt_3d_reached_l
353#endif
354
355       CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' )
356
357!
358!--    Increment time since last release
359       IF ( dt_3d_reached )  time_prel = time_prel + dt_3d
360
361!
362!--    Move Particles local to PE to a different grid cell
363       CALL lpm_move_particle
364
365!
366!--    Horizontal boundary conditions including exchange between subdmains
367       CALL lpm_exchange_horiz
368!
369!--    Pack particles (eliminate those marked for deletion),
370!--    determine new number of particles
371       CALL lpm_pack_all_arrays
372
373!
374!--    Initialize variables for the next (sub-) timestep, i.e., for marking
375!--    those particles to be deleted after the timestep
376       deleted_particles = 0
377
378       IF ( dt_3d_reached )  EXIT
379
380       first_loop_stride = .FALSE.
381    ENDDO   ! timestep loop
382
383!
384!-- Calculate the new liquid water content for each grid box
385    IF ( cloud_droplets )  CALL lpm_calc_liquid_water_content
386!
387!-- Deallocate unused memory
388    IF ( deallocate_memory  .AND.  lpm_count == step_dealloc )  THEN
389       CALL dealloc_particles_array
390       lpm_count = 0
391    ELSEIF ( deallocate_memory )  THEN
392       lpm_count = lpm_count + 1
393    ENDIF
394
395!
396!-- Set particle attributes.
397!-- Feature is not available if collision is activated, because the respective
398!-- particle attribute (class) is then used for storing the particle radius
399!-- class.
400    IF ( collision_kernel == 'none' )  CALL lpm_set_attributes
401
402!
403!-- Set particle attributes defined by the user
404    CALL user_lpm_set_attributes
405
406!
407!-- Write particle statistics (in particular the number of particles
408!-- exchanged between the subdomains) on file
409    IF ( write_particle_statistics )  CALL lpm_write_exchange_statistics
410
411    CALL cpu_log( log_point(25), 'lpm', 'stop' )
412
413 END SUBROUTINE lpm
Note: See TracBrowser for help on using the repository browser.