Home

Context Navigation

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

Last change on this file since 2946 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 123.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |