Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |