source: palm/trunk/SOURCE/lpm_init.f90 @ 2225

Last change on this file since 2225 was 2224, checked in by suehring, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 45.5 KB
Line 
1!> @file lpm_init.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
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_init.f90 2224 2017-05-15 16:38:37Z raasch $
27!
28! 2223 2017-05-15 16:38:09Z suehring
29! Add check for particle release at model top
30!
31! 2182 2017-03-17 14:27:40Z schwenkel
32! Added parameters for simplified particle initialization.
33!
34! 2122 2017-01-18 12:22:54Z hoffmann
35! Improved initialization of equilibrium aerosol radii
36! Calculation of particle ID
37!
38! 2000 2016-08-20 18:09:15Z knoop
39! Forced header and separation lines into 80 columns
40!
41! 2016-06-09 16:25:25Z suehring
42! Bugfix in determining initial particle height and grid index in case of
43! seed_follows_topography.
44! Bugfix concerning random positions, ensure that particles do not move more
45! than one grid length.
46! Bugfix logarithmic interpolation.
47! Initial setting of sgs_wf_part.
48!
49! 1890 2016-04-22 08:52:11Z hoffmann
50! Initialization of aerosol equilibrium radius not possible in supersaturated
51! environments. Therefore, a maximum supersaturation of -1 % is assumed during
52! initialization.
53!
54! 1873 2016-04-18 14:50:06Z maronga
55! Module renamed (removed _mod
56!
57! 1871 2016-04-15 11:46:09Z hoffmann
58! Initialization of aerosols added.
59!
60! 1850 2016-04-08 13:29:27Z maronga
61! Module renamed
62!
63! 1831 2016-04-07 13:15:51Z hoffmann
64! curvature_solution_effects moved to particle_attributes
65!
66! 1822 2016-04-07 07:49:42Z hoffmann
67! Unused variables removed.
68!
69! 1783 2016-03-06 18:36:17Z raasch
70! netcdf module added
71!
72! 1725 2015-11-17 13:01:51Z hoffmann
73! Bugfix: Processor-dependent seed for random function is generated before it is
74! used.
75!
76! 1691 2015-10-26 16:17:44Z maronga
77! Renamed prandtl_layer to constant_flux_layer.
78!
79! 1685 2015-10-08 07:32:13Z raasch
80! bugfix concerning vertical index offset in case of ocean
81!
82! 1682 2015-10-07 23:56:08Z knoop
83! Code annotations made doxygen readable
84!
85! 1575 2015-03-27 09:56:27Z raasch
86! initial vertical particle position is allowed to follow the topography
87!
88! 1359 2014-04-11 17:15:14Z hoffmann
89! New particle structure integrated.
90! Kind definition added to all floating point numbers.
91! lpm_init changed form a subroutine to a module.
92!
93! 1327 2014-03-21 11:00:16Z raasch
94! -netcdf_output
95!
96! 1322 2014-03-20 16:38:49Z raasch
97! REAL functions provided with KIND-attribute
98!
99! 1320 2014-03-20 08:40:49Z raasch
100! ONLY-attribute added to USE-statements,
101! kind-parameters added to all INTEGER and REAL declaration statements,
102! kinds are defined in new module kinds,
103! revision history before 2012 removed,
104! comment fields (!:) to be used for variable explanations added to
105! all variable declaration statements
106! bugfix: #if defined( __parallel ) added
107!
108! 1314 2014-03-14 18:25:17Z suehring
109! Vertical logarithmic interpolation of horizontal particle speed for particles
110! between roughness height and first vertical grid level.
111!
112! 1092 2013-02-02 11:24:22Z raasch
113! unused variables removed
114!
115! 1036 2012-10-22 13:43:42Z raasch
116! code put under GPL (PALM 3.9)
117!
118! 849 2012-03-15 10:35:09Z raasch
119! routine renamed: init_particles -> lpm_init
120! de_dx, de_dy, de_dz are allocated here (instead of automatic arrays in
121! advec_particles),
122! sort_particles renamed lpm_sort_arrays, user_init_particles renamed lpm_init
123!
124! 828 2012-02-21 12:00:36Z raasch
125! call of init_kernels, particle feature color renamed class
126!
127! 824 2012-02-17 09:09:57Z raasch
128! particle attributes speed_x|y|z_sgs renamed rvar1|2|3,
129! array particles implemented as pointer
130!
131! 667 2010-12-23 12:06:00Z suehring/gryschka
132! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
133! of arrays.
134!
135! Revision 1.1  1999/11/25 16:22:38  raasch
136! Initial revision
137!
138!
139! Description:
140! ------------
141!> This routine initializes a set of particles and their attributes (position,
142!> radius, ..) which are used by the Lagrangian particle model (see lpm).
143!------------------------------------------------------------------------------!
144 MODULE lpm_init_mod
145 
146
147    USE arrays_3d,                                                             &
148        ONLY:  de_dx, de_dy, de_dz, zu, zw, z0
149
150    USE control_parameters,                                                    &
151        ONLY:  cloud_droplets, constant_flux_layer, current_timestep_number,   &
152               dz, initializing_actions, message_string, ocean, simulated_time
153
154    USE grid_variables,                                                        &
155        ONLY:  ddx, dx, ddy, dy
156
157    USE indices,                                                               &
158        ONLY:  nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb,    &
159               nzb_w_inner, nzt
160
161    USE kinds
162
163    USE lpm_collision_kernels_mod,                                             &
164        ONLY:  init_kernels
165
166    USE netcdf_interface,                                                      &
167        ONLY:  netcdf_data_format
168
169    USE particle_attributes,                                                   &
170        ONLY:   alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,        &
171                block_offset, block_offset_def, collision_kernel,              &
172                curvature_solution_effects,                                    &
173                density_ratio, grid_particles,                                 &
174                initial_weighting_factor, ibc_par_b, ibc_par_lr, ibc_par_ns,   &
175                ibc_par_t, iran_part, log_z_z0,                                &
176                max_number_of_particle_groups, maximum_number_of_particles,    &
177                min_nr_particle, mpi_particle_type,                            &
178                number_concentration, number_particles_per_gridbox,            &
179                number_of_particles,                                           &
180                number_of_particle_groups, number_of_sublayers,                &
181                offset_ocean_nzt, offset_ocean_nzt_m1,                         &
182                particles, particle_advection_start, particle_groups,          &
183                particle_groups_type, particles_per_point,                     &
184                particle_type, pdx, pdy, pdz,                                  &
185                prt_count, psb, psl, psn, psr, pss, pst,                       &
186                radius, random_start_position, read_particles_from_restartfile,&
187                seed_follows_topography, sgs_wf_part, sort_count,              &
188                total_number_of_particles,                                     &
189                use_sgs_for_particles,                                         &
190                write_particle_statistics, uniform_particles, zero_particle,   &
191                z0_av_global
192
193    USE pegrid
194
195    USE random_function_mod,                                                   &
196        ONLY:  random_function
197
198    IMPLICIT NONE
199
200    PRIVATE
201
202    INTEGER(iwp), PARAMETER         :: PHASE_INIT    = 1  !<
203    INTEGER(iwp), PARAMETER, PUBLIC :: PHASE_RELEASE = 2  !<
204
205    INTERFACE lpm_init
206       MODULE PROCEDURE lpm_init
207    END INTERFACE lpm_init
208
209    INTERFACE lpm_create_particle
210       MODULE PROCEDURE lpm_create_particle
211    END INTERFACE lpm_create_particle
212
213    PUBLIC lpm_init, lpm_create_particle
214
215 CONTAINS
216
217!------------------------------------------------------------------------------!
218! Description:
219! ------------
220!> @todo Missing subroutine description.
221!------------------------------------------------------------------------------!
222 SUBROUTINE lpm_init
223
224    USE lpm_collision_kernels_mod,                                             &
225        ONLY:  init_kernels
226
227    IMPLICIT NONE
228
229    INTEGER(iwp) ::  i                           !<
230    INTEGER(iwp) ::  j                           !<
231    INTEGER(iwp) ::  k                           !<
232
233#if defined( __parallel )
234    INTEGER(iwp), DIMENSION(3) ::  blocklengths  !<
235    INTEGER(iwp), DIMENSION(3) ::  displacements !<
236    INTEGER(iwp), DIMENSION(3) ::  types         !<
237#endif
238
239    REAL(wp) ::  div                             !<
240    REAL(wp) ::  height_int                      !<
241    REAL(wp) ::  height_p                        !<
242    REAL(wp) ::  z_p                             !<
243    REAL(wp) ::  z0_av_local                     !<
244
245#if defined( __parallel )
246!
247!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
248!-- particle_attributes). Integer length is 4 byte, Real is 8 byte
249    blocklengths(1)  = 19;  blocklengths(2)  =   6;  blocklengths(3)  =   1
250    displacements(1) =  0;  displacements(2) = 152;  displacements(3) = 176
251
252    types(1) = MPI_REAL
253    types(2) = MPI_INTEGER
254    types(3) = MPI_UB
255    CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
256                          mpi_particle_type, ierr )
257    CALL MPI_TYPE_COMMIT( mpi_particle_type, ierr )
258#endif
259
260!
261!-- In case of oceans runs, the vertical index calculations need an offset,
262!-- because otherwise the k indices will become negative
263    IF ( ocean )  THEN
264       offset_ocean_nzt    = nzt
265       offset_ocean_nzt_m1 = nzt - 1
266    ENDIF
267
268!
269!-- Define block offsets for dividing a gridcell in 8 sub cells
270
271    block_offset(0) = block_offset_def (-1,-1,-1)
272    block_offset(1) = block_offset_def (-1,-1, 0)
273    block_offset(2) = block_offset_def (-1, 0,-1)
274    block_offset(3) = block_offset_def (-1, 0, 0)
275    block_offset(4) = block_offset_def ( 0,-1,-1)
276    block_offset(5) = block_offset_def ( 0,-1, 0)
277    block_offset(6) = block_offset_def ( 0, 0,-1)
278    block_offset(7) = block_offset_def ( 0, 0, 0)
279!
280!-- Check the number of particle groups.
281    IF ( number_of_particle_groups > max_number_of_particle_groups )  THEN
282       WRITE( message_string, * ) 'max_number_of_particle_groups =',      &
283                                  max_number_of_particle_groups ,         &
284                                  '&number_of_particle_groups reset to ', &
285                                  max_number_of_particle_groups
286       CALL message( 'lpm_init', 'PA0213', 0, 1, 0, 6, 0 )
287       number_of_particle_groups = max_number_of_particle_groups
288    ENDIF
289
290!
291!-- Set default start positions, if necessary
292    IF ( psl(1) == 9999999.9_wp )  psl(1) = -0.5_wp * dx
293    IF ( psr(1) == 9999999.9_wp )  psr(1) = ( nx + 0.5_wp ) * dx
294    IF ( pss(1) == 9999999.9_wp )  pss(1) = -0.5_wp * dy
295    IF ( psn(1) == 9999999.9_wp )  psn(1) = ( ny + 0.5_wp ) * dy
296    IF ( psb(1) == 9999999.9_wp )  psb(1) = zu(nz/2)
297    IF ( pst(1) == 9999999.9_wp )  pst(1) = psb(1)
298
299    IF ( pdx(1) == 9999999.9_wp  .OR.  pdx(1) == 0.0_wp )  pdx(1) = dx
300    IF ( pdy(1) == 9999999.9_wp  .OR.  pdy(1) == 0.0_wp )  pdy(1) = dy
301    IF ( pdz(1) == 9999999.9_wp  .OR.  pdz(1) == 0.0_wp )  pdz(1) = zu(2) - zu(1)
302
303!
304!-- If number_particles_per_gridbox is set, the parametres pdx, pdy and pdz are
305!-- calculated diagnostically. Therfore an isotropic distribution is prescribed.
306    IF ( number_particles_per_gridbox /= -1 .AND.   & 
307         number_particles_per_gridbox >= 1 )    THEN
308       pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) /  & 
309             REAL(number_particles_per_gridbox))**0.3333333_wp
310!
311!--    Ensure a smooth value (two significant digits) of distance between
312!--    particles (pdx, pdy, pdz).
313       div = 1000.0_wp
314       DO  WHILE ( pdx(1) < div )
315          div = div / 10.0_wp
316       ENDDO
317       pdx(1) = NINT( pdx(1) * 100.0_wp / div ) * div / 100.0_wp
318       pdy(1) = pdx(1)
319       pdz(1) = pdx(1)
320
321    ENDIF
322
323    DO  j = 2, number_of_particle_groups
324       IF ( psl(j) == 9999999.9_wp )  psl(j) = psl(j-1)
325       IF ( psr(j) == 9999999.9_wp )  psr(j) = psr(j-1)
326       IF ( pss(j) == 9999999.9_wp )  pss(j) = pss(j-1)
327       IF ( psn(j) == 9999999.9_wp )  psn(j) = psn(j-1)
328       IF ( psb(j) == 9999999.9_wp )  psb(j) = psb(j-1)
329       IF ( pst(j) == 9999999.9_wp )  pst(j) = pst(j-1)
330       IF ( pdx(j) == 9999999.9_wp  .OR.  pdx(j) == 0.0_wp )  pdx(j) = pdx(j-1)
331       IF ( pdy(j) == 9999999.9_wp  .OR.  pdy(j) == 0.0_wp )  pdy(j) = pdy(j-1)
332       IF ( pdz(j) == 9999999.9_wp  .OR.  pdz(j) == 0.0_wp )  pdz(j) = pdz(j-1)
333    ENDDO
334
335!
336!-- Allocate arrays required for calculating particle SGS velocities.
337!-- Initialize prefactor required for stoachastic Weil equation.
338    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
339       ALLOCATE( de_dx(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
340                 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
341                 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
342
343       sgs_wf_part = 1.0_wp / 3.0_wp   
344    ENDIF
345
346!
347!-- Allocate array required for logarithmic vertical interpolation of
348!-- horizontal particle velocities between the surface and the first vertical
349!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
350!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
351!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
352!-- To obtain exact height levels of particles, linear interpolation is applied
353!-- (see lpm_advec.f90).
354    IF ( constant_flux_layer )  THEN
355       
356       ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 
357       z_p         = zu(nzb+1) - zw(nzb)
358
359!
360!--    Calculate horizontal mean value of z0 used for logartihmic
361!--    interpolation. Note: this is not exact for heterogeneous z0.
362!--    However, sensitivity studies showed that the effect is
363!--    negligible.
364       z0_av_local  = SUM( z0(nys:nyn,nxl:nxr) )
365       z0_av_global = 0.0_wp
366
367#if defined( __parallel )
368       CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, &
369                          comm2d, ierr )
370#else
371       z0_av_global = z0_av_local
372#endif
373
374       z0_av_global = z0_av_global  / ( ( ny + 1 ) * ( nx + 1 ) )
375!
376!--    Horizontal wind speed is zero below and at z0
377       log_z_z0(0) = 0.0_wp
378!
379!--    Calculate vertical depth of the sublayers
380       height_int  = ( z_p - z0_av_global ) / REAL( number_of_sublayers, KIND=wp )
381!
382!--    Precalculate LOG(z/z0)
383       height_p    = z0_av_global
384       DO  k = 1, number_of_sublayers
385
386          height_p    = height_p + height_int
387          log_z_z0(k) = LOG( height_p / z0_av_global )
388
389       ENDDO
390
391    ENDIF
392
393!
394!-- Check boundary condition and set internal variables
395    SELECT CASE ( bc_par_b )
396   
397       CASE ( 'absorb' )
398          ibc_par_b = 1
399
400       CASE ( 'reflect' )
401          ibc_par_b = 2
402         
403       CASE DEFAULT
404          WRITE( message_string, * )  'unknown boundary condition ',   &
405                                       'bc_par_b = "', TRIM( bc_par_b ), '"'
406          CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
407         
408    END SELECT
409    SELECT CASE ( bc_par_t )
410   
411       CASE ( 'absorb' )
412          ibc_par_t = 1
413
414       CASE ( 'reflect' )
415          ibc_par_t = 2
416         
417       CASE DEFAULT
418          WRITE( message_string, * ) 'unknown boundary condition ',   &
419                                     'bc_par_t = "', TRIM( bc_par_t ), '"'
420          CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
421         
422    END SELECT
423    SELECT CASE ( bc_par_lr )
424
425       CASE ( 'cyclic' )
426          ibc_par_lr = 0
427
428       CASE ( 'absorb' )
429          ibc_par_lr = 1
430
431       CASE ( 'reflect' )
432          ibc_par_lr = 2
433         
434       CASE DEFAULT
435          WRITE( message_string, * ) 'unknown boundary condition ',   &
436                                     'bc_par_lr = "', TRIM( bc_par_lr ), '"'
437          CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
438         
439    END SELECT
440    SELECT CASE ( bc_par_ns )
441
442       CASE ( 'cyclic' )
443          ibc_par_ns = 0
444
445       CASE ( 'absorb' )
446          ibc_par_ns = 1
447
448       CASE ( 'reflect' )
449          ibc_par_ns = 2
450         
451       CASE DEFAULT
452          WRITE( message_string, * ) 'unknown boundary condition ',   &
453                                     'bc_par_ns = "', TRIM( bc_par_ns ), '"'
454          CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
455         
456    END SELECT
457
458!
459!-- Initialize collision kernels
460    IF ( collision_kernel /= 'none' )  CALL init_kernels
461
462!
463!-- For the first model run of a possible job chain initialize the
464!-- particles, otherwise read the particle data from restart file.
465    IF ( TRIM( initializing_actions ) == 'read_restart_data'  &
466         .AND.  read_particles_from_restartfile )  THEN
467
468       CALL lpm_read_restart_file
469
470    ELSE
471
472!
473!--    Allocate particle arrays and set attributes of the initial set of
474!--    particles, which can be also periodically released at later times.
475       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
476                 grid_particles(nzb+1:nzt,nys:nyn,nxl:nxr) )
477
478       maximum_number_of_particles = 0
479       number_of_particles         = 0
480
481       sort_count = 0
482       prt_count  = 0
483
484!
485!--    initialize counter for particle IDs
486       grid_particles%id_counter = 0
487
488!
489!--    Initialize all particles with dummy values (otherwise errors may
490!--    occur within restart runs). The reason for this is still not clear
491!--    and may be presumably caused by errors in the respective user-interface.
492       zero_particle = particle_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
493                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
494                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
495                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
496                                      0, 0, 0, 0, .FALSE., -1 )
497
498       particle_groups = particle_groups_type( 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
499
500!
501!--    Set values for the density ratio and radius for all particle
502!--    groups, if necessary
503       IF ( density_ratio(1) == 9999999.9_wp )  density_ratio(1) = 0.0_wp
504       IF ( radius(1)        == 9999999.9_wp )  radius(1) = 0.0_wp
505       DO  i = 2, number_of_particle_groups
506          IF ( density_ratio(i) == 9999999.9_wp )  THEN
507             density_ratio(i) = density_ratio(i-1)
508          ENDIF
509          IF ( radius(i) == 9999999.9_wp )  radius(i) = radius(i-1)
510       ENDDO
511
512       DO  i = 1, number_of_particle_groups
513          IF ( density_ratio(i) /= 0.0_wp  .AND.  radius(i) == 0 )  THEN
514             WRITE( message_string, * ) 'particle group #', i, 'has a', &
515                                        'density ratio /= 0 but radius = 0'
516             CALL message( 'lpm_init', 'PA0215', 1, 2, 0, 6, 0 )
517          ENDIF
518          particle_groups(i)%density_ratio = density_ratio(i)
519          particle_groups(i)%radius        = radius(i)
520       ENDDO
521
522!
523!--    Set a seed value for the random number generator to be exclusively
524!--    used for the particle code. The generated random numbers should be
525!--    different on the different PEs.
526       iran_part = iran_part + myid
527
528       CALL lpm_create_particle (PHASE_INIT)
529!
530!--    User modification of initial particles
531       CALL user_lpm_init
532
533!
534!--    Open file for statistical informations about particle conditions
535       IF ( write_particle_statistics )  THEN
536          CALL check_open( 80 )
537          WRITE ( 80, 8000 )  current_timestep_number, simulated_time,         &
538                              number_of_particles,                             &
539                              maximum_number_of_particles
540          CALL close_file( 80 )
541       ENDIF
542
543    ENDIF
544
545!
546!-- To avoid programm abort, assign particles array to the local version of
547!-- first grid cell
548    number_of_particles = prt_count(nzb+1,nys,nxl)
549    particles => grid_particles(nzb+1,nys,nxl)%particles(1:number_of_particles)
550!
551!-- Formats
5528000 FORMAT (I6,1X,F7.2,4X,I10,71X,I10)
553
554 END SUBROUTINE lpm_init
555
556!------------------------------------------------------------------------------!
557! Description:
558! ------------
559!> @todo Missing subroutine description.
560!------------------------------------------------------------------------------!
561 SUBROUTINE lpm_create_particle (phase)
562
563    USE lpm_exchange_horiz_mod,                                                &
564        ONLY: lpm_exchange_horiz, lpm_move_particle, realloc_particles_array
565
566    USE lpm_pack_arrays_mod,                                                   &
567        ONLY: lpm_pack_all_arrays
568
569    USE particle_attributes,                                                   &
570        ONLY: deleted_particles, monodisperse_aerosols
571
572    IMPLICIT  NONE
573
574    INTEGER(iwp)               ::  alloc_size  !< relative increase of allocated memory for particles
575    INTEGER(iwp)               ::  i           !< loop variable ( particle groups )
576    INTEGER(iwp)               ::  ip          !< index variable along x
577    INTEGER(iwp)               ::  j           !< loop variable ( particles per point )
578    INTEGER(iwp)               ::  jp          !< index variable along y
579    INTEGER(iwp)               ::  kp          !< index variable along z
580    INTEGER(iwp)               ::  loop_stride !< loop variable for initialization
581    INTEGER(iwp)               ::  n           !< loop variable ( number of particles )
582    INTEGER(iwp)               ::  new_size    !< new size of allocated memory for particles
583
584    INTEGER(iwp), INTENT(IN)   ::  phase       !< mode of inititialization
585
586    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_count !< start address of new particle
587    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  local_start !< start address of new particle
588
589    LOGICAL                    ::  first_stride !< flag for initialization
590
591    REAL(wp)                   ::  pos_x      !< increment for particle position in x     
592    REAL(wp)                   ::  pos_y      !< increment for particle position in y 
593    REAL(wp)                   ::  pos_z      !< increment for particle position in z     
594    REAL(wp)                   ::  rand_contr !< dummy argument for random position
595
596    TYPE(particle_type),TARGET ::  tmp_particle !< temporary particle used for initialization
597
598!
599!-- Calculate particle positions and store particle attributes, if
600!-- particle is situated on this PE
601    DO  loop_stride = 1, 2
602       first_stride = (loop_stride == 1)
603       IF ( first_stride )   THEN
604          local_count = 0           ! count number of particles
605       ELSE
606          local_count = prt_count   ! Start address of new particles
607       ENDIF
608
609!
610!--    Calculate initial_weighting_factor diagnostically
611       IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN
612          initial_weighting_factor =  number_concentration * 1.0E6_wp *             & 
613                                      pdx(1) * pdy(1) * pdz(1) 
614       END IF
615
616       n = 0
617       DO  i = 1, number_of_particle_groups
618
619          pos_z = psb(i)
620
621          DO WHILE ( pos_z <= pst(i) )
622
623             IF ( pos_z >= 0.0_wp  .AND.  pos_z < zw(nzt) )  THEN
624
625
626                pos_y = pss(i)
627
628                DO WHILE ( pos_y <= psn(i) )
629
630                   IF ( pos_y >= ( nys - 0.5_wp ) * dy  .AND.                  &
631                        pos_y <  ( nyn + 0.5_wp ) * dy )  THEN
632
633                      pos_x = psl(i)
634
635               xloop: DO WHILE ( pos_x <= psr(i) )
636
637                         IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            & 
638                              pos_x <  ( nxr + 0.5_wp ) * dx )  THEN
639
640                            DO  j = 1, particles_per_point
641
642                               n = n + 1
643                               tmp_particle%x             = pos_x
644                               tmp_particle%y             = pos_y
645                               tmp_particle%z             = pos_z
646                               tmp_particle%age           = 0.0_wp
647                               tmp_particle%age_m         = 0.0_wp
648                               tmp_particle%dt_sum        = 0.0_wp
649                               tmp_particle%user          = 0.0_wp !unused, free for the user
650                               tmp_particle%e_m           = 0.0_wp
651                               IF ( curvature_solution_effects )  THEN
652!
653!--                               Initial values (internal timesteps, derivative)
654!--                               for Rosenbrock method
655                                  tmp_particle%rvar1      = 1.0E-6_wp     !last Rosenbrock timestep
656                                  tmp_particle%rvar2      = 0.1E-6_wp     !dry aerosol radius
657                                  tmp_particle%rvar3      = -9999999.9_wp !unused in this configuration
658                               ELSE
659!
660!--                               Initial values for SGS velocities
661                                  tmp_particle%rvar1      = 0.0_wp
662                                  tmp_particle%rvar2      = 0.0_wp
663                                  tmp_particle%rvar3      = 0.0_wp
664                               ENDIF
665                               tmp_particle%speed_x       = 0.0_wp
666                               tmp_particle%speed_y       = 0.0_wp
667                               tmp_particle%speed_z       = 0.0_wp
668                               tmp_particle%origin_x      = pos_x
669                               tmp_particle%origin_y      = pos_y
670                               tmp_particle%origin_z      = pos_z
671                               tmp_particle%radius        = particle_groups(i)%radius
672                               tmp_particle%weight_factor = initial_weighting_factor
673                               tmp_particle%class         = 1
674                               tmp_particle%group         = i
675                               tmp_particle%id1           = 0 
676                               tmp_particle%id2           = 0 
677                               tmp_particle%particle_mask = .TRUE.
678                               tmp_particle%block_nr      = -1
679!
680!--                            Determine the grid indices of the particle position
681                              ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
682                               jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
683                               kp = tmp_particle%z / dz + 1 + offset_ocean_nzt
684
685                               IF ( seed_follows_topography )  THEN
686!
687!--                               Particle height is given relative to topography
688                                  kp = kp + nzb_w_inner(jp,ip)
689                                  tmp_particle%z = tmp_particle%z +            &
690                                                         zw(nzb_w_inner(jp,ip))
691                                  IF ( kp > nzt )  THEN
692                                     pos_x = pos_x + pdx(i)
693                                     CYCLE xloop
694                                  ENDIF
695                               ELSEIF ( .NOT. seed_follows_topography .AND.    &
696                                        tmp_particle%z <= zw(nzb_w_inner(jp,ip)) )  THEN
697                                  pos_x = pos_x + pdx(i)
698                                  CYCLE xloop                               
699                               ENDIF
700
701                               local_count(kp,jp,ip) = local_count(kp,jp,ip) + 1
702
703                               IF ( .NOT. first_stride )  THEN
704                                  IF ( ip < nxl  .OR.  jp < nys  .OR.  kp < nzb+1 )  THEN
705                                     write(6,*) 'xl ',ip,jp,kp,nxl,nys,nzb+1
706                                  ENDIF
707                                  IF ( ip > nxr  .OR.  jp > nyn  .OR.  kp > nzt )  THEN
708                                     write(6,*) 'xu ',ip,jp,kp,nxr,nyn,nzt
709                                  ENDIF
710                                  grid_particles(kp,jp,ip)%particles(local_count(kp,jp,ip)) = tmp_particle
711
712                               ENDIF
713                            ENDDO
714
715                         ENDIF
716
717                         pos_x = pos_x + pdx(i)
718
719                      ENDDO xloop
720
721                   ENDIF
722
723                   pos_y = pos_y + pdy(i)
724
725                ENDDO
726
727             ENDIF
728
729             pos_z = pos_z + pdz(i)
730
731          ENDDO
732
733       ENDDO
734
735       IF ( first_stride )  THEN
736          DO  ip = nxl, nxr
737             DO  jp = nys, nyn
738                DO  kp = nzb+1, nzt
739                   IF ( phase == PHASE_INIT )  THEN
740                      IF ( local_count(kp,jp,ip) > 0 )  THEN
741                         alloc_size = MAX( INT( local_count(kp,jp,ip) *        &
742                            ( 1.0_wp + alloc_factor / 100.0_wp ) ),            &
743                            min_nr_particle )
744                      ELSE
745                         alloc_size = min_nr_particle
746                      ENDIF
747                      ALLOCATE(grid_particles(kp,jp,ip)%particles(1:alloc_size))
748                      DO  n = 1, alloc_size
749                         grid_particles(kp,jp,ip)%particles(n) = zero_particle
750                      ENDDO
751                   ELSEIF ( phase == PHASE_RELEASE )  THEN
752                      IF ( local_count(kp,jp,ip) > 0 )  THEN
753                         new_size   = local_count(kp,jp,ip) + prt_count(kp,jp,ip)
754                         alloc_size = MAX( INT( new_size * ( 1.0_wp +          &
755                            alloc_factor / 100.0_wp ) ), min_nr_particle )
756                         IF( alloc_size > SIZE( grid_particles(kp,jp,ip)%particles) )  THEN
757                            CALL realloc_particles_array(ip,jp,kp,alloc_size)
758                         ENDIF
759                      ENDIF
760                   ENDIF
761                ENDDO
762             ENDDO
763          ENDDO
764       ENDIF
765
766    ENDDO
767
768
769
770    local_start = prt_count+1
771    prt_count   = local_count
772
773!
774!-- Calculate particle IDs
775    DO  ip = nxl, nxr
776       DO  jp = nys, nyn
777          DO  kp = nzb+1, nzt
778             number_of_particles = prt_count(kp,jp,ip)
779             IF ( number_of_particles <= 0 )  CYCLE
780             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
781
782             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
783
784                particles(n)%id1 = 10000_iwp * grid_particles(kp,jp,ip)%id_counter + kp
785                particles(n)%id2 = 10000_iwp * jp + ip
786
787                grid_particles(kp,jp,ip)%id_counter =                          &
788                                         grid_particles(kp,jp,ip)%id_counter + 1
789
790             ENDDO
791
792          ENDDO
793       ENDDO
794    ENDDO
795
796!
797!-- Initialize aerosol background spectrum
798    IF ( curvature_solution_effects  .AND.  .NOT. monodisperse_aerosols )  THEN
799       CALL lpm_init_aerosols(local_start)
800    ENDIF
801
802!
803!-- Add random fluctuation to particle positions.
804    IF ( random_start_position )  THEN
805       DO  ip = nxl, nxr
806          DO  jp = nys, nyn
807             DO  kp = nzb+1, nzt
808                number_of_particles = prt_count(kp,jp,ip)
809                IF ( number_of_particles <= 0 )  CYCLE
810                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
811!
812!--             Move only new particles. Moreover, limit random fluctuation
813!--             in order to prevent that particles move more than one grid box,
814!--             which would lead to problems concerning particle exchange
815!--             between processors in case pdx/pdy are larger than dx/dy,
816!--             respectively. 
817                DO  n = local_start(kp,jp,ip), number_of_particles
818                   IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
819                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
820                                     pdx(particles(n)%group)
821                      particles(n)%x = particles(n)%x +                        &
822                              MERGE( rand_contr, SIGN( dx, rand_contr ), &
823                                     ABS( rand_contr ) < dx                    &
824                                   ) 
825                   ENDIF
826                   IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
827                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
828                                     pdy(particles(n)%group)
829                      particles(n)%y = particles(n)%y +                        &
830                              MERGE( rand_contr, SIGN( dy, rand_contr ), &
831                                     ABS( rand_contr ) < dy                    &
832                                   ) 
833                   ENDIF
834                   IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
835                      rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
836                                     pdz(particles(n)%group)
837                      particles(n)%z = particles(n)%z +                        &
838                              MERGE( rand_contr, SIGN( dz, rand_contr ), &
839                                     ABS( rand_contr ) < dz                    &
840                                   ) 
841                   ENDIF
842                ENDDO
843!
844!--             Identify particles located outside the model domain and reflect
845!--             or absorb them if necessary.
846                CALL lpm_boundary_conds( 'bottom/top' )
847!
848!--             Furthermore, remove particles located in topography. Note, as
849!--             the particle speed is still zero at this point, wall
850!--             reflection boundary conditions will not work in this case.
851                particles =>                                                   &
852                       grid_particles(kp,jp,ip)%particles(1:number_of_particles)
853                DO  n = local_start(kp,jp,ip), number_of_particles
854                   i = ( particles(n)%x + 0.5_wp * dx ) * ddx
855                   j = ( particles(n)%y + 0.5_wp * dy ) * ddy
856                   IF ( particles(n)%z <= zw(nzb_w_inner(j,i)) )  THEN
857                      particles(n)%particle_mask = .FALSE.
858                      deleted_particles = deleted_particles + 1
859                   ENDIF
860                ENDDO
861             ENDDO
862          ENDDO
863       ENDDO
864!
865!--    Exchange particles between grid cells and processors
866       CALL lpm_move_particle
867       CALL lpm_exchange_horiz
868
869    ENDIF
870!
871!-- In case of random_start_position, delete particles identified by
872!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks,
873!-- which is needed for a fast interpolation of the LES fields on the particle
874!-- position.
875    CALL lpm_pack_all_arrays
876
877!
878!-- Determine maximum number of particles (i.e., all possible particles that
879!-- have been allocated) and the current number of particles
880    DO  ip = nxl, nxr
881       DO  jp = nys, nyn
882          DO  kp = nzb+1, nzt
883             maximum_number_of_particles = maximum_number_of_particles         &
884                                           + SIZE(grid_particles(kp,jp,ip)%particles)
885             number_of_particles         = number_of_particles                 &
886                                           + prt_count(kp,jp,ip)
887          ENDDO
888       ENDDO
889    ENDDO
890!
891!-- Calculate the number of particles of the total domain
892#if defined( __parallel )
893    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
894    CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
895    MPI_INTEGER, MPI_SUM, comm2d, ierr )
896#else
897    total_number_of_particles = number_of_particles
898#endif
899
900    RETURN
901
902 END SUBROUTINE lpm_create_particle
903
904 SUBROUTINE lpm_init_aerosols(local_start)
905
906    USE arrays_3d,                                                             &
907        ONLY: hyp, pt, q 
908
909    USE cloud_parameters,                                                      &
910        ONLY: l_d_rv, rho_l, r_v
911
912    USE constants,                                                             &
913        ONLY: pi
914
915    USE kinds
916
917    USE particle_attributes,                                                   &
918        ONLY: init_aerosol_probabilistic, molecular_weight_of_solute,          &
919              molecular_weight_of_water, n1, n2, n3, rho_s, rm1, rm2, rm3,     &
920              s1, s2, s3, vanthoff
921
922    IMPLICIT NONE
923
924    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cdf     !< CDF of aerosol spectrum
925    REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_temp  !< dry aerosol radius spectrum
926
927    REAL(wp)  :: afactor            !< curvature effects
928    REAL(wp)  :: bfactor            !< solute effects
929    REAL(wp)  :: dr                 !< width of radius bin
930    REAL(wp)  :: e_a                !< vapor pressure
931    REAL(wp)  :: e_s                !< saturation vapor pressure
932    REAL(wp)  :: n_init             !< sum of all aerosol concentrations
933    REAL(wp)  :: pdf                !< PDF of aerosol spectrum
934    REAL(wp)  :: rmin = 1.0e-8_wp   !< minimum aerosol radius
935    REAL(wp)  :: rmax = 1.0e-6_wp   !< maximum aerosol radius
936    REAL(wp)  :: rs_rand            !< random number
937    REAL(wp)  :: r_mid              !< mean radius
938    REAL(wp)  :: sigma              !< surface tension
939    REAL(wp)  :: t_int              !< temperature
940    REAL(wp)  :: weight_sum         !< sum of all weighting factors
941
942    INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  local_start !<
943
944    INTEGER(iwp)  :: n              !<
945    INTEGER(iwp)  :: nn             !<
946    INTEGER(iwp)  :: no_bins = 999  !< number of bins
947    INTEGER(iwp)  :: ip             !<
948    INTEGER(iwp)  :: jp             !<
949    INTEGER(iwp)  :: kp             !<
950
951    LOGICAL ::  new_pdf = .FALSE.   !< check if aerosol PDF has to be recalculated
952
953!
954!-- Compute aerosol background distribution
955    IF ( init_aerosol_probabilistic )  THEN
956       ALLOCATE( cdf(0:no_bins), r_temp(0:no_bins) )
957       DO n = 0, no_bins
958          r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) /            &
959                           REAL(no_bins, KIND=wp) * REAL(n, KIND=wp) )
960
961          cdf(n) = 0.0_wp
962          n_init = n1 + n2 + n3
963          IF ( n1 > 0.0_wp )  THEN
964             cdf(n) = cdf(n) + n1 / n_init * ( 0.5_wp + 0.5_wp *        &
965                                  ERF( LOG( r_temp(n) / rm1 ) /         &
966                                       ( SQRT(2.0_wp) * LOG(s1) )       &
967                                     ) )
968          ENDIF
969          IF ( n2 > 0.0_wp )  THEN
970             cdf(n) = cdf(n) + n2 / n_init * ( 0.5_wp + 0.5_wp *        &
971                                  ERF( LOG( r_temp(n) / rm2 ) /         &
972                                       ( SQRT(2.0_wp) * LOG(s2) )       &
973                                     ) )
974          ENDIF
975          IF ( n3 > 0.0_wp )  THEN
976             cdf(n) = cdf(n) + n3 / n_init * ( 0.5_wp + 0.5_wp *        &
977                                  ERF( LOG( r_temp(n) / rm3 ) /         &
978                                       ( SQRT(2.0_wp) * LOG(s3) )       &
979                                     ) )
980          ENDIF
981
982       ENDDO
983    ENDIF
984
985    DO  ip = nxl, nxr
986       DO  jp = nys, nyn
987          DO  kp = nzb+1, nzt
988
989             number_of_particles = prt_count(kp,jp,ip)
990             IF ( number_of_particles <= 0 )  CYCLE
991             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
992!
993!--          Initialize the aerosols with a predefined spectral distribution
994!--          of the dry radius (logarithmically increasing bins) and a varying
995!--          weighting factor
996             IF ( .NOT. init_aerosol_probabilistic )  THEN
997
998                new_pdf = .FALSE.
999                IF ( .NOT. ALLOCATED( r_temp ) )  THEN
1000                   new_pdf = .TRUE.
1001                ELSE
1002                   IF ( SIZE( r_temp ) .NE. &
1003                        number_of_particles - local_start(kp,jp,ip) + 2 )  THEN
1004                      new_pdf = .TRUE.
1005                      DEALLOCATE( r_temp )
1006                   ENDIF
1007                ENDIF
1008
1009                IF ( new_pdf )  THEN
1010
1011                   no_bins = number_of_particles + 1 - local_start(kp,jp,ip)
1012                   ALLOCATE( r_temp(0:no_bins) )
1013
1014                   DO n = 0, no_bins
1015                      r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / &
1016                                       REAL(no_bins, KIND=wp) *                 &
1017                                       REAL(n, KIND=wp) )
1018                   ENDDO
1019
1020                ENDIF
1021
1022!
1023!--             Calculate radius and concentration of each aerosol
1024                DO n = local_start(kp,jp,ip), number_of_particles
1025
1026                   nn = n - local_start(kp,jp,ip)
1027
1028                   r_mid = SQRT( r_temp(nn) * r_temp(nn+1) )
1029                   dr    = r_temp(nn+1) - r_temp(nn)
1030
1031                   pdf    = 0.0_wp
1032                   n_init = n1 + n2 + n3
1033                   IF ( n1 > 0.0_wp )  THEN
1034                      pdf = pdf + n1 / n_init * ( 1.0_wp / ( r_mid * LOG(s1) *      &
1035                                                             SQRT( 2.0_wp * pi )    &
1036                                                           ) *                      &
1037                                                  EXP( -( LOG( r_mid / rm1 ) )**2 / &
1038                                                       ( 2.0_wp * LOG(s1)**2 )      &
1039                                                     )                              &
1040                                                )
1041                   ENDIF
1042                   IF ( n2 > 0.0_wp )  THEN
1043                      pdf = pdf + n2 / n_init * ( 1.0_wp / ( r_mid * LOG(s2) *      &
1044                                                             SQRT( 2.0_wp * pi )    &
1045                                                           ) *                      &
1046                                                  EXP( -( LOG( r_mid / rm2 ) )**2 / &
1047                                                       ( 2.0_wp * LOG(s2)**2 )      &
1048                                                     )                              &
1049                                                )
1050                   ENDIF
1051                   IF ( n3 > 0.0_wp )  THEN
1052                      pdf = pdf + n3 / n_init * ( 1.0_wp / ( r_mid * LOG(s3) *      &
1053                                                             SQRT( 2.0_wp * pi )    &
1054                                                           ) *                      &
1055                                                  EXP( -( LOG( r_mid / rm3 ) )**2 / &
1056                                                       ( 2.0_wp * LOG(s3)**2 )      &
1057                                                     )                              &
1058                                                )
1059                   ENDIF
1060
1061                   particles(n)%rvar2         = r_mid
1062                   particles(n)%weight_factor = pdf * dr
1063
1064                END DO
1065!
1066!--             Adjust weighting factors to initialize the same number of aerosols
1067!--             in every grid box
1068                weight_sum = SUM(particles(local_start(kp,jp,ip):number_of_particles)%weight_factor)
1069
1070                particles(local_start(kp,jp,ip):number_of_particles)%weight_factor =     &
1071                   particles(local_start(kp,jp,ip):number_of_particles)%weight_factor /  &
1072                   weight_sum * initial_weighting_factor * ( no_bins + 1 )
1073
1074             ENDIF
1075!
1076!--          Initialize the aerosols with a predefined weighting factor but
1077!--          a randomly choosen dry radius
1078             IF ( init_aerosol_probabilistic )  THEN
1079
1080                DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
1081
1082                   rs_rand = -1.0_wp
1083                   DO WHILE ( rs_rand .LT. cdf(0)  .OR.  rs_rand .GE. cdf(no_bins)  )
1084                      rs_rand = random_function( iran_part )
1085                   ENDDO
1086!
1087!--                Determine aerosol dry radius by a random number generator
1088                   DO nn = 0, no_bins-1
1089                      IF ( cdf(nn) .LE. rs_rand  .AND.  cdf(nn+1) .GT. rs_rand )  THEN
1090                         particles(n)%rvar2 = r_temp(nn) + ( r_temp(nn+1) - r_temp(nn) ) / &
1091                                              ( cdf(nn+1) - cdf(nn) ) * ( rs_rand - cdf(nn) )
1092                         EXIT
1093                      ENDIF
1094                   ENDDO
1095
1096                ENDDO
1097
1098             ENDIF
1099
1100!
1101!--          Set particle radius to equilibrium radius based on the environmental
1102!--          supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids
1103!--          the sometimes lengthy growth toward their equilibrium radius within
1104!--          the simulation.
1105             t_int  = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
1106
1107             e_s = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) )
1108             e_a = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp )
1109
1110             sigma   = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
1111             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
1112
1113             bfactor = vanthoff * molecular_weight_of_water *    &
1114                       rho_s / ( molecular_weight_of_solute * rho_l )
1115!
1116!--          The formula is only valid for subsaturated environments. For
1117!--          supersaturations higher than -5 %, the supersaturation is set to -5%.
1118             IF ( e_a / e_s >= 0.95_wp )  e_a = 0.95_wp * e_s
1119
1120             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
1121!
1122!--             For details on this equation, see Eq. (14) of Khvorostyanov and
1123!--             Curry (2007, JGR)
1124                particles(n)%radius = bfactor**0.3333333_wp *                  &
1125                   particles(n)%rvar2 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
1126                   ( 1.0_wp + ( afactor / ( 3.0_wp * bfactor**0.3333333_wp *   &
1127                     particles(n)%rvar2 ) ) /                                  &
1128                     ( 1.0_wp - e_a / e_s )**0.6666666_wp                      &
1129                   )
1130
1131             ENDDO
1132
1133          ENDDO
1134       ENDDO
1135    ENDDO
1136!
1137!-- Deallocate used arrays
1138    IF ( ALLOCATED(r_temp) )  DEALLOCATE( r_temp )
1139    IF ( ALLOCATED(cdf) )     DEALLOCATE( cdf )
1140
1141 END SUBROUTINE lpm_init_aerosols
1142
1143END MODULE lpm_init_mod
Note: See TracBrowser for help on using the repository browser.