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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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