Index: palm/trunk/SOURCE/bulk_cloud_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 4495)
+++ palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 4502)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
 ! 
@@ -116,18 +119,25 @@
                precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,   &
                qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,         &
-               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc, &
-               flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr, &
-               flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr, &
-               diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr, &
-               diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr
+               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc,            &
+               flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr,                     &
+               flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr,                     &
+               diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr,                     &
+               diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr,                     &
+               ni, ni_1, ni_2, ni_3, ni_p, tni_m,                              &
+               qi, qi_1, qi_2, qi_3, qi_p, tqi_m,                              &
+               flux_l_qi, flux_l_ni, flux_s_qi, flux_s_ni,                     &
+               diss_l_qi, diss_l_ni, diss_s_qi, diss_s_ni
+
 
     USE averaging,                                                             &
-        ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av
+        ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av, ni_av, qi_av
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,&
+        ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, magnus_ice,              &
+               molecular_weight_of_solute,                                     &
                molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,&
                exner_function, exner_function_invers, ideal_gas_law_rho,       &
-               ideal_gas_law_rho_pt, barometric_formula, rd_d_rv
+               ideal_gas_law_rho_pt, barometric_formula, rd_d_rv, l_s,         &
+               ls_d_cp
 
     USE control_parameters,                                                    &
@@ -143,7 +153,9 @@
                dt_3d, dt_do2d_xy, intermediate_timestep_count,                 &
                intermediate_timestep_count_max, large_scale_forcing,           &
-               lsf_surf, pt_surface, restart_data_format_output, rho_surface, surface_pressure,   &
+               lsf_surf, pt_surface, restart_data_format_output, rho_surface,  &
+               surface_pressure,                                               &
                time_do2d_xy, message_string, initializing_actions,             &
-               ws_scheme_sca, scalar_advec, timestep_scheme, tsc, loop_optimization
+               ws_scheme_sca, scalar_advec, timestep_scheme, tsc,              &
+               loop_optimization, simulated_time
 
     USE cpulog,                                                                &
@@ -167,16 +179,19 @@
         ONLY:  threads_per_task
 
-    USE restart_data_mpi_io_mod,                                                                   &
+    USE restart_data_mpi_io_mod,                                               &
         ONLY:  rrd_mpi_io, wrd_mpi_io
 
     USE statistics,                                                            &
         ONLY:  weight_pres, weight_substep, sums_wsncs_ws_l, sums_wsnrs_ws_l,  &
-               sums_wsqcs_ws_l, sums_wsqrs_ws_l
+               sums_wsqcs_ws_l, sums_wsqrs_ws_l,                               &
+               sums_wsqis_ws_l, sums_wsnis_ws_l
 
     USE surface_mod,                                                           &
         ONLY :  bc_h,                                                          &
                 surf_bulk_cloud_model,                                         &
-                surf_microphysics_morrison, surf_microphysics_seifert, &
-                surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
+                surf_microphysics_morrison, surf_microphysics_seifert,         &
+                surf_microphysics_ice_extension,                               &
+                surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
+                surf_usm_v
 
     IMPLICIT NONE
@@ -193,4 +208,5 @@
     LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
     LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
+    LOGICAL ::  ice_crystal_sedimentation = .FALSE.        !< flag for ice crystal sedimentation
     LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
     LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
@@ -198,9 +214,10 @@
 
     LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.      !< namelist parameter
+    LOGICAL ::  microphysics_ice_extension = .FALSE.              !< use ice microphysics scheme
     LOGICAL ::  microphysics_sat_adjust = .FALSE.                !< use saturation adjust bulk scheme?
     LOGICAL ::  microphysics_kessler = .FALSE.                   !< use kessler bulk scheme?
     LOGICAL ::  microphysics_morrison = .FALSE.                  !< use 2-moment Morrison (add. prog. eq. for nc and qc)
     LOGICAL ::  microphysics_seifert = .FALSE.                   !< use 2-moment Seifert and Beheng scheme
-    LOGICAL ::  microphysics_morrison_no_rain = .FALSE.          !< use 2-moment Morrison     
+    LOGICAL ::  microphysics_morrison_no_rain = .FALSE.          !< use 2-moment Morrison
     LOGICAL ::  precipitation = .FALSE.                          !< namelist parameter
 
@@ -240,5 +257,5 @@
     REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
     REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
-!    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05Âµm)
+    REAL(wp) ::  ximin = 4.42E-14_wp       !< minimum mass of ice crystal (kg) (D~10Âµm)
     REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1Âµm)
     REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
@@ -249,4 +266,5 @@
     REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
     REAL(wp) ::  dt_micro                        !< microphysics time step
+    REAL(wp) ::  in_init = 1000.0_wp             !< initial number of potential ice nucleii
     REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
     REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
@@ -254,10 +272,14 @@
     REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
     REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
-    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0;
+    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0
+    REAL(wp) ::  start_ice_microphysics = 0.0_wp !< time from which on ice microhysics are executed
 
     REAL(wp) ::  e_s     !< saturation water vapor pressure
+    REAL(wp) ::  e_si    !< saturation water vapor pressure over ice
     REAL(wp) ::  q_s     !< saturation mixing ratio
+    REAL(wp) ::  q_si    !< saturation mixing ratio over ice
     REAL(wp) ::  sat     !< supersaturation
-    REAL(wp) ::  t_l     !< actual temperature
+    REAL(wp) ::  sat_ice !< supersaturation over ice
+    REAL(wp) ::  t_l     !< liquid-(ice) water temperature
 
     SAVE
@@ -294,12 +316,15 @@
            dt_precipitation, &
            microphysics_morrison, &
-           microphysics_morrison_no_rain, &           
+           microphysics_morrison_no_rain, &
            microphysics_sat_adjust, &
            microphysics_seifert, &
+           microphysics_ice_extension, &
            na_init, &
            nc_const, &
            precipitation, &
-           sigma_gc
-
+           sigma_gc, &
+           start_ice_microphysics, &
+           ice_crystal_sedimentation, &
+           in_init
 
     INTERFACE bcm_parin
@@ -420,5 +445,9 @@
           nc_const, &
           sigma_bulk, &
-          ventilation_effect
+          ventilation_effect, &
+          ice_crystal_sedimentation, &
+          microphysics_ice_extension, &
+          start_ice_microphysics, &
+          in_init
 
        line = ' '
@@ -517,6 +546,6 @@
           microphysics_kessler    = .FALSE.
           microphysics_morrison   = .TRUE.
-          microphysics_morrison_no_rain = .TRUE.          
-          precipitation           = .FALSE.     
+          microphysics_morrison_no_rain = .TRUE.
+          precipitation           = .FALSE.
        ELSE
           message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
@@ -546,5 +575,5 @@
        surf_microphysics_morrison = microphysics_morrison
        surf_microphysics_seifert = microphysics_seifert
-
+       surf_microphysics_ice_extension = microphysics_ice_extension
 !
 !--    Check aerosol
@@ -592,4 +621,13 @@
              unit = '1/m3'
 
+          CASE ( 'ni' )
+             IF ( .NOT.  microphysics_ice_extension )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" ' //       &
+                                 'requires ' //                                &
+                                 'microphysics_ice_extension = ".TRUE."'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = '1/m3'
+
           CASE ( 'nr' )
              IF ( .NOT.  microphysics_seifert )  THEN
@@ -611,4 +649,13 @@
 
           CASE ( 'qc' )
+             unit = 'kg/kg'
+
+          CASE ( 'qi' )
+             IF ( .NOT.  microphysics_ice_extension ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" ' //       &
+                                 'requires ' //                                &
+                                 'microphysics_ice_extension = ".TRUE."'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
              unit = 'kg/kg'
 
@@ -698,4 +745,18 @@
              ENDIF
              pr_index = 123
+             dopr_index(var_count) = pr_index
+             dopr_unit     = '1/m3'
+             unit = dopr_unit
+             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'ni' )
+             IF ( .NOT.  microphysics_ice_extension )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) //          &
+                                 ' is not implemented for' //                  &
+                                 ' microphysics_ice_extension = ".F."'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ENDIF
+             pr_index = 124
              dopr_index(var_count) = pr_index
              dopr_unit     = '1/m3'
@@ -730,6 +791,21 @@
              unit = dopr_unit
              hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
+
           CASE ( 'qc' )
              pr_index = 75
+             dopr_index(var_count) = pr_index
+             dopr_unit     = 'kg/kg'
+             unit = dopr_unit
+             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'qi' )
+             IF ( .NOT.  microphysics_ice_extension )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(var_count) ) //          &
+                                 ' is not implemented for' //                  &
+                                 ' microphysics_ice_extension = ".F."'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ENDIF
+             pr_index = 125
              dopr_index(var_count) = pr_index
              dopr_unit     = 'kg/kg'
@@ -792,9 +868,9 @@
 !
 !--       3D-cloud drop water content, cloud drop concentration arrays
-          ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
+          ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
                     qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        ENDIF
@@ -803,14 +879,24 @@
 !
 !--       3D-rain water content, rain drop concentration arrays
-          ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
-                    qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
+          ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
                     qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        ENDIF
 
+       IF ( microphysics_ice_extension )  THEN
+!
+!--       3D-cloud drop water content, cloud drop concentration arrays
+          ALLOCATE( ni_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    ni_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    ni_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    qi_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    qi_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
+                    qi_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       ENDIF
+
        IF ( ws_scheme_sca )  THEN
-
           IF ( microphysics_morrison )  THEN
              ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
@@ -819,5 +905,4 @@
              sums_wsncs_ws_l = 0.0_wp
           ENDIF
-
           IF ( microphysics_seifert )  THEN
              ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
@@ -826,5 +911,10 @@
              sums_wsnrs_ws_l = 0.0_wp
           ENDIF
-
+          IF ( microphysics_ice_extension )  THEN
+             ALLOCATE( sums_wsqis_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             ALLOCATE( sums_wsnis_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_wsqis_ws_l = 0.0_wp
+             sums_wsnis_ws_l = 0.0_wp
+          ENDIF
        ENDIF
 
@@ -835,7 +925,5 @@
 !--    advection routines.
        IF ( loop_optimization /= 'vector' )  THEN
-
           IF ( ws_scheme_sca )  THEN
-
              IF ( microphysics_morrison )  THEN
                 ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
@@ -848,5 +936,4 @@
                           diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
-
              IF ( microphysics_seifert )  THEN
                 ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
@@ -859,7 +946,15 @@
                           diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
-
-          ENDIF
-
+             IF ( microphysics_ice_extension )  THEN
+                ALLOCATE( flux_s_qi(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_qi(nzb+1:nzt,0:threads_per_task-1),           &
+                          flux_s_ni(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_ni(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          flux_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
+             ENDIF
+          ENDIF
        ENDIF
 
@@ -878,4 +973,9 @@
           nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
        ENDIF
+       IF ( microphysics_ice_extension )  THEN
+          qi => qi_1;  qi_p  => qi_2;  tqi_m  => qi_3
+          ni => ni_1;  ni_p  => ni_2;  tni_m  => ni_3
+       ENDIF
+
 
 
@@ -919,4 +1019,14 @@
              ENDIF
 !
+!--          Initialize the remaining quantities
+             IF ( microphysics_ice_extension )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qi(:,j,i) = 0.0_wp
+                      ni(:,j,i) = 0.0_wp
+                   ENDDO
+                ENDDO
+             ENDIF
+!
 !--          Liquid water content and precipitation amount
 !--          are zero at beginning of the simulation
@@ -938,4 +1048,10 @@
                 qr_p  = qr
                 nr_p  = nr
+             ENDIF
+             IF ( microphysics_ice_extension )  THEN
+                tqi_m = 0.0_wp
+                tni_m = 0.0_wp
+                qi_p  = qi
+                ni_p  = ni
              ENDIF
           ENDIF ! Only if not read_restart_data
@@ -1080,4 +1196,8 @@
                 sums_wsnrs_ws_l = 0.0_wp
              ENDIF
+             IF ( microphysics_ice_extension )  THEN
+                sums_wsqis_ws_l = 0.0_wp
+                sums_wsnis_ws_l = 0.0_wp
+             ENDIF
 
           ENDIF
@@ -1096,5 +1216,4 @@
 !> Control of microphysics for grid points i,j
 !------------------------------------------------------------------------------!
-
     SUBROUTINE bcm_actions_ij( i, j, location )
 
@@ -1121,4 +1240,9 @@
                 sums_wsnrs_ws_l = 0.0_wp
              ENDIF
+             IF ( microphysics_ice_extension )  THEN
+                sums_wsqis_ws_l = 0.0_wp
+                sums_wsnis_ws_l = 0.0_wp
+             ENDIF
+
 
           ENDIF
@@ -1143,7 +1267,7 @@
        CALL cpu_log( log_point(51), 'microphysics', 'start' )
 
-       IF ( .NOT. microphysics_sat_adjust  .AND.         &
-            ( intermediate_timestep_count == 1  .OR.                              &
-              call_microphysics_at_all_substeps ) )                               &
+       IF ( .NOT. microphysics_sat_adjust  .AND.                               &
+            ( intermediate_timestep_count == 1  .OR.                           &
+              call_microphysics_at_all_substeps ) )                            &
        THEN
 
@@ -1151,11 +1275,15 @@
 !
 !--          Calculate vertical profile of the hydrostatic pressure (hyp)
-             hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
+             hyp    = barometric_formula(zu, pt_surface *                      &
+                      exner_function(surface_pressure * 100.0_wp),             &
+                      surface_pressure * 100.0_wp)
              d_exner = exner_function_invers(hyp)
              exner = 1.0_wp / exner_function_invers(hyp)
-             hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
+             hyrho  = ideal_gas_law_rho_pt(hyp, pt(:, nys, nxl) )
 !
 !--          Compute reference density
-             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
+             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,      &
+                           pt_surface *                                        &
+                           exner_function(surface_pressure * 100.0_wp))
           ENDIF
 
@@ -1178,5 +1306,5 @@
              CALL autoconversion_kessler
              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
-             
+
 !
 !--       Here the seifert beheng scheme is used. Cloud concentration is assumed to
@@ -1184,4 +1312,8 @@
           ELSEIF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )  THEN
              CALL adjust_cloud
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL ice_nucleation
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL ice_deposition
              CALL autoconversion
              CALL accretion
@@ -1190,15 +1322,25 @@
              CALL sedimentation_rain
              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
-             
-!
-!--       Here the morrison scheme is used. No rain processes are considered and qr and nr 
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL adjust_ice
+             IF ( ice_crystal_sedimentation .AND. microphysics_ice_extension &
+                  .AND. simulated_time > start_ice_microphysics ) CALL sedimentation_ice
+
+!
+!--       Here the morrison scheme is used. No rain processes are considered and qr and nr
 !--       are not allocated
           ELSEIF ( microphysics_morrison_no_rain  .AND.  .NOT. microphysics_seifert )  THEN
              CALL activation
              CALL condensation
-             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud   
-             
-!
-!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are 
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL adjust_ice
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL ice_nucleation
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL ice_deposition
+             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
+
+!
+!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
 !--       included
           ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
@@ -1206,4 +1348,10 @@
              CALL activation
              CALL condensation
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL adjust_ice
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL ice_nucleation
+             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
+             CALL ice_deposition
              CALL autoconversion
              CALL accretion
@@ -1211,5 +1359,5 @@
              CALL evaporation_rain
              CALL sedimentation_rain
-             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud            
+             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
 
           ENDIF
@@ -1229,5 +1377,4 @@
 !> Control of microphysics for grid points i,j
 !------------------------------------------------------------------------------!
-
     SUBROUTINE bcm_non_advective_processes_ij( i, j )
 
@@ -1236,7 +1383,7 @@
        INTEGER(iwp) ::  j                 !<
 
-       IF ( .NOT. microphysics_sat_adjust  .AND.         &
-            ( intermediate_timestep_count == 1  .OR.                              &
-              call_microphysics_at_all_substeps ) )                               &
+       IF ( .NOT. microphysics_sat_adjust  .AND.                               &
+            ( intermediate_timestep_count == 1  .OR.                           &
+              call_microphysics_at_all_substeps ) )                            &
        THEN
 
@@ -1244,11 +1391,15 @@
 !
 !--          Calculate vertical profile of the hydrostatic pressure (hyp)
-             hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
+             hyp    = barometric_formula(zu, pt_surface *                      &
+                      exner_function(surface_pressure * 100.0_wp),             &
+                      surface_pressure * 100.0_wp)
              d_exner = exner_function_invers(hyp)
              exner = 1.0_wp / exner_function_invers(hyp)
-             hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
+             hyrho  = ideal_gas_law_rho_pt(hyp, pt(:, nys, nxl) )
 !
 !--          Compute reference density
-             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
+             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,      &
+                           pt_surface *                                        &
+                           exner_function(surface_pressure * 100.0_wp))
           ENDIF
 
@@ -1273,7 +1424,11 @@
 !
 !--       Here the seifert beheng scheme is used. Cloud concentration is assumed to
-!--       a constant value an qc a diagnostic value.             
+!--       a constant value an qc a diagnostic value.
           ELSEIF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )  THEN
              CALL adjust_cloud_ij( i,j )
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL ice_nucleation_ij( i,j )
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL ice_deposition_ij( i,j )
              CALL autoconversion_ij( i,j )
              CALL accretion_ij( i,j )
@@ -1282,20 +1437,35 @@
              CALL sedimentation_rain_ij( i,j )
              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
-             
-!
-!--       Here the morrison scheme is used. No rain processes are considered and qr and nr 
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL adjust_ice_ij ( i,j )
+             IF ( ice_crystal_sedimentation .AND. microphysics_ice_extension &
+                  .AND. simulated_time > start_ice_microphysics )  CALL sedimentation_ice_ij ( i,j )
+!
+!--       Here the morrison scheme is used. No rain processes are considered and qr and nr
 !--       are not allocated
           ELSEIF ( microphysics_morrison_no_rain  .AND.  .NOT. microphysics_seifert )  THEN
              CALL activation_ij( i,j )
              CALL condensation_ij( i,j )
-             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j ) 
-             
-!
-!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are 
-!--       included             
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL adjust_ice_ij ( i,j )
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL ice_nucleation_ij( i,j )
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL ice_deposition_ij( i,j )
+             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
+
+!
+!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
+!--       included
           ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
              CALL adjust_cloud_ij( i,j )
              CALL activation_ij( i,j )
              CALL condensation_ij( i,j )
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL adjust_ice_ij ( i,j )
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL ice_nucleation_ij( i,j )
+             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
+                 CALL ice_deposition_ij( i,j )
              CALL autoconversion_ij( i,j )
              CALL accretion_ij( i,j )
@@ -1303,5 +1473,5 @@
              CALL evaporation_rain_ij( i,j )
              CALL sedimentation_rain_ij( i,j )
-             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )           
+             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
 
           ENDIF
@@ -1331,5 +1501,5 @@
           IF ( microphysics_morrison )  THEN
              CALL exchange_horiz( nc, nbgp )
-             CALL exchange_horiz( qc, nbgp )          
+             CALL exchange_horiz( qc, nbgp )
           ENDIF
           IF ( microphysics_seifert ) THEN
@@ -1337,11 +1507,15 @@
              CALL exchange_horiz( nr, nbgp )
           ENDIF
+          IF ( microphysics_ice_extension ) THEN
+             CALL exchange_horiz( qi, nbgp )
+             CALL exchange_horiz( ni, nbgp )
+          ENDIF
           CALL exchange_horiz( q, nbgp )
-          CALL exchange_horiz( pt, nbgp )          
+          CALL exchange_horiz( pt, nbgp )
        ENDIF
 
 
     END SUBROUTINE bcm_exchange_horiz
-    
+
 
 
@@ -1425,5 +1599,5 @@
                                              )                                 &
                                     * MERGE( 1.0_wp, 0.0_wp,                   &
-                                       BTEST( wall_flags_total_0(k,j,i), 0 )   &
+                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
                                           )
                    IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
@@ -1517,5 +1691,5 @@
                                              )                                 &
                                    * MERGE( 1.0_wp, 0.0_wp,                    &
-                                       BTEST( wall_flags_total_0(k,j,i), 0 )   &
+                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
                                           )
                    IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
@@ -1551,4 +1725,197 @@
 
        ENDIF
+
+!
+!--    If required, calculate prognostic equations for ice crystal content
+!--    and ice crystal concentration
+       IF ( microphysics_ice_extension )  THEN
+
+          CALL cpu_log( log_point(70), 'qi-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for ice crystal content
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( qi, 'qi' )
+
+          ENDIF
+
+!
+!--       qi-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( advc_flags_s, qi, 'qi',                    &
+                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
+                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
+                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
+                                    bc_dirichlet_s  .OR.  bc_radiation_s )
+                ELSE
+                   CALL advec_s_pw( qi )
+                ENDIF
+             ELSE
+                CALL advec_s_up( qi )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( qi,                                                &
+                            surf_def_h(0)%qisws, surf_def_h(1)%qisws,          &
+                            surf_def_h(2)%qisws,                               &
+                            surf_lsm_h%qisws,    surf_usm_h%qisws,             &
+                            surf_def_v(0)%qisws, surf_def_v(1)%qisws,          &
+                            surf_def_v(2)%qisws, surf_def_v(3)%qisws,          &
+                            surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,          &
+                            surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,          &
+                            surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,          &
+                            surf_usm_v(2)%qisws, surf_usm_v(3)%qisws )
+
+!
+!--       Prognostic equation for ice crystal mixing ratio
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                !following directive is required to vectorize on Intel19
+                !DIR$ IVDEP
+                DO  k = nzb+1, nzt
+                   qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
+                                                      tsc(3) * tqi_m(k,j,i) )  &
+                                                    - tsc(5) * rdf_sc(k) *     &
+                                                               qi(k,j,i)       &
+                                             )                                 &
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
+                                          )
+                   IF ( qi_p(k,j,i) < 0.0_wp )  qi_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tqi_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
+                                         + 5.3125_wp * tqi_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(70), 'qi-equation', 'stop' )
+
+          CALL cpu_log( log_point(69), 'ni-equation', 'start' )
+!
+!--       Calculate prognostic equation for ice crystal concentration
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( ni, 'ni' )
+
+          ENDIF
+
+!
+!--       ni-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( advc_flags_s, ni, 'ni',                    &
+                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
+                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
+                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
+                                    bc_dirichlet_s  .OR.  bc_radiation_s )
+                ELSE
+                   CALL advec_s_pw( ni )
+                ENDIF
+             ELSE
+                CALL advec_s_up( ni )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( ni,                                                &
+                            surf_def_h(0)%nisws, surf_def_h(1)%nisws,          &
+                            surf_def_h(2)%nisws,                               &
+                            surf_lsm_h%nisws,    surf_usm_h%nisws,             &
+                            surf_def_v(0)%nisws, surf_def_v(1)%nisws,          &
+                            surf_def_v(2)%nisws, surf_def_v(3)%nisws,          &
+                            surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,          &
+                            surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,          &
+                            surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,          &
+                            surf_usm_v(2)%nisws, surf_usm_v(3)%nisws )
+
+!
+!--       Prognostic equation for ice crystal concentration
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                !following directive is required to vectorize on Intel19
+                !DIR$ IVDEP
+                DO  k = nzb+1, nzt
+                   ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
+                                                      tsc(3) * tni_m(k,j,i) )  &
+                                                    - tsc(5) * rdf_sc(k) *     &
+                                                               ni(k,j,i)       &
+                                             )                                 &
+                                   * MERGE( 1.0_wp, 0.0_wp,                    &
+                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
+                                          )
+                   IF ( ni_p(k,j,i) < 0.0_wp )  ni_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tni_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tni_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
+                                         + 5.3125_wp * tni_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(69), 'ni-equation', 'stop' )
+
+       ENDIF
+
 !
 !--    If required, calculate prognostic equations for rain water content
@@ -1616,5 +1983,5 @@
                                              )                                 &
                                     * MERGE( 1.0_wp, 0.0_wp,                   &
-                                       BTEST( wall_flags_total_0(k,j,i), 0 )   &
+                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
                                           )
                    IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
@@ -1708,5 +2075,5 @@
                                              )                                 &
                                    * MERGE( 1.0_wp, 0.0_wp,                    &
-                                       BTEST( wall_flags_total_0(k,j,i), 0 )   &
+                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
                                           )
                    IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
@@ -1751,5 +2118,4 @@
 !> Control of microphysics for grid points i,j
 !------------------------------------------------------------------------------!
-
     SUBROUTINE bcm_prognostic_equations_ij( i, j, i_omp_start, tn )
 
@@ -1805,5 +2171,5 @@
                                        )                                 &
                                  * MERGE( 1.0_wp, 0.0_wp,                &
-                                    BTEST( wall_flags_total_0(k,j,i), 0 )&
+                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
                                         )
              IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
@@ -1864,5 +2230,5 @@
                                        )                                 &
                                  * MERGE( 1.0_wp, 0.0_wp,                &
-                                    BTEST( wall_flags_total_0(k,j,i), 0 )&
+                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
                                         )
              IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
@@ -1885,4 +2251,130 @@
 
        ENDIF
+
+!
+!--    If required, calculate prognostic equations for ice crystal mixing ratio
+!--    and ice crystal concentration
+       IF ( microphysics_ice_extension )  THEN
+!
+!--       Calculate prognostic equation for ice crystal mixing ratio
+          tend(:,j,i) = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' ) &
+          THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( advc_flags_s, i, j, qi, 'qi', flux_s_qi,      &
+                                 diss_s_qi, flux_l_qi, diss_l_qi,              &
+                                 i_omp_start, tn,                              &
+                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
+                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
+                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
+                                 bc_dirichlet_s  .OR.  bc_radiation_s )
+             ELSE
+                CALL advec_s_pw( i, j, qi )
+             ENDIF
+          ELSE
+             CALL advec_s_up( i, j, qi )
+          ENDIF
+          CALL diffusion_s( i, j, qi,                                   &
+                            surf_def_h(0)%qisws, surf_def_h(1)%qisws,   &
+                            surf_def_h(2)%qisws,                        &
+                            surf_lsm_h%qisws,    surf_usm_h%qisws,      &
+                            surf_def_v(0)%qisws, surf_def_v(1)%qisws,   &
+                            surf_def_v(2)%qisws, surf_def_v(3)%qisws,   &
+                            surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,   &
+                            surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,   &
+                            surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,   &
+                            surf_usm_v(2)%qisws, surf_usm_v(3)%qisws )
+
+!
+!--       Prognostic equation for ice crystal mixing ratio
+          DO  k = nzb+1, nzt
+             qi_p(k,j,i) = qi(k,j,i) + ( dt_3d *                         &
+                                                ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tqi_m(k,j,i) )&
+                                                - tsc(5) * rdf_sc(k)     &
+                                                         * qi(k,j,i)     &
+                                       )                                 &
+                                 * MERGE( 1.0_wp, 0.0_wp,                &
+                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
+                                        )
+             IF ( qi_p(k,j,i) < 0.0_wp )  qi_p(k,j,i) = 0.0_wp
+          ENDDO
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb+1, nzt
+                   tqi_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb+1, nzt
+                   tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                     5.3125_wp * tqi_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+!
+!--       Calculate prognostic equation for ice crystal concentration.
+          tend(:,j,i) = 0.0_wp
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( advc_flags_s, i, j, ni, 'ni', flux_s_ni,      &
+                                 diss_s_ni, flux_l_ni, diss_l_ni,              &
+                                 i_omp_start, tn,                              &
+                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
+                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
+                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
+                                 bc_dirichlet_s  .OR.  bc_radiation_s )
+             ELSE
+                CALL advec_s_pw( i, j, ni )
+             ENDIF
+          ELSE
+             CALL advec_s_up( i, j, ni )
+          ENDIF
+          CALL diffusion_s( i, j, ni,                                    &
+                            surf_def_h(0)%nisws, surf_def_h(1)%nisws,    &
+                            surf_def_h(2)%nisws,                         &
+                            surf_lsm_h%nisws,    surf_usm_h%nisws,       &
+                            surf_def_v(0)%nisws, surf_def_v(1)%nisws,    &
+                            surf_def_v(2)%nisws, surf_def_v(3)%nisws,    &
+                            surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,    &
+                            surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,    &
+                            surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,    &
+                            surf_usm_v(2)%nisws, surf_usm_v(3)%nisws )
+
+!
+!--       Prognostic equation for ice crystal concentration
+          DO  k = nzb+1, nzt
+             ni_p(k,j,i) = ni(k,j,i) + ( dt_3d *                         &
+                                                ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tni_m(k,j,i) )&
+                                                - tsc(5) * rdf_sc(k)     &
+                                                         * ni(k,j,i)     &
+                                       )                                 &
+                                 * MERGE( 1.0_wp, 0.0_wp,                &
+                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
+                                        )
+             IF ( ni_p(k,j,i) < 0.0_wp )  ni_p(k,j,i) = 0.0_wp
+          ENDDO
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  k = nzb+1, nzt
+                   tni_m(k,j,i) = tend(k,j,i)
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  k = nzb+1, nzt
+                   tni_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                     5.3125_wp * tni_m(k,j,i)
+                ENDDO
+             ENDIF
+          ENDIF
+
+       ENDIF
+
 !
 !--    If required, calculate prognostic equations for rain water content
@@ -1929,5 +2421,5 @@
                                        )                                 &
                                  * MERGE( 1.0_wp, 0.0_wp,                &
-                                    BTEST( wall_flags_total_0(k,j,i), 0 )&
+                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
                                         )
              IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
@@ -1988,5 +2480,5 @@
                                        )                                 &
                                  * MERGE( 1.0_wp, 0.0_wp,                &
-                                    BTEST( wall_flags_total_0(k,j,i), 0 )&
+                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
                                         )
              IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
@@ -2038,4 +2530,8 @@
                    nr => nr_1;    nr_p => nr_2
                 ENDIF
+                IF ( microphysics_ice_extension )  THEN
+                   qi => qi_1;    qi_p => qi_2
+                   ni => ni_1;    ni_p => ni_2
+                ENDIF
 
              CASE ( 1 )
@@ -2049,4 +2545,9 @@
                    nr => nr_2;    nr_p => nr_1
                 ENDIF
+                IF ( microphysics_ice_extension )  THEN
+                   qi => qi_2;    qi_p => qi_1
+                   ni => ni_2;    ni_p => ni_1
+                ENDIF
+
 
           END SELECT
@@ -2093,4 +2594,28 @@
        ENDIF
 
+       IF ( microphysics_ice_extension )  THEN
+!
+!--       Surface conditions ice crysral (Dirichlet)
+!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
+!--       the k coordinate belongs to the atmospheric grid point, therefore, set
+!--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
+          DO  l = 0, 1
+          !$OMP PARALLEL DO PRIVATE( i, j, k )
+             DO  m = 1, bc_h(l)%ns
+                i = bc_h(l)%i(m)
+                j = bc_h(l)%j(m)
+                k = bc_h(l)%k(m)
+                qi_p(k+bc_h(l)%koff,j,i) = 0.0_wp
+                ni_p(k+bc_h(l)%koff,j,i) = 0.0_wp
+             ENDDO
+          ENDDO
+!
+!--       Top boundary condition for ice crystal (Dirichlet)
+          qi_p(nzt+1,:,:) = 0.0_wp
+          ni_p(nzt+1,:,:) = 0.0_wp
+
+       ENDIF
+
+
        IF ( microphysics_seifert )  THEN
 !
@@ -2129,4 +2654,8 @@
              nr_p(:,nys-1,:) = nr_p(:,nys,:)
           ENDIF
+          IF ( microphysics_ice_extension )  THEN
+             qi_p(:,nys-1,:) = qi_p(:,nys,:)
+             ni_p(:,nys-1,:) = ni_p(:,nys,:)
+          ENDIF
        ELSEIF ( bc_radiation_n )  THEN
           IF ( microphysics_morrison )  THEN
@@ -2138,4 +2667,8 @@
              nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
           ENDIF
+          IF ( microphysics_ice_extension )  THEN
+             qi_p(:,nyn+1,:) = qi_p(:,nyn,:)
+             ni_p(:,nyn+1,:) = ni_p(:,nyn,:)
+          ENDIF
        ELSEIF ( bc_radiation_l )  THEN
           IF ( microphysics_morrison )  THEN
@@ -2147,4 +2680,8 @@
              nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
           ENDIF
+          IF ( microphysics_ice_extension )  THEN
+             qi_p(:,:,nxl-1) = qi_p(:,:,nxl)
+             ni_p(:,:,nxl-1) = ni_p(:,:,nxl)
+          ENDIF
        ELSEIF ( bc_radiation_r )  THEN
           IF ( microphysics_morrison )  THEN
@@ -2155,4 +2692,8 @@
              qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
              nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
+          ENDIF
+          IF ( microphysics_ice_extension )  THEN
+             qi_p(:,:,nxr+1) = qi_p(:,:,nxr)
+             ni_p(:,:,nxr+1) = ni_p(:,:,nxr)
           ENDIF
        ENDIF
@@ -2433,4 +2974,16 @@
              ENDIF
              to_be_resorted => nc_av
+          ENDIF
+          IF ( mode == 'xy' ) grid = 'zu'
+
+       CASE ( 'ni_xy', 'ni_xz', 'ni_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ni
+          ELSE
+             IF ( .NOT. ALLOCATED( ni_av ) ) THEN
+                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ni_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ni_av
           ENDIF
           IF ( mode == 'xy' ) grid = 'zu'
@@ -2499,4 +3052,16 @@
           IF ( mode == 'xy' ) grid = 'zu'
 
+       CASE ( 'qi_xy', 'qi_xz', 'qi_yz' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => qi
+          ELSE
+             IF ( .NOT. ALLOCATED( qi_av ) ) THEN
+                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                qi_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => qi_av
+          ENDIF
+          IF ( mode == 'xy' ) grid = 'zu'
+
        CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
           IF ( av == 0 )  THEN
@@ -2534,7 +3099,7 @@
              DO  k = nzb_do, nzt_do
                 local_pf(i,j,k) = MERGE( &
-                                     to_be_resorted(k,j,i),                      &
-                                     REAL( fill_value, KIND = wp ),              &
-                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) &
+                                     to_be_resorted(k,j,i),                    &
+                                     REAL( fill_value, KIND = wp ),            &
+                                     BTEST( wall_flags_total_0(k,j,i), flag_nr )     &
                                   )
              ENDDO
@@ -2595,4 +3160,15 @@
              ENDIF
              to_be_resorted => nc_av
+          ENDIF
+
+       CASE ( 'ni' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => ni
+          ELSE
+             IF ( .NOT. ALLOCATED( ni_av ) ) THEN
+                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ni_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => ni_av
           ENDIF
 
@@ -2643,4 +3219,15 @@
           ENDIF
 
+       CASE ( 'qi' )
+          IF ( av == 0 )  THEN
+             to_be_resorted => qi
+          ELSE
+             IF ( .NOT. ALLOCATED( qi_av ) ) THEN
+                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                qi_av = REAL( fill_value, KIND = wp )
+             ENDIF
+             to_be_resorted => qi_av
+          ENDIF
+
        CASE ( 'ql' )
           IF ( av == 0 )  THEN
@@ -2675,8 +3262,8 @@
           DO  j = nys, nyn
              DO  k = nzb_do, nzt_do
-                local_pf(i,j,k) = MERGE(                                         &
-                                     to_be_resorted(k,j,i),                      &
-                                     REAL( fill_value, KIND = wp ),              &
-                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) &
+                local_pf(i,j,k) = MERGE(                                       &
+                                     to_be_resorted(k,j,i),                    &
+                                     REAL( fill_value, KIND = wp ),            &
+                                     BTEST( wall_flags_total_0(k,j,i), flag_nr )     &
                                   )
              ENDDO
@@ -2750,4 +3337,16 @@
           CASE ( 'curvature_solution_effects_bulk' )
              READ ( 13 )  curvature_solution_effects_bulk
+
+          CASE ( 'microphysics_ice_extension' )
+             READ ( 13 )  microphysics_ice_extension
+
+          CASE ( 'ice_crystal_sedimentation' )
+             READ ( 13 )  ice_crystal_sedimentation
+
+          CASE ( 'in_init' )
+             READ ( 13 )  in_init
+
+          CASE ( 'start_ice_microphysics' )
+             READ ( 13 )  start_ice_microphysics
 
           CASE DEFAULT
@@ -2782,4 +3381,10 @@
        CALL rrd_mpi_io( 'aerosol_bulk', aerosol_bulk )
        CALL rrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk )
+       CALL rrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics )
+       CALL rrd_mpi_io( 'microphysics_ice_extension', microphysics_ice_extension )
+       CALL rrd_mpi_io( 'in_init', in_init )
+       CALL rrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation )
+
+
 
     END SUBROUTINE bcm_rrd_global_mpi
@@ -2846,4 +3451,17 @@
                 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+          CASE ( 'ni' )
+             IF ( k == 1 )  READ ( 13 )  tmp_3d
+             ni(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
+                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+          CASE ( 'ni_av' )
+             IF ( .NOT. ALLOCATED( ni_av ) )  THEN
+                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             IF ( k == 1 )  READ ( 13 )  tmp_3d
+             ni_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
+                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
           CASE ( 'nr' )
              IF ( k == 1 )  READ ( 13 )  tmp_3d
@@ -2886,4 +3504,17 @@
                 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+          CASE ( 'qi' )
+             IF ( k == 1 )  READ ( 13 )  tmp_3d
+             qi(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
+                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+          CASE ( 'qi_av' )
+             IF ( .NOT. ALLOCATED( qi_av ) )  THEN
+                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             IF ( k == 1 )  READ ( 13 )  tmp_3d
+             qi_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
+                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
           CASE ( 'qc_av' )
              IF ( .NOT. ALLOCATED( qc_av ) )  THEN
@@ -2984,4 +3615,16 @@
           CALL wrd_write_string( 'curvature_solution_effects_bulk' )
           WRITE ( 14 )  curvature_solution_effects_bulk
+
+          CALL wrd_write_string( 'start_ice_microphysics' )
+          WRITE ( 14 )  start_ice_microphysics
+
+          CALL wrd_write_string( 'microphysics_ice_extension' )
+          WRITE ( 14 )  microphysics_ice_extension
+
+          CALL wrd_write_string( 'in_init' )
+          WRITE ( 14 )  in_init
+
+          CALL wrd_write_string( 'ice_crystal_sedimentation' )
+          WRITE ( 14 )  ice_crystal_sedimentation
 
        ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' )  THEN
@@ -3001,4 +3644,8 @@
           CALL wrd_mpi_io( 'aerosol_bulk', aerosol_bulk )
           CALL wrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk )
+          CALL wrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics )
+          CALL wrd_mpi_io( 'microphysics_ice_extension', microphysics_ice_extension )
+          CALL wrd_mpi_io( 'in_init', in_init )
+          CALL wrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation )
 
        ENDIF
@@ -3062,4 +3709,25 @@
 
           ENDIF
+
+          IF ( microphysics_ice_extension )  THEN
+
+             CALL wrd_write_string( 'ni' )
+             WRITE ( 14 )  ni
+
+             IF ( ALLOCATED( ni_av ) )  THEN
+                CALL wrd_write_string( 'ni_av' )
+                WRITE ( 14 )  ni_av
+             ENDIF
+
+             CALL wrd_write_string( 'qi' )
+             WRITE ( 14 )  qi
+
+             IF ( ALLOCATED( qi_av ) )  THEN
+                CALL wrd_write_string( 'qi_av' )
+                WRITE ( 14 )  qi_av
+             ENDIF
+
+          ENDIF
+
 
           IF ( microphysics_seifert )  THEN
@@ -3104,4 +3772,10 @@
              IF ( ALLOCATED( qr_av ) )  CALL wrd_mpi_io( 'qr_av', qr_av )
           ENDIF
+          IF ( microphysics_ice_extension )  THEN
+             CALL wrd_mpi_io( 'ni', ni )
+             IF ( ALLOCATED( ni_av ) )  CALL wrd_mpi_io( 'ni_av', ni_av )
+             CALL wrd_mpi_io( 'qi', qi )
+             IF ( ALLOCATED( qi_av ) )  CALL wrd_mpi_io( 'qi_av', qi_av )
+          ENDIF
 
        ENDIF
@@ -3134,13 +3808,14 @@
 !--             Predetermine flag to mask topography
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-
-                IF ( qr(k,j,i) <= eps_sb )  THEN
-                   qr(k,j,i) = 0.0_wp
-                   nr(k,j,i) = 0.0_wp
-                ELSE
-                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
-                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
-                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
-                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
+                IF ( .NOT. microphysics_morrison_no_rain )  THEN
+                   IF ( qr(k,j,i) <= eps_sb )  THEN
+                      qr(k,j,i) = 0.0_wp
+                      nr(k,j,i) = 0.0_wp
+                   ELSE
+                      IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
+                         nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
+                      ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
+                         nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
+                      ENDIF
                    ENDIF
                 ENDIF
@@ -3189,18 +3864,19 @@
           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
 
-          IF ( qr(k,j,i) <= eps_sb )  THEN
-             qr(k,j,i) = 0.0_wp
-             nr(k,j,i) = 0.0_wp
-          ELSE
-!
-!--          Adjust number of raindrops to avoid nonlinear effects in
-!--          sedimentation and evaporation of rain drops due to too small or
-!--          too big weights of rain drops (Stevens and Seifert, 2008).
-             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
-                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
-             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
-                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
-             ENDIF
-
+          IF ( .NOT.  microphysics_morrison_no_rain )  THEN
+             IF ( qr(k,j,i) <= eps_sb )  THEN
+                qr(k,j,i) = 0.0_wp
+                nr(k,j,i) = 0.0_wp
+             ELSE
+!
+!--             Adjust number of raindrops to avoid nonlinear effects in
+!--             sedimentation and evaporation of rain drops due to too small or
+!--             too big weights of rain drops (Stevens and Seifert, 2008).
+                IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
+                   nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
+                ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
+                   nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
+                ENDIF
+             ENDIF
           ENDIF
 
@@ -3219,4 +3895,81 @@
 
     END SUBROUTINE adjust_cloud_ij
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Adjust number of ice crystal to avoid nonlinear effects in sedimentation and
+!> evaporation of ice crytals due to too small or too big weights
+!> of ice crytals (Stevens and Seifert, 2008).
+!------------------------------------------------------------------------------!
+    SUBROUTINE adjust_ice
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp) ::  flag                  !< flag to indicate first grid level above
+
+       CALL cpu_log( log_point_s(97), 'adjust_ice', 'start' )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+                IF ( qi(k,j,i) <= ximin )  THEN
+                   qi(k,j,i) = 0.0_wp
+                   ni(k,j,i) = 0.0_wp
+                ELSE
+                   IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) )  THEN
+                      ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag
+                   ENDIF
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(97), 'adjust_ice', 'stop' )
+
+    END SUBROUTINE adjust_ice
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Adjust number of of ice crystal to avoid nonlinear effects in
+!> sedimentation and evaporation of ice crystals due to too small or
+!> too big weights of ice crytals (Stevens and Seifert, 2008).
+!> The same procedure is applied to cloud droplets if they are determined
+!> prognostically. Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE adjust_ice_ij( i, j )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+          IF ( qi(k,j,i) <= ximin )  THEN
+             qi(k,j,i) = 0.0_wp
+             ni(k,j,i) = 0.0_wp
+          ELSE
+             IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) )  THEN
+                ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag
+             ENDIF
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE adjust_ice_ij
+
 
 !------------------------------------------------------------------------------!
@@ -3423,4 +4176,136 @@
     END SUBROUTINE activation_ij
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate ice nucleation by applying the deposition-condensation formula as
+!> given by Meyers et al 1992 and as described in Seifert and Beheng 2006
+!------------------------------------------------------------------------------!
+    SUBROUTINE ice_nucleation
+
+       INTEGER(iwp) ::  i             !< loop index
+       INTEGER(iwp) ::  j             !< loop index
+       INTEGER(iwp) ::  k             !< loop index
+
+       LOGICAL :: isdac = .TRUE.
+
+       REAL(wp) ::  a_m92 = -0.639_wp !< parameter for nucleation
+       REAL(wp) ::  b_m92 = 12.96_wp  !< parameter for nucleation
+       REAL(wp) ::  flag              !< flag to indicate first grid level
+       REAL(wp) ::  n_in              !< number of ice nucleii
+       REAL(wp) ::  nucle = 0.0_wp    !< nucleation rate
+
+       CALL cpu_log( log_point_s(93), 'ice_nucleation', 'start' )
+
+       IF  ( isdac )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+!
+!--                Predetermine flag to mask topography
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--                Call calculation of supersaturation located in subroutine
+                   CALL supersaturation_ice ( i, j, k )
+                   nucle = 0.0_wp
+                   !IF ( zu(k) >= 1500.0_wp ) CYCLE
+                   IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp  ) THEN
+!
+!--                   Calculate ice nucleation
+                      nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp )
+                      !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin
+                      ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init)
+                   ENDIF
+
+                ENDDO
+             ENDDO
+          ENDDO
+       ELSE
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+!
+!--                Predetermine flag to mask topography
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--                Call calculation of supersaturation located in subroutine
+                   CALL supersaturation_ice ( i, j, k )
+                   nucle = 0.0_wp
+                   IF ( sat_ice > 0.0 ) THEN
+!
+!--                   Calculate ice nucleation
+                      n_in = in_init * EXP( a_m92 + b_m92 * sat_ice )
+                      nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp )
+                   ENDIF
+                   ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp )
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+
+       CALL cpu_log( log_point_s(93), 'ice_nucleation', 'stop' )
+
+    END SUBROUTINE ice_nucleation
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate ice nucleation by applying the deposition-condensation formula as
+!> given by Meyers et al 1992 and as described in Seifert and Beheng 2006
+!------------------------------------------------------------------------------!
+    SUBROUTINE ice_nucleation_ij( i, j )
+
+       INTEGER(iwp) ::  i             !< loop index
+       INTEGER(iwp) ::  j             !< loop index
+       INTEGER(iwp) ::  k             !< loop index
+
+       LOGICAL :: isdac = .TRUE.
+
+       REAL(wp) ::  a_m92 = -0.639_wp !< parameter for nucleation
+       REAL(wp) ::  b_m92 = 12.96_wp  !< parameter for nucleation
+       REAL(wp) ::  flag              !< flag to indicate first grid level
+       REAL(wp) ::  n_in              !< number of ice nucleii
+       REAL(wp) ::  nucle = 0.0_wp    !< nucleation rate
+
+       IF ( isdac )  THEN
+          DO  k = nzb+1, nzt
+!
+!--          Predetermine flag to mask topography
+             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--          Call calculation of supersaturation located in subroutine
+             CALL supersaturation_ice ( i, j, k )
+             nucle = 0.0_wp
+             !IF ( zu(k) >= 1500.0_wp ) CYCLE
+             IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp ) THEN
+!
+!--             Calculate ice nucleation
+                nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp )
+                !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin
+                ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init)
+             ENDIF
+          ENDDO
+       ELSE
+          DO  k = nzb+1, nzt
+!
+!--          Predetermine flag to mask topography
+             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--          Call calculation of supersaturation located in subroutine
+             CALL supersaturation_ice ( i, j, k )
+             nucle = 0.0_wp
+             IF ( sat_ice > 0.0 ) THEN
+!
+!--             Calculate ice nucleation
+                n_in = in_init * EXP( a_m92 + b_m92 * sat_ice )
+                nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp )
+             ENDIF
+             ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp )
+          ENDDO
+       ENDIF
+
+
+    END SUBROUTINE ice_nucleation_ij
 
 !------------------------------------------------------------------------------!
@@ -3434,18 +4319,17 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp)     ::  cond              !<
-       REAL(wp)     ::  cond_max          !<
-       REAL(wp)     ::  dc                !<
-       REAL(wp)     ::  evap              !<
-       REAL(wp)     ::  g_fac             !<
-       REAL(wp)     ::  nc_0              !<
-       REAL(wp)     ::  temp              !<
-       REAL(wp)     ::  xc                !<
-
-       REAL(wp) ::  flag                  !< flag to indicate first grid level above
+       INTEGER(iwp) ::  i                 !< loop index
+       INTEGER(iwp) ::  j                 !< loop index
+       INTEGER(iwp) ::  k                 !< loop index
+
+       REAL(wp)     ::  cond              !< condensation rate
+       REAL(wp)     ::  cond_max          !< maximum condensation rate
+       REAL(wp)     ::  dc                !< weight avageraed diameter
+       REAL(wp)     ::  evap              !< evaporation rate
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
+       REAL(wp)     ::  g_fac             !< factor 1 / Fk + Fd
+       REAL(wp)     ::  nc_0              !< integral diameter
+       REAL(wp)     ::  temp              !< actual temperature
+       REAL(wp)     ::  xc                !< mean mass droplets
 
        CALL cpu_log( log_point_s(66), 'condensation', 'start' )
@@ -3461,10 +4345,11 @@
                 CALL supersaturation ( i, j, k )
 !
-!--             Actual temperature:
-                IF ( microphysics_seifert ) THEN
-                   temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
-                ELSEIF ( microphysics_morrison_no_rain ) THEN
-                   temp = t_l + lv_d_cp * qc(k,j,i) 
-                ENDIF  
+!--             Actual temperature, t_l is calculated directly before
+!--             in supersaturation
+                IF ( microphysics_ice_extension ) THEN
+                   temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
+                ELSE
+                   temp = t_l + lv_d_cp * ql(k,j,i)
+                ENDIF
 
                 g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
@@ -3525,18 +4410,17 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                 !<
-       INTEGER(iwp) ::  j                 !<
-       INTEGER(iwp) ::  k                 !<
-
-       REAL(wp)     ::  cond              !<
-       REAL(wp)     ::  cond_max          !<
-       REAL(wp)     ::  dc                !<
-       REAL(wp)     ::  evap              !<
+       INTEGER(iwp) ::  i                 !< loop index
+       INTEGER(iwp) ::  j                 !< loop index
+       INTEGER(iwp) ::  k                 !< loop index
+
+       REAL(wp)     ::  cond              !< condensation rate
+       REAL(wp)     ::  cond_max          !< maximum condensation rate
+       REAL(wp)     ::  dc                !< weight avageraed diameter
+       REAL(wp)     ::  evap              !< evaporation rate
        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
-       REAL(wp)     ::  g_fac             !<
-       REAL(wp)     ::  nc_0              !<
-       REAL(wp)     ::  temp              !<
-       REAL(wp)     ::  xc                !<
-
+       REAL(wp)     ::  g_fac             !< factor 1 / Fk + Fd
+       REAL(wp)     ::  nc_0              !< integral diameter
+       REAL(wp)     ::  temp              !< actual temperature
+       REAL(wp)     ::  xc                !< mean mass droplets
 
        DO  k = nzb+1, nzt
@@ -3548,14 +4432,15 @@
           CALL supersaturation ( i, j, k )
 !
-!--       Actual temperature:
-          IF ( microphysics_seifert ) THEN
-             temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
-          ELSEIF ( microphysics_morrison_no_rain ) THEN
-             temp = t_l + lv_d_cp * qc(k,j,i) 
-          ENDIF  
-
-          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
-                              l_v / ( thermal_conductivity_l * temp )    &
-                              + r_v * temp / ( diff_coeff_l * e_s )      &
+!--       Actual temperature, t_l is calculated directly before
+!--       in supersaturation
+          IF ( microphysics_ice_extension ) THEN
+             temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
+          ELSE
+             temp = t_l + lv_d_cp * ql(k,j,i)
+          ENDIF
+
+          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *              &
+                              l_v / ( thermal_conductivity_l * temp )          &
+                              + r_v * temp / ( diff_coeff_l * e_s )            &
                             )
 !
@@ -3594,4 +4479,182 @@
 
     END SUBROUTINE condensation_ij
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate the growth of ice particles by water vapor deposition (after
+!> Seifert and Beheng, 2006).
+!------------------------------------------------------------------------------!
+    SUBROUTINE ice_deposition
+
+       INTEGER(iwp) ::  i                 !< loop index
+       INTEGER(iwp) ::  j                 !< loop index
+       INTEGER(iwp) ::  k                 !< loop index
+
+       REAL(wp) ::  ac = 0.09_wp        !< parameter for ice capacitance
+       REAL(wp) ::  bc = 0.33_wp        !< parameter for ice capacitance
+       REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral
+                                        !< distribution, here following gamma
+                                        !< size distribution with Âµ =1/3 and nu=0
+       REAL(wp) ::  deposition_rate     !< depositions rate
+       REAL(wp) ::  deposition_rate_max !< maximum deposition rate
+       REAL(wp) ::  sublimation_rate    !< sublimations rate
+       REAL(wp) ::  gfac_dep            !< factor
+       REAL(wp) ::  temp                !< actual temperature
+       REAL(wp) ::  xi                  !< mean mass of ice crystal
+       REAL(wp) ::  flag                !< flag to indicate first grid level above
+
+       CALL cpu_log( log_point_s(95), 'ice deposition', 'start' )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--             Call calculation of supersaturation over a plane ice surface
+                CALL supersaturation_ice ( i, j, k )
+!
+!--             Actual temperature:
+                temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
+
+                IF ( temp >= 273.15_wp ) CYCLE
+!
+!--             calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g.
+!--             Rogers and Yau, 1989
+                gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *     &
+                                    l_s / ( thermal_conductivity_l * temp )    &
+                                    + r_v * temp / ( diff_coeff_l * e_si )     &
+                                    )
+!
+!--             If there is nothing nucleated, than there is also no
+!--             deposition (above -38Â°C)
+                IF ( ni(k,j,i) <= 0.0_wp ) CYCLE
+!
+!--             calculate mean mass of ice crystal
+                xi = 0.0_wp
+                xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin )
+!
+!--             Condensation needs only to be calculated in supersaturated
+!--             regions (regarding ice)
+                IF ( sat_ice > 0.0_wp )  THEN
+!
+!--                Calculate deposition rate assuming ice crystal shape as
+!--                prescribed in Ovchinnikov et al., 2014 and a gamma size
+!--                distribution according to Seifert and Beheng with to
+!--                Âµ =1/3 and nu=0
+                   deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep *       &
+                                      fac_gamma * ac * xi**bc * ni(k,j,i)
+                   IF ( microphysics_seifert ) THEN
+                      deposition_rate_max  = q(k,j,i) -                        &
+                                             q_si - qr(k,j,i) - qi(k,j,i)
+                   ELSEIF ( microphysics_morrison_no_rain ) THEN
+                      deposition_rate_max  = q(k,j,i) -                        &
+                                             q_si - qc(k,j,i) - qi(k,j,i)
+                   ENDIF
+                   deposition_rate = MIN( deposition_rate,                     &
+                                          deposition_rate_max / dt_micro )
+
+                   qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag
+                ELSEIF ( sat_ice < 0.0_wp ) THEN
+                   sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep *       &
+                                      fac_gamma * ac * xi**bc * ni(k,j,i)
+                   sublimation_rate = MAX( sublimation_rate,                   &
+                                           -qi(k,j,i) / dt_micro )
+                   qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(95), 'ice deposition', 'stop' )
+
+    END SUBROUTINE ice_deposition
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and
+!> Kogan, 2000).
+!------------------------------------------------------------------------------!
+    SUBROUTINE ice_deposition_ij( i, j )
+
+       INTEGER(iwp) ::  i                 !< loop index
+       INTEGER(iwp) ::  j                 !< loop index
+       INTEGER(iwp) ::  k                 !< loop index
+
+       REAL(wp) ::  ac = 0.09_wp        !< parameter for ice capacitance
+       REAL(wp) ::  bc = 0.33_wp        !< parameter for ice capacitance
+       REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral
+                                        !< distribution, here following gamma
+                                        !< size distribution with nu=1, v=0
+       REAL(wp) ::  deposition_rate     !< depositions rate
+       REAL(wp) ::  deposition_rate_max !< maximum deposition rate
+       REAL(wp) ::  sublimation_rate    !< sublimations rate
+       REAL(wp) ::  gfac_dep            !< factor
+       REAL(wp) ::  temp                !< actual temperature
+       REAL(wp) ::  xi                  !< mean mass of ice crystal
+       REAL(wp) ::  flag                !< flag to indicate first grid level above
+
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--       Call calculation of supersaturation over a plane ice surface
+          CALL supersaturation_ice ( i, j, k )
+!
+!--       Actual temperature:
+          temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
+
+          IF ( temp >= 273.15_wp ) CYCLE
+!
+!--       calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g.
+!--       Rogers and Yau, 1989
+          gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *           &
+                              l_s / ( thermal_conductivity_l * temp )          &
+                              + r_v * temp / ( diff_coeff_l * e_si )           &
+                              )
+!
+!--       If there is nothing nucleated, than there is also no
+!--       deposition (above -38Â°C)
+          IF ( ni(k,j,i) <= 0.0_wp ) CYCLE
+!
+!--       calculate mean mass of ice crystal
+          xi = 0.0_wp
+          xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin )
+!
+!--       Condensation needs only to be calculated in supersaturated
+!--       regions (regarding ice)
+          IF ( sat_ice > 0.0_wp )  THEN
+!
+!--          Calculate deposition rate assuming ice crystal shape as
+!--          prescribed in Ovchinnikov et al., 2014 and a gamma size
+!--          distribution according to Seifert and Beheng with to
+!--          Âµ =1/3 and nu=0
+             deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep *             &
+                                fac_gamma * ac * xi**bc * ni(k,j,i)
+             IF ( microphysics_seifert ) THEN
+                deposition_rate_max  = q(k,j,i) -                              &
+                                       q_si - qr(k,j,i) - qi(k,j,i)
+             ELSEIF ( microphysics_morrison_no_rain ) THEN
+                deposition_rate_max  = q(k,j,i) -                              &
+                                       q_si - qc(k,j,i) - qi(k,j,i)
+             ENDIF
+             deposition_rate = MIN( deposition_rate,                           &
+                                    deposition_rate_max / dt_micro )
+
+             qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag
+          ELSEIF ( sat_ice < 0.0_wp ) THEN
+             sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep *             &
+                                fac_gamma * ac * xi**bc * ni(k,j,i)
+             sublimation_rate = MAX( sublimation_rate,                         &
+                                     -qi(k,j,i) / dt_micro )
+             qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE ice_deposition_ij
 
 
@@ -3823,5 +4886,5 @@
 !--          Autoconversion rate (Seifert and Beheng, 2006):
              autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
-                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
+                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *           &
                        ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
                        rho_surface
@@ -4055,5 +5118,5 @@
           ENDIF
 
-          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
+          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.    &
                ( nc_accr > eps_mr ) )  THEN
 !
@@ -4082,5 +5145,5 @@
 !
 !--          Accretion rate (Seifert and Beheng, 2006):
-             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
+             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                    &
                     SQRT( rho_surface * hyrho(k) )
              accr = MIN( accr, qc(k,j,i) / dt_micro )
@@ -4089,5 +5152,5 @@
              qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
              IF ( microphysics_morrison )  THEN
-                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
+                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *            &
                                              hyrho(k) * dt_micro * flag        &
                                          )
@@ -4587,5 +5650,5 @@
              IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
                 sed_nc(k) = sed_qc_const *                                     &
-                            ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
+                            ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *    &
                             ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
              ELSE
@@ -4594,8 +5657,8 @@
 
              sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
-                              nc(k,j,i) / dt_micro + sed_nc(k+1)                &
+                              nc(k,j,i) / dt_micro + sed_nc(k+1)               &
                             ) * flag
 
-             nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *               &
+             nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *             &
                                          ddzu(k+1) / hyrho(k) * dt_micro * flag
           ENDIF
@@ -4608,13 +5671,13 @@
           ENDIF
 
-          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /          &
+          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /         &
                                       dt_micro + sed_qc(k+1)                   &
                          ) * flag
 
-          q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
+          q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /    &
                                 hyrho(k) * dt_micro * flag
-          qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
+          qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /    &
                                 hyrho(k) * dt_micro * flag
-          pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
+          pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /    &
                                 hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
                                                                 * flag
@@ -4632,4 +5695,184 @@
 
     END SUBROUTINE sedimentation_cloud_ij
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Sedimentation of ice crystals
+!------------------------------------------------------------------------------!
+    SUBROUTINE sedimentation_ice
+
+       INTEGER(iwp) ::  i             !< loop index
+       INTEGER(iwp) ::  j             !< loop index
+       INTEGER(iwp) ::  k             !< loop index
+
+       REAL(wp) ::  flag              !< flag to indicate first grid level
+       REAL(wp) ::  av = 6.39_wp      !< parameter for calculating fall speed
+       REAL(wp) ::  bv = 0.1666_wp    !< parameter (1/6)
+       REAL(wp) ::  xi = 0.0_wp       !< mean mass of ice crystal
+       REAL(wp) ::  vi = 0.0_wp       !< mean fall speed of ice crystal
+       REAL(wp) ::  factor_sed_gamma_k0 = 0.76_wp !< factor for zeroth moment and
+                                                  !< Âµ =1/3 and nu=0
+       REAL(wp) ::  factor_sed_gamma_k1 = 1.61_wp !< factor for first moment and
+                                                  !< Âµ =1/3 and nu=0
+
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_ni  !< sedimentation rate zeroth moment
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qi  !< sedimentation rate fist moment
+
+       CALL cpu_log( log_point_s(96), 'sed_ice', 'start' )
+
+       sed_qi(nzt+1) = 0.0_wp
+       sed_ni(nzt+1) = 0.0_wp
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzt, nzb+1, -1
+
+                IF ( ni(k,j,i) <= 0.0_wp ) THEN
+                   xi = 0.0_wp
+                ELSE
+!--                Calculate mean mass of ice crystal
+                   xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin)
+                ENDIF
+!
+!--             calculate fall speed of ice crystal
+                vi = av * xi**bv
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--             Calculate sedimentation rate for each grid box, factors are
+!--             calculated using
+                IF ( qi(k,j,i) > eps_sb  .AND.  ni(k,j,i) >= 0.0_wp )  THEN
+                   sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag
+                   sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag
+                ELSE
+                   sed_qi(k) = 0.0_wp
+                   sed_ni(k) = 0.0_wp
+                ENDIF
+!
+!--             Calculate sedimentation: divergence of sedimentation flux
+                sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
+                                            dt_micro + sed_qi(k+1)             &
+                               ) * flag
+
+                q(k,j,i)  = q(k,j,i)  + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)&
+                                      / hyrho(k) * dt_micro * flag
+                qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)&
+                                      / hyrho(k) * dt_micro * flag
+                ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1)&
+                                      / hyrho(k) * dt_micro * flag
+
+                pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)&
+                                      / hyrho(k) * l_s / c_p * d_exner(k) *    &
+                                        dt_micro * flag
+!
+!--             Compute the precipitation rate of cloud (fog) droplets
+                IF ( call_microphysics_at_all_substeps )  THEN
+                   prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) *            &
+                                weight_substep(intermediate_timestep_count) *  &
+                                flag
+                ELSE
+                   prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(96), 'sed_ice', 'stop' )
+
+    END SUBROUTINE sedimentation_ice
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Sedimentation of ice crystals
+!------------------------------------------------------------------------------!
+    SUBROUTINE sedimentation_ice_ij( i, j )
+
+       INTEGER(iwp) ::  i             !<
+       INTEGER(iwp) ::  j             !<
+       INTEGER(iwp) ::  k             !<
+
+       REAL(wp) ::  flag           !< flag to indicate first grid level
+       REAL(wp) ::  av = 6.39_wp   !< parameter for calculating fall speed
+       REAL(wp) ::  bv = 0.1666_wp !<
+       REAL(wp) ::  xi = 0.0_wp    !< mean mass of ice crystal
+       REAL(wp) ::  vi = 0.0_wp    !< mean fall speed of ice crystal
+       REAL(wp) ::  factor_sed_gamma_k0 = 0.76_wp
+       REAL(wp) ::  factor_sed_gamma_k1 = 1.61_wp
+
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_ni  !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qi  !<
+
+       sed_qi(nzt+1) = 0.0_wp
+       sed_ni(nzt+1) = 0.0_wp
+
+       DO  k = nzt, nzb+1, -1
+          IF ( ni(k,j,i) <= 0.0_wp ) THEN
+             xi = 0.0_wp
+          ELSE
+!--          Calculate mean mass of ice crystal
+             xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin)
+          ENDIF
+!
+!--       calculate fall speed of ice crystal
+          vi = av * xi**bv
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+!
+!--       Calculate sedimentation rate for each grid box, factors are
+!--       calculated using
+          IF ( qi(k,j,i) > eps_sb  .AND.  ni(k,j,i) >= 0.0_wp )  THEN
+             sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag
+             sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag
+          ELSE
+             sed_qi(k) = 0.0_wp
+             sed_ni(k) = 0.0_wp
+          ENDIF
+!
+!--       calculate fall speed of ice crystal
+          vi = av * xi**bv
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
+
+          IF ( qi(k,j,i) > eps_sb  .AND.  ni(k,j,i) >= 0.0_wp )  THEN
+             sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag
+             sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag
+          ELSE
+             sed_qi(k) = 0.0_wp
+             sed_ni(k) = 0.0_wp
+          ENDIF
+
+          sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) /         &
+                                      dt_micro + sed_qi(k+1)                   &
+                         ) * flag
+
+          q(k,j,i)  = q(k,j,i)  + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) /    &
+                                hyrho(k) * dt_micro * flag
+          qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) /    &
+                                hyrho(k) * dt_micro * flag
+          ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1) /    &
+                                hyrho(k) * dt_micro * flag
+          pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) /    &
+                                hyrho(k) * l_s / c_p * d_exner(k) * dt_micro   &
+                                                                  * flag
+!
+!--       Compute the precipitation rate of cloud (fog) droplets
+          IF ( call_microphysics_at_all_substeps )  THEN
+             prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) *                  &
+                              weight_substep(intermediate_timestep_count) * flag
+          ELSE
+             prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE sedimentation_ice_ij
+
 
 
@@ -5053,5 +6296,5 @@
 
           sed_nr(k) = flux / dt_micro * flag
-          nr(k,j,i)  = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
+          nr(k,j,i)  = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /   &
                                     hyrho(k) * dt_micro * flag
 !
@@ -5080,9 +6323,9 @@
           sed_qr(k) = flux / dt_micro * flag
 
-          qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
+          qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /    &
                                 hyrho(k) * dt_micro * flag
-          q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
+          q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /    &
                                 hyrho(k) * dt_micro * flag
-          pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
+          pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /    &
                                 hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
                                                                 * flag
@@ -5202,15 +6445,59 @@
 !--    (see: Cuijpers + Duynkerke, 1993, JAS, 23)
        q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
+
+       IF ( .NOT. microphysics_ice_extension ) THEN
+          IF ( microphysics_seifert ) THEN
+             sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
+          ELSEIF ( microphysics_morrison_no_rain ) THEN
+             sat = ( q(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
+          ENDIF
+       ELSE
+          IF ( microphysics_seifert ) THEN
+             sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp
+          ELSEIF ( microphysics_morrison_no_rain ) THEN
+             sat = ( q(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp
+          ENDIF
+       ENDIF
+
+    END SUBROUTINE supersaturation
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Computation of the diagnostic supersaturation sat, actual liquid water
+!< temperature t_l and saturation water vapor mixing ratio q_s
+!------------------------------------------------------------------------------!
+    SUBROUTINE supersaturation_ice ( i, j, k )
+
+       INTEGER(iwp) ::  i   !< running index
+       INTEGER(iwp) ::  j   !< running index
+       INTEGER(iwp) ::  k   !< running index
+
+       REAL(wp) ::  e_a     !< water vapor pressure
+       REAL(wp) ::  temp    !< actual temperature
+!
+!--    Actual liquid water temperature:
+       t_l = exner(k) * pt(k,j,i)
+!
+!--    Actual temperature:
+       temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
+!
+!--    Calculate water vapor saturation pressure with formular using actual
+!--    temperature
+       e_si = magnus_ice( temp )
+!
+!--    Computation of ice saturation mixing ratio:
+       q_si = rd_d_rv * e_si / ( hyp(k) - e_si )
+!
+!--    Current vapor pressure
+       e_a = (   q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) * hyp(k) /     &
+             ( ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) + rd_d_rv )
 !
 !--    Supersaturation:
 !--    Not in case of microphysics_kessler or microphysics_sat_adjust
 !--    since qr is unallocated
-       IF ( microphysics_seifert ) THEN
-          sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
-       ELSEIF ( microphysics_morrison_no_rain ) THEN
-          sat = ( q(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
-       ENDIF
-
-    END SUBROUTINE supersaturation
+       sat_ice = e_a / e_si - 1.0_wp
+
+    END SUBROUTINE supersaturation_ice
 
 
@@ -5224,8 +6511,4 @@
     SUBROUTINE calc_liquid_water_content
 
-
-
-       IMPLICIT NONE
-
        INTEGER(iwp) ::  i !<
        INTEGER(iwp) ::  j !<
@@ -5234,5 +6517,4 @@
        REAL(wp)     ::  flag !< flag to indicate first grid level above surface
 
-
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -5241,32 +6523,77 @@
 !--             Predetermine flag to mask topography
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
-
 !
 !--             Call calculation of supersaturation located
                 CALL supersaturation( i, j, k )
-
 !
 !--             Compute the liquid water content
-                IF ( microphysics_seifert  .AND. .NOT. microphysics_morrison ) &
-                THEN
-                   IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
-                      qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
-                      ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
+                IF ( .NOT. microphysics_ice_extension ) THEN
+                   IF ( microphysics_seifert  .AND. .NOT.                      &
+                        microphysics_morrison )  THEN
+!--                   Seifert and Beheng scheme: saturation adjustment
+                      IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
+                         qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
+                         ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
+                      ELSE
+                         IF ( q(k,j,i) < qr(k,j,i) )  q(k,j,i) = qr(k,j,i)
+                         qc(k,j,i) = 0.0_wp
+                         ql(k,j,i) = qr(k,j,i) * flag
+                      ENDIF
+!
+!--                Morrison scheme: explicit condensation (see above)
+                   ELSEIF ( microphysics_morrison  .AND.                       &
+                            microphysics_seifert )  THEN
+                      ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
+!--                Morrison without rain: explicit condensation
+                   ELSEIF ( microphysics_morrison  .AND.                       &
+                            .NOT. microphysics_seifert )  THEN
+                      ql(k,j,i) = qc(k,j,i)
+!--                Kessler and saturation adjustment scheme
                    ELSE
-                      IF ( q(k,j,i) < qr(k,j,i) )  q(k,j,i) = qr(k,j,i)
-                      qc(k,j,i) = 0.0_wp
-                      ql(k,j,i) = qr(k,j,i) * flag
+                      IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
+                         qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
+                         ql(k,j,i) = qc(k,j,i) * flag
+                      ELSE
+                         qc(k,j,i) = 0.0_wp
+                         ql(k,j,i) = 0.0_wp
+                      ENDIF
                    ENDIF
-                ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
-                   ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
-                ELSEIF ( microphysics_morrison  .AND.  .NOT. microphysics_seifert )  THEN
-                   ql(k,j,i) = qc(k,j,i)                    
+!--             Calculations of liquid water content in case of mixed-phase
+!--             cloud microphyics
                 ELSE
-                   IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
-                      qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
-                      ql(k,j,i) = qc(k,j,i) * flag
+                   IF ( microphysics_seifert  .AND.                            &
+                        .NOT. microphysics_morrison ) &
+                   THEN
+!
+!--                   Seifert and Beheng scheme: saturation adjustment
+                      IF ( ( q(k,j,i)                                          &
+                             - q_s - qr(k,j,i) - qi(k,j,i) ) > 0.0_wp )  THEN
+                         qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) - qi(k,j,i) )&
+                                       * flag
+                         ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
+                      ELSE
+                         IF ( q(k,j,i) < ( qr(k,j,i) + qi(k,j,i) ) )  THEN
+                            q(k,j,i) = qr(k,j,i) + qi(k,j,i)
+                         ENDIF
+                         qc(k,j,i) = 0.0_wp
+                         ql(k,j,i) = qr(k,j,i) * flag
+                      ENDIF
+!--                Morrison scheme: explicit condensation (see above)
+                   ELSEIF ( microphysics_morrison  .AND.                       &
+                            microphysics_seifert )  THEN
+                      ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
+!--                Morrison without rain: explicit condensation
+                   ELSEIF ( microphysics_morrison  .AND.                       &
+                            .NOT. microphysics_seifert )  THEN
+                      ql(k,j,i) = qc(k,j,i)
+!--                Kessler and saturation adjustment scheme
                    ELSE
-                      qc(k,j,i) = 0.0_wp
-                      ql(k,j,i) = 0.0_wp
+                      IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
+                         qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
+                         ql(k,j,i) = qc(k,j,i) * flag
+                      ELSE
+                         qc(k,j,i) = 0.0_wp
+                         ql(k,j,i) = 0.0_wp
+                      ENDIF
                    ENDIF
                 ENDIF
