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

Last change on this file since 77 was 77, checked in by raasch, 14 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

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