Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |