Index: palm/trunk/SOURCE/microphysics_mod.f90
===================================================================
--- palm/trunk/SOURCE/microphysics_mod.f90	(revision 2292)
+++ palm/trunk/SOURCE/microphysics_mod.f90	(revision 2312)
@@ -20,26 +20,29 @@
 ! 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,5 +210,5 @@
            limiter_sedimentation, na_init, nc_const, sigma_gc,                 &
            ventilation_effect
-           
+
 
     INTERFACE microphysics_control
@@ -463,9 +466,9 @@
                    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,5 +481,5 @@
                       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,5 +524,5 @@
        USE particle_attributes,                                                &
           ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
-              s1, s2, s3, vanthoff
+              log_sigma, vanthoff
 
        IMPLICIT NONE
@@ -578,9 +581,9 @@
 !--             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,12 +595,12 @@
 !
 !--                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,5 +612,5 @@
 
 !
-!--                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,5 +619,5 @@
                    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,9 +625,9 @@
 
 !
-!--                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,5 +647,5 @@
 ! Description:
 ! ------------
-!> Calculate condensation rate for cloud water content (after Khairoutdinov and  
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and
 !> Kogan, 2000).
 !------------------------------------------------------------------------------!
@@ -724,5 +727,5 @@
 !--             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,5 +733,5 @@
 !
 !--             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,5 +744,5 @@
                    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,5 +894,5 @@
                                                               * 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,5 +1036,5 @@
 !
 !--                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,5 +1324,5 @@
 
        USE control_parameters,                                                 &
-           ONLY:  call_microphysics_at_all_substeps,                           & 
+           ONLY:  call_microphysics_at_all_substeps,                           &
                   intermediate_timestep_count, microphysics_morrison
 
@@ -1366,5 +1369,5 @@
                 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,5 +1452,5 @@
        USE surface_mod,                                                        &
            ONLY :  bc_h
-       
+
        IMPLICIT NONE
 
@@ -1455,5 +1458,5 @@
        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,7 +1531,7 @@
              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,5 +1542,5 @@
 !
 !--          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,5 +1610,5 @@
                 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,8 +1733,8 @@
        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,5 +1964,5 @@
        USE particle_attributes,                                                &
           ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
-              s1, s2, s3, vanthoff
+              log_sigma, vanthoff
 
        IMPLICIT NONE
@@ -2013,9 +2016,9 @@
 !--       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,8 +2030,8 @@
 !
 !--          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,5 +2039,5 @@
           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,5 +2049,5 @@
 
 !
-!--          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,5 +2056,5 @@
              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,12 +2062,12 @@
 
 !
-!--          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,5 +2079,5 @@
 ! Description:
 ! ------------
-!> Calculate condensation rate for cloud water content (after Khairoutdinov and  
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and
 !> Kogan, 2000).
 !------------------------------------------------------------------------------!
@@ -2154,5 +2157,5 @@
 !--       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,5 +2163,5 @@
 !
 !--       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,5 +2174,5 @@
              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,5 +2305,5 @@
              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,5 +2428,5 @@
 !
 !--          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,5 +2719,5 @@
           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,5 +2793,5 @@
        USE surface_mod,                                                        &
            ONLY :  bc_h
-       
+
        IMPLICIT NONE
 
@@ -2859,7 +2862,7 @@
        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,5 +2873,5 @@
 !
 !--    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,5 +2938,5 @@
           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,5 +3047,5 @@
        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
