Home

Context Navigation

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

Last change on this file since 3445 was 3445, checked in by schwenkel, 5 years ago
Minor bugfix and use of subroutine for supersaturation calculation in case of cache version
Property svn:keywords set to `Id`
File size: 161.3 KB

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