Home

Context Navigation

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

Last change on this file since 3875 was 3874, checked in by knoop, 5 years ago
Implemented non_transport_physics module interfaces
Property svn:keywords set to `Id`
File size: 192.4 KB

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