Changeset 2312 for palm/trunk/SOURCE

Timestamp:

Jul 14, 2017 8:26:51 PM (7 years ago)

Author:

hoffmann

Message:

various improvements of the LCM

Location:

palm/trunk/SOURCE

Files:

• check_parameters.f90 (modified) (69 diffs)
• data_output_ptseries.f90 (modified) (8 diffs)
• lpm_droplet_collision.f90 (modified) (10 diffs)
• lpm_droplet_condensation.f90 (modified) (12 diffs)
• lpm_init.f90 (modified) (49 diffs)
• lpm_read_restart_file.f90 (modified) (5 diffs)
• microphysics_mod.f90 (modified) (43 diffs)
• mod_particle_attributes.f90 (modified) (9 diffs)
• package_parin.f90 (modified) (12 diffs)

palm/trunk/SOURCE/check_parameters.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! PA0349 and PA0420 removed.
+!
+! 2300 2017-06-29 13:31:14Z raasch
 ! host-specific settings and checks removed
-! 
+!
 ! 2292 2017-06-20 09:51:42Z schwenkel
-! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
-! includes two more prognostic equations for cloud drop concentration (nc)  
-! and cloud water content (qc). 
-! 
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
+! includes two more prognostic equations for cloud drop concentration (nc)
+! and cloud water content (qc).
+!
 ! 2274 2017-06-09 13:27:48Z Giersch
-! Changed error messages 
-! 
+! Changed error messages
+!
 ! 2271 2017-06-09 12:34:55Z sward
 ! roughness-length check altered
 ! Error messages fixed
-! 
+!
 ! 2270 2017-06-09 12:18:47Z maronga
 ! Revised numbering (removed 2 timeseries)
-! 
+!
 ! 2259 2017-06-08 09:09:11Z gronemeier
 ! Implemented synthetic turbulence generator
@@ -50 +53 @@
 ! Doxygen comment added
 !
-! 2248 2017-06-06 13:52:54Z sward 
+! 2248 2017-06-06 13:52:54Z sward
 ! Error message fixed
 !
 ! 2209 2017-04-19 09:34:46Z kanani
 ! Check for plant canopy model output
-! 
+!
 ! 2200 2017-04-11 11:37:51Z suehring
 ! monotonic_adjustment removed
-! 
+!
 ! 2178 2017-03-17 11:07:39Z hellstea
 ! Index limits for perturbations are now set also in case of nested
@@ -65 +68 @@
 ! 2169 2017-03-06 18:16:35Z suehring
 ! Bugfix, move setting for topography grid convention to init_grid
-! 
+!
 ! 2118 2017-01-17 16:38:49Z raasch
 ! OpenACC related parts of code removed
-! 
+!
 ! 2088 2016-12-19 16:30:25Z suehring
 ! Bugfix in initial salinity profile
-! 
+!
 ! 2084 2016-12-09 15:59:42Z knoop
 ! Removed anelastic + multigrid error message
-! 
+!
 ! 2052 2016-11-08 15:14:59Z gronemeier
 ! Bugfix: remove setting of default value for recycling_width
-! 
+!
 ! 2050 2016-11-08 15:00:55Z gronemeier
 ! Implement turbulent outflow condition
-! 
+!
 ! 2044 2016-11-02 16:44:25Z knoop
 ! Added error code for anelastic approximation
-! 
+!
 ! 2042 2016-11-02 13:47:31Z suehring
 ! Additional checks for wall_heatflux, wall_humidityflux and wall_scalarflux.
 ! Bugfix, check for constant_scalarflux.
-! 
+!
 ! 2040 2016-10-26 16:58:09Z gronemeier
 ! Removed check for statistic_regions > 9.
-! 
+!
 ! 2037 2016-10-26 11:15:40Z knoop
 ! Anelastic approximation implemented
-! 
+!
 ! 2031 2016-10-21 15:11:58Z knoop
 ! renamed variable rho to rho_ocean
-! 
+!
 ! 2026 2016-10-18 10:27:02Z suehring
 ! Bugfix, enable output of s*2
-! 
+!
 ! 2024 2016-10-12 16:42:37Z kanani
 ! Added missing CASE, error message and unit for ssws*,
 ! increased number of possible output quantities to 500.
-! 
+!
 ! 2011 2016-09-19 17:29:57Z kanani
 ! Flag urban_surface is now defined in module control_parameters,
 ! changed prefix for urban surface model output to "usm_",
-! added flag lsf_exception (inipar-Namelist parameter) to allow explicit 
+! added flag lsf_exception (inipar-Namelist parameter) to allow explicit
 ! enabling of large scale forcing together with buildings on flat terrain,
 ! introduced control parameter varnamelength for LEN of var.
-! 
+!
 ! 2007 2016-08-24 15:47:17Z kanani
 ! Added checks for the urban surface model,
 ! increased counter in DO WHILE loop over data_output (for urban surface output)
-! 
+!
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 1994 2016-08-15 09:52:21Z suehring
 ! Add missing check for cloud_physics and cloud_droplets
-! 
+!
 ! 1992 2016-08-12 15:14:59Z suehring
 ! New checks for top_scalarflux
 ! Bugfixes concerning data output of profiles in case of passive_scalar
-! 
+!
 ! 1984 2016-08-01 14:48:05Z suehring
 ! Bugfix: checking for bottom and top boundary condition for humidity and scalars
-! 
+!
 ! 1972 2016-07-26 07:52:02Z maronga
 ! Removed check of lai* (done in land surface module)
-! 
+!
 ! 1970 2016-07-18 14:27:12Z suehring
 ! Bugfix, check for ibc_q_b and constant_waterflux in case of land-surface scheme
-! 
+!
 ! 1962 2016-07-13 09:16:44Z suehring
 ! typo removed
-! 
+!
 ! 1960 2016-07-12 16:34:24Z suehring
 ! Separate checks and calculations for humidity and passive scalar. Therefore,
 ! a subroutine to calculate vertical gradients is introduced, in order  to reduce
-! redundance. 
+! redundance.
 ! Additional check - large-scale forcing is not implemented for passive scalars
 ! so far.
-! 
+!
 ! 1931 2016-06-10 12:06:59Z suehring
 ! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
@@ -158 +161 @@
 ! 1841 2016-04-07 19:14:06Z raasch
 ! redundant location message removed
-! 
+!
 ! 1833 2016-04-07 14:23:03Z raasch
 ! check of spectra quantities moved to spectra_mod
-! 
+!
 ! 1829 2016-04-07 12:16:29Z maronga
 ! Bugfix: output of user defined data required reset of the string 'unit'
-! 
+!
 ! 1826 2016-04-07 12:01:39Z maronga
 ! Moved checks for radiation model to the respective module.
@@ -172 +175 @@
 ! 1824 2016-04-07 09:55:29Z gronemeier
 ! Check if roughness_length < dz/2. Moved location_message(finished) to the end.
-! 
+!
 ! 1822 2016-04-07 07:49:42Z hoffmann
 ! PALM collision kernel deleted. Collision algorithms are checked for correct
@@ -186 +189 @@
 ! 1817 2016-04-06 15:44:20Z maronga
 ! Moved checks for land_surface model to the respective module
-! 
+!
 ! 1806 2016-04-05 18:55:35Z gronemeier
 ! Check for recycling_yshift
@@ -200 +203 @@
 ! surface presribed in the land surface model. Added output of z0q.
 ! Syntax layout improved.
-! 
+!
 ! 1786 2016-03-08 05:49:27Z raasch
 ! cpp-direktives for spectra removed, check of spectra level removed
@@ -217 +220 @@
 ! 1745 2016-02-05 13:06:51Z gronemeier
 ! Bugfix: check data output intervals to be /= 0.0 in case of parallel NetCDF4
-! 
+!
 ! 1701 2015-11-02 07:43:04Z maronga
 ! Bugfix: definition of rad_net timeseries was missing
-! 
+!
 ! 1691 2015-10-26 16:17:44Z maronga
-! Added output of Obukhov length (ol) and radiative heating rates for RRTMG. 
+! Added output of Obukhov length (ol) and radiative heating rates for RRTMG.
 ! Added checks for use of radiation / lsm with topography.
-! 
+!
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
+! Code annotations made doxygen readable
+!
 ! 1606 2015-06-29 10:43:37Z maronga
 ! Added check for use of RRTMG without netCDF.
-! 
+!
 ! 1587 2015-05-04 14:19:01Z maronga
 ! Added check for set of albedos when using RRTMG
-! 
+!
 ! 1585 2015-04-30 07:05:52Z maronga
 ! Added support for RRTMG
@@ -248 +251 @@
 ! 1555 2015-03-04 17:44:27Z maronga
 ! Added output of r_a and r_s. Renumbering of LSM PA-messages.
-! 
+!
 ! 1553 2015-03-03 17:33:54Z maronga
 ! Removed check for missing soil temperature profile as default values were added.
-! 
+!
 ! 1551 2015-03-03 14:18:16Z maronga
 ! Added various checks for land surface and radiation model. In the course of this
 ! action, the length of the variable var has to be increased
-! 
+!
 ! 1504 2014-12-04 09:23:49Z maronga
 ! Bugfix: check for large_scale forcing got lost.
-! 
+!
 ! 1500 2014-12-03 17:42:41Z maronga
 ! Boundary conditions changed to dirichlet for land surface model
@@ -264 +267 @@
 ! 1496 2014-12-02 17:25:50Z maronga
 ! Added checks for the land surface model
-! 
+!
 ! 1484 2014-10-21 10:53:05Z kanani
 ! Changes due to new module structure of the plant canopy model:
 !   module plant_canopy_model_mod added,
-!   checks regarding usage of new method for leaf area density profile 
+!   checks regarding usage of new method for leaf area density profile
 !   construction added,
-!   lad-profile construction moved to new subroutine init_plant_canopy within 
+!   lad-profile construction moved to new subroutine init_plant_canopy within
 !   the module plant_canopy_model_mod,
 !   drag_coefficient renamed to canopy_drag_coeff.
 ! Missing KIND-attribute for REAL constant added
-! 
+!
 ! 1455 2014-08-29 10:47:47Z heinze
-! empty time records in volume, cross-section and masked data output prevented  
+! empty time records in volume, cross-section and masked data output prevented
 ! in case of non-parallel netcdf-output in restart runs
-! 
+!
 ! 1429 2014-07-15 12:53:45Z knoop
 ! run_description_header exended to provide ensemble_member_nr if specified
-! 
+!
 ! 1425 2014-07-05 10:57:53Z knoop
 ! bugfix: perturbation domain modified for parallel random number generator
-! 
+!
 ! 1402 2014-05-09 14:25:13Z raasch
 ! location messages modified
-! 
+!
 ! 1400 2014-05-09 14:03:54Z knoop
 ! Check random generator extended by option random-parallel
-! 
+!
 ! 1384 2014-05-02 14:31:06Z raasch
 ! location messages added
-! 
+!
 ! 1365 2014-04-22 15:03:56Z boeske
 ! Usage of large scale forcing for pt and q enabled:
-! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q), 
-! large scale subsidence (td_sub_lpt, td_sub_q) 
+! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
+! large scale subsidence (td_sub_lpt, td_sub_q)
 ! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
 ! check if use_subsidence_tendencies is used correctly
-! 
+!
 ! 1361 2014-04-16 15:17:48Z hoffmann
 ! PA0363 removed
-! PA0362 changed 
+! PA0362 changed
 !
 ! 1359 2014-04-11 17:15:14Z hoffmann
@@ -310 +313 @@
 !
 ! PA0084 not necessary for new particle structure
-! 
+!
 ! 1353 2014-04-08 15:21:23Z heinze
 ! REAL constants provided with KIND-attribute
-! 
+!
 ! 1330 2014-03-24 17:29:32Z suehring
-! In case of SGS-particle velocity advection of TKE is also allowed with 
+! In case of SGS-particle velocity advection of TKE is also allowed with
 ! dissipative 5th-order scheme.
-! 
+!
 ! 1327 2014-03-21 11:00:16Z raasch
 ! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
@@ -348 +351 @@
 ! 1241 2013-10-30 11:36:58Z heinze
 ! output for profiles of ug and vg added
-! set surface_heatflux and surface_waterflux also in dependence on 
+! set surface_heatflux and surface_waterflux also in dependence on
 ! large_scale_forcing
 ! checks for nudging and large scale forcing from external file
@@ -395 +398 @@
 !
 ! 1065 2012-11-22 17:42:36Z hoffmann
-! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without 
+! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
 !         precipitation in order to save computational resources.
 !
@@ -405 +408 @@
 ! - check cloud physics scheme (Kessler or Seifert and Beheng)
 ! - plant_canopy is not allowed
-! - currently, only cache loop_optimization is allowed 
+! - currently, only cache loop_optimization is allowed
 ! - initial profiles of nr, qr
 ! - boundary condition of nr, qr
@@ -483 +486 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE check_parameters
- 
+
 
     USE arrays_3d
@@ -507 +510 @@
     USE plant_canopy_model_mod,                                                &
         ONLY:  pcm_check_data_output, pcm_check_parameters, plant_canopy
-        
+
     USE pmc_interface,                                                         &
         ONLY:  cpl_id, nested_run
@@ -531 +534 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=1)   ::  sq                       !< 
-    CHARACTER (LEN=varnamelength)  ::  var           !< 
+    CHARACTER (LEN=1)   ::  sq                       !<
+    CHARACTER (LEN=varnamelength)  ::  var           !<
     CHARACTER (LEN=7)   ::  unit                     !<
-    CHARACTER (LEN=8)   ::  date                     !< 
-    CHARACTER (LEN=10)  ::  time                     !< 
+    CHARACTER (LEN=8)   ::  date                     !<
+    CHARACTER (LEN=10)  ::  time                     !<
     CHARACTER (LEN=20)  ::  ensemble_string          !<
     CHARACTER (LEN=15)  ::  nest_string              !<
-    CHARACTER (LEN=40)  ::  coupling_string          !< 
-    CHARACTER (LEN=100) ::  action                   !< 
-
-    INTEGER(iwp) ::  i                               !< 
+    CHARACTER (LEN=40)  ::  coupling_string          !<
+    CHARACTER (LEN=100) ::  action                   !<
+
+    INTEGER(iwp) ::  i                               !<
     INTEGER(iwp) ::  ilen                            !<
-    INTEGER(iwp) ::  j                               !< 
-    INTEGER(iwp) ::  k                               !< 
-    INTEGER(iwp) ::  kk                              !< 
-    INTEGER(iwp) ::  netcdf_data_format_save         !< 
+    INTEGER(iwp) ::  j                               !<
+    INTEGER(iwp) ::  k                               !<
+    INTEGER(iwp) ::  kk                              !<
+    INTEGER(iwp) ::  netcdf_data_format_save         !<
     INTEGER(iwp) ::  position                        !<
-    
-    LOGICAL     ::  found                            !< 
-    
-    REAL(wp)    ::  dum                              !< 
-    REAL(wp)    ::  gradient                         !< 
-    REAL(wp)    ::  remote = 0.0_wp                  !< 
-    REAL(wp)    ::  simulation_time_since_reference  !< 
+
+    LOGICAL     ::  found                            !<
+
+    REAL(wp)    ::  dum                              !<
+    REAL(wp)    ::  gradient                         !<
+    REAL(wp)    ::  remote = 0.0_wp                  !<
+    REAL(wp)    ::  simulation_time_since_reference  !<
 
 
@@ -568 +571 @@
 !-- Check the coupling mode
 !> @todo Check if any queries for other coupling modes (e.g. precursor_ocean) are missing
-    IF ( coupling_mode /= 'uncoupled'            .AND.  &  
-         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &  
-         coupling_mode /= 'ocean_to_atmosphere' )  THEN      
+    IF ( coupling_mode /= 'uncoupled'            .AND.  &
+         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
+         coupling_mode /= 'ocean_to_atmosphere' )  THEN
        message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
        CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
@@ -590 +593 @@
 !--    NOTE: coupled runs have not been implemented in the check_namelist_files
 !--    program.
-!--    check_namelist_files will need the following information of the other 
+!--    check_namelist_files will need the following information of the other
 !--    model (atmosphere/ocean).
 !       dt_coupling = remote
@@ -607 +610 @@
        ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
- 
+
        IF ( dt_coupling /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -620 +623 @@
              CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
                             status, ierr )
-          ENDIF   
+          ENDIF
           CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
-      
+
           dt_coupling = MAX( dt_max, remote )
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -635 +638 @@
           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
                          status, ierr )
-       ENDIF 
+       ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
-  
+
        IF ( restart_time /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -650 +653 @@
           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
                          status, ierr )
-       ENDIF   
+       ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
-    
+
        IF ( dt_restart /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -666 +669 @@
                          14, comm_inter, ierr )
           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
-                         status, ierr )    
-       ENDIF 
+                         status, ierr )
+       ENDIF
        CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
-    
+
        IF ( simulation_time_since_reference /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -730 +733 @@
              WRITE( message_string, * ) 'coupling mode "',                     &
                    TRIM( coupling_mode ),                                      &
-             '": nx+1 in ocean is not divisible by nx+1 in', & 
+             '": nx+1 in ocean is not divisible by nx+1 in', &
              ' atmosphere without remainder'
              CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
@@ -738 +741 @@
              WRITE( message_string, * ) 'coupling mode "',                     &
                    TRIM( coupling_mode ),                                      &
-             '": ny+1 in ocean is not divisible by ny+1 in', & 
+             '": ny+1 in ocean is not divisible by ny+1 in', &
              ' atmosphere without remainder'
 
@@ -763 +766 @@
     ENDIF
     CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
-    
+
 #endif
 
@@ -932 +935 @@
        flux_output_mode = 'dynamic'
     ENDIF
-    
+
 !
 !-- set the flux output units according to flux_output_mode
@@ -1024 +1027 @@
        CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
     ENDIF
-    
-    IF( momentum_advec == 'ws-scheme' .AND.                                    & 
+
+    IF( momentum_advec == 'ws-scheme' .AND.                                    &
         .NOT. call_psolver_at_all_substeps  ) THEN
         message_string = 'psolver must be called at each RK3 substep when "'// &
@@ -1074 +1077 @@
     ENDIF
 
-    IF ( use_sgs_for_particles  .AND.  curvature_solution_effects )  THEN
-       message_string = 'use_sgs_for_particles = .TRUE. not allowed with ' //  &
-                        'curvature_solution_effects = .TRUE.'
-       CALL message( 'check_parameters', 'PA0349', 1, 2, 0, 6, 0 )
-    ENDIF
-
 !
 !-- Set LOGICAL switches to enhance performance
@@ -1114 +1111 @@
     ENDIF
 !
-!-- Check for proper settings for microphysics 
+!-- Check for proper settings for microphysics
     IF ( cloud_physics  .AND.  cloud_droplets )  THEN
        message_string = 'cloud_physics = .TRUE. is not allowed with ' //  &
@@ -1141 +1138 @@
     END SELECT
     IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
-
-!
-!-- Collision algorithms:
-    SELECT CASE ( TRIM( collision_algorithm ) )
-
-       CASE ( 'all_or_nothing' )
-          all_or_nothing = .TRUE.
-
-       CASE ( 'average_impact' )
-          average_impact = .TRUE.
-
-       CASE DEFAULT
-          message_string = 'unknown collision algorithm: collision_algorithm = "' // &
-                           TRIM( collision_algorithm ) // '"'
-          CALL message( 'check_parameters', 'PA0420', 1, 2, 0, 6, 0 )
-
-       END SELECT
 
     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
@@ -1186 +1166 @@
                         ' is not allowed with conserve_volume_flow = .T.'
        CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
-    ENDIF       
+    ENDIF
 
 
@@ -1229 +1209 @@
 !-- When plant canopy model is used, peform addtional checks
     IF ( plant_canopy )  CALL pcm_check_parameters
-    
+
 !
 !-- Additional checks for spectra
@@ -1255 +1235 @@
                         'loop_optimization = "cache" or "vector"'
        CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
-    ENDIF 
+    ENDIF
 
 !
@@ -1288 +1268 @@
                    i = i + 1
                 ENDIF
-             ENDIF        
+             ENDIF
              IF ( gradient /= 0.0_wp )  THEN
                 IF ( k /= 1 )  THEN
@@ -1331 +1311 @@
        IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
           ug_vertical_gradient_level(1) = 0.0_wp
-       ENDIF  
+       ENDIF
 
 !
@@ -1345 +1325 @@
           vg(0) = vg_surface
           DO  k = 1, nzt+1
-             IF ( i < 11 )  THEN 
+             IF ( i < 11 )  THEN
                 IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
                      vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
@@ -1482 +1462 @@
        ENDIF
 
-         
+
     ENDIF
 
@@ -1522 +1502 @@
           CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
         ENDIF
-    ENDIF    
+    ENDIF
 
 !
@@ -1616 +1596 @@
        IF ( use_ug_for_galilei_tr                    .AND.                     &
             ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
-            ug_vertical_gradient(1) == 0.0_wp        .AND.                     & 
+            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
             vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
             vg_vertical_gradient(1) == 0.0_wp )  THEN
@@ -1913 +1893 @@
 
     ENDIF
-    
+
     IF ( passive_scalar )  THEN
 
@@ -1934 +1914 @@
 !
 !--    A fixed scalar concentration at the top implies Dirichlet boundary
-!--    condition for scalar. Hence, in this case specification of a constant 
+!--    condition for scalar. Hence, in this case specification of a constant
 !--    scalar flux is forbidden.
        IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux &
@@ -2652 +2632 @@
              dopr_unit(i)  = TRIM ( heatflux_output_unit )
              hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
-             
+
           CASE ( 'v*pt*' )
              dopr_index(i) = 62
@@ -2976 +2956 @@
              ENDIF
 
- 
+
 
           CASE DEFAULT
@@ -2987 +2967 @@
                                                      unit, dopr_unit(i) )
              ENDIF
-             
+
              IF ( unit == 'illegal' )  THEN
                 unit = ''
@@ -3251 +3231 @@
                                  'res passive_scalar = .TRUE.'
                 CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
-             ENDIF             
+             ENDIF
 
              IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
@@ -3259 +3239 @@
              IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
              IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
-             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'      
+             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
              IF ( TRIM( var ) == 't*'     )  unit = 'K'
              IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
              IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
-             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm' 
-             
+             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
+
           CASE DEFAULT
 
              CALL lsm_check_data_output ( var, unit, i, ilen, k )
- 
+
              IF ( unit == 'illegal' )  THEN
                 CALL radiation_check_data_output( var, unit, i, ilen, k )
@@ -3284 +3264 @@
                  CALL pcm_check_data_output( var, unit )
              ENDIF
-             
+
              IF ( unit == 'illegal' )  THEN
                 unit = ''
@@ -3408 +3388 @@
        ENDIF
     ENDDO
-    
+
     IF ( masks < 0  .OR.  masks > max_masks )  THEN
        WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
@@ -3426 +3406 @@
 !
 !--    Generate masks for masked data output
-!--    Parallel netcdf output is not tested so far for masked data, hence 
+!--    Parallel netcdf output is not tested so far for masked data, hence
 !--    netcdf_data_format is switched back to non-paralell output.
        netcdf_data_format_save = netcdf_data_format
@@ -3436 +3416 @@
                            '6 (parallel netCDF 4 Classic model) '//            &
                            '&are currently not supported (not yet tested) ' // &
-                           'for masked data.&Using respective non-parallel' // &  
+                           'for masked data.&Using respective non-parallel' // &
                            ' output for masked data.'
           CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
@@ -3629 +3609 @@
 !-- near the inflow and the perturbation area is further limited to ...(1)
 !-- after the initial phase of the flow.
-    
+
     IF ( bc_lr /= 'cyclic' )  THEN
        IF ( inflow_disturbance_begin == -1 )  THEN
@@ -3666 +3646 @@
     ENDIF
 
-    IF ( random_generator == 'random-parallel' )  THEN 
+    IF ( random_generator == 'random-parallel' )  THEN
        dist_nxl = nxl;  dist_nxr = nxr
        dist_nys = nys;  dist_nyn = nyn
@@ -3921 +3901 @@
 
     IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
-       message_string = 'The usage of large scale forcing from external &'//   & 
+       message_string = 'The usage of large scale forcing from external &'//   &
                         'file LSF_DATA requires humidity = .T..'
        CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
@@ -3927 +3907 @@
 
     IF ( large_scale_forcing  .AND.  passive_scalar )  THEN
-       message_string = 'The usage of large scale forcing from external &'//   & 
+       message_string = 'The usage of large scale forcing from external &'//   &
                         'file LSF_DATA is not implemented for passive scalars'
        CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
@@ -3934 +3914 @@
     IF ( large_scale_forcing  .AND.  topography /= 'flat'                      &
                               .AND.  .NOT.  lsf_exception )  THEN
-       message_string = 'The usage of large scale forcing from external &'//   & 
+       message_string = 'The usage of large scale forcing from external &'//   &
                         'file LSF_DATA is not implemented for non-flat topography'
        CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
@@ -3940 +3920 @@
 
     IF ( large_scale_forcing  .AND.  ocean )  THEN
-       message_string = 'The usage of large scale forcing from external &'//   & 
+       message_string = 'The usage of large scale forcing from external &'//   &
                         'file LSF_DATA is not implemented for ocean runs'
        CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
@@ -3955 +3935 @@
           do2d_yz_time_count = 0
           domask_time_count  = 0
-       ENDIF 
+       ENDIF
     ENDIF
 
@@ -3989 +3969 @@
 !> Check the length of data output intervals. In case of parallel NetCDF output
 !> the time levels of the output files need to be fixed. Therefore setting the
-!> output interval to 0.0s (usually used to output each timestep) is not 
+!> output interval to 0.0s (usually used to output each timestep) is not
 !> possible as long as a non-fixed timestep is used.
 !------------------------------------------------------------------------------!
@@ -4018 +3998 @@
 
     END SUBROUTINE check_dt_do
-    
-    
-    
-    
+
+
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4034 +4014 @@
 
 
-       IMPLICIT NONE    
-    
+       IMPLICIT NONE
+
        INTEGER(iwp) ::  i     !< counter
        INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
-       
-       REAL(wp)     ::  bc_t_val      !< model top gradient       
-       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity 
+
+       REAL(wp)     ::  bc_t_val      !< model top gradient
+       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity
        REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
 
-       
+
        REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
        REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
        REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
-       
+
        i = 1
        gradient = 0.0_wp
@@ -4109 +4089 @@
 
        ENDIF
-       
+
 !
 !--    In case of no given gradients, choose zero gradient conditions
@@ -4118 +4098 @@
 !--    Store gradient at the top boundary for possible Neumann boundary condition
        bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
-    
+
     END SUBROUTINE init_vertical_profiles
-    
-    
-     
+
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4129 +4109 @@
 !------------------------------------------------------------------------------!
 
-    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t ) 
-
-
-       IMPLICIT NONE    
-    
-       CHARACTER (LEN=1)   ::  sq                       !< 
+    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
+
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=1)   ::  sq                       !<
        CHARACTER (LEN=*)   ::  bc_b
        CHARACTER (LEN=*)   ::  bc_t
@@ -4171 +4151 @@
           CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
        ENDIF
-       
-    
-    END SUBROUTINE set_bc_scalars   
+
+
+    END SUBROUTINE set_bc_scalars
 
 
@@ -4180 +4160 @@
 ! Description:
 ! ------------
-!> Check for consistent settings of bottom boundary conditions for humidity 
+!> Check for consistent settings of bottom boundary conditions for humidity
 !> and scalars.
 !------------------------------------------------------------------------------!
@@ -4189 +4169 @@
 
 
-       IMPLICIT NONE    
-    
-       CHARACTER (LEN=1)   ::  sq                       !< 
+       IMPLICIT NONE
+
+       CHARACTER (LEN=1)   ::  sq                       !<
        CHARACTER (LEN=*)   ::  bc_b
        CHARACTER (LEN=*)   ::  err_nr_1
        CHARACTER (LEN=*)   ::  err_nr_2
-       
+
        INTEGER(iwp)        ::  ibc_b
-       
+
        LOGICAL             ::  constant_flux
-       
+
        REAL(wp)            ::  surface_initial_change
- 
+
 !
 !--    A given surface value implies Dirichlet boundary condition for
@@ -4220 +4200 @@
                  surface_initial_change
           CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
-       ENDIF  
-       
-    
-    END SUBROUTINE check_bc_scalars   
-    
-    
+       ENDIF
+
+
+    END SUBROUTINE check_bc_scalars
+
+
 
  END SUBROUTINE check_parameters

palm/trunk/SOURCE/data_output_ptseries.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! SGS velocities also possible for curvature_solution_effects = .TRUE. 
 !
 ! 2000 2016-08-20 18:09:15Z knoop
@@ -36 +37 @@
 !
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
+! Code annotations made doxygen readable
+!
 ! 1359 2014-04-11 17:15:14Z hoffmann
-! New particle structure integrated. 
-! 
+! New particle structure integrated.
+!
 ! 1353 2014-04-08 15:21:23Z heinze
 ! REAL constants provided with KIND-attribute
-! 
+!
 ! 1327 2014-03-21 11:00:16Z raasch
 ! -netcdf output queries
@@ -75 +76 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE data_output_ptseries
- 
+
 
     USE control_parameters,                                                    &
@@ -97 +98 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  curvature_solution_effects, grid_particles, number_of_particles,&
-               number_of_particle_groups, particles, prt_count
+        ONLY:  grid_particles, number_of_particles, number_of_particle_groups, &
+               particles, prt_count
 
     USE pegrid
@@ -164 +165 @@
                    pts_value_l(0,7)  = pts_value_l(0,7) + particles(n)%speed_y  ! mean v
                    pts_value_l(0,8)  = pts_value_l(0,8) + particles(n)%speed_z  ! mean w
-                   IF ( .NOT. curvature_solution_effects )  THEN
-                      pts_value_l(0,9)  = pts_value_l(0,9)  + particles(n)%rvar1 ! mean sgsu
-                      pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
-                      pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
-                   ENDIF
+                   pts_value_l(0,9)  = pts_value_l(0,9)  + particles(n)%rvar1 ! mean sgsu
+                   pts_value_l(0,10) = pts_value_l(0,10) + particles(n)%rvar2 ! mean sgsv
+                   pts_value_l(0,11) = pts_value_l(0,11) + particles(n)%rvar3 ! mean sgsw
                    IF ( particles(n)%speed_z > 0.0_wp )  THEN
                       pts_value_l(0,12) = pts_value_l(0,12) + 1.0_wp  ! # of upward moving prts
@@ -198 +197 @@
                       pts_value_l(jg,7)  = pts_value_l(jg,7) + particles(n)%speed_y
                       pts_value_l(jg,8)  = pts_value_l(jg,8) + particles(n)%speed_z
-                      IF ( .NOT. curvature_solution_effects )  THEN
-                         pts_value_l(jg,9)  = pts_value_l(jg,9)  + particles(n)%rvar1
-                         pts_value_l(jg,10) = pts_value_l(jg,10) + particles(n)%rvar2
-                         pts_value_l(jg,11) = pts_value_l(jg,11) + particles(n)%rvar3
-                      ENDIF
+                      pts_value_l(jg,9)  = pts_value_l(jg,9)  + particles(n)%rvar1
+                      pts_value_l(jg,10) = pts_value_l(jg,10) + particles(n)%rvar2
+                      pts_value_l(jg,11) = pts_value_l(jg,11) + particles(n)%rvar3
                       IF ( particles(n)%speed_z > 0.0_wp )  THEN
                          pts_value_l(jg,12) = pts_value_l(jg,12) + 1.0_wp
@@ -299 +296 @@
                 pts_value_l(0,25) = pts_value_l(0,25) + ( particles(n)%speed_z - &
                                                           pts_value(0,8) )**2   ! w*2
-                IF ( .NOT. curvature_solution_effects )  THEN
-                   pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
-                                                             pts_value(0,9) )**2   ! u"2
-                   pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
-                                                             pts_value(0,10) )**2  ! v"2
-                   pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
-                                                             pts_value(0,11) )**2  ! w"2
-                ENDIF
+                pts_value_l(0,26) = pts_value_l(0,26) + ( particles(n)%rvar1 - &
+                                                          pts_value(0,9) )**2   ! u"2
+                pts_value_l(0,27) = pts_value_l(0,27) + ( particles(n)%rvar2 - &
+                                                          pts_value(0,10) )**2  ! v"2
+                pts_value_l(0,28) = pts_value_l(0,28) + ( particles(n)%rvar3 - &
+                                                          pts_value(0,11) )**2  ! w"2
 !
 !--             Repeat the same for the respective particle group
@@ -324 +319 @@
                    pts_value_l(jg,25) = pts_value_l(jg,25) + ( particles(n)%speed_z - &
                                                              pts_value(jg,8) )**2
-                   IF ( .NOT. curvature_solution_effects )  THEN
-                      pts_value_l(jg,26) = pts_value_l(jg,26) + ( particles(n)%rvar1 - &
-                                                                pts_value(jg,9) )**2
-                      pts_value_l(jg,27) = pts_value_l(jg,27) + ( particles(n)%rvar2 - &
-                                                                pts_value(jg,10) )**2
-                      pts_value_l(jg,28) = pts_value_l(jg,28) + ( particles(n)%rvar3 - &
-                                                                pts_value(jg,11) )**2
-                   ENDIF
+                   pts_value_l(jg,26) = pts_value_l(jg,26) + ( particles(n)%rvar1 - &
+                                                             pts_value(jg,9) )**2
+                   pts_value_l(jg,27) = pts_value_l(jg,27) + ( particles(n)%rvar2 - &
+                                                             pts_value(jg,10) )**2
+                   pts_value_l(jg,28) = pts_value_l(jg,28) + ( particles(n)%rvar3 - &
+                                                             pts_value(jg,11) )**2
                 ENDIF
 

palm/trunk/SOURCE/lpm_droplet_collision.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Consideration of aerosol mass during collision. Average impact algorithm has
+! been removed.
+!
+! 2274 2017-06-09 13:27:48Z Giersch
 ! Changed error messages
 !
@@ -39 +43 @@
 !
 ! 1860 2016-04-13 13:21:28Z hoffmann
-! Interpolation of dissipation rate adjusted to more reasonable values. 
+! Interpolation of dissipation rate adjusted to more reasonable values.
 !
 ! 1822 2016-04-07 07:49:42Z hoffmann
@@ -49 +53 @@
 !
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
+! Code annotations made doxygen readable
+!
 ! 1359 2014-04-11 17:15:14Z hoffmann
-! New particle structure integrated. 
+! New particle structure integrated.
 ! Kind definition added to all floating point numbers.
 !
@@ -96 +100 @@
 !>              the modification of the relative velocity between the droplets,
 !>              the effect of preferential concentration, and the enhancement of
-!>              collision efficiencies. 
+!>              collision efficiencies.
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_droplet_collision (i,j,k)
- 
-
 
     USE arrays_3d,                                                             &
@@ -126 +128 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  all_or_nothing, average_impact, dissipation_classes,            &
-               hall_kernel, iran_part, number_of_particles, particles,         &
-               particle_type, prt_count, use_kernel_tables, wang_kernel
+        ONLY:  curvature_solution_effects, dissipation_classes, hall_kernel,   &
+               iran_part, number_of_particles, particles, particle_type,       &
+               prt_count, rho_s, use_kernel_tables, wang_kernel
 
     USE random_function_mod,                                                   &
@@ -150 +152 @@
     REAL(wp) ::  epsilon                 !< dissipation rate
     REAL(wp) ::  factor_volume_to_mass   !< 4.0 / 3.0 * pi * rho_l
-    REAL(wp) ::  xm                      !< mean mass of droplet m
-    REAL(wp) ::  xn                      !< mean mass of droplet n
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight  !< weighting factor
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mass    !< total mass of super droplet
+    REAL(wp) ::  xm                      !< droplet mass of super-droplet m
+    REAL(wp) ::  xn                      !< droplet mass of super-droplet n
+    REAL(wp) ::  xsm                     !< aerosol mass of super-droplet m
+    REAL(wp) ::  xsn                     !< aerosol mass of super-droplet n
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight    !< weighting factor
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mass      !< total mass of super droplet
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  aero_mass !< total aerosol mass of super droplet
 
     CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
 
     number_of_particles   = prt_count(k,j,i)
-    factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l 
-    ddV                   = 1 / ( dx * dy * dz )
+    factor_volume_to_mass = 4.0_wp / 3.0_wp * pi * rho_l
+    ddV                   = 1.0_wp / ( dx * dy * dz )
 !
 !-- Collision requires at least one super droplet inside the box
     IF ( number_of_particles > 0 )  THEN
 
-!
-!--    Now apply the different kernels
        IF ( use_kernel_tables )  THEN
 !
 !--       Fast method with pre-calculated collection kernels for
 !--       discrete radius- and dissipation-classes.
-!--
-!--       Determine dissipation class index of this gridbox
           IF ( wang_kernel )  THEN
              eclass = INT( diss(k,j,i) * 1.0E4_wp / 600.0_wp * &
@@ -180 +181 @@
              epsilon = 0.0_wp
           ENDIF
+
           IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
              eclass = 0   ! Hall kernel is used
@@ -186 +188 @@
           ENDIF
 
-!
-!--       Droplet collision are calculated using collision-coalescence
-!--       formulation proposed by Wang (see PALM documentation)
-!--       Temporary fields for total mass of super-droplet and weighting factors 
-!--       are allocated.
-          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
-
-          mass(1:number_of_particles)   = particles(1:number_of_particles)%weight_factor * &
-                                          particles(1:number_of_particles)%radius**3     * &
-                                          factor_volume_to_mass
-          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
-
-          IF ( average_impact )  THEN  ! select collision algorithm
-
-             DO  n = 1, number_of_particles
-
-                rclass_l = particles(n)%class
-                xn       = mass(n) / weight(n)
-
-                DO  m = n, number_of_particles
-
-                   rclass_s = particles(m)%class
-                   xm       = mass(m) / weight(m)
-
-                   IF ( xm .LT. xn )  THEN
-                      
-!
-!--                   Particle n collects smaller particle m
-                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
-                                               weight(n) * ddV * dt_3d
-
-                      mass(n)   = mass(n)   + mass(m)   * collection_probability
-                      weight(m) = weight(m) - weight(m) * collection_probability
-                      mass(m)   = mass(m)   - mass(m)   * collection_probability
-                   ELSEIF ( xm .GT. xn )  THEN 
-!
-!--                   Particle m collects smaller particle n
-                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
-                                               weight(m) * ddV * dt_3d
-
-                      mass(m)   = mass(m)   + mass(n)   * collection_probability
-                      weight(n) = weight(n) - weight(n) * collection_probability
-                      mass(n)   = mass(n)   - mass(n)   * collection_probability
-                   ELSE
-!
-!--                   Same-size collections. If n = m, weight is reduced by the
-!--                   number of possible same-size collections; the total mass 
-!--                   is not changed during same-size collection.
-!--                   Same-size collections of different 
-!--                   particles ( n /= m ) are treated as same-size collections 
-!--                   of ONE partilce with weight = weight(n) + weight(m) and 
-!--                   mass = mass(n) + mass(m).
-!--                   Accordingly, each particle loses the same number of 
-!--                   droplets to the other particle, but this has no effect on 
-!--                   total mass mass, since the exchanged droplets have the 
-!--                   same radius. 
-
-!--                   Note: For m = n this equation is an approximation only 
-!--                   valid for weight >> 1 (which is usually the case). The
-!--                   approximation is weight(n)-1 = weight(n).
-                      weight(n) = weight(n) - 0.5_wp * weight(n) *                &
-                                              ckernel(rclass_l,rclass_s,eclass) * &
-                                              weight(m) * ddV * dt_3d
-                      IF ( n .NE. m )  THEN 
-                         weight(m) = weight(m) - 0.5_wp * weight(m) *                &
-                                                 ckernel(rclass_l,rclass_s,eclass) * &
-                                                 weight(n) * ddV * dt_3d
-                      ENDIF
-                   ENDIF
-
-                ENDDO
-
-                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
-
-             ENDDO
-
-          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
-
-             DO  n = 1, number_of_particles
-
-                rclass_l = particles(n)%class
-                xn       = mass(n) / weight(n) ! mean mass of droplet n
-
-                DO  m = n, number_of_particles
-
-                   rclass_s = particles(m)%class
-                   xm = mass(m) / weight(m) ! mean mass of droplet m
-
-                   IF ( weight(n) .LT. weight(m) )  THEN
-!
-!--                   Particle n collects weight(n) droplets of particle m  
-                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
-                                               weight(m) * ddV * dt_3d
-
-                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
-                         mass(n)   = mass(n)   + weight(n) * xm
-                         weight(m) = weight(m) - weight(n)
-                         mass(m)   = mass(m)   - weight(n) * xm
-                      ENDIF
-
-                   ELSEIF ( weight(m) .LT. weight(n) )  THEN 
-!
-!--                   Particle m collects weight(m) droplets of particle n  
-                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
-                                               weight(n) * ddV * dt_3d
-
-                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
-                         mass(m)   = mass(m)   + weight(m) * xn
-                         weight(n) = weight(n) - weight(m)
-                         mass(n)   = mass(n)   - weight(m) * xn
-                      ENDIF
-                   ELSE
-!
-!--                   Collisions of particles of the same weighting factor. 
-!--                   Particle n collects 1/2 weight(n) droplets of particle m,
-!--                   particle m collects 1/2 weight(m) droplets of particle n.
-!--                   The total mass mass changes accordingly.
-!--                   If n = m, the first half of the droplets coalesces with the 
-!--                   second half of the droplets; mass is unchanged because
-!--                   xm = xn for n = m.
-
-!--                   Note: For m = n this equation is an approximation only 
-!--                   valid for weight >> 1 (which is usually the case). The
-!--                   approximation is weight(n)-1 = weight(n).
-                      collection_probability = ckernel(rclass_l,rclass_s,eclass) * &
-                                               weight(n) * ddV * dt_3d
-
-                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
-                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
-                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
-                         weight(n) = weight(n) - 0.5_wp * weight(m)
-                         weight(m) = weight(n)
-                      ENDIF
-                   ENDIF
-
-                ENDDO
-
-                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
-
-             ENDDO
-
-          ENDIF
-
-
-
-
-          IF ( ANY(weight < 0.0_wp) )  THEN
-                WRITE( message_string, * ) 'negative weighting factor'
-                CALL message( 'lpm_droplet_collision', 'PA0028',      &
-                               2, 2, -1, 6, 1 )
-          ENDIF
-
-          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
-                                                      ( weight(1:number_of_particles) &
-                                                        * factor_volume_to_mass       &
-                                                      )                               &
-                                                    )**0.33333333333333_wp
-
-          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
-
-          DEALLOCATE(weight, mass)
-
-       ELSEIF ( .NOT. use_kernel_tables )  THEN
-!
-!--       Collection kernels are calculated for every new
+       ELSE
+!
+!--       Collection kernels are re-calculated for every new
 !--       grid box. First, allocate memory for kernel table.
 !--       Third dimension is 1, because table is re-calculated for
@@ -359 +199 @@
 !--       the kernel is based on the previous time step
           CALL recalculate_kernel( i, j, k )
-!
-!--       Droplet collision are calculated using collision-coalescence
-!--       formulation proposed by Wang (see PALM documentation)
-!--       Temporary fields for total mass of super-droplet and weighting factors 
-!--       are allocated.
-          ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
-
-          mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
-                                        particles(1:number_of_particles)%radius**3     * &
-                                        factor_volume_to_mass
-
-          weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
-
-          IF ( average_impact )  THEN  ! select collision algorithm
-
-             DO  n = 1, number_of_particles
-
-                xn = mass(n) / weight(n) ! mean mass of droplet n
-
-                DO  m = n, number_of_particles
-
-                   xm = mass(m) / weight(m) !mean mass of droplet m
-
-                   IF ( xm .LT. xn )  THEN
-!
-!--                   Particle n collects smaller particle m 
-                      collection_probability = ckernel(n,m,1) * weight(n) *    &
-                                               ddV * dt_3d
-
-                      mass(n)   = mass(n)   + mass(m)   * collection_probability
-                      weight(m) = weight(m) - weight(m) * collection_probability
-                      mass(m)   = mass(m)   - mass(m)   * collection_probability
-                   ELSEIF ( xm .GT. xn )  THEN 
-!
-!--                   Particle m collects smaller particle n
-                      collection_probability = ckernel(n,m,1) * weight(m) *    &
-                                               ddV * dt_3d
-
-                      mass(m)   = mass(m)   + mass(n)   * collection_probability
-                      weight(n) = weight(n) - weight(n) * collection_probability
-                      mass(n)   = mass(n)   - mass(n)   * collection_probability
-                   ELSE
-!
-!--                   Same-size collections. If n = m, weight is reduced by the
-!--                   number of possible same-size collections; the total mass 
-!--                   mass is not changed during same-size collection. 
-!--                   Same-size collections of different 
-!--                   particles ( n /= m ) are treated as same-size collections 
-!--                   of ONE partilce with weight = weight(n) + weight(m) and 
-!--                   mass = mass(n) + mass(m).
-!--                   Accordingly, each particle loses the same number of 
-!--                   droplets to the other particle, but this has no effect on 
-!--                   total mass mass, since the exchanged droplets have the 
-!--                   same radius. 
+
+       ENDIF
+!
+!--    Temporary fields for total mass of super-droplet, aerosol mass, and
+!--    weighting factor are allocated.
+       ALLOCATE(mass(1:number_of_particles), weight(1:number_of_particles))
+       IF ( curvature_solution_effects )  ALLOCATE(aero_mass(1:number_of_particles))
+
+       mass(1:number_of_particles)   = particles(1:number_of_particles)%weight_factor * &
+                                       particles(1:number_of_particles)%radius**3     * &
+                                       factor_volume_to_mass
+
+       weight(1:number_of_particles) = particles(1:number_of_particles)%weight_factor
+
+       IF ( curvature_solution_effects )  THEN
+          aero_mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
+                                             particles(1:number_of_particles)%aux1**3       * &
+                                             4.0 / 3.0 * pi * rho_s
+       ENDIF
+!
+!--    Calculate collision/coalescence
+       DO  n = 1, number_of_particles
+
+          DO  m = n, number_of_particles
+!
+!--          For collisions, the weighting factor of at least one super-droplet
+!--          needs to be larger or equal to one.
+             IF ( MIN( weight(n), weight(m) ) .LT. 1.0 )  CYCLE
+!
+!--          Get mass of individual droplets (aerosols)
+             xn = mass(n) / weight(n)
+             xm = mass(m) / weight(m)
+             IF ( curvature_solution_effects )  THEN
+                xsn = aero_mass(n) / weight(n)
+                xsm = aero_mass(m) / weight(m)
+             ENDIF
+!
+!--          Probability that the necessary collisions take place
+             IF ( use_kernel_tables )  THEN
+                rclass_l = particles(n)%class
+                rclass_s = particles(m)%class
+
+                collection_probability  = MAX( weight(n), weight(m) ) *     &
+                                          ckernel(rclass_l,rclass_s,eclass) * ddV * dt_3d
+             ELSE
+                collection_probability  = MAX( weight(n), weight(m) ) *     &
+                                          ckernel(n,m,1) * ddV * dt_3d
+             ENDIF
+!
+!--          Calculate the number of collections and consider multiple collections.
+!--          (Accordingly, p_crit will be 0.0, 1.0, 2.0, ...)
+             IF ( collection_probability - FLOOR(collection_probability)    &
+                  .GT. random_function( iran_part ) )  THEN
+                collection_probability = FLOOR(collection_probability) + 1.0_wp
+             ELSE
+                collection_probability = FLOOR(collection_probability)
+             ENDIF
+
+             IF ( collection_probability .GT. 0.0 )  THEN
+!
+!--             Super-droplet n collects droplets of super-droplet m
+                IF ( weight(n) .LT. weight(m) )  THEN
+
+                   mass(n)   = mass(n)   + weight(n) * xm * collection_probability
+                   weight(m) = weight(m) - weight(n)      * collection_probability
+                   mass(m)   = mass(m)   - weight(n) * xm * collection_probability
+                   IF ( curvature_solution_effects )  THEN
+                      aero_mass(n) = aero_mass(n) + weight(n) * xsm * collection_probability
+                      aero_mass(m) = aero_mass(m) - weight(n) * xsm * collection_probability
+                   ENDIF
+
+                ELSEIF ( weight(m) .LT. weight(n) )  THEN
+
+                   mass(m)   = mass(m)   + weight(m) * xn * collection_probability
+                   weight(n) = weight(n) - weight(m)      * collection_probability
+                   mass(n)   = mass(n)   - weight(m) * xn * collection_probability
+                   IF ( curvature_solution_effects )  THEN
+                      aero_mass(m) = aero_mass(m) + weight(m) * xsn * collection_probability
+                      aero_mass(n) = aero_mass(n) - weight(m) * xsn * collection_probability
+                   ENDIF
+
+                ELSE
+!
+!--                Collisions of particles of the same weighting factor.
+!--                Particle n collects 1/2 weight(n) droplets of particle m,
+!--                particle m collects 1/2 weight(m) droplets of particle n.
+!--                The total mass mass changes accordingly.
+!--                If n = m, the first half of the droplets coalesces with the
+!--                second half of the droplets; mass is unchanged because
+!--                xm = xn for n = m.
 !--
-!--                   Note: For m = n this equation is an approximation only 
-!--                   valid for weight >> 1 (which is usually the case). The
-!--                   approximation is weight(n)-1 = weight(n).
-                      weight(n) = weight(n) - 0.5_wp * weight(n) *             &
-                                              ckernel(n,m,1) * weight(m) *     &
-                                              ddV * dt_3d
-                      IF ( n .NE. m )  THEN
-                         weight(m) = weight(m) - 0.5_wp * weight(m) *          &
-                                                 ckernel(n,m,1) * weight(n) *  &
-                                                 ddV * dt_3d
-                      ENDIF
+!--                Note: For m = n this equation is an approximation only
+!--                valid for weight >> 1 (which is usually the case). The
+!--                approximation is weight(n)-1 = weight(n).
+                   mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
+                   mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
+                   IF ( curvature_solution_effects )  THEN
+                      aero_mass(n) = aero_mass(n) + 0.5_wp * weight(n) * ( xsm - xsn )
+                      aero_mass(m) = aero_mass(m) + 0.5_wp * weight(m) * ( xsn - xsm )
                    ENDIF
-
-
-                ENDDO
-
-                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
-
-             ENDDO
-
-          ELSEIF ( all_or_nothing )  THEN  ! select collision algorithm
-
-             DO  n = 1, number_of_particles
-
-                xn = mass(n) / weight(n) ! mean mass of droplet n
-
-                DO  m = n, number_of_particles
-
-                   xm = mass(m) / weight(m) !mean mass of droplet m
-
-                   IF ( weight(n) .LT. weight(m) )  THEN
-!
-!--                   Particle n collects smaller particle m 
-                      collection_probability = ckernel(n,m,1) * weight(m) *    &
-                                               ddV * dt_3d
-
-                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
-                         mass(n)   = mass(n)   + weight(n) * xm  
-                         weight(m) = weight(m) - weight(n)
-                         mass(m)   = mass(m)   - weight(n) * xm
-                      ENDIF
-
-                   ELSEIF ( weight(m) .LT. weight(n) )  THEN
-!
-!--                   Particle m collects smaller particle n
-                      collection_probability = ckernel(n,m,1) * weight(n) *    &
-                                               ddV * dt_3d
-
-                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
-                         mass(m)   = mass(m)   + weight(m) * xn
-                         weight(n) = weight(n) - weight(m)
-                         mass(n)   = mass(n)   - weight(m) * xn
-                      ENDIF
-                   ELSE
-!
-!--                   Collisions of particles of the same weighting factor. 
-!--                   Particle n collects 1/2 weight(n) droplets of particle m,
-!--                   particle m collects 1/2 weight(m) droplets of particle n.
-!--                   The total mass mass changes accordingly.
-!--                   If n = m, the first half of the droplets coalesces with the 
-!--                   second half of the droplets; mass is unchanged because
-!--                   xm = xn for n = m.
-!--
-!--                   Note: For m = n this equation is an approximation only 
-!--                   valid for weight >> 1 (which is usually the case). The
-!--                   approximation is weight(n)-1 = weight(n).
-                      collection_probability = ckernel(n,m,1) * weight(n) *    &
-                                               ddV * dt_3d
-
-                      IF ( collection_probability .GT. random_function( iran_part ) )  THEN
-                         mass(n)   = mass(n)   + 0.5_wp * weight(n) * ( xm - xn )
-                         mass(m)   = mass(m)   + 0.5_wp * weight(m) * ( xn - xm )
-                         weight(n) = weight(n) - 0.5_wp * weight(m)
-                         weight(m) = weight(n)
-                      ENDIF
-                   ENDIF
-
-
-                ENDDO
-
-                ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
-
-             ENDDO
-
-          ENDIF
-
-          IF ( ANY(weight < 0.0_wp) )  THEN
-                WRITE( message_string, * ) 'negative weighting factor'
-                CALL message( 'lpm_droplet_collision', 'PA0028',      &
-                               2, 2, -1, 6, 1 )
-          ENDIF
-
-          particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
-                                                      ( weight(1:number_of_particles) &
-                                                        * factor_volume_to_mass       &
-                                                      )                               &
-                                                    )**0.33333333333333_wp
-
-          particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
-
-          DEALLOCATE( weight, mass, ckernel )
-
-       ENDIF
- 
+                   weight(n) = weight(n) - 0.5_wp * weight(m)
+                   weight(m) = weight(n)
+
+                ENDIF
+
+             ENDIF
+
+          ENDDO
+
+          ql_vp(k,j,i) = ql_vp(k,j,i) + mass(n) / factor_volume_to_mass
+
+       ENDDO
+
+       IF ( ANY(weight < 0.0_wp) )  THEN
+             WRITE( message_string, * ) 'negative weighting factor'
+             CALL message( 'lpm_droplet_collision', 'PA0028',      &
+                            2, 2, -1, 6, 1 )
+       ENDIF
+
+       particles(1:number_of_particles)%radius = ( mass(1:number_of_particles) /   &
+                                                   ( weight(1:number_of_particles) &
+                                                     * factor_volume_to_mass       &
+                                                   )                               &
+                                                 )**0.33333333333333_wp
+
+       IF ( curvature_solution_effects )  THEN
+          particles(1:number_of_particles)%aux1 = ( aero_mass(1:number_of_particles) / &
+                                                    ( weight(1:number_of_particles)    &
+                                                      * 4.0_wp / 3.0_wp * pi * rho_s   &
+                                                    )                                  &
+                                                  )**0.33333333333333_wp
+       ENDIF
+
+       particles(1:number_of_particles)%weight_factor = weight(1:number_of_particles)
+
+       DEALLOCATE( weight, mass )
+       IF ( curvature_solution_effects )  DEALLOCATE( aero_mass )
+       IF ( .NOT. use_kernel_tables )  DEALLOCATE( ckernel )
+
 !
 !--    Check if LWC is conserved during collision process
@@ -532 +352 @@
 
     ENDIF
- 
+
     CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' )
 

palm/trunk/SOURCE/lpm_droplet_condensation.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Rosenbrock scheme improved. Gas-kinetic effect added.
 !
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 1890 2016-04-22 08:52:11Z hoffmann
 ! Some improvements of the Rosenbrock method. If the Rosenbrock method needs more
 ! than 40 iterations to find a sufficient time setp, the model is not aborted.
-! This might lead to small erros. But they can be assumend as negligible, since 
-! the maximum timesetp of the Rosenbrock method after 40 iterations will be 
-! smaller than 10^-16 s.  
-! 
+! This might lead to small erros. But they can be assumend as negligible, since
+! the maximum timesetp of the Rosenbrock method after 40 iterations will be
+! smaller than 10^-16 s.
+!
 ! 1871 2016-04-15 11:46:09Z hoffmann
 ! Initialization of aerosols added.
@@ -53 +54 @@
 !
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
+! Code annotations made doxygen readable
+!
 ! 1359 2014-04-11 17:15:14Z hoffmann
-! New particle structure integrated. 
+! New particle structure integrated.
 ! Kind definition added to all floating point numbers.
 !
 ! 1346 2014-03-27 13:18:20Z heinze
-! Bugfix: REAL constants provided with KIND-attribute especially in call of 
+! Bugfix: REAL constants provided with KIND-attribute especially in call of
 ! intrinsic function like MAX, MIN, SIGN
 !
@@ -108 +109 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_droplet_condensation (ip,jp,kp)
- 
+
 
     USE arrays_3d,                                                             &
-        ONLY:  hyp, pt, q,  ql_c, ql_v
+        ONLY:  hyp, pt, q, ql_c, ql_v
 
     USE cloud_parameters,                                                      &
@@ -139 +140 @@
                use_kernel_tables, vanthoff, wang_kernel
 
-
     IMPLICIT NONE
 
     INTEGER(iwp) :: ip                         !<
-    INTEGER(iwp) :: internal_timestep_count    !<
     INTEGER(iwp) :: jp                         !<
-    INTEGER(iwp) :: jtry                       !<
     INTEGER(iwp) :: kp                         !<
     INTEGER(iwp) :: n                          !<
-    INTEGER(iwp) :: ros_count                  !<
- 
-    INTEGER(iwp), PARAMETER ::  maxtry = 40    !<
-
-    LOGICAL ::  repeat                         !<
-
-    REAL(wp) ::  aa                            !<
+
     REAL(wp) ::  afactor                       !< curvature effects
     REAL(wp) ::  arg                           !<
@@ -161 +153 @@
     REAL(wp) ::  delta_r                       !<
     REAL(wp) ::  diameter                      !< diameter of cloud droplets
-    REAL(wp) ::  diff_coeff_v                  !< diffusivity for water vapor
+    REAL(wp) ::  diff_coeff                    !< diffusivity for water vapor
     REAL(wp) ::  drdt                          !<
-    REAL(wp) ::  drdt_ini                      !<
     REAL(wp) ::  dt_ros                        !<
-    REAL(wp) ::  dt_ros_next                   !<
     REAL(wp) ::  dt_ros_sum                    !<
-    REAL(wp) ::  dt_ros_sum_ini                !<
     REAL(wp) ::  d2rdtdr                       !<
-    REAL(wp) ::  errmax                        !<
     REAL(wp) ::  e_a                           !< current vapor pressure
     REAL(wp) ::  e_s                           !< current saturation vapor pressure
-    REAL(wp) ::  err_ros                       !<
-    REAL(wp) ::  g1                            !<
-    REAL(wp) ::  g2                            !<
-    REAL(wp) ::  g3                            !<
-    REAL(wp) ::  g4                            !<
-    REAL(wp) ::  r_ros                         !<
-    REAL(wp) ::  r_ros_ini                     !<
-    REAL(wp) ::  sigma                         !<
-    REAL(wp) ::  thermal_conductivity_v        !< thermal conductivity for water
+    REAL(wp) ::  error                         !< local truncation error in Rosenbrock
+    REAL(wp) ::  k1                            !<
+    REAL(wp) ::  k2                            !<
+    REAL(wp) ::  r_err                         !< First order estimate of Rosenbrock radius
+    REAL(wp) ::  r_ros                         !< Rosenbrock radius
+    REAL(wp) ::  r_ros_ini                     !< initial Rosenbrock radius
+    REAL(wp) ::  r0                            !< gas-kinetic lengthscale
+    REAL(wp) ::  sigma                         !< surface tension of water
+    REAL(wp) ::  thermal_conductivity          !< thermal conductivity for water
     REAL(wp) ::  t_int                         !< temperature
     REAL(wp) ::  w_s                           !< terminal velocity of droplets
-    REAL(wp) ::  re_p                          !<
-
-!
-!-- Parameters for Rosenbrock method
-    REAL(wp), PARAMETER ::  a21 = 2.0_wp               !<
-    REAL(wp), PARAMETER ::  a31 = 48.0_wp / 25.0_wp    !<
-    REAL(wp), PARAMETER ::  a32 = 6.0_wp / 25.0_wp     !<
-    REAL(wp), PARAMETER ::  b1 = 19.0_wp / 9.0_wp      !<
-    REAL(wp), PARAMETER ::  b2 = 0.5_wp                !<
-    REAL(wp), PARAMETER ::  b3 = 25.0_wp / 108.0_wp    !<
-    REAL(wp), PARAMETER ::  b4 = 125.0_wp / 108.0_wp   !<
-    REAL(wp), PARAMETER ::  c21 = -8.0_wp              !<
-    REAL(wp), PARAMETER ::  c31 = 372.0_wp / 25.0_wp   !<
-    REAL(wp), PARAMETER ::  c32 = 12.0_wp / 5.0_wp     !<
-    REAL(wp), PARAMETER ::  c41 = -112.0_wp / 125.0_wp !<
-    REAL(wp), PARAMETER ::  c42 = -54.0_wp / 125.0_wp  !<
-    REAL(wp), PARAMETER ::  c43 = -2.0_wp / 5.0_wp     !<
-    REAL(wp), PARAMETER ::  errcon = 0.1296_wp         !<
-    REAL(wp), PARAMETER ::  e1 = 17.0_wp / 54.0_wp     !<
-    REAL(wp), PARAMETER ::  e2 = 7.0_wp / 36.0_wp      !<
-    REAL(wp), PARAMETER ::  e3 = 0.0_wp                !<
-    REAL(wp), PARAMETER ::  e4 = 125.0_wp / 108.0_wp   !<
-    REAL(wp), PARAMETER ::  eps_ros = 1.0E-3_wp        !< accuracy of Rosenbrock method
-    REAL(wp), PARAMETER ::  gam = 0.5_wp               !<
-    REAL(wp), PARAMETER ::  grow = 1.5_wp              !<
-    REAL(wp), PARAMETER ::  pgrow = -0.25_wp           !<
-    REAL(wp), PARAMETER ::  pshrnk = -1.0_wp /3.0_wp   !<
-    REAL(wp), PARAMETER ::  shrnk = 0.5_wp             !<
-
+    REAL(wp) ::  re_p                          !< particle Reynolds number
+!
+!-- Parameters for Rosenbrock method (see Verwer et al., 1999)
+    REAL(wp), PARAMETER :: prec = 1.0E-3_wp     !< precision of Rosenbrock solution
+    REAL(wp), PARAMETER :: q_increase = 1.5_wp  !< increase factor in timestep
+    REAL(wp), PARAMETER :: q_decrease = 0.9_wp  !< decrease factor in timestep
+    REAL(wp), PARAMETER :: gamma = 0.292893218814_wp !< = 1.0 - 1.0 / SQRT(2.0)
 !
 !-- Parameters for terminal velocity
@@ -224 +190 @@
     REAL(wp), DIMENSION(number_of_particles) ::  new_r                  !<
 
-
-
     CALL cpu_log( log_point_s(42), 'lpm_droplet_condens', 'start' )
 
 !
-!-- Calculate temperature, saturation vapor pressure and current vapor pressure
+!-- Absolute temperature
     t_int = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
-    e_s   = 611.0_wp * EXP( l_d_rv * ( 3.6609E-3_wp - 1.0_wp / t_int ) )
-    e_a   = q(kp,jp,ip) * hyp(kp) / ( 0.378_wp * q(kp,jp,ip) + 0.622_wp )
-!
-!-- Thermal conductivity for water (from Rogers and Yau, Table 7.1),
-!-- diffusivity for water vapor (after Hall und Pruppacher, 1976)
-    thermal_conductivity_v = 7.94048E-05_wp * t_int + 0.00227011_wp
-    diff_coeff_v           = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * &
-                             ( 101325.0_wp / hyp(kp) )
-!
-!-- Calculate effects of heat conductivity and diffusion of water vapor on the 
-!-- condensation/evaporation process (typically known as 1.0 / (F_k + F_d) )
-    ddenom  = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff_v ) +        &
+!
+!-- Saturation vapor pressure (Eq. 10 in Bolton, 1980)
+    e_s = 611.2_wp * EXP( 17.62_wp * ( t_int - 273.15_wp ) / ( t_int - 29.65_wp ) )
+!
+!-- Current vapor pressure
+    e_a = q(kp,jp,ip) * hyp(kp) / ( q(kp,jp,ip) + 0.622_wp )
+!
+!-- Thermal conductivity for water (from Rogers and Yau, Table 7.1)
+    thermal_conductivity = 7.94048E-05_wp * t_int + 0.00227011_wp
+!
+!-- Moldecular diffusivity of water vapor in air (Hall und Pruppacher, 1976)
+    diff_coeff           = 0.211E-4_wp * ( t_int / 273.15_wp )**1.94_wp * &
+                           ( 101325.0_wp / hyp(kp) )
+!
+!-- Lengthscale for gas-kinetic effects (from Mordy, 1959, p. 23):
+    r0 = diff_coeff / 0.036_wp * SQRT( 2.0_wp * pi / ( r_v * t_int ) )
+!
+!-- Calculate effects of heat conductivity and diffusion of water vapor on the
+!-- diffusional growth process (usually known as 1.0 / (F_k + F_d) )
+    ddenom  = 1.0_wp / ( rho_l * r_v * t_int / ( e_s * diff_coeff ) +          &
                          ( l_v / ( r_v * t_int ) - 1.0_wp ) * rho_l *          &
-                         l_v / ( thermal_conductivity_v * t_int )              &
+                         l_v / ( thermal_conductivity * t_int )                &
                        )
-
     new_r = 0.0_wp
-
 !
 !-- Determine ventilation effect on evaporation of large drops
@@ -264 +234 @@
           ENDIF
 !
-!--       First calculate droplet's Reynolds number
+!--       Calculate droplet's Reynolds number
           re_p = 2.0_wp * particles(n)%radius * w_s / molecular_viscosity
 !
@@ -275 +245 @@
        ELSE
 !
-!--       For small droplets or in supersaturated environments, the ventilation 
+!--       For small droplets or in supersaturated environments, the ventilation
 !--       effect does not play a role
           ventilation_effect(n) = 1.0_wp
@@ -281 +251 @@
     ENDDO
 
-!
-!-- Use analytic model for condensational growth
     IF( .NOT. curvature_solution_effects ) then
+!
+!--    Use analytic model for diffusional growth including gas-kinetic
+!--    effects (Mordy, 1959) but without the impact of aerosols.
        DO  n = 1, number_of_particles
-          arg = particles(n)%radius**2 + 2.0_wp * dt_3d * ddenom *             &
-                                         ventilation_effect(n) *               &
-                                         ( e_a / e_s - 1.0_wp )
-          arg = MAX( arg, 1.0E-16_wp )
-          new_r(n) = SQRT( arg )
+          arg      = ( particles(n)%radius + r0 )**2 + 2.0_wp * dt_3d * ddenom * &
+                                                       ventilation_effect(n) *   &
+                                                       ( e_a / e_s - 1.0_wp )
+          arg      = MAX( arg, ( 0.01E-6 + r0 )**2 )
+          new_r(n) = SQRT( arg ) - r0
        ENDDO
-    ENDIF
-
-!
-!-- If selected, use numerical solution of the condensational growth
-!-- equation (e.g., for studying the activation of aerosols).
-!-- Curvature and solutions effects are included in growth equation.
-!-- Change in Radius is calculated with a 4th-order Rosenbrock method
-!-- for stiff o.d.e's with monitoring local truncation error to adjust
-!-- stepsize (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
-    DO  n = 1, number_of_particles
-       IF ( curvature_solution_effects )  THEN
-
-          ros_count = 0
-          repeat = .TRUE.
-!
-!--       Carry out the Rosenbrock algorithm. In case of unreasonable results
-!--       the switch "repeat" will be set true and the algorithm will be carried
-!--       out again with the internal time step set to its initial (small) value.
-!--       Unreasonable results may occur if the external conditions, especially 
-!--       the supersaturation, has significantly changed compared to the last 
-!--       PALM timestep.
-          DO WHILE ( repeat )
-
-             repeat = .FALSE.
-!
-!--          Curvature effect (afactor) with surface tension parameterization
-!--          by Straka (2009)
-             sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
-             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
-!
-!--          Solute effect (bfactor)
-             bfactor = vanthoff * rho_s * particles(n)%rvar2**3 *              &
-                       molecular_weight_of_water /                             &
-                       ( rho_l * molecular_weight_of_solute )
-
-             r_ros = particles(n)%radius
-             dt_ros_sum  = 0.0_wp      ! internal integrated time (s)
-             internal_timestep_count = 0
-!
-!--          Take internal time step values from the end of last PALM time step
-             dt_ros_next = particles(n)%rvar1
-
-!
-!--          Internal time step should not be > 1.0E-2 and < 1.0E-6 
-             IF ( dt_ros_next > 1.0E-2_wp )  THEN
-                dt_ros_next = 1.0E-2_wp
-             ELSEIF ( dt_ros_next < 1.0E-6_wp )  THEN
-                dt_ros_next = 1.0E-6_wp
-             ENDIF
-
-!
-!--          If calculation of Rosenbrock method is repeated due to unreasonalble
-!--          results during previous try the initial internal time step has to be
-!--          reduced
-             IF ( ros_count > 1 )  THEN
-                dt_ros_next = dt_ros_next * 0.1_wp
-             ELSEIF ( ros_count > 5 )  THEN
-!
-!--             Prevent creation of infinite loop
-                message_string = 'ros_count > 5 in Rosenbrock method'
-                CALL message( 'lpm_droplet_condensation', 'PA0018', 2, 2, &
-                               0, 6, 0 )
-             ENDIF
-
-!
-!--          Internal time step must not be larger than PALM time step
-             dt_ros = MIN( dt_ros_next, dt_3d )
-
-!
-!--          Integrate growth equation in time unless PALM time step is reached
-             DO WHILE ( dt_ros_sum < dt_3d )
-
-                internal_timestep_count = internal_timestep_count + 1
-
-!
-!--             Derivative at starting value
-                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
+
+    ELSE
+!
+!--    Integrate the diffusional growth including gas-kinetic (Mordy, 1959),
+!--    as well as curvature and solute effects (e.g., Köhler, 1936).
+!
+!--    Curvature effect (afactor) with surface tension (sigma) by Straka (2009)
+       sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
+!
+!--    Solute effect (afactor)
+       afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
+
+       DO  n = 1, number_of_particles
+!
+!--       Solute effect (bfactor)
+          bfactor = vanthoff * rho_s * particles(n)%aux1**3 *                    &
+                    molecular_weight_of_water / ( rho_l * molecular_weight_of_solute )
+
+          dt_ros     = particles(n)%aux2  ! use previously stored Rosenbrock timestep
+          dt_ros_sum = 0.0_wp
+
+          r_ros     = particles(n)%radius  ! initialize Rosenbrock particle radius
+          r_ros_ini = r_ros
+!
+!--       Integrate growth equation using a 2nd-order Rosenbrock method
+!--       (see Verwer et al., 1999, Eq. (3.2)). The Rosenbrock method adjusts
+!--       its with internal timestep to minimize the local truncation error.
+          DO WHILE ( dt_ros_sum < dt_3d )
+
+             dt_ros = MIN( dt_ros, dt_3d - dt_ros_sum )
+
+             DO
+
+                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -    &
                                                           afactor / r_ros +    &
                                                           bfactor / r_ros**3   &
-                                                        ) / r_ros
-
-                drdt_ini       = drdt
-                dt_ros_sum_ini = dt_ros_sum
-                r_ros_ini      = r_ros
-
-!
-!--             Calculate radial derivative of dr/dt
-                d2rdtdr = ddenom * ventilation_effect(n) *                     &
-                                       ( ( 1.0_wp - e_a / e_s ) / r_ros**2 +   &
-                                         2.0_wp * afactor / r_ros**3 -         &
-                                         4.0_wp * bfactor / r_ros**5           &
-                                       )
-!
-!--             Adjust stepsize unless required accuracy is reached
-                DO  jtry = 1, maxtry+1
-
-                   IF ( jtry == maxtry+1 )  THEN
-                      message_string = 'maxtry > 40 in Rosenbrock method'
-                      CALL message( 'lpm_droplet_condensation', 'PA0347', 0,   &
-                                    1, 0, 6, 0 )
-                   ENDIF
-
-                   aa    = 1.0_wp / ( gam * dt_ros ) - d2rdtdr
-                   g1    = drdt_ini / aa
-                   r_ros = r_ros_ini + a21 * g1
-                   drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
-                                                              afactor / r_ros +    &
-                                                              bfactor / r_ros**3   &
-                                                            ) / r_ros
-
-                   g2    = ( drdt + c21 * g1 / dt_ros )&
-                           / aa
-                   r_ros = r_ros_ini + a31 * g1 + a32 * g2
-                   drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
-                                                              afactor / r_ros +    &
-                                                              bfactor / r_ros**3   &
-                                                            ) / r_ros
-
-                   g3    = ( drdt +  &
-                             ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
-                   g4    = ( drdt +  &
-                             ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
-                   r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
-
-                   dt_ros_sum = dt_ros_sum_ini + dt_ros
-
-                   IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
-                      message_string = 'zero stepsize in Rosenbrock method'
-                      CALL message( 'lpm_droplet_condensation', 'PA0348', 2,   &
-                                    2, 0, 6, 0 )
-                   ENDIF
-!
-!--                Calculate error
-                   err_ros = e1 * g1 + e2 * g2 + e3 * g3 + e4 * g4
-                   errmax  = 0.0_wp
-                   errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
-!
-!--                Leave loop if accuracy is sufficient, otherwise try again
-!--                with a reduced stepsize
-                   IF ( errmax <= 1.0_wp )  THEN
-                      EXIT
-                   ELSE
-                      dt_ros = MAX( ABS( 0.9_wp * dt_ros * errmax**pshrnk ),   &
-                                    shrnk * ABS( dt_ros ) )
-                   ENDIF
-
-                ENDDO  ! loop for stepsize adjustment
-
-!
-!--             Calculate next internal time step
-                IF ( errmax > errcon )  THEN
-                   dt_ros_next = 0.9_wp * dt_ros * errmax**pgrow
+                                                        ) / ( r_ros + r0 )
+
+                d2rdtdr = -ddenom * ventilation_effect(n) * (                  &
+                                                (e_a / e_s - 1.0) * r_ros**4 - &
+                                                afactor * r0 * r_ros**2 -      &
+                                                2.0 * afactor * r_ros**3 +     &
+                                                3.0 * bfactor * r0 +           &
+                                                4.0 * bfactor * r_ros          &
+                                                            )                  &
+                          / ( r_ros**4 * ( r_ros + r0 )**2 )
+
+                k1    = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr )
+
+                r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1)
+                r_err = r_ros
+
+                drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -   &
+                                                           afactor / r_ros +   &
+                                                           bfactor / r_ros**3  &
+                                                         ) / ( r_ros + r0 )
+
+                k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / &
+                     ( 1.0 - dt_ros * gamma * d2rdtdr )
+
+                r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5 * k2), particles(n)%aux1)
+   !
+   !--          Check error of the solution, and reduce dt_ros if necessary.
+                error = ABS(r_err - r_ros) / r_ros
+                IF ( error .GT. prec )  THEN
+                   dt_ros = SQRT( q_decrease * prec / error ) * dt_ros
+                   r_ros  = r_ros_ini
                 ELSE
-                   dt_ros_next = grow * dt_ros
+                   dt_ros_sum = dt_ros_sum + dt_ros
+                   dt_ros     = q_increase * dt_ros
+                   r_ros_ini  = r_ros
+                   EXIT
                 ENDIF
 
-!
-!--             Estimated time step is reduced if the PALM time step is exceeded
-                IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
-                   dt_ros = dt_3d - dt_ros_sum
-                ELSE
-                   dt_ros = dt_ros_next
-                ENDIF
-
-             ENDDO
-!
-!--          Store internal time step value for next PALM step
-             particles(n)%rvar1 = dt_ros_next
-
-!
-!--          Radius should not fall below 1E-8 because Rosenbrock method may
-!--          lead to errors otherwise
-             new_r(n) = MAX( r_ros, particles(n)%rvar2 )
-!
-!--          Check if calculated droplet radius change is reasonable since in
-!--          case of droplet evaporation the Rosenbrock method may lead to
-!--          secondary solutions which are physically unlikely.
-!--          Due to the solution effect the droplets may grow for relative
-!--          humidities below 100%, but change of radius should not be too 
-!--          large. In case of unreasonable droplet growth the Rosenbrock 
-!--          method is recalculated using a smaller initial time step.
-!--          Limiting values are tested for droplets down to 1.0E-7
-             IF ( new_r(n) - particles(n)%radius >= 3.0E-7_wp  .AND.  &
-                  e_a / e_s < 0.97_wp )  THEN
-                ros_count = ros_count + 1
-                repeat = .TRUE.
-             ENDIF
-
-          ENDDO    ! Rosenbrock method
-
-       ENDIF
-
-       delta_r = new_r(n) - particles(n)%radius
-
-!
-!--    Sum up the change in volume of liquid water for the respective grid
-!--    volume (this is needed later in lpm_calc_liquid_water_content for
-!--    calculating the release of latent heat)
+             END DO
+
+          END DO !Rosenbrock loop
+!
+!--       Store new particle radius
+          new_r(n) = r_ros
+!
+!--       Store internal time step value for next PALM step
+          particles(n)%aux2 = dt_ros
+
+       ENDDO !Particle loop
+
+    ENDIF
+
+    DO  n = 1, number_of_particles
+!
+!--    Sum up the change in liquid water for the respective grid
+!--    box for the computation of the release/depletion of water vapor
+!--    and heat.
        ql_c(kp,jp,ip) = ql_c(kp,jp,ip) + particles(n)%weight_factor *          &
                                    rho_l * 1.33333333_wp * pi *                &
                                    ( new_r(n)**3 - particles(n)%radius**3 ) /  &
                                    ( rho_surface * dx * dy * dz )
+!
+!--    Check if the increase in liqid water is not too big. If this is the case,
+!--    the model timestep might be too long.
        IF ( ql_c(kp,jp,ip) > 100.0_wp )  THEN
           WRITE( message_string, * ) 'k=',kp,' j=',jp,' i=',ip,      &
@@ -499 +368 @@
           CALL message( 'lpm_droplet_condensation', 'PA0143', 2, 2, -1, 6, 1 )
        ENDIF
-
-!
-!--    Change the droplet radius
-       IF ( ( new_r(n) - particles(n)%radius ) < 0.0_wp  .AND.        &
-            new_r(n) < 0.0_wp )  THEN
+!
+!--    Check if the change in the droplet radius is not too big. If this is the
+!--    case, the model timestep might be too long.
+       delta_r = new_r(n) - particles(n)%radius
+       IF ( delta_r < 0.0_wp  .AND. new_r(n) < 0.0_wp )  THEN
           WRITE( message_string, * ) '#1 k=',kp,' j=',jp,' i=',ip,    &
                        ' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int,      &
@@ -510 +379 @@
           CALL message( 'lpm_droplet_condensation', 'PA0144', 2, 2, -1, 6, 1 )
        ENDIF
-
 !
 !--    Sum up the total volume of liquid water (needed below for
 !--    re-calculating the weighting factors)
        ql_v(kp,jp,ip) = ql_v(kp,jp,ip) + particles(n)%weight_factor * new_r(n)**3
-
-       particles(n)%radius = new_r(n)
-
 !
 !--    Determine radius class of the particle needed for collision
-       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  use_kernel_tables )     &
-       THEN
+       IF ( use_kernel_tables )  THEN
           particles(n)%class = ( LOG( new_r(n) ) - rclass_lbound ) /           &
                                ( rclass_ubound - rclass_lbound ) *             &
@@ -528 +392 @@
           particles(n)%class = MAX( particles(n)%class, 1 )
        ENDIF
+ !
+ !--   Store new radius to particle features
+       particles(n)%radius = new_r(n)
 
     ENDDO

palm/trunk/SOURCE/lpm_init.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Extended particle data type. Aerosol initialization improved.
+!
+! 2305 2017-07-06 11:18:47Z hoffmann
 ! Improved calculation of particle IDs.
-! 
+!
 ! 2274 2017-06-09 13:27:48Z Giersch
 !  Changed error messages
-! 
+!
 ! 2265 2017-06-08 16:58:28Z schwenkel
 ! Unused variables removed.
-! 
+!
 ! 2263 2017-06-08 14:59:01Z schwenkel
 ! Implemented splitting and merging algorithm
-! 
+!
 ! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
 ! Adjustments according to new topography realization
-! 
-! 
+!
+!
 ! 2223 2017-05-15 16:38:09Z suehring
 ! Add check for particle release at model top
-! 
+!
 ! 2182 2017-03-17 14:27:40Z schwenkel
 ! Added parameters for simplified particle initialization.
 !
 ! 2122 2017-01-18 12:22:54Z hoffmann
-! Improved initialization of equilibrium aerosol radii 
+! Improved initialization of equilibrium aerosol radii
 ! Calculation of particle ID
 !
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 2016-06-09 16:25:25Z suehring
-! Bugfix in determining initial particle height and grid index in case of 
+! Bugfix in determining initial particle height and grid index in case of
 ! seed_follows_topography.
-! Bugfix concerning random positions, ensure that particles do not move more 
+! Bugfix concerning random positions, ensure that particles do not move more
 ! than one grid length.
 ! Bugfix logarithmic interpolation.
@@ -64 +67 @@
 !
 ! 1890 2016-04-22 08:52:11Z hoffmann
-! Initialization of aerosol equilibrium radius not possible in supersaturated 
+! Initialization of aerosol equilibrium radius not possible in supersaturated
 ! environments. Therefore, a maximum supersaturation of -1 % is assumed during
 ! initialization.
@@ -70 +73 @@
 ! 1873 2016-04-18 14:50:06Z maronga
 ! Module renamed (removed _mod
-! 
+!
 ! 1871 2016-04-15 11:46:09Z hoffmann
 ! Initialization of aerosols added.
@@ -86 +89 @@
 ! netcdf module added
 !
-! 1725 2015-11-17 13:01:51Z hoffmann 
-! Bugfix: Processor-dependent seed for random function is generated before it is 
+! 1725 2015-11-17 13:01:51Z hoffmann
+! Bugfix: Processor-dependent seed for random function is generated before it is
 ! used.
 !
@@ -97 +100 @@
 !
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
+! Code annotations made doxygen readable
 !
 ! 1575 2015-03-27 09:56:27Z raasch
@@ -103 +106 @@
 !
 ! 1359 2014-04-11 17:15:14Z hoffmann
-! New particle structure integrated. 
+! New particle structure integrated.
 ! Kind definition added to all floating point numbers.
 ! lpm_init changed form a subroutine to a module.
-! 
+!
 ! 1327 2014-03-21 11:00:16Z raasch
 ! -netcdf_output
@@ -123 +126 @@
 !
 ! 1314 2014-03-14 18:25:17Z suehring
-! Vertical logarithmic interpolation of horizontal particle speed for particles 
+! Vertical logarithmic interpolation of horizontal particle speed for particles
 ! between roughness height and first vertical grid level.
 !
@@ -146 +149 @@
 !
 ! 667 2010-12-23 12:06:00Z suehring/gryschka
-! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
 ! of arrays.
 !
@@ -159 +162 @@
 !------------------------------------------------------------------------------!
  MODULE lpm_init_mod
- 
+
     USE, INTRINSIC ::  ISO_C_BINDING
 
@@ -167 +170 @@
     USE control_parameters,                                                    &
         ONLY:  cloud_droplets, constant_flux_layer, current_timestep_number,   &
-               dz, initializing_actions, message_string, ocean, simulated_time
+               dt_3d, dz, initializing_actions, message_string, ocean,         &
+               simulated_time
 
     USE grid_variables,                                                        &
@@ -197 +201 @@
                 particles, particle_advection_start, particle_groups,          &
                 particle_groups_type, particles_per_point,                     &
-                particle_type, pdx, pdy, pdz,  prt_count, psb, psl, psn, psr,  & 
-                pss, pst, radius, random_start_position,                       &  
+                particle_type, pdx, pdy, pdz,  prt_count, psb, psl, psn, psr,  &
+                pss, pst, radius, random_start_position,                       &
                 read_particles_from_restartfile, seed_follows_topography,      &
                 sgs_wf_part, sort_count, splitting_function, splitting_mode,   &
@@ -247 +251 @@
     INTEGER(iwp) ::  k                           !<
 
-    REAL(wp) ::  div                             !< 
+    REAL(wp) ::  div                             !<
     REAL(wp) ::  height_int                      !<
     REAL(wp) ::  height_p                        !<
@@ -286 +290 @@
 !-- Check if downward-facing walls exist. This case, reflection boundary
 !-- conditions (as well as subgrid-scale velocities) may do not work
-!-- propably (not realized so far). 
+!-- propably (not realized so far).
     IF ( surf_def_h(1)%ns >= 1 )  THEN
-       WRITE( message_string, * ) 'Overhanging topography do not work '// &   
-                                  'with particles' 
+       WRITE( message_string, * ) 'Overhanging topography do not work '// &
+                                  'with particles'
        CALL message( 'lpm_init', 'PA0212', 0, 1, 0, 6, 0 )
 
-    ENDIF 
+    ENDIF
 
 !
@@ -310 +314 @@
 !-- If number_particles_per_gridbox is set, the parametres pdx, pdy and pdz are
 !-- calculated diagnostically. Therfore an isotropic distribution is prescribed.
-    IF ( number_particles_per_gridbox /= -1 .AND.   & 
+    IF ( number_particles_per_gridbox /= -1 .AND.   &
          number_particles_per_gridbox >= 1 )    THEN
-       pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) /  & 
+       pdx(1) = (( dx * dy * ( zu(2) - zu(1) ) ) /  &
              REAL(number_particles_per_gridbox))**0.3333333_wp
 !
-!--    Ensure a smooth value (two significant digits) of distance between 
+!--    Ensure a smooth value (two significant digits) of distance between
 !--    particles (pdx, pdy, pdz).
        div = 1000.0_wp
@@ -340 +344 @@
 
 !
-!-- Allocate arrays required for calculating particle SGS velocities. 
+!-- Allocate arrays required for calculating particle SGS velocities.
 !-- Initialize prefactor required for stoachastic Weil equation.
     IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
@@ -347 +351 @@
                  de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
-       sgs_wf_part = 1.0_wp / 3.0_wp    
+       sgs_wf_part = 1.0_wp / 3.0_wp
     ENDIF
 
 !
-!-- Allocate array required for logarithmic vertical interpolation of 
+!-- Allocate array required for logarithmic vertical interpolation of
 !-- horizontal particle velocities between the surface and the first vertical
-!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of 
-!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for 
-!-- several heights. Splitting into 20 sublayers turned out to be sufficient. 
+!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
+!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
+!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
 !-- To obtain exact height levels of particles, linear interpolation is applied
 !-- (see lpm_advec.f90).
     IF ( constant_flux_layer )  THEN
-       
-       ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 
+
+       ALLOCATE ( log_z_z0(0:number_of_sublayers) )
        z_p = zu(nzb+1) - zw(nzb)
 
 !
 !--    Calculate horizontal mean value of z0 used for logartihmic
-!--    interpolation. Note: this is not exact for heterogeneous z0. 
-!--    However, sensitivity studies showed that the effect is 
-!--    negligible. 
+!--    interpolation. Note: this is not exact for heterogeneous z0.
+!--    However, sensitivity studies showed that the effect is
+!--    negligible.
        z0_av_local  = SUM( surf_def_h(0)%z0 ) + SUM( surf_lsm_h%z0 ) +         &
                       SUM( surf_usm_h%z0 )
@@ -401 +405 @@
 !-- Check boundary condition and set internal variables
     SELECT CASE ( bc_par_b )
-    
+
        CASE ( 'absorb' )
           ibc_par_b = 1
@@ -407 +411 @@
        CASE ( 'reflect' )
           ibc_par_b = 2
-          
+
        CASE DEFAULT
           WRITE( message_string, * )  'unknown boundary condition ',   &
                                        'bc_par_b = "', TRIM( bc_par_b ), '"'
           CALL message( 'lpm_init', 'PA0217', 1, 2, 0, 6, 0 )
-          
+
     END SELECT
     SELECT CASE ( bc_par_t )
-    
+
        CASE ( 'absorb' )
           ibc_par_t = 1
@@ -421 +425 @@
        CASE ( 'reflect' )
           ibc_par_t = 2
-          
+
        CASE DEFAULT
           WRITE( message_string, * ) 'unknown boundary condition ',   &
                                      'bc_par_t = "', TRIM( bc_par_t ), '"'
           CALL message( 'lpm_init', 'PA0218', 1, 2, 0, 6, 0 )
-          
+
     END SELECT
     SELECT CASE ( bc_par_lr )
@@ -438 +442 @@
        CASE ( 'reflect' )
           ibc_par_lr = 2
-          
+
        CASE DEFAULT
           WRITE( message_string, * ) 'unknown boundary condition ',   &
                                      'bc_par_lr = "', TRIM( bc_par_lr ), '"'
           CALL message( 'lpm_init', 'PA0219', 1, 2, 0, 6, 0 )
-          
+
     END SELECT
     SELECT CASE ( bc_par_ns )
@@ -455 +459 @@
        CASE ( 'reflect' )
           ibc_par_ns = 2
-          
+
        CASE DEFAULT
           WRITE( message_string, * ) 'unknown boundary condition ',   &
                                      'bc_par_ns = "', TRIM( bc_par_ns ), '"'
           CALL message( 'lpm_init', 'PA0220', 1, 2, 0, 6, 0 )
-          
+
     END SELECT
     SELECT CASE ( splitting_mode )
-    
+
        CASE ( 'const' )
           i_splitting_mode = 1
@@ -472 +476 @@
        CASE ( 'gb_av' )
           i_splitting_mode = 3
-          
+
        CASE DEFAULT
           WRITE( message_string, * )  'unknown splitting condition ',   &
                                        'splitting_mode = "', TRIM( splitting_mode ), '"'
           CALL message( 'lpm_init', 'PA0146', 1, 2, 0, 6, 0 )
-          
+
     END SELECT
     SELECT CASE ( splitting_function )
-    
+
        CASE ( 'gamma' )
           isf = 1
@@ -489 +493 @@
        CASE ( 'exp' )
           isf = 3
-          
+
        CASE DEFAULT
           WRITE( message_string, * )  'unknown splitting function ',   &
                                        'splitting_function = "', TRIM( splitting_function ), '"'
           CALL message( 'lpm_init', 'PA0147', 1, 2, 0, 6, 0 )
-          
+
     END SELECT
-    
+
 
 !
@@ -524 +528 @@
 
 !
-!--    initialize counter for particle IDs 
+!--    initialize counter for particle IDs
        grid_particles%id_counter = 1
 
@@ -534 +538 @@
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                       0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
-                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
+                                      0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                       0, 0, 0_idp, .FALSE., -1 )
 
@@ -577 +581 @@
           CALL check_open( 80 )
           WRITE ( 80, 8000 )  current_timestep_number, simulated_time,         &
-                              number_of_particles                             
+                              number_of_particles
           CALL close_file( 80 )
        ENDIF
@@ -584 +588 @@
 
 !
-!-- To avoid programm abort, assign particles array to the local version of 
+!-- To avoid programm abort, assign particles array to the local version of
 !-- first grid cell
     number_of_particles = prt_count(nzb+1,nys,nxl)
@@ -608 +612 @@
 
     USE particle_attributes,                                                   &
-        ONLY: deleted_particles, monodisperse_aerosols
+        ONLY: deleted_particles
 
     IMPLICIT  NONE
@@ -631 +635 @@
     LOGICAL                    ::  first_stride !< flag for initialization
 
-    REAL(wp)                   ::  pos_x      !< increment for particle position in x      
-    REAL(wp)                   ::  pos_y      !< increment for particle position in y  
-    REAL(wp)                   ::  pos_z      !< increment for particle position in z     
+    REAL(wp)                   ::  pos_x      !< increment for particle position in x
+    REAL(wp)                   ::  pos_y      !< increment for particle position in y
+    REAL(wp)                   ::  pos_z      !< increment for particle position in z
     REAL(wp)                   ::  rand_contr !< dummy argument for random position
 
@@ -652 +656 @@
 !--    Calculate initial_weighting_factor diagnostically
        IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN
-          initial_weighting_factor =  number_concentration * 1.0E6_wp *             & 
-                                      pdx(1) * pdy(1) * pdz(1) 
+          initial_weighting_factor =  number_concentration * 1.0E6_wp *             &
+                                      pdx(1) * pdy(1) * pdz(1)
        END IF
 
@@ -677 +681 @@
                xloop: DO WHILE ( pos_x <= psr(i) )
 
-                         IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            & 
+                         IF ( pos_x >= ( nxl - 0.5_wp ) * dx  .AND.            &
                               pos_x <  ( nxr + 0.5_wp ) * dx )  THEN
 
@@ -690 +694 @@
                                tmp_particle%age_m         = 0.0_wp
                                tmp_particle%dt_sum        = 0.0_wp
-                               tmp_particle%user          = 0.0_wp !unused, free for the user
                                tmp_particle%e_m           = 0.0_wp
-                               IF ( curvature_solution_effects )  THEN
-!
-!--                               Initial values (internal timesteps, derivative)
-!--                               for Rosenbrock method
-                                  tmp_particle%rvar1      = 1.0E-6_wp     !last Rosenbrock timestep
-                                  tmp_particle%rvar2      = 0.1E-6_wp     !dry aerosol radius
-                                  tmp_particle%rvar3      = -9999999.9_wp !unused in this configuration
-                               ELSE
-!
-!--                               Initial values for SGS velocities
-                                  tmp_particle%rvar1      = 0.0_wp
-                                  tmp_particle%rvar2      = 0.0_wp
-                                  tmp_particle%rvar3      = 0.0_wp
-                               ENDIF
+                               tmp_particle%rvar1         = 0.0_wp
+                               tmp_particle%rvar2         = 0.0_wp
+                               tmp_particle%rvar3         = 0.0_wp
                                tmp_particle%speed_x       = 0.0_wp
                                tmp_particle%speed_y       = 0.0_wp
@@ -712 +704 @@
                                tmp_particle%origin_y      = pos_y
                                tmp_particle%origin_z      = pos_z
+                               IF ( curvature_solution_effects )  THEN
+                                  tmp_particle%aux1      = 0.0_wp    ! dry aerosol radius
+                                  tmp_particle%aux2      = dt_3d     ! last Rosenbrock timestep
+                               ELSE
+                                  tmp_particle%aux1      = 0.0_wp    ! free to use
+                                  tmp_particle%aux2      = 0.0_wp    ! free to use
+                               ENDIF
                                tmp_particle%radius        = particle_groups(i)%radius
                                tmp_particle%weight_factor = initial_weighting_factor
@@ -726 +725 @@
 !
 !--                            Determine surface level. Therefore, check for
-!--                            upward-facing wall on w-grid. MAXLOC will return 
+!--                            upward-facing wall on w-grid. MAXLOC will return
 !--                            the index of the lowest upward-facing wall.
                                k_surf = MAXLOC(                                &
@@ -748 +747 @@
                                   ENDIF
 !
-!--                            Skip particle release if particle position is 
+!--                            Skip particle release if particle position is
 !--                            below surface, or within topography in case
 !--                            of overhanging structures.
@@ -756 +755 @@
                                THEN
                                   pos_x = pos_x + pdx(i)
-                                  CYCLE xloop                               
+                                  CYCLE xloop
                                ENDIF
 
@@ -840 +839 @@
              particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
 
-             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
+             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
 
                 particles(n)%id = 10000_idp**3 * grid_particles(kp,jp,ip)%id_counter + &
@@ -857 +856 @@
 !
 !-- Initialize aerosol background spectrum
-    IF ( curvature_solution_effects  .AND.  .NOT. monodisperse_aerosols )  THEN
+    IF ( curvature_solution_effects )  THEN
        CALL lpm_init_aerosols(local_start)
     ENDIF
@@ -871 +870 @@
                 particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
 !
-!--             Move only new particles. Moreover, limit random fluctuation 
-!--             in order to prevent that particles move more than one grid box, 
-!--             which would lead to problems concerning particle exchange 
-!--             between processors in case pdx/pdy are larger than dx/dy, 
-!--             respectively.  
+!--             Move only new particles. Moreover, limit random fluctuation
+!--             in order to prevent that particles move more than one grid box,
+!--             which would lead to problems concerning particle exchange
+!--             between processors in case pdx/pdy are larger than dx/dy,
+!--             respectively.
                 DO  n = local_start(kp,jp,ip), number_of_particles
                    IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
@@ -883 +882 @@
                               MERGE( rand_contr, SIGN( dx, rand_contr ),       &
                                      ABS( rand_contr ) < dx                    &
-                                   ) 
+                                   )
                    ENDIF
                    IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
@@ -891 +890 @@
                               MERGE( rand_contr, SIGN( dy, rand_contr ),       &
                                      ABS( rand_contr ) < dy                    &
-                                   ) 
+                                   )
                    ENDIF
                    IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
@@ -899 +898 @@
                               MERGE( rand_contr, SIGN( dz, rand_contr ),       &
                                      ABS( rand_contr ) < dz                    &
-                                   ) 
+                                   )
                    ENDIF
                 ENDDO
@@ -908 +907 @@
 !
 !--             Furthermore, remove particles located in topography. Note, as
-!--             the particle speed is still zero at this point, wall 
+!--             the particle speed is still zero at this point, wall
 !--             reflection boundary conditions will not work in this case.
                 particles =>                                                   &
@@ -934 +933 @@
     ENDIF
 !
-!-- In case of random_start_position, delete particles identified by 
-!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks, 
-!-- which is needed for a fast interpolation of the LES fields on the particle 
+!-- In case of random_start_position, delete particles identified by
+!-- lpm_exchange_horiz and lpm_boundary_conds. Then sort particles into blocks,
+!-- which is needed for a fast interpolation of the LES fields on the particle
 !-- position.
     CALL lpm_pack_all_arrays
@@ -967 +966 @@
 
     USE arrays_3d,                                                             &
-        ONLY: hyp, pt, q  
+        ONLY: hyp, pt, q
 
     USE cloud_parameters,                                                      &
@@ -978 +977 @@
 
     USE particle_attributes,                                                   &
-        ONLY: init_aerosol_probabilistic, molecular_weight_of_solute,          &
-              molecular_weight_of_water, n1, n2, n3, rho_s, rm1, rm2, rm3,     &
-              s1, s2, s3, vanthoff
+        ONLY: aero_type, aero_weight, log_sigma, molecular_weight_of_solute,   &
+              molecular_weight_of_water, na, rho_s, rm, vanthoff
 
     IMPLICIT NONE
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cdf     !< CDF of aerosol spectrum 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_temp  !< dry aerosol radius spectrum
 
     REAL(wp)  :: afactor            !< curvature effects
     REAL(wp)  :: bfactor            !< solute effects
-    REAL(wp)  :: dr                 !< width of radius bin
+    REAL(wp)  :: dlogr              !< logarithmic width of radius bin
     REAL(wp)  :: e_a                !< vapor pressure
     REAL(wp)  :: e_s                !< saturation vapor pressure
-    REAL(wp)  :: n_init             !< sum of all aerosol concentrations
-    REAL(wp)  :: pdf                !< PDF of aerosol spectrum
-    REAL(wp)  :: rmin = 1.0e-8_wp   !< minimum aerosol radius
-    REAL(wp)  :: rmax = 1.0e-6_wp   !< maximum aerosol radius
-    REAL(wp)  :: rs_rand            !< random number
-    REAL(wp)  :: r_mid              !< mean radius
+    REAL(wp)  :: rmin = 0.005e-6_wp !< minimum aerosol radius
+    REAL(wp)  :: rmax = 10.0e-6_wp  !< maximum aerosol radius
+    REAL(wp)  :: r_mid              !< mean radius of bin
+    REAL(wp)  :: r_l                !< left radius of bin
+    REAL(wp)  :: r_r                !< right radius of bin
     REAL(wp)  :: sigma              !< surface tension
     REAL(wp)  :: t_int              !< temperature
-    REAL(wp)  :: weight_sum         !< sum of all weighting factors
 
     INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  local_start !<
 
     INTEGER(iwp)  :: n              !<
-    INTEGER(iwp)  :: nn             !<
-    INTEGER(iwp)  :: no_bins = 999  !< number of bins
     INTEGER(iwp)  :: ip             !<
     INTEGER(iwp)  :: jp             !<
     INTEGER(iwp)  :: kp             !<
 
-    LOGICAL ::  new_pdf = .FALSE.   !< check if aerosol PDF has to be recalculated
-
-!
-!-- Compute aerosol background distribution
-    IF ( init_aerosol_probabilistic )  THEN
-       ALLOCATE( cdf(0:no_bins), r_temp(0:no_bins) )
-       DO n = 0, no_bins
-          r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) /            &
-                           REAL(no_bins, KIND=wp) * REAL(n, KIND=wp) )
-
-          cdf(n) = 0.0_wp
-          n_init = n1 + n2 + n3
-          IF ( n1 > 0.0_wp )  THEN 
-             cdf(n) = cdf(n) + n1 / n_init * ( 0.5_wp + 0.5_wp *        &
-                                  ERF( LOG( r_temp(n) / rm1 ) /         &
-                                       ( SQRT(2.0_wp) * LOG(s1) )       &
-                                     ) )
-          ENDIF
-          IF ( n2 > 0.0_wp )  THEN 
-             cdf(n) = cdf(n) + n2 / n_init * ( 0.5_wp + 0.5_wp *        &
-                                  ERF( LOG( r_temp(n) / rm2 ) /         &
-                                       ( SQRT(2.0_wp) * LOG(s2) )       &
-                                     ) )
-          ENDIF
-          IF ( n3 > 0.0_wp )  THEN 
-             cdf(n) = cdf(n) + n3 / n_init * ( 0.5_wp + 0.5_wp *        &
-                                  ERF( LOG( r_temp(n) / rm3 ) /         &
-                                       ( SQRT(2.0_wp) * LOG(s3) )       &
-                                     ) )
-          ENDIF
-
-       ENDDO
+!
+!-- The following typical aerosol spectra are taken from Jaenicke (1993):
+!-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
+    IF ( TRIM(aero_type) .EQ. 'polar' )  THEN
+       na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6
+       rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6
+       log_sigma = (/ 0.245, 0.300, 0.291 /)
+    ELSEIF ( TRIM(aero_type) .EQ. 'background' )  THEN
+       na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6
+       rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6
+       log_sigma = (/ 0.645, 0.253, 0.425 /)
+    ELSEIF ( TRIM(aero_type) .EQ. 'maritime' )  THEN
+       na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6
+       rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6
+       log_sigma = (/ 0.657, 0.210, 0.396 /)
+    ELSEIF ( TRIM(aero_type) .EQ. 'continental' )  THEN
+       na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6
+       rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6
+       log_sigma = (/ 0.161, 0.217, 0.380 /)
+    ELSEIF ( TRIM(aero_type) .EQ. 'desert' )  THEN
+       na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6
+       rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6
+       log_sigma = (/ 0.247, 0.770, 0.438 /)
+    ELSEIF ( TRIM(aero_type) .EQ. 'rural' )  THEN
+       na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6
+       rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6
+       log_sigma = (/ 0.225, 0.557, 0.266 /)
+    ELSEIF ( TRIM(aero_type) .EQ. 'urban' )  THEN
+       na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6
+       rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6
+       log_sigma = (/ 0.245, 0.666, 0.337 /)
+    ELSEIF ( TRIM(aero_type) .EQ. 'user' )  THEN
+       CONTINUE
+    ELSE
+       WRITE( message_string, * ) 'unknown aerosol type ',   &
+                                'aero_type = "', TRIM( aero_type ), '"'
+       CALL message( 'lpm_init', 'PA0459', 1, 2, 0, 6, 0 )
     ENDIF
 
@@ -1052 +1048 @@
              IF ( number_of_particles <= 0 )  CYCLE
              particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+
+             dlogr   = ( LOG10(rmax) - LOG10(rmin) ) / ( number_of_particles - local_start(kp,jp,ip) + 1 )
 !
 !--          Initialize the aerosols with a predefined spectral distribution
-!--          of the dry radius (logarithmically increasing bins) and a varying 
+!--          of the dry radius (logarithmically increasing bins) and a varying
 !--          weighting factor
-             IF ( .NOT. init_aerosol_probabilistic )  THEN
-
-                new_pdf = .FALSE.
-                IF ( .NOT. ALLOCATED( r_temp ) )  THEN
-                   new_pdf = .TRUE.
+             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
+
+                r_l   = 10.0**( LOG10( rmin ) + (n-1) * dlogr )
+                r_r   = 10.0**( LOG10( rmin ) + n * dlogr )
+                r_mid = SQRT( r_l * r_r )
+
+                particles(n)%aux1          = r_mid
+                particles(n)%weight_factor =                                           &
+                   ( na(1) / ( SQRT( 2.0 * pi ) * log_sigma(1) ) *                     &
+                     EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0 * log_sigma(1)**2 ) ) +  &
+                     na(2) / ( SQRT( 2.0 * pi ) * log_sigma(2) ) *                     &
+                     EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0 * log_sigma(2)**2 ) ) +  &
+                     na(3) / ( SQRT( 2.0 * pi ) * log_sigma(3) ) *                     &
+                     EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0 * log_sigma(3)**2 ) )    &
+                   ) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dz )
+
+!
+!--             Multiply weight_factor with the namelist parameter aero_weight
+!--             to increase or decrease the number of simulated aerosols
+                particles(n)%weight_factor = particles(n)%weight_factor * aero_weight
+
+                IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) &
+                     .GT. random_function( iran_part ) )  THEN
+                   particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0
                 ELSE
-                   IF ( SIZE( r_temp ) .NE. &
-                        number_of_particles - local_start(kp,jp,ip) + 2 )  THEN
-                      new_pdf = .TRUE.
-                      DEALLOCATE( r_temp )
-                   ENDIF
+                   particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp)
                 ENDIF
-
-                IF ( new_pdf )  THEN
-
-                   no_bins = number_of_particles + 1 - local_start(kp,jp,ip)
-                   ALLOCATE( r_temp(0:no_bins) )
-
-                   DO n = 0, no_bins
-                      r_temp(n) = EXP( LOG(rmin) + ( LOG(rmax) - LOG(rmin ) ) / &
-                                       REAL(no_bins, KIND=wp) *                 &
-                                       REAL(n, KIND=wp) )
-                   ENDDO
-
-                ENDIF
-
-!
-!--             Calculate radius and concentration of each aerosol
-                DO n = local_start(kp,jp,ip), number_of_particles
-
-                   nn = n - local_start(kp,jp,ip)
-
-                   r_mid = SQRT( r_temp(nn) * r_temp(nn+1) )
-                   dr    = r_temp(nn+1) - r_temp(nn)
-
-                   pdf    = 0.0_wp
-                   n_init = n1 + n2 + n3
-                   IF ( n1 > 0.0_wp )  THEN 
-                      pdf = pdf + n1 / n_init * ( 1.0_wp / ( r_mid * LOG(s1) *      &
-                                                             SQRT( 2.0_wp * pi )    &
-                                                           ) *                      &
-                                                  EXP( -( LOG( r_mid / rm1 ) )**2 / &
-                                                       ( 2.0_wp * LOG(s1)**2 )      &
-                                                     )                              &
-                                                )
-                   ENDIF
-                   IF ( n2 > 0.0_wp )  THEN
-                      pdf = pdf + n2 / n_init * ( 1.0_wp / ( r_mid * LOG(s2) *      &
-                                                             SQRT( 2.0_wp * pi )    &
-                                                           ) *                      &
-                                                  EXP( -( LOG( r_mid / rm2 ) )**2 / &
-                                                       ( 2.0_wp * LOG(s2)**2 )      &
-                                                     )                              &
-                                                )
-                   ENDIF
-                   IF ( n3 > 0.0_wp )  THEN
-                      pdf = pdf + n3 / n_init * ( 1.0_wp / ( r_mid * LOG(s3) *      &
-                                                             SQRT( 2.0_wp * pi )    &
-                                                           ) *                      &
-                                                  EXP( -( LOG( r_mid / rm3 ) )**2 / &
-                                                       ( 2.0_wp * LOG(s3)**2 )      &
-                                                     )                              &
-                                                )
-                   ENDIF
-
-                   particles(n)%rvar2         = r_mid
-                   particles(n)%weight_factor = pdf * dr
-
-                END DO
-!
-!--             Adjust weighting factors to initialize the same number of aerosols
-!--             in every grid box
-                weight_sum = SUM(particles(local_start(kp,jp,ip):number_of_particles)%weight_factor)
-
-                particles(local_start(kp,jp,ip):number_of_particles)%weight_factor =     &
-                   particles(local_start(kp,jp,ip):number_of_particles)%weight_factor /  &
-                   weight_sum * initial_weighting_factor * ( no_bins + 1 )
-
-             ENDIF
-!
-!--          Initialize the aerosols with a predefined weighting factor but 
-!--          a randomly choosen dry radius
-             IF ( init_aerosol_probabilistic )  THEN
-
-                DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
-
-                   rs_rand = -1.0_wp
-                   DO WHILE ( rs_rand .LT. cdf(0)  .OR.  rs_rand .GE. cdf(no_bins)  )
-                      rs_rand = random_function( iran_part )
-                   ENDDO
-!
-!--                Determine aerosol dry radius by a random number generator
-                   DO nn = 0, no_bins-1
-                      IF ( cdf(nn) .LE. rs_rand  .AND.  cdf(nn+1) .GT. rs_rand )  THEN
-                         particles(n)%rvar2 = r_temp(nn) + ( r_temp(nn+1) - r_temp(nn) ) / &
-                                              ( cdf(nn+1) - cdf(nn) ) * ( rs_rand - cdf(nn) )
-                         EXIT
-                      ENDIF
-                   ENDDO
-
-                ENDDO
-
-             ENDIF
-
+!
+!--             Unnecessary particles will be deleted
+                IF ( particles(n)%weight_factor .LE. 0.0 )  particles(n)%particle_mask = .FALSE.
+
+             ENDDO
 !
 !--          Set particle radius to equilibrium radius based on the environmental
 !--          supersaturation (Khvorostyanov and Curry, 2007, JGR). This avoids
-!--          the sometimes lengthy growth toward their equilibrium radius within 
+!--          the sometimes lengthy growth toward their equilibrium radius within
 !--          the simulation.
              t_int  = pt(kp,jp,ip) * ( hyp(kp) / 100000.0_wp )**0.286_wp
@@ -1176 +1102 @@
                        rho_s / ( molecular_weight_of_solute * rho_l )
 !
-!--          The formula is only valid for subsaturated environments. For 
+!--          The formula is only valid for subsaturated environments. For
 !--          supersaturations higher than -5 %, the supersaturation is set to -5%.
              IF ( e_a / e_s >= 0.95_wp )  e_a = 0.95_wp * e_s
 
-             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
-!
-!--             For details on this equation, see Eq. (14) of Khvorostyanov and 
-!--             Curry (2007, JGR) 
+             DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles
+!
+!--             For details on this equation, see Eq. (14) of Khvorostyanov and
+!--             Curry (2007, JGR)
                 particles(n)%radius = bfactor**0.3333333_wp *                  &
-                   particles(n)%rvar2 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
+                   particles(n)%aux1 / ( 1.0_wp - e_a / e_s )**0.3333333_wp / &
                    ( 1.0_wp + ( afactor / ( 3.0_wp * bfactor**0.3333333_wp *   &
-                     particles(n)%rvar2 ) ) /                                  &
+                     particles(n)%aux1 ) ) /                                  &
                      ( 1.0_wp - e_a / e_s )**0.6666666_wp                      &
                    )
@@ -1196 +1122 @@
        ENDDO
     ENDDO
-!
-!-- Deallocate used arrays
-    IF ( ALLOCATED(r_temp) )  DEALLOCATE( r_temp )
-    IF ( ALLOCATED(cdf) )     DEALLOCATE( cdf )
 
  END SUBROUTINE lpm_init_aerosols

palm/trunk/SOURCE/lpm_read_restart_file.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Extended particle data type.
+!
+! 2305 2017-07-06 11:18:47Z hoffmann
 ! Improved calculation of particle IDs.
-! 
+!
 ! 2265 2017-06-08 16:58:28Z schwenkel
 ! Unused variables removed.
-! 
+!
 ! 2101 2017-01-05 16:42:31Z suehring
 !
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 1822 2016-04-07 07:49:42Z hoffmann
 ! Tails removed. Unused variables removed.
 !
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
+! Code annotations made doxygen readable
+!
 ! 1359 2014-04-11 17:15:14Z hoffmann
-! New particle structure integrated. 
+! New particle structure integrated.
 !
 ! 1320 2014-03-20 08:40:49Z raasch
@@ -61 +64 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_read_restart_file
- 
+
 
     USE control_parameters,                                                    &
@@ -119 +122 @@
 
 !
-!-- If less particles are stored on the restart file than prescribed by 
+!-- If less particles are stored on the restart file than prescribed by
 !-- min_nr_particle, the remainder is initialized by zero_particle to avoid
 !-- errors.
@@ -125 +128 @@
                                    0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                    0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
-                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,          &
+                                   0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp,  &
                                    0, 0, 0_idp, .FALSE., -1 )
 !
@@ -173 +176 @@
     CLOSE ( 90 )
 !
-!-- Must be called to sort particles into blocks, which is needed for a fast 
+!-- Must be called to sort particles into blocks, which is needed for a fast
 !-- interpolation of the LES fields on the particle position.
     CALL lpm_pack_all_arrays

palm/trunk/SOURCE/microphysics_mod.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
-! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
-! includes two more prognostic equations for cloud drop concentration (nc)  
-! and cloud water content (qc). 
-! - The process of activation is parameterized with a simple Twomey 
-!   activion scheme or with considering solution and curvature 
+! s1 changed to log_sigma
+!
+! 2292 2017-06-20 09:51:42Z schwenkel
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
+! includes two more prognostic equations for cloud drop concentration (nc)
+! and cloud water content (qc).
+! - The process of activation is parameterized with a simple Twomey
+!   activion scheme or with considering solution and curvature
 !   effects (Khvorostyanov and Curry ,2006).
 ! - The saturation adjustment scheme is replaced by the parameterization
 !   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
 ! - All other microphysical processes of Seifert and Beheng are used.
-!   Additionally, in those processes the reduction of cloud number concentration 
-!   is considered. 
-! 
+!   Additionally, in those processes the reduction of cloud number concentration
+!   is considered.
+!
 ! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
 ! Adjustments to new topography and surface concept
-! 
+!
 ! 2155 2017-02-21 09:57:40Z hoffmann
 ! Bugfix in the calculation of microphysical quantities on ghost points.
@@ -207 +210 @@
            limiter_sedimentation, na_init, nc_const, sigma_gc,                 &
            ventilation_effect
-           
+
 
     INTERFACE microphysics_control
@@ -463 +466 @@
                    IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
                       nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin *               &
-                                       MERGE( 1.0_wp, 0.0_wp,                  &           
+                                       MERGE( 1.0_wp, 0.0_wp,                  &
                                               BTEST( wall_flags_0(k,j,i), 0 ) )
                    ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
                       nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax *               &
-                                       MERGE( 1.0_wp, 0.0_wp,                  &           
+                                       MERGE( 1.0_wp, 0.0_wp,                  &
                                               BTEST( wall_flags_0(k,j,i), 0 ) )
                    ENDIF
@@ -478 +481 @@
                       IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
                          nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin *            &
-                                          MERGE( 1.0_wp, 0.0_wp,               &           
+                                          MERGE( 1.0_wp, 0.0_wp,               &
                                               BTEST( wall_flags_0(k,j,i), 0 ) )
                       ENDIF
@@ -521 +524 @@
        USE particle_attributes,                                                &
           ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
-              s1, s2, s3, vanthoff
+              log_sigma, vanthoff
 
        IMPLICIT NONE
@@ -578 +581 @@
 !--             Prescribe parameters for activation
 !--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
-                k_act  = 0.7_wp 
+                k_act  = 0.7_wp
 
                 IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
 !
-!--                Compute the number of activated Aerosols 
+!--                Compute the number of activated Aerosols
 !--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
                    n_act     = na_init * sat**k_act
@@ -592 +595 @@
 !
 !--                Compute activation rate after Khairoutdinov and Kogan
-!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
-                   sat_max = 1.0_wp / 100.0_wp                  
-                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       & 
-                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    & 
+!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
+                   sat_max = 1.0_wp / 100.0_wp
+                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
+                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
                        dt_micro
                 ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN
 !
-!--                Curvature effect (afactor) with surface tension 
+!--                Curvature effect (afactor) with surface tension
 !--                parameterization by Straka (2009)
                    sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
@@ -609 +612 @@
 
 !
-!--                Prescribe power index that describes the soluble fraction 
+!--                Prescribe power index that describes the soluble fraction
 !--                of an aerosol particle (beta) and mean geometric radius of
 !--                dry aerosol spectrum
@@ -616 +619 @@
                    rd0      = 0.05E-6_wp
 !
-!--                Calculate mean geometric supersaturation (s_0) with 
+!--                Calculate mean geometric supersaturation (s_0) with
 !--                parameterization by Khvorostyanov and Curry (2006)
                    s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                  &
@@ -622 +625 @@
 
 !
-!--                Calculate number of activated CCN as a function of 
+!--                Calculate number of activated CCN as a function of
 !--                supersaturation and taking Koehler theory into account
-!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)       
-                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              & 
-                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )      
+!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
+                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
+                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * log_sigma(1) ) ) )
                    activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
                 ENDIF
@@ -644 +647 @@
 ! Description:
 ! ------------
-!> Calculate condensation rate for cloud water content (after Khairoutdinov and  
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and
 !> Kogan, 2000).
 !------------------------------------------------------------------------------!
@@ -724 +727 @@
 !--             Mean weight of cloud drops
                 IF ( nc(k,j,i) <= 0.0_wp) CYCLE
-                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)                
+                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
 !
 !--             Weight averaged diameter of cloud drops:
@@ -730 +733 @@
 !
 !--             Integral diameter of cloud drops
-                nc_0 = nc(k,j,i) * dc                
+                nc_0 = nc(k,j,i) * dc
 !
 !--             Condensation needs only to be calculated in supersaturated regions
@@ -741 +744 @@
                    cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
                    cond      = MIN( cond, cond_max / dt_micro )
-               
+
                    qc(k,j,i) = qc(k,j,i) + cond * dt_micro
                 ELSEIF ( sat < 0.0_wp ) THEN
@@ -891 +894 @@
                                                               * flag
                    IF ( microphysics_morrison ) THEN
-                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         & 
+                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
                                     autocon / x0 * hyrho(k) * dt_micro * flag )
                    ENDIF
@@ -1033 +1036 @@
 !
 !--                Mean weight of cloud drops
-                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)      
+                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
 !
 !--                Parameterized turbulence effects on autoconversion (Seifert,
@@ -1321 +1324 @@
 
        USE control_parameters,                                                 &
-           ONLY:  call_microphysics_at_all_substeps,                           & 
+           ONLY:  call_microphysics_at_all_substeps,                           &
                   intermediate_timestep_count, microphysics_morrison
 
@@ -1366 +1369 @@
                 IF ( microphysics_morrison ) THEN
                    IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
-                      sed_nc(k) = sed_qc_const *                               & 
+                      sed_nc(k) = sed_qc_const *                               &
                               ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
                               ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
@@ -1449 +1452 @@
        USE surface_mod,                                                        &
            ONLY :  bc_h
-       
+
        IMPLICIT NONE
 
@@ -1455 +1458 @@
        INTEGER(iwp) ::  j             !< running index y direction
        INTEGER(iwp) ::  k             !< running index z direction
-       INTEGER(iwp) ::  k_run         !< 
+       INTEGER(iwp) ::  k_run         !<
        INTEGER(iwp) ::  l             !< running index of surface type
        INTEGER(iwp) ::  m             !< running index surface elements
@@ -1528 +1531 @@
              ENDDO
 !
-!--          Adjust boundary values using surface data type. 
+!--          Adjust boundary values using surface data type.
 !--          Upward-facing
-             surf_s = bc_h(0)%start_index(j,i)   
+             surf_s = bc_h(0)%start_index(j,i)
              surf_e = bc_h(0)%end_index(j,i)
              DO  m = surf_s, surf_e
@@ -1539 +1542 @@
 !
 !--          Downward-facing
-             surf_s = bc_h(1)%start_index(j,i)   
+             surf_s = bc_h(1)%start_index(j,i)
              surf_e = bc_h(1)%end_index(j,i)
              DO  m = surf_s, surf_e
@@ -1607 +1610 @@
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 !
-!--             Sum up all rain drop number densities which contribute to the flux 
+!--             Sum up all rain drop number densities which contribute to the flux
 !--             through k-1/2
                 flux  = 0.0_wp
@@ -1730 +1733 @@
        THEN
 !
-!--       Run over all upward-facing surface elements, i.e. non-natural, 
+!--       Run over all upward-facing surface elements, i.e. non-natural,
 !--       natural and urban
           DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)            
+             i = bc_h(0)%i(m)
              j = bc_h(0)%j(m)
              k = bc_h(0)%k(m)
@@ -1961 +1964 @@
        USE particle_attributes,                                                &
           ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
-              s1, s2, s3, vanthoff
+              log_sigma, vanthoff
 
        IMPLICIT NONE
@@ -2013 +2016 @@
 !--       Prescribe parameters for activation
 !--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
-          k_act  = 0.7_wp 
+          k_act  = 0.7_wp
 
           IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
 !
-!--          Compute the number of activated Aerosols 
+!--          Compute the number of activated Aerosols
 !--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
              n_act     = na_init * sat**k_act
@@ -2027 +2030 @@
 !
 !--          Compute activation rate after Khairoutdinov and Kogan
-!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
-             sat_max = 0.8_wp / 100.0_wp                  
-             activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              & 
-                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          & 
+!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
+             sat_max = 0.8_wp / 100.0_wp
+             activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              &
+                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          &
                   dt_micro
 
@@ -2036 +2039 @@
           ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk )  THEN
 !
-!--          Curvature effect (afactor) with surface tension 
+!--          Curvature effect (afactor) with surface tension
 !--          parameterization by Straka (2009)
              sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
@@ -2046 +2049 @@
 
 !
-!--          Prescribe power index that describes the soluble fraction 
+!--          Prescribe power index that describes the soluble fraction
 !--          of an aerosol particle (beta) and mean geometric radius of
 !--          dry aerosol spectrum
@@ -2053 +2056 @@
              rd0      = 0.05E-6_wp
 !
-!--          Calculate mean geometric supersaturation (s_0) with 
+!--          Calculate mean geometric supersaturation (s_0) with
 !--          parameterization by Khvorostyanov and Curry (2006)
              s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                        &
@@ -2059 +2062 @@
 
 !
-!--          Calculate number of activated CCN as a function of 
+!--          Calculate number of activated CCN as a function of
 !--          supersaturation and taking Koehler theory into account
 !--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
-             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                     & 
-                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )      
+             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                     &
+                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * log_sigma(1) ) ) )
              activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
-                                      
-             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)                   
+
+             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
           ENDIF
 
@@ -2076 +2079 @@
 ! Description:
 ! ------------
-!> Calculate condensation rate for cloud water content (after Khairoutdinov and  
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and
 !> Kogan, 2000).
 !------------------------------------------------------------------------------!
@@ -2154 +2157 @@
 !--       Mean weight of cloud drops
           IF ( nc_1d(k) <= 0.0_wp) CYCLE
-          xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)                
+          xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)
 !
 !--       Weight averaged diameter of cloud drops:
@@ -2160 +2163 @@
 !
 !--       Integral diameter of cloud drops
-          nc_0 = nc_1d(k) * dc                
+          nc_0 = nc_1d(k) * dc
 !
 !--       Condensation needs only to be calculated in supersaturated regions
@@ -2171 +2174 @@
              cond_max  = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
              cond      = MIN( cond, cond_max / dt_micro )
-               
+
              qc_1d(k) = qc_1d(k) + cond * dt_micro
           ELSEIF ( sat < 0.0_wp ) THEN
@@ -2302 +2305 @@
              nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
              IF ( microphysics_morrison ) THEN
-                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                & 
+                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                &
                                     autocon / x0 * hyrho(k) * dt_micro * flag )
              ENDIF
@@ -2425 +2428 @@
 !
 !--          Mean weight of cloud drops
-             xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)      
+             xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)
 !
 !--          Parameterized turbulence effects on autoconversion (Seifert,
@@ -2716 +2719 @@
           IF ( microphysics_morrison ) THEN
              IF ( qc_1d(k) > eps_sb  .AND.  nc_1d(k) > eps_mr )  THEN
-                sed_nc(k) = sed_qc_const *                                     & 
+                sed_nc(k) = sed_qc_const *                                     &
                             ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
                             ( nc_1d(k) )**( 1.0_wp / 3.0_wp )
@@ -2790 +2793 @@
        USE surface_mod,                                                        &
            ONLY :  bc_h
-       
+
        IMPLICIT NONE
 
@@ -2859 +2862 @@
        ENDDO
 !
-!--    Adjust boundary values using surface data type. 
+!--    Adjust boundary values using surface data type.
 !--    Upward facing non-natural
-       surf_s = bc_h(0)%start_index(j,i)   
+       surf_s = bc_h(0)%start_index(j,i)
        surf_e = bc_h(0)%end_index(j,i)
        DO  m = surf_s, surf_e
@@ -2870 +2873 @@
 !
 !--    Downward facing non-natural
-       surf_s = bc_h(1)%start_index(j,i)   
+       surf_s = bc_h(1)%start_index(j,i)
        surf_e = bc_h(1)%end_index(j,i)
        DO  m = surf_s, surf_e
@@ -2935 +2938 @@
           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 !
-!--       Sum up all rain drop number densities which contribute to the flux 
+!--       Sum up all rain drop number densities which contribute to the flux
 !--       through k-1/2
           flux  = 0.0_wp
@@ -3044 +3047 @@
        THEN
 
-          surf_s = bc_h(0)%start_index(j,i)   
+          surf_s = bc_h(0)%start_index(j,i)
           surf_e = bc_h(0)%end_index(j,i)
           DO  m = surf_s, surf_e

palm/trunk/SOURCE/mod_particle_attributes.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Aerosol initialization improved.
+!
+! 2305 2017-07-06 11:18:47Z hoffmann
 ! Improved calculation of particle IDs.
-! 
+!
 ! 2278 2017-06-12 13:08:18Z schwenkel
 ! Added comments
-! 
+!
 ! 2265 2017-06-08 16:58:28Z schwenkel
 ! Unused variables removed.
-! 
+!
 ! 2263 2017-06-08 14:59:01Z schwenkel
 ! Implemented splitting and merging algorithm
-! 
+!
 ! 2183 2017-03-17 14:29:15Z schwenkel
 !
 ! 2182 2017-03-17 14:27:40Z schwenkel
 ! Added parameters for simplified particle initialization.
-! 
+!
 ! 2122 2017-01-18 12:22:54Z hoffmann
 ! Calculation of particle ID
@@ -48 +51 @@
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 1936 2016-06-13 13:37:44Z suehring
 ! +deallocate_memory, step_dealloc
@@ -71 +74 @@
 !
 ! 1727 2015-11-20 07:22:02Z knoop
-! Bugfix: Cause of syntax warning gfortran preprocessor removed 
-! 
+! Bugfix: Cause of syntax warning gfortran preprocessor removed
+!
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
+! Code annotations made doxygen readable
 !
 ! 1575 2015-03-27 09:56:27Z raasch
@@ -88 +91 @@
 !------------------------------------------------------------------------------!
 MODULE particle_attributes
- 
+
 
     USE kinds
 
-    CHARACTER(LEN=15) ::  bc_par_lr = 'cyclic'                    !< left/right boundary condition
-    CHARACTER(LEN=15) ::  bc_par_ns = 'cyclic'                    !< north/south boundary condition
-    CHARACTER(LEN=15) ::  bc_par_b  = 'reflect'                   !< bottom boundary condition
-    CHARACTER(LEN=15) ::  bc_par_t  = 'absorb'                    !< top boundary condition
-    CHARACTER(LEN=15) ::  collision_algorithm = 'all_or_nothing'  !< collision algorithm
-    CHARACTER(LEN=15) ::  collision_kernel = 'none'               !< collision kernel
-    CHARACTER(LEN=5) ::  splitting_function = 'gamma'             !< function for calculation critical weighting factor
-    CHARACTER(LEN=5) ::  splitting_mode = 'const'                 !< splitting mode
-
-
-    INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step 
+    CHARACTER(LEN=15) ::  aero_type = 'maritime'                   !< aerosol type
+    CHARACTER(LEN=15) ::  bc_par_lr = 'cyclic'                     !< left/right boundary condition
+    CHARACTER(LEN=15) ::  bc_par_ns = 'cyclic'                     !< north/south boundary condition
+    CHARACTER(LEN=15) ::  bc_par_b  = 'reflect'                    !< bottom boundary condition
+    CHARACTER(LEN=15) ::  bc_par_t  = 'absorb'                     !< top boundary condition
+    CHARACTER(LEN=15) ::  collision_kernel = 'none'                !< collision kernel
+    CHARACTER(LEN=5)  ::  splitting_function = 'gamma'             !< function for calculation critical weighting factor
+    CHARACTER(LEN=5)  ::  splitting_mode = 'const'                 !< splitting mode
+
+    INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step
     INTEGER(iwp) ::  dissipation_classes = 10                     !< namelist parameter (see documentation)
     INTEGER(iwp) ::  ibc_par_b                                    !< particle bottom boundary condition dummy
@@ -113 +115 @@
     INTEGER(iwp) ::  max_number_particles_per_gridbox = 100       !< namelist parameter (see documentation)
     INTEGER(iwp) ::  merge_drp = 0                                !< number of merged droplets
-    INTEGER(iwp) ::  min_nr_particle = 50                         !< namelist parameter (see documentation)         
+    INTEGER(iwp) ::  min_nr_particle = 50                         !< namelist parameter (see documentation)
     INTEGER(iwp) ::  new_particles = 0                            !< number of new particles
-    INTEGER(iwp) ::  n_max = 100                                  !< number of radii bin for splitting functions     
-    INTEGER(iwp) ::  number_of_particles = 0                      !< number of particles for each grid box (3d array is saved on prt_count)           
-    INTEGER(iwp) ::  number_of_particle_groups = 1                !< namelist parameter (see documentation)            
-    INTEGER(iwp) ::  number_of_sublayers = 20                     !< number of sublayers for particle velocities betwenn surface and first grid level         
-    INTEGER(iwp) ::  number_particles_per_gridbox = -1            !< namelist parameter (see documentation)          
-    INTEGER(iwp) ::  offset_ocean_nzt = 0                         !< in case of oceans runs, the vertical index calculations need an offset         
+    INTEGER(iwp) ::  n_max = 100                                  !< number of radii bin for splitting functions
+    INTEGER(iwp) ::  number_of_particles = 0                      !< number of particles for each grid box (3d array is saved on prt_count)
+    INTEGER(iwp) ::  number_of_particle_groups = 1                !< namelist parameter (see documentation)
+    INTEGER(iwp) ::  number_of_sublayers = 20                     !< number of sublayers for particle velocities betwenn surface and first grid level
+    INTEGER(iwp) ::  number_particles_per_gridbox = -1            !< namelist parameter (see documentation)
+    INTEGER(iwp) ::  offset_ocean_nzt = 0                         !< in case of oceans runs, the vertical index calculations need an offset
     INTEGER(iwp) ::  offset_ocean_nzt_m1 = 0                      !< in case of oceans runs, the vertical index calculations need an offset
     INTEGER(iwp) ::  particles_per_point = 1                      !< namelist parameter (see documentation)
@@ -130 +132 @@
     INTEGER(iwp) ::  sum_merge_drp = 0                            !< sum of merged super droplets
     INTEGER(iwp) ::  sum_new_particles = 0                        !< sum of created particles (in splitting algorithm)
-    INTEGER(iwp) ::  total_number_of_particles                    !< total number of particles in the whole model domain           
+    INTEGER(iwp) ::  total_number_of_particles                    !< total number of particles in the whole model domain
     INTEGER(iwp) ::  trlp_count_sum                               !< parameter for particle exchange of PEs
     INTEGER(iwp) ::  trlp_count_recv_sum                          !< parameter for particle exchange of PEs
@@ -144 +146 @@
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  prt_count  !< 3d array of number of particles of every grid box
 
-    LOGICAL ::  all_or_nothing = .FALSE.                  !< flag for collision algorithm
-    LOGICAL ::  average_impact = .FALSE.                  !< flag for collision algortihm
-    LOGICAL ::  curvature_solution_effects = .FALSE.      !< namelist parameter (see documentation)                     
-    LOGICAL ::  deallocate_memory = .TRUE.                !< namelist parameter (see documentation)                  
+    LOGICAL ::  curvature_solution_effects = .FALSE.      !< namelist parameter (see documentation)
+    LOGICAL ::  deallocate_memory = .TRUE.                !< namelist parameter (see documentation)
     LOGICAL ::  hall_kernel = .FALSE.                     !< flag for collision kernel
-    LOGICAL ::  init_aerosol_probabilistic = .FALSE.      !< namelist parameter (see documentation)
     LOGICAL ::  merging = .FALSE.                         !< namelist parameter (see documentation)
-    LOGICAL ::  monodisperse_aerosols = .FALSE.           !< namelist parameter (see documentation)
     LOGICAL ::  particle_advection = .FALSE.              !< parameter to steer the advection of particles
-    LOGICAL ::  random_start_position = .FALSE.           !< namelist parameter (see documentation)                   
-    LOGICAL ::  read_particles_from_restartfile = .TRUE.  !< namelist parameter (see documentation)                   
+    LOGICAL ::  random_start_position = .FALSE.           !< namelist parameter (see documentation)
+    LOGICAL ::  read_particles_from_restartfile = .TRUE.  !< namelist parameter (see documentation)
     LOGICAL ::  seed_follows_topography = .FALSE.         !< namelist parameter (see documentation)
     LOGICAL ::  splitting = .FALSE.                       !< namelist parameter (see documentation)
     LOGICAL ::  use_kernel_tables = .FALSE.               !< parameter, which turns on the use of precalculated collision kernels
-    LOGICAL ::  use_sgs_for_particles = .FALSE.           !< namelist parameter (see documentation) 
+    LOGICAL ::  use_sgs_for_particles = .FALSE.           !< namelist parameter (see documentation)
     LOGICAL ::  wang_kernel = .FALSE.                     !< flag for collision kernel
     LOGICAL ::  write_particle_statistics = .FALSE.       !< namelist parameter (see documentation)
-                
+
     LOGICAL, DIMENSION(max_number_of_particle_groups) ::                       &
                 vertical_particle_advection = .TRUE.              !< Switch on/off vertical particle transport
 
+    REAL(wp) ::  aero_weight = 1.0_wp                      !< namelist parameter (see documentation)
     REAL(wp) ::  alloc_factor = 20.0_wp                    !< namelist parameter (see documentation)
     REAL(wp) ::  c_0 = 3.0_wp                              !< parameter for lagrangian timescale
     REAL(wp) ::  dt_min_part = 0.0002_wp                   !< minimum particle time step when SGS velocities are used (s)
     REAL(wp) ::  dt_prel = 9999999.9_wp                    !< namelist parameter (see documentation)
-    REAL(wp) ::  dt_write_particle_data = 9999999.9_wp     !< namelist parameter (see documentation)            
-    REAL(wp) ::  end_time_prel = 9999999.9_wp              !< namelist parameter (see documentation)         
+    REAL(wp) ::  dt_write_particle_data = 9999999.9_wp     !< namelist parameter (see documentation)
+    REAL(wp) ::  end_time_prel = 9999999.9_wp              !< namelist parameter (see documentation)
     REAL(wp) ::  initial_weighting_factor = 1.0_wp         !< namelist parameter (see documentation)
+    REAL(wp) ::  log_sigma(3) = 1.0_wp                     !< namelist parameter (see documentation)
     REAL(wp) ::  molecular_weight_of_solute = 0.05844_wp   !< mol. m. NaCl (kg mol-1)
     REAL(wp) ::  molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1)
-    REAL(wp) ::  n1 = 100.0_wp                             !< namelist parameter (see documentation)
-    REAL(wp) ::  n2 = 0.0_wp                               !< namelist parameter (see documentation)
-    REAL(wp) ::  n3 = 0.0_wp                               !< namelist parameter (see documentation)
+    REAL(wp) ::  na(3) = 0.0_wp                            !< namelist parameter (see documentation)
     REAL(wp) ::  number_concentration = -1.0_wp            !< namelist parameter (see documentation)
     REAL(wp) ::  particle_advection_start = 0.0_wp         !< namelist parameter (see documentation)
     REAL(wp) ::  radius_merge = 1.0E-7_wp                  !< namelist parameter (see documentation)
     REAL(wp) ::  radius_split = 40.0E-6_wp                 !< namelist parameter (see documentation)
-    REAL(wp) ::  rho_s = 2165.0_wp                         !< density of NaCl (kg m-3) 
-    REAL(wp) ::  rm1 = 0.05E-6_wp                          !< namelist parameter (see documentation)
-    REAL(wp) ::  rm2 = 0.05E-6_wp                          !< namelist parameter (see documentation)
-    REAL(wp) ::  rm3 = 0.05E-6_wp                          !< namelist parameter (see documentation)
-    REAL(wp) ::  s1 = 2.0_wp                               !< namelist parameter (see documentation)
-    REAL(wp) ::  s2 = 2.0_wp                               !< namelist parameter (see documentation) 
-    REAL(wp) ::  s3 = 2.0_wp                               !< namelist parameter (see documentation)
-    REAL(wp) ::  sgs_wf_part                               !< parameter for sgs 
+    REAL(wp) ::  rho_s = 2165.0_wp                         !< density of NaCl (kg m-3)
+    REAL(wp) ::  rm(3) = 1.0E-6_wp                         !< namelist parameter (see documentation)
+    REAL(wp) ::  sgs_wf_part                               !< parameter for sgs
     REAL(wp) ::  time_prel = 0.0_wp                        !< time for particle release
     REAL(wp) ::  time_write_particle_data = 0.0_wp         !< write particle data at current time on file
@@ -210 +203 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE     ::  log_z_z0   !< Precalculate LOG(z/z0)
 
-    
+
     TYPE particle_type
+        REAL(wp)     ::  aux1          !< auxiliary multi-purpose feature
+        REAL(wp)     ::  aux2          !< auxiliary multi-purpose feature
         REAL(wp)     ::  radius        !< radius of particle
         REAL(wp)     ::  age           !< age of particle
-        REAL(wp)     ::  age_m         !< 
+        REAL(wp)     ::  age_m         !<
         REAL(wp)     ::  dt_sum        !<
-        REAL(wp)     ::  user          !< varible for user
-        REAL(wp)     ::  e_m           !< interpolated sgs tke      
+        REAL(wp)     ::  e_m           !< interpolated sgs tke
         REAL(wp)     ::  origin_x      !< origin x-position of particle (changed cyclic bc)
         REAL(wp)     ::  origin_y      !< origin y-position of particle (changed cyclic bc)
         REAL(wp)     ::  origin_z      !< origin z-position of particle (changed cyclic bc)
-        REAL(wp)     ::  rvar1         !< 
+        REAL(wp)     ::  rvar1         !<
         REAL(wp)     ::  rvar2         !<
-        REAL(wp)     ::  rvar3         !< 
+        REAL(wp)     ::  rvar3         !<
         REAL(wp)     ::  speed_x       !< speed of particle in x
         REAL(wp)     ::  speed_y       !< speed of particle in y
         REAL(wp)     ::  speed_z       !< speed of particle in z
-        REAL(wp)     ::  weight_factor !< weighting factor  
-        REAL(wp)     ::  x             !< x-position 
-        REAL(wp)     ::  y             !< y-position 
-        REAL(wp)     ::  z             !< z-position 
+        REAL(wp)     ::  weight_factor !< weighting factor
+        REAL(wp)     ::  x             !< x-position
+        REAL(wp)     ::  y             !< y-position
+        REAL(wp)     ::  z             !< z-position
         INTEGER(iwp) ::  class         !< radius class needed for collision
         INTEGER(iwp) ::  group         !< number of particle group
@@ -274 +268 @@
 
 END MODULE particle_attributes
-

palm/trunk/SOURCE/package_parin.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+!
+!
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Aerosol initialization improved.
+!
+! 2263 2017-06-08 14:59:01Z schwenkel
 ! Implemented splitting and merging algorithm
-! 
+!
 ! 2183 2017-03-17 14:29:15Z schwenkel
 !
 ! 2182 2017-03-17 14:27:40Z schwenkel
 ! Added parameters for simplified particle initialization.
-! 
+!
 ! 2000 2016-08-20 18:09:15Z knoop
 ! Forced header and separation lines into 80 columns
-! 
+!
 ! 1936 2016-06-13 13:37:44Z suehring
 ! +deallocate_memory, step_dealloc
-! 
+!
 ! 1871 2016-04-15 11:46:09Z hoffmann
 ! Initialization of aerosols added.
-! 
+!
 ! 1833 2016-04-07 14:23:03Z raasch
 ! reading of spectra_par moved to spectra_mod
-! 
+!
 ! 1831 2016-04-07 13:15:51Z hoffmann
 ! curvature_solution_effects added
@@ -50 +53 @@
 ! Reading of &radiation_par moved to radiation_model_mod.
 ! Reading of &canopy_par moved to plant_canopy_model_mod.
-! 
+!
 ! 822 2016-04-07 07:49:42Z hoffmann
 ! +collision_algorithm
@@ -57 +60 @@
 ! 1817 2016-04-06 15:44:20Z maronga
 ! Reading of &lsm_par moved to land_surface_model_mod.
-! 
+!
 ! 1788 2016-03-10 11:01:04Z maronga
 ! Parameter dewfall removed.
@@ -66 +69 @@
 ! 1757 2016-02-22 15:49:32Z maronga
 ! Added parameter unscheduled_radiation_calls
-! 
+!
 ! 1691 2015-10-26 16:17:44Z maronga
 ! Added skip_time_do_lsm, skip_time_do_radiation, and emissivity
-! 
+!
 ! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-! 
+! Code annotations made doxygen readable
+!
 ! 1585 2015-04-30 07:05:52Z maronga
 ! Added several radiation_par parameters
@@ -81 +84 @@
 ! 1553 2015-03-03 17:33:54Z maronga
 ! Resorting of lsm_par
-! 
+!
 ! 1551 2015-03-03 14:18:16Z maronga
 ! Several changes in the liste for land surface model and radiation model
-! 
+!
 ! 1496 2014-12-02 17:25:50Z maronga
 ! Added support for the land surface model and radiation scheme
-! 
+!
 ! 1484 2014-10-21 10:53:05Z kanani
 ! Changes due to new module structure of the plant canopy model:
 !   module plant_canopy_model_mod added,
-!   new package/namelist canopy_par added, i.e. the canopy model is no longer 
+!   new package/namelist canopy_par added, i.e. the canopy model is no longer
 !   steered over the inipar-namelist,
 !   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
 !   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff.
 ! Changed statement tags in CONTINUE-statement
-! 
+!
 ! 1367 2014-04-23 15:18:30Z witha
 ! Bugfix: module kinds must be used
@@ -103 +106 @@
 ! +alloc_factor, + min_nr_particle
 ! -dt_sort_particles, -maximum_number_of_particles
-! 
+!
 ! 1340 2014-03-25 19:45:13Z kanani
 ! REAL constants defined as wp-kinds
@@ -138 +141 @@
 !> software packages which are used optionally in the run.
 !>
-!> @todo Perform all actions in the respective submodules and remove 
+!> @todo Perform all actions in the respective submodules and remove
 !>       package_parin
 !------------------------------------------------------------------------------!
  SUBROUTINE package_parin
- 
+
 
     USE control_parameters,                                                    &
@@ -165 +168 @@
 
     USE particle_attributes,                                                   &
-        ONLY:  alloc_factor, bc_par_b, bc_par_lr, bc_par_ns, bc_par_t,         &
-               collision_algorithm, collision_kernel,                          &
+        ONLY:  aero_type, aero_weight, alloc_factor, bc_par_b, bc_par_lr,      &
+               bc_par_ns, bc_par_t, collision_kernel,                          &
                curvature_solution_effects, deallocate_memory, density_ratio,   &
                dissipation_classes, dt_min_part, dt_prel,                      &
                dt_write_particle_data, end_time_prel, initial_weighting_factor,&
-               init_aerosol_probabilistic, max_number_particles_per_gridbox,   & 
-               merging, min_nr_particle, monodisperse_aerosols, n1, n2, n3,    &
+               log_sigma, max_number_particles_per_gridbox,                    &
+               merging, min_nr_particle, na,                                   &
                number_concentration, number_of_particle_groups,                &
                number_particles_per_gridbox, particles_per_point,              &
@@ -177 +180 @@
                psb, psl, psn, psr, pss, pst, radius, radius_classes,           &
                radius_merge, radius_split, random_start_position,              &
-               read_particles_from_restartfile, rm1, rm2, rm3,                 &
-               seed_follows_topography, splitting, splitting_factor,           & 
+               read_particles_from_restartfile, rm,                            &
+               seed_follows_topography, splitting, splitting_factor,           &
                splitting_factor_max, splitting_function, splitting_mode,       &
-               step_dealloc, s1, s2, s3, use_sgs_for_particles,                &
+               step_dealloc, use_sgs_for_particles,                            &
                vertical_particle_advection, weight_factor_merge,               &
-               weight_factor_split, write_particle_statistics                
+               weight_factor_split, write_particle_statistics
 
     IMPLICIT NONE
@@ -206 +209 @@
                                   vc_size_y, vc_size_z
 
-    NAMELIST /particles_par/      alloc_factor, bc_par_b, bc_par_lr,           &
-                                  bc_par_ns, bc_par_t, collision_algorithm,    &
+    NAMELIST /particles_par/      aero_type, aero_weight, alloc_factor,        &
+                                  bc_par_b, bc_par_lr,                         &
+                                  bc_par_ns, bc_par_t,                         &
                                   collision_kernel, curvature_solution_effects,&
                                   deallocate_memory, density_ratio,            &
@@ -214 +218 @@
                                   dt_write_particle_data,                      &
                                   end_time_prel, initial_weighting_factor,     &
-                                  init_aerosol_probabilistic,                  &
+                                  log_sigma,                                   &
                                   max_number_particles_per_gridbox, merging,   &
-                                  min_nr_particle, monodisperse_aerosols,      &
-                                  n1, n2, n3, number_concentration,            &
+                                  min_nr_particle,                             &
+                                  na, number_concentration,                    &
                                   number_of_particle_groups,                   &
                                   number_particles_per_gridbox,                &
@@ -224 +228 @@
                                   particle_maximum_age, pdx, pdy, pdz, psb,    &
                                   psl, psn, psr, pss, pst, radius,             &
-                                  radius_classes, radius_merge, radius_split,  &   
+                                  radius_classes, radius_merge, radius_split,  &
                                   random_start_position,                       &
-                                  read_particles_from_restartfile,             &
-                                  rm1, rm2, rm3,                               &
+                                  read_particles_from_restartfile, rm,         &
                                   seed_follows_topography,                     &
                                   splitting, splitting_factor,                 &
                                   splitting_factor_max, splitting_function,    &
-                                  splitting_mode, step_dealloc, s1, s2, s3,    &
+                                  splitting_mode, step_dealloc,                &
                                   use_sgs_for_particles,                       &
                                   vertical_particle_advection,                 &