source: palm/tags/release-4.0/SOURCE/lpm.f90 @ 3800

Last change on this file since 3800 was 1417, checked in by suehring, 10 years ago

last commit documented

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