Home

Context Navigation

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

Last change on this file since 3869 was 3869, checked in by knoop, 5 years ago
moving the furniture in the BCM around ;-)
Property svn:keywords set to `Id`
File size: 159.0 KB

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