Home

Context Navigation

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

Last change on this file since 3452 was 3452, checked in by schwenkel, 5 years ago
Bugfix for profiles output
Property svn:keywords set to `Id`
File size: 161.4 KB

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