Home

Context Navigation

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

Last change on this file since 3766 was 3724, checked in by kanani, 6 years ago
Correct double-used log_point_s units (bulk_cloud_model_mod, time_integration, turbulence_closure_mod)
Property svn:keywords set to `Id`
File size: 160.6 KB

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