Home

Context Navigation

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

Last change on this file since 3700 was 3700, checked in by knoop, 6 years ago
Moved user_define_netdf_grid into user_module.f90 Added module interface for the definition of additional timeseries
Property svn:keywords set to `Id`
File size: 160.5 KB

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