Home

Context Navigation

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

Last change on this file since 3814 was 3786, checked in by raasch, 6 years ago
unused variables removed, interoperable C datatypes introduced in particle type to avoid compiler warnings
Property svn:keywords set to `Id`
File size: 160.6 KB

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