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

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

optimization of two-moments cloud physics

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