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

Last change on this file since 1359 was 1359, checked in by hoffmann, 7 years ago

new Lagrangian particle structure integrated

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