Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |