source: palm/trunk/SOURCE/bulk_cloud_model_mod.f90 @ 3636

Last change on this file since 3636 was 3636, checked in by raasch, 3 years ago

nopointer option removed

  • Property svn:keywords set to Id
File size: 160.8 KB
Line 
1!> @file bulk_cloud_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: bulk_cloud_model_mod.f90 3636 2018-12-19 13:48:34Z raasch $
27! nopointer option removed
28!
29! 3622 2018-12-12 09:52:53Z schwenkel
30! Important bugfix in case of restart runs.
31!
32! 3507 2018-11-08 14:28:46Z schwenkel
33! Minor bugfixes for bcm cache version and initializing prr
34!
35! 3452 2018-10-30 13:13:34Z schwenkel
36! Bugfix for profiles output
37!
38! 3445 2018-10-29 12:23:02Z schwenkel
39! Minor bugfix and use of subroutine for supersaturation calculation in case
40! of cache version
41!
42! 3383 2018-10-19 14:22:58Z knoop
43! Modularization of all bulk cloud physics code components
44!
45! unused variables removed
46!
47! 3026 2018-05-22 10:30:53Z schwenkel
48! Changed the name specific humidity to mixing ratio, since we are computing
49! mixing ratios.
50!
51! 2718 2018-01-02 08:49:38Z maronga
52! Corrected "Former revisions" section
53!
54! 2701 2017-12-15 15:40:50Z suehring
55! Changes from last commit documented
56!
57! 2698 2017-12-14 18:46:24Z suehring
58! Bugfix in get_topography_top_index
59!
60! 2696 2017-12-14 17:12:51Z kanani
61! Change in file header (GPL part)
62!
63! 2608 2017-11-13 14:04:26Z schwenkel
64! Calculation of supersaturation in external module (diagnostic_quantities_mod).
65! Change: correct calculation of saturation specific humidity to saturation
66! mixing ratio (the factor of 0.378 vanishes).
67!
68! 2522 2017-10-05 14:20:37Z schwenkel
69! Minor bugfix
70!
71! 2375 2017-08-29 14:10:28Z schwenkel
72! Improved aerosol initilization and some minor bugfixes
73! for droplet sedimenation
74!
75! 2318 2017-07-20 17:27:44Z suehring
76! Get topography top index via Function call
77!
78! 2317 2017-07-20 17:27:19Z suehring
79! s1 changed to log_sigma
80!
81! 2292 2017-06-20 09:51:42Z schwenkel
82! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
83! includes two more prognostic equations for cloud drop concentration (nc)
84! and cloud water content (qc).
85! - The process of activation is parameterized with a simple Twomey
86!   activion scheme or with considering solution and curvature
87!   effects (Khvorostyanov and Curry ,2006).
88! - The saturation adjustment scheme is replaced by the parameterization
89!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
90! - All other microphysical processes of Seifert and Beheng are used.
91!   Additionally, in those processes the reduction of cloud number concentration
92!   is considered.
93!
94! 2233 2017-05-30 18:08:54Z suehring
95!
96! 2232 2017-05-30 17:47:52Z suehring
97! Adjustments to new topography and surface concept
98!
99! 2155 2017-02-21 09:57:40Z hoffmann
100! Bugfix in the calculation of microphysical quantities on ghost points.
101!
102! 2031 2016-10-21 15:11:58Z knoop
103! renamed variable rho to rho_ocean
104!
105! 2000 2016-08-20 18:09:15Z knoop
106! Forced header and separation lines into 80 columns
107!
108! 1850 2016-04-08 13:29:27Z maronga
109! Module renamed
110! Adapted for modularization of microphysics.
111!
112! 1845 2016-04-08 08:29:13Z raasch
113! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
114! point (instead of one point above surface)
115!
116! 1831 2016-04-07 13:15:51Z hoffmann
117! turbulence renamed collision_turbulence,
118! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
119! avaialble for microphysics_kessler.
120!
121! 1822 2016-04-07 07:49:42Z hoffmann
122! Unused variables removed.
123! Kessler scheme integrated.
124!
125! 1691 2015-10-26 16:17:44Z maronga
126! Added new routine calc_precipitation_amount. The routine now allows to account
127! for precipitation due to sedimenation of cloud (fog) droplets
128!
129! 1682 2015-10-07 23:56:08Z knoop
130! Code annotations made doxygen readable
131!
132! 1646 2015-09-02 16:00:10Z hoffmann
133! Bugfix: Wrong computation of d_mean.
134!
135! 1361 2014-04-16 15:17:48Z hoffmann
136! Bugfix in sedimentation_rain: Index corrected.
137! Vectorized version of adjust_cloud added.
138! Little reformatting of the code.
139!
140! 1353 2014-04-08 15:21:23Z heinze
141! REAL constants provided with KIND-attribute
142!
143! 1346 2014-03-27 13:18:20Z heinze
144! Bugfix: REAL constants provided with KIND-attribute especially in call of
145! intrinsic function like MAX, MIN, SIGN
146!
147! 1334 2014-03-25 12:21:40Z heinze
148! Bugfix: REAL constants provided with KIND-attribute
149!
150! 1322 2014-03-20 16:38:49Z raasch
151! REAL constants defined as wp-kind
152!
153! 1320 2014-03-20 08:40:49Z raasch
154! ONLY-attribute added to USE-statements,
155! kind-parameters added to all INTEGER and REAL declaration statements,
156! kinds are defined in new module kinds,
157! comment fields (!:) to be used for variable explanations added to
158! all variable declaration statements
159!
160! 1241 2013-10-30 11:36:58Z heinze
161! hyp and rho_ocean have to be calculated at each time step if data from external
162! file LSF_DATA are used
163!
164! 1115 2013-03-26 18:16:16Z hoffmann
165! microphyical tendencies are calculated in bcm_actions in an optimized
166! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
167!
168! 1106 2013-03-04 05:31:38Z raasch
169! small changes in code formatting
170!
171! 1092 2013-02-02 11:24:22Z raasch
172! unused variables removed
173! file put under GPL
174!
175! 1065 2012-11-22 17:42:36Z hoffmann
176! Sedimentation process implemented according to Stevens and Seifert (2008).
177! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
178! and Stevens, 2010).
179!
180! 1053 2012-11-13 17:11:03Z hoffmann
181! initial revision
182!
183! Description:
184! ------------
185!> Calculate bulk cloud microphysics.
186!------------------------------------------------------------------------------!
187 MODULE bulk_cloud_model_mod
188
189    USE arrays_3d,                                                             &
190        ONLY:  ddzu, diss, dzu, dzw, hyp, hyrho,                               &
191               nc, nc_1, nc_2, nc_3, nc_p, nr, nr_1, nr_2, nr_3, nr_p,         &
192               precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1,   &
193               qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p,         &
194               exner, zu, tnc_m, tnr_m, tqc_m, tqr_m
195
196    USE averaging,                                                             &
197        ONLY:  nc_av, nr_av, prr_av, qc_av, ql_av, qr_av
198
199    USE basic_constants_and_equations_mod,                                     &
200        ONLY:  c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,&
201               molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,&
202               exner_function, exner_function_invers, ideal_gas_law_rho,       &
203               ideal_gas_law_rho_pt, barometric_formula, rd_d_rv
204
205    USE control_parameters,                                                    &
206        ONLY:  dt_3d, dt_do2d_xy, intermediate_timestep_count,                 &
207               intermediate_timestep_count_max, large_scale_forcing,           &
208               lsf_surf, pt_surface, rho_surface, surface_pressure,            &
209               time_do2d_xy, message_string, initializing_actions
210
211    USE cpulog,                                                                &
212        ONLY:  cpu_log, log_point_s
213
214    USE grid_variables,                                                        &
215        ONLY:  dx, dy
216
217    USE indices,                                                               &
218        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt,           &
219               wall_flags_0
220
221    USE  kinds
222
223    USE statistics,                                                            &
224        ONLY:  weight_pres, weight_substep
225
226    USE surface_mod,                                                           &
227        ONLY :  bc_h, get_topography_top_index_ji, surf_bulk_cloud_model,      &
228                surf_microphysics_morrison, surf_microphysics_seifert
229
230    IMPLICIT NONE
231
232    CHARACTER (LEN=20)   ::  aerosol_bulk = 'nacl'                        !< namelist parameter
233    CHARACTER (LEN=20)   ::  cloud_scheme = 'saturation_adjust'           !< namelist parameter
234
235    LOGICAL ::  aerosol_nacl =.TRUE.                             !< nacl aerosol for bulk scheme
236    LOGICAL ::  aerosol_c3h4o4 =.FALSE.                          !< malonic acid aerosol for bulk scheme
237    LOGICAL ::  aerosol_nh4no3 =.FALSE.                          !< malonic acid aerosol for bulk scheme
238
239    LOGICAL ::  bulk_cloud_model = .FALSE.                       !< namelist parameter
240
241    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
242    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
243    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
244    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
245    LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
246
247    LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.      !< namelist parameter
248    LOGICAL ::  microphysics_sat_adjust = .FALSE.                !< use saturation adjust bulk scheme?
249    LOGICAL ::  microphysics_kessler = .FALSE.                   !< use kessler bulk scheme?
250    LOGICAL ::  microphysics_morrison = .FALSE.                  !< use 2-moment Morrison (add. prog. eq. for nc and qc)
251    LOGICAL ::  microphysics_seifert = .FALSE.                   !< use 2-moment Seifert and Beheng scheme
252    LOGICAL ::  precipitation = .FALSE.                          !< namelist parameter
253
254    REAL(wp) ::  precipitation_amount_interval = 9999999.9_wp    !< namelist parameter
255
256    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
257    REAL(wp) ::  a_2 = -7.38E-5_wp         !< coef. in turb. parametrization (cm-2 s3)
258    REAL(wp) ::  a_3 = -1.40E-2_wp         !< coef. in turb. parametrization
259    REAL(wp) ::  a_term = 9.65_wp          !< coef. for terminal velocity (m s-1)
260    REAL(wp) ::  a_vent = 0.78_wp          !< coef. for ventilation effect
261    REAL(wp) ::  b_1 = 11.45E-6_wp         !< coef. in turb. parametrization (m)
262    REAL(wp) ::  b_2 = 9.68E-6_wp          !< coef. in turb. parametrization (m)
263    REAL(wp) ::  b_3 = 0.62_wp             !< coef. in turb. parametrization
264    REAL(wp) ::  b_term = 9.8_wp           !< coef. for terminal velocity (m s-1)
265    REAL(wp) ::  b_vent = 0.308_wp         !< coef. for ventilation effect
266    REAL(wp) ::  beta_cc = 3.09E-4_wp      !< coef. in turb. parametrization (cm-2 s3)
267    REAL(wp) ::  c_1 = 4.82E-6_wp          !< coef. in turb. parametrization (m)
268    REAL(wp) ::  c_2 = 4.8E-6_wp           !< coef. in turb. parametrization (m)
269    REAL(wp) ::  c_3 = 0.76_wp             !< coef. in turb. parametrization
270    REAL(wp) ::  c_const = 0.93_wp         !< const. in Taylor-microscale Reynolds number
271    REAL(wp) ::  c_evap = 0.7_wp           !< constant in evaporation
272    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
273    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
274    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
275    REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
276    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
277    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
278    REAL(wp) ::  k_rr = 7.12_wp            !< const. rain-rain kernel (m3 kg-1 s-1)
279    REAL(wp) ::  k_br = 1000.0_wp          !< const. in breakup parametrization (m-1)
280    REAL(wp) ::  k_st = 1.2E8_wp           !< const. in drizzle parametrization (m-1 s-1)
281    REAL(wp) ::  kin_vis_air = 1.4086E-5_wp  !< kin. viscosity of air (m2 s-1)
282    REAL(wp) ::  prec_time_const = 0.001_wp  !< coef. in Kessler scheme (s-1)
283    REAL(wp) ::  ql_crit = 0.0005_wp       !< coef. in Kessler scheme (kg kg-1)
284    REAL(wp) ::  schmidt_p_1d3=0.8921121_wp  !< Schmidt number**0.33333, 0.71**0.33333
285    REAL(wp) ::  sigma_gc = 1.3_wp         !< geometric standard deviation cloud droplets
286    REAL(wp) ::  thermal_conductivity_l = 2.43E-2_wp  !< therm. cond. air (J m-1 s-1 K-1)
287    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
288    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
289    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
290    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
291    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
292    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
293
294    REAL(wp) ::  c_sedimentation = 2.0_wp        !< Courant number of sedimentation process
295    REAL(wp) ::  dpirho_l                        !< 6.0 / ( pi * rho_l )
296    REAL(wp) ::  dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
297    REAL(wp) ::  dt_micro                        !< microphysics time step
298    REAL(wp) ::  sigma_bulk = 2.0_wp             !< width of aerosol spectrum
299    REAL(wp) ::  na_init = 100.0E6_wp            !< Total particle/aerosol concentration (cm-3)
300    REAL(wp) ::  nc_const = 70.0E6_wp            !< cloud droplet concentration
301    REAL(wp) ::  dt_precipitation = 100.0_wp     !< timestep precipitation (s)
302    REAL(wp) ::  sed_qc_const                    !< const. for sedimentation of cloud water
303    REAL(wp) ::  pirho_l                         !< pi * rho_l / 6.0;
304
305    REAL(wp) ::  e_s     !< saturation water vapor pressure
306    REAL(wp) ::  q_s     !< saturation mixing ratio
307    REAL(wp) ::  sat     !< supersaturation
308    REAL(wp) ::  t_l     !< actual temperature
309
310    SAVE
311
312    PRIVATE
313
314    PUBLIC bcm_parin, &
315           bcm_check_parameters, &
316           bcm_check_data_output, &
317           bcm_check_data_output_pr, &
318           bcm_header, &
319           bcm_init_arrays, &
320           bcm_init, &
321           bcm_3d_data_averaging, &
322           bcm_data_output_2d, &
323           bcm_data_output_3d, &
324           bcm_swap_timelevel, &
325           bcm_rrd_global, &
326           bcm_rrd_local, &
327           bcm_wrd_global, &
328           bcm_wrd_local, &
329           bcm_actions, &
330           calc_liquid_water_content
331
332    PUBLIC call_microphysics_at_all_substeps, &
333           cloud_water_sedimentation, &
334           bulk_cloud_model, &
335           cloud_scheme, &
336           collision_turbulence, &
337           dt_precipitation, &
338           microphysics_morrison, &
339           microphysics_sat_adjust, &
340           microphysics_seifert, &
341           na_init, &
342           nc_const, &
343           precipitation, &
344           sigma_gc
345
346
347    INTERFACE bcm_parin
348       MODULE PROCEDURE bcm_parin
349    END INTERFACE bcm_parin
350
351    INTERFACE bcm_check_parameters
352       MODULE PROCEDURE bcm_check_parameters
353    END INTERFACE bcm_check_parameters
354
355    INTERFACE bcm_check_data_output
356       MODULE PROCEDURE bcm_check_data_output
357    END INTERFACE bcm_check_data_output
358
359    INTERFACE bcm_check_data_output_pr
360       MODULE PROCEDURE bcm_check_data_output_pr
361    END INTERFACE bcm_check_data_output_pr
362
363    INTERFACE bcm_header
364       MODULE PROCEDURE bcm_header
365    END INTERFACE bcm_header
366
367    INTERFACE bcm_init_arrays
368       MODULE PROCEDURE bcm_init_arrays
369    END INTERFACE bcm_init_arrays
370
371    INTERFACE bcm_init
372       MODULE PROCEDURE bcm_init
373    END INTERFACE bcm_init
374
375    INTERFACE bcm_3d_data_averaging
376       MODULE PROCEDURE bcm_3d_data_averaging
377    END INTERFACE bcm_3d_data_averaging
378
379    INTERFACE bcm_data_output_2d
380       MODULE PROCEDURE bcm_data_output_2d
381    END INTERFACE bcm_data_output_2d
382
383    INTERFACE bcm_data_output_3d
384       MODULE PROCEDURE bcm_data_output_3d
385    END INTERFACE bcm_data_output_3d
386
387    INTERFACE bcm_swap_timelevel
388       MODULE PROCEDURE bcm_swap_timelevel
389    END INTERFACE bcm_swap_timelevel
390
391    INTERFACE bcm_rrd_global
392       MODULE PROCEDURE bcm_rrd_global
393    END INTERFACE bcm_rrd_global
394
395    INTERFACE bcm_rrd_local
396       MODULE PROCEDURE bcm_rrd_local
397    END INTERFACE bcm_rrd_local
398
399    INTERFACE bcm_wrd_global
400       MODULE PROCEDURE bcm_wrd_global
401    END INTERFACE bcm_wrd_global
402
403    INTERFACE bcm_wrd_local
404       MODULE PROCEDURE bcm_wrd_local
405    END INTERFACE bcm_wrd_local
406
407    INTERFACE bcm_actions
408       MODULE PROCEDURE bcm_actions
409       MODULE PROCEDURE bcm_actions_ij
410    END INTERFACE bcm_actions
411
412    INTERFACE adjust_cloud
413       MODULE PROCEDURE adjust_cloud
414       MODULE PROCEDURE adjust_cloud_ij
415    END INTERFACE adjust_cloud
416
417    INTERFACE activation
418       MODULE PROCEDURE activation
419       MODULE PROCEDURE activation_ij
420    END INTERFACE activation
421
422    INTERFACE condensation
423       MODULE PROCEDURE condensation
424       MODULE PROCEDURE condensation_ij
425    END INTERFACE condensation
426
427    INTERFACE autoconversion
428       MODULE PROCEDURE autoconversion
429       MODULE PROCEDURE autoconversion_ij
430    END INTERFACE autoconversion
431
432    INTERFACE autoconversion_kessler
433       MODULE PROCEDURE autoconversion_kessler
434       MODULE PROCEDURE autoconversion_kessler_ij
435    END INTERFACE autoconversion_kessler
436
437    INTERFACE accretion
438       MODULE PROCEDURE accretion
439       MODULE PROCEDURE accretion_ij
440    END INTERFACE accretion
441
442    INTERFACE selfcollection_breakup
443       MODULE PROCEDURE selfcollection_breakup
444       MODULE PROCEDURE selfcollection_breakup_ij
445    END INTERFACE selfcollection_breakup
446
447    INTERFACE evaporation_rain
448       MODULE PROCEDURE evaporation_rain
449       MODULE PROCEDURE evaporation_rain_ij
450    END INTERFACE evaporation_rain
451
452    INTERFACE sedimentation_cloud
453       MODULE PROCEDURE sedimentation_cloud
454       MODULE PROCEDURE sedimentation_cloud_ij
455    END INTERFACE sedimentation_cloud
456
457    INTERFACE sedimentation_rain
458       MODULE PROCEDURE sedimentation_rain
459       MODULE PROCEDURE sedimentation_rain_ij
460    END INTERFACE sedimentation_rain
461
462    INTERFACE calc_precipitation_amount
463       MODULE PROCEDURE calc_precipitation_amount
464       MODULE PROCEDURE calc_precipitation_amount_ij
465    END INTERFACE calc_precipitation_amount
466
467    INTERFACE supersaturation
468       MODULE PROCEDURE supersaturation
469    END INTERFACE supersaturation
470
471    INTERFACE calc_liquid_water_content
472       MODULE PROCEDURE calc_liquid_water_content
473    END INTERFACE calc_liquid_water_content
474
475 CONTAINS
476
477
478!------------------------------------------------------------------------------!
479! Description:
480! ------------
481!> Parin for &bulk_cloud_parameters for the bulk cloud module
482!------------------------------------------------------------------------------!
483    SUBROUTINE bcm_parin
484
485
486       IMPLICIT NONE
487
488       CHARACTER (LEN=80)  ::  line  !< dummy string that contains the current line of the parameter file
489
490       NAMELIST /bulk_cloud_parameters/  &
491          aerosol_bulk, &
492          c_sedimentation, &
493          call_microphysics_at_all_substeps, &
494          bulk_cloud_model, &
495          cloud_scheme, &
496          cloud_water_sedimentation, &
497          collision_turbulence, &
498          curvature_solution_effects_bulk, &
499          dry_aerosol_radius, &
500          limiter_sedimentation, &
501          na_init, &
502          nc_const, &
503          precipitation, &
504          precipitation_amount_interval, &
505          sigma_bulk, &
506          ventilation_effect
507
508       line = ' '
509!
510!--    Try to find bulk cloud module namelist
511       REWIND ( 11 )
512       line = ' '
513       DO   WHILE ( INDEX( line, '&bulk_cloud_parameters' ) == 0 )
514          READ ( 11, '(A)', END=10 )  line
515       ENDDO
516       BACKSPACE ( 11 )
517!
518!--    Read user-defined namelist
519       READ ( 11, bulk_cloud_parameters )
520!
521!--    Set flag that indicates that the bulk cloud module is switched on
522       !bulk_cloud_model = .TRUE.
523
52410     CONTINUE
525
526
527    END SUBROUTINE bcm_parin
528
529
530!------------------------------------------------------------------------------!
531! Description:
532! ------------
533!> Check parameters routine for bulk cloud module
534!------------------------------------------------------------------------------!
535    SUBROUTINE bcm_check_parameters
536
537
538       IMPLICIT NONE
539!
540!--    Check cloud scheme
541       IF ( cloud_scheme == 'saturation_adjust' )  THEN
542          microphysics_sat_adjust = .TRUE.
543          microphysics_seifert    = .FALSE.
544          microphysics_kessler    = .FALSE.
545          precipitation           = .FALSE.
546       ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
547          microphysics_sat_adjust = .FALSE.
548          microphysics_seifert    = .TRUE.
549          microphysics_kessler    = .FALSE.
550          microphysics_morrison  = .FALSE.
551          precipitation           = .TRUE.
552       ELSEIF ( cloud_scheme == 'kessler' )  THEN
553          microphysics_sat_adjust = .FALSE.
554          microphysics_seifert    = .FALSE.
555          microphysics_kessler    = .TRUE.
556          microphysics_morrison   = .FALSE.
557          precipitation           = .TRUE.
558       ELSEIF ( cloud_scheme == 'morrison' )  THEN
559          microphysics_sat_adjust = .FALSE.
560          microphysics_seifert    = .TRUE.
561          microphysics_kessler    = .FALSE.
562          microphysics_morrison   = .TRUE.
563          precipitation           = .TRUE.
564       ELSE
565          message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
566                           TRIM( cloud_scheme ) // '"'
567          CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
568       ENDIF
569
570
571
572!
573!--    Set the default value for the integration interval of precipitation amount
574       IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
575          IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
576             precipitation_amount_interval = dt_do2d_xy
577          ELSE
578             IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
579                WRITE( message_string, * )  'precipitation_amount_interval = ',   &
580                    precipitation_amount_interval, ' must not be larger than ',   &
581                    'dt_do2d_xy = ', dt_do2d_xy
582                CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
583             ENDIF
584          ENDIF
585       ENDIF
586
587       ! TODO: find better sollution for circular dependency problem
588       surf_bulk_cloud_model = bulk_cloud_model
589       surf_microphysics_morrison = microphysics_morrison
590       surf_microphysics_seifert = microphysics_seifert
591
592!
593!--    Check aerosol
594       IF ( aerosol_bulk == 'nacl' )  THEN
595          aerosol_nacl   = .TRUE.
596          aerosol_c3h4o4 = .FALSE.
597          aerosol_nh4no3 = .FALSE.
598       ELSEIF ( aerosol_bulk == 'c3h4o4' )  THEN
599          aerosol_nacl   = .FALSE.
600          aerosol_c3h4o4 = .TRUE.
601          aerosol_nh4no3 = .FALSE.
602       ELSEIF ( aerosol_bulk == 'nh4no3' )  THEN
603          aerosol_nacl   = .FALSE.
604          aerosol_c3h4o4 = .FALSE.
605          aerosol_nh4no3 = .TRUE.
606       ELSE
607          message_string = 'unknown aerosol = "' // TRIM( aerosol_bulk ) // '"'
608          CALL message( 'check_parameters', 'PA0469', 1, 2, 0, 6, 0 )
609       ENDIF
610
611
612    END SUBROUTINE bcm_check_parameters
613
614!------------------------------------------------------------------------------!
615! Description:
616! ------------
617!> Check data output for bulk cloud module
618!------------------------------------------------------------------------------!
619    SUBROUTINE bcm_check_data_output( var, unit )
620
621       IMPLICIT NONE
622
623       CHARACTER (LEN=*) ::  unit  !<
624       CHARACTER (LEN=*) ::  var   !<
625
626       SELECT CASE ( TRIM( var ) )
627
628          CASE ( 'nc' )
629             IF ( .NOT.  microphysics_morrison )  THEN
630                message_string = 'output of "' // TRIM( var ) // '" ' //       &
631                                 'requires ' //                                &
632                                 'cloud_scheme = "morrison"'
633                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
634             ENDIF
635             unit = '1/m3'
636
637          CASE ( 'nr' )
638             IF ( .NOT.  microphysics_seifert )  THEN
639                message_string = 'output of "' // TRIM( var ) // '" ' //       &
640                                 'requires ' //                                &
641                                 'cloud_scheme = "seifert_beheng"'
642                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
643             ENDIF
644             unit = '1/m3'
645
646          CASE ( 'prr' )
647             IF ( microphysics_sat_adjust )  THEN
648                message_string = 'output of "' // TRIM( var ) // '" ' //       &
649                                 'is not available for ' //                    &
650                                 'cloud_scheme = "saturation_adjust"'
651                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
652             ENDIF
653             unit = 'kg/kg m/s'
654
655          CASE ( 'qc' )
656             unit = 'kg/kg'
657
658          CASE ( 'qr' )
659             IF ( .NOT.  microphysics_seifert ) THEN
660                message_string = 'output of "' // TRIM( var ) // '" ' //       &
661                                 'requires ' //                                &
662                                 'cloud_scheme = "seifert_beheng"'
663                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
664             ENDIF
665             unit = 'kg/kg'
666
667          CASE ( 'pra*' )
668             IF ( .NOT. microphysics_kessler  .AND.                            &
669                  .NOT. microphysics_seifert )  THEN
670                message_string = 'output of "' // TRIM( var ) // '" ' //       &
671                                 'requires ' //                                &
672                                 'cloud_scheme = "kessler" or "seifert_beheng"'
673                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
674             ENDIF
675! TODO: find sollution (maybe connected to flow_statistics redesign?)
676!              IF ( j == 1 )  THEN
677!                 message_string = 'temporal averaging of precipitation ' //     &
678!                           'amount "' // TRIM( var ) // '" is not possible'
679!                 CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
680!              ENDIF
681             unit = 'mm'
682
683          CASE ( 'prr*' )
684             IF ( .NOT. microphysics_kessler  .AND.                            &
685                  .NOT. microphysics_seifert )  THEN
686                message_string = 'output of "' // TRIM( var ) // '"' //        &
687                         ' requires' //                                        &
688                         ' cloud_scheme = "kessler" or "seifert_beheng"'
689                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
690             ENDIF
691             unit = 'mm/s'
692
693          CASE DEFAULT
694             unit = 'illegal'
695
696       END SELECT
697
698
699    END SUBROUTINE bcm_check_data_output
700
701
702!------------------------------------------------------------------------------!
703! Description:
704! ------------
705!> Check data output of profiles for bulk cloud module
706!------------------------------------------------------------------------------!
707    SUBROUTINE bcm_check_data_output_pr( variable, var_count, unit, dopr_unit )
708
709       USE arrays_3d,                                                          &
710           ONLY: zu
711
712       USE control_parameters,                                                 &
713           ONLY: data_output_pr
714
715       USE profil_parameter,                                                   &
716           ONLY: dopr_index
717
718       USE statistics,                                                         &
719           ONLY: hom, statistic_regions, pr_palm
720
721       IMPLICIT NONE
722
723       CHARACTER (LEN=*) ::  unit      !<
724       CHARACTER (LEN=*) ::  variable  !<
725       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
726
727       INTEGER(iwp) ::  var_count      !<
728       INTEGER(iwp) ::  pr_index       !<
729
730       SELECT CASE ( TRIM( variable ) )
731
732! TODO: make index generic: pr_index = pr_palm+1
733
734          CASE ( 'nc' )
735             IF ( .NOT.  microphysics_morrison )  THEN
736                message_string = 'data_output_pr = ' //                        &
737                                 TRIM( data_output_pr(var_count) ) //          &
738                                 ' is not implemented for' //                  &
739                                 ' cloud_scheme /= morrison'
740                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
741             ENDIF
742             pr_index = 123
743             dopr_index(var_count) = pr_index
744             dopr_unit     = '1/m3'
745             unit = dopr_unit
746             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
747
748          CASE ( 'nr' )
749             IF ( .NOT.  microphysics_seifert )  THEN
750                message_string = 'data_output_pr = ' //                        &
751                                 TRIM( data_output_pr(var_count) ) //          &
752                                 ' is not implemented for' //                  &
753                                 ' cloud_scheme /= seifert_beheng'
754                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
755             ENDIF
756             pr_index = 73
757             dopr_index(var_count) = pr_index
758             dopr_unit     = '1/m3'
759             unit = dopr_unit
760             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
761
762          CASE ( 'prr' )
763             IF ( microphysics_sat_adjust )  THEN
764                message_string = 'data_output_pr = ' //                        &
765                                 TRIM( data_output_pr(var_count) ) //          &
766                                 ' is not available for' //                    &
767                                 ' cloud_scheme = saturation_adjust'
768                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
769             ENDIF
770             pr_index = 76
771             dopr_index(var_count) = pr_index
772             dopr_unit     = 'kg/kg m/s'
773             unit = dopr_unit
774             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
775          CASE ( 'qc' )
776             pr_index = 75
777             dopr_index(var_count) = pr_index
778             dopr_unit     = 'kg/kg'
779             unit = dopr_unit
780             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
781
782          CASE ( 'qr' )
783             IF ( .NOT.  microphysics_seifert )  THEN
784                message_string = 'data_output_pr = ' //                        &
785                                 TRIM( data_output_pr(var_count) ) //          &
786                                 ' is not implemented for' //                  &
787                                 ' cloud_scheme /= seifert_beheng'
788                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
789             ENDIF
790             pr_index = 74
791             dopr_index(var_count) = pr_index
792             dopr_unit     = 'kg/kg'
793             unit = dopr_unit
794             hom(:,2,pr_index,:)  = SPREAD( zu, 2, statistic_regions+1 )
795
796          CASE DEFAULT
797             unit = 'illegal'
798
799       END SELECT
800
801    END SUBROUTINE bcm_check_data_output_pr
802
803
804!------------------------------------------------------------------------------!
805! Description:
806! ------------
807!> Allocate bulk cloud module arrays and define pointers
808!------------------------------------------------------------------------------!
809    SUBROUTINE bcm_init_arrays
810
811       USE indices,                                                            &
812           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
813
814
815       IMPLICIT NONE
816
817       INTEGER(iwp) ::  i !<
818       INTEGER(iwp) ::  j !<
819       INTEGER(iwp) ::  k !<
820!
821!--    Liquid water content
822       ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
823
824!
825!--    3D-cloud water content
826       IF ( .NOT. microphysics_morrison )  THEN
827          ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
828       ENDIF
829!
830!--    Precipitation amount and rate (only needed if output is switched)
831       ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg) )
832
833!
834!--    3d-precipitation rate
835       ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
836
837       IF ( microphysics_morrison )  THEN
838!
839!--       3D-cloud drop water content, cloud drop concentration arrays
840          ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
841                    nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
842                    nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
843                    qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
844                    qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
845                    qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
846       ENDIF
847
848       IF ( microphysics_seifert )  THEN
849!
850!--       3D-rain water content, rain drop concentration arrays
851          ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
852                    nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
853                    nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
854                    qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
855                    qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                    &
856                    qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
857       ENDIF
858
859!
860!--    Initial assignment of the pointers
861       ql => ql_1
862       IF ( .NOT. microphysics_morrison )  THEN
863          qc => qc_1
864       ENDIF
865       IF ( microphysics_morrison )  THEN
866          qc => qc_1;  qc_p  => qc_2;  tqc_m  => qc_3
867          nc => nc_1;  nc_p  => nc_2;  tnc_m  => nc_3
868       ENDIF
869       IF ( microphysics_seifert )  THEN
870          qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
871          nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
872       ENDIF
873
874
875    END SUBROUTINE bcm_init_arrays
876
877
878!------------------------------------------------------------------------------!
879! Description:
880! ------------
881!> Initialization of the bulk cloud module
882!------------------------------------------------------------------------------!
883    SUBROUTINE bcm_init !( dots_label, dots_unit, dots_num, dots_max )
884
885       IMPLICIT NONE
886
887       INTEGER(iwp) ::  i !<
888       INTEGER(iwp) ::  j !<
889       INTEGER(iwp) ::  k !<
890
891!        INTEGER(iwp) ::  dots_num
892!        INTEGER(iwp) ::  dots_max
893!        CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit
894!        CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label
895
896       CALL location_message( 'initializing bulk cloud module', .FALSE. )
897
898       IF ( bulk_cloud_model )  THEN
899          IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
900!
901!--          Initialize the remaining quantities
902             IF ( microphysics_morrison )  THEN
903                DO  i = nxlg, nxrg
904                   DO  j = nysg, nyng
905                      qc(:,j,i) = 0.0_wp
906                      nc(:,j,i) = 0.0_wp
907                   ENDDO
908                ENDDO
909             ENDIF
910
911             IF ( microphysics_seifert )  THEN
912                DO  i = nxlg, nxrg
913                   DO  j = nysg, nyng
914                      qr(:,j,i) = 0.0_wp
915                      nr(:,j,i) = 0.0_wp
916                   ENDDO
917                ENDDO
918             ENDIF
919!
920!--          Liquid water content and precipitation amount
921!--          are zero at beginning of the simulation
922             ql = 0.0_wp
923             qc = 0.0_wp
924             precipitation_amount = 0.0_wp
925             prr = 0.0_wp
926!
927!--          Initialize old and new time levels.
928             IF ( microphysics_morrison )  THEN
929                tqc_m = 0.0_wp
930                tnc_m = 0.0_wp
931                qc_p  = qc
932                nc_p  = nc
933             ENDIF
934             IF ( microphysics_seifert )  THEN
935                tqr_m = 0.0_wp
936                tnr_m = 0.0_wp
937                qr_p  = qr
938                nr_p  = nr
939             ENDIF
940          ENDIF ! Only if not read_restart_data
941!
942!--       constant for the sedimentation of cloud water (2-moment cloud physics)
943          sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l )                &
944                             )**( 2.0_wp / 3.0_wp ) *                             &
945                         EXP( 5.0_wp * LOG( sigma_gc )**2 )
946
947!
948!--       Calculate timestep according to precipitation
949          IF ( microphysics_seifert )  THEN
950             dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /      &
951                                w_precipitation
952          ENDIF
953
954!
955!--       Set constants for certain aerosol type
956          IF ( microphysics_morrison )  THEN
957             IF ( aerosol_nacl ) THEN
958                molecular_weight_of_solute = 0.05844_wp
959                rho_s                      = 2165.0_wp
960                vanthoff                   = 2.0_wp
961             ELSEIF ( aerosol_c3h4o4 ) THEN
962                molecular_weight_of_solute = 0.10406_wp
963                rho_s                      = 1600.0_wp
964                vanthoff                   = 1.37_wp
965             ELSEIF ( aerosol_nh4no3 ) THEN
966                molecular_weight_of_solute = 0.08004_wp
967                rho_s                      = 1720.0_wp
968                vanthoff                   = 2.31_wp
969             ENDIF
970          ENDIF
971
972!
973!--       Pre-calculate frequently calculated fractions of pi and rho_l
974          pirho_l  = pi * rho_l / 6.0_wp
975          dpirho_l = 1.0_wp / pirho_l
976
977          CALL location_message( 'finished', .TRUE. )
978
979       ELSE
980
981          CALL location_message( 'skipped', .TRUE. )
982
983       ENDIF
984
985    END SUBROUTINE bcm_init
986
987
988!------------------------------------------------------------------------------!
989! Description:
990! ------------
991!> Swapping of timelevels
992!------------------------------------------------------------------------------!
993    SUBROUTINE bcm_swap_timelevel ( mod_count )
994
995       IMPLICIT NONE
996
997       INTEGER, INTENT(IN) :: mod_count
998
999       IF ( bulk_cloud_model )  THEN
1000
1001          SELECT CASE ( mod_count )
1002
1003             CASE ( 0 )
1004
1005                IF ( microphysics_morrison )  THEN
1006                   qc => qc_1;    qc_p => qc_2
1007                   nc => nc_1;    nc_p => nc_2
1008                ENDIF
1009                IF ( microphysics_seifert )  THEN
1010                   qr => qr_1;    qr_p => qr_2
1011                   nr => nr_1;    nr_p => nr_2
1012                ENDIF
1013
1014             CASE ( 1 )
1015
1016                IF ( microphysics_morrison )  THEN
1017                   qc => qc_2;    qc_p => qc_1
1018                   nc => nc_2;    nc_p => nc_1
1019                ENDIF
1020                IF ( microphysics_seifert )  THEN
1021                   qr => qr_2;    qr_p => qr_1
1022                   nr => nr_2;    nr_p => nr_1
1023                ENDIF
1024
1025          END SELECT
1026
1027       ENDIF
1028
1029    END SUBROUTINE bcm_swap_timelevel
1030
1031
1032!------------------------------------------------------------------------------!
1033! Description:
1034! ------------
1035!> Header output for bulk cloud module
1036!------------------------------------------------------------------------------!
1037    SUBROUTINE bcm_header ( io )
1038
1039
1040       IMPLICIT NONE
1041
1042       INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
1043
1044!
1045!--    Write bulk cloud module header
1046       WRITE ( io, 1 )
1047
1048       WRITE ( io, 2 )
1049       WRITE ( io, 3 )
1050
1051       IF ( microphysics_kessler )  THEN
1052          WRITE ( io, 4 ) 'Kessler-Scheme'
1053       ENDIF
1054
1055       IF ( microphysics_seifert )  THEN
1056          WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
1057          IF ( cloud_water_sedimentation )  WRITE ( io, 5 )
1058          IF ( collision_turbulence )  WRITE ( io, 6 )
1059          IF ( ventilation_effect )  WRITE ( io, 7 )
1060          IF ( limiter_sedimentation )  WRITE ( io, 8 )
1061       ENDIF
1062
1063       WRITE ( io, 20 )
1064       WRITE ( io, 21 ) surface_pressure
1065       WRITE ( io, 22 ) r_d
1066       WRITE ( io, 23 ) rho_surface
1067       WRITE ( io, 24 ) c_p
1068       WRITE ( io, 25 ) l_v
1069
1070       IF ( microphysics_seifert )  THEN
1071          WRITE ( io, 26 ) 1.0E-6_wp * nc_const
1072          WRITE ( io, 27 ) c_sedimentation
1073       ENDIF
1074
1075
10761   FORMAT ( //' Bulk cloud module information:'/ &
1077               ' ------------------------------------------'/ )
10782   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
1079             '    total water content is used.' )
10803   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
10814   FORMAT ( '--> Precipitation parameterization via ', A )
1082
10835   FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
10846   FORMAT ( '--> Turbulence effects on precipitation process' )
10857   FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
10868   FORMAT ( '--> Slope limiter used for sedimentation process' )
1087
108820  FORMAT ( '--> Essential parameters:' )
108921  FORMAT ( '       Surface pressure             :   p_0   = ', F7.2, ' hPa')
109022  FORMAT ( '       Gas constant                 :   R     = ', F5.1, ' J/(kg K)')
109123  FORMAT ( '       Density of air               :   rho_0 = ', F6.3, ' kg/m**3')
109224  FORMAT ( '       Specific heat cap.           :   c_p   = ', F6.1, ' J/(kg K)')
109325  FORMAT ( '       Vapourization heat           :   L_v   = ', E9.2, ' J/kg')
109426  FORMAT ( '       Droplet density              :   N_c   = ', F6.1, ' 1/cm**3' )
109527  FORMAT ( '       Sedimentation Courant number :   C_s   = ', F4.1 )
1096
1097
1098    END SUBROUTINE bcm_header
1099
1100
1101!------------------------------------------------------------------------------!
1102!
1103! Description:
1104! ------------
1105!> Subroutine for averaging 3D data
1106!------------------------------------------------------------------------------!
1107    SUBROUTINE bcm_3d_data_averaging( mode, variable )
1108
1109       USE control_parameters,                                                 &
1110           ONLY:  average_count_3d
1111
1112       IMPLICIT NONE
1113
1114       CHARACTER (LEN=*) ::  mode     !<
1115       CHARACTER (LEN=*) ::  variable !<
1116
1117       INTEGER(iwp) ::  i       !< local index
1118       INTEGER(iwp) ::  j       !< local index
1119       INTEGER(iwp) ::  k       !< local index
1120
1121       IF ( mode == 'allocate' )  THEN
1122
1123          SELECT CASE ( TRIM( variable ) )
1124
1125             CASE ( 'nc' )
1126                IF ( .NOT. ALLOCATED( nc_av ) )  THEN
1127                   ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1128                ENDIF
1129                nc_av = 0.0_wp
1130
1131             CASE ( 'nr' )
1132                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
1133                   ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1134                ENDIF
1135                nr_av = 0.0_wp
1136
1137             CASE ( 'prr' )
1138                IF ( .NOT. ALLOCATED( prr_av ) )  THEN
1139                   ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1140                ENDIF
1141                prr_av = 0.0_wp
1142
1143             CASE ( 'qc' )
1144                IF ( .NOT. ALLOCATED( qc_av ) )  THEN
1145                   ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1146                ENDIF
1147                qc_av = 0.0_wp
1148
1149             CASE ( 'ql' )
1150                IF ( .NOT. ALLOCATED( ql_av ) )  THEN
1151                   ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1152                ENDIF
1153                ql_av = 0.0_wp
1154
1155             CASE ( 'qr' )
1156                IF ( .NOT. ALLOCATED( qr_av ) )  THEN
1157                   ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1158                ENDIF
1159                qr_av = 0.0_wp
1160
1161             CASE DEFAULT
1162                CONTINUE
1163
1164          END SELECT
1165
1166       ELSEIF ( mode == 'sum' )  THEN
1167
1168          SELECT CASE ( TRIM( variable ) )
1169
1170             CASE ( 'nc' )
1171                IF ( ALLOCATED( nc_av ) ) THEN
1172                   DO  i = nxlg, nxrg
1173                      DO  j = nysg, nyng
1174                         DO  k = nzb, nzt+1
1175                            nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i)
1176                         ENDDO
1177                      ENDDO
1178                   ENDDO
1179                ENDIF
1180
1181             CASE ( 'nr' )
1182                IF ( ALLOCATED( nr_av ) ) THEN
1183                   DO  i = nxlg, nxrg
1184                      DO  j = nysg, nyng
1185                         DO  k = nzb, nzt+1
1186                            nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
1187                         ENDDO
1188                      ENDDO
1189                   ENDDO
1190                ENDIF
1191
1192             CASE ( 'prr' )
1193                IF ( ALLOCATED( prr_av ) ) THEN
1194                   DO  i = nxlg, nxrg
1195                      DO  j = nysg, nyng
1196                         DO  k = nzb, nzt+1
1197                            prr_av(k,j,i) = prr_av(k,j,i) + prr(k,j,i)
1198                         ENDDO
1199                      ENDDO
1200                   ENDDO
1201                ENDIF
1202
1203             CASE ( 'qc' )
1204                IF ( ALLOCATED( qc_av ) ) THEN
1205                   DO  i = nxlg, nxrg
1206                      DO  j = nysg, nyng
1207                         DO  k = nzb, nzt+1
1208                            qc_av(k,j,i) = qc_av(k,j,i) + qc(k,j,i)
1209                         ENDDO
1210                      ENDDO
1211                   ENDDO
1212                ENDIF
1213
1214             CASE ( 'ql' )
1215                IF ( ALLOCATED( ql_av ) ) THEN
1216                   DO  i = nxlg, nxrg
1217                      DO  j = nysg, nyng
1218                         DO  k = nzb, nzt+1
1219                            ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
1220                         ENDDO
1221                      ENDDO
1222                   ENDDO
1223                ENDIF
1224
1225             CASE ( 'qr' )
1226                IF ( ALLOCATED( qr_av ) ) THEN
1227                   DO  i = nxlg, nxrg
1228                      DO  j = nysg, nyng
1229                         DO  k = nzb, nzt+1
1230                            qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i)
1231                         ENDDO
1232                      ENDDO
1233                   ENDDO
1234                ENDIF
1235
1236             CASE DEFAULT
1237                CONTINUE
1238
1239          END SELECT
1240
1241       ELSEIF ( mode == 'average' )  THEN
1242
1243          SELECT CASE ( TRIM( variable ) )
1244
1245             CASE ( 'nc' )
1246                IF ( ALLOCATED( nc_av ) ) THEN
1247                   DO  i = nxlg, nxrg
1248                      DO  j = nysg, nyng
1249                         DO  k = nzb, nzt+1
1250                            nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1251                         ENDDO
1252                      ENDDO
1253                   ENDDO
1254                ENDIF
1255
1256             CASE ( 'nr' )
1257                IF ( ALLOCATED( nr_av ) ) THEN
1258                   DO  i = nxlg, nxrg
1259                      DO  j = nysg, nyng
1260                         DO  k = nzb, nzt+1
1261                            nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1262                         ENDDO
1263                      ENDDO
1264                   ENDDO
1265                ENDIF
1266
1267             CASE ( 'prr' )
1268                IF ( ALLOCATED( prr_av ) ) THEN
1269                   DO  i = nxlg, nxrg
1270                      DO  j = nysg, nyng
1271                         DO  k = nzb, nzt+1
1272                            prr_av(k,j,i) = prr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1273                         ENDDO
1274                      ENDDO
1275                   ENDDO
1276                ENDIF
1277
1278             CASE ( 'qc' )
1279                IF ( ALLOCATED( qc_av ) ) THEN
1280                   DO  i = nxlg, nxrg
1281                      DO  j = nysg, nyng
1282                         DO  k = nzb, nzt+1
1283                            qc_av(k,j,i) = qc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1284                         ENDDO
1285                      ENDDO
1286                   ENDDO
1287                ENDIF
1288
1289             CASE ( 'ql' )
1290                IF ( ALLOCATED( ql_av ) ) THEN
1291                   DO  i = nxlg, nxrg
1292                      DO  j = nysg, nyng
1293                         DO  k = nzb, nzt+1
1294                            ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1295                         ENDDO
1296                      ENDDO
1297                   ENDDO
1298                ENDIF
1299
1300             CASE ( 'qr' )
1301                IF ( ALLOCATED( qr_av ) ) THEN
1302                   DO  i = nxlg, nxrg
1303                      DO  j = nysg, nyng
1304                         DO  k = nzb, nzt+1
1305                            qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1306                         ENDDO
1307                      ENDDO
1308                   ENDDO
1309                ENDIF
1310
1311             CASE DEFAULT
1312                CONTINUE
1313
1314          END SELECT
1315
1316       ENDIF
1317
1318    END SUBROUTINE bcm_3d_data_averaging
1319
1320
1321!------------------------------------------------------------------------------!
1322! Description:
1323! ------------
1324!> Define 2D output variables.
1325!------------------------------------------------------------------------------!
1326 SUBROUTINE bcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
1327                                two_d, nzb_do, nzt_do )
1328
1329
1330    IMPLICIT NONE
1331
1332    CHARACTER (LEN=*), INTENT(INOUT) ::  grid       !< name of vertical grid
1333    CHARACTER (LEN=*), INTENT(IN) ::  mode       !< either 'xy', 'xz' or 'yz'
1334    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
1335
1336    INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
1337    INTEGER(iwp), INTENT(IN) ::  nzb_do    !< vertical output index (bottom)
1338    INTEGER(iwp), INTENT(IN) ::  nzt_do    !< vertical output index (top)
1339
1340    INTEGER(iwp) ::  flag_nr   !< number of masking flag
1341
1342    INTEGER(iwp) ::  i         !< loop index along x-direction
1343    INTEGER(iwp) ::  j         !< loop index along y-direction
1344    INTEGER(iwp) ::  k         !< loop index along z-direction
1345
1346    LOGICAL, INTENT(INOUT) ::  found   !< flag if output variable is found
1347    LOGICAL, INTENT(INOUT) ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1348    LOGICAL ::  resorted  !< flag if output is already resorted
1349
1350    REAL(wp), PARAMETER ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
1351
1352    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf !< local
1353       !< array to which output data is resorted to
1354
1355    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
1356
1357    found = .TRUE.
1358    resorted = .FALSE.
1359!
1360!-- Set masking flag for topography for not resorted arrays
1361    flag_nr = 0    ! 0 = scalar, 1 = u, 2 = v, 3 = w
1362
1363    SELECT CASE ( TRIM( variable ) )
1364
1365       CASE ( 'nc_xy', 'nc_xz', 'nc_yz' )
1366          IF ( av == 0 )  THEN
1367             to_be_resorted => nc
1368          ELSE
1369             IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1370                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1371                nc_av = REAL( fill_value, KIND = wp )
1372             ENDIF
1373             to_be_resorted => nc_av
1374          ENDIF
1375          IF ( mode == 'xy' ) grid = 'zu'
1376
1377       CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
1378          IF ( av == 0 )  THEN
1379             to_be_resorted => nr
1380          ELSE
1381             IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1382                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1383                nr_av = REAL( fill_value, KIND = wp )
1384             ENDIF
1385             to_be_resorted => nr_av
1386          ENDIF
1387          IF ( mode == 'xy' ) grid = 'zu'
1388
1389       CASE ( 'pra*_xy' )        ! 2d-array / integral quantity => no av
1390!                CALL exchange_horiz_2d( precipitation_amount )
1391          DO  i = nxl, nxr
1392             DO  j = nys, nyn
1393                local_pf(i,j,nzb+1) =  precipitation_amount(j,i)
1394             ENDDO
1395          ENDDO
1396          precipitation_amount = 0.0_wp   ! reset for next integ. interval
1397          resorted = .TRUE.
1398          two_d = .TRUE.
1399          IF ( mode == 'xy' ) grid = 'zu1'
1400
1401       CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
1402          IF ( av == 0 )  THEN
1403!                   CALL exchange_horiz( prr, nbgp )
1404             DO  i = nxl, nxr
1405                DO  j = nys, nyn
1406                   DO  k = nzb, nzt+1
1407                      local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1)
1408                   ENDDO
1409                ENDDO
1410             ENDDO
1411          ELSE
1412             IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1413                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1414                prr_av = REAL( fill_value, KIND = wp )
1415             ENDIF
1416!                   CALL exchange_horiz( prr_av, nbgp )
1417             DO  i = nxl, nxr
1418                DO  j = nys, nyn
1419                   DO  k = nzb, nzt+1
1420                      local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1)
1421                   ENDDO
1422                ENDDO
1423             ENDDO
1424          ENDIF
1425          resorted = .TRUE.
1426          IF ( mode == 'xy' ) grid = 'zu'
1427
1428       CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
1429          IF ( av == 0 )  THEN
1430             to_be_resorted => qc
1431          ELSE
1432             IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1433                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1434                qc_av = REAL( fill_value, KIND = wp )
1435             ENDIF
1436             to_be_resorted => qc_av
1437          ENDIF
1438          IF ( mode == 'xy' ) grid = 'zu'
1439
1440       CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
1441          IF ( av == 0 )  THEN
1442             to_be_resorted => ql
1443          ELSE
1444             IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1445                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1446                ql_av = REAL( fill_value, KIND = wp )
1447             ENDIF
1448             to_be_resorted => ql_av
1449          ENDIF
1450          IF ( mode == 'xy' ) grid = 'zu'
1451
1452       CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
1453          IF ( av == 0 )  THEN
1454             to_be_resorted => qr
1455          ELSE
1456             IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1457                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1458                qr_av = REAL( fill_value, KIND = wp )
1459             ENDIF
1460             to_be_resorted => qr_av
1461          ENDIF
1462          IF ( mode == 'xy' ) grid = 'zu'
1463
1464       CASE DEFAULT
1465          found = .FALSE.
1466          grid  = 'none'
1467
1468    END SELECT
1469
1470    IF ( found .AND. .NOT. resorted )  THEN
1471       DO  i = nxl, nxr
1472          DO  j = nys, nyn
1473             DO  k = nzb_do, nzt_do
1474                local_pf(i,j,k) = MERGE( &
1475                                     to_be_resorted(k,j,i),                    &
1476                                     REAL( fill_value, KIND = wp ),            &
1477                                     BTEST( wall_flags_0(k,j,i), flag_nr )     &
1478                                  )
1479             ENDDO
1480          ENDDO
1481       ENDDO
1482    ENDIF
1483
1484 END SUBROUTINE bcm_data_output_2d
1485
1486
1487!------------------------------------------------------------------------------!
1488! Description:
1489! ------------
1490!> Define 3D output variables.
1491!------------------------------------------------------------------------------!
1492 SUBROUTINE bcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1493
1494
1495    IMPLICIT NONE
1496
1497    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
1498
1499    INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
1500    INTEGER(iwp), INTENT(IN) ::  nzb_do    !< lower limit of the data output (usually 0)
1501    INTEGER(iwp), INTENT(IN) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
1502
1503    INTEGER(iwp) ::  flag_nr   !< number of masking flag
1504
1505    INTEGER(iwp) ::  i         !< loop index along x-direction
1506    INTEGER(iwp) ::  j         !< loop index along y-direction
1507    INTEGER(iwp) ::  k         !< loop index along z-direction
1508
1509    LOGICAL, INTENT(INOUT) ::  found     !< flag if output variable is found
1510    LOGICAL ::  resorted  !< flag if output is already resorted
1511
1512    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
1513
1514    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf   !< local
1515       !< array to which output data is resorted to
1516
1517    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
1518
1519    found = .TRUE.
1520    resorted = .FALSE.
1521!
1522!-- Set masking flag for topography for not resorted arrays
1523    flag_nr = 0    ! 0 = scalar, 1 = u, 2 = v, 3 = w
1524
1525    SELECT CASE ( TRIM( variable ) )
1526
1527       CASE ( 'nc' )
1528          IF ( av == 0 )  THEN
1529             to_be_resorted => nc
1530          ELSE
1531             IF ( .NOT. ALLOCATED( nc_av ) ) THEN
1532                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1533                nc_av = REAL( fill_value, KIND = wp )
1534             ENDIF
1535             to_be_resorted => nc_av
1536          ENDIF
1537
1538       CASE ( 'nr' )
1539          IF ( av == 0 )  THEN
1540             to_be_resorted => nr
1541          ELSE
1542             IF ( .NOT. ALLOCATED( nr_av ) ) THEN
1543                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1544                nr_av = REAL( fill_value, KIND = wp )
1545             ENDIF
1546             to_be_resorted => nr_av
1547          ENDIF
1548
1549       CASE ( 'prr' )
1550          IF ( av == 0 )  THEN
1551             DO  i = nxl, nxr
1552                DO  j = nys, nyn
1553                   DO  k = nzb_do, nzt_do
1554                      local_pf(i,j,k) = prr(k,j,i)
1555                   ENDDO
1556                ENDDO
1557             ENDDO
1558          ELSE
1559             IF ( .NOT. ALLOCATED( prr_av ) ) THEN
1560                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1561                prr_av = REAL( fill_value, KIND = wp )
1562             ENDIF
1563             DO  i = nxl, nxr
1564                DO  j = nys, nyn
1565                   DO  k = nzb_do, nzt_do
1566                      local_pf(i,j,k) = prr_av(k,j,i)
1567                   ENDDO
1568                ENDDO
1569             ENDDO
1570          ENDIF
1571          resorted = .TRUE.
1572
1573       CASE ( 'qc' )
1574          IF ( av == 0 )  THEN
1575             to_be_resorted => qc
1576          ELSE
1577             IF ( .NOT. ALLOCATED( qc_av ) ) THEN
1578                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1579                qc_av = REAL( fill_value, KIND = wp )
1580             ENDIF
1581             to_be_resorted => qc_av
1582          ENDIF
1583
1584       CASE ( 'ql' )
1585          IF ( av == 0 )  THEN
1586             to_be_resorted => ql
1587          ELSE
1588             IF ( .NOT. ALLOCATED( ql_av ) ) THEN
1589                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1590                ql_av = REAL( fill_value, KIND = wp )
1591             ENDIF
1592             to_be_resorted => ql_av
1593          ENDIF
1594
1595       CASE ( 'qr' )
1596          IF ( av == 0 )  THEN
1597             to_be_resorted => qr
1598          ELSE
1599             IF ( .NOT. ALLOCATED( qr_av ) ) THEN
1600                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1601                qr_av = REAL( fill_value, KIND = wp )
1602             ENDIF
1603             to_be_resorted => qr_av
1604          ENDIF
1605
1606       CASE DEFAULT
1607          found = .FALSE.
1608
1609    END SELECT
1610
1611
1612    IF ( found .AND. .NOT. resorted )  THEN
1613       DO  i = nxl, nxr
1614          DO  j = nys, nyn
1615             DO  k = nzb_do, nzt_do
1616                local_pf(i,j,k) = MERGE(                                       &
1617                                     to_be_resorted(k,j,i),                    &
1618                                     REAL( fill_value, KIND = wp ),            &
1619                                     BTEST( wall_flags_0(k,j,i), flag_nr )     &
1620                                  )
1621             ENDDO
1622          ENDDO
1623       ENDDO
1624    ENDIF
1625
1626 END SUBROUTINE bcm_data_output_3d
1627
1628
1629!------------------------------------------------------------------------------!
1630! Description:
1631! ------------
1632!> This routine reads the respective restart data for the bulk cloud module.
1633!------------------------------------------------------------------------------!
1634    SUBROUTINE bcm_rrd_global( found )
1635
1636
1637       USE control_parameters,                                                 &
1638           ONLY: length, restart_string
1639
1640
1641       IMPLICIT NONE
1642
1643       LOGICAL, INTENT(OUT)  ::  found
1644
1645
1646       found = .TRUE.
1647
1648       SELECT CASE ( restart_string(1:length) )
1649
1650          CASE ( 'c_sedimentation' )
1651             READ ( 13 )  c_sedimentation
1652
1653          CASE ( 'bulk_cloud_model' )
1654             READ ( 13 )  bulk_cloud_model
1655
1656          CASE ( 'cloud_scheme' )
1657             READ ( 13 )  cloud_scheme
1658
1659          CASE ( 'cloud_water_sedimentation' )
1660             READ ( 13 )  cloud_water_sedimentation
1661
1662          CASE ( 'collision_turbulence' )
1663             READ ( 13 )  collision_turbulence
1664
1665          CASE ( 'limiter_sedimentation' )
1666             READ ( 13 )  limiter_sedimentation
1667
1668          CASE ( 'nc_const' )
1669             READ ( 13 )  nc_const
1670
1671          CASE ( 'precipitation' )
1672             READ ( 13 ) precipitation
1673
1674          CASE ( 'ventilation_effect' )
1675             READ ( 13 )  ventilation_effect
1676
1677          CASE ( 'na_init' )
1678             READ ( 13 )  na_init
1679
1680          CASE ( 'dry_aerosol_radius' )
1681             READ ( 13 )  dry_aerosol_radius
1682
1683          CASE ( 'sigma_bulk' )
1684             READ ( 13 )  sigma_bulk
1685
1686          CASE ( 'aerosol_bulk' )
1687             READ ( 13 )  aerosol_bulk
1688
1689          CASE ( 'curvature_solution_effects_bulk' )
1690             READ ( 13 )  curvature_solution_effects_bulk
1691
1692
1693!          CASE ( 'global_paramter' )
1694!             READ ( 13 )  global_parameter
1695!          CASE ( 'global_array' )
1696!             IF ( .NOT. ALLOCATED( global_array ) )  ALLOCATE( global_array(1:10) )
1697!             READ ( 13 )  global_array
1698
1699          CASE DEFAULT
1700
1701             found = .FALSE.
1702
1703       END SELECT
1704
1705
1706    END SUBROUTINE bcm_rrd_global
1707
1708
1709!------------------------------------------------------------------------------!
1710! Description:
1711! ------------
1712!> This routine reads the respective restart data for the bulk cloud module.
1713!------------------------------------------------------------------------------!
1714    SUBROUTINE bcm_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,       &
1715                              nxr_on_file, nynf, nync, nyn_on_file, nysf,      &
1716                              nysc, nys_on_file, tmp_2d, tmp_3d, found )
1717
1718
1719       USE control_parameters
1720
1721       USE indices
1722
1723       USE pegrid
1724
1725
1726       IMPLICIT NONE
1727
1728       INTEGER(iwp) ::  i               !<
1729       INTEGER(iwp) ::  k               !<
1730       INTEGER(iwp) ::  nxlc            !<
1731       INTEGER(iwp) ::  nxlf            !<
1732       INTEGER(iwp) ::  nxl_on_file     !<
1733       INTEGER(iwp) ::  nxrc            !<
1734       INTEGER(iwp) ::  nxrf            !<
1735       INTEGER(iwp) ::  nxr_on_file     !<
1736       INTEGER(iwp) ::  nync            !<
1737       INTEGER(iwp) ::  nynf            !<
1738       INTEGER(iwp) ::  nyn_on_file     !<
1739       INTEGER(iwp) ::  nysc            !<
1740       INTEGER(iwp) ::  nysf            !<
1741       INTEGER(iwp) ::  nys_on_file     !<
1742
1743       LOGICAL, INTENT(OUT)  ::  found
1744
1745       REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
1746       REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
1747
1748!
1749!-- Here the reading of user-defined restart data follows:
1750!-- Sample for user-defined output
1751
1752
1753       found = .TRUE.
1754
1755       SELECT CASE ( restart_string(1:length) )
1756
1757          CASE ( 'nc' )
1758             IF ( k == 1 )  READ ( 13 )  tmp_3d
1759             nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1760                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1761
1762          CASE ( 'nc_av' )
1763             IF ( .NOT. ALLOCATED( nc_av ) )  THEN
1764                ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1765             ENDIF
1766             IF ( k == 1 )  READ ( 13 )  tmp_3d
1767             nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1768                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1769
1770          CASE ( 'nr' )
1771             IF ( k == 1 )  READ ( 13 )  tmp_3d
1772             nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1773                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1774
1775          CASE ( 'nr_av' )
1776             IF ( .NOT. ALLOCATED( nr_av ) )  THEN
1777                ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1778             ENDIF
1779             IF ( k == 1 )  READ ( 13 )  tmp_3d
1780             nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1781                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1782
1783          CASE ( 'precipitation_amount' )
1784             IF ( k == 1 )  READ ( 13 )  tmp_2d
1785             precipitation_amount(nysc-nbgp:nync+nbgp,                   &
1786                                  nxlc-nbgp:nxrc+nbgp)  =                &
1787                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1788
1789          CASE ( 'prr' )
1790             IF ( .NOT. ALLOCATED( prr ) )  THEN
1791                ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1792             ENDIF
1793             IF ( k == 1 )  READ ( 13 )  tmp_3d
1794             prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =            &
1795                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1796
1797          CASE ( 'prr_av' )
1798             IF ( .NOT. ALLOCATED( prr_av ) )  THEN
1799                ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1800             ENDIF
1801             IF ( k == 1 )  READ ( 13 )  tmp_3d
1802             prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
1803                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1804
1805          CASE ( 'qc' )
1806             IF ( k == 1 )  READ ( 13 )  tmp_3d
1807             qc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1808                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1809
1810          CASE ( 'qc_av' )
1811             IF ( .NOT. ALLOCATED( qc_av ) )  THEN
1812                ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1813             ENDIF
1814             IF ( k == 1 )  READ ( 13 )  tmp_3d
1815             qc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1816                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1817
1818          CASE ( 'ql' )
1819             IF ( k == 1 )  READ ( 13 )  tmp_3d
1820             ql(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1821                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1822
1823          CASE ( 'ql_av' )
1824             IF ( .NOT. ALLOCATED( ql_av ) )  THEN
1825                ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1826             ENDIF
1827             IF ( k == 1 )  READ ( 13 )  tmp_3d
1828             ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1829                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1830
1831          CASE ( 'qr' )
1832             IF ( k == 1 )  READ ( 13 )  tmp_3d
1833             qr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =             &
1834                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1835
1836          CASE ( 'qr_av' )
1837             IF ( .NOT. ALLOCATED( qr_av ) )  THEN
1838                ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1839             ENDIF
1840             IF ( k == 1 )  READ ( 13 )  tmp_3d
1841             qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =          &
1842                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1843!
1844          CASE DEFAULT
1845
1846             found = .FALSE.
1847
1848          END SELECT
1849
1850
1851    END SUBROUTINE bcm_rrd_local
1852
1853
1854!------------------------------------------------------------------------------!
1855! Description:
1856! ------------
1857!> This routine writes the respective restart data for the bulk cloud module.
1858!------------------------------------------------------------------------------!
1859    SUBROUTINE bcm_wrd_global
1860
1861
1862       IMPLICIT NONE
1863
1864       CALL wrd_write_string( 'c_sedimentation' )
1865       WRITE ( 14 )  c_sedimentation
1866
1867       CALL wrd_write_string( 'bulk_cloud_model' )
1868       WRITE ( 14 )  bulk_cloud_model
1869
1870       CALL wrd_write_string( 'cloud_scheme' )
1871       WRITE ( 14 )  cloud_scheme
1872
1873       CALL wrd_write_string( 'cloud_water_sedimentation' )
1874       WRITE ( 14 )  cloud_water_sedimentation
1875
1876       CALL wrd_write_string( 'collision_turbulence' )
1877       WRITE ( 14 )  collision_turbulence
1878
1879       CALL wrd_write_string( 'limiter_sedimentation' )
1880       WRITE ( 14 )  limiter_sedimentation
1881
1882       CALL wrd_write_string( 'nc_const' )
1883       WRITE ( 14 )  nc_const
1884
1885       CALL wrd_write_string( 'precipitation' )
1886       WRITE ( 14 )  precipitation
1887
1888       CALL wrd_write_string( 'ventilation_effect' )
1889       WRITE ( 14 )  ventilation_effect
1890
1891       CALL wrd_write_string( 'na_init' )
1892       WRITE ( 14 )  na_init
1893
1894       CALL wrd_write_string( 'dry_aerosol_radius' )
1895       WRITE ( 14 )  dry_aerosol_radius
1896
1897       CALL wrd_write_string( 'sigma_bulk' )
1898       WRITE ( 14 )  sigma_bulk
1899
1900       CALL wrd_write_string( 'aerosol_bulk' )
1901       WRITE ( 14 )  aerosol_bulk
1902
1903       CALL wrd_write_string( 'curvature_solution_effects_bulk' )
1904       WRITE ( 14 )  curvature_solution_effects_bulk
1905
1906
1907! needs preceeding allocation if array
1908!       CALL wrd_write_string( 'global_parameter' )
1909!       WRITE ( 14 )  global_parameter
1910
1911!       IF ( ALLOCATED( inflow_damping_factor ) )  THEN
1912!          CALL wrd_write_string( 'inflow_damping_factor' )
1913!          WRITE ( 14 )  inflow_damping_factor
1914!       ENDIF
1915
1916
1917    END SUBROUTINE bcm_wrd_global
1918
1919
1920!------------------------------------------------------------------------------!
1921! Description:
1922! ------------
1923!> This routine writes the respective restart data for the bulk cloud module.
1924!------------------------------------------------------------------------------!
1925    SUBROUTINE bcm_wrd_local
1926
1927
1928       IMPLICIT NONE
1929
1930       IF ( ALLOCATED( prr ) )  THEN
1931          CALL wrd_write_string( 'prr' )
1932          WRITE ( 14 )  prr
1933       ENDIF
1934
1935       IF ( ALLOCATED( prr_av ) )  THEN
1936          CALL wrd_write_string( 'prr_av' )
1937          WRITE ( 14 )  prr_av
1938       ENDIF
1939
1940       IF ( ALLOCATED( precipitation_amount ) )  THEN
1941          CALL wrd_write_string( 'precipitation_amount' )
1942          WRITE ( 14 )  precipitation_amount
1943       ENDIF
1944
1945       CALL wrd_write_string( 'ql' )
1946       WRITE ( 14 )  ql
1947
1948       IF ( ALLOCATED( ql_av ) )  THEN
1949          CALL wrd_write_string( 'ql_av' )
1950          WRITE ( 14 )  ql_av
1951       ENDIF
1952
1953       CALL wrd_write_string( 'qc' )
1954       WRITE ( 14 )  qc
1955
1956       IF ( ALLOCATED( qc_av ) )  THEN
1957          CALL wrd_write_string( 'qc_av' )
1958          WRITE ( 14 )  qc_av
1959       ENDIF
1960
1961       IF ( microphysics_morrison )  THEN
1962
1963          CALL wrd_write_string( 'nc' )
1964          WRITE ( 14 )  nc
1965
1966          IF ( ALLOCATED( nc_av ) )  THEN
1967             CALL wrd_write_string( 'nc_av' )
1968             WRITE ( 14 )  nc_av
1969          ENDIF
1970
1971       ENDIF
1972
1973       IF ( microphysics_seifert )  THEN
1974
1975          CALL wrd_write_string( 'nr' )
1976          WRITE ( 14 )  nr
1977
1978          IF ( ALLOCATED( nr_av ) )  THEN
1979             CALL wrd_write_string( 'nr_av' )
1980             WRITE ( 14 )  nr_av
1981          ENDIF
1982
1983          CALL wrd_write_string( 'qr' )
1984          WRITE ( 14 )  qr
1985
1986          IF ( ALLOCATED( qr_av ) )  THEN
1987             CALL wrd_write_string( 'qr_av' )
1988             WRITE ( 14 )  qr_av
1989          ENDIF
1990
1991       ENDIF
1992
1993
1994    END SUBROUTINE bcm_wrd_local
1995
1996
1997!------------------------------------------------------------------------------!
1998! Description:
1999! ------------
2000!> Control of microphysics for all grid points
2001!------------------------------------------------------------------------------!
2002    SUBROUTINE bcm_actions
2003
2004       IMPLICIT NONE
2005
2006       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
2007!
2008!--       Calculate vertical profile of the hydrostatic pressure (hyp)
2009          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
2010          d_exner = exner_function_invers(hyp)
2011          exner = 1.0_wp / exner_function_invers(hyp)
2012          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
2013!
2014!--       Compute reference density
2015          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
2016       ENDIF
2017
2018!
2019!--    Compute length of time step
2020       IF ( call_microphysics_at_all_substeps )  THEN
2021          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
2022       ELSE
2023          dt_micro = dt_3d
2024       ENDIF
2025
2026!
2027!--    Reset precipitation rate
2028       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
2029
2030!
2031!--    Compute cloud physics
2032       IF ( microphysics_kessler )  THEN
2033
2034          CALL autoconversion_kessler
2035          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
2036
2037       ELSEIF ( microphysics_seifert )  THEN
2038
2039          CALL adjust_cloud
2040          IF ( microphysics_morrison )  CALL activation
2041          IF ( microphysics_morrison )  CALL condensation
2042          CALL autoconversion
2043          CALL accretion
2044          CALL selfcollection_breakup
2045          CALL evaporation_rain
2046          CALL sedimentation_rain
2047          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
2048
2049       ENDIF
2050
2051       CALL calc_precipitation_amount
2052
2053    END SUBROUTINE bcm_actions
2054
2055!------------------------------------------------------------------------------!
2056! Description:
2057! ------------
2058!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
2059!> evaporation of rain drops due to too small or too big weights
2060!> of rain drops (Stevens and Seifert, 2008).
2061!------------------------------------------------------------------------------!
2062    SUBROUTINE adjust_cloud
2063
2064       IMPLICIT NONE
2065
2066       INTEGER(iwp) ::  i                 !<
2067       INTEGER(iwp) ::  j                 !<
2068       INTEGER(iwp) ::  k                 !<
2069
2070       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2071
2072       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
2073
2074       DO  i = nxlg, nxrg
2075          DO  j = nysg, nyng
2076             DO  k = nzb+1, nzt
2077!
2078!--             Predetermine flag to mask topography
2079                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2080
2081                IF ( qr(k,j,i) <= eps_sb )  THEN
2082                   qr(k,j,i) = 0.0_wp
2083                   nr(k,j,i) = 0.0_wp
2084                ELSE
2085                   IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
2086                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
2087                   ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
2088                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
2089                   ENDIF
2090                ENDIF
2091
2092                IF ( microphysics_morrison ) THEN
2093                   IF ( qc(k,j,i) <= eps_sb )  THEN
2094                      qc(k,j,i) = 0.0_wp
2095                      nc(k,j,i) = 0.0_wp
2096                   ELSE
2097                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
2098                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * flag
2099                      ENDIF
2100                   ENDIF
2101                ENDIF
2102
2103             ENDDO
2104          ENDDO
2105       ENDDO
2106
2107       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
2108
2109    END SUBROUTINE adjust_cloud
2110
2111!------------------------------------------------------------------------------!
2112! Description:
2113! ------------
2114!> Calculate number of activated condensation nucleii after simple activation
2115!> scheme of Twomey, 1959.
2116!------------------------------------------------------------------------------!
2117    SUBROUTINE activation
2118
2119       IMPLICIT NONE
2120
2121       INTEGER(iwp) ::  i                 !<
2122       INTEGER(iwp) ::  j                 !<
2123       INTEGER(iwp) ::  k                 !<
2124
2125       REAL(wp)     ::  activ             !<
2126       REAL(wp)     ::  afactor           !<
2127       REAL(wp)     ::  beta_act          !<
2128       REAL(wp)     ::  bfactor           !<
2129       REAL(wp)     ::  k_act             !<
2130       REAL(wp)     ::  n_act             !<
2131       REAL(wp)     ::  n_ccn             !<
2132       REAL(wp)     ::  s_0               !<
2133       REAL(wp)     ::  sat_max           !<
2134       REAL(wp)     ::  sigma             !<
2135       REAL(wp)     ::  sigma_act         !<
2136
2137       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2138
2139       CALL cpu_log( log_point_s(65), 'activation', 'start' )
2140
2141       DO  i = nxlg, nxrg
2142          DO  j = nysg, nyng
2143             DO  k = nzb+1, nzt
2144!
2145!--             Predetermine flag to mask topography
2146                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2147
2148!
2149!--             Call calculation of supersaturation located in subroutine
2150                CALL supersaturation ( i, j, k )
2151!
2152!--             Prescribe parameters for activation
2153!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
2154                k_act  = 0.7_wp
2155                activ  = 0.0_wp
2156
2157
2158                IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
2159!
2160!--                Compute the number of activated Aerosols
2161!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
2162                   n_act     = na_init * sat**k_act
2163!
2164!--                Compute the number of cloud droplets
2165!--                (see: Morrison + Grabowski, 2007, JAS, 64)
2166!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
2167
2168!
2169!--                Compute activation rate after Khairoutdinov and Kogan
2170!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
2171                   sat_max = 1.0_wp / 100.0_wp
2172                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
2173                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
2174                       dt_micro
2175                ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
2176!
2177!--                Curvature effect (afactor) with surface tension
2178!--                parameterization by Straka (2009)
2179                   sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
2180                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
2181!
2182!--                Solute effect (bfactor)
2183                   bfactor = vanthoff * molecular_weight_of_water *            &
2184                       rho_s / ( molecular_weight_of_solute * rho_l )
2185
2186!
2187!--                Prescribe power index that describes the soluble fraction
2188!--                of an aerosol particle (beta)
2189!--                (see: Morrison + Grabowski, 2007, JAS, 64)
2190                   beta_act  = 0.5_wp
2191                   sigma_act = sigma_bulk**( 1.0_wp + beta_act )
2192!
2193!--                Calculate mean geometric supersaturation (s_0) with
2194!--                parameterization by Khvorostyanov and Curry (2006)
2195                   s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *    &
2196                       ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
2197
2198!
2199!--                Calculate number of activated CCN as a function of
2200!--                supersaturation and taking Koehler theory into account
2201!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
2202                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
2203                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
2204                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
2205                ENDIF
2206
2207                nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
2208
2209             ENDDO
2210          ENDDO
2211       ENDDO
2212
2213       CALL cpu_log( log_point_s(65), 'activation', 'stop' )
2214
2215    END SUBROUTINE activation
2216
2217
2218!------------------------------------------------------------------------------!
2219! Description:
2220! ------------
2221!> Calculate condensation rate for cloud water content (after Khairoutdinov and
2222!> Kogan, 2000).
2223!------------------------------------------------------------------------------!
2224    SUBROUTINE condensation
2225
2226       IMPLICIT NONE
2227
2228       INTEGER(iwp) ::  i                 !<
2229       INTEGER(iwp) ::  j                 !<
2230       INTEGER(iwp) ::  k                 !<
2231
2232       REAL(wp)     ::  cond              !<
2233       REAL(wp)     ::  cond_max          !<
2234       REAL(wp)     ::  dc                !<
2235       REAL(wp)     ::  evap              !<
2236       REAL(wp)     ::  g_fac             !<
2237       REAL(wp)     ::  nc_0              !<
2238       REAL(wp)     ::  temp              !<
2239       REAL(wp)     ::  xc                !<
2240
2241       REAL(wp) ::  flag                  !< flag to indicate first grid level above
2242
2243       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
2244
2245       DO  i = nxlg, nxrg
2246          DO  j = nysg, nyng
2247             DO  k = nzb+1, nzt
2248!
2249!--             Predetermine flag to mask topography
2250                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2251!
2252!--             Call calculation of supersaturation
2253                CALL supersaturation ( i, j, k )
2254!
2255!--             Actual temperature:
2256                temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2257
2258                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
2259                                    l_v / ( thermal_conductivity_l * temp )    &
2260                                    + r_v * temp / ( diff_coeff_l * e_s )      &
2261                                    )
2262!
2263!--             Mean weight of cloud drops
2264                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
2265                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
2266!
2267!--             Weight averaged diameter of cloud drops:
2268                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2269!
2270!--             Integral diameter of cloud drops
2271                nc_0 = nc(k,j,i) * dc
2272!
2273!--             Condensation needs only to be calculated in supersaturated regions
2274                IF ( sat > 0.0_wp )  THEN
2275!
2276!--                Condensation rate of cloud water content
2277!--                after KK scheme.
2278!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
2279                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2280                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
2281                   cond      = MIN( cond, cond_max / dt_micro )
2282
2283                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
2284                ELSEIF ( sat < 0.0_wp ) THEN
2285                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
2286                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
2287
2288                   qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
2289                ENDIF
2290                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
2291                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
2292                ENDIF
2293             ENDDO
2294          ENDDO
2295       ENDDO
2296
2297       CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
2298
2299    END SUBROUTINE condensation
2300
2301
2302!------------------------------------------------------------------------------!
2303! Description:
2304! ------------
2305!> Autoconversion rate (Seifert and Beheng, 2006).
2306!------------------------------------------------------------------------------!
2307    SUBROUTINE autoconversion
2308
2309       IMPLICIT NONE
2310
2311       INTEGER(iwp) ::  i                 !<
2312       INTEGER(iwp) ::  j                 !<
2313       INTEGER(iwp) ::  k                 !<
2314
2315       REAL(wp)     ::  alpha_cc          !<
2316       REAL(wp)     ::  autocon           !<
2317       REAL(wp)     ::  dissipation       !<
2318       REAL(wp)     ::  flag              !< flag to mask topography grid points
2319       REAL(wp)     ::  k_au              !<
2320       REAL(wp)     ::  l_mix             !<
2321       REAL(wp)     ::  nc_auto           !<
2322       REAL(wp)     ::  nu_c              !<
2323       REAL(wp)     ::  phi_au            !<
2324       REAL(wp)     ::  r_cc              !<
2325       REAL(wp)     ::  rc                !<
2326       REAL(wp)     ::  re_lambda         !<
2327       REAL(wp)     ::  sigma_cc          !<
2328       REAL(wp)     ::  tau_cloud         !<
2329       REAL(wp)     ::  xc                !<
2330
2331       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
2332
2333       DO  i = nxlg, nxrg
2334          DO  j = nysg, nyng
2335             DO  k = nzb+1, nzt
2336!
2337!--             Predetermine flag to mask topography
2338                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2339
2340                IF ( microphysics_morrison ) THEN
2341                   nc_auto = nc(k,j,i)
2342                ELSE
2343                   nc_auto = nc_const
2344                ENDIF
2345
2346                IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
2347
2348                   k_au = k_cc / ( 20.0_wp * x0 )
2349!
2350!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
2351!--                (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
2352                   tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) +         &
2353                                    qc(k,j,i) ), 0.0_wp )
2354!
2355!--                Universal function for autoconversion process
2356!--                (Seifert and Beheng, 2006):
2357                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
2358                            ( 1.0_wp - tau_cloud**0.68_wp )**3
2359!
2360!--                Shape parameter of gamma distribution (Geoffroy et al., 2010):
2361!--                (Use constant nu_c = 1.0_wp instead?)
2362                   nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
2363!
2364!--                Mean weight of cloud droplets:
2365                   xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin) 
2366!
2367!--                Parameterized turbulence effects on autoconversion (Seifert,
2368!--                Nuijens and Stevens, 2010)
2369                   IF ( collision_turbulence )  THEN
2370!
2371!--                   Weight averaged radius of cloud droplets:
2372                      rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
2373
2374                      alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
2375                      r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
2376                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
2377!
2378!--                   Mixing length (neglecting distance to ground and
2379!--                   stratification)
2380                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
2381!
2382!--                   Limit dissipation rate according to Seifert, Nuijens and
2383!--                   Stevens (2010)
2384                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
2385!
2386!--                   Compute Taylor-microscale Reynolds number:
2387                      re_lambda = 6.0_wp / 11.0_wp *                           &
2388                                  ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *   &
2389                                  SQRT( 15.0_wp / kin_vis_air ) *              &
2390                                  dissipation**( 1.0_wp / 6.0_wp )
2391!
2392!--                   The factor of 1.0E4 is needed to convert the dissipation
2393!--                   rate from m2 s-3 to cm2 s-3.
2394                      k_au = k_au * ( 1.0_wp +                                 &
2395                                      dissipation * 1.0E4_wp *                 &
2396                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
2397                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
2398                                                                      r_cc ) / &
2399                                                        sigma_cc )**2          &
2400                                                      ) + beta_cc              &
2401                                      )                                        &
2402                                    )
2403                   ENDIF
2404!
2405!--                Autoconversion rate (Seifert and Beheng, 2006):
2406                   autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /    &
2407                             ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *     &
2408                             ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &
2409                             rho_surface
2410                   autocon = MIN( autocon, qc(k,j,i) / dt_micro )
2411
2412                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
2413                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
2414                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
2415                                                              * flag
2416                   IF ( microphysics_morrison ) THEN
2417                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
2418                                    autocon / x0 * hyrho(k) * dt_micro * flag )
2419                   ENDIF
2420
2421                ENDIF
2422
2423             ENDDO
2424          ENDDO
2425       ENDDO
2426
2427       CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
2428
2429    END SUBROUTINE autoconversion
2430
2431
2432!------------------------------------------------------------------------------!
2433! Description:
2434! ------------
2435!> Autoconversion process (Kessler, 1969).
2436!------------------------------------------------------------------------------!
2437    SUBROUTINE autoconversion_kessler
2438
2439
2440       IMPLICIT NONE
2441
2442       INTEGER(iwp) ::  i      !<
2443       INTEGER(iwp) ::  j      !<
2444       INTEGER(iwp) ::  k      !<
2445       INTEGER(iwp) ::  k_wall !< topgraphy top index
2446
2447       REAL(wp)    ::  dqdt_precip !<
2448       REAL(wp)    ::  flag        !< flag to mask topography grid points
2449
2450       DO  i = nxlg, nxrg
2451          DO  j = nysg, nyng
2452!
2453!--          Determine vertical index of topography top
2454             k_wall = get_topography_top_index_ji( j, i, 's' )
2455             DO  k = nzb+1, nzt
2456!
2457!--             Predetermine flag to mask topography
2458                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2459
2460                IF ( qc(k,j,i) > ql_crit )  THEN
2461                   dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
2462                ELSE
2463                   dqdt_precip = 0.0_wp
2464                ENDIF
2465
2466                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
2467                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
2468                pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp *     &
2469                                        d_exner(k)              * flag
2470
2471!
2472!--             Compute the rain rate (stored on surface grid point)
2473                prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
2474
2475             ENDDO
2476          ENDDO
2477       ENDDO
2478
2479   END SUBROUTINE autoconversion_kessler
2480
2481
2482!------------------------------------------------------------------------------!
2483! Description:
2484! ------------
2485!> Accretion rate (Seifert and Beheng, 2006).
2486!------------------------------------------------------------------------------!
2487    SUBROUTINE accretion
2488
2489       IMPLICIT NONE
2490
2491       INTEGER(iwp) ::  i                 !<
2492       INTEGER(iwp) ::  j                 !<
2493       INTEGER(iwp) ::  k                 !<
2494
2495       REAL(wp)     ::  accr              !<
2496       REAL(wp)     ::  flag              !< flag to mask topography grid points
2497       REAL(wp)     ::  k_cr              !<
2498       REAL(wp)     ::  nc_accr           !<
2499       REAL(wp)     ::  phi_ac            !<
2500       REAL(wp)     ::  tau_cloud         !<
2501       REAL(wp)     ::  xc                !<
2502
2503
2504       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
2505
2506       DO  i = nxlg, nxrg
2507          DO  j = nysg, nyng
2508             DO  k = nzb+1, nzt
2509!
2510!--             Predetermine flag to mask topography
2511                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2512
2513                IF ( microphysics_morrison ) THEN
2514                   nc_accr = nc(k,j,i)
2515                ELSE
2516                   nc_accr = nc_const
2517                ENDIF
2518
2519                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
2520                                             .AND.  ( nc_accr > eps_mr ) ) THEN
2521!
2522!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
2523                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
2524!
2525!--                Universal function for accretion process (Seifert and
2526!--                Beheng, 2001):
2527                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
2528
2529!
2530!--                Mean weight of cloud drops
2531                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
2532!
2533!--                Parameterized turbulence effects on autoconversion (Seifert,
2534!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
2535!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
2536                   IF ( collision_turbulence )  THEN
2537                      k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                      &
2538                                       MIN( 600.0_wp,                          &
2539                                            diss(k,j,i) * 1.0E4_wp )**0.25_wp  &
2540                                     )
2541                   ELSE
2542                      k_cr = k_cr0
2543                   ENDIF
2544!
2545!--                Accretion rate (Seifert and Beheng, 2006):
2546                   accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *              &
2547                          SQRT( rho_surface * hyrho(k) )
2548                   accr = MIN( accr, qc(k,j,i) / dt_micro )
2549
2550                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
2551                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
2552                   IF ( microphysics_morrison )  THEN
2553                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i),                  &
2554                                         accr / xc * hyrho(k) * dt_micro * flag)
2555                   ENDIF
2556
2557                ENDIF
2558
2559             ENDDO
2560          ENDDO
2561       ENDDO
2562
2563       CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
2564
2565    END SUBROUTINE accretion
2566
2567
2568!------------------------------------------------------------------------------!
2569! Description:
2570! ------------
2571!> Collisional breakup rate (Seifert, 2008).
2572!------------------------------------------------------------------------------!
2573    SUBROUTINE selfcollection_breakup
2574
2575       IMPLICIT NONE
2576
2577       INTEGER(iwp) ::  i                 !<
2578       INTEGER(iwp) ::  j                 !<
2579       INTEGER(iwp) ::  k                 !<
2580
2581       REAL(wp)     ::  breakup           !<
2582       REAL(wp)     ::  dr                !<
2583       REAL(wp)     ::  flag              !< flag to mask topography grid points
2584       REAL(wp)     ::  phi_br            !<
2585       REAL(wp)     ::  selfcoll          !<
2586
2587       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
2588
2589       DO  i = nxlg, nxrg
2590          DO  j = nysg, nyng
2591             DO  k = nzb+1, nzt
2592!
2593!--             Predetermine flag to mask topography
2594                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2595
2596                IF ( qr(k,j,i) > eps_sb )  THEN
2597!
2598!--                Selfcollection rate (Seifert and Beheng, 2001):
2599                   selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) *                   &
2600                              SQRT( hyrho(k) * rho_surface )
2601!
2602!--                Weight averaged diameter of rain drops:
2603                   dr = ( hyrho(k) * qr(k,j,i) /                               &
2604                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2605!
2606!--                Collisional breakup rate (Seifert, 2008):
2607                   IF ( dr >= 0.3E-3_wp )  THEN
2608                      phi_br  = k_br * ( dr - 1.1E-3_wp )
2609                      breakup = selfcoll * ( phi_br + 1.0_wp )
2610                   ELSE
2611                      breakup = 0.0_wp
2612                   ENDIF
2613
2614                   selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
2615                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
2616
2617                ENDIF
2618             ENDDO
2619          ENDDO
2620       ENDDO
2621
2622       CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
2623
2624    END SUBROUTINE selfcollection_breakup
2625
2626
2627!------------------------------------------------------------------------------!
2628! Description:
2629! ------------
2630!> Evaporation of precipitable water. Condensation is neglected for
2631!> precipitable water.
2632!------------------------------------------------------------------------------!
2633    SUBROUTINE evaporation_rain
2634
2635       IMPLICIT NONE
2636
2637       INTEGER(iwp) ::  i                 !<
2638       INTEGER(iwp) ::  j                 !<
2639       INTEGER(iwp) ::  k                 !<
2640
2641       REAL(wp)     ::  dr                !<
2642       REAL(wp)     ::  evap              !<
2643       REAL(wp)     ::  evap_nr           !<
2644       REAL(wp)     ::  f_vent            !<
2645       REAL(wp)     ::  flag              !< flag to mask topography grid points
2646       REAL(wp)     ::  g_evap            !<
2647       REAL(wp)     ::  lambda_r          !<
2648       REAL(wp)     ::  mu_r              !<
2649       REAL(wp)     ::  mu_r_2            !<
2650       REAL(wp)     ::  mu_r_5d2          !<
2651       REAL(wp)     ::  nr_0              !<
2652       REAL(wp)     ::  temp              !<
2653       REAL(wp)     ::  xr                !<
2654
2655       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
2656
2657       DO  i = nxlg, nxrg
2658          DO  j = nysg, nyng
2659             DO  k = nzb+1, nzt
2660!
2661!--             Predetermine flag to mask topography
2662                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2663
2664                IF ( qr(k,j,i) > eps_sb )  THEN
2665
2666!
2667!--                Call calculation of supersaturation 
2668                   CALL supersaturation ( i, j, k )
2669!
2670!--                Evaporation needs only to be calculated in subsaturated regions
2671                   IF ( sat < 0.0_wp )  THEN
2672!
2673!--                   Actual temperature:
2674                      temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
2675
2676                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
2677                                          l_v / ( thermal_conductivity_l * temp ) &
2678                                          + r_v * temp / ( diff_coeff_l * e_s )   &
2679                                        )
2680!
2681!--                   Mean weight of rain drops
2682                      xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
2683!
2684!--                   Weight averaged diameter of rain drops:
2685                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
2686!
2687!--                   Compute ventilation factor and intercept parameter
2688!--                   (Seifert and Beheng, 2006; Seifert, 2008):
2689                      IF ( ventilation_effect )  THEN
2690!
2691!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
2692!--                      2005; Stevens and Seifert, 2008):
2693                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
2694                                                          ( dr - 1.4E-3_wp ) ) )
2695!
2696!--                      Slope parameter of gamma distribution (Seifert, 2008):
2697                         lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *  &
2698                                      ( mu_r + 1.0_wp )                        &
2699                                    )**( 1.0_wp / 3.0_wp ) / dr
2700
2701                         mu_r_2   = mu_r + 2.0_wp
2702                         mu_r_5d2 = mu_r + 2.5_wp
2703
2704                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
2705                                  lambda_r**( -mu_r_2 ) + b_vent *              &
2706                                  schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
2707                                  gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) *  &
2708                                  ( 1.0_wp -                                    &
2709                                    0.5_wp * ( b_term / a_term ) *              &
2710                                    ( lambda_r / ( c_term + lambda_r )          &
2711                                    )**mu_r_5d2 -                               &
2712                                    0.125_wp * ( b_term / a_term )**2 *         &
2713                                    ( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
2714                                    )**mu_r_5d2 -                               &
2715                                    0.0625_wp * ( b_term / a_term )**3 *        &
2716                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
2717                                    )**mu_r_5d2 -                               &
2718                                    0.0390625_wp * ( b_term / a_term )**4 *     &
2719                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
2720                                    )**mu_r_5d2                                 &
2721                                  )
2722
2723                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
2724                                  gamm( mu_r + 1.0_wp )
2725                      ELSE
2726                         f_vent = 1.0_wp
2727                         nr_0   = nr(k,j,i) * dr
2728                      ENDIF
2729   !
2730   !--                Evaporation rate of rain water content (Seifert and
2731   !--                Beheng, 2006):
2732                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
2733                                hyrho(k)
2734                      evap    = MAX( evap, -qr(k,j,i) / dt_micro )
2735                      evap_nr = MAX( c_evap * evap / xr * hyrho(k),            &
2736                                     -nr(k,j,i) / dt_micro )
2737
2738                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
2739                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
2740
2741                   ENDIF
2742                ENDIF
2743
2744             ENDDO
2745          ENDDO
2746       ENDDO
2747
2748       CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
2749
2750    END SUBROUTINE evaporation_rain
2751
2752
2753!------------------------------------------------------------------------------!
2754! Description:
2755! ------------
2756!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
2757!------------------------------------------------------------------------------!
2758    SUBROUTINE sedimentation_cloud
2759
2760
2761       IMPLICIT NONE
2762
2763       INTEGER(iwp) ::  i             !<
2764       INTEGER(iwp) ::  j             !<
2765       INTEGER(iwp) ::  k             !<
2766
2767       REAL(wp) ::  flag              !< flag to mask topography grid points
2768       REAL(wp) ::  nc_sedi           !<
2769
2770       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
2771       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
2772
2773
2774       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
2775
2776       sed_qc(nzt+1) = 0.0_wp
2777       sed_nc(nzt+1) = 0.0_wp
2778
2779       DO  i = nxlg, nxrg
2780          DO  j = nysg, nyng
2781             DO  k = nzt, nzb+1, -1
2782!
2783!--             Predetermine flag to mask topography
2784                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2785
2786                IF ( microphysics_morrison ) THEN
2787                   nc_sedi = nc(k,j,i)
2788                ELSE
2789                   nc_sedi = nc_const
2790                ENDIF
2791
2792!
2793!--             Sedimentation fluxes for number concentration are only calculated
2794!--             for cloud_scheme = 'morrison'
2795                IF ( microphysics_morrison ) THEN
2796                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
2797                      sed_nc(k) = sed_qc_const *                               &
2798                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
2799                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
2800                   ELSE
2801                      sed_nc(k) = 0.0_wp
2802                   ENDIF
2803
2804                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
2805                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
2806                                  ) * flag
2807
2808                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
2809                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
2810                ENDIF
2811
2812                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
2813                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
2814                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
2815                                                                           flag
2816                ELSE
2817                   sed_qc(k) = 0.0_wp
2818                ENDIF
2819
2820                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
2821                                            dt_micro + sed_qc(k+1)             &
2822                               ) * flag
2823
2824                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
2825                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
2826                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
2827                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
2828                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
2829                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
2830                                        d_exner(k) * dt_micro            * flag
2831
2832!
2833!--             Compute the precipitation rate due to cloud (fog) droplets
2834                IF ( call_microphysics_at_all_substeps )  THEN
2835                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
2836                                * weight_substep(intermediate_timestep_count)  &
2837                                * flag
2838                ELSE
2839                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
2840                ENDIF
2841
2842             ENDDO
2843          ENDDO
2844       ENDDO
2845
2846       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
2847
2848    END SUBROUTINE sedimentation_cloud
2849
2850
2851!------------------------------------------------------------------------------!
2852! Description:
2853! ------------
2854!> Computation of sedimentation flux. Implementation according to Stevens
2855!> and Seifert (2008). Code is based on UCLA-LES.
2856!------------------------------------------------------------------------------!
2857    SUBROUTINE sedimentation_rain
2858
2859       IMPLICIT NONE
2860
2861       INTEGER(iwp) ::  i             !< running index x direction
2862       INTEGER(iwp) ::  j             !< running index y direction
2863       INTEGER(iwp) ::  k             !< running index z direction
2864       INTEGER(iwp) ::  k_run         !<
2865       INTEGER(iwp) ::  m             !< running index surface elements
2866       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
2867       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
2868
2869       REAL(wp)     ::  c_run                      !<
2870       REAL(wp)     ::  d_max                      !<
2871       REAL(wp)     ::  d_mean                     !<
2872       REAL(wp)     ::  d_min                      !<
2873       REAL(wp)     ::  dr                         !<
2874       REAL(wp)     ::  flux                       !<
2875       REAL(wp)     ::  flag                       !< flag to mask topography grid points
2876       REAL(wp)     ::  lambda_r                   !<
2877       REAL(wp)     ::  mu_r                       !<
2878       REAL(wp)     ::  z_run                      !<
2879
2880       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
2881       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
2882       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
2883       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
2884       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
2885       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
2886       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
2887       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
2888
2889       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
2890
2891!
2892!--    Compute velocities
2893       DO  i = nxlg, nxrg
2894          DO  j = nysg, nyng
2895             DO  k = nzb+1, nzt
2896!
2897!--             Predetermine flag to mask topography
2898                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2899
2900                IF ( qr(k,j,i) > eps_sb )  THEN
2901!
2902!--                Weight averaged diameter of rain drops:
2903                   dr = ( hyrho(k) * qr(k,j,i) /                               &
2904                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
2905!
2906!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
2907!--                Stevens and Seifert, 2008):
2908                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
2909                                                     ( dr - 1.4E-3_wp ) ) )
2910!
2911!--                Slope parameter of gamma distribution (Seifert, 2008):
2912                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
2913                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
2914
2915                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
2916                                               a_term - b_term * ( 1.0_wp +    &
2917                                                  c_term /                     &
2918                                                  lambda_r )**( -1.0_wp *      &
2919                                                  ( mu_r + 1.0_wp ) )          &
2920                                              )                                &
2921                                ) * flag
2922
2923                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
2924                                               a_term - b_term * ( 1.0_wp +    &
2925                                                  c_term /                     &
2926                                                  lambda_r )**( -1.0_wp *      &
2927                                                  ( mu_r + 4.0_wp ) )          &
2928                                             )                                 &
2929                                ) * flag
2930                ELSE
2931                   w_nr(k) = 0.0_wp
2932                   w_qr(k) = 0.0_wp
2933                ENDIF
2934             ENDDO
2935!
2936!--          Adjust boundary values using surface data type.
2937!--          Upward-facing
2938             surf_s = bc_h(0)%start_index(j,i)
2939             surf_e = bc_h(0)%end_index(j,i)
2940             DO  m = surf_s, surf_e
2941                k         = bc_h(0)%k(m)
2942                w_nr(k-1) = w_nr(k)
2943                w_qr(k-1) = w_qr(k)
2944             ENDDO
2945!
2946!--          Downward-facing
2947             surf_s = bc_h(1)%start_index(j,i)
2948             surf_e = bc_h(1)%end_index(j,i)
2949             DO  m = surf_s, surf_e
2950                k         = bc_h(1)%k(m)
2951                w_nr(k+1) = w_nr(k)
2952                w_qr(k+1) = w_qr(k)
2953             ENDDO
2954!
2955!--          Model top boundary value
2956             w_nr(nzt+1) = 0.0_wp
2957             w_qr(nzt+1) = 0.0_wp
2958!
2959!--          Compute Courant number
2960             DO  k = nzb+1, nzt
2961!
2962!--             Predetermine flag to mask topography
2963                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2964
2965                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
2966                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
2967                          dt_micro * ddzu(k) * flag
2968                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
2969                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
2970                          dt_micro * ddzu(k) * flag
2971             ENDDO
2972!
2973!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2974             IF ( limiter_sedimentation )  THEN
2975
2976                DO k = nzb+1, nzt
2977!
2978!--                Predetermine flag to mask topography
2979                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2980
2981                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
2982                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
2983                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
2984
2985                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
2986                                                              2.0_wp * d_max,  &
2987                                                              ABS( d_mean ) )  &
2988                                                      * flag
2989
2990                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
2991                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
2992                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
2993
2994                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
2995                                                              2.0_wp * d_max,  &
2996                                                              ABS( d_mean ) )
2997                ENDDO
2998
2999             ELSE
3000
3001                nr_slope = 0.0_wp
3002                qr_slope = 0.0_wp
3003
3004             ENDIF
3005
3006             sed_nr(nzt+1) = 0.0_wp
3007             sed_qr(nzt+1) = 0.0_wp
3008!
3009!--          Compute sedimentation flux
3010             DO  k = nzt, nzb+1, -1
3011!
3012!--             Predetermine flag to mask topography
3013                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3014!
3015!--             Sum up all rain drop number densities which contribute to the flux
3016!--             through k-1/2
3017                flux  = 0.0_wp
3018                z_run = 0.0_wp ! height above z(k)
3019                k_run = k
3020                c_run = MIN( 1.0_wp, c_nr(k) )
3021                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
3022                   flux  = flux + hyrho(k_run) *                               &
3023                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
3024                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
3025                                              * flag
3026                   z_run = z_run + dzu(k_run) * flag
3027                   k_run = k_run + 1          * flag
3028                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
3029                                              * flag
3030                ENDDO
3031!
3032!--             It is not allowed to sediment more rain drop number density than
3033!--             available
3034                flux = MIN( flux,                                              &
3035                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
3036                            dt_micro                                           &
3037                          )
3038
3039                sed_nr(k) = flux / dt_micro * flag
3040                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
3041                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3042!
3043!--             Sum up all rain water content which contributes to the flux
3044!--             through k-1/2
3045                flux  = 0.0_wp
3046                z_run = 0.0_wp ! height above z(k)
3047                k_run = k
3048                c_run = MIN( 1.0_wp, c_qr(k) )
3049
3050                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
3051
3052                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
3053                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
3054                                  0.5_wp ) * c_run * dzu(k_run) * flag
3055                   z_run = z_run + dzu(k_run)                   * flag
3056                   k_run = k_run + 1                            * flag
3057                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
3058                                                                * flag
3059
3060                ENDDO
3061!
3062!--             It is not allowed to sediment more rain water content than
3063!--             available
3064                flux = MIN( flux,                                              &
3065                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
3066                            dt_micro                                           &
3067                          )
3068
3069                sed_qr(k) = flux / dt_micro * flag
3070
3071                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
3072                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3073                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
3074                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
3075                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
3076                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
3077                                        d_exner(k) * dt_micro            * flag
3078!
3079!--             Compute the rain rate
3080                IF ( call_microphysics_at_all_substeps )  THEN
3081                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
3082                                * weight_substep(intermediate_timestep_count) &
3083                                * flag
3084                ELSE
3085                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
3086                ENDIF
3087
3088             ENDDO
3089          ENDDO
3090       ENDDO
3091
3092       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
3093
3094    END SUBROUTINE sedimentation_rain
3095
3096
3097!------------------------------------------------------------------------------!
3098! Description:
3099! ------------
3100!> Computation of the precipitation amount due to gravitational settling of
3101!> rain and cloud (fog) droplets
3102!------------------------------------------------------------------------------!
3103    SUBROUTINE calc_precipitation_amount
3104
3105       IMPLICIT NONE
3106
3107       INTEGER(iwp) ::  i             !< running index x direction
3108       INTEGER(iwp) ::  j             !< running index y direction
3109       INTEGER(iwp) ::  k             !< running index y direction
3110       INTEGER(iwp) ::  m             !< running index surface elements
3111
3112       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
3113            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
3114            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
3115       THEN
3116!
3117!--       Run over all upward-facing surface elements, i.e. non-natural,
3118!--       natural and urban
3119          DO  m = 1, bc_h(0)%ns
3120             i = bc_h(0)%i(m)
3121             j = bc_h(0)%j(m)
3122             k = bc_h(0)%k(m)
3123             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
3124                                               prr(k,j,i) * hyrho(k) * dt_3d
3125          ENDDO
3126
3127       ENDIF
3128
3129    END SUBROUTINE calc_precipitation_amount
3130
3131
3132!------------------------------------------------------------------------------!
3133! Description:
3134! ------------
3135!> Control of microphysics for grid points i,j
3136!------------------------------------------------------------------------------!
3137
3138    SUBROUTINE bcm_actions_ij( i, j )
3139
3140       IMPLICIT NONE
3141
3142       INTEGER(iwp) ::  i                 !<
3143       INTEGER(iwp) ::  j                 !<
3144
3145       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
3146!
3147!--       Calculate vertical profile of the hydrostatic pressure (hyp)
3148          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
3149          d_exner = exner_function_invers(hyp)
3150          exner = 1.0_wp / exner_function_invers(hyp)
3151          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
3152!
3153!--       Compute reference density
3154          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
3155       ENDIF
3156
3157!
3158!--    Compute length of time step
3159       IF ( call_microphysics_at_all_substeps )  THEN
3160          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
3161       ELSE
3162          dt_micro = dt_3d
3163       ENDIF
3164!
3165!--    Reset precipitation rate
3166       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
3167
3168!
3169!--    Compute cloud physics
3170       IF( microphysics_kessler )  THEN
3171
3172          CALL autoconversion_kessler( i,j )
3173          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
3174
3175       ELSEIF ( microphysics_seifert )  THEN
3176
3177          CALL adjust_cloud( i,j )
3178          IF ( microphysics_morrison ) CALL activation( i,j )
3179          IF ( microphysics_morrison ) CALL condensation( i,j )
3180          CALL autoconversion( i,j )
3181          CALL accretion( i,j )
3182          CALL selfcollection_breakup( i,j )
3183          CALL evaporation_rain( i,j )
3184          CALL sedimentation_rain( i,j )
3185          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
3186
3187       ENDIF
3188
3189       CALL calc_precipitation_amount( i,j )
3190
3191    END SUBROUTINE bcm_actions_ij
3192
3193!------------------------------------------------------------------------------!
3194! Description:
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!> The same procedure is applied to cloud droplets if they are determined
3200!> prognostically. Call for grid point i,j
3201!------------------------------------------------------------------------------!
3202    SUBROUTINE adjust_cloud_ij( i, j )
3203
3204       IMPLICIT NONE
3205
3206       INTEGER(iwp) ::  i                 !<
3207       INTEGER(iwp) ::  j                 !<
3208       INTEGER(iwp) ::  k                 !<
3209
3210       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
3211
3212       DO  k = nzb+1, nzt
3213!
3214!--       Predetermine flag to mask topography
3215          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3216
3217          IF ( qr(k,j,i) <= eps_sb )  THEN
3218             qr(k,j,i) = 0.0_wp
3219             nr(k,j,i) = 0.0_wp
3220          ELSE
3221!
3222!--          Adjust number of raindrops to avoid nonlinear effects in
3223!--          sedimentation and evaporation of rain drops due to too small or
3224!--          too big weights of rain drops (Stevens and Seifert, 2008).
3225             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
3226                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
3227             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
3228                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
3229             ENDIF
3230
3231          ENDIF
3232
3233          IF ( microphysics_morrison ) THEN
3234             IF ( qc(k,j,i) <= eps_sb )  THEN
3235                qc(k,j,i) = 0.0_wp
3236                nc(k,j,i) = 0.0_wp
3237             ELSE
3238                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
3239                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
3240                ENDIF
3241             ENDIF
3242          ENDIF
3243
3244       ENDDO
3245
3246    END SUBROUTINE adjust_cloud_ij
3247
3248!------------------------------------------------------------------------------!
3249! Description:
3250! ------------
3251!> Calculate number of activated condensation nucleii after simple activation
3252!> scheme of Twomey, 1959.
3253!------------------------------------------------------------------------------!
3254    SUBROUTINE activation_ij( i, j )
3255
3256       IMPLICIT NONE
3257
3258       INTEGER(iwp) ::  i                 !<
3259       INTEGER(iwp) ::  j                 !<
3260       INTEGER(iwp) ::  k                 !<
3261
3262       REAL(wp)     ::  activ             !<
3263       REAL(wp)     ::  afactor           !<
3264       REAL(wp)     ::  beta_act          !<
3265       REAL(wp)     ::  bfactor           !<
3266       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3267       REAL(wp)     ::  k_act             !<
3268       REAL(wp)     ::  n_act             !<
3269       REAL(wp)     ::  n_ccn             !<
3270       REAL(wp)     ::  s_0               !<
3271       REAL(wp)     ::  sat_max           !<
3272       REAL(wp)     ::  sigma             !<
3273       REAL(wp)     ::  sigma_act         !<
3274
3275       DO  k = nzb+1, nzt
3276!
3277!--       Predetermine flag to mask topography
3278          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3279!
3280!--       Call calculation of supersaturation 
3281          CALL supersaturation ( i, j, k )
3282!
3283!--       Prescribe parameters for activation
3284!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
3285          k_act  = 0.7_wp
3286          activ  = 0.0_wp
3287
3288          IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
3289!
3290!--          Compute the number of activated Aerosols
3291!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
3292             n_act     = na_init * sat**k_act
3293!
3294!--          Compute the number of cloud droplets
3295!--          (see: Morrison + Grabowski, 2007, JAS, 64)
3296!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
3297
3298!
3299!--          Compute activation rate after Khairoutdinov and Kogan
3300!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
3301             sat_max = 0.8_wp / 100.0_wp
3302             activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
3303                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
3304                  dt_micro
3305
3306             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
3307          ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
3308!
3309!--          Curvature effect (afactor) with surface tension
3310!--          parameterization by Straka (2009)
3311             sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
3312             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
3313!
3314!--          Solute effect (bfactor)
3315             bfactor = vanthoff * molecular_weight_of_water *                  &
3316                  rho_s / ( molecular_weight_of_solute * rho_l )
3317
3318!
3319!--          Prescribe power index that describes the soluble fraction
3320!--          of an aerosol particle (beta).
3321!--          (see: Morrison + Grabowski, 2007, JAS, 64)
3322             beta_act  = 0.5_wp
3323             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
3324!
3325!--          Calculate mean geometric supersaturation (s_0) with
3326!--          parameterization by Khvorostyanov and Curry (2006)
3327             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
3328               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
3329
3330!
3331!--          Calculate number of activated CCN as a function of
3332!--          supersaturation and taking Koehler theory into account
3333!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
3334             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
3335                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
3336             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
3337
3338             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
3339          ENDIF
3340
3341       ENDDO
3342
3343    END SUBROUTINE activation_ij
3344
3345!------------------------------------------------------------------------------!
3346! Description:
3347! ------------
3348!> Calculate condensation rate for cloud water content (after Khairoutdinov and
3349!> Kogan, 2000).
3350!------------------------------------------------------------------------------!
3351    SUBROUTINE condensation_ij( i, j )
3352
3353       IMPLICIT NONE
3354
3355       INTEGER(iwp) ::  i                 !<
3356       INTEGER(iwp) ::  j                 !<
3357       INTEGER(iwp) ::  k                 !<
3358
3359       REAL(wp)     ::  cond              !<
3360       REAL(wp)     ::  cond_max          !<
3361       REAL(wp)     ::  dc                !<
3362       REAL(wp)     ::  evap              !<
3363       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3364       REAL(wp)     ::  g_fac             !<
3365       REAL(wp)     ::  nc_0              !<
3366       REAL(wp)     ::  temp              !<
3367       REAL(wp)     ::  xc                !<
3368
3369
3370       DO  k = nzb+1, nzt
3371!
3372!--       Predetermine flag to mask topography
3373          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3374!
3375!--       Call calculation of supersaturation 
3376          CALL supersaturation ( i, j, k )
3377!
3378!--       Actual temperature:
3379          temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3380
3381          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
3382                              l_v / ( thermal_conductivity_l * temp )    &
3383                              + r_v * temp / ( diff_coeff_l * e_s )      &
3384                            )
3385!
3386!--       Mean weight of cloud drops
3387          IF ( nc(k,j,i) <= 0.0_wp) CYCLE
3388          xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
3389!
3390!--       Weight averaged diameter of cloud drops:
3391          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3392!
3393!--       Integral diameter of cloud drops
3394          nc_0 = nc(k,j,i) * dc
3395!
3396!--       Condensation needs only to be calculated in supersaturated regions
3397          IF ( sat > 0.0_wp )  THEN
3398!
3399!--          Condensation rate of cloud water content
3400!--          after KK scheme.
3401!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
3402             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3403             cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
3404             cond      = MIN( cond, cond_max / dt_micro )
3405
3406             qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
3407          ELSEIF ( sat < 0.0_wp ) THEN
3408             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
3409             evap      = MAX( evap, -qc(k,j,i) / dt_micro )
3410
3411             qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
3412          ENDIF
3413       ENDDO
3414
3415    END SUBROUTINE condensation_ij
3416
3417
3418!------------------------------------------------------------------------------!
3419! Description:
3420! ------------
3421!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
3422!------------------------------------------------------------------------------!
3423    SUBROUTINE autoconversion_ij( i, j )
3424
3425       IMPLICIT NONE
3426
3427       INTEGER(iwp) ::  i                 !<
3428       INTEGER(iwp) ::  j                 !<
3429       INTEGER(iwp) ::  k                 !<
3430
3431       REAL(wp)     ::  alpha_cc          !<
3432       REAL(wp)     ::  autocon           !<
3433       REAL(wp)     ::  dissipation       !<
3434       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3435       REAL(wp)     ::  k_au              !<
3436       REAL(wp)     ::  l_mix             !<
3437       REAL(wp)     ::  nc_auto           !<
3438       REAL(wp)     ::  nu_c              !<
3439       REAL(wp)     ::  phi_au            !<
3440       REAL(wp)     ::  r_cc              !<
3441       REAL(wp)     ::  rc                !<
3442       REAL(wp)     ::  re_lambda         !<
3443       REAL(wp)     ::  sigma_cc          !<
3444       REAL(wp)     ::  tau_cloud         !<
3445       REAL(wp)     ::  xc                !<
3446
3447       DO  k = nzb+1, nzt
3448!
3449!--       Predetermine flag to mask topography
3450          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3451          IF ( microphysics_morrison ) THEN
3452             nc_auto = nc(k,j,i)
3453          ELSE
3454             nc_auto = nc_const
3455          ENDIF
3456
3457          IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
3458
3459             k_au = k_cc / ( 20.0_wp * x0 )
3460!
3461!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
3462!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
3463             tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
3464                              0.0_wp )
3465!
3466!--          Universal function for autoconversion process
3467!--          (Seifert and Beheng, 2006):
3468             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
3469!
3470!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
3471!--          (Use constant nu_c = 1.0_wp instead?)
3472             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
3473!
3474!--          Mean weight of cloud droplets:
3475             xc = hyrho(k) * qc(k,j,i) / nc_auto
3476!
3477!--          Parameterized turbulence effects on autoconversion (Seifert,
3478!--          Nuijens and Stevens, 2010)
3479             IF ( collision_turbulence )  THEN
3480!
3481!--             Weight averaged radius of cloud droplets:
3482                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
3483
3484                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
3485                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
3486                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
3487!
3488!--             Mixing length (neglecting distance to ground and stratification)
3489                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
3490!
3491!--             Limit dissipation rate according to Seifert, Nuijens and
3492!--             Stevens (2010)
3493                dissipation = MIN( 0.06_wp, diss(k,j,i) )
3494!
3495!--             Compute Taylor-microscale Reynolds number:
3496                re_lambda = 6.0_wp / 11.0_wp *                                 &
3497                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
3498                            SQRT( 15.0_wp / kin_vis_air ) *                    &
3499                            dissipation**( 1.0_wp / 6.0_wp )
3500!
3501!--             The factor of 1.0E4 is needed to convert the dissipation rate
3502!--             from m2 s-3 to cm2 s-3.
3503                k_au = k_au * ( 1.0_wp +                                       &
3504                                dissipation * 1.0E4_wp *                       &
3505                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
3506                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
3507                                                  sigma_cc )**2                &
3508                                                ) + beta_cc                    &
3509                                )                                              &
3510                              )
3511             ENDIF
3512!
3513!--          Autoconversion rate (Seifert and Beheng, 2006):
3514             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
3515                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
3516                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
3517                       rho_surface
3518             autocon = MIN( autocon, qc(k,j,i) / dt_micro )
3519
3520             qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
3521             qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
3522             nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
3523             IF ( microphysics_morrison ) THEN
3524                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
3525                                    autocon / x0 * hyrho(k) * dt_micro * flag )
3526             ENDIF
3527
3528          ENDIF
3529
3530       ENDDO
3531
3532    END SUBROUTINE autoconversion_ij
3533
3534!------------------------------------------------------------------------------!
3535! Description:
3536! ------------
3537!> Autoconversion process (Kessler, 1969).
3538!------------------------------------------------------------------------------!
3539    SUBROUTINE autoconversion_kessler_ij( i, j )
3540
3541
3542       IMPLICIT NONE
3543
3544       INTEGER(iwp) ::  i      !<
3545       INTEGER(iwp) ::  j      !<
3546       INTEGER(iwp) ::  k      !<
3547       INTEGER(iwp) ::  k_wall !< topography top index
3548
3549       REAL(wp)    ::  dqdt_precip !<
3550       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
3551
3552!
3553!--    Determine vertical index of topography top
3554       k_wall = get_topography_top_index_ji( j, i, 's' )
3555       DO  k = nzb+1, nzt
3556!
3557!--       Predetermine flag to mask topography
3558          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3559
3560          IF ( qc(k,j,i) > ql_crit )  THEN
3561             dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
3562          ELSE
3563             dqdt_precip = 0.0_wp
3564          ENDIF
3565
3566          qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
3567          q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
3568          pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
3569                                                                 * flag
3570
3571!
3572!--       Compute the rain rate (stored on surface grid point)
3573          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
3574
3575       ENDDO
3576
3577    END SUBROUTINE autoconversion_kessler_ij
3578
3579!------------------------------------------------------------------------------!
3580! Description:
3581! ------------
3582!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
3583!------------------------------------------------------------------------------!
3584    SUBROUTINE accretion_ij( i, j )
3585
3586       IMPLICIT NONE
3587
3588       INTEGER(iwp) ::  i                 !<
3589       INTEGER(iwp) ::  j                 !<
3590       INTEGER(iwp) ::  k                 !<
3591
3592       REAL(wp)     ::  accr              !<
3593       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3594       REAL(wp)     ::  k_cr              !<
3595       REAL(wp)     ::  nc_accr           !<
3596       REAL(wp)     ::  phi_ac            !<
3597       REAL(wp)     ::  tau_cloud         !<
3598       REAL(wp)     ::  xc                !<
3599
3600
3601       DO  k = nzb+1, nzt
3602!
3603!--       Predetermine flag to mask topography
3604          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3605          IF ( microphysics_morrison ) THEN
3606             nc_accr = nc(k,j,i)
3607          ELSE
3608             nc_accr = nc_const
3609          ENDIF
3610
3611          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
3612               ( nc_accr > eps_mr ) )  THEN
3613!
3614!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
3615             tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
3616!
3617!--          Universal function for accretion process
3618!--          (Seifert and Beheng, 2001):
3619             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
3620
3621!
3622!--          Mean weight of cloud drops
3623             xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
3624!
3625!--          Parameterized turbulence effects on autoconversion (Seifert,
3626!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
3627!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
3628             IF ( collision_turbulence )  THEN
3629                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
3630                                 MIN( 600.0_wp,                                &
3631                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
3632                               )
3633             ELSE
3634                k_cr = k_cr0
3635             ENDIF
3636!
3637!--          Accretion rate (Seifert and Beheng, 2006):
3638             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
3639                    SQRT( rho_surface * hyrho(k) )
3640             accr = MIN( accr, qc(k,j,i) / dt_micro )
3641
3642             qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
3643             qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
3644             IF ( microphysics_morrison )  THEN
3645                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
3646                                             hyrho(k) * dt_micro * flag        &
3647                                         )
3648             ENDIF
3649
3650
3651          ENDIF
3652
3653       ENDDO
3654
3655    END SUBROUTINE accretion_ij
3656
3657
3658!------------------------------------------------------------------------------!
3659! Description:
3660! ------------
3661!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
3662!------------------------------------------------------------------------------!
3663    SUBROUTINE selfcollection_breakup_ij( i, j )
3664
3665       IMPLICIT NONE
3666
3667       INTEGER(iwp) ::  i                 !<
3668       INTEGER(iwp) ::  j                 !<
3669       INTEGER(iwp) ::  k                 !<
3670
3671       REAL(wp)     ::  breakup           !<
3672       REAL(wp)     ::  dr                !<
3673       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3674       REAL(wp)     ::  phi_br            !<
3675       REAL(wp)     ::  selfcoll          !<
3676
3677       DO  k = nzb+1, nzt
3678!
3679!--       Predetermine flag to mask topography
3680          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3681
3682          IF ( qr(k,j,i) > eps_sb )  THEN
3683!
3684!--          Selfcollection rate (Seifert and Beheng, 2001):
3685             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
3686!
3687!--          Weight averaged diameter of rain drops:
3688             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
3689!
3690!--          Collisional breakup rate (Seifert, 2008):
3691             IF ( dr >= 0.3E-3_wp )  THEN
3692                phi_br  = k_br * ( dr - 1.1E-3_wp )
3693                breakup = selfcoll * ( phi_br + 1.0_wp )
3694             ELSE
3695                breakup = 0.0_wp
3696             ENDIF
3697
3698             selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
3699             nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
3700
3701          ENDIF
3702       ENDDO
3703
3704    END SUBROUTINE selfcollection_breakup_ij
3705
3706
3707!------------------------------------------------------------------------------!
3708! Description:
3709! ------------
3710!> Evaporation of precipitable water. Condensation is neglected for
3711!> precipitable water. Call for grid point i,j
3712!------------------------------------------------------------------------------!
3713    SUBROUTINE evaporation_rain_ij( i, j )
3714
3715       IMPLICIT NONE
3716
3717       INTEGER(iwp) ::  i                 !<
3718       INTEGER(iwp) ::  j                 !<
3719       INTEGER(iwp) ::  k                 !<
3720
3721       REAL(wp)     ::  dr                !<
3722       REAL(wp)     ::  evap              !<
3723       REAL(wp)     ::  evap_nr           !<
3724       REAL(wp)     ::  f_vent            !<
3725       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
3726       REAL(wp)     ::  g_evap            !<
3727       REAL(wp)     ::  lambda_r          !<
3728       REAL(wp)     ::  mu_r              !<
3729       REAL(wp)     ::  mu_r_2            !<
3730       REAL(wp)     ::  mu_r_5d2          !<
3731       REAL(wp)     ::  nr_0              !<
3732       REAL(wp)     ::  temp              !<
3733       REAL(wp)     ::  xr                !<
3734
3735       DO  k = nzb+1, nzt
3736!
3737!--       Predetermine flag to mask topography
3738          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3739
3740          IF ( qr(k,j,i) > eps_sb )  THEN
3741!
3742!--          Call calculation of supersaturation 
3743             CALL supersaturation ( i, j, k )
3744!
3745!--          Evaporation needs only to be calculated in subsaturated regions
3746             IF ( sat < 0.0_wp )  THEN
3747!
3748!--             Actual temperature:
3749                temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
3750
3751                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
3752                                    ( thermal_conductivity_l * temp ) +        &
3753                                    r_v * temp / ( diff_coeff_l * e_s )        &
3754                                  )
3755!
3756!--             Mean weight of rain drops
3757                xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
3758!
3759!--             Weight averaged diameter of rain drops:
3760                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
3761!
3762!--             Compute ventilation factor and intercept parameter
3763!--             (Seifert and Beheng, 2006; Seifert, 2008):
3764                IF ( ventilation_effect )  THEN
3765!
3766!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
3767!--                Stevens and Seifert, 2008):
3768                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
3769!
3770!--                Slope parameter of gamma distribution (Seifert, 2008):
3771                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
3772                                ( mu_r + 1.0_wp )                              &
3773                              )**( 1.0_wp / 3.0_wp ) / dr
3774
3775                   mu_r_2   = mu_r + 2.0_wp
3776                   mu_r_5d2 = mu_r + 2.5_wp
3777
3778                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
3779                            b_vent * schmidt_p_1d3 *                           &
3780                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
3781                            lambda_r**( -mu_r_5d2 ) *                          &
3782                            ( 1.0_wp -                                         &
3783                              0.5_wp * ( b_term / a_term ) *                   &
3784                              ( lambda_r / ( c_term + lambda_r )               &
3785                              )**mu_r_5d2 -                                    &
3786                              0.125_wp * ( b_term / a_term )**2 *              &
3787                              ( lambda_r / ( 2.0_wp * c_term + lambda_r )      &
3788                              )**mu_r_5d2 -                                    &
3789                              0.0625_wp * ( b_term / a_term )**3 *             &
3790                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
3791                              )**mu_r_5d2 -                                    &
3792                              0.0390625_wp * ( b_term / a_term )**4 *          &
3793                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
3794                              )**mu_r_5d2                                      &
3795                            )
3796
3797                   nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /           &
3798                            gamm( mu_r + 1.0_wp )
3799                ELSE
3800                   f_vent = 1.0_wp
3801                   nr_0   = nr(k,j,i) * dr
3802                ENDIF
3803!
3804!--             Evaporation rate of rain water content (Seifert and Beheng, 2006):
3805                evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
3806                evap    = MAX( evap, -qr(k,j,i) / dt_micro )
3807                evap_nr = MAX( c_evap * evap / xr * hyrho(k),                  &
3808                               -nr(k,j,i) / dt_micro )
3809
3810                qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
3811                nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
3812
3813             ENDIF
3814          ENDIF
3815
3816       ENDDO
3817
3818    END SUBROUTINE evaporation_rain_ij
3819
3820
3821!------------------------------------------------------------------------------!
3822! Description:
3823! ------------
3824!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
3825!> Call for grid point i,j
3826!------------------------------------------------------------------------------!
3827    SUBROUTINE sedimentation_cloud_ij( i, j )
3828
3829       IMPLICIT NONE
3830
3831       INTEGER(iwp) ::  i             !<
3832       INTEGER(iwp) ::  j             !<
3833       INTEGER(iwp) ::  k             !<
3834
3835       REAL(wp)     ::  flag    !< flag to indicate first grid level above surface
3836       REAL(wp)     ::  nc_sedi !<
3837
3838       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc  !<
3839       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
3840
3841       sed_qc(nzt+1) = 0.0_wp
3842       sed_nc(nzt+1) = 0.0_wp
3843
3844
3845       DO  k = nzt, nzb+1, -1
3846!
3847!--       Predetermine flag to mask topography
3848          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3849          IF ( microphysics_morrison ) THEN
3850             nc_sedi = nc(k,j,i)
3851          ELSE
3852             nc_sedi = nc_const
3853          ENDIF
3854!
3855!--       Sedimentation fluxes for number concentration are only calculated
3856!--       for cloud_scheme = 'morrison'
3857          IF ( microphysics_morrison ) THEN
3858             IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
3859                sed_nc(k) = sed_qc_const *                                     &
3860                            ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
3861                            ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
3862             ELSE
3863                sed_nc(k) = 0.0_wp
3864             ENDIF
3865
3866             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
3867                              nc(k,j,i) / dt_micro + sed_nc(k+1)                &
3868                            ) * flag
3869
3870             nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *               &
3871                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
3872          ENDIF
3873
3874          IF ( qc(k,j,i) > eps_sb  .AND.  nc_sedi > eps_mr )  THEN
3875             sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *        &
3876                         ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
3877          ELSE
3878             sed_qc(k) = 0.0_wp
3879          ENDIF
3880
3881          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /          &
3882                                      dt_micro + sed_qc(k+1)                   &
3883                         ) * flag
3884
3885          q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3886                                hyrho(k) * dt_micro * flag
3887          qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3888                                hyrho(k) * dt_micro * flag
3889          pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
3890                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
3891                                                                * flag
3892
3893!
3894!--       Compute the precipitation rate of cloud (fog) droplets
3895          IF ( call_microphysics_at_all_substeps )  THEN
3896             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *                  &
3897                              weight_substep(intermediate_timestep_count) * flag
3898          ELSE
3899             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
3900          ENDIF
3901
3902       ENDDO
3903
3904    END SUBROUTINE sedimentation_cloud_ij
3905
3906
3907!------------------------------------------------------------------------------!
3908! Description:
3909! ------------
3910!> Computation of sedimentation flux. Implementation according to Stevens
3911!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
3912!------------------------------------------------------------------------------!
3913    SUBROUTINE sedimentation_rain_ij( i, j )
3914
3915       IMPLICIT NONE
3916
3917       INTEGER(iwp) ::  i             !< running index x direction
3918       INTEGER(iwp) ::  j             !< running index y direction
3919       INTEGER(iwp) ::  k             !< running index z direction
3920       INTEGER(iwp) ::  k_run         !<
3921       INTEGER(iwp) ::  m             !< running index surface elements
3922       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
3923       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
3924
3925       REAL(wp)     ::  c_run                      !<
3926       REAL(wp)     ::  d_max                      !<
3927       REAL(wp)     ::  d_mean                     !<
3928       REAL(wp)     ::  d_min                      !<
3929       REAL(wp)     ::  dr                         !<
3930       REAL(wp)     ::  flux                       !<
3931       REAL(wp)     ::  flag                       !< flag to indicate first grid level above surface
3932       REAL(wp)     ::  lambda_r                   !<
3933       REAL(wp)     ::  mu_r                       !<
3934       REAL(wp)     ::  z_run                      !<
3935
3936       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
3937       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
3938       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
3939       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
3940       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
3941       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
3942       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
3943       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
3944
3945!
3946!--    Compute velocities
3947       DO  k = nzb+1, nzt
3948!
3949!--       Predetermine flag to mask topography
3950          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3951
3952          IF ( qr(k,j,i) > eps_sb )  THEN
3953!
3954!--          Weight averaged diameter of rain drops:
3955             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
3956!
3957!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
3958!--          Stevens and Seifert, 2008):
3959             mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
3960!
3961!--          Slope parameter of gamma distribution (Seifert, 2008):
3962             lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *              &
3963                          ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
3964
3965             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
3966                                         a_term - b_term * ( 1.0_wp +          &
3967                                            c_term / lambda_r )**( -1.0_wp *   &
3968                                            ( mu_r + 1.0_wp ) )                &
3969                                        )                                      &
3970                          ) * flag
3971             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
3972                                         a_term - b_term * ( 1.0_wp +          &
3973                                            c_term / lambda_r )**( -1.0_wp *   &
3974                                            ( mu_r + 4.0_wp ) )                &
3975                                       )                                       &
3976                          ) * flag
3977          ELSE
3978             w_nr(k) = 0.0_wp
3979             w_qr(k) = 0.0_wp
3980          ENDIF
3981       ENDDO
3982!
3983!--    Adjust boundary values using surface data type.
3984!--    Upward facing non-natural
3985       surf_s = bc_h(0)%start_index(j,i)
3986       surf_e = bc_h(0)%end_index(j,i)
3987       DO  m = surf_s, surf_e
3988          k         = bc_h(0)%k(m)
3989          w_nr(k-1) = w_nr(k)
3990          w_qr(k-1) = w_qr(k)
3991       ENDDO
3992!
3993!--    Downward facing non-natural
3994       surf_s = bc_h(1)%start_index(j,i)
3995       surf_e = bc_h(1)%end_index(j,i)
3996       DO  m = surf_s, surf_e
3997          k         = bc_h(1)%k(m)
3998          w_nr(k+1) = w_nr(k)
3999          w_qr(k+1) = w_qr(k)
4000       ENDDO
4001!
4002!--    Neumann boundary condition at model top
4003       w_nr(nzt+1) = 0.0_wp
4004       w_qr(nzt+1) = 0.0_wp
4005!
4006!--    Compute Courant number
4007       DO  k = nzb+1, nzt
4008!
4009!--       Predetermine flag to mask topography
4010          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4011
4012          c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
4013                    dt_micro * ddzu(k) * flag
4014          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
4015                    dt_micro * ddzu(k) * flag
4016       ENDDO
4017!
4018!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
4019       IF ( limiter_sedimentation )  THEN
4020
4021          DO k = nzb+1, nzt
4022!
4023!--          Predetermine flag to mask topography
4024             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4025
4026             d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
4027             d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
4028             d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
4029
4030             qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
4031                                                        2.0_wp * d_max,        &
4032                                                        ABS( d_mean ) ) * flag
4033
4034             d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
4035             d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
4036             d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
4037
4038             nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
4039                                                        2.0_wp * d_max,        &
4040                                                        ABS( d_mean ) ) * flag
4041          ENDDO
4042
4043       ELSE
4044
4045          nr_slope = 0.0_wp
4046          qr_slope = 0.0_wp
4047
4048       ENDIF
4049
4050       sed_nr(nzt+1) = 0.0_wp
4051       sed_qr(nzt+1) = 0.0_wp
4052!
4053!--    Compute sedimentation flux
4054       DO  k = nzt, nzb+1, -1
4055!
4056!--       Predetermine flag to mask topography
4057          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4058!
4059!--       Sum up all rain drop number densities which contribute to the flux
4060!--       through k-1/2
4061          flux  = 0.0_wp
4062          z_run = 0.0_wp ! height above z(k)
4063          k_run = k
4064          c_run = MIN( 1.0_wp, c_nr(k) )
4065          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
4066             flux  = flux + hyrho(k_run) *                                     &
4067                     ( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0_wp - c_run ) *  &
4068                     0.5_wp ) * c_run * dzu(k_run) * flag
4069             z_run = z_run + dzu(k_run)            * flag
4070             k_run = k_run + 1                     * flag
4071             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
4072          ENDDO
4073!
4074!--       It is not allowed to sediment more rain drop number density than
4075!--       available
4076          flux = MIN( flux,                                                    &
4077                      hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * dt_micro )
4078
4079          sed_nr(k) = flux / dt_micro * flag
4080          nr(k,j,i)  = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
4081                                    hyrho(k) * dt_micro * flag
4082!
4083!--       Sum up all rain water content which contributes to the flux
4084!--       through k-1/2
4085          flux  = 0.0_wp
4086          z_run = 0.0_wp ! height above z(k)
4087          k_run = k
4088          c_run = MIN( 1.0_wp, c_qr(k) )
4089
4090          DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
4091
4092             flux  = flux + hyrho(k_run) *                                     &
4093                     ( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0_wp - c_run ) *  &
4094                     0.5_wp ) * c_run * dzu(k_run) * flag
4095             z_run = z_run + dzu(k_run)            * flag
4096             k_run = k_run + 1                     * flag
4097             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
4098
4099          ENDDO
4100!
4101!--       It is not allowed to sediment more rain water content than available
4102          flux = MIN( flux,                                                    &
4103                      hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_micro )
4104
4105          sed_qr(k) = flux / dt_micro * flag
4106
4107          qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4108                                hyrho(k) * dt_micro * flag
4109          q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4110                                hyrho(k) * dt_micro * flag
4111          pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
4112                                hyrho(k) * lv_d_cp * d_exner(k) * dt_micro     &
4113                                                                * flag
4114!
4115!--       Compute the rain rate
4116          IF ( call_microphysics_at_all_substeps )  THEN
4117             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
4118                          * weight_substep(intermediate_timestep_count) * flag
4119          ELSE
4120             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
4121          ENDIF
4122
4123       ENDDO
4124
4125    END SUBROUTINE sedimentation_rain_ij
4126
4127
4128!------------------------------------------------------------------------------!
4129! Description:
4130! ------------
4131!> This subroutine computes the precipitation amount due to gravitational
4132!> settling of rain and cloud (fog) droplets
4133!------------------------------------------------------------------------------!
4134    SUBROUTINE calc_precipitation_amount_ij( i, j )
4135
4136       IMPLICIT NONE
4137
4138       INTEGER(iwp) ::  i             !< running index x direction
4139       INTEGER(iwp) ::  j             !< running index y direction
4140       INTEGER(iwp) ::  k             !< running index z direction
4141       INTEGER(iwp) ::  m             !< running index surface elements
4142       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
4143       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
4144
4145       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
4146            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
4147            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
4148       THEN
4149
4150          surf_s = bc_h(0)%start_index(j,i)
4151          surf_e = bc_h(0)%end_index(j,i)
4152          DO  m = surf_s, surf_e
4153             k                         = bc_h(0)%k(m)
4154             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
4155                                               prr(k,j,i) * hyrho(k) * dt_3d
4156          ENDDO
4157
4158       ENDIF
4159
4160    END SUBROUTINE calc_precipitation_amount_ij
4161
4162
4163!------------------------------------------------------------------------------!
4164! Description:
4165! ------------
4166!> Computation of the diagnostic supersaturation sat, actual liquid water
4167!< temperature t_l and saturation water vapor mixing ratio q_s
4168!------------------------------------------------------------------------------!
4169    SUBROUTINE supersaturation ( i,j,k )
4170
4171       IMPLICIT NONE
4172
4173       INTEGER(iwp) ::  i   !< running index
4174       INTEGER(iwp) ::  j   !< running index
4175       INTEGER(iwp) ::  k   !< running index
4176
4177       REAL(wp) ::  alpha   !< correction factor
4178!
4179!--    Actual liquid water temperature:
4180       t_l = exner(k) * pt(k,j,i)
4181!
4182!--    Calculate water vapor saturation pressure
4183       e_s = magnus( t_l )
4184!
4185!--    Computation of saturation mixing ratio:
4186       q_s   = rd_d_rv * e_s / ( hyp(k) - e_s )
4187!
4188!--    Correction factor
4189       alpha = rd_d_rv * lv_d_rd * lv_d_cp / ( t_l * t_l )
4190!
4191!--    Correction of the approximated value
4192!--    (see: Cuijpers + Duynkerke, 1993, JAS, 23)
4193       q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s )
4194!
4195!--    Supersaturation:
4196!--    Not in case of microphysics_kessler or microphysics_sat_adjust
4197!--    since qr is unallocated
4198       IF ( .NOT. microphysics_kessler  .AND.                                  &
4199            .NOT.  microphysics_sat_adjust )  THEN
4200          sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
4201       ENDIF
4202
4203    END SUBROUTINE supersaturation
4204
4205
4206!------------------------------------------------------------------------------!
4207! Description:
4208! ------------
4209!> Calculation of the liquid water content (0%-or-100%-scheme). This scheme is
4210!> used by the one and the two moment cloud physics scheme. Using the two moment
4211!> scheme, this calculation results in the cloud water content.
4212!------------------------------------------------------------------------------!
4213    SUBROUTINE calc_liquid_water_content
4214
4215
4216
4217       IMPLICIT NONE
4218
4219       INTEGER(iwp) ::  i !<
4220       INTEGER(iwp) ::  j !<
4221       INTEGER(iwp) ::  k !<
4222
4223       REAL(wp)     ::  flag !< flag to indicate first grid level above surface
4224
4225
4226       DO  i = nxlg, nxrg
4227          DO  j = nysg, nyng
4228             DO  k = nzb+1, nzt
4229!
4230!--             Predetermine flag to mask topography
4231                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4232
4233   !
4234   !--          Call calculation of supersaturation located
4235                CALL supersaturation( i, j, k )
4236
4237   !
4238   !--          Compute the liquid water content
4239                IF ( microphysics_seifert  .AND. .NOT. microphysics_morrison ) &
4240                THEN
4241                   IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
4242                      qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag
4243                      ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag
4244                   ELSE
4245                      IF ( q(k,j,i) < qr(k,j,i) )  q(k,j,i) = qr(k,j,i)
4246                      qc(k,j,i) = 0.0_wp
4247                      ql(k,j,i) = qr(k,j,i) * flag
4248                   ENDIF
4249                ELSEIF ( microphysics_morrison )  THEN
4250                   ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag
4251                ELSE
4252                   IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
4253                      qc(k,j,i) = ( q(k,j,i) - q_s ) * flag
4254                      ql(k,j,i) = qc(k,j,i) * flag
4255                   ELSE
4256                      qc(k,j,i) = 0.0_wp
4257                      ql(k,j,i) = 0.0_wp
4258                   ENDIF
4259                ENDIF
4260             ENDDO
4261          ENDDO
4262       ENDDO
4263
4264    END SUBROUTINE calc_liquid_water_content
4265
4266!------------------------------------------------------------------------------!
4267! Description:
4268! ------------
4269!> This function computes the gamma function (Press et al., 1992).
4270!> The gamma function is needed for the calculation of the evaporation
4271!> of rain drops.
4272!------------------------------------------------------------------------------!
4273    FUNCTION gamm( xx )
4274
4275       IMPLICIT NONE
4276
4277       INTEGER(iwp) ::  j            !<
4278
4279       REAL(wp)     ::  gamm         !<
4280       REAL(wp)     ::  ser          !<
4281       REAL(wp)     ::  tmp          !<
4282       REAL(wp)     ::  x_gamm       !<
4283       REAL(wp)     ::  xx           !<
4284       REAL(wp)     ::  y_gamm       !<
4285
4286
4287       REAL(wp), PARAMETER  ::  stp = 2.5066282746310005_wp               !<
4288       REAL(wp), PARAMETER  ::  cof(6) = (/ 76.18009172947146_wp,      &
4289                                           -86.50532032941677_wp,      &
4290                                            24.01409824083091_wp,      &
4291                                            -1.231739572450155_wp,     &
4292                                             0.1208650973866179E-2_wp, &
4293                                            -0.5395239384953E-5_wp /)     !<
4294
4295       x_gamm = xx
4296       y_gamm = x_gamm
4297       tmp = x_gamm + 5.5_wp
4298       tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
4299       ser = 1.000000000190015_wp
4300
4301       DO  j = 1, 6
4302          y_gamm = y_gamm + 1.0_wp
4303          ser    = ser + cof( j ) / y_gamm
4304       ENDDO
4305
4306!
4307!--    Until this point the algorithm computes the logarithm of the gamma
4308!--    function. Hence, the exponential function is used.
4309!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
4310       gamm = EXP( tmp ) * stp * ser / x_gamm
4311
4312       RETURN
4313
4314    END FUNCTION gamm
4315
4316 END MODULE bulk_cloud_model_mod
Note: See TracBrowser for help on using the repository browser.