Home

Context Navigation

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

Last change on this file since 2297 was 2292, checked in by schwenkel, 7 years ago
implementation of new bulk microphysics scheme
Property svn:keywords set to `Id`
File size: 123.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |