Home

Context Navigation

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

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