Home

Context Navigation

microphysics.f90 @ 1849

Last change on this file since 1849 was 1849, checked in by hoffmann, 8 years ago
lpm_droplet_condensation improved, microphysics partially modularized
Property svn:keywords set to `Id`
File size: 86.6 KB

Line
1	!> @file microphysics.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms
6	! of the GNU General Public License as published by the Free Software Foundation,
7	! either version 3 of the License, or (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with
14	! PALM. If not, see <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2016 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! ------------------
21	! Adapted for modularization of microphysics.
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: microphysics.f90 1849 2016-04-08 11:33:18Z hoffmann $
26	!
27	! 1845 2016-04-08 08:29:13Z raasch
28	! nzb_2d replaced by nzb_s_inner, Kessler precipitation is stored at surface
29	! point (instead of one point above surface)
30	!
31	! 1831 2016-04-07 13:15:51Z hoffmann
32	! turbulence renamed collision_turbulence,
33	! drizzle renamed cloud_water_sedimentation. cloud_water_sedimentation also
34	! avaialble for microphysics_kessler.
35	!
36	! 1822 2016-04-07 07:49:42Z hoffmann
37	! Unused variables removed.
38	! Kessler scheme integrated.
39	!
40	! 1691 2015-10-26 16:17:44Z maronga
41	! Added new routine calc_precipitation_amount. The routine now allows to account
42	! for precipitation due to sedimenation of cloud (fog) droplets
43	!
44	! 1682 2015-10-07 23:56:08Z knoop
45	! Code annotations made doxygen readable
46	!
47	! 1646 2015-09-02 16:00:10Z hoffmann
48	! Bugfix: Wrong computation of d_mean.
49	!
50	! 1361 2014-04-16 15:17:48Z hoffmann
51	! Bugfix in sedimentation_rain: Index corrected.
52	! Vectorized version of adjust_cloud added.
53	! Little reformatting of the code.
54	!
55	! 1353 2014-04-08 15:21:23Z heinze
56	! REAL constants provided with KIND-attribute
57	!
58	! 1346 2014-03-27 13:18:20Z heinze
59	! Bugfix: REAL constants provided with KIND-attribute especially in call of
60	! intrinsic function like MAX, MIN, SIGN
61	!
62	! 1334 2014-03-25 12:21:40Z heinze
63	! Bugfix: REAL constants provided with KIND-attribute
64	!
65	! 1322 2014-03-20 16:38:49Z raasch
66	! REAL constants defined as wp-kind
67	!
68	! 1320 2014-03-20 08:40:49Z raasch
69	! ONLY-attribute added to USE-statements,
70	! kind-parameters added to all INTEGER and REAL declaration statements,
71	! kinds are defined in new module kinds,
72	! comment fields (!:) to be used for variable explanations added to
73	! all variable declaration statements
74	!
75	! 1241 2013-10-30 11:36:58Z heinze
76	! hyp and rho have to be calculated at each time step if data from external
77	! file LSF_DATA are used
78	!
79	! 1115 2013-03-26 18:16:16Z hoffmann
80	! microphyical tendencies are calculated in microphysics_control in an optimized
81	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
82	!
83	! 1106 2013-03-04 05:31:38Z raasch
84	! small changes in code formatting
85	!
86	! 1092 2013-02-02 11:24:22Z raasch
87	! unused variables removed
88	! file put under GPL
89	!
90	! 1065 2012-11-22 17:42:36Z hoffmann
91	! Sedimentation process implemented according to Stevens and Seifert (2008).
92	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
93	! and Stevens, 2010).
94	!
95	! 1053 2012-11-13 17:11:03Z hoffmann
96	! initial revision
97	!
98	! Description:
99	! ------------
100	!> Calculate bilk cloud microphysics.
101	!------------------------------------------------------------------------------!
102	MODULE microphysics_mod
103
104	USE kinds
105
106	IMPLICIT NONE
107
108	LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation
109	LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter
110	LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects
111	LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect
112
113	REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3)
114	REAL(wp) :: a_2 = -7.38E-5_wp !< coef. in turb. parametrization (cm-2 s3)
115	REAL(wp) :: a_3 = -1.40E-2_wp !< coef. in turb. parametrization
116	REAL(wp) :: a_term = 9.65_wp !< coef. for terminal velocity (m s-1)
117	REAL(wp) :: a_vent = 0.78_wp !< coef. for ventilation effect
118	REAL(wp) :: b_1 = 11.45E-6_wp !< coef. in turb. parametrization (m)
119	REAL(wp) :: b_2 = 9.68E-6_wp !< coef. in turb. parametrization (m)
120	REAL(wp) :: b_3 = 0.62_wp !< coef. in turb. parametrization
121	REAL(wp) :: b_term = 9.8_wp !< coef. for terminal velocity (m s-1)
122	REAL(wp) :: b_vent = 0.308_wp !< coef. for ventilation effect
123	REAL(wp) :: beta_cc = 3.09E-4_wp !< coef. in turb. parametrization (cm-2 s3)
124	REAL(wp) :: c_1 = 4.82E-6_wp !< coef. in turb. parametrization (m)
125	REAL(wp) :: c_2 = 4.8E-6_wp !< coef. in turb. parametrization (m)
126	REAL(wp) :: c_3 = 0.76_wp !< coef. in turb. parametrization
127	REAL(wp) :: c_const = 0.93_wp !< const. in Taylor-microscale Reynolds number
128	REAL(wp) :: c_evap = 0.7_wp !< constant in evaporation
129	REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1)
130	REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1)
131	REAL(wp) :: eps_sb = 1.0E-20_wp !< threshold in two-moments scheme
132	REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1)
133	REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1)
134	REAL(wp) :: k_rr = 7.12_wp !< const. rain-rain kernel (m3 kg-1 s-1)
135	REAL(wp) :: k_br = 1000.0_wp !< const. in breakup parametrization (m-1)
136	REAL(wp) :: k_st = 1.2E8_wp !< const. in drizzle parametrization (m-1 s-1)
137	REAL(wp) :: kappa_rr = 60.7_wp !< const. in collision kernel (kg-1/3)
138	REAL(wp) :: kin_vis_air = 1.4086E-5_wp !< kin. viscosity of air (m2 s-1)
139	REAL(wp) :: prec_time_const = 0.001_wp !< coef. in Kessler scheme (s-1)
140	REAL(wp) :: ql_crit = 0.0005_wp !< coef. in Kessler scheme (kg kg-1)
141	REAL(wp) :: schmidt_p_1d3=0.8921121_wp !< Schmidt number0.33333, 0.710.33333
142	REAL(wp) :: sigma_gc = 1.3_wp !< geometric standard deviation cloud droplets
143	REAL(wp) :: thermal_conductivity_l = 2.43E-2_wp !< therm. cond. air (J m-1 s-1 K-1)
144	REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1)
145	REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg)
146	REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg)
147	REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg)
148
149	REAL(wp) :: c_sedimentation = 2.0_wp !< Courant number of sedimentation process
150	REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l )
151	REAL(wp) :: dt_micro !< microphysics time step
152	REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration
153	REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s)
154	REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water
155	REAL(wp) :: pirho_l !< pi * rho_l / 6.0;
156
157	REAL(wp), DIMENSION(:), ALLOCATABLE :: nc_1d !<
158	REAL(wp), DIMENSION(:), ALLOCATABLE :: nr_1d !<
159	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_1d !<
160	REAL(wp), DIMENSION(:), ALLOCATABLE :: qc_1d !<
161	REAL(wp), DIMENSION(:), ALLOCATABLE :: qr_1d !<
162	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_1d !<
163
164	SAVE
165
166	PRIVATE
167	PUBLIC microphysics_control, microphysics_init
168
169	PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation, &
170	dt_precipitation, limiter_sedimentation, nc_const, sigma_gc, &
171	ventilation_effect
172
173	INTERFACE microphysics_control
174	MODULE PROCEDURE microphysics_control
175	MODULE PROCEDURE microphysics_control_ij
176	END INTERFACE microphysics_control
177
178	INTERFACE adjust_cloud
179	MODULE PROCEDURE adjust_cloud
180	MODULE PROCEDURE adjust_cloud_ij
181	END INTERFACE adjust_cloud
182
183	INTERFACE autoconversion
184	MODULE PROCEDURE autoconversion
185	MODULE PROCEDURE autoconversion_ij
186	END INTERFACE autoconversion
187
188	INTERFACE autoconversion_kessler
189	MODULE PROCEDURE autoconversion_kessler
190	MODULE PROCEDURE autoconversion_kessler_ij
191	END INTERFACE autoconversion_kessler
192
193	INTERFACE accretion
194	MODULE PROCEDURE accretion
195	MODULE PROCEDURE accretion_ij
196	END INTERFACE accretion
197
198	INTERFACE selfcollection_breakup
199	MODULE PROCEDURE selfcollection_breakup
200	MODULE PROCEDURE selfcollection_breakup_ij
201	END INTERFACE selfcollection_breakup
202
203	INTERFACE evaporation_rain
204	MODULE PROCEDURE evaporation_rain
205	MODULE PROCEDURE evaporation_rain_ij
206	END INTERFACE evaporation_rain
207
208	INTERFACE sedimentation_cloud
209	MODULE PROCEDURE sedimentation_cloud
210	MODULE PROCEDURE sedimentation_cloud_ij
211	END INTERFACE sedimentation_cloud
212
213	INTERFACE sedimentation_rain
214	MODULE PROCEDURE sedimentation_rain
215	MODULE PROCEDURE sedimentation_rain_ij
216	END INTERFACE sedimentation_rain
217
218	INTERFACE calc_precipitation_amount
219	MODULE PROCEDURE calc_precipitation_amount
220	MODULE PROCEDURE calc_precipitation_amount_ij
221	END INTERFACE calc_precipitation_amount
222
223	CONTAINS
224	!------------------------------------------------------------------------------!
225	! Description:
226	! ------------
227	!> Initialization of bulk microphysics
228	!------------------------------------------------------------------------------!
229	SUBROUTINE microphysics_init
230
231	USE arrays_3d, &
232	ONLY: dzu
233
234	USE constants, &
235	ONLY: pi
236
237	USE cloud_parameters, &
238	ONLY: rho_l
239
240	USE control_parameters, &
241	ONLY: microphysics_seifert
242
243	USE indices, &
244	ONLY: nzb, nzt
245
246	IMPLICIT NONE
247
248	!
249	!-- constant for the sedimentation of cloud water (2-moment cloud physics)
250	sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) &
251	)*( 2.0_wp / 3.0_wp ) &
252	EXP( 5.0_wp * LOG( sigma_gc )**2 )
253
254	!
255	!-- Calculate timestep according to precipitation
256	IF ( microphysics_seifert ) THEN
257	dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) / &
258	w_precipitation
259	ENDIF
260
261	!
262	!-- Pre-calculate frequently calculated fractions of pi and rho_l
263	pirho_l = pi * rho_l / 6.0_wp
264	dpirho_l = 1.0_wp / pirho_l
265
266	!
267	!-- Allocate 1D microphysics arrays
268	ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1), &
269	qc_1d(nzb:nzt+1) )
270
271	IF ( microphysics_seifert ) THEN
272	ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) )
273	ENDIF
274
275	!
276	!-- Initialze nc_1d with nc_const
277	nc_1d = nc_const
278
279	END SUBROUTINE microphysics_init
280
281
282	!------------------------------------------------------------------------------!
283	! Description:
284	! ------------
285	!> Control of microphysics for all grid points
286	!------------------------------------------------------------------------------!
287	SUBROUTINE microphysics_control
288
289	USE arrays_3d, &
290	ONLY: hyp, pt_init, prr, zu
291
292	USE cloud_parameters, &
293	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
294
295	USE control_parameters, &
296	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
297	intermediate_timestep_count, large_scale_forcing, &
298	lsf_surf, microphysics_kessler, microphysics_seifert, &
299	pt_surface, rho_surface,surface_pressure
300
301	USE indices, &
302	ONLY: nzb, nzt
303
304	USE kinds
305
306	USE statistics, &
307	ONLY: weight_pres
308
309	IMPLICIT NONE
310
311	INTEGER(iwp) :: k !<
312
313	REAL(wp) :: t_surface !<
314
315	IF ( large_scale_forcing .AND. lsf_surf ) THEN
316	!
317	!-- Calculate:
318	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
319	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
320	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
321	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
322	DO k = nzb, nzt+1
323	hyp(k) = surface_pressure * 100.0_wp * &
324	( ( t_surface - g / cp * zu(k) ) / &
325	t_surface )**(1.0_wp / 0.286_wp)
326	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
327	t_d_pt(k) = 1.0_wp / pt_d_t(k)
328	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
329	ENDDO
330
331	!
332	!-- Compute reference density
333	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
334	ENDIF
335
336	!
337	!-- Compute length of time step
338	IF ( call_microphysics_at_all_substeps ) THEN
339	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
340	ELSE
341	dt_micro = dt_3d
342	ENDIF
343
344	!
345	!-- Reset precipitation rate
346	IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp
347
348	!
349	!-- Compute cloud physics
350	IF ( microphysics_kessler ) THEN
351
352	CALL autoconversion_kessler
353	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
354
355	ELSEIF ( microphysics_seifert ) THEN
356
357	CALL adjust_cloud
358	CALL autoconversion
359	CALL accretion
360	CALL selfcollection_breakup
361	CALL evaporation_rain
362	CALL sedimentation_rain
363	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud
364
365	ENDIF
366
367	CALL calc_precipitation_amount
368
369	END SUBROUTINE microphysics_control
370
371	!------------------------------------------------------------------------------!
372	! Description:
373	! ------------
374	!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
375	!> evaporation of rain drops due to too small or too big weights
376	!> of rain drops (Stevens and Seifert, 2008).
377	!------------------------------------------------------------------------------!
378	SUBROUTINE adjust_cloud
379
380	USE arrays_3d, &
381	ONLY: qr, nr
382
383	USE cloud_parameters, &
384	ONLY: hyrho
385
386	USE cpulog, &
387	ONLY: cpu_log, log_point_s
388
389	USE indices, &
390	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
391
392	USE kinds
393
394	IMPLICIT NONE
395
396	INTEGER(iwp) :: i !<
397	INTEGER(iwp) :: j !<
398	INTEGER(iwp) :: k !<
399
400	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
401
402	DO i = nxl, nxr
403	DO j = nys, nyn
404	DO k = nzb_s_inner(j,i)+1, nzt
405	IF ( qr(k,j,i) <= eps_sb ) THEN
406	qr(k,j,i) = 0.0_wp
407	nr(k,j,i) = 0.0_wp
408	ELSE
409	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
410	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
411	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
412	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
413	ENDIF
414	ENDIF
415	ENDDO
416	ENDDO
417	ENDDO
418
419	CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' )
420
421	END SUBROUTINE adjust_cloud
422
423
424	!------------------------------------------------------------------------------!
425	! Description:
426	! ------------
427	!> Autoconversion rate (Seifert and Beheng, 2006).
428	!------------------------------------------------------------------------------!
429	SUBROUTINE autoconversion
430
431	USE arrays_3d, &
432	ONLY: diss, dzu, nr, qc, qr
433
434	USE cloud_parameters, &
435	ONLY: hyrho
436
437	USE control_parameters, &
438	ONLY: rho_surface
439
440	USE cpulog, &
441	ONLY: cpu_log, log_point_s
442
443	USE grid_variables, &
444	ONLY: dx, dy
445
446	USE indices, &
447	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
448
449	USE kinds
450
451	IMPLICIT NONE
452
453	INTEGER(iwp) :: i !<
454	INTEGER(iwp) :: j !<
455	INTEGER(iwp) :: k !<
456
457	REAL(wp) :: alpha_cc !<
458	REAL(wp) :: autocon !<
459	REAL(wp) :: dissipation !<
460	REAL(wp) :: k_au !<
461	REAL(wp) :: l_mix !<
462	REAL(wp) :: nu_c !<
463	REAL(wp) :: phi_au !<
464	REAL(wp) :: r_cc !<
465	REAL(wp) :: rc !<
466	REAL(wp) :: re_lambda !<
467	REAL(wp) :: sigma_cc !<
468	REAL(wp) :: tau_cloud !<
469	REAL(wp) :: xc !<
470
471	CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
472
473	DO i = nxl, nxr
474	DO j = nys, nyn
475	DO k = nzb_s_inner(j,i)+1, nzt
476
477	IF ( qc(k,j,i) > eps_sb ) THEN
478
479	k_au = k_cc / ( 20.0_wp * x0 )
480	!
481	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
482	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
483	tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
484	!
485	!-- Universal function for autoconversion process
486	!-- (Seifert and Beheng, 2006):
487	phi_au = 600.0_wp * tau_cloud*0.68_wp &
488	( 1.0_wp - tau_cloud0.68_wp )3
489	!
490	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
491	!-- (Use constant nu_c = 1.0_wp instead?)
492	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
493	!
494	!-- Mean weight of cloud droplets:
495	xc = hyrho(k) * qc(k,j,i) / nc_const
496	!
497	!-- Parameterized turbulence effects on autoconversion (Seifert,
498	!-- Nuijens and Stevens, 2010)
499	IF ( collision_turbulence ) THEN
500	!
501	!-- Weight averaged radius of cloud droplets:
502	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
503
504	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
505	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
506	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
507	!
508	!-- Mixing length (neglecting distance to ground and
509	!-- stratification)
510	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
511	!
512	!-- Limit dissipation rate according to Seifert, Nuijens and
513	!-- Stevens (2010)
514	dissipation = MIN( 0.06_wp, diss(k,j,i) )
515	!
516	!-- Compute Taylor-microscale Reynolds number:
517	re_lambda = 6.0_wp / 11.0_wp * &
518	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
519	SQRT( 15.0_wp / kin_vis_air ) * &
520	dissipation**( 1.0_wp / 6.0_wp )
521	!
522	!-- The factor of 1.0E4 is needed to convert the dissipation
523	!-- rate from m2 s-3 to cm2 s-3.
524	k_au = k_au * ( 1.0_wp + &
525	dissipation * 1.0E4_wp * &
526	( re_lambda * 1.0E-3_wp )*0.25_wp &
527	( alpha_cc * EXP( -1.0_wp * ( ( rc - &
528	r_cc ) / &
529	sigma_cc )**2 &
530	) + beta_cc &
531	) &
532	)
533	ENDIF
534	!
535	!-- Autoconversion rate (Seifert and Beheng, 2006):
536	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
537	( nu_c + 1.0_wp )*2 qc(k,j,i)*2 xc*2 &
538	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
539	rho_surface
540	autocon = MIN( autocon, qc(k,j,i) / dt_micro )
541
542	qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
543	qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
544	nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
545
546	ENDIF
547
548	ENDDO
549	ENDDO
550	ENDDO
551
552	CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' )
553
554	END SUBROUTINE autoconversion
555
556
557	!------------------------------------------------------------------------------!
558	! Description:
559	! ------------
560	!> Autoconversion process (Kessler, 1969).
561	!------------------------------------------------------------------------------!
562	SUBROUTINE autoconversion_kessler
563
564	USE arrays_3d, &
565	ONLY: dzw, pt, prr, q, qc
566
567	USE cloud_parameters, &
568	ONLY: l_d_cp, pt_d_t
569
570	USE indices, &
571	ONLY: nxl, nxr, nyn, nys, nzb_s_inner, nzt
572
573	USE kinds
574
575
576	IMPLICIT NONE
577
578	INTEGER(iwp) :: i !<
579	INTEGER(iwp) :: j !<
580	INTEGER(iwp) :: k !<
581
582	REAL(wp) :: dqdt_precip !<
583
584	DO i = nxl, nxr
585	DO j = nys, nyn
586	DO k = nzb_s_inner(j,i)+1, nzt
587
588	IF ( qc(k,j,i) > ql_crit ) THEN
589	dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
590	ELSE
591	dqdt_precip = 0.0_wp
592	ENDIF
593
594	qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro
595	q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro
596	pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * &
597	pt_d_t(k)
598
599	!
600	!-- Compute the rain rate (stored on surface grid point)
601	prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) + &
602	dqdt_precip * dzw(k)
603
604	ENDDO
605	ENDDO
606	ENDDO
607
608	END SUBROUTINE autoconversion_kessler
609
610
611	!------------------------------------------------------------------------------!
612	! Description:
613	! ------------
614	!> Accretion rate (Seifert and Beheng, 2006).
615	!------------------------------------------------------------------------------!
616	SUBROUTINE accretion
617
618	USE arrays_3d, &
619	ONLY: diss, qc, qr
620
621	USE cloud_parameters, &
622	ONLY: hyrho
623
624	USE control_parameters, &
625	ONLY: rho_surface
626
627	USE cpulog, &
628	ONLY: cpu_log, log_point_s
629
630	USE indices, &
631	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
632
633	USE kinds
634
635	IMPLICIT NONE
636
637	INTEGER(iwp) :: i !<
638	INTEGER(iwp) :: j !<
639	INTEGER(iwp) :: k !<
640
641	REAL(wp) :: accr !<
642	REAL(wp) :: k_cr !<
643	REAL(wp) :: phi_ac !<
644	REAL(wp) :: tau_cloud !<
645
646	CALL cpu_log( log_point_s(56), 'accretion', 'start' )
647
648	DO i = nxl, nxr
649	DO j = nys, nyn
650	DO k = nzb_s_inner(j,i)+1, nzt
651
652	IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN
653	!
654	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
655	tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
656	!
657	!-- Universal function for accretion process (Seifert and
658	!-- Beheng, 2001):
659	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
660	!
661	!-- Parameterized turbulence effects on autoconversion (Seifert,
662	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
663	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
664	IF ( collision_turbulence ) THEN
665	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
666	MIN( 600.0_wp, &
667	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
668	)
669	ELSE
670	k_cr = k_cr0
671	ENDIF
672	!
673	!-- Accretion rate (Seifert and Beheng, 2006):
674	accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &
675	SQRT( rho_surface * hyrho(k) )
676	accr = MIN( accr, qc(k,j,i) / dt_micro )
677
678	qr(k,j,i) = qr(k,j,i) + accr * dt_micro
679	qc(k,j,i) = qc(k,j,i) - accr * dt_micro
680
681	ENDIF
682
683	ENDDO
684	ENDDO
685	ENDDO
686
687	CALL cpu_log( log_point_s(56), 'accretion', 'stop' )
688
689	END SUBROUTINE accretion
690
691
692	!------------------------------------------------------------------------------!
693	! Description:
694	! ------------
695	!> Collisional breakup rate (Seifert, 2008).
696	!------------------------------------------------------------------------------!
697	SUBROUTINE selfcollection_breakup
698
699	USE arrays_3d, &
700	ONLY: nr, qr
701
702	USE cloud_parameters, &
703	ONLY: hyrho
704
705	USE control_parameters, &
706	ONLY: rho_surface
707
708	USE cpulog, &
709	ONLY: cpu_log, log_point_s
710
711	USE indices, &
712	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
713
714	USE kinds
715
716	IMPLICIT NONE
717
718	INTEGER(iwp) :: i !<
719	INTEGER(iwp) :: j !<
720	INTEGER(iwp) :: k !<
721
722	REAL(wp) :: breakup !<
723	REAL(wp) :: dr !<
724	REAL(wp) :: phi_br !<
725	REAL(wp) :: selfcoll !<
726
727	CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
728
729	DO i = nxl, nxr
730	DO j = nys, nyn
731	DO k = nzb_s_inner(j,i)+1, nzt
732	IF ( qr(k,j,i) > eps_sb ) THEN
733	!
734	!-- Selfcollection rate (Seifert and Beheng, 2001):
735	selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * &
736	SQRT( hyrho(k) * rho_surface )
737	!
738	!-- Weight averaged diameter of rain drops:
739	dr = ( hyrho(k) * qr(k,j,i) / &
740	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
741	!
742	!-- Collisional breakup rate (Seifert, 2008):
743	IF ( dr >= 0.3E-3_wp ) THEN
744	phi_br = k_br * ( dr - 1.1E-3_wp )
745	breakup = selfcoll * ( phi_br + 1.0_wp )
746	ELSE
747	breakup = 0.0_wp
748	ENDIF
749
750	selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
751	nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
752
753	ENDIF
754	ENDDO
755	ENDDO
756	ENDDO
757
758	CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' )
759
760	END SUBROUTINE selfcollection_breakup
761
762
763	!------------------------------------------------------------------------------!
764	! Description:
765	! ------------
766	!> Evaporation of precipitable water. Condensation is neglected for
767	!> precipitable water.
768	!------------------------------------------------------------------------------!
769	SUBROUTINE evaporation_rain
770
771	USE arrays_3d, &
772	ONLY: hyp, nr, pt, q, qc, qr
773
774	USE cloud_parameters, &
775	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
776
777	USE constants, &
778	ONLY: pi
779
780	USE cpulog, &
781	ONLY: cpu_log, log_point_s
782
783	USE indices, &
784	ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt
785
786	USE kinds
787
788	IMPLICIT NONE
789
790	INTEGER(iwp) :: i !<
791	INTEGER(iwp) :: j !<
792	INTEGER(iwp) :: k !<
793
794	REAL(wp) :: alpha !<
795	REAL(wp) :: dr !<
796	REAL(wp) :: e_s !<
797	REAL(wp) :: evap !<
798	REAL(wp) :: evap_nr !<
799	REAL(wp) :: f_vent !<
800	REAL(wp) :: g_evap !<
801	REAL(wp) :: lambda_r !<
802	REAL(wp) :: mu_r !<
803	REAL(wp) :: mu_r_2 !<
804	REAL(wp) :: mu_r_5d2 !<
805	REAL(wp) :: nr_0 !<
806	REAL(wp) :: q_s !<
807	REAL(wp) :: sat !<
808	REAL(wp) :: t_l !<
809	REAL(wp) :: temp !<
810	REAL(wp) :: xr !<
811
812	CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
813
814	DO i = nxl, nxr
815	DO j = nys, nyn
816	DO k = nzb_s_inner(j,i)+1, nzt
817	IF ( qr(k,j,i) > eps_sb ) THEN
818	!
819	!-- Actual liquid water temperature:
820	t_l = t_d_pt(k) * pt(k,j,i)
821	!
822	!-- Saturation vapor pressure at t_l:
823	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
824	( t_l - 35.86_wp ) &
825	)
826	!
827	!-- Computation of saturation humidity:
828	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
829	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
830	q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / &
831	( 1.0_wp + alpha * q_s )
832	!
833	!-- Supersaturation:
834	sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
835	!
836	!-- Evaporation needs only to be calculated in subsaturated regions
837	IF ( sat < 0.0_wp ) THEN
838	!
839	!-- Actual temperature:
840	temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
841
842	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &
843	l_v / ( thermal_conductivity_l * temp ) &
844	+ r_v * temp / ( diff_coeff_l * e_s ) &
845	)
846	!
847	!-- Mean weight of rain drops
848	xr = hyrho(k) * qr(k,j,i) / nr(k,j,i)
849	!
850	!-- Weight averaged diameter of rain drops:
851	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
852	!
853	!-- Compute ventilation factor and intercept parameter
854	!-- (Seifert and Beheng, 2006; Seifert, 2008):
855	IF ( ventilation_effect ) THEN
856	!
857	!-- Shape parameter of gamma distribution (Milbrandt and Yau,
858	!-- 2005; Stevens and Seifert, 2008):
859	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
860	( dr - 1.4E-3_wp ) ) )
861	!
862	!-- Slope parameter of gamma distribution (Seifert, 2008):
863	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
864	( mu_r + 1.0_wp ) &
865	)**( 1.0_wp / 3.0_wp ) / dr
866
867	mu_r_2 = mu_r + 2.0_wp
868	mu_r_5d2 = mu_r + 2.5_wp
869
870	f_vent = a_vent * gamm( mu_r_2 ) * &
871	lambda_r*( -mu_r_2 ) + b_vent &
872	schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *&
873	gamm( mu_r_5d2 ) * lambda_r*( -mu_r_5d2 ) &
874	( 1.0_wp - &
875	0.5_wp * ( b_term / a_term ) * &
876	( lambda_r / ( c_term + lambda_r ) &
877	)**mu_r_5d2 - &
878	0.125_wp * ( b_term / a_term )*2 &
879	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
880	)**mu_r_5d2 - &
881	0.0625_wp * ( b_term / a_term )*3 &
882	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
883	)**mu_r_5d2 - &
884	0.0390625_wp * ( b_term / a_term )*4 &
885	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
886	)**mu_r_5d2 &
887	)
888
889	nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / &
890	gamm( mu_r + 1.0_wp )
891	ELSE
892	f_vent = 1.0_wp
893	nr_0 = nr(k,j,i) * dr
894	ENDIF
895	!
896	!-- Evaporation rate of rain water content (Seifert and
897	!-- Beheng, 2006):
898	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / &
899	hyrho(k)
900	evap = MAX( evap, -qr(k,j,i) / dt_micro )
901	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
902	-nr(k,j,i) / dt_micro )
903
904	qr(k,j,i) = qr(k,j,i) + evap * dt_micro
905	nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
906
907	ENDIF
908	ENDIF
909
910	ENDDO
911	ENDDO
912	ENDDO
913
914	CALL cpu_log( log_point_s(58), 'evaporation', 'stop' )
915
916	END SUBROUTINE evaporation_rain
917
918
919	!------------------------------------------------------------------------------!
920	! Description:
921	! ------------
922	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
923	!------------------------------------------------------------------------------!
924	SUBROUTINE sedimentation_cloud
925
926	USE arrays_3d, &
927	ONLY: ddzu, dzu, pt, prr, q, qc
928
929	USE cloud_parameters, &
930	ONLY: hyrho, l_d_cp, pt_d_t
931
932	USE control_parameters, &
933	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
934
935	USE cpulog, &
936	ONLY: cpu_log, log_point_s
937
938	USE indices, &
939	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
940
941	USE kinds
942
943	USE statistics, &
944	ONLY: weight_substep
945
946
947	IMPLICIT NONE
948
949	INTEGER(iwp) :: i !<
950	INTEGER(iwp) :: j !<
951	INTEGER(iwp) :: k !<
952
953	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
954
955	CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
956
957	sed_qc(nzt+1) = 0.0_wp
958
959	DO i = nxl, nxr
960	DO j = nys, nyn
961	DO k = nzt, nzb_s_inner(j,i)+1, -1
962
963	IF ( qc(k,j,i) > eps_sb ) THEN
964	sed_qc(k) = sed_qc_const * nc_const*( -2.0_wp / 3.0_wp ) &
965	( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
966	ELSE
967	sed_qc(k) = 0.0_wp
968	ENDIF
969
970	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &
971	dt_micro + sed_qc(k+1) &
972	)
973
974	q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
975	ddzu(k+1) / hyrho(k) * dt_micro
976	qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &
977	ddzu(k+1) / hyrho(k) * dt_micro
978	pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &
979	ddzu(k+1) / hyrho(k) * l_d_cp * &
980	pt_d_t(k) * dt_micro
981
982	!
983	!-- Compute the precipitation rate due to cloud (fog) droplets
984	IF ( call_microphysics_at_all_substeps ) THEN
985	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &
986	* weight_substep(intermediate_timestep_count)
987	ELSE
988	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
989	ENDIF
990
991	ENDDO
992	ENDDO
993	ENDDO
994
995	CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
996
997	END SUBROUTINE sedimentation_cloud
998
999
1000	!------------------------------------------------------------------------------!
1001	! Description:
1002	! ------------
1003	!> Computation of sedimentation flux. Implementation according to Stevens
1004	!> and Seifert (2008). Code is based on UCLA-LES.
1005	!------------------------------------------------------------------------------!
1006	SUBROUTINE sedimentation_rain
1007
1008	USE arrays_3d, &
1009	ONLY: ddzu, dzu, nr, pt, prr, q, qr
1010
1011	USE cloud_parameters, &
1012	ONLY: hyrho, l_d_cp, pt_d_t
1013
1014	USE control_parameters, &
1015	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
1016	USE cpulog, &
1017	ONLY: cpu_log, log_point_s
1018
1019	USE indices, &
1020	ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
1021
1022	USE kinds
1023
1024	USE statistics, &
1025	ONLY: weight_substep
1026
1027	IMPLICIT NONE
1028
1029	INTEGER(iwp) :: i !<
1030	INTEGER(iwp) :: j !<
1031	INTEGER(iwp) :: k !<
1032	INTEGER(iwp) :: k_run !<
1033
1034	REAL(wp) :: c_run !<
1035	REAL(wp) :: d_max !<
1036	REAL(wp) :: d_mean !<
1037	REAL(wp) :: d_min !<
1038	REAL(wp) :: dr !<
1039	REAL(wp) :: flux !<
1040	REAL(wp) :: lambda_r !<
1041	REAL(wp) :: mu_r !<
1042	REAL(wp) :: z_run !<
1043
1044	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
1045	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
1046	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
1047	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
1048	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
1049	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
1050	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
1051	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
1052
1053	CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
1054
1055	!
1056	!-- Compute velocities
1057	DO i = nxl, nxr
1058	DO j = nys, nyn
1059	DO k = nzb_s_inner(j,i)+1, nzt
1060	IF ( qr(k,j,i) > eps_sb ) THEN
1061	!
1062	!-- Weight averaged diameter of rain drops:
1063	dr = ( hyrho(k) * qr(k,j,i) / &
1064	nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
1065	!
1066	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1067	!-- Stevens and Seifert, 2008):
1068	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &
1069	( dr - 1.4E-3_wp ) ) )
1070	!
1071	!-- Slope parameter of gamma distribution (Seifert, 2008):
1072	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
1073	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1074
1075	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
1076	a_term - b_term * ( 1.0_wp + &
1077	c_term / &
1078	lambda_r )*( -1.0_wp &
1079	( mu_r + 1.0_wp ) ) &
1080	) &
1081	)
1082
1083	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
1084	a_term - b_term * ( 1.0_wp + &
1085	c_term / &
1086	lambda_r )*( -1.0_wp &
1087	( mu_r + 4.0_wp ) ) &
1088	) &
1089	)
1090	ELSE
1091	w_nr(k) = 0.0_wp
1092	w_qr(k) = 0.0_wp
1093	ENDIF
1094	ENDDO
1095	!
1096	!-- Adjust boundary values
1097	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
1098	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
1099	w_nr(nzt+1) = 0.0_wp
1100	w_qr(nzt+1) = 0.0_wp
1101	!
1102	!-- Compute Courant number
1103	DO k = nzb_s_inner(j,i)+1, nzt
1104	c_nr(k) = 0.25_wp * ( w_nr(k-1) + &
1105	2.0_wp * w_nr(k) + w_nr(k+1) ) * &
1106	dt_micro * ddzu(k)
1107	c_qr(k) = 0.25_wp * ( w_qr(k-1) + &
1108	2.0_wp * w_qr(k) + w_qr(k+1) ) * &
1109	dt_micro * ddzu(k)
1110	ENDDO
1111	!
1112	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
1113	IF ( limiter_sedimentation ) THEN
1114
1115	DO k = nzb_s_inner(j,i)+1, nzt
1116	d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
1117	d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
1118	d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
1119
1120	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
1121	2.0_wp * d_max, &
1122	ABS( d_mean ) )
1123
1124	d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
1125	d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
1126	d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
1127
1128	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
1129	2.0_wp * d_max, &
1130	ABS( d_mean ) )
1131	ENDDO
1132
1133	ELSE
1134
1135	nr_slope = 0.0_wp
1136	qr_slope = 0.0_wp
1137
1138	ENDIF
1139
1140	sed_nr(nzt+1) = 0.0_wp
1141	sed_qr(nzt+1) = 0.0_wp
1142	!
1143	!-- Compute sedimentation flux
1144	DO k = nzt, nzb_s_inner(j,i)+1, -1
1145	!
1146	!-- Sum up all rain drop number densities which contribute to the flux
1147	!-- through k-1/2
1148	flux = 0.0_wp
1149	z_run = 0.0_wp ! height above z(k)
1150	k_run = k
1151	c_run = MIN( 1.0_wp, c_nr(k) )
1152	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
1153	flux = flux + hyrho(k_run) * &
1154	( nr(k_run,j,i) + nr_slope(k_run) * &
1155	( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
1156	z_run = z_run + dzu(k_run)
1157	k_run = k_run + 1
1158	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
1159	ENDDO
1160	!
1161	!-- It is not allowed to sediment more rain drop number density than
1162	!-- available
1163	flux = MIN( flux, &
1164	hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &
1165	dt_micro &
1166	)
1167
1168	sed_nr(k) = flux / dt_micro
1169	nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &
1170	ddzu(k+1) / hyrho(k) * dt_micro
1171	!
1172	!-- Sum up all rain water content which contributes to the flux
1173	!-- through k-1/2
1174	flux = 0.0_wp
1175	z_run = 0.0_wp ! height above z(k)
1176	k_run = k
1177	c_run = MIN( 1.0_wp, c_qr(k) )
1178
1179	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
1180
1181	flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &
1182	qr_slope(k_run) * ( 1.0_wp - c_run ) * &
1183	0.5_wp ) * c_run * dzu(k_run)
1184	z_run = z_run + dzu(k_run)
1185	k_run = k_run + 1
1186	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
1187
1188	ENDDO
1189	!
1190	!-- It is not allowed to sediment more rain water content than
1191	!-- available
1192	flux = MIN( flux, &
1193	hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &
1194	dt_micro &
1195	)
1196
1197	sed_qr(k) = flux / dt_micro
1198
1199	qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
1200	ddzu(k+1) / hyrho(k) * dt_micro
1201	q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &
1202	ddzu(k+1) / hyrho(k) * dt_micro
1203	pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &
1204	ddzu(k+1) / hyrho(k) * l_d_cp * &
1205	pt_d_t(k) * dt_micro
1206	!
1207	!-- Compute the rain rate
1208	IF ( call_microphysics_at_all_substeps ) THEN
1209	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
1210	* weight_substep(intermediate_timestep_count)
1211	ELSE
1212	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
1213	ENDIF
1214
1215	ENDDO
1216	ENDDO
1217	ENDDO
1218
1219	CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
1220
1221	END SUBROUTINE sedimentation_rain
1222
1223
1224	!------------------------------------------------------------------------------!
1225	! Description:
1226	! ------------
1227	!> Computation of the precipitation amount due to gravitational settling of
1228	!> rain and cloud (fog) droplets
1229	!------------------------------------------------------------------------------!
1230	SUBROUTINE calc_precipitation_amount
1231
1232	USE arrays_3d, &
1233	ONLY: precipitation_amount, prr
1234
1235	USE cloud_parameters, &
1236	ONLY: hyrho
1237
1238	USE control_parameters, &
1239	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
1240	intermediate_timestep_count, intermediate_timestep_count_max,&
1241	precipitation_amount_interval, time_do2d_xy
1242
1243	USE indices, &
1244	ONLY: nxl, nxr, nys, nyn, nzb_s_inner
1245
1246	USE kinds
1247
1248	IMPLICIT NONE
1249
1250	INTEGER(iwp) :: i !:
1251	INTEGER(iwp) :: j !:
1252
1253
1254	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
1255	( .NOT. call_microphysics_at_all_substeps .OR. &
1256	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
1257	THEN
1258
1259	DO i = nxl, nxr
1260	DO j = nys, nyn
1261
1262	precipitation_amount(j,i) = precipitation_amount(j,i) + &
1263	prr(nzb_s_inner(j,i)+1,j,i) * &
1264	hyrho(nzb_s_inner(j,i)+1) * dt_3d
1265
1266	ENDDO
1267	ENDDO
1268	ENDIF
1269
1270	END SUBROUTINE calc_precipitation_amount
1271
1272
1273	!------------------------------------------------------------------------------!
1274	! Description:
1275	! ------------
1276	!> Control of microphysics for grid points i,j
1277	!------------------------------------------------------------------------------!
1278
1279	SUBROUTINE microphysics_control_ij( i, j )
1280
1281	USE arrays_3d, &
1282	ONLY: hyp, nr, pt, pt_init, prr, q, qc, qr, zu
1283
1284	USE cloud_parameters, &
1285	ONLY: cp, hyrho, pt_d_t, r_d, t_d_pt
1286
1287	USE control_parameters, &
1288	ONLY: call_microphysics_at_all_substeps, dt_3d, g, &
1289	intermediate_timestep_count, large_scale_forcing, &
1290	lsf_surf, microphysics_seifert, microphysics_kessler, &
1291	pt_surface, rho_surface, surface_pressure
1292
1293	USE indices, &
1294	ONLY: nzb, nzt
1295
1296	USE kinds
1297
1298	USE statistics, &
1299	ONLY: weight_pres
1300
1301	IMPLICIT NONE
1302
1303	INTEGER(iwp) :: i !<
1304	INTEGER(iwp) :: j !<
1305	INTEGER(iwp) :: k !<
1306
1307	REAL(wp) :: t_surface !<
1308
1309	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1310	!
1311	!-- Calculate:
1312	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
1313	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
1314	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
1315	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1316	DO k = nzb, nzt+1
1317	hyp(k) = surface_pressure * 100.0_wp * &
1318	( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp)
1319	pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
1320	t_d_pt(k) = 1.0_wp / pt_d_t(k)
1321	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
1322	ENDDO
1323	!
1324	!-- Compute reference density
1325	rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface )
1326	ENDIF
1327
1328	!
1329	!-- Compute length of time step
1330	IF ( call_microphysics_at_all_substeps ) THEN
1331	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
1332	ELSE
1333	dt_micro = dt_3d
1334	ENDIF
1335
1336	!
1337	!-- Use 1d arrays
1338	q_1d(:) = q(:,j,i)
1339	pt_1d(:) = pt(:,j,i)
1340	qc_1d(:) = qc(:,j,i)
1341	nc_1d(:) = nc_const
1342	IF ( microphysics_seifert ) THEN
1343	qr_1d(:) = qr(:,j,i)
1344	nr_1d(:) = nr(:,j,i)
1345	ENDIF
1346
1347	!
1348	!-- Reset precipitation rate
1349	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp
1350
1351	!
1352	!-- Compute cloud physics
1353	IF( microphysics_kessler ) THEN
1354
1355	CALL autoconversion_kessler( i,j )
1356	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
1357
1358	ELSEIF ( microphysics_seifert ) THEN
1359
1360	CALL adjust_cloud( i,j )
1361	CALL autoconversion( i,j )
1362	CALL accretion( i,j )
1363	CALL selfcollection_breakup( i,j )
1364	CALL evaporation_rain( i,j )
1365	CALL sedimentation_rain( i,j )
1366	IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )
1367
1368	ENDIF
1369
1370	CALL calc_precipitation_amount( i,j )
1371
1372	!
1373	!-- Store results on the 3d arrays
1374	q(:,j,i) = q_1d(:)
1375	pt(:,j,i) = pt_1d(:)
1376	IF ( microphysics_seifert ) THEN
1377	qr(:,j,i) = qr_1d(:)
1378	nr(:,j,i) = nr_1d(:)
1379	ENDIF
1380
1381	END SUBROUTINE microphysics_control_ij
1382
1383	!------------------------------------------------------------------------------!
1384	! Description:
1385	! ------------
1386	!> Adjust number of raindrops to avoid nonlinear effects in
1387	!> sedimentation and evaporation of rain drops due to too small or
1388	!> too big weights of rain drops (Stevens and Seifert, 2008).
1389	!> The same procedure is applied to cloud droplets if they are determined
1390	!> prognostically. Call for grid point i,j
1391	!------------------------------------------------------------------------------!
1392	SUBROUTINE adjust_cloud_ij( i, j )
1393
1394	USE cloud_parameters, &
1395	ONLY: hyrho
1396
1397	USE indices, &
1398	ONLY: nzb_s_inner, nzt
1399
1400	USE kinds
1401
1402	IMPLICIT NONE
1403
1404	INTEGER(iwp) :: i !<
1405	INTEGER(iwp) :: j !<
1406	INTEGER(iwp) :: k !<
1407
1408	DO k = nzb_s_inner(j,i)+1, nzt
1409
1410	IF ( qr_1d(k) <= eps_sb ) THEN
1411	qr_1d(k) = 0.0_wp
1412	nr_1d(k) = 0.0_wp
1413	ELSE
1414	!
1415	!-- Adjust number of raindrops to avoid nonlinear effects in
1416	!-- sedimentation and evaporation of rain drops due to too small or
1417	!-- too big weights of rain drops (Stevens and Seifert, 2008).
1418	IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN
1419	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
1420	ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN
1421	nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
1422	ENDIF
1423
1424	ENDIF
1425
1426	ENDDO
1427
1428	END SUBROUTINE adjust_cloud_ij
1429
1430
1431	!------------------------------------------------------------------------------!
1432	! Description:
1433	! ------------
1434	!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
1435	!------------------------------------------------------------------------------!
1436	SUBROUTINE autoconversion_ij( i, j )
1437
1438	USE arrays_3d, &
1439	ONLY: diss, dzu
1440
1441	USE cloud_parameters, &
1442	ONLY: hyrho
1443
1444	USE control_parameters, &
1445	ONLY: rho_surface
1446
1447	USE grid_variables, &
1448	ONLY: dx, dy
1449
1450	USE indices, &
1451	ONLY: nzb_s_inner, nzt
1452
1453	USE kinds
1454
1455	IMPLICIT NONE
1456
1457	INTEGER(iwp) :: i !<
1458	INTEGER(iwp) :: j !<
1459	INTEGER(iwp) :: k !<
1460
1461	REAL(wp) :: alpha_cc !<
1462	REAL(wp) :: autocon !<
1463	REAL(wp) :: dissipation !<
1464	REAL(wp) :: k_au !<
1465	REAL(wp) :: l_mix !<
1466	REAL(wp) :: nu_c !<
1467	REAL(wp) :: phi_au !<
1468	REAL(wp) :: r_cc !<
1469	REAL(wp) :: rc !<
1470	REAL(wp) :: re_lambda !<
1471	REAL(wp) :: sigma_cc !<
1472	REAL(wp) :: tau_cloud !<
1473	REAL(wp) :: xc !<
1474
1475	DO k = nzb_s_inner(j,i)+1, nzt
1476
1477	IF ( qc_1d(k) > eps_sb ) THEN
1478
1479	k_au = k_cc / ( 20.0_wp * x0 )
1480	!
1481	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
1482	!-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
1483	tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
1484	!
1485	!-- Universal function for autoconversion process
1486	!-- (Seifert and Beheng, 2006):
1487	phi_au = 600.0_wp * tau_cloud*0.68_wp ( 1.0_wp - tau_cloud0.68_wp )3
1488	!
1489	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
1490	!-- (Use constant nu_c = 1.0_wp instead?)
1491	nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp )
1492	!
1493	!-- Mean weight of cloud droplets:
1494	xc = hyrho(k) * qc_1d(k) / nc_1d(k)
1495	!
1496	!-- Parameterized turbulence effects on autoconversion (Seifert,
1497	!-- Nuijens and Stevens, 2010)
1498	IF ( collision_turbulence ) THEN
1499	!
1500	!-- Weight averaged radius of cloud droplets:
1501	rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
1502
1503	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
1504	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
1505	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
1506	!
1507	!-- Mixing length (neglecting distance to ground and stratification)
1508	l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
1509	!
1510	!-- Limit dissipation rate according to Seifert, Nuijens and
1511	!-- Stevens (2010)
1512	dissipation = MIN( 0.06_wp, diss(k,j,i) )
1513	!
1514	!-- Compute Taylor-microscale Reynolds number:
1515	re_lambda = 6.0_wp / 11.0_wp * &
1516	( l_mix / c_const )*( 2.0_wp / 3.0_wp ) &
1517	SQRT( 15.0_wp / kin_vis_air ) * &
1518	dissipation**( 1.0_wp / 6.0_wp )
1519	!
1520	!-- The factor of 1.0E4 is needed to convert the dissipation rate
1521	!-- from m2 s-3 to cm2 s-3.
1522	k_au = k_au * ( 1.0_wp + &
1523	dissipation * 1.0E4_wp * &
1524	( re_lambda * 1.0E-3_wp )*0.25_wp &
1525	( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &
1526	sigma_cc )**2 &
1527	) + beta_cc &
1528	) &
1529	)
1530	ENDIF
1531	!
1532	!-- Autoconversion rate (Seifert and Beheng, 2006):
1533	autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &
1534	( nu_c + 1.0_wp )*2 qc_1d(k)*2 xc*2 &
1535	( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )*2 ) &
1536	rho_surface
1537	autocon = MIN( autocon, qc_1d(k) / dt_micro )
1538
1539	qr_1d(k) = qr_1d(k) + autocon * dt_micro
1540	qc_1d(k) = qc_1d(k) - autocon * dt_micro
1541	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
1542
1543	ENDIF
1544
1545	ENDDO
1546
1547	END SUBROUTINE autoconversion_ij
1548
1549	!------------------------------------------------------------------------------!
1550	! Description:
1551	! ------------
1552	!> Autoconversion process (Kessler, 1969).
1553	!------------------------------------------------------------------------------!
1554	SUBROUTINE autoconversion_kessler_ij( i, j )
1555
1556	USE arrays_3d, &
1557	ONLY: dzw, prr
1558
1559	USE cloud_parameters, &
1560	ONLY: l_d_cp, pt_d_t
1561
1562	USE indices, &
1563	ONLY: nzb_s_inner, nzt
1564
1565	USE kinds
1566
1567
1568	IMPLICIT NONE
1569
1570	INTEGER(iwp) :: i !<
1571	INTEGER(iwp) :: j !<
1572	INTEGER(iwp) :: k !<
1573
1574	REAL(wp) :: dqdt_precip !<
1575
1576	DO k = nzb_s_inner(j,i)+1, nzt
1577
1578	IF ( qc_1d(k) > ql_crit ) THEN
1579	dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit )
1580	ELSE
1581	dqdt_precip = 0.0_wp
1582	ENDIF
1583
1584	qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro
1585	q_1d(k) = q_1d(k) - dqdt_precip * dt_micro
1586	pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
1587
1588	!
1589	!-- Compute the rain rate (stored on surface grid point)
1590	prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) + &
1591	dqdt_precip * dzw(k)
1592
1593	ENDDO
1594
1595	END SUBROUTINE autoconversion_kessler_ij
1596
1597	!------------------------------------------------------------------------------!
1598	! Description:
1599	! ------------
1600	!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
1601	!------------------------------------------------------------------------------!
1602	SUBROUTINE accretion_ij( i, j )
1603
1604	USE arrays_3d, &
1605	ONLY: diss
1606
1607	USE cloud_parameters, &
1608	ONLY: hyrho
1609
1610	USE control_parameters, &
1611	ONLY: rho_surface
1612
1613	USE indices, &
1614	ONLY: nzb_s_inner, nzt
1615
1616	USE kinds
1617
1618	IMPLICIT NONE
1619
1620	INTEGER(iwp) :: i !<
1621	INTEGER(iwp) :: j !<
1622	INTEGER(iwp) :: k !<
1623
1624	REAL(wp) :: accr !<
1625	REAL(wp) :: k_cr !<
1626	REAL(wp) :: phi_ac !<
1627	REAL(wp) :: tau_cloud !<
1628
1629	DO k = nzb_s_inner(j,i)+1, nzt
1630	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN
1631	!
1632	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
1633	tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
1634	!
1635	!-- Universal function for accretion process
1636	!-- (Seifert and Beheng, 2001):
1637	phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
1638	!
1639	!-- Parameterized turbulence effects on autoconversion (Seifert,
1640	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
1641	!-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
1642	IF ( collision_turbulence ) THEN
1643	k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &
1644	MIN( 600.0_wp, &
1645	diss(k,j,i) * 1.0E4_wp )**0.25_wp &
1646	)
1647	ELSE
1648	k_cr = k_cr0
1649	ENDIF
1650	!
1651	!-- Accretion rate (Seifert and Beheng, 2006):
1652	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
1653	accr = MIN( accr, qc_1d(k) / dt_micro )
1654
1655	qr_1d(k) = qr_1d(k) + accr * dt_micro
1656	qc_1d(k) = qc_1d(k) - accr * dt_micro
1657
1658	ENDIF
1659
1660	ENDDO
1661
1662	END SUBROUTINE accretion_ij
1663
1664
1665	!------------------------------------------------------------------------------!
1666	! Description:
1667	! ------------
1668	!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
1669	!------------------------------------------------------------------------------!
1670	SUBROUTINE selfcollection_breakup_ij( i, j )
1671
1672	USE cloud_parameters, &
1673	ONLY: hyrho
1674
1675	USE control_parameters, &
1676	ONLY: rho_surface
1677
1678	USE indices, &
1679	ONLY: nzb_s_inner, nzt
1680
1681	USE kinds
1682
1683	IMPLICIT NONE
1684
1685	INTEGER(iwp) :: i !<
1686	INTEGER(iwp) :: j !<
1687	INTEGER(iwp) :: k !<
1688
1689	REAL(wp) :: breakup !<
1690	REAL(wp) :: dr !<
1691	REAL(wp) :: phi_br !<
1692	REAL(wp) :: selfcoll !<
1693
1694	DO k = nzb_s_inner(j,i)+1, nzt
1695	IF ( qr_1d(k) > eps_sb ) THEN
1696	!
1697	!-- Selfcollection rate (Seifert and Beheng, 2001):
1698	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface )
1699	!
1700	!-- Weight averaged diameter of rain drops:
1701	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1702	!
1703	!-- Collisional breakup rate (Seifert, 2008):
1704	IF ( dr >= 0.3E-3_wp ) THEN
1705	phi_br = k_br * ( dr - 1.1E-3_wp )
1706	breakup = selfcoll * ( phi_br + 1.0_wp )
1707	ELSE
1708	breakup = 0.0_wp
1709	ENDIF
1710
1711	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
1712	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
1713
1714	ENDIF
1715	ENDDO
1716
1717	END SUBROUTINE selfcollection_breakup_ij
1718
1719
1720	!------------------------------------------------------------------------------!
1721	! Description:
1722	! ------------
1723	!> Evaporation of precipitable water. Condensation is neglected for
1724	!> precipitable water. Call for grid point i,j
1725	!------------------------------------------------------------------------------!
1726	SUBROUTINE evaporation_rain_ij( i, j )
1727
1728	USE arrays_3d, &
1729	ONLY: hyp
1730
1731	USE cloud_parameters, &
1732	ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
1733
1734	USE constants, &
1735	ONLY: pi
1736
1737	USE indices, &
1738	ONLY: nzb_s_inner, nzt
1739
1740	USE kinds
1741
1742	IMPLICIT NONE
1743
1744	INTEGER(iwp) :: i !<
1745	INTEGER(iwp) :: j !<
1746	INTEGER(iwp) :: k !<
1747
1748	REAL(wp) :: alpha !<
1749	REAL(wp) :: dr !<
1750	REAL(wp) :: e_s !<
1751	REAL(wp) :: evap !<
1752	REAL(wp) :: evap_nr !<
1753	REAL(wp) :: f_vent !<
1754	REAL(wp) :: g_evap !<
1755	REAL(wp) :: lambda_r !<
1756	REAL(wp) :: mu_r !<
1757	REAL(wp) :: mu_r_2 !<
1758	REAL(wp) :: mu_r_5d2 !<
1759	REAL(wp) :: nr_0 !<
1760	REAL(wp) :: q_s !<
1761	REAL(wp) :: sat !<
1762	REAL(wp) :: t_l !<
1763	REAL(wp) :: temp !<
1764	REAL(wp) :: xr !<
1765
1766	DO k = nzb_s_inner(j,i)+1, nzt
1767	IF ( qr_1d(k) > eps_sb ) THEN
1768	!
1769	!-- Actual liquid water temperature:
1770	t_l = t_d_pt(k) * pt_1d(k)
1771	!
1772	!-- Saturation vapor pressure at t_l:
1773	e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / &
1774	( t_l - 35.86_wp ) &
1775	)
1776	!
1777	!-- Computation of saturation humidity:
1778	q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
1779	alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
1780	q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s )
1781	!
1782	!-- Supersaturation:
1783	sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
1784	!
1785	!-- Evaporation needs only to be calculated in subsaturated regions
1786	IF ( sat < 0.0_wp ) THEN
1787	!
1788	!-- Actual temperature:
1789	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
1790
1791	g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / &
1792	( thermal_conductivity_l * temp ) + &
1793	r_v * temp / ( diff_coeff_l * e_s ) &
1794	)
1795	!
1796	!-- Mean weight of rain drops
1797	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
1798	!
1799	!-- Weight averaged diameter of rain drops:
1800	dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
1801	!
1802	!-- Compute ventilation factor and intercept parameter
1803	!-- (Seifert and Beheng, 2006; Seifert, 2008):
1804	IF ( ventilation_effect ) THEN
1805	!
1806	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1807	!-- Stevens and Seifert, 2008):
1808	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1809	!
1810	!-- Slope parameter of gamma distribution (Seifert, 2008):
1811	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
1812	( mu_r + 1.0_wp ) &
1813	)**( 1.0_wp / 3.0_wp ) / dr
1814
1815	mu_r_2 = mu_r + 2.0_wp
1816	mu_r_5d2 = mu_r + 2.5_wp
1817
1818	f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + &
1819	b_vent * schmidt_p_1d3 * &
1820	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
1821	lambda_r*( -mu_r_5d2 ) &
1822	( 1.0_wp - &
1823	0.5_wp * ( b_term / a_term ) * &
1824	( lambda_r / ( c_term + lambda_r ) &
1825	)**mu_r_5d2 - &
1826	0.125_wp * ( b_term / a_term )*2 &
1827	( lambda_r / ( 2.0_wp * c_term + lambda_r ) &
1828	)**mu_r_5d2 - &
1829	0.0625_wp * ( b_term / a_term )*3 &
1830	( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
1831	)**mu_r_5d2 - &
1832	0.0390625_wp * ( b_term / a_term )*4 &
1833	( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
1834	)**mu_r_5d2 &
1835	)
1836
1837	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / &
1838	gamm( mu_r + 1.0_wp )
1839	ELSE
1840	f_vent = 1.0_wp
1841	nr_0 = nr_1d(k) * dr
1842	ENDIF
1843	!
1844	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
1845	evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k)
1846	evap = MAX( evap, -qr_1d(k) / dt_micro )
1847	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
1848	-nr_1d(k) / dt_micro )
1849
1850	qr_1d(k) = qr_1d(k) + evap * dt_micro
1851	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
1852
1853	ENDIF
1854	ENDIF
1855
1856	ENDDO
1857
1858	END SUBROUTINE evaporation_rain_ij
1859
1860
1861	!------------------------------------------------------------------------------!
1862	! Description:
1863	! ------------
1864	!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
1865	!> Call for grid point i,j
1866	!------------------------------------------------------------------------------!
1867	SUBROUTINE sedimentation_cloud_ij( i, j )
1868
1869	USE arrays_3d, &
1870	ONLY: ddzu, dzu, prr
1871
1872	USE cloud_parameters, &
1873	ONLY: hyrho, l_d_cp, pt_d_t
1874
1875	USE control_parameters, &
1876	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
1877
1878	USE indices, &
1879	ONLY: nzb, nzb_s_inner, nzt
1880
1881	USE kinds
1882
1883	USE statistics, &
1884	ONLY: weight_substep
1885
1886	IMPLICIT NONE
1887
1888	INTEGER(iwp) :: i !<
1889	INTEGER(iwp) :: j !<
1890	INTEGER(iwp) :: k !<
1891
1892	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
1893
1894	sed_qc(nzt+1) = 0.0_wp
1895
1896	DO k = nzt, nzb_s_inner(j,i)+1, -1
1897	IF ( qc_1d(k) > eps_sb ) THEN
1898	sed_qc(k) = sed_qc_const * nc_1d(k)*( -2.0_wp / 3.0_wp ) &
1899	( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
1900	ELSE
1901	sed_qc(k) = 0.0_wp
1902	ENDIF
1903
1904	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
1905	dt_micro + sed_qc(k+1) &
1906	)
1907
1908	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
1909	hyrho(k) * dt_micro
1910	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
1911	hyrho(k) * dt_micro
1912	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
1913	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
1914
1915	!
1916	!-- Compute the precipitation rate of cloud (fog) droplets
1917	IF ( call_microphysics_at_all_substeps ) THEN
1918	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * &
1919	weight_substep(intermediate_timestep_count)
1920	ELSE
1921	prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
1922	ENDIF
1923
1924	ENDDO
1925
1926	END SUBROUTINE sedimentation_cloud_ij
1927
1928
1929	!------------------------------------------------------------------------------!
1930	! Description:
1931	! ------------
1932	!> Computation of sedimentation flux. Implementation according to Stevens
1933	!> and Seifert (2008). Code is based on UCLA-LES. Call for grid point i,j
1934	!------------------------------------------------------------------------------!
1935	SUBROUTINE sedimentation_rain_ij( i, j )
1936
1937	USE arrays_3d, &
1938	ONLY: ddzu, dzu, prr
1939
1940	USE cloud_parameters, &
1941	ONLY: hyrho, l_d_cp, pt_d_t
1942
1943	USE control_parameters, &
1944	ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count
1945
1946	USE indices, &
1947	ONLY: nzb, nzb_s_inner, nzt
1948
1949	USE kinds
1950
1951	USE statistics, &
1952	ONLY: weight_substep
1953
1954	IMPLICIT NONE
1955
1956	INTEGER(iwp) :: i !<
1957	INTEGER(iwp) :: j !<
1958	INTEGER(iwp) :: k !<
1959	INTEGER(iwp) :: k_run !<
1960
1961	REAL(wp) :: c_run !<
1962	REAL(wp) :: d_max !<
1963	REAL(wp) :: d_mean !<
1964	REAL(wp) :: d_min !<
1965	REAL(wp) :: dr !<
1966	REAL(wp) :: flux !<
1967	REAL(wp) :: lambda_r !<
1968	REAL(wp) :: mu_r !<
1969	REAL(wp) :: z_run !<
1970
1971	REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<
1972	REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<
1973	REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<
1974	REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<
1975	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<
1976	REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<
1977	REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<
1978	REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<
1979
1980	!
1981	!-- Compute velocities
1982	DO k = nzb_s_inner(j,i)+1, nzt
1983	IF ( qr_1d(k) > eps_sb ) THEN
1984	!
1985	!-- Weight averaged diameter of rain drops:
1986	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
1987	!
1988	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
1989	!-- Stevens and Seifert, 2008):
1990	mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) )
1991	!
1992	!-- Slope parameter of gamma distribution (Seifert, 2008):
1993	lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &
1994	( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
1995
1996	w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
1997	a_term - b_term * ( 1.0_wp + &
1998	c_term / lambda_r )*( -1.0_wp &
1999	( mu_r + 1.0_wp ) ) &
2000	) &
2001	)
2002	w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &
2003	a_term - b_term * ( 1.0_wp + &
2004	c_term / lambda_r )*( -1.0_wp &
2005	( mu_r + 4.0_wp ) ) &
2006	) &
2007	)
2008	ELSE
2009	w_nr(k) = 0.0_wp
2010	w_qr(k) = 0.0_wp
2011	ENDIF
2012	ENDDO
2013	!
2014	!-- Adjust boundary values
2015	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
2016	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
2017	w_nr(nzt+1) = 0.0_wp
2018	w_qr(nzt+1) = 0.0_wp
2019	!
2020	!-- Compute Courant number
2021	DO k = nzb_s_inner(j,i)+1, nzt
2022	c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &
2023	dt_micro * ddzu(k)
2024	c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &
2025	dt_micro * ddzu(k)
2026	ENDDO
2027	!
2028	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
2029	IF ( limiter_sedimentation ) THEN
2030
2031	DO k = nzb_s_inner(j,i)+1, nzt
2032	d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
2033	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
2034	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
2035
2036	qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
2037	2.0_wp * d_max, &
2038	ABS( d_mean ) )
2039
2040	d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
2041	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
2042	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
2043
2044	nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &
2045	2.0_wp * d_max, &
2046	ABS( d_mean ) )
2047	ENDDO
2048
2049	ELSE
2050
2051	nr_slope = 0.0_wp
2052	qr_slope = 0.0_wp
2053
2054	ENDIF
2055
2056	sed_nr(nzt+1) = 0.0_wp
2057	sed_qr(nzt+1) = 0.0_wp
2058	!
2059	!-- Compute sedimentation flux
2060	DO k = nzt, nzb_s_inner(j,i)+1, -1
2061	!
2062	!-- Sum up all rain drop number densities which contribute to the flux
2063	!-- through k-1/2
2064	flux = 0.0_wp
2065	z_run = 0.0_wp ! height above z(k)
2066	k_run = k
2067	c_run = MIN( 1.0_wp, c_nr(k) )
2068	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
2069	flux = flux + hyrho(k_run) * &
2070	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * &
2071	0.5_wp ) * c_run * dzu(k_run)
2072	z_run = z_run + dzu(k_run)
2073	k_run = k_run + 1
2074	c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
2075	ENDDO
2076	!
2077	!-- It is not allowed to sediment more rain drop number density than
2078	!-- available
2079	flux = MIN( flux, &
2080	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
2081
2082	sed_nr(k) = flux / dt_micro
2083	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
2084	hyrho(k) * dt_micro
2085	!
2086	!-- Sum up all rain water content which contributes to the flux
2087	!-- through k-1/2
2088	flux = 0.0_wp
2089	z_run = 0.0_wp ! height above z(k)
2090	k_run = k
2091	c_run = MIN( 1.0_wp, c_qr(k) )
2092
2093	DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )
2094
2095	flux = flux + hyrho(k_run) * &
2096	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * &
2097	0.5_wp ) * c_run * dzu(k_run)
2098	z_run = z_run + dzu(k_run)
2099	k_run = k_run + 1
2100	c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
2101
2102	ENDDO
2103	!
2104	!-- It is not allowed to sediment more rain water content than available
2105	flux = MIN( flux, &
2106	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
2107
2108	sed_qr(k) = flux / dt_micro
2109
2110	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
2111	hyrho(k) * dt_micro
2112	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
2113	hyrho(k) * dt_micro
2114	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
2115	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
2116	!
2117	!-- Compute the rain rate
2118	IF ( call_microphysics_at_all_substeps ) THEN
2119	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &
2120	* weight_substep(intermediate_timestep_count)
2121	ELSE
2122	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
2123	ENDIF
2124
2125	ENDDO
2126
2127	END SUBROUTINE sedimentation_rain_ij
2128
2129
2130	!------------------------------------------------------------------------------!
2131	! Description:
2132	! ------------
2133	!> This subroutine computes the precipitation amount due to gravitational
2134	!> settling of rain and cloud (fog) droplets
2135	!------------------------------------------------------------------------------!
2136	SUBROUTINE calc_precipitation_amount_ij( i, j )
2137
2138	USE arrays_3d, &
2139	ONLY: precipitation_amount, prr
2140
2141	USE cloud_parameters, &
2142	ONLY: hyrho
2143
2144	USE control_parameters, &
2145	ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, &
2146	intermediate_timestep_count, intermediate_timestep_count_max,&
2147	precipitation_amount_interval, time_do2d_xy
2148
2149	USE indices, &
2150	ONLY: nzb_s_inner
2151
2152	USE kinds
2153
2154	IMPLICIT NONE
2155
2156	INTEGER(iwp) :: i !:
2157	INTEGER(iwp) :: j !:
2158
2159
2160	IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
2161	( .NOT. call_microphysics_at_all_substeps .OR. &
2162	intermediate_timestep_count == intermediate_timestep_count_max ) ) &
2163	THEN
2164
2165	precipitation_amount(j,i) = precipitation_amount(j,i) + &
2166	prr(nzb_s_inner(j,i)+1,j,i) * &
2167	hyrho(nzb_s_inner(j,i)+1) * dt_3d
2168	ENDIF
2169
2170	END SUBROUTINE calc_precipitation_amount_ij
2171
2172	!------------------------------------------------------------------------------!
2173	! Description:
2174	! ------------
2175	!> This function computes the gamma function (Press et al., 1992).
2176	!> The gamma function is needed for the calculation of the evaporation
2177	!> of rain drops.
2178	!------------------------------------------------------------------------------!
2179	FUNCTION gamm( xx )
2180
2181	USE kinds
2182
2183	IMPLICIT NONE
2184
2185	INTEGER(iwp) :: j !<
2186
2187	REAL(wp) :: gamm !<
2188	REAL(wp) :: ser !<
2189	REAL(wp) :: tmp !<
2190	REAL(wp) :: x_gamm !<
2191	REAL(wp) :: xx !<
2192	REAL(wp) :: y_gamm !<
2193
2194
2195	REAL(wp), PARAMETER :: stp = 2.5066282746310005_wp !<
2196	REAL(wp), PARAMETER :: cof(6) = (/ 76.18009172947146_wp, &
2197	-86.50532032941677_wp, &
2198	24.01409824083091_wp, &
2199	-1.231739572450155_wp, &
2200	0.1208650973866179E-2_wp, &
2201	-0.5395239384953E-5_wp /) !<
2202
2203	x_gamm = xx
2204	y_gamm = x_gamm
2205	tmp = x_gamm + 5.5_wp
2206	tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp
2207	ser = 1.000000000190015_wp
2208
2209	DO j = 1, 6
2210	y_gamm = y_gamm + 1.0_wp
2211	ser = ser + cof( j ) / y_gamm
2212	ENDDO
2213
2214	!
2215	!-- Until this point the algorithm computes the logarithm of the gamma
2216	!-- function. Hence, the exponential function is used.
2217	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
2218	gamm = EXP( tmp ) * stp * ser / x_gamm
2219
2220	RETURN
2221
2222	END FUNCTION gamm
2223
2224	END MODULE microphysics_mod

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

Context Navigation

source: palm/trunk/SOURCE/microphysics.f90 @ 1849

Download in other formats: