Home

Context Navigation

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

Last change on this file since 3424 was 3383, checked in by knoop, 6 years ago
renamed string "microphysics module" to "bulk cloud module"
Property svn:keywords set to `Id`
File size: 163.9 KB

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