Home

Context Navigation

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

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