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

Last change on this file since 1066 was 1066, checked in by hoffmann, 9 years ago

last commit documented

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