Changeset 4502 for palm

Timestamp:

Apr 17, 2020 4:14:16 PM (11 months ago)

Author:

schwenkel

Message:

Implementation of ice microphysics

Location:

palm/trunk/SOURCE

Files:

• advec_s_bc.f90 (modified) (5 diffs)
• advec_ws.f90 (modified) (7 diffs)
• basic_constants_and_equations_mod.f90 (modified) (5 diffs)
• bulk_cloud_model_mod.f90 (modified) (100 diffs)
• compute_vpt.f90 (modified) (4 diffs)
• flow_statistics.f90 (modified) (6 diffs)
• init_masks.f90 (modified) (4 diffs)
• modules.f90 (modified) (13 diffs)
• netcdf_interface_mod.f90 (modified) (7 diffs)
• surface_data_output_mod.f90 (modified) (7 diffs)
• surface_mod.f90 (modified) (211 diffs)
• time_integration.f90 (modified) (9 diffs)

palm/trunk/SOURCE/advec_s_bc.f90

@@ -24 +24 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4488 2020-04-03 11:34:29Z raasch
 ! file re-formatted to follow the PALM coding standard
 !
@@ -957 +960 @@
           ENDDO
 
-       ELSEIF ( sk_char == 'qr' )  THEN
-
-!
-!--       Rain water content boundary condition at the bottom boundary: Dirichlet (fixed surface
-!--       rain water content).
+       ELSEIF ( sk_char == 'qi' )  THEN
+
+!
+!--       Ice crystal content boundary condition at the bottom boundary:
+!--       Dirichlet (fixed surface rain water content).
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -971 +974 @@
 
 !
-!--       Rain water content boundary condition at the top boundary: Dirichlet
+!--       Ice crystal content boundary condition at the top boundary: Dirichlet
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -979 +982 @@
           ENDDO
 
-       ELSEIF ( sk_char == 'nc' )  THEN
-
-!
-!--       Cloud drop concentration boundary condition at the bottom boundary: Dirichlet (fixed
-!--       surface cloud drop concentration).
+       ELSEIF ( sk_char == 'qr' )  THEN
+
+!
+!--       Rain water content boundary condition at the bottom boundary: Dirichlet (fixed surface
+!--       rain water content).
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -993 +996 @@
 
 !
+!--       Rain water content boundary condition at the top boundary: Dirichlet
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+             ENDDO
+          ENDDO
+
+       ELSEIF ( sk_char == 'nc' )  THEN
+
+!
+!--       Cloud drop concentration boundary condition at the bottom boundary: Dirichlet (fixed
+!--       surface cloud drop concentration).
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
+                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+             ENDDO
+          ENDDO
+
+!
 !--       Cloud drop concentration boundary condition at the top boundary: Dirichlet
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+             ENDDO
+          ENDDO
+
+       ELSEIF ( sk_char == 'ni' )  THEN
+
+!
+!--       Ice crystal concentration boundary condition at the bottom boundary:
+!--       Dirichlet (fixed surface cloud drop concentration).
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
+                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+             ENDDO
+          ENDDO
+
+!
+!--       Ice crystal concentration boundary condition at the top boundary: Dirichlet
           DO  i = nxl, nxr
              DO  j = nys, nyn

palm/trunk/SOURCE/advec_ws.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4469 2020-03-23 14:31:00Z suehring
 ! fix mistakenly committed version
 !
@@ -240 +243 @@
                sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l,                &
                sums_wsqcs_ws_l, sums_wsqrs_ws_l,                               &
+               sums_wsqis_ws_l, sums_wsnis_ws_l,                               &
                sums_wsncs_ws_l, sums_wsnrs_ws_l,                               &
                hom, weight_substep
@@ -1938 +1942 @@
              ENDDO
 
+          CASE ( 'qi' )
+
+             DO  k = nzb, nzt
+                sums_wsqis_ws_l(k,tn)  = sums_wsqis_ws_l(k,tn) +               &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
+             ENDDO
 
           CASE ( 'qr' )
@@ -1954 +1968 @@
              DO  k = nzb, nzt
                 sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn) +               &
+                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
+                    ) * weight_substep(intermediate_timestep_count)
+             ENDDO
+
+          CASE ( 'ni' )
+
+             DO  k = nzb, nzt
+                sums_wsnis_ws_l(k,tn)  = sums_wsnis_ws_l(k,tn) +               &
                     ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
                                 * ( w(k,j,i) - hom(k,1,3,0)                 )  &
@@ -3798 +3823 @@
            CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
               sk_num = 9
+           CASE ( 'ni' )
+              sk_num = 10
+           CASE ( 'qi' )
+              sk_num = 11
            CASE DEFAULT
               sk_num = 0
@@ -3824 +3853 @@
        !$ACC PRESENT(sums_wsqrs_ws_l, sums_wsncs_ws_l) &
        !$ACC PRESENT(sums_wsnrs_ws_l, sums_wsss_ws_l) &
+       !$ACC PRESENT(sums_wsnis_ws_l, sums_wsqis_ws_l) &
        !$ACC PRESENT(sums_salsa_ws_l)
        DO  i = nxl, nxr
@@ -4517 +4547 @@
                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )   &
                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 10 )
+                       !$ACC ATOMIC
+                       sums_wsnis_ws_l(k,tn)  = sums_wsnis_ws_l(k,tn)          &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
+                    CASE ( 11 )
+                       !$ACC ATOMIC
+                       sums_wsqis_ws_l(k,tn)  = sums_wsqis_ws_l(k,tn)          &
+                          + ( flux_t(k)                                        &
+                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
+                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
+                            + diss_t(k)                                        &
+                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
+                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
+                            ) * weight_substep(intermediate_timestep_count)
 
                 END SELECT

palm/trunk/SOURCE/basic_constants_and_equations_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4400 2020-02-10 20:32:41Z suehring
 ! Move routine to transform coordinates from netcdf_interface_mod to 
 ! basic_constants_and_equations_mod
@@ -92 +95 @@
     REAL(wp), PARAMETER ::  g_d_cp  = g   / c_p   !< precomputed g / c_p
     REAL(wp), PARAMETER ::  lv_d_cp = l_v / c_p   !< precomputed l_v / c_p
+    REAL(wp), PARAMETER ::  ls_d_cp = l_s / c_p   !< precomputed l_s / c_p
     REAL(wp), PARAMETER ::  lv_d_rd = l_v / r_d   !< precomputed l_v / r_d
     REAL(wp), PARAMETER ::  rd_d_rv = r_d / r_v   !< precomputed r_d / r_v
@@ -106 +110 @@
     PRIVATE magnus_0d, &
             magnus_1d, &
+            magnus_tl_0d, &
+            magnus_tl_1d, &
+            magnus_0d_ice, &
+            magnus_1d_ice, &
             ideal_gas_law_rho_0d, &
             ideal_gas_law_rho_1d, &
@@ -126 +134 @@
        MODULE PROCEDURE magnus_1d
     END INTERFACE magnus
+
+    INTERFACE magnus_tl
+       MODULE PROCEDURE magnus_tl_0d
+       MODULE PROCEDURE magnus_tl_1d
+    END INTERFACE magnus_tl
+
+    INTERFACE magnus_ice
+       MODULE PROCEDURE magnus_0d_ice
+       MODULE PROCEDURE magnus_1d_ice
+    END INTERFACE magnus_ice
 
     INTERFACE ideal_gas_law_rho
@@ -337 +355 @@
 ! Description:
 ! ------------
+!> This function computes the magnus formula (Press et al., 1992) using the
+!> (ice-) liquid water potential temperature.
+!> The magnus formula is needed to calculate the saturation vapor pressure over
+!> a plane liquid water surface
+!------------------------------------------------------------------------------!
+    FUNCTION magnus_tl_0d( t_l )
+
+       IMPLICIT NONE
+
+       REAL(wp), INTENT(IN) ::  t_l  !< liquid water temperature (K)
+
+       REAL(wp) ::  magnus_tl_0d
+!
+!--    Saturation vapor pressure for a specific temperature:
+       magnus_tl_0d =  610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /             &
+                                                    ( t_l - 35.86_wp  ) )
+
+    END FUNCTION magnus_tl_0d
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This function computes the magnus formula (Press et al., 1992) using the
+!> (ice-) liquid water potential temperature.
+!> The magnus formula is needed to calculate the saturation vapor pressure over
+!> a plane liquid water surface
+!------------------------------------------------------------------------------!
+    FUNCTION magnus_tl_1d( t_l )
+
+       IMPLICIT NONE
+
+       REAL(wp), INTENT(IN), DIMENSION(:) ::  t_l  !< liquid water temperature (K)
+
+       REAL(wp), DIMENSION(size(t_l)) ::  magnus_tl_1d
+!
+!--    Saturation vapor pressure for a specific temperature:
+       magnus_tl_1d =  610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /      &
+                                                    ( t_l - 35.86_wp  ) )
+
+    END FUNCTION magnus_tl_1d
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This function computes the magnus formula (Press et al., 1992).
+!> The magnus formula is needed to calculate the saturation vapor pressure over
+!> a plane ice surface
+!------------------------------------------------------------------------------!
+    FUNCTION magnus_0d_ice( t )
+
+       IMPLICIT NONE
+
+       REAL(wp), INTENT(IN) ::  t  !< temperature (K)
+
+       REAL(wp) ::  magnus_0d_ice
+!
+!--    Saturation vapor pressure for a specific temperature:
+       magnus_0d_ice =  611.2_wp * EXP( 22.46_wp * ( t - degc_to_k ) /         &
+                                                   ( t - 0.53_wp  ) )
+
+    END FUNCTION magnus_0d_ice
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This function computes the magnus formula (Press et al., 1992).
+!> The magnus formula is needed to calculate the saturation vapor pressure over
+!> a plane ice surface
+!------------------------------------------------------------------------------!
+    FUNCTION magnus_1d_ice( t )
+
+       IMPLICIT NONE
+
+       REAL(wp), INTENT(IN), DIMENSION(:) ::  t  !< temperature (K)
+
+       REAL(wp), DIMENSION(size(t)) ::  magnus_1d_ice
+!
+!--    Saturation vapor pressure for a specific temperature:
+       magnus_1d_ice =  611.2_wp * EXP( 22.46_wp * ( t - degc_to_k ) /          &
+                                                   ( t - 0.53_wp  ) )
+
+    END FUNCTION magnus_1d_ice
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Compute the ideal gas law for scalar arguments.
 !------------------------------------------------------------------------------!

palm/trunk/SOURCE/bulk_cloud_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
 ! 
@@ -116 +119 @@
                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 +153 @@
                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 +179 @@
         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 +208 @@
     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 +214 @@
 
     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 +257 @@
     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 +266 @@
     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 +272 @@
     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 +316 @@
            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 +445 @@
           nc_const, &
           sigma_bulk, &
-          ventilation_effect
+          ventilation_effect, &
+          ice_crystal_sedimentation, &
+          microphysics_ice_extension, &
+          start_ice_microphysics, &
+          in_init
 
        line = ' '
@@ -517 +546 @@
           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 +575 @@
        surf_microphysics_morrison = microphysics_morrison
        surf_microphysics_seifert = microphysics_seifert
-
+       surf_microphysics_ice_extension = microphysics_ice_extension
 !
 !--    Check aerosol
@@ -592 +621 @@
              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 +649 @@
 
           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 +745 @@
              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 +791 @@
              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 +868 @@
 !
 !--       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 +879 @@
 !
 !--       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 +905 @@
              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 +911 @@
              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 +925 @@
 !--    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 +936 @@
                           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 +946 @@
                           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 +973 @@
           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 +1019 @@
              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 +1048 @@
                 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 +1196 @@
                 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 +1216 @@
 !> Control of microphysics for grid points i,j
 !------------------------------------------------------------------------------!
-
     SUBROUTINE bcm_actions_ij( i, j, location )
 
@@ -1121 +1240 @@
                 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 +1267 @@
        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 +1275 @@
 !
 !--          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 +1306 @@
              CALL autoconversion_kessler
              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
-             
+
 !
 !--       Here the seifert beheng scheme is used. Cloud concentration is assumed to
@@ -1184 +1312 @@
           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 +1322 @@
              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 +1348 @@
              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 +1359 @@
              CALL evaporation_rain
              CALL sedimentation_rain
-             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud            
+             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
 
           ENDIF
@@ -1229 +1377 @@
 !> Control of microphysics for grid points i,j
 !------------------------------------------------------------------------------!
-
     SUBROUTINE bcm_non_advective_processes_ij( i, j )
 
@@ -1236 +1383 @@
        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 +1391 @@
 !
 !--          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 +1424 @@
 !
 !--       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 +1437 @@
              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 +1473 @@
              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 +1501 @@
           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 +1507 @@
              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 +1599 @@
                                              )                                 &
                                     * 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 +1691 @@
                                              )                                 &
                                    * 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 +1725 @@
 
        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 +1983 @@
                                              )                                 &
                                     * 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 +2075 @@
                                              )                                 &
                                    * 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 +2118 @@
 !> Control of microphysics for grid points i,j
 !------------------------------------------------------------------------------!
-
     SUBROUTINE bcm_prognostic_equations_ij( i, j, i_omp_start, tn )
 
@@ -1805 +2171 @@
                                        )                                 &
                                  * 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 +2230 @@
                                        )                                 &
                                  * 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 +2251 @@
 
        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 +2421 @@
                                        )                                 &
                                  * 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 +2480 @@
                                        )                                 &
                                  * 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 +2530 @@
                    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 +2545 @@
                    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 +2594 @@
        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 +2654 @@
              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 +2667 @@
              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 +2680 @@
              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 +2692 @@
              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 +2974 @@
              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 +3052 @@
           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 +3099 @@
              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 +3160 @@
              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 +3219 @@
           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 +3262 @@
           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 +3337 @@
           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 +3381 @@
        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 +3451 @@
                 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 +3504 @@
                 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 +3615 @@
           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 +3644 @@
           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 +3709 @@
 
           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 +3772 @@
              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 +3808 @@
 !--             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 +3864 @@
           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 +3895 @@
 
     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 +4176 @@
     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 +4319 @@
        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 +4345 @@
                 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 +4410 @@
        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 +4432 @@
           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 +4479 @@
 
     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 +4886 @@
 !--          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 +5118 @@
           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 +5145 @@
 !
 !--          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 +5152 @@
              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 +5650 @@
              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 +5657 @@
 
              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 +5671 @@
           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 +5695 @@
 
     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 +6296 @@
 
           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 +6323 @@
           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 +6445 @@
 !--    (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 +6511 @@
     SUBROUTINE calc_liquid_water_content
 
-
-
-       IMPLICIT NONE
-
        INTEGER(iwp) ::  i !<
        INTEGER(iwp) ::  j !<
@@ -5234 +6517 @@
        REAL(wp)     ::  flag !< flag to indicate first grid level above surface
 
-
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
@@ -5241 +6523 @@
 !--             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

palm/trunk/SOURCE/compute_vpt.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -42 +45 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  pt, q, ql, vpt, d_exner
+        ONLY:  pt, q, ql, vpt, d_exner, qi
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  lv_d_cp
+        ONLY:  lv_d_cp, ls_d_cp
 
     USE control_parameters,                                                    &
@@ -56 +59 @@
 
     USE bulk_cloud_model_mod,                                                  &
-        ONLY:  bulk_cloud_model
+        ONLY:  bulk_cloud_model, microphysics_ice_extension
 
     IMPLICIT NONE
@@ -64 +67 @@
     IF ( .NOT. bulk_cloud_model  .AND.  .NOT. cloud_droplets )  THEN
        vpt = pt * ( 1.0_wp + 0.61_wp * q )
-    ELSE IF (bulk_cloud_model)  THEN
+    ELSEIF ( bulk_cloud_model  .AND.  .NOT. microphysics_ice_extension )  THEN
        DO  k = nzb, nzt+1
-          vpt(k,:,:) = ( pt(k,:,:) + d_exner(k) * lv_d_cp * ql(k,:,:) ) *      &
-                       ( 1.0_wp + 0.61_wp * q(k,:,:) - 1.61_wp * ql(k,:,:) ) 
+              vpt(k,:,:) = ( pt(k,:,:) + d_exner(k) * lv_d_cp * ql(k,:,:) ) *  &
+                       ( 1.0_wp + 0.61_wp * q(k,:,:) - 1.61_wp *  ql(k,:,:)  )
+       ENDDO
+    ELSEIF ( bulk_cloud_model  .AND.  microphysics_ice_extension ) THEN
+       DO  k = nzb, nzt+1
+          vpt(k,:,:) = ( pt(k,:,:) + d_exner(k) * lv_d_cp * ql(k,:,:)   +      &
+                                     d_exner(k) * ls_d_cp * qi(k,:,:) ) *      &
+                       ( 1.0_wp + 0.61_wp * q(k,:,:) - 1.61_wp          *      &
+                                           ( ql(k,:,:) + qi(k,:,:) ) )
        ENDDO
     ELSE

palm/trunk/SOURCE/flow_statistics.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4472 2020-03-24 12:21:00Z Giersch
 ! Calculations of the Kolmogorov lengt scale eta implemented
 !
@@ -84 +87 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
-               momentumflux_output_conversion, nc, nr, p, prho, prr, pt, q,    &
-               qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
+               momentumflux_output_conversion, nc, ni, nr, p, prho, prr, pt, q,&
+               qc, qi, ql, qr, rho_air, rho_air_zw, rho_ocean, s,              &
                sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion,  &
                zw, d_exner
@@ -93 +96 @@
 
     USE bulk_cloud_model_mod,                                                  &
-        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
+        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert,   &
+              microphysics_ice_extension
 
     USE chem_modules,                                                          &
@@ -1414 +1418 @@
                                                                 flag
                          ENDIF
+                         IF ( microphysics_ice_extension )  THEN
+                            sums_l(k,124,tn) = sums_l(k,124,tn) + ni(k,j,i) *  &
+                                                                rmask(j,i,sr) *&
+                                                                flag
+                            sums_l(k,125,tn) = sums_l(k,125,tn) + qi(k,j,i) *  &
+                                                                rmask(j,i,sr) *&
+                                                                flag
+                         ENDIF
+
                          IF ( microphysics_seifert )  THEN
                             sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
@@ -1915 +1928 @@
              sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
              sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
-             sums(k,123)           = sums(k,123) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
+             sums(k,123:125)       = sums(k,123:125) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
           ENDIF
        ENDDO
@@ -2029 +2042 @@
        hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
        hom(:,1,123,sr) = sums(:,123)   ! nc
+       hom(:,1,124,sr) = sums(:,124)   ! ni
+       hom(:,1,125,sr) = sums(:,125)   ! qi
        hom(:,1,73,sr) = sums(:,73)     ! nr
        hom(:,1,74,sr) = sums(:,74)     ! qr

palm/trunk/SOURCE/init_masks.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4444 2020-03-05 15:59:50Z raasch
 ! bugfix: cpp-directives for serial mode added
 ! 
@@ -58 +61 @@
 
     USE bulk_cloud_model_mod,                                                  &
-        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
+        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert,   &
+              microphysics_ice_extension
 
     USE control_parameters,                                                    &
@@ -268 +272 @@
                 unit = '1/m3'
 
+             CASE ( 'ni' )
+                IF ( .NOT. bulk_cloud_model )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
+                        '" requires bulk_cloud_model = .TRUE.'
+                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
+                 ELSEIF ( .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. bulk_cloud_model )  THEN
@@ -331 +347 @@
                         '" requires bulk_cloud_model = .TRUE.'
                    CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'kg/kg'
+
+             CASE ( 'qi' )
+                IF ( .NOT. bulk_cloud_model )  THEN
+                   message_string = 'output of "' // TRIM( var ) // '" ' //    &
+                            'requires bulk_cloud_model = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
+                ELSEIF ( .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'

palm/trunk/SOURCE/modules.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! +restart_data_format, restart_data_format_input|output, include_total_domain_boundaries
 ! 
@@ -229 +232 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_e              !< artificial numerical dissipation flux at south face of grid box - subgrid-scale TKE
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nc             !< artificial numerical dissipation flux at south face of grid box - clouddrop-number concentration
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_ni             !< artificial numerical dissipation flux at south face of grid box - ice crystal-number concentration
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nr             !< artificial numerical dissipation flux at south face of grid box - raindrop-number concentration
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_pt             !< artificial numerical dissipation flux at south face of grid box - potential temperature
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_q              !< artificial numerical dissipation flux at south face of grid box - mixing ratio
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qc             !< artificial numerical dissipation flux at south face of grid box - cloudwater
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qi             !< artificial numerical dissipation flux at south face of grid box - ice crystal mixing ratio
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qr             !< artificial numerical dissipation flux at south face of grid box - rainwater
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_s              !< artificial numerical dissipation flux at south face of grid box - passive scalar
@@ -244 +249 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_e              !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nc             !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_ni             !< 6th-order advective flux at south face of grid box - icecrystal-number concentration
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nr             !< 6th-order advective flux at south face of grid box - raindrop-number concentration
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_pt             !< 6th-order advective flux at south face of grid box - potential temperature
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_q              !< 6th-order advective flux at south face of grid box - mixing ratio
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qc             !< 6th-order advective flux at south face of grid box - cloudwater
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qi             !< 6th-order advective flux at south face of grid box - ice crystal
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qr             !< 6th-order advective flux at south face of grid box - rainwater
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_s              !< 6th-order advective flux at south face of grid box - passive scalar
@@ -271 +278 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_e    !< artificial numerical dissipation flux at left face of grid box - subgrid-scale TKE
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nc   !< artificial numerical dissipation flux at left face of grid box - clouddrop-number concentration
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_ni   !< artificial numerical dissipation flux at left face of grid box - ice crystal-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nr   !< artificial numerical dissipation flux at left face of grid box - raindrop-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_pt   !< artificial numerical dissipation flux at left face of grid box - potential temperature
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_q    !< artificial numerical dissipation flux at left face of grid box - mixing ratio
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qc   !< artificial numerical dissipation flux at left face of grid box - cloudwater
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qi   !< artificial numerical dissipation flux at left face of grid box - ice crystal
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qr   !< artificial numerical dissipation flux at left face of grid box - rainwater
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_s    !< artificial numerical dissipation flux at left face of grid box - passive scalar
@@ -284 +293 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_e    !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nc   !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_ni   !< 6th-order advective flux at south face of grid box - ice crystal-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nr   !< 6th-order advective flux at south face of grid box - raindrop-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_pt   !< 6th-order advective flux at south face of grid box - potential temperature
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_q    !< 6th-order advective flux at south face of grid box - mixing ratio
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qc   !< 6th-order advective flux at south face of grid box - cloudwater
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qi   !< 6th-order advective flux at south face of grid box - ice crystal
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qr   !< 6th-order advective flux at south face of grid box - rainwater
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_s    !< 6th-order advective flux at south face of grid box - passive scalar
@@ -324 +335 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_2    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_3    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_1    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_2    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_3    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_1    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_2    !< pointer for swapping of timelevels for respective quantity
@@ -336 +350 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_2    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_3    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_1    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_2    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_3    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_v    !< pointer: volume of liquid water
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_vp   !< pointer: liquid water weighting factor
@@ -367 +384 @@
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nc         !< pointer: cloud drop number density
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nc_p       !< pointer: prognostic value of cloud drop number density
+    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ni         !< pointer: ice crystal number density
+    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ni_p       !< pointer: prognostic value of ice crystal number density
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nr         !< pointer: rain drop number density
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nr_p       !< pointer: prognostic value of rain drop number density
@@ -375 +394 @@
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  q_p        !< pointer: prognostic value of mixing ratio
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qc         !< pointer: cloud water content
-    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qc_p       !< pointer: cloud water content
+    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qc_p       !< pointer: prognostic value cloud water content
+    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qi         !< pointer: ice crystal content
+    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qi_p       !< pointer: prognostic value ice crystal content
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ql         !< pointer: liquid water content
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ql_c       !< pointer: change in liquid water content due to
@@ -389 +410 @@
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  te_m       !< pointer: weighted tendency of e for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tnc_m      !< pointer: weighted tendency of nc for previous sub-timestep (Runge-Kutta)
+    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tni_m      !< pointer: weighted tendency of ni for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tnr_m      !< pointer: weighted tendency of nr for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tpt_m      !< pointer: weighted tendency of pt for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tq_m       !< pointer: weighted tendency of q for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tqc_m      !< pointer: weighted tendency of qc for previous sub-timestep (Runge-Kutta)
+    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tqi_m      !< pointer: weighted tendency of qi for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tqr_m      !< pointer: weighted tendency of qr for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ts_m       !< pointer: weighted tendency of s for previous sub-timestep (Runge-Kutta)
@@ -462 +485 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  lpt_av        !< avg. liquid water potential temperature
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_av         !< avg. cloud drop number density
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_av         !< avg. ice crystal number density
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_av         !< avg. rain drop number density
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  p_av          !< avg. perturbation pressure
@@ -471 +495 @@
                                                                       !< (or total water content with active cloud physics)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_av         !< avg. cloud water content
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_av         !< avg. ice crystal content
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_av         !< avg. liquid water content
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_c_av       !< avg. change in liquid water content due to
@@ -1423 +1448 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_ws2_ws_l    !< subdomain sum of vertical momentum flux w'w' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsncs_ws_l  !< subdomain sum of vertical clouddrop-number concentration flux w'nc' (5th-order advection scheme only)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsnis_ws_l  !< subdomain sum of vertical ice crystal concentration flux w'ni' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsnrs_ws_l  !< subdomain sum of vertical raindrop-number concentration flux w'nr' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wspts_ws_l  !< subdomain sum of vertical sensible heat flux w'pt' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqs_ws_l   !< subdomain sum of vertical latent heat flux w'q' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqcs_ws_l  !< subdomain sum of vertical cloudwater flux w'qc' (5th-order advection scheme only)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqis_ws_l  !< subdomain sum of vertical ice crystal flux w'qi' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqrs_ws_l  !< subdomain sum of vertical rainwater flux w'qr' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wssas_ws_l  !< subdomain sum of vertical salinity flux w'sa' (5th-order advection scheme only)

palm/trunk/SOURCE/netcdf_interface_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4455 2020-03-11 12:20:29Z Giersch
 ! Axis attribute added to netcdf output
 ! 
@@ -880 +883 @@
                 CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       's', 'theta', 'thetal', 'thetav' )
+                       's', 'theta', 'thetal', 'thetav', 'qi', 'ni' )
 
                    grid_x = 'x'
@@ -1644 +1647 @@
                 CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       's', 'theta', 'thetal', 'thetav' )
+                       's', 'theta', 'thetal', 'thetav', 'qi', 'ni' )
 
                    grid_x = 'x'
@@ -2659 +2662 @@
 !
 !--                   Most variables are defined on the zu grid
-                      CASE ( 'e_xy', 'nc_xy', 'nr_xy', 'p_xy',                 &
+                      CASE ( 'e_xy', 'nc_xy', 'ni_xy', 'nr_xy', 'p_xy',        &
                              'pc_xy', 'pr_xy', 'prr_xy', 'q_xy',               &
-                             'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
+                             'qc_xy', 'qi_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',  &
                              'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
                              's_xy',                                           &
@@ -3583 +3586 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz',           &
-                          'pr_xz', 'prr_xz', 'q_xz', 'qc_xz',                  &
+                   CASE ( 'e_xz', 'nc_xz', 'ni_xz', 'nr_xz', 'p_xz', 'pc_xz',  &
+                          'pr_xz', 'prr_xz', 'q_xz', 'qc_xz', 'qi_xz',         &
                           'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
                           'qv_xz', 's_xz',                                     &
@@ -4460 +4463 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz',           &
-                          'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'ql_yz',          &
+                   CASE ( 'e_yz', 'nc_yz', 'ni_yz', 'nr_yz', 'p_yz', 'pc_yz',  &
+                          'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'qi_yz', 'ql_yz', &
                           'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
                           's_yz',                                              &
@@ -5853 +5856 @@
 !
 !--             Most variables are defined on the zu levels
-                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
+                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'ni', 'qi',    &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
                        'rho_sea_water', 's', 'sa', &

palm/trunk/SOURCE/surface_data_output_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4500 2020-04-17 10:12:45Z suehring
 ! - Correct output of ground/wall heat flux at USM surfaces
 ! - Add conversion factor to heat and momentum-flux output
@@ -1695 +1698 @@
                ENDIF
 
+            CASE ( 'qis' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_data_output_collect( surf_def_h(0)%qis,         &
+                                               surf_def_h(1)%qis,              &
+                                               surf_lsm_h%qis,                 &
+                                               surf_usm_h%qis,                 &
+                                               surf_def_v(0)%qis,              &
+                                               surf_lsm_v(0)%qis,              &
+                                               surf_usm_v(0)%qis,              &
+                                               surf_def_v(1)%qis,              &
+                                               surf_lsm_v(1)%qis,              &
+                                               surf_usm_v(1)%qis,              &
+                                               surf_def_v(2)%qis,              &
+                                               surf_lsm_v(2)%qis,              &
+                                               surf_usm_v(2)%qis,              &
+                                               surf_def_v(3)%qis,              &
+                                               surf_lsm_v(3)%qis,              &
+                                               surf_usm_v(3)%qis )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+
+               ENDIF
+
+            CASE ( 'nis' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_data_output_collect( surf_def_h(0)%nis,         &
+                                               surf_def_h(1)%nis,              &
+                                               surf_lsm_h%nis,                 &
+                                               surf_usm_h%nis,                 &
+                                               surf_def_v(0)%nis,              &
+                                               surf_lsm_v(0)%nis,              &
+                                               surf_usm_v(0)%nis,              &
+                                               surf_def_v(1)%nis,              &
+                                               surf_lsm_v(1)%nis,              &
+                                               surf_usm_v(1)%nis,              &
+                                               surf_def_v(2)%nis,              &
+                                               surf_lsm_v(2)%nis,              &
+                                               surf_usm_v(2)%nis,              &
+                                               surf_def_v(3)%nis,              &
+                                               surf_lsm_v(3)%nis,              &
+                                               surf_usm_v(3)%nis )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+
+               ENDIF
+
             CASE ( 'qrs' )
 !
@@ -2153 +2214 @@
                                                surf_lsm_v(3)%ncsws,            &
                                                surf_usm_v(3)%ncsws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+
+               ENDIF
+
+
+            CASE ( 'qisws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_data_output_collect( surf_def_h(0)%qisws,       &
+                                               surf_def_h(1)%qisws,            &
+                                               surf_lsm_h%qisws,               &
+                                               surf_usm_h%qisws,               &
+                                               surf_def_v(0)%qisws,            &
+                                               surf_lsm_v(0)%qisws,            &
+                                               surf_usm_v(0)%qisws,            &
+                                               surf_def_v(1)%qisws,            &
+                                               surf_lsm_v(1)%qisws,            &
+                                               surf_usm_v(1)%qisws,            &
+                                               surf_def_v(2)%qisws,            &
+                                               surf_lsm_v(2)%qisws,            &
+                                               surf_usm_v(2)%qisws,            &
+                                               surf_def_v(3)%qisws,            &
+                                               surf_lsm_v(3)%qisws,            &
+                                               surf_usm_v(3)%qisws )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+
+               ENDIF
+
+            CASE ( 'nisws' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_data_output_collect( surf_def_h(0)%nisws,       &
+                                               surf_def_h(1)%nisws,            &
+                                               surf_lsm_h%nisws,               &
+                                               surf_usm_h%nisws,               &
+                                               surf_def_v(0)%nisws,            &
+                                               surf_lsm_v(0)%nisws,            &
+                                               surf_usm_v(0)%nisws,            &
+                                               surf_def_v(1)%nisws,            &
+                                               surf_lsm_v(1)%nisws,            &
+                                               surf_usm_v(1)%nisws,            &
+                                               surf_def_v(2)%nisws,            &
+                                               surf_lsm_v(2)%nisws,            &
+                                               surf_usm_v(2)%nisws,            &
+                                               surf_def_v(3)%nisws,            &
+                                               surf_lsm_v(3)%nisws,            &
+                                               surf_usm_v(3)%nisws )
                ELSE
 !
@@ -3120 +3240 @@
                                            surf_usm_v(3)%ncs, n_out )
 
+            CASE ( 'qis' )
+               CALL surface_data_output_sum_up( surf_def_h(0)%qis,             &
+                                           surf_def_h(1)%qis,                  &
+                                           surf_lsm_h%qis,                     &
+                                           surf_usm_h%qis,                     &
+                                           surf_def_v(0)%qis,                  &
+                                           surf_lsm_v(0)%qis,                  &
+                                           surf_usm_v(0)%qis,                  &
+                                           surf_def_v(1)%qis,                  &
+                                           surf_lsm_v(1)%qis,                  &
+                                           surf_usm_v(1)%qis,                  &
+                                           surf_def_v(2)%qis,                  &
+                                           surf_lsm_v(2)%qis,                  &
+                                           surf_usm_v(2)%qis,                  &
+                                           surf_def_v(3)%qis,                  &
+                                           surf_lsm_v(3)%qis,                  &
+                                           surf_usm_v(3)%qrs, n_out )
+
+            CASE ( 'nis' )
+               CALL surface_data_output_sum_up( surf_def_h(0)%nis,             &
+                                           surf_def_h(1)%nis,                  &
+                                           surf_lsm_h%nis,                     &
+                                           surf_usm_h%nis,                     &
+                                           surf_def_v(0)%nis,                  &
+                                           surf_lsm_v(0)%nis,                  &
+                                           surf_usm_v(0)%nis,                  &
+                                           surf_def_v(1)%nis,                  &
+                                           surf_lsm_v(1)%nis,                  &
+                                           surf_usm_v(1)%nis,                  &
+                                           surf_def_v(2)%nis,                  &
+                                           surf_lsm_v(2)%nis,                  &
+                                           surf_usm_v(2)%nis,                  &
+                                           surf_def_v(3)%nis,                  &
+                                           surf_lsm_v(3)%nis,                  &
+                                           surf_usm_v(3)%nis, n_out )
+
             CASE ( 'qrs' )
                CALL surface_data_output_sum_up( surf_def_h(0)%qrs,             &
@@ -3411 +3567 @@
                                            surf_lsm_v(3)%ncsws,                &
                                            surf_usm_v(3)%ncsws, n_out )
+
+            CASE ( 'qisws' )
+               CALL surface_data_output_sum_up( surf_def_h(0)%qisws,           &
+                                           surf_def_h(1)%qisws,                &
+                                           surf_lsm_h%qisws,                   &
+                                           surf_usm_h%qisws,                   &
+                                           surf_def_v(0)%qisws,                &
+                                           surf_lsm_v(0)%qisws,                &
+                                           surf_usm_v(0)%qisws,                &
+                                           surf_def_v(1)%qisws,                &
+                                           surf_lsm_v(1)%qisws,                &
+                                           surf_usm_v(1)%qisws,                &
+                                           surf_def_v(2)%qisws,                &
+                                           surf_lsm_v(2)%qisws,                &
+                                           surf_usm_v(2)%qisws,                &
+                                           surf_def_v(3)%qisws,                &
+                                           surf_lsm_v(3)%qisws,                &
+                                           surf_usm_v(3)%qisws, n_out )
+
+            CASE ( 'nisws' )
+               CALL surface_data_output_sum_up( surf_def_h(0)%nisws,           &
+                                           surf_def_h(1)%nisws,                &
+                                           surf_lsm_h%nisws,                   &
+                                           surf_usm_h%nisws,                   &
+                                           surf_def_v(0)%nisws,                &
+                                           surf_lsm_v(0)%nisws,                &
+                                           surf_usm_v(0)%nisws,                &
+                                           surf_def_v(1)%nisws,                &
+                                           surf_lsm_v(1)%nisws,                &
+                                           surf_usm_v(1)%nisws,                &
+                                           surf_def_v(2)%nisws,                &
+                                           surf_lsm_v(2)%nisws,                &
+                                           surf_usm_v(2)%nisws,                &
+                                           surf_def_v(3)%nisws,                &
+                                           surf_lsm_v(3)%nisws,                &
+                                           surf_usm_v(3)%nisws, n_out )
 
             CASE ( 'qrsws' )
@@ -4533 +4725 @@
                unit = 'm/s'
 
-            CASE ( 'ss', 'qcs', 'ncs', 'qrs', 'nrs' )
+            CASE ( 'ss', 'qcs', 'ncs', 'qis', 'nis', 'qrs', 'nrs' )
                unit = '1'
 
@@ -4545 +4737 @@
                unit = 'm2/s2'
 
-            CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' )
+            CASE ( 'qcsws', 'ncsws', 'qisws', 'nisws', 'qrsws', 'nrsws', 'sasws' )
 
             CASE ( 'shf' )

palm/trunk/SOURCE/surface_mod.f90

@@ -26 +26 @@
 ! -----------------
 ! $Id$
+! Implementation of ice microphysics
+! 
+! 4495 2020-04-13 20:11:20Z raasch
 ! restart data handling with MPI-IO added
 ! 
@@ -148 +151 @@
 
 !
-!-- Data type used to identify grid-points where horizontal boundary conditions 
-!-- are applied 
+!-- Data type used to identify grid-points where horizontal boundary conditions
+!-- are applied
     TYPE bc_type
        INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element
@@ -156 +159 @@
        INTEGER(iwp) :: ns   !< number of surface elements on the PE
 
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i !< x-index linking to the PALM 3D-grid  
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j !< y-index linking to the PALM 3D-grid    
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k !< z-index linking to the PALM 3D-grid   
-
-       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i) 
-       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< end index within surface data type for given (j,i)  
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i !< x-index linking to the PALM 3D-grid
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j !< y-index linking to the PALM 3D-grid
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k !< z-index linking to the PALM 3D-grid
+
+       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
+       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< end index within surface data type for given (j,i)
 
     END TYPE bc_type
 !
-!-- Data type used to identify and treat surface-bounded grid points  
+!-- Data type used to identify and treat surface-bounded grid points
     TYPE surf_type
 
        LOGICAL ::  albedo_from_ascii = .FALSE. !< flag indicating that albedo for urban surfaces is input via ASCII format (just for a workaround)
-    
+
        INTEGER(iwp) :: ioff                                !< offset value in x-direction, used to determine index of surface element
        INTEGER(iwp) :: joff                                !< offset value in y-direction, used to determine index of surface element
@@ -175 +178 @@
        INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
 
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid  
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid    
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid       
-
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  facing  !< Bit indicating surface orientation 
-     
-       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i) 
-       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< End index within surface data type for given (j,i)  
-
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_mo      !< surface-layer height 
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid
+
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  facing  !< Bit indicating surface orientation
+
+       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i)
+       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< End index within surface data type for given (j,i)
+
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_mo      !< surface-layer height
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  uvw_abs   !< absolute surface-parallel velocity
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  us        !< friction velocity
@@ -192 +195 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcs       !< scaling parameter qc
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncs       !< scaling parameter nc
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qis       !< scaling parameter qi
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nis       !< scaling parameter ni
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrs       !< scaling parameter qr
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrs       !< scaling parameter nr
@@ -205 +210 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< mixing ratio at first grid level
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  vpt1      !< virtual potential temperature at first grid level
-       
+
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  css     !< scaling parameter chemical species
 !
 !--    Define arrays for surface fluxes
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  usws      !< vertical momentum flux for u-component at horizontal surfaces
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vsws      !< vertical momentum flux for v-component at horizontal surfaces 
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vsws      !< vertical momentum flux for v-component at horizontal surfaces
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf       !< surface flux sensible heat
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws      !< surface flux latent heat 
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws      !< surface flux latent heat
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  ssws      !< surface flux passive scalar
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcsws     !< surface flux qc
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncsws     !< surface flux nc
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qisws     !< surface flux qi
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nisws     !< surface flux ni
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrsws     !< surface flux qr
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrsws     !< surface flux nr
@@ -224 +231 @@
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  amsws   !< surface flux aerosol mass:   dim 1: flux, dim 2: bin
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gtsws   !< surface flux gesous tracers: dim 1: flux, dim 2: gas
-       
+
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  cssws   !< surface flux chemical species
 !
@@ -240 +247 @@
        CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  vegetation_type_name  !< water type at name surface element
        CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  water_type_name       !< water type at name surface element
-       
+
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nzt_pavement     !< top index for pavement in soil
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  building_type    !< building type at surface element
@@ -246 +253 @@
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  vegetation_type  !< vegetation type at surface element
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  water_type       !< water type at surface element
-       
+
        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  albedo_type   !< albedo type, for each fraction (wall,green,window or vegetation,pavement water)
 
@@ -254 +261 @@
        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  water_surface       !< flag parameter for water surfaces
        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  vegetation_surface  !< flag parameter for natural land surfaces
-       
+
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  albedo            !< broadband albedo for each surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  emissivity        !< emissivity of the surface, for each fraction  (LSM: vegetation, water, pavement; USM: wall, green, window)
@@ -271 +278 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pt_surface        !< skin-surface temperature
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  vpt_surface       !< skin-surface virtual temperature
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net           !< net radiation 
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net           !< net radiation
        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l         !< net radiation, used in USM
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil/ wall (W/m/K) 
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_green    !< heat conductivity of green soil (W/m/K) 
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_window   !< heat conductivity of windows (W/m/K) 
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)   
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil/ wall (W/m/K)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_green    !< heat conductivity of green soil (W/m/K)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_window   !< heat conductivity of windows (W/m/K)
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_in           !< incoming longwave radiation
@@ -302 +309 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg            !< surface flux of latent heat (vegetation portion)
 
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a                 !< aerodynamic resistance 
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a                 !< aerodynamic resistance
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a_green           !< aerodynamic resistance at green fraction
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a_window          !< aerodynamic resistance at window fraction
@@ -310 +317 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_s                 !< total surface resistance (combination of r_soil and r_canopy)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy_min        !< minimum canopy (stomatal) resistance
-       
+
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm             !< near surface air potential temperature at distance 10 cm from the surface (K)
-       
+
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  alpha_vg          !< coef. of Van Genuchten
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_w          !< hydraulic diffusivity of soil (?)
@@ -322 +329 @@
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_sat             !< saturation soil moisture (m3/m3)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_wilt            !< soil moisture at permanent wilting point (m3/m3)
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_vg              !< coef. Van Genuchten  
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_vg              !< coef. Van Genuchten
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total_def   !< default volumetric heat capacity of the (soil) layer (J/m3/K)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total       !< volumetric heat capacity of the actual soil matrix (J/m3/K)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  root_fr           !< root fraction within the soil layers
-       
+
 !--    Indoor model variables
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  waste_heat          !< waste heat 
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  waste_heat          !< waste heat
 !
 !--    Urban surface variables
@@ -334 +341 @@
 
        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  isroof_surf          !< flag indicating roof surfaces
-       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  ground_level         !< flag indicating ground floor level surfaces 
+       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  ground_level         !< flag indicating ground floor level surfaces
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_summer  !< indoor target temperature summer
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_winter  !< indoor target temperature summer        
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_winter  !< indoor target temperature summer
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface           !< heat capacity of the wall surface skin (J/m2/K)
@@ -373 +380 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw            !< longwave radiation falling to local surface including radiation from reflections
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw           !< total longwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
-       
+
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  n_vg_green          !< vangenuchten parameters
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  alpha_vg_green      !< vangenuchten parameters
@@ -384 +391 @@
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_wall_stag      !< wall grid spacing (edge-edge)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_wall_stag     !< 1/dz_wall_stag
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_wall_m         !< t_wall prognostic array 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_wall_m         !< t_wall prognostic array
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw                !< wall layer depths (m)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_window      !< volumetric heat capacity of the window material ( J m-3 K-1 ) (= 2.19E6)
@@ -391 +398 @@
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_window_stag    !< window grid spacing (edge-edge)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_window_stag   !< 1/dz_window_stag
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_window_m       !< t_window prognostic array 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_window_m       !< t_window prognostic array
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_window         !< window layer depths (m)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_green       !< volumetric heat capacity of the green material ( J m-3 K-1 ) (= 2.19E6)
@@ -399 +406 @@
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_green_stag     !< green grid spacing (edge-edge)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_green_stag    !< 1/dz_green_stag
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_green_m        !< t_green prognostic array 
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_green_m        !< t_green prognostic array
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_green          !< green layer depths (m)
 
@@ -421 +428 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
-       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface  
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_av       !< average of wghf_eb
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_window_av  !< average of wghf_eb window
@@ -437 +444 @@
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm_av       !< average of theta_10cm (K)
-       
+
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_wall_av      !< Average of t_wall
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_window_av    !< Average of t_window
@@ -455 +462 @@
     TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_usm_v  !< vertical urban surfaces (North, South, East, West)
 
-    INTEGER(iwp), PARAMETER ::  ind_veg_wall  = 0            !< index for vegetation / wall-surface fraction, used for access of albedo, emissivity, etc., for each surface type   
+    INTEGER(iwp), PARAMETER ::  ind_veg_wall  = 0            !< index for vegetation / wall-surface fraction, used for access of albedo, emissivity, etc., for each surface type
     INTEGER(iwp), PARAMETER ::  ind_pav_green = 1            !< index for pavement / green-wall surface fraction, used for access of albedo, emissivity, etc., for each surface type
     INTEGER(iwp), PARAMETER ::  ind_wat_win   = 2            !< index for water / window-surface fraction, used for access of albedo, emissivity, etc., for each surface type
 
-    INTEGER(iwp) ::  ns_h_on_file(0:2)                       !< total number of horizontal surfaces with the same facing, required for writing restart data 
-    INTEGER(iwp) ::  ns_v_on_file(0:3)                       !< total number of vertical surfaces with the same facing, required for writing restart data 
-
-    LOGICAL ::  vertical_surfaces_exist = .FALSE.   !< flag indicating that there are vertical urban/land surfaces 
+    INTEGER(iwp) ::  ns_h_on_file(0:2)                       !< total number of horizontal surfaces with the same facing, required for writing restart data
+    INTEGER(iwp) ::  ns_v_on_file(0:3)                       !< total number of vertical surfaces with the same facing, required for writing restart data
+
+    LOGICAL ::  vertical_surfaces_exist = .FALSE.   !< flag indicating that there are vertical urban/land surfaces
                                                     !< in the domain (required to activiate RTM)
 
@@ -468 +475 @@
     LOGICAL ::  surf_microphysics_morrison = .FALSE.   !< use 2-moment Morrison (add. prog. eq. for nc and qc)
     LOGICAL ::  surf_microphysics_seifert = .FALSE.    !< use 2-moment Seifert and Beheng scheme
+    LOGICAL ::  surf_microphysics_ice_extension = .FALSE. !< use 2-moment Seifert and Beheng scheme
 
 
@@ -473 +481 @@
 
     PRIVATE
-    
+
     INTERFACE init_bc
        MODULE PROCEDURE init_bc
@@ -481 +489 @@
        MODULE PROCEDURE init_single_surface_properties
     END INTERFACE init_single_surface_properties
-    
+
     INTERFACE init_surfaces
        MODULE PROCEDURE init_surfaces
@@ -522 +530 @@
            surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,      &
            vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison,             &
-           surf_microphysics_seifert
+           surf_microphysics_seifert, surf_microphysics_ice_extension
 !
 !-- Public subroutines and functions
@@ -544 +552 @@
 ! Description:
 ! ------------
-!> Initialize data type for setting boundary conditions at horizontal and  
-!> vertical surfaces. 
+!> Initialize data type for setting boundary conditions at horizontal and
+!> vertical surfaces.
 !------------------------------------------------------------------------------!
     SUBROUTINE init_bc
@@ -559 +567 @@
        INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of horizontal surfaces at (j,i)-grid point
        INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index of horizontal surface elements
-       
+
        INTEGER(iwp), DIMENSION(0:3) ::  num_v         !< number of vertical surfaces on subdomain
        INTEGER(iwp), DIMENSION(0:3) ::  num_v_kji     !< number of vertical surfaces at (j,i)-grid point
        INTEGER(iwp), DIMENSION(0:3) ::  start_index_v !< local start index of vertical surface elements
 !
-!--    Set offset indices, i.e. index difference between surface element and 
+!--    Set offset indices, i.e. index difference between surface element and
 !--    surface-bounded grid point.
 !--    Horizontal surfaces - no horizontal offsets
@@ -596 +604 @@
 !
 !--    Initialize data structure for horizontal surfaces, i.e. count the number
-!--    of surface elements, allocate and initialize the respective index arrays, 
-!--    and set the respective start and end indices at each (j,i)-location. 
-!--    The index space is defined also over the ghost points, so that e.g. 
-!--    boundary conditions for diagnostic quanitities can be set on ghost  
-!--    points so that no exchange is required any more. 
+!--    of surface elements, allocate and initialize the respective index arrays,
+!--    and set the respective start and end indices at each (j,i)-location.
+!--    The index space is defined also over the ghost points, so that e.g.
+!--    boundary conditions for diagnostic quanitities can be set on ghost
+!--    points so that no exchange is required any more.
        DO  l = 0, 1
 !
@@ -608 +616 @@
              DO  j = nysg, nyng
                 DO  k = nzb+1, nzt
-!         
+!
 !--                Check if current gridpoint belongs to the atmosphere
                    IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
@@ -619 +627 @@
              ENDDO
           ENDDO
-!         
+!
 !--       Save the number of horizontal surface elements
           bc_h(l)%ns = num_h(l)
@@ -631 +639 @@
           bc_h(l)%start_index = 1
           bc_h(l)%end_index   = 0
-          
+
           num_h(l)         = 1
           start_index_h(l) = 1
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-          
+
                 num_h_kji(l) = 0
                 DO  k = nzb+1, nzt
-!         
+!
 !--                Check if current gridpoint belongs to the atmosphere
                    IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
-!         
+!
 !--                   Upward-facing
                       IF ( .NOT. BTEST( wall_flags_total_0(k+bc_h(l)%koff,     &
@@ -666 +674 @@
 !
 !--    Initialize data structure for vertical surfaces, i.e. count the number
-!--    of surface elements, allocate and initialize the respective index arrays, 
+!--    of surface elements, allocate and initialize the respective index arrays,
 !--    and set the respective start and end indices at each (j,i)-location.
        DO  l = 0, 3
@@ -675 +683 @@
              DO  j = nys, nyn
                 DO  k = nzb+1, nzt
-!         
+!
 !--                Check if current gridpoint belongs to the atmosphere
                    IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
@@ -686 +694 @@
              ENDDO
           ENDDO
-!         
+!
 !--       Save the number of horizontal surface elements
           bc_v(l)%ns = num_v(l)
 !
 !--       ALLOCATE arrays for horizontal surfaces. In contrast to the
-!--       horizontal surfaces, the index space is not defined over the 
-!--       ghost points. 
+!--       horizontal surfaces, the index space is not defined over the
+!--       ghost points.
           ALLOCATE( bc_v(l)%i(1:bc_v(l)%ns) )
           ALLOCATE( bc_v(l)%j(1:bc_v(l)%ns) )
@@ -700 +708 @@
           bc_v(l)%start_index = 1
           bc_v(l)%end_index   = 0
-          
+
           num_v(l)         = 1
           start_index_v(l) = 1
           DO  i = nxl, nxr
              DO  j = nys, nyn
-          
+
                 num_v_kji(l) = 0
                 DO  k = nzb+1, nzt
-!         
+!
 !--                Check if current gridpoint belongs to the atmosphere
                    IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
-!         
+!
 !--                   Upward-facing
                       IF ( .NOT. BTEST( wall_flags_total_0(k+bc_v(l)%koff,     &
@@ -741 +749 @@
 ! ------------
 !> Initialize horizontal and vertical surfaces. Counts the number of default-,
-!> natural and urban surfaces and allocates memory, respectively. 
+!> natural and urban surfaces and allocates memory, respectively.
 !------------------------------------------------------------------------------!
     SUBROUTINE init_surface_arrays
@@ -751 +759 @@
        IMPLICIT NONE
 
-       INTEGER(iwp)                 ::  i         !< running index x-direction 
+       INTEGER(iwp)                 ::  i         !< running index x-direction
        INTEGER(iwp)                 ::  j         !< running index y-direction
        INTEGER(iwp)                 ::  k         !< running index z-direction
        INTEGER(iwp)                 ::  l         !< index variable for surface facing
-       INTEGER(iwp)                 ::  num_lsm_h !< number of horizontally-aligned natural surfaces 
-       INTEGER(iwp)                 ::  num_usm_h !< number of horizontally-aligned urban surfaces 
-
-       INTEGER(iwp), DIMENSION(0:2) ::  num_def_h !< number of horizontally-aligned default surfaces 
-       INTEGER(iwp), DIMENSION(0:3) ::  num_def_v !< number of vertically-aligned default surfaces 
-       INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces 
-       INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces 
+       INTEGER(iwp)                 ::  num_lsm_h !< number of horizontally-aligned natural surfaces
+       INTEGER(iwp)                 ::  num_usm_h !< number of horizontally-aligned urban surfaces
+
+       INTEGER(iwp), DIMENSION(0:2) ::  num_def_h !< number of horizontally-aligned default surfaces
+       INTEGER(iwp), DIMENSION(0:3) ::  num_def_v !< number of vertically-aligned default surfaces
+       INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces
+       INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces
 
        INTEGER(iwp)              ::  num_surf_v_l !< number of vertically-aligned local urban/land surfaces
@@ -768 +776 @@
        LOGICAL ::  building                       !< flag indicating building grid point
        LOGICAL ::  terrain                        !< flag indicating natural terrain grid point
-       LOGICAL ::  unresolved_building            !< flag indicating a grid point where actually a building is 
-                                                  !< defined but not resolved by the vertical grid 
+       LOGICAL ::  unresolved_building            !< flag indicating a grid point where actually a building is
+                                                  !< defined but not resolved by the vertical grid
 
        num_def_h = 0
@@ -780 +788 @@
 !--    Surfaces are classified according to the input data read from static
 !--    input file. If no input file is present, all surfaces are classified
-!--    either as natural, urban, or default, depending on the setting of 
+!--    either as natural, urban, or default, depending on the setting of
 !--    land_surface and urban_surface. To control this, use the control
 !--    flag topo_no_distinct
 !
-!--    Count number of horizontal surfaces on local domain 
+!--    Count number of horizontal surfaces on local domain
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -798 +806 @@
 !
 !--                   Determine flags indicating a terrain surface, a building
-!--                   surface, 
+!--                   surface,
                       terrain  = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .OR.       &
                                  topo_no_distinct
@@ -804 +812 @@
                                  topo_no_distinct
 !
-!--                   unresolved_building indicates a surface with equal height 
+!--                   unresolved_building indicates a surface with equal height
 !--                   as terrain but with a non-grid resolved building on top.
 !--                   These surfaces will be flagged as urban surfaces.
@@ -813 +821 @@
                       IF ( land_surface  .AND.  terrain  .AND.                 &
                            .NOT. unresolved_building )  THEN
-                         num_lsm_h    = num_lsm_h    + 1 
+                         num_lsm_h    = num_lsm_h    + 1
 !
 !--                   Urban surface tpye
                       ELSEIF ( urban_surface  .AND.  building )  THEN
-                         num_usm_h    = num_usm_h    + 1 
+                         num_usm_h    = num_usm_h    + 1
 !
 !--                   Default-surface type
                       ELSEIF ( .NOT. land_surface    .AND.                     &
                                .NOT. urban_surface )  THEN
-                               
+
                          num_def_h(0) = num_def_h(0) + 1
 !
 !--                   Unclassifified surface-grid point. Give error message.
-                      ELSE 
+                      ELSE
                          WRITE( message_string, * )                           &
                                           'Unclassified upward-facing ' //    &
@@ -840 +848 @@
                       num_def_h(2) = num_def_h(2) + 1
 !
-!--                Check for any other downward-facing surface. So far only for 
+!--                Check for any other downward-facing surface. So far only for
 !--                default surface type.
                    ELSEIF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
                       num_def_h(1) = num_def_h(1) + 1
-                   ENDIF 
+                   ENDIF
 
                 ENDIF
@@ -851 +859 @@
        ENDDO
 !
-!--    Count number of vertical surfaces on local domain 
+!--    Count number of vertical surfaces on local domain
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -869 +877 @@
                       unresolved_building = BTEST( wall_flags_total_0(k,j-1,i), 5 )  &
                                      .AND.  BTEST( wall_flags_total_0(k,j-1,i), 6 )
-                                     
+
                       IF (  land_surface  .AND.  terrain  .AND.                &
                            .NOT. unresolved_building )  THEN
-                         num_lsm_v(0) = num_lsm_v(0) + 1 
+                         num_lsm_v(0) = num_lsm_v(0) + 1
                       ELSEIF ( urban_surface  .AND.  building )  THEN
-                         num_usm_v(0) = num_usm_v(0) + 1 
+                         num_usm_v(0) = num_usm_v(0) + 1
 !
 !--                   Default-surface type
                       ELSEIF ( .NOT. land_surface    .AND.                     &
                                .NOT. urban_surface )  THEN
-                         num_def_v(0) = num_def_v(0) + 1 
+                         num_def_v(0) = num_def_v(0) + 1
 !
 !--                   Unclassifified surface-grid point. Give error message.
-                      ELSE 
+                      ELSE
                          WRITE( message_string, * )                            &
                                           'Unclassified northward-facing ' //  &
@@ -900 +908 @@
                       building = BTEST( wall_flags_total_0(k,j+1,i), 6 )  .OR.       &
                                  topo_no_distinct
-                                 
+
                       unresolved_building = BTEST( wall_flags_total_0(k,j+1,i), 5 )  &
-                                     .AND.  BTEST( wall_flags_total_0(k,j+1,i), 6 ) 
-                                
+                                     .AND.  BTEST( wall_flags_total_0(k,j+1,i), 6 )
+
                       IF (  land_surface  .AND.  terrain  .AND.                &
                            .NOT. unresolved_building )  THEN
-                         num_lsm_v(1) = num_lsm_v(1) + 1 
+                         num_lsm_v(1) = num_lsm_v(1) + 1
                       ELSEIF ( urban_surface  .AND.  building )  THEN
-                         num_usm_v(1) = num_usm_v(1) + 1 
+                         num_usm_v(1) = num_usm_v(1) + 1
 !
 !--                   Default-surface type
                       ELSEIF ( .NOT. land_surface    .AND.                     &
                                .NOT. urban_surface )  THEN
-                         num_def_v(1) = num_def_v(1) + 1 
+                         num_def_v(1) = num_def_v(1) + 1
 !
 !--                   Unclassifified surface-grid point. Give error message.
-                      ELSE 
+                      ELSE
                          WRITE( message_string, * )                            &
                                           'Unclassified southward-facing ' //  &
@@ -934 +942 @@
                       building = BTEST( wall_flags_total_0(k,j,i-1), 6 )  .OR.       &
                                  topo_no_distinct
-                                 
+
                       unresolved_building = BTEST( wall_flags_total_0(k,j,i-1), 5 )  &
                                      .AND.  BTEST( wall_flags_total_0(k,j,i-1), 6 )
-                                     
+
                       IF (  land_surface  .AND.  terrain  .AND.                &
                            .NOT. unresolved_building )  THEN
-                         num_lsm_v(2) = num_lsm_v(2) + 1 
+                         num_lsm_v(2) = num_lsm_v(2) + 1
                       ELSEIF ( urban_surface  .AND.  building )  THEN
-                         num_usm_v(2) = num_usm_v(2) + 1 
+                         num_usm_v(2) = num_usm_v(2) + 1
 !
 !--                   Default-surface type
                       ELSEIF ( .NOT. land_surface    .AND.                     &
                                .NOT. urban_surface )  THEN
-                         num_def_v(2) = num_def_v(2) + 1 
+                         num_def_v(2) = num_def_v(2) + 1
 !
 !--                   Unclassifified surface-grid point. Give error message.
-                      ELSE 
+                      ELSE
                          WRITE( message_string, * )                            &
                                           'Unclassified eastward-facing ' //   &
@@ -968 +976 @@
                       building = BTEST( wall_flags_total_0(k,j,i+1), 6 )  .OR.       &
                                  topo_no_distinct
-                                 
+
                       unresolved_building = BTEST( wall_flags_total_0(k,j,i+1), 5 )  &
                                      .AND.  BTEST( wall_flags_total_0(k,j,i+1), 6 )
-                                 
+
                       IF (  land_surface  .AND.  terrain  .AND.                &
                            .NOT. unresolved_building )  THEN
-                         num_lsm_v(3) = num_lsm_v(3) + 1 
+                         num_lsm_v(3) = num_lsm_v(3) + 1
                       ELSEIF ( urban_surface  .AND.  building )  THEN
-                         num_usm_v(3) = num_usm_v(3) + 1 
+                         num_usm_v(3) = num_usm_v(3) + 1
 !
 !--                   Default-surface type
                       ELSEIF ( .NOT. land_surface    .AND.                     &
                                .NOT. urban_surface )  THEN
-                         num_def_v(3) = num_def_v(3) + 1 
+                         num_def_v(3) = num_def_v(3) + 1
 !
 !--                   Unclassifified surface-grid point. Give error message.
-                      ELSE 
+                      ELSE
                          WRITE( message_string, * )                            &
                                           'Unclassified westward-facing ' //   &
@@ -1010 +1018 @@
 !
 !--    Horizontal surface, natural type, so far only upward-facing
-       surf_lsm_h%ns    = num_lsm_h  
+       surf_lsm_h%ns    = num_lsm_h
 !
 !--    Horizontal surface, urban type, so far only upward-facing
-       surf_usm_h%ns    = num_usm_h    
+       surf_usm_h%ns    = num_usm_h
 !
 !--    Vertical surface, default type, northward facing
@@ -1051 +1059 @@
        surf_usm_v(3)%ns = num_usm_v(3)
 !
-!--    Allocate required attributes for horizontal surfaces - default type. 
+!--    Allocate required attributes for horizontal surfaces - default type.
 !--    Upward-facing (l=0) and downward-facing (l=1).
        DO  l = 0, 1
@@ -1060 +1068 @@
        CALL allocate_surface_attributes_h_top ( surf_def_h(2), nys, nyn, nxl, nxr )
 !
-!--    Allocate required attributes for horizontal surfaces - natural type. 
+!--    Allocate required attributes for horizontal surfaces - natural type.
        CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
 !
-!--    Allocate required attributes for horizontal surfaces - urban type. 
+!--    Allocate required attributes for horizontal surfaces - urban type.
        CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
 
 !
-!--    Allocate required attributes for vertical surfaces. 
+!--    Allocate required attributes for vertical surfaces.
 !--    Northward-facing (l=0), southward-facing (l=1), eastward-facing (l=2)
 !--    and westward-facing (l=3).
@@ -1101 +1109 @@
 #endif
        IF ( num_surf_v > 0 )  vertical_surfaces_exist = .TRUE.
-        
+
 
     END SUBROUTINE init_surface_arrays
@@ -1117 +1125 @@
 
        INTEGER(iwp) ::  l     !<
-       
+
        !$ACC ENTER DATA &
        !$ACC COPYIN(surf_def_h(0:2)) &
@@ -1162 +1170 @@
 
        INTEGER(iwp) ::  l     !<
-       
+
        ! Delete data in surf_def_h(0:2)
        DO  l = 0, 1
@@ -1199 +1207 @@
 ! Description:
 ! ------------
-!> Deallocating memory for upward and downward-facing horizontal surface types, 
-!> except for top fluxes. 
+!> Deallocating memory for upward and downward-facing horizontal surface types,
+!> except for top fluxes.
 !------------------------------------------------------------------------------!
     SUBROUTINE deallocate_surface_attributes_h( surfaces )
@@ -1242 +1250 @@
 !
 !--    Vertical momentum fluxes of u and v
-       DEALLOCATE ( surfaces%usws )  
-       DEALLOCATE ( surfaces%vsws )  
+       DEALLOCATE ( surfaces%usws )
+       DEALLOCATE ( surfaces%vsws )
 !
 !--    Required in production_e
-       IF ( .NOT. constant_diffusion )  THEN    
-          DEALLOCATE ( surfaces%u_0 )  
+       IF ( .NOT. constant_diffusion )  THEN
+          DEALLOCATE ( surfaces%u_0 )
           DEALLOCATE ( surfaces%v_0 )
-       ENDIF 
+       ENDIF
 !
 !--    Characteristic temperature and surface flux of sensible heat
-       DEALLOCATE ( surfaces%ts )    
+       DEALLOCATE ( surfaces%ts )
        DEALLOCATE ( surfaces%shf )
 !
 !--    surface temperature
-       DEALLOCATE ( surfaces%pt_surface ) 
+       DEALLOCATE ( surfaces%pt_surface )
 !
 !--    Characteristic humidity and surface flux of latent heat
-       IF ( humidity )  THEN          
-          DEALLOCATE ( surfaces%qs ) 
-          DEALLOCATE ( surfaces%qsws )  
+       IF ( humidity )  THEN
+          DEALLOCATE ( surfaces%qs )
+          DEALLOCATE ( surfaces%qsws )
           DEALLOCATE ( surfaces%q_surface   )
           DEALLOCATE ( surfaces%vpt_surface )
-       ENDIF 
+       ENDIF
 !
 !--    Characteristic scalar and surface flux of scalar
        IF ( passive_scalar )  THEN
-          DEALLOCATE ( surfaces%ss )   
-          DEALLOCATE ( surfaces%ssws ) 
-       ENDIF 
+          DEALLOCATE ( surfaces%ss )
+          DEALLOCATE ( surfaces%ssws )
+       ENDIF
 !
 !--    Scaling parameter (cs*) and surface flux of chemical species
        IF ( air_chemistry )  THEN
-          DEALLOCATE ( surfaces%css )   
-          DEALLOCATE ( surfaces%cssws ) 
-       ENDIF 
+          DEALLOCATE ( surfaces%css )
+          DEALLOCATE ( surfaces%cssws )
+       ENDIF
 !
 !--    Arrays for storing potential temperature and
@@ -1283 +1291 @@
        DEALLOCATE ( surfaces%qv1 )
        DEALLOCATE ( surfaces%vpt1 )
-       
-!
-!--       
+
+!
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
           DEALLOCATE ( surfaces%qcs )
@@ -1293 +1301 @@
        ENDIF
 !
-!--       
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
           DEALLOCATE ( surfaces%qrs )
@@ -1301 +1309 @@
        ENDIF
 !
+!--
+       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
+          DEALLOCATE ( surfaces%qis )
+          DEALLOCATE ( surfaces%nis )
+          DEALLOCATE ( surfaces%qisws )
+          DEALLOCATE ( surfaces%nisws )
+       ENDIF
+!
 !--    Salinity surface flux
        IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
@@ -1310 +1326 @@
 ! Description:
 ! ------------
-!> Allocating memory for upward and downward-facing horizontal surface types, 
-!> except for top fluxes. 
+!> Allocating memory for upward and downward-facing horizontal surface types,
+!> except for top fluxes.
 !------------------------------------------------------------------------------!
     SUBROUTINE allocate_surface_attributes_h( surfaces,                        &
@@ -1326 +1342 @@
 
 !
-!--    Allocate arrays for start and end index of horizontal surface type 
+!--    Allocate arrays for start and end index of horizontal surface type
 !--    for each (j,i)-grid point. This is required e.g. in diffion_x, which is
-!--    called for each (j,i). In order to find the location where the 
-!--    respective flux is store within the surface-type, start- and end- 
+!--    called for each (j,i). In order to find the location where the
+!--    respective flux is store within the surface-type, start- and end-
 !--    index are stored for each (j,i). For example, each (j,i) can have
 !--    several entries where fluxes for horizontal surfaces might be stored,
@@ -1370 +1386 @@
 !
 !--    Vertical momentum fluxes of u and v
-       ALLOCATE ( surfaces%usws(1:surfaces%ns) )  
-       ALLOCATE ( surfaces%vsws(1:surfaces%ns) )  
+       ALLOCATE ( surfaces%usws(1:surfaces%ns) )
+       ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
 !
 !--    Required in production_e
-       IF ( .NOT. constant_diffusion )  THEN    
-          ALLOCATE ( surfaces%u_0(1:surfaces%ns) )  
+       IF ( .NOT. constant_diffusion )  THEN
+          ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
           ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
-       ENDIF 
+       ENDIF
 !
 !--    Characteristic temperature and surface flux of sensible heat
-       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )    
+       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )
        ALLOCATE ( surfaces%shf(1:surfaces%ns) )
 !
 !--    Surface temperature
-       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) ) 
+       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
 !
 !--    Characteristic humidity, surface flux of latent heat, and surface virtual potential temperature
        IF ( humidity )  THEN
-          ALLOCATE ( surfaces%qs(1:surfaces%ns)   ) 
-          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )      
-          ALLOCATE ( surfaces%q_surface(1:surfaces%ns)   ) 
-          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) ) 
-       ENDIF 
+          ALLOCATE ( surfaces%qs(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
+          ALLOCATE ( surfaces%q_surface(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )
+       ENDIF
 
 !
 !--    Characteristic scalar and surface flux of scalar
        IF ( passive_scalar )  THEN
-          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )   
-          ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
-       ENDIF 
+          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
+       ENDIF
 !
 !--    Scaling parameter (cs*) and surface flux of chemical species
        IF ( air_chemistry )  THEN
-          ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )   
-          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
-       ENDIF 
+          ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )
+          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
+       ENDIF
 !
 !--    Arrays for storing potential temperature and
@@ -1413 +1429 @@
        ALLOCATE ( surfaces%vpt1(1:surfaces%ns) )
 !
-!--       
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
           ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
@@ -1421 +1437 @@
        ENDIF
 !
-!--       
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
           ALLOCATE ( surfaces%qrs(1:surfaces%ns)   )
@@ -1428 +1444 @@
           ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
        ENDIF
+
+!
+!--
+       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
+          ALLOCATE ( surfaces%qis(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%nis(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
+          ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
+       ENDIF
+
 !
 !--    Salinity surface flux
@@ -1445 +1471 @@
 
        IMPLICIT NONE
-   
+
        TYPE(surf_type) ::  surfaces  !< respective surface type
-   
+
        !$ACC EXIT DATA &
        !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
@@ -1473 +1499 @@
           !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
        ENDIF
-   
+
     END SUBROUTINE exit_surface_attributes_h
 #endif
@@ -1522 +1548 @@
 ! Description:
 ! ------------
-!> Deallocating memory for model-top fluxes  
+!> Deallocating memory for model-top fluxes
 !------------------------------------------------------------------------------!
     SUBROUTINE deallocate_surface_attributes_h_top( surfaces )
@@ -1539 +1565 @@
        DEALLOCATE ( surfaces%k )
 
-       IF ( .NOT. constant_diffusion )  THEN    
-          DEALLOCATE ( surfaces%u_0 )  
+       IF ( .NOT. constant_diffusion )  THEN
+          DEALLOCATE ( surfaces%u_0 )
           DEALLOCATE ( surfaces%v_0 )
-       ENDIF 
+       ENDIF
 !
 !--    Vertical momentum fluxes of u and v
-       DEALLOCATE ( surfaces%usws )  
-       DEALLOCATE ( surfaces%vsws )  
+       DEALLOCATE ( surfaces%usws )
+       DEALLOCATE ( surfaces%vsws )
 !
 !--    Sensible heat flux
@@ -1553 +1579 @@
 !--    Latent heat flux
        IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere')  THEN
-          DEALLOCATE ( surfaces%qsws )      
-       ENDIF 
+          DEALLOCATE ( surfaces%qsws )
+       ENDIF
 !
 !--    Scalar flux
        IF ( passive_scalar )  THEN
-          DEALLOCATE ( surfaces%ssws ) 
-       ENDIF 
+          DEALLOCATE ( surfaces%ssws )
+       ENDIF
 !
 !--    Chemical species flux
        IF ( air_chemistry )  THEN
-          DEALLOCATE ( surfaces%cssws ) 
-       ENDIF 
-!
-!--       
+          DEALLOCATE ( surfaces%cssws )
+       ENDIF
+!
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
           DEALLOCATE ( surfaces%qcsws )
@@ -1572 +1598 @@
        ENDIF
 !
-!--       
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
           DEALLOCATE ( surfaces%qrsws )
           DEALLOCATE ( surfaces%nrsws )
        ENDIF
+
+!
+!--
+       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
+          DEALLOCATE ( surfaces%qisws )
+          DEALLOCATE ( surfaces%nisws )
+       ENDIF
 !
 !--    Salinity flux
@@ -1587 +1620 @@
 ! Description:
 ! ------------
-!> Allocating memory for model-top fluxes  
+!> Allocating memory for model-top fluxes
 !------------------------------------------------------------------------------!
     SUBROUTINE allocate_surface_attributes_h_top( surfaces,                    &
@@ -1611 +1644 @@
        ALLOCATE ( surfaces%k(1:surfaces%ns)  )
 
-       IF ( .NOT. constant_diffusion )  THEN    
-          ALLOCATE ( surfaces%u_0(1:surfaces%ns) )  
+       IF ( .NOT. constant_diffusion )  THEN
+          ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
           ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
-       ENDIF 
+       ENDIF
 !
 !--    Vertical momentum fluxes of u and v
-       ALLOCATE ( surfaces%usws(1:surfaces%ns) )  
-       ALLOCATE ( surfaces%vsws(1:surfaces%ns) )  
+       ALLOCATE ( surfaces%usws(1:surfaces%ns) )
+       ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
 !
 !--    Sensible heat flux
@@ -1625 +1658 @@
 !--    Latent heat flux
        IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere')  THEN
-          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )      
-       ENDIF 
+          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
+       ENDIF
 !
 !--    Scalar flux
        IF ( passive_scalar )  THEN
-          ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
-       ENDIF 
+          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
+       ENDIF
 !
 !--    Chemical species flux
        IF ( air_chemistry )  THEN
-          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
-       ENDIF 
-!
-!--       
+          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
+       ENDIF
+!
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
           ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
@@ -1644 +1677 @@
        ENDIF
 !
-!--       
+!--
        IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
           ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
           ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
+       ENDIF
+
+!
+!--
+       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
+          ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
+          ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
        ENDIF
 !
@@ -1665 +1705 @@
 
        IMPLICIT NONE
-   
+
        TYPE(surf_type) ::  surfaces  !< respective surface type
-   
+
        !$ACC EXIT DATA &
        !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
@@ -1683 +1723 @@
           !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
        ENDIF
-   
+
     END SUBROUTINE exit_surface_attributes_h_top
 #endif
@@ -1721 +1761 @@
 ! Description:
 ! ------------
-!> Deallocating memory for vertical surface types. 
+!> Deallocating memory for vertical surface types.
 !------------------------------------------------------------------------------!
     SUBROUTINE deallocate_surface_attributes_v( surfaces )
@@ -1731 +1771 @@
 
 !
-!--    Allocate arrays for start and end index of vertical surface type 
+!--    Allocate arrays for start and end index of vertical surface type
 !--    for each (j,i)-grid point. This is required in diffion_x, which is
-!--    called for each (j,i). In order to find the location where the 
-!--    respective flux is store within the surface-type, start- and end- 
+!--    called for each (j,i). In order to find the location where the
+!--    respective flux is store within the surface-type, start- and end-
 !--    index are stored for each (j,i). For example, each (j,i) can have
-!--    several entries where fluxes for vertical surfaces might be stored.  
-!--    In the flat case, where no vertical walls exit, set indicies such 
-!--    that loop in diffusion routines will not be entered. 
+!--    several entries where fluxes for vertical surfaces might be stored.
+!--    In the flat case, where no vertical walls exit, set indicies such
+!--    that loop in diffusion routines will not be entered.
        DEALLOCATE ( surfaces%start_index )
        DEALLOCATE ( surfaces%end_index )
@@ -1766 +1806 @@
 !
 !--    Allocate Obukhov length and bulk Richardson number. Actually, at
-!--    vertical surfaces these are only required for natural surfaces.  
+!--    vertical surfaces these are only required for natural surfaces.
 !--    for natural land surfaces
-       DEALLOCATE( surfaces%ol ) 
-       DEALLOCATE( surfaces%rib ) 
-!
-!--    Allocate arrays for surface momentum fluxes for u and v. For u at north- 
+       DEALLOCATE( surfaces%ol )
+       DEALLOCATE( surfaces%rib )
+!
+!--    Allocate arrays for surface momentum fluxes for u and v. For u at north-
 !--    and south-facing surfaces, for v at east- and west-facing surfaces.
        DEALLOCATE ( surfaces%mom_flux_uv )
 !
 !--    Allocate array for surface momentum flux for w - wsus and wsvs
-       DEALLOCATE ( surfaces%mom_flux_w ) 
-!
-!--    Allocate array for surface momentum flux for subgrid-scale tke wsus and 
+       DEALLOCATE ( surfaces%mom_flux_w )
+!
+!--    Allocate array for surface momentum flux for subgrid-scale tke wsus and
 !--    wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
 !--    on surface.
-       DEALLOCATE ( surfaces%mom_flux_tke )  
+       DEALLOCATE ( surfaces%mom_flux_tke )
 !
 !--    Characteristic temperature and surface flux of sensible heat
-       DEALLOCATE ( surfaces%ts )    
+       DEALLOCATE ( surfaces%ts )
        DEALLOCATE ( surfaces%shf )
 !
 !--    surface temperature
-       DEALLOCATE ( surfaces%pt_surface ) 
+       DEALLOCATE ( surfaces%pt_surface )
 !
 !--    Characteristic humidity and surface flux of latent heat
        IF ( humidity )  THEN
-          DEALLOCATE ( surfaces%qs ) 
-          DEALLOCATE ( surfaces%qsws )  
+          DEALLOCATE ( surfaces%qs )
+          DEALLOCATE ( surfaces%qsws )
           DEALLOCATE ( surfaces%q_surface   )
           DEALLOCATE ( surfaces%vpt_surface )
-       ENDIF 
+       ENDIF
 !
 !--    Characteristic scalar and surface flux of scalar
        IF ( passive_scalar )  THEN
-          DEALLOCATE ( surfaces%ss )   
-          DEALLOCATE ( surfaces%ssws ) 
+          DEALLOCATE ( surfaces%ss )
+          DEALLOCATE ( surfaces%ssws )
        ENDIF
 !
 !--    Scaling parameter (cs*) and surface flux of chemical species
        IF ( air_chemistry )  THEN
-             DEALLOCATE ( surfaces%css )   
-             DEALLOCATE ( surfaces%cssws ) 
-       ENDIF 
+             DEALLOCATE ( surfaces%css )
+             DEALLOCATE ( surfaces%cssws )
+       ENDIF
 !
 !--    Arrays for storing potential temperature and
@@ -1829 +1869 @@
           DEALLOCATE ( surfaces%nrsws )
        ENDIF
+
+       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
+          DEALLOCATE ( surfaces%qis )
+          DEALLOCATE ( surfaces%nis )
+          DEALLOCATE ( surfaces%qisws )
+          DEALLOCATE ( surfaces%nisws )
+       ENDIF
+
 !
 !--    Salinity surface flux
@@ -1839 +1887 @@
 ! Description:
 ! ------------
-!> Allocating memory for vertical surface types. 
+!> Allocating memory for vertical surface types.
 !------------------------------------------------------------------------------!
     SUBROUTINE allocate_surface_attributes_v( surfaces,                        &
@@ -1854 +1902 @@
 
 !
-!--    Allocate arrays for start and end index of vertical surface type 
+!--    Allocate arrays for start and end index of vertical surface type
 !--    for each (j,i)-grid point. This is required in diffion_x, which is
-!--    called for each (j,i). In order to find the location where the 
-!--    respective flux is store within the surface-type, start- and end- 
+!--    called for each (j,i). In order to find the location where the
+!--    respective flux is store within the surface-type, start- and end-
 !--    index are stored for each (j,i). For example, each (j,i) can have
-!--    several entries where fluxes for vertical surfaces might be stored.  
-!--    In the flat case, where no vertical walls exit, set indicies such 
-!--    that loop in diffusion routines will not be entered. 
+!--    several entries where fluxes for vertical surfaces might be stored.
+!--    In the flat case, where no vertical walls exit, set indicies such
+!--    that loop in diffusion routines will not be entered.
        ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
        ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l)   )
@@ -1891 +1939 @@
 !
 !--    Allocate Obukhov length and bulk Richardson number. Actually, at
-!--    vertical surfaces these are only required for natural surfaces.  
+!--    vertical surfaces these are only required for natural surfaces.
 !--    for natural land surfaces
-       ALLOCATE( surfaces%ol(1:surfaces%ns)  ) 
-       ALLOCATE( surfaces%rib(1:surfaces%ns) ) 
-!
-!--    Allocate arrays for surface momentum fluxes for u and v. For u at north- 
+       ALLOCATE( surfaces%ol(1:surfaces%ns)  )
+       ALLOCATE( surfaces%rib(1:surfaces%ns) )
+!
+!--    Allocate arrays for surface momentum fluxes for u and v. For u at north-
 !--    and south-facing surfaces, for v at east- and west-facing surfaces.
        ALLOCATE ( surfaces%mom_flux_uv(1:surfaces%ns) )
 !
 !--    Allocate array for surface momentum flux for w - wsus and wsvs
-       ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) ) 
-!
-!--    Allocate array for surface momentum flux for subgrid-scale tke wsus and 
+       ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) )
+!
+!--    Allocate array for surface momentum flux for subgrid-scale tke wsus and
 !--    wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
 !--    on surface.
-       ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) )  
+       ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) )
 !
 !--    Characteristic temperature and surface flux of sensible heat
-       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )    
+       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )
        ALLOCATE ( surfaces%shf(1:surfaces%ns) )
 !
 !--    surface temperature
-       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) ) 
+       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
 !
 !--    Characteristic humidity and surface flux of latent heat
        IF ( humidity )  THEN
-          ALLOCATE ( surfaces%qs(1:surfaces%ns)          ) 
-          ALLOCATE ( surfaces%qsws(1:surfaces%ns)        )   
+          ALLOCATE ( surfaces%qs(1:surfaces%ns)          )
+          ALLOCATE ( surfaces%qsws(1:surfaces%ns)        )
           ALLOCATE ( surfaces%q_surface(1:surfaces%ns)   )
-          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )            
-       ENDIF 
+          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )
+       ENDIF
 !
 !--    Characteristic scalar and surface flux of scalar
        IF ( passive_scalar )  THEN
-          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )   
-          ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
+          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
        ENDIF
 !
 !--    Scaling parameter (cs*) and surface flux of chemical species
        IF ( air_chemistry )  THEN
-             ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )   
-             ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
-       ENDIF 
+             ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )
+             ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
+       ENDIF
 !
 !--    Arrays for storing potential temperature and
@@ -1954 +2002 @@
           ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
        ENDIF
+
+       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
+          ALLOCATE ( surfaces%qis(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%nis(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
+          ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
+       ENDIF
 !
 !--    Salinity surface flux
@@ -1964 +2019 @@
 ! Description:
 ! ------------
-!> Exit memory for vertical surface types. 
+!> Exit memory for vertical surface types.
 !------------------------------------------------------------------------------!
 #if defined( _OPENACC )
@@ -1996 +2051 @@
 ! Description:
 ! ------------
-!> Enter memory for vertical surface types. 
+!> Enter memory for vertical surface types.
 !------------------------------------------------------------------------------!
 #if defined( _OPENACC )
     SUBROUTINE enter_surface_attributes_v( surfaces )
-   
+
        IMPLICIT NONE
-   
+
        TYPE(surf_type) ::  surfaces  !< respective surface type
-   
+
        !$ACC ENTER DATA &
        !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
@@ -2021 +2076 @@
        !$ACC COPYIN(surfaces%pt1(1:surfaces%ns)) &
        !$ACC COPYIN(surfaces%qv1(1:surfaces%ns))
-   
+
     END SUBROUTINE enter_surface_attributes_v
 #endif
@@ -2028 +2083 @@
 ! Description:
 ! ------------
-!> Initialize surface elements, i.e. set initial values for surface fluxes, 
-!> friction velocity, calcuation of start/end indices, etc. .