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

Last change on this file since 1114 was 1107, checked in by raasch, 11 years ago

last commit documented

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