Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |