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

Last change on this file since 1092 was 1092, checked in by raasch, 9 years ago

unused variables remove from several routines

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