Index: palm/trunk/SOURCE/lpm_init_mod.f90
===================================================================
--- palm/trunk/SOURCE/lpm_init_mod.f90	(revision 1851)
+++ palm/trunk/SOURCE/lpm_init_mod.f90	(revision 1871)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Initialization of aerosols added.
 ! 
 ! Former revisions:
@@ -509,4 +509,7 @@
         ONLY: lpm_pack_all_arrays
 
+    USE particle_attributes,                                                   &
+        ONLY: monodisperse_aerosols
+
     IMPLICIT  NONE
 
@@ -581,7 +584,7 @@
 !--                            Initial values (internal timesteps, derivative)
 !--                            for Rosenbrock method
-                               tmp_particle%rvar1      = 1.0E-12_wp
-                               tmp_particle%rvar2      = 1.0E-3_wp
-                               tmp_particle%rvar3      = -9999999.9_wp
+                               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
                             ELSE
 !
@@ -686,4 +689,11 @@
     local_start = prt_count+1
     prt_count   = local_count
+
+!
+!-- Initialize aerosol background spectrum
+    IF ( curvature_solution_effects  .AND.  .NOT. monodisperse_aerosols )  THEN
+       CALL lpm_init_aerosols(local_start)
+    ENDIF
+
 !
 !-- Add random fluctuation to particle positions
@@ -759,3 +769,235 @@
  END SUBROUTINE lpm_create_particle
 
+ SUBROUTINE lpm_init_aerosols(local_start)
+
+    USE arrays_3d,                                                             &
+        ONLY: hyp, pt, q  
+
+    USE cloud_parameters,                                                      &
+        ONLY: l_d_rv, rho_l
+
+    USE constants,                                                             &
+        ONLY: pi
+
+    USE kinds
+
+    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
+
+    IMPLICIT NONE
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cdf     !< CDF of aerosol spectrum 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_temp  !< dry aerosol radius spectrum
+
+    REAL(wp)  :: bfactor            !< solute effects
+    REAL(wp)  :: dr                 !< 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)  :: t_int              !< temperature
+    REAL(wp)  :: weight_sum         !< sum of all weighting factors
+
+    INTEGER(iwp), DIMENSION(:,:,:), 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
+    ENDIF
+
+    DO  ip = nxl, nxr
+       DO  jp = nys, nyn
+          DO  kp = nzb+1, nzt
+
+             number_of_particles = prt_count(kp,jp,ip)
+             IF ( number_of_particles <= 0 )  CYCLE
+             particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+!
+!--          Initialize the aerosols with a predefined spectral distribution
+!--          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.
+                ELSE
+                   IF ( SIZE( r_temp ) .NE. &
+                        number_of_particles - local_start(kp,jp,ip) + 2 )  THEN
+                      new_pdf = .TRUE.
+                      DEALLOCATE( r_temp )
+                   ENDIF
+                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
+
+!
+!--          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 simulation.
+             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 )
+
+!
+!--          The formula is only valid for subsaturated environments. In (super-)
+!--          saturated air, the inital radius is used.
+             IF ( e_a / e_s < 1.0_wp )  THEN
+
+                DO  n = local_start(kp,jp,ip), number_of_particles  !only new particles 
+
+                   bfactor             = vanthoff * molecular_weight_of_water *    &
+                                         rho_s * particles(n)%rvar2**3 /           &
+                                         ( molecular_weight_of_solute * rho_l )
+                   particles(n)%radius = particles(n)%rvar2 * ( bfactor /          &
+                                         particles(n)%rvar2**3 )**(1.0_wp/3.0_wp) *&
+                                         ( 1.0_wp - e_a / e_s )**(-1.0_wp/3.0_wp)
+
+                ENDDO
+
+             ENDIF
+
+          ENDDO
+       ENDDO
+    ENDDO
+!
+!-- Deallocate used arrays
+    IF ( ALLOCATED(r_temp) )  DEALLOCATE( r_temp )
+    IF ( ALLOCATED(cdf) )     DEALLOCATE( cdf )
+
+ END SUBROUTINE lpm_init_aerosols
+
 END MODULE lpm_init_mod
