Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |