Changeset 825


Ignore:
Timestamp:
Feb 19, 2012 3:03:44 AM (10 years ago)
Author:
raasch
Message:

New:
---

droplet growth by condensation may include curvature and solution effects.
Steered by new inipar-parameter curvature_solution_effects.
(advec_particles, check_parameters, header, init_cloud_physics, init_particles, modules, parin, read_var_list, write_var_list)

mean/min/max particle radius added as output quantity. (data_output_ptseries, modules)

Changed:


Initialisation of temporary particle array for resorting removed.
(advec_particles)

particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
(advec_particles, data_output_ptseries, modules, init_particles, particle_boundary_conds)

routine wang_kernel and respective module renamed lpm_collision_kernels.
Package (particle) parameters wang_collision_kernel and turbulence_effects_on_collision
replaced by parameter collision_kernel.
(Makefile, advec_particles, check_parameters, diffusion_e, init_3d_model, modules, package_parin, time_integration, new: lpm_collision_kernels, deleted: wang_kernel)

Errors:


Location:
palm/trunk/SOURCE
Files:
10 edited
1 moved

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r791 r825  
    44# Current revisions:
    55# ------------------
    6 #
     6# wang_kernel renamed lpm_collision_kernels
    77#
    88# Former revisions:
     
    8383        init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \
    8484        init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
    85         interaction_droplets_ptq.f90 local_flush.f90 \
    86         local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
    87         local_tremain_ini.f90 message.f90 modules.f90 netcdf.f90 \
     85        interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
     86        local_stop.f90 local_system.f90 local_tremain.f90 local_tremain_ini.f90 \
     87        lpm_collision_kernels.f90 message.f90 modules.f90 netcdf.f90 \
    8888        package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \
    8989        plant_canopy_model.f90 poisfft.f90 \
     
    108108        user_particle_attributes.f90 user_read_restart_data.f90 \
    109109        user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
    110         wang_kernel.f90 write_3d_binary.f90 write_compressed.f90 \
     110        write_3d_binary.f90 write_compressed.f90 \
    111111        write_var_list.f90
    112112
     
    124124        data_output_3d.o diffusion_e.o diffusion_s.o diffusion_u.o \
    125125        diffusion_v.o diffusion_w.o diffusivities.o disturb_field.o \
    126         disturb_heatflux.o eqn_state_seawater.o exchange_horiz.o exchange_horiz_2d.o fft_xy.o \
    127         flow_statistics.o global_min_max.o header.o impact_of_latent_heat.o \
    128         inflow_turbulence.o init_1d_model.o init_3d_model.o init_advec.o init_cloud_physics.o \
    129         init_coupling.o init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o init_pegrid.o \
    130         init_pt_anomaly.o init_rankine.o init_slope.o \
     126        disturb_heatflux.o eqn_state_seawater.o exchange_horiz.o \
     127        exchange_horiz_2d.o fft_xy.o flow_statistics.o global_min_max.o \
     128        header.o impact_of_latent_heat.o inflow_turbulence.o init_1d_model.o \
     129        init_3d_model.o init_advec.o init_cloud_physics.o init_coupling.o \
     130        init_dvrp.o init_grid.o init_masks.o init_ocean.o init_particles.o \
     131        init_pegrid.o init_pt_anomaly.o init_rankine.o init_slope.o \
    131132        interaction_droplets_ptq.o local_flush.o local_getenv.o local_stop.o \
    132         local_system.o local_tremain.o local_tremain_ini.o message.o modules.o \
     133        local_system.o local_tremain.o local_tremain_ini.o \
     134        lpm_collision_kernels.f90 message.o modules.o \
    133135        netcdf.o package_parin.o palm.o parin.o particle_boundary_conds.o \
    134136        plant_canopy_model.o poisfft.o \
     
    151153        user_module.o user_parin.o user_particle_attributes.o \
    152154        user_read_restart_data.o user_spectra.o user_statistics.o \
    153         wall_fluxes.o wang_kernel.o write_3d_binary.o write_compressed.o \
     155        wall_fluxes.o write_3d_binary.o write_compressed.o \
    154156        write_var_list.o
    155157
     
    181183
    182184
    183 advec_particles.o: modules.o random_function.o wang_kernel.o
     185advec_particles.o: modules.o lpm_collision_kernels.o random_function.o
    184186advec_s_bc.o: modules.o
    185187advec_s_pw.o: modules.o
     
    257259local_tremain.o: modules.o
    258260local_tremain_ini.o: modules.o
     261lpm_collision_kernels.o: modules.o user_module.o
    259262message.o: modules.o
    260263modules.o: modules.f90
     
    328331user_statistics.o: modules.o user_module.o
    329332wall_fluxes.o: modules.o
    330 wang_kernel.o: modules.o user_module.o
    331333write_3d_binary.o: modules.o random_function.o
    332334write_compressed.o: modules.o
  • palm/trunk/SOURCE/advec_particles.f90

    r824 r825  
    44! Current revisions:
    55! ------------------
    6 ! droplet growth by condensation may include curvature and solution effects.
    7 ! initialisation of temporary particle array for resorting removed.
    8 ! particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
     6! droplet growth by condensation may include curvature and solution effects,
     7! initialisation of temporary particle array for resorting removed,
     8! particle attributes speed_x|y|z_sgs renamed rvar1|2|3,
     9! module wang_kernel_mod renamed lpm_collision_kernels_mod,
     10! wang_collision_kernel renamed wang_kernel
    911!
    1012! Former revisions:
     
    130132    USE indices
    131133    USE interfaces
     134    USE lpm_collision_kernels_mod
    132135    USE netcdf_control
    133136    USE particle_attributes
     
    135138    USE random_function_mod
    136139    USE statistics
    137     USE wang_kernel_mod
    138140
    139141    IMPLICIT NONE
     
    200202    CALL cpu_log( log_point(25), 'advec_particles', 'start' )
    201203
    202 !    IF ( number_of_particles /= number_of_tails )  THEN
    203 !       WRITE (9,*) '--- advec_particles: #1'
    204 !       WRITE (9,*) '    #of p=',number_of_particles,' #of t=',number_of_tails
    205 !       CALL local_flush( 9 )
    206 !    ENDIF
     204
    207205!
    208206!-- Write particle data on file for later analysis.
     
    270268    CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' )
    271269
    272 !    WRITE ( 9, * ) '*** advec_particles: ##0.3'
    273 !    CALL local_flush( 9 )
    274 !    nd = 0
    275 !    DO  n = 1, number_of_particles
    276 !       IF ( .NOT. particle_mask(n) )  nd = nd + 1
    277 !    ENDDO
    278 !    IF ( nd /= deleted_particles ) THEN
    279 !       WRITE (9,*) '*** nd=',nd,' deleted_particles=',deleted_particles
    280 !       CALL local_flush( 9 )
    281 !       CALL MPI_ABORT( comm2d, 9999, ierr )
    282 !    ENDIF
    283 
    284270!
    285271!-- Particle (droplet) growth by condensation/evaporation and collision
     
    357343!
    358344!--       Change in radius by condensation/evaporation
    359           IF ( curvature_solution_effects .AND. particles(n)%radius < 1.0E-6 ) &
    360           THEN
    361 !
    362 !--          Curvature and solutions effects are included in growth equation.
    363 !--          Change in Radius is calculated with 4th-order Rosenbrock method
    364 !--          for stiff o.d.e's with monitoring local truncation error to adjust
    365 !--          stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
    366 !--          For larger radii the simple analytic method (see ELSE) gives
    367 !--          almost the same results.
    368 !
    369 !--          Surface tension after (Straka, 2009)
    370              sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
    371 
    372              r_ros = particles(n)%radius
    373              dt_ros_sum  = 0.0      ! internal integrated time (s)
    374              internal_timestep_count = 0
    375 
    376              ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
    377                                ( l_v / ( r_v * t_int ) - 1.0 ) *               &
    378                                rho_l * l_v / ( thermal_conductivity_l * t_int )&
    379                              )
    380              
    381              afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
    382 
    383              IF ( particles(n)%rvar3 == -9999999.9 )  THEN
    384 !
    385 !--             First particle timestep. Derivative has to be calculated.
    386                 drdt_m  = ddenom / r_ros * ( e_a / e_s - 1.0 -    &
    387                                              afactor / r_ros +    &
    388                                              bfactor / r_ros**3 )
    389              ELSE
    390 !
    391 !--             Take value from last PALM timestep
    392                 drdt_m = particles(n)%rvar3
    393              ENDIF
    394 !
    395 !--          Take internal timestep values from the end of last PALM timestep
    396              dt_ros_last = particles(n)%rvar1
    397              dt_ros_next = particles(n)%rvar2
    398 !
    399 !--          Internal timestep must not be larger than PALM timestep
    400              dt_ros = MIN( dt_ros_next, dt_3d )
    401 !
    402 !--          Integrate growth equation in time unless PALM timestep is reached
    403              DO WHILE ( dt_ros_sum < dt_3d )
    404 
    405                 internal_timestep_count = internal_timestep_count + 1
    406 
    407 !
    408 !--             Derivative at starting value
    409                 drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
    410                                           bfactor / r_ros**3 )
    411                 drdt_ini       = drdt
    412                 dt_ros_sum_ini = dt_ros_sum
    413                 r_ros_ini      = r_ros
    414 
    415 !
    416 !--             Calculate time derivative of dr/dt
    417                 d2rdt2 = ( drdt - drdt_m ) / dt_ros_last
    418 
    419 !
    420 !--             Calculate radial derivative of dr/dt
    421                 d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
    422                                      2.0 * afactor / r_ros**3 -     &
    423                                      4.0 * bfactor / r_ros**5 )
    424 !
    425 !--             Adjust stepsize unless required accuracy is reached
    426                 DO  jtry = 1, maxtry+1
    427 
    428                    IF ( jtry == maxtry+1 )  THEN
    429                       message_string = 'maxtry > 40 in Rosenbrock method'
    430                       CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 )
    431                    ENDIF
    432 
    433                    aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
    434                    g1    = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
    435                    r_ros = r_ros_ini + a21 * g1
    436                    drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
    437                                               afactor / r_ros + &
    438                                               bfactor / r_ros**3 )
    439 
    440                    g2    = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
    441                            / aa
    442                    r_ros = r_ros_ini + a31 * g1 + a32 * g2
    443                    drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
    444                                               afactor / r_ros + &
    445                                               bfactor / r_ros**3 )
    446 
    447                    g3    = ( drdt + dt_ros * c3x * d2rdt2 + &
    448                              ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
    449                    g4    = ( drdt + dt_ros * c4x * d2rdt2 + &
    450                              ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
    451                    r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
    452 
    453                    dt_ros_sum = dt_ros_sum_ini + dt_ros
    454 
    455                    IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
    456                       message_string = 'zero stepsize in Rosenbrock method'
    457                       CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 )
    458                    ENDIF
    459 !
    460 !--                Calculate error
    461                    err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4
    462                    errmax  = 0.0
    463                    errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
    464 !
    465 !--                Leave loop if accuracy is sufficient, otherwise try again
    466 !--                with a reduced stepsize
    467                    IF ( errmax <= 1.0 )  THEN
    468                       EXIT
    469                    ELSE
    470                       dt_ros = SIGN(                                           &
    471                                     MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
    472                                          shrnk * ABS( dt_ros ) ),              &
    473                                     dt_ros                                     &
    474                                    )
    475                    ENDIF
    476 
    477                 ENDDO  ! loop for stepsize adjustment
    478 
    479 !
    480 !--             Calculate next internal timestep
    481                 IF ( errmax > errcon )  THEN
    482                    dt_ros_next = 0.9 * dt_ros * errmax**pgrow
    483                 ELSE
    484                    dt_ros_next = grow * dt_ros
    485                 ENDIF
    486 
    487 !
    488 !--             Estimated timestep is reduced if the PALM time step is exceeded
    489                 dt_ros_last = dt_ros
    490                 IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
    491                    dt_ros = dt_3d - dt_ros_sum
    492                 ELSE
    493                    dt_ros = dt_ros_next
    494                 ENDIF
    495 
    496                 drdt_m = drdt
    497 
    498              ENDDO
    499 !
    500 !--          Store derivative and internal timestep values for next PALM step
    501              particles(n)%rvar1 = dt_ros_last
    502              particles(n)%rvar2 = dt_ros_next
    503              particles(n)%rvar3 = drdt_m
    504 
    505              new_r = r_ros
    506 !
    507 !--          Radius should not fall below 1E-8 because Rosenbrock method may
    508 !--          lead to errors otherwise
    509              new_r = MAX( new_r, 1.0E-8 )
    510 
    511           ELSE
     345          IF ( particles(n)%radius >= 1.0E-6  .OR.  &
     346               .NOT. curvature_solution_effects )  THEN
    512347!
    513348!--          Approximation for large radii, where curvature and solution
     
    523358                new_r = SQRT( arg )
    524359             ENDIF
     360          ENDIF
     361
     362          IF ( curvature_solution_effects  .AND.                              &
     363               ( ( particles(n)%radius < 1.0E-6 ) .OR. ( new_r < 1.0E-6 ) ) ) &
     364          THEN
     365!
     366!--          Curvature and solutions effects are included in growth equation.
     367!--          Change in Radius is calculated with 4th-order Rosenbrock method
     368!--          for stiff o.d.e's with monitoring local truncation error to adjust
     369!--          stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
     370!--          For larger radii the simple analytic method (see ELSE) gives
     371!--          almost the same results.
     372!
     373!--          Surface tension after (Straka, 2009)
     374             sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
     375
     376             r_ros = particles(n)%radius
     377             dt_ros_sum  = 0.0      ! internal integrated time (s)
     378             internal_timestep_count = 0
     379
     380             ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
     381                               ( l_v / ( r_v * t_int ) - 1.0 ) *               &
     382                               rho_l * l_v / ( thermal_conductivity_l * t_int )&
     383                             )
     384             
     385             afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
     386
     387             IF ( particles(n)%rvar3 == -9999999.9 )  THEN
     388!
     389!--             First particle timestep. Derivative has to be calculated.
     390                drdt_m  = ddenom / r_ros * ( e_a / e_s - 1.0 -    &
     391                                             afactor / r_ros +    &
     392                                             bfactor / r_ros**3 )
     393             ELSE
     394!
     395!--             Take value from last PALM timestep
     396                drdt_m = particles(n)%rvar3
     397             ENDIF
     398!
     399!--          Take internal timestep values from the end of last PALM timestep
     400             dt_ros_last = particles(n)%rvar1
     401             dt_ros_next = particles(n)%rvar2
     402!
     403!--          Internal timestep must not be larger than PALM timestep
     404             dt_ros = MIN( dt_ros_next, dt_3d )
     405!
     406!--          Integrate growth equation in time unless PALM timestep is reached
     407             DO WHILE ( dt_ros_sum < dt_3d )
     408
     409                internal_timestep_count = internal_timestep_count + 1
     410
     411!
     412!--             Derivative at starting value
     413                drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
     414                                          bfactor / r_ros**3 )
     415                drdt_ini       = drdt
     416                dt_ros_sum_ini = dt_ros_sum
     417                r_ros_ini      = r_ros
     418
     419!
     420!--             Calculate time derivative of dr/dt
     421                d2rdt2 = ( drdt - drdt_m ) / dt_ros_last
     422
     423!
     424!--             Calculate radial derivative of dr/dt
     425                d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
     426                                     2.0 * afactor / r_ros**3 -     &
     427                                     4.0 * bfactor / r_ros**5 )
     428!
     429!--             Adjust stepsize unless required accuracy is reached
     430                DO  jtry = 1, maxtry+1
     431
     432                   IF ( jtry == maxtry+1 )  THEN
     433                      message_string = 'maxtry > 40 in Rosenbrock method'
     434                      CALL message( 'advec_particles', 'PA0347', 2, 2, 0, 6, 0 )
     435                   ENDIF
     436
     437                   aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
     438                   g1    = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
     439                   r_ros = r_ros_ini + a21 * g1
     440                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
     441                                              afactor / r_ros + &
     442                                              bfactor / r_ros**3 )
     443
     444                   g2    = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
     445                           / aa
     446                   r_ros = r_ros_ini + a31 * g1 + a32 * g2
     447                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
     448                                              afactor / r_ros + &
     449                                              bfactor / r_ros**3 )
     450
     451                   g3    = ( drdt + dt_ros * c3x * d2rdt2 + &
     452                             ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
     453                   g4    = ( drdt + dt_ros * c4x * d2rdt2 + &
     454                             ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
     455                   r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
     456
     457                   dt_ros_sum = dt_ros_sum_ini + dt_ros
     458
     459                   IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
     460                      message_string = 'zero stepsize in Rosenbrock method'
     461                      CALL message( 'advec_particles', 'PA0348', 2, 2, 0, 6, 0 )
     462                   ENDIF
     463!
     464!--                Calculate error
     465                   err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4
     466                   errmax  = 0.0
     467                   errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
     468!
     469!--                Leave loop if accuracy is sufficient, otherwise try again
     470!--                with a reduced stepsize
     471                   IF ( errmax <= 1.0 )  THEN
     472                      EXIT
     473                   ELSE
     474                      dt_ros = SIGN(                                           &
     475                                    MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
     476                                         shrnk * ABS( dt_ros ) ),              &
     477                                    dt_ros                                     &
     478                                   )
     479                   ENDIF
     480
     481                ENDDO  ! loop for stepsize adjustment
     482
     483!
     484!--             Calculate next internal timestep
     485                IF ( errmax > errcon )  THEN
     486                   dt_ros_next = 0.9 * dt_ros * errmax**pgrow
     487                ELSE
     488                   dt_ros_next = grow * dt_ros
     489                ENDIF
     490
     491!
     492!--             Estimated timestep is reduced if the PALM time step is exceeded
     493                dt_ros_last = dt_ros
     494                IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
     495                   dt_ros = dt_3d - dt_ros_sum
     496                ELSE
     497                   dt_ros = dt_ros_next
     498                ENDIF
     499
     500                drdt_m = drdt
     501
     502             ENDDO
     503!
     504!--          Store derivative and internal timestep values for next PALM step
     505             particles(n)%rvar1 = dt_ros_last
     506             particles(n)%rvar2 = dt_ros_next
     507             particles(n)%rvar3 = drdt_m
     508
     509             new_r = r_ros
     510!
     511!--          Radius should not fall below 1E-8 because Rosenbrock method may
     512!--          lead to errors otherwise
     513             new_r = MAX( new_r, 1.0E-8 )
    525514
    526515          ENDIF
     
    570559!
    571560!--    Particle growth by collision
     561       IF ( collision_kernel /= 'none' )  THEN
     562
    572563       CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
    573 
    574        IF ( wang_collision_kernel ) THEN
    575564
    576565       DO  i = nxl, nxr
     
    609598                   pse = psi + prt_count(k,j,i)-1
    610599
    611                    ALLOCATE( kern(psi:pse,psi:pse) )
    612 
    613 !
    614 !--                Calculate collision kernel for all particles in grid box
    615                    CALL colker( i, j, k, kern )
    616 !
    617 !--                collison kernel is calculated in cm**3/s but needed in m**3/s
    618                    kern = kern * 1.0E-6
    619 
    620                    DO  n = pse, psi+1, -1
    621 
    622                       integral = 0.0
    623                       lw_max   = 0.0
    624 
    625 !
    626 !--                   Calculate growth of collector particle radius using all
    627 !--                   droplets smaller than current droplet
    628                       DO  is = psi, n-1
    629 
    630                          integral = integral +                               &
    631                                     ( particles(is)%radius**3 * kern(n,is) * &
    632                                       particles(is)%weight_factor )
    633 !
    634 !--                      Calculate volume of liquid water of the collected
    635 !--                      droplets which is the maximum liquid water available
    636 !--                      for droplet growth   
    637                          lw_max = lw_max + ( particles(is)%radius**3 * &
    638                                              particles(is)%weight_factor )
     600                   IF ( hall_kernel  .OR.  wang_kernel )  THEN
     601
     602                      ALLOCATE( kern(psi:pse,psi:pse) )
     603
     604!
     605!--                   Calculate collision kernel for all particles in grid box
     606                      CALL colker( i, j, k, kern )
     607!
     608!--                   Collison kernel is calculated in cm**3/s but needed
     609!--                   in m**3/s
     610                      kern = kern * 1.0E-6   ! to be moved to colker
     611
     612                      DO  n = pse, psi+1, -1
     613
     614                         integral = 0.0
     615                         lw_max   = 0.0
     616
     617!
     618!--                      Calculate growth of collector particle radius using all
     619!--                      droplets smaller than current droplet
     620                         DO  is = psi, n-1
     621
     622                            integral = integral +                              &
     623                                       ( particles(is)%radius**3 * kern(n,is) *&
     624                                         particles(is)%weight_factor )
     625!
     626!--                         Calculate volume of liquid water of the collected
     627!--                         droplets which is the maximum liquid water available
     628!--                         for droplet growth   
     629                            lw_max = lw_max + ( particles(is)%radius**3 * &
     630                                                particles(is)%weight_factor )
     631
     632                         ENDDO
     633
     634!
     635!--                      Change in radius of collector droplet due to collision
     636                         delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
     637                                   integral * dt_3d * ddx * ddy / dz
     638
     639!
     640!--                      Change in volume of collector droplet due to collision
     641                         delta_v = particles(n)%weight_factor                  &
     642                                   * ( ( particles(n)%radius + delta_r )**3    &
     643                                       - particles(n)%radius**3 )
     644
     645                         IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
     646!-- replace by message call
     647                            PRINT*, 'Particle has grown to much because the',  &
     648                                    ' change of volume of particle is larger', &
     649                                    ' than liquid water available!'
     650
     651                         ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
     652
     653                            DO  is = psi, n-1
     654
     655                               particles(is)%weight_factor = 0.0
     656                               particle_mask(is)  = .FALSE.
     657                               deleted_particles = deleted_particles + 1
     658
     659                            ENDDO
     660
     661                         ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
     662!
     663!--                         Calculate new weighting factor of collected droplets
     664                            DO  is = psi, n-1
     665
     666                               particles(is)%weight_factor = &
     667                                  particles(is)%weight_factor - &
     668                                  ( ( kern(n,is) * particles(is)%weight_factor &
     669                                    / integral ) * delta_v )
     670
     671                               IF ( particles(is)%weight_factor < 0.0 )  THEN
     672                                     WRITE( message_string, * ) 'negative ',   &
     673                                                        'weighting factor: ',  &
     674                                                     particles(is)%weight_factor
     675                                     CALL message( 'advec_particles', '', 2, 2,&
     676                                                   -1, 6, 1 )
     677
     678                               ELSEIF ( particles(is)%weight_factor == 0.0 )  &
     679                               THEN
     680
     681                                   particles(is)%weight_factor = 0.0
     682                                   particle_mask(is)  = .FALSE.
     683                                   deleted_particles = deleted_particles + 1
     684
     685                               ENDIF
     686
     687                            ENDDO
     688
     689                         ENDIF
     690
     691                         particles(n)%radius = particles(n)%radius + delta_r
     692                         ql_vp(k,j,i) = ql_vp(k,j,i) + &
     693                                                    particles(n)%weight_factor &
     694                                                    * particles(n)%radius**3
    639695
    640696                      ENDDO
    641697
    642 !
    643 !--                   Change in radius of collector droplet due to collision
    644                       delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
    645                                 integral * dt_3d * ddx * ddy / dz
    646 
    647 !
    648 !--                   Change in volume of collector droplet due to collision
    649                       delta_v = particles(n)%weight_factor                     &
    650                                 * ( ( particles(n)%radius + delta_r )**3       &
    651                                     - particles(n)%radius**3 )
    652 
    653                       IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
    654 
    655                          PRINT*, 'Particle has grown to much because the       &
    656                                   change of volume of particle is larger than  &
    657                                   liquid water available!'
    658 
    659                       ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 )  THEN
    660 
    661                          DO  is = psi, n-1
    662 
    663                             particles(is)%weight_factor = 0.0
    664                             particle_mask(is)  = .FALSE.
    665                             deleted_particles = deleted_particles + 1
    666 
     698                      DEALLOCATE( kern )
     699
     700                   ELSEIF ( palm_kernel )  THEN
     701!
     702!--                   PALM collision kernel
     703!
     704!--                   Calculate the mean radius of all those particles which
     705!--                   are of smaller size than the current particle and
     706!--                   use this radius for calculating the collision efficiency
     707                      DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
     708
     709                         sl_r3 = 0.0
     710                         sl_r4 = 0.0
     711
     712                         DO is = n-1, psi, -1
     713                            IF ( particles(is)%radius < particles(n)%radius )  &
     714                            THEN
     715                               sl_r3 = sl_r3 + particles(is)%weight_factor     &
     716                                             *( particles(is)%radius**3 )
     717                               sl_r4 = sl_r4 + particles(is)%weight_factor     &
     718                                             *( particles(is)%radius**4 )
     719                            ENDIF
    667720                         ENDDO
    668721
    669                       ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
    670 !
    671 !--                      Calculate new weighting factor of collected droplets
    672                          DO  is = psi, n-1
    673 
    674                             particles(is)%weight_factor = &
    675                                particles(is)%weight_factor - &
    676                                ( ( kern(n,is) * particles(is)%weight_factor / &
    677                                    integral ) * delta_v )
    678 
    679                             IF ( particles(is)%weight_factor < 0.0 )  THEN
    680                                   WRITE( message_string, * ) 'negative ',      &
     722                         IF ( ( sl_r3 ) > 0.0 )  THEN
     723                            mean_r = ( sl_r4 ) / ( sl_r3 )
     724
     725                            CALL collision_efficiency( mean_r,                 &
     726                                                    particles(n)%radius,       &
     727                                                    effective_coll_efficiency )
     728
     729                         ELSE
     730                            effective_coll_efficiency = 0.0
     731                         ENDIF
     732
     733                         IF ( effective_coll_efficiency > 1.0  .OR.            &
     734                              effective_coll_efficiency < 0.0 )                &
     735                         THEN
     736                            WRITE( message_string, * )  'collision_efficien' , &
     737                                   'cy out of range:' ,effective_coll_efficiency
     738                            CALL message( 'advec_particles', 'PA0145', 2, 2,   &
     739                                          -1, 6, 1 )
     740                         ENDIF
     741
     742!
     743!--                      Interpolation of ...
     744                         ii = particles(n)%x * ddx
     745                         jj = particles(n)%y * ddy
     746                         kk = ( particles(n)%z + 0.5 * dz ) / dz
     747
     748                         x  = particles(n)%x - ii * dx
     749                         y  = particles(n)%y - jj * dy
     750                         aa = x**2          + y**2
     751                         bb = ( dx - x )**2 + y**2
     752                         cc = x**2          + ( dy - y )**2
     753                         dd = ( dx - x )**2 + ( dy - y )**2
     754                         gg = aa + bb + cc + dd
     755
     756                         ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
     757                                                              ql(kk,jj,ii+1)   &
     758                                    + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
     759                                                              ql(kk,jj+1,ii+1) &
     760                                    ) / ( 3.0 * gg )
     761
     762                         ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
     763                                                            ql(kk+1,jj,ii+1)   &
     764                                    + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
     765                                                            ql(kk+1,jj+1,ii+1) &
     766                                    ) / ( 3.0 * gg )
     767
     768                         ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
     769                                             ( ql_int_u - ql_int_l )
     770
     771!
     772!--                      Interpolate u velocity-component
     773                         ii = ( particles(n)%x + 0.5 * dx ) * ddx
     774                         jj =   particles(n)%y * ddy
     775                         kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist
     776
     777                         IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1
     778
     779                         x  = particles(n)%x + ( 0.5 - ii ) * dx
     780                         y  = particles(n)%y - jj * dy
     781                         aa = x**2          + y**2
     782                         bb = ( dx - x )**2 + y**2
     783                         cc = x**2          + ( dy - y )**2
     784                         dd = ( dx - x )**2 + ( dy - y )**2
     785                         gg = aa + bb + cc + dd
     786
     787                         u_int_l = ( (gg-aa) * u(kk,jj,ii)   + (gg-bb) *       &
     788                                                              u(kk,jj,ii+1)    &
     789                                   + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
     790                                                              u(kk,jj+1,ii+1)  &
     791                                   ) / ( 3.0 * gg ) - u_gtrans
     792                         IF ( kk+1 == nzt+1 )  THEN
     793                            u_int = u_int_l
     794                         ELSE
     795                            u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
     796                                                             u(kk+1,jj,ii+1)   &
     797                                      + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
     798                                                             u(kk+1,jj+1,ii+1) &
     799                                      ) / ( 3.0 * gg ) - u_gtrans
     800                            u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
     801                                              * ( u_int_u - u_int_l )
     802                         ENDIF
     803
     804!
     805!--                      Same procedure for interpolation of the v velocity-com-
     806!--                      ponent (adopt index k from u velocity-component)
     807                         ii =   particles(n)%x * ddx
     808                         jj = ( particles(n)%y + 0.5 * dy ) * ddy
     809
     810                         x  = particles(n)%x - ii * dx
     811                         y  = particles(n)%y + ( 0.5 - jj ) * dy
     812                         aa = x**2          + y**2
     813                         bb = ( dx - x )**2 + y**2
     814                         cc = x**2          + ( dy - y )**2
     815                         dd = ( dx - x )**2 + ( dy - y )**2
     816                         gg = aa + bb + cc + dd
     817
     818                         v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *   &
     819                                                               v(kk,jj,ii+1)   &
     820                                   + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
     821                                                               v(kk,jj+1,ii+1) &
     822                                   ) / ( 3.0 * gg ) - v_gtrans
     823                         IF ( kk+1 == nzt+1 )  THEN
     824                            v_int = v_int_l
     825                         ELSE
     826                            v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
     827                                                               v(kk+1,jj,ii+1) &
     828                                      + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
     829                                                             v(kk+1,jj+1,ii+1) &
     830                                      ) / ( 3.0 * gg ) - v_gtrans
     831                            v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
     832                                              * ( v_int_u - v_int_l )
     833                         ENDIF
     834
     835!
     836!--                      Same procedure for interpolation of the w velocity-com-
     837!--                      ponent (adopt index i from v velocity-component)
     838                         jj = particles(n)%y * ddy
     839                         kk = particles(n)%z / dz
     840
     841                         x  = particles(n)%x - ii * dx
     842                         y  = particles(n)%y - jj * dy
     843                         aa = x**2          + y**2
     844                         bb = ( dx - x )**2 + y**2
     845                         cc = x**2          + ( dy - y )**2
     846                         dd = ( dx - x )**2 + ( dy - y )**2
     847                         gg = aa + bb + cc + dd
     848
     849                         w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *   &
     850                                                                 w(kk,jj,ii+1) &
     851                                   + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
     852                                                               w(kk,jj+1,ii+1) &
     853                                   ) / ( 3.0 * gg )
     854                         IF ( kk+1 == nzt+1 )  THEN
     855                            w_int = w_int_l
     856                         ELSE
     857                            w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
     858                                                               w(kk+1,jj,ii+1) &
     859                                      + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
     860                                                             w(kk+1,jj+1,ii+1) &
     861                                      ) / ( 3.0 * gg )
     862                            w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
     863                                              * ( w_int_u - w_int_l )
     864                         ENDIF
     865
     866!
     867!--                      Change in radius due to collision
     868                         delta_r = effective_coll_efficiency / 3.0             &
     869                                   * pi * sl_r3 * ddx * ddy / dz               &
     870                                   * SQRT( ( u_int - particles(n)%speed_x )**2 &
     871                                         + ( v_int - particles(n)%speed_y )**2 &
     872                                         + ( w_int - particles(n)%speed_z )**2 &
     873                                         ) * dt_3d
     874!
     875!--                      Change in volume due to collision
     876                         delta_v = particles(n)%weight_factor                  &
     877                                   * ( ( particles(n)%radius + delta_r )**3    &
     878                                       - particles(n)%radius**3 )
     879
     880!
     881!--                      Check if collected particles provide enough LWC for
     882!--                      volume change of collector particle
     883                         IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0 )  THEN
     884
     885                            delta_r = ( ( sl_r3/particles(n)%weight_factor )   &
     886                                        + particles(n)%radius**3 )**( 1./3. )  &
     887                                        - particles(n)%radius
     888
     889                            DO  is = n-1, psi, -1
     890                               IF ( particles(is)%radius < &
     891                                    particles(n)%radius )  THEN
     892                                  particles(is)%weight_factor = 0.0
     893                                  particle_mask(is) = .FALSE.
     894                                  deleted_particles = deleted_particles + 1
     895                               ENDIF
     896                            ENDDO
     897
     898                         ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0 )  THEN
     899
     900                            DO  is = n-1, psi, -1
     901                               IF ( particles(is)%radius < particles(n)%radius &
     902                                    .AND.  sl_r3 > 0.0 )  THEN
     903                                  particles(is)%weight_factor =                &
     904                                               ( ( particles(is)%weight_factor &
     905                                               * ( particles(is)%radius**3 ) ) &
     906                                               - ( delta_v                     &
     907                                               * particles(is)%weight_factor   &
     908                                               * ( particles(is)%radius**3 )   &
     909                                               / sl_r3 ) )                     &
     910                                               / ( particles(is)%radius**3 )
     911
     912                                  IF ( particles(is)%weight_factor < 0.0 )  THEN
     913                                     WRITE( message_string, * ) 'negative ',   &
    681914                                                     'weighting factor: ',     &
    682915                                                     particles(is)%weight_factor
    683                                   CALL message( 'advec_particles', '', 2, 2,   &
    684                                                 -1, 6, 1 )
    685 
    686                             ELSEIF ( particles(is)%weight_factor == 0.0 )  THEN
    687 
    688                                 particles(is)%weight_factor = 0.0
    689                                 particle_mask(is)  = .FALSE.
    690                                 deleted_particles = deleted_particles + 1
    691 
    692                             ENDIF
    693 
    694                          ENDDO
    695 
    696                       ENDIF
    697 
    698                       particles(n)%radius = particles(n)%radius + delta_r
    699                       ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor &
    700                                                     * particles(n)%radius**3
    701 
    702                    ENDDO
     916                                     CALL message( 'advec_particles', '', 2,   &
     917                                                   2, -1, 6, 1 )
     918                                  ENDIF
     919                               ENDIF
     920                            ENDDO
     921                         ENDIF
     922
     923                         particles(n)%radius = particles(n)%radius + delta_r
     924                         ql_vp(k,j,i) = ql_vp(k,j,i) +               &
     925                                        particles(n)%weight_factor * &
     926                                        ( particles(n)%radius**3 )
     927                      ENDDO
     928
     929                   ENDIF  ! collision kernel
    703930
    704931                   ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
    705932                                                 * particles(psi)%radius**3
    706933
    707                    DEALLOCATE( kern )
    708934
    709935                ELSE IF ( prt_count(k,j,i) == 1 )  THEN
     
    735961       ENDDO
    736962
    737        ELSE
    738 
    739        DO  i = nxl, nxr
    740           DO  j = nys, nyn
    741              DO  k = nzb+1, nzt
    742 !
    743 !--             Collision requires at least two particles in the box
    744                 IF ( prt_count(k,j,i) > 1 )  THEN
    745 !
    746 !--                First, sort particles within the gridbox by their size,
    747 !--                using Shell's method (see Numerical Recipes)
    748 !--                NOTE: In case of using particle tails, the re-sorting of
    749 !--                ----  tails would have to be included here!
    750                    psi = prt_start_index(k,j,i) - 1
    751                    inc = 1
    752                    DO WHILE ( inc <= prt_count(k,j,i) )
    753                       inc = 3 * inc + 1
    754                    ENDDO
    755 
    756                    DO WHILE ( inc > 1 )
    757                       inc = inc / 3
    758                       DO  is = inc+1, prt_count(k,j,i)
    759                          tmp_particle = particles(psi+is)
    760                          js = is
    761                          DO WHILE ( particles(psi+js-inc)%radius >             &
    762                                     tmp_particle%radius )
    763                             particles(psi+js) = particles(psi+js-inc)
    764                             js = js - inc
    765                             IF ( js <= inc )  EXIT
    766                          ENDDO
    767                          particles(psi+js) = tmp_particle
    768                       ENDDO
    769                    ENDDO
    770 
    771 !
    772 !--                Calculate the mean radius of all those particles which
    773 !--                are of smaller size than the current particle
    774 !--                and use this radius for calculating the collision efficiency
    775                    psi = prt_start_index(k,j,i)
    776 
    777                    DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
    778                       sl_r3 = 0.0
    779                       sl_r4 = 0.0
    780 
    781                       DO is = n-1, psi, -1
    782                          IF ( particles(is)%radius < particles(n)%radius ) THEN
    783                             sl_r3 = sl_r3 + particles(is)%weight_factor        &
    784                                           *( particles(is)%radius**3 )
    785                             sl_r4 = sl_r4 + particles(is)%weight_factor        &
    786                                           *( particles(is)%radius**4 )
    787                          ENDIF
    788                       ENDDO
    789 
    790                       IF ( ( sl_r3 ) > 0.0 )  THEN
    791                          mean_r = ( sl_r4 ) / ( sl_r3 )
    792 
    793                          CALL collision_efficiency( mean_r,                    &
    794                                                     particles(n)%radius,       &
    795                                                     effective_coll_efficiency )
    796 
    797                       ELSE
    798                          effective_coll_efficiency = 0.0
    799                       ENDIF
    800 
    801                       IF ( effective_coll_efficiency > 1.0  .OR.               &
    802                            effective_coll_efficiency < 0.0 )                   &
    803                       THEN
    804                          WRITE( message_string, * )  'collision_efficiency ' , &
    805                                    'out of range:' ,effective_coll_efficiency
    806                          CALL message( 'advec_particles', 'PA0145', 2, 2, -1,  &
    807                                        6, 1 )
    808                       ENDIF
    809 
    810 !
    811 !--                   Interpolation of ...
    812                       ii = particles(n)%x * ddx
    813                       jj = particles(n)%y * ddy
    814                       kk = ( particles(n)%z + 0.5 * dz ) / dz
    815 
    816                       x  = particles(n)%x - ii * dx
    817                       y  = particles(n)%y - jj * dy
    818                       aa = x**2          + y**2
    819                       bb = ( dx - x )**2 + y**2
    820                       cc = x**2          + ( dy - y )**2
    821                       dd = ( dx - x )**2 + ( dy - y )**2
    822                       gg = aa + bb + cc + dd
    823 
    824                       ql_int_l = ( ( gg-aa ) * ql(kk,jj,ii)   + ( gg-bb ) *    &
    825                                                               ql(kk,jj,ii+1)   &
    826                                  + ( gg-cc ) * ql(kk,jj+1,ii) + ( gg-dd ) *    &
    827                                                               ql(kk,jj+1,ii+1) &
    828                                  ) / ( 3.0 * gg )
    829 
    830                       ql_int_u = ( ( gg-aa ) * ql(kk+1,jj,ii)   + ( gg-bb ) *  &
    831                                                             ql(kk+1,jj,ii+1)   &
    832                                  + ( gg-cc ) * ql(kk+1,jj+1,ii) + ( gg-dd ) *  &
    833                                                             ql(kk+1,jj+1,ii+1) &
    834                                  ) / ( 3.0 * gg )
    835 
    836                       ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *   &
    837                                           ( ql_int_u - ql_int_l )
    838 
    839 !
    840 !--                   Interpolate u velocity-component
    841                       ii = ( particles(n)%x + 0.5 * dx ) * ddx
    842                       jj =   particles(n)%y * ddy
    843                       kk = ( particles(n)%z + 0.5 * dz ) / dz  ! only if eq.dist
    844 
    845                       IF ( ( particles(n)%z - zu(kk) ) > ( 0.5*dz ) )  kk = kk+1
    846 
    847                       x  = particles(n)%x + ( 0.5 - ii ) * dx
    848                       y  = particles(n)%y - jj * dy
    849                       aa = x**2          + y**2
    850                       bb = ( dx - x )**2 + y**2
    851                       cc = x**2          + ( dy - y )**2
    852                       dd = ( dx - x )**2 + ( dy - y )**2
    853                       gg = aa + bb + cc + dd
    854 
    855                       u_int_l = ( ( gg-aa ) * u(kk,jj,ii)   + ( gg-bb ) *      &
    856                                                               u(kk,jj,ii+1)    &
    857                                 + ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) *      &
    858                                                               u(kk,jj+1,ii+1)  &
    859                                 ) / ( 3.0 * gg ) - u_gtrans
    860                       IF ( kk+1 == nzt+1 )  THEN
    861                          u_int = u_int_l
    862                       ELSE
    863                          u_int_u = ( ( gg-aa ) * u(kk+1,jj,ii)   + ( gg-bb ) * &
    864                                                              u(kk+1,jj,ii+1)   &
    865                                    + ( gg-cc ) * u(kk+1,jj+1,ii) + ( gg-dd ) * &
    866                                                              u(kk+1,jj+1,ii+1) &
    867                                    ) / ( 3.0 * gg ) - u_gtrans
    868                          u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz *  &
    869                                            ( u_int_u - u_int_l )
    870                       ENDIF
    871 
    872 !
    873 !--                   Same procedure for interpolation of the v velocity-compo-
    874 !--                   nent (adopt index k from u velocity-component)
    875                       ii =   particles(n)%x * ddx
    876                       jj = ( particles(n)%y + 0.5 * dy ) * ddy
    877 
    878                       x  = particles(n)%x - ii * dx
    879                       y  = particles(n)%y + ( 0.5 - jj ) * dy
    880                       aa = x**2          + y**2
    881                       bb = ( dx - x )**2 + y**2
    882                       cc = x**2          + ( dy - y )**2
    883                       dd = ( dx - x )**2 + ( dy - y )**2
    884                       gg = aa + bb + cc + dd
    885 
    886                       v_int_l = ( ( gg-aa ) * v(kk,jj,ii)   + ( gg-bb ) *      &
    887                                                                v(kk,jj,ii+1)   &
    888                                 + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *      &
    889                                                                v(kk,jj+1,ii+1) &
    890                                 ) / ( 3.0 * gg ) - v_gtrans
    891                       IF ( kk+1 == nzt+1 )  THEN
    892                          v_int = v_int_l
    893                       ELSE
    894                          v_int_u = ( ( gg-aa ) * v(kk+1,jj,ii)   + ( gg-bb ) * &
    895                                                                v(kk+1,jj,ii+1) &
    896                                    + ( gg-cc ) * v(kk+1,jj+1,ii) + ( gg-dd ) * &
    897                                                              v(kk+1,jj+1,ii+1) &
    898                                 ) / ( 3.0 * gg ) - v_gtrans
    899                          v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz * &
    900                                            ( v_int_u - v_int_l )
    901                       ENDIF
    902 
    903 !
    904 !--                   Same procedure for interpolation of the w velocity-compo-
    905 !--                   nent (adopt index i from v velocity-component)
    906                       jj = particles(n)%y * ddy
    907                       kk = particles(n)%z / dz
    908 
    909                       x  = particles(n)%x - ii * dx
    910                       y  = particles(n)%y - jj * dy
    911                       aa = x**2          + y**2
    912                       bb = ( dx - x )**2 + y**2
    913                       cc = x**2          + ( dy - y )**2
    914                       dd = ( dx - x )**2 + ( dy - y )**2
    915                       gg = aa + bb + cc + dd
    916 
    917                       w_int_l = ( ( gg-aa ) * w(kk,jj,ii)   + ( gg-bb ) *      &
    918                                                                  w(kk,jj,ii+1) &
    919                                 + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *      &
    920                                                                w(kk,jj+1,ii+1) &
    921                                 ) / ( 3.0 * gg )
    922                       IF ( kk+1 == nzt+1 )  THEN
    923                          w_int = w_int_l
    924                       ELSE
    925                          w_int_u = ( ( gg-aa ) * w(kk+1,jj,ii)   + ( gg-bb ) * &
    926                                                                w(kk+1,jj,ii+1) &
    927                                    + ( gg-cc ) * w(kk+1,jj+1,ii) + ( gg-dd ) * &
    928                                                              w(kk+1,jj+1,ii+1) &
    929                                    ) / ( 3.0 * gg )
    930                          w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz * &
    931                                            ( w_int_u - w_int_l )
    932                       ENDIF
    933 
    934 !
    935 !--                Change in radius due to collision
    936                       delta_r = effective_coll_efficiency / 3.0                &
    937                                 * pi * sl_r3 * ddx * ddy / dz                  &
    938                                 * SQRT( ( u_int - particles(n)%speed_x )**2    &
    939                                       + ( v_int - particles(n)%speed_y )**2    &
    940                                       + ( w_int - particles(n)%speed_z )**2    &
    941                                       ) * dt_3d
    942 !
    943 !--                Change in volume due to collision
    944                       delta_v = particles(n)%weight_factor                     &
    945                                 * ( ( particles(n)%radius + delta_r )**3       &
    946                                     - particles(n)%radius**3 )
    947 
    948 !
    949 !--                Check if collected particles provide enough LWC for volume
    950 !--                change of collector particle
    951                       IF ( delta_v >= sl_r3 .and. sl_r3 > 0.0 ) THEN
    952 
    953                          delta_r = ( ( sl_r3/particles(n)%weight_factor )      &
    954                                      + particles(n)%radius**3 )**( 1./3. )     &
    955                                      - particles(n)%radius
    956 
    957                          DO is = n-1, psi, -1
    958                             IF ( particles(is)%radius < particles(n)%radius )  &
    959                             THEN
    960                                  particles(is)%weight_factor = 0.0
    961                                  particle_mask(is)  = .FALSE.
    962                                  deleted_particles = deleted_particles + 1
    963                             ENDIF
    964                          ENDDO
    965 
    966                       ELSE IF ( delta_v < sl_r3 .and. sl_r3 > 0.0 ) THEN
    967 
    968                          DO is = n-1, psi, -1
    969                             IF ( particles(is)%radius < particles(n)%radius    &
    970                                  .and. sl_r3 > 0.0 ) THEN
    971                                  particles(is)%weight_factor =                 &
    972                                                ( ( particles(is)%weight_factor &
    973                                                * ( particles(is)%radius**3 ) ) &
    974                                                - ( delta_v                     &
    975                                                * particles(is)%weight_factor   &
    976                                                * ( particles(is)%radius**3 )   &
    977                                                / sl_r3 ) )                     &
    978                                                / ( particles(is)%radius**3 )
    979 
    980                                IF (particles(is)%weight_factor < 0.0) THEN
    981                                   WRITE( message_string, * ) 'negative ',      &
    982                                                      'weighting factor: ',     &
    983                                                      particles(is)%weight_factor
    984                                   CALL message( 'advec_particles', '', 2, 2,   &
    985                                                 -1, 6, 1 )
    986                                ENDIF
    987                             ENDIF
    988                          ENDDO
    989                       ENDIF
    990 
    991                       particles(n)%radius = particles(n)%radius + delta_r
    992                       ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor &
    993                                       *( particles(n)%radius**3 )
    994                    ENDDO
    995 
    996                    ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
    997                                   *( particles(psi)%radius**3 )
    998 
    999                 ELSE IF ( prt_count(k,j,i) == 1 )  THEN
    1000                    psi = prt_start_index(k,j,i)
    1001                    ql_vp(k,j,i) =  particles(psi)%weight_factor                &
    1002                                   *( particles(psi)%radius**3 )
    1003                 ENDIF
    1004 !
    1005 !--            Check if condensation of LWC was conserved
    1006 !--            during collision process
    1007                 IF (ql_v(k,j,i) /= 0.0 ) THEN
    1008                    IF (ql_vp(k,j,i)/ql_v(k,j,i) >= 1.00001 .or.                &
    1009                        ql_vp(k,j,i)/ql_v(k,j,i) <= 0.99999  ) THEN
    1010                       WRITE( message_string, * ) 'LWC is not conserved during',&
    1011                                                  ' collision! ',               &
    1012                                                  'LWC after condensation: ',   &
    1013                                                  ql_v(k,j,i),                  &
    1014                                                  ' LWC after collision: ',     &
    1015                                                  ql_vp(k,j,i)
    1016                       CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 )
    1017                    ENDIF
    1018                 ENDIF
    1019 
    1020              ENDDO
    1021           ENDDO
    1022        ENDDO
    1023 
    1024        ENDIF
     963       ENDIF   ! collision handling
    1025964
    1026965       CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
    1027966
    1028     ENDIF
     967    ENDIF   ! cloud droplet handling
    1029968
    1030969
     
    22722211             particles(n)%speed_x = particles(n)%speed_x * exp_term + &
    22732212                                    u_int * ( 1.0 - exp_term )
    2274               particles(n)%speed_y = particles(n)%speed_y * exp_term + &
     2213             particles(n)%speed_y = particles(n)%speed_y * exp_term + &
    22752214                                    v_int * ( 1.0 - exp_term )
    2276               particles(n)%speed_z = particles(n)%speed_z * exp_term +       &
    2277                               ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
    2278                                     * ( 1.0 - exp_term )
     2215             particles(n)%speed_z = particles(n)%speed_z * exp_term +       &
     2216                             ( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
     2217                           * ( 1.0 - exp_term )
    22792218          ENDIF
    22802219
  • palm/trunk/SOURCE/check_parameters.f90

    r824 r825  
    44! Current revisions:
    55! -----------------
    6 ! check for curvature_solution_effects
     6! check for collision_kernel and curvature_solution_effects
    77!
    88! Former revisions:
     
    709709    ENDIF
    710710
     711!
     712!-- Collision kernels:
     713    SELECT CASE ( TRIM( collision_kernel ) )
     714
     715       CASE ( 'hall' )
     716          hall_kernel = .TRUE.
     717
     718       CASE ( 'palm' )
     719          palm_kernel = .TRUE.
     720
     721       CASE ( 'wang' )
     722          wang_kernel = .TRUE.
     723
     724       CASE ( 'none' )
     725
     726
     727       CASE DEFAULT
     728          message_string = 'unknown collision kernel: collision_kernel = "' // &
     729                           TRIM( collision_kernel ) // '"'
     730          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
     731
     732    END SELECT
     733
     734
    711735    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
    712736         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
  • palm/trunk/SOURCE/data_output_ptseries.f90

    r824 r825  
    44! Current revisions:
    55! -----------------
    6 ! mean particle radius added as output quantity.
     6! mean/minimum/maximum particle radius added as output quantity.
    77! particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
    88!
     
    3535!------------------------------------------------------------------------------!
    3636
     37    USE cloud_parameters
    3738    USE control_parameters
    3839    USE cpulog
     
    4849    INTEGER ::  i, inum, j, n
    4950
    50     REAL, DIMENSION(0:number_of_particle_groups,30) ::  pts_value, pts_value_l
    51 
     51    REAL, DIMENSION(:,:), ALLOCATABLE ::  pts_value, pts_value_l
    5252
    5353
     
    7070    ENDIF
    7171
     72    ALLOCATE( pts_value(0:number_of_particle_groups,dopts_num), &
     73              pts_value_l(0:number_of_particle_groups,dopts_num) )
     74
    7275    pts_value_l = 0.0
     76    pts_value_l(:,16) = 9999999.9    ! for calculation of minimum radius
    7377
    7478!
     
    8892       pts_value_l(0,7)  = pts_value_l(0,7) + particles(n)%speed_y     ! mean v
    8993       pts_value_l(0,8)  = pts_value_l(0,8) + particles(n)%speed_z     ! mean w
    90        pts_value_l(0,9)  = pts_value_l(0,9) + particles(n)%rvar1     ! mean sgsu
    91        pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2    ! mean sgsv
    92        pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3    ! mean sgsw
     94       IF ( .NOT. curvature_solution_effects )  THEN
     95          pts_value_l(0,9)  = pts_value_l(0,9) + particles(n)%rvar1  ! mean sgsu
     96          pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
     97          pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
     98       ENDIF
    9399       IF ( particles(n)%speed_z > 0.0 )  THEN
    94100          pts_value_l(0,12) = pts_value_l(0,12) + 1.0  ! # of upward moving prts
     
    99105                                            particles(n)%speed_z ! mean w down
    100106       ENDIF
    101        pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius
    102        pts_value_l(0,16) = number_of_particles
    103        pts_value_l(0,17) = number_of_particles
     107       pts_value_l(0,15) = pts_value_l(0,15) + particles(n)%radius ! mean rad
     108       pts_value_l(0,16) = MIN( pts_value_l(0,16), particles(n)%radius ) ! minrad
     109       pts_value_l(0,17) = MAX( pts_value_l(0,17), particles(n)%radius ) ! maxrad
     110       pts_value_l(0,18) = number_of_particles
     111       pts_value_l(0,19) = number_of_particles
    104112
    105113!
     
    119127          pts_value_l(j,7)  = pts_value_l(j,7)  + particles(n)%speed_y
    120128          pts_value_l(j,8)  = pts_value_l(j,8)  + particles(n)%speed_z
    121           pts_value_l(j,9)  = pts_value_l(j,9)  + particles(n)%rvar1
    122           pts_value_l(j,10) = pts_value_l(j,10) + particles(n)%rvar2
    123           pts_value_l(j,11) = pts_value_l(j,11) + particles(n)%rvar3
     129          IF ( .NOT. curvature_solution_effects )  THEN
     130             pts_value_l(j,9)  = pts_value_l(j,9)  + particles(n)%rvar1
     131             pts_value_l(j,10) = pts_value_l(j,10) + particles(n)%rvar2
     132             pts_value_l(j,11) = pts_value_l(j,11) + particles(n)%rvar3
     133          ENDIF
    124134          IF ( particles(n)%speed_z > 0.0 )  THEN
    125135             pts_value_l(j,12) = pts_value_l(j,12) + 1.0
     
    129139          ENDIF
    130140          pts_value_l(j,15) = pts_value_l(j,15) + particles(n)%radius
    131           pts_value_l(j,16) = pts_value_l(j,16) + 1.0
    132           pts_value_l(j,17) = pts_value_l(j,17) + 1.0
     141          pts_value_l(j,16) = MIN( pts_value(j,16), particles(n)%radius )
     142          pts_value_l(j,17) = MAX( pts_value(j,17), particles(n)%radius )
     143          pts_value_l(j,18) = pts_value_l(j,18) + 1.0
     144          pts_value_l(j,19) = pts_value_l(j,19) + 1.0
    133145
    134146       ENDIF
     
    146158    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    147159    CALL MPI_ALLREDUCE( pts_value_l(0,16), pts_value(0,16), inum, MPI_REAL, &
     160                        MPI_MIN, comm2d, ierr )
     161    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     162    CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, &
    148163                        MPI_MAX, comm2d, ierr )
    149164    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    150     CALL MPI_ALLREDUCE( pts_value_l(0,17), pts_value(0,17), inum, MPI_REAL, &
     165    CALL MPI_ALLREDUCE( pts_value_l(0,18), pts_value(0,18), inum, MPI_REAL, &
     166                        MPI_MAX, comm2d, ierr )
     167    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     168    CALL MPI_ALLREDUCE( pts_value_l(0,19), pts_value(0,19), inum, MPI_REAL, &
    151169                        MPI_MIN, comm2d, ierr )
    152170#else
    153     pts_value(:,1:17) = pts_value_l(:,1:17)
     171    pts_value(:,1:19) = pts_value_l(:,1:19)
    154172#endif
    155173
    156174!
    157 !-- Normalize the above calculated quantities with the total number of
    158 !-- particles
     175!-- Normalize the above calculated quantities (except min/max values) with the
     176!-- total number of particles
    159177    IF ( number_of_particle_groups > 1 )  THEN
    160178       inum = number_of_particle_groups
     
    186204    DO  n = 1, number_of_particles
    187205
    188        pts_value_l(0,18) = pts_value_l(0,18) + ( particles(n)%x - &
     206       pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%x - &
    189207                           particles(n)%origin_x - pts_value(0,2) )**2 ! x*2
    190        pts_value_l(0,19) = pts_value_l(0,19) + ( particles(n)%y - &
     208       pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%y - &
    191209                           particles(n)%origin_y - pts_value(0,3) )**2 ! y*2
    192        pts_value_l(0,20) = pts_value_l(0,20) + ( particles(n)%z - &
     210       pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%z - &
    193211                           particles(n)%origin_z - pts_value(0,4) )**2 ! z*2
    194        pts_value_l(0,21) = pts_value_l(0,21) + ( particles(n)%speed_x - &
    195                                                pts_value(0,6) )**2     ! u*2
    196        pts_value_l(0,22) = pts_value_l(0,22) + ( particles(n)%speed_y - &
    197                                                pts_value(0,7) )**2     ! v*2
    198        pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_z - &
    199                                                pts_value(0,8) )**2     ! w*2
    200        pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%rvar1 - &
    201                                                pts_value(0,9) )**2     ! u"2
    202        pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%rvar2 - &
    203                                                pts_value(0,10) )**2    ! v"2
    204        pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar3 - &
    205                                                pts_value(0,11) )**2    ! w"2
     212       pts_value_l(0,23) = pts_value_l(0,23) + ( particles(n)%speed_x - &
     213                                                pts_value(0,6) )**2   ! u*2
     214       pts_value_l(0,24) = pts_value_l(0,24) + ( particles(n)%speed_y - &
     215                                                 pts_value(0,7) )**2   ! v*2
     216       pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
     217                                                 pts_value(0,8) )**2   ! w*2
     218       IF ( .NOT. curvature_solution_effects )  THEN
     219          pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
     220                                                    pts_value(0,9) )**2   ! u"2
     221          pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
     222                                                    pts_value(0,10) )**2  ! v"2
     223          pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
     224                                                    pts_value(0,11) )**2  ! w"2
     225       ENDIF
    206226!
    207227!--    Repeat the same for the respective particle group
     
    209229          j = particles(n)%group
    210230
    211           pts_value_l(j,18) = pts_value_l(j,18) + ( particles(n)%x - &
     231          pts_value_l(j,20) = pts_value_l(j,20) + ( particles(n)%x - &
    212232                              particles(n)%origin_x - pts_value(j,2) )**2
    213           pts_value_l(j,19) = pts_value_l(j,19) + ( particles(n)%y - &
     233          pts_value_l(j,21) = pts_value_l(j,21) + ( particles(n)%y - &
    214234                              particles(n)%origin_y - pts_value(j,3) )**2
    215           pts_value_l(j,20) = pts_value_l(j,20) + ( particles(n)%z - &
     235          pts_value_l(j,22) = pts_value_l(j,22) + ( particles(n)%z - &
    216236                              particles(n)%origin_z - pts_value(j,4) )**2
    217           pts_value_l(j,21) = pts_value_l(j,21) + ( particles(n)%speed_x - &
    218                                                   pts_value(j,6) )**2
    219           pts_value_l(j,22) = pts_value_l(j,22) + ( particles(n)%speed_y - &
    220                                                   pts_value(j,7) )**2
    221           pts_value_l(j,23) = pts_value_l(j,23) + ( particles(n)%speed_z - &
    222                                                   pts_value(j,8) )**2
    223           pts_value_l(j,24) = pts_value_l(j,24) + ( particles(n)%rvar1 -   &
    224                                                   pts_value(j,9) )**2
    225           pts_value_l(j,25) = pts_value_l(j,25) + ( particles(n)%rvar2 -   &
    226                                                   pts_value(j,10) )**2
    227           pts_value_l(j,26) = pts_value_l(j,26) + ( particles(n)%rvar3 -   &
    228                                                   pts_value(j,11) )**2
     237          pts_value_l(j,23) = pts_value_l(j,23) + ( particles(n)%speed_x - &
     238                                                    pts_value(j,6) )**2
     239          pts_value_l(j,24) = pts_value_l(j,24) + ( particles(n)%speed_y - &
     240                                                    pts_value(j,7) )**2
     241          pts_value_l(j,25) = pts_value_l(j,25) + ( particles(n)%speed_z - &
     242                                                    pts_value(j,8) )**2
     243          IF ( .NOT. curvature_solution_effects )  THEN
     244             pts_value_l(j,26) = pts_value_l(j,26) + ( particles(n)%rvar1 - &
     245                                                       pts_value(j,9) )**2
     246             pts_value_l(j,27) = pts_value_l(j,27) + ( particles(n)%rvar2 - &
     247                                                       pts_value(j,10) )**2
     248             pts_value_l(j,28) = pts_value_l(j,28) + ( particles(n)%rvar3 - &
     249                                                       pts_value(j,11) )**2
     250          ENDIF
    229251       ENDIF
    230252
    231253    ENDDO
    232254
    233     pts_value_l(0,27) = ( number_of_particles - pts_value(0,1) / numprocs )**2
     255    pts_value_l(0,29) = ( number_of_particles - pts_value(0,1) / numprocs )**2
    234256                                                 ! variance of particle numbers
    235257    IF ( number_of_particle_groups > 1 )  THEN
    236258       DO  j = 1, number_of_particle_groups
    237           pts_value_l(j,27) = ( pts_value_l(j,1) - &
     259          pts_value_l(j,29) = ( pts_value_l(j,1) - &
    238260                                pts_value(j,1) / numprocs )**2
    239261       ENDDO
     
    246268
    247269    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    248     CALL MPI_ALLREDUCE( pts_value_l(0,18), pts_value(0,18), inum*10, MPI_REAL, &
     270    CALL MPI_ALLREDUCE( pts_value_l(0,20), pts_value(0,20), inum*10, MPI_REAL, &
    249271                        MPI_SUM, comm2d, ierr )
    250272#else
    251     pts_value(:,18:27) = pts_value_l(:,18:27)
     273    pts_value(:,20:29) = pts_value_l(:,20:29)
    252274#endif
    253275
     
    264286
    265287       IF ( pts_value(j,1) > 0.0 )  THEN
    266           pts_value(j,18:26) = pts_value(j,18:26) / pts_value(j,1)
    267        ENDIF
    268        pts_value(j,27) = pts_value(j,27) / numprocs
     288          pts_value(j,20:28) = pts_value(j,20:28) / pts_value(j,1)
     289       ENDIF
     290       pts_value(j,29) = pts_value(j,29) / numprocs
    269291
    270292    ENDDO
     
    286308#endif
    287309
     310    DEALLOCATE( pts_value, pts_value_l )
     311
    288312    CALL cpu_log( log_point(36), 'data_output_ptseries','stop', 'nobarrier' )
    289313
  • palm/trunk/SOURCE/diffusion_e.f90

    r791 r825  
    44! Current revisions:
    55! -----------------
    6 !
     6! wang_collision_kernel renamed wang_kernel
    77!
    88! Former revisions:
     
    164164!--          Store dissipation if needed for calculating the sgs particle
    165165!--          velocities
    166              IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
     166             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
    167167                DO  j = nys, nyn
    168168                   DO  k = nzb_s_inner(j,i)+1, nzt
     
    254254!--          Store dissipation if needed for calculating the sgs particle
    255255!--          velocities
    256              IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
     256             IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
    257257                DO  j = nys, nyn
    258258                   DO  k = nzb_s_inner(j,i)+1, nzt
     
    268268!
    269269!--    Boundary condition for dissipation
    270        IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
     270       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
    271271          DO  i = nxl, nxr
    272272             DO  j = nys, nyn
     
    372372!
    373373!--    Store dissipation if needed for calculating the sgs particle velocities
    374        IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
     374       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
    375375          DO  k = nzb_s_inner(j,i)+1, nzt
    376376             diss(k,j,i) = dissipation(k)
  • palm/trunk/SOURCE/header.f90

    r824 r825  
    14781478       WRITE ( io, 433 )
    14791479       IF ( curvature_solution_effects )  WRITE ( io, 434 )
     1480       IF ( collision_kernel /= 'none' )  THEN
     1481          WRITE ( io, 435 )  TRIM( collision_kernel )
     1482       ELSE
     1483          WRITE ( io, 436 )
     1484       ENDIF
    14801485    ENDIF
    14811486
     
    19551960434 FORMAT ('    Curvature and solution effecs are considered for growth of', &
    19561961                 ' droplets < 1.0E-6 m')
     1962435 FORMAT ('    Droplet collision is handled by ',A,'-kernel')
     1963436 FORMAT ('    Droplet collision is switched off')
    19571964450 FORMAT (//' LES / Turbulence quantities:'/ &
    19581965              ' ---------------------------'/)
  • palm/trunk/SOURCE/init_3d_model.f90

    r791 r825  
    77! Current revisions:
    88! ------------------
    9 !
     9! wang_collision_kernel renamed wang_kernel
    1010!
    1111! Former revisions:
     
    334334!-- 3D-array for storing the dissipation, needed for calculating the sgs
    335335!-- particle velocities
    336     IF ( use_sgs_for_particles  .OR.  wang_collision_kernel )  THEN
     336    IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
    337337       ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    338338    ELSE
  • palm/trunk/SOURCE/lpm_collision_kernels.f90

    r824 r825  
    1 MODULE wang_kernel_mod
     1MODULE lpm_collision_kernels_mod
    22
    33!------------------------------------------------------------------------------!
    44! Current revisions:
    55! -----------------
    6 !
     6! routine renamed from wang_kernel to lpm_collision_kernels,
     7! turbulence_effects on collision replaced by wang_kernel
    78!
    89! Former revisions:
     
    106107       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
    107108
    108        IF ( turbulence_effects_on_collision )  THEN
     109       IF ( wang_kernel )  THEN
    109110
    110111          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
     
    756757   END SUBROUTINE turb_enhance_eff
    757758
    758  END MODULE wang_kernel_mod
     759 END MODULE lpm_collision_kernels_mod
  • palm/trunk/SOURCE/modules.f90

    r824 r825  
    66! +bfactor, curvature_solution_effects, eps_ros, molecular_weight_of_water,
    77! vanthoff, -b_cond in cloud_parameters,
    8 ! dopts_num increased to 27, particle attributes speed_x|y|z_sgs renamed
     8! dopts_num increased to 29, particle attributes speed_x|y|z_sgs renamed
    99! rvar1|2|3
     10! wang_collision_kernel and turbulence_effects_on_collision replaced by
     11! collision_kernel, hall_kernel, palm_kernel, wang_kernel
    1012!
    1113! Former revisions:
     
    10471049#endif
    10481050
    1049     INTEGER, PARAMETER ::  dopr_norm_num = 7, dopts_num = 27, dots_max = 100
     1051    INTEGER, PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100
    10501052
    10511053    INTEGER ::  dots_num = 23
     
    10621064          (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
    10631065             'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
    1064              'w_up   ', 'w_down ', 'radius ', 'npt_max', 'npt_min', 'x*2    ', &
    1065              'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', 'w*2    ', 'u"2    ', &
    1066              'v"2    ', 'w"2    ', 'npt*2  ' /)
     1066             'w_up   ', 'w_down ', 'radius ', 'r_min  ', 'r_max  ', 'npt_max', &
     1067             'npt_min', 'x*2    ', 'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', &
     1068             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
    10671069
    10681070    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit = &
    10691071          (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
    10701072             'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
    1071              'm/s    ', 'm/s    ', 'm      ', 'number ', 'number ', 'm2    ', &
    1072              'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', &
    1073              'm2/s2  ', 'm2/s2  ', 'number2' /)
     1073             'm/s    ', 'm/s    ', 'm      ', 'm      ', 'm      ', 'number ', &
     1074             'number ', 'm2     ', 'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', &
     1075             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
    10741076
    10751077    CHARACTER (LEN=7), DIMENSION(dots_max) :: dots_label = &
     
    11841186
    11851187    CHARACTER (LEN=15)  ::  bc_par_lr = 'cyclic',  bc_par_ns = 'cyclic', &
    1186                             bc_par_b  = 'reflect', bc_par_t  = 'absorb'
     1188                            bc_par_b  = 'reflect', bc_par_t  = 'absorb', &
     1189                            collision_kernel = 'none'
    11871190
    11881191#if defined( __parallel )
     
    12061209    INTEGER, DIMENSION(:,:,:), ALLOCATABLE ::  prt_count, prt_start_index
    12071210
    1208     LOGICAL ::  turbulence_effects_on_collision = .FALSE.,                     &
     1211    LOGICAL ::  hall_kernel = .FALSE., palm_kernel = .FALSE.,                  &
    12091212                particle_advection = .FALSE., random_start_position = .FALSE., &
    12101213                read_particles_from_restartfile = .TRUE.,                      &
    12111214                uniform_particles = .TRUE., use_particle_tails = .FALSE.,      &
    12121215                use_sgs_for_particles = .FALSE.,                               &
    1213                 wang_collision_kernel = .FALSE.,                               &
    1214                 write_particle_statistics = .FALSE.
     1216                wang_kernel = .FALSE., write_particle_statistics = .FALSE.
    12151217
    12161218    LOGICAL, DIMENSION(max_number_of_particle_groups) ::                       &
  • palm/trunk/SOURCE/package_parin.f90

    r791 r825  
    44! Current revisions:
    55! -----------------
    6 !
     6! wang_collision_kernel and turbulence_effects_on_collision in particles_par
     7! replaced by collision_kernel
    78!
    89! Former revisions:
     
    7172
    7273    NAMELIST /particles_par/      bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,    &
     74                                  collision_kernel,                            &
    7375                                  density_ratio, radius, dt_dopts,             &
    7476                                  dt_min_part, dt_prel, dt_sort_particles,     &
     
    8789                                  read_particles_from_restartfile,             &
    8890                                  skip_particles_for_tail,                     &
    89                                   turbulence_effects_on_collision,             &
    9091                                  use_particle_tails, use_sgs_for_particles,   &
    9192                                  vertical_particle_advection,                 &
    92                                   wang_collision_kernel,                       &
    9393                                  write_particle_statistics
    9494    NAMELIST /spectra_par/        averaging_interval_sp, comp_spectra_level,   &
  • palm/trunk/SOURCE/time_integration.f90

    r791 r825  
    44! Current revisions:
    55! -----------------
    6 !
     6! wang_collision_kernel renamed wang_kernel
    77!
    88! Former revisions:
     
    218218             CALL exchange_horiz( ql_vp, nbgp )
    219219          ENDIF
    220           IF ( wang_collision_kernel )  CALL exchange_horiz( diss, nbgp )
     220          IF ( wang_kernel )  CALL exchange_horiz( diss, nbgp )
    221221
    222222          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
Note: See TracChangeset for help on using the changeset viewer.