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

Last change on this file since 2320 was 2263, checked in by schwenkel, 7 years ago

Implemented splitting and merging algorithm

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