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

Last change on this file since 1815 was 1683, checked in by knoop, 9 years ago

last commit documented

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