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

Last change on this file since 1207 was 1116, checked in by hoffmann, 12 years ago

last commit documented

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