Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |