Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |