source: palm/tags/release-6.0/SOURCE/lpm.f90 @ 3984

Last change on this file since 3984 was 2801, checked in by thiele, 6 years ago

Introduce particle transfer in nested models

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