Home

Context Navigation

bulk_cloud_model_mod.f90 @ 4268

Last change on this file since 4268 was 4268, checked in by schwenkel, 4 years ago
Introducing module interface for boundary conditions and move module specific boundary conditions into their modules
Property svn:keywords set to `Id`
File size: 202.0 KB

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