Home

Context Navigation

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

Last change on this file since 1263 was 1242, checked in by heinze, 11 years ago
Last commmit documented
Property svn:keywords set to `Id`
File size: 28.5 KB

Line
1	MODULE microphysics_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: microphysics.f90 1242 2013-10-30 11:50:11Z fricke $
27	!
28	! 1241 2013-10-30 11:36:58Z heinze
29	! hyp and rho have to be calculated at each time step if data from external
30	! file LSF_DATA are used
31	!
32	! 1115 2013-03-26 18:16:16Z hoffmann
33	! microphyical tendencies are calculated in microphysics_control in an optimized
34	! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
35	!
36	! 1106 2013-03-04 05:31:38Z raasch
37	! small changes in code formatting
38	!
39	! 1092 2013-02-02 11:24:22Z raasch
40	! unused variables removed
41	! file put under GPL
42	!
43	! 1065 2012-11-22 17:42:36Z hoffmann
44	! Sedimentation process implemented according to Stevens and Seifert (2008).
45	! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
46	! and Stevens, 2010).
47	!
48	! 1053 2012-11-13 17:11:03Z hoffmann
49	! initial revision
50	!
51	! Description:
52	! ------------
53	! Calculate cloud microphysics according to the two moment bulk
54	! scheme by Seifert and Beheng (2006).
55	!------------------------------------------------------------------------------!
56
57	PRIVATE
58	PUBLIC microphysics_control
59
60	INTERFACE microphysics_control
61	MODULE PROCEDURE microphysics_control
62	MODULE PROCEDURE microphysics_control_ij
63	END INTERFACE microphysics_control
64
65	INTERFACE adjust_cloud
66	MODULE PROCEDURE adjust_cloud
67	MODULE PROCEDURE adjust_cloud_ij
68	END INTERFACE adjust_cloud
69
70	INTERFACE autoconversion
71	MODULE PROCEDURE autoconversion
72	MODULE PROCEDURE autoconversion_ij
73	END INTERFACE autoconversion
74
75	INTERFACE accretion
76	MODULE PROCEDURE accretion
77	MODULE PROCEDURE accretion_ij
78	END INTERFACE accretion
79
80	INTERFACE selfcollection_breakup
81	MODULE PROCEDURE selfcollection_breakup
82	MODULE PROCEDURE selfcollection_breakup_ij
83	END INTERFACE selfcollection_breakup
84
85	INTERFACE evaporation_rain
86	MODULE PROCEDURE evaporation_rain
87	MODULE PROCEDURE evaporation_rain_ij
88	END INTERFACE evaporation_rain
89
90	INTERFACE sedimentation_cloud
91	MODULE PROCEDURE sedimentation_cloud
92	MODULE PROCEDURE sedimentation_cloud_ij
93	END INTERFACE sedimentation_cloud
94
95	INTERFACE sedimentation_rain
96	MODULE PROCEDURE sedimentation_rain
97	MODULE PROCEDURE sedimentation_rain_ij
98	END INTERFACE sedimentation_rain
99
100	CONTAINS
101
102
103	!------------------------------------------------------------------------------!
104	! Call for all grid points
105	!------------------------------------------------------------------------------!
106	SUBROUTINE microphysics_control
107
108	USE arrays_3d
109	USE cloud_parameters
110	USE control_parameters
111	USE grid_variables
112	USE indices
113	USE statistics
114
115	IMPLICIT NONE
116
117	INTEGER :: i, j, k
118
119
120	DO i = nxl, nxr
121	DO j = nys, nyn
122	DO k = nzb_s_inner(j,i)+1, nzt
123
124	ENDDO
125	ENDDO
126	ENDDO
127
128	END SUBROUTINE microphysics_control
129
130	SUBROUTINE adjust_cloud
131
132	USE arrays_3d
133	USE cloud_parameters
134	USE indices
135
136	IMPLICIT NONE
137
138	INTEGER :: i, j, k
139
140
141	DO i = nxl, nxr
142	DO j = nys, nyn
143	DO k = nzb_s_inner(j,i)+1, nzt
144
145	ENDDO
146	ENDDO
147	ENDDO
148
149	END SUBROUTINE adjust_cloud
150
151
152	SUBROUTINE autoconversion
153
154	USE arrays_3d
155	USE cloud_parameters
156	USE control_parameters
157	USE grid_variables
158	USE indices
159
160	IMPLICIT NONE
161
162	INTEGER :: i, j, k
163
164
165	DO i = nxl, nxr
166	DO j = nys, nyn
167	DO k = nzb_s_inner(j,i)+1, nzt
168
169	ENDDO
170	ENDDO
171	ENDDO
172
173	END SUBROUTINE autoconversion
174
175
176	SUBROUTINE accretion
177
178	USE arrays_3d
179	USE cloud_parameters
180	USE control_parameters
181	USE indices
182
183	IMPLICIT NONE
184
185	INTEGER :: i, j, k
186
187
188	DO i = nxl, nxr
189	DO j = nys, nyn
190	DO k = nzb_s_inner(j,i)+1, nzt
191
192	ENDDO
193	ENDDO
194	ENDDO
195
196	END SUBROUTINE accretion
197
198
199	SUBROUTINE selfcollection_breakup
200
201	USE arrays_3d
202	USE cloud_parameters
203	USE control_parameters
204	USE indices
205
206	IMPLICIT NONE
207
208	INTEGER :: i, j, k
209
210
211	DO i = nxl, nxr
212	DO j = nys, nyn
213	DO k = nzb_s_inner(j,i)+1, nzt
214
215	ENDDO
216	ENDDO
217	ENDDO
218
219	END SUBROUTINE selfcollection_breakup
220
221
222	SUBROUTINE evaporation_rain
223
224	USE arrays_3d
225	USE cloud_parameters
226	USE constants
227	USE control_parameters
228	USE indices
229
230	IMPLICIT NONE
231
232	INTEGER :: i, j, k
233
234
235	DO i = nxl, nxr
236	DO j = nys, nyn
237	DO k = nzb_s_inner(j,i)+1, nzt
238
239	ENDDO
240	ENDDO
241	ENDDO
242
243	END SUBROUTINE evaporation_rain
244
245
246	SUBROUTINE sedimentation_cloud
247
248	USE arrays_3d
249	USE cloud_parameters
250	USE constants
251	USE control_parameters
252	USE indices
253
254	IMPLICIT NONE
255
256	INTEGER :: i, j, k
257
258
259	DO i = nxl, nxr
260	DO j = nys, nyn
261	DO k = nzb_s_inner(j,i)+1, nzt
262
263	ENDDO
264	ENDDO
265	ENDDO
266
267	END SUBROUTINE sedimentation_cloud
268
269
270	SUBROUTINE sedimentation_rain
271
272	USE arrays_3d
273	USE cloud_parameters
274	USE constants
275	USE control_parameters
276	USE indices
277	USE statistics
278
279	IMPLICIT NONE
280
281	INTEGER :: i, j, k
282
283
284	DO i = nxl, nxr
285	DO j = nys, nyn
286	DO k = nzb_s_inner(j,i)+1, nzt
287
288	ENDDO
289	ENDDO
290	ENDDO
291
292	END SUBROUTINE sedimentation_rain
293
294
295	!------------------------------------------------------------------------------!
296	! Call for grid point i,j
297	!------------------------------------------------------------------------------!
298
299	SUBROUTINE microphysics_control_ij( i, j )
300
301	USE arrays_3d
302	USE cloud_parameters
303	USE control_parameters
304	USE grid_variables
305	USE indices
306	USE statistics
307
308	IMPLICIT NONE
309
310	INTEGER :: i, j, k
311	REAL :: t_surface
312
313	IF ( large_scale_forcing .AND. lsf_surf ) THEN
314	!
315	!-- Calculate:
316	!-- pt / t : ratio of potential and actual temperature (pt_d_t)
317	!-- t / pt : ratio of actual and potential temperature (t_d_pt)
318	!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
319	t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286
320	DO k = nzb, nzt+1
321	hyp(k) = surface_pressure * 100.0 * &
322	( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
323	pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
324	t_d_pt(k) = 1.0 / pt_d_t(k)
325	hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
326	ENDDO
327	!
328	!-- Compute reference density
329	rho_surface = surface_pressure * 100.0 / ( r_d * t_surface )
330	ENDIF
331
332
333	dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
334	!
335	!-- Adjust unrealistic values
336	IF ( precipitation ) CALL adjust_cloud( i,j )
337	!
338	!-- Use 1-d arrays
339	q_1d(:) = q(:,j,i)
340	pt_1d(:) = pt(:,j,i)
341	qc_1d(:) = qc(:,j,i)
342	nc_1d(:) = nc_const
343	IF ( precipitation ) THEN
344	qr_1d(:) = qr(:,j,i)
345	nr_1d(:) = nr(:,j,i)
346	ENDIF
347	!
348	!-- Compute cloud physics
349	IF ( precipitation ) THEN
350	CALL autoconversion( i,j )
351	CALL accretion( i,j )
352	CALL selfcollection_breakup( i,j )
353	CALL evaporation_rain( i,j )
354	CALL sedimentation_rain( i,j )
355	ENDIF
356
357	IF ( drizzle ) CALL sedimentation_cloud( i,j )
358	!
359	!-- Derive tendencies
360	tend_q(:,j,i) = ( q_1d(:) - q(:,j,i) ) / dt_micro
361	tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro
362	IF ( precipitation ) THEN
363	tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro
364	tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro
365	ENDIF
366
367	END SUBROUTINE microphysics_control_ij
368
369	SUBROUTINE adjust_cloud_ij( i, j )
370
371	USE arrays_3d
372	USE cloud_parameters
373	USE indices
374
375	IMPLICIT NONE
376
377	INTEGER :: i, j, k
378	!
379	!-- Adjust number of raindrops to avoid nonlinear effects in
380	!-- sedimentation and evaporation of rain drops due to too small or
381	!-- too big weights of rain drops (Stevens and Seifert, 2008).
382	!-- The same procedure is applied to cloud droplets if they are determined
383	!-- prognostically.
384	DO k = nzb_s_inner(j,i)+1, nzt
385
386	IF ( qr(k,j,i) <= eps_sb ) THEN
387	qr(k,j,i) = 0.0
388	nr(k,j,i) = 0.0
389	ELSE
390	!
391	!-- Adjust number of raindrops to avoid nonlinear effects in
392	!-- sedimentation and evaporation of rain drops due to too small or
393	!-- too big weights of rain drops (Stevens and Seifert, 2008).
394	IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN
395	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
396	ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN
397	nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
398	ENDIF
399
400	ENDIF
401
402	ENDDO
403
404	END SUBROUTINE adjust_cloud_ij
405
406
407	SUBROUTINE autoconversion_ij( i, j )
408
409	USE arrays_3d
410	USE cloud_parameters
411	USE control_parameters
412	USE grid_variables
413	USE indices
414
415	IMPLICIT NONE
416
417	INTEGER :: i, j, k
418	REAL :: alpha_cc, autocon, epsilon, k_au, l_mix, nu_c, phi_au, &
419	r_cc, rc, re_lambda, selfcoll, sigma_cc, tau_cloud, xc
420
421
422	k_au = k_cc / ( 20.0 * x0 )
423
424	DO k = nzb_s_inner(j,i)+1, nzt
425
426	IF ( qc_1d(k) > eps_sb ) THEN
427	!
428	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
429	!-- (1.0 - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
430	tau_cloud = 1.0 - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
431	!
432	!-- Universal function for autoconversion process
433	!-- (Seifert and Beheng, 2006):
434	phi_au = 600.0 * tau_cloud*0.68 ( 1.0 - tau_cloud0.68 )3
435	!
436	!-- Shape parameter of gamma distribution (Geoffroy et al., 2010):
437	!-- (Use constant nu_c = 1.0 instead?)
438	nu_c = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * qc(k,j,i) - 0.28 )
439	!
440	!-- Mean weight of cloud droplets:
441	xc = hyrho(k) * qc_1d(k) / nc_1d(k)
442	!
443	!-- Parameterized turbulence effects on autoconversion (Seifert,
444	!-- Nuijens and Stevens, 2010)
445	IF ( turbulence ) THEN
446	!
447	!-- Weight averaged radius of cloud droplets:
448	rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0 )
449
450	alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
451	r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0 + b_3 * nu_c )
452	sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0 + c_3 * nu_c )
453	!
454	!-- Mixing length (neglecting distance to ground and stratification)
455	l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0 )
456	!
457	!-- Limit dissipation rate according to Seifert, Nuijens and
458	!-- Stevens (2010)
459	epsilon = MIN( 0.06, diss(k,j,i) )
460	!
461	!-- Compute Taylor-microscale Reynolds number:
462	re_lambda = 6.0 / 11.0 * ( l_mix / c_const )*( 2.0 / 3.0 ) &
463	SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0 )
464	!
465	!-- The factor of 1.0E4 is needed to convert the dissipation rate
466	!-- from m2 s-3 to cm2 s-3.
467	k_au = k_au * ( 1.0 + &
468	epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )*0.25 &
469	( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) / &
470	sigma_cc )**2 ) + beta_cc ) )
471	ENDIF
472	!
473	!-- Autoconversion rate (Seifert and Beheng, 2006):
474	autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) / &
475	( nu_c + 1.0 )*2 qc_1d(k)*2 xc*2 &
476	( 1.0 + phi_au / ( 1.0 - tau_cloud )*2 ) &
477	rho_surface
478	autocon = MIN( autocon, qc_1d(k) / dt_micro )
479
480	qr_1d(k) = qr_1d(k) + autocon * dt_micro
481	qc_1d(k) = qc_1d(k) - autocon * dt_micro
482	nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
483
484	ENDIF
485
486	ENDDO
487
488	END SUBROUTINE autoconversion_ij
489
490
491	SUBROUTINE accretion_ij( i, j )
492
493	USE arrays_3d
494	USE cloud_parameters
495	USE control_parameters
496	USE indices
497
498	IMPLICIT NONE
499
500	INTEGER :: i, j, k
501	REAL :: accr, k_cr, phi_ac, tau_cloud, xc
502
503	DO k = nzb_s_inner(j,i)+1, nzt
504	IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN
505	!
506	!-- Intern time scale of coagulation (Seifert and Beheng, 2006):
507	tau_cloud = 1.0 - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
508	!
509	!-- Universal function for accretion process
510	!-- (Seifert and Beheng, 2001):
511	phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
512	phi_ac = ( phi_ac2 )2
513	!
514	!-- Parameterized turbulence effects on autoconversion (Seifert,
515	!-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
516	!-- convert the dissipation (diss) from m2 s-3 to cm2 s-3.
517	IF ( turbulence ) THEN
518	k_cr = k_cr0 * ( 1.0 + 0.05 * &
519	MIN( 600.0, diss(k,j,i) * 1.0E4 )**0.25 )
520	ELSE
521	k_cr = k_cr0
522	ENDIF
523	!
524	!-- Accretion rate (Seifert and Beheng, 2006):
525	accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * &
526	SQRT( rho_surface * hyrho(k) )
527	accr = MIN( accr, qc_1d(k) / dt_micro )
528
529	qr_1d(k) = qr_1d(k) + accr * dt_micro
530	qc_1d(k) = qc_1d(k) - accr * dt_micro
531
532	ENDIF
533
534	ENDDO
535
536	END SUBROUTINE accretion_ij
537
538
539	SUBROUTINE selfcollection_breakup_ij( i, j )
540
541	USE arrays_3d
542	USE cloud_parameters
543	USE control_parameters
544	USE indices
545
546	IMPLICIT NONE
547
548	INTEGER :: i, j, k
549	REAL :: breakup, dr, phi_br, selfcoll
550
551	DO k = nzb_s_inner(j,i)+1, nzt
552	IF ( qr_1d(k) > eps_sb ) THEN
553	!
554	!-- Selfcollection rate (Seifert and Beheng, 2001):
555	selfcoll = k_rr * nr_1d(k) * qr_1d(k) * &
556	SQRT( hyrho(k) * rho_surface )
557	!
558	!-- Weight averaged diameter of rain drops:
559	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 )
560	!
561	!-- Collisional breakup rate (Seifert, 2008):
562	IF ( dr >= 0.3E-3 ) THEN
563	phi_br = k_br * ( dr - 1.1E-3 )
564	breakup = selfcoll * ( phi_br + 1.0 )
565	ELSE
566	breakup = 0.0
567	ENDIF
568
569	selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
570	nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
571
572	ENDIF
573	ENDDO
574
575	END SUBROUTINE selfcollection_breakup_ij
576
577
578	SUBROUTINE evaporation_rain_ij( i, j )
579	!
580	!-- Evaporation of precipitable water. Condensation is neglected for
581	!-- precipitable water.
582
583	USE arrays_3d
584	USE cloud_parameters
585	USE constants
586	USE control_parameters
587	USE indices
588
589	IMPLICIT NONE
590
591	INTEGER :: i, j, k
592	REAL :: alpha, dr, e_s, evap, evap_nr, f_vent, g_evap, lambda_r, &
593	mu_r, mu_r_2, mu_r_5d2, nr_0, q_s, sat, t_l, temp, xr
594
595	DO k = nzb_s_inner(j,i)+1, nzt
596	IF ( qr_1d(k) > eps_sb ) THEN
597	!
598	!-- Actual liquid water temperature:
599	t_l = t_d_pt(k) * pt_1d(k)
600	!
601	!-- Saturation vapor pressure at t_l:
602	e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86 ) )
603	!
604	!-- Computation of saturation humidity:
605	q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
606	alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )
607	q_s = q_s * ( 1.0 + alpha * q_1d(k) ) / ( 1.0 + alpha * q_s )
608	!
609	!-- Supersaturation:
610	sat = MIN( 0.0, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 )
611	!
612	!-- Actual temperature:
613	temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
614
615	g_evap = 1.0 / ( ( l_v / ( r_v * temp ) - 1.0 ) * l_v / &
616	( thermal_conductivity_l * temp ) + r_v * temp / &
617	( diff_coeff_l * e_s ) )
618	!
619	!-- Mean weight of rain drops
620	xr = hyrho(k) * qr_1d(k) / nr_1d(k)
621	!
622	!-- Weight averaged diameter of rain drops:
623	dr = ( xr * dpirho_l )**( 1.0 / 3.0 )
624	!
625	!-- Compute ventilation factor and intercept parameter
626	!-- (Seifert and Beheng, 2006; Seifert, 2008):
627	IF ( ventilation_effect ) THEN
628	!
629	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
630	!-- Stevens and Seifert, 2008):
631	mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) )
632	!
633	!-- Slope parameter of gamma distribution (Seifert, 2008):
634	lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) * &
635	( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr
636
637	mu_r_2 = mu_r + 2.0
638	mu_r_5d2 = mu_r + 2.5
639	f_vent = a_vent * gamm( mu_r_2 ) * &
640	lambda_r**( -mu_r_2 ) + &
641	b_vent * schmidt_p_1d3 * &
642	SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * &
643	lambda_r*( -mu_r_5d2 ) &
644	( 1.0 - 0.5 * ( b_term / a_term ) * &
645	( lambda_r / &
646	( c_term + lambda_r ) )**mu_r_5d2 - &
647	0.125 * ( b_term / a_term )*2 &
648	( lambda_r / &
649	( 2.0 * c_term + lambda_r ) )**mu_r_5d2 - &
650	0.0625 * ( b_term / a_term )*3 &
651	( lambda_r / &
652	( 3.0 * c_term + lambda_r ) )**mu_r_5d2 - &
653	0.0390625 * ( b_term / a_term )*4 &
654	( lambda_r / &
655	( 4.0 * c_term + lambda_r ) )**mu_r_5d2 )
656	nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0 ) / &
657	gamm( mu_r + 1.0 )
658	ELSE
659	f_vent = 1.0
660	nr_0 = nr_1d(k) * dr
661	ENDIF
662	!
663	!-- Evaporation rate of rain water content (Seifert and Beheng, 2006):
664	evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / &
665	hyrho(k)
666
667	evap = MAX( evap, -qr_1d(k) / dt_micro )
668	evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
669	-nr_1d(k) / dt_micro )
670
671	qr_1d(k) = qr_1d(k) + evap * dt_micro
672	nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
673	ENDIF
674
675	ENDDO
676
677	END SUBROUTINE evaporation_rain_ij
678
679
680	SUBROUTINE sedimentation_cloud_ij( i, j )
681
682	USE arrays_3d
683	USE cloud_parameters
684	USE constants
685	USE control_parameters
686	USE indices
687
688	IMPLICIT NONE
689
690	INTEGER :: i, j, k
691	REAL :: sed_qc_const
692
693	REAL, DIMENSION(nzb:nzt+1) :: sed_qc
694
695	!
696	!-- Sedimentation of cloud droplets (Heus et al., 2010):
697	sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))*( 2.0 / 3.0 ) &
698	EXP( 5.0 * LOG( sigma_gc )**2 )
699
700	sed_qc(nzt+1) = 0.0
701
702	DO k = nzt, nzb_s_inner(j,i)+1, -1
703	IF ( qc_1d(k) > eps_sb ) THEN
704	sed_qc(k) = sed_qc_const * nc_1d(k)*( -2.0 / 3.0 ) &
705	( qc_1d(k) * hyrho(k) )**( 5.0 / 3.0 )
706	ELSE
707	sed_qc(k) = 0.0
708	ENDIF
709
710	sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / &
711	dt_micro + sed_qc(k+1) )
712
713	q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
714	hyrho(k) * dt_micro
715	qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
716	hyrho(k) * dt_micro
717	pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
718	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
719
720	ENDDO
721
722	END SUBROUTINE sedimentation_cloud_ij
723
724
725	SUBROUTINE sedimentation_rain_ij( i, j )
726
727	USE arrays_3d
728	USE cloud_parameters
729	USE constants
730	USE control_parameters
731	USE indices
732	USE statistics
733
734	IMPLICIT NONE
735
736	INTEGER :: i, j, k, k_run
737	REAL :: c_run, d_max, d_mean, d_min, dr, dt_sedi, flux, lambda_r, &
738	mu_r, z_run
739
740	REAL, DIMENSION(nzb:nzt+1) :: c_nr, c_qr, d_nr, d_qr, nr_slope, &
741	qr_slope, sed_nr, sed_qr, w_nr, w_qr
742	!
743	!-- Computation of sedimentation flux. Implementation according to Stevens
744	!-- and Seifert (2008).
745	IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0
746	!
747	!-- Compute velocities
748	DO k = nzb_s_inner(j,i)+1, nzt
749	IF ( qr_1d(k) > eps_sb ) THEN
750	!
751	!-- Weight averaged diameter of rain drops:
752	dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 )
753	!
754	!-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
755	!-- Stevens and Seifert, 2008):
756	mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) )
757	!
758	!-- Slope parameter of gamma distribution (Seifert, 2008):
759	lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) * &
760	( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr
761
762	w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + &
763	c_term / lambda_r )*( -1.0 ( mu_r + 1.0 ) ) ) )
764	w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + &
765	c_term / lambda_r )*( -1.0 ( mu_r + 4.0 ) ) ) )
766	ELSE
767	w_nr(k) = 0.0
768	w_qr(k) = 0.0
769	ENDIF
770	ENDDO
771	!
772	!-- Adjust boundary values
773	w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
774	w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
775	w_nr(nzt+1) = 0.0
776	w_qr(nzt+1) = 0.0
777	!
778	!-- Compute Courant number
779	DO k = nzb_s_inner(j,i)+1, nzt
780	c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * &
781	dt_micro * ddzu(k)
782	c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * &
783	dt_micro * ddzu(k)
784	ENDDO
785	!
786	!-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
787	IF ( limiter_sedimentation ) THEN
788
789	DO k = nzb_s_inner(j,i)+1, nzt
790	d_mean = 0.5 * ( qr_1d(k+1) + qr_1d(k-1) )
791	d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
792	d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)
793
794	qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
795	ABS( d_mean ) )
796
797	d_mean = 0.5 * ( nr_1d(k+1) + nr_1d(k-1) )
798	d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
799	d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)
800
801	nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
802	ABS( d_mean ) )
803	ENDDO
804
805	ELSE
806
807	nr_slope = 0.0
808	qr_slope = 0.0
809
810	ENDIF
811
812	sed_nr(nzt+1) = 0.0
813	sed_qr(nzt+1) = 0.0
814	!
815	!-- Compute sedimentation flux
816	DO k = nzt, nzb_s_inner(j,i)+1, -1
817	!
818	!-- Sum up all rain drop number densities which contribute to the flux
819	!-- through k-1/2
820	flux = 0.0
821	z_run = 0.0 ! height above z(k)
822	k_run = k
823	c_run = MIN( 1.0, c_nr(k) )
824	DO WHILE ( c_run > 0.0 .AND. k_run <= nzt )
825	flux = flux + hyrho(k_run) * &
826	( nr_1d(k_run) + nr_slope(k_run) * ( 1.0 - c_run ) * &
827	0.5 ) * c_run * dzu(k_run)
828	z_run = z_run + dzu(k_run)
829	k_run = k_run + 1
830	c_run = MIN( 1.0, c_nr(k_run) - z_run * ddzu(k_run) )
831	ENDDO
832	!
833	!-- It is not allowed to sediment more rain drop number density than
834	!-- available
835	flux = MIN( flux, &
836	hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
837
838	sed_nr(k) = flux / dt_micro
839	nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &
840	hyrho(k) * dt_micro
841	!
842	!-- Sum up all rain water content which contributes to the flux
843	!-- through k-1/2
844	flux = 0.0
845	z_run = 0.0 ! height above z(k)
846	k_run = k
847	c_run = MIN( 1.0, c_qr(k) )
848
849	DO WHILE ( c_run > 0.0 .AND. k_run <= nzt-1 )
850
851	flux = flux + hyrho(k_run) * &
852	( qr_1d(k_run) + qr_slope(k_run) * ( 1.0 - c_run ) * &
853	0.5 ) * c_run * dzu(k_run)
854	z_run = z_run + dzu(k_run)
855	k_run = k_run + 1
856	c_run = MIN( 1.0, c_qr(k_run) - z_run * ddzu(k_run) )
857
858	ENDDO
859	!
860	!-- It is not allowed to sediment more rain water content than available
861	flux = MIN( flux, &
862	hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
863
864	sed_qr(k) = flux / dt_micro
865
866	qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
867	hyrho(k) * dt_micro
868	q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
869	hyrho(k) * dt_micro
870	pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
871	hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
872	!
873	!-- Compute the rain rate
874	prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * &
875	weight_substep(intermediate_timestep_count)
876	ENDDO
877
878	!
879	!-- Precipitation amount
880	IF ( intermediate_timestep_count == intermediate_timestep_count_max &
881	.AND. ( dt_do2d_xy - time_do2d_xy ) < &
882	precipitation_amount_interval ) THEN
883
884	precipitation_amount(j,i) = precipitation_amount(j,i) + &
885	prr(nzb_s_inner(j,i)+1,j,i) * &
886	hyrho(nzb_s_inner(j,i)+1) * dt_3d
887	ENDIF
888
889	END SUBROUTINE sedimentation_rain_ij
890
891
892	!
893	!-- This function computes the gamma function (Press et al., 1992).
894	!-- The gamma function is needed for the calculation of the evaporation
895	!-- of rain drops.
896	FUNCTION gamm( xx )
897
898	USE cloud_parameters
899
900	IMPLICIT NONE
901
902	REAL :: gamm, ser, tmp, x_gamm, xx, y_gamm
903	INTEGER :: j
904
905
906	x_gamm = xx
907	y_gamm = x_gamm
908	tmp = x_gamm + 5.5
909	tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp
910	ser = 1.000000000190015
911
912	DO j = 1, 6
913	y_gamm = y_gamm + 1.0
914	ser = ser + cof( j ) / y_gamm
915	ENDDO
916
917	!
918	!-- Until this point the algorithm computes the logarithm of the gamma
919	!-- function. Hence, the exponential function is used.
920	! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
921	gamm = EXP( tmp ) * stp * ser / x_gamm
922
923	RETURN
924
925	END FUNCTION gamm
926
927	END MODULE microphysics_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |