Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |