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

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

last commit documented

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