Changeset 4502


Ignore:
Timestamp:
Apr 17, 2020 4:14:16 PM (12 months ago)
Author:
schwenkel
Message:

Implementation of ice microphysics

Location:
palm/trunk/SOURCE
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_s_bc.f90

    r4488 r4502  
    2424! -----------------
    2525! $Id$
     26! Implementation of ice microphysics
     27!
     28! 4488 2020-04-03 11:34:29Z raasch
    2629! file re-formatted to follow the PALM coding standard
    2730!
     
    957960          ENDDO
    958961
    959        ELSEIF ( sk_char == 'qr' )  THEN
    960 
    961 !
    962 !--       Rain water content boundary condition at the bottom boundary: Dirichlet (fixed surface
    963 !--       rain water content).
     962       ELSEIF ( sk_char == 'qi' )  THEN
     963
     964!
     965!--       Ice crystal content boundary condition at the bottom boundary:
     966!--       Dirichlet (fixed surface rain water content).
    964967          DO  i = nxl, nxr
    965968             DO  j = nys, nyn
     
    971974
    972975!
    973 !--       Rain water content boundary condition at the top boundary: Dirichlet
     976!--       Ice crystal content boundary condition at the top boundary: Dirichlet
    974977          DO  i = nxl, nxr
    975978             DO  j = nys, nyn
     
    979982          ENDDO
    980983
    981        ELSEIF ( sk_char == 'nc' )  THEN
    982 
    983 !
    984 !--       Cloud drop concentration boundary condition at the bottom boundary: Dirichlet (fixed
    985 !--       surface cloud drop concentration).
     984       ELSEIF ( sk_char == 'qr' )  THEN
     985
     986!
     987!--       Rain water content boundary condition at the bottom boundary: Dirichlet (fixed surface
     988!--       rain water content).
    986989          DO  i = nxl, nxr
    987990             DO  j = nys, nyn
     
    993996
    994997!
     998!--       Rain water content boundary condition at the top boundary: Dirichlet
     999          DO  i = nxl, nxr
     1000             DO  j = nys, nyn
     1001                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     1002                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     1003             ENDDO
     1004          ENDDO
     1005
     1006       ELSEIF ( sk_char == 'nc' )  THEN
     1007
     1008!
     1009!--       Cloud drop concentration boundary condition at the bottom boundary: Dirichlet (fixed
     1010!--       surface cloud drop concentration).
     1011          DO  i = nxl, nxr
     1012             DO  j = nys, nyn
     1013                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     1014                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     1015                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     1016             ENDDO
     1017          ENDDO
     1018
     1019!
    9951020!--       Cloud drop concentration boundary condition at the top boundary: Dirichlet
     1021          DO  i = nxl, nxr
     1022             DO  j = nys, nyn
     1023                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     1024                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     1025             ENDDO
     1026          ENDDO
     1027
     1028       ELSEIF ( sk_char == 'ni' )  THEN
     1029
     1030!
     1031!--       Ice crystal concentration boundary condition at the bottom boundary:
     1032!--       Dirichlet (fixed surface cloud drop concentration).
     1033          DO  i = nxl, nxr
     1034             DO  j = nys, nyn
     1035                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     1036                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     1037                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     1038             ENDDO
     1039          ENDDO
     1040
     1041!
     1042!--       Ice crystal concentration boundary condition at the top boundary: Dirichlet
    9961043          DO  i = nxl, nxr
    9971044             DO  j = nys, nyn
  • palm/trunk/SOURCE/advec_ws.f90

    r4469 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4469 2020-03-23 14:31:00Z suehring
    2730! fix mistakenly committed version
    2831!
     
    240243               sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l,                &
    241244               sums_wsqcs_ws_l, sums_wsqrs_ws_l,                               &
     245               sums_wsqis_ws_l, sums_wsnis_ws_l,                               &
    242246               sums_wsncs_ws_l, sums_wsnrs_ws_l,                               &
    243247               hom, weight_substep
     
    19381942             ENDDO
    19391943
     1944          CASE ( 'qi' )
     1945
     1946             DO  k = nzb, nzt
     1947                sums_wsqis_ws_l(k,tn)  = sums_wsqis_ws_l(k,tn) +               &
     1948                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     1949                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     1950                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     1951                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
     1952                    ) * weight_substep(intermediate_timestep_count)
     1953             ENDDO
    19401954
    19411955          CASE ( 'qr' )
     
    19541968             DO  k = nzb, nzt
    19551969                sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn) +               &
     1970                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     1971                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     1972                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     1973                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
     1974                    ) * weight_substep(intermediate_timestep_count)
     1975             ENDDO
     1976
     1977          CASE ( 'ni' )
     1978
     1979             DO  k = nzb, nzt
     1980                sums_wsnis_ws_l(k,tn)  = sums_wsnis_ws_l(k,tn) +               &
    19561981                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    19571982                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     
    37983823           CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' )
    37993824              sk_num = 9
     3825           CASE ( 'ni' )
     3826              sk_num = 10
     3827           CASE ( 'qi' )
     3828              sk_num = 11
    38003829           CASE DEFAULT
    38013830              sk_num = 0
     
    38243853       !$ACC PRESENT(sums_wsqrs_ws_l, sums_wsncs_ws_l) &
    38253854       !$ACC PRESENT(sums_wsnrs_ws_l, sums_wsss_ws_l) &
     3855       !$ACC PRESENT(sums_wsnis_ws_l, sums_wsqis_ws_l) &
    38263856       !$ACC PRESENT(sums_salsa_ws_l)
    38273857       DO  i = nxl, nxr
     
    45174547                               *   ABS(w(k,j,i) - hom(k,1,3,0)             )   &
    45184548                           ) * weight_substep(intermediate_timestep_count)
     4549                    CASE ( 10 )
     4550                       !$ACC ATOMIC
     4551                       sums_wsnis_ws_l(k,tn)  = sums_wsnis_ws_l(k,tn)          &
     4552                          + ( flux_t(k)                                        &
     4553                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     4554                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     4555                            + diss_t(k)                                        &
     4556                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     4557                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
     4558                            ) * weight_substep(intermediate_timestep_count)
     4559                    CASE ( 11 )
     4560                       !$ACC ATOMIC
     4561                       sums_wsqis_ws_l(k,tn)  = sums_wsqis_ws_l(k,tn)          &
     4562                          + ( flux_t(k)                                        &
     4563                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     4564                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     4565                            + diss_t(k)                                        &
     4566                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     4567                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
     4568                            ) * weight_substep(intermediate_timestep_count)
    45194569
    45204570                END SELECT
  • palm/trunk/SOURCE/basic_constants_and_equations_mod.f90

    r4400 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4400 2020-02-10 20:32:41Z suehring
    2730! Move routine to transform coordinates from netcdf_interface_mod to
    2831! basic_constants_and_equations_mod
     
    9295    REAL(wp), PARAMETER ::  g_d_cp  = g   / c_p   !< precomputed g / c_p
    9396    REAL(wp), PARAMETER ::  lv_d_cp = l_v / c_p   !< precomputed l_v / c_p
     97    REAL(wp), PARAMETER ::  ls_d_cp = l_s / c_p   !< precomputed l_s / c_p
    9498    REAL(wp), PARAMETER ::  lv_d_rd = l_v / r_d   !< precomputed l_v / r_d
    9599    REAL(wp), PARAMETER ::  rd_d_rv = r_d / r_v   !< precomputed r_d / r_v
     
    106110    PRIVATE magnus_0d, &
    107111            magnus_1d, &
     112            magnus_tl_0d, &
     113            magnus_tl_1d, &
     114            magnus_0d_ice, &
     115            magnus_1d_ice, &
    108116            ideal_gas_law_rho_0d, &
    109117            ideal_gas_law_rho_1d, &
     
    126134       MODULE PROCEDURE magnus_1d
    127135    END INTERFACE magnus
     136
     137    INTERFACE magnus_tl
     138       MODULE PROCEDURE magnus_tl_0d
     139       MODULE PROCEDURE magnus_tl_1d
     140    END INTERFACE magnus_tl
     141
     142    INTERFACE magnus_ice
     143       MODULE PROCEDURE magnus_0d_ice
     144       MODULE PROCEDURE magnus_1d_ice
     145    END INTERFACE magnus_ice
    128146
    129147    INTERFACE ideal_gas_law_rho
     
    337355! Description:
    338356! ------------
     357!> This function computes the magnus formula (Press et al., 1992) using the
     358!> (ice-) liquid water potential temperature.
     359!> The magnus formula is needed to calculate the saturation vapor pressure over
     360!> a plane liquid water surface
     361!------------------------------------------------------------------------------!
     362    FUNCTION magnus_tl_0d( t_l )
     363
     364       IMPLICIT NONE
     365
     366       REAL(wp), INTENT(IN) ::  t_l  !< liquid water temperature (K)
     367
     368       REAL(wp) ::  magnus_tl_0d
     369!
     370!--    Saturation vapor pressure for a specific temperature:
     371       magnus_tl_0d =  610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /             &
     372                                                    ( t_l - 35.86_wp  ) )
     373
     374    END FUNCTION magnus_tl_0d
     375
     376!------------------------------------------------------------------------------!
     377! Description:
     378! ------------
     379!> This function computes the magnus formula (Press et al., 1992) using the
     380!> (ice-) liquid water potential temperature.
     381!> The magnus formula is needed to calculate the saturation vapor pressure over
     382!> a plane liquid water surface
     383!------------------------------------------------------------------------------!
     384    FUNCTION magnus_tl_1d( t_l )
     385
     386       IMPLICIT NONE
     387
     388       REAL(wp), INTENT(IN), DIMENSION(:) ::  t_l  !< liquid water temperature (K)
     389
     390       REAL(wp), DIMENSION(size(t_l)) ::  magnus_tl_1d
     391!
     392!--    Saturation vapor pressure for a specific temperature:
     393       magnus_tl_1d =  610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /      &
     394                                                    ( t_l - 35.86_wp  ) )
     395
     396    END FUNCTION magnus_tl_1d
     397
     398!------------------------------------------------------------------------------!
     399! Description:
     400! ------------
     401!> This function computes the magnus formula (Press et al., 1992).
     402!> The magnus formula is needed to calculate the saturation vapor pressure over
     403!> a plane ice surface
     404!------------------------------------------------------------------------------!
     405    FUNCTION magnus_0d_ice( t )
     406
     407       IMPLICIT NONE
     408
     409       REAL(wp), INTENT(IN) ::  t  !< temperature (K)
     410
     411       REAL(wp) ::  magnus_0d_ice
     412!
     413!--    Saturation vapor pressure for a specific temperature:
     414       magnus_0d_ice =  611.2_wp * EXP( 22.46_wp * ( t - degc_to_k ) /         &
     415                                                   ( t - 0.53_wp  ) )
     416
     417    END FUNCTION magnus_0d_ice
     418
     419!------------------------------------------------------------------------------!
     420! Description:
     421! ------------
     422!> This function computes the magnus formula (Press et al., 1992).
     423!> The magnus formula is needed to calculate the saturation vapor pressure over
     424!> a plane ice surface
     425!------------------------------------------------------------------------------!
     426    FUNCTION magnus_1d_ice( t )
     427
     428       IMPLICIT NONE
     429
     430       REAL(wp), INTENT(IN), DIMENSION(:) ::  t  !< temperature (K)
     431
     432       REAL(wp), DIMENSION(size(t)) ::  magnus_1d_ice
     433!
     434!--    Saturation vapor pressure for a specific temperature:
     435       magnus_1d_ice =  611.2_wp * EXP( 22.46_wp * ( t - degc_to_k ) /          &
     436                                                   ( t - 0.53_wp  ) )
     437
     438    END FUNCTION magnus_1d_ice
     439
     440!------------------------------------------------------------------------------!
     441! Description:
     442! ------------
    339443!> Compute the ideal gas law for scalar arguments.
    340444!------------------------------------------------------------------------------!
  • palm/trunk/SOURCE/bulk_cloud_model_mod.f90

    r4495 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4495 2020-04-13 20:11:20Z raasch
    2730! restart data handling with MPI-IO added
    2831!
     
    116119               precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,   &
    117120               qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,         &
    118                exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc, &
    119                flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr, &
    120                flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr, &
    121                diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr, &
    122                diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr
     121               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc,            &
     122               flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr,                     &
     123               flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr,                     &
     124               diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr,                     &
     125               diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr,                     &
     126               ni, ni_1, ni_2, ni_3, ni_p, tni_m,                              &
     127               qi, qi_1, qi_2, qi_3, qi_p, tqi_m,                              &
     128               flux_l_qi, flux_l_ni, flux_s_qi, flux_s_ni,                     &
     129               diss_l_qi, diss_l_ni, diss_s_qi, diss_s_ni
     130
    123131
    124132    USE averaging,                                                             &
    125         ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av
     133        ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av, ni_av, qi_av
    126134
    127135    USE basic_constants_and_equations_mod,                                     &
    128         ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,&
     136        ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, magnus_ice,              &
     137               molecular_weight_of_solute,                                     &
    129138               molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,&
    130139               exner_function, exner_function_invers, ideal_gas_law_rho,       &
    131                ideal_gas_law_rho_pt, barometric_formula, rd_d_rv
     140               ideal_gas_law_rho_pt, barometric_formula, rd_d_rv, l_s,         &
     141               ls_d_cp
    132142
    133143    USE control_parameters,                                                    &
     
    143153               dt_3d, dt_do2d_xy, intermediate_timestep_count,                 &
    144154               intermediate_timestep_count_max, large_scale_forcing,           &
    145                lsf_surf, pt_surface, restart_data_format_output, rho_surface, surface_pressure,   &
     155               lsf_surf, pt_surface, restart_data_format_output, rho_surface,  &
     156               surface_pressure,                                               &
    146157               time_do2d_xy, message_string, initializing_actions,             &
    147                ws_scheme_sca, scalar_advec, timestep_scheme, tsc, loop_optimization
     158               ws_scheme_sca, scalar_advec, timestep_scheme, tsc,              &
     159               loop_optimization, simulated_time
    148160
    149161    USE cpulog,                                                                &
     
    167179        ONLY:  threads_per_task
    168180
    169     USE restart_data_mpi_io_mod,                                                                   &
     181    USE restart_data_mpi_io_mod,                                               &
    170182        ONLY:  rrd_mpi_io, wrd_mpi_io
    171183
    172184    USE statistics,                                                            &
    173185        ONLY:  weight_pres, weight_substep, sums_wsncs_ws_l, sums_wsnrs_ws_l,  &
    174                sums_wsqcs_ws_l, sums_wsqrs_ws_l
     186               sums_wsqcs_ws_l, sums_wsqrs_ws_l,                               &
     187               sums_wsqis_ws_l, sums_wsnis_ws_l
    175188
    176189    USE surface_mod,                                                           &
    177190        ONLY :  bc_h,                                                          &
    178191                surf_bulk_cloud_model,                                         &
    179                 surf_microphysics_morrison, surf_microphysics_seifert, &
    180                 surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
     192                surf_microphysics_morrison, surf_microphysics_seifert,         &
     193                surf_microphysics_ice_extension,                               &
     194                surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
     195                surf_usm_v
    181196
    182197    IMPLICIT NONE
     
    193208    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
    194209    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
     210    LOGICAL ::  ice_crystal_sedimentation = .FALSE.        !< flag for ice crystal sedimentation
    195211    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
    196212    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
     
    198214
    199215    LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.      !< namelist parameter
     216    LOGICAL ::  microphysics_ice_extension = .FALSE.              !< use ice microphysics scheme
    200217    LOGICAL ::  microphysics_sat_adjust = .FALSE.                !< use saturation adjust bulk scheme?
    201218    LOGICAL ::  microphysics_kessler = .FALSE.                   !< use kessler bulk scheme?
    202219    LOGICAL ::  microphysics_morrison = .FALSE.                  !< use 2-moment Morrison (add. prog. eq. for nc and qc)
    203220    LOGICAL ::  microphysics_seifert = .FALSE.                   !< use 2-moment Seifert and Beheng scheme
    204     LOGICAL ::  microphysics_morrison_no_rain = .FALSE.          !< use 2-moment Morrison     
     221    LOGICAL ::  microphysics_morrison_no_rain = .FALSE.          !< use 2-moment Morrison
    205222    LOGICAL ::  precipitation = .FALSE.                          !< namelist parameter
    206223
     
    240257    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
    241258    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
    242 !    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
     259    REAL(wp) ::  ximin = 4.42E-14_wp       !< minimum mass of ice crystal (kg) (D~10µm)
    243260    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
    244261    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
     
    249266    REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
    250267    REAL(wp) ::  dt_micro                        !< microphysics time step
     268    REAL(wp) ::  in_init = 1000.0_wp             !< initial number of potential ice nucleii
    251269    REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
    252270    REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
     
    254272    REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
    255273    REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
    256     REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0;
     274    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0
     275    REAL(wp) ::  start_ice_microphysics = 0.0_wp !< time from which on ice microhysics are executed
    257276
    258277    REAL(wp) ::  e_s     !< saturation water vapor pressure
     278    REAL(wp) ::  e_si    !< saturation water vapor pressure over ice
    259279    REAL(wp) ::  q_s     !< saturation mixing ratio
     280    REAL(wp) ::  q_si    !< saturation mixing ratio over ice
    260281    REAL(wp) ::  sat     !< supersaturation
    261     REAL(wp) ::  t_l     !< actual temperature
     282    REAL(wp) ::  sat_ice !< supersaturation over ice
     283    REAL(wp) ::  t_l     !< liquid-(ice) water temperature
    262284
    263285    SAVE
     
    294316           dt_precipitation, &
    295317           microphysics_morrison, &
    296            microphysics_morrison_no_rain, &           
     318           microphysics_morrison_no_rain, &
    297319           microphysics_sat_adjust, &
    298320           microphysics_seifert, &
     321           microphysics_ice_extension, &
    299322           na_init, &
    300323           nc_const, &
    301324           precipitation, &
    302            sigma_gc
    303 
     325           sigma_gc, &
     326           start_ice_microphysics, &
     327           ice_crystal_sedimentation, &
     328           in_init
    304329
    305330    INTERFACE bcm_parin
     
    420445          nc_const, &
    421446          sigma_bulk, &
    422           ventilation_effect
     447          ventilation_effect, &
     448          ice_crystal_sedimentation, &
     449          microphysics_ice_extension, &
     450          start_ice_microphysics, &
     451          in_init
    423452
    424453       line = ' '
     
    517546          microphysics_kessler    = .FALSE.
    518547          microphysics_morrison   = .TRUE.
    519           microphysics_morrison_no_rain = .TRUE.         
    520           precipitation           = .FALSE.     
     548          microphysics_morrison_no_rain = .TRUE.
     549          precipitation           = .FALSE.
    521550       ELSE
    522551          message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
     
    546575       surf_microphysics_morrison = microphysics_morrison
    547576       surf_microphysics_seifert = microphysics_seifert
    548 
     577       surf_microphysics_ice_extension = microphysics_ice_extension
    549578!
    550579!--    Check aerosol
     
    592621             unit = '1/m3'
    593622
     623          CASE ( 'ni' )
     624             IF ( .NOT.  microphysics_ice_extension )  THEN
     625                message_string = 'output of "' // TRIM( var ) // '" ' //       &
     626                                 'requires ' //                                &
     627                                 'microphysics_ice_extension = ".TRUE."'
     628                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     629             ENDIF
     630             unit = '1/m3'
     631
    594632          CASE ( 'nr' )
    595633             IF ( .NOT.  microphysics_seifert )  THEN
     
    611649
    612650          CASE ( 'qc' )
     651             unit = 'kg/kg'
     652
     653          CASE ( 'qi' )
     654             IF ( .NOT.  microphysics_ice_extension ) THEN
     655                message_string = 'output of "' // TRIM( var ) // '" ' //       &
     656                                 'requires ' //                                &
     657                                 'microphysics_ice_extension = ".TRUE."'
     658                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     659             ENDIF
    613660             unit = 'kg/kg'
    614661
     
    698745             ENDIF
    699746             pr_index = 123
     747             dopr_index(var_count) = pr_index
     748             dopr_unit     = '1/m3'
     749             unit = dopr_unit
     750             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
     751
     752          CASE ( 'ni' )
     753             IF ( .NOT.  microphysics_ice_extension )  THEN
     754                message_string = 'data_output_pr = ' //                        &
     755                                 TRIM( data_output_pr(var_count) ) //          &
     756                                 ' is not implemented for' //                  &
     757                                 ' microphysics_ice_extension = ".F."'
     758                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
     759             ENDIF
     760             pr_index = 124
    700761             dopr_index(var_count) = pr_index
    701762             dopr_unit     = '1/m3'
     
    730791             unit = dopr_unit
    731792             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
     793
    732794          CASE ( 'qc' )
    733795             pr_index = 75
     796             dopr_index(var_count) = pr_index
     797             dopr_unit     = 'kg/kg'
     798             unit = dopr_unit
     799             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
     800
     801          CASE ( 'qi' )
     802             IF ( .NOT.  microphysics_ice_extension )  THEN
     803                message_string = 'data_output_pr = ' //                        &
     804                                 TRIM( data_output_pr(var_count) ) //          &
     805                                 ' is not implemented for' //                  &
     806                                 ' microphysics_ice_extension = ".F."'
     807                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
     808             ENDIF
     809             pr_index = 125
    734810             dopr_index(var_count) = pr_index
    735811             dopr_unit     = 'kg/kg'
     
    792868!
    793869!--       3D-cloud drop water content, cloud drop concentration arrays
    794           ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    795                     nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    796                     nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    797                     qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    798                     qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
     870          ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     871                    nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     872                    nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     873                    qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     874                    qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    799875                    qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    800876       ENDIF
     
    803879!
    804880!--       3D-rain water content, rain drop concentration arrays
    805           ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    806                     nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    807                     nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    808                     qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
    809                     qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
     881          ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     882                    nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     883                    nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     884                    qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     885                    qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
    810886                    qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    811887       ENDIF
    812888
     889       IF ( microphysics_ice_extension )  THEN
     890!
     891!--       3D-cloud drop water content, cloud drop concentration arrays
     892          ALLOCATE( ni_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     893                    ni_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     894                    ni_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     895                    qi_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     896                    qi_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                       &
     897                    qi_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     898       ENDIF
     899
    813900       IF ( ws_scheme_sca )  THEN
    814 
    815901          IF ( microphysics_morrison )  THEN
    816902             ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     
    819905             sums_wsncs_ws_l = 0.0_wp
    820906          ENDIF
    821 
    822907          IF ( microphysics_seifert )  THEN
    823908             ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     
    826911             sums_wsnrs_ws_l = 0.0_wp
    827912          ENDIF
    828 
     913          IF ( microphysics_ice_extension )  THEN
     914             ALLOCATE( sums_wsqis_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     915             ALLOCATE( sums_wsnis_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     916             sums_wsqis_ws_l = 0.0_wp
     917             sums_wsnis_ws_l = 0.0_wp
     918          ENDIF
    829919       ENDIF
    830920
     
    835925!--    advection routines.
    836926       IF ( loop_optimization /= 'vector' )  THEN
    837 
    838927          IF ( ws_scheme_sca )  THEN
    839 
    840928             IF ( microphysics_morrison )  THEN
    841929                ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
     
    848936                          diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    849937             ENDIF
    850 
    851938             IF ( microphysics_seifert )  THEN
    852939                ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),           &
     
    859946                          diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    860947             ENDIF
    861 
    862           ENDIF
    863 
     948             IF ( microphysics_ice_extension )  THEN
     949                ALLOCATE( flux_s_qi(nzb+1:nzt,0:threads_per_task-1),           &
     950                          diss_s_qi(nzb+1:nzt,0:threads_per_task-1),           &
     951                          flux_s_ni(nzb+1:nzt,0:threads_per_task-1),           &
     952                          diss_s_ni(nzb+1:nzt,0:threads_per_task-1) )
     953                ALLOCATE( flux_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     954                          diss_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     955                          flux_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     956                          diss_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
     957             ENDIF
     958          ENDIF
    864959       ENDIF
    865960
     
    878973          nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
    879974       ENDIF
     975       IF ( microphysics_ice_extension )  THEN
     976          qi => qi_1;  qi_p  => qi_2;  tqi_m  => qi_3
     977          ni => ni_1;  ni_p  => ni_2;  tni_m  => ni_3
     978       ENDIF
     979
    880980
    881981
     
    9191019             ENDIF
    9201020!
     1021!--          Initialize the remaining quantities
     1022             IF ( microphysics_ice_extension )  THEN
     1023                DO  i = nxlg, nxrg
     1024                   DO  j = nysg, nyng
     1025                      qi(:,j,i) = 0.0_wp
     1026                      ni(:,j,i) = 0.0_wp
     1027                   ENDDO
     1028                ENDDO
     1029             ENDIF
     1030!
    9211031!--          Liquid water content and precipitation amount
    9221032!--          are zero at beginning of the simulation
     
    9381048                qr_p  = qr
    9391049                nr_p  = nr
     1050             ENDIF
     1051             IF ( microphysics_ice_extension )  THEN
     1052                tqi_m = 0.0_wp
     1053                tni_m = 0.0_wp
     1054                qi_p  = qi
     1055                ni_p  = ni
    9401056             ENDIF
    9411057          ENDIF ! Only if not read_restart_data
     
    10801196                sums_wsnrs_ws_l = 0.0_wp
    10811197             ENDIF
     1198             IF ( microphysics_ice_extension )  THEN
     1199                sums_wsqis_ws_l = 0.0_wp
     1200                sums_wsnis_ws_l = 0.0_wp
     1201             ENDIF
    10821202
    10831203          ENDIF
     
    10961216!> Control of microphysics for grid points i,j
    10971217!------------------------------------------------------------------------------!
    1098 
    10991218    SUBROUTINE bcm_actions_ij( i, j, location )
    11001219
     
    11211240                sums_wsnrs_ws_l = 0.0_wp
    11221241             ENDIF
     1242             IF ( microphysics_ice_extension )  THEN
     1243                sums_wsqis_ws_l = 0.0_wp
     1244                sums_wsnis_ws_l = 0.0_wp
     1245             ENDIF
     1246
    11231247
    11241248          ENDIF
     
    11431267       CALL cpu_log( log_point(51), 'microphysics', 'start' )
    11441268
    1145        IF ( .NOT. microphysics_sat_adjust  .AND.         &
    1146             ( intermediate_timestep_count == 1  .OR.                              &
    1147               call_microphysics_at_all_substeps ) )                               &
     1269       IF ( .NOT. microphysics_sat_adjust  .AND.                               &
     1270            ( intermediate_timestep_count == 1  .OR.                           &
     1271              call_microphysics_at_all_substeps ) )                            &
    11481272       THEN
    11491273
     
    11511275!
    11521276!--          Calculate vertical profile of the hydrostatic pressure (hyp)
    1153              hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
     1277             hyp    = barometric_formula(zu, pt_surface *                      &
     1278                      exner_function(surface_pressure * 100.0_wp),             &
     1279                      surface_pressure * 100.0_wp)
    11541280             d_exner = exner_function_invers(hyp)
    11551281             exner = 1.0_wp / exner_function_invers(hyp)
    1156              hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
     1282             hyrho  = ideal_gas_law_rho_pt(hyp, pt(:, nys, nxl) )
    11571283!
    11581284!--          Compute reference density
    1159              rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
     1285             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,      &
     1286                           pt_surface *                                        &
     1287                           exner_function(surface_pressure * 100.0_wp))
    11601288          ENDIF
    11611289
     
    11781306             CALL autoconversion_kessler
    11791307             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
    1180              
     1308
    11811309!
    11821310!--       Here the seifert beheng scheme is used. Cloud concentration is assumed to
     
    11841312          ELSEIF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )  THEN
    11851313             CALL adjust_cloud
     1314             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1315             CALL ice_nucleation
     1316             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1317             CALL ice_deposition
    11861318             CALL autoconversion
    11871319             CALL accretion
     
    11901322             CALL sedimentation_rain
    11911323             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
    1192              
    1193 !
    1194 !--       Here the morrison scheme is used. No rain processes are considered and qr and nr
     1324             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1325             CALL adjust_ice
     1326             IF ( ice_crystal_sedimentation .AND. microphysics_ice_extension &
     1327                  .AND. simulated_time > start_ice_microphysics ) CALL sedimentation_ice
     1328
     1329!
     1330!--       Here the morrison scheme is used. No rain processes are considered and qr and nr
    11951331!--       are not allocated
    11961332          ELSEIF ( microphysics_morrison_no_rain  .AND.  .NOT. microphysics_seifert )  THEN
    11971333             CALL activation
    11981334             CALL condensation
    1199              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud   
    1200              
    1201 !
    1202 !--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
     1335             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1336             CALL adjust_ice
     1337             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1338             CALL ice_nucleation
     1339             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1340             CALL ice_deposition
     1341             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
     1342
     1343!
     1344!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
    12031345!--       included
    12041346          ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
     
    12061348             CALL activation
    12071349             CALL condensation
     1350             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1351             CALL adjust_ice
     1352             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1353             CALL ice_nucleation
     1354             IF ( microphysics_ice_extension  .AND.  simulated_time > start_ice_microphysics )     &
     1355             CALL ice_deposition
    12081356             CALL autoconversion
    12091357             CALL accretion
     
    12111359             CALL evaporation_rain
    12121360             CALL sedimentation_rain
    1213              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud           
     1361             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
    12141362
    12151363          ENDIF
     
    12291377!> Control of microphysics for grid points i,j
    12301378!------------------------------------------------------------------------------!
    1231 
    12321379    SUBROUTINE bcm_non_advective_processes_ij( i, j )
    12331380
     
    12361383       INTEGER(iwp) ::  j                 !<
    12371384
    1238        IF ( .NOT. microphysics_sat_adjust  .AND.         &
    1239             ( intermediate_timestep_count == 1  .OR.                              &
    1240               call_microphysics_at_all_substeps ) )                               &
     1385       IF ( .NOT. microphysics_sat_adjust  .AND.                               &
     1386            ( intermediate_timestep_count == 1  .OR.                           &
     1387              call_microphysics_at_all_substeps ) )                            &
    12411388       THEN
    12421389
     
    12441391!
    12451392!--          Calculate vertical profile of the hydrostatic pressure (hyp)
    1246              hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
     1393             hyp    = barometric_formula(zu, pt_surface *                      &
     1394                      exner_function(surface_pressure * 100.0_wp),             &
     1395                      surface_pressure * 100.0_wp)
    12471396             d_exner = exner_function_invers(hyp)
    12481397             exner = 1.0_wp / exner_function_invers(hyp)
    1249              hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
     1398             hyrho  = ideal_gas_law_rho_pt(hyp, pt(:, nys, nxl) )
    12501399!
    12511400!--          Compute reference density
    1252              rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
     1401             rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp,      &
     1402                           pt_surface *                                        &
     1403                           exner_function(surface_pressure * 100.0_wp))
    12531404          ENDIF
    12541405
     
    12731424!
    12741425!--       Here the seifert beheng scheme is used. Cloud concentration is assumed to
    1275 !--       a constant value an qc a diagnostic value.             
     1426!--       a constant value an qc a diagnostic value.
    12761427          ELSEIF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )  THEN
    12771428             CALL adjust_cloud_ij( i,j )
     1429             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1430                 CALL ice_nucleation_ij( i,j )
     1431             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1432                 CALL ice_deposition_ij( i,j )
    12781433             CALL autoconversion_ij( i,j )
    12791434             CALL accretion_ij( i,j )
     
    12821437             CALL sedimentation_rain_ij( i,j )
    12831438             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
    1284              
    1285 !
    1286 !--       Here the morrison scheme is used. No rain processes are considered and qr and nr
     1439             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1440                 CALL adjust_ice_ij ( i,j )
     1441             IF ( ice_crystal_sedimentation .AND. microphysics_ice_extension &
     1442                  .AND. simulated_time > start_ice_microphysics )  CALL sedimentation_ice_ij ( i,j )
     1443!
     1444!--       Here the morrison scheme is used. No rain processes are considered and qr and nr
    12871445!--       are not allocated
    12881446          ELSEIF ( microphysics_morrison_no_rain  .AND.  .NOT. microphysics_seifert )  THEN
    12891447             CALL activation_ij( i,j )
    12901448             CALL condensation_ij( i,j )
    1291              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
    1292              
    1293 !
    1294 !--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
    1295 !--       included             
     1449             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1450                 CALL adjust_ice_ij ( i,j )
     1451             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1452                 CALL ice_nucleation_ij( i,j )
     1453             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1454                 CALL ice_deposition_ij( i,j )
     1455             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
     1456
     1457!
     1458!--       Here the full morrison scheme is used and all processes of Seifert and Beheng are
     1459!--       included
    12961460          ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
    12971461             CALL adjust_cloud_ij( i,j )
    12981462             CALL activation_ij( i,j )
    12991463             CALL condensation_ij( i,j )
     1464             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1465                 CALL adjust_ice_ij ( i,j )
     1466             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1467                 CALL ice_nucleation_ij( i,j )
     1468             IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics )       &
     1469                 CALL ice_deposition_ij( i,j )
    13001470             CALL autoconversion_ij( i,j )
    13011471             CALL accretion_ij( i,j )
     
    13031473             CALL evaporation_rain_ij( i,j )
    13041474             CALL sedimentation_rain_ij( i,j )
    1305              IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )           
     1475             IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
    13061476
    13071477          ENDIF
     
    13311501          IF ( microphysics_morrison )  THEN
    13321502             CALL exchange_horiz( nc, nbgp )
    1333              CALL exchange_horiz( qc, nbgp )         
     1503             CALL exchange_horiz( qc, nbgp )
    13341504          ENDIF
    13351505          IF ( microphysics_seifert ) THEN
     
    13371507             CALL exchange_horiz( nr, nbgp )
    13381508          ENDIF
     1509          IF ( microphysics_ice_extension ) THEN
     1510             CALL exchange_horiz( qi, nbgp )
     1511             CALL exchange_horiz( ni, nbgp )
     1512          ENDIF
    13391513          CALL exchange_horiz( q, nbgp )
    1340           CALL exchange_horiz( pt, nbgp )         
     1514          CALL exchange_horiz( pt, nbgp )
    13411515       ENDIF
    13421516
    13431517
    13441518    END SUBROUTINE bcm_exchange_horiz
    1345    
     1519
    13461520
    13471521
     
    14251599                                             )                                 &
    14261600                                    * MERGE( 1.0_wp, 0.0_wp,                   &
    1427                                        BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1601                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
    14281602                                          )
    14291603                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
     
    15171691                                             )                                 &
    15181692                                   * MERGE( 1.0_wp, 0.0_wp,                    &
    1519                                        BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1693                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
    15201694                                          )
    15211695                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
     
    15511725
    15521726       ENDIF
     1727
     1728!
     1729!--    If required, calculate prognostic equations for ice crystal content
     1730!--    and ice crystal concentration
     1731       IF ( microphysics_ice_extension )  THEN
     1732
     1733          CALL cpu_log( log_point(70), 'qi-equation', 'start' )
     1734
     1735!
     1736!--       Calculate prognostic equation for ice crystal content
     1737          sbt = tsc(2)
     1738          IF ( scalar_advec == 'bc-scheme' )  THEN
     1739
     1740             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1741!
     1742!--             Bott-Chlond scheme always uses Euler time step. Thus:
     1743                sbt = 1.0_wp
     1744             ENDIF
     1745             tend = 0.0_wp
     1746             CALL advec_s_bc( qi, 'qi' )
     1747
     1748          ENDIF
     1749
     1750!
     1751!--       qi-tendency terms with no communication
     1752          IF ( scalar_advec /= 'bc-scheme' )  THEN
     1753             tend = 0.0_wp
     1754             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1755                IF ( ws_scheme_sca )  THEN
     1756                   CALL advec_s_ws( advc_flags_s, qi, 'qi',                    &
     1757                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
     1758                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
     1759                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
     1760                                    bc_dirichlet_s  .OR.  bc_radiation_s )
     1761                ELSE
     1762                   CALL advec_s_pw( qi )
     1763                ENDIF
     1764             ELSE
     1765                CALL advec_s_up( qi )
     1766             ENDIF
     1767          ENDIF
     1768
     1769          CALL diffusion_s( qi,                                                &
     1770                            surf_def_h(0)%qisws, surf_def_h(1)%qisws,          &
     1771                            surf_def_h(2)%qisws,                               &
     1772                            surf_lsm_h%qisws,    surf_usm_h%qisws,             &
     1773                            surf_def_v(0)%qisws, surf_def_v(1)%qisws,          &
     1774                            surf_def_v(2)%qisws, surf_def_v(3)%qisws,          &
     1775                            surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,          &
     1776                            surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,          &
     1777                            surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,          &
     1778                            surf_usm_v(2)%qisws, surf_usm_v(3)%qisws )
     1779
     1780!
     1781!--       Prognostic equation for ice crystal mixing ratio
     1782          DO  i = nxl, nxr
     1783             DO  j = nys, nyn
     1784                !following directive is required to vectorize on Intel19
     1785                !DIR$ IVDEP
     1786                DO  k = nzb+1, nzt
     1787                   qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
     1788                                                      tsc(3) * tqi_m(k,j,i) )  &
     1789                                                    - tsc(5) * rdf_sc(k) *     &
     1790                                                               qi(k,j,i)       &
     1791                                             )                                 &
     1792                                    * MERGE( 1.0_wp, 0.0_wp,                   &
     1793                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1794                                          )
     1795                   IF ( qi_p(k,j,i) < 0.0_wp )  qi_p(k,j,i) = 0.0_wp
     1796                ENDDO
     1797             ENDDO
     1798          ENDDO
     1799
     1800!
     1801!--       Calculate tendencies for the next Runge-Kutta step
     1802          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1803             IF ( intermediate_timestep_count == 1 )  THEN
     1804                DO  i = nxl, nxr
     1805                   DO  j = nys, nyn
     1806                      DO  k = nzb+1, nzt
     1807                         tqi_m(k,j,i) = tend(k,j,i)
     1808                      ENDDO
     1809                   ENDDO
     1810                ENDDO
     1811             ELSEIF ( intermediate_timestep_count < &
     1812                      intermediate_timestep_count_max )  THEN
     1813                DO  i = nxl, nxr
     1814                   DO  j = nys, nyn
     1815                      DO  k = nzb+1, nzt
     1816                         tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
     1817                                         + 5.3125_wp * tqi_m(k,j,i)
     1818                      ENDDO
     1819                   ENDDO
     1820                ENDDO
     1821             ENDIF
     1822          ENDIF
     1823
     1824          CALL cpu_log( log_point(70), 'qi-equation', 'stop' )
     1825
     1826          CALL cpu_log( log_point(69), 'ni-equation', 'start' )
     1827!
     1828!--       Calculate prognostic equation for ice crystal concentration
     1829          sbt = tsc(2)
     1830          IF ( scalar_advec == 'bc-scheme' )  THEN
     1831
     1832             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1833!
     1834!--             Bott-Chlond scheme always uses Euler time step. Thus:
     1835                sbt = 1.0_wp
     1836             ENDIF
     1837             tend = 0.0_wp
     1838             CALL advec_s_bc( ni, 'ni' )
     1839
     1840          ENDIF
     1841
     1842!
     1843!--       ni-tendency terms with no communication
     1844          IF ( scalar_advec /= 'bc-scheme' )  THEN
     1845             tend = 0.0_wp
     1846             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1847                IF ( ws_scheme_sca )  THEN
     1848                   CALL advec_s_ws( advc_flags_s, ni, 'ni',                    &
     1849                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
     1850                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
     1851                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
     1852                                    bc_dirichlet_s  .OR.  bc_radiation_s )
     1853                ELSE
     1854                   CALL advec_s_pw( ni )
     1855                ENDIF
     1856             ELSE
     1857                CALL advec_s_up( ni )
     1858             ENDIF
     1859          ENDIF
     1860
     1861          CALL diffusion_s( ni,                                                &
     1862                            surf_def_h(0)%nisws, surf_def_h(1)%nisws,          &
     1863                            surf_def_h(2)%nisws,                               &
     1864                            surf_lsm_h%nisws,    surf_usm_h%nisws,             &
     1865                            surf_def_v(0)%nisws, surf_def_v(1)%nisws,          &
     1866                            surf_def_v(2)%nisws, surf_def_v(3)%nisws,          &
     1867                            surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,          &
     1868                            surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,          &
     1869                            surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,          &
     1870                            surf_usm_v(2)%nisws, surf_usm_v(3)%nisws )
     1871
     1872!
     1873!--       Prognostic equation for ice crystal concentration
     1874          DO  i = nxl, nxr
     1875             DO  j = nys, nyn
     1876                !following directive is required to vectorize on Intel19
     1877                !DIR$ IVDEP
     1878                DO  k = nzb+1, nzt
     1879                   ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
     1880                                                      tsc(3) * tni_m(k,j,i) )  &
     1881                                                    - tsc(5) * rdf_sc(k) *     &
     1882                                                               ni(k,j,i)       &
     1883                                             )                                 &
     1884                                   * MERGE( 1.0_wp, 0.0_wp,                    &
     1885                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1886                                          )
     1887                   IF ( ni_p(k,j,i) < 0.0_wp )  ni_p(k,j,i) = 0.0_wp
     1888                ENDDO
     1889             ENDDO
     1890          ENDDO
     1891
     1892!
     1893!--       Calculate tendencies for the next Runge-Kutta step
     1894          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1895             IF ( intermediate_timestep_count == 1 )  THEN
     1896                DO  i = nxl, nxr
     1897                   DO  j = nys, nyn
     1898                      DO  k = nzb+1, nzt
     1899                         tni_m(k,j,i) = tend(k,j,i)
     1900                      ENDDO
     1901                   ENDDO
     1902                ENDDO
     1903             ELSEIF ( intermediate_timestep_count < &
     1904                      intermediate_timestep_count_max )  THEN
     1905                DO  i = nxl, nxr
     1906                   DO  j = nys, nyn
     1907                      DO  k = nzb+1, nzt
     1908                         tni_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
     1909                                         + 5.3125_wp * tni_m(k,j,i)
     1910                      ENDDO
     1911                   ENDDO
     1912                ENDDO
     1913             ENDIF
     1914          ENDIF
     1915
     1916          CALL cpu_log( log_point(69), 'ni-equation', 'stop' )
     1917
     1918       ENDIF
     1919
    15531920!
    15541921!--    If required, calculate prognostic equations for rain water content
     
    16161983                                             )                                 &
    16171984                                    * MERGE( 1.0_wp, 0.0_wp,                   &
    1618                                        BTEST( wall_flags_total_0(k,j,i), 0 )   &
     1985                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
    16191986                                          )
    16201987                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
     
    17082075                                             )                                 &
    17092076                                   * MERGE( 1.0_wp, 0.0_wp,                    &
    1710                                        BTEST( wall_flags_total_0(k,j,i), 0 )   &
     2077                                             BTEST( wall_flags_total_0(k,j,i), 0 )   &
    17112078                                          )
    17122079                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
     
    17512118!> Control of microphysics for grid points i,j
    17522119!------------------------------------------------------------------------------!
    1753 
    17542120    SUBROUTINE bcm_prognostic_equations_ij( i, j, i_omp_start, tn )
    17552121
     
    18052171                                       )                                 &
    18062172                                 * MERGE( 1.0_wp, 0.0_wp,                &
    1807                                     BTEST( wall_flags_total_0(k,j,i), 0 )&
     2173                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
    18082174                                        )
    18092175             IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
     
    18642230                                       )                                 &
    18652231                                 * MERGE( 1.0_wp, 0.0_wp,                &
    1866                                     BTEST( wall_flags_total_0(k,j,i), 0 )&
     2232                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
    18672233                                        )
    18682234             IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
     
    18852251
    18862252       ENDIF
     2253
     2254!
     2255!--    If required, calculate prognostic equations for ice crystal mixing ratio
     2256!--    and ice crystal concentration
     2257       IF ( microphysics_ice_extension )  THEN
     2258!
     2259!--       Calculate prognostic equation for ice crystal mixing ratio
     2260          tend(:,j,i) = 0.0_wp
     2261          IF ( timestep_scheme(1:5) == 'runge' ) &
     2262          THEN
     2263             IF ( ws_scheme_sca )  THEN
     2264                CALL advec_s_ws( advc_flags_s, i, j, qi, 'qi', flux_s_qi,      &
     2265                                 diss_s_qi, flux_l_qi, diss_l_qi,              &
     2266                                 i_omp_start, tn,                              &
     2267                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
     2268                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
     2269                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2270                                 bc_dirichlet_s  .OR.  bc_radiation_s )
     2271             ELSE
     2272                CALL advec_s_pw( i, j, qi )
     2273             ENDIF
     2274          ELSE
     2275             CALL advec_s_up( i, j, qi )
     2276          ENDIF
     2277          CALL diffusion_s( i, j, qi,                                   &
     2278                            surf_def_h(0)%qisws, surf_def_h(1)%qisws,   &
     2279                            surf_def_h(2)%qisws,                        &
     2280                            surf_lsm_h%qisws,    surf_usm_h%qisws,      &
     2281                            surf_def_v(0)%qisws, surf_def_v(1)%qisws,   &
     2282                            surf_def_v(2)%qisws, surf_def_v(3)%qisws,   &
     2283                            surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws,   &
     2284                            surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws,   &
     2285                            surf_usm_v(0)%qisws, surf_usm_v(1)%qisws,   &
     2286                            surf_usm_v(2)%qisws, surf_usm_v(3)%qisws )
     2287
     2288!
     2289!--       Prognostic equation for ice crystal mixing ratio
     2290          DO  k = nzb+1, nzt
     2291             qi_p(k,j,i) = qi(k,j,i) + ( dt_3d *                         &
     2292                                                ( tsc(2) * tend(k,j,i) + &
     2293                                                  tsc(3) * tqi_m(k,j,i) )&
     2294                                                - tsc(5) * rdf_sc(k)     &
     2295                                                         * qi(k,j,i)     &
     2296                                       )                                 &
     2297                                 * MERGE( 1.0_wp, 0.0_wp,                &
     2298                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
     2299                                        )
     2300             IF ( qi_p(k,j,i) < 0.0_wp )  qi_p(k,j,i) = 0.0_wp
     2301          ENDDO
     2302!
     2303!--       Calculate tendencies for the next Runge-Kutta step
     2304          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2305             IF ( intermediate_timestep_count == 1 )  THEN
     2306                DO  k = nzb+1, nzt
     2307                   tqi_m(k,j,i) = tend(k,j,i)
     2308                ENDDO
     2309             ELSEIF ( intermediate_timestep_count < &
     2310                      intermediate_timestep_count_max )  THEN
     2311                DO  k = nzb+1, nzt
     2312                   tqi_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
     2313                                     5.3125_wp * tqi_m(k,j,i)
     2314                ENDDO
     2315             ENDIF
     2316          ENDIF
     2317
     2318!
     2319!--       Calculate prognostic equation for ice crystal concentration.
     2320          tend(:,j,i) = 0.0_wp
     2321          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2322             IF ( ws_scheme_sca )  THEN
     2323                CALL advec_s_ws( advc_flags_s, i, j, ni, 'ni', flux_s_ni,      &
     2324                                 diss_s_ni, flux_l_ni, diss_l_ni,              &
     2325                                 i_omp_start, tn,                              &
     2326                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
     2327                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
     2328                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
     2329                                 bc_dirichlet_s  .OR.  bc_radiation_s )
     2330             ELSE
     2331                CALL advec_s_pw( i, j, ni )
     2332             ENDIF
     2333          ELSE
     2334             CALL advec_s_up( i, j, ni )
     2335          ENDIF
     2336          CALL diffusion_s( i, j, ni,                                    &
     2337                            surf_def_h(0)%nisws, surf_def_h(1)%nisws,    &
     2338                            surf_def_h(2)%nisws,                         &
     2339                            surf_lsm_h%nisws,    surf_usm_h%nisws,       &
     2340                            surf_def_v(0)%nisws, surf_def_v(1)%nisws,    &
     2341                            surf_def_v(2)%nisws, surf_def_v(3)%nisws,    &
     2342                            surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws,    &
     2343                            surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws,    &
     2344                            surf_usm_v(0)%nisws, surf_usm_v(1)%nisws,    &
     2345                            surf_usm_v(2)%nisws, surf_usm_v(3)%nisws )
     2346
     2347!
     2348!--       Prognostic equation for ice crystal concentration
     2349          DO  k = nzb+1, nzt
     2350             ni_p(k,j,i) = ni(k,j,i) + ( dt_3d *                         &
     2351                                                ( tsc(2) * tend(k,j,i) + &
     2352                                                  tsc(3) * tni_m(k,j,i) )&
     2353                                                - tsc(5) * rdf_sc(k)     &
     2354                                                         * ni(k,j,i)     &
     2355                                       )                                 &
     2356                                 * MERGE( 1.0_wp, 0.0_wp,                &
     2357                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
     2358                                        )
     2359             IF ( ni_p(k,j,i) < 0.0_wp )  ni_p(k,j,i) = 0.0_wp
     2360          ENDDO
     2361!
     2362!--       Calculate tendencies for the next Runge-Kutta step
     2363          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2364             IF ( intermediate_timestep_count == 1 )  THEN
     2365                DO  k = nzb+1, nzt
     2366                   tni_m(k,j,i) = tend(k,j,i)
     2367                ENDDO
     2368             ELSEIF ( intermediate_timestep_count < &
     2369                      intermediate_timestep_count_max )  THEN
     2370                DO  k = nzb+1, nzt
     2371                   tni_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
     2372                                     5.3125_wp * tni_m(k,j,i)
     2373                ENDDO
     2374             ENDIF
     2375          ENDIF
     2376
     2377       ENDIF
     2378
    18872379!
    18882380!--    If required, calculate prognostic equations for rain water content
     
    19292421                                       )                                 &
    19302422                                 * MERGE( 1.0_wp, 0.0_wp,                &
    1931                                     BTEST( wall_flags_total_0(k,j,i), 0 )&
     2423                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
    19322424                                        )
    19332425             IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
     
    19882480                                       )                                 &
    19892481                                 * MERGE( 1.0_wp, 0.0_wp,                &
    1990                                     BTEST( wall_flags_total_0(k,j,i), 0 )&
     2482                                          BTEST( wall_flags_total_0(k,j,i), 0 )&
    19912483                                        )
    19922484             IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
     
    20382530                   nr => nr_1;    nr_p => nr_2
    20392531                ENDIF
     2532                IF ( microphysics_ice_extension )  THEN
     2533                   qi => qi_1;    qi_p => qi_2
     2534                   ni => ni_1;    ni_p => ni_2
     2535                ENDIF
    20402536
    20412537             CASE ( 1 )
     
    20492545                   nr => nr_2;    nr_p => nr_1
    20502546                ENDIF
     2547                IF ( microphysics_ice_extension )  THEN
     2548                   qi => qi_2;    qi_p => qi_1
     2549                   ni => ni_2;    ni_p => ni_1
     2550                ENDIF
     2551
    20512552
    20522553          END SELECT
     
    20932594       ENDIF
    20942595
     2596       IF ( microphysics_ice_extension )  THEN
     2597!
     2598!--       Surface conditions ice crysral (Dirichlet)
     2599!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
     2600!--       the k coordinate belongs to the atmospheric grid point, therefore, set
     2601!--       qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls
     2602          DO  l = 0, 1
     2603          !$OMP PARALLEL DO PRIVATE( i, j, k )
     2604             DO  m = 1, bc_h(l)%ns
     2605                i = bc_h(l)%i(m)
     2606                j = bc_h(l)%j(m)
     2607                k = bc_h(l)%k(m)
     2608                qi_p(k+bc_h(l)%koff,j,i) = 0.0_wp
     2609                ni_p(k+bc_h(l)%koff,j,i) = 0.0_wp
     2610             ENDDO
     2611          ENDDO
     2612!
     2613!--       Top boundary condition for ice crystal (Dirichlet)
     2614          qi_p(nzt+1,:,:) = 0.0_wp
     2615          ni_p(nzt+1,:,:) = 0.0_wp
     2616
     2617       ENDIF
     2618
     2619
    20952620       IF ( microphysics_seifert )  THEN
    20962621!
     
    21292654             nr_p(:,nys-1,:) = nr_p(:,nys,:)
    21302655          ENDIF
     2656          IF ( microphysics_ice_extension )  THEN
     2657             qi_p(:,nys-1,:) = qi_p(:,nys,:)
     2658             ni_p(:,nys-1,:) = ni_p(:,nys,:)
     2659          ENDIF
    21312660       ELSEIF ( bc_radiation_n )  THEN
    21322661          IF ( microphysics_morrison )  THEN
     
    21382667             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
    21392668          ENDIF
     2669          IF ( microphysics_ice_extension )  THEN
     2670             qi_p(:,nyn+1,:) = qi_p(:,nyn,:)
     2671             ni_p(:,nyn+1,:) = ni_p(:,nyn,:)
     2672          ENDIF
    21402673       ELSEIF ( bc_radiation_l )  THEN
    21412674          IF ( microphysics_morrison )  THEN
     
    21472680             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
    21482681          ENDIF
     2682          IF ( microphysics_ice_extension )  THEN
     2683             qi_p(:,:,nxl-1) = qi_p(:,:,nxl)
     2684             ni_p(:,:,nxl-1) = ni_p(:,:,nxl)
     2685          ENDIF
    21492686       ELSEIF ( bc_radiation_r )  THEN
    21502687          IF ( microphysics_morrison )  THEN
     
    21552692             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
    21562693             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
     2694          ENDIF
     2695          IF ( microphysics_ice_extension )  THEN
     2696             qi_p(:,:,nxr+1) = qi_p(:,:,nxr)
     2697             ni_p(:,:,nxr+1) = ni_p(:,:,nxr)
    21572698          ENDIF
    21582699       ENDIF
     
    24332974             ENDIF
    24342975             to_be_resorted => nc_av
     2976          ENDIF
     2977          IF ( mode == 'xy' ) grid = 'zu'
     2978
     2979       CASE ( 'ni_xy', 'ni_xz', 'ni_yz' )
     2980          IF ( av == 0 )  THEN
     2981             to_be_resorted => ni
     2982          ELSE
     2983             IF ( .NOT. ALLOCATED( ni_av ) ) THEN
     2984                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2985                ni_av = REAL( fill_value, KIND = wp )
     2986             ENDIF
     2987             to_be_resorted => ni_av
    24352988          ENDIF
    24362989          IF ( mode == 'xy' ) grid = 'zu'
     
    24993052          IF ( mode == 'xy' ) grid = 'zu'
    25003053
     3054       CASE ( 'qi_xy', 'qi_xz', 'qi_yz' )
     3055          IF ( av == 0 )  THEN
     3056             to_be_resorted => qi
     3057          ELSE
     3058             IF ( .NOT. ALLOCATED( qi_av ) ) THEN
     3059                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3060                qi_av = REAL( fill_value, KIND = wp )
     3061             ENDIF
     3062             to_be_resorted => qi_av
     3063          ENDIF
     3064          IF ( mode == 'xy' ) grid = 'zu'
     3065
    25013066       CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
    25023067          IF ( av == 0 )  THEN
     
    25343099             DO  k = nzb_do, nzt_do
    25353100                local_pf(i,j,k) = MERGE( &
    2536                                      to_be_resorted(k,j,i),                      &
    2537                                      REAL( fill_value, KIND = wp ),              &
    2538                                      BTEST( wall_flags_total_0(k,j,i), flag_nr ) &
     3101                                     to_be_resorted(k,j,i),                    &
     3102                                     REAL( fill_value, KIND = wp ),            &
     3103                                     BTEST( wall_flags_total_0(k,j,i), flag_nr )     &
    25393104                                  )
    25403105             ENDDO
     
    25953160             ENDIF
    25963161             to_be_resorted => nc_av
     3162          ENDIF
     3163
     3164       CASE ( 'ni' )
     3165          IF ( av == 0 )  THEN
     3166             to_be_resorted => ni
     3167          ELSE
     3168             IF ( .NOT. ALLOCATED( ni_av ) ) THEN
     3169                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3170                ni_av = REAL( fill_value, KIND = wp )
     3171             ENDIF
     3172             to_be_resorted => ni_av
    25973173          ENDIF
    25983174
     
    26433219          ENDIF
    26443220
     3221       CASE ( 'qi' )
     3222          IF ( av == 0 )  THEN
     3223             to_be_resorted => qi
     3224          ELSE
     3225             IF ( .NOT. ALLOCATED( qi_av ) ) THEN
     3226                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3227                qi_av = REAL( fill_value, KIND = wp )
     3228             ENDIF
     3229             to_be_resorted => qi_av
     3230          ENDIF
     3231
    26453232       CASE ( 'ql' )
    26463233          IF ( av == 0 )  THEN
     
    26753262          DO  j = nys, nyn
    26763263             DO  k = nzb_do, nzt_do
    2677                 local_pf(i,j,k) = MERGE(                                         &
    2678                                      to_be_resorted(k,j,i),                      &
    2679                                      REAL( fill_value, KIND = wp ),              &
    2680                                      BTEST( wall_flags_total_0(k,j,i), flag_nr ) &
     3264                local_pf(i,j,k) = MERGE(                                       &
     3265                                     to_be_resorted(k,j,i),                    &
     3266                                     REAL( fill_value, KIND = wp ),            &
     3267                                     BTEST( wall_flags_total_0(k,j,i), flag_nr )     &
    26813268                                  )
    26823269             ENDDO
     
    27503337          CASE ( 'curvature_solution_effects_bulk' )
    27513338             READ ( 13 )  curvature_solution_effects_bulk
     3339
     3340          CASE ( 'microphysics_ice_extension' )
     3341             READ ( 13 )  microphysics_ice_extension
     3342
     3343          CASE ( 'ice_crystal_sedimentation' )
     3344             READ ( 13 )  ice_crystal_sedimentation
     3345
     3346          CASE ( 'in_init' )
     3347             READ ( 13 )  in_init
     3348
     3349          CASE ( 'start_ice_microphysics' )
     3350             READ ( 13 )  start_ice_microphysics
    27523351
    27533352          CASE DEFAULT
     
    27823381       CALL rrd_mpi_io( 'aerosol_bulk', aerosol_bulk )
    27833382       CALL rrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk )
     3383       CALL rrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics )
     3384       CALL rrd_mpi_io( 'microphysics_ice_extension', microphysics_ice_extension )
     3385       CALL rrd_mpi_io( 'in_init', in_init )
     3386       CALL rrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation )
     3387
     3388
    27843389
    27853390    END SUBROUTINE bcm_rrd_global_mpi
     
    28463451                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    28473452
     3453          CASE ( 'ni' )
     3454             IF ( k == 1 )  READ ( 13 )  tmp_3d
     3455             ni(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3456                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3457
     3458          CASE ( 'ni_av' )
     3459             IF ( .NOT. ALLOCATED( ni_av ) )  THEN
     3460                ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3461             ENDIF
     3462             IF ( k == 1 )  READ ( 13 )  tmp_3d
     3463             ni_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3464                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3465
    28483466          CASE ( 'nr' )
    28493467             IF ( k == 1 )  READ ( 13 )  tmp_3d
     
    28863504                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    28873505
     3506          CASE ( 'qi' )
     3507             IF ( k == 1 )  READ ( 13 )  tmp_3d
     3508             qi(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
     3509                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3510
     3511          CASE ( 'qi_av' )
     3512             IF ( .NOT. ALLOCATED( qi_av ) )  THEN
     3513                ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3514             ENDIF
     3515             IF ( k == 1 )  READ ( 13 )  tmp_3d
     3516             qi_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
     3517                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3518
    28883519          CASE ( 'qc_av' )
    28893520             IF ( .NOT. ALLOCATED( qc_av ) )  THEN
     
    29843615          CALL wrd_write_string( 'curvature_solution_effects_bulk' )
    29853616          WRITE ( 14 )  curvature_solution_effects_bulk
     3617
     3618          CALL wrd_write_string( 'start_ice_microphysics' )
     3619          WRITE ( 14 )  start_ice_microphysics
     3620
     3621          CALL wrd_write_string( 'microphysics_ice_extension' )
     3622          WRITE ( 14 )  microphysics_ice_extension
     3623
     3624          CALL wrd_write_string( 'in_init' )
     3625          WRITE ( 14 )  in_init
     3626
     3627          CALL wrd_write_string( 'ice_crystal_sedimentation' )
     3628          WRITE ( 14 )  ice_crystal_sedimentation
    29863629
    29873630       ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' )  THEN
     
    30013644          CALL wrd_mpi_io( 'aerosol_bulk', aerosol_bulk )
    30023645          CALL wrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk )
     3646          CALL wrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics )
     3647          CALL wrd_mpi_io( 'microphysics_ice_extension', microphysics_ice_extension )
     3648          CALL wrd_mpi_io( 'in_init', in_init )
     3649          CALL wrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation )
    30033650
    30043651       ENDIF
     
    30623709
    30633710          ENDIF
     3711
     3712          IF ( microphysics_ice_extension )  THEN
     3713
     3714             CALL wrd_write_string( 'ni' )
     3715             WRITE ( 14 )  ni
     3716
     3717             IF ( ALLOCATED( ni_av ) )  THEN
     3718                CALL wrd_write_string( 'ni_av' )
     3719                WRITE ( 14 )  ni_av
     3720             ENDIF
     3721
     3722             CALL wrd_write_string( 'qi' )
     3723             WRITE ( 14 )  qi
     3724
     3725             IF ( ALLOCATED( qi_av ) )  THEN
     3726                CALL wrd_write_string( 'qi_av' )
     3727                WRITE ( 14 )  qi_av
     3728             ENDIF
     3729
     3730          ENDIF
     3731
    30643732
    30653733          IF ( microphysics_seifert )  THEN
     
    31043772             IF ( ALLOCATED( qr_av ) )  CALL wrd_mpi_io( 'qr_av', qr_av )
    31053773          ENDIF
     3774          IF ( microphysics_ice_extension )  THEN
     3775             CALL wrd_mpi_io( 'ni', ni )
     3776             IF ( ALLOCATED( ni_av ) )  CALL wrd_mpi_io( 'ni_av', ni_av )
     3777             CALL wrd_mpi_io( 'qi', qi )
     3778             IF ( ALLOCATED( qi_av ) )  CALL wrd_mpi_io( 'qi_av', qi_av )
     3779          ENDIF
    31063780
    31073781       ENDIF
     
    31343808!--             Predetermine flag to mask topography
    31353809                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    3136 
    3137                 IF ( qr(k,j,i) <= eps_sb )  THEN
    3138                    qr(k,j,i) = 0.0_wp
    3139                    nr(k,j,i) = 0.0_wp
    3140                 ELSE
    3141                    IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
    3142                       nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
    3143                    ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
    3144                       nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
     3810                IF ( .NOT. microphysics_morrison_no_rain )  THEN
     3811                   IF ( qr(k,j,i) <= eps_sb )  THEN
     3812                      qr(k,j,i) = 0.0_wp
     3813                      nr(k,j,i) = 0.0_wp
     3814                   ELSE
     3815                      IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
     3816                         nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
     3817                      ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
     3818                         nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
     3819                      ENDIF
    31453820                   ENDIF
    31463821                ENDIF
     
    31893864          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    31903865
    3191           IF ( qr(k,j,i) <= eps_sb )  THEN
    3192              qr(k,j,i) = 0.0_wp
    3193              nr(k,j,i) = 0.0_wp
    3194           ELSE
    3195 !
    3196 !--          Adjust number of raindrops to avoid nonlinear effects in
    3197 !--          sedimentation and evaporation of rain drops due to too small or
    3198 !--          too big weights of rain drops (Stevens and Seifert, 2008).
    3199              IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
    3200                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
    3201              ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
    3202                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
    3203              ENDIF
    3204 
     3866          IF ( .NOT.  microphysics_morrison_no_rain )  THEN
     3867             IF ( qr(k,j,i) <= eps_sb )  THEN
     3868                qr(k,j,i) = 0.0_wp
     3869                nr(k,j,i) = 0.0_wp
     3870             ELSE
     3871!
     3872!--             Adjust number of raindrops to avoid nonlinear effects in
     3873!--             sedimentation and evaporation of rain drops due to too small or
     3874!--             too big weights of rain drops (Stevens and Seifert, 2008).
     3875                IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
     3876                   nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
     3877                ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
     3878                   nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
     3879                ENDIF
     3880             ENDIF
    32053881          ENDIF
    32063882
     
    32193895
    32203896    END SUBROUTINE adjust_cloud_ij
     3897
     3898!------------------------------------------------------------------------------!
     3899! Description:
     3900! ------------
     3901!> Adjust number of ice crystal to avoid nonlinear effects in sedimentation and
     3902!> evaporation of ice crytals due to too small or too big weights
     3903!> of ice crytals (Stevens and Seifert, 2008).
     3904!------------------------------------------------------------------------------!
     3905    SUBROUTINE adjust_ice
     3906
     3907       IMPLICIT NONE
     3908
     3909       INTEGER(iwp) ::  i                 !<
     3910       INTEGER(iwp) ::  j                 !<
     3911       INTEGER(iwp) ::  k                 !<
     3912
     3913       REAL(wp) ::  flag                  !< flag to indicate first grid level above
     3914
     3915       CALL cpu_log( log_point_s(97), 'adjust_ice', 'start' )
     3916
     3917       DO  i = nxl, nxr
     3918          DO  j = nys, nyn
     3919             DO  k = nzb+1, nzt
     3920!
     3921!--             Predetermine flag to mask topography
     3922                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3923                IF ( qi(k,j,i) <= ximin )  THEN
     3924                   qi(k,j,i) = 0.0_wp
     3925                   ni(k,j,i) = 0.0_wp
     3926                ELSE
     3927                   IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) )  THEN
     3928                      ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag
     3929                   ENDIF
     3930                ENDIF
     3931             ENDDO
     3932          ENDDO
     3933       ENDDO
     3934
     3935       CALL cpu_log( log_point_s(97), 'adjust_ice', 'stop' )
     3936
     3937    END SUBROUTINE adjust_ice
     3938
     3939!------------------------------------------------------------------------------!
     3940! Description:
     3941! ------------
     3942!> Adjust number of of ice crystal to avoid nonlinear effects in
     3943!> sedimentation and evaporation of ice crystals due to too small or
     3944!> too big weights of ice crytals (Stevens and Seifert, 2008).
     3945!> The same procedure is applied to cloud droplets if they are determined
     3946!> prognostically. Call for grid point i,j
     3947!------------------------------------------------------------------------------!
     3948    SUBROUTINE adjust_ice_ij( i, j )
     3949
     3950       IMPLICIT NONE
     3951
     3952       INTEGER(iwp) ::  i                 !<
     3953       INTEGER(iwp) ::  j                 !<
     3954       INTEGER(iwp) ::  k                 !<
     3955
     3956       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
     3957
     3958       DO  k = nzb+1, nzt
     3959!
     3960!--       Predetermine flag to mask topography
     3961          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3962          IF ( qi(k,j,i) <= ximin )  THEN
     3963             qi(k,j,i) = 0.0_wp
     3964             ni(k,j,i) = 0.0_wp
     3965          ELSE
     3966             IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) )  THEN
     3967                ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag
     3968             ENDIF
     3969          ENDIF
     3970       ENDDO
     3971
     3972    END SUBROUTINE adjust_ice_ij
     3973
    32213974
    32223975!------------------------------------------------------------------------------!
     
    34234176    END SUBROUTINE activation_ij
    34244177
     4178!------------------------------------------------------------------------------!
     4179! Description:
     4180! ------------
     4181!> Calculate ice nucleation by applying the deposition-condensation formula as
     4182!> given by Meyers et al 1992 and as described in Seifert and Beheng 2006
     4183!------------------------------------------------------------------------------!
     4184    SUBROUTINE ice_nucleation
     4185
     4186       INTEGER(iwp) ::  i             !< loop index
     4187       INTEGER(iwp) ::  j             !< loop index
     4188       INTEGER(iwp) ::  k             !< loop index
     4189
     4190       LOGICAL :: isdac = .TRUE.
     4191
     4192       REAL(wp) ::  a_m92 = -0.639_wp !< parameter for nucleation
     4193       REAL(wp) ::  b_m92 = 12.96_wp  !< parameter for nucleation
     4194       REAL(wp) ::  flag              !< flag to indicate first grid level
     4195       REAL(wp) ::  n_in              !< number of ice nucleii
     4196       REAL(wp) ::  nucle = 0.0_wp    !< nucleation rate
     4197
     4198       CALL cpu_log( log_point_s(93), 'ice_nucleation', 'start' )
     4199
     4200       IF  ( isdac )  THEN
     4201          DO  i = nxl, nxr
     4202             DO  j = nys, nyn
     4203                DO  k = nzb+1, nzt
     4204!
     4205!--                Predetermine flag to mask topography
     4206                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4207!
     4208!--                Call calculation of supersaturation located in subroutine
     4209                   CALL supersaturation_ice ( i, j, k )
     4210                   nucle = 0.0_wp
     4211                   !IF ( zu(k) >= 1500.0_wp ) CYCLE
     4212                   IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp  ) THEN
     4213!
     4214!--                   Calculate ice nucleation
     4215                      nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp )
     4216                      !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin
     4217                      ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init)
     4218                   ENDIF
     4219
     4220                ENDDO
     4221             ENDDO
     4222          ENDDO
     4223       ELSE
     4224          DO  i = nxl, nxr
     4225             DO  j = nys, nyn
     4226                DO  k = nzb+1, nzt
     4227!
     4228!--                Predetermine flag to mask topography
     4229                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4230!
     4231!--                Call calculation of supersaturation located in subroutine
     4232                   CALL supersaturation_ice ( i, j, k )
     4233                   nucle = 0.0_wp
     4234                   IF ( sat_ice > 0.0 ) THEN
     4235!
     4236!--                   Calculate ice nucleation
     4237                      n_in = in_init * EXP( a_m92 + b_m92 * sat_ice )
     4238                      nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp )
     4239                   ENDIF
     4240                   ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp )
     4241                ENDDO
     4242             ENDDO
     4243          ENDDO
     4244       ENDIF
     4245
     4246       CALL cpu_log( log_point_s(93), 'ice_nucleation', 'stop' )
     4247
     4248    END SUBROUTINE ice_nucleation
     4249
     4250
     4251!------------------------------------------------------------------------------!
     4252! Description:
     4253! ------------
     4254!> Calculate ice nucleation by applying the deposition-condensation formula as
     4255!> given by Meyers et al 1992 and as described in Seifert and Beheng 2006
     4256!------------------------------------------------------------------------------!
     4257    SUBROUTINE ice_nucleation_ij( i, j )
     4258
     4259       INTEGER(iwp) ::  i             !< loop index
     4260       INTEGER(iwp) ::  j             !< loop index
     4261       INTEGER(iwp) ::  k             !< loop index
     4262
     4263       LOGICAL :: isdac = .TRUE.
     4264
     4265       REAL(wp) ::  a_m92 = -0.639_wp !< parameter for nucleation
     4266       REAL(wp) ::  b_m92 = 12.96_wp  !< parameter for nucleation
     4267       REAL(wp) ::  flag              !< flag to indicate first grid level
     4268       REAL(wp) ::  n_in              !< number of ice nucleii
     4269       REAL(wp) ::  nucle = 0.0_wp    !< nucleation rate
     4270
     4271       IF ( isdac )  THEN
     4272          DO  k = nzb+1, nzt
     4273!
     4274!--          Predetermine flag to mask topography
     4275             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4276!
     4277!--          Call calculation of supersaturation located in subroutine
     4278             CALL supersaturation_ice ( i, j, k )
     4279             nucle = 0.0_wp
     4280             !IF ( zu(k) >= 1500.0_wp ) CYCLE
     4281             IF ( sat_ice >= 0.05_wp  .OR.  ql(k,j,i) >= 0.001E-3_wp ) THEN
     4282!
     4283!--             Calculate ice nucleation
     4284                nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp )
     4285                !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin
     4286                ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init)
     4287             ENDIF
     4288          ENDDO
     4289       ELSE
     4290          DO  k = nzb+1, nzt
     4291!
     4292!--          Predetermine flag to mask topography
     4293             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4294!
     4295!--          Call calculation of supersaturation located in subroutine
     4296             CALL supersaturation_ice ( i, j, k )
     4297             nucle = 0.0_wp
     4298             IF ( sat_ice > 0.0 ) THEN
     4299!
     4300!--             Calculate ice nucleation
     4301                n_in = in_init * EXP( a_m92 + b_m92 * sat_ice )
     4302                nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp )
     4303             ENDIF
     4304             ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp )
     4305          ENDDO
     4306       ENDIF
     4307
     4308
     4309    END SUBROUTINE ice_nucleation_ij
    34254310
    34264311!------------------------------------------------------------------------------!
     
    34344319       IMPLICIT NONE
    34354320
    3436        INTEGER(iwp) ::  i                 !<
    3437        INTEGER(iwp) ::  j                 !<
    3438        INTEGER(iwp) ::  k                 !<
    3439 
    3440        REAL(wp)     ::  cond              !<
    3441        REAL(wp)     ::  cond_max          !<
    3442        REAL(wp)     ::  dc                !<
    3443        REAL(wp)     ::  evap              !<
    3444        REAL(wp)     ::  g_fac             !<
    3445        REAL(wp)     ::  nc_0              !<
    3446        REAL(wp)     ::  temp              !<
    3447        REAL(wp)     ::  xc                !<
    3448 
    3449        REAL(wp) ::  flag                  !< flag to indicate first grid level above
     4321       INTEGER(iwp) ::  i                 !< loop index
     4322       INTEGER(iwp) ::  j                 !< loop index
     4323       INTEGER(iwp) ::  k                 !< loop index
     4324
     4325       REAL(wp)     ::  cond              !< condensation rate
     4326       REAL(wp)     ::  cond_max          !< maximum condensation rate
     4327       REAL(wp)     ::  dc                !< weight avageraed diameter
     4328       REAL(wp)     ::  evap              !< evaporation rate
     4329       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
     4330       REAL(wp)     ::  g_fac             !< factor 1 / Fk + Fd
     4331       REAL(wp)     ::  nc_0              !< integral diameter
     4332       REAL(wp)     ::  temp              !< actual temperature
     4333       REAL(wp)     ::  xc                !< mean mass droplets
    34504334
    34514335       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
     
    34614345                CALL supersaturation ( i, j, k )
    34624346!
    3463 !--             Actual temperature:
    3464                 IF ( microphysics_seifert ) THEN
    3465                    temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
    3466                 ELSEIF ( microphysics_morrison_no_rain ) THEN
    3467                    temp = t_l + lv_d_cp * qc(k,j,i)
    3468                 ENDIF 
     4347!--             Actual temperature, t_l is calculated directly before
     4348!--             in supersaturation
     4349                IF ( microphysics_ice_extension ) THEN
     4350                   temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
     4351                ELSE
     4352                   temp = t_l + lv_d_cp * ql(k,j,i)
     4353                ENDIF
    34694354
    34704355                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
     
    35254410       IMPLICIT NONE
    35264411
    3527        INTEGER(iwp) ::  i                 !<
    3528        INTEGER(iwp) ::  j                 !<
    3529        INTEGER(iwp) ::  k                 !<
    3530 
    3531        REAL(wp)     ::  cond              !<
    3532        REAL(wp)     ::  cond_max          !<
    3533        REAL(wp)     ::  dc                !<
    3534        REAL(wp)     ::  evap              !<
     4412       INTEGER(iwp) ::  i                 !< loop index
     4413       INTEGER(iwp) ::  j                 !< loop index
     4414       INTEGER(iwp) ::  k                 !< loop index
     4415
     4416       REAL(wp)     ::  cond              !< condensation rate
     4417       REAL(wp)     ::  cond_max          !< maximum condensation rate
     4418       REAL(wp)     ::  dc                !< weight avageraed diameter
     4419       REAL(wp)     ::  evap              !< evaporation rate
    35354420       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
    3536        REAL(wp)     ::  g_fac             !<
    3537        REAL(wp)     ::  nc_0              !<
    3538        REAL(wp)     ::  temp              !<
    3539        REAL(wp)     ::  xc                !<
    3540 
     4421       REAL(wp)     ::  g_fac             !< factor 1 / Fk + Fd
     4422       REAL(wp)     ::  nc_0              !< integral diameter
     4423       REAL(wp)     ::  temp              !< actual temperature
     4424       REAL(wp)     ::  xc                !< mean mass droplets
    35414425
    35424426       DO  k = nzb+1, nzt
     
    35484432          CALL supersaturation ( i, j, k )
    35494433!
    3550 !--       Actual temperature:
    3551           IF ( microphysics_seifert ) THEN
    3552              temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
    3553           ELSEIF ( microphysics_morrison_no_rain ) THEN
    3554              temp = t_l + lv_d_cp * qc(k,j,i)
    3555           ENDIF 
    3556 
    3557           g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
    3558                               l_v / ( thermal_conductivity_l * temp )    &
    3559                               + r_v * temp / ( diff_coeff_l * e_s )      &
     4434!--       Actual temperature, t_l is calculated directly before
     4435!--       in supersaturation
     4436          IF ( microphysics_ice_extension ) THEN
     4437             temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
     4438          ELSE
     4439             temp = t_l + lv_d_cp * ql(k,j,i)
     4440          ENDIF
     4441
     4442          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *              &
     4443                              l_v / ( thermal_conductivity_l * temp )          &
     4444                              + r_v * temp / ( diff_coeff_l * e_s )            &
    35604445                            )
    35614446!
     
    35944479
    35954480    END SUBROUTINE condensation_ij
     4481
     4482!------------------------------------------------------------------------------!
     4483! Description:
     4484! ------------
     4485!> Calculate the growth of ice particles by water vapor deposition (after
     4486!> Seifert and Beheng, 2006).
     4487!------------------------------------------------------------------------------!
     4488    SUBROUTINE ice_deposition
     4489
     4490       INTEGER(iwp) ::  i                 !< loop index
     4491       INTEGER(iwp) ::  j                 !< loop index
     4492       INTEGER(iwp) ::  k                 !< loop index
     4493
     4494       REAL(wp) ::  ac = 0.09_wp        !< parameter for ice capacitance
     4495       REAL(wp) ::  bc = 0.33_wp        !< parameter for ice capacitance
     4496       REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral
     4497                                        !< distribution, here following gamma
     4498                                        !< size distribution with µ =1/3 and nu=0
     4499       REAL(wp) ::  deposition_rate     !< depositions rate
     4500       REAL(wp) ::  deposition_rate_max !< maximum deposition rate
     4501       REAL(wp) ::  sublimation_rate    !< sublimations rate
     4502       REAL(wp) ::  gfac_dep            !< factor
     4503       REAL(wp) ::  temp                !< actual temperature
     4504       REAL(wp) ::  xi                  !< mean mass of ice crystal
     4505       REAL(wp) ::  flag                !< flag to indicate first grid level above
     4506
     4507       CALL cpu_log( log_point_s(95), 'ice deposition', 'start' )
     4508
     4509       DO  i = nxl, nxr
     4510          DO  j = nys, nyn
     4511             DO  k = nzb+1, nzt
     4512!
     4513!--             Predetermine flag to mask topography
     4514                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4515!
     4516!--             Call calculation of supersaturation over a plane ice surface
     4517                CALL supersaturation_ice ( i, j, k )
     4518!
     4519!--             Actual temperature:
     4520                temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
     4521
     4522                IF ( temp >= 273.15_wp ) CYCLE
     4523!
     4524!--             calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g.
     4525!--             Rogers and Yau, 1989
     4526                gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *     &
     4527                                    l_s / ( thermal_conductivity_l * temp )    &
     4528                                    + r_v * temp / ( diff_coeff_l * e_si )     &
     4529                                    )
     4530!
     4531!--             If there is nothing nucleated, than there is also no
     4532!--             deposition (above -38°C)
     4533                IF ( ni(k,j,i) <= 0.0_wp ) CYCLE
     4534!
     4535!--             calculate mean mass of ice crystal
     4536                xi = 0.0_wp
     4537                xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin )
     4538!
     4539!--             Condensation needs only to be calculated in supersaturated
     4540!--             regions (regarding ice)
     4541                IF ( sat_ice > 0.0_wp )  THEN
     4542!
     4543!--                Calculate deposition rate assuming ice crystal shape as
     4544!--                prescribed in Ovchinnikov et al., 2014 and a gamma size
     4545!--                distribution according to Seifert and Beheng with to
     4546!--                µ =1/3 and nu=0
     4547                   deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep *       &
     4548                                      fac_gamma * ac * xi**bc * ni(k,j,i)
     4549                   IF ( microphysics_seifert ) THEN
     4550                      deposition_rate_max  = q(k,j,i) -                        &
     4551                                             q_si - qr(k,j,i) - qi(k,j,i)
     4552                   ELSEIF ( microphysics_morrison_no_rain ) THEN
     4553                      deposition_rate_max  = q(k,j,i) -                        &
     4554                                             q_si - qc(k,j,i) - qi(k,j,i)
     4555                   ENDIF
     4556                   deposition_rate = MIN( deposition_rate,                     &
     4557                                          deposition_rate_max / dt_micro )
     4558
     4559                   qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag
     4560                ELSEIF ( sat_ice < 0.0_wp ) THEN
     4561                   sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep *       &
     4562                                      fac_gamma * ac * xi**bc * ni(k,j,i)
     4563                   sublimation_rate = MAX( sublimation_rate,                   &
     4564                                           -qi(k,j,i) / dt_micro )
     4565                   qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag
     4566                ENDIF
     4567             ENDDO
     4568          ENDDO
     4569       ENDDO
     4570
     4571       CALL cpu_log( log_point_s(95), 'ice deposition', 'stop' )
     4572
     4573    END SUBROUTINE ice_deposition
     4574
     4575!------------------------------------------------------------------------------!
     4576! Description:
     4577! ------------
     4578!> Calculate condensation rate for cloud water content (after Khairoutdinov and
     4579!> Kogan, 2000).
     4580!------------------------------------------------------------------------------!
     4581    SUBROUTINE ice_deposition_ij( i, j )
     4582
     4583       INTEGER(iwp) ::  i                 !< loop index
     4584       INTEGER(iwp) ::  j                 !< loop index
     4585       INTEGER(iwp) ::  k                 !< loop index
     4586
     4587       REAL(wp) ::  ac = 0.09_wp        !< parameter for ice capacitance
     4588       REAL(wp) ::  bc = 0.33_wp        !< parameter for ice capacitance
     4589       REAL(wp) ::  fac_gamma = 0.76_wp !< parameter to describe spectral
     4590                                        !< distribution, here following gamma
     4591                                        !< size distribution with nu=1, v=0
     4592       REAL(wp) ::  deposition_rate     !< depositions rate
     4593       REAL(wp) ::  deposition_rate_max !< maximum deposition rate
     4594       REAL(wp) ::  sublimation_rate    !< sublimations rate
     4595       REAL(wp) ::  gfac_dep            !< factor
     4596       REAL(wp) ::  temp                !< actual temperature
     4597       REAL(wp) ::  xi                  !< mean mass of ice crystal
     4598       REAL(wp) ::  flag                !< flag to indicate first grid level above
     4599
     4600       DO  k = nzb+1, nzt
     4601!
     4602!--       Predetermine flag to mask topography
     4603          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4604!
     4605!--       Call calculation of supersaturation over a plane ice surface
     4606          CALL supersaturation_ice ( i, j, k )
     4607!
     4608!--       Actual temperature:
     4609          temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
     4610
     4611          IF ( temp >= 273.15_wp ) CYCLE
     4612!
     4613!--       calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g.
     4614!--       Rogers and Yau, 1989
     4615          gfac_dep  = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) *           &
     4616                              l_s / ( thermal_conductivity_l * temp )          &
     4617                              + r_v * temp / ( diff_coeff_l * e_si )           &
     4618                              )
     4619!
     4620!--       If there is nothing nucleated, than there is also no
     4621!--       deposition (above -38°C)
     4622          IF ( ni(k,j,i) <= 0.0_wp ) CYCLE
     4623!
     4624!--       calculate mean mass of ice crystal
     4625          xi = 0.0_wp
     4626          xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin )
     4627!
     4628!--       Condensation needs only to be calculated in supersaturated
     4629!--       regions (regarding ice)
     4630          IF ( sat_ice > 0.0_wp )  THEN
     4631!
     4632!--          Calculate deposition rate assuming ice crystal shape as
     4633!--          prescribed in Ovchinnikov et al., 2014 and a gamma size
     4634!--          distribution according to Seifert and Beheng with to
     4635!--          µ =1/3 and nu=0
     4636             deposition_rate  = 4.0_wp * pi * sat_ice * gfac_dep *             &
     4637                                fac_gamma * ac * xi**bc * ni(k,j,i)
     4638             IF ( microphysics_seifert ) THEN
     4639                deposition_rate_max  = q(k,j,i) -                              &
     4640                                       q_si - qr(k,j,i) - qi(k,j,i)
     4641             ELSEIF ( microphysics_morrison_no_rain ) THEN
     4642                deposition_rate_max  = q(k,j,i) -                              &
     4643                                       q_si - qc(k,j,i) - qi(k,j,i)
     4644             ENDIF
     4645             deposition_rate = MIN( deposition_rate,                           &
     4646                                    deposition_rate_max / dt_micro )
     4647
     4648             qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag
     4649          ELSEIF ( sat_ice < 0.0_wp ) THEN
     4650             sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep *             &
     4651                                fac_gamma * ac * xi**bc * ni(k,j,i)
     4652             sublimation_rate = MAX( sublimation_rate,                         &
     4653                                     -qi(k,j,i) / dt_micro )
     4654             qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag
     4655          ENDIF
     4656       ENDDO
     4657
     4658    END SUBROUTINE ice_deposition_ij
    35964659
    35974660
     
    38234886!--          Autoconversion rate (Seifert and Beheng, 2006):
    38244887             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
    3825                        ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
     4888                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *           &
    38264889                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
    38274890                       rho_surface
     
    40555118          ENDIF
    40565119
    4057           IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
     5120          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.    &
    40585121               ( nc_accr > eps_mr ) )  THEN
    40595122!
     
    40825145!
    40835146!--          Accretion rate (Seifert and Beheng, 2006):
    4084              accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
     5147             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                    &
    40855148                    SQRT( rho_surface * hyrho(k) )
    40865149             accr = MIN( accr, qc(k,j,i) / dt_micro )
     
    40895152             qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
    40905153             IF ( microphysics_morrison )  THEN
    4091                 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
     5154                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *            &
    40925155                                             hyrho(k) * dt_micro * flag        &
    40935156                                         )
     
    45875650             IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
    45885651                sed_nc(k) = sed_qc_const *                                     &
    4589                             ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
     5652                            ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *    &
    45905653                            ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
    45915654             ELSE
     
    45945657
    45955658             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
    4596                               nc(k,j,i) / dt_micro + sed_nc(k+1)                &
     5659                              nc(k,j,i) / dt_micro + sed_nc(k+1)               &
    45975660                            ) * flag
    45985661
    4599              nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *               &
     5662             nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *             &
    46005663                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
    46015664          ENDIF
     
    46085671          ENDIF
    46095672
    4610           sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /          &
     5673          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /         &
    46115674                                      dt_micro + sed_qc(k+1)                   &
    46125675                         ) * flag
    46135676
    4614           q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
     5677          q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /    &
    46155678                                hyrho(k) * dt_micro * flag
    4616           qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
     5679          qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /    &
    46175680                                hyrho(k) * dt_micro * flag
    4618           pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
     5681          pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /    &
    46195682                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
    46205683                                                                * flag
     
    46325695
    46335696    END SUBROUTINE sedimentation_cloud_ij
     5697
     5698!------------------------------------------------------------------------------!
     5699! Description:
     5700! ------------
     5701!> Sedimentation of ice crystals
     5702!------------------------------------------------------------------------------!
     5703    SUBROUTINE sedimentation_ice
     5704
     5705       INTEGER(iwp) ::  i             !< loop index
     5706       INTEGER(iwp) ::  j             !< loop index
     5707       INTEGER(iwp) ::  k             !< loop index
     5708
     5709       REAL(wp) ::  flag              !< flag to indicate first grid level
     5710       REAL(wp) ::  av = 6.39_wp      !< parameter for calculating fall speed
     5711       REAL(wp) ::  bv = 0.1666_wp    !< parameter (1/6)
     5712       REAL(wp) ::  xi = 0.0_wp       !< mean mass of ice crystal
     5713       REAL(wp) ::  vi = 0.0_wp       !< mean fall speed of ice crystal
     5714       REAL(wp) ::  factor_sed_gamma_k0 = 0.76_wp !< factor for zeroth moment and
     5715                                                  !< µ =1/3 and nu=0
     5716       REAL(wp) ::  factor_sed_gamma_k1 = 1.61_wp !< factor for first moment and
     5717                                                  !< µ =1/3 and nu=0
     5718
     5719       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_ni  !< sedimentation rate zeroth moment
     5720       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qi  !< sedimentation rate fist moment
     5721
     5722       CALL cpu_log( log_point_s(96), 'sed_ice', 'start' )
     5723
     5724       sed_qi(nzt+1) = 0.0_wp
     5725       sed_ni(nzt+1) = 0.0_wp
     5726
     5727       DO  i = nxl, nxr
     5728          DO  j = nys, nyn
     5729             DO  k = nzt, nzb+1, -1
     5730
     5731                IF ( ni(k,j,i) <= 0.0_wp ) THEN
     5732                   xi = 0.0_wp
     5733                ELSE
     5734!--                Calculate mean mass of ice crystal
     5735                   xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin)
     5736                ENDIF
     5737!
     5738!--             calculate fall speed of ice crystal
     5739                vi = av * xi**bv
     5740!
     5741!--             Predetermine flag to mask topography
     5742                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     5743!
     5744!--             Calculate sedimentation rate for each grid box, factors are
     5745!--             calculated using
     5746                IF ( qi(k,j,i) > eps_sb  .AND.  ni(k,j,i) >= 0.0_wp )  THEN
     5747                   sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag
     5748                   sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag
     5749                ELSE
     5750                   sed_qi(k) = 0.0_wp
     5751                   sed_ni(k) = 0.0_wp
     5752                ENDIF
     5753!
     5754!--             Calculate sedimentation: divergence of sedimentation flux
     5755                sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
     5756                                            dt_micro + sed_qi(k+1)             &
     5757                               ) * flag
     5758
     5759                q(k,j,i)  = q(k,j,i)  + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)&
     5760                                      / hyrho(k) * dt_micro * flag
     5761                qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)&
     5762                                      / hyrho(k) * dt_micro * flag
     5763                ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1)&
     5764                                      / hyrho(k) * dt_micro * flag
     5765
     5766                pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)&
     5767                                      / hyrho(k) * l_s / c_p * d_exner(k) *    &
     5768                                        dt_micro * flag
     5769!
     5770!--             Compute the precipitation rate of cloud (fog) droplets
     5771                IF ( call_microphysics_at_all_substeps )  THEN
     5772                   prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) *            &
     5773                                weight_substep(intermediate_timestep_count) *  &
     5774                                flag
     5775                ELSE
     5776                   prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag
     5777                ENDIF
     5778
     5779             ENDDO
     5780          ENDDO
     5781       ENDDO
     5782
     5783       CALL cpu_log( log_point_s(96), 'sed_ice', 'stop' )
     5784
     5785    END SUBROUTINE sedimentation_ice
     5786
     5787
     5788!------------------------------------------------------------------------------!
     5789! Description:
     5790! ------------
     5791!> Sedimentation of ice crystals
     5792!------------------------------------------------------------------------------!
     5793    SUBROUTINE sedimentation_ice_ij( i, j )
     5794
     5795       INTEGER(iwp) ::  i             !<
     5796       INTEGER(iwp) ::  j             !<
     5797       INTEGER(iwp) ::  k             !<
     5798
     5799       REAL(wp) ::  flag           !< flag to indicate first grid level
     5800       REAL(wp) ::  av = 6.39_wp   !< parameter for calculating fall speed
     5801       REAL(wp) ::  bv = 0.1666_wp !<
     5802       REAL(wp) ::  xi = 0.0_wp    !< mean mass of ice crystal
     5803       REAL(wp) ::  vi = 0.0_wp    !< mean fall speed of ice crystal
     5804       REAL(wp) ::  factor_sed_gamma_k0 = 0.76_wp
     5805       REAL(wp) ::  factor_sed_gamma_k1 = 1.61_wp
     5806
     5807       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_ni  !<
     5808       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qi  !<
     5809
     5810       sed_qi(nzt+1) = 0.0_wp
     5811       sed_ni(nzt+1) = 0.0_wp
     5812
     5813       DO  k = nzt, nzb+1, -1
     5814          IF ( ni(k,j,i) <= 0.0_wp ) THEN
     5815             xi = 0.0_wp
     5816          ELSE
     5817!--          Calculate mean mass of ice crystal
     5818             xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin)
     5819          ENDIF
     5820!
     5821!--       calculate fall speed of ice crystal
     5822          vi = av * xi**bv
     5823!
     5824!--       Predetermine flag to mask topography
     5825          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     5826!
     5827!--       Calculate sedimentation rate for each grid box, factors are
     5828!--       calculated using
     5829          IF ( qi(k,j,i) > eps_sb  .AND.  ni(k,j,i) >= 0.0_wp )  THEN
     5830             sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag
     5831             sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag
     5832          ELSE
     5833             sed_qi(k) = 0.0_wp
     5834             sed_ni(k) = 0.0_wp
     5835          ENDIF
     5836!
     5837!--       calculate fall speed of ice crystal
     5838          vi = av * xi**bv
     5839!
     5840!--       Predetermine flag to mask topography
     5841          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     5842
     5843          IF ( qi(k,j,i) > eps_sb  .AND.  ni(k,j,i) >= 0.0_wp )  THEN
     5844             sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag
     5845             sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag
     5846          ELSE
     5847             sed_qi(k) = 0.0_wp
     5848             sed_ni(k) = 0.0_wp
     5849          ENDIF
     5850
     5851          sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) /         &
     5852                                      dt_micro + sed_qi(k+1)                   &
     5853                         ) * flag
     5854
     5855          q(k,j,i)  = q(k,j,i)  + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) /    &
     5856                                hyrho(k) * dt_micro * flag
     5857          qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) /    &
     5858                                hyrho(k) * dt_micro * flag
     5859          ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1) /    &
     5860                                hyrho(k) * dt_micro * flag
     5861          pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) /    &
     5862                                hyrho(k) * l_s / c_p * d_exner(k) * dt_micro   &
     5863                                                                  * flag
     5864!
     5865!--       Compute the precipitation rate of cloud (fog) droplets
     5866          IF ( call_microphysics_at_all_substeps )  THEN
     5867             prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) *                  &
     5868                              weight_substep(intermediate_timestep_count) * flag
     5869          ELSE
     5870             prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag
     5871          ENDIF
     5872
     5873       ENDDO
     5874
     5875    END SUBROUTINE sedimentation_ice_ij
     5876
    46345877
    46355878
     
    50536296
    50546297          sed_nr(k) = flux / dt_micro * flag
    5055           nr(k,j,i)  = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
     6298          nr(k,j,i)  = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /   &
    50566299                                    hyrho(k) * dt_micro * flag
    50576300!
     
    50806323          sed_qr(k) = flux / dt_micro * flag
    50816324
    5082           qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
     6325          qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /    &
    50836326                                hyrho(k) * dt_micro * flag
    5084           q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
     6327          q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /    &
    50856328                                hyrho(k) * dt_micro * flag
    5086           pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
     6329          pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /    &
    50876330                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
    50886331                                                                * flag
     
    52026445!--    (see: Cuijpers + Duynkerke, 1993, JAS, 23)
    52036446       q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
     6447
     6448       IF ( .NOT. microphysics_ice_extension ) THEN
     6449          IF ( microphysics_seifert ) THEN
     6450             sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
     6451          ELSEIF ( microphysics_morrison_no_rain ) THEN
     6452             sat = ( q(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
     6453          ENDIF
     6454       ELSE
     6455          IF ( microphysics_seifert ) THEN
     6456             sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp
     6457          ELSEIF ( microphysics_morrison_no_rain ) THEN
     6458             sat = ( q(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp
     6459          ENDIF
     6460       ENDIF
     6461
     6462    END SUBROUTINE supersaturation
     6463
     6464!------------------------------------------------------------------------------!
     6465! Description:
     6466! ------------
     6467!> Computation of the diagnostic supersaturation sat, actual liquid water
     6468!< temperature t_l and saturation water vapor mixing ratio q_s
     6469!------------------------------------------------------------------------------!
     6470    SUBROUTINE supersaturation_ice ( i, j, k )
     6471
     6472       INTEGER(iwp) ::  i   !< running index
     6473       INTEGER(iwp) ::  j   !< running index
     6474       INTEGER(iwp) ::  k   !< running index
     6475
     6476       REAL(wp) ::  e_a     !< water vapor pressure
     6477       REAL(wp) ::  temp    !< actual temperature
     6478!
     6479!--    Actual liquid water temperature:
     6480       t_l = exner(k) * pt(k,j,i)
     6481!
     6482!--    Actual temperature:
     6483       temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i)
     6484!
     6485!--    Calculate water vapor saturation pressure with formular using actual
     6486!--    temperature
     6487       e_si = magnus_ice( temp )
     6488!
     6489!--    Computation of ice saturation mixing ratio:
     6490       q_si = rd_d_rv * e_si / ( hyp(k) - e_si )
     6491!
     6492!--    Current vapor pressure
     6493       e_a = (   q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) * hyp(k) /     &
     6494             ( ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) + rd_d_rv )
    52046495!
    52056496!--    Supersaturation:
    52066497!--    Not in case of microphysics_kessler or microphysics_sat_adjust
    52076498!--    since qr is unallocated
    5208        IF ( microphysics_seifert ) THEN
    5209           sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
    5210        ELSEIF ( microphysics_morrison_no_rain ) THEN
    5211           sat = ( q(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
    5212        ENDIF
    5213 
    5214     END SUBROUTINE supersaturation
     6499       sat_ice = e_a / e_si - 1.0_wp
     6500
     6501    END SUBROUTINE supersaturation_ice
    52156502
    52166503
     
    52246511    SUBROUTINE calc_liquid_water_content
    52256512
    5226 
    5227 
    5228        IMPLICIT NONE
    5229 
    52306513       INTEGER(iwp) ::  i !<
    52316514       INTEGER(iwp) ::  j !<
     
    52346517       REAL(wp)     ::  flag !< flag to indicate first grid level above surface
    52356518
    5236 
    52376519       DO  i = nxlg, nxrg
    52386520          DO  j = nysg, nyng
     
    52416523!--             Predetermine flag to mask topography
    52426524                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    5243 
    52446525!
    52456526!--             Call calculation of supersaturation located
    52466527                CALL supersaturation( i, j, k )
    5247 
    52486528!
    52496529!--             Compute the liquid water content
    5250                 IF ( microphysics_seifert  .AND. .NOT. microphysics_morrison ) &
    5251                 THEN
    5252                    IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
    5253                       qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
    5254                       ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
     6530                IF ( .NOT. microphysics_ice_extension ) THEN
     6531                   IF ( microphysics_seifert  .AND. .NOT.                      &
     6532                        microphysics_morrison )  THEN
     6533!--                   Seifert and Beheng scheme: saturation adjustment
     6534                      IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
     6535                         qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
     6536                         ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
     6537                      ELSE
     6538                         IF ( q(k,j,i) < qr(k,j,i) )  q(k,j,i) = qr(k,j,i)
     6539                         qc(k,j,i) = 0.0_wp
     6540                         ql(k,j,i) = qr(k,j,i) * flag
     6541                      ENDIF
     6542!
     6543!--                Morrison scheme: explicit condensation (see above)
     6544                   ELSEIF ( microphysics_morrison  .AND.                       &
     6545                            microphysics_seifert )  THEN
     6546                      ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
     6547!--                Morrison without rain: explicit condensation
     6548                   ELSEIF ( microphysics_morrison  .AND.                       &
     6549                            .NOT. microphysics_seifert )  THEN
     6550                      ql(k,j,i) = qc(k,j,i)
     6551!--                Kessler and saturation adjustment scheme
    52556552                   ELSE
    5256                       IF ( q(k,j,i) < qr(k,j,i) )  q(k,j,i) = qr(k,j,i)
    5257                       qc(k,j,i) = 0.0_wp
    5258                       ql(k,j,i) = qr(k,j,i) * flag
     6553                      IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
     6554                         qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
     6555                         ql(k,j,i) = qc(k,j,i) * flag
     6556                      ELSE
     6557                         qc(k,j,i) = 0.0_wp
     6558                         ql(k,j,i) = 0.0_wp
     6559                      ENDIF
    52596560                   ENDIF
    5260                 ELSEIF ( microphysics_morrison  .AND.  microphysics_seifert )  THEN
    5261                    ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
    5262                 ELSEIF ( microphysics_morrison  .AND.  .NOT. microphysics_seifert )  THEN
    5263                    ql(k,j,i) = qc(k,j,i)                   
     6561!--             Calculations of liquid water content in case of mixed-phase
     6562!--             cloud microphyics
    52646563                ELSE
    5265                    IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
    5266                       qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
    5267                       ql(k,j,i) = qc(k,j,i) * flag
     6564                   IF ( microphysics_seifert  .AND.                            &
     6565                        .NOT. microphysics_morrison ) &
     6566                   THEN
     6567!
     6568!--                   Seifert and Beheng scheme: saturation adjustment
     6569                      IF ( ( q(k,j,i)                                          &
     6570                             - q_s - qr(k,j,i) - qi(k,j,i) ) > 0.0_wp )  THEN
     6571                         qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) - qi(k,j,i) )&
     6572                                       * flag
     6573                         ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
     6574                      ELSE
     6575                         IF ( q(k,j,i) < ( qr(k,j,i) + qi(k,j,i) ) )  THEN
     6576                            q(k,j,i) = qr(k,j,i) + qi(k,j,i)
     6577                         ENDIF
     6578                         qc(k,j,i) = 0.0_wp
     6579                         ql(k,j,i) = qr(k,j,i) * flag
     6580                      ENDIF
     6581!--                Morrison scheme: explicit condensation (see above)
     6582                   ELSEIF ( microphysics_morrison  .AND.                       &
     6583                            microphysics_seifert )  THEN
     6584                      ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
     6585!--                Morrison without rain: explicit condensation
     6586                   ELSEIF ( microphysics_morrison  .AND.                       &
     6587                            .NOT. microphysics_seifert )  THEN
     6588                      ql(k,j,i) = qc(k,j,i)
     6589!--                Kessler and saturation adjustment scheme
    52686590                   ELSE
    5269                       qc(k,j,i) = 0.0_wp
    5270                       ql(k,j,i) = 0.0_wp
     6591                      IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
     6592                         qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
     6593                         ql(k,j,i) = qc(k,j,i) * flag
     6594                      ELSE
     6595                         qc(k,j,i) = 0.0_wp
     6596                         ql(k,j,i) = 0.0_wp
     6597                      ENDIF
    52716598                   ENDIF
    52726599                ENDIF
  • palm/trunk/SOURCE/compute_vpt.f90

    r4360 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4360 2020-01-07 11:25:50Z suehring
    2730! Corrected "Former revisions" section
    2831!
     
    4245
    4346    USE arrays_3d,                                                             &
    44         ONLY:  pt, q, ql, vpt, d_exner
     47        ONLY:  pt, q, ql, vpt, d_exner, qi
    4548
    4649    USE basic_constants_and_equations_mod,                                     &
    47         ONLY:  lv_d_cp
     50        ONLY:  lv_d_cp, ls_d_cp
    4851
    4952    USE control_parameters,                                                    &
     
    5659
    5760    USE bulk_cloud_model_mod,                                                  &
    58         ONLY:  bulk_cloud_model
     61        ONLY:  bulk_cloud_model, microphysics_ice_extension
    5962
    6063    IMPLICIT NONE
     
    6467    IF ( .NOT. bulk_cloud_model  .AND.  .NOT. cloud_droplets )  THEN
    6568       vpt = pt * ( 1.0_wp + 0.61_wp * q )
    66     ELSE IF (bulk_cloud_model)  THEN
     69    ELSEIF ( bulk_cloud_model  .AND.  .NOT. microphysics_ice_extension )  THEN
    6770       DO  k = nzb, nzt+1
    68           vpt(k,:,:) = ( pt(k,:,:) + d_exner(k) * lv_d_cp * ql(k,:,:) ) *      &
    69                        ( 1.0_wp + 0.61_wp * q(k,:,:) - 1.61_wp * ql(k,:,:) )
     71              vpt(k,:,:) = ( pt(k,:,:) + d_exner(k) * lv_d_cp * ql(k,:,:) ) *  &
     72                       ( 1.0_wp + 0.61_wp * q(k,:,:) - 1.61_wp *  ql(k,:,:)  )
     73       ENDDO
     74    ELSEIF ( bulk_cloud_model  .AND.  microphysics_ice_extension ) THEN
     75       DO  k = nzb, nzt+1
     76          vpt(k,:,:) = ( pt(k,:,:) + d_exner(k) * lv_d_cp * ql(k,:,:)   +      &
     77                                     d_exner(k) * ls_d_cp * qi(k,:,:) ) *      &
     78                       ( 1.0_wp + 0.61_wp * q(k,:,:) - 1.61_wp          *      &
     79                                           ( ql(k,:,:) + qi(k,:,:) ) )
    7080       ENDDO
    7181    ELSE
  • palm/trunk/SOURCE/flow_statistics.f90

    r4472 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4472 2020-03-24 12:21:00Z Giersch
    2730! Calculations of the Kolmogorov lengt scale eta implemented
    2831!
     
    8487    USE arrays_3d,                                                             &
    8588        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
    86                momentumflux_output_conversion, nc, nr, p, prho, prr, pt, q,    &
    87                qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
     89               momentumflux_output_conversion, nc, ni, nr, p, prho, prr, pt, q,&
     90               qc, qi, ql, qr, rho_air, rho_air_zw, rho_ocean, s,              &
    8891               sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion,  &
    8992               zw, d_exner
     
    9396
    9497    USE bulk_cloud_model_mod,                                                  &
    95         ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
     98        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert,   &
     99              microphysics_ice_extension
    96100
    97101    USE chem_modules,                                                          &
     
    14141418                                                                flag
    14151419                         ENDIF
     1420                         IF ( microphysics_ice_extension )  THEN
     1421                            sums_l(k,124,tn) = sums_l(k,124,tn) + ni(k,j,i) *  &
     1422                                                                rmask(j,i,sr) *&
     1423                                                                flag
     1424                            sums_l(k,125,tn) = sums_l(k,125,tn) + qi(k,j,i) *  &
     1425                                                                rmask(j,i,sr) *&
     1426                                                                flag
     1427                         ENDIF
     1428
    14161429                         IF ( microphysics_seifert )  THEN
    14171430                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
     
    19151928             sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
    19161929             sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
    1917              sums(k,123)           = sums(k,123) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
     1930             sums(k,123:125)       = sums(k,123:125) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
    19181931          ENDIF
    19191932       ENDDO
     
    20292042       hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
    20302043       hom(:,1,123,sr) = sums(:,123)   ! nc
     2044       hom(:,1,124,sr) = sums(:,124)   ! ni
     2045       hom(:,1,125,sr) = sums(:,125)   ! qi
    20312046       hom(:,1,73,sr) = sums(:,73)     ! nr
    20322047       hom(:,1,74,sr) = sums(:,74)     ! qr
  • palm/trunk/SOURCE/init_masks.f90

    r4444 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4444 2020-03-05 15:59:50Z raasch
    2730! bugfix: cpp-directives for serial mode added
    2831!
     
    5861
    5962    USE bulk_cloud_model_mod,                                                  &
    60         ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
     63        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert,   &
     64              microphysics_ice_extension
    6165
    6266    USE control_parameters,                                                    &
     
    268272                unit = '1/m3'
    269273
     274             CASE ( 'ni' )
     275                IF ( .NOT. bulk_cloud_model )  THEN
     276                   WRITE ( message_string, * ) 'output of "', TRIM( var ),     &
     277                        '" requires bulk_cloud_model = .TRUE.'
     278                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
     279                 ELSEIF ( .NOT. microphysics_ice_extension ) THEN
     280                   message_string = 'output of "' // TRIM( var ) // '" ' //    &
     281                         'requires  microphysics_ice_extension = .TRUE.'
     282                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     283                ENDIF
     284                unit = '1/m3'
     285
    270286             CASE ( 'nr' )
    271287                IF ( .NOT. bulk_cloud_model )  THEN
     
    331347                        '" requires bulk_cloud_model = .TRUE.'
    332348                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
     349                ENDIF
     350                unit = 'kg/kg'
     351
     352             CASE ( 'qi' )
     353                IF ( .NOT. bulk_cloud_model )  THEN
     354                   message_string = 'output of "' // TRIM( var ) // '" ' //    &
     355                            'requires bulk_cloud_model = .TRUE.'
     356                   CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
     357                ELSEIF ( .NOT. microphysics_ice_extension ) THEN
     358                   message_string = 'output of "' // TRIM( var ) // '" ' //    &
     359                            'requires microphysics_ice_extension = .TRUE.'
     360                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
    333361                ENDIF
    334362                unit = 'kg/kg'
  • palm/trunk/SOURCE/modules.f90

    r4495 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4495 2020-04-13 20:11:20Z raasch
    2730! +restart_data_format, restart_data_format_input|output, include_total_domain_boundaries
    2831!
     
    229232    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_e              !< artificial numerical dissipation flux at south face of grid box - subgrid-scale TKE
    230233    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nc             !< artificial numerical dissipation flux at south face of grid box - clouddrop-number concentration
     234    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_ni             !< artificial numerical dissipation flux at south face of grid box - ice crystal-number concentration
    231235    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nr             !< artificial numerical dissipation flux at south face of grid box - raindrop-number concentration
    232236    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_pt             !< artificial numerical dissipation flux at south face of grid box - potential temperature
    233237    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_q              !< artificial numerical dissipation flux at south face of grid box - mixing ratio
    234238    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qc             !< artificial numerical dissipation flux at south face of grid box - cloudwater
     239    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qi             !< artificial numerical dissipation flux at south face of grid box - ice crystal mixing ratio
    235240    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qr             !< artificial numerical dissipation flux at south face of grid box - rainwater
    236241    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_s              !< artificial numerical dissipation flux at south face of grid box - passive scalar
     
    244249    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_e              !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
    245250    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nc             !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
     251    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_ni             !< 6th-order advective flux at south face of grid box - icecrystal-number concentration
    246252    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nr             !< 6th-order advective flux at south face of grid box - raindrop-number concentration
    247253    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_pt             !< 6th-order advective flux at south face of grid box - potential temperature
    248254    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_q              !< 6th-order advective flux at south face of grid box - mixing ratio
    249255    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qc             !< 6th-order advective flux at south face of grid box - cloudwater
     256    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qi             !< 6th-order advective flux at south face of grid box - ice crystal
    250257    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qr             !< 6th-order advective flux at south face of grid box - rainwater
    251258    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_s              !< 6th-order advective flux at south face of grid box - passive scalar
     
    271278    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_e    !< artificial numerical dissipation flux at left face of grid box - subgrid-scale TKE
    272279    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nc   !< artificial numerical dissipation flux at left face of grid box - clouddrop-number concentration
     280    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_ni   !< artificial numerical dissipation flux at left face of grid box - ice crystal-number concentration
    273281    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nr   !< artificial numerical dissipation flux at left face of grid box - raindrop-number concentration
    274282    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_pt   !< artificial numerical dissipation flux at left face of grid box - potential temperature
    275283    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_q    !< artificial numerical dissipation flux at left face of grid box - mixing ratio
    276284    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qc   !< artificial numerical dissipation flux at left face of grid box - cloudwater
     285    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qi   !< artificial numerical dissipation flux at left face of grid box - ice crystal
    277286    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qr   !< artificial numerical dissipation flux at left face of grid box - rainwater
    278287    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_s    !< artificial numerical dissipation flux at left face of grid box - passive scalar
     
    284293    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_e    !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
    285294    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nc   !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
     295    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_ni   !< 6th-order advective flux at south face of grid box - ice crystal-number concentration
    286296    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nr   !< 6th-order advective flux at south face of grid box - raindrop-number concentration
    287297    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_pt   !< 6th-order advective flux at south face of grid box - potential temperature
    288298    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_q    !< 6th-order advective flux at south face of grid box - mixing ratio
    289299    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qc   !< 6th-order advective flux at south face of grid box - cloudwater
     300    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qi   !< 6th-order advective flux at south face of grid box - ice crystal
    290301    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qr   !< 6th-order advective flux at south face of grid box - rainwater
    291302    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_s    !< 6th-order advective flux at south face of grid box - passive scalar
     
    324335    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_2    !< pointer for swapping of timelevels for respective quantity
    325336    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_3    !< pointer for swapping of timelevels for respective quantity
     337    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_1    !< pointer for swapping of timelevels for respective quantity
     338    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_2    !< pointer for swapping of timelevels for respective quantity
     339    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_3    !< pointer for swapping of timelevels for respective quantity
    326340    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_1    !< pointer for swapping of timelevels for respective quantity
    327341    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_2    !< pointer for swapping of timelevels for respective quantity
     
    336350    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_2    !< pointer for swapping of timelevels for respective quantity
    337351    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_3    !< pointer for swapping of timelevels for respective quantity
     352    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_1    !< pointer for swapping of timelevels for respective quantity
     353    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_2    !< pointer for swapping of timelevels for respective quantity
     354    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_3    !< pointer for swapping of timelevels for respective quantity
    338355    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_v    !< pointer: volume of liquid water
    339356    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_vp   !< pointer: liquid water weighting factor
     
    367384    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nc         !< pointer: cloud drop number density
    368385    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nc_p       !< pointer: prognostic value of cloud drop number density
     386    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ni         !< pointer: ice crystal number density
     387    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ni_p       !< pointer: prognostic value of ice crystal number density
    369388    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nr         !< pointer: rain drop number density
    370389    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  nr_p       !< pointer: prognostic value of rain drop number density
     
    375394    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  q_p        !< pointer: prognostic value of mixing ratio
    376395    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qc         !< pointer: cloud water content
    377     REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qc_p       !< pointer: cloud water content
     396    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qc_p       !< pointer: prognostic value cloud water content
     397    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qi         !< pointer: ice crystal content
     398    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  qi_p       !< pointer: prognostic value ice crystal content
    378399    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ql         !< pointer: liquid water content
    379400    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ql_c       !< pointer: change in liquid water content due to
     
    389410    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  te_m       !< pointer: weighted tendency of e for previous sub-timestep (Runge-Kutta)
    390411    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tnc_m      !< pointer: weighted tendency of nc for previous sub-timestep (Runge-Kutta)
     412    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tni_m      !< pointer: weighted tendency of ni for previous sub-timestep (Runge-Kutta)
    391413    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tnr_m      !< pointer: weighted tendency of nr for previous sub-timestep (Runge-Kutta)
    392414    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tpt_m      !< pointer: weighted tendency of pt for previous sub-timestep (Runge-Kutta)
    393415    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tq_m       !< pointer: weighted tendency of q for previous sub-timestep (Runge-Kutta)
    394416    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tqc_m      !< pointer: weighted tendency of qc for previous sub-timestep (Runge-Kutta)
     417    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tqi_m      !< pointer: weighted tendency of qi for previous sub-timestep (Runge-Kutta)
    395418    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  tqr_m      !< pointer: weighted tendency of qr for previous sub-timestep (Runge-Kutta)
    396419    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  ts_m       !< pointer: weighted tendency of s for previous sub-timestep (Runge-Kutta)
     
    462485    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  lpt_av        !< avg. liquid water potential temperature
    463486    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_av         !< avg. cloud drop number density
     487    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ni_av         !< avg. ice crystal number density
    464488    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_av         !< avg. rain drop number density
    465489    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  p_av          !< avg. perturbation pressure
     
    471495                                                                      !< (or total water content with active cloud physics)
    472496    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_av         !< avg. cloud water content
     497    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qi_av         !< avg. ice crystal content
    473498    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_av         !< avg. liquid water content
    474499    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_c_av       !< avg. change in liquid water content due to
     
    14231448    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_ws2_ws_l    !< subdomain sum of vertical momentum flux w'w' (5th-order advection scheme only)
    14241449    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsncs_ws_l  !< subdomain sum of vertical clouddrop-number concentration flux w'nc' (5th-order advection scheme only)
     1450    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsnis_ws_l  !< subdomain sum of vertical ice crystal concentration flux w'ni' (5th-order advection scheme only)
    14251451    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsnrs_ws_l  !< subdomain sum of vertical raindrop-number concentration flux w'nr' (5th-order advection scheme only)
    14261452    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wspts_ws_l  !< subdomain sum of vertical sensible heat flux w'pt' (5th-order advection scheme only)
    14271453    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqs_ws_l   !< subdomain sum of vertical latent heat flux w'q' (5th-order advection scheme only)
    14281454    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqcs_ws_l  !< subdomain sum of vertical cloudwater flux w'qc' (5th-order advection scheme only)
     1455    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqis_ws_l  !< subdomain sum of vertical ice crystal flux w'qi' (5th-order advection scheme only)
    14291456    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqrs_ws_l  !< subdomain sum of vertical rainwater flux w'qr' (5th-order advection scheme only)
    14301457    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

    r4455 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4455 2020-03-11 12:20:29Z Giersch
    2730! Axis attribute added to netcdf output
    2831!
     
    880883                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
    881884                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
    882                        's', 'theta', 'thetal', 'thetav' )
     885                       's', 'theta', 'thetal', 'thetav', 'qi', 'ni' )
    883886
    884887                   grid_x = 'x'
     
    16441647                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
    16451648                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
    1646                        's', 'theta', 'thetal', 'thetav' )
     1649                       's', 'theta', 'thetal', 'thetav', 'qi', 'ni' )
    16471650
    16481651                   grid_x = 'x'
     
    26592662!
    26602663!--                   Most variables are defined on the zu grid
    2661                       CASE ( 'e_xy', 'nc_xy', 'nr_xy', 'p_xy',                 &
     2664                      CASE ( 'e_xy', 'nc_xy', 'ni_xy', 'nr_xy', 'p_xy',        &
    26622665                             'pc_xy', 'pr_xy', 'prr_xy', 'q_xy',               &
    2663                              'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
     2666                             'qc_xy', 'qi_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',  &
    26642667                             'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
    26652668                             's_xy',                                           &
     
    35833586!
    35843587!--                Most variables are defined on the zu grid
    3585                    CASE ( 'e_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz',           &
    3586                           'pr_xz', 'prr_xz', 'q_xz', 'qc_xz',                  &
     3588                   CASE ( 'e_xz', 'nc_xz', 'ni_xz', 'nr_xz', 'p_xz', 'pc_xz',  &
     3589                          'pr_xz', 'prr_xz', 'q_xz', 'qc_xz', 'qi_xz',         &
    35873590                          'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
    35883591                          'qv_xz', 's_xz',                                     &
     
    44604463!
    44614464!--                Most variables are defined on the zu grid
    4462                    CASE ( 'e_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz',           &
    4463                           'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'ql_yz',          &
     4465                   CASE ( 'e_yz', 'nc_yz', 'ni_yz', 'nr_yz', 'p_yz', 'pc_yz',  &
     4466                          'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'qi_yz', 'ql_yz', &
    44644467                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
    44654468                          's_yz',                                              &
     
    58535856!
    58545857!--             Most variables are defined on the zu levels
    5855                 CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
     5858                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'ni', 'qi',    &
    58565859                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
    58575860                       'rho_sea_water', 's', 'sa', &
  • palm/trunk/SOURCE/surface_data_output_mod.f90

    r4500 r4502  
    2525! -----------------
    2626! $Id$
     27! Implementation of ice microphysics
     28!
     29! 4500 2020-04-17 10:12:45Z suehring
    2730! - Correct output of ground/wall heat flux at USM surfaces
    2831! - Add conversion factor to heat and momentum-flux output
     
    16951698               ENDIF
    16961699
     1700            CASE ( 'qis' )
     1701!
     1702!--            Output of instantaneous data
     1703               IF ( av == 0 )  THEN
     1704                  CALL surface_data_output_collect( surf_def_h(0)%qis,         &
     1705                                               surf_def_h(1)%qis,              &
     1706                                               surf_lsm_h%qis,                 &
     1707                                               surf_usm_h%qis,                 &
     1708                                               surf_def_v(0)%qis,              &
     1709                                               surf_lsm_v(0)%qis,              &
     1710                                               surf_usm_v(0)%qis,              &
     1711                                               surf_def_v(1)%qis,              &
     1712                                               surf_lsm_v(1)%qis,              &
     1713                                               surf_usm_v(1)%qis,              &
     1714                                               surf_def_v(2)%qis,              &
     1715                                               surf_lsm_v(2)%qis,              &
     1716                                               surf_usm_v(2)%qis,              &
     1717                                               surf_def_v(3)%qis,              &
     1718                                               surf_lsm_v(3)%qis,              &
     1719                                               surf_usm_v(3)%qis )
     1720               ELSE
     1721!
     1722!--               Output of averaged data
     1723                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
     1724                                        REAL( average_count_surf, KIND=wp )
     1725                  surfaces%var_av(:,n_out) = 0.0_wp
     1726
     1727               ENDIF
     1728
     1729            CASE ( 'nis' )
     1730!
     1731!--            Output of instantaneous data
     1732               IF ( av == 0 )  THEN
     1733                  CALL surface_data_output_collect( surf_def_h(0)%nis,         &
     1734                                               surf_def_h(1)%nis,              &
     1735                                               surf_lsm_h%nis,                 &
     1736                                               surf_usm_h%nis,                 &
     1737                                               surf_def_v(0)%nis,              &
     1738                                               surf_lsm_v(0)%nis,              &
     1739                                               surf_usm_v(0)%nis,              &
     1740                                               surf_def_v(1)%nis,              &
     1741                                               surf_lsm_v(1)%nis,              &
     1742                                               surf_usm_v(1)%nis,              &
     1743                                               surf_def_v(2)%nis,              &
     1744                                               surf_lsm_v(2)%nis,              &
     1745                                               surf_usm_v(2)%nis,              &
     1746                                               surf_def_v(3)%nis,              &
     1747                                               surf_lsm_v(3)%nis,              &
     1748                                               surf_usm_v(3)%nis )
     1749               ELSE
     1750!
     1751!--               Output of averaged data
     1752                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
     1753                                        REAL( average_count_surf, KIND=wp )
     1754                  surfaces%var_av(:,n_out) = 0.0_wp
     1755
     1756               ENDIF
     1757
    16971758            CASE ( 'qrs' )
    16981759!
     
    21532214                                               surf_lsm_v(3)%ncsws,            &
    21542215                                               surf_usm_v(3)%ncsws )
     2216               ELSE
     2217!
     2218!--               Output of averaged data
     2219                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
     2220                                        REAL( average_count_surf, KIND=wp )
     2221                  surfaces%var_av(:,n_out) = 0.0_wp
     2222
     2223               ENDIF
     2224
     2225
     2226            CASE ( 'qisws' )
     2227!
     2228!--            Output of instantaneous data
     2229               IF ( av == 0 )  THEN
     2230                  CALL surface_data_output_collect( surf_def_h(0)%qisws,       &
     2231                                               surf_def_h(1)%qisws,            &
     2232                                               surf_lsm_h%qisws,               &
     2233                                               surf_usm_h%qisws,               &
     2234                                               surf_def_v(0)%qisws,            &
     2235                                               surf_lsm_v(0)%qisws,            &
     2236                                               surf_usm_v(0)%qisws,            &
     2237                                               surf_def_v(1)%qisws,            &
     2238                                               surf_lsm_v(1)%qisws,            &
     2239                                               surf_usm_v(1)%qisws,            &
     2240                                               surf_def_v(2)%qisws,            &
     2241                                               surf_lsm_v(2)%qisws,            &
     2242                                               surf_usm_v(2)%qisws,            &
     2243                                               surf_def_v(3)%qisws,            &
     2244                                               surf_lsm_v(3)%qisws,            &
     2245                                               surf_usm_v(3)%qisws )
     2246               ELSE
     2247!
     2248!--               Output of averaged data
     2249                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
     2250                                        REAL( average_count_surf, KIND=wp )
     2251                  surfaces%var_av(:,n_out) = 0.0_wp
     2252
     2253               ENDIF
     2254
     2255            CASE ( 'nisws' )
     2256!
     2257!--            Output of instantaneous data
     2258               IF ( av == 0 )  THEN
     2259                  CALL surface_data_output_collect( surf_def_h(0)%nisws,       &
     2260                                               surf_def_h(1)%nisws,            &
     2261                                               surf_lsm_h%nisws,               &
     2262                                               surf_usm_h%nisws,               &
     2263                                               surf_def_v(0)%nisws,            &
     2264                                               surf_lsm_v(0)%nisws,            &
     2265                                               surf_usm_v(0)%nisws,            &
     2266                                               surf_def_v(1)%nisws,            &
     2267                                               surf_lsm_v(1)%nisws,            &
     2268                                               surf_usm_v(1)%nisws,            &
     2269                                               surf_def_v(2)%nisws,            &
     2270                                               surf_lsm_v(2)%nisws,            &
     2271                                               surf_usm_v(2)%nisws,            &
     2272                                               surf_def_v(3)%nisws,            &
     2273                                               surf_lsm_v(3)%nisws,            &
     2274                                               surf_usm_v(3)%nisws )
    21552275               ELSE
    21562276!
     
    31203240                                           surf_usm_v(3)%ncs, n_out )
    31213241
     3242            CASE ( 'qis' )
     3243               CALL surface_data_output_sum_up( surf_def_h(0)%qis,             &
     3244                                           surf_def_h(1)%qis,                  &
     3245                                           surf_lsm_h%qis,                     &
     3246                                           surf_usm_h%qis,                     &
     3247                                           surf_def_v(0)%qis,                  &
     3248                                           surf_lsm_v(0)%qis,                  &
     3249                                           surf_usm_v(0)%qis,                  &
     3250                                           surf_def_v(1)%qis,                  &
     3251                                           surf_lsm_v(1)%qis,                  &
     3252                                           surf_usm_v(1)%qis,                  &
     3253                                           surf_def_v(2)%qis,                  &
     3254                                           surf_lsm_v(2)%qis,                  &
     3255                                           surf_usm_v(2)%qis,                  &
     3256                                           surf_def_v(3)%qis,                  &
     3257                                           surf_lsm_v(3)%qis,                  &
     3258                                           surf_usm_v(3)%qrs, n_out )
     3259
     3260            CASE ( 'nis' )
     3261               CALL surface_data_output_sum_up( surf_def_h(0)%nis,             &
     3262                                           surf_def_h(1)%nis,                  &
     3263                                           surf_lsm_h%nis,                     &
     3264                                           surf_usm_h%nis,                     &
     3265                                           surf_def_v(0)%nis,                  &
     3266                                           surf_lsm_v(0)%nis,                  &
     3267                                           surf_usm_v(0)%nis,                  &
     3268                                           surf_def_v(1)%nis,                  &
     3269                                           surf_lsm_v(1)%nis,                  &
     3270                                           surf_usm_v(1)%nis,                  &
     3271                                           surf_def_v(2)%nis,                  &
     3272                                           surf_lsm_v(2)%nis,                  &
     3273                                           surf_usm_v(2)%nis,                  &
     3274                                           surf_def_v(3)%nis,                  &
     3275                                           surf_lsm_v(3)%nis,                  &
     3276                                           surf_usm_v(3)%nis, n_out )
     3277
    31223278            CASE ( 'qrs' )
    31233279               CALL surface_data_output_sum_up( surf_def_h(0)%qrs,             &
     
    34113567                                           surf_lsm_v(3)%ncsws,                &
    34123568                                           surf_usm_v(3)%ncsws, n_out )
     3569
     3570            CASE ( 'qisws' )
     3571               CALL surface_data_output_sum_up( surf_def_h(0)%qisws,           &
     3572                                           surf_def_h(1)%qisws,                &
     3573                                           surf_lsm_h%qisws,                   &
     3574                                           surf_usm_h%qisws,                   &
     3575                                           surf_def_v(0)%qisws,                &
     3576                                           surf_lsm_v(0)%qisws,                &
     3577                                           surf_usm_v(0)%qisws,                &
     3578                                           surf_def_v(1)%qisws,                &
     3579                                           surf_lsm_v(1)%qisws,                &
     3580                                           surf_usm_v(1)%qisws,                &
     3581                                           surf_def_v(2)%qisws,                &
     3582                                           surf_lsm_v(2)%qisws,                &
     3583                                           surf_usm_v(2)%qisws,                &
     3584                                           surf_def_v(3)%qisws,                &
     3585                                           surf_lsm_v(3)%qisws,                &
     3586                                           surf_usm_v(3)%qisws, n_out )
     3587
     3588            CASE ( 'nisws' )
     3589               CALL surface_data_output_sum_up( surf_def_h(0)%nisws,           &
     3590                                           surf_def_h(1)%nisws,                &
     3591                                           surf_lsm_h%nisws,                   &
     3592                                           surf_usm_h%nisws,                   &
     3593                                           surf_def_v(0)%nisws,                &
     3594                                           surf_lsm_v(0)%nisws,                &
     3595                                           surf_usm_v(0)%nisws,                &
     3596                                           surf_def_v(1)%nisws,                &
     3597                                           surf_lsm_v(1)%nisws,                &
     3598                                           surf_usm_v(1)%nisws,                &
     3599                                           surf_def_v(2)%nisws,                &
     3600                                           surf_lsm_v(2)%nisws,                &
     3601                                           surf_usm_v(2)%nisws,                &
     3602                                           surf_def_v(3)%nisws,                &
     3603                                           surf_lsm_v(3)%nisws,                &
     3604                                           surf_usm_v(3)%nisws, n_out )
    34133605
    34143606            CASE ( 'qrsws' )
     
    45334725               unit = 'm/s'
    45344726
    4535             CASE ( 'ss', 'qcs', 'ncs', 'qrs', 'nrs' )
     4727            CASE ( 'ss', 'qcs', 'ncs', 'qis', 'nis', 'qrs', 'nrs' )
    45364728               unit = '1'
    45374729
     
    45454737               unit = 'm2/s2'
    45464738
    4547             CASE ( 'qcsws', 'ncsws', 'qrsws', 'nrsws', 'sasws' )
     4739            CASE ( 'qcsws', 'ncsws', 'qisws', 'nisws', 'qrsws', 'nrsws', 'sasws' )
    45484740
    45494741            CASE ( 'shf' )
  • palm/trunk/SOURCE/surface_mod.f90

    r4495 r4502  
    2626! -----------------
    2727! $Id$
     28! Implementation of ice microphysics
     29!
     30! 4495 2020-04-13 20:11:20Z raasch
    2831! restart data handling with MPI-IO added
    2932!
     
    148151
    149152!
    150 !-- Data type used to identify grid-points where horizontal boundary conditions 
    151 !-- are applied 
     153!-- Data type used to identify grid-points where horizontal boundary conditions
     154!-- are applied
    152155    TYPE bc_type
    153156       INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element
     
    156159       INTEGER(iwp) :: ns   !< number of surface elements on the PE
    157160
    158        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i !< x-index linking to the PALM 3D-grid 
    159        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j !< y-index linking to the PALM 3D-grid   
    160        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k !< z-index linking to the PALM 3D-grid   
    161 
    162        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i) 
    163        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< end index within surface data type for given (j,i) 
     161       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i !< x-index linking to the PALM 3D-grid
     162       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j !< y-index linking to the PALM 3D-grid
     163       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k !< z-index linking to the PALM 3D-grid
     164
     165       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
     166       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< end index within surface data type for given (j,i)
    164167
    165168    END TYPE bc_type
    166169!
    167 !-- Data type used to identify and treat surface-bounded grid points 
     170!-- Data type used to identify and treat surface-bounded grid points
    168171    TYPE surf_type
    169172
    170173       LOGICAL ::  albedo_from_ascii = .FALSE. !< flag indicating that albedo for urban surfaces is input via ASCII format (just for a workaround)
    171    
     174
    172175       INTEGER(iwp) :: ioff                                !< offset value in x-direction, used to determine index of surface element
    173176       INTEGER(iwp) :: joff                                !< offset value in y-direction, used to determine index of surface element
     
    175178       INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
    176179
    177        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid 
    178        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid   
    179        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid       
    180 
    181        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  facing  !< Bit indicating surface orientation 
    182      
    183        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i) 
    184        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< End index within surface data type for given (j,i) 
    185 
    186        REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_mo      !< surface-layer height 
     180       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid
     181       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid
     182       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid
     183
     184       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  facing  !< Bit indicating surface orientation
     185
     186       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< Start index within surface data type for given (j,i)
     187       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index   !< End index within surface data type for given (j,i)
     188
     189       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_mo      !< surface-layer height
    187190       REAL(wp), DIMENSION(:), ALLOCATABLE ::  uvw_abs   !< absolute surface-parallel velocity
    188191       REAL(wp), DIMENSION(:), ALLOCATABLE ::  us        !< friction velocity
     
    192195       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcs       !< scaling parameter qc
    193196       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncs       !< scaling parameter nc
     197       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qis       !< scaling parameter qi
     198       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nis       !< scaling parameter ni
    194199       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrs       !< scaling parameter qr
    195200       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrs       !< scaling parameter nr
     
    205210       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qv1       !< mixing ratio at first grid level
    206211       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vpt1      !< virtual potential temperature at first grid level
    207        
     212
    208213       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  css     !< scaling parameter chemical species
    209214!
    210215!--    Define arrays for surface fluxes
    211216       REAL(wp), DIMENSION(:), ALLOCATABLE ::  usws      !< vertical momentum flux for u-component at horizontal surfaces
    212        REAL(wp), DIMENSION(:), ALLOCATABLE ::  vsws      !< vertical momentum flux for v-component at horizontal surfaces 
     217       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vsws      !< vertical momentum flux for v-component at horizontal surfaces
    213218
    214219       REAL(wp), DIMENSION(:), ALLOCATABLE ::  shf       !< surface flux sensible heat
    215        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws      !< surface flux latent heat 
     220       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws      !< surface flux latent heat
    216221       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ssws      !< surface flux passive scalar
    217222       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcsws     !< surface flux qc
    218223       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncsws     !< surface flux nc
     224       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qisws     !< surface flux qi
     225       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nisws     !< surface flux ni
    219226       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrsws     !< surface flux qr
    220227       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrsws     !< surface flux nr
     
    224231       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  amsws   !< surface flux aerosol mass:   dim 1: flux, dim 2: bin
    225232       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  gtsws   !< surface flux gesous tracers: dim 1: flux, dim 2: gas
    226        
     233
    227234       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  cssws   !< surface flux chemical species
    228235!
     
    240247       CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  vegetation_type_name  !< water type at name surface element
    241248       CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE ::  water_type_name       !< water type at name surface element
    242        
     249
    243250       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nzt_pavement     !< top index for pavement in soil
    244251       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  building_type    !< building type at surface element
     
    246253       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  vegetation_type  !< vegetation type at surface element
    247254       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  water_type       !< water type at surface element
    248        
     255
    249256       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  albedo_type   !< albedo type, for each fraction (wall,green,window or vegetation,pavement water)
    250257
     
    254261       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  water_surface       !< flag parameter for water surfaces
    255262       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  vegetation_surface  !< flag parameter for natural land surfaces
    256        
     263
    257264       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  albedo            !< broadband albedo for each surface fraction (LSM: vegetation, water, pavement; USM: wall, green, window)
    258265       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  emissivity        !< emissivity of the surface, for each fraction  (LSM: vegetation, water, pavement; USM: wall, green, window)
     
    271278       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  pt_surface        !< skin-surface temperature
    272279       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  vpt_surface       !< skin-surface virtual temperature
    273        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net           !< net radiation 
     280       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net           !< net radiation
    274281       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  rad_net_l         !< net radiation, used in USM
    275        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil/ wall (W/m/K) 
    276        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_green    !< heat conductivity of green soil (W/m/K) 
    277        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_window   !< heat conductivity of windows (W/m/K) 
    278        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)   
     282       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h          !< heat conductivity of soil/ wall (W/m/K)
     283       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_green    !< heat conductivity of green soil (W/m/K)
     284       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_window   !< heat conductivity of windows (W/m/K)
     285       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_h_def      !< default heat conductivity of soil (W/m/K)
    279286
    280287       REAL(wp), DIMENSION(:), ALLOCATABLE ::  rad_lw_in           !< incoming longwave radiation
     
    302309       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws_veg            !< surface flux of latent heat (vegetation portion)
    303310
    304        REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a                 !< aerodynamic resistance 
     311       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a                 !< aerodynamic resistance
    305312       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a_green           !< aerodynamic resistance at green fraction
    306313       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_a_window          !< aerodynamic resistance at window fraction
     
    310317       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_s                 !< total surface resistance (combination of r_soil and r_canopy)
    311318       REAL(wp), DIMENSION(:), ALLOCATABLE ::  r_canopy_min        !< minimum canopy (stomatal) resistance
    312        
     319
    313320       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm             !< near surface air potential temperature at distance 10 cm from the surface (K)
    314        
     321
    315322       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  alpha_vg          !< coef. of Van Genuchten
    316323       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lambda_w          !< hydraulic diffusivity of soil (?)
     
    322329       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_sat             !< saturation soil moisture (m3/m3)
    323330       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  m_wilt            !< soil moisture at permanent wilting point (m3/m3)
    324        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_vg              !< coef. Van Genuchten 
     331       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  n_vg              !< coef. Van Genuchten
    325332       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total_def   !< default volumetric heat capacity of the (soil) layer (J/m3/K)
    326333       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total       !< volumetric heat capacity of the actual soil matrix (J/m3/K)
    327334       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  root_fr           !< root fraction within the soil layers
    328        
     335
    329336!--    Indoor model variables
    330        REAL(wp), DIMENSION(:), ALLOCATABLE ::  waste_heat          !< waste heat 
     337       REAL(wp), DIMENSION(:), ALLOCATABLE ::  waste_heat          !< waste heat
    331338!
    332339!--    Urban surface variables
     
    334341
    335342       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  isroof_surf          !< flag indicating roof surfaces
    336        LOGICAL, DIMENSION(:), ALLOCATABLE  ::  ground_level         !< flag indicating ground floor level surfaces 
     343       LOGICAL, DIMENSION(:), ALLOCATABLE  ::  ground_level         !< flag indicating ground floor level surfaces
    337344
    338345       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_summer  !< indoor target temperature summer
    339        REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_winter  !< indoor target temperature summer       
     346       REAL(wp), DIMENSION(:), ALLOCATABLE ::  target_temp_winter  !< indoor target temperature summer
    340347
    341348       REAL(wp), DIMENSION(:), ALLOCATABLE ::  c_surface           !< heat capacity of the wall surface skin (J/m2/K)
     
    373380       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinlw            !< longwave radiation falling to local surface including radiation from reflections
    374381       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfoutlw           !< total longwave radiation outgoing from nonvirtual surfaces surfaces after all reflection
    375        
     382
    376383       REAL(wp), DIMENSION(:), ALLOCATABLE ::  n_vg_green          !< vangenuchten parameters
    377384       REAL(wp), DIMENSION(:), ALLOCATABLE ::  alpha_vg_green      !< vangenuchten parameters
     
    384391       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_wall_stag      !< wall grid spacing (edge-edge)
    385392       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_wall_stag     !< 1/dz_wall_stag
    386        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_wall_m         !< t_wall prognostic array 
     393       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_wall_m         !< t_wall prognostic array
    387394       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw                !< wall layer depths (m)
    388395       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_window      !< volumetric heat capacity of the window material ( J m-3 K-1 ) (= 2.19E6)
     
    391398       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_window_stag    !< window grid spacing (edge-edge)
    392399       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_window_stag   !< 1/dz_window_stag
    393        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_window_m       !< t_window prognostic array 
     400       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_window_m       !< t_window prognostic array
    394401       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_window         !< window layer depths (m)
    395402       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_green       !< volumetric heat capacity of the green material ( J m-3 K-1 ) (= 2.19E6)
     
    399406       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dz_green_stag     !< green grid spacing (edge-edge)
    400407       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddz_green_stag    !< 1/dz_green_stag
    401        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_green_m        !< t_green prognostic array 
     408       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tt_green_m        !< t_green prognostic array
    402409       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  zw_green          !< green layer depths (m)
    403410
     
    421428       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
    422429       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
    423        REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface 
     430       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface
    424431       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_av       !< average of wghf_eb
    425432       REAL(wp), DIMENSION(:), ALLOCATABLE ::  wghf_eb_window_av  !< average of wghf_eb window
     
    437444
    438445       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_10cm_av       !< average of theta_10cm (K)
    439        
     446
    440447       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_wall_av      !< Average of t_wall
    441448       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  t_window_av    !< Average of t_window
     
    455462    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_usm_v  !< vertical urban surfaces (North, South, East, West)
    456463
    457     INTEGER(iwp), PARAMETER ::  ind_veg_wall  = 0            !< index for vegetation / wall-surface fraction, used for access of albedo, emissivity, etc., for each surface type   
     464    INTEGER(iwp), PARAMETER ::  ind_veg_wall  = 0            !< index for vegetation / wall-surface fraction, used for access of albedo, emissivity, etc., for each surface type
    458465    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
    459466    INTEGER(iwp), PARAMETER ::  ind_wat_win   = 2            !< index for water / window-surface fraction, used for access of albedo, emissivity, etc., for each surface type
    460467
    461     INTEGER(iwp) ::  ns_h_on_file(0:2)                       !< total number of horizontal surfaces with the same facing, required for writing restart data 
    462     INTEGER(iwp) ::  ns_v_on_file(0:3)                       !< total number of vertical surfaces with the same facing, required for writing restart data 
    463 
    464     LOGICAL ::  vertical_surfaces_exist = .FALSE.   !< flag indicating that there are vertical urban/land surfaces 
     468    INTEGER(iwp) ::  ns_h_on_file(0:2)                       !< total number of horizontal surfaces with the same facing, required for writing restart data
     469    INTEGER(iwp) ::  ns_v_on_file(0:3)                       !< total number of vertical surfaces with the same facing, required for writing restart data
     470
     471    LOGICAL ::  vertical_surfaces_exist = .FALSE.   !< flag indicating that there are vertical urban/land surfaces
    465472                                                    !< in the domain (required to activiate RTM)
    466473
     
    468475    LOGICAL ::  surf_microphysics_morrison = .FALSE.   !< use 2-moment Morrison (add. prog. eq. for nc and qc)
    469476    LOGICAL ::  surf_microphysics_seifert = .FALSE.    !< use 2-moment Seifert and Beheng scheme
     477    LOGICAL ::  surf_microphysics_ice_extension = .FALSE. !< use 2-moment Seifert and Beheng scheme
    470478
    471479
     
    473481
    474482    PRIVATE
    475    
     483
    476484    INTERFACE init_bc
    477485       MODULE PROCEDURE init_bc
     
    481489       MODULE PROCEDURE init_single_surface_properties
    482490    END INTERFACE init_single_surface_properties
    483    
     491
    484492    INTERFACE init_surfaces
    485493       MODULE PROCEDURE init_surfaces
     
    522530           surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,      &
    523531           vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison,             &
    524            surf_microphysics_seifert
     532           surf_microphysics_seifert, surf_microphysics_ice_extension
    525533!
    526534!-- Public subroutines and functions
     
    544552! Description:
    545553! ------------
    546 !> Initialize data type for setting boundary conditions at horizontal and 
    547 !> vertical surfaces. 
     554!> Initialize data type for setting boundary conditions at horizontal and
     555!> vertical surfaces.
    548556!------------------------------------------------------------------------------!
    549557    SUBROUTINE init_bc
     
    559567       INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of horizontal surfaces at (j,i)-grid point
    560568       INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index of horizontal surface elements
    561        
     569
    562570       INTEGER(iwp), DIMENSION(0:3) ::  num_v         !< number of vertical surfaces on subdomain
    563571       INTEGER(iwp), DIMENSION(0:3) ::  num_v_kji     !< number of vertical surfaces at (j,i)-grid point
    564572       INTEGER(iwp), DIMENSION(0:3) ::  start_index_v !< local start index of vertical surface elements
    565573!
    566 !--    Set offset indices, i.e. index difference between surface element and 
     574!--    Set offset indices, i.e. index difference between surface element and
    567575!--    surface-bounded grid point.
    568576!--    Horizontal surfaces - no horizontal offsets
     
    596604!
    597605!--    Initialize data structure for horizontal surfaces, i.e. count the number
    598 !--    of surface elements, allocate and initialize the respective index arrays, 
    599 !--    and set the respective start and end indices at each (j,i)-location. 
    600 !--    The index space is defined also over the ghost points, so that e.g. 
    601 !--    boundary conditions for diagnostic quanitities can be set on ghost 
    602 !--    points so that no exchange is required any more. 
     606!--    of surface elements, allocate and initialize the respective index arrays,
     607!--    and set the respective start and end indices at each (j,i)-location.
     608!--    The index space is defined also over the ghost points, so that e.g.
     609!--    boundary conditions for diagnostic quanitities can be set on ghost
     610!--    points so that no exchange is required any more.
    603611       DO  l = 0, 1
    604612!
     
    608616             DO  j = nysg, nyng
    609617                DO  k = nzb+1, nzt
    610 !         
     618!
    611619!--                Check if current gridpoint belongs to the atmosphere
    612620                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
     
    619627             ENDDO
    620628          ENDDO
    621 !         
     629!
    622630!--       Save the number of horizontal surface elements
    623631          bc_h(l)%ns = num_h(l)
     
    631639          bc_h(l)%start_index = 1
    632640          bc_h(l)%end_index   = 0
    633          
     641
    634642          num_h(l)         = 1
    635643          start_index_h(l) = 1
    636644          DO  i = nxlg, nxrg
    637645             DO  j = nysg, nyng
    638          
     646
    639647                num_h_kji(l) = 0
    640648                DO  k = nzb+1, nzt
    641 !         
     649!
    642650!--                Check if current gridpoint belongs to the atmosphere
    643651                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
    644 !         
     652!
    645653!--                   Upward-facing
    646654                      IF ( .NOT. BTEST( wall_flags_total_0(k+bc_h(l)%koff,     &
     
    666674!
    667675!--    Initialize data structure for vertical surfaces, i.e. count the number
    668 !--    of surface elements, allocate and initialize the respective index arrays, 
     676!--    of surface elements, allocate and initialize the respective index arrays,
    669677!--    and set the respective start and end indices at each (j,i)-location.
    670678       DO  l = 0, 3
     
    675683             DO  j = nys, nyn
    676684                DO  k = nzb+1, nzt
    677 !         
     685!
    678686!--                Check if current gridpoint belongs to the atmosphere
    679687                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
     
    686694             ENDDO
    687695          ENDDO
    688 !         
     696!
    689697!--       Save the number of horizontal surface elements
    690698          bc_v(l)%ns = num_v(l)
    691699!
    692700!--       ALLOCATE arrays for horizontal surfaces. In contrast to the
    693 !--       horizontal surfaces, the index space is not defined over the 
    694 !--       ghost points. 
     701!--       horizontal surfaces, the index space is not defined over the
     702!--       ghost points.
    695703          ALLOCATE( bc_v(l)%i(1:bc_v(l)%ns) )
    696704          ALLOCATE( bc_v(l)%j(1:bc_v(l)%ns) )
     
    700708          bc_v(l)%start_index = 1
    701709          bc_v(l)%end_index   = 0
    702          
     710
    703711          num_v(l)         = 1
    704712          start_index_v(l) = 1
    705713          DO  i = nxl, nxr
    706714             DO  j = nys, nyn
    707          
     715
    708716                num_v_kji(l) = 0
    709717                DO  k = nzb+1, nzt
    710 !         
     718!
    711719!--                Check if current gridpoint belongs to the atmosphere
    712720                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
    713 !         
     721!
    714722!--                   Upward-facing
    715723                      IF ( .NOT. BTEST( wall_flags_total_0(k+bc_v(l)%koff,     &
     
    741749! ------------
    742750!> Initialize horizontal and vertical surfaces. Counts the number of default-,
    743 !> natural and urban surfaces and allocates memory, respectively. 
     751!> natural and urban surfaces and allocates memory, respectively.
    744752!------------------------------------------------------------------------------!
    745753    SUBROUTINE init_surface_arrays
     
    751759       IMPLICIT NONE
    752760
    753        INTEGER(iwp)                 ::  i         !< running index x-direction 
     761       INTEGER(iwp)                 ::  i         !< running index x-direction
    754762       INTEGER(iwp)                 ::  j         !< running index y-direction
    755763       INTEGER(iwp)                 ::  k         !< running index z-direction
    756764       INTEGER(iwp)                 ::  l         !< index variable for surface facing
    757        INTEGER(iwp)                 ::  num_lsm_h !< number of horizontally-aligned natural surfaces 
    758        INTEGER(iwp)                 ::  num_usm_h !< number of horizontally-aligned urban surfaces 
    759 
    760        INTEGER(iwp), DIMENSION(0:2) ::  num_def_h !< number of horizontally-aligned default surfaces 
    761        INTEGER(iwp), DIMENSION(0:3) ::  num_def_v !< number of vertically-aligned default surfaces 
    762        INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces 
    763        INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces 
     765       INTEGER(iwp)                 ::  num_lsm_h !< number of horizontally-aligned natural surfaces
     766       INTEGER(iwp)                 ::  num_usm_h !< number of horizontally-aligned urban surfaces
     767
     768       INTEGER(iwp), DIMENSION(0:2) ::  num_def_h !< number of horizontally-aligned default surfaces
     769       INTEGER(iwp), DIMENSION(0:3) ::  num_def_v !< number of vertically-aligned default surfaces
     770       INTEGER(iwp), DIMENSION(0:3) ::  num_lsm_v !< number of vertically-aligned natural surfaces
     771       INTEGER(iwp), DIMENSION(0:3) ::  num_usm_v !< number of vertically-aligned urban surfaces
    764772
    765773       INTEGER(iwp)              ::  num_surf_v_l !< number of vertically-aligned local urban/land surfaces
     
    768776       LOGICAL ::  building                       !< flag indicating building grid point
    769777       LOGICAL ::  terrain                        !< flag indicating natural terrain grid point
    770        LOGICAL ::  unresolved_building            !< flag indicating a grid point where actually a building is 
    771                                                   !< defined but not resolved by the vertical grid 
     778       LOGICAL ::  unresolved_building            !< flag indicating a grid point where actually a building is
     779                                                  !< defined but not resolved by the vertical grid
    772780
    773781       num_def_h = 0
     
    780788!--    Surfaces are classified according to the input data read from static
    781789!--    input file. If no input file is present, all surfaces are classified
    782 !--    either as natural, urban, or default, depending on the setting of 
     790!--    either as natural, urban, or default, depending on the setting of
    783791!--    land_surface and urban_surface. To control this, use the control
    784792!--    flag topo_no_distinct
    785793!
    786 !--    Count number of horizontal surfaces on local domain 
     794!--    Count number of horizontal surfaces on local domain
    787795       DO  i = nxl, nxr
    788796          DO  j = nys, nyn
     
    798806!
    799807!--                   Determine flags indicating a terrain surface, a building
    800 !--                   surface, 
     808!--                   surface,
    801809                      terrain  = BTEST( wall_flags_total_0(k-1,j,i), 5 )  .OR.       &
    802810                                 topo_no_distinct
     
    804812                                 topo_no_distinct
    805813!
    806 !--                   unresolved_building indicates a surface with equal height 
     814!--                   unresolved_building indicates a surface with equal height
    807815!--                   as terrain but with a non-grid resolved building on top.
    808816!--                   These surfaces will be flagged as urban surfaces.
     
    813821                      IF ( land_surface  .AND.  terrain  .AND.                 &
    814822                           .NOT. unresolved_building )  THEN
    815                          num_lsm_h    = num_lsm_h    + 1 
     823                         num_lsm_h    = num_lsm_h    + 1
    816824!
    817825!--                   Urban surface tpye
    818826                      ELSEIF ( urban_surface  .AND.  building )  THEN
    819                          num_usm_h    = num_usm_h    + 1 
     827                         num_usm_h    = num_usm_h    + 1
    820828!
    821829!--                   Default-surface type
    822830                      ELSEIF ( .NOT. land_surface    .AND.                     &
    823831                               .NOT. urban_surface )  THEN
    824                                
     832
    825833                         num_def_h(0) = num_def_h(0) + 1
    826834!
    827835!--                   Unclassifified surface-grid point. Give error message.
    828                       ELSE 
     836                      ELSE
    829837                         WRITE( message_string, * )                           &
    830838                                          'Unclassified upward-facing ' //    &
     
    840848                      num_def_h(2) = num_def_h(2) + 1
    841849!
    842 !--                Check for any other downward-facing surface. So far only for 
     850!--                Check for any other downward-facing surface. So far only for
    843851!--                default surface type.
    844852                   ELSEIF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
    845853                      num_def_h(1) = num_def_h(1) + 1
    846                    ENDIF 
     854                   ENDIF
    847855
    848856                ENDIF
     
    851859       ENDDO
    852860!
    853 !--    Count number of vertical surfaces on local domain 
     861!--    Count number of vertical surfaces on local domain
    854862       DO  i = nxl, nxr
    855863          DO  j = nys, nyn
     
    869877                      unresolved_building = BTEST( wall_flags_total_0(k,j-1,i), 5 )  &
    870878                                     .AND.  BTEST( wall_flags_total_0(k,j-1,i), 6 )
    871                                      
     879
    872880                      IF (  land_surface  .AND.  terrain  .AND.                &
    873881                           .NOT. unresolved_building )  THEN
    874                          num_lsm_v(0) = num_lsm_v(0) + 1 
     882                         num_lsm_v(0) = num_lsm_v(0) + 1
    875883                      ELSEIF ( urban_surface  .AND.  building )  THEN
    876                          num_usm_v(0) = num_usm_v(0) + 1 
     884                         num_usm_v(0) = num_usm_v(0) + 1
    877885!
    878886!--                   Default-surface type
    879887                      ELSEIF ( .NOT. land_surface    .AND.                     &
    880888                               .NOT. urban_surface )  THEN
    881                          num_def_v(0) = num_def_v(0) + 1 
     889                         num_def_v(0) = num_def_v(0) + 1
    882890!
    883891!--                   Unclassifified surface-grid point. Give error message.
    884                       ELSE 
     892                      ELSE
    885893                         WRITE( message_string, * )                            &
    886894                                          'Unclassified northward-facing ' //  &
     
    900908                      building = BTEST( wall_flags_total_0(k,j+1,i), 6 )  .OR.       &
    901909                                 topo_no_distinct
    902                                  
     910
    903911                      unresolved_building = BTEST( wall_flags_total_0(k,j+1,i), 5 )  &
    904                                      .AND.  BTEST( wall_flags_total_0(k,j+1,i), 6 ) 
    905                                
     912                                     .AND.  BTEST( wall_flags_total_0(k,j+1,i), 6 )
     913
    906914                      IF (  land_surface  .AND.  terrain  .AND.                &
    907915                           .NOT. unresolved_building )  THEN
    908                          num_lsm_v(1) = num_lsm_v(1) + 1 
     916                         num_lsm_v(1) = num_lsm_v(1) + 1
    909917                      ELSEIF ( urban_surface  .AND.  building )  THEN
    910                          num_usm_v(1) = num_usm_v(1) + 1 
     918                         num_usm_v(1) = num_usm_v(1) + 1
    911919!
    912920!--                   Default-surface type
    913921                      ELSEIF ( .NOT. land_surface    .AND.                     &
    914922                               .NOT. urban_surface )  THEN
    915                          num_def_v(1) = num_def_v(1) + 1 
     923                         num_def_v(1) = num_def_v(1) + 1
    916924!
    917925!--                   Unclassifified surface-grid point. Give error message.
    918                       ELSE 
     926                      ELSE
    919927                         WRITE( message_string, * )                            &
    920928                                          'Unclassified southward-facing ' //  &
     
    934942                      building = BTEST( wall_flags_total_0(k,j,i-1), 6 )  .OR.       &
    935943                                 topo_no_distinct
    936                                  
     944
    937945                      unresolved_building = BTEST( wall_flags_total_0(k,j,i-1), 5 )  &
    938946                                     .AND.  BTEST( wall_flags_total_0(k,j,i-1), 6 )
    939                                      
     947
    940948                      IF (  land_surface  .AND.  terrain  .AND.                &
    941949                           .NOT. unresolved_building )  THEN
    942                          num_lsm_v(2) = num_lsm_v(2) + 1 
     950                         num_lsm_v(2) = num_lsm_v(2) + 1
    943951                      ELSEIF ( urban_surface  .AND.  building )  THEN
    944                          num_usm_v(2) = num_usm_v(2) + 1 
     952                         num_usm_v(2) = num_usm_v(2) + 1
    945953!
    946954!--                   Default-surface type
    947955                      ELSEIF ( .NOT. land_surface    .AND.                     &
    948956                               .NOT. urban_surface )  THEN
    949                          num_def_v(2) = num_def_v(2) + 1 
     957                         num_def_v(2) = num_def_v(2) + 1
    950958!
    951959!--                   Unclassifified surface-grid point. Give error message.
    952                       ELSE 
     960                      ELSE
    953961                         WRITE( message_string, * )                            &
    954962                                          'Unclassified eastward-facing ' //   &
     
    968976                      building = BTEST( wall_flags_total_0(k,j,i+1), 6 )  .OR.       &
    969977                                 topo_no_distinct
    970                                  
     978
    971979                      unresolved_building = BTEST( wall_flags_total_0(k,j,i+1), 5 )  &
    972980                                     .AND.  BTEST( wall_flags_total_0(k,j,i+1), 6 )
    973                                  
     981
    974982                      IF (  land_surface  .AND.  terrain  .AND.                &
    975983                           .NOT. unresolved_building )  THEN
    976                          num_lsm_v(3) = num_lsm_v(3) + 1 
     984                         num_lsm_v(3) = num_lsm_v(3) + 1
    977985                      ELSEIF ( urban_surface  .AND.  building )  THEN
    978                          num_usm_v(3) = num_usm_v(3) + 1 
     986                         num_usm_v(3) = num_usm_v(3) + 1
    979987!
    980988!--                   Default-surface type
    981989                      ELSEIF ( .NOT. land_surface    .AND.                     &
    982990                               .NOT. urban_surface )  THEN
    983                          num_def_v(3) = num_def_v(3) + 1 
     991                         num_def_v(3) = num_def_v(3) + 1
    984992!
    985993!--                   Unclassifified surface-grid point. Give error message.
    986                       ELSE 
     994                      ELSE
    987995                         WRITE( message_string, * )                            &
    988996                                          'Unclassified westward-facing ' //   &
     
    10101018!
    10111019!--    Horizontal surface, natural type, so far only upward-facing
    1012        surf_lsm_h%ns    = num_lsm_h 
     1020       surf_lsm_h%ns    = num_lsm_h
    10131021!
    10141022!--    Horizontal surface, urban type, so far only upward-facing
    1015        surf_usm_h%ns    = num_usm_h   
     1023       surf_usm_h%ns    = num_usm_h
    10161024!
    10171025!--    Vertical surface, default type, northward facing
     
    10511059       surf_usm_v(3)%ns = num_usm_v(3)
    10521060!
    1053 !--    Allocate required attributes for horizontal surfaces - default type. 
     1061!--    Allocate required attributes for horizontal surfaces - default type.
    10541062!--    Upward-facing (l=0) and downward-facing (l=1).
    10551063       DO  l = 0, 1
     
    10601068       CALL allocate_surface_attributes_h_top ( surf_def_h(2), nys, nyn, nxl, nxr )
    10611069!
    1062 !--    Allocate required attributes for horizontal surfaces - natural type. 
     1070!--    Allocate required attributes for horizontal surfaces - natural type.
    10631071       CALL allocate_surface_attributes_h ( surf_lsm_h, nys, nyn, nxl, nxr )
    10641072!
    1065 !--    Allocate required attributes for horizontal surfaces - urban type. 
     1073!--    Allocate required attributes for horizontal surfaces - urban type.
    10661074       CALL allocate_surface_attributes_h ( surf_usm_h, nys, nyn, nxl, nxr )
    10671075
    10681076!
    1069 !--    Allocate required attributes for vertical surfaces. 
     1077!--    Allocate required attributes for vertical surfaces.
    10701078!--    Northward-facing (l=0), southward-facing (l=1), eastward-facing (l=2)
    10711079!--    and westward-facing (l=3).
     
    11011109#endif
    11021110       IF ( num_surf_v > 0 )  vertical_surfaces_exist = .TRUE.
    1103        
     1111
    11041112
    11051113    END SUBROUTINE init_surface_arrays
     
    11171125
    11181126       INTEGER(iwp) ::  l     !<
    1119        
     1127
    11201128       !$ACC ENTER DATA &
    11211129       !$ACC COPYIN(surf_def_h(0:2)) &
     
    11621170
    11631171       INTEGER(iwp) ::  l     !<
    1164        
     1172
    11651173       ! Delete data in surf_def_h(0:2)
    11661174       DO  l = 0, 1
     
    11991207! Description:
    12001208! ------------
    1201 !> Deallocating memory for upward and downward-facing horizontal surface types, 
    1202 !> except for top fluxes. 
     1209!> Deallocating memory for upward and downward-facing horizontal surface types,
     1210!> except for top fluxes.
    12031211!------------------------------------------------------------------------------!
    12041212    SUBROUTINE deallocate_surface_attributes_h( surfaces )
     
    12421250!
    12431251!--    Vertical momentum fluxes of u and v
    1244        DEALLOCATE ( surfaces%usws ) 
    1245        DEALLOCATE ( surfaces%vsws ) 
     1252       DEALLOCATE ( surfaces%usws )
     1253       DEALLOCATE ( surfaces%vsws )
    12461254!
    12471255!--    Required in production_e
    1248        IF ( .NOT. constant_diffusion )  THEN   
    1249           DEALLOCATE ( surfaces%u_0 ) 
     1256       IF ( .NOT. constant_diffusion )  THEN
     1257          DEALLOCATE ( surfaces%u_0 )
    12501258          DEALLOCATE ( surfaces%v_0 )
    1251        ENDIF 
     1259       ENDIF
    12521260!
    12531261!--    Characteristic temperature and surface flux of sensible heat
    1254        DEALLOCATE ( surfaces%ts )   
     1262       DEALLOCATE ( surfaces%ts )
    12551263       DEALLOCATE ( surfaces%shf )
    12561264!
    12571265!--    surface temperature
    1258        DEALLOCATE ( surfaces%pt_surface ) 
     1266       DEALLOCATE ( surfaces%pt_surface )
    12591267!
    12601268!--    Characteristic humidity and surface flux of latent heat
    1261        IF ( humidity )  THEN         
    1262           DEALLOCATE ( surfaces%qs ) 
    1263           DEALLOCATE ( surfaces%qsws ) 
     1269       IF ( humidity )  THEN
     1270          DEALLOCATE ( surfaces%qs )
     1271          DEALLOCATE ( surfaces%qsws )
    12641272          DEALLOCATE ( surfaces%q_surface   )
    12651273          DEALLOCATE ( surfaces%vpt_surface )
    1266        ENDIF 
     1274       ENDIF
    12671275!
    12681276!--    Characteristic scalar and surface flux of scalar
    12691277       IF ( passive_scalar )  THEN
    1270           DEALLOCATE ( surfaces%ss )   
    1271           DEALLOCATE ( surfaces%ssws ) 
    1272        ENDIF 
     1278          DEALLOCATE ( surfaces%ss )
     1279          DEALLOCATE ( surfaces%ssws )
     1280       ENDIF
    12731281!
    12741282!--    Scaling parameter (cs*) and surface flux of chemical species
    12751283       IF ( air_chemistry )  THEN
    1276           DEALLOCATE ( surfaces%css )   
    1277           DEALLOCATE ( surfaces%cssws ) 
    1278        ENDIF 
     1284          DEALLOCATE ( surfaces%css )
     1285          DEALLOCATE ( surfaces%cssws )
     1286       ENDIF
    12791287!
    12801288!--    Arrays for storing potential temperature and
     
    12831291       DEALLOCATE ( surfaces%qv1 )
    12841292       DEALLOCATE ( surfaces%vpt1 )
    1285        
    1286 !
    1287 !--       
     1293
     1294!
     1295!--
    12881296       IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
    12891297          DEALLOCATE ( surfaces%qcs )
     
    12931301       ENDIF
    12941302!
    1295 !--       
     1303!--
    12961304       IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
    12971305          DEALLOCATE ( surfaces%qrs )
     
    13011309       ENDIF
    13021310!
     1311!--
     1312       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
     1313          DEALLOCATE ( surfaces%qis )
     1314          DEALLOCATE ( surfaces%nis )
     1315          DEALLOCATE ( surfaces%qisws )
     1316          DEALLOCATE ( surfaces%nisws )
     1317       ENDIF
     1318!
    13031319!--    Salinity surface flux
    13041320       IF ( ocean_mode )  DEALLOCATE ( surfaces%sasws )
     
    13101326! Description:
    13111327! ------------
    1312 !> Allocating memory for upward and downward-facing horizontal surface types, 
    1313 !> except for top fluxes. 
     1328!> Allocating memory for upward and downward-facing horizontal surface types,
     1329!> except for top fluxes.
    13141330!------------------------------------------------------------------------------!
    13151331    SUBROUTINE allocate_surface_attributes_h( surfaces,                        &
     
    13261342
    13271343!
    1328 !--    Allocate arrays for start and end index of horizontal surface type 
     1344!--    Allocate arrays for start and end index of horizontal surface type
    13291345!--    for each (j,i)-grid point. This is required e.g. in diffion_x, which is
    1330 !--    called for each (j,i). In order to find the location where the 
    1331 !--    respective flux is store within the surface-type, start- and end- 
     1346!--    called for each (j,i). In order to find the location where the
     1347!--    respective flux is store within the surface-type, start- and end-
    13321348!--    index are stored for each (j,i). For example, each (j,i) can have
    13331349!--    several entries where fluxes for horizontal surfaces might be stored,
     
    13701386!
    13711387!--    Vertical momentum fluxes of u and v
    1372        ALLOCATE ( surfaces%usws(1:surfaces%ns) ) 
    1373        ALLOCATE ( surfaces%vsws(1:surfaces%ns) ) 
     1388       ALLOCATE ( surfaces%usws(1:surfaces%ns) )
     1389       ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
    13741390!
    13751391!--    Required in production_e
    1376        IF ( .NOT. constant_diffusion )  THEN   
    1377           ALLOCATE ( surfaces%u_0(1:surfaces%ns) ) 
     1392       IF ( .NOT. constant_diffusion )  THEN
     1393          ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
    13781394          ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
    1379        ENDIF 
     1395       ENDIF
    13801396!
    13811397!--    Characteristic temperature and surface flux of sensible heat
    1382        ALLOCATE ( surfaces%ts(1:surfaces%ns)  )   
     1398       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )
    13831399       ALLOCATE ( surfaces%shf(1:surfaces%ns) )
    13841400!
    13851401!--    Surface temperature
    1386        ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) ) 
     1402       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
    13871403!
    13881404!--    Characteristic humidity, surface flux of latent heat, and surface virtual potential temperature
    13891405       IF ( humidity )  THEN
    1390           ALLOCATE ( surfaces%qs(1:surfaces%ns)   ) 
    1391           ALLOCATE ( surfaces%qsws(1:surfaces%ns) )     
    1392           ALLOCATE ( surfaces%q_surface(1:surfaces%ns)   ) 
    1393           ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) ) 
    1394        ENDIF 
     1406          ALLOCATE ( surfaces%qs(1:surfaces%ns)   )
     1407          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
     1408          ALLOCATE ( surfaces%q_surface(1:surfaces%ns)   )
     1409          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )
     1410       ENDIF
    13951411
    13961412!
    13971413!--    Characteristic scalar and surface flux of scalar
    13981414       IF ( passive_scalar )  THEN
    1399           ALLOCATE ( surfaces%ss(1:surfaces%ns)   )   
    1400           ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
    1401        ENDIF 
     1415          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )
     1416          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
     1417       ENDIF
    14021418!
    14031419!--    Scaling parameter (cs*) and surface flux of chemical species
    14041420       IF ( air_chemistry )  THEN
    1405           ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )   
    1406           ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
    1407        ENDIF 
     1421          ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )
     1422          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
     1423       ENDIF
    14081424!
    14091425!--    Arrays for storing potential temperature and
     
    14131429       ALLOCATE ( surfaces%vpt1(1:surfaces%ns) )
    14141430!
    1415 !--       
     1431!--
    14161432       IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
    14171433          ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
     
    14211437       ENDIF
    14221438!
    1423 !--       
     1439!--
    14241440       IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
    14251441          ALLOCATE ( surfaces%qrs(1:surfaces%ns)   )
     
    14281444          ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
    14291445       ENDIF
     1446
     1447!
     1448!--
     1449       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
     1450          ALLOCATE ( surfaces%qis(1:surfaces%ns)   )
     1451          ALLOCATE ( surfaces%nis(1:surfaces%ns)   )
     1452          ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
     1453          ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
     1454       ENDIF
     1455
    14301456!
    14311457!--    Salinity surface flux
     
    14451471
    14461472       IMPLICIT NONE
    1447    
     1473
    14481474       TYPE(surf_type) ::  surfaces  !< respective surface type
    1449    
     1475
    14501476       !$ACC EXIT DATA &
    14511477       !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
     
    14731499          !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
    14741500       ENDIF
    1475    
     1501
    14761502    END SUBROUTINE exit_surface_attributes_h
    14771503#endif
     
    15221548! Description:
    15231549! ------------
    1524 !> Deallocating memory for model-top fluxes 
     1550!> Deallocating memory for model-top fluxes
    15251551!------------------------------------------------------------------------------!
    15261552    SUBROUTINE deallocate_surface_attributes_h_top( surfaces )
     
    15391565       DEALLOCATE ( surfaces%k )
    15401566
    1541        IF ( .NOT. constant_diffusion )  THEN   
    1542           DEALLOCATE ( surfaces%u_0 ) 
     1567       IF ( .NOT. constant_diffusion )  THEN
     1568          DEALLOCATE ( surfaces%u_0 )
    15431569          DEALLOCATE ( surfaces%v_0 )
    1544        ENDIF 
     1570       ENDIF
    15451571!
    15461572!--    Vertical momentum fluxes of u and v
    1547        DEALLOCATE ( surfaces%usws ) 
    1548        DEALLOCATE ( surfaces%vsws ) 
     1573       DEALLOCATE ( surfaces%usws )
     1574       DEALLOCATE ( surfaces%vsws )
    15491575!
    15501576!--    Sensible heat flux
     
    15531579!--    Latent heat flux
    15541580       IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere')  THEN
    1555           DEALLOCATE ( surfaces%qsws )     
    1556        ENDIF 
     1581          DEALLOCATE ( surfaces%qsws )
     1582       ENDIF
    15571583!
    15581584!--    Scalar flux
    15591585       IF ( passive_scalar )  THEN
    1560           DEALLOCATE ( surfaces%ssws ) 
    1561        ENDIF 
     1586          DEALLOCATE ( surfaces%ssws )
     1587       ENDIF
    15621588!
    15631589!--    Chemical species flux
    15641590       IF ( air_chemistry )  THEN
    1565           DEALLOCATE ( surfaces%cssws ) 
    1566        ENDIF 
    1567 !
    1568 !--       
     1591          DEALLOCATE ( surfaces%cssws )
     1592       ENDIF
     1593!
     1594!--
    15691595       IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
    15701596          DEALLOCATE ( surfaces%qcsws )
     
    15721598       ENDIF
    15731599!
    1574 !--       
     1600!--
    15751601       IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
    15761602          DEALLOCATE ( surfaces%qrsws )
    15771603          DEALLOCATE ( surfaces%nrsws )
    15781604       ENDIF
     1605
     1606!
     1607!--
     1608       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
     1609          DEALLOCATE ( surfaces%qisws )
     1610          DEALLOCATE ( surfaces%nisws )
     1611       ENDIF
    15791612!
    15801613!--    Salinity flux
     
    15871620! Description:
    15881621! ------------
    1589 !> Allocating memory for model-top fluxes 
     1622!> Allocating memory for model-top fluxes
    15901623!------------------------------------------------------------------------------!
    15911624    SUBROUTINE allocate_surface_attributes_h_top( surfaces,                    &
     
    16111644       ALLOCATE ( surfaces%k(1:surfaces%ns)  )
    16121645
    1613        IF ( .NOT. constant_diffusion )  THEN   
    1614           ALLOCATE ( surfaces%u_0(1:surfaces%ns) ) 
     1646       IF ( .NOT. constant_diffusion )  THEN
     1647          ALLOCATE ( surfaces%u_0(1:surfaces%ns) )
    16151648          ALLOCATE ( surfaces%v_0(1:surfaces%ns) )
    1616        ENDIF 
     1649       ENDIF
    16171650!
    16181651!--    Vertical momentum fluxes of u and v
    1619        ALLOCATE ( surfaces%usws(1:surfaces%ns) ) 
    1620        ALLOCATE ( surfaces%vsws(1:surfaces%ns) ) 
     1652       ALLOCATE ( surfaces%usws(1:surfaces%ns) )
     1653       ALLOCATE ( surfaces%vsws(1:surfaces%ns) )
    16211654!
    16221655!--    Sensible heat flux
     
    16251658!--    Latent heat flux
    16261659       IF ( humidity .OR. coupling_mode == 'ocean_to_atmosphere')  THEN
    1627           ALLOCATE ( surfaces%qsws(1:surfaces%ns) )     
    1628        ENDIF 
     1660          ALLOCATE ( surfaces%qsws(1:surfaces%ns) )
     1661       ENDIF
    16291662!
    16301663!--    Scalar flux
    16311664       IF ( passive_scalar )  THEN
    1632           ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
    1633        ENDIF 
     1665          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
     1666       ENDIF
    16341667!
    16351668!--    Chemical species flux
    16361669       IF ( air_chemistry )  THEN
    1637           ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
    1638        ENDIF 
    1639 !
    1640 !--       
     1670          ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
     1671       ENDIF
     1672!
     1673!--
    16411674       IF ( surf_bulk_cloud_model .AND. surf_microphysics_morrison)  THEN
    16421675          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
     
    16441677       ENDIF
    16451678!
    1646 !--       
     1679!--
    16471680       IF ( surf_bulk_cloud_model .AND. surf_microphysics_seifert)  THEN
    16481681          ALLOCATE ( surfaces%qrsws(1:surfaces%ns) )
    16491682          ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
     1683       ENDIF
     1684
     1685!
     1686!--
     1687       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
     1688          ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
     1689          ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
    16501690       ENDIF
    16511691!
     
    16651705
    16661706       IMPLICIT NONE
    1667    
     1707
    16681708       TYPE(surf_type) ::  surfaces  !< respective surface type
    1669    
     1709
    16701710       !$ACC EXIT DATA &
    16711711       !$ACC DELETE(surfaces%start_index(nys:nyn,nxl:nxr)) &
     
    16831723          !$ACC DELETE(surfaces%v_0(1:surfaces%ns))
    16841724       ENDIF
    1685    
     1725
    16861726    END SUBROUTINE exit_surface_attributes_h_top
    16871727#endif
     
    17211761! Description:
    17221762! ------------
    1723 !> Deallocating memory for vertical surface types. 
     1763!> Deallocating memory for vertical surface types.
    17241764!------------------------------------------------------------------------------!
    17251765    SUBROUTINE deallocate_surface_attributes_v( surfaces )
     
    17311771
    17321772!
    1733 !--    Allocate arrays for start and end index of vertical surface type 
     1773!--    Allocate arrays for start and end index of vertical surface type
    17341774!--    for each (j,i)-grid point. This is required in diffion_x, which is
    1735 !--    called for each (j,i). In order to find the location where the 
    1736 !--    respective flux is store within the surface-type, start- and end- 
     1775!--    called for each (j,i). In order to find the location where the
     1776!--    respective flux is store within the surface-type, start- and end-
    17371777!--    index are stored for each (j,i). For example, each (j,i) can have
    1738 !--    several entries where fluxes for vertical surfaces might be stored. 
    1739 !--    In the flat case, where no vertical walls exit, set indicies such 
    1740 !--    that loop in diffusion routines will not be entered. 
     1778!--    several entries where fluxes for vertical surfaces might be stored.
     1779!--    In the flat case, where no vertical walls exit, set indicies such
     1780!--    that loop in diffusion routines will not be entered.
    17411781       DEALLOCATE ( surfaces%start_index )
    17421782       DEALLOCATE ( surfaces%end_index )
     
    17661806!
    17671807!--    Allocate Obukhov length and bulk Richardson number. Actually, at
    1768 !--    vertical surfaces these are only required for natural surfaces. 
     1808!--    vertical surfaces these are only required for natural surfaces.
    17691809!--    for natural land surfaces
    1770        DEALLOCATE( surfaces%ol ) 
    1771        DEALLOCATE( surfaces%rib ) 
    1772 !
    1773 !--    Allocate arrays for surface momentum fluxes for u and v. For u at north- 
     1810       DEALLOCATE( surfaces%ol )
     1811       DEALLOCATE( surfaces%rib )
     1812!
     1813!--    Allocate arrays for surface momentum fluxes for u and v. For u at north-
    17741814!--    and south-facing surfaces, for v at east- and west-facing surfaces.
    17751815       DEALLOCATE ( surfaces%mom_flux_uv )
    17761816!
    17771817!--    Allocate array for surface momentum flux for w - wsus and wsvs
    1778        DEALLOCATE ( surfaces%mom_flux_w ) 
    1779 !
    1780 !--    Allocate array for surface momentum flux for subgrid-scale tke wsus and 
     1818       DEALLOCATE ( surfaces%mom_flux_w )
     1819!
     1820!--    Allocate array for surface momentum flux for subgrid-scale tke wsus and
    17811821!--    wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
    17821822!--    on surface.
    1783        DEALLOCATE ( surfaces%mom_flux_tke ) 
     1823       DEALLOCATE ( surfaces%mom_flux_tke )
    17841824!
    17851825!--    Characteristic temperature and surface flux of sensible heat
    1786        DEALLOCATE ( surfaces%ts )   
     1826       DEALLOCATE ( surfaces%ts )
    17871827       DEALLOCATE ( surfaces%shf )
    17881828!
    17891829!--    surface temperature
    1790        DEALLOCATE ( surfaces%pt_surface ) 
     1830       DEALLOCATE ( surfaces%pt_surface )
    17911831!
    17921832!--    Characteristic humidity and surface flux of latent heat
    17931833       IF ( humidity )  THEN
    1794           DEALLOCATE ( surfaces%qs ) 
    1795           DEALLOCATE ( surfaces%qsws ) 
     1834          DEALLOCATE ( surfaces%qs )
     1835          DEALLOCATE ( surfaces%qsws )
    17961836          DEALLOCATE ( surfaces%q_surface   )
    17971837          DEALLOCATE ( surfaces%vpt_surface )
    1798        ENDIF 
     1838       ENDIF
    17991839!
    18001840!--    Characteristic scalar and surface flux of scalar
    18011841       IF ( passive_scalar )  THEN
    1802           DEALLOCATE ( surfaces%ss )   
    1803           DEALLOCATE ( surfaces%ssws ) 
     1842          DEALLOCATE ( surfaces%ss )
     1843          DEALLOCATE ( surfaces%ssws )
    18041844       ENDIF
    18051845!
    18061846!--    Scaling parameter (cs*) and surface flux of chemical species
    18071847       IF ( air_chemistry )  THEN
    1808              DEALLOCATE ( surfaces%css )   
    1809              DEALLOCATE ( surfaces%cssws ) 
    1810        ENDIF 
     1848             DEALLOCATE ( surfaces%css )
     1849             DEALLOCATE ( surfaces%cssws )
     1850       ENDIF
    18111851!
    18121852!--    Arrays for storing potential temperature and
     
    18291869          DEALLOCATE ( surfaces%nrsws )
    18301870       ENDIF
     1871
     1872       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
     1873          DEALLOCATE ( surfaces%qis )
     1874          DEALLOCATE ( surfaces%nis )
     1875          DEALLOCATE ( surfaces%qisws )
     1876          DEALLOCATE ( surfaces%nisws )
     1877       ENDIF
     1878
    18311879!
    18321880!--    Salinity surface flux
     
    18391887! Description:
    18401888! ------------
    1841 !> Allocating memory for vertical surface types. 
     1889!> Allocating memory for vertical surface types.
    18421890!------------------------------------------------------------------------------!
    18431891    SUBROUTINE allocate_surface_attributes_v( surfaces,                        &
     
    18541902
    18551903!
    1856 !--    Allocate arrays for start and end index of vertical surface type 
     1904!--    Allocate arrays for start and end index of vertical surface type
    18571905!--    for each (j,i)-grid point. This is required in diffion_x, which is
    1858 !--    called for each (j,i). In order to find the location where the 
    1859 !--    respective flux is store within the surface-type, start- and end- 
     1906!--    called for each (j,i). In order to find the location where the
     1907!--    respective flux is store within the surface-type, start- and end-
    18601908!--    index are stored for each (j,i). For example, each (j,i) can have
    1861 !--    several entries where fluxes for vertical surfaces might be stored. 
    1862 !--    In the flat case, where no vertical walls exit, set indicies such 
    1863 !--    that loop in diffusion routines will not be entered. 
     1909!--    several entries where fluxes for vertical surfaces might be stored.
     1910!--    In the flat case, where no vertical walls exit, set indicies such
     1911!--    that loop in diffusion routines will not be entered.
    18641912       ALLOCATE ( surfaces%start_index(nys_l:nyn_l,nxl_l:nxr_l) )
    18651913       ALLOCATE ( surfaces%end_index(nys_l:nyn_l,nxl_l:nxr_l)   )
     
    18911939!
    18921940!--    Allocate Obukhov length and bulk Richardson number. Actually, at
    1893 !--    vertical surfaces these are only required for natural surfaces. 
     1941!--    vertical surfaces these are only required for natural surfaces.
    18941942!--    for natural land surfaces
    1895        ALLOCATE( surfaces%ol(1:surfaces%ns)  ) 
    1896        ALLOCATE( surfaces%rib(1:surfaces%ns) ) 
    1897 !
    1898 !--    Allocate arrays for surface momentum fluxes for u and v. For u at north- 
     1943       ALLOCATE( surfaces%ol(1:surfaces%ns)  )
     1944       ALLOCATE( surfaces%rib(1:surfaces%ns) )
     1945!
     1946!--    Allocate arrays for surface momentum fluxes for u and v. For u at north-
    18991947!--    and south-facing surfaces, for v at east- and west-facing surfaces.
    19001948       ALLOCATE ( surfaces%mom_flux_uv(1:surfaces%ns) )
    19011949!
    19021950!--    Allocate array for surface momentum flux for w - wsus and wsvs
    1903        ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) ) 
    1904 !
    1905 !--    Allocate array for surface momentum flux for subgrid-scale tke wsus and 
     1951       ALLOCATE ( surfaces%mom_flux_w(1:surfaces%ns) )
     1952!
     1953!--    Allocate array for surface momentum flux for subgrid-scale tke wsus and
    19061954!--    wsvs; first index usvs or vsws, second index for wsus or wsvs, depending
    19071955!--    on surface.
    1908        ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) ) 
     1956       ALLOCATE ( surfaces%mom_flux_tke(0:1,1:surfaces%ns) )
    19091957!
    19101958!--    Characteristic temperature and surface flux of sensible heat
    1911        ALLOCATE ( surfaces%ts(1:surfaces%ns)  )   
     1959       ALLOCATE ( surfaces%ts(1:surfaces%ns)  )
    19121960       ALLOCATE ( surfaces%shf(1:surfaces%ns) )
    19131961!
    19141962!--    surface temperature
    1915        ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) ) 
     1963       ALLOCATE ( surfaces%pt_surface(1:surfaces%ns) )
    19161964!
    19171965!--    Characteristic humidity and surface flux of latent heat
    19181966       IF ( humidity )  THEN
    1919           ALLOCATE ( surfaces%qs(1:surfaces%ns)          ) 
    1920           ALLOCATE ( surfaces%qsws(1:surfaces%ns)        )   
     1967          ALLOCATE ( surfaces%qs(1:surfaces%ns)          )
     1968          ALLOCATE ( surfaces%qsws(1:surfaces%ns)        )
    19211969          ALLOCATE ( surfaces%q_surface(1:surfaces%ns)   )
    1922           ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )           
    1923        ENDIF 
     1970          ALLOCATE ( surfaces%vpt_surface(1:surfaces%ns) )
     1971       ENDIF
    19241972!
    19251973!--    Characteristic scalar and surface flux of scalar
    19261974       IF ( passive_scalar )  THEN
    1927           ALLOCATE ( surfaces%ss(1:surfaces%ns)   )   
    1928           ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
     1975          ALLOCATE ( surfaces%ss(1:surfaces%ns)   )
     1976          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
    19291977       ENDIF
    19301978!
    19311979!--    Scaling parameter (cs*) and surface flux of chemical species
    19321980       IF ( air_chemistry )  THEN
    1933              ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )   
    1934              ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) ) 
    1935        ENDIF 
     1981             ALLOCATE ( surfaces%css(1:nvar,1:surfaces%ns)   )
     1982             ALLOCATE ( surfaces%cssws(1:nvar,1:surfaces%ns) )
     1983       ENDIF
    19361984!
    19371985!--    Arrays for storing potential temperature and
     
    19542002          ALLOCATE ( surfaces%nrsws(1:surfaces%ns) )
    19552003       ENDIF
     2004
     2005       IF ( surf_bulk_cloud_model .AND. surf_microphysics_ice_extension)  THEN
     2006          ALLOCATE ( surfaces%qis(1:surfaces%ns)   )
     2007          ALLOCATE ( surfaces%nis(1:surfaces%ns)   )
     2008          ALLOCATE ( surfaces%qisws(1:surfaces%ns) )
     2009          ALLOCATE ( surfaces%nisws(1:surfaces%ns) )
     2010       ENDIF
    19562011!
    19572012!--    Salinity surface flux
     
    19642019! Description:
    19652020! ------------
    1966 !> Exit memory for vertical surface types. 
     2021!> Exit memory for vertical surface types.
    19672022!------------------------------------------------------------------------------!
    19682023#if defined( _OPENACC )
     
    19962051! Description:
    19972052! ------------
    1998 !> Enter memory for vertical surface types. 
     2053!> Enter memory for vertical surface types.
    19992054!------------------------------------------------------------------------------!
    20002055#if defined( _OPENACC )
    20012056    SUBROUTINE enter_surface_attributes_v( surfaces )
    2002    
     2057
    20032058       IMPLICIT NONE
    2004    
     2059
    20052060       TYPE(surf_type) ::  surfaces  !< respective surface type
    2006    
     2061
    20072062       !$ACC ENTER DATA &
    20082063       !$ACC COPYIN(surfaces%start_index(nys:nyn,nxl:nxr)) &
     
    20212076       !$ACC COPYIN(surfaces%pt1(1:surfaces%ns)) &
    20222077       !$ACC COPYIN(surfaces%qv1(1:surfaces%ns))
    2023    
     2078
    20242079    END SUBROUTINE enter_surface_attributes_v
    20252080#endif
     
    20282083! Description:
    20292084! ------------
    2030 !> Initialize surface elements, i.e. set initial values for surface fluxes, 
    2031 !> friction velocity, calcuation of start/end indices, etc. .