Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |