Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/SOURCE/init_particles.f90 @ 513

Last change on this file since 513 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list
  • Property svn:keywords set to Id
File size: 21.8 KB
Line 
1 SUBROUTINE init_particles
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: init_particles.f90 392 2009-09-24 10:39:14Z heinze $
11!
12! 336 2009-06-10 11:19:35Z raasch
13! Maximum number of tails is calculated from maximum number of particles and
14! skip_particles_for_tail,
15! output of messages replaced by message handling routine
16! Bugfix: arrays for tails are allocated with a minimum size of 10 tails if
17! there is no tail initially
18!
19! 150 2008-02-29 08:19:58Z raasch
20! Setting offset_ocean_* needed for calculating vertical indices within ocean
21! runs
22!
23! 117 2007-10-11 03:27:59Z raasch
24! Sorting of particles only in case of cloud droplets
25!
26! 106 2007-08-16 14:30:26Z raasch
27! variable iran replaced by iran_part
28!
29! 82 2007-04-16 15:40:52Z raasch
30! Preprocessor directives for old systems removed
31!
32! 70 2007-03-18 23:46:30Z raasch
33! displacements for mpi_particle_type changed, age_m initialized,
34! particles-package is now part of the default code
35!
36! 16 2007-02-15 13:16:47Z raasch
37! Bugfix: MPI_REAL in MPI_ALLREDUCE replaced by MPI_INTEGER
38!
39! r4 | raasch | 2007-02-13 12:33:16 +0100 (Tue, 13 Feb 2007)
40! RCS Log replace by Id keyword, revision history cleaned up
41!
42! Revision 1.24  2007/02/11 13:00:17  raasch
43! Bugfix: allocation of tail_mask and new_tail_id in case of restart-runs
44! Bugfix: __ was missing in a cpp-directive
45!
46! Revision 1.1  1999/11/25 16:22:38  raasch
47! Initial revision
48!
49!
50! Description:
51! ------------
52! This routine initializes a set of particles and their attributes (position,
53! radius, ..). Advection of these particles is carried out by advec_particles,
54! plotting is done in data_output_dvrp.
55!------------------------------------------------------------------------------!
56
57    USE arrays_3d
58    USE control_parameters
59    USE dvrp_variables
60    USE grid_variables
61    USE indices
62    USE particle_attributes
63    USE pegrid
64    USE random_function_mod
65
66
67    IMPLICIT NONE
68
69    CHARACTER (LEN=10) ::  particle_binary_version, version_on_file
70
71    INTEGER ::  i, j, n, nn
72#if defined( __parallel )
73    INTEGER, DIMENSION(3) ::  blocklengths, displacements, types
74#endif
75    LOGICAL ::  uniform_particles_l
76    REAL    ::  factor, pos_x, pos_y, pos_z, value
77
78
79#if defined( __parallel )
80!
81!-- Define MPI derived datatype for FORTRAN datatype particle_type (see module
82!-- particle_attributes). Integer length is 4 byte, Real is 8 byte
83    blocklengths(1)  = 19;  blocklengths(2)  =   4;  blocklengths(3)  =   1
84    displacements(1) =  0;  displacements(2) = 152;  displacements(3) = 168
85
86    types(1) = MPI_REAL
87    types(2) = MPI_INTEGER
88    types(3) = MPI_UB
89    CALL MPI_TYPE_STRUCT( 3, blocklengths, displacements, types, &
90                          mpi_particle_type, ierr )
91    CALL MPI_TYPE_COMMIT( mpi_particle_type, ierr )
92#endif
93
94!
95!-- In case of oceans runs, the vertical index calculations need an offset,
96!-- because otherwise the k indices will become negative
97    IF ( ocean )  THEN
98       offset_ocean_nzt    = nzt
99       offset_ocean_nzt_m1 = nzt - 1
100    ENDIF
101
102
103!
104!-- Check the number of particle groups.
105    IF ( number_of_particle_groups > max_number_of_particle_groups )  THEN
106       WRITE( message_string, * ) 'max_number_of_particle_groups =',      &
107                                  max_number_of_particle_groups ,         &
108                                  '&number_of_particle_groups reset to ', &
109                                  max_number_of_particle_groups
110       CALL message( 'init_particles', 'PA0213', 0, 1, 0, 6, 0 )
111       number_of_particle_groups = max_number_of_particle_groups
112    ENDIF
113
114!
115!-- Set default start positions, if necessary
116    IF ( psl(1) == 9999999.9 )  psl(1) = -0.5 * dx
117    IF ( psr(1) == 9999999.9 )  psr(1) = ( nx + 0.5 ) * dx
118    IF ( pss(1) == 9999999.9 )  pss(1) = -0.5 * dy
119    IF ( psn(1) == 9999999.9 )  psn(1) = ( ny + 0.5 ) * dy
120    IF ( psb(1) == 9999999.9 )  psb(1) = zu(nz/2)
121    IF ( pst(1) == 9999999.9 )  pst(1) = psb(1)
122
123    IF ( pdx(1) == 9999999.9  .OR.  pdx(1) == 0.0 )  pdx(1) = dx
124    IF ( pdy(1) == 9999999.9  .OR.  pdy(1) == 0.0 )  pdy(1) = dy
125    IF ( pdz(1) == 9999999.9  .OR.  pdz(1) == 0.0 )  pdz(1) = zu(2) - zu(1)
126
127    DO  j = 2, number_of_particle_groups
128       IF ( psl(j) == 9999999.9 )  psl(j) = psl(j-1)
129       IF ( psr(j) == 9999999.9 )  psr(j) = psr(j-1)
130       IF ( pss(j) == 9999999.9 )  pss(j) = pss(j-1)
131       IF ( psn(j) == 9999999.9 )  psn(j) = psn(j-1)
132       IF ( psb(j) == 9999999.9 )  psb(j) = psb(j-1)
133       IF ( pst(j) == 9999999.9 )  pst(j) = pst(j-1)
134       IF ( pdx(j) == 9999999.9  .OR.  pdx(j) == 0.0 )  pdx(j) = pdx(j-1)
135       IF ( pdy(j) == 9999999.9  .OR.  pdy(j) == 0.0 )  pdy(j) = pdy(j-1)
136       IF ( pdz(j) == 9999999.9  .OR.  pdz(j) == 0.0 )  pdz(j) = pdz(j-1)
137    ENDDO
138
139!
140!-- For the first model run of a possible job chain initialize the
141!-- particles, otherwise read the particle data from file.
142    IF ( TRIM( initializing_actions ) == 'read_restart_data'  &
143         .AND.  read_particles_from_restartfile )  THEN
144
145!
146!--    Read particle data from previous model run.
147!--    First open the input unit.
148       IF ( myid_char == '' )  THEN
149          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN'//myid_char, &
150                     FORM='UNFORMATTED' )
151       ELSE
152          OPEN ( 90, FILE='PARTICLE_RESTART_DATA_IN/'//myid_char, &
153                     FORM='UNFORMATTED' )
154       ENDIF
155
156!
157!--    First compare the version numbers
158       READ ( 90 )  version_on_file
159       particle_binary_version = '3.0'
160       IF ( TRIM( version_on_file ) /= TRIM( particle_binary_version ) )  THEN
161          message_string = 'version mismatch concerning data from prior ' // &
162                           'run &version on file    = "' //                  &
163                                         TRIM( version_on_file ) //          &
164                           '&version in program = "' //                      &
165                                         TRIM( particle_binary_version ) // '"'
166          CALL message( 'init_particles', 'PA0214', 1, 2, 0, 6, 0 )
167       ENDIF
168
169!
170!--    Read some particle parameters and the size of the particle arrays,
171!--    allocate them and read their contents.
172       READ ( 90 )  bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,                  &
173                    maximum_number_of_particles, maximum_number_of_tailpoints, &
174                    maximum_number_of_tails, number_of_initial_particles,      &
175                    number_of_particles, number_of_particle_groups,            &
176                    number_of_tails, particle_groups, time_prel,               &
177                    time_write_particle_data, uniform_particles
178
179       IF ( number_of_initial_particles /= 0 )  THEN
180          ALLOCATE( initial_particles(1:number_of_initial_particles) )
181          READ ( 90 )  initial_particles
182       ENDIF
183
184       ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
185                 prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
186                 particle_mask(maximum_number_of_particles),         &
187                 particles(maximum_number_of_particles) )
188
189       READ ( 90 )  prt_count, prt_start_index
190       READ ( 90 )  particles
191
192       IF ( use_particle_tails )  THEN
193          ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
194                    maximum_number_of_tails),                                 &
195                    new_tail_id(maximum_number_of_tails),                     &
196                    tail_mask(maximum_number_of_tails) )
197          READ ( 90 )  particle_tail_coordinates
198       ENDIF
199
200       CLOSE ( 90 )
201
202    ELSE
203
204!
205!--    Allocate particle arrays and set attributes of the initial set of
206!--    particles, which can be also periodically released at later times.
207!--    Also allocate array for particle tail coordinates, if needed.
208       ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
209                 prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
210                 particle_mask(maximum_number_of_particles),         &
211                 particles(maximum_number_of_particles) )
212
213!
214!--    Initialize all particles with dummy values (otherwise errors may
215!--    occur within restart runs). The reason for this is still not clear
216!--    and may be presumably caused by errors in the respective user-interface.
217       particles = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
218                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
219                                  0.0, 0, 0, 0, 0 )
220       particle_groups = particle_groups_type( 0.0, 0.0, 0.0, 0.0 )
221
222!
223!--    Set the default particle size used for dvrp plots
224       IF ( dvrp_psize == 9999999.9 )  dvrp_psize = 0.2 * dx
225
226!
227!--    Set values for the density ratio and radius for all particle
228!--    groups, if necessary
229       IF ( density_ratio(1) == 9999999.9 )  density_ratio(1) = 0.0
230       IF ( radius(1)        == 9999999.9 )  radius(1) = 0.0
231       DO  i = 2, number_of_particle_groups
232          IF ( density_ratio(i) == 9999999.9 )  THEN
233             density_ratio(i) = density_ratio(i-1)
234          ENDIF
235          IF ( radius(i) == 9999999.9 )  radius(i) = radius(i-1)
236       ENDDO
237
238       DO  i = 1, number_of_particle_groups
239          IF ( density_ratio(i) /= 0.0  .AND.  radius(i) == 0 )  THEN
240             WRITE( message_string, * ) 'particle group #', i, 'has a', &
241                                        'density ratio /= 0 but radius = 0'
242             CALL message( 'init_particles', 'PA0215', 1, 2, 0, 6, 0 )
243          ENDIF
244          particle_groups(i)%density_ratio = density_ratio(i)
245          particle_groups(i)%radius        = radius(i)
246       ENDDO
247
248!
249!--    Calculate particle positions and store particle attributes, if
250!--    particle is situated on this PE
251       n = 0
252
253       DO  i = 1, number_of_particle_groups
254
255          pos_z = psb(i)
256
257          DO WHILE ( pos_z <= pst(i) )
258
259             pos_y = pss(i)
260
261             DO WHILE ( pos_y <= psn(i) )
262
263                IF ( pos_y >= ( nys - 0.5 ) * dy  .AND.  &
264                     pos_y <  ( nyn + 0.5 ) * dy )  THEN
265
266                   pos_x = psl(i)
267
268                   DO WHILE ( pos_x <= psr(i) )
269
270                      IF ( pos_x >= ( nxl - 0.5 ) * dx  .AND.  &
271                           pos_x <  ( nxr + 0.5 ) * dx )  THEN
272
273                         DO  j = 1, particles_per_point
274
275                            n = n + 1
276                            IF ( n > maximum_number_of_particles )  THEN
277                               WRITE( message_string, * ) 'number of initial', &
278                                      'particles (', n, ') exceeds',           &
279                                      '&maximum_number_of_particles (',        &
280                                      maximum_number_of_particles, ') on PE ', &
281                                             myid
282                               CALL message( 'init_particles', 'PA0216', &
283                                                                 2, 2, -1, 6, 1 )
284                            ENDIF
285                            particles(n)%x             = pos_x
286                            particles(n)%y             = pos_y
287                            particles(n)%z             = pos_z
288                            particles(n)%age           = 0.0
289                            particles(n)%age_m         = 0.0
290                            particles(n)%dt_sum        = 0.0
291                            particles(n)%dvrp_psize    = dvrp_psize
292                            particles(n)%e_m           = 0.0
293                            particles(n)%speed_x       = 0.0
294                            particles(n)%speed_x_sgs   = 0.0
295                            particles(n)%speed_y       = 0.0
296                            particles(n)%speed_y_sgs   = 0.0
297                            particles(n)%speed_z       = 0.0
298                            particles(n)%speed_z_sgs   = 0.0
299                            particles(n)%origin_x      = pos_x
300                            particles(n)%origin_y      = pos_y
301                            particles(n)%origin_z      = pos_z
302                            particles(n)%radius      = particle_groups(i)%radius
303                            particles(n)%weight_factor =initial_weighting_factor
304                            particles(n)%color         = 1
305                            particles(n)%group         = i
306                            particles(n)%tailpoints    = 0
307                            IF ( use_particle_tails  .AND. &
308                                 MOD( n, skip_particles_for_tail ) == 0 )  THEN
309                               number_of_tails         = number_of_tails + 1
310!
311!--                            This is a temporary provisional setting (see
312!--                            further below!)
313                               particles(n)%tail_id    = number_of_tails
314                            ELSE
315                               particles(n)%tail_id    = 0
316                            ENDIF
317
318                         ENDDO
319
320                      ENDIF
321
322                      pos_x = pos_x + pdx(i)
323
324                   ENDDO
325
326                ENDIF
327
328                pos_y = pos_y + pdy(i)
329
330             ENDDO
331
332             pos_z = pos_z + pdz(i)
333
334          ENDDO
335
336       ENDDO
337
338       number_of_initial_particles = n
339       number_of_particles         = n
340
341!
342!--    Calculate the number of particles and tails of the total domain
343#if defined( __parallel )
344       CALL MPI_ALLREDUCE( number_of_particles, total_number_of_particles, 1, &
345                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
346       CALL MPI_ALLREDUCE( number_of_tails, total_number_of_tails, 1, &
347                           MPI_INTEGER, MPI_SUM, comm2d, ierr )
348#else
349       total_number_of_particles = number_of_particles
350       total_number_of_tails     = number_of_tails
351#endif
352
353!
354!--    Set a seed value for the random number generator to be exclusively
355!--    used for the particle code. The generated random numbers should be
356!--    different on the different PEs.
357       iran_part = iran_part + myid
358
359!
360!--    User modification of initial particles
361       CALL user_init_particles
362
363!
364!--    Store the initial set of particles for release at later times
365       IF ( number_of_initial_particles /= 0 )  THEN
366          ALLOCATE( initial_particles(1:number_of_initial_particles) )
367          initial_particles(1:number_of_initial_particles) = &
368                                        particles(1:number_of_initial_particles)
369       ENDIF
370
371!
372!--    Add random fluctuation to particle positions
373       IF ( random_start_position )  THEN
374
375          DO  n = 1, number_of_initial_particles
376             IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
377                particles(n)%x = particles(n)%x + &
378                                 ( random_function( iran_part ) - 0.5 ) * &
379                                 pdx(particles(n)%group)
380                IF ( particles(n)%x  <=  ( nxl - 0.5 ) * dx )  THEN
381                   particles(n)%x = ( nxl - 0.4999999999 ) * dx
382                ELSEIF ( particles(n)%x  >=  ( nxr + 0.5 ) * dx )  THEN
383                   particles(n)%x = ( nxr + 0.4999999999 ) * dx
384                ENDIF
385             ENDIF
386             IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
387                particles(n)%y = particles(n)%y + &
388                                 ( random_function( iran_part ) - 0.5 ) * &
389                                 pdy(particles(n)%group)
390                IF ( particles(n)%y  <=  ( nys - 0.5 ) * dy )  THEN
391                   particles(n)%y = ( nys - 0.4999999999 ) * dy
392                ELSEIF ( particles(n)%y  >=  ( nyn + 0.5 ) * dy )  THEN
393                   particles(n)%y = ( nyn + 0.4999999999 ) * dy
394                ENDIF
395             ENDIF
396             IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
397                particles(n)%z = particles(n)%z + &
398                                 ( random_function( iran_part ) - 0.5 ) * &
399                                 pdz(particles(n)%group)
400             ENDIF
401          ENDDO
402       ENDIF
403
404!
405!--    Sort particles in the sequence the gridboxes are stored in the memory.
406!--    Only required if cloud droplets are used.
407       IF ( cloud_droplets )  CALL sort_particles
408
409!
410!--    Open file for statistical informations about particle conditions
411       IF ( write_particle_statistics )  THEN
412          CALL check_open( 80 )
413          WRITE ( 80, 8000 )  current_timestep_number, simulated_time, &
414                              number_of_initial_particles,             &
415                              maximum_number_of_particles
416          CALL close_file( 80 )
417       ENDIF
418
419!
420!--    Check if particles are really uniform in color and radius (dvrp_size)
421!--    (uniform_particles is preset TRUE)
422       IF ( uniform_particles )  THEN
423          IF ( number_of_initial_particles == 0 )  THEN
424             uniform_particles_l = .TRUE.
425          ELSE
426             n = number_of_initial_particles
427             IF ( MINVAL( particles(1:n)%dvrp_psize  ) ==     &
428                  MAXVAL( particles(1:n)%dvrp_psize  )  .AND. &
429                  MINVAL( particles(1:n)%color ) ==     &
430                  MAXVAL( particles(1:n)%color ) )  THEN
431                uniform_particles_l = .TRUE.
432             ELSE
433                uniform_particles_l = .FALSE.
434             ENDIF
435          ENDIF
436
437#if defined( __parallel )
438          CALL MPI_ALLREDUCE( uniform_particles_l, uniform_particles, 1, &
439                              MPI_LOGICAL, MPI_LAND, comm2d, ierr )
440#else
441          uniform_particles = uniform_particles_l
442#endif
443
444       ENDIF
445
446!
447!--    Particles will probably become none-uniform, if their size and color
448!--    will be determined by flow variables
449       IF ( particle_color /= 'none'  .OR.  particle_dvrpsize /= 'none' )  THEN
450          uniform_particles = .FALSE.
451       ENDIF
452
453!
454!--    Set the beginning of the particle tails and their age
455       IF ( use_particle_tails )  THEN
456!
457!--       Choose the maximum number of tails with respect to the maximum number
458!--       of particles and skip_particles_for_tail
459          maximum_number_of_tails = maximum_number_of_particles / &
460                                    skip_particles_for_tail
461
462!
463!--       Create a minimum number of tails in case that there is no tail
464!--       initially (otherwise, index errors will occur when adressing the
465!--       arrays below)
466          IF ( maximum_number_of_tails == 0 )  maximum_number_of_tails = 10
467
468          ALLOCATE( particle_tail_coordinates(maximum_number_of_tailpoints,5, &
469                    maximum_number_of_tails),                                 &
470                    new_tail_id(maximum_number_of_tails),                     &
471                    tail_mask(maximum_number_of_tails) )
472
473          particle_tail_coordinates  = 0.0
474          minimum_tailpoint_distance = minimum_tailpoint_distance**2
475          number_of_initial_tails    = number_of_tails
476
477          nn = 0
478          DO  n = 1, number_of_particles
479!
480!--          Only for those particles marked above with a provisional tail_id
481!--          tails will be created. Particles now get their final tail_id.
482             IF ( particles(n)%tail_id /= 0 )  THEN
483
484                nn = nn + 1
485                particles(n)%tail_id = nn
486
487                particle_tail_coordinates(1,1,nn) = particles(n)%x
488                particle_tail_coordinates(1,2,nn) = particles(n)%y
489                particle_tail_coordinates(1,3,nn) = particles(n)%z
490                particle_tail_coordinates(1,4,nn) = particles(n)%color
491                particles(n)%tailpoints = 1
492                IF ( minimum_tailpoint_distance /= 0.0 )  THEN
493                   particle_tail_coordinates(2,1,nn) = particles(n)%x
494                   particle_tail_coordinates(2,2,nn) = particles(n)%y
495                   particle_tail_coordinates(2,3,nn) = particles(n)%z
496                   particle_tail_coordinates(2,4,nn) = particles(n)%color
497                   particle_tail_coordinates(1:2,5,nn) = 0.0
498                   particles(n)%tailpoints = 2
499                ENDIF
500
501             ENDIF
502          ENDDO
503       ENDIF
504
505!
506!--    Plot initial positions of particles (only if particle advection is
507!--    switched on from the beginning of the simulation (t=0))
508       IF ( particle_advection_start == 0.0 )  CALL data_output_dvrp
509
510    ENDIF
511
512!
513!-- Check boundary condition and set internal variables
514    SELECT CASE ( bc_par_b )
515   
516       CASE ( 'absorb' )
517          ibc_par_b = 1
518
519       CASE ( 'reflect' )
520          ibc_par_b = 2
521         
522       CASE DEFAULT
523          WRITE( message_string, * )  'unknown boundary condition ',   &
524                                       'bc_par_b = "', TRIM( bc_par_b ), '"'
525          CALL message( 'init_particles', 'PA0217', 1, 2, 0, 6, 0 )
526         
527    END SELECT
528    SELECT CASE ( bc_par_t )
529   
530       CASE ( 'absorb' )
531          ibc_par_t = 1
532
533       CASE ( 'reflect' )
534          ibc_par_t = 2
535         
536       CASE DEFAULT
537          WRITE( message_string, * ) 'unknown boundary condition ',   &
538                                     'bc_par_t = "', TRIM( bc_par_t ), '"'
539          CALL message( 'init_particles', 'PA0218', 1, 2, 0, 6, 0 )
540         
541    END SELECT
542    SELECT CASE ( bc_par_lr )
543
544       CASE ( 'cyclic' )
545          ibc_par_lr = 0
546
547       CASE ( 'absorb' )
548          ibc_par_lr = 1
549
550       CASE ( 'reflect' )
551          ibc_par_lr = 2
552         
553       CASE DEFAULT
554          WRITE( message_string, * ) 'unknown boundary condition ',   &
555                                     'bc_par_lr = "', TRIM( bc_par_lr ), '"'
556          CALL message( 'init_particles', 'PA0219', 1, 2, 0, 6, 0 )
557         
558    END SELECT
559    SELECT CASE ( bc_par_ns )
560
561       CASE ( 'cyclic' )
562          ibc_par_ns = 0
563
564       CASE ( 'absorb' )
565          ibc_par_ns = 1
566
567       CASE ( 'reflect' )
568          ibc_par_ns = 2
569         
570       CASE DEFAULT
571          WRITE( message_string, * ) 'unknown boundary condition ',   &
572                                     'bc_par_ns = "', TRIM( bc_par_ns ), '"'
573          CALL message( 'init_particles', 'PA0220', 1, 2, 0, 6, 0 )
574         
575    END SELECT
576!
577!-- Formats
5788000 FORMAT (I6,1X,F7.2,4X,I6,71X,I6)
579
580 END SUBROUTINE init_particles
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |