Home

Context Navigation

source: palm/trunk/SOURCE/microphysics_mod.f90 @ 2579

Last change on this file since 2579 was 2522, checked in by schwenkel, 7 years ago
minor bugfix
Property svn:keywords set to `Id`
File size: 124.6 KB

Rev	Line
[1850]	1	!> @file microphysics_mod.f90
[2000]	2	!------------------------------------------------------------------------------!
[1093]	3	! This file is part of PALM.
	4	!
[2000]	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.
[1093]	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	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1093]	19	!
[1000]	20	! Current revisions:
[1092]	21	! ------------------
[2318]	22	!
	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: microphysics_mod.f90 2522 2017-10-05 14:20:37Z knoop $
[2522]	27	! Minor bugfix
	28	!
	29	! 2375 2017-08-29 14:10:28Z schwenkel
[2375]	30	! Improved aerosol initilization and some minor bugfixes
	31	! for droplet sedimenation
	32	!
	33	! 2318 2017-07-20 17:27:44Z suehring
[2318]	34	! Get topography top index via Function call
	35	!
	36	! 2317 2017-07-20 17:27:19Z suehring
[2312]	37	! s1 changed to log_sigma
	38	!
	39	! 2292 2017-06-20 09:51:42Z schwenkel
	40	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
	41	! includes two more prognostic equations for cloud drop concentration (nc)
	42	! and cloud water content (qc).
	43	! - The process of activation is parameterized with a simple Twomey
	44	! activion scheme or with considering solution and curvature
[2292]	45	! effects (Khvorostyanov and Curry ,2006).
	46	! - The saturation adjustment scheme is replaced by the parameterization
	47	! of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
	48	! - All other microphysical processes of Seifert and Beheng are used.
[2312]	49	! Additionally, in those processes the reduction of cloud number concentration
	50	! is considered.
	51	!
[2292]	52	! 2233 2017-05-30 18:08:54Z suehring
[1321]	53	!
[2233]	54	! 2232 2017-05-30 17:47:52Z suehring
	55	! Adjustments to new topography and surface concept
[2312]	56	!
[2156]	57	! 2155 2017-02-21 09:57:40Z hoffmann
	58	! Bugfix in the calculation of microphysical quantities on ghost points.
	59	!
[2032]	60	! 2031 2016-10-21 15:11:58Z knoop
	61	! renamed variable rho to rho_ocean
[2155]	62	!
[2001]	63	! 2000 2016-08-20 18:09:15Z knoop
	64	! Forced header and separation lines into 80 columns
[2155]	65	!
[1851]	66	! 1850 2016-04-08 13:29:27Z maronga
	67	! Module renamed
	68	! Adapted for modularization of microphysics.
	69	!
[1846]	70	! 1845 2016-04-08 08:29:13Z raasch
	71	! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
	72	! point (instead of one point above surface)
	73	!
[1832]	74	! 1831 2016-04-07 13:15:51Z hoffmann
	75	! turbulence renamed collision_turbulence,
	76	! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
	77	! avaialble for microphysics_kessler.
	78	!
[1823]	79	! 1822 2016-04-07 07:49:42Z hoffmann
	80	! Unused variables removed.
	81	! Kessler scheme integrated.
	82	!
[1692]	83	! 1691 2015-10-26 16:17:44Z maronga
	84	! Added new routine calc_precipitation_amount. The routine now allows to account
	85	! for precipitation due to sedimenation of cloud (fog) droplets
[2155]	86	!
[1683]	87	! 1682 2015-10-07 23:56:08Z knoop
[2155]	88	! Code annotations made doxygen readable
[1683]	89	!
[1647]	90	! 1646 2015-09-02 16:00:10Z hoffmann
	91	! Bugfix: Wrong computation of d_mean.
	92	!
[1362]	93	! 1361 2014-04-16 15:17:48Z hoffmann
	94	! Bugfix in sedimentation_rain: Index corrected.
	95	! Vectorized version of adjust_cloud added.
	96	! Little reformatting of the code.
[2155]	97	!
[1354]	98	! 1353 2014-04-08 15:21:23Z heinze
	99	! REAL constants provided with KIND-attribute
[2155]	100	!
[1347]	101	! 1346 2014-03-27 13:18:20Z heinze
[2155]	102	! Bugfix: REAL constants provided with KIND-attribute especially in call of
[1347]	103	! intrinsic function like MAX, MIN, SIGN
[2155]	104	!
[1335]	105	! 1334 2014-03-25 12:21:40Z heinze
	106	! Bugfix: REAL constants provided with KIND-attribute
	107	!
[1323]	108	! 1322 2014-03-20 16:38:49Z raasch
	109	! REAL constants defined as wp-kind
	110	!
[1321]	111	! 1320 2014-03-20 08:40:49Z raasch
[1320]	112	! ONLY-attribute added to USE-statements,
	113	! kind-parameters added to all INTEGER and REAL declaration statements,
	114	! kinds are defined in new module kinds,
	115	! comment fields (!:) to be used for variable explanations added to
	116	! all variable declaration statements
[1000]	117	!
[1242]	118	! 1241 2013-10-30 11:36:58Z heinze
[2031]	119	! hyp and rho_ocean have to be calculated at each time step if data from external
[1242]	120	! file LSF_DATA are used
	121	!
[1116]	122	! 1115 2013-03-26 18:16:16Z hoffmann
	123	! microphyical tendencies are calculated in microphysics_control in an optimized
	124	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
	125	!
[1107]	126	! 1106 2013-03-04 05:31:38Z raasch
	127	! small changes in code formatting
	128	!
[1093]	129	! 1092 2013-02-02 11:24:22Z raasch
	130	! unused variables removed
	131	! file put under GPL
	132	!
[1066]	133	! 1065 2012-11-22 17:42:36Z hoffmann
	134	! Sedimentation process implemented according to Stevens and Seifert (2008).
[1115]	135	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
[1066]	136	! and Stevens, 2010).
	137	!
[1054]	138	! 1053 2012-11-13 17:11:03Z hoffmann
	139	! initial revision
[2155]	140	!
[1000]	141	! Description:
	142	! ------------
[1849]	143	!> Calculate bilk cloud microphysics.
[1000]	144	!------------------------------------------------------------------------------!
[1682]	145	MODULE microphysics_mod
[1000]	146
[1849]	147	USE kinds
	148
	149	IMPLICIT NONE
	150
[2292]	151	LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation
	152	LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
	153	LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter
	154	LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects
	155	LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect
[1849]	156
	157	REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3)
	158	REAL(wp) :: a_2 = -7.38E-5_wp !< coef. in turb. parametrization (cm-2 s3)
	159	REAL(wp) :: a_3 = -1.40E-2_wp !< coef. in turb. parametrization
	160	REAL(wp) :: a_term = 9.65_wp !< coef. for terminal velocity (m s-1)
	161	REAL(wp) :: a_vent = 0.78_wp !< coef. for ventilation effect
	162	REAL(wp) :: b_1 = 11.45E-6_wp !< coef. in turb. parametrization (m)
	163	REAL(wp) :: b_2 = 9.68E-6_wp !< coef. in turb. parametrization (m)
	164	REAL(wp) :: b_3 = 0.62_wp !< coef. in turb. parametrization
	165	REAL(wp) :: b_term = 9.8_wp !< coef. for terminal velocity (m s-1)
	166	REAL(wp) :: b_vent = 0.308_wp !< coef. for ventilation effect
	167	REAL(wp) :: beta_cc = 3.09E-4_wp !< coef. in turb. parametrization (cm-2 s3)
	168	REAL(wp) :: c_1 = 4.82E-6_wp !< coef. in turb. parametrization (m)
	169	REAL(wp) :: c_2 = 4.8E-6_wp !< coef. in turb. parametrization (m)
	170	REAL(wp) :: c_3 = 0.76_wp !< coef. in turb. parametrization
	171	REAL(wp) :: c_const = 0.93_wp !< const. in Taylor-microscale Reynolds number
	172	REAL(wp) :: c_evap = 0.7_wp !< constant in evaporation
	173	REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1)
	174	REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1)
[2155]	175	REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme
[2292]	176	REAL(wp) :: eps_mr = 0.0_wp !< threshold for morrison scheme
[1849]	177	REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1)
	178	REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1)
	179	REAL(wp) :: k_rr = 7.12_wp !< const. rain-rain kernel (m3 kg-1 s-1)
	180	REAL(wp) :: k_br = 1000.0_wp !< const. in breakup parametrization (m-1)
	181	REAL(wp) :: k_st = 1.2E8_wp !< const. in drizzle parametrization (m-1 s-1)
	182	REAL(wp) :: kappa_rr = 60.7_wp !< const. in collision kernel (kg-1/3)
	183	REAL(wp) :: kin_vis_air = 1.4086E-5_wp !< kin. viscosity of air (m2 s-1)
	184	REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1)
	185	REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1)
	186	REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number0.33333, 0.710.33333
	187	REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets
	188	REAL(wp) :: thermal_conductivity_l = 2.43E-2_wp !< therm. cond. air (J m-1 s-1 K-1)
	189	REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1)
	190	REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg)
[2292]	191	REAL(wp) :: xamin = 5.24E-19_wp !< average aerosol mass (kg) (~ 0.05Âµm)
	192	REAL(wp) :: xcmin = 4.18E-15_wp !< minimum cloud drop size (kg) (~ 1Âµm)
	193	REAL(wp) :: xcmax = 2.6E-10_wp !< maximum cloud drop size (kg) (~ 40Âµm)
[1849]	194	REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg)
	195	REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg)
	196
[2375]	197	REAL(wp) :: c_sedimentation = 2.0_wp !< Courant number of sedimentation process
	198	REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l )
	199	REAL(wp) :: dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius
	200	REAL(wp) :: dt_micro !< microphysics time step
	201	REAL(wp) :: sigma_bulk = 2.0_wp !< width of aerosol spectrum
	202	REAL(wp) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3)
	203	REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration
	204	REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s)
	205	REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water
	206	REAL(wp) :: pirho_l !< pi * rho_l / 6.0;
[1849]	207
	208	REAL(wp), DIMENSION(:), ALLOCATABLE :: nc_1d !<
	209	REAL(wp), DIMENSION(:), ALLOCATABLE :: nr_1d !<
	210	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_1d !<
	211	REAL(wp), DIMENSION(:), ALLOCATABLE :: qc_1d !<
	212	REAL(wp), DIMENSION(:), ALLOCATABLE :: qr_1d !<
	213	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_1d !<
	214
	215	SAVE
	216
[1000]	217	PRIVATE
[1849]	218	PUBLIC microphysics_control, microphysics_init
[1000]	219
[2292]	220	PUBLIC cloud_water_sedimentation, collision_turbulence, &
[2375]	221	curvature_solution_effects_bulk, c_sedimentation, &
	222	dry_aerosol_radius, dt_precipitation, &
	223	limiter_sedimentation, na_init, nc_const, sigma_gc, sigma_bulk, &
[1849]	224	ventilation_effect
	225
[2312]	226
[1115]	227	INTERFACE microphysics_control
	228	MODULE PROCEDURE microphysics_control
	229	MODULE PROCEDURE microphysics_control_ij
	230	END INTERFACE microphysics_control
[1022]	231
[1115]	232	INTERFACE adjust_cloud
	233	MODULE PROCEDURE adjust_cloud
	234	MODULE PROCEDURE adjust_cloud_ij
	235	END INTERFACE adjust_cloud
	236
[2292]	237	INTERFACE activation
	238	MODULE PROCEDURE activation
	239	MODULE PROCEDURE activation_ij
	240	END INTERFACE activation
	241
	242	INTERFACE condensation
	243	MODULE PROCEDURE condensation
	244	MODULE PROCEDURE condensation_ij
	245	END INTERFACE condensation
	246
[1000]	247	INTERFACE autoconversion
	248	MODULE PROCEDURE autoconversion
	249	MODULE PROCEDURE autoconversion_ij
	250	END INTERFACE autoconversion
	251
[1822]	252	INTERFACE autoconversion_kessler
	253	MODULE PROCEDURE autoconversion_kessler
	254	MODULE PROCEDURE autoconversion_kessler_ij
	255	END INTERFACE autoconversion_kessler
	256
[1000]	257	INTERFACE accretion
	258	MODULE PROCEDURE accretion
	259	MODULE PROCEDURE accretion_ij
	260	END INTERFACE accretion
[1005]	261
	262	INTERFACE selfcollection_breakup
	263	MODULE PROCEDURE selfcollection_breakup
	264	MODULE PROCEDURE selfcollection_breakup_ij
	265	END INTERFACE selfcollection_breakup
[1012]	266
	267	INTERFACE evaporation_rain
	268	MODULE PROCEDURE evaporation_rain
	269	MODULE PROCEDURE evaporation_rain_ij
	270	END INTERFACE evaporation_rain
	271
	272	INTERFACE sedimentation_cloud
	273	MODULE PROCEDURE sedimentation_cloud
	274	MODULE PROCEDURE sedimentation_cloud_ij
	275	END INTERFACE sedimentation_cloud
[2155]	276
[1012]	277	INTERFACE sedimentation_rain
	278	MODULE PROCEDURE sedimentation_rain
	279	MODULE PROCEDURE sedimentation_rain_ij
	280	END INTERFACE sedimentation_rain
	281
[1691]	282	INTERFACE calc_precipitation_amount
	283	MODULE PROCEDURE calc_precipitation_amount
	284	MODULE PROCEDURE calc_precipitation_amount_ij
	285	END INTERFACE calc_precipitation_amount
	286
[1000]	287	CONTAINS
[1849]	288	!------------------------------------------------------------------------------!
	289	! Description:
	290	! ------------
	291	!> Initialization of bulk microphysics
	292	!------------------------------------------------------------------------------!
	293	SUBROUTINE microphysics_init
[1000]	294
[1849]	295	USE arrays_3d, &
	296	ONLY: dzu
[1000]	297
[1849]	298	USE constants, &
	299	ONLY: pi
	300
	301	USE cloud_parameters, &
[2375]	302	ONLY: molecular_weight_of_solute, rho_l, rho_s, vanthoff
[1849]	303
	304	USE control_parameters, &
[2375]	305	ONLY: aerosol_nacl, aerosol_c3h4o4 , aerosol_nh4no3, &
	306	microphysics_morrison, microphysics_seifert
[1849]	307
	308	USE indices, &
	309	ONLY: nzb, nzt
	310
	311	IMPLICIT NONE
	312
	313	!
	314	!-- constant for the sedimentation of cloud water (2-moment cloud physics)
	315	sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) &
	316	)*( 2.0_wp / 3.0_wp ) &
	317	EXP( 5.0_wp * LOG( sigma_gc )**2 )
	318
	319	!
	320	!-- Calculate timestep according to precipitation
	321	IF ( microphysics_seifert ) THEN
	322	dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / &
	323	w_precipitation
	324	ENDIF
	325
	326	!
[2375]	327	!-- Set constants for certain aerosol type
	328	IF ( microphysics_morrison ) THEN
	329	IF ( aerosol_nacl ) THEN
	330	molecular_weight_of_solute = 0.05844_wp
	331	rho_s = 2165.0_wp
	332	vanthoff = 2.0_wp
	333	ELSEIF ( aerosol_c3h4o4 ) THEN
	334	molecular_weight_of_solute = 0.10406_wp
	335	rho_s = 1600.0_wp
	336	vanthoff = 1.37_wp
	337	ELSEIF ( aerosol_nh4no3 ) THEN
	338	molecular_weight_of_solute = 0.08004_wp
	339	rho_s = 1720.0_wp
	340	vanthoff = 2.31_wp
	341	ENDIF
	342	ENDIF
	343
	344	!
[1849]	345	!-- Pre-calculate frequently calculated fractions of pi and rho_l
	346	pirho_l = pi * rho_l / 6.0_wp
	347	dpirho_l = 1.0_wp / pirho_l
	348
	349	!
	350	!-- Allocate 1D microphysics arrays
[2292]	351	ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1), &
[1849]	352	qc_1d(nzb:nzt+1) )
	353
[2292]	354	IF ( microphysics_morrison .OR. microphysics_seifert ) THEN
	355	ALLOCATE ( nc_1d(nzb:nzt+1) )
	356	ENDIF
	357
[1849]	358	IF ( microphysics_seifert ) THEN
	359	ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
	360	ENDIF
	361
	362	END SUBROUTINE microphysics_init
	363
	364
[1000]	365	!------------------------------------------------------------------------------!
[1682]	366	! Description:
	367	! ------------
[1849]	368	!> Control of microphysics for all grid points
[1000]	369	!------------------------------------------------------------------------------!
[1115]	370	SUBROUTINE microphysics_control
[1022]	371
[1361]	372	USE arrays_3d, &
[1849]	373	ONLY: hyp, pt_init, prr, zu
[1361]	374
	375	USE cloud_parameters, &
[1849]	376	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
[1361]	377
	378	USE control_parameters, &
[1849]	379	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
	380	intermediate_timestep_count, large_scale_forcing, &
[2292]	381	lsf_surf, microphysics_kessler, microphysics_morrison, &
	382	microphysics_seifert, pt_surface, &
	383	rho_surface,surface_pressure
[1361]	384
	385	USE indices, &
	386	ONLY: nzb, nzt
	387
[1320]	388	USE kinds
[1115]	389
[1361]	390	USE statistics, &
	391	ONLY: weight_pres
	392
[1115]	393	IMPLICIT NONE
	394
[1682]	395	INTEGER(iwp) :: k !<
[1115]	396
[1682]	397	REAL(wp) :: t_surface !<
[1361]	398
	399	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	400	!
	401	!-- Calculate:
	402	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	403	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	404	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
	405	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
	406	DO k = nzb, nzt+1
	407	hyp(k) = surface_pressure * 100.0_wp * &
	408	( ( t_surface - g / cp * zu(k) ) / &
	409	t_surface )**(1.0_wp / 0.286_wp)
	410	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	411	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]	412	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1115]	413	ENDDO
[1822]	414
[1361]	415	!
	416	!-- Compute reference density
	417	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
	418	ENDIF
[1115]	419
[1361]	420	!
[2155]	421	!-- Compute length of time step
[1361]	422	IF ( call_microphysics_at_all_substeps ) THEN
	423	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	424	ELSE
	425	dt_micro = dt_3d
	426	ENDIF
	427
	428	!
[1822]	429	!-- Reset precipitation rate
	430	IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp
	431
	432	!
[1361]	433	!-- Compute cloud physics
[1822]	434	IF ( microphysics_kessler ) THEN
	435
	436	CALL autoconversion_kessler
[1831]	437	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
[1822]	438
	439	ELSEIF ( microphysics_seifert ) THEN
	440
[1361]	441	CALL adjust_cloud
[2292]	442	IF ( microphysics_morrison ) CALL activation
	443	IF ( microphysics_morrison ) CALL condensation
[1361]	444	CALL autoconversion
	445	CALL accretion
	446	CALL selfcollection_breakup
	447	CALL evaporation_rain
	448	CALL sedimentation_rain
[1831]	449	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
[1361]	450
[1691]	451	ENDIF
	452
[1822]	453	CALL calc_precipitation_amount
	454
[1115]	455	END SUBROUTINE microphysics_control
	456
[1682]	457	!------------------------------------------------------------------------------!
	458	! Description:
	459	! ------------
[2155]	460	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
	461	!> evaporation of rain drops due to too small or too big weights
[1682]	462	!> of rain drops (Stevens and Seifert, 2008).
	463	!------------------------------------------------------------------------------!
[1115]	464	SUBROUTINE adjust_cloud
	465
[1361]	466	USE arrays_3d, &
[2292]	467	ONLY: qc, nc, qr, nr
[1361]	468
	469	USE cloud_parameters, &
[1849]	470	ONLY: hyrho
[1361]	471
[2292]	472	USE control_parameters, &
	473	ONLY: microphysics_morrison
	474
[1361]	475	USE cpulog, &
	476	ONLY: cpu_log, log_point_s
	477
	478	USE indices, &
[2317]	479	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	480
[1320]	481	USE kinds
[1022]	482
	483	IMPLICIT NONE
	484
[1682]	485	INTEGER(iwp) :: i !<
	486	INTEGER(iwp) :: j !<
	487	INTEGER(iwp) :: k !<
[1022]	488
[1361]	489	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
	490
[2155]	491	DO i = nxlg, nxrg
	492	DO j = nysg, nyng
[2232]	493	DO k = nzb+1, nzt
[1361]	494	IF ( qr(k,j,i) <= eps_sb ) THEN
	495	qr(k,j,i) = 0.0_wp
	496	nr(k,j,i) = 0.0_wp
	497	ELSE
	498	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
[2232]	499	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * &
[2312]	500	MERGE( 1.0_wp, 0.0_wp, &
[2232]	501	BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]	502	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
[2232]	503	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * &
[2312]	504	MERGE( 1.0_wp, 0.0_wp, &
[2232]	505	BTEST( wall_flags_0(k,j,i), 0 ) )
[1361]	506	ENDIF
	507	ENDIF
[2292]	508	IF ( microphysics_morrison ) THEN
	509	IF ( qc(k,j,i) <= eps_sb ) THEN
	510	qc(k,j,i) = 0.0_wp
	511	nc(k,j,i) = 0.0_wp
	512	ELSE
	513	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
	514	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * &
[2312]	515	MERGE( 1.0_wp, 0.0_wp, &
[2292]	516	BTEST( wall_flags_0(k,j,i), 0 ) )
	517	ENDIF
	518	ENDIF
	519	ENDIF
[1022]	520	ENDDO
	521	ENDDO
	522	ENDDO
	523
[1361]	524	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
	525
[1115]	526	END SUBROUTINE adjust_cloud
[1022]	527
[2292]	528	!------------------------------------------------------------------------------!
	529	! Description:
	530	! ------------
	531	!> Calculate number of activated condensation nucleii after simple activation
	532	!> scheme of Twomey, 1959.
	533	!------------------------------------------------------------------------------!
	534	SUBROUTINE activation
[1106]	535
[2292]	536	USE arrays_3d, &
	537	ONLY: hyp, nc, nr, pt, q, qc, qr
	538
	539	USE cloud_parameters, &
[2375]	540	ONLY: hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute, &
	541	molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt, &
	542	vanthoff
[2292]	543
	544	USE constants, &
	545	ONLY: pi
	546
	547	USE cpulog, &
	548	ONLY: cpu_log, log_point_s
	549
	550	USE indices, &
	551	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	552
	553	USE kinds
	554
	555	USE control_parameters, &
	556	ONLY: simulated_time
	557
	558	IMPLICIT NONE
	559
	560	INTEGER(iwp) :: i !<
	561	INTEGER(iwp) :: j !<
	562	INTEGER(iwp) :: k !<
	563
	564	REAL(wp) :: activ !<
	565	REAL(wp) :: afactor !<
	566	REAL(wp) :: alpha !<
	567	REAL(wp) :: beta_act !<
	568	REAL(wp) :: bfactor !<
	569	REAL(wp) :: e_s !<
	570	REAL(wp) :: k_act !<
	571	REAL(wp) :: n_act !<
	572	REAL(wp) :: n_ccn !<
	573	REAL(wp) :: q_s !<
	574	REAL(wp) :: s_0 !<
	575	REAL(wp) :: sat !<
	576	REAL(wp) :: sat_max !<
	577	REAL(wp) :: sigma !<
[2375]	578	REAL(wp) :: sigma_act !<
[2292]	579	REAL(wp) :: t_int !<
	580	REAL(wp) :: t_l !<
	581
	582	CALL cpu_log( log_point_s(65), 'activation', 'start' )
	583
	584	DO i = nxlg, nxrg
	585	DO j = nysg, nyng
	586	DO k = nzb+1, nzt
	587
	588	!
	589	!-- Actual liquid water temperature:
	590	t_l = t_d_pt(k) * pt(k,j,i)
	591
	592	!
	593	!-- Calculate actual temperature
	594	t_int = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
	595	!
	596	!-- Saturation vapor pressure at t_l:
	597	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	598	( t_l - 35.86_wp ) &
	599	)
	600	!
	601	!-- Computation of saturation humidity:
	602	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	603	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	604	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
	605	( 1.0_wp + alpha * q_s )
	606
	607	!-- Supersaturation:
	608	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
	609
	610	!
	611	!-- Prescribe parameters for activation
	612	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
[2312]	613	k_act = 0.7_wp
[2522]	614	activ = 0.0_wp
[2292]	615
[2522]	616
	617	IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
[2292]	618	!
[2312]	619	!-- Compute the number of activated Aerosols
[2292]	620	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
	621	n_act = na_init * sat**k_act
	622	!
	623	!-- Compute the number of cloud droplets
	624	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
	625	! activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
	626
	627	!
	628	!-- Compute activation rate after Khairoutdinov and Kogan
[2312]	629	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
	630	sat_max = 1.0_wp / 100.0_wp
	631	activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN &
	632	( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / &
[2292]	633	dt_micro
[2522]	634	ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN
[2292]	635	!
[2312]	636	!-- Curvature effect (afactor) with surface tension
[2292]	637	!-- parameterization by Straka (2009)
	638	sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
	639	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
	640	!
	641	!-- Solute effect (bfactor)
	642	bfactor = vanthoff * molecular_weight_of_water * &
	643	rho_s / ( molecular_weight_of_solute * rho_l )
	644
	645	!
[2312]	646	!-- Prescribe power index that describes the soluble fraction
[2375]	647	!-- of an aerosol particle (beta)
[2292]	648	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
[2375]	649	beta_act = 0.5_wp
	650	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
[2292]	651	!
[2312]	652	!-- Calculate mean geometric supersaturation (s_0) with
[2292]	653	!-- parameterization by Khvorostyanov and Curry (2006)
[2375]	654	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
[2292]	655	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
	656
	657	!
[2312]	658	!-- Calculate number of activated CCN as a function of
[2292]	659	!-- supersaturation and taking Koehler theory into account
[2312]	660	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
	661	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
[2375]	662	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
[2292]	663	activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
	664	ENDIF
	665
	666	nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
	667
	668	ENDDO
	669	ENDDO
	670	ENDDO
	671
	672	CALL cpu_log( log_point_s(65), 'activation', 'stop' )
	673
	674	END SUBROUTINE activation
	675
	676
[1682]	677	!------------------------------------------------------------------------------!
	678	! Description:
	679	! ------------
[2312]	680	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
[2292]	681	!> Kogan, 2000).
	682	!------------------------------------------------------------------------------!
	683	SUBROUTINE condensation
	684
	685	USE arrays_3d, &
	686	ONLY: hyp, nr, pt, q, qc, qr, nc
	687
	688	USE cloud_parameters, &
	689	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
	690
	691	USE constants, &
	692	ONLY: pi
	693
	694	USE cpulog, &
	695	ONLY: cpu_log, log_point_s
	696
	697	USE indices, &
	698	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	699
	700	USE kinds
	701
	702	USE control_parameters, &
	703	ONLY: simulated_time
	704
	705	IMPLICIT NONE
	706
	707	INTEGER(iwp) :: i !<
	708	INTEGER(iwp) :: j !<
	709	INTEGER(iwp) :: k !<
	710
	711	REAL(wp) :: alpha !<
	712	REAL(wp) :: cond !<
	713	REAL(wp) :: cond_max !<
	714	REAL(wp) :: dc !<
	715	REAL(wp) :: e_s !<
	716	REAL(wp) :: evap !<
	717	REAL(wp) :: evap_nc !<
	718	REAL(wp) :: g_fac !<
	719	REAL(wp) :: nc_0 !<
	720	REAL(wp) :: q_s !<
	721	REAL(wp) :: sat !<
	722	REAL(wp) :: t_l !<
	723	REAL(wp) :: temp !<
	724	REAL(wp) :: xc !<
	725
	726	CALL cpu_log( log_point_s(66), 'condensation', 'start' )
	727
	728	DO i = nxlg, nxrg
	729	DO j = nysg, nyng
	730	DO k = nzb+1, nzt
	731	!
	732	!-- Actual liquid water temperature:
	733	t_l = t_d_pt(k) * pt(k,j,i)
	734	!
	735	!-- Saturation vapor pressure at t_l:
	736	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	737	( t_l - 35.86_wp ) &
	738	)
	739	!
	740	!-- Computation of saturation humidity:
	741	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	742	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	743	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
	744	( 1.0_wp + alpha * q_s )
	745
	746	!-- Supersaturation:
	747	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
	748
	749	!
	750	!-- Actual temperature:
	751	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
	752
	753	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	754	l_v / ( thermal_conductivity_l * temp ) &
	755	+ r_v * temp / ( diff_coeff_l * e_s ) &
	756	)
	757	!
	758	!-- Mean weight of cloud drops
	759	IF ( nc(k,j,i) <= 0.0_wp) CYCLE
[2312]	760	xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
[2292]	761	!
	762	!-- Weight averaged diameter of cloud drops:
	763	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	764	!
	765	!-- Integral diameter of cloud drops
[2312]	766	nc_0 = nc(k,j,i) * dc
[2292]	767	!
	768	!-- Condensation needs only to be calculated in supersaturated regions
	769	IF ( sat > 0.0_wp ) THEN
	770	!
	771	!-- Condensation rate of cloud water content
	772	!-- after KK scheme.
	773	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
	774	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	775	cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
	776	cond = MIN( cond, cond_max / dt_micro )
[2312]	777
[2292]	778	qc(k,j,i) = qc(k,j,i) + cond * dt_micro
	779	ELSEIF ( sat < 0.0_wp ) THEN
	780	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	781	evap = MAX( evap, -qc(k,j,i) / dt_micro )
	782
	783	qc(k,j,i) = qc(k,j,i) + evap * dt_micro
	784	ENDIF
	785	IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN
	786	nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
	787	ENDIF
	788	ENDDO
	789	ENDDO
	790	ENDDO
	791
	792	CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
	793
	794	END SUBROUTINE condensation
	795
	796
	797	!------------------------------------------------------------------------------!
	798	! Description:
	799	! ------------
[1682]	800	!> Autoconversion rate (Seifert and Beheng, 2006).
	801	!------------------------------------------------------------------------------!
[1000]	802	SUBROUTINE autoconversion
	803
[1361]	804	USE arrays_3d, &
[2292]	805	ONLY: diss, dzu, nc, nr, qc, qr
[1361]	806
	807	USE cloud_parameters, &
[1849]	808	ONLY: hyrho
[1361]	809
	810	USE control_parameters, &
[2292]	811	ONLY: microphysics_morrison, rho_surface
[1361]	812
	813	USE cpulog, &
	814	ONLY: cpu_log, log_point_s
	815
	816	USE grid_variables, &
	817	ONLY: dx, dy
	818
	819	USE indices, &
[2317]	820	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	821
[1320]	822	USE kinds
[1000]	823
	824	IMPLICIT NONE
	825
[1682]	826	INTEGER(iwp) :: i !<
	827	INTEGER(iwp) :: j !<
	828	INTEGER(iwp) :: k !<
[1000]	829
[2155]	830	REAL(wp) :: alpha_cc !<
[1682]	831	REAL(wp) :: autocon !<
	832	REAL(wp) :: dissipation !<
[2232]	833	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	834	REAL(wp) :: k_au !<
	835	REAL(wp) :: l_mix !<
[2292]	836	REAL(wp) :: nc_auto !<
[1682]	837	REAL(wp) :: nu_c !<
	838	REAL(wp) :: phi_au !<
	839	REAL(wp) :: r_cc !<
	840	REAL(wp) :: rc !<
	841	REAL(wp) :: re_lambda !<
	842	REAL(wp) :: sigma_cc !<
	843	REAL(wp) :: tau_cloud !<
	844	REAL(wp) :: xc !<
[1361]	845
	846	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
	847
[2155]	848	DO i = nxlg, nxrg
	849	DO j = nysg, nyng
[2232]	850	DO k = nzb+1, nzt
	851	!
	852	!-- Predetermine flag to mask topography
	853	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]	854
[2292]	855	nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
	856
[2522]	857	IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN
[1361]	858
	859	k_au = k_cc / ( 20.0_wp * x0 )
	860	!
	861	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
	862	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
[2522]	863	tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + &
	864	qc(k,j,i) ), 0.0_wp )
[1361]	865	!
[2155]	866	!-- Universal function for autoconversion process
[1361]	867	!-- (Seifert and Beheng, 2006):
	868	phi_au = 600.0_wp * tau_cloud*0.68_wp &
	869	( 1.0_wp - tau_cloud0.68_wp )3
	870	!
	871	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
	872	!-- (Use constant nu_c = 1.0_wp instead?)
	873	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
	874	!
	875	!-- Mean weight of cloud droplets:
[2522]	876	xc = MAX( hyrho(k) * qc(k,j,i) / nc_auto, xcmin)
[1361]	877	!
[2155]	878	!-- Parameterized turbulence effects on autoconversion (Seifert,
[1361]	879	!-- Nuijens and Stevens, 2010)
[1831]	880	IF ( collision_turbulence ) THEN
[1361]	881	!
	882	!-- Weight averaged radius of cloud droplets:
	883	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	884
	885	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	886	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	887	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
	888	!
[2155]	889	!-- Mixing length (neglecting distance to ground and
[1361]	890	!-- stratification)
	891	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
	892	!
[2155]	893	!-- Limit dissipation rate according to Seifert, Nuijens and
[1361]	894	!-- Stevens (2010)
	895	dissipation = MIN( 0.06_wp, diss(k,j,i) )
	896	!
	897	!-- Compute Taylor-microscale Reynolds number:
	898	re_lambda = 6.0_wp / 11.0_wp * &
	899	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	900	SQRT( 15.0_wp / kin_vis_air ) * &
	901	dissipation**( 1.0_wp / 6.0_wp )
	902	!
[2155]	903	!-- The factor of 1.0E4 is needed to convert the dissipation
[1361]	904	!-- rate from m2 s-3 to cm2 s-3.
	905	k_au = k_au * ( 1.0_wp + &
	906	dissipation * 1.0E4_wp * &
	907	( re_lambda * 1.0E-3_wp )*0.25_wp &
[2155]	908	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
[1361]	909	r_cc ) / &
	910	sigma_cc )**2 &
	911	) + beta_cc &
	912	) &
	913	)
	914	ENDIF
	915	!
	916	!-- Autoconversion rate (Seifert and Beheng, 2006):
	917	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	918	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
	919	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
	920	rho_surface
	921	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
	922
[2232]	923	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
	924	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
	925	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro &
	926	* flag
[2292]	927	IF ( microphysics_morrison ) THEN
[2312]	928	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * &
[2292]	929	autocon / x0 * hyrho(k) * dt_micro * flag )
	930	ENDIF
[1361]	931
	932	ENDIF
	933
[1000]	934	ENDDO
	935	ENDDO
	936	ENDDO
	937
[1361]	938	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
	939
[1000]	940	END SUBROUTINE autoconversion
	941
[1106]	942
[1682]	943	!------------------------------------------------------------------------------!
	944	! Description:
	945	! ------------
[1822]	946	!> Autoconversion process (Kessler, 1969).
	947	!------------------------------------------------------------------------------!
	948	SUBROUTINE autoconversion_kessler
	949
	950	USE arrays_3d, &
[1849]	951	ONLY: dzw, pt, prr, q, qc
[1822]	952
	953	USE cloud_parameters, &
[1849]	954	ONLY: l_d_cp, pt_d_t
[1822]	955
	956	USE indices, &
[2317]	957	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1822]	958
	959	USE kinds
	960
[2317]	961	USE surface_mod, &
	962	ONLY: get_topography_top_index
[1822]	963
[2317]	964
[1822]	965	IMPLICIT NONE
	966
[2232]	967	INTEGER(iwp) :: i !<
	968	INTEGER(iwp) :: j !<
	969	INTEGER(iwp) :: k !<
	970	INTEGER(iwp) :: k_wall !< topgraphy top index
[1822]	971
	972	REAL(wp) :: dqdt_precip !<
[2232]	973	REAL(wp) :: flag !< flag to mask topography grid points
[1822]	974
[2155]	975	DO i = nxlg, nxrg
	976	DO j = nysg, nyng
[2232]	977	!
	978	!-- Determine vertical index of topography top
[2317]	979	k_wall = get_topography_top_index( j, i, 's' )
[2232]	980	DO k = nzb+1, nzt
	981	!
	982	!-- Predetermine flag to mask topography
	983	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]	984
	985	IF ( qc(k,j,i) > ql_crit ) THEN
	986	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
	987	ELSE
	988	dqdt_precip = 0.0_wp
	989	ENDIF
	990
[2232]	991	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
	992	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag
[1845]	993	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * &
[2232]	994	pt_d_t(k) * flag
[1822]	995
	996	!
[1845]	997	!-- Compute the rain rate (stored on surface grid point)
[2232]	998	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]	999
	1000	ENDDO
	1001	ENDDO
	1002	ENDDO
	1003
	1004	END SUBROUTINE autoconversion_kessler
	1005
	1006
	1007	!------------------------------------------------------------------------------!
	1008	! Description:
	1009	! ------------
[1682]	1010	!> Accretion rate (Seifert and Beheng, 2006).
	1011	!------------------------------------------------------------------------------!
[1005]	1012	SUBROUTINE accretion
[1000]	1013
[1361]	1014	USE arrays_3d, &
[2292]	1015	ONLY: diss, qc, qr, nc
[1361]	1016
	1017	USE cloud_parameters, &
[1849]	1018	ONLY: hyrho
[1361]	1019
	1020	USE control_parameters, &
[2292]	1021	ONLY: microphysics_morrison, rho_surface
[1361]	1022
	1023	USE cpulog, &
	1024	ONLY: cpu_log, log_point_s
	1025
	1026	USE indices, &
[2317]	1027	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1028
[1320]	1029	USE kinds
[1005]	1030
[1000]	1031	IMPLICIT NONE
	1032
[1682]	1033	INTEGER(iwp) :: i !<
	1034	INTEGER(iwp) :: j !<
	1035	INTEGER(iwp) :: k !<
[1000]	1036
[1682]	1037	REAL(wp) :: accr !<
[2232]	1038	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1039	REAL(wp) :: k_cr !<
[2292]	1040	REAL(wp) :: nc_accr !<
[1682]	1041	REAL(wp) :: phi_ac !<
	1042	REAL(wp) :: tau_cloud !<
[2292]	1043	REAL(wp) :: xc !<
[1361]	1044
[2292]	1045
[1361]	1046	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
	1047
[2155]	1048	DO i = nxlg, nxrg
	1049	DO j = nysg, nyng
[2232]	1050	DO k = nzb+1, nzt
	1051	!
	1052	!-- Predetermine flag to mask topography
	1053	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1000]	1054
[2292]	1055	nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
	1056
	1057	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) &
	1058	.AND. ( nc_accr > eps_mr ) ) THEN
[1361]	1059	!
	1060	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]	1061	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
[1361]	1062	!
[2155]	1063	!-- Universal function for accretion process (Seifert and
[1361]	1064	!-- Beheng, 2001):
	1065	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[2292]	1066
[1361]	1067	!
[2292]	1068	!-- Mean weight of cloud drops
[2312]	1069	xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
[2292]	1070	!
[2155]	1071	!-- Parameterized turbulence effects on autoconversion (Seifert,
	1072	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	1073	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]	1074	IF ( collision_turbulence ) THEN
[1361]	1075	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	1076	MIN( 600.0_wp, &
	1077	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	1078	)
	1079	ELSE
[2155]	1080	k_cr = k_cr0
[1361]	1081	ENDIF
	1082	!
	1083	!-- Accretion rate (Seifert and Beheng, 2006):
	1084	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
	1085	SQRT( rho_surface * hyrho(k) )
	1086	accr = MIN( accr, qc(k,j,i) / dt_micro )
	1087
[2232]	1088	qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
	1089	qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
[2292]	1090	IF ( microphysics_morrison ) THEN
	1091	nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), &
	1092	accr / xc * hyrho(k) * dt_micro * flag)
	1093	ENDIF
[1361]	1094
	1095	ENDIF
	1096
[1005]	1097	ENDDO
	1098	ENDDO
[1000]	1099	ENDDO
	1100
[1361]	1101	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
	1102
[1005]	1103	END SUBROUTINE accretion
[1000]	1104
[1106]	1105
[1682]	1106	!------------------------------------------------------------------------------!
	1107	! Description:
	1108	! ------------
	1109	!> Collisional breakup rate (Seifert, 2008).
	1110	!------------------------------------------------------------------------------!
[1005]	1111	SUBROUTINE selfcollection_breakup
[1000]	1112
[1361]	1113	USE arrays_3d, &
	1114	ONLY: nr, qr
	1115
	1116	USE cloud_parameters, &
[1849]	1117	ONLY: hyrho
[1361]	1118
	1119	USE control_parameters, &
[1849]	1120	ONLY: rho_surface
[1361]	1121
	1122	USE cpulog, &
	1123	ONLY: cpu_log, log_point_s
	1124
	1125	USE indices, &
[2317]	1126	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1127
[1320]	1128	USE kinds
[2155]	1129
[1000]	1130	IMPLICIT NONE
	1131
[1682]	1132	INTEGER(iwp) :: i !<
	1133	INTEGER(iwp) :: j !<
	1134	INTEGER(iwp) :: k !<
[1000]	1135
[1682]	1136	REAL(wp) :: breakup !<
	1137	REAL(wp) :: dr !<
[2232]	1138	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1139	REAL(wp) :: phi_br !<
	1140	REAL(wp) :: selfcoll !<
[1361]	1141
	1142	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
	1143
[2155]	1144	DO i = nxlg, nxrg
	1145	DO j = nysg, nyng
[2232]	1146	DO k = nzb+1, nzt
	1147	!
	1148	!-- Predetermine flag to mask topography
	1149	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1150
[1361]	1151	IF ( qr(k,j,i) > eps_sb ) THEN
	1152	!
	1153	!-- Selfcollection rate (Seifert and Beheng, 2001):
	1154	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
	1155	SQRT( hyrho(k) * rho_surface )
	1156	!
	1157	!-- Weight averaged diameter of rain drops:
	1158	dr = ( hyrho(k) * qr(k,j,i) / &
	1159	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	1160	!
	1161	!-- Collisional breakup rate (Seifert, 2008):
	1162	IF ( dr >= 0.3E-3_wp ) THEN
	1163	phi_br = k_br * ( dr - 1.1E-3_wp )
	1164	breakup = selfcoll * ( phi_br + 1.0_wp )
	1165	ELSE
	1166	breakup = 0.0_wp
	1167	ENDIF
[1000]	1168
[1361]	1169	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
[2232]	1170	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
[1361]	1171
[2155]	1172	ENDIF
[1000]	1173	ENDDO
	1174	ENDDO
	1175	ENDDO
	1176
[1361]	1177	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
	1178
[1005]	1179	END SUBROUTINE selfcollection_breakup
[1000]	1180
[1106]	1181
[1682]	1182	!------------------------------------------------------------------------------!
	1183	! Description:
	1184	! ------------
[2155]	1185	!> Evaporation of precipitable water. Condensation is neglected for
[1682]	1186	!> precipitable water.
	1187	!------------------------------------------------------------------------------!
[1012]	1188	SUBROUTINE evaporation_rain
[1000]	1189
[1361]	1190	USE arrays_3d, &
	1191	ONLY: hyp, nr, pt, q, qc, qr
	1192
	1193	USE cloud_parameters, &
[1849]	1194	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1361]	1195
	1196	USE constants, &
	1197	ONLY: pi
	1198
	1199	USE cpulog, &
	1200	ONLY: cpu_log, log_point_s
	1201
	1202	USE indices, &
[2317]	1203	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1204
[1320]	1205	USE kinds
[1012]	1206
	1207	IMPLICIT NONE
	1208
[1682]	1209	INTEGER(iwp) :: i !<
	1210	INTEGER(iwp) :: j !<
	1211	INTEGER(iwp) :: k !<
[1361]	1212
[1682]	1213	REAL(wp) :: alpha !<
	1214	REAL(wp) :: dr !<
	1215	REAL(wp) :: e_s !<
	1216	REAL(wp) :: evap !<
	1217	REAL(wp) :: evap_nr !<
	1218	REAL(wp) :: f_vent !<
[2232]	1219	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1220	REAL(wp) :: g_evap !<
	1221	REAL(wp) :: lambda_r !<
	1222	REAL(wp) :: mu_r !<
	1223	REAL(wp) :: mu_r_2 !<
	1224	REAL(wp) :: mu_r_5d2 !<
	1225	REAL(wp) :: nr_0 !<
	1226	REAL(wp) :: q_s !<
	1227	REAL(wp) :: sat !<
	1228	REAL(wp) :: t_l !<
	1229	REAL(wp) :: temp !<
	1230	REAL(wp) :: xr !<
[1361]	1231
	1232	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
	1233
[2155]	1234	DO i = nxlg, nxrg
	1235	DO j = nysg, nyng
[2232]	1236	DO k = nzb+1, nzt
	1237	!
	1238	!-- Predetermine flag to mask topography
	1239	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1240
[1361]	1241	IF ( qr(k,j,i) > eps_sb ) THEN
	1242	!
	1243	!-- Actual liquid water temperature:
	1244	t_l = t_d_pt(k) * pt(k,j,i)
	1245	!
	1246	!-- Saturation vapor pressure at t_l:
	1247	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	1248	( t_l - 35.86_wp ) &
	1249	)
	1250	!
	1251	!-- Computation of saturation humidity:
	1252	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	1253	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	1254	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
	1255	( 1.0_wp + alpha * q_s )
	1256	!
	1257	!-- Supersaturation:
	1258	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
	1259	!
	1260	!-- Evaporation needs only to be calculated in subsaturated regions
	1261	IF ( sat < 0.0_wp ) THEN
	1262	!
	1263	!-- Actual temperature:
	1264	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
[2155]	1265
[1361]	1266	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	1267	l_v / ( thermal_conductivity_l * temp ) &
	1268	+ r_v * temp / ( diff_coeff_l * e_s ) &
	1269	)
	1270	!
	1271	!-- Mean weight of rain drops
	1272	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
	1273	!
	1274	!-- Weight averaged diameter of rain drops:
	1275	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
	1276	!
[2155]	1277	!-- Compute ventilation factor and intercept parameter
[1361]	1278	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	1279	IF ( ventilation_effect ) THEN
	1280	!
[2155]	1281	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
[1361]	1282	!-- 2005; Stevens and Seifert, 2008):
	1283	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	1284	( dr - 1.4E-3_wp ) ) )
	1285	!
	1286	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1287	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1288	( mu_r + 1.0_wp ) &
	1289	)**( 1.0_wp / 3.0_wp ) / dr
[1012]	1290
[1361]	1291	mu_r_2 = mu_r + 2.0_wp
[2155]	1292	mu_r_5d2 = mu_r + 2.5_wp
[1361]	1293
	1294	f_vent = a_vent * gamm( mu_r_2 ) * &
	1295	lambda_r*( -mu_r_2 ) + b_vent &
	1296	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
	1297	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
	1298	( 1.0_wp - &
	1299	0.5_wp * ( b_term / a_term ) * &
	1300	( lambda_r / ( c_term + lambda_r ) &
	1301	)**mu_r_5d2 - &
	1302	0.125_wp * ( b_term / a_term )*2 &
	1303	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	1304	)**mu_r_5d2 - &
	1305	0.0625_wp * ( b_term / a_term )*3 &
	1306	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	1307	)**mu_r_5d2 - &
[2155]	1308	0.0390625_wp * ( b_term / a_term )*4 &
[1361]	1309	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	1310	)**mu_r_5d2 &
	1311	)
	1312
	1313	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
[2155]	1314	gamm( mu_r + 1.0_wp )
[1361]	1315	ELSE
	1316	f_vent = 1.0_wp
	1317	nr_0 = nr(k,j,i) * dr
	1318	ENDIF
	1319	!
[2155]	1320	!-- Evaporation rate of rain water content (Seifert and
[1361]	1321	!-- Beheng, 2006):
	1322	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
	1323	hyrho(k)
	1324	evap = MAX( evap, -qr(k,j,i) / dt_micro )
	1325	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	1326	-nr(k,j,i) / dt_micro )
	1327
[2232]	1328	qr(k,j,i) = qr(k,j,i) + evap * dt_micro * flag
	1329	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
[1361]	1330
	1331	ENDIF
[2155]	1332	ENDIF
[1361]	1333
[1012]	1334	ENDDO
	1335	ENDDO
	1336	ENDDO
	1337
[1361]	1338	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
	1339
[1012]	1340	END SUBROUTINE evaporation_rain
	1341
[1106]	1342
[1682]	1343	!------------------------------------------------------------------------------!
	1344	! Description:
	1345	! ------------
	1346	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
	1347	!------------------------------------------------------------------------------!
[1012]	1348	SUBROUTINE sedimentation_cloud
	1349
[1361]	1350	USE arrays_3d, &
[2292]	1351	ONLY: ddzu, dzu, nc, pt, prr, q, qc
[1361]	1352
	1353	USE cloud_parameters, &
[1849]	1354	ONLY: hyrho, l_d_cp, pt_d_t
[1361]	1355
	1356	USE control_parameters, &
[2312]	1357	ONLY: call_microphysics_at_all_substeps, &
[2292]	1358	intermediate_timestep_count, microphysics_morrison
[1361]	1359
	1360	USE cpulog, &
	1361	ONLY: cpu_log, log_point_s
	1362
	1363	USE indices, &
[2317]	1364	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1365
[1320]	1366	USE kinds
[1691]	1367
	1368	USE statistics, &
	1369	ONLY: weight_substep
	1370
[2155]	1371
[1012]	1372	IMPLICIT NONE
	1373
[1849]	1374	INTEGER(iwp) :: i !<
	1375	INTEGER(iwp) :: j !<
	1376	INTEGER(iwp) :: k !<
[1361]	1377
[2292]	1378	REAL(wp) :: flag !< flag to mask topography grid points
	1379	REAL(wp) :: nc_sedi !<
	1380
[2232]	1381	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[2292]	1382	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
[1361]	1383
[2292]	1384
[1361]	1385	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
	1386
	1387	sed_qc(nzt+1) = 0.0_wp
[2375]	1388	sed_nc(nzt+1) = 0.0_wp
[1361]	1389
[2155]	1390	DO i = nxlg, nxrg
	1391	DO j = nysg, nyng
[2232]	1392	DO k = nzt, nzb+1, -1
	1393	!
	1394	!-- Predetermine flag to mask topography
	1395	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	1396	nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison )
[1012]	1397
[2292]	1398	!
	1399	!-- Sedimentation fluxes for number concentration are only calculated
	1400	!-- for cloud_scheme = 'morrison'
	1401	IF ( microphysics_morrison ) THEN
	1402	IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN
[2312]	1403	sed_nc(k) = sed_qc_const * &
[2292]	1404	( qc(k,j,i) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
	1405	( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
	1406	ELSE
	1407	sed_nc(k) = 0.0_wp
	1408	ENDIF
	1409
	1410	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
	1411	nc(k,j,i) / dt_micro + sed_nc(k+1) &
	1412	) * flag
	1413
	1414	nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * &
	1415	ddzu(k+1) / hyrho(k) * dt_micro * flag
	1416	ENDIF
	1417
[2375]	1418	IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN
[2292]	1419	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
[2232]	1420	( qc(k,j,i) * hyrho(k) )*( 5.0_wp / 3.0_wp ) &
	1421	flag
[1361]	1422	ELSE
	1423	sed_qc(k) = 0.0_wp
	1424	ENDIF
	1425
	1426	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
	1427	dt_micro + sed_qc(k+1) &
[2232]	1428	) * flag
[1361]	1429
	1430	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
[2232]	1431	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1432	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
[2232]	1433	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1434	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
	1435	ddzu(k+1) / hyrho(k) * l_d_cp * &
[2232]	1436	pt_d_t(k) * dt_micro * flag
[1361]	1437
[1691]	1438	!
	1439	!-- Compute the precipitation rate due to cloud (fog) droplets
[1822]	1440	IF ( call_microphysics_at_all_substeps ) THEN
	1441	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &
[2232]	1442	* weight_substep(intermediate_timestep_count) &
	1443	* flag
[1822]	1444	ELSE
[2232]	1445	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]	1446	ENDIF
	1447
[1012]	1448	ENDDO
	1449	ENDDO
	1450	ENDDO
	1451
[1361]	1452	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
	1453
[1012]	1454	END SUBROUTINE sedimentation_cloud
	1455
[1106]	1456
[1682]	1457	!------------------------------------------------------------------------------!
	1458	! Description:
	1459	! ------------
	1460	!> Computation of sedimentation flux. Implementation according to Stevens
	1461	!> and Seifert (2008). Code is based on UCLA-LES.
	1462	!------------------------------------------------------------------------------!
[1012]	1463	SUBROUTINE sedimentation_rain
	1464
[1361]	1465	USE arrays_3d, &
[1849]	1466	ONLY: ddzu, dzu, nr, pt, prr, q, qr
[1361]	1467
	1468	USE cloud_parameters, &
[1849]	1469	ONLY: hyrho, l_d_cp, pt_d_t
[1361]	1470
	1471	USE control_parameters, &
[1849]	1472	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1361]	1473	USE cpulog, &
	1474	ONLY: cpu_log, log_point_s
	1475
	1476	USE indices, &
[2317]	1477	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1361]	1478
[1320]	1479	USE kinds
[1012]	1480
[1361]	1481	USE statistics, &
	1482	ONLY: weight_substep
[2155]	1483
[2232]	1484	USE surface_mod, &
	1485	ONLY : bc_h
[2312]	1486
[1012]	1487	IMPLICIT NONE
	1488
[2232]	1489	INTEGER(iwp) :: i !< running index x direction
	1490	INTEGER(iwp) :: j !< running index y direction
	1491	INTEGER(iwp) :: k !< running index z direction
[2312]	1492	INTEGER(iwp) :: k_run !<
[2232]	1493	INTEGER(iwp) :: l !< running index of surface type
	1494	INTEGER(iwp) :: m !< running index surface elements
	1495	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	1496	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1361]	1497
[1682]	1498	REAL(wp) :: c_run !<
	1499	REAL(wp) :: d_max !<
	1500	REAL(wp) :: d_mean !<
	1501	REAL(wp) :: d_min !<
	1502	REAL(wp) :: dr !<
	1503	REAL(wp) :: flux !<
[2232]	1504	REAL(wp) :: flag !< flag to mask topography grid points
[1682]	1505	REAL(wp) :: lambda_r !<
	1506	REAL(wp) :: mu_r !<
	1507	REAL(wp) :: z_run !<
[1361]	1508
[1682]	1509	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	1510	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	1511	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	1512	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	1513	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	1514	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	1515	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	1516	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1361]	1517
	1518	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
[1682]	1519
[1361]	1520	!
[2155]	1521	!-- Compute velocities
	1522	DO i = nxlg, nxrg
	1523	DO j = nysg, nyng
[2232]	1524	DO k = nzb+1, nzt
	1525	!
	1526	!-- Predetermine flag to mask topography
	1527	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1528
[1361]	1529	IF ( qr(k,j,i) > eps_sb ) THEN
	1530	!
	1531	!-- Weight averaged diameter of rain drops:
	1532	dr = ( hyrho(k) * qr(k,j,i) / &
	1533	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
	1534	!
	1535	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	1536	!-- Stevens and Seifert, 2008):
	1537	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
	1538	( dr - 1.4E-3_wp ) ) )
	1539	!
	1540	!-- Slope parameter of gamma distribution (Seifert, 2008):
	1541	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	1542	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1012]	1543
[1361]	1544	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1545	a_term - b_term * ( 1.0_wp + &
	1546	c_term / &
	1547	lambda_r )*( -1.0_wp &
	1548	( mu_r + 1.0_wp ) ) &
	1549	) &
[2232]	1550	) * flag
[1361]	1551
	1552	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	1553	a_term - b_term * ( 1.0_wp + &
	1554	c_term / &
	1555	lambda_r )*( -1.0_wp &
	1556	( mu_r + 4.0_wp ) ) &
	1557	) &
[2232]	1558	) * flag
[1361]	1559	ELSE
	1560	w_nr(k) = 0.0_wp
	1561	w_qr(k) = 0.0_wp
	1562	ENDIF
[1012]	1563	ENDDO
[1361]	1564	!
[2312]	1565	!-- Adjust boundary values using surface data type.
[2232]	1566	!-- Upward-facing
[2312]	1567	surf_s = bc_h(0)%start_index(j,i)
[2232]	1568	surf_e = bc_h(0)%end_index(j,i)
	1569	DO m = surf_s, surf_e
	1570	k = bc_h(0)%k(m)
	1571	w_nr(k-1) = w_nr(k)
	1572	w_qr(k-1) = w_qr(k)
	1573	ENDDO
	1574	!
	1575	!-- Downward-facing
[2312]	1576	surf_s = bc_h(1)%start_index(j,i)
[2232]	1577	surf_e = bc_h(1)%end_index(j,i)
	1578	DO m = surf_s, surf_e
	1579	k = bc_h(1)%k(m)
	1580	w_nr(k+1) = w_nr(k)
	1581	w_qr(k+1) = w_qr(k)
	1582	ENDDO
	1583	!
	1584	!-- Model top boundary value
[1361]	1585	w_nr(nzt+1) = 0.0_wp
	1586	w_qr(nzt+1) = 0.0_wp
	1587	!
	1588	!-- Compute Courant number
[2232]	1589	DO k = nzb+1, nzt
	1590	!
	1591	!-- Predetermine flag to mask topography
	1592	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1593
[1361]	1594	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
	1595	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[2232]	1596	dt_micro * ddzu(k) * flag
[1361]	1597	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
	1598	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[2232]	1599	dt_micro * ddzu(k) * flag
[2155]	1600	ENDDO
[1361]	1601	!
	1602	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	1603	IF ( limiter_sedimentation ) THEN
	1604
[2232]	1605	DO k = nzb+1, nzt
	1606	!
	1607	!-- Predetermine flag to mask topography
	1608	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1609
[1646]	1610	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
[1361]	1611	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
	1612	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
	1613
	1614	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1615	2.0_wp * d_max, &
[2232]	1616	ABS( d_mean ) ) &
	1617	* flag
[1361]	1618
[1646]	1619	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
[1361]	1620	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
	1621	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
	1622
	1623	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	1624	2.0_wp * d_max, &
	1625	ABS( d_mean ) )
	1626	ENDDO
	1627
	1628	ELSE
	1629
	1630	nr_slope = 0.0_wp
	1631	qr_slope = 0.0_wp
	1632
	1633	ENDIF
	1634
	1635	sed_nr(nzt+1) = 0.0_wp
	1636	sed_qr(nzt+1) = 0.0_wp
	1637	!
	1638	!-- Compute sedimentation flux
[2232]	1639	DO k = nzt, nzb+1, -1
[1361]	1640	!
[2232]	1641	!-- Predetermine flag to mask topography
	1642	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1643	!
[2312]	1644	!-- Sum up all rain drop number densities which contribute to the flux
[1361]	1645	!-- through k-1/2
	1646	flux = 0.0_wp
	1647	z_run = 0.0_wp ! height above z(k)
	1648	k_run = k
	1649	c_run = MIN( 1.0_wp, c_nr(k) )
	1650	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1651	flux = flux + hyrho(k_run) * &
	1652	( nr(k_run,j,i) + nr_slope(k_run) * &
[2232]	1653	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) &
	1654	* flag
	1655	z_run = z_run + dzu(k_run) * flag
	1656	k_run = k_run + 1 * flag
	1657	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) &
	1658	* flag
[1361]	1659	ENDDO
	1660	!
[2155]	1661	!-- It is not allowed to sediment more rain drop number density than
[1361]	1662	!-- available
	1663	flux = MIN( flux, &
	1664	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
	1665	dt_micro &
	1666	)
	1667
[2232]	1668	sed_nr(k) = flux / dt_micro * flag
[1361]	1669	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
[2232]	1670	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1671	!
[2155]	1672	!-- Sum up all rain water content which contributes to the flux
[1361]	1673	!-- through k-1/2
	1674	flux = 0.0_wp
	1675	z_run = 0.0_wp ! height above z(k)
	1676	k_run = k
	1677	c_run = MIN( 1.0_wp, c_qr(k) )
	1678
	1679	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
	1680
	1681	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
	1682	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	1683	0.5_wp ) * c_run * dzu(k_run) * flag
	1684	z_run = z_run + dzu(k_run) * flag
	1685	k_run = k_run + 1 * flag
	1686	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) &
	1687	* flag
[1361]	1688
	1689	ENDDO
	1690	!
[2155]	1691	!-- It is not allowed to sediment more rain water content than
[1361]	1692	!-- available
	1693	flux = MIN( flux, &
	1694	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
	1695	dt_micro &
	1696	)
	1697
[2232]	1698	sed_qr(k) = flux / dt_micro * flag
[1361]	1699
	1700	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
[2232]	1701	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1702	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
[2232]	1703	ddzu(k+1) / hyrho(k) * dt_micro * flag
[1361]	1704	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
	1705	ddzu(k+1) / hyrho(k) * l_d_cp * &
[2232]	1706	pt_d_t(k) * dt_micro * flag
[1361]	1707	!
	1708	!-- Compute the rain rate
	1709	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	1710	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
[2232]	1711	* weight_substep(intermediate_timestep_count) &
	1712	* flag
[1361]	1713	ELSE
[2232]	1714	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]	1715	ENDIF
	1716
	1717	ENDDO
[1012]	1718	ENDDO
	1719	ENDDO
	1720
[1691]	1721	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
	1722
	1723	END SUBROUTINE sedimentation_rain
	1724
	1725
	1726	!------------------------------------------------------------------------------!
	1727	! Description:
	1728	! ------------
	1729	!> Computation of the precipitation amount due to gravitational settling of
	1730	!> rain and cloud (fog) droplets
	1731	!------------------------------------------------------------------------------!
	1732	SUBROUTINE calc_precipitation_amount
	1733
[1849]	1734	USE arrays_3d, &
	1735	ONLY: precipitation_amount, prr
	1736
[1691]	1737	USE cloud_parameters, &
[1849]	1738	ONLY: hyrho
[1691]	1739
	1740	USE control_parameters, &
	1741	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	1742	intermediate_timestep_count, intermediate_timestep_count_max,&
	1743	precipitation_amount_interval, time_do2d_xy
	1744
	1745	USE indices, &
[2232]	1746	ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
[1691]	1747
	1748	USE kinds
	1749
[2232]	1750	USE surface_mod, &
	1751	ONLY : bc_h
	1752
[1691]	1753	IMPLICIT NONE
	1754
[2232]	1755	INTEGER(iwp) :: i !< running index x direction
	1756	INTEGER(iwp) :: j !< running index y direction
	1757	INTEGER(iwp) :: k !< running index y direction
	1758	INTEGER(iwp) :: m !< running index surface elements
	1759	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	1760	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1691]	1761
	1762	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	1763	( .NOT. call_microphysics_at_all_substeps .OR. &
	1764	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	1765	THEN
[2232]	1766	!
[2312]	1767	!-- Run over all upward-facing surface elements, i.e. non-natural,
[2232]	1768	!-- natural and urban
	1769	DO m = 1, bc_h(0)%ns
[2312]	1770	i = bc_h(0)%i(m)
[2232]	1771	j = bc_h(0)%j(m)
	1772	k = bc_h(0)%k(m)
	1773	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	1774	prr(k,j,i) * hyrho(k) * dt_3d
	1775	ENDDO
[1691]	1776
[1361]	1777	ENDIF
	1778
[1691]	1779	END SUBROUTINE calc_precipitation_amount
[1361]	1780
[1012]	1781
[1000]	1782	!------------------------------------------------------------------------------!
[1682]	1783	! Description:
	1784	! ------------
[1849]	1785	!> Control of microphysics for grid points i,j
[1000]	1786	!------------------------------------------------------------------------------!
[1022]	1787
[1115]	1788	SUBROUTINE microphysics_control_ij( i, j )
	1789
[1320]	1790	USE arrays_3d, &
[2292]	1791	ONLY: hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu
[1115]	1792
[1320]	1793	USE cloud_parameters, &
[1849]	1794	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
[1320]	1795
	1796	USE control_parameters, &
[1849]	1797	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
	1798	intermediate_timestep_count, large_scale_forcing, &
[2292]	1799	lsf_surf, microphysics_morrison, microphysics_seifert, &
	1800	microphysics_kessler, pt_surface, rho_surface, &
	1801	surface_pressure
[1320]	1802
	1803	USE indices, &
	1804	ONLY: nzb, nzt
	1805
	1806	USE kinds
	1807
	1808	USE statistics, &
	1809	ONLY: weight_pres
	1810
[1022]	1811	IMPLICIT NONE
	1812
[1682]	1813	INTEGER(iwp) :: i !<
	1814	INTEGER(iwp) :: j !<
	1815	INTEGER(iwp) :: k !<
[1115]	1816
[1682]	1817	REAL(wp) :: t_surface !<
[1320]	1818
[1361]	1819	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[1241]	1820	!
	1821	!-- Calculate:
	1822	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
	1823	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
	1824	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
[1353]	1825	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1241]	1826	DO k = nzb, nzt+1
[1353]	1827	hyp(k) = surface_pressure * 100.0_wp * &
[1361]	1828	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
[1353]	1829	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
	1830	t_d_pt(k) = 1.0_wp / pt_d_t(k)
[2155]	1831	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
[1241]	1832	ENDDO
	1833	!
	1834	!-- Compute reference density
[1353]	1835	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
[1241]	1836	ENDIF
	1837
[1361]	1838	!
[2155]	1839	!-- Compute length of time step
[1361]	1840	IF ( call_microphysics_at_all_substeps ) THEN
	1841	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
	1842	ELSE
	1843	dt_micro = dt_3d
	1844	ENDIF
[1241]	1845
[1115]	1846	!
[1361]	1847	!-- Use 1d arrays
[1115]	1848	q_1d(:) = q(:,j,i)
	1849	pt_1d(:) = pt(:,j,i)
	1850	qc_1d(:) = qc(:,j,i)
[1822]	1851	IF ( microphysics_seifert ) THEN
[1115]	1852	qr_1d(:) = qr(:,j,i)
	1853	nr_1d(:) = nr(:,j,i)
	1854	ENDIF
[2292]	1855	IF ( microphysics_morrison ) THEN
	1856	nc_1d(:) = nc(:,j,i)
	1857	ENDIF
[1361]	1858
[2292]	1859
[1115]	1860	!
[1822]	1861	!-- Reset precipitation rate
	1862	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
	1863
	1864	!
[1115]	1865	!-- Compute cloud physics
[1822]	1866	IF( microphysics_kessler ) THEN
	1867
	1868	CALL autoconversion_kessler( i,j )
[1831]	1869	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
[1822]	1870
	1871	ELSEIF ( microphysics_seifert ) THEN
	1872
	1873	CALL adjust_cloud( i,j )
[2292]	1874	IF ( microphysics_morrison ) CALL activation( i,j )
	1875	IF ( microphysics_morrison ) CALL condensation( i,j )
[1115]	1876	CALL autoconversion( i,j )
	1877	CALL accretion( i,j )
	1878	CALL selfcollection_breakup( i,j )
	1879	CALL evaporation_rain( i,j )
	1880	CALL sedimentation_rain( i,j )
[1831]	1881	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
[1115]	1882
[1691]	1883	ENDIF
	1884
[1822]	1885	CALL calc_precipitation_amount( i,j )
	1886
[1115]	1887	!
[1361]	1888	!-- Store results on the 3d arrays
	1889	q(:,j,i) = q_1d(:)
	1890	pt(:,j,i) = pt_1d(:)
[2292]	1891	IF ( microphysics_morrison ) THEN
	1892	qc(:,j,i) = qc_1d(:)
	1893	nc(:,j,i) = nc_1d(:)
	1894	ENDIF
[1822]	1895	IF ( microphysics_seifert ) THEN
[1361]	1896	qr(:,j,i) = qr_1d(:)
	1897	nr(:,j,i) = nr_1d(:)
[1115]	1898	ENDIF
	1899
	1900	END SUBROUTINE microphysics_control_ij
	1901
[1682]	1902	!------------------------------------------------------------------------------!
	1903	! Description:
	1904	! ------------
[2155]	1905	!> Adjust number of raindrops to avoid nonlinear effects in
[1682]	1906	!> sedimentation and evaporation of rain drops due to too small or
	1907	!> too big weights of rain drops (Stevens and Seifert, 2008).
	1908	!> The same procedure is applied to cloud droplets if they are determined
	1909	!> prognostically. Call for grid point i,j
	1910	!------------------------------------------------------------------------------!
[1115]	1911	SUBROUTINE adjust_cloud_ij( i, j )
	1912
[1320]	1913	USE cloud_parameters, &
[1849]	1914	ONLY: hyrho
[1320]	1915
[2292]	1916	USE control_parameters, &
	1917	ONLY: microphysics_morrison
	1918
[1320]	1919	USE indices, &
[2232]	1920	ONLY: nzb, nzt, wall_flags_0
[1320]	1921
	1922	USE kinds
	1923
[1115]	1924	IMPLICIT NONE
	1925
[1682]	1926	INTEGER(iwp) :: i !<
	1927	INTEGER(iwp) :: j !<
	1928	INTEGER(iwp) :: k !<
	1929
[2232]	1930	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1022]	1931
[2232]	1932	DO k = nzb+1, nzt
	1933	!
	1934	!-- Predetermine flag to mask topography
	1935	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	1936
[1361]	1937	IF ( qr_1d(k) <= eps_sb ) THEN
	1938	qr_1d(k) = 0.0_wp
	1939	nr_1d(k) = 0.0_wp
[1065]	1940	ELSE
[1022]	1941	!
[2155]	1942	!-- Adjust number of raindrops to avoid nonlinear effects in
[1048]	1943	!-- sedimentation and evaporation of rain drops due to too small or
[1065]	1944	!-- too big weights of rain drops (Stevens and Seifert, 2008).
[1361]	1945	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
[2232]	1946	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
[1361]	1947	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
[2232]	1948	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
[1048]	1949	ENDIF
[1115]	1950
[1022]	1951	ENDIF
[1115]	1952
[2292]	1953	IF ( microphysics_morrison ) THEN
	1954	IF ( qc_1d(k) <= eps_sb ) THEN
	1955	qc_1d(k) = 0.0_wp
	1956	nc_1d(k) = 0.0_wp
	1957	ELSE
	1958	IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) ) THEN
	1959	nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag
	1960	ENDIF
	1961	ENDIF
	1962	ENDIF
	1963
[1022]	1964	ENDDO
	1965
[1115]	1966	END SUBROUTINE adjust_cloud_ij
[1022]	1967
[2292]	1968	!------------------------------------------------------------------------------!
	1969	! Description:
	1970	! ------------
	1971	!> Calculate number of activated condensation nucleii after simple activation
	1972	!> scheme of Twomey, 1959.
	1973	!------------------------------------------------------------------------------!
	1974	SUBROUTINE activation_ij( i, j )
[1106]	1975
[2292]	1976	USE arrays_3d, &
	1977	ONLY: hyp, nr, pt, q, qc, qr, nc
	1978
	1979	USE cloud_parameters, &
[2375]	1980	ONLY: hyrho, l_d_cp, l_d_r, l_v, molecular_weight_of_solute, &
	1981	molecular_weight_of_water, rho_l, rho_s, r_v, t_d_pt, &
	1982	vanthoff
[2292]	1983
	1984	USE constants, &
	1985	ONLY: pi
	1986
	1987	USE cpulog, &
	1988	ONLY: cpu_log, log_point_s
	1989
	1990	USE indices, &
	1991	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	1992
	1993	USE kinds
	1994
	1995	USE control_parameters, &
	1996	ONLY: simulated_time
	1997
	1998	IMPLICIT NONE
	1999
	2000	INTEGER(iwp) :: i !<
	2001	INTEGER(iwp) :: j !<
	2002	INTEGER(iwp) :: k !<
	2003
	2004	REAL(wp) :: activ !<
	2005	REAL(wp) :: afactor !<
	2006	REAL(wp) :: alpha !<
	2007	REAL(wp) :: beta_act !<
	2008	REAL(wp) :: bfactor !<
	2009	REAL(wp) :: e_s !<
	2010	REAL(wp) :: k_act !<
	2011	REAL(wp) :: n_act !<
	2012	REAL(wp) :: n_ccn !<
	2013	REAL(wp) :: q_s !<
	2014	REAL(wp) :: s_0 !<
	2015	REAL(wp) :: sat !<
	2016	REAL(wp) :: sat_max !<
	2017	REAL(wp) :: sigma !<
[2375]	2018	REAL(wp) :: sigma_act !<
[2292]	2019	REAL(wp) :: t_int !<
	2020	REAL(wp) :: t_l !<
	2021
	2022	DO k = nzb+1, nzt
	2023	!
	2024	!-- Actual liquid water temperature:
	2025	t_l = t_d_pt(k) * pt_1d(k)
	2026
	2027	!
	2028	!-- Calculate actual temperature
	2029	t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp
	2030	!
	2031	!-- Saturation vapor pressure at t_l:
	2032	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	2033	( t_l - 35.86_wp ) &
	2034	)
	2035	!
	2036	!-- Computation of saturation humidity:
	2037	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	2038	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	2039	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / &
	2040	( 1.0_wp + alpha * q_s )
	2041
	2042	!-- Supersaturation:
	2043	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
	2044
	2045	!
	2046	!-- Prescribe parameters for activation
	2047	!-- (see: Bott + Trautmann, 2002, Atm. Res., 64)
[2312]	2048	k_act = 0.7_wp
[2522]	2049	activ = 0.0_wp
[2292]	2050
	2051	IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
	2052	!
[2312]	2053	!-- Compute the number of activated Aerosols
[2292]	2054	!-- (see: Twomey, 1959, Pure and applied Geophysics, 43)
	2055	n_act = na_init * sat**k_act
	2056	!
	2057	!-- Compute the number of cloud droplets
	2058	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
	2059	! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
	2060
	2061	!
	2062	!-- Compute activation rate after Khairoutdinov and Kogan
[2312]	2063	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
	2064	sat_max = 0.8_wp / 100.0_wp
	2065	activ = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN &
	2066	( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) / &
[2292]	2067	dt_micro
	2068
	2069	nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
	2070	ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN
	2071	!
[2312]	2072	!-- Curvature effect (afactor) with surface tension
[2292]	2073	!-- parameterization by Straka (2009)
	2074	sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
	2075	afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
	2076	!
	2077	!-- Solute effect (bfactor)
	2078	bfactor = vanthoff * molecular_weight_of_water * &
	2079	rho_s / ( molecular_weight_of_solute * rho_l )
	2080
	2081	!
[2312]	2082	!-- Prescribe power index that describes the soluble fraction
[2375]	2083	!-- of an aerosol particle (beta).
[2292]	2084	!-- (see: Morrison + Grabowski, 2007, JAS, 64)
[2375]	2085	beta_act = 0.5_wp
	2086	sigma_act = sigma_bulk**( 1.0_wp + beta_act )
[2292]	2087	!
[2312]	2088	!-- Calculate mean geometric supersaturation (s_0) with
[2292]	2089	!-- parameterization by Khvorostyanov and Curry (2006)
[2375]	2090	s_0 = dry_aerosol_radius *(- ( 1.0_wp + beta_act ) ) &
[2292]	2091	( 4.0_wp * afactor*3 / ( 27.0_wp bfactor ) )**0.5_wp
	2092
	2093	!
[2312]	2094	!-- Calculate number of activated CCN as a function of
[2292]	2095	!-- supersaturation and taking Koehler theory into account
	2096	!-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
[2375]	2097	n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &
	2098	LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
[2292]	2099	activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
[2312]	2100
	2101	nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
[2292]	2102	ENDIF
	2103
	2104	ENDDO
	2105
	2106	END SUBROUTINE activation_ij
	2107
[1682]	2108	!------------------------------------------------------------------------------!
	2109	! Description:
	2110	! ------------
[2312]	2111	!> Calculate condensation rate for cloud water content (after Khairoutdinov and
[2292]	2112	!> Kogan, 2000).
	2113	!------------------------------------------------------------------------------!
	2114	SUBROUTINE condensation_ij( i, j )
	2115
	2116	USE arrays_3d, &
	2117	ONLY: hyp, nr, pt, q, qc, qr, nc
	2118
	2119	USE cloud_parameters, &
	2120	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
	2121
	2122	USE constants, &
	2123	ONLY: pi
	2124
	2125	USE cpulog, &
	2126	ONLY: cpu_log, log_point_s
	2127
	2128	USE indices, &
	2129	ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt
	2130
	2131	USE kinds
	2132
	2133	USE control_parameters, &
	2134	ONLY: simulated_time
	2135
	2136	IMPLICIT NONE
	2137
	2138	INTEGER(iwp) :: i !<
	2139	INTEGER(iwp) :: j !<
	2140	INTEGER(iwp) :: k !<
	2141
	2142	REAL(wp) :: alpha !<
	2143	REAL(wp) :: cond !<
	2144	REAL(wp) :: cond_max !<
	2145	REAL(wp) :: dc !<
	2146	REAL(wp) :: e_s !<
	2147	REAL(wp) :: evap !<
	2148	REAL(wp) :: evap_nc !<
	2149	REAL(wp) :: g_fac !<
	2150	REAL(wp) :: nc_0 !<
	2151	REAL(wp) :: q_s !<
	2152	REAL(wp) :: sat !<
	2153	REAL(wp) :: t_l !<
	2154	REAL(wp) :: temp !<
	2155	REAL(wp) :: xc !<
	2156
	2157
	2158	DO k = nzb+1, nzt
	2159	!
	2160	!-- Actual liquid water temperature:
	2161	t_l = t_d_pt(k) * pt_1d(k)
	2162	!
	2163	!-- Saturation vapor pressure at t_l:
	2164	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	2165	( t_l - 35.86_wp ) &
	2166	)
	2167	!
	2168	!-- Computation of saturation humidity:
	2169	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
	2170	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
	2171	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / &
	2172	( 1.0_wp + alpha * q_s )
	2173
	2174	!-- Supersaturation:
	2175	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
	2176
	2177
	2178	!
	2179	!-- Actual temperature:
	2180	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
	2181
	2182	g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
	2183	l_v / ( thermal_conductivity_l * temp ) &
	2184	+ r_v * temp / ( diff_coeff_l * e_s ) &
	2185	)
	2186	!
	2187	!-- Mean weight of cloud drops
	2188	IF ( nc_1d(k) <= 0.0_wp) CYCLE
[2312]	2189	xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)
[2292]	2190	!
	2191	!-- Weight averaged diameter of cloud drops:
	2192	dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
	2193	!
	2194	!-- Integral diameter of cloud drops
[2312]	2195	nc_0 = nc_1d(k) * dc
[2292]	2196	!
	2197	!-- Condensation needs only to be calculated in supersaturated regions
	2198	IF ( sat > 0.0_wp ) THEN
	2199	!
	2200	!-- Condensation rate of cloud water content
	2201	!-- after KK scheme.
	2202	!-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
	2203	cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	2204	cond_max = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
	2205	cond = MIN( cond, cond_max / dt_micro )
[2312]	2206
[2292]	2207	qc_1d(k) = qc_1d(k) + cond * dt_micro
	2208	ELSEIF ( sat < 0.0_wp ) THEN
	2209	evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
	2210	evap = MAX( evap, -qc_1d(k) / dt_micro )
	2211
	2212	qc_1d(k) = qc_1d(k) + evap * dt_micro
	2213	ENDIF
	2214	ENDDO
	2215
	2216	END SUBROUTINE condensation_ij
	2217
	2218
	2219	!------------------------------------------------------------------------------!
	2220	! Description:
	2221	! ------------
[1682]	2222	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
	2223	!------------------------------------------------------------------------------!
[1005]	2224	SUBROUTINE autoconversion_ij( i, j )
[1000]	2225
[1320]	2226	USE arrays_3d, &
[1849]	2227	ONLY: diss, dzu
[1115]	2228
[1320]	2229	USE cloud_parameters, &
[1849]	2230	ONLY: hyrho
[1320]	2231
	2232	USE control_parameters, &
[2292]	2233	ONLY: microphysics_morrison, rho_surface
[1320]	2234
	2235	USE grid_variables, &
	2236	ONLY: dx, dy
	2237
	2238	USE indices, &
[2232]	2239	ONLY: nzb, nzt, wall_flags_0
[1320]	2240
	2241	USE kinds
	2242
[1000]	2243	IMPLICIT NONE
	2244
[1682]	2245	INTEGER(iwp) :: i !<
	2246	INTEGER(iwp) :: j !<
	2247	INTEGER(iwp) :: k !<
[1000]	2248
[2155]	2249	REAL(wp) :: alpha_cc !<
[1682]	2250	REAL(wp) :: autocon !<
	2251	REAL(wp) :: dissipation !<
[2232]	2252	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2253	REAL(wp) :: k_au !<
	2254	REAL(wp) :: l_mix !<
[2292]	2255	REAL(wp) :: nc_auto !<
[1682]	2256	REAL(wp) :: nu_c !<
	2257	REAL(wp) :: phi_au !<
	2258	REAL(wp) :: r_cc !<
	2259	REAL(wp) :: rc !<
	2260	REAL(wp) :: re_lambda !<
	2261	REAL(wp) :: sigma_cc !<
	2262	REAL(wp) :: tau_cloud !<
	2263	REAL(wp) :: xc !<
[1106]	2264
[2232]	2265	DO k = nzb+1, nzt
	2266	!
	2267	!-- Predetermine flag to mask topography
	2268	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	2269	nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
[1000]	2270
[2522]	2271	IF ( qc_1d(k) > eps_sb .AND. nc_auto > eps_mr ) THEN
[1361]	2272
	2273	k_au = k_cc / ( 20.0_wp * x0 )
[1012]	2274	!
[1048]	2275	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[1353]	2276	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
[2522]	2277	tau_cloud = MAX( 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ), &
	2278	0.0_wp )
[1012]	2279	!
[2155]	2280	!-- Universal function for autoconversion process
[1012]	2281	!-- (Seifert and Beheng, 2006):
[1361]	2282	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
[1012]	2283	!
	2284	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
[1353]	2285	!-- (Use constant nu_c = 1.0_wp instead?)
[1361]	2286	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
[1012]	2287	!
	2288	!-- Mean weight of cloud droplets:
[2292]	2289	xc = hyrho(k) * qc_1d(k) / nc_auto
[1012]	2290	!
[2155]	2291	!-- Parameterized turbulence effects on autoconversion (Seifert,
[1065]	2292	!-- Nuijens and Stevens, 2010)
[1831]	2293	IF ( collision_turbulence ) THEN
[1065]	2294	!
	2295	!-- Weight averaged radius of cloud droplets:
[1353]	2296	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
[1065]	2297
[1353]	2298	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
	2299	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
	2300	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
[1065]	2301	!
	2302	!-- Mixing length (neglecting distance to ground and stratification)
[1334]	2303	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
[1065]	2304	!
[2155]	2305	!-- Limit dissipation rate according to Seifert, Nuijens and
[1065]	2306	!-- Stevens (2010)
[1361]	2307	dissipation = MIN( 0.06_wp, diss(k,j,i) )
[1065]	2308	!
	2309	!-- Compute Taylor-microscale Reynolds number:
[1361]	2310	re_lambda = 6.0_wp / 11.0_wp * &
	2311	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
	2312	SQRT( 15.0_wp / kin_vis_air ) * &
	2313	dissipation**( 1.0_wp / 6.0_wp )
[1065]	2314	!
	2315	!-- The factor of 1.0E4 is needed to convert the dissipation rate
	2316	!-- from m2 s-3 to cm2 s-3.
[1361]	2317	k_au = k_au * ( 1.0_wp + &
	2318	dissipation * 1.0E4_wp * &
	2319	( re_lambda * 1.0E-3_wp )*0.25_wp &
	2320	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
	2321	sigma_cc )**2 &
	2322	) + beta_cc &
	2323	) &
	2324	)
[1065]	2325	ENDIF
	2326	!
[1012]	2327	!-- Autoconversion rate (Seifert and Beheng, 2006):
[1361]	2328	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
	2329	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
	2330	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
[1115]	2331	rho_surface
	2332	autocon = MIN( autocon, qc_1d(k) / dt_micro )
[1106]	2333
[2232]	2334	qr_1d(k) = qr_1d(k) + autocon * dt_micro * flag
	2335	qc_1d(k) = qc_1d(k) - autocon * dt_micro * flag
	2336	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
[2292]	2337	IF ( microphysics_morrison ) THEN
[2312]	2338	nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp * &
[2292]	2339	autocon / x0 * hyrho(k) * dt_micro * flag )
	2340	ENDIF
[1115]	2341
[1005]	2342	ENDIF
[1000]	2343
	2344	ENDDO
	2345
[1005]	2346	END SUBROUTINE autoconversion_ij
	2347
[1822]	2348	!------------------------------------------------------------------------------!
	2349	! Description:
	2350	! ------------
	2351	!> Autoconversion process (Kessler, 1969).
	2352	!------------------------------------------------------------------------------!
	2353	SUBROUTINE autoconversion_kessler_ij( i, j )
[1106]	2354
[1822]	2355	USE arrays_3d, &
[1849]	2356	ONLY: dzw, prr
[1822]	2357
	2358	USE cloud_parameters, &
[1849]	2359	ONLY: l_d_cp, pt_d_t
[1822]	2360
	2361	USE indices, &
[2317]	2362	ONLY: nzb, nzt, wall_flags_0
[1822]	2363
	2364	USE kinds
	2365
[2317]	2366	USE surface_mod, &
	2367	ONLY: get_topography_top_index
[1822]	2368
[2317]	2369
[1822]	2370	IMPLICIT NONE
	2371
[2232]	2372	INTEGER(iwp) :: i !<
	2373	INTEGER(iwp) :: j !<
	2374	INTEGER(iwp) :: k !<
	2375	INTEGER(iwp) :: k_wall !< topography top index
[1822]	2376
	2377	REAL(wp) :: dqdt_precip !<
[2232]	2378	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1822]	2379
[2232]	2380	!
	2381	!-- Determine vertical index of topography top
[2317]	2382	k_wall = get_topography_top_index( j, i, 's' )
[2232]	2383	DO k = nzb+1, nzt
	2384	!
	2385	!-- Predetermine flag to mask topography
	2386	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1822]	2387
	2388	IF ( qc_1d(k) > ql_crit ) THEN
	2389	dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
	2390	ELSE
	2391	dqdt_precip = 0.0_wp
	2392	ENDIF
	2393
[2232]	2394	qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
	2395	q_1d(k) = q_1d(k) - dqdt_precip * dt_micro * flag
	2396	pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
[1822]	2397
	2398	!
[1845]	2399	!-- Compute the rain rate (stored on surface grid point)
[2232]	2400	prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
[1822]	2401
	2402	ENDDO
	2403
	2404	END SUBROUTINE autoconversion_kessler_ij
	2405
[1682]	2406	!------------------------------------------------------------------------------!
	2407	! Description:
	2408	! ------------
	2409	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
	2410	!------------------------------------------------------------------------------!
[1005]	2411	SUBROUTINE accretion_ij( i, j )
	2412
[1320]	2413	USE arrays_3d, &
[1849]	2414	ONLY: diss
[1115]	2415
[1320]	2416	USE cloud_parameters, &
[1849]	2417	ONLY: hyrho
[1320]	2418
	2419	USE control_parameters, &
[2292]	2420	ONLY: microphysics_morrison, rho_surface
[1320]	2421
	2422	USE indices, &
[2232]	2423	ONLY: nzb, nzt, wall_flags_0
[1320]	2424
	2425	USE kinds
	2426
[1005]	2427	IMPLICIT NONE
	2428
[1682]	2429	INTEGER(iwp) :: i !<
	2430	INTEGER(iwp) :: j !<
	2431	INTEGER(iwp) :: k !<
[1005]	2432
[1682]	2433	REAL(wp) :: accr !<
[2232]	2434	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2435	REAL(wp) :: k_cr !<
[2292]	2436	REAL(wp) :: nc_accr !<
[1682]	2437	REAL(wp) :: phi_ac !<
	2438	REAL(wp) :: tau_cloud !<
[2292]	2439	REAL(wp) :: xc !<
[1320]	2440
[2292]	2441
[2232]	2442	DO k = nzb+1, nzt
	2443	!
	2444	!-- Predetermine flag to mask topography
	2445	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	2446	nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
[2232]	2447
[2292]	2448	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) .AND. &
	2449	( nc_accr > eps_mr ) ) THEN
[1012]	2450	!
[1048]	2451	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
[2155]	2452	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
[1012]	2453	!
[2155]	2454	!-- Universal function for accretion process
[1048]	2455	!-- (Seifert and Beheng, 2001):
[1361]	2456	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
[2292]	2457
[1012]	2458	!
[2292]	2459	!-- Mean weight of cloud drops
[2312]	2460	xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)
[2292]	2461	!
[2155]	2462	!-- Parameterized turbulence effects on autoconversion (Seifert,
	2463	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
[1361]	2464	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
[1831]	2465	IF ( collision_turbulence ) THEN
[1361]	2466	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
	2467	MIN( 600.0_wp, &
	2468	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
	2469	)
[1065]	2470	ELSE
[2155]	2471	k_cr = k_cr0
[1065]	2472	ENDIF
	2473	!
[1012]	2474	!-- Accretion rate (Seifert and Beheng, 2006):
[2292]	2475	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * &
	2476	SQRT( rho_surface * hyrho(k) )
[1115]	2477	accr = MIN( accr, qc_1d(k) / dt_micro )
[1106]	2478
[2232]	2479	qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
	2480	qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
[2292]	2481	IF ( microphysics_morrison ) THEN
	2482	nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc * &
	2483	hyrho(k) * dt_micro * flag &
	2484	)
	2485	ENDIF
[1115]	2486
[2292]	2487
[1005]	2488	ENDIF
[1106]	2489
[1005]	2490	ENDDO
	2491
[1000]	2492	END SUBROUTINE accretion_ij
	2493
[1005]	2494
[1682]	2495	!------------------------------------------------------------------------------!
	2496	! Description:
	2497	! ------------
	2498	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
	2499	!------------------------------------------------------------------------------!
[1005]	2500	SUBROUTINE selfcollection_breakup_ij( i, j )
	2501
[1320]	2502	USE cloud_parameters, &
[1849]	2503	ONLY: hyrho
[1320]	2504
	2505	USE control_parameters, &
[1849]	2506	ONLY: rho_surface
[1320]	2507
	2508	USE indices, &
[2232]	2509	ONLY: nzb, nzt, wall_flags_0
[1320]	2510
	2511	USE kinds
[2155]	2512
[1005]	2513	IMPLICIT NONE
	2514
[1682]	2515	INTEGER(iwp) :: i !<
	2516	INTEGER(iwp) :: j !<
	2517	INTEGER(iwp) :: k !<
[1005]	2518
[1682]	2519	REAL(wp) :: breakup !<
	2520	REAL(wp) :: dr !<
[2232]	2521	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2522	REAL(wp) :: phi_br !<
	2523	REAL(wp) :: selfcoll !<
[1320]	2524
[2232]	2525	DO k = nzb+1, nzt
	2526	!
	2527	!-- Predetermine flag to mask topography
	2528	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2529
[1115]	2530	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	2531	!
[1115]	2532	!-- Selfcollection rate (Seifert and Beheng, 2001):
[1361]	2533	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
[1012]	2534	!
[1115]	2535	!-- Weight averaged diameter of rain drops:
[1334]	2536	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	2537	!
[1048]	2538	!-- Collisional breakup rate (Seifert, 2008):
[1353]	2539	IF ( dr >= 0.3E-3_wp ) THEN
	2540	phi_br = k_br * ( dr - 1.1E-3_wp )
	2541	breakup = selfcoll * ( phi_br + 1.0_wp )
[1005]	2542	ELSE
[1353]	2543	breakup = 0.0_wp
[1005]	2544	ENDIF
[1048]	2545
[1115]	2546	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
[2232]	2547	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
[1106]	2548
[2155]	2549	ENDIF
[1005]	2550	ENDDO
	2551
	2552	END SUBROUTINE selfcollection_breakup_ij
	2553
[1106]	2554
[1682]	2555	!------------------------------------------------------------------------------!
	2556	! Description:
	2557	! ------------
[2155]	2558	!> Evaporation of precipitable water. Condensation is neglected for
[1682]	2559	!> precipitable water. Call for grid point i,j
	2560	!------------------------------------------------------------------------------!
[1012]	2561	SUBROUTINE evaporation_rain_ij( i, j )
	2562
[1320]	2563	USE arrays_3d, &
[1849]	2564	ONLY: hyp
[1048]	2565
[1320]	2566	USE cloud_parameters, &
[1849]	2567	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
[1320]	2568
	2569	USE constants, &
	2570	ONLY: pi
	2571
	2572	USE indices, &
[2232]	2573	ONLY: nzb, nzt, wall_flags_0
[1320]	2574
	2575	USE kinds
	2576
[1012]	2577	IMPLICIT NONE
	2578
[1682]	2579	INTEGER(iwp) :: i !<
	2580	INTEGER(iwp) :: j !<
	2581	INTEGER(iwp) :: k !<
[1012]	2582
[1682]	2583	REAL(wp) :: alpha !<
	2584	REAL(wp) :: dr !<
	2585	REAL(wp) :: e_s !<
	2586	REAL(wp) :: evap !<
	2587	REAL(wp) :: evap_nr !<
	2588	REAL(wp) :: f_vent !<
[2232]	2589	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2590	REAL(wp) :: g_evap !<
	2591	REAL(wp) :: lambda_r !<
	2592	REAL(wp) :: mu_r !<
	2593	REAL(wp) :: mu_r_2 !<
	2594	REAL(wp) :: mu_r_5d2 !<
	2595	REAL(wp) :: nr_0 !<
	2596	REAL(wp) :: q_s !<
	2597	REAL(wp) :: sat !<
	2598	REAL(wp) :: t_l !<
	2599	REAL(wp) :: temp !<
	2600	REAL(wp) :: xr !<
[1320]	2601
[2232]	2602	DO k = nzb+1, nzt
	2603	!
	2604	!-- Predetermine flag to mask topography
	2605	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2606
[1115]	2607	IF ( qr_1d(k) > eps_sb ) THEN
[1012]	2608	!
	2609	!-- Actual liquid water temperature:
[1115]	2610	t_l = t_d_pt(k) * pt_1d(k)
[1012]	2611	!
	2612	!-- Saturation vapor pressure at t_l:
[1361]	2613	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
	2614	( t_l - 35.86_wp ) &
	2615	)
[1012]	2616	!
	2617	!-- Computation of saturation humidity:
[1361]	2618	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
[1353]	2619	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
[1361]	2620	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
[1012]	2621	!
[1106]	2622	!-- Supersaturation:
[1361]	2623	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
[1012]	2624	!
[1361]	2625	!-- Evaporation needs only to be calculated in subsaturated regions
	2626	IF ( sat < 0.0_wp ) THEN
[1012]	2627	!
[1361]	2628	!-- Actual temperature:
	2629	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
[2155]	2630
[1361]	2631	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
	2632	( thermal_conductivity_l * temp ) + &
	2633	r_v * temp / ( diff_coeff_l * e_s ) &
	2634	)
[1012]	2635	!
[1361]	2636	!-- Mean weight of rain drops
	2637	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
[1115]	2638	!
[1361]	2639	!-- Weight averaged diameter of rain drops:
	2640	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	2641	!
[2155]	2642	!-- Compute ventilation factor and intercept parameter
[1361]	2643	!-- (Seifert and Beheng, 2006; Seifert, 2008):
	2644	IF ( ventilation_effect ) THEN
[1115]	2645	!
[1361]	2646	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	2647	!-- Stevens and Seifert, 2008):
	2648	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
	2649	!
	2650	!-- Slope parameter of gamma distribution (Seifert, 2008):
	2651	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	2652	( mu_r + 1.0_wp ) &
	2653	)**( 1.0_wp / 3.0_wp ) / dr
[1115]	2654
[1361]	2655	mu_r_2 = mu_r + 2.0_wp
[2155]	2656	mu_r_5d2 = mu_r + 2.5_wp
[1361]	2657
[2155]	2658	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
[1361]	2659	b_vent * schmidt_p_1d3 * &
	2660	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
	2661	lambda_r*( -mu_r_5d2 ) &
	2662	( 1.0_wp - &
	2663	0.5_wp * ( b_term / a_term ) * &
	2664	( lambda_r / ( c_term + lambda_r ) &
	2665	)**mu_r_5d2 - &
	2666	0.125_wp * ( b_term / a_term )*2 &
	2667	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
	2668	)**mu_r_5d2 - &
	2669	0.0625_wp * ( b_term / a_term )*3 &
	2670	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
	2671	)**mu_r_5d2 - &
[2155]	2672	0.0390625_wp * ( b_term / a_term )*4 &
[1361]	2673	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
	2674	)**mu_r_5d2 &
	2675	)
	2676
	2677	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
[2155]	2678	gamm( mu_r + 1.0_wp )
[1361]	2679	ELSE
	2680	f_vent = 1.0_wp
	2681	nr_0 = nr_1d(k) * dr
	2682	ENDIF
[1012]	2683	!
[1361]	2684	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
	2685	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
	2686	evap = MAX( evap, -qr_1d(k) / dt_micro )
	2687	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
	2688	-nr_1d(k) / dt_micro )
[1106]	2689
[2232]	2690	qr_1d(k) = qr_1d(k) + evap * dt_micro * flag
	2691	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
[1115]	2692
[1361]	2693	ENDIF
[2155]	2694	ENDIF
[1106]	2695
[1012]	2696	ENDDO
	2697
	2698	END SUBROUTINE evaporation_rain_ij
	2699
[1106]	2700
[1682]	2701	!------------------------------------------------------------------------------!
	2702	! Description:
	2703	! ------------
[2155]	2704	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
[1682]	2705	!> Call for grid point i,j
	2706	!------------------------------------------------------------------------------!
[1012]	2707	SUBROUTINE sedimentation_cloud_ij( i, j )
	2708
[1320]	2709	USE arrays_3d, &
[1849]	2710	ONLY: ddzu, dzu, prr
[1320]	2711
	2712	USE cloud_parameters, &
[1849]	2713	ONLY: hyrho, l_d_cp, pt_d_t
[1320]	2714
	2715	USE control_parameters, &
[2292]	2716	ONLY: call_microphysics_at_all_substeps, &
	2717	intermediate_timestep_count, microphysics_morrison
[1320]	2718
	2719	USE indices, &
[2232]	2720	ONLY: nzb, nzb, nzt, wall_flags_0
[1320]	2721
	2722	USE kinds
[2155]	2723
[1691]	2724	USE statistics, &
	2725	ONLY: weight_substep
	2726
[1012]	2727	IMPLICIT NONE
	2728
[1849]	2729	INTEGER(iwp) :: i !<
	2730	INTEGER(iwp) :: j !<
	2731	INTEGER(iwp) :: k !<
[1106]	2732
[2292]	2733	REAL(wp) :: flag !< flag to indicate first grid level above surface
	2734	REAL(wp) :: nc_sedi !<
	2735
	2736	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<
[2232]	2737	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
[1115]	2738
[1353]	2739	sed_qc(nzt+1) = 0.0_wp
[2375]	2740	sed_nc(nzt+1) = 0.0_wp
[1012]	2741
[2375]	2742
[2232]	2743	DO k = nzt, nzb+1, -1
	2744	!
	2745	!-- Predetermine flag to mask topography
	2746	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[2292]	2747	nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison )
	2748	!
	2749	!-- Sedimentation fluxes for number concentration are only calculated
	2750	!-- for cloud_scheme = 'morrison'
	2751	IF ( microphysics_morrison ) THEN
	2752	IF ( qc_1d(k) > eps_sb .AND. nc_1d(k) > eps_mr ) THEN
[2312]	2753	sed_nc(k) = sed_qc_const * &
[2292]	2754	( qc_1d(k) * hyrho(k) )*( 2.0_wp / 3.0_wp ) &
	2755	( nc_1d(k) )**( 1.0_wp / 3.0_wp )
	2756	ELSE
	2757	sed_nc(k) = 0.0_wp
	2758	ENDIF
[2232]	2759
[2292]	2760	sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &
	2761	nc_1d(k) / dt_micro + sed_nc(k+1) &
	2762	) * flag
	2763
	2764	nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) * &
	2765	ddzu(k+1) / hyrho(k) * dt_micro * flag
	2766	ENDIF
	2767
[2375]	2768	IF ( qc_1d(k) > eps_sb .AND. nc_sedi > eps_mr ) THEN
[2292]	2769	sed_qc(k) = sed_qc_const * nc_sedi*( -2.0_wp / 3.0_wp ) &
[2232]	2770	( qc_1d(k) * hyrho(k) )*( 5.0_wp / 3.0_wp ) flag
[1115]	2771	ELSE
[1353]	2772	sed_qc(k) = 0.0_wp
[1012]	2773	ENDIF
[1115]	2774
[1361]	2775	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
	2776	dt_micro + sed_qc(k+1) &
[2232]	2777	) * flag
[1115]	2778
[1361]	2779	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2780	hyrho(k) * dt_micro * flag
[2155]	2781	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2782	hyrho(k) * dt_micro * flag
[1361]	2783	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
[2232]	2784	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1115]	2785
[1691]	2786	!
	2787	!-- Compute the precipitation rate of cloud (fog) droplets
[1822]	2788	IF ( call_microphysics_at_all_substeps ) THEN
[2232]	2789	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * &
	2790	weight_substep(intermediate_timestep_count) * flag
[1822]	2791	ELSE
[2232]	2792	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
[1691]	2793	ENDIF
	2794
[1012]	2795	ENDDO
	2796
	2797	END SUBROUTINE sedimentation_cloud_ij
	2798
[1106]	2799
[1682]	2800	!------------------------------------------------------------------------------!
	2801	! Description:
	2802	! ------------
	2803	!> Computation of sedimentation flux. Implementation according to Stevens
	2804	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
	2805	!------------------------------------------------------------------------------!
[1012]	2806	SUBROUTINE sedimentation_rain_ij( i, j )
	2807
[1320]	2808	USE arrays_3d, &
[1849]	2809	ONLY: ddzu, dzu, prr
[1320]	2810
	2811	USE cloud_parameters, &
[1849]	2812	ONLY: hyrho, l_d_cp, pt_d_t
[1320]	2813
	2814	USE control_parameters, &
[1849]	2815	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
[1320]	2816
	2817	USE indices, &
[2232]	2818	ONLY: nzb, nzb, nzt, wall_flags_0
[1320]	2819
	2820	USE kinds
	2821
	2822	USE statistics, &
	2823	ONLY: weight_substep
[2155]	2824
[2232]	2825	USE surface_mod, &
	2826	ONLY : bc_h
[2312]	2827
[1012]	2828	IMPLICIT NONE
	2829
[2232]	2830	INTEGER(iwp) :: i !< running index x direction
	2831	INTEGER(iwp) :: j !< running index y direction
	2832	INTEGER(iwp) :: k !< running index z direction
	2833	INTEGER(iwp) :: k_run !<
	2834	INTEGER(iwp) :: m !< running index surface elements
	2835	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	2836	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1012]	2837
[1682]	2838	REAL(wp) :: c_run !<
	2839	REAL(wp) :: d_max !<
	2840	REAL(wp) :: d_mean !<
	2841	REAL(wp) :: d_min !<
	2842	REAL(wp) :: dr !<
	2843	REAL(wp) :: flux !<
[2232]	2844	REAL(wp) :: flag !< flag to indicate first grid level above surface
[1682]	2845	REAL(wp) :: lambda_r !<
	2846	REAL(wp) :: mu_r !<
	2847	REAL(wp) :: z_run !<
[1320]	2848
[1682]	2849	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
	2850	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
	2851	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
	2852	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
	2853	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
	2854	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
	2855	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
	2856	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
[1320]	2857
[1012]	2858	!
[2155]	2859	!-- Compute velocities
[2232]	2860	DO k = nzb+1, nzt
	2861	!
	2862	!-- Predetermine flag to mask topography
	2863	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2864
[1115]	2865	IF ( qr_1d(k) > eps_sb ) THEN
	2866	!
	2867	!-- Weight averaged diameter of rain drops:
[1334]	2868	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
[1115]	2869	!
	2870	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
	2871	!-- Stevens and Seifert, 2008):
[1353]	2872	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
[1115]	2873	!
	2874	!-- Slope parameter of gamma distribution (Seifert, 2008):
[1361]	2875	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
	2876	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
[1115]	2877
[1361]	2878	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	2879	a_term - b_term * ( 1.0_wp + &
	2880	c_term / lambda_r )*( -1.0_wp &
	2881	( mu_r + 1.0_wp ) ) &
	2882	) &
[2232]	2883	) * flag
[1361]	2884	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
	2885	a_term - b_term * ( 1.0_wp + &
	2886	c_term / lambda_r )*( -1.0_wp &
	2887	( mu_r + 4.0_wp ) ) &
	2888	) &
[2232]	2889	) * flag
[1065]	2890	ELSE
[1353]	2891	w_nr(k) = 0.0_wp
	2892	w_qr(k) = 0.0_wp
[1065]	2893	ENDIF
	2894	ENDDO
[1048]	2895	!
[2312]	2896	!-- Adjust boundary values using surface data type.
[2232]	2897	!-- Upward facing non-natural
[2312]	2898	surf_s = bc_h(0)%start_index(j,i)
[2232]	2899	surf_e = bc_h(0)%end_index(j,i)
	2900	DO m = surf_s, surf_e
	2901	k = bc_h(0)%k(m)
	2902	w_nr(k-1) = w_nr(k)
	2903	w_qr(k-1) = w_qr(k)
	2904	ENDDO
	2905	!
	2906	!-- Downward facing non-natural
[2312]	2907	surf_s = bc_h(1)%start_index(j,i)
[2232]	2908	surf_e = bc_h(1)%end_index(j,i)
	2909	DO m = surf_s, surf_e
	2910	k = bc_h(1)%k(m)
	2911	w_nr(k+1) = w_nr(k)
	2912	w_qr(k+1) = w_qr(k)
	2913	ENDDO
	2914	!
	2915	!-- Neumann boundary condition at model top
[1353]	2916	w_nr(nzt+1) = 0.0_wp
	2917	w_qr(nzt+1) = 0.0_wp
[1065]	2918	!
	2919	!-- Compute Courant number
[2232]	2920	DO k = nzb+1, nzt
	2921	!
	2922	!-- Predetermine flag to mask topography
	2923	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2924
[1361]	2925	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
[2232]	2926	dt_micro * ddzu(k) * flag
[1361]	2927	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
[2232]	2928	dt_micro * ddzu(k) * flag
[2155]	2929	ENDDO
[1065]	2930	!
	2931	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
	2932	IF ( limiter_sedimentation ) THEN
	2933
[2232]	2934	DO k = nzb+1, nzt
	2935	!
	2936	!-- Predetermine flag to mask topography
	2937	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2938
[1646]	2939	d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
[1115]	2940	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
	2941	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
[1065]	2942
[1361]	2943	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	2944	2.0_wp * d_max, &
[2232]	2945	ABS( d_mean ) ) * flag
[1065]	2946
[1646]	2947	d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
[1115]	2948	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
	2949	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
[1065]	2950
[1361]	2951	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
	2952	2.0_wp * d_max, &
[2232]	2953	ABS( d_mean ) ) * flag
[1022]	2954	ENDDO
[1048]	2955
[1065]	2956	ELSE
[1106]	2957
[1353]	2958	nr_slope = 0.0_wp
	2959	qr_slope = 0.0_wp
[1106]	2960
[1065]	2961	ENDIF
[1115]	2962
[1353]	2963	sed_nr(nzt+1) = 0.0_wp
	2964	sed_qr(nzt+1) = 0.0_wp
[1065]	2965	!
	2966	!-- Compute sedimentation flux
[2232]	2967	DO k = nzt, nzb+1, -1
[1065]	2968	!
[2232]	2969	!-- Predetermine flag to mask topography
	2970	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	2971	!
[2312]	2972	!-- Sum up all rain drop number densities which contribute to the flux
[1065]	2973	!-- through k-1/2
[1353]	2974	flux = 0.0_wp
	2975	z_run = 0.0_wp ! height above z(k)
[1065]	2976	k_run = k
[1346]	2977	c_run = MIN( 1.0_wp, c_nr(k) )
[1353]	2978	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1361]	2979	flux = flux + hyrho(k_run) * &
	2980	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	2981	0.5_wp ) * c_run * dzu(k_run) * flag
	2982	z_run = z_run + dzu(k_run) * flag
	2983	k_run = k_run + 1 * flag
	2984	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
[1022]	2985	ENDDO
	2986	!
[2155]	2987	!-- It is not allowed to sediment more rain drop number density than
[1065]	2988	!-- available
[1361]	2989	flux = MIN( flux, &
[1115]	2990	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
[1065]	2991
[2232]	2992	sed_nr(k) = flux / dt_micro * flag
[1361]	2993	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
[2232]	2994	hyrho(k) * dt_micro * flag
[1065]	2995	!
[2155]	2996	!-- Sum up all rain water content which contributes to the flux
[1065]	2997	!-- through k-1/2
[1353]	2998	flux = 0.0_wp
	2999	z_run = 0.0_wp ! height above z(k)
[1065]	3000	k_run = k
[1346]	3001	c_run = MIN( 1.0_wp, c_qr(k) )
[1106]	3002
[1361]	3003	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
[1106]	3004
[1361]	3005	flux = flux + hyrho(k_run) * &
	3006	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
[2232]	3007	0.5_wp ) * c_run * dzu(k_run) * flag
	3008	z_run = z_run + dzu(k_run) * flag
	3009	k_run = k_run + 1 * flag
	3010	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
[1106]	3011
[1065]	3012	ENDDO
	3013	!
	3014	!-- It is not allowed to sediment more rain water content than available
[1361]	3015	flux = MIN( flux, &
[1115]	3016	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
[1065]	3017
[2232]	3018	sed_qr(k) = flux / dt_micro * flag
[1115]	3019
[1361]	3020	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	3021	hyrho(k) * dt_micro * flag
[1361]	3022	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	3023	hyrho(k) * dt_micro * flag
[1361]	3024	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
[2232]	3025	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
[1065]	3026	!
	3027	!-- Compute the rain rate
[1361]	3028	IF ( call_microphysics_at_all_substeps ) THEN
[1691]	3029	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
[2232]	3030	* weight_substep(intermediate_timestep_count) * flag
[1361]	3031	ELSE
[2232]	3032	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
[1361]	3033	ENDIF
	3034
[1065]	3035	ENDDO
[1115]	3036
[1691]	3037	END SUBROUTINE sedimentation_rain_ij
[1012]	3038
[1691]	3039
	3040	!------------------------------------------------------------------------------!
	3041	! Description:
	3042	! ------------
	3043	!> This subroutine computes the precipitation amount due to gravitational
	3044	!> settling of rain and cloud (fog) droplets
	3045	!------------------------------------------------------------------------------!
	3046	SUBROUTINE calc_precipitation_amount_ij( i, j )
	3047
[1849]	3048	USE arrays_3d, &
	3049	ONLY: precipitation_amount, prr
	3050
[1691]	3051	USE cloud_parameters, &
[1849]	3052	ONLY: hyrho
[1691]	3053
	3054	USE control_parameters, &
	3055	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
	3056	intermediate_timestep_count, intermediate_timestep_count_max,&
[1822]	3057	precipitation_amount_interval, time_do2d_xy
[1691]	3058
	3059	USE indices, &
[2232]	3060	ONLY: nzb, nzt, wall_flags_0
[1691]	3061
	3062	USE kinds
	3063
[2232]	3064	USE surface_mod, &
	3065	ONLY : bc_h
	3066
[1691]	3067	IMPLICIT NONE
	3068
[2232]	3069	INTEGER(iwp) :: i !< running index x direction
	3070	INTEGER(iwp) :: j !< running index y direction
	3071	INTEGER(iwp) :: k !< running index z direction
	3072	INTEGER(iwp) :: m !< running index surface elements
	3073	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	3074	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1691]	3075
	3076	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
	3077	( .NOT. call_microphysics_at_all_substeps .OR. &
	3078	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
	3079	THEN
	3080
[2312]	3081	surf_s = bc_h(0)%start_index(j,i)
[2232]	3082	surf_e = bc_h(0)%end_index(j,i)
	3083	DO m = surf_s, surf_e
	3084	k = bc_h(0)%k(m)
	3085	precipitation_amount(j,i) = precipitation_amount(j,i) + &
	3086	prr(k,j,i) * hyrho(k) * dt_3d
	3087	ENDDO
	3088
[1048]	3089	ENDIF
	3090
[1691]	3091	END SUBROUTINE calc_precipitation_amount_ij
[1012]	3092
[1361]	3093	!------------------------------------------------------------------------------!
[1682]	3094	! Description:
	3095	! ------------
	3096	!> This function computes the gamma function (Press et al., 1992).
[2155]	3097	!> The gamma function is needed for the calculation of the evaporation
[1682]	3098	!> of rain drops.
[1361]	3099	!------------------------------------------------------------------------------!
[2155]	3100	FUNCTION gamm( xx )
[1320]	3101
	3102	USE kinds
	3103
[2155]	3104	IMPLICIT NONE
[1106]	3105
[1682]	3106	INTEGER(iwp) :: j !<
[1320]	3107
[1682]	3108	REAL(wp) :: gamm !<
	3109	REAL(wp) :: ser !<
	3110	REAL(wp) :: tmp !<
	3111	REAL(wp) :: x_gamm !<
	3112	REAL(wp) :: xx !<
	3113	REAL(wp) :: y_gamm !<
[1320]	3114
[1849]	3115
	3116	REAL(wp), PARAMETER :: stp = 2.5066282746310005_wp !<
	3117	REAL(wp), PARAMETER :: cof(6) = (/ 76.18009172947146_wp, &
	3118	-86.50532032941677_wp, &
	3119	24.01409824083091_wp, &
	3120	-1.231739572450155_wp, &
	3121	0.1208650973866179E-2_wp, &
	3122	-0.5395239384953E-5_wp /) !<
	3123
	3124	x_gamm = xx
	3125	y_gamm = x_gamm
[1353]	3126	tmp = x_gamm + 5.5_wp
[1849]	3127	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
[1334]	3128	ser = 1.000000000190015_wp
[1106]	3129
[2155]	3130	DO j = 1, 6
	3131	y_gamm = y_gamm + 1.0_wp
	3132	ser = ser + cof( j ) / y_gamm
[1106]	3133	ENDDO
	3134
[2155]	3135	!
	3136	!-- Until this point the algorithm computes the logarithm of the gamma
	3137	!-- function. Hence, the exponential function is used.
	3138	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
	3139	gamm = EXP( tmp ) * stp * ser / x_gamm
[1106]	3140
[2155]	3141	RETURN
[1012]	3142
[2155]	3143	END FUNCTION gamm
	3144
[1012]	3145	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |