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

Last change on this file since 3633 was 3560, checked in by raasch, 6 years ago

bugfix to guarantee correct particle releases in case that the release interval is smaller than the model timestep

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