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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

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