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

Last change on this file since 1834 was 1823, checked in by hoffmann, 5 years ago

last commit documented

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